text:
module NSum

using Grid
#using Cubature
#using Base.Collections
#using Distributions
import Base: ndims, show, length, collect

export nsum

typealias Location Vector{Int}
typealias Dimension Vector{Int}
#mins :: NTuple{K,Int}

type EvalPoint
    location::Location
    vals::Vector{Float64}
end

abstract Region{K}

type RegionTree{K} <: Region{K}
    root::Location
    sub1::Region{K}
    sub2::Region{K}
    level::Int
    dims::Dimension
    cutdim::Int
end

type RegionLeaf{K} <: Region{K}
    root::Location
    evals::Vector{EvalPoint}
    level::Int
    dims::Dimension
end

function _split(root, dims, level, cut_plane=1)
    newdims = copy(dims)
    K = length(dims)
    newdims[cut_plane]/=2
    newroot = copy(root)
    newroot[cut_plane]+=newdims[cut_plane]
    #@show cut_plane, root, newroot, dims, newdims
    leaf1 = RegionLeaf{K}(root, [], level+1, newdims)
    leaf2 = RegionLeaf{K}(newroot, [], level+1, newdims)
    # This is already handled by tree_from_leaf. FIXME Delete
    #for (i,c) in enumerate(Counter(zeros(K).+2))
        ## Now generate all points for leaf1 and leaf2
        #newpt = (c .- 1) .* div(newdims,2)
        #p1 = newpt .+ leaf1.root
        #p2 = newpt .+ leaf2.root

        #push!(leaf1.evals, EvalPoint(p1,[]))
        #push!(leaf2.evals, EvalPoint(p2,[]))
    #end
    return (leaf1, leaf2)
end
split_empty_tree(r::Region, cut_plane=1) = _split(r.root, r.dims, r.level, cut_plane)

function RegionTree(dims::Dimension)
    k = length(dims)
    leaf1, leaf2 = _split(zeros(k), dims, 1)
    RegionTree{k}(zeros(k), leaf1, leaf2, 1, dims, 1)
end

ndims{K}(::Type{RegionLeaf{K}}) = K
ndims{K}(::Type{RegionTree{K}}) = K
ndims(r::Region) = ndims(typeof(r))

graball(curr,state) = push!(state, curr)
collect(r::Region) = (allnodes = Region[]; dfs(graball, allnodes, r); allnodes)

function show(io::IO, r::Region)
    if isa(r, RegionTree)
        allnodes = collect(r)
    else
        allnodes = [r]
    end

    for i=r.level:maxdepth(r)
        print(io, "Level $i: ")
        levelnodes = filter(allnodes) do x
            x.level == i
        end
        for n in sort(levelnodes, by=x->x.root, lt=lexless)
            print(io, n.root)
            if isa(n, RegionLeaf) && !isempty(n.evals[1].vals)
                print(io, n.evals[1].vals)
            end
        end
        println(io, "")
    end
end

function tree_from_leaf{K}(r::RegionLeaf{K}, cut_plane)
    leaf1, leaf2 = _split(r.root, r.dims, r.level, cut_plane)
    # Copy evals from r into leaf1 and leaf2
    #leaf1_buf = EvalPoints[]
    #leaf2_buf = EvalPoints[]
    #for (i,c) in enumerate(Counter(zeros(K).+2))
        ## first save all pre-computed points
        #if c[cut_plane] == 1 # r.evals[i] point goes to leaf1
            #push!(leaf1_buf, r.evals[i])
        #else
            #push!(leaf2_buf, r.evals[i])
        #end
    #end
    newdims = leaf1.dims
    for (i,c) in enumerate(Counter(zeros(K).+2))
        # Now generate all points for leaf1 and leaf2
        newpt = (c .- 1) .* div(newdims,2)
        p1 = newpt .+ leaf1.root
        p2 = newpt .+ leaf2.root

        p1prev = p2prev = false
        # Have we already computed p1?
        for pt in r.evals
            if p1 == pt.location
                p1prev = true
                push!(leaf1.evals, pt)
            elseif p2 == pt.location
                p2prev = true
                push!(leaf2.evals, pt)
            end
        end
        p1prev || push!(leaf1.evals, EvalPoint(p1,[]))
        p2prev || push!(leaf2.evals, EvalPoint(p2,[]))
    end
    RegionTree{K}(r.root, leaf1, leaf2, r.level, r.dims, cut_plane)
end

function pick_cut_plane(r::RegionLeaf)
    K = ndims(r)
    assert(any(r.dims .> 2))
    candidates = find(x->x>2, r.dims) # only consider directions with enough space
    if rand()<0.2 || length(r.evals) == 0 || any(x->isempty(x.vals), r.evals) # no choice but to pick randomly
        return -Inf, candidates[rand(1:length(candidates))]
    else
        discrep = -Inf
        direc = -1
        #@show r.root, length(r.evals), r.evals
        for eval in r.evals[2:end-1]
            poss_direc = findfirst(eval.location .- r.evals[1].location)
            #@show poss_direc, eval
            r.dims[poss_direc] == 2 && continue
            grad = r.evals[1].vals[end] - eval.vals[end]
            if grad > discrep
                discrep = grad
                direc = poss_direc
            end
        end
        return discrep, direc
    end
end

function cut!(r::RegionTree)
    if isa(r.sub1, RegionTree)
        cut!(r.sub1)
    elseif any(r.sub1.dims .> 2) #RegionLeaf, so we have to cut from above
        _, cut_plane = pick_cut_plane(r.sub1)
        r.sub1 = tree_from_leaf(r.sub1, cut_plane)
    end

    # argh.. I know about DRY, but how?
    if isa(r.sub2, RegionTree)
        cut!(r.sub2)
    elseif any(r.sub2.dims .> 2) #RegionLeaf, so we have to cut from above
        _, cut_plane = pick_cut_plane(r.sub2)
        r.sub2 = tree_from_leaf(r.sub2, cut_plane)
    end
    nothing
end

function dfs(f::Function, state, r::RegionLeaf, onlyleaves=false)
    f(r, state)
end

function dfs(f::Function, state, r::RegionTree, onlyleaves=false)
    !onlyleaves && f(r, state)
    tovisit = Region[r.sub1, r.sub2]
    while length(tovisit) > 0
        curr = pop!(tovisit)
        if !onlyleaves # run on all
            f(curr, state)
        elseif isa(curr,RegionLeaf) # or only on leaves
            f(curr, state)
        end
        isa(curr, RegionTree) && push!(tovisit, curr.sub1, curr.sub2)
    end
end

function recursive_dfs(f::Function, state, r::RegionTree, onlyleaves=false)
    assert(isa(state, Dict))
    if !haskey(state, :history)
        state[:history] = Region[r]
    else
        push!(state[:history],r)
    end
    !onlyleaves && f(r, state)
    recursive_dfs(f, state, r.sub1, onlyleaves)
    recursive_dfs(f, state, r.sub2, onlyleaves)
    pop!(state[:history]) # clean up after ourselves
end

function recursive_dfs(f::Function, state, r::RegionLeaf, onlyleaves=false)
    push!(state[:history],r)
    f(r, state)
    pop!(state[:history]) # clean up after ourselves
end

function maxdepth_(r::Region, state)
    if r.level > state[1]
        state[1] = r.level
    end
end

function minlen_(r::Region, state)
    if r.len < state[1]
        state[1] = r.len
    end
end

length_(r::Region, state) = state[1] += 1
length(r::Region) = (state = {0}; dfs(length_, state, r); state[1])

maxdepth(r::Region) = (state = {0}; dfs(maxdepth_, state, r); state[1])
minlen(r::Region) = (state = {typemax(Int)}; dfs(minlen_, state, r); state[1])

function treesum(r::Region, fdim)
    s = zeros(fdim)
    dfs(treesum_, s, r, true)
    return s
end

treesum_(r::RegionTree, state) = nothing
function treesum_(r::RegionLeaf, state)
    cellvolume = prod(r.dims) / 2^ndims(r)
    for eval in r.evals
        state[:] += eval.vals * cellvolume
    end
end

addifdirty(r::RegionTree, state) = nothing
function addifdirty(r::RegionLeaf, state)
    assert(length(r.evals) > 0)
    for eval in r.evals
        if length(eval.vals)==0 # Don't want to add already computed ones
            push!(state, eval)
        end
    end
end

function max_(r::RegionLeaf, state)
    if all(r.dims .== 2) # Ignore fully explored regions
        return
    end
    k = ndims(r)
    buf = zeros(state[:fdim])
    treesum_(r, buf)
    if !haskey(state,:maxhistory) || buf[end] > state[:maxvalue]
        state[:maxhistory] = copy(state[:history])
        state[:maxvalue] = buf[end]
    end
end

function find_max(r::Region, fdim)
    s = {:fdim=>fdim}
    recursive_dfs(max_, s, r, true)
    !haskey(s,:maxhistory) && return r # We've evaluated everything
    return s[:maxhistory]
end

function find_min_containing(r::Region, pt)
    @assert all(0 .<= pt .- r.root .< r.dims)
    while true
        r = pt[r.cutdim] .< r.dims[r.cutdim] ? r.sub1 : r.sub2
        isa(r, RegionLeaf) && return r
    end
end

function find_uneven_branch(r::RegionTree)
    # we are still high, so we need to recurse downwards
    # return maximum discrepancy branch:
    path1 = find_uneven_branch(r.sub1)
    path2 = find_uneven_branch(r.sub2)
    maxpath = path1[1] > path2[1] ? path1 : path2
    push!(maxpath[2], r)
    return maxpath
end

function find_uneven_branch(r::RegionLeaf)
    if all(r.dims .== 2) # but we can't go lower
        return(-Inf, Region[r])
    else
        # Get a measure of discrepancy:
        discrep,direc = pick_cut_plane(r)
        return (discrep, Region[r])
    end
end

function nsum(fdim::Integer, f::Function, maxs; abstol=0.01, maxevals=30)
    #pseudocode:
    #create octree using closest larger point to maxs
    #subdivide a few times, adding points to evaluate to dirty list
    #evaluate in batch
    #update region values
    lens = map(nextpow2,iround(maxs))
    r = RegionTree(lens)
    for i=1:3
        cut!(r)
    end

    dirtylist = {}
    dfs(addifdirty, dirtylist, r, true)
    f(dirtylist)
    tots = zeros(fdim)
    dfs(treesum_, tots, r, true)
    #println("current estimate: $tots")

    count = 0
    findtype = 1
    totnum = fdim==1 ? 1 : fdim - 1
    while any(abs(tots[1:totnum].-1) .> abstol) && count < maxevals
        count += 1

        if findtype == 1
            maxr = find_max(r, fdim)
            r === maxr[end] && break # Exhausted everything
            maxr = maxr[end-1]
            #println("found max: $(maxr.root) w/ dims: $(maxr.dims)")

            #get regions in other 2^D-1 directions from max
            #others = Array(Any, 2^D-1)
            #resize!(others, 0)
            #for c in Counter(zeros(D).+2)
                #newloc = (c .- 1) .* 2 .-1 .+ maxr.mins
                #any(newloc .< 0) && continue
                #push!(others, newloc)
            #end
            ##println(others)

            ##subdivide these other regions
            #rs = map(others) do x
                #find_min_containing(r, x)
            #end
            cut!(maxr)
            cut!(maxr)
            findtype $= 1
        else
            dis, maxrs = find_uneven_branch(r)
            dis === -Inf && break # Exhausted everything
            maxr = maxrs[2] # 1 is the leaf
            #println("found uneven point: $(maxr.root) w/ dims: $(maxr.dims)")
            cut!(maxr)
            cut!(maxr)
            findtype $= 1
        end

        dirtylist = {}
        dfs(addifdirty, dirtylist, r, true)

        if length(dirtylist)>0
            f(dirtylist)
        end
        fill!(tots, 0.0)
        dfs(treesum_, tots, r, true)
        #println("Tots estimate: $tots")
    end
    #println("Finished $count iters")
    tots, r
