text:
module MakeGlassBox

using HDF5
using StaticArrays
using Statistics
using Random
using LinearAlgebra
using SpecialFunctions
using Printf
using Parameters
#using .Threads
using PyPlot

using OctreeBH

export generate_ICfile_glass_box, plot_convergence
export Params

function vec2svec(vec::Matrix{T}) where {T}
    ndim = size(vec,1)
    svec = [SVector{ndim,T}(vec[:,i]) for i in 1:size(vec,2)]
end
function mat2smat(mat::Array{T,3}) where {T}
    ndim = size(mat,1)
    smat = [SMatrix{ndim,ndim,T}(mat[:,:,i]) for i in 1:size(mat,3)]
end
const KERNELCONST = (16.0/pi)
@inline ramp(x) = max(0, x);
@inline function kernel_cubic(x::T) where {T}
    return (KERNELCONST * (ramp(1.0 - x)^3 - 4.0 * ramp(0.5 - x)^3))
end

const SOLARMASS = 1.989e+43
const PROTONMASS = 1.6726e-24
const XH = 0.71
const BOLTZMANN = 1.38e-16
const GAMMA = 5.0 / 3

@with_kw struct Params{NDIM,T}

    boxsize::T

    nH::T = 0.
    T_mu::T = 0.

    UnitLength_in_cm::T = 3.08568e+21
    UnitMass_in_g::T = 1.989e+43
    UnitVelocity_in_cm_per_s::T = 1e5
    u0::T = T_mu * BOLTZMANN / (GAMMA-1.0) / PROTONMASS / UnitVelocity_in_cm_per_s^2
    Mgas_tot::T = nH * PROTONMASS * (boxsize*UnitLength_in_cm)^3 / XH / UnitMass_in_g

    seed::Int64 = 0

    #Ngas is always known s.t. we can calculate h0 & Lint below
    #OTOH, mgas can be zero if Ngas is user input
    mgas::T = 0.0
    Ngas::Int64 = (mgas > 0.) ? round(Int64, Mgas_tot / mgas) : Ngas

    filename::String = ""

    ms::T = 1.0
    flag_plot_part::Bool = false
    flag_plot_conv::Bool = false

    #NDIM::Int64 = 3
    geo_fac::T = pi^(0.5*NDIM)/gamma(0.5*NDIM+1)

    Nngb::T = (100.0)^(NDIM/3.0)

    unit_boxsize = 1.0
    Lint::T = unit_boxsize / Ngas^(1.0/NDIM) #inter-particle spacing
    h0::T = Lint * (Nngb / geo_fac)^(1.0/NDIM)

    Nstep::Int64 = 30
    Njump::Int64 = 3
    Nidx::Int64 = NDIM-1
    A0::T = 1.0
    alpha::T = 0.1 / sqrt(Lint)
end

function calc_force!(force,Nngb_true,Nngb_weight,X,V,tree,par::Params{NDIM,T}) where {NDIM,T}

    @unpack boxsize, mgas, geo_fac, Lint, h0, Nidx, A0, alpha = par
    #NDIM=3
    unit_boxsize = 1.0
    boxsizes = SVector{NDIM}(ones(NDIM)) .* unit_boxsize

    Threads.@threads for i in eachindex(X)
        force[i] = SVector{NDIM}(zeros(NDIM))
        #idx_ngbs = get_scatter_ngb_tree(X[i], tree, boxsizes)
        idx_ngbs = get_gather_ngb_tree(X[i], h0, tree, boxsizes)
        Nngb_weight[i] = 0.0
        for k in eachindex(idx_ngbs)
            j = idx_ngbs[k]
            if i==j continue end
            dx = nearest.(X[i] - X[j], boxsizes) ./ Lint #from j to i, positive = respulsive!
            dr = norm(dx)
            force[i] += (A0/dr^Nidx .* (dx./dr) - alpha .* V[i])
            Wij = kernel_cubic(dr*Lint/h0)
            Nngb_weight[i] += Wij
        end
        Nngb_weight[i] *= geo_fac
        Nngb_true[i] = length(idx_ngbs)
    end
end

function make_glass_box(X, V, par::Params{NDIM,T}) where {NDIM,T}

    @unpack boxsize, Ngas, h0, Lint, ms, flag_plot_part, flag_plot_conv, geo_fac, Nngb, Nstep, Njump, Nidx = par
    #NDIM=3
    mass = mass_H2 = mass_CO = zeros(Ngas)

    unit_boxsize = T(1)
    boxsizes = SVector{NDIM,T}(ones(NDIM)) .* unit_boxsize #use boxsize=1 here!!! otherwise double-counting!!!!!!
    center = boxsizes .* T(0.5)

    hsml = ones(Ngas) .* h0

    part = [Data{NDIM,T}(X[i], i, hsml[i], mass[i]) for i in 1:Ngas]
    tree = buildtree(part, center, boxsizes);

    Nngb_true = zeros(Ngas)
    Nngb_weight = zeros(Ngas)
    force = [SVector{NDIM}(zeros(NDIM)) for _ in eachindex(X)]
    fac_dt = T(0.3)

    std_dNngb_true = zeros(Nstep+1)
    std_dNngb_weight = zeros(Nstep+1)
    calc_force!(force,Nngb_true,Nngb_weight,X,V,tree,par)

    std_dNngb_true[1] = std(Nngb_true .- Nngb) / Nngb
    std_dNngb_weight[1] = std(Nngb_weight .- Nngb) / Nngb

    dt = fac_dt * minimum( @. sqrt(Lint / norm(force)) )

    @. V = V + 0.5*dt * force
    for i in 1:Nstep
        println("istep = ", i, "  Nstep = ", Nstep)
        if flag_plot_part && i%Njump==0
            fig, ax = subplots(1, 1, figsize=(6, 6))
            ax.plot(getindex.(X,1), getindex.(X,2), ".", ms=ms)
            ax.set_aspect("equal")
            ax.axis([0,1,0,1])
            fig.tight_layout()
            fig.savefig("ic_relax_"*string(div(i,Njump))*".png")
            clf()
        end
        @. X = X + dt * V
        for j in eachindex(X)
            X[j] -= convert.(T, X[j] .> boxsizes) .* boxsizes
            X[j] += convert.(T, X[j] .< 0.0) .* boxsizes
        end

        part = [Data{NDIM,T}(X[i], i, hsml[i], mass[i]) for i in 1:Ngas]
        tree = buildtree(part, center, boxsizes);

        calc_force!(force,Nngb_true,Nngb_weight,X,V,tree,par)
        @. V = V + dt * force
        dt = fac_dt * minimum( @. sqrt(Lint / norm(force)) )

        std_dNngb_true[i+1] = std(Nngb_true .- Nngb) / Nngb
        std_dNngb_weight[i+1] = std(Nngb_weight .- Nngb) / Nngb
    end
    #return Nngb_true, std_dNngb_true, Nngb_weight, std_dNngb_weight
    if flag_plot_conv == true
        plot_convergence(std_dNngb_true, std_dNngb_weight, Nngb_true, Nngb_weight, Nngb)
    end
end

function generate_ICfile_glass_box(par::Params{NDIM,T}) where{NDIM,T}

    @unpack nH, T_mu, boxsize, Ngas, mgas, Mgas_tot, u0, filename, seed = par

    mgas = Mgas_tot / Ngas

    @show mgas, Ngas, boxsize
    id_gas = collect(Int32, 1:Ngas)
    mass = ones(Ngas) .* mgas
    u = ones(Ngas) .* u0;
    #hsml = ones(Ngas) .* h0

    ########## prepare glass configuration ##########
    unit_boxsize = 1.0
    #pos = rand(MersenneTwister(seed), NDIM, Ngas) .* unit_boxsize;
    pos = zeros(NDIM, Ngas);
    vel = zeros(NDIM, Ngas);

    Random.seed!(seed);
    X = [SVector{NDIM,T}(rand(NDIM)) for _ in 1:Ngas] .* unit_boxsize
    V = [SVector{NDIM,T}(zeros(NDIM)) for _ in 1:Ngas]
    #Nngb_true, std_dNngb_true, Nngb_weight, std_dNngb_weight =
    @time make_glass_box(X,V,par);

    #save back to pos and scale to the actual boxsize
    for i in eachindex(X)
        pos[:,i] = X[i] .* boxsize;  #do this only to positions; velocities remain zero
    end

    if filename == ""
        filename = "ics"
        if nH > 0.
            filename *= "_nH" * @sprintf("%.0e", nH) * "_T" * @sprintf("%.0e", T_mu)
        else
            filename *= "_Mtot" * @sprintf("%.0e", Mgas_tot) * "_u" * @sprintf("%.0e", u0)
        end
        filename *= "_boxsize" * @sprintf("%.2e", boxsize) *
                    "_Ngas" * @sprintf("%.0e", Ngas) *
                    "_mgas" * @sprintf("%.1e", mgas) * ".hdf5"
    end

    ########## write to file ##########
    save_gadget_ics(filename, pos, vel, id_gas, mass, u, Ngas, boxsize)

    return pos, vel, id_gas, mass, u, Ngas, boxsize
end

function save_gadget_ics(filename, pos, vel, id_gas, mass, u, Ngas, boxsize)
    T = Float32

    fid=h5open(filename,"w")

    grp_head = create_group(fid,"Header");
    attributes(fid["Header"])["NumPart_ThisFile"]       = Int32[Ngas, 0, 0, 0, 0, 0]
    attributes(fid["Header"])["MassTable"]              = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    attributes(fid["Header"])["Time"]                   = 0.0
    attributes(fid["Header"])["Redshift"]               = 0.0
    attributes(fid["Header"])["Flag_Sfr"]               = 0
    attributes(fid["Header"])["Flag_Feedback"]          = 0
    attributes(fid["Header"])["NumPart_Total"]          = Int32[Ngas, 0, 0, 0, 0, 0]
    attributes(fid["Header"])["Flag_Cooling"]           = 0
    attributes(fid["Header"])["NumFilesPerSnapshot"]    = 1
    attributes(fid["Header"])["BoxSize"]                = boxsize
    attributes(fid["Header"])["Omega0"]                 = 0.27
    attributes(fid["Header"])["OmegaLambda"]            = 0.73
    attributes(fid["Header"])["HubbleParam"]            = 1.0
    attributes(fid["Header"])["Flag_StellarAge"]        = 0
    attributes(fid["Header"])["Flag_Metals"]            = 0
    attributes(fid["Header"])["NumPart_Total_HighWord"] = UInt32[0,0,0,0,0,0]
    attributes(fid["Header"])["flag_entropy_instead_u"] = 0
    attributes(fid["Header"])["Flag_DoublePrecision"]   = 0
    attributes(fid["Header"])["Flag_IC_Info"]           = 0
    #attributes(fid["Header"])["lpt_scalingfactor"] =

    grp_part = create_group(fid,"PartType0");
    h5write(filename, "PartType0/Coordinates"   , T.(pos))
    h5write(filename, "PartType0/Velocities"    , T.(vel))
    h5write(filename, "PartType0/ParticleIDs"   , id_gas)
    h5write(filename, "PartType0/Masses"        , T.(mass))
    h5write(filename, "PartType0/InternalEnergy", T.(u))

    close(fid)
end

function plot_convergence(std_dNngb_true::Vector{T}, std_dNngb_weight::Vector{T}, Nngb_true::Vector{T}, Nngb_weight::Vector{T}, Nngb::T) where {T}
    ##### convergence test
    fig = figure("convergence test",figsize=(9,4))
    subplot(121)
    plot(collect(Int32,1:length(std_dNngb_true)), std_dNngb_true, "-", label="true")
    plot(collect(Int32,1:length(std_dNngb_weight)), std_dNngb_weight, "-", label="weighted")
    ylabel("std(dNngb / Nngb)")
    ylim(0,0.22)
    legend(loc="best", frameon=false)

    subplot(122)
    hist(mean( (Nngb ./ Nngb_true ) )   .* Nngb_true  ./Nngb, bins=30, label="true"    , histtype="step", normed=true, range=[0.85,1.15])
    hist(mean( (Nngb ./ Nngb_weight ) ) .* Nngb_weight./Nngb, bins=30, label="weighted", histtype="step", normed=true, range=[0.85,1.15])
    ylabel("Nngb_true / Nngb (final)")
    legend(loc="best", frameon=false)

    tight_layout()
    Ngas = length(Nngb_true)
    savefig("dNngb_converge_"*string(Ngas)*".png")
    0
