text:
# Make modification to the content clear the cache
function Base.setproperty!(o::Org, name::Symbol, value)
    if name == :contents
        o.cache = OrgCache(o)
        setfield!(o, name, value)
    elseif name == :settings
        setfield!(o, name, value)
    else # Anything other than .contents or .settings is a cache reference
        setfield!(getfield(o, :cache), name, value)
    end
end

function Base.getproperty(o::Org, name::Symbol)
    if name == :contents || name == :settings || name == :cache
        getfield(o, name)
    else
        getproperty(getfield(o, :cache), name)
    end
end

function Base.getproperty(c::OrgCache, name::Symbol)
    if isnothing(getfield(c, name))
        setfield!(c, name, gencache(c, Val(name)))
    end
    getfield(c, name)
end

function gencache(c::OrgCache, ::Val{:elements})
    Vector{Union{Heading, Section, GreaterElement, Element}}(
        OrgElementIterator(c.doc) |> collect)
end

function gencache(c::OrgCache, ::Val{:components})
    Vector{OrgComponent}(OrgIterator(c.doc) |> collect)
end

function gencache(c::OrgCache, ::Val{:headings})
    Vector{Heading}(filter(h -> h isa Heading, c.elements))
end

function gencache(c::OrgCache, ::Val{:keywords})
    kwdict = Dict{AbstractString, Vector{AbstractString}}()
    for kw in filter(k -> k isa Keyword, c.elements)
        if haskey(kwdict, kw.key)
            push!(kwdict[kw.key], kw.value)
        else
            kwdict[kw.key] = [kw.value]
        end
    end
    kwdict
end

function gencache(c::OrgCache, ::Val{:macros})
    macros = Dict{AbstractString, Function}()
    if haskey(c.keywords, "macro")
        for (name, replacement) in split.(c.keywords["macro"], r"[ \t]+", limit=2)
            if startswith(replacement, "(eval ")
                macros[name] = _ -> nothing
            else
                n = getproperty.(eachmatch(r"\$\d+", replacement), :match) |> unique |> length
                macros[name] = arguments -> if length(arguments) == n
                    replace(replacement, r"\$\d+" => i -> arguments[parse(Int, i[2:end])])
                end
            end
        end
    end
    macros
end

function gencache(c::OrgCache, ::Val{:footnotes})
    footnotes = Dict{Union{AbstractString, FootnoteReference}, Tuple{Int, Union{FootnoteReference, FootnoteDefinition}}}()
    i = 1
    for f in filter(f -> f isa FootnoteDefinition || f isa FootnoteReference, c.components)
        if !isnothing(f.definition)
            footnotes[something(f.label, f)] = (i, f)
            i += 1
        end
    end
    footnotes
end


text:
module TestUtilities

# using Revise
using Test
using MLJBase

end # module
true

ext: jl


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
$(TYPEDEF)

Real symmetric or complex Hermitian positive semidefinite cone of dimension
`dim` in svec format.

    $(FUNCTIONNAME){T, R}(dim::Int)
"""
mutable struct PosSemidefTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    dim::Int
    side::Int
    is_complex::Bool
    rt2::T

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}

    mat::Matrix{R}
    mat2::Matrix{R}
    mat3::Matrix{R}
    mat4::Matrix{R}
    inv_mat::Matrix{R}
    fact_mat::Cholesky{R}

    function PosSemidefTri{T, R}(
        dim::Int,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert dim >= 1
        cone = new{T, R}()
        cone.dim = dim
        cone.rt2 = sqrt(T(2))
        cone.is_complex = (R <: Complex)
        cone.side = svec_side(R, dim)
        return cone
    end
end

use_dual_barrier(::PosSemidefTri) = false

reset_data(cone::PosSemidefTri) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = false)

use_sqrt_hess_oracles(::Int, cone::PosSemidefTri) = true

function setup_extra_data!(
    cone::PosSemidefTri{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    cone.mat = zeros(R, cone.side, cone.side)
    cone.mat2 = zero(cone.mat)
    cone.mat3 = zero(cone.mat)
    cone.mat4 = zero(cone.mat)
    cone.inv_mat = zero(cone.mat)
    return cone
end

get_nu(cone::PosSemidefTri) = cone.side

function set_initial_point!(arr::AbstractVector, cone::PosSemidefTri)
    incr = (cone.is_complex ? 2 : 1)
    arr .= 0
    k = 1
    @inbounds for i in 1:cone.side
        arr[k] = 1
        k += incr * i + 1
    end
    return arr
end

function update_feas(cone::PosSemidefTri)
    @assert !cone.feas_updated

    svec_to_smat!(cone.mat, cone.point, cone.rt2)
    copyto!(cone.mat2, cone.mat)
    cone.fact_mat = cholesky!(Hermitian(cone.mat2, :U), check = false)
    cone.is_feas = isposdef(cone.fact_mat)

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::PosSemidefTri)
    svec_to_smat!(cone.mat3, cone.dual_point, cone.rt2)
    return isposdef(cholesky!(Hermitian(cone.mat3, :U), check = false))
end

function update_grad(cone::PosSemidefTri)
    @assert cone.is_feas

    inv_fact!(cone.inv_mat, cone.fact_mat)
    smat_to_svec!(cone.grad, cone.inv_mat, cone.rt2)
    cone.grad .*= -1
    copytri!(cone.mat, 'U', true)

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::PosSemidefTri)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    copytri!(cone.inv_mat, 'U', true)
    symm_kron!(cone.hess.data, cone.inv_mat, cone.rt2)
    cone.hess_updated = true
    return cone.hess
end

function update_inv_hess(cone::PosSemidefTri)
    @assert is_feas(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    symm_kron!(cone.inv_hess.data, cone.mat, cone.rt2)
    cone.inv_hess_updated = true
    return cone.inv_hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::PosSemidefTri,
    )
    @assert is_feas(cone)

    @inbounds for i in 1:size(arr, 2)
        svec_to_smat!(cone.mat4, view(arr, :, i), cone.rt2)
        copytri!(cone.mat4, 'U', true)
        rdiv!(cone.mat4, cone.fact_mat)
        ldiv!(cone.fact_mat, cone.mat4)
        smat_to_svec!(view(prod, :, i), cone.mat4, cone.rt2)
    end

    return prod
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::PosSemidefTri,
    )
    @assert is_feas(cone)

    @inbounds for i in 1:size(arr, 2)
        svec_to_smat!(cone.mat4, view(arr, :, i), cone.rt2)
        mul!(cone.mat3, Hermitian(cone.mat4, :U), cone.mat)
        mul!(cone.mat4, Hermitian(cone.mat, :U), cone.mat3)
        smat_to_svec!(view(prod, :, i), cone.mat4, cone.rt2)
    end

    return prod
end

function sqrt_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::PosSemidefTri,
    )
    @assert is_feas(cone)

    @inbounds for i in 1:size(arr, 2)
        svec_to_smat!(cone.mat4, view(arr, :, i), cone.rt2)
        copytri!(cone.mat4, 'U', true)
        rdiv!(cone.mat4, cone.fact_mat.U)
        ldiv!(cone.fact_mat.U', cone.mat4)
        smat_to_svec!(view(prod, :, i), cone.mat4, cone.rt2)
    end

    return prod
end

function inv_sqrt_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::PosSemidefTri,
    )
    @assert is_feas(cone)

    @inbounds for i in 1:size(arr, 2)
        svec_to_smat!(cone.mat4, view(arr, :, i), cone.rt2)
        copytri!(cone.mat4, 'U', true)
        rmul!(cone.mat4, cone.fact_mat.U')
        lmul!(cone.fact_mat.U, cone.mat4)
        smat_to_svec!(view(prod, :, i), cone.mat4, cone.rt2)
    end

    return prod
end

function dder3(cone::PosSemidefTri, dir::AbstractVector)
    @assert cone.grad_updated

    S = copytri!(svec_to_smat!(cone.mat4, dir, cone.rt2), 'U', true)
    ldiv!(cone.fact_mat, S)
    rdiv!(S, cone.fact_mat.U)
    mul!(cone.mat3, S, S') # TODO use outer prod function
    smat_to_svec!(cone.dder3, cone.mat3, cone.rt2)

    return cone.dder3
end


text:
# This file is a part of EncodedArrays.jl, licensed under the MIT License (MIT).

using EncodedArrays
using Test

using BitOperations

@testset "varlen_diff_codec" begin
    a = zigzagdec(UInt64(0x4c6b5f94758))
    data = [a, a + 3, a - 2]
    codec = VarlenDiffArrayCodec()

    encoded = Vector{UInt8}()
    EncodedArrays.encode_data!(encoded, codec, data)
    @test encoded == [0xd8, 0x8e, 0xe5, 0xaf, 0xeb, 0x98, 0x01, 0x06, 0x09]

    data_dec = view(similar(data), :) # Use view to test that decoder doesn't resize
    @test EncodedArrays.decode_data!(data_dec, codec, encoded) === data_dec
    @test data_dec == data

    data_dec = Vector{Int64}()
    @test EncodedArrays.decode_data!(data_dec, codec, encoded) === data_dec
    @test data_dec == data

    data_dec = Vector{Int64}(undef, 2)
    @test EncodedArrays.decode_data!(data_dec, codec, encoded) === data_dec
    @test data_dec == data

    data_dec = Vector{UInt64}()
    @test EncodedArrays.decode_data!(data_dec, codec, encoded) === data_dec
    @test data_dec == data

    b = unsigned(a)
    data = [b, b + 3, b - 2]

    encoded = Vector{UInt8}()
    EncodedArrays.encode_data!(encoded, codec, data)
    @test encoded == [0xd8, 0x8e, 0xe5, 0xaf, 0xeb, 0x98, 0x01, 0x06, 0x09]

    data_dec = view(similar(data), :) # Use view to test that decoder doesn't resize
    @test EncodedArrays.decode_data!(data_dec, codec, encoded) === data_dec
    @test data_dec == data

    data_dec = Vector{Int64}()
    @test EncodedArrays.decode_data!(data_dec, codec, encoded) === data_dec
    @test data_dec == data
end # testset


text:
using NonlinearEigenproblemsTest
using BenchmarkTools
using Statistics
using Printf
using HTTP
using Dates
using TimeZones
using JSON

include("pretty_formatting.jl")

import Base.+
import Base.zero

+(t1::BenchmarkTools.TrialEstimate, t2::BenchmarkTools.TrialEstimate) =
    BenchmarkTools.TrialEstimate(
        t1.params,
        t1.time + t2.time,
        t1.gctime + t2.gctime,
        t1.memory + t2.memory,
        t1.allocs + t2.allocs)

zero(t::BenchmarkTools.TrialEstimate) = BenchmarkTools.TrialEstimate(t.params, 0, 0, 0, 0)

const SEPARATOR = "_~/" # this string should not occur in benchmark names

recursive_mapreduce(f, op, rec, root) =
    mapreduce(x -> recursive_mapreduce(f, op, rec, x), op, rec(root); init = f(root))

mutable struct AggregatedTrial
    name
    key
    trial
end

"Hierarchically aggregate benchmark trials"
function aggregate_trials(aggregated_trials, test, chain)
    chain = [chain; test]
    if haskey(NonlinearEigenproblemsTest.benchmark_results, test)
        trial = minimum(NonlinearEigenproblemsTest.benchmark_results[test])
        for i in 1:length(chain)
            t = chain[i]
            name = " "^((i-1)*2) * t.description
            key = i == 1 ? "" : join([x.description for x in chain[2:i]], SEPARATOR)
            aggr = get!(aggregated_trials, t, AggregatedTrial(name, key, zero(trial)))
            aggr.trial += trial
        end
    end
    foreach(t -> aggregate_trials(aggregated_trials, t, chain), test.results)
    return aggregated_trials
end

"Save collected benchmark results to a JSON file on disk"
function save_benchmark(test_results, runtime_seconds, file_name)
    aggregated_trials = aggregate_trials(Dict(), test_results, [])
    benchmark = BenchmarkGroup()
    foreach(t -> benchmark[t.key] = t.trial, values(aggregated_trials))

    dict = Dict()

    dict["time"] = Dates.now(Dates.UTC)
    dict["runtime"] = runtime_seconds
    dict["config"] = @sprintf("seconds=%.1f", NonlinearEigenproblemsTest.benchmark_duration_seconds)
    dict["julia-version"] = string(VERSION)
    dict["julia-commit"] = Base.GIT_VERSION_INFO.commit_short
    dict["cpu"] = @sprintf("%s (%s), %d threads, %.0f GB memory",
        Sys.cpu_info()[1].model,
        Sys.CPU_NAME,
        Sys.CPU_THREADS,
        Sys.total_memory() / 1024^3)

    git_commit = safe_read_command(`git log -1 --format="%h"`)
    git_branch = safe_read_command(`git rev-parse --abbrev-ref HEAD`)
    git_branch == "HEAD" && (git_branch = safe_read_command(pipeline(pipeline(`git describe --all --exact-match`, stderr=devnull), `sed 's=.*/=='`)))
    isempty(git_branch) && (git_branch = "N/A")
    git_details = safe_read_command(`git log -1 --format="%ai, %s"`)
    dict["git"] = @sprintf("%s (%s), %s", git_commit, git_branch, git_details)

    io = IOBuffer()
    BenchmarkTools.save(io, benchmark)
    dict["benchmark"] = JSONText(String(io.data))

    open(file_name, "w") do f
        JSON.print(f, dict)
    end

    println("Benchmark saved to file $(abspath(file_name))")
end

"Load benchmark results from a JSON file on disk"
load_benchmark_from_file(file_name) = load_benchmark(JSON.parsefile(file_name))

"Load benchmark results from a JSON string"
load_benchmark_from_string(str) = load_benchmark(JSON.parse(str))

function load_benchmark(dict)
    # convert UTC time to formatted local time
    dt = DateTime(dict["time"])
    zdt = astimezone(ZonedDateTime(dt, tz"UTC"), localzone())
    dict["time"] = Dates.format(zdt, "yyyy-mm-dd HH:MM:SS")

    io = IOBuffer()
    JSON.print(io, dict["benchmark"])
    dict["benchmark"] = BenchmarkTools.load(IOBuffer(io.data))

    return dict
end

function safe_read_command(cmd)
    try
        chomp(read(cmd, String))
    catch
        ""
    end
end

struct hierarchical_benchmark
    trial
    indented_name
    children

    hierarchical_benchmark(trial, name) = new(trial, name, Dict())
end

components(x) = x[1] == "" ? [] : split(x[1], SEPARATOR)

"Reconstruct the testset hierarchy from deserialized benchmark data"
function get_testset_hierarchy(benchmark)
    tests = sort(collect(benchmark), by = x -> length(components(x)))
    root = hierarchical_benchmark(tests[1], "All tests")
    for t in tests[2:end]
        comps = components(t)
        parent = reduce((parent, child) -> parent.children[child], comps[1:end-1]; init = root)
        name = " "^(length(comps)*2) * comps[end]
        parent.children[comps[end]] = hierarchical_benchmark(t, name)
    end
    return root
end

"Hierarchically sort and iterate over benchmark trials"
function iterate_hierarchy(f, test)
    f(test)
    for child in sort([x for x in values(test.children)]; by = t -> -t.trial[2].time)
        iterate_hierarchy(f, child)
    end
end

const WHITE = 15

color(c, text) = @sprintf("\e[38;5;%dm%s\e[0m", c, text)

const TIME_TOLERANCE_PERCENT = [
    (-10, 10)
    ( -5, 114)
    (  5, 7)
    ( 10, 214)
    (Inf, 9)]

time_color(diff) =
    TIME_TOLERANCE_PERCENT[findfirst(t -> (isnan(diff) ? 0 : diff) < t[1], TIME_TOLERANCE_PERCENT)][2]

const MEMORY_TOLERANCE_PERCENT = [
    ( -5, 10)
    ( -1, 114)
    (  1, 7)
    (  5, 214)
    (Inf, 9)]

memory_color(diff) =
    MEMORY_TOLERANCE_PERCENT[findfirst(t -> (isnan(diff) ? 0 : diff) < t[1], MEMORY_TOLERANCE_PERCENT)][2]

function print_table_header(test_hierarchy, time_memory_width, time_padding, memory_padding)
    column_width = recursive_mapreduce(x -> length(x.indented_name), max, parent -> values(parent.children), test_hierarchy)
    column_title = "Test name"
    divider = color(WHITE, "─"^(column_width + time_memory_width)) * "\n"
    header = color(WHITE, @sprintf("%s%sTime%sMemory", column_title,
        " "^(column_width - length(column_title) + time_padding), " "^memory_padding))
    @printf("%s%s\n%s", divider, header, divider)
    return column_width, divider
end

function print_benchmark_details(title, d)
    runstr = d["time"]
    haskey(d, "runtime") && (runstr *= @sprintf("; runtime=%.0fs", d["runtime"]))
    haskey(d, "config") && (runstr *= "; $(d["config"])")
    if haskey(d, "julia-version")
        runstr *= "; Julia $(d["julia-version"])"
        haskey(d, "julia-commit") && (runstr *= " ($(d["julia-commit"]))")
    end

    println("\n",
        color(WHITE, title), "\n",
        color(WHITE, "Run:"), " $runstr\n",
        color(WHITE, "CPU:"), " $(d["cpu"])\n",
        color(WHITE, "Git:"), " $(d["git"])")
end

"Prints a single JSON file benchmark to stdout."
function print_benchmark(baseline_file)
    base_dict = load_benchmark_from_file(baseline_file)
    base = base_dict["benchmark"][1]

    root = get_testset_hierarchy(base)
    column_width, divider = print_table_header(root, 23, 6, 7)

    iterate_hierarchy(root) do test
        trial = test.trial[2]
        @printf("%s%s    %s     %s\n",
            test.indented_name,
            " "^(column_width - length(test.indented_name)),
            prettytime(trial.time, true), prettymemory(trial.memory, true))
    end

    print(divider)

    print_benchmark_details("Details", base_dict)
end

"Prints the comparison of two JSON file benchmarks to stdout."
function print_benchmark_comparison(baseline_file, alternative_file)
    base_dict = load_benchmark_from_file(baseline_file)
    base = base_dict["benchmark"][1]

    alt_dict = load_benchmark_from_file(alternative_file)
    alt = alt_dict["benchmark"][1]

    root = get_testset_hierarchy(alt)
    column_width, divider = print_table_header(root, 39, 9, 15)

    iterate_hierarchy(root) do test
        key = test.trial[1]
        trial = test.trial[2]

        if haskey(base, key)
            j = judge(trial, base[key])
            time_diff = 100 * (j.ratio.time - 1)
            memory_diff = 100 * (j.ratio.memory - 1)
        else
            time_diff = NaN
            memory_diff = NaN
        end

        @printf("%s%s    %s     %s\n",
            test.indented_name,
            " "^(column_width - length(test.indented_name)),
            color(time_color(time_diff), prettytime(trial.time, true) * " (" * prettypercent(time_diff, true) * ")"),
            color(memory_color(memory_diff), prettymemory(trial.memory, true) * " (" * prettypercent(memory_diff, true) * ")"))
    end

    print(divider)

    print_benchmark_details("Baseline", base_dict)
    print_benchmark_details("Alternative", alt_dict)
end


text:
"""
Transfer matrix for ZEUS experiment
high-x data.