end

end


text:
# -*- coding: utf-8 -*-
using GeophysicalFlows, Plots, Printf, Random

 dev = CPU()

n = 256
stepper = "FilteredETDRK4"
dt = 0.03
tf =60
nsteps = Int(tf / dt)
nsubs = round(Int, nsteps/100)

L = 2π
ν = 1e-19
nν = 4

prob = SurfaceQG.Problem(dev; nx=n, Lx=L, dt=dt,
    stepper=stepper, ν = ν , nν = nν )

sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y

X, Y = gridpoints(grid)
b = @. exp(-((X+L/4)^2 + 4*(Y+L/4)^2))
b .+= @. exp(-((X+L/4)^2 + 4*(Y-L/4)^2))
b .-= @. exp(-((X-L/4)^2 + 4*(Y+L/4)^2))
b .-= @. exp(-((X-L/4)^2 + 4*(Y-L/4)^2))

SurfaceQG.set_b!(prob, b )

heatmap(x, y, vars.b', aspectratio = 1,
               c = :deep,
            clim = (-1, 1),
           xlims = (-grid.Lx/2, grid.Lx/2),
           ylims = (-grid.Ly/2, grid.Ly/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
ylabel = "y",
title = "buoyancy b ",
      framestyle = :box)

B  = Diagnostic(SurfaceQG.buoyancy_variance, prob; nsteps=nsteps)
KE = Diagnostic(SurfaceQG.kinetic_energy, prob; nsteps=nsteps)
Dᵇ = Diagnostic(SurfaceQG.buoyancy_dissipation, prob; nsteps=nsteps)
diags = [B, KE, Dᵇ] # A list of Diagnostics types passed to `stepforward!`. Diagnostics are updated every timestep.

# +
base_filename = string("four_vortices_", n)

datapath = "./"
plotpath = "./"

dataname = joinpath(datapath, base_filename)
plotname = joinpath(plotpath, base_filename)

if !isdir(plotpath); mkdir(plotpath); end
if !isdir(datapath); mkdir(datapath); end
# -

get_sol(prob) = sol
get_u(prob) = irfft(im * grid.l .* sqrt.(grid.invKrsq) .* sol, grid.nx)
out = Output(prob, dataname, (:sol, get_sol), (:u, get_u))

function plot_output(prob)
  b = prob.vars.b

  pb = heatmap(x, y, Array(b'),
       aspectratio = 1,
                 c = :deep,
              clim = (-1, 1),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "buoyancy bₛ",
        framestyle = :box)

  pKE = plot(1,
             label = "kinetic energy ∫½(uₛ²+vₛ²)dxdy/L²",
         linewidth = 2,
            legend = :bottomright,
             alpha = 0.7,
             xlims = (0, tf),
             ylims = (0, 1e-2),
            xlabel = "t")

  pb² = plot(1,
             label = "buoyancy variance ∫bₛ²dxdy/L²",
         linecolor = :red,
            legend = :bottomright,
         linewidth = 2,
             alpha = 0.7,
             xlims = (0, tf),
             ylims = (0, 2e-2),
            xlabel = "t")

  layout = @layout [a{0.5w} Plots.grid(2, 1)]
  p = plot(pb, pKE, pb², layout=layout, size = (900, 500))

  return p
end

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:round(Int, nsteps/nsubs)
  if j % (500 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log1 = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min",
          clock.step, clock.t, cfl, (time()-startwalltime)/60)

    log2 = @sprintf("buoyancy variance: %.2e, buoyancy variance dissipation: %.2e",
              B.data[B.i], Dᵇ.data[Dᵇ.i])

    println(log1)

    println(log2)
  end

  p[1][1][:z] = vars.b
  p[1][:title] = "buoyancy, t=" * @sprintf("%.2f", clock.t)
  push!(p[2][1], KE.t[KE.i], KE.data[KE.i])
  push!(p[3][1], B.t[B.i], B.data[B.i])

  stepforward!(prob, diags, nsubs)
  SurfaceQG.updatevars!(prob)
end

mp4(anim, "four_vortices.mp4", fps=14)

pu = heatmap(x, y, Array(vars.u'),
     aspectratio = 1,
               c = :balance,
            clim = (-maximum(abs.(vars.u)), maximum(abs.(vars.u))),
           xlims = (-L/2, L/2),
           ylims = (-L/2, L/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")",
      framestyle = :box)

pv = heatmap(x, y, Array(vars.v'),
     aspectratio = 1,
               c = :balance,
            clim = (-maximum(abs.(vars.v)), maximum(abs.(vars.v))),
           xlims = (-L/2, L/2),
           ylims = (-L/2, L/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")",
      framestyle = :box)

pb = heatmap(x, y, Array(vars.b'),
     aspectratio = 1,
               c = :deep,
            clim = (-1, 1),
           xlims = (-L/2, L/2),
           ylims = (-L/2, L/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "bₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")",
      framestyle = :box)

layout = @layout [a{0.5h}; b{0.5w} c{0.5w}]

plot_final = plot(pb, pu, pv, layout=layout, size = (800, 800))

ext: jl


text:
"""
Drawings and calculations useful for research in the theory of
Iterated Functions Systems (IFSs) in the euclidean and complex planes.
"""
module SDDIFS

using Reexport

@reexport using SDDCore, SDDGeometry, SDDGraphics

include("IFSs.jl")

export
    AbstractIFS,
    AbstractIFSComplex,
    IFSComplex,
    IFSComplexAffine,
    #WIFSComplex,
    #WIFSComplexAffine,
    #AbstractIFSR2,
    #IFSR2,
    #IFSR2Affine,
    #WIFSR2,
    #WIFSR2Affine,
    #AbstractIFSR3,
    #IFSR3,
    #IFSR3Affine,
    #WIFSR3,
    #WIFSR3Affine,
    contractions,
    size

include("attractors.jl")

export
    drawattractorC
    #drawattractorR2
    #drawattractorR3

#include("densistymaps.jl")
#include("trappedpoints.jl")
#include("fractaldimensions.jl")

end # module


text:
import Base.Sort: Forward, Ordering, Algorithm, lt

using StatsBase

function getmedians(x, n)
    q,r = divrem(length(x), n+1)

    if q == 0
        return x
    end
    buckets = [q for _ in 1:n+1]
    for i in 1:r
        buckets[i] += 1
    end
    pop!(buckets)
    x[cumsum(buckets)]
end

function sortandsample_array(ord, xs, nsamples, presorted=false)
    r = sample(1:length(xs), min(length(xs), nsamples),
               replace=false, ordered=true)
    if !presorted
        sorted = sort(xs, order=ord)
        chunk = tochunk(sorted)
        samples = sorted[r]
    else
        chunk = nothing # avoid communicating metadata if already sorted
        samples = chunk[r]
    end
    (chunk, samples)
end

function batchedsplitmerge(chunks, splitters, batchsize, start_proc=1; merge=merge_sorted, by=identity, sub=getindex, order=default_ord)
    if batchsize >= length(chunks)
        if isempty(splitters)
            return [collect_merge(merge, chunks)]
        else
            return splitmerge(chunks, splitters, merge, by, sub, order)
        end
    end

    # group chunks into batches:
    q, r = divrem(length(chunks), batchsize)
    b = [batchsize for _ in 1:q]
    r != 0 && push!(b, r)
    batch_ranges = map(UnitRange, cumsum([1, b[1:end-1];]), cumsum(b))
    batches = map(x->chunks[x], batch_ranges)

    # splitmerge each batch
    topsplits, lowersplits = splitter_levels(order, splitters, length(chunks), batchsize)

    sorted_batches = map(batches) do b
        isempty(topsplits) ?
            [collect_merge(merge, b)] :
            splitmerge(b, topsplits, merge, by, sub, order)
    end

    range_groups = transpose_vecvec(sorted_batches)

    chunks = []
    p = start_proc
    for i = 1:length(range_groups)
        s = lowersplits[i]
        group = range_groups[i]
        if !isempty(s)
            cs = batchedsplitmerge(group, s, batchsize; merge = merge, by=by, sub=sub, order=order)
            append!(chunks, cs)
        else
            push!(chunks, collect_merge(merge, group))
        end
    end
    return chunks
end

function collect_merge(merge, group)
    #delayed((xs...) -> treereduce(merge, Any[xs...]))(group...)
    t = treereduce(delayed(merge), group)
end

# Given sorted chunks, splits each chunk according to splitters
# then merges corresponding splits together to form length(splitters) + 1 sorted chunks
# these chunks will be in turn sorted
function splitmerge(chunks, splitters, merge, by, sub, ord)
    c1 = map(c->splitchunk(c, splitters, by, sub, ord), chunks)
    map(cs->collect_merge(merge, cs), transpose_vecvec(c1))
end

function splitchunk(c, splitters, by, sub, ord)
    function getbetween(xs, lo, hi, ord)
        i = searchsortedlast(xs, lo, order=ord)
        j = searchsortedlast(xs, hi, order=ord)
        (i+1):j
    end

    function getgt(xs, lo, ord)
        i = searchsortedlast(xs, lo, order=ord)
        (i+1):length(xs)
    end

    function getlt(xs, hi, ord)
        j = searchsortedlast(xs, hi, order=ord)
        1:j
    end

    between = map((hi, lo) -> delayed(c->sub(c, getbetween(by(c), hi, lo, ord)))(c),
                  splitters[1:end-1], splitters[2:end])
    hi = splitters[1]
    lo = splitters[end]
    [delayed(c->sub(c, getlt(by(c), hi, ord)))(c);
     between; delayed(c->sub(c, getgt(by(c), lo, ord)))(c)]
end

# transpose a vector of vectors
function transpose_vecvec(xs)
    map(1:length(xs[1])) do i
        map(x->x[i], xs)
    end
end

function _promote_array{T,S}(x::AbstractArray{T}, y::AbstractArray{S})
    Q = promote_type(T,S)
    samehost = Distributed.check_same_host(procs())
    ok = (isa(x, Array) || isa(x, SharedArray)) && (isa(y, Array) || isa(y, SharedArray))
    if samehost && ok && isbits(Q)
        return Array{Q}(length(x)+length(y)) #, pids=procs())
    else
        return similar(x, Q, length(x)+length(y))
    end
end

function merge_sorted(ord::Ordering, x::AbstractArray, y::AbstractArray)
    z = _promote_array(x, y)
    _merge_sorted(ord, z, x, y)
end

function _merge_sorted{T, S}(ord::Ordering, z::AbstractArray{T}, x::AbstractArray{T}, y::AbstractArray{S})
    n = length(x) + length(y)
    i = 1; j = 1; k = 1
    len_x = length(x)
    len_y = length(y)
    @inbounds while i <= len_x && j <= len_y
        if lt(ord, x[i], y[j])
            z[k] = x[i]
            i += 1
        else
            z[k] = y[j]
            j += 1
        end
        k += 1
    end
    remaining, m = i <= len_x ? (x, i) : (y, j)
    @inbounds while k <= n
        z[k] = remaining[m]
        k += 1
        m += 1
    end
    z
end

function splitter_levels(ord, splitters, nchunks, batchsize)
    # final number of chunks
    noutchunks = length(splitters) + 1
    # chunks per batch
    perbatch = ceil(Int, nchunks / batchsize)
    root = getmedians(splitters, perbatch-1)

    subsplits = []
    i = 1
    for c in root
        j = findlast(x->lt(ord, x, c), splitters)
        push!(subsplits, splitters[i:j])
        i = j+2
    end
    push!(subsplits, splitters[i:end])
    root, subsplits
end

arrayorvcat(x::AbstractArray,y::AbstractArray) = vcat(x,y)
arrayorvcat(x,y) = [x,y]

const default_ord = Base.Sort.ord(isless, identity, false, Forward)

function dsort_chunks(cs, nchunks=length(cs), nsamples=2000;
                      merge = merge_sorted,
                      by=identity,
                      sub=getindex,
                      order=default_ord,
                      batchsize=max(2, nworkers()),
                      splitters=nothing,
                      chunks_presorted=false,
                      sortandsample = (x,ns, presorted)->sortandsample_array(order, x,ns, presorted),
                      affinities=workers(),
                     )
    if splitters !== nothing
        nsamples = 0 # no samples needed
    end

    cs1 = map(c->delayed(sortandsample)(c, nsamples, chunks_presorted), cs)
    xs = collect(treereduce(delayed(vcat), cs1))
    if length(cs1) == 1
        xs = [xs]
    end

    if splitters === nothing
        samples = reduce((a,b)->merge_sorted(order, a, b), map(x->x[2], xs))
        splitters = getmedians(samples, nchunks-1)
    end

    cs = batchedsplitmerge(map((x,c) -> first(x) === nothing ? c : first(x), xs, cs), splitters, batchsize; merge=merge, by=by, sub=sub, order=order)
    for (w, c) in zip(Iterators.cycle(affinities), cs)
        propagate_affinity!(c, Dagger.OSProc(w) => 1)
    end
    cs
end

function propagate_affinity!(c, aff)
    if !isa(c, Thunk)
        return
    end
    if !isnull(c.affinity)
        push!(get(c.affinity), aff)
    else
        c.affinity = [aff]
    end

    for t in c.inputs
        propagate_affinity!(t, aff)
    end
end

function Base.sort(v::ArrayOp;
               lt=Base.isless,
               by=identity,
               nchunks=nothing,
               rev::Bool=false,
               batchsize=max(2, nworkers()),
               nsamples=2000,
               order::Ordering=default_ord)
    v1 = compute(v)
    ord = Base.Sort.ord(lt,by,rev,order)
    nchunks = nchunks === nothing ? length(v1.chunks) : nchunks
    cs = dsort_chunks(v1.chunks, nchunks, nsamples,
                      order=ord, merge=(x,y)->merge_sorted(ord, x,y))
    t=delayed((xs...)->[xs...]; meta=true)(cs...)
    chunks = compute(t)
    dmn = ArrayDomain((1:sum(length(domain(c)) for c in chunks),))
    DArray(eltype(v1), dmn, map(domain, chunks), chunks)
end


text:
"""
	AbstractKnotPoint{Nx,Nu,V,T}

Stores the states, controls, time, and time step at a single point along a trajectory of a
forced dynamical system with `n` states and `m` controls.

# Interface
All instances of `AbstractKnotPoint` should support the following methods:

	state_dim(z)::Integer         # state vector dimension
	control_dim(z)::Integer       # control vector dimension
	getparams(z)                  # arbitrary parameters, usually the time and time step
	getdata(z)::V                 # the underlying vector
"""
abstract type AbstractKnotPoint{Nx,Nu,V,T} <: AbstractVector{T} end

dims(z::AbstractKnotPoint) = (state_dim(z), control_dim(z))
getstate(z::AbstractKnotPoint, v) = view(v, 1:state_dim(z))
getcontrol(z::AbstractKnotPoint, v) = begin n,m = dims(z); view(v, n+1:n+m) end
getstate(z::AbstractKnotPoint{Nx,Nu,<:SVector}, v) where {Nx,Nu} = v[SVector{Nx}(1:Nx)]
getcontrol(z::AbstractKnotPoint{Nx,Nu,<:SVector}, v) where {Nx,Nu} = v[SVector{Nu}(Nx+1:Nx+Nu)]

state(z::AbstractKnotPoint) = getstate(z, getdata(z))
control(z::AbstractKnotPoint) = getcontrol(z, getdata(z))

setdata!(z::AbstractKnotPoint, v) = z.z .= v
setdata!(z::AbstractKnotPoint{<:Any,<:Any,<:SVector}, v) = z.z = v

setstate!(z::AbstractKnotPoint, x) = state(z) .= x
setcontrol!(z::AbstractKnotPoint, u) = control(z) .= u
setstate!(z::AbstractKnotPoint{<:Any,<:Any,<:SVector}, x) = setdata!(z, [x; control(z)])
setcontrol!(z::AbstractKnotPoint{<:Any,<:Any,<:SVector}, u) = setdata!(z, [state(z); u])

time(z::AbstractKnotPoint) = getparams(z).t
timestep(z::AbstractKnotPoint) = getparams(z).dt

is_terminal(z::AbstractKnotPoint) = getparams(z).dt ≈ 0.0

# Array interface
Base.size(z::AbstractKnotPoint) = (state_dim(z) + control_dim(z), )
Base.getindex(z::AbstractKnotPoint, i::Int) = Base.getindex(getdata(z), i)
Base.setindex!(z::AbstractKnotPoint, v, i::Integer) = Base.setindex!(getdata(z), v, i)
Base.IndexStyle(z::AbstractKnotPoint) = IndexLinear()

for (name,mutable) in [(:KnotPoint, true), (:StaticKnotPoint, false)]
    expr = quote
        struct $name{Nx,Nu,V,T} <: AbstractKnotPoint{Nx,Nu,V,T}
            z::V
            t::Float64
            dt::Float64
            n::Int
            m::Int
            function $name{Nx,Nu}(z::V,t,dt) where {Nx,Nu,V}
                @assert Nx > 0
                @assert Nu > 0
                new{Nx,Nu,V,eltype(V)}(z, t, dt, Nx, Nu)
            end
            function $name{Nx,Nu}(x::V, u::V, t,dt) where {Nx,Nu,V}
                @assert Nx > 0
                @assert Nu > 0
                new{Nx,Nu,V,eltype(V)}([x;u], t, dt, Nx, Nu)
            end
            function $name(n::Integer, m::Integer, z::V, t, dt) where V
                @assert n > 0
                @assert m > 0
                new{Any,Any,V,eltype(V)}(z, t, dt, n, m)
            end
            function $name{Nx,Nu}(n::Integer, m::Integer, z::V, t, dt) where {Nx,Nu,V}
                Nx != Any && @assert Nx == n
                Nu != Any && @assert Nu == m
                new{Nx,Nu,V,eltype(V)}(z, t, dt, n, m)
            end
        end
        function $name(x::StaticVector{Nx}, u::StaticVector{Nu}, t, dt) where {Nx, Nu}
            $name{Nx,Nu}([x;u], t, dt)
        end
        function $name(x::AbstractVector, u::AbstractVector, t, dt)
            $name(length(x), length(u), [x;u], t, dt)
        end
        function Base.copy(z::$name{Nx,Nu}) where {Nx,Nu}
            $name{Nx,Nu}(z.n, z.m, z.z, z.t, z.dt)
        end

        state_dim(z::$name{Any}) = z.n
        control_dim(z::$name{<:Any,Any}) = z.m
        state_dim(z::$name{Nx}) where Nx = Nx
        control_dim(z::$name{<:Any,Nu}) where Nu = Nu

        @inline getparams(z::$name) = (t=z.t, dt=z.dt)
        @inline getdata(z::$name) = z.z

        # Array interface
        Base.similar(z::$name{Nx,Nu,V,T}, ::Type{S}, dims::Dims) where {Nx,Nu,V,T,S} =
            KnotPoint{Nx,Nu}(z.n, z.m, similar(z.z, S, dims), z.t, z.dt)
        Base.similar(z::$name{Nx,Nu,<:SVector,T}, ::Type{S}, dims::Dims) where {Nx,Nu,T,S} =
            KnotPoint{Nx,Nu}(z.n, z.m, similar(MVector{Nx+Nu,S}))
    end
    # Set struct mutability
    if mutable
        struct_expr = expr.args[2]
        @assert struct_expr.head == :struct
        struct_expr.args[1] = mutable
    end
    eval(expr)
end

function (::Type{<:KnotPoint{Nx,Nu}})(z, t, dt) where {Nx,Nu}
    KnotPoint{Nx,Nu}(z, t, dt)
end

function StaticKnotPoint(z::AbstractKnotPoint{Nx,Nu}, v::AbstractVector) where {Nx,Nu}
    StaticKnotPoint{Nx,Nu}(state_dim(z), control_dim(z), v, time(z), timestep(z))
end
@inline StaticKnotPoint(z::AbstractKnotPoint) = StaticKnotPoint(z, getdata(z))
@inline setdata(z::StaticKnotPoint, v) = StaticKnotPoint(z, v)
setstate(z::StaticKnotPoint, x) = setdata(z, [x; control(z)])
setcontrol(z::StaticKnotPoint, u) = setdata(z, [state(z); u])

function Base.:*(c::Real, z::KP) where {KP<:AbstractKnotPoint}
    KP(getdata(z)*c, getparams(z)...)
end


text:
# ---
# title: 1652. Defuse the Bomb
# id: problem1652
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/defuse-the-bomb/description/>
# hidden: true
# ---
#
# You have a bomb to defuse, and your time is running out! Your informer will
# provide you with a **circular** array `code` of length of `n` and a key `k`.
#
# To decrypt the code, you must replace every number. All the numbers are
# replaced **simultaneously**.
#
#   * If `k > 0`, replace the `ith` number with the sum of the **next** `k` numbers.
#   * If `k < 0`, replace the `ith` number with the sum of the **previous** `k` numbers.
#   * If `k == 0`, replace the `ith` number with `0`.
#
# As `code` is circular, the next element of `code[n-1]` is `code[0]`, and the
# previous element of `code[0]` is `code[n-1]`.
#
# Given the **circular** array `code` and an integer key `k`, return _the
# decrypted code to defuse the bomb_!
#
#
#
# **Example 1:**
#
#
#
#     Input: code = [5,7,1,4], k = 3
#     Output: [12,10,16,13]
#     Explanation: Each number is replaced by the sum of the next 3 numbers. The decrypted code is [7+1+4, 1+4+5, 4+5+7, 5+7+1]. Notice that the numbers wrap around.
#
#
# **Example 2:**
#
#
#
#     Input: code = [1,2,3,4], k = 0
#     Output: [0,0,0,0]
#     Explanation: When k is zero, the numbers are replaced by 0.
#
#
# **Example 3:**
#
#
#
#     Input: code = [2,4,9,3], k = -2
#     Output: [12,5,6,13]
#     Explanation: The decrypted code is [3+9, 2+3, 4+2, 9+4]. Notice that the numbers wrap around again. If k is negative, the sum is of the **previous** numbers.
#
#
#
#
# **Constraints:**
#
#   * `n == code.length`
#   * `1 <= n <= 100`
#   * `1 <= code[i] <= 100`
#   * `-(n - 1) <= k <= n - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using DistributionsDiff
using Test
using Distributions
using ChainRulesCore
import ForwardDiff

@testset "gradlogpdf" begin
    @testset "Exponential" begin
        x = 2.0
        θ = 3.0
        E, E_pullback = rrule(Exponential, θ)
        lp, lp_pullback = rrule(logpdf, E, x)

        vm = 1
        (_, res1, res2) = lp_pullback(vm)

        grad = ForwardDiff.gradient(pair -> logpdf(Exponential(pair[2]), pair[1]), [x, θ])
        @test extern(res1) ≈ grad[1]
        @test extern(res2) ≈ grad[2]
    end

    @testset "Normal" begin
        x = 2.0
        μ = 3.0
        σ = 1.0
        N, N_pullback = rrule(Normal, μ, σ)
        lp, lp_pullback = rrule(logpdf, N, x)

        vm = 1
        (_, res1, res2, res3) = lp_pullback(vm)

        grad = ForwardDiff.gradient(pair -> logpdf(Normal(pair[2], pair[3]), pair[1]), [x, μ, σ])
        @test extern(res1) ≈ grad[1]
        @test extern(res2) ≈ grad[2]
        @test extern(res3) ≈ grad[3]
    end

    # @testset "Arcsine" begin
    #     x = 4.0
    #     a = 3.0
    #     b = 5.0
    #     N, N_pullback = rrule(Normal, μ, σ)
    #     lp, lp_pullback = rrule(logpdf, N, x)

    #     vm = 1
    #     (_, res1, res2, res3) = lp_pullback(vm)

    #     grad = ForwardDiff.gradient(pair -> logpdf(Normal(pair[2], pair[3]), pair[1]), [x, μ, σ])
    #     @test extern(res1) ≈ grad[1]
    #     @test extern(res2) ≈ grad[2]
    #     @test extern(res3) ≈ grad[3]
    # end

    @testset "Beta" begin
        x = 0.4
        α = 3.0
        β = 5.0
        B, B_pullback = rrule(Beta, α, β)
        lp, lp_pullback = rrule(logpdf, B, x)

        vm = 1
        (_, res1, res2, res3) = lp_pullback(vm)

        grad = ForwardDiff.gradient(pair -> logpdf(Beta(pair[2], pair[3]), pair[1]), [x, α, β])
        @test extern(res1) ≈ grad[1]
        @test extern(res2) ≈ grad[2]
        @test extern(res3) ≈ grad[3]
    end

end


text:
#InspectDR: Rendering rectangular axes & ticks
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#==Main types
===============================================================================#

#==Rendering base plot elements
===============================================================================#

#Render main plot annotation (titles, axis labels, ...)
#-------------------------------------------------------------------------------
function render_baseannotation(ctx::CairoContext, rplot::RPlot2D, lyt::PlotLayout, a::Annotation)
	usemarkup=false #want user to be able to make fancy labels - but can cause crashes.
	TIMESTAMP_OFFSET = 3 #WANTCONST
	bb =  rplot.bb #WANTCONST
	databb =  rplot.databb #WANTCONST

	#Title
#	xcenter = (bb.xmin+bb.xmax)/2 #Entire plot BB.
	xcenter = (databb.xmin+databb.xmax)/2 #Data-area BB only (looks less odd).
	pt = Point2D(xcenter, databb.ymin-lyt.voffset_title)
	render(ctx, a.title, pt, lyt.font_title, align=ALIGN_HCENTER|ALIGN_BOTTOM, markup=usemarkup)

	#X-axis label
	xcenter = (databb.xmin+databb.xmax)/2
	pt = Point2D(xcenter, databb.ymax+lyt.voffset_xaxislabel)
	render(ctx, a.xlabel, pt, lyt.font_axislabel, align=ALIGN_HCENTER|ALIGN_TOP, markup=usemarkup)

	#Strip annotation
	nstrips = min(length(a.ylabels), length(rplot.strips))
	for i in 1:nstrips
		rstrip = rplot.strips[i]
		sbb = rstrip.bb
		ycenter = (sbb.ymin+sbb.ymax)/2
		xcenter = (sbb.xmin+sbb.xmax)/2

		#Strip label
		if i <= length(a.ystriplabels)
			lbl = a.ystriplabels[i]
			pt = Point2D(xcenter, sbb.ymin-lyt.voffset_ystriplabel)
			render(ctx, lbl, pt, lyt.font_striplabel, align=ALIGN_HCENTER|ALIGN_BOTTOM, markup=usemarkup)
		end

		#Y-axis label
		if i <= length(a.ylabels)
			lbl = a.ylabels[i]
			pt = Point2D(sbb.xmin-lyt.hoffset_yaxislabel, ycenter)
			render(ctx, lbl, pt, lyt.font_axislabel, align=ALIGN_HCENTER|ALIGN_BOTTOM, angle=-π/2, markup=usemarkup)
		end
	end

	#Time stamp
	if lyt.enable_timestamp
		pt = Point2D(bb.xmax-TIMESTAMP_OFFSET, bb.ymax-TIMESTAMP_OFFSET)
		render(ctx, a.timestamp, pt, lyt.font_time, align=ALIGN_RIGHT|ALIGN_BOTTOM, markup=usemarkup)
	end
end

#Render frame around graph
#-------------------------------------------------------------------------------
function render_graphframe(ctx::CairoContext, rstrip::RStrip2D, aa::AreaAttributes)
Cairo.save(ctx)
	setlinestyle(ctx, LineStyle(aa.line))
	Cairo.rectangle(ctx, rstrip.bb)
	Cairo.stroke(ctx)
Cairo.restore(ctx)
	return
end

#==Render grid
===============================================================================#
render_vlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, xlines::AbstractGridLines) = nothing
function render_vlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, xlines::GridLines)
	if xlines.displaymajor
		setlinestyle(ctx, lyt.line_gridmajor)
		for xline in xlines.major
			x = apply(xf, Point2D(xline, 0)).x
			drawline(ctx, Point2D(x, graphbb.ymin), Point2D(x, graphbb.ymax))
		end
	end
	if xlines.displayminor
		setlinestyle(ctx, lyt.line_gridminor)
		for xline in xlines.minor
			x = apply(xf, Point2D(xline, 0)).x
			drawline(ctx, Point2D(x, graphbb.ymin), Point2D(x, graphbb.ymax))
		end
	end
end
render_hlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, ylines::AbstractGridLines) = nothing
function render_hlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, ylines::GridLines)
	if ylines.displaymajor
		setlinestyle(ctx, lyt.line_gridmajor)
		for yline in ylines.major
			y = apply(xf, Point2D(0, yline)).y
			drawline(ctx, Point2D(graphbb.xmin, y), Point2D(graphbb.xmax, y))
		end
	end
	if ylines.displayminor
		setlinestyle(ctx, lyt.line_gridminor)
		for yline in ylines.minor
			y = apply(xf, Point2D(0, yline)).y
			drawline(ctx, Point2D(graphbb.xmin, y), Point2D(graphbb.xmax, y))
		end
	end
end

function render_grid(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, grid::GridRect)
	Cairo.save(ctx) #-----
	render_vlines(ctx, rstrip.bb, rstrip.xf, lyt, grid.xlines)
	render_hlines(ctx, rstrip.bb, rstrip.xf, lyt, grid.ylines)
	Cairo.restore(ctx) #-----
end

#==Render ticks
===============================================================================#

#Default label:
function render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, ::AxisScale)
	tstr = ""
	if false #Old display system.
		#TODO: deprecate
		tstr = "$val"
		if length(tstr) > 7 #HACK!
			tstr = @sprintf("%0.1e", val)
		end
	else
		tstr = formatted(val, fmt.fmt, showexp=!fmt.splitexp)
	end
	render(ctx, tstr, pt, font, align=align)
end

render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, scale::LogScale{:e}) =
	render_power(ctx, "e", val, pt, font, align)
render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, scale::LogScale{2}) =
	render_power(ctx, "2", val, pt, font, align)