end

end #module MakeGlassBox


text:
export URI, URIAuthority, @uri_str

const uriregex = r"^(([^:/?#]+):)?(//([^/?#]*))?([^?#]*)(\?([^#]*))?(#(.*))?"
const authorityregex = r"^(([^/?#\[\]]*)@)?(([^:/?#\[\]@]+)|(\[([^/?#@\[\]]+)\]))(:([0-9]+))?$"
const ipv4regex = r"^[0-9]{1,3}\.[0-9]{1,3}\.[0-9]{1,3}\.[0-9]{1,3}$"
struct URIAuthority
    userinfo::Union{AbstractString,Nothing}
    host::Union{AbstractString,IPv4,IPv6,Nothing}
    port::Union{UInt,Nothing}
end # struct URIAuthority

struct URI
    scheme::Union{AbstractString,Nothing}
    authority::Union{URIAuthority,Nothing}
    path::Union{AbstractString,Nothing}
    query::Union{AbstractString,Nothing}
    fragment::Union{AbstractString,Nothing}
end # struct URI

function URIAuthority(authority::AbstractString)
    parts = match(authorityregex, authority)
    if parts == nothing
        return URIAuthority(nothing, nothing, nothing)
    end

    userinfo = parts[2]
    host = if parts[4] != nothing
        if occursin(ipv4regex, parts[4])
            try
                IPv4(parts[4])
            catch e
                parts[4]
            end
        else
            parts[4]
        end
    elseif parts[6] != nothing
        try
            IPv6(parts[6])
        catch e
            "[$(parts[4])]"
        end
    else
        nothing
    end

    port = if parts[8] == nothing
        nothing
    else
        parse(UInt, parts[8])
    end

    URIAuthority(userinfo, host, port)
end

function URI(uri::AbstractString)
    parts = match(uriregex, uri)
    scheme = parts[2]
    authority = if parts[4] == nothing nothing else URIAuthority(parts[4]) end
    path = parts[5]
    query = parts[7]
    fragment = parts[9]
    URI(scheme, authority, path, query, fragment)
end

macro uri_str(uri)
    URI(uri)
end

function Base.print(io::IO, authority::URIAuthority)
    if authority.userinfo != nothing
        print(io, authority.userinfo)
        print(io, '@')
    end

    if authority.host isa AbstractString || authority.host isa IPv4
        print(io, authority.host)
    elseif authority.host isa IPv6
        print(io, '[')
        print(io, authority.host)
        print(io, ']')
    end

    if authority.port != nothing
        print(io, ':')
        print(io, authority.port)
    end
end

function Base.print(io::IO, uri::URI)
    if uri.scheme != nothing
        print(io, uri.scheme)
        print(io, ':')
    end
    if uri.authority != nothing
        print(io, "//")
        print(io, uri.authority)
    end
    if uri.path != nothing
        print(io, uri.path)
    end
    if uri.query != nothing
        print(io, '?')
        print(io, uri.query)
    end
    if (uri.fragment != nothing)
        print(io, '#')
        print(io, uri.fragment)
    end
end

ext: jl


text:
using Test

#--- initial data (smooth, support compact). Zero velocity for compatibility between models.
init     = Init(x->exp.(-1/8*x.^4),x-> 0. * x )

#--- Tests on shallow water models

para = ( ϵ  = 0.01, μ = 0.01)  # physical parameters
paraX= ( N  = 2^6, L  = 4)   # mesh with 64 collocation points on [-4,4]
paraT= ( T  = 1e-1, dt = 1e-2) # timegrid with 10 instants: t=[0.0:1.0:10.0]/100
param = merge(para,paraX)  # used to construct models

models=AbstractModel[];precisions=Real[]

# Build shallow water models
push!(models,SerreGreenNaghdi(param;
            dealias = 0,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = true,
            restart	= 15,
            maxiter	= 25,
            label	= "Green-Naghdi with GMRES"
			) )
push!(precisions,para.μ^2)

push!(models,SerreGreenNaghdi(param::NamedTuple;
			dealias = 1,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "Green-Naghdi with LU"
			) )
push!(precisions,para.μ^2)

push!(models,WhithamGreenNaghdi(param::NamedTuple;
			SGN=true,
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "(Whitham-)Green-Naghdi with LU"
			) )
push!(precisions,para.μ^2)

push!(models,WhithamGreenNaghdi(param;
			SGN=true,
            dealias = false,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 15,
            maxiter	= 25,
            label	= "(Whitham-)Green-Naghdi with GMRES"
			) )
push!(precisions,para.μ^2)

push!(models,WhithamGreenNaghdi(param;
            dealias = 0,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 5,
            maxiter	= 15,
            label	= "Whitham-Green-Naghdi with GMRES") )
push!(precisions,para.μ^2*param.ϵ)

push!(models,WhithamGreenNaghdi(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "Whitham-Green-Naghdi with LU") )
push!(precisions,para.μ^2*param.ϵ)

push!(models,IsobeKakinuma(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "Isobe-Kakinuma with LU") )
push!(precisions,para.μ^3)

push!(models,IsobeKakinuma(param;
            dealias = 0,
            ktol	= 1e-14,
            iterate = true,
            gtol	= 1e-14,
            precond = true,
            restart	= 10,
            maxiter	= 20,
            label	= "Isobe-Kakinuma with GMRES") )
push!(precisions,para.μ^3)

push!(models,NonHydrostatic(param;
            dealias = 2,
            ktol	= 1e-14,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "non-hydrostatic with GMRES") )
push!(precisions,para.μ/5)

push!(models,NonHydrostatic(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "non-hydrostatic with LU") )
push!(precisions,para.μ/5)

push!(models,SquareRootDepth(param;
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 15,
            maxiter	= 20,
            label	= "square-root depth with GMRES") )
push!(precisions,max(para.μ*para.ϵ,para.μ^2))

push!(models,SquareRootDepth(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "square-root depth with LU"
			) )
push!(precisions,max(para.μ*para.ϵ,para.μ^2))

push!(models,SaintVenant(param;
                    dealias=0,ktol=1e-12,
                    label="Saint-Venant"
					) )
push!(precisions,para.μ)

push!(models,Boussinesq(param;
                    a=-1//2,b=5//12,
                    dealias=0,ktol=1e-12,
                    label="Boussinesq"
					) )
push!(precisions,max(para.μ*para.ϵ/4,para.μ^2))

push!(models,WhithamBoussinesq(param;
				Boussinesq=true,
				α = 1, a = -1/2, b = 5/12,
				dealias = 0,
				ktol	= 1e-14,
				label 	= "(Whitham-)Boussinesq"
				) )
push!(precisions,max(para.μ*para.ϵ/4,para.μ^2))

push!(models,WhithamBoussinesq(param;
				α = 1.5, a = -1, b = 3,
				dealias = 0,
				ktol	= 1e-12,
				label 	= "Whitham-Boussinesq"
				) )
push!(precisions,para.μ*para.ϵ/4)

# Build reference problem (water waves)
modelWW =  WaterWaves(merge(param,(ν=1,));
					IL	    = false,
					method  = 3,
					tol	    = 1e-8,
					maxiter = 5,
					dealias	= 0,
					ktol	= 1e-10,
					label	= "water waves (shallow water)",
					verbose	= false)
pbWW = Problem( modelWW, init, paraT )
solve!(pbWW;verbose=false)
ηWW,vWW,xWW=solution(pbWW)

for i in 1:length(models)
    # build the initial-value problem
    problem = Problem( models[i], init, paraT )
    # solve the initial-value problem
    solve!(problem ; verbose=false)
    # check the order of magnitude of the precision is correct (from above and from below)
    @testset "Initial value problem with the model $(models[i].label)" begin
		@test isapprox(solution(problem,x=xWW)[1] , ηWW, rtol = precisions[i]*paraT.T )
		@test !isapprox(solution(problem,x=xWW)[1] , ηWW, rtol = precisions[i]*paraT.T/10 )
    end
end

#--- Tests on deep water models
models=AbstractModel[];precisions=Real[]

para = ( ϵ  = 0.1, μ = 100)  # physical parameters
paraX= ( N  = 2^6, L  = 4)   # mesh with 64 collocation points on [-4,4]
paraT= ( T  = 1e-1, dt = 1e-2) # timegrid with 10 instants: t=[0.0:1.0:10.0]/100
param = merge(para,paraX)  # used to construct models

# Build deep layer models
push!(models, modifiedMatsuno(merge(param,(ν=1/√para.μ,));
			ktol=01e-14,dealias=0,
			label="modified Matsuno"
			) )
push!(precisions, para.ϵ.^2*para.μ)

push!(models,Airy(param;
                    label="Airy"
					) )
push!(precisions,para.ϵ*√para.μ)

push!(models, WWn(param;
			n		= 1,
			δ		= 0.1,
			m		= (-1,2),
			ktol	= 1e-12,
			dealias	= 1,
			label	= "WW1"
			) )
push!(precisions, para.ϵ*√para.μ)

push!(models, WWn(merge(param,(ν=1/√para.μ,));
			n		= 7,
			δ		= 0.1,
			m		= -Inf,
			ktol	= 1e-12,
			dealias	= 0,
			label	= "WW2"
			) )
push!(precisions, para.ϵ.^2*para.μ)

push!(models, WWn(param;
			n		= 3,
			δ		= 0.01,
			ktol	= 1e-14,
			dealias	= 0,
			label	= "WW3"
			) )
push!(precisions, para.ϵ.^3*para.μ^(3/2))

# Build reference problem (water waves)
modelWW = WaterWaves(param;
					IL	    = false,
					method  = 2,
					tol	    = 1e-12,
					maxiter = 25,
					dealias	= 0,
					ktol	= 1e-10,
					label	= "water waves (deep water)",
					verbose	= false)
pbWW = Problem( modelWW, init, paraT )
solve!(pbWW;verbose=false)
ηWW,vWW,xWW=solution(pbWW)

for i in 1:length(models)
    # build the initial-value problem
    problem = Problem( models[i], init, paraT )
    # solve the initial-value problem
    solve!(problem ; verbose=false)
    # check the order of magnitude of the precision is correct (from above and from below)
    @testset "Initial value problem with the model $(models[i].label)" begin
		@test isapprox(solution(problem,x=xWW)[1] , ηWW, rtol = precisions[i]*paraT.T/50 )
		@test !isapprox(solution(problem,x=xWW)[1] , ηWW, rtol = precisions[i]*paraT.T/500 )

    end
end

#--- Tests on infinite layer models
models=AbstractModel[];precisions=Real[]

para = ( ϵ  = 0.1, μ = Inf)  # physical parameters
paraX= ( N  = 2^6, L  = 4)   # mesh with 64 collocation points on [-4,4]
paraT= ( T  = 1e-1, dt = 1e-2) # timegrid with 10 instants: t=[0.0:1.0:10.0]/100
param = merge(para,paraX)  # used to construct models

# Build infinite layer models
push!(models, modifiedMatsuno(param;
			ktol=01e-14,dealias=0,
			label="modified Matsuno"
			) )
push!(precisions, para.ϵ.^2)

push!(models, WWn(param;
			n		= 1,
			δ		= 0.1,
			m		= (-1,2),
			ktol	= 1e-12,
			dealias	= 1,
			label	= "WW1") )
push!(precisions, para.ϵ)

push!(models, WWn(merge(param,(ν=0,));
			n		= 7,
			δ		= 0.1,
			m		= -Inf,
			ktol	= 1e-12,
			dealias	= 0,
			label	= "WW2") )
push!(precisions, para.ϵ.^2)

push!(models, WWn(param;
			IL		= true,
			n		= 3,
			δ		= 0.01,
			ktol	= 0,
			dealias	= 0,
			label	= "WW3") )
push!(precisions, para.ϵ.^3)

push!(models,AkersNicholls_fast( param;
                    dealias=true, label="fast deep quadratic") )
push!(precisions, para.ϵ.^2)

push!(models,AkersNicholls( param;
                    dealias=false, label="naive deep quadratic") )
push!(precisions, para.ϵ.^2)

push!(models,Matsuno_fast( param;
                    dealias=true, label="fast Matsuno") )
push!(precisions, para.ϵ.^2)

push!(models,Matsuno( param;
                    dealias=false, label="naive Matsuno") )
push!(precisions, para.ϵ.^2)

# Build reference problem (water waves)
modelWW =  WaterWaves(param;
					IL	    = true,
					method  = 1,
					tol	    = 1e-12,
					maxiter = 5,
					dealias	= 0,
					ktol	= 1e-14,
					label	= "water waves (infinite layer)",
					verbose	= false)