N = T K S
T = Transfer matrix (including data luminosity)
K - Radiative corrections
S = Born-level integrated sigma

Based on C++ code by R. Aggarwal.
"""

include("../data/zeus_transfer_matrix.jl")
include("../data/zeus_bin_edges.jl")

export get_pred_N
export get_data_events
export get_TM_elements
export get_K_elements
export get_L_data
export get_bin_info

export q2bins_M_begin
export q2bins_M_end
export xbins_M_begin
export xbins_M_end

function get_data_events(eMPp)
    if eMPp == 0
        return Data_Events_ePp
    else
        return Data_Events_eMp
    end
end

function get_K_elements(eMPp)
    if eMPp == 0
        return Kii_ePp
    else
        return Kii_eMp
    end
end

function get_L_data(eMPp)
    if eMPp == 0
        return Ld_ePp
    else
        return Ld_eMp
    end
end

function get_TM_elements(eMPp)

    lumi_data = get_L_data(eMPp)

    if eMPp == 1
        TM_elements = copy(TM_Elements_eMp)
    elseif eMPp == 0
        TM_elements = copy(TM_Elements_ePp)
    end

    for i in 1:429
        for j in 1:153
         TM_elements[i,j] = TM_elements[i,j] * lumi_data
        end
    end

    return TM_elements
end

function get_pred_N(Integ_xsec, eMPp)

    TM = get_TM_elements(eMPp);
    K = get_K_elements(eMPp)

    xsec_pred = zeros(153)
    for j in 1:153
        temp=0.
        for i in 1:429

        temp = temp +(TM[i,j]./K[i]) * Integ_xsec[i];
        end
        xsec_pred[j] = temp
    end

    return xsec_pred;

end

"""
    get_bin_info(n)

