text:
using LaTeXStrings
import AlgebraicMultigrid
import DPFEHM
import IterativeSolvers
import JLD
import LinearAlgebra
import PyPlot
import Statistics

function solve_cg_amg(A, b)
	Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.smoothed_aggregation(A))
	x, ch = IterativeSolvers.cg(A, b; maxiter=10000, reltol=1e-12, log=true, Pl=Pl)
	return x, ch.data[:resnorm]
end

function solve_cholesky(A, b)
	Af = LinearAlgebra.cholesky(A)
	x = Af \ b
	return x, [sqrt(sum((A * x .- b) .^ 2) / sum(b .^ 2))]
end

function setup_problem(n)
	mins = [0, 0, 0]#meters
	maxs = [100, 100, 100]#meters
	ns = [n, n, n]
	coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid3d(mins, maxs, ns)
	Qs = zeros(size(coords, 2))
	dirichletnodes = Int[]
	dirichleths = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if coords[1, i] == mins[1]
			push!(dirichletnodes, i)
			dirichleths[i] = 1.0
		elseif coords[1, i] == maxs[1]
			push!(dirichletnodes, i)
			dirichleths[i] = 0.0
		end
	end
	Ks = ones(length(neighbors))
	args = (zeros(length(Qs) - length(dirichletnodes)), Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, ones(length(Qs)), ones(length(Qs)))
	b = -DPFEHM.groundwater_residuals(args...)
	A = DPFEHM.groundwater_h(args...)
	return A, b, (neighbors, areasoverlengths, Ks, Qs, dirichletnodes, dirichleths)
end

for i = 1:10#this loop gets compilation out of the way
	A, b, _ = setup_problem(2 ^ 2)
	solve_cholesky(A, b)
	solve_cg_amg(A, b)
end
if !isfile("porous_media_results.jld")
	ns = 2 .^ (2:7)
	num_solves = 10
	cholesky_times = zeros(length(ns), num_solves)
	cg_times = zeros(length(ns), num_solves)
	for (j, n) in enumerate(ns)
		A, b = setup_problem(n)
		these_cholesky_times = Float64[]
		for i = 1:num_solves
			@show n, i
			t = @elapsed solve_cholesky(A, b)
			@show t
			cholesky_times[j, i] = t
			t = @elapsed solve_cg_amg(A, b)
			@show t
			cg_times[j, i] = t
		end
	end
	JLD.save("porous_media_results.jld", "cholesky_times", cholesky_times, "cg_times", cg_times, "ns", ns)
end

function plotit(cholesky_times, cg_times, ns, ax)
	ax.loglog(ns, Statistics.median(cholesky_times; dims=2), ".", label="Cholesky", color="C0", ms=10, alpha=0.75)
	ax.loglog(ns, Statistics.median(cg_times; dims=2), ".", label="CG+AMG", color="C1", ms=10, alpha=0.75)
	for i = 1:size(cholesky_times, 1)
		ax.plot([ns[i], ns[i]], sort(cholesky_times[i, :])[[2, end - 1]], color="C0", alpha=0.75)
		ax.plot([ns[i], ns[i]], sort(cg_times[i, :])[[2, end - 1]], color="C1", alpha=0.75)
	end
	ax.set(xlabel="Degrees of freedom", ylabel="Time [s]")
	ax.legend()
end

cholesky_times, cg_times, ns = JLD.load("porous_media_results.jld", "cholesky_times", "cg_times", "ns")
fig, axs = PyPlot.subplots(1, 3; figsize=(12, 4))
plotit(cholesky_times, cg_times, ns .^ 3, axs[3])

function setup_problem_2d(n)
	mins = [0, 0]#meters
	maxs = [100, 100]#meters
	ns = [n, n]
	coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid2d(mins, maxs, ns, 1.0)
	Qs = zeros(size(coords, 2))
	dirichletnodes = Int[]
	dirichleths = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if coords[1, i] == mins[1]
			push!(dirichletnodes, i)
			dirichleths[i] = 1.0
		elseif coords[1, i] == maxs[1]
			push!(dirichletnodes, i)
			dirichleths[i] = 0.0
		end
	end
	Ks = ones(length(neighbors))
	args = (zeros(length(Qs) - length(dirichletnodes)), Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, ones(length(Qs)), ones(length(Qs)))
	b = -DPFEHM.groundwater_residuals(args...)
	A = DPFEHM.groundwater_h(args...)
	return A, b, (neighbors, areasoverlengths, Ks, Qs, dirichletnodes, dirichleths)
end

for i = 1:10#this loop gets compilation out of the way
	A, b, _ = setup_problem_2d(2 ^ 2)
	solve_cholesky(A, b)
	solve_cg_amg(A, b)
end
if !isfile("porous_media_results_2d.jld")
	ns = 2 .^ (3:13)
	num_solves = 10
	cholesky_times = zeros(length(ns), num_solves)
	cg_times = zeros(length(ns), num_solves)
	for (j, n) in enumerate(ns)
		A, b = setup_problem_2d(n)
		these_cholesky_times = Float64[]
		for i = 1:num_solves
			@show n, i
			t = @elapsed solve_cholesky(A, b)
			@show t
			cholesky_times[j, i] = t
			t = @elapsed solve_cg_amg(A, b)
			@show t
			cg_times[j, i] = t
		end
	end
	JLD.save("porous_media_results_2d.jld", "cholesky_times", cholesky_times, "cg_times", cg_times, "ns", ns)
end

cholesky_times, cg_times, ns = JLD.load("porous_media_results_2d.jld", "cholesky_times", "cg_times", "ns")
plotit(cholesky_times, cg_times, ns .^ 2, axs[2])

function setup_problem_1d(n)
	mins = [0]#meters
	maxs = [100]#meters
	ns = [n]
	coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid1d(mins, maxs, ns, 1.0, 1.0)
	Qs = zeros(size(coords, 2))
	dirichletnodes = Int[]
	dirichleths = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if coords[1, i] == mins[1]
			push!(dirichletnodes, i)
			dirichleths[i] = 1.0
		elseif coords[1, i] == maxs[1]
			push!(dirichletnodes, i)
			dirichleths[i] = 0.0
		end
	end
	Ks = ones(length(neighbors))
	args = (zeros(length(Qs) - length(dirichletnodes)), Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, ones(length(Qs)), ones(length(Qs)))
	b = -DPFEHM.groundwater_residuals(args...)
	A = DPFEHM.groundwater_h(args...)
	return A, b, (neighbors, areasoverlengths, Ks, Qs, dirichletnodes, dirichleths)
end

for i = 1:10#this loop gets compilation out of the way
	A, b, _ = setup_problem_1d(2 ^ 2)
	solve_cholesky(A, b)
	solve_cg_amg(A, b)
end
if !isfile("porous_media_results_1d.jld")
	ns = 10 .^ (1:8)
	num_solves = 10
	cholesky_times = zeros(length(ns), num_solves)
	cg_times = zeros(length(ns), num_solves)
	for (j, n) in enumerate(ns)
		A, b = setup_problem_1d(n)
		these_cholesky_times = Float64[]
		for i = 1:num_solves
			@show n, i
			t = @elapsed solve_cholesky(A, b)
			@show t
			cholesky_times[j, i] = t
			t = @elapsed solve_cg_amg(A, b)
			@show t
			cg_times[j, i] = t
		end
	end
	JLD.save("porous_media_results_1d.jld", "cholesky_times", cholesky_times, "cg_times", cg_times, "ns", ns)
end

cholesky_times, cg_times, ns = JLD.load("porous_media_results_1d.jld", "cholesky_times", "cg_times", "ns")
plotit(cholesky_times, cg_times, ns, axs[1])
for i = 1:3
	axs[i].set_title("$(i)D")
end
fig.tight_layout()
display(fig)
println()
fig.savefig("porous_media_comparison.pdf")
PyPlot.close(fig)


text:
sidebar = sidebar_layout(
    side_panel(
        h1("hi there"),
        slider("slider","Slider",-100,11),
        text_input("textinput", "Add a number")
    ),
    main_panel(
        text_output("foo"),
        text_output("bar"),
        h1(text_output("baz")),
        p("just some more text")
    )
)

@test typeof(sidebar) <: Tuple
@test typeof(sidebar[1]) <: Matte.UIElement

tabs = tabs_layout(
    tab_panel("Tab 1", (
        h1("hi there"),
        slider("slider","Slider",-100,11),
        text_input("textinput", "Add a number")
    )),
    tab_panel("Tab 2", (
        text_output("foo"),
        text_output("bar"),
        h1(text_output("baz")),
        p("just some more text")
    ))
)

@test typeof(tabs) <: Tuple
@test typeof(tabs[1]) <: Matte.UIElement

module TestTablayout
using Matte
const title = "Test Tab Layout"
function ui()
tabs_layout(
  tab_panel(
    "Tab 1",
    (
      h1("Content for tab 1 goes here"),
      br(),
      slider("slider1", "Number 1", 0, 100)
    )
  ),
  tab_panel(
    "Tab 2",
    sidebar_layout(
      side_panel(h1("You can even inset sidebar layouts inside tab!")),
      main_panel(plots_output("my_plot"))
    )
  ),
  tab_panel(
    "Tab 3",
    "You can have as many tabs as you like..."
  )
)
end
module Server
end

end

run_app(TestTablayout, async = true)

# Test that we can actually fetch the javascript that powers Matte
res = HTTP.request("GET", "http://localhost:8000/")
@test res.status == 200
stop_app()

custom = custom_grid_layout((
    custom_grid_row(
        custom_grid_column(
            (
                h1(text_output("baz")),
                p("just some more text")
            ),
            width = 6
        ),
        custom_grid_column(
            (
                h1("hi there"),
                slider("slider","Slider",-100,11),
                text_input("textinput", "Add a number")
            ),
            width = 6
        )
    )
))

@test typeof(custom) <: Tuple
@test typeof(custom[1]) <: Matte.UIElement

ext: jl


text:
using FountainCodes, LinearAlgebra, SparseArrays, Test

for i in 1:100
    rng = MersenneTwister(i)
    A = sprand(rng, 100, 100, 0.04)
    nonzeros(A) .= 1
    enc = zeros(size(A, 1))
    try
        src = decode(A, enc)
        @test count(isnan, src) == 0
    catch e
        if e isa RankDeficiencyException
            r = rank(A)
            @test r < size(A, 1)
        else
            rethrow()
        end
    end
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
function sortkey(key_values, a)
    indperm=sortperm(key_values)
    return a[indperm]
end

function sortkey(f::Function, a)
    sortkey(f(a), a)
end


text:
## Julia 1.1.0

# Initialize
using Plots
plotlyjs() # Backend PlotlyJS

angles = 200 * 2π
Θ = range(0, stop = angles, length = 10000)

x(n) = cos.(Θ) + (1.0/n) .* cos.(n .* Θ)
y(n) = sin.(Θ) - (1.0/n) .* sin.(n .* Θ)

v = exp(0.0625)
plot(x(v), y(v), size=(960, 960))
savefig("./png/round-spiral.png")

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
#    Problem : GROUP A
#    *********
#    The ENGVAL1 problem.
#    This problem is a sum of 2n-2 groups, n-1 of which contain 2 nonlinear
#    elements
#
#    Original SIF Source: problem 31 in
#    Ph.L. Toint,
#    "Test problems for partially separable optimization and results
#    for the routine PSPMIN",
#    Report 83/4, Department of Mathematics, FUNDP (Namur, B), 1983.
#
#    See also Buckley#172 (p. 52)
#    SIF input: Ph. Toint and N. Gould, Dec 1989.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    ENGVAL1.SIF classification OUR2-AN-V-0
#
#    N is the number of variables
#
# Daniel Henderson, 08/2021

f = (x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)-1
        γ = x[i]^2 + x[i+1]^2
        fx += γ^2 + (3.0 - 4.0x[i])
    end
    return fx
end