render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, scale::LogScale{10}) =
	render_power(ctx, "10", val, pt, font, align)

function render_axisscalelabel(ctx::CairoContext, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, ::AxisScale)
	tstr = formatted_exp(fmt.fmt)
	render(ctx, tstr, pt, font, align=align)
end

#Render ticks: Well-defined GridLines
#-------------------------------------------------------------------------------
function render_xticks(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, xlines::GridLines, ticklabels::Bool)
	graphbb = rstrip.bb
	xs = rstrip.xscale

	#TODO: Fix LineAttributes to have concrete type:
	line_frame = LineStyle(lyt.frame_data.line) #Convert/do some validation
	line_major = LineStyle(line_frame.style, lyt.line_gridmajor.width, line_frame.color)
	line_minor = LineStyle(line_frame.style, lyt.line_gridminor.width, line_frame.color)
	tframe = line_frame.width #WANTCONST
	fmt = TickLabelFormatting(lyt.format_xtick, xlines.rnginfo)
	yframe = graphbb.ymax
	ylabel = graphbb.ymax + lyt.voffset_xticklabel

Cairo.save(ctx) #-----
	setlinestyle(ctx, line_major)
	for xtick in xlines.major
		x = apply(rstrip.xf, Point2D(xtick, 0)).x
		if ticklabels
			render_ticklabel(ctx, xtick, Point2D(x, ylabel), lyt.font_ticklabel, ALIGN_TOP|ALIGN_HCENTER, fmt, xs)
		end
		drawline(ctx, Point2D(x, yframe), Point2D(x, yframe-lyt.length_tickmajor))
	end
	setlinestyle(ctx, line_minor)
	for xtick in xlines.minor
		x = apply(rstrip.xf, Point2D(xtick, 0)).x
		drawline(ctx, Point2D(x, yframe), Point2D(x, yframe-lyt.length_tickminor))
	end
	if fmt.splitexp && ticklabels
		xlabel = graphbb.xmax + tframe
		render_axisscalelabel(ctx, Point2D(xlabel, yframe), lyt.font_ticklabel, ALIGN_BOTTOM|ALIGN_LEFT, fmt, xs)
	end