pbWW = Problem( modelWW, init, paraT )
solve!(pbWW;verbose=false)
ηWW,vWW,xWW=solution(pbWW)

for i in 1:length(models)
    # build the initial-value problem
    problem = Problem( models[i], init, paraT )
    # solve the initial-value problem
    solve!(problem ; verbose=false)
    # check the order of magnitude of the precision is correct (from above and from below)
    @testset "Initial value problem with the model $(models[i].label)" begin
		@test isapprox(solution(problem,x=xWW)[1] , ηWW, rtol = precisions[i]*paraT.T/50 )
		@test !isapprox(solution(problem,x=xWW)[1] , ηWW, rtol = precisions[i]*paraT.T/500 )

    end
end


text:
using Test
using OrdinaryDiffEq

import ClimaCore:
    ClimaCore, Domains, Topologies, Meshes, Geometry, Fields, Spaces
using LinearAlgebra, IntervalSets

@testset "DiffEq Solvers" begin
    domain = Domains.RectangleDomain(
        Geometry.XPoint(-2π) .. Geometry.XPoint(2π),
        Geometry.YPoint(-2π) .. Geometry.YPoint(2π),
        x1periodic = true,
        x2periodic = true,
    )

    n1, n2 = 3, 4
    Nq = 4
    mesh = Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = Topologies.Topology2D(mesh)
    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    # y = [cos(x1+x2+t), sin(x1+x2+t)]
    # dy/dt = [-y[2], y[1]]

    f(x, t) = (c = cos(x.x + x.y + t), s = sin(x.x + x.y + t))
    function dfdt(y, _, t)
        broadcast(y -> (c = -y.s, s = y.c), y)
    end
    function dfdt!(dydt, y, _, t)
        broadcast!(y -> (c = -y.s, s = y.c), dydt, y)
    end

    y0 = f.(Fields.coordinate_field(space), 0.0)
    y1 = f.(Fields.coordinate_field(space), 1.0)

    # Solve the ODE operator
    prob = ODEProblem(dfdt!, y0, (0.0, 1.0))
    sol = solve(prob, Tsit5(), reltol = 1e-6)

    @test norm(sol(1.0) .- y1) <= 1e-6 * norm(y1)
end

# implicit solvers
# require use of FieldVector
@testset "DiffEq Implicit Solvers" begin
    domain = Domains.IntervalDomain(
        Geometry.ZPoint(0.0),
        Geometry.ZPoint(pi);
        boundary_tags = (:left, :right),
    )
    mesh = Meshes.IntervalMesh(domain; nelems = 16)

    center_space = Spaces.CenterFiniteDifferenceSpace(mesh)
    # y(z, t) = z*exp(α*t)
    # dydt(z, t) = α*y(t)
    zcenters = Fields.coordinate_field(center_space).z
    y = Fields.FieldVector(z = zcenters)

    function f!(dydt, y, α, t)
        dydt.z .= α .* y.z
    end

    function f_jac!(J, y, α, t)
        copyto!(J, α * LinearAlgebra.I)
    end

    prob = ODEProblem(ODEFunction(f!), copy(y), (0.0, 1.0), 0.1)
    sol = solve(prob, SSPRK22(), dt = 0.01)

    sol = solve(prob, ImplicitEuler(), reltol = 1e-6)
    y1 = similar(y)
    y1 .= y .* exp(0.1)

    @test norm(sol[end] .- y1) <= 1e-4 * norm(y1)
end


text:
struct Node
    x::Float64
    y::Float64
    id::String
end

mutable struct Way
    id::String
    visible::Bool
    version::Int
    changeset::String
    timestamp::String
    user::String
    uid::String
    nodes::Array{Node}
    tags::Dict
    Way() = new()
end

mutable struct Relation
    id::String
    visible::Bool
    version::Int
    changeset::String
    timestamp::String
    user::String
    uid::String
    members::Array{Any}
    tags::Dict
    Relation() = new()
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## work with AbstractVector{UInt8} via I/O primitives ##

# Stateful string
mutable struct GenericIOBuffer{T<:AbstractVector{UInt8}} <: IO
    data::T # T should support: getindex, setindex!, length, copyto!, and resize!
    readable::Bool
    writable::Bool
    seekable::Bool # if not seekable, implementation is free to destroy (compact) past read data
    append::Bool # add data at end instead of at pointer
    size::Int # end pointer (and write pointer if append == true)
    maxsize::Int # fixed array size (typically pre-allocated)
    ptr::Int # read (and maybe write) pointer
    mark::Int # reset mark location for ptr (or <0 for no mark)

    function GenericIOBuffer{T}(data::T, readable::Bool, writable::Bool, seekable::Bool, append::Bool,
                                maxsize::Integer) where T<:AbstractVector{UInt8}
        require_one_based_indexing(data)
        new(data,readable,writable,seekable,append,length(data),maxsize,1,-1)
    end
end
const IOBuffer = GenericIOBuffer{Vector{UInt8}}

function GenericIOBuffer(data::T, readable::Bool, writable::Bool, seekable::Bool, append::Bool,
                         maxsize::Integer) where T<:AbstractVector{UInt8}
    GenericIOBuffer{T}(data, readable, writable, seekable, append, maxsize)
end

# allocate Vector{UInt8}s for IOBuffer storage that can efficiently become Strings
StringVector(n::Integer) = unsafe_wrap(Vector{UInt8}, _string_n(n))

# IOBuffers behave like Files. They are typically readable and writable. They are seekable. (They can be appendable).

"""
    IOBuffer([data::AbstractVector{UInt8}]; keywords...) -> IOBuffer

Create an in-memory I/O stream, which may optionally operate on a pre-existing array.

It may take optional keyword arguments:
- `read`, `write`, `append`: restricts operations to the buffer; see `open` for details.
- `truncate`: truncates the buffer size to zero length.
- `maxsize`: specifies a size beyond which the buffer may not be grown.
- `sizehint`: suggests a capacity of the buffer (`data` must implement `sizehint!(data, size)`).

When `data` is not given, the buffer will be both readable and writable by default.

# Examples
```jldoctest
julia> io = IOBuffer();

julia> write(io, "JuliaLang is a GitHub organization.", " It has many members.")
56

julia> String(take!(io))
"JuliaLang is a GitHub organization. It has many members."

julia> io = IOBuffer(b"JuliaLang is a GitHub organization.")
IOBuffer(data=UInt8[...], readable=true, writable=false, seekable=true, append=false, size=35, maxsize=Inf, ptr=1, mark=-1)

julia> read(io, String)
"JuliaLang is a GitHub organization."

julia> write(io, "This isn't writable.")
ERROR: ArgumentError: ensureroom failed, IOBuffer is not writeable

julia> io = IOBuffer(UInt8[], read=true, write=true, maxsize=34)
IOBuffer(data=UInt8[...], readable=true, writable=true, seekable=true, append=false, size=0, maxsize=34, ptr=1, mark=-1)

julia> write(io, "JuliaLang is a GitHub organization.")
34

julia> String(take!(io))
"JuliaLang is a GitHub organization"

julia> length(read(IOBuffer(b"data", read=true, truncate=false)))
4

julia> length(read(IOBuffer(b"data", read=true, truncate=true)))
0
```
"""
function IOBuffer(
        data::AbstractVector{UInt8};
        read::Union{Bool,Nothing}=nothing,
        write::Union{Bool,Nothing}=nothing,
        append::Union{Bool,Nothing}=nothing,
        truncate::Union{Bool,Nothing}=nothing,
        maxsize::Integer=typemax(Int),
        sizehint::Union{Integer,Nothing}=nothing)
    if maxsize < 0
        throw(ArgumentError("negative maxsize"))
    end
    if sizehint !== nothing
        sizehint!(data, sizehint)
    end
    flags = open_flags(read=read, write=write, append=append, truncate=truncate)
    buf = GenericIOBuffer(data, flags.read, flags.write, true, flags.append, Int(maxsize))
    if flags.truncate
        buf.size = 0
    end
    return buf
end

function IOBuffer(;
        read::Union{Bool,Nothing}=true,
        write::Union{Bool,Nothing}=true,
        append::Union{Bool,Nothing}=nothing,
        truncate::Union{Bool,Nothing}=true,
        maxsize::Integer=typemax(Int),
        sizehint::Union{Integer,Nothing}=nothing)
    size = sizehint !== nothing ? Int(sizehint) : maxsize != typemax(Int) ? Int(maxsize) : 32
    flags = open_flags(read=read, write=write, append=append, truncate=truncate)
    buf = IOBuffer(
        StringVector(size),
        read=flags.read,
        write=flags.write,
        append=flags.append,
        truncate=flags.truncate,
        maxsize=maxsize)
    fill!(buf.data, 0)
    return buf
end

# PipeBuffers behave like Unix Pipes. They are typically readable and writable, they act appendable, and are not seekable.

"""
    PipeBuffer(data::Vector{UInt8}=UInt8[]; maxsize::Integer = typemax(Int))

An [`IOBuffer`](@ref) that allows reading and performs writes by appending.
Seeking and truncating are not supported.
See [`IOBuffer`](@ref) for the available constructors.
If `data` is given, creates a `PipeBuffer` to operate on a data vector,
optionally specifying a size beyond which the underlying `Array` may not be grown.
"""
PipeBuffer(data::Vector{UInt8}=UInt8[]; maxsize::Int = typemax(Int)) =
    GenericIOBuffer(data,true,true,false,true,maxsize)
PipeBuffer(maxsize::Integer) = (x = PipeBuffer(StringVector(maxsize), maxsize = maxsize); x.size=0; x)

function copy(b::GenericIOBuffer)
    ret = typeof(b)(b.writable ? copy(b.data) : b.data,
                    b.readable, b.writable, b.seekable, b.append, b.maxsize)
    ret.size = b.size
    ret.ptr  = b.ptr
    return ret
end

show(io::IO, b::GenericIOBuffer) = print(io, "IOBuffer(data=UInt8[...], ",
                                      "readable=", b.readable, ", ",
                                      "writable=", b.writable, ", ",
                                      "seekable=", b.seekable, ", ",
                                      "append=",   b.append, ", ",
                                      "size=",     b.size, ", ",
                                      "maxsize=",  b.maxsize == typemax(Int) ? "Inf" : b.maxsize, ", ",
                                      "ptr=",      b.ptr, ", ",
                                      "mark=",     b.mark, ")")

@noinline function _throw_not_readable()
    # See https://github.com/JuliaLang/julia/issues/29688.
    throw(ArgumentError("read failed, IOBuffer is not readable"))
end

function unsafe_read(from::GenericIOBuffer, p::Ptr{UInt8}, nb::UInt)
    from.readable || _throw_not_readable()
    avail = bytesavailable(from)
    adv = min(avail, nb)
    GC.@preserve from unsafe_copyto!(p, pointer(from.data, from.ptr), adv)
    from.ptr += adv
    if nb > avail
        throw(EOFError())
    end
    nothing
end

function peek(from::GenericIOBuffer, T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}})
    from.readable || _throw_not_readable()
    avail = bytesavailable(from)
    nb = sizeof(T)
    if nb > avail
        throw(EOFError())
    end
    GC.@preserve from begin
        ptr::Ptr{T} = pointer(from.data, from.ptr)
        x = unsafe_load(ptr)
    end
    return x
end

function read(from::GenericIOBuffer, T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}})
    x = peek(from, T)
    from.ptr += sizeof(T)
    return x
end

function read_sub(from::GenericIOBuffer, a::AbstractArray{T}, offs, nel) where T
    require_one_based_indexing(a)
    from.readable || _throw_not_readable()
    if offs+nel-1 > length(a) || offs < 1 || nel < 0
        throw(BoundsError())
    end
    if isbitstype(T) && isa(a,Array)
        nb = UInt(nel * sizeof(T))
        GC.@preserve a unsafe_read(from, pointer(a, offs), nb)
    else
        for i = offs:offs+nel-1
            a[i] = read(to, T)
        end
    end
    return a
end

@inline function read(from::GenericIOBuffer, ::Type{UInt8})
    from.readable || _throw_not_readable()
    ptr = from.ptr
    size = from.size
    if ptr > size
        throw(EOFError())
    end
    @inbounds byte = from.data[ptr]::UInt8
    from.ptr = ptr + 1
    return byte
end