g! = (g, x) -> begin
    for i in firstindex(x):lastindex(x)-1
        γ = x[i]^2 + x[i+1]^2
        g[i] += 4.0γ*x[i] - 4.0
        g[i+1] += 4.0γ*x[i+1]
    end
    return g
end

fg! = (g, x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)-1
        γ = x[i]^2 + x[i+1]^2
        fx += γ^2 + (3.0 - 4.0x[i])
        g[i] += 4.0γ*x[i] - 4.0
        g[i+1] += 4.0γ*x[i+1]
    end
    return fx, g
end

init = (n::Int=5000) -> begin
	n < 2 && @warn("ENGVAL1: number of variables must be ≥ 2")
	n = max(n, 2)

    return n, 2.0*ones(n)
end

TestSet["ENGVAL1"] = UncProgram("ENGVAL1",  f, g!, fg!, init)


text:
using ReachabilityAnalysis, Plots

@taylorize function seir2!(du,u,p,t)
  S, E, I, R, α, β, γ = u

  βIS = β * (I * S)
  αE = α*E
  γI = γ*I

  du[1] = -βIS      # dS
  du[2] = βIS - αE  # dE
  du[3] = -γI + αE  # dI
  du[4] = γI        # dR

  #uncertain parameters
  local zerou = zero(u[1])
  du[5] = zerou
  du[6] = zerou
  du[7] = zerou
end

E₀ = 1e-4
u₀ = [1-E₀, E₀, 0, 0]
param_error = 0.01
α = 0.2 ± param_error
β = 1.0 ± 0.0
γ = 0.5 ± param_error
p = [α, β, γ]
X0 = IntervalBox(vcat(u₀, p));
prob = @ivp(x' = seir2!(x), dim=7, x(0) ∈ X0);

sol = solve(prob, tspan=(0.0, 200.0), alg=TMJets21a(orderT=7, orderQ=1));
solz = overapproximate(sol, Zonotope);

LazySets.set_ztol(Float64, 1e-13)
fig = plot(legend=:outerright)
plot!(fig, solz, vars=(0, 1), color=:blue, lw=0.0, lab="S")
plot!(fig, solz, vars=(0, 2), color=:green, lw=0.0, lab="E")
plot!(fig, solz, vars=(0, 3), color=:red, lw=0.0, lab="I")
plot!(solz, vars=(0, 4), color=:grey, lw=0.0, lab="R")
fig


text:
#=
    Sorting networks that process 17,..25 values are given.

    These sorting networks should minimize conditional swaps.
    There may exist one or more comparable alternative solutions;
      those may perform marginally better for some kinds of hardware.

    Functions are written with groups of internally paralleli1zable
    statements given together and sequential steps are separated.

    `swapsortr(x1::T, .., xn::T)::NTuple{n,T}`
    `swapsortr(tup::NTuple{n,T})::NTuple{n,T}`

    networks were selected using software by John Gamble
        http://pages.ripco.net/~jgamble/nw.html
        http://search.cpan.org/dist/Algorithm-Networksort/
=#

const MORE_ARGS_MIN = 17
const MORE_ARGS_MAX = 25

#=
   sort NTuples of 17..25 values
=#

for N in collect(MORE_ARGS_MIN:MORE_ARGS_MAX)
    @eval swapsortr(x::NTuple{$N, T}) where T = swapsortr(x...)
end

#=
    sorting networks for 17..25 values
    using known compact realizations
=#

#    sort 17 values with 74 maxmins in 12 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
                f::T, g::T, h::T, i::T, j::T,
                k::T, l::T, m::T, n::T, o::T,
                p::T, q::T) where {T}

    a, q = max_min(a, q)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)

    i, q = max_min(i, q)
    a, e = max_min(a, e)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    h, l = max_min(h, l)
    b, d = max_min(b, d)
    n, p = max_min(n, p)

    i, m = max_min(i, m)
    e, q = max_min(e, q)
    a, c = max_min(a, c)
    f, h = max_min(f, h)
    j, l = max_min(j, l)

    e, i = max_min(e, i)
    m, q = max_min(m, q)
    d, j = max_min(d, j)
    h, n = max_min(h, n)
    a, b = max_min(a, b)

    e, g = max_min(e, g)
    i, k = max_min(i, k)
    m, o = max_min(m, o)
    c, q = max_min(c, q)
    d, f = max_min(d, f)
    h, j = max_min(h, j)
    l, n = max_min(l, n)

    c, i = max_min(c, i)
    g, m = max_min(g, m)
    k, q = max_min(k, q)

    c, e = max_min(c, e)
    g, i = max_min(g, i)
    k, m = max_min(k, m)
    o, q = max_min(o, q)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    b, q = max_min(b, q)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    h, o = max_min(h, o)
    j, q = max_min(j, q)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q
end

#    sort 18 values with 82 maxmins in 13 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, g = max_min(c, g)
    d, h = max_min(d, h)
    k, o = max_min(k, o)
    l, p = max_min(l, p)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    g, k = max_min(g, k)
    h, l = max_min(h, l)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    e, q = max_min(e, q)
    f, r = max_min(f, r)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    m, q = max_min(m, q)
    n, r = max_min(n, r)
    a, b = max_min(a, b)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    c, q = max_min(c, q)
    d, r = max_min(d, r)

    c, i = max_min(c, i)
    d, j = max_min(d, j)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    k, q = max_min(k, q)
    l, r = max_min(l, r)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    g, i = max_min(g, i)
    h, j = max_min(h, j)
    k, m = max_min(k, m)
    l, n = max_min(l, n)
    o, q = max_min(o, q)
    p, r = max_min(p, r)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)

    b, q = max_min(b, q)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    h, o = max_min(h, o)

    b, i = max_min(b, i)
    j, q = max_min(j, q)
    d, g = max_min(d, g)
    h, k = max_min(h, k)
    l, o = max_min(l, o)

    b, e = max_min(b, e)
    f, i = max_min(f, i)
    j, m = max_min(j, m)
    n, q = max_min(n, q)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r
end

#    sort 19 values with 91 maxmins in 14 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, h = max_min(d, h)
    l, p = max_min(l, p)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    h, l = max_min(h, l)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    e, q = max_min(e, q)
    f, r = max_min(f, r)
    g, s = max_min(g, s)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    m, q = max_min(m, q)
    n, r = max_min(n, r)
    o, s = max_min(o, s)
    a, b = max_min(a, b)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    q, s = max_min(q, s)
    d, r = max_min(d, r)

    c, q = max_min(c, q)
    d, j = max_min(d, j)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    l, r = max_min(l, r)

    c, i = max_min(c, i)
    k, q = max_min(k, q)
    d, f = max_min(d, f)
    h, j = max_min(h, j)
    l, n = max_min(l, n)
    p, r = max_min(p, r)

    c, e = max_min(c, e)
    g, i = max_min(g, i)
    k, m = max_min(k, m)
    o, q = max_min(o, q)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, m = max_min(f, m)
    h, o = max_min(h, o)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    j, q = max_min(j, q)
    l, s = max_min(l, s)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s
end

#    sort 20 values with 97 maxmins in 14 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, t = max_min(d, t)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    l, t = max_min(l, t)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)
    e, q = max_min(e, q)
    f, r = max_min(f, r)
    g, s = max_min(g, s)
    h, t = max_min(h, t)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    h, l = max_min(h, l)
    m, q = max_min(m, q)
    n, r = max_min(n, r)
    o, s = max_min(o, s)
    p, t = max_min(p, t)
    a, b = max_min(a, b)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    q, s = max_min(q, s)
    r, t = max_min(r, t)

    c, q = max_min(c, q)
    d, r = max_min(d, r)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    s, t = max_min(s, t)

    c, i = max_min(c, i)
    d, j = max_min(d, j)
    k, q = max_min(k, q)
    l, r = max_min(l, r)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    g, i = max_min(g, i)
    h, j = max_min(h, j)
    k, m = max_min(k, m)
    l, n = max_min(l, n)
    o, q = max_min(o, q)
    p, r = max_min(p, r)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, m = max_min(f, m)
    h, o = max_min(h, o)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    j, q = max_min(j, q)
    l, s = max_min(l, s)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t
end

#    sort 21 values with 107 maxmins in 15 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, t = max_min(d, t)
    e, u = max_min(e, u)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    l, t = max_min(l, t)
    m, u = max_min(m, u)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)
    q, u = max_min(q, u)
    f, r = max_min(f, r)
    g, s = max_min(g, s)
    h, t = max_min(h, t)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, q = max_min(e, q)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    h, l = max_min(h, l)
    n, r = max_min(n, r)
    o, s = max_min(o, s)
    p, t = max_min(p, t)
    a, b = max_min(a, b)

    e, i = max_min(e, i)
    m, q = max_min(m, q)
    f, h = max_min(f, h)
    j, l = max_min(j, l)
    n, p = max_min(n, p)
    r, t = max_min(r, t)

    e, g = max_min(e, g)
    i, k = max_min(i, k)
    m, o = max_min(m, o)
    q, s = max_min(q, s)
    d, r = max_min(d, r)
    h, n = max_min(h, n)

    c, q = max_min(c, q)
    g, u = max_min(g, u)
    d, j = max_min(d, j)
    l, r = max_min(l, r)

    c, i = max_min(c, i)
    g, m = max_min(g, m)
    k, q = max_min(k, q)
    o, u = max_min(o, u)
    d, f = max_min(d, f)
    h, j = max_min(h, j)
    l, n = max_min(l, n)
    p, r = max_min(p, r)

    c, e = max_min(c, e)
    g, i = max_min(g, i)
    k, m = max_min(k, m)
    o, q = max_min(o, q)
    s, u = max_min(s, u)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)
    s, t = max_min(s, t)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, u = max_min(f, u)
    h, o = max_min(h, o)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    j, q = max_min(j, q)
    l, s = max_min(l, s)
    n, u = max_min(n, u)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)
    r, u = max_min(r, u)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)
    t, u = max_min(t, u)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u
end

#    sort 22 values with 114 maxmins in 15 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, t = max_min(d, t)
    e, u = max_min(e, u)
    f, v = max_min(f, v)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    l, t = max_min(l, t)
    m, u = max_min(m, u)
    n, v = max_min(n, v)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)
    q, u = max_min(q, u)
    r, v = max_min(r, v)
    g, s = max_min(g, s)
    h, t = max_min(h, t)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, q = max_min(e, q)
    f, r = max_min(f, r)
    g, k = max_min(g, k)
    h, l = max_min(h, l)
    o, s = max_min(o, s)
    p, t = max_min(p, t)
    a, b = max_min(a, b)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    m, q = max_min(m, q)
    n, r = max_min(n, r)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    q, s = max_min(q, s)
    r, t = max_min(r, t)

    c, q = max_min(c, q)
    d, r = max_min(d, r)
    g, u = max_min(g, u)
    h, v = max_min(h, v)

    c, i = max_min(c, i)
    d, j = max_min(d, j)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    k, q = max_min(k, q)
    l, r = max_min(l, r)
    o, u = max_min(o, u)
    p, v = max_min(p, v)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    g, i = max_min(g, i)
    h, j = max_min(h, j)
    k, m = max_min(k, m)
    l, n = max_min(l, n)
    o, q = max_min(o, q)
    p, r = max_min(p, r)
    s, u = max_min(s, u)
    t, v = max_min(t, v)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)
    s, t = max_min(s, t)
    u, v = max_min(u, v)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, u = max_min(f, u)
    h, o = max_min(h, o)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    j, q = max_min(j, q)
    l, s = max_min(l, s)
    n, u = max_min(n, u)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)
    r, u = max_min(r, u)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)
    t, u = max_min(t, u)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v
end