Cairo.restore(ctx) #-----
end
function render_yticks(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, ylines::GridLines)
	#TODO: Fix LineAttributes to have concrete type:
	line_frame = LineStyle(lyt.frame_data.line) #Convert/do some validation
	line_major = LineStyle(line_frame.style, lyt.line_gridmajor.width, line_frame.color)
	line_minor = LineStyle(line_frame.style, lyt.line_gridminor.width, line_frame.color)
	tframe = line_frame.width #WANTCONST
	graphbb = rstrip.bb
	ys = rstrip.yscale
	fmt = TickLabelFormatting(lyt.format_ytick, ylines.rnginfo)
	xframe = graphbb.xmin
	xlabel = graphbb.xmin - lyt.hoffset_yticklabel
Cairo.save(ctx) #-----
	setlinestyle(ctx, line_major)
	for ytick in ylines.major
		y = apply(rstrip.xf, Point2D(0, ytick)).y
		render_ticklabel(ctx, ytick, Point2D(xlabel, y), lyt.font_ticklabel, ALIGN_RIGHT|ALIGN_VCENTER, fmt, ys)
		drawline(ctx, Point2D(xframe, y), Point2D(xframe+lyt.length_tickmajor, y))
	end
	setlinestyle(ctx, line_minor)
	for ytick in ylines.minor
		y = apply(rstrip.xf, Point2D(0, ytick)).y
		drawline(ctx, Point2D(xframe, y), Point2D(xframe+lyt.length_tickminor, y))
	end
	if fmt.splitexp
		ylabel = graphbb.ymin - tframe
		render_axisscalelabel(ctx, Point2D(xframe, ylabel), lyt.font_ticklabel, ALIGN_BOTTOM|ALIGN_LEFT, fmt, ys)
	end