function peek(from::GenericIOBuffer, ::Type{UInt8})
    from.readable || _throw_not_readable()
    if from.ptr > from.size
        throw(EOFError())
    end
    return from.data[from.ptr]::UInt8
end

read(from::GenericIOBuffer, ::Type{Ptr{T}}) where {T} = convert(Ptr{T}, read(from, UInt))

isreadable(io::GenericIOBuffer) = io.readable
iswritable(io::GenericIOBuffer) = io.writable

# TODO: GenericIOBuffer is not iterable, so doesn't really have a length.
# This should maybe be sizeof() instead.
#length(io::GenericIOBuffer) = (io.seekable ? io.size : bytesavailable(io))
bytesavailable(io::GenericIOBuffer) = io.size - io.ptr + 1
position(io::GenericIOBuffer) = io.ptr-1

function skip(io::GenericIOBuffer, n::Integer)
    seekto = io.ptr + n
    n < 0 && return seek(io, seekto-1) # Does error checking
    io.ptr = min(seekto, io.size+1)
    return io
end

function seek(io::GenericIOBuffer, n::Integer)
    if !io.seekable
        ismarked(io) || throw(ArgumentError("seek failed, IOBuffer is not seekable and is not marked"))
        n == io.mark || throw(ArgumentError("seek failed, IOBuffer is not seekable and n != mark"))
    end
    # TODO: REPL.jl relies on the fact that this does not throw (by seeking past the beginning or end
    #       of an GenericIOBuffer), so that would need to be fixed in order to throw an error here
    #(n < 0 || n > io.size) && throw(ArgumentError("Attempted to seek outside IOBuffer boundaries."))
    #io.ptr = n+1
    io.ptr = max(min(n+1, io.size+1), 1)
    return io
end

function seekend(io::GenericIOBuffer)
    io.ptr = io.size+1
    return io
end

function truncate(io::GenericIOBuffer, n::Integer)
    io.writable || throw(ArgumentError("truncate failed, IOBuffer is not writeable"))
    io.seekable || throw(ArgumentError("truncate failed, IOBuffer is not seekable"))
    n < 0 && throw(ArgumentError("truncate failed, n bytes must be ≥ 0, got $n"))
    n > io.maxsize && throw(ArgumentError("truncate failed, $(n) bytes is exceeds IOBuffer maxsize $(io.maxsize)"))
    if n > length(io.data)
        resize!(io.data, n)
    end
    io.data[io.size+1:n] .= 0
    io.size = n
    io.ptr = min(io.ptr, n+1)
    ismarked(io) && io.mark > n && unmark(io)
    return io
end

function compact(io::GenericIOBuffer)
    io.writable || throw(ArgumentError("compact failed, IOBuffer is not writeable"))
    io.seekable && throw(ArgumentError("compact failed, IOBuffer is seekable"))
    local ptr::Int, bytes_to_move::Int
    if ismarked(io) && io.mark < io.ptr
        if io.mark == 0 return end
        ptr = io.mark
        bytes_to_move = bytesavailable(io) + (io.ptr-io.mark)
    else
        ptr = io.ptr
        bytes_to_move = bytesavailable(io)
    end
    copyto!(io.data, 1, io.data, ptr, bytes_to_move)
    io.size -= ptr - 1
    io.ptr -= ptr - 1
    io.mark -= ptr - 1
    return io
end

@noinline function ensureroom_slowpath(io::GenericIOBuffer, nshort::UInt)
    io.writable || throw(ArgumentError("ensureroom failed, IOBuffer is not writeable"))
    if !io.seekable
        if !ismarked(io) && io.ptr > 1 && io.size <= io.ptr - 1
            io.ptr = 1
            io.size = 0
        else
            datastart = ismarked(io) ? io.mark : io.ptr
            if (io.size+nshort > io.maxsize) ||
                (datastart > 4096 && datastart > io.size - io.ptr) ||
                (datastart > 262144)
                # apply somewhat arbitrary heuristics to decide when to destroy
                # old, read data to make more room for new data
                compact(io)
            end
        end
    end
    return
end

@inline ensureroom(io::GenericIOBuffer, nshort::Int) = ensureroom(io, UInt(nshort))
@inline function ensureroom(io::GenericIOBuffer, nshort::UInt)
    if !io.writable || (!io.seekable && io.ptr > 1)
        ensureroom_slowpath(io, nshort)
    end
    n = min((nshort % Int) + (io.append ? io.size : io.ptr-1), io.maxsize)
    l = length(io.data)
    if n > l
        _growend!(io.data, (n - l) % UInt)
    end
    return io
end

eof(io::GenericIOBuffer) = (io.ptr-1 == io.size)

@noinline function close(io::GenericIOBuffer{T}) where T
    io.readable = false
    io.writable = false
    io.seekable = false
    io.size = 0
    io.maxsize = 0
    io.ptr = 1
    io.mark = -1
    if io.writable
        resize!(io.data, 0)
    end
    nothing
end

isopen(io::GenericIOBuffer) = io.readable || io.writable || io.seekable || bytesavailable(io) > 0

"""
    take!(b::IOBuffer)

Obtain the contents of an `IOBuffer` as an array, without copying. Afterwards, the
`IOBuffer` is reset to its initial state.

# Examples
```jldoctest
julia> io = IOBuffer();

julia> write(io, "JuliaLang is a GitHub organization.", " It has many members.")
56

julia> String(take!(io))
"JuliaLang is a GitHub organization. It has many members."
```
"""
function take!(io::GenericIOBuffer)
    ismarked(io) && unmark(io)
    if io.seekable
        nbytes = io.size
        data = copyto!(StringVector(nbytes), 1, io.data, 1, nbytes)
    else
        nbytes = bytesavailable(io)
        data = read!(io,StringVector(nbytes))
    end
    if io.writable
        io.ptr = 1
        io.size = 0
    end
    return data
end
function take!(io::IOBuffer)
    ismarked(io) && unmark(io)
    if io.seekable
        data = io.data
        if io.writable
            maxsize = (io.maxsize == typemax(Int) ? 0 : min(length(io.data),io.maxsize))
            io.data = StringVector(maxsize)
        else
            data = copy(data)
        end
        resize!(data,io.size)
    else
        nbytes = bytesavailable(io)
        a = StringVector(nbytes)
        data = read!(io, a)
    end
    if io.writable
        io.ptr = 1
        io.size = 0
    end
    return data
end

function write(to::GenericIOBuffer, from::GenericIOBuffer)
    if to === from
        from.ptr = from.size + 1
        return 0
    end
    written::Int = write_sub(to, from.data, from.ptr, bytesavailable(from))
    from.ptr += written
    return written
end

function unsafe_write(to::GenericIOBuffer, p::Ptr{UInt8}, nb::UInt)
    ensureroom(to, nb)
    ptr = (to.append ? to.size+1 : to.ptr)
    written = Int(min(nb, length(to.data) - ptr + 1))
    towrite = written
    d = to.data
    while towrite > 0
        @inbounds d[ptr] = unsafe_load(p)
        ptr += 1
        p += 1
        towrite -= 1
    end
    to.size = max(to.size, ptr - 1)
    if !to.append
        to.ptr += written
    end
    return written
end

function write_sub(to::GenericIOBuffer, a::AbstractArray{UInt8}, offs, nel)
    require_one_based_indexing(a)
    if offs+nel-1 > length(a) || offs < 1 || nel < 0
        throw(BoundsError())
    end
    GC.@preserve a unsafe_write(to, pointer(a, offs), UInt(nel))
end

@inline function write(to::GenericIOBuffer, a::UInt8)
    ensureroom(to, UInt(1))
    ptr = (to.append ? to.size+1 : to.ptr)
    if ptr > to.maxsize
        return 0
    else
        to.data[ptr] = a
    end
    to.size = max(to.size, ptr)
    if !to.append
        to.ptr += 1
    end
    return sizeof(UInt8)
end

readbytes!(io::GenericIOBuffer, b::Array{UInt8}, nb=length(b)) = readbytes!(io, b, Int(nb))
function readbytes!(io::GenericIOBuffer, b::Array{UInt8}, nb::Int)
    nr = min(nb, bytesavailable(io))
    if length(b) < nr
        resize!(b, nr)
    end
    read_sub(io, b, 1, nr)
    return nr
end
read(io::GenericIOBuffer) = read!(io,StringVector(bytesavailable(io)))
readavailable(io::GenericIOBuffer) = read(io)
read(io::GenericIOBuffer, nb::Integer) = read!(io,StringVector(min(nb, bytesavailable(io))))

function occursin(delim::UInt8, buf::IOBuffer)
    p = pointer(buf.data, buf.ptr)
    q = GC.@preserve buf ccall(:memchr,Ptr{UInt8},(Ptr{UInt8},Int32,Csize_t),p,delim,bytesavailable(buf))
    return q != C_NULL
end

function occursin(delim::UInt8, buf::GenericIOBuffer)
    data = buf.data
    for i = buf.ptr:buf.size
        @inbounds b = data[i]
        b == delim && return true
    end
    return false
end

function readuntil(io::GenericIOBuffer, delim::UInt8; keep::Bool=false)
    lb = 70
    A = StringVector(lb)
    nread = 0
    nout = 0
    data = io.data
    for i = io.ptr : io.size
        @inbounds b = data[i]
        nread += 1
        if keep || b != delim
            nout += 1
            if nout > lb
                lb = nout*2
                resize!(A, lb)
            end
            @inbounds A[nout] = b
        end
        if b == delim
            break
        end
    end
    io.ptr += nread
    if lb != nout
        resize!(A, nout)
    end
    A
end

# copy-free crc32c of IOBuffer:
function _crc32c(io::IOBuffer, nb::Integer, crc::UInt32=0x00000000)
    nb < 0 && throw(ArgumentError("number of bytes to checksum must be ≥ 0, got $nb"))
    io.readable || _throw_not_readable()
    n = min(nb, bytesavailable(io))
    n == 0 && return crc
    crc = GC.@preserve io unsafe_crc32c(pointer(io.data, io.ptr), n, crc)
    io.ptr += n
    return crc
end
_crc32c(io::IOBuffer, crc::UInt32=0x00000000) = _crc32c(io, bytesavailable(io), crc)

max_stars_repo_name: ma-laforge/julia

max_issues_repo_name: ma-laforge/julia

max_forks_repo_name: ma-laforge/julia


text:
### A Pluto.jl notebook ###
# v0.18.0

using Markdown
using InteractiveUtils

# ╔═╡ d0379e1e-8a6b-11ec-0ae6-3506c99b8004
import PlutoLinks, PlutoHooks

# ╔═╡ 33e53084-ddca-4147-ae09-fd9c0deec242
"""
Run a code block using debounce, i.e. only run it when variable list has stabilized through `timeout`.

```julia
@use_debounce([var], 1) do
	sleep(2.)
	@info "stable variable \$var"
end
```
"""
macro use_debounce(f, deps, timeout=1)
	quote
		let
			value, set_value = @PlutoHooks.use_state(nothing)

			@PlutoLinks.use_task($(esc(deps))) do
				sleep($(esc(timeout)))
				set_value($(esc(f))())
			end
			value
		end
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoHooks = "0ff47ea0-7a50-410d-8455-4348d5de0774"
PlutoLinks = "0ff47ea0-7a50-410d-8455-4348d5de0420"

