text:
using PyPlot, HDF5, MRIReco, LinearAlgebra

filename = "rawdata_brain_radial_96proj_12ch.h5"
data = permutedims(h5read(filename, "rawdata"),[3,2,1,4])
traj = permutedims(h5read(filename, "trajectory"),[3,2,1])
N = 300
Nc = 12

#############################################
# load data and form Acquisition data object
##############################################
tr = Trajectory(reshape(traj[1:2,:,:],2,:) ./ N, 96, 512, circular=false)
dat = Array{Array{Complex{Float64},2},3}(undef,1,1,1)
dat[1,1,1] = 1.e8.*reshape(data,:,12)
acqData = AcquisitionData(tr, dat, encodingSize=[N,N,1])

################################
# generate coil sensitivity maps
################################
@info "Espirit"
acqDataCart = regrid2d(acqData, (N,N); cgnr_iter=3)
sensitivity = espirit(acqDataCart, (6,6), 30, eigThresh_1=0.02, eigThresh_2=0.98)

##########################
# reference reconstruction
##########################
@info "reference reco"
params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = (N,N)
params[:regularization] = "L2"
params[:λ] = 1.e-2
params[:iterations] = 100
params[:solver] = "cgnr"
params[:senseMaps] = reshape(sensitivity, N, N, 1, Nc)

img_ref = reconstruction(acqData, params).data

##############################
# undersampled reconstructions
##############################
@info "undersampled reco"
rf = [1,2,3,4]
Δ = Vector{Vector{Float64}}(undef,4)
δ = Vector{Vector{Float64}}(undef,4)
img_cg = Vector{Array{ComplexF64,5}}(undef,4)
img_cg1 = Vector{Array{ComplexF64,5}}(undef,4)
img_sc = Vector{Array{ComplexF64,5}}(undef,4)
params[:iterations] = 20
params[:relTol] = 1.e-3
params[:λ] = 1.e-2
for r in rf
  @info "r=$r"
  # undersample profiles
  acqDataSub = convertUndersampledData(sample_kspace(acqData, Float64.(r), "regular"))
  # SENSE reconstruction while monitoring error
  params[:reco] = "multiCoil"
  params[:solverInfo] = SolverInfo(ComplexF64,store_solutions=true)
  img_cg[r] = reconstruction(acqDataSub, params).data
  img_cg1[r] = reshape(params[:solverInfo].x_iter[2],N,N,1,1,1)
  Δ[r] = [norm(vec(img_ref).-params[:solverInfo].x_iter[i])/norm(vec(img_ref)) for i=1:length(params[:solverInfo].x_iter)]
  δ[r] = params[:solverInfo].convMeas ./ params[:solverInfo].convMeas[1]
  # single-channel reconstruction
  params[:reco] = "direct"
  img_sc[r] = reconstruction(acqDataSub, params).data
end

##############
# plot figures
##############
figure("brain reconstruction error",figsize=(8,4))
subplot(1,2,1)
plot(log.(Δ[1]),label="r=1")
plot(log.(Δ[2]),label="r=2")
plot(log.(Δ[3]),label="r=3")
plot(log.(Δ[4]),label="r=4")
xlabel("iteration numbers")
ylabel("log(Δ)")
xticks(collect(0:4:14))
subplot(1,2,2)
plot(log.(δ[1]),label="r=1")
plot(log.(δ[2]),label="r=2")
plot(log.(δ[3]),label="r=3")
plot(log.(δ[4]),label="r=4")
xlabel("iteration numbers")
ylabel("log(δ)")
xticks(collect(0:4:14))
legend()
subplots_adjust(wspace=0.2,hspace=0.0,left=0.07,bottom=0.13,right=0.99,top=0.95)
savefig("Fig4.png",dpi=300)

figure("brain reconstructions",figsize=(4.5,6))
subplot(4,3,1); title("Single coil")
subplot(4,3,2); title("Initial")
subplot(4,3,3); title("Final")
for r in rf
  # plots
  subplot(4,3,(r-1)*3+1)
  imshow(abs.(img_sc[r][end:-1:1,end:-1:1,1,1,1]),cmap="gray")
  ylabel("r=$r")
  xticks([]);yticks([])
  subplot(4,3,(r-1)*3+2)
  imshow(abs.(img_cg1[r][end:-1:1,end:-1:1,1,1,1]),cmap="gray")
  xticks([]);yticks([])
  subplot(4,3,(r-1)*3+3)
  imshow(abs.(img_cg[r][end:-1:1,end:-1:1,1,1,1]),cmap="gray")
  xticks([]);yticks([])
end
subplots_adjust(wspace=0.05,hspace=0.05,left=0.05,bottom=0.0,right=1.0,top=0.95)
savefig("Fig5.png",dpi=300)


text:
using Viznet: canvas
import Viznet
using Compose: CurvePrimitive, Form
using Yao
using BitBasis

export CircuitStyles, CircuitGrid, circuit_canvas, vizcircuit, pt, cm, darktheme!, lighttheme!

module CircuitStyles
    using Compose
    import Viznet
    const r = Ref(0.2)
    const lw = Ref(1pt)
    const textsize = Ref(16pt)
    const paramtextsize = Ref(10pt)
    const fontfamily = Ref("monospace")
    const linecolor = Ref("#000000")
    const gate_bgcolor = Ref("transparent")
    const textcolor = Ref("#000000")
    const scale = Ref(1.0)

    # formats
    struct ComposeSVG end
    struct Tikz end

    # Context is gate dependent information
    Base.@kwdef struct Context
        hspace::Float64
        wspace::Float64
    end

    abstract type Gadget end
    struct Box{FT} <: Gadget
        height::FT
        width::FT
    end
    struct Cross <: Gadget end
    struct Dot <: Gadget end
    struct NDot <: Gadget end
    struct OPlus <: Gadget end
    struct MeasureBox <: Gadget end
    struct Text{FT} <: Gadget
        fontsize::FT
    end
    struct Line <: Gadget end
    get_width(::Cross) = 0.0
    get_width(::Dot) = r[]/3
    get_width(::NDot) = r[]/3
    get_width(::OPlus) = r[]*1.4
    boxsize(b::Gadget, params) = (w = get_width(b); (w, w))
    function boxsize(b::Box, params)
        hspace, wspace = params.hspace, params.wspace
        return b.width + wspace, b.height + hspace
    end
    function boxsize(::MeasureBox, params)
        return 2 * r[], 2 * r[]
    end

    function render(::ComposeSVG, b::Box, params)
        WIDTH, HEIGHT = boxsize(b, params)
        compose(context(), rectangle(-WIDTH/2, -HEIGHT/2, WIDTH, HEIGHT), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(lw[]))
    end

    #G() = compose(context(), rectangle(-r[], -r[], 2*r[], 2*r[]), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(lw[]))
    render(::ComposeSVG, d::Dot, params) = compose(context(), circle(0.0, 0.0, get_width(d)/2), fill(linecolor[]), linewidth(0))
    render(::ComposeSVG, d::NDot, params) = compose(context(), circle(0.0, 0.0, get_width(d)/2), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(lw[]))
    render(::ComposeSVG, d::Cross, params) = compose(context(), xgon(0.0, 0.0, r[], 4), fill(linecolor[]), linewidth(0))
    function render(::ComposeSVG, d::OPlus, params)
        w = get_width(d) / 2
        compose(context(),
                    (context(), circle(0.0, 0.0, w), stroke(linecolor[]), linewidth(lw[]), fill("transparent")),
                    (context(), polygon([(-w, 0.0), (w, 0.0)]), stroke(linecolor[]), linewidth(lw[])),
                    (context(), polygon([(0.0, -w), (0.0, w)]), stroke(linecolor[]), linewidth(lw[]))
               )
    end
    render(::ComposeSVG, ::MeasureBox, params) = compose(context(),
                        rectangle(-r[], -r[], 2*r[], 2*r[]), fill(gate_bgcolor[]), stroke(linecolor[]), linewidth(lw[]),
                        compose(context(), curve((-0.8*r[], 0.5*r[]), (-0.8*r[], -0.6*r[]), (0.8*r[], -0.6*r[]), (0.8*r[], 0.5*r[])), stroke(linecolor[]), linewidth(lw[])),
                        compose(context(), line([(0.0, 0.5*r[]), (0.7*r[], -0.4*r[])]), stroke(linecolor[]), linewidth(lw[])),
        begin
            ns = Viznet.nodestyle(:triangle, fill(linecolor[]); r=0.1*r[], θ=atan(0.7, 0.9))
            Viznet.inner_most_containers(ns) do c
                Viznet.update_locs!(c.form_children, [(0.7*r[], -0.4*r[])])
            end
            ns
        end
        )

    render(::ComposeSVG, ::Line, params) = compose(context(), line(), stroke(linecolor[]), linewidth(lw[]))
    render(::ComposeSVG, t::Text, params) = compose(context(), text(0.0, 0.0, "", hcenter, vcenter), fontsize(t.fontsize*scale[]), fill(textcolor[]), font(fontfamily[]))

    Base.@kwdef struct GateStyles
        g = Box(2*r[], 2*r[])
        c = Dot()
        x = Cross()
        nc = NDot()
        not = OPlus()
        measure = MeasureBox()

        # other styles
        line = Line()
        text = Text(textsize[])
        smalltext = Text(paramtextsize[])
    end
end

struct CircuitGrid
    backend
    frontier::Vector{Int}
    w_depth::Float64
    w_line::Float64
    gatestyles::CircuitStyles.GateStyles
end

nline(c::CircuitGrid) = length(c.frontier)
depth(c::CircuitGrid) = frontier(c, 1, nline(c))
Base.getindex(c::CircuitGrid, i, j) = (c.w_depth*i, c.w_line*j)
Base.typed_vcat(c::CircuitGrid, ij1, ij2) = (c[ij1...], c[ij2...])

function CircuitGrid(nline::Int; w_depth=1.0, w_line=1.0, gatestyles=CircuitStyles.GateStyles(), backend=CircuitStyles.ComposeSVG())
    CircuitGrid(backend, zeros(Int, nline), w_depth, w_line, gatestyles)
end

function frontier(c::CircuitGrid, args...)
    maximum(i->c.frontier[i], min(args..., nline(c)):max(args..., 1))
end

function _draw!(c::CircuitGrid, loc_brush_texts)
    isempty(loc_brush_texts) && return
    # a loc can be a integer, or a range
    loc_brush_texts = sort(loc_brush_texts, by=x->maximum(x[1]))
    locs = Iterators.flatten(getindex.(loc_brush_texts, 1)) |> collect

    # get the gate width and the circuit depth to draw
    boxwidths, boxheights = Float64[], Float64[]
    for (j, b, txt) in loc_brush_texts
        length(j) == 0 && continue
        wspace, fontsize = text_width_and_size(txt)
        hspace = (maximum(j)-minimum(j)) * c.w_line
        context = CircuitStyles.Context(; wspace=wspace,
                                        hspace=hspace)
        boxwidth, boxheight = CircuitStyles.boxsize(b, context)
        push!(boxwidths, boxwidth)
        push!(boxheights, boxheight)
    end
    max_width = maximum(boxwidths)
    ncolumn = max(1, ceil(Int, max_width/c.w_depth + 0.2))  # 0.1 is the minimum gap between two columns

    ipre = frontier(c, minimum(locs):maximum(locs)...)
    i = ipre + ncolumn/2
    local jpre
    for (k, ((j, b, txt), boxheight)) in enumerate(zip(loc_brush_texts, boxheights))
        length(j) == 0 && continue
        wspace, fontsize = text_width_and_size(txt)
        hspace = (maximum(j)-minimum(j)) * c.w_line
        jmid = (minimum(j)+maximum(j))/2
        context = CircuitStyles.Context(; wspace=wspace,
                                        hspace=hspace)
        CircuitStyles.render(c.backend, b, context) >> c[i, jmid]
        CircuitStyles.render(c.backend, CircuitStyles.Text(fontsize), context) >> (c[i, jmid], txt)
        # use line to connect blocks in the same gate
        if k!=1
            CircuitStyles.render(c.backend, c.gatestyles.line, context) >> c[(i, jmid-boxheight/2/c.w_line); (i, jpre)]
        end
        jpre = jmid + boxheight/2/c.w_line
    end

    #jmin, jmax = min(locs..., nline(c)), max(locs..., 1)
    # connect horizontal lines
    for (width, (j, b, txt)) in zip(boxwidths, loc_brush_texts)
        for jj in j
            CircuitStyles.render(CircuitStyles.ComposeSVG(), c.gatestyles.line, nothing) >> c[(c.frontier[jj], jj); (i-width/2/c.w_depth, jj)]
            CircuitStyles.render(CircuitStyles.ComposeSVG(), c.gatestyles.line, nothing) >> c[(i+width/2/c.w_depth, jj); (ipre+ncolumn, jj)]
            c.frontier[jj] = ipre + ncolumn
        end
    end
    for j in setdiff(minimum(locs):maximum(locs), locs)
        CircuitStyles.render(CircuitStyles.ComposeSVG(), c.gatestyles.line, nothing) >> c[(c.frontier[j], j); (ipre+ncolumn, j)]
        c.frontier[j] = ipre + ncolumn
    end
end

function text_width_and_size(text)
    W = maximum(x->textwidth(x), split(text, "\n"))
    fontsize = W > 3 ? CircuitStyles.paramtextsize[] : CircuitStyles.textsize[]
    # -2 because the gate has a default size
    width = max(W - 4, 0) * fontsize.value * 0.025  # mm to cm
    return width, fontsize
end

function initialize!(c::CircuitGrid; starting_texts, starting_offset)
    starting_texts !== nothing && for j=1:nline(c)
        CircuitStyles.render(c.backend, c.gatestyles.text, nothing) >> (c[starting_offset, j], string(starting_texts[j]))
    end
end

function finalize!(c::CircuitGrid; show_ending_bar, ending_offset, ending_texts)
    i = frontier(c, 1, nline(c))
    for j=1:nline(c)
        show_ending_bar && c.gatestyles.line >> c[(i, j-0.2); (i, j+0.2)]
        CircuitStyles.render(c.backend, c.gatestyles.line, nothing) >> c[(i, j); (c.frontier[j], j)]
        ending_texts !== nothing && c.gatestyles.text >> (c[i+ending_offset, j], string(ending_texts[j]))
    end