#    sort 23 values with 122 maxmins in 15 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T, w::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, t = max_min(d, t)
    e, u = max_min(e, u)
    f, v = max_min(f, v)
    g, w = max_min(g, w)
    h, p = max_min(h, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    l, t = max_min(l, t)
    m, u = max_min(m, u)
    n, v = max_min(n, v)
    o, w = max_min(o, w)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)
    q, u = max_min(q, u)
    r, v = max_min(r, v)
    s, w = max_min(s, w)
    h, t = max_min(h, t)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, q = max_min(e, q)
    f, r = max_min(f, r)
    g, s = max_min(g, s)
    h, l = max_min(h, l)
    p, t = max_min(p, t)
    u, w = max_min(u, w)
    a, b = max_min(a, b)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    m, q = max_min(m, q)
    n, r = max_min(n, r)
    o, s = max_min(o, s)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    q, s = max_min(q, s)
    r, t = max_min(r, t)

    c, q = max_min(c, q)
    d, r = max_min(d, r)
    g, u = max_min(g, u)
    h, v = max_min(h, v)

    c, i = max_min(c, i)
    d, j = max_min(d, j)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    k, q = max_min(k, q)
    l, r = max_min(l, r)
    o, u = max_min(o, u)
    p, v = max_min(p, v)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    g, i = max_min(g, i)
    h, j = max_min(h, j)
    k, m = max_min(k, m)
    l, n = max_min(l, n)
    o, q = max_min(o, q)
    p, r = max_min(p, r)
    s, u = max_min(s, u)
    t, v = max_min(t, v)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)
    s, t = max_min(s, t)
    u, v = max_min(u, v)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, u = max_min(f, u)
    h, w = max_min(h, w)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    h, o = max_min(h, o)
    j, q = max_min(j, q)
    l, s = max_min(l, s)
    n, u = max_min(n, u)
    p, w = max_min(p, w)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)
    r, u = max_min(r, u)
    t, w = max_min(t, w)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)
    t, u = max_min(t, u)
    v, w = max_min(v, w)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w
end

#    sort 24 values with 127 maxmins in 15 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T, w::T, x::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, t = max_min(d, t)
    e, u = max_min(e, u)
    f, v = max_min(f, v)
    g, w = max_min(g, w)
    h, x = max_min(h, x)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    l, t = max_min(l, t)
    m, u = max_min(m, u)
    n, v = max_min(n, v)
    o, w = max_min(o, w)
    p, x = max_min(p, x)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)
    q, u = max_min(q, u)
    r, v = max_min(r, v)
    s, w = max_min(s, w)
    t, x = max_min(t, x)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, q = max_min(e, q)
    f, r = max_min(f, r)
    g, s = max_min(g, s)
    h, t = max_min(h, t)
    u, w = max_min(u, w)
    v, x = max_min(v, x)
    a, b = max_min(a, b)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    h, l = max_min(h, l)
    m, q = max_min(m, q)
    n, r = max_min(n, r)
    o, s = max_min(o, s)
    p, t = max_min(p, t)
    w, x = max_min(w, x)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    q, s = max_min(q, s)
    r, t = max_min(r, t)

    c, q = max_min(c, q)
    d, r = max_min(d, r)
    g, u = max_min(g, u)
    h, v = max_min(h, v)

    c, i = max_min(c, i)
    d, j = max_min(d, j)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    k, q = max_min(k, q)
    l, r = max_min(l, r)
    o, u = max_min(o, u)
    p, v = max_min(p, v)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    g, i = max_min(g, i)
    h, j = max_min(h, j)
    k, m = max_min(k, m)
    l, n = max_min(l, n)
    o, q = max_min(o, q)
    p, r = max_min(p, r)
    s, u = max_min(s, u)
    t, v = max_min(t, v)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)
    s, t = max_min(s, t)
    u, v = max_min(u, v)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, u = max_min(f, u)
    h, w = max_min(h, w)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    h, o = max_min(h, o)
    j, q = max_min(j, q)
    l, s = max_min(l, s)
    n, u = max_min(n, u)
    p, w = max_min(p, w)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)
    r, u = max_min(r, u)
    t, w = max_min(t, w)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)
    t, u = max_min(t, u)
    v, w = max_min(v, w)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x
end

#    sort 25 values with 138 maxmins in 15 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T, w::T, x::T, y::T) where {T}

    a, q = max_min(a, q)
    b, r = max_min(b, r)
    c, s = max_min(c, s)
    d, t = max_min(d, t)
    e, u = max_min(e, u)
    f, v = max_min(f, v)
    g, w = max_min(g, w)
    h, x = max_min(h, x)
    i, y = max_min(i, y)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)
    q, y = max_min(q, y)

    i, q = max_min(i, q)
    j, r = max_min(j, r)
    k, s = max_min(k, s)
    l, t = max_min(l, t)
    m, u = max_min(m, u)
    n, v = max_min(n, v)
    o, w = max_min(o, w)
    p, x = max_min(p, x)
    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    i, m = max_min(i, m)
    j, n = max_min(j, n)
    k, o = max_min(k, o)
    l, p = max_min(l, p)
    q, u = max_min(q, u)
    r, v = max_min(r, v)
    s, w = max_min(s, w)
    t, x = max_min(t, x)
    a, c = max_min(a, c)
    b, d = max_min(b, d)

    e, q = max_min(e, q)
    f, r = max_min(f, r)
    g, s = max_min(g, s)
    h, t = max_min(h, t)
    m, y = max_min(m, y)
    v, x = max_min(v, x)
    a, b = max_min(a, b)

    e, i = max_min(e, i)
    f, j = max_min(f, j)
    g, k = max_min(g, k)
    h, l = max_min(h, l)
    m, q = max_min(m, q)
    n, r = max_min(n, r)
    o, s = max_min(o, s)
    p, t = max_min(p, t)
    u, y = max_min(u, y)

    e, g = max_min(e, g)
    f, h = max_min(f, h)
    i, k = max_min(i, k)
    j, l = max_min(j, l)
    m, o = max_min(m, o)
    n, p = max_min(n, p)
    q, s = max_min(q, s)
    r, t = max_min(r, t)
    u, w = max_min(u, w)

    c, q = max_min(c, q)
    d, r = max_min(d, r)
    g, u = max_min(g, u)
    h, v = max_min(h, v)
    k, y = max_min(k, y)

    c, i = max_min(c, i)
    d, j = max_min(d, j)
    g, m = max_min(g, m)
    h, n = max_min(h, n)
    k, q = max_min(k, q)
    l, r = max_min(l, r)
    o, u = max_min(o, u)
    p, v = max_min(p, v)
    s, y = max_min(s, y)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    g, i = max_min(g, i)
    h, j = max_min(h, j)
    k, m = max_min(k, m)
    l, n = max_min(l, n)
    o, q = max_min(o, q)
    p, r = max_min(p, r)
    s, u = max_min(s, u)
    t, v = max_min(t, v)
    w, y = max_min(w, y)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)
    q, r = max_min(q, r)
    s, t = max_min(s, t)
    u, v = max_min(u, v)
    w, x = max_min(w, x)

    b, q = max_min(b, q)
    d, s = max_min(d, s)
    f, u = max_min(f, u)
    h, w = max_min(h, w)
    j, y = max_min(j, y)

    b, i = max_min(b, i)
    d, k = max_min(d, k)
    f, m = max_min(f, m)
    h, o = max_min(h, o)
    j, q = max_min(j, q)
    l, s = max_min(l, s)
    n, u = max_min(n, u)
    p, w = max_min(p, w)
    r, y = max_min(r, y)

    b, e = max_min(b, e)
    d, g = max_min(d, g)
    f, i = max_min(f, i)
    h, k = max_min(h, k)
    j, m = max_min(j, m)
    l, o = max_min(l, o)
    n, q = max_min(n, q)
    p, s = max_min(p, s)
    r, u = max_min(r, u)
    t, w = max_min(t, w)
    v, y = max_min(v, y)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)
    n, o = max_min(n, o)
    p, q = max_min(p, q)
    r, s = max_min(r, s)
    t, u = max_min(t, u)
    v, w = max_min(v, w)
    x, y = max_min(x, y)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y
end


text:
module GraphicalLinearAlgebra
export LinearFunctions, FreeLinearFunctions, LinearRelations,
  FreeLinearRelations, LinearMapDom, LinearMap,
  LinearOpDom, LinearOperator,
  Ob, Hom, dom, codom, compose, ⋅, ∘, id, oplus, ⊕, mzero, braid,
  dagger, dunit, docunit, mcopy, Δ, delete, ◊, mmerge, ∇, create, □,
  pplus, zero, coplus, cozero, plus, +, meet, top, join, bottom,
  scalar, antipode, antipode, adjoint, evaluate

import Base: +
using AutoHashEquals
using LinearMaps
import LinearMaps: adjoint
const LMs = LinearMaps
using LinearOperators
import LinearOperators:
  adjoint, opEye, opExtension, opRestriction, opZeros

using ...Catlab, ...Doctrines
import ...Doctrines:
  Ob, Hom, dom, codom, compose, ⋅, ∘, id, oplus, ⊕, mzero, braid,
  dagger, dunit, dcounit, mcopy, Δ, delete, ◊, mmerge, ∇, create, □,
  plus, zero, coplus, cozero, meet, top, join, bottom
using ...Programs
import ...Programs: evaluate_hom

# Doctrines
###########

""" Doctrine of *linear functions*, aka linear maps

Functional fragment of graphical linear algebra.
"""
@theory AdditiveSymmetricMonoidalCategory(Ob,Hom) => LinearFunctions(Ob,Hom) begin
  # Copying and deleting maps.
  mcopy(A::Ob)::(A → (A ⊕ A))
  @op (Δ) := mcopy
  delete(A::Ob)::(A → mzero())
  @op (◊) := delete

  # Addition and zero maps.
  plus(A::Ob)::((A ⊕ A) → A)
  @op (+) := plus
  zero(A::Ob)::(mzero() → A)

  plus(f::(A → B), g::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)
  adjoint(f::(A → B))::(B → A) ⊣ (A::Ob, B::Ob)

  scalar(A::Ob, c::Number)::(A → A)
  antipode(A::Ob)::(A → A)

  # Axioms
  antipode(A) == scalar(A, -1) ⊣ (A::Ob)

  Δ(A) == Δ(A) ⋅ σ(A,A) ⊣ (A::Ob)
  Δ(A) ⋅ (Δ(A) ⊕ id(A)) == Δ(A) ⋅ (id(A) ⊕ Δ(A)) ⊣ (A::Ob)
  Δ(A) ⋅ (◊(A) ⊕ id(A)) == id(A) ⊣ (A::Ob)
  plus(A) == σ(A,A) ⋅ plus(A) ⊣ (A::Ob)
  (plus(A) ⊕ id(A)) ⋅ plus(A) == (id(A) ⊕ plus(A)) ⋅ plus(A) ⊣ (A::Ob)
  (zero(A) ⊕ id(A)) ⋅ plus(A) == id(A) ⊣ (A::Ob)
  plus(A) ⋅ Δ(A) == ((Δ(A) ⊕ Δ(A)) ⋅ (id(A) ⊕ (σ(A, A) ⊕ id(A)))) ⋅ (plus(A) ⊕ plus(A)) ⊣ (A::Ob)
  plus(A) ⋅ ◊(A) == ◊(A) ⊕ ◊(A) ⊣ (A::Ob)
  zero(A) ⋅ Δ(A) == zero(A) ⊕ zero(A) ⊣ (A::Ob)
  zero(A) ⋅ ◊(A) == id(mzero()) ⊣ (A::Ob)
  scalar(A, a) ⋅ scalar(A, b) == scalar(A, a*b) ⊣ (A::Ob, a::Number, b::Number)
  scalar(A, 1) == id(A) ⊣ (A::Ob)
  scalar(A, a) ⋅ Δ(A) == Δ(A) ⋅ (scalar(A, a) ⊕ scalar(A, a)) ⊣ (A::Ob, a::Number)
  scalar(A, a) ⋅ ◊(A) == ◊(A) ⊣ (A::Ob, a::Number)
  (Δ(A) ⋅ (scalar(A, a) ⊕ scalar(A, b))) ⋅ plus(A) == scalar(A, a+b) ⊣ (A::Ob, a::Number, b::Number)
  scalar(A, 0) == ◊(A) ⋅ zero(A) ⊣ (A::Ob)
  zero(A) ⋅ scalar(A, a) == zero(A) ⊣ (A::Ob, a::Number)

  plus(A) ⋅ f == (f ⊕ f) ⋅ plus(B) ⊣ (A::Ob, B::Ob, f::(A → B))
  scalar(A, c) ⋅ f == f ⋅ scalar(B, c) ⊣ (A::Ob, B::Ob, c::Number, f::(A → B))