[compat]
PlutoHooks = "~0.0.4"
PlutoLinks = "~0.1.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "b55aae9a2bf436fc797d9c253a900913e0e90178"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.9.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "6b0440822974cab904c8b14d79743565140567f6"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.2.1"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoHooks]]
deps = ["InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "bda062fe28bab89e96281ba3047afa515b06e8e1"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.4"

[[deps.PlutoLinks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "PlutoHooks", "Revise", "UUIDs"]
git-tree-sha1 = "f2a6fcc4cd926a9f97ddf102e97afcc2a29b720a"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0420"
version = "0.1.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "2f9d4d6679b5f0394c52731db3794166f49d5131"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.3.1"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═d0379e1e-8a6b-11ec-0ae6-3506c99b8004
# ╠═33e53084-ddca-4147-ae09-fd9c0deec242
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
macro linklibrary_add(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_add_e))

  posittype = Symbol(:UInt, n)

  esc(:(padd = (a, b) -> ccall( ($fnname, $path), $posittype, ($posittype, $posittype), a, b )))
end

#spot unit testing
@linklibrary_add(8,0)
@test padd(0x00, 0x00) == 0x00  #0 + 0 == 0
@test padd(0x00, 0x40) == 0x40  #0 + 1 == 1
@test padd(0x40, 0x40) == 0x60  #1 + 1 == 2
@test padd(0x40, 0xc0) == 0x00  #1 + -1 == 0
@test padd(0x80, 0xc0) == 0x80  #Inf + 1 == Inf

#comprehensive integration testing
@linklibrary_add(8,0)
for x = 0x00:0xFF, y = 0x00:0xFF
  if (x != 0x80 || y != 0x80) #ignore the one corner case that doesn't work.
    @test (x,y,padd(x, y)) == (x,y,reinterpret(UInt8, (reinterpret(Posit{8,0}, x) + reinterpret(Posit{8,0}, y))))
  end
end
@linklibrary_add(8,1)
for x = 0x00:0xFF, y = 0x00:0xFF
  if (x != 0x80 || y != 0x80) #ignore the one corner case that doesn't work.
    @test (x,y,padd(x, y)) == (x,y,reinterpret(UInt8, (reinterpret(Posit{8,1}, x) + reinterpret(Posit{8,1}, y))))
  end
end

@linklibrary_add(8,2)
for x = 0x00:0xFF, y = 0x00:0xFF
  if (x != 0x80 || y != 0x80) #ignore the one corner case that doesn't work.
    @test (x,y,padd(x, y)) == (x,y,reinterpret(UInt8, (reinterpret(Posit{8,2}, x) + reinterpret(Posit{8,2}, y))))
  end
end

# testing 16-bit number types.

macro test_big_bittype_add(N, ES)
  quote
    @linklibrary_add($N, $ES)

    println("testing addition on p$($N)e$($ES)")

    UType = UIntLookup[$N]

    for idx = 1:10000
      x = rand(UType)
      y = rand(UType)
      if (x != topbits($N) || y != topbits($N)) #ignore the one corner case that doesn't work.
        @test (x,y,padd(x, y)) == (x,y,reinterpret(UType, (reinterpret(Posit{$N,$ES}, x) + reinterpret(Posit{$N,$ES}, y))))
      end
    end
  end
end

@test_big_bittype_add(16,0)
@test_big_bittype_add(16,1)
@test_big_bittype_add(16,2)

@test_big_bittype_add(32,0)
@test_big_bittype_add(32,1)
@test_big_bittype_add(32,2)
@test_big_bittype_add(32,3)


text:
using Tokens
using Test

function testlexer()
end


text:
    ```
    Chlorine, a gas first used by the Germans during the war, as it is an extremely irritating and poisonous gas.
    However, it can also sterilize germ-infected water, rendering it safe for consumption.
    It's made in a lab by heating Manganese dioxide and concentrated hydrochloric acid.
    ```

    """
    Given that we use enough hydrochloric acid,
    how much chlorine in Grams (default)  or (L)iters at NTP do we get from x grams of Manganese dioxide.
    """
    function magnese_dioxide_to_chlorine(magneseGrams,unit = "grams")
        if lowercase(unit) == "grams"
            return (magneseGrams/86.94) * 71
        elseif lowercase(unit) == "litre" || unit == "liter"
            return (magneseGrams/86.94) * 22.4
        end
    end

    """
    How much hydrochloric acid do we need to convert X amount of Manganese dioxide into chlorine.
    """
    function hydrochloric_acid_needed(magneseGrams)
        hydroc = (36.5 * 4) * magneseGrams / 86.94
        return hydroc
    end


text:
#
#
#

function run_dsatur_heuristc(infile_name::ASCIIString)
	# Read the graph from the disk (as adjacency list)
	G = simpleGraph(infile_name)

	# Apply the dsatur method
	number_color, vertex_color = dsatur(G.Graph)

	# Improve the coloring using Iterated Greedy techniques

	# Improve the coloring using min-conflicts local search

	# Return the coloring

end

function run_optdsatur_heuristc()
end

function dsatur(G::Array{Array{Int64},1})
	l = length(G)
	vertex_color = zeros(Int,l)
	degrees = Array{Tuple{Int,Int}}(l)
	satur = zeros(Int,l)
	uncolored = Set(collect(1:l))
	color_counter = 1

	for i in 1:l
		degrees[i] = (length(G[i]), i)
	end

	#coloring max degree vertex

	max_deg_vertex = maximum(degrees)[2]
	vertex_color[max_deg_vertex] = color_counter

	delete!(uncolored, max_deg_vertex)

	#Increase saturation of max_degree_vertex neighbours
	for j in G[max_deg_vertex]
		satur[j]+=1
	end

	coloring = 0
	while !isempty(uncolored)

		max_satur_degree_vertices = []
		aux = []
		for i in uncolored
			push!(aux, (satur[i], i))
		end

		max_sat = findmax(aux)[1][1]
		for i in aux
			if i[1] == max_sat
				push!(max_satur_degree_vertices, i[2])
			end
		end

		coloring = max_satur_degree_vertices[1]

		# If there are more than 1 vertex with max_satur, untie by max_degree
		if length(max_satur_degree_vertices) > 1
			max_degree = -1
			node_index = -1
			for i in max_satur_degree_vertices
				degree, node_index = degrees[i]
				if degree > max_degree
					coloring = node_index
					max_degree = degree
				end
			end
		end

		# Closure to return how many vertex have some color.
		function amount_color(node_indexes, color_number, vertex_color)
			color_counter = 0
			for i in node_indexes
				if vertex_color[i] == color_number
					color_counter +=1
				end
			end
			return color_counter
		end

		# Coloring
		for number_color in 1:color_counter
			if amount_color(G[coloring], number_color, vertex_color) == 0
				vertex_color[coloring] = number_color
				break
			end
		end

		#In case it wasnt colored
		if vertex_color[coloring] == 0
			color_counter+=1
			vertex_color[coloring] = color_counter
		end

		delete!(uncolored, coloring)

		for j in G[coloring]
			if amount_color(G[j], vertex_color[coloring], vertex_color) == 1
				satur[j] +=1
			end
		end
	end
	return number_color, vertex_color
end

function optdsatur()
end

function iteratedGreedy()

end

function localSearch()

end

type colorGraph
end


text:
"""
    update_partial_residuals!(r, y, x, perm, b, k)

A subroutine to compute the partial residuals `r = Y - X*b` in-place based on a permutation vector `perm` that indexes the nonzeroes in `b`.
"""
function update_partial_residuals!{T <: Float}(
    r    :: DenseVector{T},
    y    :: DenseVector{T},
    x    :: DenseMatrix{T},
    perm :: DenseArray{Int,1},
    b    :: DenseVector{T},
    k    :: Int
)
    copy!(r, y)
    @inbounds for j = 1:k
        idx = perm[j]
        @inbounds @simd for i in eachindex(y)
            r[i] += -b[idx]*x[i,idx]
        end
    end
    return nothing
end

"""
    exchange_leastsq!(v, x, y, k) -> b

This function minimizes the residual sum of squares

    0.5*sumabs2( y - x*b )

subject to `b` having no more than `k` nonzero components. The function will compute a `b` for a given value of `k`.
For optimal accuracy, this function should be run for multiple values of `k` over a path of model sizes.
For optimal performance, allocate `v` once and then reuse it between different runs of `exchange_leastsq!`.

Arguments:

- `v` is the `ELSQVariables` object housing all temporary arrays, including `b`
- `x` is the n x p statistical design matrix.
- `y` is the n-dimensional response vector.
- `k` is the desired number of nonzero components in `b`.

Optional Arguments:

- `n` and `p` are the dimensions of `x`; the former defaults to `length(y)` while the latter defaults to `size(x,2)`.
- `window` is an `Int` to dictate the dimension of the search window for potentially exchanging predictors.
   Defaults to `k`, which tells `exchange_leastsq!` to search all active predictors.
   Decreasing this quantity tells the algorithm to search through fewer current active predictors,
   which can decrease compute time but can also degrade model recovery performance.
- `max_iter` is the maximum permissible number of iterations. Defaults to `100`.
- `tol` is the convergence tolerance. Defaults to `1e-6`.
- `quiet` is a `Bool` to control output. Defaults to `false` (full output).
"""
function exchange_leastsq!{T <: Float}(
    v        :: ELSQVariables{T},
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    k        :: Int;
    n        :: Int  = length(y),
    p        :: Int  = size(x,2),
    window   :: Int  = k,         # <--- search *all* active predictors
    max_iter :: Int  = 100,
    tol      :: T    = convert(T, 1e-6),
    quiet    :: Bool = true
)
    # initial value for previous residual sum of squares
    old_rss = oftype(tol,Inf)

    # obtain top r components of bvec in magnitude
    selectperm!(v.perm, v.b, k, by=abs, rev=true, initialized=true)

    # compute partial residuals based on top r components of perm vector
    update_partial_residuals!(v.r, y, x, v.perm, v.b, k)

    # save value of RSS before starting algorithm
    rss = sumabs2(v.r) / 2

    # compute inner products of X and residuals
    # this is basically the negative gradient
    BLAS.gemv!('T', one(T), x, v.r, zero(T), v.df)

    # outer loop controls number of total iterations for algorithm run on one r
    for iter = 1:(max_iter)

        # output algorithm progress to console
        quiet || println("\titer = ", iter, ", RSS = ", rss)

        # middle loop tests each of top r parameters (by magnitude?)
        for i = abs(k-window+1):k

            # save information for current value of i
            l     = v.perm[i]
            betal = v.b[l]
            #update_col!(v.tempn, x, l)  # tempn now holds X[:,l]
            copy!(v.tempn, view(x, :, l))

            # if necessary, compute inner products of current predictor against all other predictors
            # store this information in Dict inner
            # for current index, hold dot products in memory for duration of inner loop
            get_inner_product!(v.dotprods, v.tempn, v, x, l)

            # subroutine compares current predictor i against all predictors k+1, k+2, ..., p
            # these predictors are candidates for inclusion in set
            # _exlstsq_innerloop! find best new predictor r
            a, b, r, adb = _exlstsq_innerloop!(v, k, i, p, tol)

            # now want to update residuals with current best predictor
            m = update_current_best_predictor!(v, x, betal, adb, r)

            # if necessary, compute inner product of current predictor against all other predictors
            # save in our Dict for future reference
            get_inner_product!(v.tempp, v.tempn2, v, x, m)

            # also update df
            axpymbz!(v.df, betal, v.dotprods, adb, v.tempp)

            # now swap best predictor with current predictor
            _swap_predictors!(v, i, r, m, adb)

        end # end middle loop over predictors

        # update residual sum of squares
        rss = sumabs2(v.r) / 2

        # test for descent failure
        # if no descent failure, then test for convergence
        # if not converged, then save RSS and check finiteness
        # if not converged and still finite, then save RSS and continue
        ascent    = rss > old_rss + tol
        converged = abs(old_rss - rss) / abs(old_rss + 1) < tol
        old_rss   = rss
        ascent && print_descent_error(iter, old_rss, rss)
        converged && return nothing

        check_finiteness(rss)

    end # end outer iteration loop

    # at this point, maximum iterations reached
    # warn and return
    print_maxiter(max_iter, rss)
    return nothing

end # end exchange_leastsq!

"""
    exlstsq(x, y) -> ELSQResults

For a statistical model `b`, `exlstsq` minimizes the residual sum of squares

    0.5*sumabs2( y - x*b )

subject to `b` having no more than `k` nonzero components. `exlstsq` will compute a `b` for several model sizes `k`.

Arguments:

- `x` is the n x p statistical design matrix.
- `y` is the n-dimensional response vector.

Optional Arguments:

- `v` is the `ELSQVariables` object housing all temporary arrays, including `b`.
- `models` is the integer vector of model sizes to test. It defaults to `collect(1:p)`, where `p = min(20, size(x,2))`.
- `window` is an `Int` to dictate the maximum size of the search window for potentially exchanging predictors.
   Defaults to `max(20, min(maximum(models), size(x,2)))` or the greater of `20` and `r`, where `r` is the lesser of the maximum model size and the number of predictors.
   Decreasing this quantity tells the algorithm to search through fewer current active predictors,
   which can decrease compute time but can also degrade model recovery performance.
- `max_iter` is the maximum permissible number of iterations. Defaults to `100`.
- `tol` is the convergence tolerance. Defaults to `1e-6`.
- `quiet` is a `Bool` to control output. Defaults to `false` (full output).
"""
function exlstsq{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T};
    v        :: ELSQVariables{T} = ELSQVariables(x, y),
    models   :: DenseVector{Int} = collect(1:min(20, size(x,2))),
    window   :: Int  = max(20, min(maximum(models), size(x,2))),
    max_iter :: Int  = 100,
    tol      :: T    = convert(T, 1e-6),
    quiet    :: Bool = true
)
    # dimensions of problem
    nmodels = length(models)
    n,p = size(x)

    # error checking
    errorcheck(x, y, tol, max_iter, window, p)

    # initialize sparse matrix to return
    betas = spzeros(T, p, nmodels)

    # loop through models
    for i in models

        # monitor output
        quiet || println("Testing model size $i.")
        exchange_leastsq!(v, x, y, i, window=min(window,i), max_iter=max_iter, tol=tol, quiet=quiet, n=n, p=p)
        betas[:,i] = sparse(sdata(v.b))
    end

    # return matrix of betas
    return betas