end

# elementary
function draw!(c::CircuitGrid, p::PrimitiveBlock, address, controls)
    bts = length(controls)>=1 ? get_cbrush_texts(c, p) : get_brush_texts(c, p)
    _draw!(c, [controls..., (getindex.(Ref(address), occupied_locs(p)), bts[1], bts[2])])
end
function draw!(c::CircuitGrid, p::Daggered{<:PrimitiveBlock}, address, controls)
    bts = length(controls)>=1 ? get_cbrush_texts(c, content(p)) : get_brush_texts(c, content(p))
    _draw!(c, [controls..., (getindex.(Ref(address), occupied_locs(p)), bts[1], bts[2]*"'")])
end

function draw!(c::CircuitGrid, p::Scale, address, controls)
    fp = YaoBlocks.factor(p)
    if !(abs(fp) ≈ 1)
        error("can not visualize non-phase factor.")
    end
    draw!(c, YaoBlocks.phase(angle(fp)), [first(address)], controls)
    draw!(c, p.content, address, controls)
end
# Special primitive gates
function draw!(c::CircuitGrid, ::I2Gate, address, controls)
    return
end
function draw!(c::CircuitGrid, ::IdentityGate, address, controls)
    return
end
function draw!(c::CircuitGrid, p::ConstGate.SWAPGate, address, controls)
    bts = [(c.gatestyles.x, ""), (c.gatestyles.x, "")]
    _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...])
end
function draw!(c::CircuitGrid, p::ConstGate.CNOTGate, address, controls)
    bts = [(c.gatestyles.c, ""), (c.gatestyles.x, "")]
    _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...])
end
function draw!(c::CircuitGrid, p::ConstGate.CZGate, address, controls)
    bts = [(c.gatestyles.c, ""), (c.gatestyles.c, "")]
    _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...])
end
function draw!(c::CircuitGrid, p::ConstGate.ToffoliGate, address, controls)
    bts = [(c.gatestyles.c, ""), (c.gatestyles.c, ""), (c.gatestyles.x, "")]
    _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...])
end

# composite
function draw!(c::CircuitGrid, p::ChainBlock, address, controls)
    draw!.(Ref(c), subblocks(p), Ref(address), Ref(controls))
end

function draw!(c::CircuitGrid, p::PutBlock, address, controls)
    locs = [address[i] for i in p.locs]
    draw!(c, p.content, locs, controls)
end

function draw!(c::CircuitGrid, m::YaoBlocks.Measure, address, controls)
    @assert length(controls) == 0
    if m.postprocess isa RemoveMeasured
        error("can not visualize post-processing: `RemoveMeasured`.")
    end
    if !(m.operator isa ComputationalBasis)
        error("can not visualize measure blocks for operators")
    end
    locs = m.locations isa AllLocs ? collect(1:nqudits(m)) : [address[i] for i in m.locations]
    for (i, loc) in enumerate(locs)
        _draw!(c, [(loc, c.gatestyles.measure, "")])
        if m.postprocess isa ResetTo
            val = readbit(m.postprocess.x, i)
            _draw!(c, [(loc, c.gatestyles.g, val == 1 ? "P₁" : "P₀")])
        end
    end
end

function draw!(c::CircuitGrid, cb::ControlBlock{GT,C}, address, controls) where {GT,C}
    ctrl_locs = [address[i] for i in cb.ctrl_locs]
    locs = [address[i] for i in cb.locs]
    mycontrols = [(loc, (bit == 1 ? c.gatestyles.c : c.gatestyles.nc), "") for (loc, bit)=zip(ctrl_locs, cb.ctrl_config)]
    draw!(c, cb.content, locs, [controls..., mycontrols...])
end

for B in [:LabelBlock, :GeneralMatrixBlock]
    @eval function draw!(c::CircuitGrid, cb::$B, address, controls)
        length(address) == 0 && return
        #is_continuous_chunk(address) || error("address not continuous in a block marked as continous.")
        _draw!(c, [controls..., (address, c.gatestyles.g, string(cb))])
    end
end
for GT in [:KronBlock, :RepeatedBlock, :CachedBlock, :Subroutine, :(YaoBlocks.AD.NoParams)]
    @eval function draw!(c::CircuitGrid, p::$GT, address, controls)
        draw!(c, YaoBlocks.Optimise.to_basictypes(p), address, controls)
    end
end
# function draw!(c::CircuitGrid, b::GeneralMatrixBlock{GT}, address, controls) where {GT}
#     length(address) == 0 && return
#     nline = maximum(address)-minimum(address)
#     g = CircuitStyles.MULTIGATE(nline*c.w_line, get_textwidth(b.tag))
#     _draw!(c, [controls..., (address, g, b.tag)])
# end

for (GATE, SYM) in [(:XGate, :Rx), (:YGate, :Ry), (:ZGate, :Rz)]
    @eval get_brush_texts(c, b::RotationGate{D,T,<:$GATE}) where {D,T} = (c.gatestyles.g, "$($(SYM))($(pretty_angle(b.theta)))")
end
get_brush_texts(c, b::EasyBuild.FSimGate) = (c.gatestyles.g, "FSim($(pretty_angle(b.theta)), $(pretty_angle(b.phi)))")
get_brush_texts(c, ::EasyBuild.SqrtWGate) = (c.gatestyles.g, "√W")
get_brush_texts(c, ::EasyBuild.SqrtXGate) = (c.gatestyles.g, "√X")
get_brush_texts(c, ::EasyBuild.SqrtYGate) = (c.gatestyles.g, "√Y")

pretty_angle(theta) = string(theta)
function pretty_angle(theta::AbstractFloat)
    c = ZXCalculus.continued_fraction(theta/π, 10)
    if c.den < 100
        res = if c.num == 1
            "π"
        elseif c.num==0
            "0"
        elseif c.num==-1
            "-π"
        else
            "$(c.num)π"
        end
        if c.den != 1
            res *= "/$(c.den)"
        end
        res
    else
        "$(round(theta; digits=2))"
    end
end

get_brush_texts(c, ::ConstGate.SdagGate) = (c.gatestyles.g, "S'")
get_brush_texts(c, ::ConstGate.TdagGate) = (c.gatestyles.g, "T'")
get_brush_texts(c, ::ConstGate.PuGate) = (c.gatestyles.g, "P+")
get_brush_texts(c, ::ConstGate.PdGate) = (c.gatestyles.g, "P-")
get_brush_texts(c, ::ConstGate.P0Gate) = (c.gatestyles.g, "P₀")
get_brush_texts(c, ::ConstGate.P1Gate) = (c.gatestyles.g, "P₁")
get_brush_texts(c, b::PrimitiveBlock) = (c.gatestyles.g, string(b))
get_brush_texts(c, b::TimeEvolution) = (c.gatestyles.g, string(b))
get_brush_texts(c, b::ShiftGate) = (c.gatestyles.g, "ϕ($(pretty_angle(b.theta)))")
get_brush_texts(c, b::PhaseGate) = (c.gatestyles.c, "    $(pretty_angle(b.theta))\n")
function get_brush_texts(c, b::T) where T<:ConstantGate
    namestr = string(T.name.name)
    if endswith(namestr, "Gate")
        namestr = namestr[1:end-4]
    end
    # Fix!
    (c.gatestyles.g, namestr)
end

get_cbrush_texts(c, b::PrimitiveBlock) = get_brush_texts(c, b)
get_cbrush_texts(c, ::XGate) = (c.gatestyles.not, "")
get_cbrush_texts(c, ::ZGate) = (c.gatestyles.c, "")

# front end
plot(blk::AbstractBlock; kwargs...) = vizcircuit(blk; kwargs...)
function vizcircuit(blk::AbstractBlock; w_depth=0.85, w_line=0.75, scale=1.0, show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, ending_texts=nothing, ending_offset=0.3, graphsize=1.0, gatestyles=CircuitStyles.GateStyles(), backend=CircuitStyles.ComposeSVG())
    CircuitStyles.scale[] = scale
    img = circuit_canvas(nqubits(blk); w_depth, w_line, show_ending_bar, starting_texts, starting_offset, ending_texts, ending_offset, graphsize, gatestyles, backend) do c
        blk >> c
    end
    CircuitStyles.scale[] = 1.0
    return img |> rescale(scale)
end

function circuit_canvas(f, nline::Int; backend=CircuitStyles.ComposeSVG(), w_depth=0.85, w_line=0.75, show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, ending_texts=nothing, ending_offset=0.3, graphsize=1.0, gatestyles=CircuitStyles.GateStyles())
    c = CircuitGrid(nline; w_depth, w_line, gatestyles, backend)
    Viznet.empty_cache!()
    initialize!(c; starting_texts, starting_offset)
    f(c)
    finalize!(c; show_ending_bar, ending_texts, ending_offset)
    # flush!
    g = compose(context(),
        Viznet.flush!(Viznet.TEXT_CACHE)...,
        Viznet.flush!(Viznet.EDGE_CACHE)...,
        Viznet.flush!(Viznet.NODE_CACHE)...,
        )
    a, b = (depth(c)+1)*w_depth, nline*w_line
    Compose.set_default_graphic_size(a*2.5*graphsize*cm, b*2.5*graphsize*cm)
    compose(context(0.5/a, -0.5/b, 1/a, 1/b), g)
end

Base.:>>(blk::AbstractBlock, c::CircuitGrid) = draw!(c, blk, collect(1:nqudits(blk)), [])
Base.:>>(blk::Function, c::CircuitGrid) = blk(nline(c)) >> c

function rescale(factor)
    a, b = Compose.default_graphic_width, Compose.default_graphic_height
    Compose.set_default_graphic_size(a*factor, b*factor)
    graph -> compose(context(), graph)
end

vizcircuit(; kwargs...) = c->vizcircuit(c; kwargs...)

function darktheme!()
    const CircuitStyles.linecolor[] = "#FFFFFF"
    const CircuitStyles.textcolor[] = "#FFFFFF"
end

function lighttheme!()
    const CircuitStyles.linecolor[] = "#000000"
    const CircuitStyles.textcolor[] = "#000000"
end

ext: jl


text:
export @asmcall

@generated function _asmcall(::Val{asm}, ::Val{constraints}, ::Val{side_effects},
                             ::Val{rettyp}, ::Val{argtyp}, args...) where
                            {asm, constraints, side_effects, rettyp, argtyp}
    # create a function
    llvm_rettyp = convert(LLVMType, rettyp)
    llvm_argtyp = LLVMType[convert.(LLVMType, [argtyp.parameters...])...]
    llvm_f, llvm_ft = create_function(llvm_rettyp, llvm_argtyp)

    inline_asm = InlineAsm(llvm_ft, String(asm), String(constraints), side_effects)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        val = call!(builder, inline_asm, collect(parameters(llvm_f)))
        if rettyp == Nothing
            ret!(builder)
        else
            ret!(builder, val)
        end
    end

    call_function(llvm_f, rettyp, argtyp,
                  Expr(:tuple, (:(args[$i]) for i in 1:length(args))...))
end

"""
    @asmcall asm::String [constraints::String] [side_effects::Bool=false]
             rettyp=Nothing argtyp=Tuple{} args...

Call some inline assembly `asm`, optionally constrained by `constraints` and denoting other
side effects in `side_effects`, specifying the return type in `rettyp` and types of
arguments as a tuple-type in `argtyp`.
"""
:(@asmcall)

macro asmcall(asm::String, constraints::String, side_effects::Bool,
              rettyp::Union{Expr,Symbol}=:(Nothing), argtyp::Expr=:(Tuple{}), args...)
    asm_val = Val{Symbol(asm)}()
    constraints_val = Val{Symbol(constraints)}()
    return esc(:(LLVM.Interop._asmcall($asm_val, $constraints_val,
                                       Val{$side_effects}(), Val{$rettyp}(), Val{$argtyp}(),
                                       $(args...))))
end

# shorthand: no side_effects
macro asmcall(asm::String, constraints::String,
              rettyp::Union{Expr,Symbol}=:(Nothing), argtyp::Expr=:(Tuple{}), args...)
    esc(:(LLVM.Interop.@asmcall $asm $constraints false $rettyp $argtyp $(args...)))
end

# shorthand: no side_effects or constraints
macro asmcall(asm::String,
              rettyp::Union{Expr,Symbol}=:(Nothing), argtyp::Expr=:(Tuple{}), args...)
    esc(:(LLVM.Interop.@asmcall $asm "" $rettyp $argtyp $(args...)))
end


text:
import Pkg

if get(ENV, "CI", nothing) === nothing
    Pkg.activate("..")
end

using Documenter, Echidna

DocMeta.setdocmeta!(Echidna, :DocTestSetup,
                    :(using Echidna); recursive=true)