Get the bin egdes of the ZEUS
detector space for a given
bin number, `n`.
"""
function get_bin_info(n::Integer)

    if n < 1 || n > 153
        @error "Bin number n should be [1, 153]"
    end

    @info "ZEUS detector bin" n BinQ2low[n] BinQ2high[n] Binxlow[n] Binxhigh[n]

    return ([BinQ2low[n], BinQ2high[n]], [Binxlow[n], Binxhigh[n]])
end

q2bins_M_begin = Float64.([400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 11000, 11000, 11000, 20000, 20000]);

q2bins_M_end = Float64.([650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 650, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 950, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1300, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1500, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2100, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2400, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 2800, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3200, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 3800, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 4500, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 6000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 11000, 11000, 11000, 11000, 11000, 11000, 11000, 11000, 11000, 11000, 11000, 11000, 20000, 20000, 20000, 30000, 30000]);

xbins_M_begin = Float64.([0.0035, 0.007, 0.01, 0.012, 0.015, 0.017, 0.019, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.041, 0.044, 0.047, 0.05, 0.054, 0.059, 0.062, 0.065, 0.0675, 0.07, 0.073, 0.076, 0.079, 0.082, 0.085, 0.09, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.24, 0.26, 0.28, 0.3, 0.35, 0.4, 0.42, 0.435, 0.45, 0.5, 0.6, 0.7, 0.0035, 0.007, 0.01, 0.012, 0.015, 0.017, 0.019, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.041, 0.044, 0.046, 0.05, 0.054, 0.06, 0.07, 0.079, 0.088, 0.098, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.21, 0.23, 0.26, 0.28, 0.31, 0.35, 0.4, 0.42, 0.435, 0.45, 0.5, 0.6, 0.007, 0.01, 0.012, 0.015, 0.018, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.041, 0.044, 0.046, 0.05, 0.059, 0.067, 0.076, 0.085, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.17, 0.19, 0.21, 0.23, 0.25, 0.28, 0.31, 0.35, 0.4, 0.425, 0.45, 0.55, 0.007, 0.01, 0.012, 0.015, 0.018, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.044, 0.05, 0.059, 0.067, 0.076, 0.085, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.205, 0.22, 0.24, 0.26, 0.28, 0.3, 0.34, 0.4, 0.5, 0.01, 0.015, 0.018, 0.021, 0.024, 0.029, 0.033, 0.038, 0.042, 0.048, 0.056, 0.063, 0.072, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.22, 0.24, 0.27, 0.3, 0.34, 0.4, 0.5, 0.01, 0.015, 0.017, 0.019, 0.022, 0.026, 0.03, 0.034, 0.038, 0.042, 0.048, 0.054, 0.062, 0.07, 0.078, 0.088, 0.1, 0.105, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.25, 0.3, 0.36, 0.5, 0.01, 0.018, 0.022, 0.025, 0.028, 0.032, 0.038, 0.044, 0.05, 0.059, 0.066, 0.076, 0.086, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.21, 0.23, 0.25, 0.28, 0.31, 0.34, 0.39, 0.45, 0.01, 0.022, 0.025, 0.028, 0.032, 0.038, 0.044, 0.054, 0.06, 0.07, 0.078, 0.088, 0.11, 0.14, 0.17, 0.21, 0.26, 0.31, 0.43, 0.01, 0.025, 0.029, 0.033, 0.039, 0.048, 0.058, 0.063, 0.08, 0.1, 0.13, 0.16, 0.2, 0.24, 0.29, 0.34, 0.46, 0.025, 0.028, 0.032, 0.038, 0.044, 0.05, 0.058, 0.065, 0.076, 0.086, 0.096, 0.12, 0.15, 0.18, 0.22, 0.27, 0.32, 0.38, 0.5, 0.025, 0.032, 0.038, 0.044, 0.05, 0.062, 0.072, 0.088, 0.11, 0.14, 0.17, 0.21, 0.25, 0.3, 0.36, 0.42, 0.54, 0.035, 0.042, 0.05, 0.06, 0.07, 0.085, 0.1, 0.115, 0.13, 0.16, 0.2, 0.24, 0.29, 0.34, 0.4, 0.46, 0.58, 0.035, 0.05, 0.07, 0.096, 0.12, 0.15, 0.19, 0.23, 0.28, 0.33, 0.39, 0.45, 0.51, 0.63, 0.035, 0.05, 0.07, 0.096, 0.12, 0.15, 0.18, 0.22, 0.27, 0.32, 0.38, 0.44, 0.5, 0.56, 0.69, 0.065, 0.1, 0.13, 0.16, 0.2, 0.24, 0.29, 0.35, 0.41, 0.47, 0.53, 0.59, 0.73, 0.105, 0.15, 0.19, 0.23, 0.28, 0.33, 0.39, 0.45, 0.51, 0.57, 0.64, 0.78, 0.13, 0.256, 0.6, 0.25, 0.6]);

xbins_M_end = Float64.([0.007, 0.01, 0.012, 0.015, 0.017, 0.019, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.041, 0.044, 0.047, 0.05, 0.054, 0.059, 0.062, 0.065, 0.0675, 0.07, 0.073, 0.076, 0.079, 0.082, 0.085, 0.09, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.24, 0.26, 0.28, 0.3, 0.35, 0.4, 0.42, 0.435, 0.45, 0.5, 0.6, 0.7, 0.999999, 0.007, 0.01, 0.012, 0.015, 0.017, 0.019, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.041, 0.044, 0.046, 0.05, 0.054, 0.06, 0.07, 0.079, 0.088, 0.098, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.21, 0.23, 0.26, 0.28, 0.31, 0.35, 0.4, 0.42, 0.435, 0.45, 0.5, 0.6, 0.999999, 0.01, 0.012, 0.015, 0.018, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.041, 0.044, 0.046, 0.05, 0.059, 0.067, 0.076, 0.085, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.17, 0.19, 0.21, 0.23, 0.25, 0.28, 0.31, 0.35, 0.4, 0.425, 0.45, 0.55, 0.99999, 0.01, 0.012, 0.015, 0.018, 0.022, 0.024, 0.026, 0.028, 0.03, 0.032, 0.034, 0.036, 0.038, 0.044, 0.05, 0.059, 0.067, 0.076, 0.085, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.205, 0.22, 0.24, 0.26, 0.28, 0.3, 0.34, 0.4, 0.5, 0.999999, 0.015, 0.018, 0.021, 0.024, 0.029, 0.033, 0.038, 0.042, 0.048, 0.056, 0.063, 0.072, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.22, 0.24, 0.27, 0.3, 0.34, 0.4, 0.5, 0.999999, 0.015, 0.017, 0.019, 0.022, 0.026, 0.03, 0.034, 0.038, 0.042, 0.048, 0.054, 0.062, 0.07, 0.078, 0.088, 0.1, 0.105, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.25, 0.3, 0.36, 0.5, 0.999999, 0.018, 0.022, 0.025, 0.028, 0.032, 0.038, 0.044, 0.05, 0.059, 0.066, 0.076, 0.086, 0.096, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.21, 0.23, 0.25, 0.28, 0.31, 0.34, 0.39, 0.45, 0.999999, 0.022, 0.025, 0.028, 0.032, 0.038, 0.044, 0.054, 0.06, 0.07, 0.078, 0.088, 0.11, 0.14, 0.17, 0.21, 0.26, 0.31, 0.43, 0.999999, 0.025, 0.029, 0.033, 0.039, 0.048, 0.058, 0.063, 0.08, 0.1, 0.13, 0.16, 0.2, 0.24, 0.29, 0.34, 0.46, 0.999999, 0.028, 0.032, 0.038, 0.044, 0.05, 0.058, 0.065, 0.076, 0.086, 0.096, 0.12, 0.15, 0.18, 0.22, 0.27, 0.32, 0.38, 0.5, 0.999999, 0.032, 0.038, 0.044, 0.05, 0.062, 0.072, 0.088, 0.11, 0.14, 0.17, 0.21, 0.25, 0.3, 0.36, 0.42, 0.54, 0.999999, 0.042, 0.05, 0.06, 0.07, 0.085, 0.1, 0.115, 0.13, 0.16, 0.2, 0.24, 0.29, 0.34, 0.4, 0.46, 0.58, 0.999999, 0.05, 0.07, 0.096, 0.12, 0.15, 0.19, 0.23, 0.28, 0.33, 0.39, 0.45, 0.51, 0.63, 0.999999, 0.05, 0.07, 0.096, 0.12, 0.15, 0.18, 0.22, 0.27, 0.32, 0.38, 0.44, 0.5, 0.56, 0.69, 0.99999, 0.1, 0.13, 0.16, 0.2, 0.24, 0.29, 0.35, 0.41, 0.47, 0.53, 0.59, 0.73, 1, 0.15, 0.19, 0.23, 0.28, 0.33, 0.39, 0.45, 0.51, 0.57, 0.64, 0.78, 0.99999, 0.256, 0.6, 0.99999, 0.6, 0.99999]);


text:
# VE execution support

export @veda, vedaconvert, vefunction

## high-level @cuda interface

"""
    @veda [kwargs...] func(args...)

High-level interface for executing code on a Vector Engine. The `@veda` macro should prefix a call,
with `func` a callable function or object that should return nothing. It will be compiled to
a VE function upon first use, and to a certain extent arguments will be converted and
managed automatically using `vedaconvert`. Finally, a call to `vedacall` is
performed, scheduling a kernel launch on the current VEDA context.

Several keyword arguments are supported that influence the behavior of `@veda`.
- `launch`: whether to launch this kernel, defaults to `true`. If `false` the returned
  kernel object should be launched by calling it and passing arguments again.
- `dynamic`: use dynamic parallelism to launch device-side kernels, defaults to `false`.
- arguments that influence kernel compilation: see [`vefunction`](@ref) and
  [`dynamic_vefunction`](@ref)
- arguments that influence kernel launch: see [`VEDA.HostKernel`](@ref) and
  [`VEDA.DeviceKernel`](@ref)
"""
macro veda(ex...)
    # destructure the `@veda` expression
    call = ex[end]
    kwargs = ex[1:end-1]

    # destructure the kernel call
    Meta.isexpr(call, :call) || throw(ArgumentError("second argument to @veda should be a function call"))
    f = call.args[1]
    args = call.args[2:end]

    code = quote end
    vars, var_exprs = assign_args!(code, args)

    # group keyword argument
    macro_kwargs, compiler_kwargs, call_kwargs, other_kwargs =
        split_kwargs(kwargs,
                     [:launch],
                     [:name, :global_hooks],
                     [:stream])
    if !isempty(other_kwargs)
        key,val = first(other_kwargs).args
        throw(ArgumentError("Unsupported keyword argument '$key'"))
    end

    # handle keyword arguments that influence the macro's behavior
    launch = true
    for kwarg in macro_kwargs
        key,val = kwarg.args
        if key == :launch
            isa(val, Bool) || throw(ArgumentError("`launch` keyword argument to @veda should be a constant value"))
            launch = val::Bool
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end
    if !launch && !isempty(call_kwargs)
        error("@veda with launch=false does not support launch-time keyword arguments; use them when calling the kernel")
    end

    # FIXME: macro hygiene wrt. escaping kwarg values (this broke with 1.5)
    #        we esc() the whole thing now, necessitating gensyms...
    @gensym f_var kernel_f kernel_args kernel_tt kernel
    # regular, host-side kernel launch
    #
    # convert the function, its arguments, call the compiler and launch the kernel
    # while keeping the original arguments alive
    push!(code.args,
        quote
            $f_var = $f
            GC.@preserve $(vars...) $f_var begin
                local $kernel_f = $vedaconvert($f_var)
                local $kernel_args = map($vedaconvert, ($(var_exprs...),))
                local $kernel_tt = Tuple{map(Core.Typeof, $kernel_args)...}
                local $kernel = $vefunction($kernel_f, $kernel_tt; $(compiler_kwargs...))
                if $launch
                    $kernel($(var_exprs...); $(call_kwargs...))
                end
                $kernel
            end
         end)
    return esc(code)
end

## host to device value conversion

struct Adaptor end

# convert host pointers to device pointers
Adapt.adapt_storage(to::Adaptor, p::VEPtr{T}) where {T} = reinterpret(Ptr{T}, p)

## Base.RefValue isn't VE compatible, so provide a compatible alternative
#struct VERefValue{T} <: Ref{T}
#  x::T
#end
#Base.getindex(r::VERefValue) = r.x
#Adapt.adapt_structure(to::Adaptor, r::Base.RefValue) = VERefValue(adapt(to, r[]))

#Adapt.adapt_structure(to::Adaptor, x::Base.RefValue{T}) where {T} =
#    Base.unsafe_convert(VERef{T}, x)

## broadcast sometimes passes a ref(type), resulting in a GPU-incompatible DataType box.
## avoid that by using a special kind of ref that knows about the boxed type.
#struct VERefType{T} <: Ref{DataType} end
#Base.getindex(r::VERefType{T}) where T = T
#Adapt.adapt_structure(to::Adaptor, r::Base.RefValue{<:Union{DataType,Type}}) = VERefType{r[]}()

"""
    vedaconvert(x)

This function is called for every argument to be passed to a kernel, allowing it to be
converted to a device-friendly format. By default, the function does nothing and returns the
input object `x` as-is.

Do not add methods to this function, but instead extend the underlying Adapt.jl package and
register methods for the the `VectorEngine.Adaptor` type.
"""
vedaconvert(arg) = adapt(Adaptor(), arg)

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

## abstract kernel functionality

abstract type AbstractKernel{F,TT} end

struct HostKernel{F,TT} <: AbstractKernel{F,TT}
    #ctx::VEContext
    mod::VEModule
    fun::VEFunction
end

"""
    (::HostKernel)(args...; kwargs...)
    (::DeviceKernel)(args...; kwargs...)

Low-level interface to call a compiled kernel, passing GPU-compatible arguments in `args`.
For a higher-level interface, use [`@veda`](@ref).

The following keyword arguments are supported:
- `stream` (defaults to the default stream)
"""
AbstractKernel

@generated function call(kernel::AbstractKernel{F,TT}, args...; call_kwargs...) where {F,TT}
    sig = Base.signature_type(F, TT)
    args = (:F, (:( args[$i] ) for i in 1:length(args))...)

    # filter out ghost arguments that shouldn't be passed
    #predicate = dt -> isghosttype(dt) || Core.Compiler.isconstType(dt)
    #to_pass = map(!predicate, sig.parameters)
    to_pass = map(!isghosttype, sig.parameters)
    call_t =                  Type[x[1] for x in zip(sig.parameters,  to_pass) if x[2]]
    call_args = Union{Expr,Symbol}[x[1] for x in zip(args, to_pass)            if x[2]]
    @debug("call_t=$call_t")
    @debug("call_args=$call_args")
    # replace non-isbits arguments (they should be unused, or compilation would have failed)
    # alternatively, allow `launch` with non-isbits arguments.
    # make sure the type is none of those passed on stack by VE.
    for (i,dt) in enumerate(call_t)
        if !isbitstype(dt) && !isstructtype(dt) && !(dt <: Base.RefValue)
            @debug("!isbitstype detected: argument $i $dt")
            call_t[i] = Ptr{Any}
            call_args[i] = :C_NULL
        end
    end

    # finalize types
    call_tt = Base.to_tuple_type(call_t)

    quote
        Base.@_inline_meta

        err, veargs = vecall(kernel, $call_tt, $(call_args...); call_kwargs...)
        return veargs
    end
end

using LLVM, LLVM.Interop # for the correct definition of isghosttype

"""
    vefunction(f, tt=Tuple{}; kwargs...)

Low-level interface to compile a function invocation for the currently-active GPU, returning
a callable kernel object.

The following keyword arguments are supported:
- `name`: overrides the name that the kernel will have in the generated code
- `global_hooks`: specifies maps from global variable name to initializer hook