Cairo.restore(ctx) #-----
end
#Render ticks for colorscale (z-values)
function render_zticks(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, csbb::BoundingBox,
		zext::PExtents1D, zlines::GridLines, zs::AxisScale)
	tframe = DReal(lyt.frame_colorscale.line.width) #WANTCONST TODO: Fix LineAttributes to have concrete type
	fmt = TickLabelFormatting(lyt.format_ztick, zlines.rnginfo)
	xf = LinearTransform(zext, csbb.ymax, csbb.ymin) #ymax corresponds to min z value
	xframe = csbb.xmax
	xlabel = xframe + lyt.hoffset_yticklabel
	for ztick in zlines.major
		y = apply(xf, ztick)
		render_ticklabel(ctx, ztick, Point2D(xlabel, y), lyt.font_ticklabel, ALIGN_LEFT|ALIGN_VCENTER, fmt, zs)
		drawline(ctx, Point2D(xframe, y), Point2D(xframe-lyt.length_tickmajor, y))
	end
	for ztick in zlines.minor
		y = apply(xf, ztick)
		drawline(ctx, Point2D(xframe, y), Point2D(xframe-lyt.length_tickminor, y))
	end
	if fmt.splitexp
		ylabel = csbb.ymin - tframe
		render_axisscalelabel(ctx, Point2D(xframe, ylabel), lyt.font_ticklabel, ALIGN_BOTTOM|ALIGN_RIGHT, fmt, zs)
	end
end

#Render ticks: UndefinedGridLines
#-------------------------------------------------------------------------------
function render_xticks(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout,
		xlines::UndefinedGridLines, ticklabels::Bool)
	graphbb = rstrip.bb
	fmt = TickLabelFormatting(NoRangeDisplayInfo())
	yframe = graphbb.ymax
	ylabel = graphbb.ymax + lyt.voffset_xticklabel
	for (x, xlabel) in [(graphbb.xmin, xlines.minline), (graphbb.xmax, xlines.maxline)]
		if ticklabels
			render_ticklabel(ctx, xlabel, Point2D(x, ylabel), lyt.font_ticklabel, ALIGN_TOP|ALIGN_HCENTER, fmt, LinScale())
		end
		drawline(ctx, Point2D(x, yframe), Point2D(x, yframe-lyt.length_tickmajor))
	end
end
function render_yticks(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, ylines::UndefinedGridLines)
	fmt = TickLabelFormatting(NoRangeDisplayInfo())
	graphbb = rstrip.bb
	xframe = graphbb.xmin
	xlabel = graphbb.xmin - lyt.hoffset_yticklabel
	for (y, ylabel) in [(graphbb.ymax, ylines.minline), (graphbb.ymin, ylines.maxline)]
		render_ticklabel(ctx, ylabel, Point2D(xlabel, y), lyt.font_ticklabel, ALIGN_RIGHT|ALIGN_VCENTER, fmt, LinScale())
		drawline(ctx, Point2D(xframe, y), Point2D(xframe+lyt.length_tickmajor, y))
	end
end

function render_ticks(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, xticklabels::Bool)
	#NOTE: Use coordinate grid (cgrid) to render ticks
	render_xticks(ctx, rstrip, lyt, rstrip.cgrid.xlines, xticklabels)
	render_yticks(ctx, rstrip, lyt, rstrip.cgrid.ylines)
end

#==High-level rendering
===============================================================================#

#Render axis labels, ticks, ...
#-------------------------------------------------------------------------------
function render_axes(ctx::CairoContext, rstrip::RStrip2D, lyt::PlotLayout, xticklabels::Bool)
	render_graphframe(ctx, rstrip, lyt.frame_data)
	render_ticks(ctx, rstrip, lyt, xticklabels)
end

#Last line


text:
"""
`solve(P::Problem, O::Options=DefaultOptions)::Tuple{Float64, Float64, Vector}`

Solve the problem `P` using options `O`.

See `Problem` and `Options` for more info about how to state a problem and change the
default options.
"""
function solve(P::Problem, O::Options=DefaultOptions)::Tuple{Float64, Float64, Vector}

    # Unpack the data structures
    @unpack aᵢⱼ = P
    @unpack n, iterₘₐₓ, ϵₚ, ϵⱼ, get_x₀ = O

    # Get the matrix
    A = Symmetric(Array{Float64}(undef, n, n))
    for i = 1:n, j = i:n
        A.data[i,j] = aᵢⱼ(i, j)
    end

    # Get a vector of initial values
    x₀ = get_x₀(A)

    # Find the maximum modulo eigenvalue

    xₙ = copy(x₀)
    λₘₐₓ = 0
    λₚᵣₑᵥ = 1

    while abs(1 - λₘₐₓ / λₚᵣₑᵥ) > ϵₚ
        λₚᵣₑᵥ = λₘₐₓ
        x̃ₙ₊₁ = A * xₙ
        λₘₐₓ = x̃ₙ₊₁[1]
        xₙ = x̃ₙ₊₁ / λₘₐₓ
    end

    # Find the minimum modulo eigenvalue

    xₙ = copy(x₀)
    λ⁻¹ₘᵢₙ = 0
    λ⁻¹ₚᵣₑᵥ = 1

    while abs(1 - λ⁻¹ₘᵢₙ / λ⁻¹ₚᵣₑᵥ) > ϵₚ
        λ⁻¹ₚᵣₑᵥ = λ⁻¹ₘᵢₙ
        x̃ₙ₊₁ = A \ xₙ
        λ⁻¹ₘᵢₙ = x̃ₙ₊₁[1]
        xₙ = x̃ₙ₊₁ / λ⁻¹ₘᵢₙ
    end

    # Use the Jacobi eigenvalue algorithm to find all eigenvalues

    λ = Vector{Float64}(undef, n)

    for _ in 1:iterₘₐₓ

        aₘₐₓ = k = l = 0

        for i in 1:n-1, j in i+1:n
            if abs(A[i,j]) ≥ aₘₐₓ
                aₘₐₓ = abs(A[i,j])
                k = i; l = j
            end
        end

        if aₘₐₓ < ϵⱼ
            λ = Diagonal(A).diag
            break
        end

        Δ = A[l,l] - A[k,k]

        t = if abs(A[k,l]) < abs(Δ) * 1e-36
            A[k,l] / Δ
        else
            ϕ = Δ / (2 * A[k,l])
            if ϕ < 0
                -1 / (abs(ϕ) + √(ϕ^2 + 1))
            else
                1 / (abs(ϕ) + √(ϕ^2 + 1))
            end
        end

        c = 1 / √(t^2 + 1)
        s = t * c
        τ = s / (1 + c)

        aₜₑₘₚ = A[k,l]
        A.data[k,l] = 0
        A.data[k,k] -= t * aₜₑₘₚ
        A.data[l,l] += t * aₜₑₘₚ

        for i in 1:k-1
            aₜₑₘₚ = A[i,k]
            A.data[i,k] = aₜₑₘₚ - s * (A[i,l] + τ * aₜₑₘₚ)
            A.data[i,l] += s * (aₜₑₘₚ - τ * A[i,l])
        end

        for i in k+1:l-1
            aₜₑₘₚ = A[k,i]
            A.data[k,i] = aₜₑₘₚ - s * (A[i,l] + τ * A[k,i])
            A.data[i,l] += s * (aₜₑₘₚ - τ * A[i,l])
        end

        for i in l+1:n
            aₜₑₘₚ = A[k,i]
            A.data[k,i] = aₜₑₘₚ - s * (A[l,i] + τ * aₜₑₘₚ)
            A.data[l,i] += s * (aₜₑₘₚ - τ * A[l,i])
        end

    end

    return λₘₐₓ, 1 / λ⁻¹ₘᵢₙ, λ

end


text:
"""
    FunctionizeBridge{T}

The `FunctionizeBridge` converts a `SingleVariable` objective into a
`ScalarAffineFunction{T}` objective.
"""
struct FunctionizeBridge{T} <: AbstractBridge end
function bridge_objective(::Type{FunctionizeBridge{T}}, model::MOI.ModelLike,
                           func::MOI.SingleVariable) where T
    F = MOI.ScalarAffineFunction{T}
    MOI.set(model, MOI.ObjectiveFunction{F}(), convert(F, func))
    return FunctionizeBridge{T}()
end