end


text:
module ImageView

using FixedPointNumbers, Colors, ColorVectorSpace, Images, StatsBase
using MappedArrays, RoundingIntegers
using Gtk.ShortNames, GtkReactive, Graphics
using Gtk.GConstants.GtkAlign: GTK_ALIGN_START, GTK_ALIGN_END, GTK_ALIGN_FILL
import AxisArrays
using AxisArrays: Axis, axisnames, axisvalues

import Images: scaleminmax

export AnnotationText, AnnotationPoint, AnnotationPoints,
       AnnotationLine, AnnotationLines, AnnotationBox
export CLim, annotate!, canvasgrid, imshow, imshow!, imshow_gui, imlink,
       roi, scalebar, slice2d

const AbstractGray{T} = Color{T,1}
const GrayLike = Union{AbstractGray,Number}
const FixedColorant{T<:FixedPoint} = Colorant{T}

include("slicing.jl")

"""
    CLim(cmin, cmax)

Specify contrast limits where `x <= cmin` will be rendered as black,
and `x >= cmax` will be rendered as white.
"""
struct CLim{T}
    min::T
    max::T
end
CLim(min, max) = CLim(promote(min, max)...)
Base.convert(::Type{CLim{T}}, clim::CLim) where {T} = CLim(convert(T, clim.min),
                                                    convert(T, clim.max))
Base.eltype(::CLim{T}) where {T} = T

"""
    closeall()

Closes all windows opened by ImageView2.
"""
function closeall()
    for (w, _) in window_wrefs
        destroy(w)
    end
    nothing
end

const window_wrefs = WeakKeyDict{Gtk.GtkWindowLeaf,Nothing}()

"""
    imshow()

Choose an image to display via a file dialog.
"""
imshow() = imshow(load(open_dialog("Pick an image to display")))

"""
    imshow!(canvas, img) -> drawsignal
    imshow!(canvas, img::Signal, zr::Signal{ZoomRegion}) -> drawsignal
    imshow!(frame::Frame, canvas, img::Signal, zr::Signal{ZoomRegion}) -> drawsignal
    imshow!(..., annotations=Signal(Dict{UInt,Any}()))

Display the image `img`, in the specified `canvas`. Use the version
with `zr` if you have already turned on rubber-banding or other
pan/zoom interactivity for `canvas`. Returns the Reactive `drawsignal`
used for updating the canvas.

If you supply `frame`, then the pixel aspect ratio will be set to that
of `pixelspacing(img)`.

With any of these forms, you may optionally supply `annotations`.

This only creates the `draw` method for `canvas`; mouse- or key-based
interactivity can be set up via [`imshow`](@ref) or, at a lower level,
using GtkReactive's tools:

- `init_zoom_rubberband`
- `init_zoom_scroll`
- `init_pan_scroll`
- `init_pan_drag`

# Example

```julia
using ImageView, GtkReactive, Gtk.ShortNames, TestImages
# Create a window with a canvas in it
win = Window()
c = canvas(UserUnit)
push!(win, c)
Gtk.showall(win)
# Load images
mri = testimage("mri")
# Display the image
imshow!(c, mri[:,:,1])
# Update with a different image
imshow!(c, mri[:,:,8])
"""
function imshow!(canvas::GtkReactive.Canvas{UserUnit},
                 imgsig::Signal,
                 zr::Signal{ZoomRegion{T}},
                 annotations::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}())) where T<:RInteger
    draw(canvas, imgsig, annotations) do cnvs, image, anns
        copy!(cnvs, image)
        set_coordinates(cnvs, value(zr))
        draw_annotations(cnvs, anns)
    end
end

function imshow!(frame::Frame,
                 canvas::GtkReactive.Canvas{UserUnit},
                 imgsig::Signal,
                 zr::Signal{ZoomRegion{T}},
                 annotations::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}())) where T<:RInteger
    draw(canvas, imgsig, annotations) do cnvs, image, anns
        copy!(cnvs, image)
        set_coordinates(cnvs, value(zr))
        set_aspect!(frame, image)
        draw_annotations(cnvs, anns)
    end
end

# Without a ZoomRegion, there's no risk that the apsect ratio needs to
# change dynamically, so it can be set once and left. Consequently we
# don't need `frame` variants of the remaining methods.
function imshow!(canvas::GtkReactive.Canvas,
                 imgsig::Signal,
                 annotations::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}()))
    draw(canvas, imgsig, annotations) do cnvs, image, anns
        copy!(cnvs, image)
        set_coordinates(cnvs, axes(image))
        draw_annotations(cnvs, anns)
    end
end

# Simple non-interactive image display
function imshow!(canvas::GtkReactive.Canvas,
                 img::AbstractMatrix,
                 annotations::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}()))
    draw(canvas, annotations) do cnvs, anns
        copy!(cnvs, img)
        set_coordinates(cnvs, axes(img))
        draw_annotations(cnvs, anns)
    end
    nothing
end

"""
    imshow(img; axes=(1,2), name="ImageView") -> guidict
    imshow(img, clim; axes=(1,2), name="ImageView") -> guidict
    imshow(img, clim, zoomregion, slicedata, annotations; axes=(1,2), name="ImageView") -> guidict

Display the image `img` in a new window titled with `name`, returning
a dictionary `guidict` containing any Reactive signals or GtkReactive
widgets. If the image is 3 or 4 dimensional, GUI controls will be
added for slicing along "extra" axes. By default the two-dimensional
slice containing axes 1 and 2 are shown, but that can be changed by
passing a different setting for `axes`.

If the image is grayscale, by default contrast is set by a
`scaleminmax` object whose end-points can be modified by
right-clicking on the image. If `clim == nothing`, the image's own
native contrast is used (`clamp01nan`).  You may also pass a custom
contrast function.

Finally, you may specify [`GtkReactive.ZoomRegion`](@ref) and
[`SliceData`](@ref) signals. See also [`roi`](@ref), as well as any
`annotations` that you wish to apply.
"""
function imshow(img::AbstractArray;
                axes=default_axes(img), name="ImageView", scalei=identity, aspect=:auto,
                kwargs...)
    imgmapped = kwhandler(_mappedarray(scalei, img), axes; kwargs...)
    zr, sd = roi(imgmapped, axes)
    v = slice2d(imgmapped, value(zr), sd)
    imshow(imgmapped, default_clim(v), zr, sd; name=name, aspect=aspect)
end

imshow(img::AbstractVector; kwargs...) = imshow(reshape(img, :, 1); kwargs...)

function imshow(c::GtkReactive.Canvas, img::AbstractMatrix, anns=Signal(Dict{UInt,Any}());
                kwargs...)
    f = parent(widget(c))
    imshow(f, c, img, default_clim(img), roi(img, default_axes(img))..., anns; kwargs...)
end

function imshow(img::AbstractArray, clim;
                axes = default_axes(img), name="ImageView", aspect=:auto)
    imshow(img, clim, roi(img, axes)...; name=name, aspect=aspect)
end

function imshow(img::AbstractArray, clim,
                zr::Signal{ZoomRegion{T}}, sd::SliceData,
                anns=Signal(Dict{UInt,Any}());
                name="ImageView", aspect=:auto) where T
    v = slice2d(img, value(zr), sd)
    ps = map(abs, pixelspacing(v))
    csz = default_canvas_size(fullsize(value(zr)), ps[2]/ps[1])
    guidict = imshow_gui(csz, sd; name=name, aspect=aspect)
    guidict["hoverinfo"] = map(guidict["canvas"].mouse.motion; name="hoverinfo") do btn
        hoverinfo(guidict["status"], btn, img, sd)
    end

    roidict = imshow(guidict["frame"], guidict["canvas"], img,
                     wrap_signal(clim), zr, sd, anns)

    win = guidict["window"]
    Gtk.showall(win)
    dct = Dict("gui"=>guidict, "clim"=>clim, "roi"=>roidict, "annotations"=>anns)
    GtkReactive.gc_preserve(win, dct)
    return dct
end

function imshow(frame::Gtk.GtkFrame, canvas::GtkReactive.Canvas,
                img::AbstractArray, clim::Union{Nothing,Signal{<:CLim}},
                zr::Signal{ZoomRegion{T}}, sd::SliceData,
                anns::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}())) where T
    imgsig = map(zr, sd.signals...; name="imgsig") do r, s...
        while length(s) < 2
            s = (s..., 1)
        end
        for (h, ann) in value(anns)
            setvalid!(ann, s...)
        end
        slice2d(img, r, sd)
    end
    set_aspect!(frame, value(imgsig))
    imgc = prep_contrast(canvas, imgsig, clim)
    GtkReactive.gc_preserve(frame, imgc)

    roidict = imshow(frame, canvas, imgc, zr, anns)
    roidict["slicedata"] = sd
    roidict
end

# For things that are not AbstractArrays, we don't offer the clim
# option.  We also don't display hoverinfo, as there is no guarantee
# that one can quickly compute intensities at a point.
function imshow(img;
                axes = default_axes(img), name="ImageView", aspect=:auto)
    zr, sd = roi(img, axes)
    imshow(img, zr, sd; name=name, aspect=aspect)
end

function imshow(img,
                zr::Signal{ZoomRegion{T}}, sd::SliceData,
                anns=Signal(Dict{UInt,Any}());
                name="ImageView", aspect=:auto) where T
    v = slice2d(img, value(zr), sd)
    ps = map(abs, pixelspacing(v))
    csz = default_canvas_size(fullsize(value(zr)), ps[2]/ps[1])
    guidict = imshow_gui(csz, sd; name=name, aspect=aspect)

    roidict = imshow(guidict["frame"], guidict["canvas"], img, zr, sd, anns)

    win = guidict["window"]
    Gtk.showall(win)
    dct = Dict("gui"=>guidict, "roi"=>roidict)
    GtkReactive.gc_preserve(win, dct)
    return dct
end

function imshow(frame::Gtk.GtkFrame, canvas::GtkReactive.Canvas,
                img, zr::Signal{ZoomRegion{T}}, sd::SliceData,
                anns::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}())) where T
    imgsig = map(zr, sd.signals...; name="imgsig") do r, s...
        slice2d(img, r, sd)
    end
    set_aspect!(frame, value(imgsig))
    GtkReactive.gc_preserve(frame, imgsig)

    roidict = imshow(frame, canvas, imgsig, zr, anns)
    roidict["slicedata"] = sd
    GtkReactive.gc_preserve(frame, roidict)
    roidict
end

"""
    guidict = imshow_gui(canvassize; name="ImageView", aspect=:auto)
    guidict = imshow_gui(canvassize, slicedata, gridsize=(1,1); name="ImageView", aspect=:auto)

Create an image-viewer GUI. By default creates a single canvas, but
with custom `gridsize` you can create a grid of canvases. `canvassize
= (szx, szy)` describes the desired size of the (or each)
canvas. `slicedata` is an object created by [`roi`](@ref) that encodes
the necessary information for creating player widgets for viewing
multidimensional images.
"""
function imshow_gui(canvassize::Tuple{Int,Int},
                    sd::SliceData=SliceData{false}(),
                    gridsize::Tuple{Int,Int} = (1,1);
                    name = "ImageView", aspect=:auto)
    winsize = canvas_size(screen_size(), map(*, canvassize, gridsize))
    win = Window(name, winsize...)
    window_wrefs[win] = nothing
    signal_connect(win, :destroy) do w
        delete!(window_wrefs, win)
    end
    vbox = Box(:v)
    push!(win, vbox)
    if gridsize == (1,1)
        frames, canvases = frame_canvas(aspect)
        g = frames
    else
        g, frames, canvases = canvasgrid(gridsize, aspect)
    end
    push!(vbox, g)
    status = Label("")
    set_gtk_property!(status, :halign, Gtk.GConstants.GtkAlign.START)
    push!(vbox, status)

    guidict = Dict("window"=>win, "vbox"=>vbox, "frame"=>frames, "status"=>status,
                   "canvas"=>canvases)

    # Add the player controls
    if !isempty(sd)
        players = [player(sd.signals[i], axisvalues(sd.axs[i])[1]; id=i) for i = 1:length(sd)]
        guidict["players"] = players
        hbox = Box(:h)
        for p in players
            push!(hbox, frame(p))
        end
        push!(guidict["vbox"], hbox)
    end

    guidict