end

@syntax FreeLinearFunctions(ObExpr,HomExpr) LinearFunctions begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = new(f,g; strict=true) # No normalization!
end

""" Doctrine of *linear relations*

The full relational language of graphical linear algebra. This is an abelian
bicategory of relations (`AbelianBicategoryRelations`), written additively.
"""
@signature LinearFunctions(Ob,Hom) => LinearRelations(Ob,Hom) begin
  # Dagger category.
  dagger(f::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)

  # Self-dual compact closed category.
  dunit(A::Ob)::(mzero() → (A ⊕ A))
  dcounit(A::Ob)::((A ⊕ A) → mzero())

  # Merging and creating relations (converses of copying and deleting maps).
  mmerge(A::Ob)::((A ⊕ A) → A)
  @op (∇) := mmerge
  create(A::Ob)::(mzero() → A)
  @op (□) := create

  # Co-addition and co-zero relations (converses of addition and zero maps)
  coplus(A::Ob)::(A → (A ⊕ A))
  cozero(A::Ob)::(A → mzero())

  # Lattice of linear relations.
  meet(f::(A → B), g::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)
  top(A::Ob, B::Ob)::(A → B)
  join(f::(A → B), g::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)
  bottom(A::Ob, B::Ob)::(A → B)
end

@syntax FreeLinearRelations(ObExpr,HomExpr) LinearRelations begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = new(f,g; strict=true) # No normalization!
end

# Evaluation
############

# LinearMaps instance
#--------------------

@auto_hash_equals struct LinearMapDom
  N::Int
end

@instance LinearFunctions(LinearMapDom, LinearMap) begin
  @import adjoint, +

  dom(f::LinearMap) = LinearMapDom(size(f,2))
  codom(f::LinearMap) = LinearMapDom(size(f,1))

  compose(f::LinearMap, g::LinearMap) = g*f
  id(V::LinearMapDom) = LMs.UniformScalingMap(1, V.N)

  oplus(V::LinearMapDom, W::LinearMapDom) = LinearMapDom(V.N + W.N)
  oplus(f::LinearMap, g::LinearMap) = LMs.BlockDiagonalMap(f, g)
  mzero(::Type{LinearMapDom}) = LinearMapDom(0)
  braid(V::LinearMapDom, W::LinearMapDom) =
    LinearMap(braid_lm(V.N), braid_lm(W.N), W.N+V.N, V.N+W.N)

  mcopy(V::LinearMapDom) = LinearMap(mcopy_lm, plus_lm, 2*V.N, V.N)
  delete(V::LinearMapDom) = LinearMap(delete_lm, zero_lm(V.N), 0, V.N)
  plus(V::LinearMapDom) = LinearMap(plus_lm, mcopy_lm, V.N, 2*V.N)
  zero(V::LinearMapDom) = LinearMap(zero_lm(V.N), delete_lm, V.N, 0)

  plus(f::LinearMap, g::LinearMap) = f+g
  scalar(V::LinearMapDom, c::Number) = LMs.UniformScalingMap(c, V.N)
  antipode(V::LinearMapDom) = LMs.UniformScalingMap(-1, V.N)
end

braid_lm(n::Int) = x::AbstractVector -> vcat(x[n+1:end], x[1:n])
mcopy_lm(x::AbstractVector) = vcat(x, x)
delete_lm(x::AbstractVector) = eltype(x)[]
plus_lm(x::AbstractVector) = begin
  n = length(x) ÷ 2
  x[1:n] + x[n+1:end]
end
zero_lm(n::Int) = x::AbstractVector -> zeros(eltype(x), n)

# LinearOps instance
#-------------------

@auto_hash_equals struct LinearOpDom
  N::Int
end

@instance LinearFunctions(LinearOpDom, LinearOperator) begin
  @import adjoint, +

  dom(f::LinearOperator) = LinearOpDom(size(f,2))
  codom(f::LinearOperator) = LinearOpDom(size(f,1))

  compose(f::LinearOperator, g::LinearOperator) = g*f
  id(V::LinearOpDom) = opEye(V.N)

  oplus(V::LinearOpDom, W::LinearOpDom) = LinearOpDom(V.N + W.N)
  oplus(f::LinearOperator, g::LinearOperator) = begin
    dom_total   = size(f,2) + size(g,2)
    codom_total = size(f,1) + size(g,1)
    dom_f   = 1:size(f,2)
    codom_f = 1:size(f,1)
    dom_g   = (size(f,2)+1):dom_total
    codom_g = (size(f,1)+1):codom_total
    fOp = opExtension(codom_f, codom_total)*f*opRestriction(dom_f,dom_total)
    gOp = opExtension(codom_g, codom_total)*g*opRestriction(dom_g,dom_total)
    fOp + gOp
  end
  mzero(::Type{LinearOpDom}) = LinearOpDom(0)
  braid(V::LinearOpDom, W::LinearOpDom) =
    opExtension(1:W.N, V.N+W.N) * opRestriction((V.N+1):(V.N+W.N),V.N+W.N) +
    opExtension((W.N+1):(V.N+W.N), V.N+W.N) * opRestriction(1:V.N,V.N+W.N)
  mcopy(V::LinearOpDom) =
    opExtension(1:V.N, 2*V.N)+opExtension((V.N+1):(2*V.N), 2*V.N)
  delete(V::LinearOpDom) = opZeros(0, V.N)
  plus(V::LinearOpDom) =
    opRestriction(1:V.N, 2*V.N)+opRestriction((V.N+1):(2*V.N), 2*V.N)
  zero(V::LinearOpDom) = opZeros(V.N, 0)

  plus(f::LinearOperator, g::LinearOperator) = f+g
  scalar(V::LinearOpDom, c::Number) = opEye(typeof(c),V.N)*c
  antipode(V::LinearOpDom) = scalar(V,-1)
end

# Catlab evaluate
#----------------

function evaluate_hom(f::FreeLinearFunctions.Hom{:generator}, xs::Vector;
                      generators::AbstractDict=Dict())
  M = generators[f]
  x = reduce(vcat, xs; init=eltype(M)[])
  [ M*x ]
end

function evaluate_hom(f::FreeLinearFunctions.Hom{:plus}, xs::Vector; kw...)
  if first(f) isa ObExpr
    # Addition map.
    [ reduce(+, xs) ]
  else
    # Sum of linear maps.
    mapreduce(+, args(f)) do g
      evaluate_hom(g, xs; kw...)
    end
  end
end
function evaluate_hom(f::FreeLinearFunctions.Hom{:zero}, xs::Vector;
                      generators::AbstractDict=Dict())
  map(collect(codom(f))) do A
    dims = generators[A]
    zeros(dims...)
  end
end

evaluate_hom(f::FreeLinearFunctions.Hom{:scalar}, xs::Vector; kw...) = last(f) .* xs
evaluate_hom(f::FreeLinearFunctions.Hom{:antipode}, xs::Vector; kw...) = -1 .* xs

end


text:
module Day04

using AdventOfCode2019

function day04(input::String = readInput(joinpath(@__DIR__, "..", "data", "day04.txt")))
    low, high = parse.(Int, split(input, "-"))
    !((low >= 100_000) && (high <= 999_999)) && throw(AssertionError("Invalid input range"))

    digs = Array{UInt8,1}(undef, 6)
    reverse!(digits!(digs, low))
    first_valid!(digs)

    count1, count2 = 0, 0
    hdigs = Array{UInt8,1}(undef, 6)
    reverse!(digits!(hdigs, high))
    while true
        if has_pair(digs)
            count1 += 1
            if has_exact_pair(digs)
                count2 += 1
            end
        end
        !next!(digs, hdigs) && break
    end
    return [count1, count2]
end

"""
    has_pair(digits::Array{UInt8,1})

Return `true` if the sequence `digits` contains a pair,
otherwise `false`.
"""
function has_pair(digits::Array{UInt8,1})
    for i = 1:length(digits)-1
        digits[i] == digits[i+1] && return true
    end
    return false
end

"""
    has_exact_pair(digits::Array{UInt8,1})

Return `true` if the sequense `digits` contains an exact pair,
otherwise `false`.
"""
function has_exact_pair(digits::Array{UInt8,1})
    digits[1] == digits[2] != digits[3] && return true
    for i = 2:length(digits)-2
        digits[i-1] != digits[i] == digits[i+1] != digits[i+2] && return true
    end
    digits[end] == digits[end-1] != digits[end-2] && return true
    return false
end

"""
    first_valid!(digits::Array{UInt8,1})

Compute the first valid non-decreasing array of digits.
"""
function first_valid!(digits::Array{UInt8,1})
    for i = 1:length(digits)-1
        if digits[i] > digits[i+1]
            for j = i+1:length(digits)
                digits[j] = digits[i]
            end
            break
        end
    end
end

"""
    next!(digits::Array{UInt8,1}, high::Array{UInt8})

Compute the next non-decreasing array of digits. Return `true`
if the value represented by that next array of digits is less
or equal the value represented by the upper bound `high`.
"""
function next!(digits::Array{UInt8,1}, high::Array{UInt8,1})
    i = length(digits)
    while digits[i] == 9
        i -= 1
        if i == 0
            return false
        end
    end
    digits[i] += 1
    for j = i+1:length(digits)
        digits[j] = digits[i]
    end
    for i=1:length(digits)-1
        if digits[i] < high[i]
            return true
        elseif digits[i] > high[i]
            return false
        end
    end
    return digits[end] <= high[end]
end

end # module


text:
a = 1

while a < 10
a = a + 1
end

println(a)


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

import ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq:
    OrdinaryDiffEq,
    ODEProblem,
    ODEFunction,
    solve,
    SSPRK33,
    Rosenbrock23,
    ImplicitEuler

using Logging: global_logger
using TerminalLoggers: TerminalLogger
using LinearAlgebra
global_logger(TerminalLogger())

const FT = Float64

# https://github.com/CliMA/CLIMAParameters.jl/blob/master/src/Planet/planet_parameters.jl#L5
const MSLP = 1e5 # mean sea level pressure
const grav = 9.8 # gravitational constant
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const γ = 1.4 # heat capacity ratio
const C_p = R_d * γ / (γ - 1) # heat capacity at constant pressure
const C_v = R_d / (γ - 1) # heat capacit at constant volume
const R_m = R_d # moist R, assumed to be dry

domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(0.0),
    Geometry.ZPoint{FT}(30e3),
    boundary_tags = (:bottom, :top),
)
#mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
mesh = Meshes.IntervalMesh(domain; nelems = 30)

cspace = Spaces.CenterFiniteDifferenceSpace(mesh)
fspace = Spaces.FaceFiniteDifferenceSpace(cspace)