function supports_objective_function(
    ::Type{<:FunctionizeBridge}, ::Type{MOI.SingleVariable})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:FunctionizeBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:FunctionizeBridge})
    return Tuple{DataType, DataType}[]
end
function MOIB.set_objective_function_type(::Type{FunctionizeBridge{T}}) where T
    return MOI.ScalarAffineFunction{T}
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::FunctionizeBridge, ::MOI.NumberOfVariables)
    return 0
end
function MOI.get(bridge::FunctionizeBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

# No variables or constraints are created in this bridge so there is nothing to
# delete.
function MOI.delete(model::MOI.ModelLike, bridge::FunctionizeBridge) end

function MOI.set(::MOI.ModelLike, ::MOI.ObjectiveSense,
                 ::FunctionizeBridge, ::MOI.OptimizationSense)
    # `FunctionizeBridge` is sense agnostic, therefore, we don't need to change
    # anything.
end
function MOI.get(model::MOI.ModelLike,
                 attr::MOIB.ObjectiveFunctionValue{MOI.SingleVariable},
                 bridge::FunctionizeBridge{T}) where T
    F = MOI.ScalarAffineFunction{T}
    return MOI.get(model, MOIB.ObjectiveFunctionValue{F}(attr.result_index))
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ObjectiveFunction{MOI.SingleVariable},
                 bridge::FunctionizeBridge{T}) where T
    F = MOI.ScalarAffineFunction{T}
    func = MOI.get(model, MOI.ObjectiveFunction{F}())
    return convert(MOI.SingleVariable, func)
end


text:
using ECharts

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]

y = [11, 11, 15, 13, 12, 13, 10]

bm = bar(x, hcat(0.95 .* y, 1.25 .* y, y))

# https://github.com/randyzwitch/ECharts.jl

max_stars_repo_name: huangzhen/aquila

max_issues_repo_name: huangzhen/aquila

max_forks_repo_name: huangzhen/aquila


text:
import StaticArrays

struct Vec3 <: StaticArrays.FieldVector{3, Float32}
    x::Float32
    y::Float32
    z::Float32
end

@inline shfl_sync(mask, val::Vec3, src) = Vec3(
    CUDA.shfl_sync(mask, val.x, src),
    CUDA.shfl_sync(mask, val.y, src),
    CUDA.shfl_sync(mask, val.z, src)
)


text: println("Happy New Year 1396")

max_stars_repo_name: AliMD/happy-new-year

max_issues_repo_name: AliMD/happy-new-year

max_forks_repo_name: AliMD/happy-new-year


text:
export PatchedSystem, PatchedSystemCache

"""
    PatchedSystem(F::AbstractSystem, patch, v::PVector)

Augment the system `F` with the given patch `v`. This results in the system `[F(x,t); v ⋅ x - 1]`
"""
struct PatchedSystem{F<:AbstractSystem, PS<:AbstractAffinePatchState} <: AbstractSystem
    system::F
    patch::PS
end

function PatchedSystem(F::AbstractSystem, p::AbstractAffinePatch, x)
    PatchedSystem(F, state(p, x))
end

function Base.size(F::PatchedSystem)
    m, n = size(F.system)
    (m + nequations(F.patch), n)
end

struct PatchedSystemCache{FC, T} <: AbstractSystemCache
    cache::FC
    A::Matrix{T} # intermediate storage of the jacobian
    b::Vector{T} # intermediate storage for the evaluation
end

function cache(ph::PatchedSystem, x)
    F = ph.system
    c = cache(F, x)
    b = evaluate(F, x, c)
    A = similar(b, size(F))
    PatchedSystemCache(c, A, b)
end

function evaluate!(u, F::PatchedSystem, x, c::PatchedSystemCache)
    M, N = size(F.system)

    # F(x)
    evaluate!(c.b, F.system, x, c.cache)
    @inbounds for i=1:M
        u[i] = c.b[i]
    end
    evaluate!(u, F.patch, x)
    u
end

function jacobian!(U, F::PatchedSystem, x::PVector, c::PatchedSystemCache)
    M, N = size(F.system)
    # J_F(x)
    jacobian!(c.A, F.system, x, c.cache)
    @inbounds for j=1:N, i=1:M
        U[i, j] = c.A[i, j]
    end
    jacobian!(U, F.patch, x)
    U
end

function evaluate_and_jacobian!(u, U, F::PatchedSystem, x::PVector, c::PatchedSystemCache)
    M, N = size(F.system)
    evaluate_and_jacobian!(c.b, c.A, F.system, x, c.cache)

    @inbounds for j=1:N, i=1:M
        U[i, j] = c.A[i, j]
    end
    jacobian!(U, F.patch, x)

    @inbounds for i=1:M
        u[i] = c.b[i]
    end
    evaluate!(u, F.patch, x)

    nothing
end

function evaluate(F::PatchedSystem, x, c::PatchedSystemCache)
    M, N = size(F)
    evaluate!(similar(c.b, M), F, x, c)
end

function jacobian(F::PatchedSystem, x, c::PatchedSystemCache)
    M, N = size(F)
    jacobian!(similar(c.A, M, N), F, x, c)
end


text:
TB=TightBinding

@testset "NRLTB" begin

println("Testing NRLTBModel: ")

if TESTDEPTH == 0
   Orbitals = ()
elseif TESTDEPTH == 1
   Orbitals = (:sp,)
else
   Orbitals = (:sp, :spd)
end

for orbitals in Orbitals
   println("Test NRLTB with $(orbitals)")
   at = (1,2,2) * bulk(:Si, pbc=(true, true, true), cubic=true)
   @show length(at)
   β, fixed_eF = 10.0, true
   tbm = TB.NRLTB.NRLTBModel(:Si, TB.GrandPotential(β, 0.0),
                           #  TB.FermiDiracSmearing(β, fixed_eF=fixed_eF),
                           orbitals = orbitals,
                           #   bzquad=TB.MPGrid(at, (4,0,0)) )
                             bzquad = TB.GammaPoint() )
   print("Test setting Nel ... ")
   TB.set_δNel!(tbm, at)
   println("ok.")
   @test true
   @show TB.get_eF(tbm.potential)

   # perturb positions a bit
   X = positions(at)
   X[2] += 0.1 * rand(JVecF)
   set_positions!(at, X)

   print("check that hamiltonian evaluates ... ")
   H, M = hamiltonian(tbm, at)
   println("ok.")
   print("check that `densitymatrix` evaluates ... ")
   Γ = densitymatrix(tbm, at)
   println("ok.")
   print("check that `energy` evaluates ... ")
   E = energy(tbm, at)
   println("ok : E = ", E, ".")
   print("check that `forces` evaluates ... ")
   frc = forces(tbm, at) |> mat
   println("ok : |f|∞ = ", vecnorm(frc, Inf), ".")
   @test true

   println("-------------------------------------------")
   println("  Finite-difference test with NRLTBModel:  ")
   println("-------------------------------------------")
   @test fdtest(tbm, at, verbose=true)
   println("-------------------------------------------")
end

println("Test sparse hamiltonian assembly")
for N = 2:6
   at = (N,N,N) * bulk(:Si, pbc=false, cubic=true)
   tbm = TightBinding.NRLTB.NRLTBModel(:Si, FermiDiracSmearing(1.0))
   H, M = hamiltonian(tbm, at)
   Hf, Mf =  TightBinding.evaluate(tbm.H, at, zeros(3), T = full)
   @show typeof(H)
   @test issparse(H) == (N > 3)
   @show vecnorm(full(H) - Hf, Inf)
   @test vecnorm(full(H) - Hf, Inf) == 0.0
end

end


text:
struct ImageGenerator
  model
  img
  target_class::Int
end

"""
    ImageGenerator(model, target_class)

Function to initialize the ImageGenerator for the `model` and `target_class`. Returns a `struct` with fields `model`, `img` and `target`. Here, the image is randomly generated which is then to be further updated.
"""
function ImageGenerator(model, target_class::Int)
  target_class > 1000 && error("Invalid Target Class")
  img = param(rand(224, 224, 3, 1)) |> gpu
  ImageGenerator(model.layers[1:end-1] |> gpu, img, target_class)
end

"""
    ImageGenerator(save_dir; img_name = "generated", niters::Int = 150, lr = 6.0)

Generates the image for which the prediction probability is maximized for a `target class` after `niter` iterations using backpropagation.
Before using this function, the ImageGenerator must be initialized using `ImageGenerator(model, target_class)`.

# Arguments
 - `save_dir`: Directory to which the generated image is to be saved
 - `img_name="generated"`: Image name with which file is to be saved withoout file extension
 - `niters::Int=150`: number of iterations for which the generated image is updated with backpropagation
 - `lr=6.0`: Learning rate to be used while updating the generated image with respect to backpropagated gradients to maximize the specific class.

"""
function (IG::ImageGenerator)(save_dir; img_name = "generated", niters::Int = 150, lr = 6.0)
  isdir(save_dir) || error("Storage must be a directory")
  mask = Float32.(zeros(1000, 1)) |> gpu
  mask[IG.target_class, 1] = Float32(1.0)
  local mean_img = Float32.(reshape([0.485, 0.456, 0.406], 1, 1, 3, 1)) |> gpu
  IG.img .= IG.img .- mean_img
  for i in 1:niters
    set_training_false()
    outputs, back = Zygote.pullback(IG.model,IG.img)
    set_training_true()
    loss = outputs .* mask
    grads = back(loss)
    println("Loss after Iteration $i is $(sum(loss)) and Probability is $(softmax(outputs)[IG.target_class, 1])")
    IG.img .= IG.img .+ lr * grads/sqrt(mean(abs2.(grads)))
    save("$(save_dir)/$(img_name)_iteration_$i.png", generate_image(IG.img))
  end
end


text:
# codegen.jl - utils to generate code in different formats from ExGraph and EinGraph

include("codegens/vec.jl")
include("codegens/buf.jl")
include("codegens/cuda.jl")
include("codegens/cudavec.jl")
include("codegens/gpu.jl")

function autoselect_codegen(inputs)
    # assuming inputs is Base.Iterators.Pairs or Vector{Any} with pairs
    first_array_idx = findfirst(a -> isa(a, AbstractArray) && eltype(a) != Any,
                                [v for (k,v) in inputs])
    first_array_idx != nothing || return VectorCodeGen()
    eltyp = eltype(inputs[first_array_idx][2])
    # we don't want to include CuArrays as dependency, so working on strings
    if any(startswith(string(typeof(v)), "CuArray") for (k, v) in inputs)
        return CuCodeGen(eltyp)
    else
        return BufCodeGen(eltyp)
    end
end


text:
#### Subdomain statistics

compute_subdomain_statistics(m::AtmosModel, state, aux, t) =
    compute_subdomain_statistics(
        m,
        state,
        aux,
        t,
        m.turbconv.micro_phys.statistical_model,
    )