end

fullsize(zr::ZoomRegion) =
    map(i->length(UnitRange{Int}(i)), (zr.fullview.y, zr.fullview.x))

"""
    grid, frames, canvases = canvasgrid((ny, nx))

Create a grid of `ny`-by-`nx` canvases for drawing. `grid` is a
GtkGrid layout, `frames` is an `ny`-by-`nx` array of
GtkAspectRatioFrames that contain each canvas, and `canvases` is an
`ny`-by-`nx` array of canvases.
"""
function canvasgrid(gridsize::Tuple{Int,Int}, aspect=:auto)
    g = Grid()
    frames = Matrix{Any}(undef, gridsize)
    canvases = Matrix{Any}(undef, gridsize)
    for j = 1:gridsize[2], i = 1:gridsize[1]
        f, c = frame_canvas(aspect)
        g[j,i] = f
        frames[i,j] = f
        canvases[i,j] = c
    end
    return g, frames, canvases
end

function frame_canvas(aspect)
    f = aspect==:none ? Frame() : AspectFrame("", 0.5, 0.5, 1)
    set_gtk_property!(f, :expand, true)
    set_gtk_property!(f, :shadow_type, Gtk.GConstants.GtkShadowType.NONE)
    c = canvas(UserUnit)
    push!(f, widget(c))
    f, c
end

"""
    imshow(canvas, imgsig::Signal) -> guidict
    imshow(canvas, imgsig::Signal, zr::Signal{ZoomRegion}) -> guidict
    imshow(frame::Frame, canvas, imgsig::Signal, zr::Signal{ZoomRegion}) -> guidict

Display `imgsig` (a `Signal` of an image) in `canvas`, setting up
panning and zooming. Optionally include a `frame` for preserving
aspect ratio. `imgsig` must be two-dimensional (but can be a
Signal-view of a higher-dimensional object).

# Example

```julia
using ImageView, TestImages, Gtk
mri = testimage("mri");
# Create a canvas `c`. There are other approaches, like stealing one from a previous call
# to `imshow`, or using GtkReactive directly.
guidict = imshow_gui((300, 300))
c = guidict["canvas"];
# To see anything you have to call `showall` on the window (once)
Gtk.showall(guidict["window"])
# Create the image Signal
imgsig = Signal(mri[:,:,1]);
# Show it
imshow(c, imgsig)
# Now anytime you want to update, just push! a new image
push!(imgsig, mri[:,:,8])
```
"""
function imshow(canvas::GtkReactive.Canvas{UserUnit},
                imgsig::Signal,
                zr::Signal{ZoomRegion{T}}=Signal(ZoomRegion(value(imgsig))),
                anns::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}())) where T<:RInteger
    zoomrb = init_zoom_rubberband(canvas, zr)
    zooms = init_zoom_scroll(canvas, zr)
    pans = init_pan_scroll(canvas, zr)
    pand = init_pan_drag(canvas, zr)
    redraw = imshow!(canvas, imgsig, zr, anns)
    dct = Dict("image roi"=>imgsig, "zoomregion"=>zr, "zoom_rubberband"=>zoomrb,
               "zoom_scroll"=>zooms, "pan_scroll"=>pans, "pan_drag"=>pand,
               "redraw"=>redraw)
    GtkReactive.gc_preserve(widget(canvas), dct)
    dct
end

function imshow(frame::Frame,
                canvas::GtkReactive.Canvas{UserUnit},
                imgsig::Signal,
                zr::Signal{ZoomRegion{T}},
                anns::Signal{Dict{UInt,Any}}=Signal(Dict{UInt,Any}())) where T<:RInteger
    zoomrb = init_zoom_rubberband(canvas, zr)
    zooms = init_zoom_scroll(canvas, zr)
    pans = init_pan_scroll(canvas, zr)
    pand = init_pan_drag(canvas, zr)
    redraw = imshow!(frame, canvas, imgsig, zr, anns)
    dct = Dict("image roi"=>imgsig, "zoomregion"=>zr, "zoom_rubberband"=>zoomrb,
               "zoom_scroll"=>zooms, "pan_scroll"=>pans, "pan_drag"=>pand,
               "redraw"=>redraw)
    GtkReactive.gc_preserve(widget(canvas), dct)
    dct
end

"""
    imshowlabeled(img, label)

Display `img`, but showing the pixel's `label` rather than the color
value in the status bar.
"""
function imshowlabeled(img::AbstractArray, label::AbstractArray; proplist...)
    axes(img) == axes(label) || throw(DimensionMismatch("axes $(axes(label)) of label array disagree with axes $(axes(img)) of the image"))
    guidict = imshow(img; proplist...)
    gui = guidict["gui"]
    sd = guidict["roi"]["slicedata"]
    close(gui["hoverinfo"])
    gui["hoverinfo"] = map(gui["canvas"].mouse.motion; name="hoverinfo") do btn
        hoverinfo(gui["status"], btn, label, sd)
    end
    guidict
end

function hoverinfo(lbl, btn, img, sd::SliceData{transpose}) where transpose
    io = IOBuffer()
    y, x = round(Int, btn.position.y.val), round(Int, btn.position.x.val)
    axes = sliceinds(img, transpose ? (x, y) : (y, x), makeslices(sd)...)
    if checkbounds(Bool, img, axes...)
        print(io, '[', y, ',', x, "] ")
        show(IOContext(io, :compact=>true), img[axes...])
        set_gtk_property!(lbl, :label, String(take!(io)))
    else
        set_gtk_property!(lbl, :label, "")
    end
end

default_clim(img) = nothing
default_clim(img::AbstractArray{C}) where {C<:GrayLike} = _default_clim(img, eltype(C))
default_clim(img::AbstractArray{C}) where {C<:AbstractRGB} = _default_clim(img, eltype(C))
_default_clim(img, ::Type{Bool}) = nothing
_default_clim(img, ::Type{T}) where {T} = _deflt_clim(img)
function _deflt_clim(img::AbstractMatrix)
    minval = nanz(minfinite(img))
    maxval = nanz(maxfinite(img))
    if minval == maxval
        minval = zero(typeof(minval))
        maxval = one(typeof(maxval))
    end
    Signal(CLim(saferound(gray(minval)), saferound(gray(maxval))); name="CLim")
end

function _deflt_clim(img::AbstractMatrix{T}) where {T<:AbstractRGB}
    minval = RGB(0.0,0.0,0.0)
    maxval = RGB(1.0,1.0,1.0)
    Signal(CLim(minval, maxval); name="CLim")
end

saferound(x::Integer) = convert(RInteger, x)
saferound(x) = x

default_axes(::AbstractVector) = (1,)
default_axes(img) = (1, 2)
default_axes(img::AxisArray) = axisnames(img)[[1,2]]

#default_view(img) = view(img, :, :, ntuple(d->1, ndims(img)-2)...)
#default_view(img::Signal) = default_view(value(img))

# default_slices(img) = ntuple(d->PlayerInfo(Signal(1), axes(img, d+2)), ndims(img)-2)

function histsignals(enabled::Signal, defaultimg, img::Signal, clim::Signal{CLim{T}}) where {T<:GrayLike}
    return [map(filterwhen(enabled, defaultimg, img), enabled; name="histsig") do image, _  # `enabled` fixes issue #168
        cl = value(clim)
        smin = float(nanz(min(minfinite(image), cl.min)))
        smax = float(nanz(max(maxfinite(image), cl.max)))
        if smax == smin
            smax = smin+1
        end
        rng = range(smin, stop=smax, length=300)
        fit(Histogram, mappedarray(nanz, vec(channelview(image))), rng; closed=:right)
    end]
end

channel_clim(f, clim::CLim{T}) where {T<:AbstractRGB} = CLim(f(clim.min), f(clim.max))
channel_clims(clim::CLim{T}) where {T<:AbstractRGB} = map(f->channel_clim(f, clim), (red, green, blue))

function mapped_channel_clims(clim::Signal{CLim{T}}) where {T<:AbstractRGB}
    inits = channel_clims(value(clim))
    rsig = map(x->channel_clim(red, x), clim; init=inits[1])
    gsig = map(x->channel_clim(green, x), clim; init=inits[1])
    bsig = map(x->channel_clim(blue, x), clim; init=inits[1])
    return [rsig;gsig;bsig]
end

function histsignals(enabled::Signal, defaultimg, img::Signal, clim::Signal{CLim{T}}) where {T<:AbstractRGB}
    rv = map(x->mappedarray(red, x), filterwhen(enabled, defaultimg, img); name="redview")
    gv = map(x->mappedarray(green,x), filterwhen(enabled, defaultimg, img); name="greenview")
    bv = map(x->mappedarray(blue, x), filterwhen(enabled, defaultimg, img); name="blueview")
    cls = mapped_channel_clims(clim) #note currently this gets called twice, also in contrast gui creation (a bit inefficient/awkward)
    histsigs = []
    push!(histsigs, histsignals(enabled, mappedarray(red, defaultimg), rv, cls[1])[1])
    push!(histsigs, histsignals(enabled, mappedarray(green, defaultimg), gv, cls[2])[1])
    push!(histsigs, histsignals(enabled, mappedarray(blue, defaultimg), bv, cls[3])[1])
    return histsigs
end

function scaleminmax(::Type{Tout}, cmin::AbstractRGB{T}, cmax::AbstractRGB{T}) where {T,Tout}
    r = scaleminmax(T, red(cmin), red(cmax))
    g = scaleminmax(T, green(cmin), green(cmax))
    b = scaleminmax(T, blue(cmin), blue(cmax))
    return x->Tout(nanz(r(red(x))), nanz(g(green(x))), nanz(b(blue(x))))
end

function safeminmax(cmin::T, cmax::T) where {T<:GrayLike}
    if !(cmin < cmax)
        cmax = cmin+1
    end
    return cmin, cmax
end

function safeminmax(cmin::T, cmax::T) where {T<:AbstractRGB}
    rmin, rmax = safeminmax(red(cmin), red(cmax))
    gmin, gmax = safeminmax(green(cmin), green(cmax))
    bmin, bmax = safeminmax(blue(cmin), blue(cmax))
    return T(rmin, gmin, bmin), T(rmax, gmax, bmax)
end

function prep_contrast(img::Signal, clim::Signal{CLim{T}}) where {T}
    # Set up the signals to calculate the histogram of intensity
    enabled = Signal(false; name="contrast_enabled") # skip hist calculation if the contrast gui isn't open
    histsigs = histsignals(enabled, value(img), img, clim)
    # Return a signal corresponding to the scaled image
    imgc = map(img, clim; name="clim-mapped image") do image, cl
        cmin, cmax = safeminmax(cl.min, cl.max)
        smm = scaleminmax(outtype(T), cmin, cmax)
        mappedarray(smm, image)
    end
    enabled, histsigs, imgc
end

outtype(::Type{T}) where T<:GrayLike         = Gray{N0f8}
outtype(::Type{C}) where C<:Color            = RGB{N0f8}
outtype(::Type{C}) where C<:TransparentColor = RGBA{N0f8}

function prep_contrast(canvas, img::Signal, clim::Signal{CLim{T}}) where T
    enabled, histsigs, imgsig = prep_contrast(img, clim)
    # Set up the right-click to open the contrast gui
    push!(canvas.preserved, create_contrast_popup(canvas, enabled, histsigs, clim))
    imgsig
end

prep_contrast(canvas, img::Signal, f) =
    map(image->mappedarray(f, image), img; name="f-mapped image")
prep_contrast(canvas, img::Signal{A}, ::Nothing) where {A<:AbstractArray} =
    prep_contrast(canvas, img, clamp01nan)
prep_contrast(canvas, img::Signal, ::Nothing) = img

nanz(x) = ifelse(isnan(x), zero(x), x)
nanz(x::FixedPoint) = x
nanz(x::Integer) = x