makedocs(
    modules = [Echidna],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    clean = true,
    sitename = "Echidna.jl",
    authors = "Peter Hawthorne",
    checkdocs = :exports,
    pages = Any["index.md"]
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(
        repo = "github.com/phawthorne/Echidna.git"
    )
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
# this test is designed to test the ability of the package to use a state space that is not an AbstractArray.
# using Revise
using Test, Random, Setfield
using PseudoArcLengthContinuation
const PALC = PseudoArcLengthContinuation
####################################################################################################
# We start with a simple Fold problem
using LinearAlgebra
function F0(x::Vector, r)
	out = r .+  x .- x.^3
end

opt_newton0 = PALC.NewtonPar(tol = 1e-11, verbose = true)
	out0, hist, flag = @time PALC.newton(
		x -> F0(x, 1.),
		[0.8],
		opt_newton0)

opts_br0 = PALC.ContinuationPar(dsmin = 0.001, dsmax = 0.021, ds= -0.01, pMax = 4.1, pMin = -1., newtonOptions = setproperties(opt_newton0; maxIter = 70, tol = 1e-8), detectBifurcation = 0, maxSteps = 150, detectFold = true)

	br0, u1 = @time PALC.continuation(
		F0,
		out0, 1.0,
		opts_br0, printSolution = (x,p) -> x[1])

# plotBranch(br0);title!("a")

outfold, hist, flag = @time PALC.newtonFold(
	F0,
	(x, r) -> diagm(0 => 1 .- 3 .* x.^2),
	(x, r) -> diagm(0 => 1 .- 3 .* x.^2),
	(x, r, v1, v2) -> -6 .* x .* v1 .* v2,
	br0, 2, #index of the fold point
	opts_br0.newtonOptions)
		flag && printstyled(color=:red, "--> We found a Fold Point at α = ",outfold.p, ", from ", br0.foldpoint[1][3], "\n")

####################################################################################################
# Here is a more involved example

function Fb(x::BorderedArray, r, s = 1.0)
	BorderedArray(r .+  s .* x.u .- (x.u).^3, x.p - 0.0)
end

# there is no finite differences defined, so we need to provide a linearsolve
# we could also have used GMRES. We define a custom Jacobian which will be used for evaluation and jacobian inverse
struct Jacobian
	x
	r
	s
end

# We express the jacobian operator
function (J::Jacobian)(dx)
	BorderedArray( (J.s .- 3 .* ((J.x).u).^2) .* dx.u, dx.p)
end

struct linsolveBd <: PALC.AbstractBorderedLinearSolver end

function (l::linsolveBd)(J, dx)
	x = J.x
	r = J.r
	out = BorderedArray(dx.u ./ (J.s .- 3 .* (x.u).^2), dx.p)
	out, true, 1
end

sol = BorderedArray([0.8], 0.0)

opt_newton = PALC.NewtonPar(tol = 1e-11, verbose = true, linsolver = linsolveBd())
out, hist, flag = @time PALC.newton(
		x -> Fb(x, 1., 1.),
		x -> Jacobian(x, 1., 1.),
		sol,
		opt_newton)

opts_br = PALC.ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= -0.01, pMax = 4.1, pMin = -1., newtonOptions = setproperties(opt_newton; maxIter = 70, tol = 1e-8), detectBifurcation = 0, maxSteps = 150)

	br, u1 = @time PALC.continuation(
		(x, r) -> Fb(x, r, 1.),
		(x, r) -> Jacobian(x, r, 1.),
		out, 1.,
		opts_br; printSolution = (x,p) -> x.u[1])

# plotBranch(br);title!("")

outfold, hist, flag = @time PALC.newtonFold(
	(x, r) -> Fb(x, r),
	(x, r) -> Jacobian(x, r, 1.),
	(x, r) -> Jacobian(x, r, 1.),
	(x, r, v1, v2) -> BorderedArray(-6 .* x.u .* v1.u .* v2.u, 0.),
	br, 1, #index of the fold point
	opts_br.newtonOptions)
		flag && printstyled(color=:red, "--> We found a Fold Point at α = ", outfold.p, ", from ", br.foldpoint[1][3],"\n")

outfoldco, hist, flag = @time PALC.continuationFold(
	(x, r, s) -> Fb(x, r, s),
	(x, r, s) -> Jacobian(x, r, s),
	(x, r, s) -> Jacobian(x, r, s),
	s -> ((x, r, v1, v2) -> BorderedArray(-6 .* x.u .* v1.u .* v2.u, 0.)),
	br, 1,
	1.0, plot = false, opts_br)

# try with newtonDeflation
printstyled(color=:green, "--> test with Newton deflation 1")
deflationOp = DeflationOperator(2.0, (x, y) -> dot(x, y), 1.0, [zero(sol)])
soldef0 = BorderedArray([0.1], 0.0)
soldef1, _, _ = @time PALC.newton(
	x -> Fb(x, 0., 1.),
	x -> Jacobian(x, 0., 1.),
	soldef0,
	opt_newton, deflationOp)

push!(deflationOp, soldef1)

Random.seed!(1231)
printstyled(color=:green, "--> test with Newton deflation 2")
soldef2, _, _ = @time PALC.newton(
	x -> Fb(x, 0., 1.),
	x -> Jacobian(x, 0., 1.),
	rmul!(soldef0,rand()),
	opt_newton, deflationOp)

# test indexing
deflationOp[1]

# test length
length(deflationOp)

# test pop!
pop!(deflationOp)

####################################################################################################
using KrylovKit

function Fr(x::RecursiveVec, r, s = 1.0)
	out = similar(x)
	for ii=1:length(x)
		out[ii] .= r .+  s .* x[ii] .- x[ii].^3
	end
	out
end

# there is no finite differences defined, so we need to provide a linearsolve
# we could also have used GMRES. We define a custom Jacobian which will be used for evaluation and jacobian inverse
struct JacobianR
	x
	s
end

# We express the jacobian operator
function (J::JacobianR)(dx)
	out = similar(dx)
	for ii=1:length(out)
		out[ii] .= (J.s .- 3 .* (J.x[ii]).^2) .* dx[ii]
	end
	return out
end

struct linsolveBd_r <: PALC.AbstractBorderedLinearSolver end

function (l::linsolveBd_r)(J, dx)
	x = J.x
	out = similar(dx)
	for ii=1:length(out)
		out[ii] .= dx[ii] ./ (J.s .- 3 .* (x[ii]).^2)
	end
	out, true, 1
end

opt_newton0 = PALC.NewtonPar(tol = 1e-11, verbose = false, linsolver = linsolveBd_r())
	out0, hist, flag = @time PALC.newton(
		x -> Fr(x, 1.),
		x -> JacobianR(x, 1),
		RecursiveVec([1 .+ 0.1*rand(10) for _=1:2]),
		opt_newton0)

opts_br0 = PALC.ContinuationPar(dsmin = 0.001, dsmax = 0.2, ds= -0.01, pMin = -1.1, pMax = 1.1, newtonOptions = opt_newton0)

br0, u1 = @time PALC.continuation(
	(x, r) -> Fr(x, r),
	(x, r) -> JacobianR(x, 1),
	out0, 0.9, opts_br0;
	plot = false,
	printSolution = (x,p) -> x[1][1])

br0, u1 = @time PALC.continuation(
	(x, r) -> Fr(x, r),
	(x, r) -> JacobianR(x, 1),
	out0, 0.9, opts_br0;
	tangentAlgo = BorderedPred(),
	plot = false,
	printSolution = (x,p) -> x[1][1])

outfold, hist, flag = @time PALC.newtonFold(
	(x, r) -> Fr(x, r),
	(x, r) -> JacobianR(x, 1),
	(x, r) -> JacobianR(x, 1),
	(x, r, v1, v2) -> RecursiveVec([-6 .* x[ii] .* v1[ii] .* v2[ii] for ii=1:length(x)]),
	br0, 1, #index of the fold point
	opt_newton0)
		flag && printstyled(color=:red, "--> We found a Fold Point at α = ", outfold.p, ", from ", br0.foldpoint[1][3],"\n")

outfoldco, hist, flag = @time PALC.continuationFold(
	(x, r, s) -> Fr(x, r, s),
	(x, r, s) -> JacobianR(x, s),
	(x, r, s) -> JacobianR(x, s),
	s -> ((x, r, v1, v2) -> RecursiveVec([-6 .* x[ii] .* v1[ii] .* v2[ii] for ii=1:length(x)])),
	br0, 1,	1.0, opts_br0;
	tangentAlgo = SecantPred(), plot = false)

outfoldco, hist, flag = @time PALC.continuationFold(
	(x, r, s) -> Fr(x, r, s),
	(x, r, s) -> JacobianR(x, s),
	(x, r, s) -> JacobianR(x, s),
	s -> ((x, r, v1, v2) -> RecursiveVec([-6 .* x[ii] .* v1[ii] .* v2[ii] for ii=1:length(x)])),
	br0, 1,	1.0, opts_br0;
	tangentAlgo = BorderedPred(), plot = false)

# try with newtonDeflation
printstyled(color=:green, "--> test with Newton deflation 1")
deflationOp = DeflationOperator(2.0, (x, y) -> dot(x, y), 1.0, [(out0)])
soldef1, _, _ = @time PALC.newton(
	x -> Fr(x, 0.),
	x -> JacobianR(x, 0.),
	out0,
	opt_newton0, deflationOp)

push!(deflationOp, soldef1)


text:
function demo(outdir::AbstractString="results")
  if outdir == "results"
    outdir = joinpath(dirname(pathof(DCEMRI)), "..", "demo", "results")
  end
  outdir = abspath(outdir)
  isdir(outdir) || mkdir(outdir)
  println("Processing in vivo data ...")

  datafile = joinpath(dirname(pathof(DCEMRI)), "..", "demo", "invivo.mat")
  # run the model
  results = fitdata(datafile=datafile, outfile="$outdir/results.mat", models=[2])

  # plot the results
  println("Plotting results ...")
  makeplots(results; outdir=outdir)
  println("Results can be found in $outdir")
  println("Demo run complete.")
end


text:
export
    MtlCommandBuffer, commit!

const MTLCommandBuffer = Ptr{MtCommandBuffer}

"""
    MtlCommandBuffer(queue::MtlCommandQueue; unretained_references=false)

A container that stores encoded commands for the GPU to execute.

If `unretained_references=true` it doesn't hold strong references to any objects
required to execute the command buffer

[Metal Docs](https://developer.apple.com/documentation/metal/mtlcommandbuffer?language=objc)

[Unretained references](https://developer.apple.com/documentation/metal/mtlcommandqueue/1508684-commandbufferwithunretainedrefer?language=objc)
"""
mutable struct MtlCommandBuffer
    handle::MTLCommandBuffer
    queue::MtlCommandQueue
end

Base.convert(::Type{MTLCommandBuffer}, q::MtlCommandBuffer) = q.handle
Base.unsafe_convert(::Type{MTLCommandBuffer}, q::MtlCommandBuffer) = convert(MTLCommandBuffer, q.handle)

Base.:(==)(a::MtlCommandBuffer, b::MtlCommandBuffer) = a.handle == b.handle
Base.hash(q::MtlCommandBuffer, h::UInt) = hash(q.handle, h)

# Constructor
function MtlCommandBuffer(queue::MtlCommandQueue; retain_references = true)
    if retain_references
        handle = mtNewCommandBuffer(queue)
    else
        handle = mtNewCommandBufferWithUnretainedReferences(queue)
    end

    obj = MtlCommandBuffer(handle, queue)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(cmdbuf::MtlCommandBuffer)
    if cmdbuf.handle !== C_NULL
        mtCommandBufferRelease(cmdbuf)
    end
end

## Properties
device(l::MtlCommandBuffer) = MtlDevice(Bool, mtCommandBufferDevice(l))
queue(l::MtlCommandBuffer) = l.queue
retained_references(l::MtlCommandBuffer) = mtCommandBufferRetainedReferences(l)
label(l::MtlCommandBuffer) = unsafe_string_maybe(mtCommandBufferLabel(l))

##STATUS
status(q::MtlCommandBuffer) = mtCommandBufferStatus(q)
kernel_start_time(q::MtlCommandBuffer) = mtCommandBufferKernelStartTime(q)
kernel_end_time(q::MtlCommandBuffer) = mtCommandBufferKernelEndTime(q)
gpu_start_time(q::MtlCommandBuffer) = mtCommandBufferGPUStartTime(q)
gpu_end_time(q::MtlCommandBuffer) = mtCommandBufferGPUEndTime(q)
execution_error(q::MtlCommandBuffer) = NsError_maybe(mtCommandBufferError(q))

# Operations
"""
    enqueue!(q::MtlCommandBuffer)

Enqueueing a command buffer reserves a place for the command buffer on the command
queue without committing the command buffer for execution. When this command buffer
is later committed, it keeps its position in the queue. You enqueue command buffers
so that you can create multiple command buffers with a fixed order of execution without
encoding the command buffers serially. You can use other threads to encode commands
into the command buffers and those threads can complete in any order.

[enqueue](https://developer.apple.com/documentation/metal/mtlcommandbuffer/1443019-enqueue?language=objc)
"""
enqueue!(q::MtlCommandBuffer) = mtCommandBufferEnqueue(q)

"""
    waitUntilScheduled(commandBuffer)

Blocks execution of the current thread until the command buffer is scheduled.
This method returns after the command buffer has been scheduled and all code
blocks registered by addScheduledHandler: have been invoked. A command buffer
is considered scheduled after all its dependencies are resolved, and it is sent
to the GPU for execution.
"""
waitUntilScheduled(q::MtlCommandBuffer) = mtCommandBufferWaitUntilScheduled(q)

"""
    waitUntilCompleted(cmdbuf::MtlCommandBuffer)
    Base.wait(cmdbuf::MtlCommandBuffer)

Blocks execution of the current thread until execution of the command
buffer is completed.
This method returns after the command buffer is completed and all code
blocks registered by addCompletedHandler: are invoked.
"""
waitUntilCompleted(q::MtlCommandBuffer) = mtCommandBufferWaitUntilCompleted(q)
Base.wait(q::MtlCommandBuffer) = waitUntilCompleted(q)

function show(io::IO, ::MIME"text/plain", q::MtlCommandBuffer)
    println(io, "MtlCommandBuffer:")
    println(io, " handle  : ", q.handle)
    println(io, " queue   : ", q.queue)
    print(io, "  status : ", status(q))
end

## One-Api like stuff
"""
    encode_signal!(buf::MtlCommandBuffer, ev::MtlEvent, val::UInt)

Encodes a command that signals the given event, updating it to a new value.

You can't encode a signal event if the command buffer has an active command encoder.
Metal signals the event after all commands scheduled prior to this command
have finished executing. If the new event value is greater than the event's
current value, Metal updates the event's value to the new value. Commands
waiting on the event are allowed to run if the new value is equal to or
greater than the value for which they are waiting. For shared events, this
update similarly triggers notification handlers waiting on the event.
"""
append_signal!(buf::MtlCommandBuffer, ev::MtlEvent, val::UInt) =
    mtEncodeSignalEvent(buf, ev, val)

"""
    encode_wait!(buf::MtlCommandBuffer, ev::MtlEvent, val::UInt)

Encodes a command that blocks the execution of the command buffer
until the given event reaches the given value.

You can't encode a signal event if the command buffer has an active command encoder.
When the device object reaches the command for the wait event, the
device object waits until the event is signaled with a value
equal to or larger than the provided value. While waiting, the
GPU executes commands that appear earlier than the wait command,
but doesn't start any commands that appear after it. Execution continues
immediately if the event already has an equal or larger value.
"""
append_wait!(buf::MtlCommandBuffer, ev::MtlEvent, val::UInt) =
    mtEncodeWaitForEvent(buf, ev, val)