"""
    compute_subdomain_statistics(
        m::AtmosModel{FT},
        state::Vars,
        aux::Vars,
        t::Real,
        statistical_model::SubdomainMean,
    ) where {FT}

Returns a cloud fraction and cloudy and dry thermodynamic
states in the subdomain.
"""
function compute_subdomain_statistics(
    m::AtmosModel{FT},
    state::Vars,
    aux::Vars,
    t::Real,
    statistical_model::SubdomainMean,
) where {FT}
    ts_en = recover_thermo_state_en(m, state, aux)
    cloud_frac = has_condensate(ts_en) ? FT(1) : FT(0)
    dry = ts_en
    cloudy = ts_en
    return (dry = dry, cloudy = cloudy, cloud_frac = cloud_frac)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using Parallelism
using Distributed: ProcessExitedException
using Memento
using Memento.TestUtils
using Test

const LOGGER = getlogger()

@testset "Parallelism.jl" begin
   include("robust_pmap.jl")
   include("tmap.jl")
end


text:
using SymPy, Documenter, SymFEL

makedocs(
    sitename="SymFEL.jl",
    format = Documenter.HTML(mathengine=Documenter.Writers.HTMLWriter.MathJax2()),

    pages=[
        "Home" => "index.md",
        "One dimensional examples" => "examples1.md",
        "Two dimensional examples" => "examples2.md",
        "Three dimensional examples" => "examples3.md"
        ]
    )


text:
module JungleHelperGeckoTutorialTrigger

using ..Ahorn, Maple

@mapdef Trigger "JungleHelper/GeckoTutorialTrigger" GeckoTutorialTrigger(x::Integer, y::Integer, width::Integer=8, height::Integer=8, geckoId::String="geckoId", showTutorial::Bool=true)

const placements = Ahorn.PlacementDict(
    "Gecko Tutorial (Jungle Helper)" => Ahorn.EntityPlacement(
       GeckoTutorialTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
struct AlwaysAssertionError <: Exception
    msg::String
end


text:
using GraphQLParser
using GraphQLParser: Document, Operation, FragmentDefinition,
    Field, FragmentSpread, InlineFragment, SelectionSet, Variable,
    Argument, Directive, VariableDefinition, Loc
using Test

@testset "GraphQLParser.jl" begin
    include("parser_tests.jl")
    include("validation.jl")
end


text:
function _maybe_download(url, filepath)
    if !isfile(filepath)
        mkpath(dirname(filepath))
        HTTP.download(url, filepath)
    end
    filepath
end


text:
import Base: isempty, start, next, done, length, getindex, minimum, maximum, search, show

export Treap, TreapNode

type TreapNode{K}
	priority::Float64
	size::Int
	key::K
	left::TreapNode{K}
	right::TreapNode{K}

	TreapNode(key::K, priority::Float64, left::TreapNode{K}, right::TreapNode{K}) =
		new(priority, left.size + right.size + 1, key, left, right)

	TreapNode(key::K, priority::Float64) =
		new(priority, 1, key, TreapNode{K}(), TreapNode{K}())

	TreapNode() = new(Inf, 0)
end

isempty(t::TreapNode) = t.size == 0
show(io::IO, t::TreapNode) = show(io, "key:$(t.key) size:$(t.size) priority:$(t.priority)")
start(t::TreapNode) = 1
next(t::TreapNode, state) = (t[state], state + 1)
done(t::TreapNode, state) = state > length(t)
length(t::TreapNode) = t.size
key(t::TreapNode) = t.key
left(t::TreapNode) = t.left
right(t::TreapNode) = t.right

function getindex{K}(t::TreapNode{K}, index::Int)
	1 <= index <= t.size || throw(KeyError(index))
	while t.left.size != index - 1
		if index <= t.left.size
			t = t.left
		else
			index = index - t.left.size - 1
			t = t.right
		end
	end
	t.key
end

function minimum{K}(t::TreapNode{K})
	isempty(t) && error("An empty treap has no minimum.")
	while !isempty(t.left) t = t.left end
	t.key
end

function maximum{K}(t::TreapNode{K})
	isempty(t) && error("An empty treap has no maximum.")
	while !isempty(t.right) t = t.right end
	t.key
end

function add!{K}(t::TreapNode{K}, key::K)
	isempty(t) && return TreapNode{K}(key, rand())
	t.size += 1
	if key < t.key
		t.left = add!(t.left, key)
		t.left.priority < t.priority ? rotate_right!(t) : t
	else
		@assert t.key < key "A treap may not contain duplicate keys: $key, $(t.key)"
		t.right = add!(t.right, key)
		t.right.priority < t.priority ? rotate_left!(t) : t
	end
end

function merge!{K}(left::TreapNode{K}, right::TreapNode{K})
	isempty(left) && return right
	isempty(right) && return left

	if left.priority < right.priority
		result = left
		result.size += right.size
		result.right = merge!(left.right, right)
	else
		result = right
		result.size += left.size
		result.left = merge!(left, result.left)
	end
	result
end

function remove!{K}(t::TreapNode{K}, key::K)
	isempty(t) && throw(KeyError(key))

	t.size -= 1
	if key == t.key
		merge!(t.left, t.right)
	elseif key < t.key
		t.left = remove!(t.left, key)
		t.left.priority < t.priority ? rotate_right!(t) : t
	else
		t.right = remove!(t.right, key)
		t.right.priority < t.priority ? rotate_left!(t) : t
	end
end

function rotate_right!{K}(root::TreapNode{K})
	@assert !isempty(root)
	newroot = root.left

	# Rotate
	root.left = newroot.right
	newroot.right = root

	# Update sizes
	root.size -= newroot.left.size + 1
	newroot.size += root.right.size + 1
	newroot
end

function rotate_left!{K}(root::TreapNode{K})
	@assert !isempty(root)
	newroot = root.right

	# Rotate
	root.right = newroot.left
	newroot.left = root

	# Update sizes
	root.size -= newroot.right.size + 1
	newroot.size += root.left.size + 1
	newroot
end

type Treap{K}
	root::TreapNode{K}
	Treap() = new(TreapNode{K}())
end
add!{K}(t::Treap{K}, key::K) = t.root = add!(t.root, key)
remove!{K}(t::Treap{K}, key::K) = t.root = remove!(t.root, key)
length(t::Treap) = length(t.root)
getindex(t::Treap, n::Int) = getindex(t.root, n)
start(t::Treap) = start(t.root)
next(t::Treap, state) = next(t.root, state)
done(t::Treap, state) = done(t.root, state)
minimum(t::Treap) = minimum(t.root)
maximum(t::Treap) = maximum(t.root)
key(t::Treap) = key(t.root)
left(t::Treap) = left(t.root)
right(t::Treap) = right(t.root)
root(t::Treap) = t.root

function test_treap()
	n = 10000
	a = shuffle([i for i in 1:n])
	t = Treap{Int}()

	for i in 1:n
		add!(t, a[i])
	end
	sort!(a)

	@assert !isempty(t)
	@assert length(t) == length(a) "$(length(t)) != $(length(a))"
	@assert maximum(t) == maximum(a)
	@assert minimum(t) == minimum(a)

	# Check that all elements have been added to
	# the treap and are in sorted order.
	for (i, v) in enumerate(a)
		@assert t[i] == v
	end

	for v in a
		remove!(t, v)
	end

	@assert length(t) == 0
	@assert isempty(t)

	println("Treap: Test succeeded.")

end

function benchmark_treap(n)
	gc_disable()

	t = Treap{Int}()
	a = shuffle([i for i in 1:n])

	println("Timing $n insert operations.")
	@time for i in 1:n
		t = add!(t, a[i])
	end

	@assert length(t) == n

	println("Timing $n random access operations.")
	@time for i in 1:n
		t[rand(1:n)]
	end

	println("Timing $n remove operations.")
	@time for i in 1:n
		remove!(t, a[i])
	end

	@assert isempty(t)

	gc_enable()
end

# test_treap()
# benchmark_treap(100000)

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
@doc raw"""
    cost_acceleration_bezier(
        M::AbstractManifold,
        B::AbstractVector{P},
        degrees::AbstractVector{<:Integer},
        T::AbstractVector{<:AbstractFloat},
    ) where {P}

compute the value of the discrete Acceleration of the composite Bezier curve

```math
\sum_{i=1}^{N-1}\frac{d^2_2 [ B(t_{i-1}), B(t_{i}), B(t_{i+1})]}{\Delta_t^3}
```

where for this formula the `pts` along the curve are equispaced and denoted by
``t_i``, ``i=1,…,N``, and ``d_2`` refers to the second order absolute difference [`costTV2`](@ref)
(squared). Note that the Beziér-curve is given in reduces form as a point on a `PowerManifold`,
together with the `degrees` of the segments and assuming a differentiable curve, the segments
can internally be reconstructed.

This acceleration discretization was introduced in[^BergmannGousenbourger2018].

# See also

[`grad_acceleration_bezier`](@ref), [`cost_L2_acceleration_bezier`](@ref), [`grad_L2_acceleration_bezier`](@ref)

[^BergmannGousenbourger2018]:
    > Bergmann, R. and Gousenbourger, P.-Y.: A variational model for data fitting on
    > manifolds by minimizing the acceleration of a Bézier curve.
    > Frontiers in Applied Mathematics and Statistics (2018).
    > doi [10.3389/fams.2018.00059](http://dx.doi.org/10.3389/fams.2018.00059),
    > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
"""
function cost_acceleration_bezier(
    M::AbstractManifold,
    B::AbstractVector{P},
    degrees::AbstractVector{<:Integer},
    T::AbstractVector{<:AbstractFloat},
) where {P}
    Bt = get_bezier_segments(M, B, degrees, :differentiable)
    p = de_casteljau(M, Bt, T)
    n = length(T)
    f = p[[1, 3:n..., n]]
    b = p[[1, 1:(n - 2)..., n]]
    d = distance.(Ref(M), p, shortest_geodesic.(Ref(M), f, b, Ref(0.5))) .^ 2
    samplingFactor = 1 / (((max(T...) - min(T...)) / (n - 1))^3)
    return samplingFactor * sum(d)
end
@doc raw"""
    cost_L2_acceleration_bezier(M,B,pts,λ,d)

compute the value of the discrete Acceleration of the composite Bezier curve
together with a data term, i.e.

````math
\frac{λ}{2}\sum_{i=0}^{N} d_{\mathcal M}(d_i, c_B(i))^2+
\sum_{i=1}^{N-1}\frac{d^2_2 [ B(t_{i-1}), B(t_{i}), B(t_{i+1})]}{\Delta_t^3}
````

where for this formula the `pts` along the curve are equispaced and denoted by
``t_i`` and ``d_2`` refers to the second order absolute difference [`costTV2`](@ref)
(squared), the junction points are denoted by ``p_i``, and to each ``p_i`` corresponds
one data item in the manifold points given in `d`. For details on the acceleration
approximation, see [`cost_acceleration_bezier`](@ref).
Note that the Beziér-curve is given in reduces form as a point on a `PowerManifold`,
together with the `degrees` of the segments and assuming a differentiable curve, the
segments can internally be reconstructed.

# See also

[`grad_L2_acceleration_bezier`](@ref), [`cost_acceleration_bezier`](@ref), [`grad_acceleration_bezier`](@ref)
"""
function cost_L2_acceleration_bezier(
    M::AbstractManifold,
    B::AbstractVector{P},
    degrees::AbstractVector{<:Integer},
    T::AbstractVector{<:AbstractFloat},
    λ::AbstractFloat,
    d::AbstractVector{P},
) where {P}
    Bt = get_bezier_segments(M, B, degrees, :differentiable)
    p = get_bezier_junctions(M, Bt)
    return cost_acceleration_bezier(M, B, degrees, T) +
           λ / 2 * sum((distance.(Ref(M), p, d)) .^ 2)
end

@doc raw"""
    costIntrICTV12(M, f, u, v, α, β)

Compute the intrinsic infimal convolution model, where the addition is replaced
by a mid point approach and the two functions involved are [`costTV2`](@ref)
and [`costTV`](@ref). The model reads

```math
E(u,v) =
  \frac{1}{2}\sum_{i ∈ \mathcal G}
    d_{\mathcal M}\bigl(g(\frac{1}{2},v_i,w_i),f_i\bigr)
  +\alpha\bigl( β\mathrm{TV}(v) + (1-β)\mathrm{TV}_2(w) \bigr).
```
"""
function costIntrICTV12(M::AbstractManifold, f, u, v, α, β)
    IC = 1 / 2 * distance(M, shortest_geodesic(M, u, v, 0.5), f)^2
    TV12 = β * costTV(M, u) + (1 - β) * costTV2(M, v)
    return IC + α * TV12
end

@doc raw"""
    costL2TV(M, f, α, x)

compute the ``ℓ^2``-TV functional on the `PowerManifold manifold `M` for given
(fixed) data `f` (on `M`), a nonnegative weight `α`, and evaluated at `x` (on `M`),
i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \alpha \operatorname{TV}(x)
```

# See also
[`costTV`](@ref)
"""
costL2TV(M, f, α, x) = 1 / 2 * distance(M, f, x)^2 + α * costTV(M, x)

@doc raw"""
    costL2TVTV2(M, f, α, β, x)

compute the ``ℓ^2``-TV-TV2 functional on the `PowerManifold` manifold `M` for
given (fixed) data `f` (on `M`), nonnegative weight `α`, `β`, and evaluated
at `x` (on `M`), i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \alpha\operatorname{TV}(x)
  + β\operatorname{TV}_2(x)
```

# See also
[`costTV`](@ref), [`costTV2`](@ref)
"""
function costL2TVTV2(M::PowerManifold, f, α, β, x)
    return 1 / 2 * distance(M, f, x)^2 + α * costTV(M, x) + β * costTV2(M, x)
end

@doc raw"""
    costL2TV2(M, f, β, x)

compute the ``ℓ^2``-TV2 functional on the `PowerManifold` manifold `M`
for given data `f`, nonnegative parameter `β`, and evaluated at `x`, i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + β\operatorname{TV}_2(x)
```

# See also
[`costTV2`](@ref)
"""
function costL2TV2(M::PowerManifold, f, β, x)
    return 1 / 2 * distance(M, f, x)^2 + β * costTV2(M, x)
end

@doc raw"""
    costTV(M, x, p)

Compute the ``\operatorname{TV}^p`` functional for a tuple `pT` of pointss
on a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) `M`, i.e.

```math
E(x_1,x_2) = d_{\mathcal M}^p(x_1,x_2), \quad x_1,x_2 ∈ \mathcal M
```

# See also

[`grad_TV`](@ref), [`prox_TV`](@ref)
"""
function costTV(M::AbstractManifold, x::Tuple{T,T}, p::Int=1) where {T}
    return distance(M, x[1], x[2])^p
end
@doc raw"""
    costTV(M,x [,p=2,q=1])

Compute the ``\operatorname{TV}^p`` functional for data `x`on the `PowerManifold`
manifold `M`, i.e. ``\mathcal M = \mathcal N^n``, where ``n ∈ \mathbb N^k`` denotes
the dimensions of the data `x`.
Let ``\mathcal I_i`` denote the forward neighbors, i.e. with ``\mathcal G`` as all
indices from ``\mathbf{1} ∈ \mathbb N^k`` to ``n`` we have
``\mathcal I_i = \{i+e_j, j=1,…,k\}\cap \mathcal G``.
The formula reads

```math
E^q(x) = \sum_{i ∈ \mathcal G}
  \bigl( \sum_{j ∈  \mathcal I_i} d^p_{\mathcal M}(x_i,x_j) \bigr)^{q/p}.
```

# See also
[`grad_TV`](@ref), [`prox_TV`](@ref)
"""
function costTV(M::PowerManifold, x, p=1, q=1)
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R)
    cost = fill(0.0, Tuple(power_size))
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i -> (i == k) ? 1 : 0, d)) #k th unit vector
        for i in R # iterate over all pixel
            j = i + ek # compute neighbor
            if all(map(<=, j.I, maxInd.I)) # is this neighbor in range?
                cost[i] += costTV(M.manifold, (x[M, Tuple(i)...], x[M, Tuple(j)...]), p)
            end
        end
    end
    cost = (cost) .^ (1 / p)
    if q > 0
        return sum(cost .^ q)^(1 / q)
    else
        return cost
    end