function create_contrast_popup(canvas, enabled, hists, clim)
    popupmenu = Menu()
    contrast = MenuItem("Contrast...")
    push!(popupmenu, contrast)
    Gtk.showall(popupmenu)
    push!(canvas.preserved, map(canvas.mouse.buttonpress; name="open contrast GUI") do btn
        if btn.button == 3 && btn.clicktype == BUTTON_PRESS
            popup(popupmenu, btn.gtkevent)
        end
    end)
    signal_connect(contrast, :activate) do widget
        push!(enabled, true)
        contrast_gui(enabled, hists, clim)
    end
end

function map_image_roi(img, zr::Signal{ZoomRegion{T}}, slices...) where T
    map(zr, slices...; name="map_image_roi") do r, s...
        cv = r.currentview
        view(img, UnitRange{Int}(cv.y), UnitRange{Int}(cv.x), s...)
    end
end
map_image_roi(img::Signal, zr::Signal{ZoomRegion{T}}, slices...) where {T} = img

function set_aspect!(frame::AspectFrame, image)
    ps = map(abs, pixelspacing(image))
    sz = map(length, axes(image))
    r = sz[2]*ps[2]/(sz[1]*ps[1])
    set_gtk_property!(frame, :ratio, r)
    nothing
end
set_aspect!(frame, image) = nothing

"""
    default_canvas_size(imagesize, pixelaspectratio=1) -> (xsz, ysz)

Compute the canvas size for an image of size `imagesize` with the
defined `pixelaspectratio`. Note that `imagesize` is supplied in
coordinate order, i.e., (y, x) order, whereas the returned canvas size
is in Gtk order, i.e., (x, y) order.
"""
default_canvas_size(imgsz::Tuple{Integer,Integer}, pixelaspectratio::Number=1) =
    pixelaspectratio >= 1 ? (round(Int, pixelaspectratio*imgsz[2]), Int(imgsz[1])) :
        (Int(imgsz[2]), round(Int, imgsz[1]/pixelaspectratio))

"""
    canvas_size(win, requested_size) -> (xsz, ysz)
    canvas_size(screensize, requested_size) -> (xsz, ysz)

Limit the requested canvas size by the screen size. Both the output
and `screensize` are supplied in Gtk order (x, y).

When supplying a GtkWindow `win`, the canvas size is limited to 60% of
the total screen size.
"""
function canvas_size(win::Gtk.GtkWindowLeaf, requestedsize_xy; minsize=100)
    ssz = screen_size(win)
    canvas_size(map(x->0.6*x, ssz), requestedsize_xy; minsize=minsize)
end

function canvas_size(screensize_xy, requestedsize_xy; minsize=100)
    f = minimum(map(/, screensize_xy, requestedsize_xy))
    if f > 1
        fmn = maximum(map(/, (minsize,minsize), requestedsize_xy))
        f = max(1, min(f, fmn))
    end
    (round(Int, f*requestedsize_xy[1]), round(Int, f*requestedsize_xy[2]))
end

function kwhandler(img, axs; flipx=false, flipy=false, kwargs...)
    if flipx || flipy
        inds = AbstractRange[axes(img)...]
        setrange!(inds, _axisdim(img, axs[1]), flipy)
        setrange!(inds, _axisdim(img, axs[2]), flipx)
        img = view(img, inds...)
    end
    # isempty(kwargs) || error("cannot handle ", kwargs)
    img
end
function setrange!(inds, ax::Integer, flip)
    ind = inds[ax]
    inds[ax] = flip ? (last(ind):-1:first(ind)) : ind
    inds
end
_axisdim(img, ax::Integer) = ax
_axisdim(img, ax::Axis) = axisdim(img, ax)
_axisdim(img, ax) = axisdim(img, Axis{ax})

isgray(img::AbstractArray{T}) where {T<:Real} = true
isgray(img::AbstractArray{T}) where {T<:AbstractGray} = true
isgray(img) = false

_mappedarray(f, img) = mappedarray(f, img)
_mappedarray(f, img::AxisArray) = AxisArray(mappedarray(f, img.data), AxisArrays.axes(img))
_mappedarray(f, img::ImageMeta) = shareproperties(img, _mappedarray(f, data(img)))

wrap_signal(x) = Signal(x)
wrap_signal(x::Signal) = x
wrap_signal(::Nothing) = nothing

include("link.jl")
include("contrast_gui.jl")
include("annotations.jl")

end # module


text:
function fetch_catalog(rcf, stagedir)
    # note: this call to read_photoobj_files considers only primary detections.
    catalog = SDSSIO.read_photoobj_files([rcf,], stagedir)

    # we're ignoring really faint sources entirely...not even using them to
    # render the background
    isnt_faint(entry) = maximum(entry.star_fluxes) >= MIN_FLUX
    filter(isnt_faint, catalog)
end

function in_box(entry::CatalogEntry, box::BoundingBox)
    box.ramin < entry.pos[1] < box.ramax &&
        box.decmin < entry.pos[2] < box.decmax
end

function invert_rcf_array(rcfs)
    max_run = maximum([rcf.run for rcf in rcfs])
    max_camcol = maximum([rcf.camcol for rcf in rcfs])
    max_field = maximum([rcf.field for rcf in rcfs])

    rcf_to_index = Array(Int64, max_run, max_camcol, max_field)

    fill!(rcf_to_index, -1)

    # this should be really fast, each node could do it
    for n in 1:length(rcfs)
        rcf = rcfs[n]
        rcf_to_index[rcf.run, rcf.camcol, rcf.field] = n
    end

    rcf_to_index
end

function load_images(box, rcfs, stagedir)
    # TODO: set `total_threads`. We may not want to use all threads possible
    # on all nodes, since this operation hits disk. Maybe just 1 thread per node?
    total_threads = 5
    num_fields = length(rcfs)
    fields_per_thread = ceil(Int, num_fields / total_threads)

    #TODO: make `images` a global array
    # (each cell of `images` contains B=5 tiled images)
    images = Array(Vector{TiledImage}, num_fields)

    #TODO: make `catalog_offset` a global array too.
    # It stores first index of each field's sources in the catalog array.
    catalog_offset = Array(Int64, num_fields)

    #TODO: make `task_offset` a global array too.
    # It stores first index of each field's tasks in the tasks array.
    task_offset = Array(Int64, num_fields)

    #TODO: use dtree to have each thread on each node
    #  load some fields.
    #Note: this loop and the next one over threads hit disk, so we may not want
    #  to use all available threads. (The third loop over threads is the
    #  only cpu intensive loop.)
    for m in 0:(total_threads - 1)
        n0 = 1 + m * fields_per_thread
        nend = min(num_fields, (m + 1) * fields_per_thread)

        for n in n0:nend
            # first, load the field array
            rcf = rcfs[n]
            raw_images = SDSSIO.load_field_images(rcf, stagedir)
            @assert(length(raw_images) == 5)
            images[n] = [TiledImage(img) for img in raw_images]

            # second, load the `catalog_offset` and `task_count` arrays with
            # a number of sources for each field.
            # (We'll accumulate the entries later.)
            local_catalog = fetch_catalog(rcf, stagedir)

            # we'll use sources outside of the box to render the background,
            # but we won't optimize them
            local_tasks = filter(s->in_box(s, box), local_catalog)

            catalog_offset[n] = length(local_catalog)
            task_offset[n] = length(local_tasks)
        end
    end

   # folds right, converting each field's count to its offset in `catalog`
   catalog_size = 0
    for n in 1:num_fields
        catalog_size += catalog_offset[n]
        catalog_offset[n] = catalog_size
    end

    # folds right, converting each field's count to its offset in `tasks`
    num_tasks = 0
    for n in 1:num_fields
        num_tasks += task_offset[n]
        task_offset[n] = num_tasks
    end

    images, catalog_offset, task_offset
end

function load_catalog(box, rcfs, catalog_offset, task_offset, stagedir)
    # TODO: set `total_threads`. We may not want to use all threads possible
    # on all nodes, since this operation hits disk. Maybe just 1 thread per node?
    total_threads = 3
    num_fields = length(rcfs)
    fields_per_thread = ceil(Int, num_fields / total_threads)

    catalog_size = catalog_offset[end]
    num_tasks = task_offset[end]

    #TODO: make the `catalog` a GlobalArray
    catalog = Array(Tuple{CatalogEntry, RunCamcolField}, catalog_size)

    # entries in `tasks` are indexes into `catalog`
    #TODO: make the `tasks` a GlobalArray
    tasks = Array(Int64, num_tasks)

    #let's read the catalog from disk again, rather than storing it in memory
    for m in 0:(total_threads-1)
        n0 = 1 + m * fields_per_thread
        nend = min(num_fields, (m + 1) * fields_per_thread)

        for n in n0:nend
            rcf = rcfs[n]

            local_catalog = fetch_catalog(rcf, stagedir)

            local_t = 0
            for local_s in 1:length(local_catalog)
                # `s` is the global index of this source (in the source array)
                offset = n == 1 ? 0 : catalog_offset[n - 1]
                s = offset + local_s
                # Note: each field's primary detections are stored contiguously
                # in the source array
                entry = local_catalog[local_s]
                catalog[s] = (entry, rcf)

                if in_box(entry, box)
                    local_t += 1
                    offset = n == 1 ? 0 : task_offset[n - 1]
                    t = offset + local_t
                    tasks[t] = s
                end
            end
        end
    end

    catalog, tasks
end

function optimize_sources(images, catalog, tasks,
                          catalog_offset, task_offset,
                          rcf_to_index, stagedir)
    #TODO: set `nnodes` to the number of nodes, and
    #`nthreads` to the number of threads per node.
    num_nodes = 2
    threads_per_node = 3
    total_threads = num_nodes * threads_per_node

    num_tasks = task_offset[end]

    tasks_per_thread = ceil(Int, num_tasks / total_threads)

    results = Dict[]

    for m in 0:(total_threads-1)
        t0 = 1 + m * tasks_per_thread
        tend = min(num_tasks, (m + 1) * tasks_per_thread)

        for t in t0:tend
            s = tasks[t]
            entry, primary_rcf = catalog[s]
            t_box = BoundingBox(entry.pos[1] - 1e-8, entry.pos[1] + 1e-8,
                                entry.pos[2] - 1e-8, entry.pos[2] + 1e-8)
            surrounding_rcfs = get_overlapping_fields(t_box, stagedir)
            local_images = Vector{TiledImage}[]
            local_catalog = CatalogEntry[];

            for rcf in surrounding_rcfs
                n = rcf_to_index[rcf.run, rcf.camcol, rcf.field]

                push!(local_images, images[n])

                s_a = n == 1 ? 1 : catalog_offset[n - 1] + 1
                s_b = catalog_offset[n]
                for neighbor in catalog[s_a:s_b]
                    push!(local_catalog, neighbor[1])
                end
            end

            flat_images = TiledImage[]
            for img5 in local_images, img in img5
                push!(flat_images, img)
            end

            i = findfirst(local_catalog, entry)
            neighbor_indexes = Infer.find_neighbors([i,], local_catalog, flat_images)[1]
            neighbors = local_catalog[neighbor_indexes]

            vs_opt = Infer.infer_source(flat_images, neighbors, entry)
            # TODO: write vs_opt (the results) to disk, to a global array of size num_tasks,
            # or to a local array, and then flush the array to disk later.

            push!(results, Dict(
                "thing_id"=>entry.thing_id,
                "objid"=>entry.objid,
                "ra"=>entry.pos[1],
                "dec"=>entry.pos[2],
                "vs"=>vs_opt))
        end
    end

    results
end

"""
Fit the Celeste model to sources in a given ra, dec range,
based on data from specified fields
- box: a bounding box specifying a region of sky
"""
function divide_sources_and_infer(
                box::BoundingBox,
                stagedir::String;
                timing=InferTiming(),
                outdir=".")
    # read the run-camcol-field triplets for this box
    rcfs = get_overlapping_fields(box, stagedir)

    # mulithreaded operation 1 (of 3) -- loads 25TB from disk for SDSS
    images, catalog_offset, task_offset = load_images(box, rcfs, stagedir)

    # mulithreaded operation 2 (of 3) -- loads 4TB from disk for SDSS
    catalog, tasks = load_catalog(box, rcfs, catalog_offset, task_offset, stagedir)

    rcf_to_index = invert_rcf_array(rcfs)

    # mulithreaded operation 3 (of 3) -- little disk access, cpu intensive
    results = optimize_sources(images, catalog, tasks,
                               catalog_offset, task_offset,
                               rcf_to_index, stagedir)
end