##

commit!(q::MtlCommandBuffer) = mtCommandBufferCommit(q)
commit!(q::Vector{MtlCommandBuffer}) = map(mtCommandBufferCommit, q)

function commit!(f::Base.Callable, queue::MtlCommandQueue; kwargs...)
    cmdbuf = MtlCommandBuffer(f, queue; kwargs...)
    commit!(cmdbuf)
    return cmdbuf
end

function MtlCommandBuffer(f::Base.Callable, queue::MtlCommandQueue; kwargs...)
    cmdbuf = MtlCommandBuffer(queue; kwargs...)
    f(cmdbuf)
    return cmdbuf
end


text:
############################################################
###  ASE Neighborlist implementation
############################################################

"""### type ASENeighborList <: AbstractNeighborList

This makes available the functionality of the ASE neighborlist implementation.
The neighborlist itself is actually stored in `ASEAtoms.po`, but
attaching an `ASENeighborList` will indicate this.

#### Keyword arguments:

 * skin = 0.3
 * sorted
 * self_interaction
 * bothways
"""
type ASENeighborList #  <: AbstractNeighborList
    po::PyObject
    cutoffs::Vector{Float64}
end

# # default constructor from a list of cut-offs
# ASENeighborList(cutoffs::Vector{Float64}; kwargs...) =
#     ASENeighborList(ase_neiglist.NeighborList(cutoffs; kwargs...))

# constructor from an ASEAtoms object, with a single cut-off
# this generates a list of cut-offs, then the neighborulist, then
# builds the neighborlist from the ASEAtoms object
ASENeighborList(atm::ASEAtoms, cutoff::Float64; kwargs...) =
    ASENeighborList(atm, cutoff * ones(length(atm)); kwargs...)

# constructor from an ASEAtoms object, with multiple cutoffs
# this also builds the neighborlist
function ASENeighborList(atm::ASEAtoms, cutoffs::Vector{Float64};
                         bothways=true, self_interaction=false, kwargs...)
    po = ase_neiglist.NeighborList(cutoffs * 0.5;
                                   bothways=bothways,
                                   self_interaction=self_interaction,
                                   kwargs...)
    nlist = ASENeighborList(po, cutoffs)
    update!(nlist, atm)
    return nlist
end

# # regain the cutoffs vector
# _get_cutoffs_ref_(nlist::ASENeighborList) =
#     PyArray(nlist.po["cutoffs"])

"""`update!(nlist::ASENeighborList, atm::ASEAtoms)`

checks whether the atom positions have moved by more than the skin
and rebuilds the list if they have done so.
"""
update!(nlist::ASENeighborList, atm::ASEAtoms) =  nlist.po[:update](atm.po)

"""`build!(nlist::ASENeighborList, atm::ASEAtoms)`

force rebuild of the neighborlist
"""
build!(nlist::ASENeighborList, atm::ASEAtoms) = nlist.po[:build](atm.po)

"""`get_neighbors(n::Integer, neiglist::ASENeighborList) -> (indices, offsets)`

Return neighbors and offsets of atom number n. The positions of the neighbor
atoms can then calculated like this: (python code)

```
          indices, offsets = nl.get_neighbors(42)
          for i, offset in zip(indices, offsets):
              print(atoms.positions[i] + dot(offset, atoms.get_cell()))
```

If `get_neighbors(a)` gives atom b as a neighbor,
    then `get_neighbors(b)` will not return a as a neighbor, unless
    `bothways=True` was used.
"""
function get_neighbors(n::Integer, neiglist::ASENeighborList)
    indices, offset = neiglist.po[:get_neighbors](n-1)
    indices .+= 1
    return (indices, offset)
end
# "alias for `get_neighbors`"
# get_neighbors = get_neighbors
# "alias for `get_neighbors`"
# neighbors = get_neighbors
"alias for `get_neighbors`"
neighbors = get_neighbors

"""`get_neighbors(n::Integer, neiglist::ASENeighborList, atm::ASEAtoms)
              -> (indices::Vector{Int}, s::Vector{Float64}, r::Matrix{Float64})`

 * `indices`: indices of atom positions
 * `s`: scalar distances of neighbours
 * `r`: relative positions of neighbours (vectorial dist)

This is a convenience function that does some of the work of constructing the
neighborhood. This is probably fairly inefficient to use since it has
to construct a `PyArray` object for the positions every time it is called.
Instead, use the iterator. Problem is, this is not much faster, because
the PyCall conversion overhead is so horrendous.
"""
function get_neighbors(n::Integer, neiglist::ASENeighborList, atm::ASEAtoms;
                       rcut=Inf)
    (inds, offsets) = ASE.neighbors(n, neiglist)
    X = _get_positions_ref_(atm)
    cell = get_cell(atm)
    r = X[inds, :]' + cell' * offsets' .- slice(X, n, :)
    s = sqrt(sumabs2(r, 1))
    I = find(s .<= neiglist.cutoffs[n])
    return inds[I], s[I], r[:, I]
end

##################### NEIGHBOURHOOD ITERATOR ###################
# this is about twice as fast as `get_neighbours`
# which indicates that, either, the ASE neighbour list is very slow
#   or, the overhead from the python call is horrendous!
#

type ASEAtomIteratorState
    at::ASEAtoms
    neiglist::ASENeighborList
    n::Int         # iteration index
    X::Matrix{Float64}
    cell_t::Matrix{Float64}
end