# https://github.com/CliMA/Thermodynamics.jl/blob/main/src/TemperatureProfiles.jl#L115-L155
# https://clima.github.io/Thermodynamics.jl/dev/TemperatureProfiles/#DecayingTemperatureProfile
function decaying_temperature_profile(z; T_virt_surf = 280.0, T_min_ref = 230.0)
    # Scale height for surface temperature
    H_sfc = R_d * T_virt_surf / grav
    H_t = H_sfc

    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv = T_virt_surf - T_min_ref
    Tv = T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / T_virt_surf
    p =
        MSLP * exp(
            (
                -H_t *
                (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
            ) / (H_sfc * (1 - ΔTv′^2)),
        )
    ρ = p / (R_d * Tv)
    ρθ = ρ * Tv * (MSLP / p)^(R_m / C_p)
    return (ρ = ρ, ρθ = ρθ)
end

Π(ρθ) = C_p * (R_d * ρθ / MSLP)^(R_m / C_v)
Φ(z) = grav * z

function discrete_hydrostatic_balance!(
    ρ,
    w,
    ρθ,
    Δz::Float64,
    _grav::Float64,
    Π::Function,
)
    # compute θ such that
    #   I(θ)[i+1/2] = -g / ∂f(Π(ρθ))
    # discretely, then set
    #   ρ = ρθ/θ
    for i in 1:(length(ρ) - 1)
        #  ρ[i+1] = ρθ[i+1]/(-2Δz*_grav/(Π(ρθ[i+1]) - Π(ρθ[i])) - ρθ[i]/ρ[i])
        ρ[i + 1] =
            ρθ[i + 1] /
            (-2 * _grav / ((Π(ρθ[i + 1]) - Π(ρθ[i])) / Δz) - ρθ[i] / ρ[i])

        ρ[i + 1] =
            ρθ[i + 1] /
            (1 / ((-2 * _grav) * (Π(ρθ[i + 1]) - Π(ρθ[i]))Δz) - ρθ[i] / ρ[i])

        ∂Π∂z = (Π(ρθ[i + 1]) - Π(ρθ[i])) / Δz
    end
end

zc = Fields.coordinate_field(cspace)
Yc = decaying_temperature_profile.(zc.z)
w = Geometry.Cartesian3Vector.(zeros(FT, fspace))
zf = parent(Fields.coordinate_field(fspace).z)
Δz = zf[2:end] - zf[1:(end - 1)]
Y_init = copy(Yc)
w_init = copy(w)
Y = Fields.FieldVector(Yc = Yc, w = w)

function tendency!(dY, Y, _, t)
    Yc = Y.Yc
    w = Y.w
    dYc = dY.Yc
    dw = dY.w

    If = Operators.InterpolateC2F()
    ∂ = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
        top = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
    )
    ∂f = Operators.GradientC2F()
    B = Operators.SetBoundaryOperator(
        bottom = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
        top = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
    )

    @. dYc.ρ = -(∂(w * If(Yc.ρ)))
    @. dYc.ρθ = -(∂(w * If(Yc.ρθ)))
    @. dw = B(
        Geometry.CartesianVector(
            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z)),
        ),
    )
    return dY
end

dY = tendency!(similar(Y), Y, nothing, 0.0)

function jacobian!(J, Y, p, t)
    # copyto!(J, LinearAlgebra.I)
    # @info length(Y)
    # @info Y[1]

    J .= 0.0

    @info "Jacobian computation!!!!"
    # N cells
    N = div(length(Y) - 1, 3)

    ρ, ρθ, w = Y[1:N], Y[(N + 1):(2N)], Y[(2N + 1):(3N + 1)]

    # construct cell center
    ρh = [ρ[1]; (ρ[1:(N - 1)] + ρ[2:N]) / 2.0; ρ[N]]
    ρθh = [ρθ[1]; (ρθ[1:(N - 1)] + ρθ[2:N]) / 2.0; ρθ[N]]

    Πc = Π.(ρθ)
    Πh = [NaN64; (Πc[1:(N - 1)] + Πc[2:N]) / 2.0; NaN64]
    Δzh = [NaN64; (Δz[1:(N - 1)] + Δz[2:N]) / 2.0; NaN64]

    for i in 1:N
        J[i, i + 2N] = ρh[i] / Δz[i]
        J[i, i + 2N + 1] = -ρh[i + 1] / Δz[i]
    end

    for i in 1:N
        J[i + N, i + 2N] = ρθh[i] / Δz[i]
        J[i + N, i + 2N + 1] = -ρθh[i + 1] / Δz[i]
    end

    # 0 for i = 1, N+1
    for i in 2:N
        J[i + 2N, (i - 1)] = -grav / (2 * ρh[i])
        J[i + 2N, (i - 1) + 1] = -grav / (2 * ρh[i])

        J[i + 2N, (i - 1) + N] = (γ - 1) * Πh[i] ./ (ρh[i] * Δzh[i])
        J[i + 2N, (i - 1) + 1 + N] = -(γ - 1) * Πh[i] ./ (ρh[i] * Δzh[i])
    end

    return J

    # D_ρ = diagm(0=>-ρh/Δz, -1=>ρh/Δz)[1:N, 1:N-1]
    # D_Θ = diagm(0=>-ρθh/Δz, -1=>ρθh/Δz)[1:N, 1:N-1]
    # G_W = (γ - 1) * diagm(0=>Πh./ρh/Δz, 1=>-Πh./ρh/Δz)[1:N-1, 1:N]
    # A_W = diagm(0=>-ones(N-1)./ρh/2, 1=>-ones(N-1)./ρh/2)[1:N-1, 1:N]

    # P = ([zeros(N,N)     zeros(N,N)      D_ρ;
    #       zeros(N,N)     zeros(N,N)      D_Θ
    #       A_W*_grav        G_W          zeros(N+1,N+1)])

end

function linsolve!(::Type{Val{:init}}, f, u0; kwargs...)
    function _linsolve!(x, A, b, update_matrix = false; kwargs...)
        #     _A = RecursiveFactorization.lu(A)
        #   ldiv!(x,_A,b)
        #   @info "norm: ", norm(x), norm(A), norm(b)
        x = copy(b)
        N = div(length(x) - 1, 3)
        J = A[(2N + 1):(3N + 1), (2N + 1):(3N + 1)]
        J +=
            -A[(2N + 1):(3N + 1), 1:N] *
            (Diagonal(A[1:N, 1:N]) \ A[1:N, (2N + 1):(3N + 1)])
        J +=
            -A[(2N + 1):(3N + 1), (N + 1):(2N)] * (
                Diagonal(A[(N + 1):(2N), (N + 1):(2N)]) \
                A[(N + 1):(2N), (2N + 1):(3N + 1)]
            )

        x[(2N + 1):(3N + 1)] +=
            -A[(2N + 1):(3N + 1), 1:N] * (Diagonal(A[1:N, 1:N]) \ b[1:N])
        x[(2N + 1):(3N + 1)] +=
            -A[(2N + 1):(3N + 1), (N + 1):(2N)] *
            (Diagonal(A[(N + 1):(2N), (N + 1):(2N)]) \ b[(N + 1):(2N)])

        x[(2N + 1):(3N + 1)] .= Tridiagonal(J) \ x[(2N + 1):(3N + 1)]
        x[1:N] .=
            Diagonal(A[1:N, 1:N]) \
            (b[1:N] - A[1:N, (2N + 1):(3N + 1)] * x[(2N + 1):(3N + 1)])
        x[(N + 1):(2N)] .=
            Diagonal(A[(N + 1):(2N), (N + 1):(2N)]) \ (
                b[(N + 1):(2N)] -
                A[(N + 1):(2N), (2N + 1):(3N + 1)] * x[(2N + 1):(3N + 1)]
            )

    end
end

Δt = 600.0
ndays = 1.0