The output of this function is automatically cached, i.e. you can simply call `vefunction`
in a hot path without degrading performance. New code will be generated automatically, when
function definitions change, or when different types or keyword arguments are provided.
"""
function vefunction(f::Core.Function, tt::Type=Tuple{}; name=nothing, device=0,
                    global_hooks=NamedTuple(), kwargs...)
    source = FunctionSpec(f, tt, true, name)
    cache = get!(()->Dict{UInt,Any}(), vefunction_cache, device)
    #isa = default_isa(device)
    target = VECompilerTarget()
    params = VECompilerParams(device, global_hooks, rv_loop_vectorize!)
    job = CompilerJob(target, source, params)
    GPUCompiler.cached_compilation(cache, job, vefunction_compile, vefunction_link)::HostKernel{f,tt}
end

const vefunction_cache = Dict{UInt,Dict{UInt,Any}}()

function vefunction_compile(@nospecialize(job::CompilerJob))
    # compile
    method_instance, mi_meta = GPUCompiler.emit_julia(job)
    ir, ir_meta = GPUCompiler.emit_llvm(job, method_instance; ctx=JuliaContext())
    kernel = ir_meta.entry
    if haskey(ENV, "JULIA_VE_IR") && ENV["JULIA_VE_IR"] == "1"
        @show ir
    end

    obj, obj_meta = GPUCompiler.emit_asm(job, ir; format=LLVM.API.LLVMObjectFile)

    # find undefined globals and calculate sizes
    globals = map(gbl->Symbol(LLVM.name(gbl))=>llvmsize(eltype(llvmtype(gbl))),
                  filter(x->isextinit(x), collect(LLVM.globals(ir))))
    entry = LLVM.name(kernel)
    dispose(ir)

    return (;obj, entry, globals)
end

function vefunction_link(@nospecialize(job::CompilerJob), compiled)
    device = job.params.device
    global_hooks = job.params.global_hooks
    obj, entry, globals = compiled.obj, compiled.entry, compiled.globals

    # create executable and kernel
    obj = codeunits(obj)
    mod = VEModule(obj)
    fun = VEFunction(mod, entry)
    kernel = HostKernel{job.source.f,job.source.tt}(mod, fun)

    # initialize globals from hooks
    for gname in first.(globals)
        hook = nothing
        if haskey(default_global_hooks, gname)
            hook = default_global_hooks[gname]
        elseif haskey(global_hooks, gname)
            hook = global_hooks[gname]
        end
        if hook !== nothing
            @debug "Initializing global $gname"
            gbl = get_global(exe, gname)
            hook(gbl, mod, device)
        else
            @debug "Uninitialized global $gname"
            continue
        end
    end

    return kernel
end

function (kernel::HostKernel)(args...; kwargs...)
    call(kernel, map(vedaconvert, args)...; kwargs...)
end

@inline function vecall(kernel::HostKernel, tt, args...; kwargs...)
    err, veargs = VEDA.vecall(kernel.fun, tt, args...; kwargs...)
end


text:
"""
kl\\_gmres(x0, b, atv, V, eta, ptv=nothing; kl_store=nothing;
             orth = "cgs2", side="right", lmaxit=-1, pdata=nothing)

Gmres linear solver. Handles preconditioning and restarts.
Uses gmres_base which is completely oblivious to these things.

The deal is

Input:\n
x0:  initial iterate, this is usually zero for nonlinear solvers

b: right hand side (duh!)

atv:  matrix-vector product which depends on precomputed data pdta
      I expect you to use pdata most or all of the time, so it is not
      an optional argument, even if it's nothing (at least for now).
      If your mat-vec is just A*v, you have to write a function where
      A is the precomputed data.
      API for atv is ```av=atv(v,pdata)```

V:  Preallocated n x K array for the Krylov vectors. I store the initial
    normalized residual in column 1, so  you have at most K-1 iterations
    before gmres\\_base returns a failure. kl\\_gmres will handle the
    restarts and, if lmaxit > 0, keep going until you hit lmaxit GMRES
    iterations. You may allocate V in Float32 and save on storage. The
    benefit from doing this is not dramatic in terms of CPU time.

eta: Termination happens when ||b - Ax|| <= eta || b ||

ptv:  preconditioner-vector product, which will also use pdata. The
      default is nothing, which is no preconditioning at all.
      API for ptv is px=ptv(x,pdata)

Keyword arguments

kl\\_store: You have the option of giving me some room
         for the vectors gmres needs. These include copies of x0 and b,
         which I will not overwrite and a couple of vectors I use
         in the iteration. If you're only doing a linear solve, it
         does no harm to let me allocate those vectores in kl\\_gmres.
         The way to do this is ```kl_store=kstore(n,"gmres")``` where n
         is the number of unknows. I call this myself in the initialization
         phase if you don't do it ahead of me.

pdata: precomputed data. The default is nothing, but that ain't gonna
        work well for nonlinear equations.

orth: your choice of the wise default, classical Gram-Schmidt twice,
       or something slower and less stable. Those are classical once (really
       bad) or a couple variants of modified Gram-Schmidt. mgs2 is what I
       used in my old matlab codes. Not terrible, but far from great.

side: left or right preconditioning. The default is "right".

lmaxit: maximum number of linear iterations. The default is -1, which
        means that the maximum number of linear iterations is K-1, which
        is all V will allow without restarts. If lmaxit > K-1, then the
        iteration will restart until you consume lmaxit iterations or
        terminate successfully.

Other parameters on the way.

Output:\n
A named tuple (sol, reshist, lits, idid)

where

sol= final result
reshist = residual norm history
lits = number of iterations
idid = status of the iteration
       true -> converged
       false -> failed to converge

### Examples from the docstrings for kl\\_gmres

In these examples you have the matrix and use
```
function atv(x, A)
    return A * x
end
```
to compute the matvec.

#### Three dimensional problem.

Will converge in the correct three iterations
only if you orthogonalize with CGS twice.

```jldoctest
julia> function atv(x, A)
           return A * x
       end
atv (generic function with 1 method)

julia> A = [0.001 0 0; 0 0.0011 0; 0 0 1.e4];

julia> V = zeros(3, 10); b = [1.0; 1.0; 1.0]; x0 = zeros(3);

julia> gout = kl_gmres(x0, b, atv, V, 1.e-10; pdata = A);

julia> gout.reshist
4-element Array{Float64,1}:
 1.73205e+00
 1.41421e+00
 6.72673e-02
 1.97712e-34

julia> norm(b - A*gout.sol,Inf)
1.28536e-10
```

#### Integral equation. Notice that pdata has the kernel of the
operator and we do the matvec directly. Just like the previous example.
We put the grid information and, for this artifical example, the solution
in the precomputed data.

```jldoctest
julia> function integop(u, pdata)
           K = pdata.K
           return u - K * u
       end
integop (generic function with 1 method)

julia> function integopinit(n)
           h = 1 / n
           X = collect(0.5*h:h:1.0-0.5*h)
           K = [ker(x, y) for x in X, y in X]
           K .*= h
           sol = [usol(x) for x in X]
           f = sol - K * sol
           pdata = (K = K, xe = sol, f = f)
           return pdata
       end
integopinit (generic function with 1 method)

julia> function usol(x)
           return exp.(x) .* log.(2.0 * x .+ 1.0)
       end
usol (generic function with 1 method)

julia> function ker(x, y)
           ker = 0.1 * sin(x + exp(y))
       end
ker (generic function with 1 method)

julia> n=100; pdata = integopinit(n); ue = pdata.xe; f=pdata.f;

julia> u0 = zeros(size(f)); V = zeros(n, 20); V32=zeros(Float32,n,20);

julia> gout = kl_gmres(u0, f, integop, V, 1.e-10; pdata = pdata);

julia> gout32 = kl_gmres(u0, f, integop, V32, 1.e-10; pdata = pdata);

julia> [norm(gout.sol-ue,Inf) norm(gout32.sol-ue,Inf)]
1×2 Array{Float64,2}:
 4.44089e-16  2.93700e-07

julia> [gout.reshist gout32.reshist]
4×2 Array{Float64,2}:
 1.48252e+01  1.48252e+01
 5.52337e-01  5.52337e-01
 1.77741e-03  1.77742e-03
 1.29876e-19  8.73568e-11
```

"""
function kl_gmres(
    x0,
    b,
    atv,
    V,
    eta,
    ptv = nothing;
    kl_store = nothing,
    orth = "cgs2",
    side = "right",
    lmaxit = -1,
    pdata = nothing,
)

    # Build some precomputed data to inform KL_atv about
    # preconditioning ...
    # Do not overwrite the initial iterate or the right hand side.
    n = length(x0)
    # Get the vectors GMRES needs internally and make room to
    # copy the initial iterate and right side
    (kl_store !== nothing) || (kl_store=kstore(n,"gmres"))
        y0 = kl_store[1]
        y0 .= x0
        rhs = kl_store[2]
        rhs .= b
    # Two vectors for internals
        linsol = kl_store[3]
        restmp = kl_store[4]
    #
    if side == "right" || ptv == nothing
        itsleft = false
    else
        itsleft = true
        rhs .= ptv(rhs, pdata)
    end
    (n, K) = size(V)
    K > 1 || error("Must allocate for GMRES iterations. V must have
                   at least two columns")
    klmaxit = lmaxit
    lmaxit > 0 || (lmaxit = K - 1)
    #
    itvec = maxitvec(K, lmaxit)
    ip = 1
    idid = false
    Kpdata =
        (pdata = pdata, side = side, ptv = ptv, atv = atv, linsol = linsol, restmp = restmp)
    gout = []
    #
    # Restarted GMRES loop.
    #
    while ip <= length(itvec) && idid == false
        localout =
            gmres_base(y0, rhs, Katv, V, eta, Kpdata; lmaxit = itvec[ip], orth = orth)
        idid = localout.idid
        gout = outup(gout, localout, ip, klmaxit)
        reslen = length(localout.reshist)
        #
        # Update the termination criterion for the restart.
        # gmres_base overwrites y0 with the solution
        #
        idid || (eta = eta * localout.rho0 / localout.reshist[reslen])
        ip += 1
    end
    #
    # Fixup the solution if preconditioning from the right.
    #
    sol = y0
    if side == "left" || ptv == nothing
        return (sol = sol, reshist = gout.reshist, lits = gout.lits, idid = gout.idid)
    else
        sol .= ptv(sol, pdata)
        return (sol = sol, reshist = gout.reshist, lits = gout.lits, idid = gout.idid)
    end
end

"""
Katv(x,Kpdata)

Builds a matrix-vector product to hand to gmres_base or bicgstab_base.
Puts the preconditioner in there on the correct side.
"""
function Katv(x, Kpdata)
    #    y=copy(x)
    y = Kpdata.linsol
    pdata = Kpdata.pdata
    ptv = Kpdata.ptv
    atv = Kpdata.atv
    side = Kpdata.side
    sideok = (side == "left") || (side == "right")
    sideok || error(
        "Bad preconditioner side in Krylov solver, input side = ",
        side,
        ". Side must be \"left\" or \"right\" ",
    )
    if ptv == nothing
        y .= atv(x, pdata)
        return y
    elseif side == "left"
        y .= atv(x, pdata)
        return ptv(y, pdata)
    elseif side == "right"
        y .= ptv(x, pdata)
        return atv(y, pdata)
    end
end

"""
gmres_base(x0, b, atv, V, eta, pdata; orth="cgs2", lmaxit=-1)

Base GMRES solver. This is GMRES(m) with no restarts and no preconditioning.
The idea for the future is that it'll be called by kl_gmres (linear
solver) which is the backend of klgmres.