ASEAtomIteratorState(at::ASEAtoms, neiglist::ASENeighborList) =
    ASEAtomIteratorState( at, neiglist, 0,
                          get_positions(at),
                          get_cell(at)' )

import Base.start
start(I::Tuple{ASEAtoms,ASENeighborList}) =
    ASEAtomIteratorState(I...)

import Base.done
done(I::Tuple{ASEAtoms,ASENeighborList}, state::ASEAtomIteratorState) =
    (state.n == length(state.at))

import Base.next
function next(I::Tuple{ASEAtoms,ASENeighborList}, state::ASEAtomIteratorState)
    state.n += 1
    inds, offsets = neighbors(state.n, state.neiglist)
    r = state.X[:,inds] + state.cell_t * offsets' .- state.X[:,state.n]
    s = sqrt(sumabs2(r, 1))
    I = find(s .<= I[2].cutoffs[state.n])
    return (state.n, inds, s, r), state
    # now find the indices that are actually within the cut-off
    #  TODO: resolve the nasty issue of 1/2 cut-off first?!?!?
    # I = find(s .< 3.0)   # state.cutoffs[state.n])  (hard-coded for debugging)
    # return (state.n, inds[I], s[I], r[:,I]), state
end


text:
const BCACHE_COUNT = something(VectorizationBase.CACHE_COUNT[3], 1);
const BSIZE = Int(something(cache_size(Float64, Val(3)), 393216));

const BCACHE_LOCK = Atomic{UInt}(zero(UInt))

struct BCache{T<:Union{UInt,Nothing}}
    p::Ptr{Float64}
    i::T
end
BCache(i::Integer) = BCache(pointer(BCACHE)+8BSIZE*i, i % UInt)
BCache(::Nothing) = BCache(pointer(BCACHE), nothing)

@inline Base.pointer(b::BCache) = b.p
@inline Base.unsafe_convert(::Type{Ptr{T}}, b::BCache) where {T} = Base.unsafe_convert(Ptr{T}, b.p)

function _use_bcache()
    while true
        x = one(UInt)
        if BCACHE_COUNT > 1
            for i ∈ 0:BCACHE_COUNT-1
                if iszero(atomic_or!(BCACHE_LOCK, x) & x) # we've now set the flag, `i` was free
                    return BCache(i)
                end
                x <<= one(UInt)
            end
        else
            if iszero(atomic_or!(BCACHE_LOCK, x) & x) # we've now set the flag, `i` was free
                return BCache(nothing)
            end
        end
        pause()
    end
end
_free_bcache!(b::BCache{UInt}) = (atomic_xor!(BCACHE_LOCK, one(UInt) << b.i); nothing)
_free_bcache!(b::BCache{Nothing}) = reseet_bcache_lock!()#atomic_xor!(BCACHE_LOCK, one(UInt))

"""
  reset_bcache_lock!()

Currently not using try/finally in matmul routine, despite locking.
So if it errors for some reason, you may need to manually call `reset_bcache_lock!()`.
"""
reseet_bcache_lock!() = (BCACHE_LOCK[] = zero(UInt); nothing)

let ityp = "i$(8sizeof(UInt))"
    @eval begin
        @inline function _atomic_load(ptr::Ptr{UInt})
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %v = load atomic volatile $(ityp), $(ityp)* %p acquire, align $(Base.gc_alignment(UInt))
              ret $(ityp) %v
            """), UInt, Tuple{Ptr{UInt}}, ptr)
        end
        @inline function _atomic_store!(ptr::Ptr{UInt}, x::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              store atomic volatile $(ityp) %1, $(ityp)* %p release, align $(Base.gc_alignment(UInt))
              ret void
            """), Cvoid, Tuple{Ptr{UInt}, UInt}, ptr, x)
        end
        @inline function _atomic_cas_cmp!(ptr::Ptr{UInt}, cmp::UInt, newval::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %c = cmpxchg volatile i64* %p, i64 %1, i64 %2 acq_rel acquire
              %bit = extractvalue { i64, i1 } %c, 1
              %bool = zext i1 %bit to i8
              ret i8 %bool
            """), Bool, Tuple{Ptr{UInt}, UInt, UInt}, ptr, cmp, newval)
        end
        @inline function _atomic_add!(ptr::Ptr{UInt}, x::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %v = atomicrmw volatile add $(ityp)* %p, $(ityp) %1 acq_rel
              ret $(ityp) %v
            """), UInt, Tuple{Ptr{UInt}, UInt}, ptr, x)
        end
        @inline function _atomic_or!(ptr::Ptr{UInt}, x::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %v = atomicrmw volatile or $(ityp)* %p, $(ityp) %1 acq_rel
              ret $(ityp) %v
            """), UInt, Tuple{Ptr{UInt}, UInt}, ptr, x)
        end
        @inline function _atomic_xor!(ptr::Ptr{UInt}, x::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %v = atomicrmw volatile xor $(ityp)* %p, $(ityp) %1 acq_rel
              ret $(ityp) %v
            """), UInt, Tuple{Ptr{UInt}, UInt}, ptr, x)
        end
        @inline function _atomic_max!(ptr::Ptr{UInt}, x::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %v = atomicrmw volatile umax $(ityp)* %p, $(ityp) %1 acq_rel
              ret $(ityp) %v
            """), UInt, Tuple{Ptr{UInt}, UInt}, ptr, x)
        end
        @inline function _atomic_min!(ptr::Ptr{UInt}, x::UInt)
            Base.llvmcall($("""
              %p = inttoptr $(ityp) %0 to $(ityp)*
              %v = atomicrmw volatile umin $(ityp)* %p, $(ityp) %1 acq_rel
              ret $(ityp) %v
            """), UInt, Tuple{Ptr{UInt}, UInt}, ptr, x)
        end
    end
end


text:
@testset "1220.count-vowels-permutation.jl" begin
    @test count_vowel_permutation(1) == 5
    @test count_vowel_permutation(2) == 10
    @test count_vowel_permutation(5) == 68
    @test count_vowel_permutation(10^4) == 76428576
end


text:
using Jeff, Test, DelimitedFiles

LAYERS = Array{Any}(undef, (3, 4))
LAYERS[1, 1] = 0.0
LAYERS[1, 2] = 0.0
LAYERS[1, 3] = 0.0
LAYERS[1, 4] = 0.0
LAYERS[2, 1] = 100.
LAYERS[2, 2] = 3.47
LAYERS[2, 3] = 0.
LAYERS[2, 4] = 2.
LAYERS[3, 1] = 0.0
LAYERS[3, 2] = 2.07
LAYERS[3, 3] = 0.0
LAYERS[3, 4] = 3.

@testset "abeles" begin
    data = readdlm(string(pwd(), Base.Filesystem.path_separator, "theoretical.txt"))
    q_test = data[:, 1]
    r_test = data[:, 2]

    result = Jeff.abeles(q_test, LAYERS)

    @test all(isapprox.(result, r_test, atol=1e-9))
end

@testset "constant_smearing" begin
    data = readdlm(string(pwd(), Base.Filesystem.path_separator, "theoretical.txt"))
    q_test = data[:, 1]
    r_test = data[:, 2]

    result = Jeff.constant_smearing(q_test, LAYERS, 0.0, 1.0, 0.0)

    @test all(isapprox.(result, r_test, atol=1e-9))
end


text:
@time @testset "Homogeneous Black brane tests:" begin

    grid = SpecCartGrid3D(
        x_min            = -5.0,
        x_max            =  5.0,
        x_nodes          =  6,
        y_min            = -5.0,
        y_max            =  5.0,
        y_nodes          =  6,
        u_outer_min      =  0.1,
        u_outer_max      =  1.003,
        u_outer_domains  =  4,
        u_outer_nodes    =  32,
        u_inner_nodes    =  12,
    )

    id = BlackBrane()

    evoleq = AffineNull(
        phi0       = 0.0,
        potential  = ZeroPotential(),
        gaugecondition = AdS5_3_1.ConstantGauge(),
    )

    # atlas of grid configuration and respective SystemPartition
    atlas     = Atlas(grid)
    systems   = SystemPartition(grid)

    # allocate variables
    boundary  = Boundary(grid)
    gauge     = Gauge(grid)
    gauge_t   = similar(gauge)

    bulkevols      = BulkEvolvedPartition(grid)
    bulkconstrains = BulkConstrainedPartition(grid)
    bulkderivs     = BulkDerivPartition(grid)
    horizoncache   = AdS5_3_1.HorizonCache(systems[end], evoleq.gaugecondition.fd_order)

    # initial conditions
    id(bulkconstrains, bulkevols, bulkderivs, boundary, gauge,
       horizoncache, systems, evoleq)

    # function to solve the nested system, given the initial data
    nested = AdS5_3_1.Nested(systems, bulkconstrains, bulkderivs)

    # solve nested system for the constrained variables
    nested(bulkevols, boundary, gauge, evoleq)

    # analyze data

    Du_in  = systems[1].Du
    Du_out = systems[2].Du

    # inner grid
    bulk_in  = bulkconstrains[1]
    chart_in = atlas[1]
    uu_in    = chart_in.coords[1][:]

    # first outer grid
    bulk_out_1  = bulkconstrains[2]
    chart_out_1 = atlas[2]
    uu_out_1    = chart_out_1.coords[1][:]

    # take u-derivatives
    Du_A_in    = Du_in * bulk_in.A
    Du_A_out_1 = Du_out * bulk_out_1.A

    # Compare value of A and its u-derivative on the inner-outer domain interface

    xi0         = gauge.xi[1,1,1]
    A0_in       = bulk_in.A[:,1,1]
    A0_out      = bulk_out_1.A[:,1,1]
    Du_A0_in    = Du_A_in[:,1,1]
    Du_A0_out_1 = Du_A_out_1[:,1,1]

    Aout_in_end = 1/uu_in[end]^2 + 2*xi0/uu_in[end] + xi0^2 + uu_in[end]^2 * A0_in[end]

    Du_Aout_in_end = -2/uu_in[end]^3 - 2*xi0/uu_in[end]^2 + 2*uu_in[end] * A0_in[end] +
        uu_in[end]^2 * Du_A0_in[end]

    # in-out interface
    @test Aout_in_end    ≈ A0_out[1]
    @test Du_Aout_in_end ≈ Du_A0_out_1[1]

    # remaining outer grids
    for i in 2:grid.u_outer_domains
        Du_1    = systems[i].Du
        bulk_1  = bulkconstrains[i]
        chart_1 = atlas[i]
        uu_1    = chart_1.coords[1][:]

        Du_2    = systems[i+1].Du
        bulk_2  = bulkconstrains[i+1]
        chart_2 = atlas[i+1]
        uu_2    = chart_2.coords[1][:]

        # take u-derivatives
        Du_A_1    = Du_1 * bulk_1.A
        Du_A_2    = Du_2 * bulk_2.A

        # 1-2 interface
        @test all(bulk_1.A[end,:,:] .≈ bulk_2.A[1,:,:])
        @test all(Du_A_1[end,:,:]   .≈ Du_A_2[1,:,:])
    end

    # now check if the result matches the expected value. since the
    # configuration is homogeneous, let's compare at
    i = 1; j = 1

    # first let's check inner grid
    A_in    = bulk_in.A[:,i,j]
    A_u_in  = 2 .* uu_in .* bulk_in.A[:,i,j] .+ uu_in .* uu_in .* Du_A_in[:,i,j]

    a4 = -id.energy_dens/0.75
    A_in_exact = a4 ./ (1 .+ 2 .* xi0 .* uu_in .+ xi0^2 .* uu_in .* uu_in)

    A_u_in_exact = 2 .* uu_in .* a4 ./ (1 .+ xi0 .* uu_in).^3

    @test A_in    ≈ A_in_exact
    @test A_u_in  ≈ A_u_in_exact

    # now for the outer grids
    for i in 2:grid.u_outer_domains+1
        Du    = systems[i].Du
        bulk  = bulkconstrains[i]
        chart = atlas[i]
        uu    = chart.coords[1][:]

        # take u-derivative
        Du_A  = Du * bulk.A

        A     = bulk.A[:,i,j]
        A_u   = Du_A[:,i,j]

        A_exact  = 1 ./ uu.^2 .+ 2 .* xi0 ./ uu .+ xi0^2 .+
            uu.^2 .* a4 ./ (1 .+ 2 .* xi0 .* uu .+ xi0^2 .* uu.^2)

        A_u_exact = 2 .* (.- 1 .- uu .* xi0 + a4 .* uu.^4 ./ (1 .+ uu .* xi0).^3 ) ./ uu.^3

        @test A   ≈ A_exact
        @test A_u ≈ A_u_exact

    end

end


text:
const MOD = 20201227

function root(value)
    for i in 1:10^9
        if powermod(7, i, MOD) == value
            return i
        end
    end
end

# println(powermod(5764801, root(17807724), MOD) == 14897079) # test
println("P1: ", powermod(19774466, root(7290641), MOD))
println("P2: resort deposit")

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
type MIDIEvent <: TrackEvent
    dT::Int
    status::UInt8
    data::Array{UInt8,1}
end

function isstatusbyte(b::UInt8)
    (b & 0b10000000) == 0b10000000
end

function isdatabyte(b::UInt8)
    !isstatusbyte(b)
end

function isMIDIevent(b::UInt8)
    !ismetaevent(b) && !issysexevent(b)
end

function channelnumber(m::MIDIEvent)
    0x0F & m.status
end

function readMIDIevent(dT::Int64, f::IO, laststatus::UInt8)
    statusbyte = read(f, UInt8)
    highnybble = statusbyte & 0b11110000

    toread = 0
    if haskey(EVENTTYPETOLENGTH, highnybble)
        laststatus = statusbyte
        toread = EVENTTYPETOLENGTH[highnybble]
    else # Running status is in use
        toread = EVENTTYPETOLENGTH[laststatus & 0b11110000]
        statusbyte = laststatus
        skip(f, -1)
    end

    data = read(f, UInt8, toread)

    MIDIEvent(dT, statusbyte, data)
end

function writeevent(f::IO, event::MIDIEvent, writestatus::Bool)
    writevariablelength(f, event.dT)

    if writestatus
        write(f, event.status)
    end

    write(f, event.data)
end

function writeevent(f::IO, event::MIDIEvent)
    writeevent(f, event, true)
end


text:
using TransformUtils
using Base.Test

print("[TEST] constructors SO3, Quaternion, AngleAxis... ")
q = Quaternion(0)
ss = so3(0)
R = SO3(0)
A = AngleAxis(0)
T = SE3(0)
SE3(T.t, q)
SE3(T.t, R)
SE3(T.t, A)
println("[SUCCESS]")

print("[TEST] comparison functions for SO3, Quaternion, AngleAxis... ")
@test compare(Quaternion(0),Quaternion(0))
@test !compare(Quaternion(0),Quaternion(0,[1;0;0]))
@test compare(SO3(0),SO3(0))
IrotSO3 = SO3(0)
for i in 1:1000
  r = 0.1*randn(3)
  if norm(r) > 1e-10
    @test !compare(IrotSO3, convert(SO3,so3( r )) )
  end
end
@test compare(AngleAxis(0), AngleAxis(0))
println("[SUCCESS]")

print("[TEST] trivial case quaterion -> SO3 -> Euler -> quaternion... ")
q = Quaternion(0)
R = convert(SO3, q)
E = convert(Euler, R)
@test compare(q,convert(Quaternion, E))
# @show SO3(0) * Quaternion(0) * so3(0) * AngleAxis(0)
@test compare(SO3(0), SO3(0) * Quaternion(0) * so3(0) * AngleAxis(0))
println("[SUCCESS]")

println("[TEST] convert functions ")
warn("Need better coverage on convert function tests")
# print("[SUCSSESS]")

print("[TEST] compare SO3 and quaternion rotations... ")
q = Quaternion(0)
R = convert(SO3,q)
A = convert(AngleAxis, q)

for dAA in
    [AngleAxis(pi/4, [1;0;0]);
    so3(0.1*randn(3));
    so3(0.1*randn(3));
    so3(0.1*randn(3));
    convert(Quaternion, Euler(0.1*randn(3)));
    convert(SO3,Euler(0.1*randn(3)))  ]
  @show dAA
  q1 = q*dAA
  R1 = R*dAA
  AA1 = A*dAA
  @test compare( q1, convert(Quaternion, R1) )
  @test compare( R1, convert(SO3, q1) )
  @test compare( q1, convert(Quaternion,AA1) )
  @test compare( R1, convert(SO3,AA1) )
end
println("[SUCCESS]")

print("[TEST] basic SE3 mechanics... ")
# SE3 tests
a = SE3([0;0;0],SO3(0))
b = SE3([1;0;0],SO3(0))
@test compare(a.R,b.R)
c = a*b
@test compare(b,c)

ap = SE3(a.t, a.R*so3(0.1*randn(3)))
@test !compare(ap,a)

back = ap*b*SE3(zeros(3),transpose(ap.R))
@test compare(SO3(0),back.R)
println("[SUCCESS]")


text:
module NonNumeric

# Test interpolation with a multi-valued type

using Interpolations
using Compat

import Base: +, -, *, /

immutable MyPair{T}
    first::T
    second::T
end

# Here's the interface your type must define
(+)(p1::MyPair, p2::MyPair) = MyPair(p1.first+p2.first, p1.second+p2.second)
(-)(p1::MyPair, p2::MyPair) = MyPair(p1.first-p2.first, p1.second-p2.second)
(*)(n::Number, p::MyPair) = MyPair(n*p.first, n*p.second)
(*)(p::MyPair, n::Number) = n*p
(/)(p::MyPair, n::Number) = MyPair(p.first/n, p.second/n)
Base.one{T}(::Type{MyPair{T}}) = MyPair(one(T),one(T))
Base.zero{T}(::Type{MyPair{T}}) = MyPair(zero(T),zero(T))
Base.promote_rule{T1,T2<:Number}(::Type{MyPair{T1}}, ::Type{T2}) = MyPair{promote_type(T1,T2)}
Base.promote_op{T1,T2<:Number}(::typeof(@functorize(*)), ::Type{MyPair{T1}}, ::Type{T2}) = MyPair{promote_type(T1,T2)}
Base.promote_op{T1<:Number,T2}(::typeof(@functorize(*)), ::Type{T1}, ::Type{MyPair{T2}}) = MyPair{promote_type(T1,T2)}

# 1d
A = reinterpret(MyPair{Float64}, rand(2, 10), (10,))
itp = interpolate(A, BSpline(Constant()), OnGrid())
itp[3.2]
itp = interpolate(A, BSpline(Linear()), OnGrid())
itp[3.2]
itp = interpolate(A, BSpline(Quadratic(Flat())), OnGrid())
itp[3.2]

# 2d
A = reinterpret(MyPair{Float64}, rand(2, 10, 5), (10,5))
itp = interpolate(A, BSpline(Constant()), OnGrid())
itp[3.2,1.8]
itp = interpolate(A, BSpline(Linear()), OnGrid())
itp[3.2,1.8]
itp = interpolate(A, BSpline(Quadratic(Flat())), OnGrid())
itp[3.2,1.8]

end


text:
mutable struct MHSamplingTrackingData
    X::Vector{Vector{Float64}}
    X_fdens::Vector{Float64}
    X_fcost::Vector{Float64}
    X_viable::Vector{Bool}
    ratios_polytope_hits::Vector{Float64}
    fevals_density::Int64
    fevals_cost::Int64
    successes::Int64
end

function MHSamplingTrackingData()
    return MHSamplingTrackingData( Vector{Vector{Float64}}() , Vector{Float64}() , Vector{Float64}() , Vector{Bool}(), Vector{Float64}() , 0 , 0 , 0 )
end
function Base.hcat(a::MHSamplingTrackingData,b::MHSamplingTrackingData)
    return MHSamplingTrackingData( [a.X;b.X] , [a.X_fdens;b.X_fdens] , [a.X_fcost;b.X_fcost] , [a.X_viable;b.X_viable], [a.ratios_polytope_hits;b.ratios_polytope_hits] , a.fevals_density+b.fevals_density , a.fevals_cost+b.fevals_cost  , a.successes+b.successes )
end

"""
  mh_step( f , threshold , G::Array{Float64,2} , h::Array{Float64,1} , x::Array{Float64,1} , Sigma::Array{Float64,2} , max_tries::Int64  , tdata::MHSamplingTrackingData )

computes the mh step for the uniform distribution

"""
function mh_step( f , threshold::Float64 , G::Array{Float64,2} , h::Array{Float64,1} , x::Array{Float64,1} , Sigma::Array{Float64,2} , max_tries::Int64  , tdata::MHSamplingTrackingData )
    return mh_step( Void() , f , threshold , G::Array{Float64,2} , h::Array{Float64,1} , x::Array{Float64,1} , NaN , Sigma::Array{Float64,2} , max_tries::Int64  , tdata::MHSamplingTrackingData  )
end

"""
function mhstep( f , threshold , G::Array{Float64,2} , hG::Array{Float64,1} , x::Array{Float64,1} , x_density::Float64 , Sigma::Array{Float64,2} , max_tries::Integer64 )
[ x_next , x_next_cost , X_NonViableInPolytope , num_rejected_jumps , num_fevals , accepted ] = MHStep_01( f , threshold , G , h , x , Sigma , max_tries )

 performs a step of metropolis hastings, i.e. it generates proposal
 samples according to Sigma, and then accepts with probability 0 if
 outside the viable space, and accepts with probability (f_density(x_next) / x_f) if inside the
 viable space.

RETURNS: ::Tuple{ Bool , Array{Float64,1} , Float64  , Float64   ,  Integer64 }
                  succ ,       x    ,       fdens(x) ,  fcost(x) ,
"""
function mh_step( f_density , f , threshold , G::Array{Float64,2} , h::Array{Float64,1} , x::Array{Float64,1} , x_density::Float64 , Sigma::Array{Float64,2} , max_tries::Int64  , tdata::MHSamplingTrackingData )

    #dim = numel(x);
    factor_TriesPolytope = 100

    accepted = false
    n_tries  = 0

    # sample the whole batch.. :)
    x_c = mvnrnd(x,Sigma, factor_TriesPolytope * max_tries )

    # sort out the ones which are not in polytope..
    x_c_inPolytope = all( broadcast( - ,G*x_c,h) .< 0 , 1 )

    x_c = x_c[:,x_c_inPolytope[:]]
    #x_c = x_c[:,1:min(max_tries,size(x_c,2))]

    push!( tdata.ratios_polytope_hits , size(x_c,2)*1.0 / (factor_TriesPolytope * max_tries) )

    #num_rejected_jumps = 0max_tries - sum(x_c_inPolytope)

    num_fevals_fdensity = 0
    num_fevals_fcost    = 0

    zi = 1

    x_c_i      = NaN * ones( length(x) )
    x_c_i_cost = NaN

    # check if we sample non-unif.
    nonunif_sampling = !isa( f_density , Void )
    xfdens_next::Float64 = NaN # holds the value of the density function at x_next

    while( ~accepted && zi <= min( size(x_c,2) , max_tries ) )
        x_c_i = x_c[:,zi]
        # if nonunif. then first find next mcmc step
        if(nonunif_sampling)
            xfdens_next = f_density(x_c_i)
            num_fevals_fdensity += 1
            mcmc_step = rand() < ( xfdens_next / x_density )
            #print(" $(xfdens_next) / $(x_density) = $(xfdens_next/x_density) \n")
            if(!mcmc_step)
                zi = zi + 1
                continue
            end
        end

        num_fevals_fcost = num_fevals_fcost + 1
        x_c_i_cost   = f(x_c_i)

        accepted     = x_c_i_cost <= threshold

        zi = zi + 1

        push!(tdata.X   , x_c_i)
        push!(tdata.X_fdens , xfdens_next)
        push!(tdata.X_fcost , x_c_i_cost)
        push!(tdata.X_viable , accepted)
    end

    if(~accepted)
        #X_NonViableInPolytope = x_c(:,1:end)
        #num_rejected_jumps = zi-1
        tdata.fevals_cost     += num_fevals_fcost
        tdata.fevals_density  += num_fevals_fdensity

        x_next             = x
        x_next_cost        = NaN
        return (false,x_next,xfdens_next,x_next_cost,num_fevals_fdensity,num_fevals_fcost) # ::Tuple{ Bool , Array{Float64,1} , Float64 , Integer64 }
    else
        #X_NonViableInPolytope = x_c(:,1:zi-2);
        #num_rejected_jumps = num_rejected_jumps + zi - 1
        tdata.fevals_cost     += num_fevals_fcost
        tdata.fevals_density  += num_fevals_fdensity
        tdata.successes += 1

        x_next             = x_c_i
        x_next_cost        = x_c_i_cost
        return (true,x_next,xfdens_next,x_next_cost,num_fevals_fdensity,num_fevals_fcost) # ::Tuple{ Bool , Array{Float64,1} , Float64 ,  Integer64 }
    end
end


text:
module BinaryTreeImmutable

# Adopted from
# https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/binarytrees.html#binarytrees

using Base.Threads
using Printf

struct Node
    l::Union{Nothing, Node}
    r::Union{Nothing, Node}
end

function make(n::Int)
    return n === 0 ? Node(nothing, nothing) : Node(make(n-1), make(n-1))
end

function check(node::Node)
    return  1 + (node.l === nothing ? 0 : check(node.l) + check(node.r))
end

function binary_trees(io, n::Int)
    @printf io "stretch tree of depth %jd\t check: %jd\n" n+1 check(make(n+1))

    long_tree = make(n)
    minDepth = 4
    resultSize = div((n - minDepth), 2) + 1
    results = Vector{String}(undef, resultSize)
    Threads.@threads for depth in minDepth:2:n
        c = 0
        niter = 1 << (n - depth + minDepth)
        for _ in 1:niter
            c += check(make(depth))
        end
        index = div((depth - minDepth),2) + 1
        results[index] = @sprintf "%jd\t trees of depth %jd\t check: %jd\n" niter depth c
    end

    for i in results
        write(io, i)
    end

    @printf io "long lived tree of depth %jd\t check: %jd\n" n check(long_tree)
end

end #module

include("../../utils.jl")

using .BinaryTreeImmutable

@gctime BinaryTreeImmutable.binary_trees(devnull, 21)

max_stars_repo_name: vtjnash/GCBenchmarks

max_issues_repo_name: vtjnash/GCBenchmarks

max_forks_repo_name: vtjnash/GCBenchmarks


text:
function EvaluateFisherMatrixElement!(FisherMatrix::Fisherαβ, Cov::aₗₘCovariance,
    ∂Cℓα::∂Cℓ, ∂Cℓβ::∂Cℓ, Parα::String, Parβ::String)
    αMatrix = mymatmul(∂Cℓα.∂CℓArray,Cov.Covariance⁻¹)
    βMatrix = mymatmul(∂Cℓβ.∂CℓArray,Cov.Covariance⁻¹)
    Fisherℓ = mymatmul(αMatrix, βMatrix)
    fisherelement = SumℓAndTrace(Fisherℓ)
    FisherMatrix.FisherDict[Parα*"_"*Parβ] = fisherelement
end

function mymatmul(A::Array{Float64, 3}, B::Array{Float64, 3})
    C = zeros(size(A))
    @avx for l in 1:length(A[:,1,1])
        for i in 1:length(A[1,:,1])
            for j in 1:length(B[1,1,:])
                for k in 1:length(B[1,:,1])
                    C[l,i,j] += A[l,i,k]*B[l,k,j]
                end
            end
        end
    end
    return C
end

function SumℓAndTrace(Fisherℓ::Array{Float64, 3})
    fisherelement = 0
    @avx for ℓ in 1:size(Fisherℓ)[1]
        for i in 1:size(Fisherℓ)[2]
            fisherelement += Fisherℓ[ℓ, i, i]
        end
    end
    return fisherelement
end

function EvaluateFisherMatrixElement!(FisherMatrix::Fisherαβ, Cov::CℓCovariance,
    ∂Cℓα::∂Cℓ, ∂Cℓβ::∂Cℓ, Parα::String, Parβ::String)
    ℓnumber = length(∂Cℓα.∂CℓArray[:,1,1])
    inumber = length(∂Cℓα.∂CℓArray[1,:,1])
    L = EliminationMatrix(inumber)

    vecp∂CℓαᵀCov⁻¹ = zeros(1,floor(Int,inumber*0.5*(inumber+1)))
    fisher_temp = zeros(1,1)
    fisherelement = 0
    for ℓ in 1:ℓnumber
        vec∂Cℓα = vec(∂Cℓα.∂CℓArray[ℓ,:,:])
        vecp∂Cℓα = zeros(floor(Int,inumber*0.5*(inumber+1)))
        LinearAlgebra.mul!(vecp∂Cℓα, L, vec∂Cℓα)

        vec∂Cℓβ = vec(∂Cℓβ.∂CℓArray[ℓ,:,:])
        vecp∂Cℓβ = zeros(floor(Int,inumber*0.5*(inumber+1)))
        LinearAlgebra.mul!(vecp∂Cℓβ, L, vec∂Cℓβ)

        Covariance⁻¹ = Cov.Covariance⁻¹[ℓ,:,:]
        LinearAlgebra.mul!(vecp∂CℓαᵀCov⁻¹, transpose(vecp∂Cℓα), Covariance⁻¹)
        LinearAlgebra.mul!(fisher_temp, vecp∂CℓαᵀCov⁻¹, vecp∂Cℓβ)
        fisherelement += fisher_temp[1,1]
    end
    FisherMatrix.FisherDict[Parα*"_"*Parβ] = fisherelement
end


text:
using Pkg

cd("2022_Scandalist/Julia")
Pkg.activate(".")
Pkg.instantiate()

using Revise

include("../src/tweet_mining.jl")
include("../src/tweet_helpers.jl")
include("../src/tweet_graphing.jl")

using Chain
using CairoMakie
using AlgebraOfGraphics
using DataFrames
using DataFramesMeta
using Dates
using Graphs
using SparseArrays
using DataStructures
using MetaGraphs

Pepsi_data = "Data/pepsi.csv" |> read_tw_data

Redcross_data = "Data/RedCrossAU.csv" |> read_tw_data

mit_data = "Data/MIT2.csv" |> read_tw_data

g,userL = generate_author_retweet_source_graph(mit_data)
author_mg = generate_unipartite_projection(g, userL)

generate_small_user_author_plot(g, 5)
generate_small_user_author_plot(author_mg, 5)
generate_most_connected_user_author_plot(author_mg)

mit_data.created_at = Date.(ZonedDateTime.(mit_data.created_at))
sort(mit_data, order(:created_at))

everyone_pepsi = @chain Pepsi_data begin
    @subset(week.(:date) .== 16)
    get_all_username()
end

set_aog_theme!()

sort(Pepsi_data,:quote_count, rev = true)
Pepsi_data

Pepsi_data = add_mentions(Pepsi_data)

# Wrangling functions

get_volume_by(Pepsi_data, week)

get_cumulative_volume_by(Pepsi_data, week)

get_interactions_by(Pepsi_data, dayofyear, mean)

get_interactions_by(Pepsi_data, dayofyear, sum)

add_mentions(Pepsi_data)

get_mentions_by(Pepsi_data, week, mean)

get_mentions_by(Pepsi_data, week, sum)

# Plotting functions

generate_volume_by_plot(Pepsi_data, dayofyear)

generate_cumulative_volume_by_plot(Pepsi_data, dayofyear)

generate_avg_mentions_by_plot(Pepsi_data, dayofyear)

generate_total_mentions_by_plot(Pepsi_data, dayofyear)

generate_avg_retweets_by_plot(Pepsi_data, dayofyear)

generate_total_retweets_by_plot(Pepsi_data, dayofyear)

generate_avg_quotes_by_plot(Pepsi_data, dayofyear)

generate_total_quotes_by_plot(Pepsi_data, dayofyear)

# ROUGH MENTIONS NETWORK TESTING

g = generate_mentions_graph(Pepsi_data)

savegraph("mention_graph.lgz", g)
g = loadgraph("mention_graph.lgz")

cen_df = DataFrame(user_name=total_nodes, centrality=eigenvector_centrality(g))
print(first(sort(cen_df, :centrality, rev=true), 20))

c = counter(filter(!=("NA"), Pepsi_data.sourcetweet_author_id))
c_df = DataFrame(source_id=collect(keys(c)), count=collect(values(c)))
print(first(sort(c_df, :count, rev=true), 20))

# ROUGH REPLY / QUOTES NETWORK TESTING

g = generate_reply_quote_graph(Redcross_data)

savegraph("reply_quote_graph.lgz", g)
g = loadgraph("reply_quote_graph.lgz")

mg = add_meta_reply_quote_graph(Redcross_data, g)

slice_mg = take_slice(mg, Date(2020,2,3), Date(2020,2,10))

function graph_metric(x)
    return mean(eigenvector_centrality(x))
end

plot_graph_metric(mg, graph_metric, Week(1), Date.(ZonedDateTime.(Redcross_data.created_at))[1], Date.(ZonedDateTime.(Redcross_data.created_at))[size(Redcross_data,1)])

#User plot

user_info = groupby(select(Redcross_data, :user_username, :retweet_count, :like_count, :quote_count, :reply_count), :user_username)
user_info = combine(user_info, nrow => :ntweets, :retweet_count => sum => :nretweets, :like_count => sum => :nlikes, :quote_count => sum => :nquotes, :reply_count => sum => :nreplies)

user_info = sort(user_info, order(:nretweets, rev=false))

y = filter(x -> x > 100, user_info.nretweets)
x = collect(1:1:size(y,1))
users = data((; x, y)) *
visual(Lines) *
mapping(:x, :y)

draw(users)


text:
include("euler/euler.jl")
using .Calculus: fibonacci_index, fibonacci_numbers
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    index = fibonacci_index(n)
    fibonacci = fibonacci_numbers(index + 1, Int)
    last_sum, new_sum = 0, 0
    result = sum(fibonacci[1:4]) - 1
    for i ∈ 5:index
        last_sum, new_sum = new_sum, new_sum + last_sum + fibonacci[i - 3]
        result += new_sum + fibonacci[i]
    end
    return result
end

compute(13)

compute(23_416_728_348_467_685)

@benchmark compute(23_416_728_348_467_685)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
module TestImportedNames

using AutoImports: importednames
using CSTParser: CSTParser
using Test

inames(code) = importednames(CSTParser.parse(CSTParser.ParseState(code), true)[1].args[1])

@testset begin
    @test inames("using A") == ["A"]
    @test inames("using A: B, C") == ["B", "C"]
    @test inames("using .A: B, C") == ["B", "C"]
    @test inames("using A.B") == ["B"]
    @test inames("using A.B: C") == ["C"]
end

end  # module


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

#VERSION >= v"0.4.0-dev" && __precompile__()

module ParallelAccelerator

import CompilerTools.DebugMsg
DebugMsg.init()

importall CompilerTools
using CompilerTools.OptFramework

#######################################################################
# Constants and functions below this line are not exported, and should
# be used either fully qualified, or explicitly imported.
#######################################################################

const OFF_MODE = 0
const HOST_MODE = 1
const OFFLOAD1_MODE = 2
const OFFLOAD2_MODE = 3
const TASK_MODE = 4
const THREADS_MODE = 5

num_acc_allocs = 0
num_acc_parfors = 0

"""
(internal) Sets the number of allocations left in the accelerated function after optimizations.
"""
function set_num_acc_allocs(allocs::Int)
    global num_acc_allocs = allocs
    return
end

"""
(internal) Sets the number of parfors generated for the accelerated function after optimizations.
"""
function set_num_acc_parfors(parfors::Int)
    global num_acc_parfors = parfors
    return
end

"""
Returns the number of allocations left in the accelerated function after optimizations.
"""
function get_num_acc_allocs()
    return num_acc_allocs
end

"""
Returns the number of parfors generated for the accelerated function after optimizations.
"""
function get_num_acc_parfors()
    return num_acc_parfors
end

"""
Generate a file path to the directory above the one containing this source file.
This should be the root of the package.
"""
function getPackageRoot()
  joinpath(dirname(@__FILE__), "..")
end

const USE_ICC = 0
const USE_GCC = 1
const USE_MINGW = 2
const NONE = 3

use_bcpp = 0
backend_compiler = USE_ICC
mkl_lib = ""
openblas_lib = ""
sys_blas = 0
openmp_supported = 0
package_root = getPackageRoot()

#config file overrides backend_compiler variable
if isfile("$package_root/deps/generated/config.jl")
  include("$package_root/deps/generated/config.jl")
end

function getMklLib()
  return mkl_lib
end

function getOpenblasLib()
  return openblas_lib
end

function getSysBlas()
  return sys_blas
end

function getUseBcpp()
  return use_bcpp
end

function getBackendCompiler()
  return backend_compiler
end

cached_mode = nothing

"""
Return internal mode number by looking up environment variable "PROSPECT_MODE".
"""
function getPseMode()
    global cached_mode
    if cached_mode == nothing
        if haskey(ENV,"PROSPECT_MODE")
            env_mode = ENV["PROSPECT_MODE"]
            if (env_mode == "cgen" || env_mode == "offload1" || env_mode == "offload2" || env_mode == "task") && backend_compiler == NONE
                println("ParallelAccelerator backend CGen requested but no C compiler is installed...")
                println("...reverting to Julia native threading backend.")
                ENV["PROSPECT_MODE"] = "threads"
            end
        end

        if haskey(ENV,"PROSPECT_MODE")
           mode = ENV["PROSPECT_MODE"]
        else
           mode = "threads"
        end

        if mode == "none" || mode == "off"
          cached_mode = OFF_MODE
        elseif mode == "cgen"
          cached_mode = HOST_MODE
        elseif mode == "offload1"
          cached_mode = OFFLOAD1_MODE
        elseif mode == "offload2"
          cached_mode = OFFLOAD2_MODE
        elseif mode == "task"
          cached_mode = TASK_MODE
        elseif mode == "threads"
          cached_mode = THREADS_MODE
        else
          error(string("Unknown PROSPECT_MODE = ", mode))
        end
    end

    return cached_mode
end

const NO_TASK_MODE = 0
const HOST_TASK_MODE = 1
const PHI_TASK_MODE = 2
const DYNAMIC_TASK_MODE = 3

"""
Return internal mode number by looking up environment variable "PROSPECT_TASK_MODE".
If not specified, it defaults to NO_TASK_MODE, or DYNAMIC_TASK_MODE when
getPseMode() is TASK_MODE.
"""
function getTaskMode()
  if haskey(ENV,"PROSPECT_TASK_MODE")
     mode = ENV["PROSPECT_TASK_MODE"]
  else
    if getPseMode() == TASK_MODE
      mode = "dynamic"
    else
      mode = "none"
    end
  end
  if mode == "none" || mode == "off"
    NO_TASK_MODE
  elseif mode == "task"
    HOST_TASK_MODE
  elseif mode == "phi"
    PHI_TASK_MODE
  elseif mode == "dynamic"
    DYNAMIC_TASK_MODE
  else
    error(string("Unknown PROSPECT_TASK_MODE = ", mode))
  end
end

type UnsupportedFeature <: Exception
  text :: AbstractString
end

"""
Call this function if you want to embed binary-code of ParallelAccelerator into your Julia build to speed-up @acc compilation time.
It will attempt to add a userimg.jl file to your Julia distribution and then re-build Julia.
"""
function embed(julia_root)
  # See if the specified julia_root path exists.
  if !ispath(julia_root)
    println("The specified path to the Julia source code, ", julia_root, ", was not valid.")
    return nothing
  end

  # Do a check to see if this is likely to be a Julia distribution by looking for a base directory.
  base_dir = string(julia_root, "/base")
  if !ispath(base_dir)
    println("The specified path to the Julia source code, ", julia_root, ", does not seem to be a Julia source distribution; base dir not found.")
    return nothing
  end

  # The contents of the userimg.jl file.
  current_file = @__FILE__
  userimg_contents = string("""
# This is a dummy use of the ParallelAccelerator module.  It will
# cause a brief delay once, when Julia compiles, in exchange for *not*
# having a delay every time the package is used.

Base.reinit_stdio()

include("$current_file")

# This is to force the target of CGen's atexit hook to transitively compile so that
# compilation will not be attempted atexit time during Julia's build process.
ParallelAccelerator.CGen.CGen_finalize()

module __UserimgDummyModule__

using ParallelAccelerator
const dottimes=GlobalRef(ParallelAccelerator.API, :*)
const dotplus=GlobalRef(ParallelAccelerator.API, :-)
const runStencil=GlobalRef(ParallelAccelerator.API, :runStencil)

@eval __userimg_dummy_fn__(A::Array{Float64,2}, B::Array{Float64,2}) = begin
    \$runStencil((a, b) -> a[0,0] = b[0,0], A, B, 1, :oob_skip)
    \$dotplus(2, \$dottimes(A, B))
end

ParallelAccelerator.accelerate(__userimg_dummy_fn__, (Array{Float64,2},Array{Float64,2},))

end
""")

  userimg_file = string(base_dir, "/userimg.jl")
  # If there is already a userimg.jl file then just report the need for a manual merge.
  if isfile(userimg_file)
    println("The Julia source tree already seems to have a userimg.jl file.")
    println("A file ParallelAccelerator_userimg.jl has been created in ", base_dir, ".")
    println("Please manually merge this file with userimg.jl and then re-compile Julia.")

    outfile = string(base_dir, "/ParallelAccelerator_userimg.jl")
    f = open(outfile, "w")
    println(f, userimg_contents)
    close(f)
  else
    # There is no existing userimg.jl file so create one with the userimg_contents.
    f = open(userimg_file, "w")
    println(f, userimg_contents)
    close(f)

    curwd = pwd()    # Get the current directory so we can switch back to that later.
    cd(julia_root)   # Switch the the root of the julia distribution where we can run make.
    run(`make`)      # Re-build Julia.
    cd(curwd)        # Go back to the initial working directory.
  end

  return nothing
end

"""
This version of embed tries to use JULIA_HOME to find the root of the source distribution.
It then calls the version above specifying the path.
"""
function embed()
  embed(joinpath(JULIA_HOME, "..", ".."))
end

# a hack to make accelerate function and DomainIR mutually recursive.
_accelerate(function_name, signature) = accelerate(function_name, signature, 0)

# data types for declaring size array for data sources
# will be matched in CGen to proper C arrays
type H5SizeArr_t
end
type SizeArr_t
end

function show_backtrace()
    bt = backtrace()
    s = sprint(io->Base.show_backtrace(io, bt))
    println(s)
end

using Compat

include("api.jl")
include("domain-ir.jl")
include("parallel-ir.jl")
include("j2c-array.jl")
include("cgen.jl")
#include("callgraph.jl")
include("comprehension.jl")
include("driver.jl")

importall .Driver

"""
Called when the package is loaded to do initialization.
"""
function __init__()
    package_root = getPackageRoot()
    if Compat.is_linux()
        ld_env_key = "LD_LIBRARY_PATH"
    elseif Compat.is_apple()
        ld_env_key = "DYLD_LIBRARY_PATH"
    elseif Compat.is_windows()
        ld_env_key = "PATH"
    else
        throw(string("System is not linux, apple, or windows...giving up."))
    end
    prefix = ""
    # See if the LD_LIBRARY_PATH environment varible is already set.
    if haskey(ENV, ld_env_key)
        # Prepare the current value to have a new path added to it.
        prefix = ENV[ld_env_key] * ":"
    end
    # Add the bin directory off the package root to the LD_LIBRARY_PATH.
    ENV[ld_env_key] = string(prefix, package_root, "bin")

    # Must add LD_LIBRARY_PATH for MIC
    if getPseMode() == OFFLOAD1_MODE || getPseMode() == OFFLOAD2_MODE || getPseMode() == TASK_MODE
        ENV["LD_LIBRARY_PATH"]=string(CGen.generated_file_dir, ":", ENV["LD_LIBRARY_PATH"])
        #ENV["MIC_LD_LIBRARY_PATH"]=string(generated_file_dir, ":", ENV["MIC_LD_LIBRARY_PATH"])
    end

    if getPseMode() == OFF_MODE
      addOptPass(runStencilMacro, PASS_MACRO)
      #addOptPass(cleanupAPI, PASS_MACRO)
    else
      addOptPass(captureOperators, PASS_MACRO)
      addOptPass(toCartesianArray, PASS_MACRO)
      addOptPass(expandParMacro, PASS_MACRO)
      addOptPass(extractCallGraph, PASS_TYPED)
      addOptPass(toDomainIR, PASS_TYPED)
      addOptPass(toParallelIR, PASS_TYPED)
      addOptPass(toFlatParfors, PASS_TYPED)
      if getPseMode() == THREADS_MODE
        addOptPass(toJulia, PASS_TYPED)
      else
        addOptPass(toCGen, PASS_TYPED)
      end
    end
end

import .API.@par
import .API.runStencil
import .API.cartesianarray
import .API.parallel_for
export accelerate, @acc, @noacc, @par, runStencil, cartesianarray, parallel_for

end


text:
module PyRhodium

using PyCall
using PyPlot
using IterableTables
using Distributions
using DataFrames
using DataStructures

export Model, Parameter, Response, Lever, RealLever, IntegerLever, CategoricalLever,
    PermutationLever, SubsetLever, Constraint, Brush, DataSet, pandas_dataframe,
    named_tuple, named_tuples, optimize, scatter2d, scatter3d, pairs,
    parallel_coordinates, apply, evaluate, sample_lhs, set_parameters!,
    set_levers!, set_responses!, set_constraints!, set_uncertainties!,
    Prim, PrimBox, find_box, find_all, show_tradeoff, stats, limits,
    Cart, show_tree, print_tree, save, save_pdf, save_png,
    Sensitivity, SAResult, sa, oat, plot, plot_sobol, find

const rhodium = PyNULL()
const prim = PyNULL()
const pd = PyNULL()
const seaborn = PyNULL()

function __init__()
    copy!(rhodium, pyimport("rhodium"))
    copy!(prim, pyimport("prim"))
    copy!(pd, pyimport("pandas"))
    copy!(seaborn, pyimport("seaborn"))

    # TBD: see if it works to simply call __init__(function) without storing the julia function
    py"""
    from rhodium import *
    class JuliaModel(Model):

        def __init__(self, function, **kwargs):
            super(JuliaModel, self).__init__(self._evaluate)
            self.j_function = function

        def _evaluate(self, **kwargs):
            result = self.j_function(**kwargs)
            return result
    """

    # Create a Python class that can store the SAResult (a subclass of dict) in
    # an instance var so we can access it without conversion. Otherwise, PyCall
    # converts it to a Dict and we can't use it as an argument to the plot routines
    # that are methods of rhodium.SARsult.
    py"""
    import rhodium

    class SAResultContainer(object):
        def __init__(self, sa_result):
            self.sa_result = sa_result

    def my_sa(*args, **kwargs):
        sa_result = rhodium.sa(*args, **kwargs)
        return SAResultContainer(sa_result)
    """
end

include("core.jl")
include("prim.jl")
include("cart.jl")
include("sa.jl")

end # module


text:
#=

Purpose: Use the Schur method to find the bounded solution of the
         discrete Lyapunov equation:

         X = B + A*X*A'.

Inputs:  A (n*n)
         B (n*n) = Symmetric positive semidefinite

Outputs: X (n*n) = Symmetric positive semidefinite solution.

Based on Octave code written by A.S. Hodel.

=#

function dlyap{T<:AbstractFloat}(a::Array{T,2},b::Array{T,2})

  n = size(a,1)
  x = zeros(n,n)
  j = n

  (s,u) = schur(a)
  b = u'b*u

  while j > 0
    j1 = j
    if j == 1
      block = 1
    elseif s[j,j-1] > 2*eps(T)*norm(s)
      block = 2
      j -= 1
    else
      block = 1
    end
    lhs = kron(s[j:j1,j:j1],s) - I # I = eye(block*n)
    rhs = vec(b[:,j:j1])
    if j1 < n
      rhs2 = s*(x[:,(j1+1):n]*s[j:j1,(j1+1):n]')
      rhs += vec(rhs2)
    end
    w = -lhs\rhs
    x[:,j] = w[1:n]
    if block == 2
      x[:,j1] = w[(n+1):block*n]
    end
    j -= 1
  end

  x = u*x*u'

  return x

end


text:
####################################################################################
####################################################################################
# Trapezoidal rule
####################################################################################
####################################################################################

# Regular grid for 1D array
trapz(y::AbstractVector,dx::Real=1)  = trapzRegular1D(
									OffsetArrays.no_offset_view(y),dx)

# Regular grid for nD array
function trapz(y::AbstractArray,dx::Real=1;axis=1)
	T = promote_type(eltype(y), Float64)

	if axis != 1
		perm = [axis]
		append!(perm,[i for i in 1:axis-1])
		append!(perm,[i for i in axis+1:ndims(y)])
		y = permutedims(y,perm)
	end

	trapzRegular1D(OffsetArrays.no_offset_view(y),dx)
end

# Irregular grid for 1D array
trapz(y::AbstractVector,x::AbstractVector{<:Real}) = trapzIrregular1D(
												OffsetArrays.no_offset_view(y),
												OffsetArrays.no_offset_view(x))

# Irregular grid for nD array
function trapz(y::AbstractArray,x::AbstractVector{<:Real};axis=1)
	T = promote_type(eltype(y), Float64)

	if axis != 1
		perm = [axis]
		append!(perm,[i for i in 1:axis-1])
		append!(perm,[i for i in axis+1:ndims(y)])
		y = permutedims(y,perm)
	end

	trapzIrregular1D(OffsetArrays.no_offset_view(y),OffsetArrays.no_offset_view(x))
end

####################################################################################

function trapzweight(Nelems,index)
	@assert(1<=index<=Nelems,"index not in range")
	if index==1 || index==Nelems
		return 1
	end
	return 2
end

######################################################################################

# Regular grid for 1D array
function trapzRegular1D(y::AbstractVector,dx::Real=1)
	N = size(y,1)
	T = promote_type(eltype(y),Float64)
	int_y = zero(T)
	for (ind,y_i) in enumerate(y)
		w = trapzweight(N,ind)
		int_y += w*y_i
	end
	return int_y * dx/2
end

function trapzRegular1D!(y::AbstractVector,
	int_y::AbstractArray{T},dx::Real=1) where {T<:Number}

	N = size(y,1)

	fill!(int_y,zero(T))

	for (ind,y_i) in enumerate(y)
		w = trapzweight(N,ind)
		@. int_y += w * y_i * dx/2
	end

	return int_y
end

# UnitRanges can be integrated exactly
function trapzRegular1D(y::AbstractUnitRange,dx::Real=1)
	(last(y)^2 - first(y)^2) * dx/2
end

# N dimensional array, integrate along the first axis
function trapzRegular1D(y::AbstractArray,dx::Real=1)

	trailing_axes,inds = trailing_indices(y)

	T = promote_type(eltype(y),Float64)
	int_y = zeros(T,inds)

	for ind_t in trailing_axes
		int_y[ind_t] = trapzRegular1D(view(y,:,ind_t),dx)
	end
	return int_y
end

function trapzRegular1D!(y::AbstractArray,
	int_y::AbstractArray{<:Number},dx::Real=1)

	trailing_axes,_ = trailing_indices(y)

	for ind_t in trailing_axes
		int_y[ind_t] = trapzRegular1D(view(y,:,ind_t),dx)
	end

	return int_y
end

function trapzRegular1DAdd!(y::AbstractArray,
	int_y::AbstractArray{<:Number},dx::Real=1)

	trailing_axes,_ = trailing_indices(y)

	for ind_t in trailing_axes
		int_y[ind_t] += trapzRegular1D(view(y,:,ind_t),dx)
	end

	return int_y
end

####################################################################################

# Irregular grid for 1D array
function trapzIrregular1D(y::AbstractVector,x::AbstractVector{<:Real})
	first_ind = first(axes(y,1))
	last_ind = last(axes(y,1))

	first_ind_x = first(axes(x,1))
	last_ind_x = last(axes(x,1))

	T = promote_type(eltype(y),Float64)
	int_y = zero(T)

	for (ind_x,ind_y) in zip(first_ind_x:last_ind_x-1,first_ind:last_ind-1)
		dx = x[ind_x+1] - x[ind_x]
		int_y += (y[ind_y]+y[ind_y+1])/2 * dx
	end
	return int_y
end

function trapzIrregular1D!(y::AbstractVector,
	int_y::AbstractArray{T},x::AbstractVector{<:Real}) where {T<:Number}

	first_ind = first(axes(y,1))
	last_ind = last(axes(y,1))

	first_ind_x = first(axes(x,1))
	last_ind_x = last(axes(x,1))

	fill!(int_y,zero(T))

	for (ind_x,ind_y) in zip(first_ind_x:last_ind_x-1,first_ind:last_ind-1)
		dx = x[ind_x+1] - x[ind_x]
		@. int_y += (y[ind_y]+y[ind_y+1])/2 * dx
	end
	return first(int_y)
end

# N dimensional array, integrate along the first axis
function trapzIrregular1D(y::AbstractArray,x::AbstractVector{<:Real})

	trailing_axes,inds = trailing_indices(y)

	T = promote_type(eltype(y),Float64)
	int_y = zeros(T,inds)

	for ind_t in trailing_axes
		int_y[ind_t] = trapzIrregular1D(view(y,:,ind_t),x)
	end
	return int_y
end

function trapzIrregular1D!(y::AbstractArray,
	int_y::AbstractArray{<:Number},x::AbstractVector{<:Real})

	trailing_axes,_ = trailing_indices(y,axis)

	for ind_t in trailing_axes
		int_y[ind_t] = trapzIrregular1D(view(y,:,ind_t),x)
	end
	return int_y
end

function trapzIrregular1DAdd!(y::AbstractArray,
	int_y::AbstractArray{<:Number},x::AbstractVector{<:Real})

	trailing_axes,_ = trailing_indices(y,axis)

	for ind_t in trailing_axes
		int_y[ind_t] += trapzIrregular1D(view(y,:,ind_t),x)
	end

	return int_y
end

####################################################################################
# Error in trapz
#####################################################################################

function trapzerror(y::AbstractVector,dx::Number)::Number
	y_2nddiff = diff(diff(y))
	x_range = (length(y)-1)*dx
	return max(eps(),abs(x_range/12*mean(y_2nddiff)))
end

#####################################################################################
# trapz on a sphere
#####################################################################################

function trapzCircle(y::AbstractVector,
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,1));
	ϕ2π::Bool=true):: Array{<:Number,0}

	if ϕ2π
		z = trapz(y,x=ϕ)
	else

		info("Assuming that 2π is not included,"*
		" will append it to the array before integrating.\n"*
		"If 2π is included set ϕ2π=true in the function call")

		check_last_value(ϕ,2π)

		push!(y,y[1])

		z = trapz(y,x=push!(ϕ,2π))

	end

	return z
end

function trapzSpherical2D(y::AbstractArray{<:Number,2},
	θ::AbstractArray{<:Real,1}=linspace(0,π,size(y,1)),
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,2));
	θint="clenshaw_quadrature",GLweights = nothing,
	ϕ2π=true)

	if θint=="clenshaw_quadrature"

		z = clenshaw_curtis_quadrature(y) :: Vector{<:Number}

	elseif θint=="gauss_quadrature"

		if GLweights != nothing
			z = gauss_quadrature(y,GLweights) :: Vector{<:Number}
		end

	else

		z = trapz(sin.(θ).*y,x=θ) :: Vector{<:Number}

	end

	trapzCircle(z,ϕ,ϕ2π=ϕ2π) :: Array{<:Number,0}
end

function trapzSpherical3D(y::AbstractArray{<:Number,3},
	r::Union{Real,AbstractArray{<:Real,1}}=1,
	θ::AbstractArray{<:Real,1}=linspace(0,π,size(y,2)),
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,3));
	θint="clenshaw_quadrature",GLweights = nothing,
	ϕ2π=true):: Array{<:Number,0}

	if typeof(r) <: AbstractArray
		z = trapz(r.^2 .* y,x=r)
	else
		dr =r
		z = dr^2 .* trapz((1:size(y,1)).^2 .* y,dx=dr)
	end

	trapzSpherical2D(z,θ,ϕ,GLweights=GLweights,ϕ2π=ϕ2π,θint=θint)
end


text:
using XCB
using Test
using UnPack

function on_button_pressed(details::EventDetails)
    x, y = details.location
    click = details.data.button
    state = details.data.state
    buttons_pressed = pressed_buttons(state)
    printed_state = isempty(buttons_pressed) ? "" : "with $(join(string.(buttons_pressed), ", ")) button$(length(buttons_pressed) > 1 ? "s" : "") held"
    @info "$click at $x, $y $printed_state"
end

function on_key_pressed(wm::XWindowManager, details::EventDetails)
    @unpack win, data = details
    send = XCB.send(wm, win)
    km = wm.keymap
    @info keystroke_info(km, details)
    @unpack key_name, key, input, modifiers = data
    kc = KeyCombination(key, modifiers)
    set_title(win, "Random title $(rand())")
    if kc ∈ [key"q", key"ctrl+q", key"f4"]
        throw(CloseWindow(win, "Received closing request from user input"))
    elseif kc == key"s"
        curr_extent = XCB.extent(win)
        XCB.set_extent(win, curr_extent .+ 1)
    elseif kc == key"i"
        open("keymap.txt", "w") do io
            write(io, String(km))
        end
    elseif kc == key"f"
        @info "Faking input: sending key AD01 to quit (requires an english keyboard layout to be translated to the relevant symbol 'q')"
        send(key_event_from_name(wm.keymap, :AD01, KeyModifierState(), KeyPressed()))
    else
        gc = win.gc
        set_attributes(gc, [XCB.XCB_GC_FOREGROUND], [rand(1:16_777_215)])
        XCB.@flush XCB.xcb_poly_fill_rectangle(win.conn, win.id, gc.id, UInt32(1), r)
    end
end

r = Ref(XCB.xcb_rectangle_t(20, 20, 60, 60))

const is_xvfb = ENV["DISPLAY"] == ":99"

function test()
    conn = Connection()
    screen = current_screen(conn)
    win = XCBWindow(conn, screen; x=0, y=1000, border_width=50, window_title="XCB window", icon_title="XCB", attributes=[XCB.XCB_CW_BACK_PIXEL], values=[screen.black_pixel])
    ctx = GraphicsContext(win, attributes=[XCB.XCB_GC_FOREGROUND, XCB.XCB_GC_GRAPHICS_EXPOSURES], values=[screen.black_pixel, 0])
    attach_graphics_context!(win, ctx)
    wm = XWindowManager(conn, [win])

    set_callbacks!(wm, win, WindowCallbacks(;
        on_resize = x -> @info("Window size changed: $(x.data.new_dimensions)"),
        on_mouse_button_pressed = on_button_pressed,
        on_mouse_button_released = x -> @info("Released mouse button $(x.data.button)"),
        on_key_pressed = x -> on_key_pressed(wm, x),
        on_key_released = x -> @info("Released $(KeyCombination(x.data.key, x.data.modifiers))"),
        on_pointer_enter = x -> @info("Entering window at $(x.location)"),
        on_pointer_leave = x -> @info("Leaving window at $(x.location)"),
        on_pointer_move = x -> @info("Moving pointer at $(x.location)"),
        on_expose = x -> @info("Window exposed")
    ))

    send = XCB.send(wm, win)

    if is_xvfb
        @info "- Running window asynchronously"
        task = run(wm, Asynchronous(); warn_unknown=true)
        @info "- Sending fake inputs"
        send(MouseEvent(ButtonLeft(), MouseState(), ButtonPressed()))
        send(MouseEvent(ButtonLeft(), MouseState(), ButtonReleased()))
        send(PointerEntersWindowEvent())
        send(PointerMovesEvent(MouseState(), KeyModifierState()))
        send(PointerLeavesWindowEvent())
        send(key_event_from_name(wm.keymap, :AC04, KeyModifierState(), KeyReleased()))
        send(key_event_from_name(wm.keymap, :AC04, KeyModifierState(), KeyPressed()))
        @info "- Waiting for window to close"
        wait(task)
    else
        run(wm, Synchronous(); warn_unknown=true, poll=false)
    end
end

test()


text:
using Spectra
using Base.Test

# write your own tests here
println("\nMatrix size: 100")
srand(123)
n = 100
M = randn(n, n)
M = M + M'
A = DenseMatProd(M)
v0 = rand(n) - 0.5

r1 =    eigs(M; nev = 5, ncv = 12, tol = 1e-10, v0 = v0)
r2 = symeigs(A; nev = 5, ncv = 12, tol = 1e-10, v0 = v0)

@test maximum(abs(r1[1] - r2[1])) < 1e-10

println("Error: ", maximum(abs(r1[1] - r2[1])))

t1 = @elapsed(for i = 1:10    eigs(M; nev = 5, ncv = 12, tol = 1e-10, v0 = v0) end)
t2 = @elapsed(for i = 1:10 symeigs(A; nev = 5, ncv = 12, tol = 1e-10, v0 = v0) end)

println("Time of eigs(): $(t1)")
println("Time of symeigs(): $(t2)")

println("\nMatrix size: 1000")
srand(123)
n = 1000
M = randn(n, n)
M = M + M'
A = DenseMatProd(M)
v0 = rand(n) - 0.5

r1 =    eigs(M; nev = 10, ncv = 30, tol = 1e-10, v0 = v0)
r2 = symeigs(A; nev = 10, ncv = 30, tol = 1e-10, v0 = v0)

@test maximum(abs(r1[1] - r2[1])) < 1e-10

println("Error: ", maximum(abs(r1[1] - r2[1])))

t1 = @elapsed(for i = 1:10    eigs(M; nev = 10, ncv = 30, tol = 1e-10, v0 = v0) end)
t2 = @elapsed(for i = 1:10 symeigs(A; nev = 10, ncv = 30, tol = 1e-10, v0 = v0) end)

println("Time of eigs(): $(t1)")
println("Time of symeigs(): $(t2)")


text:
include("NumField/Elem.jl")
include("NumField/Field.jl")


text:
struct Taylens{N,T<:Real} <: ImplicitOp{Basis,Spin,Pix}
    # pixel remapping
    i::Matrix{Int}
    j::Matrix{Int}

    # residual displacement
    rx::Matrix{T}
    ry::Matrix{T}

    # precomputed quantities
    kα::Dict{Any,Matrix{Complex{T}}}
    xα::Dict{Any,Matrix{T}}

end

function Taylens(ϕ::FlatS0{P,T}, N) where {P,T}

    @unpack kx,ky,Nx,Ny,Δx = fieldinfo(ϕ)

    # total displacement
    d = ∇*ϕ
    dx, dy = d[1][:Ix], d[2][:Ix]

    # nearest pixel displacement
    indexwrap(ind::Int64, uplim)  = mod(ind - 1, uplim) + 1
    di, dj = (round.(Int,d/Δx) for d=(dx,dy))
    i = indexwrap.(di .+ (1:Nx)', Nx)
    j = indexwrap.(dj .+ (1:Ny) , Ny)

    # residual displacement
    rx, ry = ((d - i.*Δx) for (d,i)=[(dx,di),(dy,dj)])

    # precomputation
    kα = Dict{Any,Matrix{Complex{T}}}()
    xα = Dict{Any,Matrix{T}}()
    for n in 1:N, α₁ in 0:n
        kα[n,α₁] = im ^ n .* kx' .^ α₁ .* ky[1:Ny÷2+1] .^ (n - α₁)
        xα[n,α₁] = rx .^ α₁ .* ry .^ (n - α₁) ./ factorial(α₁) ./ factorial(n - α₁)
    end

    Taylens{N,T}(i,j,rx,ry,kα,xα)
end

# our implementation of Taylens
function *(L::Taylens{N}, f::FlatS0{P,T}) where {N,P,T}

    intlense(fx) = getindex.(Ref(fx), L.j, L.i)
    fl = f[:Il]

    # lens to the nearest whole pixel
    Lfx = intlense(f[:Ix])

    # add in Taylor series correction
    for n in 1:N, α₁ in 0:n
        Lfx .+= L.xα[n,α₁] .* intlense(fieldinfo(f).FFT \ (L.kα[n,α₁] .* fl))
    end

    FlatMap{P}(Lfx)
end

function *(L::Taylens, f::FieldTuple)
    Łf = Ł(f)
    F = typeof(Łf)
    F(map(f->L*f, Łf.fs))
end


text:
using DataValues
using Test

mutable struct TestNType{T}
    v::DataValue{T}
end

@enum TestEnum enum_val_a enum_val_b

@testset "DataValues" begin

include("scalar/test_core.jl")
include("scalar/test_operations.jl")
include("scalar/test_basederived.jl")
include("scalar/test_broadcast.jl")

include("array/test_broadcast.jl")
include("array/test_constructors.jl")
include("array/test_indexing.jl")
# include("array/test_map.jl")
include("array/test_datavaluematrix.jl")
include("array/test_datavaluevector.jl")
include("array/test_primitives.jl")
include("array/test_reduce.jl")
include("array/test_show.jl")
# include("array/test_subarray.jl")
include("array/test_typedefs.jl")

end