# Solve the ODE operator
prob = ODEProblem(
    ODEFunction(
        tendency!,
        jac = jacobian!,
        jac_prototype = zeros(length(Y), length(Y)),
        tgrad = (dT, Y, p, t) -> fill!(dT, 0),
    ),
    Y,
    (0.0, 60 * 60 * 24 * ndays),
)
# 60 * 60 * 24 * ndays
sol = solve(
    prob,
    # ImplicitEuler(),
    Rosenbrock23(linsolve = linsolve!),
    dt = Δt,
    saveat = 60 * 60, # save every hour
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "hydrostatic_implicit"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

z_centers = parent(Fields.coordinate_field(cspace))
z_faces = parent(Fields.coordinate_field(fspace))

function hydrostatic_plot(u, Yc_init, w_init; title = "", size = (1024, 600))
    sub_plt1 = Plots.plot(
        parent(Yc_init.ρ),
        z_centers,
        marker = :circle,
        xlabel = "ρ",
        label = "T=0",
    )
    sub_plt1 = Plots.plot!(sub_plt1, parent(u.Yc.ρ), z_centers, label = "T")

    sub_plt2 = Plots.plot(
        parent(w_init),
        z_faces,
        marker = :circle,
        xlim = (-0.2, 0.2),
        xlabel = "ω",
        label = "T=0",
    )
    sub_plt2 = Plots.plot!(sub_plt2, parent(u.w), z_faces, label = "T")

    sub_plt3 = Plots.plot(
        parent(Yc_init.ρθ),
        z_centers,
        marker = :circle,
        xlabel = "ρθ",
        label = "T=0",
    )
    sub_plt3 = Plots.plot!(sub_plt3, parent(u.Yc.ρθ), z_centers, label = "T")

    return Plots.plot(
        sub_plt1,
        sub_plt2,
        sub_plt3,
        title = title,
        layout = (1, 3),
        size = size,
    )
end

anim = Plots.@animate for (i, u) in enumerate(sol.u)
    hydrostatic_plot(u, Y_init, w_init, title = "Hour $(i)")
end
Plots.mp4(anim, joinpath(path, "hydrostatic.mp4"), fps = 10)

Plots.png(
    hydrostatic_plot(sol[end], Y_init, w_init),
    joinpath(path, "hydrostatic_end.png"),
)

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("output/$(dirname)/hydrostatic_end.png", "Hydrostatic End")


text:
function pairwise_knn(X::AbstractMatrix{S},
                      n_neighbors,
                      metric) where {S <: AbstractFloat, V <: AbstractVector}
    num_points = size(X, 2)
    all_dists = Array{S}(undef, num_points, num_points)
    pairwise!(all_dists, metric, X)
    # all_dists is symmetric distance matrix
    knns = [partialsortperm(view(all_dists, :, i), 2:n_neighbors+1) for i in 1:num_points]
    dists = [all_dists[:, i][knns[i]] for i in eachindex(knns)]
    return hcat(knns...), hcat(dists...)
end


text:
using R
using Test

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


text:
@static if VERSION < v"1.1"

    using Base: promote_typejoin

    # backport #30076 just for Rows
    @inline function collect_to!(destination::Rows{Item}, iterator, offset, state) where Item
        # collect to destination array, checking the type of each result. if a result does not
        # match, widen the result type and re-dispatch.
        index = offset
        while true
            result = iterate(iterator, state)
            result === nothing && break
            item, state = result
            if item isa Item || typeof(item) === Item
                @inbounds destination[index] = item::Item
                index = index + 1
            else
                new = setindex_widen_up_to(destination, item, index)
                return collect_to!(new, iterator, index + 1, state)
            end
        end
        destination
    end

    @inline function setindex_widen_up_to(destination::AbstractArray{Item}, item, index) where Item
        new = similar(destination, promote_typejoin(Item, typeof(item)))
        copyto!(new, firstindex(new), destination, firstindex(destination), index - 1)
        @inbounds new[index] = item
        return new
    end

    @inline function grow_to!(destination::Rows, iterator, state)
        Item = eltype(destination)
        result = iterate(iterator, state)
        while result !== nothing
            item, state = result
            if item isa Item || typeof(item) === Item
                push!(destination, item::Item)
            else
                new = push_widen(destination, item)
                return grow_to!(new, iterator, state)
            end
            result = iterate(iterator, state)
        end
        destination
    end

    @inline function push_widen(destination, item)
        new = sizehint!(empty(
            destination,
            promote_typejoin(eltype(destination), typeof(item))
        ), length(destination))
        if new isa AbstractSet
            # TODO: merge back these two branches when copy! is re-enabled for sets/vectors
            union!(new, destination)
        else
            append!(new, destination)
        end
        push!(new, item)
        new
    end
end


text:
using Pollen
using GraphPlot
using FluxTraining
using FilePathsBase

project = Pollen.documentationproject(FluxTraining)
Pollen.fullbuild(project, Pollen.FileBuilder(Pollen.HTML(), p"dev/"))


text:
#=
Implement an autocomplete system. That is, given a query string s and a set of all possible query strings, return all strings in the set that have s as a prefix.

For example, given the query string de and the set of strings [dog, deer, deal], return [deer, deal].

Hint: Try preprocessing the dictionary into a more efficient data structure to speed up queries.
=#

using Test
include("Solutions/problem11_autocomplete.jl")

@test autocomplete("de", ["dog", "deer", "deal"]) == ["deer", "deal"]
@test autocomplete("or", ["origami", "parents", "organisation", "original", "officer", "or", "open", "meteor", "prow", "cabbage"]) == ["origami", "organisation", "original", "or"]

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using ReinforcementLearningExperiments

run(E`JuliaRL_BasicDQN_CartPole`)
run(E`JuliaRL_BC_CartPole`)
run(E`JuliaRL_DQN_CartPole`)
run(E`JuliaRL_PrioritizedDQN_CartPole`)
run(E`JuliaRL_Rainbow_CartPole`)
run(E`JuliaRL_QRDQN_CartPole`)
run(E`JuliaRL_REMDQN_CartPole`)
run(E`JuliaRL_IQN_CartPole`)
run(E`JuliaRL_VPG_CartPole`)
run(E`JuliaRL_BasicDQN_MountainCar`)
run(E`JuliaRL_DQN_MountainCar`)
run(E`JuliaRL_A2C_CartPole`)
run(E`JuliaRL_A2CGAE_CartPole`)
run(E`JuliaRL_PPO_CartPole`)
run(E`JuliaRL_MAC_CartPole`)
run(E`JuliaRL_DDPG_Pendulum`)
run(E`JuliaRL_SAC_Pendulum`)
run(E`JuliaRL_TD3_Pendulum`)
run(E`JuliaRL_PPO_Pendulum`)

run(E`JuliaRL_BasicDQN_SingleRoomUndirected`)


text:
export air_density, absolute_pressure, vapor_pressure_over_water

"""
    air_density(P, T, ϕ)

Calculate density of humid air [kg/m^3] given absolute pressure `P` [Pa], temperature `T` [K] and relative humidity `ϕ` [1].
"""
function air_density(P, T, ϕ)
    R_d = 287.058
    R_v = 461.495

    P_sat = vapor_pressure_over_water(T - 273.15)
    P_v = ϕ * P_sat
    P_d = P - P_v

    (P_d / R_d + P_v / R_v) / T
end

"""
    absolute_pressure(p0, h)

Calculate the absolute pressure [Pa] at altitude `h` [m] for a given sea level static pressure [Pa].
"""
function absolute_pressure(P0, h)
    R = 8.3144598 # universal gas constant
    M = 0.0289644 # molar mass of Earth's air
    g0 = 9.80665 # gravitational acceleration

    if h < 11000
        Tb = 288.15
        Pb = 101325.00
        hb = 0.0
        Lb = -0.0065 # standard temperature lapse rate (K/m)
    else
        error("not implemented for altitude $(h) m")
    end

    P0 * (Tb / (Tb + Lb * (h - hb)))^(g0 * M / R / Lb)
end

"""
    vapor_pressure_over_water(T)

Compute the saturation vapor pressure over liquid water in Pa given temperature `T` in °C.

# References
- Arden L Buck (1981) - New Equations for Computing Vapor Pressure and Enhancement Factor
- Arden L Buck (1996)
"""
function vapor_pressure_over_water(Tc)
    1000 * 0.61121 * exp((18.678 - Tc / 234.5) * (Tc / (257.14 + Tc)))
end


text:
using DelimitedFiles
using ArgParse

function ArgParse.parse_item(::Array{Int64, 1}, x::AbstractString)
    a = split(x, ('[', ']', ' ', ','))
    b = []
    for i in a
        (i == "") && (continue)
        push!(b, parse(Int64, i))
    end
    return b
end

include("dnn_to_milp.jl")
include("overlap_percent.jl")
include("types.jl")
include("data_parsing.jl")
include("sensor_placement.jl")
include("writer.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--node_file"
        help = "DNN node bias data file"
        arg_type = String
        default = "./data/Bias_1.csv"
    "--edge_file"
        help = "DNN edge weight data file"
        arg_type = String
        default = "./data/Node_1.csv"
    "--range_file"
        help = "DNN file for normalizing inputs"
        arg_type = String
        default = "./data/Ranges_1.csv"
    "--pmax"
        help = "required pmax (default = 0.1)"
        arg_type = Float64
        default = 0.1
    "--num_layers"
        help = "number of layers including input and output layers in DNN"
        arg_type = Int
        default = 4
    "--num_nodes_in_layer"
        help = "number of nodes in each layer"
        arg_type = Array{Int64, 1}
        default = [6, 20, 20, 1]
    "--casefile"
        help = "case file for tunnel orientation"
        arg_type = String
        default = "./data/TestCase_1.csv"
    "--wallcasefile"
        help = "data for walls"
        arg_type = String
        default = "./data/WallWPs_1.csv"
    "--output_straight"
        help = "output filename for straight placement"
        arg_type = String
        default = "./output/straight.csv"
    "--output_optimal"
        help = "output filename for optimal placement"
        arg_type = String
        default = "./output/optimal.csv"
end

parsed_args = parse_args(ARGS, s)

node_file = parsed_args["node_file"]
edge_file = parsed_args["edge_file"]
range_file = parsed_args["range_file"]
num_layers = parsed_args["num_layers"]
num_nodes_in_layer = parsed_args["num_nodes_in_layer"]
case_file = parsed_args["casefile"]
wall_case_file = parsed_args["wallcasefile"]

case_data = readdlm(case_file)

op = get_overlap_percentage(parsed_args)

drop_distance = case_data[3] * (1 - op)

@info "overlap percentage: $(round(op, digits=4))"
@info "drop distance: $(round(drop_distance, digits=2))"

tunnel_data = get_data(parsed_args)

straight_placement = get_straight_placements(tunnel_data, drop_distance)

over_placement = get_over_placements(tunnel_data, straight_placement, drop_distance)

op_placement = get_optimal_placements(tunnel_data,
    straight_placement, drop_distance)

write_placement(op_placement, op; file = parsed_args["output_optimal"])
write_placement(straight_placement, op; file = parsed_args["output_straight"])


text:
export QRDQNLearner, quantile_huber_loss

function quantile_huber_loss(ŷ, y; κ=1.0f0)
    N, B = size(y)
    Δ = reshape(y, N, 1, B) .- reshape(ŷ, 1, N, B)
    abs_error = abs.(Δ)
    quadratic = min.(abs_error, κ)
    linear = abs_error .- quadratic
    huber_loss = 0.5f0 .* quadratic .* quadratic .+ κ .* linear

    cum_prob = Zygote.ignore() do
        send_to_device(device(y), range(0.5f0 / N; length=N, step=1.0f0 / N))
    end
    loss = Zygote.dropgrad(abs.(cum_prob .- (Δ .< 0))) .* huber_loss
    mean(sum(loss;dims=1))
end

mutable struct QRDQNLearner{Tq <: AbstractApproximator,Tt <: AbstractApproximator,Tf,R} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    min_replay_history::Int
    update_freq::Int
    update_step::Int
    target_update_freq::Int
    sampler::NStepBatchSampler
    n_quantile::Int
    loss_func::Tf
    rng::R
    # for recording
    loss::Float32
end

"""
    QRDQNLearner(;kwargs...)

See paper: [Distributional Reinforcement Learning with Quantile Regression](https://arxiv.org/pdf/1710.10044.pdf)

# Keywords

- `approximator`::[`AbstractApproximator`](@ref): used to get quantile-values of a batch of states. The output should be of size `(n_quantile, n_action)`.
- `target_approximator`::[`AbstractApproximator`](@ref): similar to `approximator`, but used to estimate the quantile values of the next state batch.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `update_freq::Int=1`: the frequency of updating the `approximator`.
- `n_quantile::Int=1`: the number of quantiles.
- `target_update_freq::Int=100`: the frequency of syncing `target_approximator`.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `traces = SARTS`, set to `SLARTSL` if you are to apply to an environment of `FULL_ACTION_SET`.
- `loss_func`=[`quantile_huber_loss`](@ref).
"""
function QRDQNLearner(;
    approximator,
    target_approximator,
    stack_size::Union{Int,Nothing}=nothing,
    γ::Float32=0.99f0,
    batch_size::Int=32,
    update_horizon::Int=1,
    min_replay_history::Int=32,
    update_freq::Int=1,
    n_quantile::Int=1,
    target_update_freq::Int=100,
    traces=SARTS,
    update_step=0,
    loss_func=quantile_huber_loss,
    rng=Random.GLOBAL_RNG
)
    copyto!(approximator, target_approximator)
    sampler = NStepBatchSampler{traces}(;
        γ=γ,
        n=update_horizon,
        stack_size=stack_size,
        batch_size=batch_size,
    )

    N = n_quantile

    QRDQNLearner(
        approximator,
        target_approximator,
        min_replay_history,
        update_freq,
        update_step,
        target_update_freq,
        sampler,
        N,
        loss_func,
        rng,
        0.0f0,
    )
end

Flux.functor(x::QRDQNLearner) = (Q = x.approximator, Qₜ = x.target_approximator),
y -> begin
    x = @set x.approximator = y.Q
    x = @set x.target_approximator = y.Qₜ
    x
end

function (learner::QRDQNLearner)(env)
    s = send_to_device(device(learner.approximator), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    q = reshape(learner.approximator(s), learner.n_quantile, :)
    vec(mean(q, dims=1)) |> send_to_host
end

function RLBase.update!(learner::QRDQNLearner, batch::NamedTuple)
    Q = learner.approximator
    Qₜ = learner.target_approximator
    γ = learner.sampler.γ
    n = learner.sampler.n
    batch_size = learner.sampler.batch_size
    N = learner.n_quantile
    D = device(Q)
    loss_func = learner.loss_func

    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = CartesianIndex.(a, 1:batch_size)

    target_quantiles = reshape(Qₜ(s′), N, :, batch_size)
    qₜ = dropdims(mean(target_quantiles; dims=1); dims=1)
    aₜ = dropdims(argmax(qₜ, dims=1); dims=1)
    @views target_quantile_aₜ = target_quantiles[:, aₜ]
    y = reshape(r, 1, batch_size) .+ γ .* reshape(1 .- t, 1, batch_size) .* target_quantile_aₜ

    gs = gradient(params(Q)) do
        q = reshape(Q(s), N, :, batch_size)
        @views ŷ = q[:, a]

        loss = loss_func(ŷ, y)

        ignore() do
            learner.loss = loss
        end
        loss
    end

    update!(Q, gs)
end


text:
using Base.Test
using Lint

# lint the code
#@lintpragma( "Ignore dimensionless array field [regv]" )
@test isempty(lintpkg( "FactorizationMachines", returnMsgs=true))


text:
# Connect

import HTTP: request
import JSON: parse

"""
    struct IEXException <: Exception

    IEX Cloud uses [HTTP response codes](https://www.restapitutorial.com/httpstatuscodes.html)
    to indicate the success or failure of an API request.

    Found at [API Version](https://iexcloud.io/docs/api/#error-codes)
"""
struct IEXException <: Exception
    msg::AbstractString

    function IEXException(code::Int16)
        msg = join(get(IEXErrorCodes, code, "Unknown Error Code"), " or ")
        new(msg)
    end
end

const IEXErrorCodes = Dict{Int, Array{String,1}}(
    400 => [
        "Invalid values were supplied for the API request",
        "No symbol provided",
        "Batch request types parameter requires a valid value"
    ],
    401 => [
        "Hashed token authorization is restricted",
        "Hashed token authorization is required",
        "The requested data is marked restricted and the account does not have access.",
        "An API key is required to access the requested endpoint.",
        "The secret key is required to access to requested endpoint.",
        "The referer in the request header is not allowed due to API token domain restrictions."
    ],
    402 => [
        "You have exceeded your allotted message quota and pay-as-you-go is not enabled.",
        "The requested endpoint is not available to free accounts.",
        "The requested data is not available to your current tier."
    ],
    403 => [
        "Hashed token authorization is invalid.",
        "The provided API token has been disabled",
        "The provided API token is not valid.",
        "A test token was used for a production endpoint.",
        "A production token was used for a sandbox endpoint.",
        "Your pay-as-you-go circuit breaker has been engaged and further requests are not allowed.",
        "Your account is currently inactive."
    ],
    404 => [
        "Unknown symbol provided",
        "Resource not found"
    ],
    413 => ["Maximum number of types values provided in a batch request."],
    429 => ["Too many requests hit the API too quickly. An exponential backoff of your requests is recommended."],
    451 => ["The requested data requires additional permission to access."],
    500 => ["Something went wrong on an IEX Cloud server."]
)

"""
    global Iex connection

    mutable struct Connection <: Any

        Fields
        ≡≡≡≡≡≡≡≡

        url     :: String
        version :: String
        token   :: String

"""
mutable struct Connection <: Any
    url::String
    version::String
    token::String

    Connection(;url::String="", version::String="", token::String="") = begin
        @assert !isempty(url) "Iex URL is missing!"
        @assert !isempty(version) "Iex version is missing!"
        @assert !isempty(token) "Iex token is missing!"

        new(url, version, token)
    end
end

"""
    show method override for Connection

"""
Base.show(io::IO, z::Connection) = print(io,
    "{\n  URL: $(z.url)\n  Version: $(z.version)\n  Token: $(z.token)\n}"
)

"""
    set_connection_url!(conn::Connection, url::String)

Modifies url for connection object conn
"""
function set_connection_url!(conn::Connection, url::String)
    conn.url = url
end

"""
    set_connection_version!(conn::Connection, version::String)

Modifies version for connection object conn
"""
function set_connection_version!(conn::Connection, version::String)
    conn.version = version
end

"""
    set_connection_token!(conn::Connection, token::String)

Modifies url for connection object conn
"""
function set_connection_token!(conn::Connection, token::String)
    conn.token = token
end

"""
    status(conn::Connection)

Get status for IEX connection
"""
function status(conn::Connection)
    url = conn.url * conn.version * "/status"
    r = try
        request("GET", url;)
    catch e
        throw(IEXException[e.status])
    end
    return parse(String(r.body))
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

t = [0:0.2:2pi;]
plot(x=cos.(t), y=sin.(t), Geom.line(preserve_order=true))
plot(x=cos.(t), y=sin.(t), Geom.path)


text:
"""
struct containing grids used in a simulation

# Examples

```jldoctest
julia> using UltraDark

julia> len = 1;

julia> resol = 16;

julia> Grids(len, resol);

```
"""
struct Grids
    "Array of x positions"
    x::Array{Float64,3}
    "Array of y positions"
    y::Array{Float64,3}
    "Array of z positions"
    z::Array{Float64,3}
    "Array of x Fourier modes"
    kx::Array{Float64,3}
    "Array of y Fourier modes"
    ky::Array{Float64,3}
    "Array of z Fourier modes"
    kz::Array{Float64,3}
    "Fourier space postition array"
    k::Array{Float64,3}
    "Array of x Fourier modes for use with `rfft`"
    rkx::Array{Float64,3}
    "Array of y Fourier modes for use with `rfft`"
    rky::Array{Float64,3}
    "Array of z Fourier modes for use with `rfft`"
    rkz::Array{Float64,3}
    "Fourier space postition array for use with `rfft`"
    rk::Array{Float64,3}
    "ψ field"
    ψx::Array{Complex{Float64},3}
    "ψ field in Fourier space"
    ψk::Array{Complex{Float64},3}
    "density field ρ"
    ρx::Array{Float64,3}
    "density field ρ in Fourier space"
    ρk::Array{Complex{Float64},3}
    "gravitational potential field Φ"
    Φx::Array{Float64,3}
    "gravitational potential field Φ in fourier space"
    Φk::Array{Complex{Float64},3}

    "FFT plan for complex-to-complex transforms"
    fft_plan
    "FFT plan for real-to-complex transforms"
    rfft_plan

    "Indices at which k==0.0"
    k_vanish_indices::Vector{CartesianIndex{3}}
    "Indices at which rk==0.0"
    rk_vanish_indices#::Vector{CartesianIndex{3}}

    function Grids(x, y, z, kx, ky, kz, k, rkx, rky, rkz, rk, ψx, ψk, ρx, ρk, Φx, Φk)
        n_dims = 3
        resol_tuple = (size(x)[1], size(y)[2], size(z)[3])
        resol_tuple_realfft = (size(x)[1] ÷ 2 + 1, size(y)[2], size(z)[3])

        for var in [x, y, z, k, rk, ψx, ψk, ρx, ρk, Φx, Φk]
            @assert(ndims(var) == n_dims)
        end
        for var in [k, ψx, ψk, ρx, Φx]
            @assert(size(var) == resol_tuple)
        end
        for var in [rk, ρk, Φk]
            @assert(size(var) == resol_tuple_realfft)
        end
        @assert(size(x) == (resol_tuple[1], 1, 1))
        @assert(size(y) == (1, resol_tuple[2], 1))
        @assert(size(z) == (1, 1, resol_tuple[3]))

        FFTW.set_num_threads(Threads.nthreads())
        fft_plan = plan_fft(
            zeros(Complex{Float64}, resol_tuple),
            flags=FFTW.MEASURE
        )
        inv(fft_plan)

        rfft_plan = plan_rfft(
            zeros(Float64, resol_tuple),
            flags=FFTW.MEASURE
        )
        inv(rfft_plan)

        k_vanish_indices = findall(x -> x==0.0, k)
        rk_vanish_indices = findall(x -> x==0.0, rk)

        new(x, y, z, kx, ky, kz, k, rkx, rky, rkz, rk, ψx, ψk, ρx, ρk, Φx, Φk, fft_plan, rfft_plan, k_vanish_indices, rk_vanish_indices)
    end
end

"""
    Grids(length, resol::Int)

Constructor for `Grids`

Create an empty grid with length `length` and resolution `resol`

# Examples

```jldoctest
julia> using UltraDark

julia> Grids(1.0, 64);

```
"""
function Grids(length, resol::Int)::Grids

    Grids((length, length, length), (resol, resol, resol))

end

"""
    Grids(length_tuple, resol_tuple::Tuple{Int, Int, Int})

Constructor for `Grids`

Create an empty `length[1]`x`length[2]`x`length[3]` grid with resolution
`resol[1]`x`resol[2]`x`resol[3]`.

# Examples

```jldoctest
julia> using UltraDark

julia> Grids((1.0, 1.0, 0.5), (64, 64, 32));

```
"""
function Grids(length_tuple, resol_tuple::Tuple{Int, Int, Int})::Grids

    resol_tuple_realfft = (resol_tuple[1] ÷ 2 + 1, resol_tuple[2], resol_tuple[3])

    @assert length_tuple[1] / resol_tuple[1] ≈ length_tuple[2] / resol_tuple[2]
    @assert length_tuple[1] / resol_tuple[1] ≈ length_tuple[3] / resol_tuple[3]

    ψx = zeros(Complex{Float64}, resol_tuple)
    ψk = zeros(Complex{Float64}, resol_tuple)

    ρx = zeros(Float64, resol_tuple)
    ρk = zeros(Complex{Float64}, resol_tuple_realfft)

    Φx = zeros(Float64, resol_tuple)
    Φk = zeros(Complex{Float64}, resol_tuple_realfft)

    x = reshape(
                range(
                      -length_tuple[1] / 2 + length_tuple[1] / 2resol_tuple[1],
                      +length_tuple[1] / 2 - length_tuple[1] / 2resol_tuple[1],
                      length=resol_tuple[1]
                     ),
                resol_tuple[1], 1, 1
    )

    y = reshape(
                range(
                      -length_tuple[2] / 2 + length_tuple[2] / 2resol_tuple[2],
                      +length_tuple[2] / 2 - length_tuple[2] / 2resol_tuple[2],
                      length=resol_tuple[2]
                     ),
                1, resol_tuple[2], 1
    )

    z = reshape(
                range(
                      -length_tuple[3] / 2 + length_tuple[3] / 2resol_tuple[3],
                      +length_tuple[3] / 2 - length_tuple[3] / 2resol_tuple[3],
                      length=resol_tuple[3]
                     ),
                1, 1, resol_tuple[3]
    )

    kvec = k_vec(length_tuple, resol_tuple)

    kx = reshape(kvec[1], size(kvec[1])[1], 1, 1)
    ky = reshape(kvec[2], 1, size(kvec[2])[1], 1)
    kz = reshape(kvec[3], 1, 1, size(kvec[3])[1])
    k_norm = (kx.^2 .+ ky.^2 .+ kz.^2).^0.5

    rkvec = rk_vec(length_tuple, resol_tuple)

    rkx = reshape(rkvec[1], size(rkvec[1])[1], 1, 1)
    rky = reshape(rkvec[2], 1, size(rkvec[2])[1], 1)
    rkz = reshape(rkvec[3], 1, 1, size(rkvec[3])[1])
    rk_norm = (rkx.^2 .+ rky.^2 .+ rkz.^2).^0.5

    Grids(
        x, y, z,
        kx, ky, kz, k_norm,
        rkx, rky, rkz, rk_norm,
        ψx,
        ψk,
        ρx,
        ρk,
        Φx,
        Φk,
    )
end

"""
    k_vec(lengths, resols)

Calculate the Fourier frequencies of a box with side lengths `lengths` and resolutions `resols`

# Examples

```jldoctest
julia> using UltraDark: k_vec

julia> kvec = k_vec((2π, 2π, 2π), (4, 4, 4));

julia> kvec[1]
4-element AbstractFFTs.Frequencies{Float64}:
  0.0
  1.0
 -2.0
 -1.0

```
"""
function k_vec(lengths, resols)
    sample_rate = 2π .* resols ./ lengths

    kx = fftfreq(resols[1], sample_rate[1])

    ky = fftfreq(resols[2], sample_rate[2])
    kz = fftfreq(resols[3], sample_rate[3])

    kvec = (kx, ky, kz)
end

function rk_vec(lengths, resols)
    sample_rate = 2π .* resols ./ lengths

    kx = rfftfreq(resols[1], sample_rate[1])

    ky = fftfreq(resols[2], sample_rate[2])
    kz = fftfreq(resols[3], sample_rate[3])

    rkvec = (kx, ky, kz)
end

function k_norm(lengths, resols)
    kvec = k_vec(lengths, resols)

    kx = reshape(kvec[1], size(kvec[1])[1], 1, 1)
    ky = reshape(kvec[2], 1, size(kvec[2])[1], 1)
    kz = reshape(kvec[3], 1, 1, size(kvec[3])[1])

    (kx.^2 .+ ky.^2 .+ kz.^2).^0.5
end

function rk_norm(lengths, resols)
    kvec = rk_vec(lengths, resols)

    kx = reshape(kvec[1], size(kvec[1])[1], 1, 1)
    ky = reshape(kvec[2], 1, size(kvec[2])[1], 1)
    kz = reshape(kvec[3], 1, 1, size(kvec[3])[1])

    (kx.^2 .+ ky.^2 .+ kz.^2).^0.5
end


text:
# --------------
# Initialization
# --------------
Pkg.add("ArrayViews")
Pkg.add("Docile")
include("../src/LLLplus.jl")

using LLLplus
using Base.Test

# --------------
# tests with small matrices
# --------------
println("tests with small matrices...")

# Matrix from http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=442
H =[1   9   1   2;
    1   8   8   3;
    7   4   5   1;
    2   6   7   1];
B1=[2   3  -2  -4
    3  -1   2   1
    1   1   6  -4
    1   3  -1   3];
T1=[0   0   1  -1
    0   1  -1   0
    0   0   0   1
    1  -3   3  -2];
(B,T) = lll(H);
@test sum(abs(B-B1))==0
@test sum(abs(T-T1))==0
println("...done")

# --------------
# Run tests with random matrices
# --------------

println("\nIn all the following tests, the first time includes the "*
        "JIT compilation; \nfor the second execution the "*
        "compilation is already done and the time\n"*
        "should be faster.\n")

# Time LLL decomposition of a 1000x1000 real matrix with randn entries
N = 1000;
H = randn(N,N);
println("Testing LLL on $(N)x$(N) real matrix...")
@time (B,T,Q,R) = lll(H);
@time (B,T,Q,R) = lll(H);

# Time LLL, Seysen, VBLAST decompositions of a 10x10 complex matrix with
# randn entries
N = 10;
println("\nTesting LLL on $(N)x$(N) complex matrix...")
H = randn(N,N) + im*randn(N,N);
@time (B,T,Q,R) = lll(H);
@time (B,T,Q,R) = lll(H);
# println("Testing Seysen on same $(N)x$(N) complex matrix...")
# @time (B,T) = seysen(H);
# @time (B,T) = seysen(H);
println("Testing VBLAST on same $(N)x$(N) complex matrix...")
@time (W,P,B) = vblast(H);
@time (W,P,B) = vblast(H);

# Test sphere decoder
Ns = 100000;
N = 2;
println("\nTesting sphere decoder on $(Ns) samples of $(N)x$(N) BPSK system...")
NN = randn(N,Ns)/sqrt(100);
H = [.6 -2.6; -.5 .95];
C = [-1,1];
Z = rand(1:2,N,Ns);
X = C[Z];
Y = H*X+NN;
@time Zt = hard_sphere(Y,H,2);
@time Zt = hard_sphere(Y,H,2);
Zh = Zt+1;
errRate = sum(abs(Z-Zh))/Ns;
println("Error Rate is $(errRate). It should be zero or very small.\n")

# --------------
# test norm for matrix from http://www.latticechallenge.org/
# --------------

println("Testing now with 200x200 matrix from latticechallenge.org.")
println("The col norm of the input should be 233, col norm of the "*
        "reduced bases should be 30.")
mat = readdlm("challenge-200.mod",Int64) #run from parent directory

nrms = zeros(200,1)
for ix=1:200
    nrms[ix] = norm(mat[:,ix])
end
println("max col-norm of input is $(maximum(nrms))")

@time (B,T) = lll(mat);
nrms = zeros(200,1)
for ix=1:200
    nrms[ix] = norm(B[:,ix])
end
maxlll = maximum(nrms)
println("max col-norm of lll-reduced basis is $(maxlll)")
@test maxlll<=30+1e-6

@time (B,T) = seysen(mat);
nrms = zeros(200,1)
for ix=1:200
    nrms[ix] = norm(B[:,ix])
end
maxSeysen = maximum(nrms)
println("max column norm of seysen-reduced basis is $(maxSeysen)")
@test maxSeysen<=30+1e-6


text:
#
# Run all unit tests in UnitulAstrodynamics.jl
#

include("TwoBody/test_twobody.jl")
include("ThreeBody/test_threebody.jl")
include("NBody/test_nbody.jl")
include("Propagators/test_propagators.jl")
include("OrbitPlots/test_plots.jl")


text:
#
# This module is a wrapper for the C, unstructured (thus serial) interface of PIPS-NLP
#

module PipsNlpSolverSerial

using StructJuMPSolverInterface
using Libdl

export createProblem, solveProblem

function __init__()
try
  sharedLib=ENV["PIPS_NLP_SHARED_LIB"]
  #explicitly check if the file exists (because dlopen sometimes does not throw an error for invalid filenames, resulting in a seg fault)
  if(!isfile(sharedLib))
    error(string("The specified shared library ([", sharedLib, "]) does not exist"))
  end
  global libparpipsnlp=Libdl.dlopen(get(ENV,"PIPS_NLP_SHARED_LIB",""))
catch
  @warn("Could not load PIPS-NLP shared library. Make sure the ENV variable 'PIPS_NLP_SHARED_LIB' points to its location, usually in the PIPS repo at PIPS/build_pips/PIPS-NLP/libpipsnlp.so")
  rethrow()
end
end

function convert_to_c_idx(indicies)
    for i in 1:length(indicies)
        indicies[i] = indicies[i] - 1
    end
end

mutable struct PipsNlpProblem
    ref::Ptr{Nothing}
    n::Int
    m::Int
    x::Vector{Float64}
    g::Vector{Float64}
    obj_val::Float64
    status::Int

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h  # Can be nothing

    #jac , hess
    nzJac::Int
    nzHess::Int

    # For MathProgBase
    sense::Symbol

    function PipsNlpProblem(ref::Ptr{Nothing}, n, m, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, nzJac, nzHess)
        prob = new(ref, n, m, zeros(Float64, n), zeros(Float64, m), 0.0, 0,
        eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, nzJac, nzHess,
        :Min)

        finalizer(freeProblem, prob)
        # Return the object we just made
        prob
    end
end

###########################################################################
# Callback wrappers
###########################################################################
# Objective (eval_f)
function eval_f_wrapper(x_ptr::Ptr{Float64}, obj_ptr::Ptr{Float64}, user_data::Ptr{Nothing})
    # println(" julia - eval_f_wrapper " );
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::PipsNlpProblem
    # Calculate the new objective
    new_obj = convert(Float64, prob.eval_f(unsafe_wrap(Array,x_ptr, prob.n)))::Float64
    # Fill out the pointer
    unsafe_store!(obj_ptr, new_obj)
    # Done
    return Int32(1)
end

# Constraints (eval_g)
function eval_g_wrapper(x_ptr::Ptr{Float64}, g_ptr::Ptr{Float64}, user_data::Ptr{Nothing})
    # println(" julia - eval_g_wrapper " );
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::PipsNlpProblem
    # Calculate the new constraint values
    new_g = unsafe_wrap(Array,g_ptr, prob.m)
    prob.eval_g(unsafe_wrap(Array,x_ptr, prob.n), new_g)
    # Done
    return Int32(1)
end

# Objective gradient (eval_grad_f)
function eval_grad_f_wrapper(x_ptr::Ptr{Float64}, grad_f_ptr::Ptr{Float64}, user_data::Ptr{Nothing})
    # println(" julia -  eval_grad_f_wrapper " );
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::PipsNlpProblem
    # Calculate the gradient
    new_grad_f = unsafe_wrap(Array,grad_f_ptr, Int(prob.n))
    prob.eval_grad_f(unsafe_wrap(Array,x_ptr, Int(prob.n)), new_grad_f)
    if prob.sense == :Max
        new_grad_f *= -1.0
    end
    # Done
    return Int32(1)
end

# Jacobian (eval_jac_g)
function eval_jac_g_wrapper(x_ptr::Ptr{Float64}, values_ptr::Ptr{Float64}, iRow::Ptr{Cint}, jCol::Ptr{Cint},  user_data::Ptr{Nothing})
    # println(" julia -  eval_jac_g_wrapper " );
    # Extract Julia the problem from the pointer
    #@show user_data
    prob = unsafe_pointer_to_objref(user_data)::PipsNlpProblem
    #@show prob
    # Determine mode
    mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
    x = unsafe_wrap(Array,x_ptr, prob.n)
    irows = unsafe_wrap(Array,iRow, Int(prob.nzJac))
    kcols = unsafe_wrap(Array,jCol, Int(prob.n+1))
    values = unsafe_wrap(Array,values_ptr, Int(prob.nzJac))
    prob.eval_jac_g(x, mode, irows, kcols, values)
    if mode == :Structure
	    convert_to_c_idx(irows)
	    convert_to_c_idx(kcols)
	end
    # Done
    return Int32(1)
end

# Hessian
function eval_h_wrapper(x_ptr::Ptr{Float64}, lambda_ptr::Ptr{Float64}, values_ptr::Ptr{Float64}, iRow::Ptr{Cint}, jCol::Ptr{Cint}, user_data::Ptr{Nothing})
    # println(" julia - eval_h_wrapper " );
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::PipsNlpProblem
    # Did the user specify a Hessian
    if prob.eval_h === nothing
        # No Hessian provided
        return Int32(0)
    else
        # Determine mode
        mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
        x = unsafe_wrap(Array,x_ptr, prob.n)
        lambda = unsafe_wrap(Array,lambda_ptr, prob.m)
        irows = unsafe_wrap(Array,iRow, Int(prob.nzHess))
        kcols = unsafe_wrap(Array,jCol, Int(prob.n+1))
        values = unsafe_wrap(Array,values_ptr, Int(prob.nzHess))
        obj_factor = 1.0
        if prob.sense == :Max
            obj_factor *= -1.0
        end
        prob.eval_h(x, mode, irows, kcols, obj_factor, lambda, values)
        if mode == :Structure
	        convert_to_c_idx(irows)
	        convert_to_c_idx(kcols)
        end
        # Done
        return Int32(1)
    end
end

###########################################################################
# C function wrappers
###########################################################################
function createProblem(n::Int,m::Int,
    x_L::Vector{Float64},x_U::Vector{Float64},
    g_L::Vector{Float64},g_U::Vector{Float64},
    nzJac::Int, nzHess::Int,
    eval_f, eval_g,eval_grad_f,eval_jac_g,eval_h)

    @assert n == length(x_L) == length(x_U)
    @assert m == length(g_L) == length(g_U)
    eval_f_cb = @cfunction(eval_f_wrapper,Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Nothing}) )
    eval_g_cb = @cfunction(eval_g_wrapper,Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Nothing}) )
    eval_grad_f_cb = @cfunction(eval_grad_f_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Nothing}) )
    eval_jac_g_cb = @cfunction(eval_jac_g_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{Nothing}))
    eval_h_cb = @cfunction(eval_h_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{Nothing}))

    ret = ccall(Libdl.dlsym(libparpipsnlp,:CreatePipsNlpProblem),Ptr{Nothing},
            (Cint, Cint,
            Ptr{Float64}, Ptr{Float64},
            Ptr{Float64}, Ptr{Float64},
            Cint, Cint,
            Ptr{Nothing},Ptr{Nothing},
            Ptr{Nothing},Ptr{Nothing}, Ptr{Nothing}),
            n,m,
            x_L,x_U,
            g_L,g_U,
            nzJac, nzHess,
            eval_f_cb, eval_g_cb,
            eval_grad_f_cb, eval_jac_g_cb, eval_h_cb
            )
    # println(" ccall CreatePipsNlpProblem done ")

    if ret == C_NULL
        error("PIPS-NLP: Failed to construct problem.")
    else
        return(PipsNlpProblem(ret, n, m, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, nzJac, nzHess))
    end
end

function solveProblem(prob::PipsNlpProblem)
    # @show "solveProblem"

    final_objval = [0.0]
    ret = ccall(Libdl.dlsym(libparpipsnlp,:PipsNlpSolve), Cint,
            (Ptr{Nothing}, Ptr{Float64}, Ptr{Float64}, Any),
            prob.ref, final_objval, prob.x, prob)
    prob.obj_val = final_objval[1]
    prob.status = Int(ret)

    return prob.status
end

function freeProblem(prob::PipsNlpProblem)
    # @show "freeProblem"
    ret = ccall(Libdl.dlsym(libparpipsnlp,:FreePipsNlpProblem),
            Nothing, (Ptr{Nothing},),
            prob.ref)
    # @show ret
    return ret
end

end