gmres_base overwrites x0 with the solution. This is one of many reasons
that you should not invoke it directly.
"""
function gmres_base(x0, b, atv, V, eta, pdata; orth = "cgs2", lmaxit = -1)

    (n, m) = size(V)
    #
    # Allocate for Givens
    #
    #    kmax = m - 1
    kmax = m
    lmaxit == -1 || (kmax = lmaxit)
    kmax > m - 1 && error("lmaxit error in gmres_base")
    r = pdata.restmp
    r .= b
    T = eltype(V)
    h = zeros(T, kmax + 1, kmax + 1)
    c = zeros(kmax + 1)
    s = zeros(kmax + 1)
    #
    # Don't do the mat-vec if the intial iterate is zero
    #
    y = pdata.linsol
    (norm(x0) == 0.0) || (r .-= atv(x0, pdata))
    #    (norm(x0) == 0.0) || (y .= atv(x0, pdata); r .-=y;)
    #
    #
    rho0 = norm(r)
    rho = rho0
    #
    # Initial residual = 0? This can't be good.
    #
    rho == 0.0 && error("Initial resdiual in kl_gmres is zero. Why?")
    #
    g = zeros(size(c))
    g[1] = rho
    errtol = eta * rho
    reshist = []
    #
    # Initialize
    #
    idid = true
    push!(reshist, rho)
    k = 0
    #
    # Showtime!
    #
    @views V[:, 1] .= r / rho
    beta = rho
    while (rho > errtol) && (k < kmax)
        k += 1
        @views V[:, k+1] .= atv(V[:, k], pdata)
        @views vv = vec(V[:, k+1])
        @views hv = vec(h[1:k+1, k])
        @views Vkm = V[:, 1:k]
        #
        # Don't mourn. Orthogonalize!
        #
        Orthogonalize!(Vkm, hv, vv, orth)
        #
        # Build information for new Givens rotations.
        #
        if k > 1
            hv = @view h[1:k, k]
            giveapp!(c[1:k-1], s[1:k-1], hv, k - 1)
        end
        nu = norm(h[k:k+1, k])
        if nu != 0
            c[k] = conj(h[k, k] / nu)
            s[k] = -h[k+1, k] / nu
            h[k, k] = c[k] * h[k, k] - s[k] * h[k+1, k]
            h[k+1, k] = 0.0
            gv = @view g[k:k+1]
            giveapp!(c[k], s[k], gv, 1)
        end
        #
        # Update the residual norm.
        #
        rho = abs(g[k+1])
        push!(reshist, rho)
    end
    #
    # At this point either k = kmax or rho < errtol.
    # It's time to compute x and check out.
    #
    y = h[1:k, 1:k] \ g[1:k]
    qmf = view(V, 1:n, 1:k)
    #    mul!(r, qmf, y)
    #    r .= qmf*y
    #    x .+= r
    #    sol = x0
    #    mul!(sol, qmf, y, 1.0, 1.0)
    mul!(x0, qmf, y, 1.0, 1.0)
    (rho <= errtol) || (idid = false)
    k > 0 || println("GMRES iteration terminates on entry.")
    return (rho0 = rho0, reshist = Float64.(reshist), lits = k, idid = idid)
end

function giveapp!(c, s, vin, k)
    for i = 1:k
        w1 = c[i] * vin[i] - s[i] * vin[i+1]
        w2 = s[i] * vin[i] + c[i] * vin[i+1]
        vin[i:i+1] .= [w1, w2]
    end
    return vin
end

#
# The functions maxitvec and outup manage the restarts.
# There is no reason to look at them or fiddle with them.
#

function maxitvec(K, lmaxit)
    levels = Int.(ceil(lmaxit / (K - 1)))
    itvec = ones(Int, levels)
    itvec[1:levels-1] .= K - 1
    remainder = lmaxit - (levels - 1) * (K - 1)
    itvec[levels] = remainder
    return itvec
end

function outup(gout, localout, ip, klmaxit)
    idid = localout.idid
    #
    # If I'm doing restarts I won't store the last residual
    # unless the iteration is successful. The reason is that
    # I will add that residual to the list when I restart.
    #
    if idid || klmaxit == -1
        lreshist = localout.reshist
    else
        lk = length(localout.reshist)
        lreshist = localout.reshist[1:lk-1]
    end
    if ip == 1
        reshist = lreshist
        lits = localout.lits
    else
        reshist = gout.reshist
        append!(reshist, lreshist)
        lits = gout.lits + localout.lits
    end
    gout = (reshist = reshist, lits = lits, idid = idid)
    return gout
end


text:
"""
    systematic_sample(N, d=Normal(0,1); permute=true)
returns a `Vector` of length `N` sampled systematically from the distribution `d`. If `permute=false`, this vector will be sorted.
"""
function systematic_sample(N, d=Normal(0,1); permute=true)
    e   = 0.5/N # rand()/N
    y   = e:1/N:1
    o = map(y) do y
        quantile(d,y)
    end
    permute && permute!(o, randperm(N))
    o
end


text:
module Dagitty

include("DAG.jl")
include("drawdag.jl")
include("is_d_separated.jl")
include("implied_conditional_independencies.jl")
include("utils.jl")
include("adj_sets.jl")

end


text:
"""
   getRecord retrieves data out of an event for a given recordSpec. Note that
   the handle that is returned must be treated as a C++ object (e.g. with `icxx`)
"""
function getRecord(he::HeistEvent, rs::HeistRecordSpec)

   # The below builds cxxt"SomeType"
   # See https://discourse.julialang.org/t/how-to-build-a-string-for-icxx/7354
   s = @eval @cxxt_str $(rs.productName)

   return icxx"return $(he.galleryEvent).getValidHandle<$s>($(rs.artInputTag));"
end


text:
const julia_getModuleAt_request_type = JSONRPC.RequestType("julia/getModuleAt", VersionedTextDocumentPositionParams, String)
const julia_getCurrentBlockRange_request_type = JSONRPC.RequestType("julia/getCurrentBlockRange", VersionedTextDocumentPositionParams, Position)
const julia_getDocAt_request_type = JSONRPC.RequestType("julia/getDocAt", VersionedTextDocumentPositionParams, String)


text:
using JunoDocs
using Test

#Pkg.checkout("Juno")


text:
using JuMP
using SCIP

mutable struct Problem
    model
    x
    y
end

function make_model(nonlin=false)::Problem
    m = Model(SCIP.Optimizer)

    @variable(m, x[1:9], Bin)
    @variable(m, -1 ≤ y[1:9] ≤ 1)

    @constraint(m, c1[i=1:9], x[i] <= 3y[i])
    if nonlin
        @NLconstraint(m, sum(y[i] for i=1:9)^2 == prod(y[i] for i=1:9))
    end

    @objective(m, Max, sum(y[i] - i*x[i] for i=1:9))

    Problem(m, x, y)
end

function add_lazycons(prob::Problem)
    # some state for the lazy callback
    counter = 1

    function lazycb(cb)
        println("### lazycb is called! ($counter)")
        @lazyconstraint(cb, prob.x[counter] == prob.y[counter])

        counter = counter % 9 + 1
    end

    addlazycallback(prob.model, lazycb)
end

println("<<< First try without nonlinear constraints:")
prob = make_model(false)
add_lazycons(prob)
solve(prob.model)

println("<<< Now try with nonlinear constraints:")
prob = make_model(true)
add_lazycons(prob)
solve(prob.model)


text:
# SPDX-License-Identifier: MIT

#======================================================================
Forward recursion
======================================================================#

function αrecursion(π::AbstractVector{SR},
                    Tᵀ::AbstractMatrix{SR},
                    lhs::AbstractMatrix{SR}) where SR <: Semiring
    S, N = length(π), size(lhs, 2)
    α = similar(lhs, SR, S, N)
    buffer = similar(lhs[:, 1])

    @views elmul!(α[:,1], π, lhs[:,1])
    @views for n in 2:N
        matmul!(buffer, Tᵀ, α[:,n-1])
        elmul!(α[:,n], buffer, lhs[:,n])
    end
    α
end

#======================================================================
Backward recursion
======================================================================#

function βrecursion(T::AbstractMatrix{SR},
                    lhs::AbstractMatrix{SR}) where SR <: Semiring
    S, N = size(T, 1), size(lhs, 2)
    β = fill!(similar(lhs, SR, S, N), one(SR))
    buffer = similar(lhs[:, 1], S)

    @views for n in N-1:-1:1
        elmul!(buffer, β[:,n+1], lhs[:,n+1])
        matmul!(β[:,n], T, buffer)
    end
    β
end

# Do the backward recursion and multiply with the α immediately.
# This avoids to allocate a second array for the β messages.
function βrecursion_mulα!(α::AbstractMatrix{SR}, T::AbstractMatrix{SR},
                    lhs::AbstractMatrix{SR}) where SR <: Semiring
    S, N = size(T, 1), size(lhs, 2)
    βₘ = fill!(similar(lhs[:, 1], S), one(SR))
    buffer = similar(lhs[:, 1], S)

    @views for n in N-1:-1:1
        elmul!(buffer, βₘ, lhs[:,n+1])
        matmul!(βₘ, T, buffer)
        elmul!(α[:,n], α[:,n], βₘ)
    end
    α
end

#======================================================================
Specialized algorithms
======================================================================#

function _expand(lhs::AbstractMatrix{T}, seqlength = size(lhs, 2)) where T
    S, N = size(lhs)
    retval = hcat(vcat(lhs, zeros(T, 1, N)), zeros(T, S+1))
    @views fill!(retval[1:end-1,seqlength+1:end], zero(T))
    @views fill!(retval[end,seqlength+1:end], one(T))
    retval
end

function _split(fsm, lhs)
	offset = 0
	retval = []
	for r in fsm.ranges
		nr = (r[1]:r[end]-1) .- offset
		push!(retval, lhs[nr, :])
		offset += 1
	end
	retval
end

function _drop_extradims(fsm, γ)
    [γ[r[1]:r[end]-1, 1:end-1] for r in fsm.ranges]
end

"""
    pdfposteriors(mfsm::MatrixFSM{SF} lhs) where SF <: Semifield
    pdfposteriors(union(mfsm1, mfsm2, ...), batch_lhs)