end
@doc raw"""
    costTV2(M,(x1,x2,x3) [,p=1])

Compute the ``\operatorname{TV}_2^p`` functional for the 3-tuple of points
`(x1,x2,x3)`on the [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) `M`. Denote by

```math
  \mathcal C = \bigl\{ c ∈  \mathcal M \ |\ g(\tfrac{1}{2};x_1,x_3) \text{ for some geodesic }g\bigr\}
```

the set of mid points between ``x_1`` and ``x_3``. Then the function reads

```math
d_2^p(x_1,x_2,x_3) = \min_{c ∈ \mathcal C} d_{\mathcal M}(c,x_2).
```

# See also
[`grad_TV2`](@ref), [`prox_TV2`](@ref)
"""
function costTV2(M::MT, x::Tuple{T,T,T}, p=1) where {MT<:AbstractManifold,T}
    # note that here mid_point returns the closest to x2 from the e midpoints between x1 x3
    return 1 / p * distance(M, mid_point(M, x[1], x[3]), x[2])^p
end
@doc raw"""
    costTV2(M,x [,p=1])

compute the ``\operatorname{TV}_2^p`` functional for data `x` on the
`PowerManifold` manifoldmanifold `M`, i.e. ``\mathcal M = \mathcal N^n``,
where ``n ∈ \mathbb N^k`` denotes the dimensions of the data `x`.
Let ``\mathcal I_i^{\pm}`` denote the forward and backward neighbors, respectively,
i.e. with ``\mathcal G`` as all indices from ``\mathbf{1} ∈ \mathbb N^k`` to ``n`` we
have ``\mathcal I^\pm_i = \{i\pm e_j, j=1,…,k\}\cap \mathcal I``.
The formula then reads

```math
E(x) = \sum_{i ∈ \mathcal I,\ j_1 ∈  \mathcal I^+_i,\ j_2 ∈  \mathcal I^-_i}
d^p_{\mathcal M}(c_i(x_{j_1},x_{j_2}), x_i),
```

where ``c_i(⋅,⋅)`` denotes the mid point between its two arguments that is
nearest to ``x_i``.

# See also
[`grad_TV2`](@ref), [`prox_TV2`](@ref)
"""
function costTV2(M::PowerManifold, x, p::Int=1, Sum::Bool=true)
    Tt = Tuple(power_dimensions(M))
    R = CartesianIndices(Tt)
    d = length(Tt)
    minInd, maxInd = first(R), last(R)
    cost = fill(0.0, Tt)
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i -> (i == k) ? 1 : 0, d)) #k th unit vector
        for i in R # iterate over all pixel
            jF = i + ek # compute forward neighbor
            jB = i - ek # compute backward neighbor
            if all(map(<=, jF.I, maxInd.I)) && all(map(>=, jB.I, minInd.I)) # are neighbors in range?
                cost[i] += costTV2(
                    M.manifold,
                    (x[M, Tuple(jB)...], x[M, Tuple(i)...], x[M, Tuple(jF)...]),
                    p,
                )
            end
        end # i in R
    end # directions
    if p != 1
        cost = (cost) .^ (1 / p)
    end
    if Sum
        return sum(cost)
    else
        return cost
    end
end


text:
include("header.jl")
using LinearAlgebra

@testset "LinearAlgebra" begin

    diagm1(x,y,z) = diagm(1=>x,0=>y,-1=>z)

    function krontest(w1,w2)
        sum(kron(w1,w2))
    end

    function qrtest(w)
        Q,R = qr(w)
        sum(Q) + 2sum(R)
    end

    function lqtest(w)
        Q,L = lq(w)
        sum(Q) + 2sum(L)
    end

    function svdtest(w)
        U,S,V = svd(w)
        sum(U) + 2sum(S) + 3sum(V)
    end

    for T in (Float32, Float64)
        w = randn(T,3,2)
        wt = randn(T,2,3)
        wsquare = I + 0.2*randn(T,3,3)
        wposdef = wsquare' * wsquare
        udg = T[1.,2.,3.]
        dg = T[4.,5.,6.]
        ldg  = T[7.,8.,9.]

        @test gradcheck(adjoint,w)
        @test gradcheck(det,wsquare)
        @test gradcheck(diag,wsquare)
        @test_skip gradcheck(diag,w)             #TODO no support for non-square matrices yet
        @test_skip gradcheck(diagm1,udg,dg,ldg)   #TODO diagm not implemented yet
        @test gradcheck(dot,udg,ldg)
        @test gradcheck(dot,w,copy(w))
        @test gradcheck(inv,wsquare)
        @test gradcheck(krontest,w,copy(w); delta=0.01)
        @test gradcheck(logabsdet,wsquare)
        @test gradcheck(logdet,wposdef)
        @test gradcheck(norm,w,1; args=1)
        @test gradcheck(norm,w,2; args=1)
        @test gradcheck(norm,w,Inf; args=1)
        @test_skip gradcheck(qrtest,w) #TODO iterator error
        @test_skip gradcheck(lqtest,w) #TODO iterator error
        @test_skip gradcheck(svdtest,w) #TODO iterator error
        @test gradcheck(tr,wsquare)
        @test gradcheck(transpose,w)
        @test gradcheck(tril,w)
        @test gradcheck(triu,w)
        @test gradcheck(*,w,wt)
        @test_skip gradcheck(^,wsquare,randn()) #TODO: real matrix powers
        @test_skip gradcheck(^,wsquare,rand(-3:3)) #TODO: integer matrix powers
    end
end