Calculate the conditional posterior of "assigning" the \$n\$th frame
to the \$i\$th pdf. The output is a tuple `γ, ttl` where `γ` is a matrix
(# pdf x # frames) and `ttl` is the total probability of the sequence.
This function can also be called in "batch-mode" by providing a union
of compiled fsm and a 3D tensor containing the per-pdf, per-frame
and per-batch values.
"""
pdfposteriors

function pdfposteriors(mfsm::MatrixFSM{SR},
                       in_lhs::AbstractMatrix{T}) where {SR<:LogSemifield,T}
    # Convert the likelihood matrix to the mfsm' semiring.
    lhs = copyto!(similar(in_lhs, SR), in_lhs)

    # Add one frame and one dimension for the final state.
    expanded_lhs = _expand(lhs)

    S = size(mfsm.C, 1)
    N = size(expanded_lhs, 2)

    # Get the per-state likelihood matrix.
    state_lhs = matmul!(similar(lhs, S, N), mfsm.C, expanded_lhs)

    α = αrecursion(mfsm.π, mfsm.Tᵀ, state_lhs)
    state_γ = βrecursion_mulα!(α, mfsm.T, state_lhs)

    # Transform the per-state γs to per-likelihoods γs.
    γ = matmul!(expanded_lhs, mfsm.Cᵀ, state_γ) # re-use `expanded_lhs` memory.

    # Re-normalize the γs.
    sums = sum(γ, dims = 1)
    eldiv!(γ, γ, sums)

    # Convert the result to the Real-semiring
    out = copyto!(similar(in_lhs), γ[1:end-1,1:end-1])
    out = map!(exp, out, out)

    out, convert(T, minimum(sums))
end

function pdfposteriors(mfsm::UnionMatrixFSM{SR},
                       in_lhs::AbstractArray{T,3},
                       seqlengths,
                      ) where {SR<:LogSemifield,T}

    # Convert the likelihood tensor to the mfsm's semiring.
    lhs_tensor = copyto!(similar(in_lhs, SR), in_lhs)

    # Reshape the 3D tensor to have a matrix of size BK x N where
    # B is the number of elements in the batch.
    lhs = vcat(
        map(t -> _expand(t...),
            zip(eachslice(lhs_tensor, dims = 3), seqlengths))...
    )

    S = size(mfsm.C, 1)       # number of states
    K = size(in_lhs, 1) + 1   # number of pdfs
    N = size(in_lhs, 2) + 1   # number of frames

    # Get the per-state likelihoods.
    state_lhs = matmul!(similar(lhs, SR, S, N), mfsm.C, lhs)

    α = αrecursion(mfsm.π, mfsm.Tᵀ, state_lhs)
    state_γ = βrecursion_mulα!(α, mfsm.T, state_lhs)

    # Transform the per-state γs to per-likelihoods γs.
    γ = matmul!(lhs, mfsm.Cᵀ, state_γ) # re-use `lhs` memory.
    γ = permutedims(reshape(γ, K, :, N), (1, 3, 2))

    # Re-normalize each element of the batch.
    sums = sum(γ, dims = 1)
    eldiv!(γ, γ, sums)

    # Convert the result to the Real-semiring
    out = copyto!(similar(in_lhs), γ[1:end-1,1:end-1,:])
    out = map!(exp, out, out)

    ttl = dropdims(minimum(sums, dims = (1, 2)), dims = (1, 2))
    out, copyto!(similar(ttl, T), ttl)
end

@deprecate stateposteriors(mfsm, lhs) pdfposteriors(mfsm, lhs)

"""
    maxstateposteriors(mfsm, lhs)

Calculate the posterior of "assigning" the \$n\$th frame
to the \$i\$th state conditioned on all other states maximizing the
likelihood of the sequence. The output is a matrix `μ` (# pdf x # frames).
"""
function maxstateposteriors(mfsm::MatrixFSM{SR},
                            in_lhs::AbstractArray) where SR<:TropicalSemiring
    # Convert the FSM and the likelihood matrix to the Log-semifield.
    lhs = copyto!(similar(in_lhs, SR), in_lhs)

    # Add one frame and one dimension for the final state.
    expanded_lhs = _expand(lhs)

    S = size(mfsm.C, 1)
    N = size(expanded_lhs, 2)

    # Expand the likelihood matrix to get the per-state likelihoods.
    state_lhs = matmul!(similar(lhs, S, N), mfsm.C, expanded_lhs)

    α = αrecursion(mfsm.π, mfsm.Tᵀ, state_lhs)
    β = βrecursion(mfsm.T, state_lhs)
    elmul!(similar(state_lhs), α, β)[1:end-1,1:end-1]
end

function findmaxstate(mfsm::MatrixFSM{SR}, μ::AbstractArray{SR},
                      prevstate::Int) where {SR,T}
    m = typemin(SR)
    s = -1
    for (i, μᵢ) in pairs(μ)
        if μᵢ >= m && mfsm.T[prevstate,i] ≠ zero(SR)
            m = μᵢ
            s = i
        end
    end
    s
end

"""
    bestpath(mfsm, μ)

Return the sequence of state with the highest value from `μ`, where
`μ` is the ouput of [`maxstateposteriors`](@ref).
"""
function bestpath(mfsm::MatrixFSM{SR}, μ::AbstractArray{SR}) where SR
    seq = Vector{Int}(undef, size(μ, 2))
    seq[1] = argmax(μ[:,1])
    for n in 2:size(μ,2)
        seq[n] = findmaxstate(mfsm, μ[:,n], seq[n-1])
    end
    seq
end


text:
#EasyPlot utilities to handle multi-dimensional datasets
#-------------------------------------------------------------------------------

#==NOTE:
This module includes tools to use in case the plotting (rendering) tool does
not support multi-dimensional datasets directly.

 - Also provides support for folded x-axis scale.
 - iparam: Parameter index (used to pick color)
==#

#==Interface to build plots on different backends.
===============================================================================#

"""
`abstract type AbstractWfrmBuilder`

EasyPlot provides:
 - `addwfrm(b::AbstractWfrmBuilder, wfrm::Waveform, wfrmidx::Int)`

Must implement:
 - addwfrm(b::T, d::DataF1, label::String, l::LineAttributes, g::GlyphAttributes, strip::Int) where T<:AbstractWfrmBuilder
 - needsfold(b::T) where T<:AbstractWfrmBuilder -> FoldedAxis() or nothing
"""
abstract type AbstractWfrmBuilder; end

#Register symbols that need to be implemented by backend interfacing code
addwfrm(b::AbstractWfrmBuilder, d::DataF1, label::String, l::LineAttributes, g::GlyphAttributes, strip::Int) =
	throw(MethodError(addwfrm, (b, d, label, l, g, strip)))
needsfold(b::AbstractWfrmBuilder) = throw(MethodError(needsfold, (b,)))

#==Catch missing implementations
===============================================================================#
_addwfrm(b::AbstractWfrmBuilder, d::T, args...; kwargs...) where T<:DataMD =
	throw(ArgumentError("Plotting if $T datasets is not supported."))

#==Helper functions
===============================================================================#
_get_crnid(v::Integer) = string(v)
_get_crnid(v::Real) = SI(v, ndigits=3) #Limit significant digits for readability

#==Adding DataF1 waveforms
===============================================================================#

#AbstractWfrmBuilder requested to fold data on axis:
function _addwfrm_fold(b::AbstractWfrmBuilder, f::FoldedAxis, d::DataF1,
   id::String, l::LineAttributes, g::GlyphAttributes, strip::Int)
	xspan = f.xmax - f.xmin
	if f.xmin != 0
		@warn("buildeye(): Does not currently support non-zero start time")
	end
	eye = buildeye(d, f.foldinterval, xspan, tstart=f.xstart)

	cur_id = id
	for segment in eye.data
		addwfrm(b, segment, cur_id, l, g, strip)
		cur_id = "" #Should avoid clutter in legend
	end
end

#Add a waveform to a typical plot:
function _addwfrm(b::AbstractWfrmBuilder, d::DataF1,
   id::String, l::LineAttributes, g::GlyphAttributes, strip::Int, iparam::Int)
	cur_id = l
	if nothing == cur_id.color
		cur_id = LineAttributes(l.style, l.width, iparam)
	end

	faxis = needsfold(b)
	if isnothing(faxis)
		addwfrm(b, d, id, cur_id, g, strip)
	else
		_addwfrm_fold(b, faxis, d, id, cur_id, g, strip)
	end
end

#==Adding DataRS waveforms
===============================================================================#

#Add collection of DataRS{DataF1} results:
function _addwfrm(b::AbstractWfrmBuilder, d::DataRS{DataF1}, crnid::String, id::String,
	l::LineAttributes, g::GlyphAttributes, strip::Int, iparam::Int)
	crnid = (""==crnid) ? crnid : "$crnid / "
	sweepname = d.sweep.id
	for i in 1:length(d.elem)
		v = d.sweep.v[i]
		curcrnid = "$crnid$sweepname=" * _get_crnid(v)
		cur_id = "$id; $curcrnid"
		wfrm = _addwfrm(b, d.elem[i], cur_id, l, g, strip, i)
	end
end

#Add collection of DataRS{Number} results:
function _addwfrm(b::AbstractWfrmBuilder, d::DataRS{T}, crnid::String, id::String,
	l::LineAttributes, g::GlyphAttributes, strip::Int, iparam::Int) where T<:Number
	cur_id = (""==crnid) ? id : "$id; $crnid"
	return _addwfrm(b, DataF1(d.sweep.v, d.elem), cur_id, l, g, strip, iparam)
end

#Add collection of DataRS{DataRS} results:
function _addwfrm(b::AbstractWfrmBuilder, d::DataRS{DataRS}, crnid::String, id::String,
	l::LineAttributes, g::GlyphAttributes, strip::Int, iparam::Int)
	crnid = (""==crnid) ? crnid : "$crnid / "
	sweepname = d.sweep.id
	for i in 1:length(d.elem)
		v = d.sweep.v[i]
		curcrnid = "$crnid$sweepname=" * _get_crnid(v)
		wfrm = _addwfrm(b, d.elem[i], curcrnid, id, l, g, strip, i)
	end
end

#If corner id !exists, use "" & relay call:
_addwfrm(b::AbstractWfrmBuilder, d::DataRS, id::String,
	l::LineAttributes, g::GlyphAttributes, strip::Int, iparam::Int) =
	_addwfrm(b, d, "", id, l, g, strip, iparam)

#==Adding DataHR waveforms
===============================================================================#

#Add waveforms from a collection of DataHR{DataF1}:
function _addwfrm(b::AbstractWfrmBuilder, d::DataHR{DataF1}, id::String,
	l::LineAttributes, g::GlyphAttributes, strip::Int, wfrmidx::Int)
	sweepnames = names(sweeps(d))

	for inds in subscripts(d)
		values = coordinates(d, inds)
		crnid=join(["$k="*_get_crnid(v) for (k,v) in zip(sweepnames,values)], " / ")
		cur_id = "$id; $crnid"
		cur_wfrmidx = inds[end]
		wfrm = _addwfrm(b, d.elem[inds...], cur_id, l, g, strip, cur_wfrmidx)
	end
end

#Add waveforms from a collection of DataHR{DataF1}:
#(Convert DataHR{Number} => DataHR{DataF1} & add):
function _addwfrm(b::AbstractWfrmBuilder, d::DataHR{T}, id::String,
	l::LineAttributes, g::GlyphAttributes, strip::Int, wfrmidx::Int) where T<:Number
	return _addwfrm(b, DataHR{DataF1}(d), id, l, g, strip, wfrmidx)
end

#==External interface (for backend implementation)
===============================================================================#
"""
    addwfrm(b::AbstractWfrmBuilder, wfrm::Waveform, wfrmidx::Int)

Have `EasyPlot` module break down a multi-dimensional `Waveform` into individual
traces, and call upon `AbstractWfrmBuilder` at the end to add them one-by-one.

Also segment data again in order to emulate a folded x-axis if `::FoldedAxis`
is returned by:
```julia-repl
needsfold(b::T) where T<:AbstractWfrmBuilder

# Notes
 - `wfrmidx` identifies order of added waveform for color selection purposes.
 - `wfrmidx` might get reset to 1 for each subplot, or might be continuous across all plots.
"""
addwfrm(b::AbstractWfrmBuilder, wfrm::Waveform, wfrmidx::Int) =
	_addwfrm(b, wfrm.data, wfrm.label, wfrm.line, wfrm.glyph, wfrm.strip, wfrmidx)

#Last line


text:
using QuasiArrays, DomainSets, Test
import QuasiArrays: ApplyQuasiMatrix, UnionVcat

@testset "Concatenation" begin
    @testset "Hcat" begin
        A = QuasiArray(randn(2,2), ([0,0.5], Base.OneTo(2)))
        B = [A A]
        @test B isa ApplyQuasiMatrix{Float64,typeof(hcat)}
        @test axes(B) == (axes(A,1), Base.OneTo(4))
        @test B[0.0,1] == B[0.0,3] == A[0.0,1]
        @test B == B
        @test B ≠ A
        @test B ≠ [A A A]
        @test_throws BoundsError B[0.1,1]
        @test_throws BoundsError B[0.0,5]
        @test QuasiArray(B) == B
    end

    @testset "Inclusion Union" begin
        a = Inclusion(0:0.5:2)
        b = Inclusion([1,3,4])
        c = union(a,b)
        @test 0.5 in c
        @test 3 in c
        @test !(2.5 in c)
        @test c[0.5] ≡ 0.5
        @test c[3] ≡ c[3.0] ≡ 3.0
        @test_throws BoundsError c[2.5]

        @test Inclusion{Float64}([1,2]) ∪ Inclusion{Float64}([3,4]) isa Inclusion{Float64}
    end

    @testset "UnionVcat" begin
        @testset "vector" begin
            a = QuasiVector(randn(5), 0:0.5:2)
            b = QuasiVector([5,6,8], [1,3,4])
            c = UnionVcat(a,b)
            @test axes(c,1) == Inclusion(UnionDomain(0:0.5:2, [1,3,4]))
            @test c[0.5] ≡ a[0.5]
            @test c[3] ≡ c[3.0] ≡ 6.0
            @test_throws BoundsError c[2.5]
        end
        @testset "matrix" begin
            A = QuasiArray(randn(5,2), (0:0.5:2, [1,4]))
            B = QuasiArray(randn(3,2), (3:0.5:4, [1,4]))
            C = UnionVcat(A, B)
            @test C[0,1] == A[0,1]
            @test C[0.5,1] == A[0.5,1]
            @test C[0.5,4] == A[0.5,4]
            @test C[3,1] == B[3,1]
            @test C[3.5,1] == B[3.5,1]
            @test C[3.5,4] == B[3.5,4]
            @test_throws BoundsError C[2.5,1]
            @test_throws BoundsError C[3,2]

            @test_throws ArgumentError UnionVcat(A, QuasiArray(randn(3,1), (3:0.5:4, [1])))
        end
    end
end


text:
module TestTimescaleDB

using Test
using TSML

function test_timescaledb()
  tsdb = TimescaleDB()
  mpipeline = tsdb
  try
    respipe = fit_transform!(mpipeline)
    @test sum(respipe.Value) ≈ 560997.0
    @test nrow(respipe) == 48349
    @test ncol(respipe) == 2
    respipe = fit_transform(mpipeline)
    @test sum(respipe.Value) ≈ 560997.0
    @test nrow(respipe) == 48349
    @test ncol(respipe) == 2
  catch e
  end
end
@testset "TimescaleDB" begin
  test_timescaledb()
end

end


text:
mutable struct SGDParams
    lr::Float64
    lr₀::Float64
    N_sgd::Int
    N_mc::Int
    N_mc0::Int
end

function SGDParams(;lr=0.1, N_sgd=1, N_mc=10^3)
    SGDParams(lr, lr, N_sgd, N_mc, N_mc)
end

function sgd(;L=2,
              P=1,
              hX = 1.0,
              θ = rand(3),
              SGD=SGDParams(),
              path=false, #"../data/2d-tfim/",
              print=false)

    num_chain = nworkers()

    θ[3:3:3P] .= zeros(P) # just for hZ=0
    h = (ZZ = 1.0,
         X = hX,
         Z = 0.0)

    M = Model(L, θ, h)
    Obs = [Thermo(P) for i in 1:num_chain]
    Grad = [Gradient(Obs[i], P) for i in 1:num_chain]
    Sto = [StochasticMatrix(Obs[i], P) for i in 1:num_chain]
    C = [Config(size(M)) for i in 1:num_chain]

    E, Eσ, mag, Mσ, grad = pmc(C, M, Obs, Grad, SGD.N_mc, Sto, num_chain)

    for i in 1:SGD.N_sgd
        M = Model(L, θ, h)
        Obs = [Thermo(P) for i in 1:num_chain]
        @time E, Eσ, mag, Mσ, grad = pmc(C, M, Obs, Grad, SGD.N_mc, Sto, num_chain)

        normalize!(grad)
        θ.+= SGD.lr * grad
        θ.= abs.(θ)
        θ[3:3:3P] .= zeros(P) # just for hZ=0

        SGD.lr = SGD.lr₀ * 0.4512 ^ floor(i / 10)
        SGD.N_mc += 10*floor(1/SGD.lr) + 100

        if print == true
            println("\nIteration :", i)
            # println("θ is :", θ[θ .> 0])
            println("N is :", SGD.N_mc)
            println("Δ is :", SGD.lr)
            println("E is :", E/L^2)
            println("M is :", mag/L^2)
            println("gradient is :", grad)
        end

        if path != false
            params = θ[θ .> 0]
            result = vcat(sum(params), params, E/L^2, Eσ/L^2, mag/L^2, Mσ/L^2)'
            open(path, "a") do io
                writedlm(io, result)
            end
        end
    end

    return E
end


text:
export decision_tree

import ..DecisionTreeModel

using mlpack._Internal.io

import mlpack_jll
const decision_treeLibrary = mlpack_jll.libmlpack_julia_decision_tree

# Call the C binding of the mlpack decision_tree binding.
function decision_tree_mlpackMain()
  success = ccall((:decision_tree, decision_treeLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module decision_tree_internal
  import ..decision_treeLibrary

import ...DecisionTreeModel

# Get the value of a model pointer parameter of type DecisionTreeModel.
function IOGetParamDecisionTreeModel(paramName::String)::DecisionTreeModel
  DecisionTreeModel(ccall((:IO_GetParamDecisionTreeModelPtr, decision_treeLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type DecisionTreeModel.
function IOSetParamDecisionTreeModel(paramName::String, model::DecisionTreeModel)
  ccall((:IO_SetParamDecisionTreeModelPtr, decision_treeLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeDecisionTreeModel(stream::IO, model::DecisionTreeModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeDecisionTreeModelPtr, decision_treeLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeDecisionTreeModel(stream::IO)::DecisionTreeModel
  buffer = read(stream)
  DecisionTreeModel(ccall((:DeserializeDecisionTreeModelPtr, decision_treeLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    decision_tree(; [input_model, labels, maximum_depth, minimum_gain_split, minimum_leaf_size, print_training_accuracy, print_training_error, test, test_labels, training, verbose, weights])

Train and evaluate using a decision tree.  Given a dataset containing numeric or
categorical features, and associated labels for each point in the dataset, this
program can train a decision tree on that data.

The training set and associated labels are specified with the `training` and
`labels` parameters, respectively.  The labels should be in the range [0,
num_classes - 1]. Optionally, if `labels` is not specified, the labels are
assumed to be the last dimension of the training dataset.

When a model is trained, the `output_model` output parameter may be used to save
the trained model.  A model may be loaded for predictions with the `input_model`
parameter.  The `input_model` parameter may not be specified when the `training`
parameter is specified.  The `minimum_leaf_size` parameter specifies the minimum
number of training points that must fall into each leaf for it to be split.  The
`minimum_gain_split` parameter specifies the minimum gain that is needed for the
node to split.  The `maximum_depth` parameter specifies the maximum depth of the
tree.  If `print_training_error` is specified, the training error will be
printed.

Test data may be specified with the `test` parameter, and if performance numbers
are desired for that test set, labels may be specified with the `test_labels`
parameter.  Predictions for each test point may be saved via the `predictions`
output parameter.  Class probabilities for each prediction may be saved with the
`probabilities` output parameter.

For example, to train a decision tree with a minimum leaf size of 20 on the
dataset contained in `data` with labels `labels`, saving the output model to
`tree` and printing the training error, one could call

```julia
julia> using CSV
julia> labels = CSV.read("labels.csv"; type=Int)
julia> tree, _, _ = decision_tree(labels=labels,
            minimum_gain_split=0.001, minimum_leaf_size=20,
            print_training_accuracy=1, training=data)
```

Then, to use that model to classify points in `test_set` and print the test
error given the labels `test_labels` using that model, while saving the
predictions for each point to `predictions`, one could call

```julia
julia> using CSV
julia> test_labels = CSV.read("test_labels.csv"; type=Int)
julia> _, predictions, _ = decision_tree(input_model=tree,
            test=test_set, test_labels=test_labels)
```

# Arguments

 - `input_model::DecisionTreeModel`: Pre-trained decision tree, to be used
      with test points.
 - `labels::Array{Int, 1}`: Training labels.
 - `maximum_depth::Int`: Maximum depth of the tree (0 means no limit).
      Default value `0`.

 - `minimum_gain_split::Float64`: Minimum gain for node splitting.
      Default value `1e-07`.

 - `minimum_leaf_size::Int`: Minimum number of points in a leaf.  Default
      value `20`.

 - `print_training_accuracy::Bool`: Print the training accuracy.  Default
      value `false`.

 - `print_training_error::Bool`: Print the training error (deprecated;
      will be removed in mlpack 4.0.0).  Default value `false`.

 - `test::Tuple{Array{Bool, 1}, Array{Float64, 2}}`: Testing dataset (may
      be categorical).
 - `test_labels::Array{Int, 1}`: Test point labels, if accuracy
      calculation is desired.
 - `training::Tuple{Array{Bool, 1}, Array{Float64, 2}}`: Training dataset
      (may be categorical).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

 - `weights::Array{Float64, 2}`: The weight of labels

# Return values

 - `output_model::DecisionTreeModel`: Output for trained decision tree.
 - `predictions::Array{Int, 1}`: Class predictions for each test point.
 - `probabilities::Array{Float64, 2}`: Class probabilities for each test
      point.

"""
function decision_tree(;
                       input_model::Union{DecisionTreeModel, Missing} = missing,
                       labels = missing,
                       maximum_depth::Union{Int, Missing} = missing,
                       minimum_gain_split::Union{Float64, Missing} = missing,
                       minimum_leaf_size::Union{Int, Missing} = missing,
                       print_training_accuracy::Union{Bool, Missing} = missing,
                       print_training_error::Union{Bool, Missing} = missing,
                       test::Union{Tuple{Array{Bool, 1}, Array{Float64, 2}}, Missing} = missing,
                       test_labels = missing,
                       training::Union{Tuple{Array{Bool, 1}, Array{Float64, 2}}, Missing} = missing,
                       verbose::Union{Bool, Missing} = missing,
                       weights = missing,
                       points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, decision_treeLibrary), Nothing, ());

  IORestoreSettings("Decision tree")

  # Process each input argument before calling mlpackMain().
  if !ismissing(input_model)
    decision_tree_internal.IOSetParamDecisionTreeModel("input_model", convert(DecisionTreeModel, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(maximum_depth)
    IOSetParam("maximum_depth", convert(Int, maximum_depth))
  end
  if !ismissing(minimum_gain_split)
    IOSetParam("minimum_gain_split", convert(Float64, minimum_gain_split))
  end
  if !ismissing(minimum_leaf_size)
    IOSetParam("minimum_leaf_size", convert(Int, minimum_leaf_size))
  end
  if !ismissing(print_training_accuracy)
    IOSetParam("print_training_accuracy", convert(Bool, print_training_accuracy))
  end
  if !ismissing(print_training_error)
    IOSetParam("print_training_error", convert(Bool, print_training_error))
  end
  if !ismissing(test)
    IOSetParam("test", convert(Tuple{Array{Bool, 1}, Array{Float64, 2}}, test), points_are_rows)
  end
  if !ismissing(test_labels)
    IOSetParamURow("test_labels", test_labels)
  end
  if !ismissing(training)
    IOSetParam("training", convert(Tuple{Array{Bool, 1}, Array{Float64, 2}}, training), points_are_rows)
  end
  if !ismissing(weights)
    IOSetParamMat("weights", weights, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  decision_tree_mlpackMain()

  return decision_tree_internal.IOGetParamDecisionTreeModel("output_model"),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
@testset "refine_root" begin
    # Define a number of problems for which we expect the method to
    # work without issues. The problems are defined as (f, enclosure,
    # answer, atol, rtol).
    problems = [
        # Problems with the root away from zero so that the relative
        # tolerance is relevant
        [
            (sin, Arblib.add_error!(root, Mag(1)), root, 0, sqrt(eps(Arb))) for
            root in Arb(π) .* [-10:-1; 1:10]
        ]...,
        (
            x -> exp(x - 1) - 1,
            Arblib.add_error!(Arb(1.5), Mag(1)),
            Arb(1),
            0,
            sqrt(eps(Arb)),
        ),
        (
            x -> exp(x + 1) - 1,
            Arblib.add_error!(Arb(-1.5), Mag(1)),
            Arb(-1),
            0,
            sqrt(eps(Arb)),
        ),

        # Problems where the root is at zero so that the absolute
        # tolerance is relevant
        (sin, Arblib.add_error!(Arb(0.1), Mag(1)), Arb(0), sqrt(eps(Arb)), 0),
        (x -> exp(x) - 1, Arblib.add_error!(Arb(0.1), Mag(1)), Arb(0), sqrt(eps(Arb)), 0),

        # ArbPoly
        (
            Arblib.fromroots(ArbPoly, [1, 2, 3]),
            Arblib.add_error!(Arb(2.001), Mag(0.01)),
            Arb(2),
            0,
            sqrt(eps(Arb)),
        ),
        (
            Arblib.fromroots(ArbPoly, [-Arb(π), Arb(π)]),
            Arblib.add_error!(Arb(π), Mag(0.1)),
            Arb(π),
            0,
            sqrt(eps(Arb)),
        ),
    ]

    for (f, enclosure, answer, atol, rtol) in problems
        root = ArbExtras.refine_root(f, enclosure; atol, rtol)
        aerror = 2Arblib.radius(Arb, root)
        rerror = aerror / root
        @test Arblib.overlaps(root, answer)
        @test aerror < atol || rerror < rtol
    end

    # Test that it returns NaN in case it can't prove we have a root unless strict is true
    # Double root
    @test isnan(ArbExtras.refine_root(x -> cos(x) - 1, Arblib.add_error!(Arb(0), Mag(1))))
    @test isequal(
        ArbExtras.refine_root(
            x -> cos(x) - 1,
            Arblib.add_error!(Arb(0), Mag(1)),
            strict = false,
        ),
        Arblib.add_error!(Arb(0), Mag(1)),
    )

    # Test that it stops correctly if the enclosure doesn't improve
    # and that it runs the minimum number of iterations still. We do
    # this by implementing a version of sin that counts the number of
    # times it's called.
    count = 0
    g = x -> begin
        count += 1
        return sin(x)
    end
    ArbExtras.refine_root(g, Arb(π), rtol = 0, min_iterations = 5)
    @test count == 10

    # Root on endpoint of interval
    # Enclosure which has 1 as an endpoint
    enclosure = Arblib.add_error!(Arb(1.125), Arblib.set_ui_2exp!(Mag(), unsigned(1), -3))
    @test isnan(ArbExtras.refine_root(sinpi, enclosure))
    root = ArbExtras.refine_root(sinpi, enclosure, strict = false)
    @test Arblib.overlaps(root, enclosure)
    @test contains(root, 1)

    # Test that it doesn't find a root if the root is just outside the
    # enclosure
    @test !contains(enclosure + Arb(1e-40), 1)
    @test isnan(ArbExtras.refine_root(sinpi, enclosure + Arb(1e-40)))

    # Check max_iterations

    # The function x -> x - 1 but with a weird implementation such
    # that for Arb it returns union(x, 1) - 1. Note that this still
    # gives a correct enclosure. This is done so that the convergence
    # for the Newton iterations is very slow.
    f(x::Arb) = begin
        return union(Arb(1), x) - 1
    end
    f(x) = x - 1

    @test Arblib.radius(ArbExtras.refine_root(f, Arblib.add_error!(Arb(1.1), Mag(1)))) <
          Arblib.radius(
        ArbExtras.refine_root(f, Arblib.add_error!(Arb(1.1), Mag(1)), atol = 1e-1),
    )
    @test Arblib.radius(
        ArbExtras.refine_root(f, Arblib.add_error!(Arb(1.1), Mag(1)), max_iterations = 40),
    ) < Arblib.radius(
        ArbExtras.refine_root(f, Arblib.add_error!(Arb(1.1), Mag(1)), max_iterations = 20),
    )
end


text:
using PkgBenchmark
using Trixi

let results = judge(Trixi,
            BenchmarkConfig(juliacmd=`$(Base.julia_cmd()) --check-bounds=no --threads=1`), # target
            BenchmarkConfig(juliacmd=`$(Base.julia_cmd()) --check-bounds=no --threads=1`, id="main") # baseline
        )

    export_markdown(joinpath(pathof(Trixi) |> dirname |> dirname, "benchmark", "results_threads1.md"), results)
end

let results = judge(Trixi,
            BenchmarkConfig(juliacmd=`$(Base.julia_cmd()) --check-bounds=no --threads=2`), # target
            BenchmarkConfig(juliacmd=`$(Base.julia_cmd()) --check-bounds=no --threads=2`, id="main") # baseline
        )

    export_markdown(joinpath(pathof(Trixi) |> dirname |> dirname, "benchmark", "results_threads2.md"), results)
end


text:
struct Transpose{T,N} <: ArrayOp{T,N}
    f::Function
    input::ArrayOp
end

function Transpose(f,x::ArrayOp)
    @assert 1 <= ndims(x) && ndims(x) <= 2
    Transpose{eltype(x), 2}(f,x)
end
function size(x::Transpose)
    sz = size(x.input)
    if length(sz) == 1
        (1, sz[1])
    else
        (sz[2], sz[1])
    end
end

transpose(x::ArrayOp) = Transpose(transpose, x)
transpose(x::Union{Chunk, Thunk}) = Thunk(transpose, x)

adjoint(x::ArrayOp) = Transpose(adjoint, x)
adjoint(x::Union{Chunk, Thunk}) = Thunk(adjoint, x)

function adjoint(x::ArrayDomain{2})
    d = indexes(x)
    ArrayDomain(d[2], d[1])
end
function adjoint(x::ArrayDomain{1})
    d = indexes(x)
    ArrayDomain(1, d[1])
end

function _ctranspose(x::AbstractArray)
    Any[delayed(adjoint)(x[j,i]) for i=1:size(x,2), j=1:size(x,1)]
end

function stage(ctx::Context, node::Transpose)
    inp = cached_stage(ctx, node.input)
    thunks = _ctranspose(chunks(inp))
    DArray(eltype(inp), domain(inp)', domainchunks(inp)', thunks)
end

import Base: *, +

struct MatMul{T, N} <: ArrayOp{T, N}
    a::ArrayOp
    b::ArrayOp
end

function mul_size(a,b)
  if ndims(b) == 1
    (size(a,1),)
  else
    (size(a,1), size(b,2))
  end
end
size(x::MatMul) = mul_size(x.a, x.b)
MatMul(a,b) =
  MatMul{promote_type(eltype(a), eltype(b)), length(mul_size(a,b))}(a,b)
(*)(a::ArrayOp, b::ArrayOp) = MatMul(a,b)
# Bonus method for matrix-vector multiplication
(*)(a::ArrayOp, b::Vector) = MatMul(a,PromotePartition(b))
(*)(a::AbstractArray, b::ArrayOp) = MatMul(PromotePartition(a), b)

function (*)(a::ArrayDomain{2}, b::ArrayDomain{2})

    if size(a, 2) != size(b, 1)
        throw(DimensionMismatch("The domains cannot be multiplied"))
    end

    ArrayDomain((indexes(a)[1], indexes(b)[2]))
end
function (*)(a::ArrayDomain{2}, b::ArrayDomain{1})
    if size(a, 2) != length(b)
        throw(DimensionMismatch("The domains cannot be multiplied"))
    end
    ArrayDomain((indexes(a)[1],))
end

function (*)(a::Blocks{2}, b::Blocks{2})
    Blocks(a.blocksize[1], b.blocksize[2])
end
(*)(a::Blocks{2}, b::Blocks{1}) =
    Blocks((a.blocksize[1],))

function (+)(a::ArrayDomain, b::ArrayDomain)
    if a == b
        DimensionMismatch("The domains cannot be added")
    end
    a
end

(*)(a::Union{Chunk, Thunk}, b::Union{Chunk, Thunk}) = Thunk(*, a,b)
(+)(a::Union{Chunk, Thunk}, b::Union{Chunk, Thunk}) = Thunk(+, a,b)

# we define our own matmat and matvec multiply
# for computing the new domains and thunks.
function _mul(a::Matrix, b::Matrix; T=eltype(a))
    c = Array{T}(undef, (size(a,1), size(b,2)))
    n = size(a, 2)
    for i=1:size(a,1)
        for j=1:size(b, 2)
            c[i,j] = treereduce(+, map(*, reshape(a[i,:], (n,)), b[:, j]))
        end
    end
    c
end

function _mul(a::Matrix, b::Vector; T=eltype(b))
    c = Array{T}(undef, size(a,1))
    n = size(a,2)
    for i=1:size(a,1)
        c[i] = treereduce(+, map(*, reshape(a[i, :], (n,)), b))
    end
    c
end

function _mul(a::Vector, b::Vector; T=eltype(b))
    @assert length(b) == 1
    [x * b[1] for x in a]
end

function promote_distribution(ctx::Context, m::MatMul, a,b)
    iscompat = try domain(a) * domain(b); true
               catch e; false end
    if iscompat
        return a,b
    end

    pa = domainchunks(a)
    pb = domainchunks(b)

    d = DomainBlocks((1,1), (pa.cumlength[2], pb.cumlength[2])) # FIXME: this is not generic
    a, cached_stage(ctx, Distribute(d, b))
end

function stage_operands(ctx::Context, m::MatMul, a, b)
    if size(a, 2) != size(b, 1)
        error(DimensionMismatch("Inputs to * have incompatible size"))
    end
    # take the row distribution of a and get b onto that.

    stg_a = cached_stage(ctx, a)
    stg_b = cached_stage(ctx, b)
    promote_distribution(ctx, m, stg_a, stg_b)
end

"An operand which should be distributed as per convenience"
function stage_operands(ctx::Context, ::MatMul, a::ArrayOp, b::PromotePartition{T,1}) where T
    stg_a = cached_stage(ctx, a)
    dmn_a = domain(stg_a)
    dchunks_a = domainchunks(stg_a)
    dmn_b = domain(b.data)
    if size(dmn_a, 2) != size(dmn_b, 1)
        throw(DimensionMismatch("Cannot promote array of domain $(dmn_b) to multiply with an array of size $(dmn_a)"))
    end
    dmn_out = DomainBlocks((1,),(dchunks_a.cumlength[2],))

    stg_a, cached_stage(ctx, Distribute(dmn_out, b.data))
end

function stage_operands(ctx::Context, ::MatMul, a::PromotePartition, b::ArrayOp)

    if size(a, 2) != size(b, 1)
        throw(DimensionMismatch("Cannot promote array of domain $(dmn_b) to multiply with an array of size $(dmn_a)"))
    end
    stg_b = cached_stage(ctx, b)

    ps = domainchunks(stg_b)
    dmn_out = DomainBlocks((1,1),([size(a.data, 1)], ps.cumlength[1],))
    cached_stage(ctx, Distribute(dmn_out, a.data)), stg_b
end

function stage(ctx::Context, mul::MatMul)
    a, b = stage_operands(ctx, mul, mul.a, mul.b)
    d = domain(a)*domain(b)
    DArray(Any, d, domainchunks(a)*domainchunks(b),
                          _mul(chunks(a), chunks(b); T=Thunk))
end

### Scale

struct Scale{T,N} <: ArrayOp{T,N}
    l::ArrayOp
    r::ArrayOp
end
Scale(l::ArrayOp{Tl}, r::ArrayOp{Tr,N}) where {Tl, Tr, N} =
  Scale{promote_type(Tl, Tr), N}(l,r)

size(s::Scale) = size(s.l)

scale(l::Number, r::ArrayOp) = BlockwiseOp(x->scale(l, x), (r,))
scale(l::Vector, r::ArrayOp) = scale(PromotePartition(l), r)
(*)(l::Diagonal, r::ArrayOp) = Scale(PromotePartition(l.diag), r)
scale(l::ArrayOp, r::ArrayOp) = Scale(l, r)

function stage_operand(ctx::Context, ::Scale, a, b::PromotePartition)
    ps = domainchunks(a)
    b_parts = DomainBlocks((1,), (ps.cumlength[1],))
    cached_stage(ctx, Distribute(b_parts, b.data))
end

function stage_operand(ctx::Context, ::Scale, a, b)
    cached_stage(ctx, b)
end

function _scale(l, r)
    res = similar(r, Any)
    for i=1:length(l)
        res[i,:] = map(x->Thunk((a,b) -> Diagonal(a)*b, l[i], x), r[i,:])
    end
    res
end

function stage(ctx::Context, scal::Scale)
    r = cached_stage(ctx, scal.r)
    l = stage_operand(ctx, scal, r, scal.l)

    @assert size(domain(r), 1) == size(domain(l), 1)

    scal_parts = _scale(chunks(l), chunks(r))
    DArray(Any, domain(r), domainchunks(r), scal_parts)
end

struct Concat{T,N} <: ArrayOp{T,N}
    axis::Int
    inputs::Tuple
end
Concat(axis::Int, inputs::Tuple) =
  Concat{promote_type(map(eltype, inputs)...),
         ndims(inputs[1])}(axis, inputs)

function size(c::Concat)
    sz = [size(c.inputs[1])...]
    sz[c.axis] = sum(map(x->size(x, c.axis), c.inputs))
    (sz...,)
end

function Base.cat(d::ArrayDomain, ds::ArrayDomain...; dims::Int)
    h = (d)
    out_idxs = [x for x in indexes(h)]
    len = sum(map(x->length(indexes(x)[dims]), (d, ds...)))
    fst = first(out_idxs[dims])
    out_idxs[dims] = fst:(fst+len-1)
    ArrayDomain(out_idxs)
end

function stage(ctx::Context, c::Concat)
    inp = Any[cached_stage(ctx, x) for x in c.inputs]

    dmns = map(domain, inp)
    dims = [[i == c.axis ? 0 : i for i in size(d)] for d in dmns]
    if !all(map(x -> x == dims[1], dims[2:end]))
        error("Inputs to cat do not have compatible dimensions.")
    end

    dmn = cat(dmns..., dims = c.axis)
    dmnchunks = cumulative_domains(cat(map(domainchunks, inp)..., dims = c.axis))
    thunks = cat(map(chunks, inp)..., dims = c.axis)
    T = promote_type(map(eltype, inp)...)
    DArray(T, dmn, dmnchunks, thunks)
end

Base.cat(x::ArrayOp, xs::ArrayOp...; dims::Int) = Concat(dims, (x, xs...))

Base.hcat(xs::ArrayOp...) = cat(xs..., dims=2)
Base.vcat(xs::ArrayOp...) = cat(xs..., dims=1)


text:
url = "http://vincentarelbundock.github.io/Rdatasets/csv/quantreg/engel.csv"
Data = CSV.read(download(url))

out_ip = qreg(@formula(foodexp~income), Data, IP())
out_irls = qreg(@formula(foodexp~income), Data, IRLS())

@test norm( out_ip.model.beta - [81.4822; 0.560181]) < 1e-4
@test norm( out_irls.model.beta - [81.4823; 0.560181]) < 1e-4

#=
...... Interior point ......
TableRegressionModel{QRegModel,Float64}:

Coefficients:
             Estimate Std.Error t value
(Intercept)   81.4822   14.6345 5.56783
income       0.560181 0.0131756 42.5164

.......... IRLS ............
TableRegressionModel{QRegModel,Float64}:

Coefficients:
             Estimate Std.Error t value
(Intercept)   81.4823   14.6345 5.56783
income       0.560181 0.0131756 42.5164
=#