text:
using CUDA: CuArray, CuDevice, CuContext, CuPtr, device, device!, synchronize
using OffsetArrays
using Oceananigans.Grids: AbstractGrid
import Base: length

const GPUVar = Union{CuArray, CuContext, CuPtr, Ptr}

###
### Multi Region Object
###

struct MultiRegionObject{R, D}
    regions :: R
    devices :: D
end

###
### Convenience structs
###

struct Reference{R}
    ref :: R
end

struct Iterate{I}
    iter :: I
end

###
### Multi region functions
###

@inline getdevice(a, i)                     = nothing
@inline getdevice(cu::GPUVar, i)            = CUDA.device(cu)
@inline getdevice(cu::OffsetArray, i)       = getdevice(cu.parent)
@inline getdevice(mo::MultiRegionObject, i) = mo.devices[i]

@inline getdevice(a)               = nothing
@inline getdevice(cu::GPUVar)      = CUDA.device(cu)
@inline getdevice(cu::OffsetArray) = getdevice(cu.parent)

@inline switch_device!(a)                        = nothing
@inline switch_device!(dev::Int)                 = CUDA.device!(dev)
@inline switch_device!(dev::CuDevice)            = CUDA.device!(dev)
@inline switch_device!(dev::Tuple, i)            = switch_device!(dev[i])
@inline switch_device!(mo::MultiRegionObject, i) = switch_device!(getdevice(mo, i))

@inline getregion(a, i) = a
@inline getregion(ref::Reference, i)        = ref.ref
@inline getregion(iter::Iterate, i)         = iter.iter[i]
@inline getregion(mo::MultiRegionObject, i) = _getregion(mo.regions[i], i)
@inline getregion(p::Pair, i)               = p.first => _getregion(p.second, i)

@inline _getregion(a, i) = a
@inline _getregion(ref::Reference, i)        = ref.ref
@inline _getregion(iter::Iterate, i)         = iter.iter[i]
@inline _getregion(mo::MultiRegionObject, i) = getregion(mo.regions[i], i)
@inline _getregion(p::Pair, i)               = p.first => getregion(p.second, i)

## The implementation of `getregion` for a Tuple forces the compiler to infer the size of the Tuple
@inline getregion(t::Tuple{}, i)                            = ()
@inline getregion(t::Tuple{<:Any}, i)                       = (_getregion(t[1], i), )
@inline getregion(t::Tuple{<:Any, <:Any}, i)                = (_getregion(t[1], i), _getregion(t[2], i))
@inline getregion(t::Tuple{<:Any, <:Any, <:Any}, i)         = (_getregion(t[1], i), _getregion(t[2], i), _getregion(t[3], i))
@inline getregion(t::Tuple, i)                              = (_getregion(t[1], i), _getregion(t[2:end], i)...)

@inline _getregion(t::Tuple{}, i)                           = ()
@inline _getregion(t::Tuple{<:Any}, i)                      = (getregion(t[1], i), )
@inline _getregion(t::Tuple{<:Any, <:Any}, i)               = (getregion(t[1], i), getregion(t[2], i))
@inline _getregion(t::Tuple{<:Any, <:Any, <:Any}, i)        = (getregion(t[1], i), getregion(t[2], i), getregion(t[3], i))
@inline _getregion(t::Tuple, i)                             = (getregion(t[1], i), getregion(t[2:end], i)...)

@inline getregion(nt::NamedTuple, i)   = NamedTuple{keys(nt)}(_getregion(Tuple(nt), i))
@inline _getregion(nt::NamedTuple, i)  = NamedTuple{keys(nt)}(getregion(Tuple(nt), i))

@inline isregional(a)                   = false
@inline isregional(::MultiRegionObject) = true

@inline isregional(t::Tuple{}) = false
@inline isregional(nt::NT) where NT<:NamedTuple{<:Any, Tuple{Tuple{}}} = false
for func in [:isregional, :devices, :switch_device!]
    @eval begin
        @inline $func(t::Union{Tuple, NamedTuple}) = $func(first(t))
    end
end

@inline devices(mo::MultiRegionObject) = mo.devices

Base.getindex(mo::MultiRegionObject, i, args...) = Base.getindex(mo.regions, i, args...)
Base.length(mo::MultiRegionObject)               = Base.length(mo.regions)

# For non-returning functions -> can we make it NON BLOCKING? This seems to be synchronous!
@inline function apply_regionally!(func!, args...; kwargs...)
    mra = isnothing(findfirst(isregional, args)) ? nothing : args[findfirst(isregional, args)]
    mrk = isnothing(findfirst(isregional, kwargs)) ? nothing : kwargs[findfirst(isregional, kwargs)]
    isnothing(mra) && isnothing(mrk) && return func!(args...; kwargs...)

    if isnothing(mra)
        devs = devices(mrk)
    else
        devs = devices(mra)
    end

    for (r, dev) in enumerate(devs)
        switch_device!(dev)
        func!((getregion(arg, r) for arg in args)...; (getregion(kwarg, r) for kwarg in kwargs)...)
    end

    sync_all_devices!(devs)
end

# For functions with return statements -> BLOCKING! (use as seldom as possible)
@inline function construct_regionally(constructor, args...; kwargs...)
    mra = isnothing(findfirst(isregional, args)) ? nothing : args[findfirst(isregional, args)]
    mrk = isnothing(findfirst(isregional, kwargs)) ? nothing : kwargs[findfirst(isregional, kwargs)]
    isnothing(mra) && isnothing(mrk) && return constructor(args...; kwargs...)

    if isnothing(mra)
        devs = devices(mrk)
    else
        devs = devices(mra)
    end

    res = Vector(undef, length(devs))
    for (r, dev) in enumerate(devs)
        switch_device!(dev)
        res[r] = constructor((getregion(arg, r) for arg in args)...; (getregion(kwarg, r) for kwarg in kwargs)...)
    end
    sync_all_devices!(devs)

    return MultiRegionObject(Tuple(res), devs)
end

@inline sync_all_devices!(grid::AbstractGrid)    = nothing
@inline sync_all_devices!(mo::MultiRegionObject) = sync_all_devices!(devices(mo))

@inline function sync_all_devices!(devices)
    for dev in devices
        switch_device!(dev)
        sync_device!(dev)
    end
end

@inline sync_device!(::CuDevice) = CUDA.device_synchronize()
@inline sync_device!(dev)        = nothing

"""
macro `@apply_regionally` to distribute locally the function calls

calls `compute_regionally` in case of a returning value and `apply_regionally!`
in case of no return
"""
macro apply_regionally(expr)
    if expr.head == :call
        func = expr.args[1]
        args = expr.args[2:end]
        multi_region = quote
            apply_regionally!($func, $(args...))
        end
        return quote
            $(esc(multi_region))
        end
    elseif expr.head == :(=)
        ret = expr.args[1]
        exp = expr.args[2]
        func = exp.args[1]
        args = exp.args[2:end]
        multi_region = quote
            $ret = construct_regionally($func, $(args...))
        end
        return quote
            $(esc(multi_region))
        end
    elseif expr.head == :block
        new_expr = deepcopy(expr)
        for (idx, arg) in enumerate(expr.args)
            if arg isa Expr && arg.head == :call
                func = arg.args[1]
                args = arg.args[2:end]
                new_expr.args[idx] = quote
                    apply_regionally!($func, $(args...))
                end
            elseif arg isa Expr && arg.head == :(=)
                ret = arg.args[1]
                exp = arg.args[2]
                func = exp.args[1]
                args = exp.args[2:end]
                new_expr.args[idx] = quote
                    $ret = construct_regionally!($func, $(args...))
                end
            end
        end
        return quote
            $(esc(new_expr))
        end
    end
end


text:
function test_modeterms(scenario)
    @unpack tx, rx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)

    ea = TEST_MODES[scenario][1]
    modeequation = PhysicalModeEquation(ea, tx.frequency, waveguide)

    groundsampler = GroundSampler(rx.distance, rx.fieldcomponent)
    sampler = Sampler(rx.distance, rx.fieldcomponent, 0.0)

    # each field
    for fc in (Fields.Ez, Fields.Ey, Fields.Ex)
        groundsampler = GroundSampler(rx.distance, fc)
        sampler = Sampler(rx.distance, fc, 0.0)

        tx1, rx1 = LMP.modeterms(modeequation, tx, sampler)
        tx2, rx2 = LMP.modeterms(modeequation, tx, groundsampler)  # specialized

        @test tx1 ≈ tx2
        @test rx1 ≈ rx2
    end

    # frequency mismatch with modeequation
    txwrong = Transmitter(15e3)
    @test_throws ArgumentError LMP.modeterms(modeequation, txwrong, sampler)
end

function test_Efield(scenario)
    @unpack tx, rx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)

    modes = TEST_MODES[scenario]

    X = LMP.distance(rx, tx)

    E1 = LMP.Efield(modes, waveguide, tx, rx)  # out-of-place

    singlerx = GroundSampler(1000e3, rx.fieldcomponent)
    E2 = LMP.Efield(modes, waveguide, tx, singlerx)  # specialized

    distidx = findfirst(x->x==1000e3, X)
    @test E2 ≈ E1[distidx]
end

@testset "modesum.jl" begin
    @info "Testing modesum"

    for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
        multiplespecies_scenario)

        test_modeterms(scn)
        test_Efield(scn)
    end
end

ext: jl


text:
using MyPkg
using Test
using Dates

@testset "MyPkg.jl" begin
    @test MyPkg.todayis() == dayname(now())
end

max_stars_repo_name: glennmoy/MyPkg

max_issues_repo_name: glennmoy/MyPkg

max_forks_repo_name: glennmoy/MyPkg


text:
##### generic methods (fallback) #####

## sampling

rand(d::UnivariateDistribution) = quantile(d, rand())

rand!(d::UnivariateDistribution, A::AbstractArray) = _rand!(sampler(d), A)
rand(d::UnivariateDistribution, n::Int) = _rand!(sampler(d), Array(eltype(d), n))
rand(d::UnivariateDistribution, shp::Dims) = _rand!(sampler(d), Array(eltype(d), shp))

## domain

isbounded(d::UnivariateDistribution) = isupperbounded(d) && islowerbounded(d)
hasfinitesupport(d::DiscreteUnivariateDistribution) = isbounded(d)
hasfinitesupport(d::ContinuousUnivariateDistribution) = false

insupport(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) && (minimum(d) <= x <= maximum(d))
insupport(d::ContinuousUnivariateDistribution, x::Real) = minimum(d) <= x <= maximum(d)

function insupport!{D<:UnivariateDistribution}(r::AbstractArray, d::Union(D,Type{D}), X::AbstractArray)
    length(r) == length(X) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    for i in 1 : length(X)
        @inbounds r[i] = insupport(d, X[i])
    end
    return r
end

insupport{D<:UnivariateDistribution}(d::Union(D,Type{D}), X::AbstractArray) =
     insupport!(BitArray(size(X)), d, X)

## statistics

std(d::UnivariateDistribution) = sqrt(var(d))
median(d::UnivariateDistribution) = quantile(d, 0.5)
modes(d::UnivariateDistribution) = [mode(d)]
entropy(d::UnivariateDistribution, b::Real) = entropy(d) / log(b)

isplatykurtic(d::UnivariateDistribution) = kurtosis(d) > 0.0
isleptokurtic(d::UnivariateDistribution) = kurtosis(d) < 0.0
ismesokurtic(d::UnivariateDistribution) = kurtosis(d) == 0.0

function kurtosis(d::Distribution, correction::Bool)
    if correction
        return kurtosis(d)
    else
        return kurtosis(d) + 3.0
    end
end

excess(d::Distribution) = kurtosis(d)
excess_kurtosis(d::Distribution) = kurtosis(d)
proper_kurtosis(d::Distribution) = kurtosis(d, false)

## pdf, cdf, and friends

logpdf(d::UnivariateDistribution, x::Number) = log(pdf(d, x))
cdf(d::DiscreteUnivariateDistribution, k::Real) = sum([pdf(d,i) for i in minimum(d):k])
ccdf(d::UnivariateDistribution, q::Real) = 1.0 - cdf(d, q)
cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p)

logcdf(d::UnivariateDistribution, q::Real) = log(cdf(d,q))
logccdf(d::UnivariateDistribution, q::Real) = log(ccdf(d,q))
invlogccdf(d::UnivariateDistribution, lp::Real) = quantile(d, -expm1(lp))
invlogcdf(d::UnivariateDistribution, lp::Real) = quantile(d, exp(lp))

# vectorized versions
for fun in [:pdf, :logpdf,
            :cdf, :logcdf,
            :ccdf, :logccdf,
            :invlogcdf, :invlogccdf,
            :quantile, :cquantile]

    _fun! = symbol(string('_', fun, '!'))
    fun! = symbol(string(fun, '!'))

    @eval begin
        function ($_fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray)
            for i in 1 : length(X)
                r[i] = ($fun)(d, X[i])
            end
            return r
        end

        function ($fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray)
            length(r) == length(X) ||
                throw(ArgumentError("Inconsistent array dimensions."))
            $(_fun!)(r, d, X)
        end

        ($fun)(d::UnivariateDistribution, X::AbstractArray) =
            $(_fun!)(Array(Float64, size(X)), d, X)
    end
end

## loglikelihood

function _loglikelihood(d::UnivariateDistribution, X::AbstractArray)
    ll = 0.0
    for i in 1:length(X)
        @inbounds ll += logpdf(d, X[i])
    end
    return ll
end

loglikelihood(d::UnivariateDistribution, X::AbstractArray) =
    _loglikelihood(d, X)

##### specific distributions #####

for finame in ["arcsine.jl",
               "bernoulli.jl",
               "beta.jl",
               "betaprime.jl",
               "binomial.jl",
               "categorical.jl",
               "cauchy.jl",
               "chi.jl",
               "chisq.jl",
               "cosine.jl",
               "discreteuniform.jl",
               "empirical.jl",
               "exponential.jl",
               "fdist.jl",
               "frechet.jl",
               "gamma.jl",
               "edgeworth.jl",
               "erlang.jl",
               "geometric.jl",
               "gumbel.jl",
               "hypergeometric.jl",
               "inversegamma.jl",
               "inversegaussian.jl",
               "kolmogorov.jl",
               "ksdist.jl",
               "ksonesided.jl",
               "laplace.jl",
               "levy.jl",
               "logistic.jl",
               "lognormal.jl",
               "negativebinomial.jl",
               "noncentralbeta.jl",
               "noncentralchisq.jl",
               "noncentralf.jl",
               "noncentralhypergeometric.jl",
               "noncentralt.jl",
               "normal.jl",
               "normalcanon.jl",
               "pareto.jl",
               "poisson.jl",
               "rayleigh.jl",
               "skellam.jl",
               "symtriangular.jl",
               "tdist.jl",
               "triangular.jl",
               "uniform.jl",
               "vonmises.jl",
               "weibull.jl"
               ]
    include(joinpath("univariate", finame))
end


text:
@with_kw struct FastGrid
    resolution::Float64 = 1.0
    #tight::Bool         = false
end

# This is the main function
function solve(solver::FastGrid, problem::Problem) #original
    result = true
    delta = solver.resolution

    center = problem.input.center
    radius = problem.input.radius[1]

    (W, b) = (problem.network.layers[1].weights, problem.network.layers[1].bias)

    np = pyimport("numpy")
    q, e, p = np.linalg.svd(W)

    input = forward_affine_map(solver, W, b, problem.input)

    lower, upper = low(input), high(input)
    n_hypers_per_dim = BigInt.(max.(ceil.(Int, (upper-lower) / delta), 1))

    k_1 = size(W, 1)
    k_0 = size(W, 2)

    hps = Array{Hyperplane, 1}(undef, k_0-k_1)

    for i in k_1+1:k_0
        hps[i-k_1] = Hyperplane(p[i,:], dot(p[i,:], center))
    end

    S = hps[1]
    if k_0-k_1 > 1
        for i in 2:k_0-k_1
            S = intersection(S, hps[i])
        end
    end

    C = vcat(W, -W)
    d = zeros(2k_1)

    #kb = p[k_1+1:k_0,:] #kernel basis
    #kc = kb * center

    #= Al = zeros(k_0, k_0)
    Au = zeros(k_0, k_0)

    dl = zeros(k_1)
    du = zeros(k_1)

    bl = zeros(k_0)
    bu = zeros(k_0)

    pl = zeros(k_0)
    pu = zeros(k_0) =#

    #count3 = BigInt(0)
    count4 = BigInt(0)
    #println("All: " * string(prod(n_hypers_per_dim)) * " - " * string(prod(n_hypers_per_dim.-2)))

    # preallocate work arrays
    local_lower, local_upper, CI = similar(lower), similar(lower), similar(lower)
    for i in 1:prod(n_hypers_per_dim)
        n = i
        for j in firstindex(CI):lastindex(CI)
            n, CI[j] = fldmod1(n, n_hypers_per_dim[j])
        end
        @. local_lower = lower + delta * (CI - 1)
        @. local_upper = min(local_lower + delta, upper)
        hyper = Hyperrectangle(low = local_lower, high = local_upper)

        d = vcat(local_upper - b, b - local_lower)

        P_i = HPolyhedron(C, d)

        inter = intersection(problem.input, P_i)

        if isempty(inter) == false
            #count3 += 1
            O_i = box_approximation(intersection(P_i, S))
            inner = true
            for j in 1:k_0
                if low(O_i)[j] ≤ low(problem.input)[j]
                    inner = false
                    break
                end

                if high(problem.input)[j] ≤ high(O_i)[j]
                    inner = false
                    break
                end
            end

            if inner == false
                reach = forward_network(solver, problem.network, hyper)
                count4 += 1
                if !issubset(reach, problem.output)
                    result = false
                end
            end
        end
    end

    println("Verified: " * string(count4))

    if result
        return BasicResult(:holds)
    end
    return BasicResult(:violated)
end

function forward_network(solver::FastGrid, nnet::Network, input::Hyperrectangle)
    layers = nnet.layers
    act = layers[1].activation
    reach = Hyperrectangle(low = act.(low(input)), high = act.(high(input)))

    for i in 2:length(layers)
        reach = forward_layer(solver, layers[i], reach)
    end
    return reach
end

# This function is called by forward_network
function forward_layer(solver::FastGrid, L::Layer, input::Hyperrectangle)
    (W, b, act) = (L.weights, L.bias, L.activation)
    center = zeros(size(W, 1))
    gamma  = zeros(size(W, 1))
    for j in 1:size(W, 1)
        node = Node(W[j,:], b[j], act)
        center[j], gamma[j] = forward_node(solver, node, input)
    end
    return Hyperrectangle(center, gamma)
end

function forward_affine_map(solver::FastGrid, W::Matrix, b::Vector, input::Hyperrectangle)
    center = W * input.center + b
    radius = abs.(W) * input.radius
    return Hyperrectangle(center, radius)
end

function forward_node(solver::FastGrid, node::Node, input::Hyperrectangle)
    output    = node.w' * input.center + node.b
    deviation = sum(abs.(node.w) .* input.radius)
    β    = node.act(output)  # TODO expert suggestion for variable name. beta? β? O? x?
    βmax = node.act(output + deviation)
    βmin = node.act(output - deviation)
    return ((βmax + βmin)/2, (βmax - βmin)/2)
    #if solver.tight
        #return ((βmax + βmin)/2, (βmax - βmin)/2)
    #else
        #return (β, max(abs(βmax - β), abs(βmin - β)))
    #end
end

#distance of point to hyperPlane cx = d
function distance(point::Vector, c::Vector, d::Real)
    if length(point) == length(c)
        return abs(sum(point .* c) - d)/sqrt(sum(abs2.(c)))
    else
        error("Dimesion dismatch for point and constraint")
    end
end


text:
using Printf
using Glob
using JLD2
using PyCall

using Oceananigans.Utils: prettytime

np = pyimport("numpy")
ma = pyimport("numpy.ma")
plt = pyimport("matplotlib.pyplot")
mticker = pyimport("matplotlib.ticker")
ccrs = pyimport("cartopy.crs")
cmocean = pyimport("cmocean")

function plot_cubed_sphere_tracer_field!(fig, ax, var, grid; add_colorbar, transform, cmap, vmin, vmax)

    Nf = grid["faces"] |> keys |> length
    Nx = grid["faces/1/grid/Nx"]
    Ny = grid["faces/1/grid/Ny"]
    Hx = grid["faces/1/grid/Hx"]
    Hy = grid["faces/1/grid/Hy"]

    for face in 1:Nf
        λᶠᶠᵃ = grid["faces/$face/grid/λᶠᶠᵃ"][1+Hx:Nx+Hx+1, 1+Hy:Ny+Hy+1]
        φᶠᶠᵃ = grid["faces/$face/grid/φᶠᶠᵃ"][1+Hx:Nx+Hx+1, 1+Hy:Ny+Hy+1]

        var_face = var[face][:, :, 1]

        # Remove very specific problematic grid cells near λ = 180° on face 6 that mess up the plot.
        # May be related to https://github.com/SciTools/cartopy/issues/1151
        # Not needed for Orthographic or NearsidePerspective projection so let's use those.
        # if face == 6
        #     for i in 1:Nx+1, j in 1:Ny+1
        #         if isapprox(λᶠᶠᵃ[i, j], -180, atol=15)
        #             var_face[min(i, Nx), min(j, Ny)] = NaN
        #         end
        #     end
        # end

        var_face_masked = ma.masked_where(np.isnan(var_face), var_face)

        im = ax.pcolormesh(λᶠᶠᵃ, φᶠᶠᵃ, var_face_masked; transform, cmap, vmin, vmax)

        # Add colorbar below all the subplots.
        if add_colorbar && face == Nf
            ax_cbar = fig.add_axes([0.25, 0.1, 0.5, 0.02])
            fig.colorbar(im, cax=ax_cbar, orientation="horizontal")
        end

        ax.set_global()
    end

    return ax
end

function animate_rossby_haurwitz(; projections=[ccrs.Robinson()])

    ## Extract data

    file = jldopen("cubed_sphere_rossby_haurwitz.jld2")

    iterations = parse.(Int, keys(file["timeseries/t"]))

    ## Makie movie of tracer field h

    for (n, i) in enumerate(iterations)
        @info "Plotting iteration $i/$(iterations[end]) (frame $n/$(length(iterations)))..."

        u = file["timeseries/u/$i"]
        v = file["timeseries/v/$i"]
        η = file["timeseries/η/$i"]

        t = prettytime(file["timeseries/t/$i"])
        plot_title = "Rossby-Haurwitz wave (mode 4): η(λ, φ) at t = $t"

        fig = plt.figure(figsize=(16, 9))
        n_subplots = length(projections)
        subplot_kwargs = (transform=ccrs.PlateCarree(), cmap=cmocean.cm.balance, vmin=8000-80, vmax=8000+80)

        for (n, projection) in enumerate(projections)
            ax = fig.add_subplot(1, n_subplots, n, projection=projection)
            plot_cubed_sphere_tracer_field!(fig, ax, η, file["grid"]; add_colorbar = (n == n_subplots), subplot_kwargs...)
            n_subplots == 1 && ax.set_title(plot_title)

            gl = ax.gridlines(color="gray", alpha=0.5, linestyle="--")
            gl.xlocator = mticker.FixedLocator(-180:30:180)
            gl.ylocator = mticker.FixedLocator(-80:20:80)
        end

        n_subplots > 1 && fig.suptitle(plot_title, y=0.85)

        filename = @sprintf("cubed_sphere_rossby_haurwitz_%04d.png", n)
        plt.savefig(filename, dpi=200, bbox_inches="tight")
        plt.close(fig)
    end

    close(file)

    filename_pattern = "cubed_sphere_rossby_haurwitz_%04d.png"
    output_filename  = "cubed_sphere_rossby_haurwitz.mp4"

    # Need extra crop video filter in case we end up with odd number of pixels in width or height.
    # See: https://stackoverflow.com/a/29582287
    run(`ffmpeg -y -i $filename_pattern -c:v libx264 -vf "fps=30, crop=trunc(iw/2)*2:trunc(ih/2)*2" -pix_fmt yuv420p $output_filename`)

    [rm(f) for f in glob("cubed_sphere_rossby_haurwitz_*.png")]

    return nothing
end


text:
function SQIL_callback(𝒟)
    𝒟[:r] .= 0
end

function SQIL(;π,
               S,
               𝒟_demo,
               normalize_demo::Bool=true,
               solver=SAC,
               log::NamedTuple=(;),
               kwargs...)

    !haskey(𝒟_demo, :r) && error("SQIL requires a reward value for the demonstrations")
    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π)))
    𝒟_demo = 𝒟_demo |> device(π)
    solver(;π=π,
            S=S,
            post_experience_callback=SQIL_callback,
            extra_buffers=[𝒟_demo],
            buffer_fractions=[1/2, 1/2],
            log=(dir="log/sqil", period=500, log...),
            kwargs...)
end


text:
## silly module to emulate D3 syntax.
module DThree

export AsIs, @I_str, D3

using JSON
using PyCall
using PySide
using Mustache

tpl = Mustache.template_from_file(Pkg.dir("PySide", "tpl", "d3.html"))

type AsIs
    x
end

macro I_str(x)
    AsIs(x)
end

## Take a function, get its args as array of symbols. There must be better way...
## Helper functions for bind callback
function get_args(li::LambdaStaticData)
    e = li.ast
    if !isa(e, Expr)
        e = Base.uncompressed_ast(li)
    end
    argnames = e.args[1]
    ## return array of symbols -- not args
    if isa(argnames[1], Expr)
        argnames = map(u -> u.args[1], argnames)
    end

    argnames
end

function get_args(m::Method)
    li = m.func.code
    get_args(li)
end

function get_args(f::Function)
    try
        get_args(f.env.defs.func)
    catch e
        get_args(f.code)
    end
end

## figure out how to make [0, I"d3.max(data)] -> [0, d3.max(data)]
## tojson(x::Any) = to_json(x)
## tojson(x::AsIs) = string(x.x)
## tojson(x::Tupe) = to_json([tojson(x) for x in x])

prep(d3, x::Any) = json(x)
prep(d3, x::AsIs) = string(x.x)
function prep(d3, f::Function)
    ## much trickier. We have julia object pushed into webpage which callback to this function
    julia = qnew_class_instance("Julia")
    qconnect(julia, :evaluate, (args...) -> f(args...))
    nm = randstring(10)
    d3.web[:page]()[:currentFrame]()[:addToJavaScriptWindowObject](nm, julia)
    args = join(map(string, get_args(f)), ", ")
    "function($args) {$nm.evaluate($args)}"
end

## D3 instances can chain commands with interface nearly same as d3.js
## d3.selectAll("p").style("color","#00f").render()
## render passes JavaScript to browser. If you want to assign to variable name, pass as argument.
## If method not defined, then can do something like d3._(:selectAll, "p"). ...
## by default this uses JavaScript d3 object as receiver. This can be changed with d3.receiver("chart"). ...
## arguments are converted via json, except:
## * functions are treated as callbacks into julia. These are asynchronous.
## * use I"x" or AsIs("x") to treat object "as is". This is needed to quote JavaScript functions
## TODO: generate these so members can be added from a list of symbols
type D3
    web
    cmd
    loaded
    receiver
    eval_js
    render
    _
    select
    selectAll
    data
    enter
    append
    style
    attr
    text
    domain
    range
    scale_linear
    svg_axis
    scale
    ticks
    orient
    call
    max
    function D3(parent::PyCall.PyObject, style::String,dataset::String)
        web = PySide.QtWebKit[:QWebView](parent)
        this = new(web, "d3", false)

        tmp = tempname() * ".html"

        io = open(tmp, "w")
        Mustache.render(io, tpl, Dict([(:style, style), (:dataset, dataset)]))
        close(io)

        qconnect(web, :loadFinished, (bool) -> this.loaded = bool)
        qinvoke(web, :load, QtCore[:QUrl]("file:///" * tmp))

        this.receiver = (value) -> begin this.cmd = value; this end
        this.eval_js = (value) -> begin
            #println(value)
            ## a bit kludgey, making sure page is loaded before evaluating
            if this.loaded
                this.web[:page]()[:currentFrame]()[:evaluateJavaScript](value)
                return(nothing)
            else
                ctr, max_steps = 1, 100
                timer = QtCore[:QTimer]()
                qconnect(timer, :timeout) do
                  ctr += 1
                  if this.loaded
                      timer[:stop]()
                      this.web[:page]()[:currentFrame]()[:evaluateJavaScript](value)
                  elseif ctr > max_steps
                      println("Failed to load page?")
                      timer[:stop]()
                  else
                      println(ctr)
                  end
                  return(nothing)
                end
                timer[:start](50)
            end
        end
        this.render = (args...) -> begin
            if length(args) > 0
                this.cmd = "var " * args[1] * " = " * this.cmd
            end
            this.eval_js(this.cmd)
            this.receiver("d3")
            nothing
        end
        this._ = (meth, args...) ->  begin
            args = map(u -> prep(this, u), args)
            args = join(args, ", ")
            cmd = this.cmd
            this.cmd = "$(this.cmd).$meth($args)"
            this
        end
        ## generate in @eval map meth by replacing _ with symbol(replace("scale_linear", "_", "."))
        this.select    = (args...) -> this._(:select, args...)
        this.selectAll = (args...) -> this._(:selectAll, args...)
        this.data      = (args...) -> this._(:data, args...)
        this.enter     = (args...) -> this._(:enter, args...)
        this.append    = (args...) -> this._(:append, args...)
        this.style     = (args...) -> this._(:style, args...)
        this.attr      = (args...) -> this._(:attr, args...)
        this.text      = (args...) -> this._(:text, args...)
        this.domain    = (args...) -> this._(:domain, args...)
        this.range     = (args...) -> this._(:range, args...)
        this.scale_linear = (args...) -> this._("scale.linear", args...)
        this.svg_axis  = (args...) -> this._("svg.axis", args...)
        this.scale     = (args...) -> this._(:scale, args...)
        this.ticks     = (args...) -> this._(:ticks, args...)
        this.orient    = (args...) -> this._(:orient, args...)
        this.call      = (args...) -> this._(:call, args...)
        this.max       = (args...) -> this._(:max, args...)

        this
    end
    D3(parent::PyCall.PyObject, dataset::String) = D3(parent, "", dataset)
end

end ## module


text:
using AQFED, Test

import AQFED.Random

function testSkip(N::Int)
   myrange = [1535, 1536, 2000, 2046, 2047, 2048, 2049, 2050, 100000, 100001]
   for i in myrange
      gen1 = AQFED.Random.Mixmax17()
      n = UInt64(i)
      AQFED.Random.skip(gen1, n)
      u1 = rand(gen1, UInt64)
      gen2 = AQFED.Random.Mixmax17()
      for j = 1:n
         rand(gen2, UInt64)
      end
      u2 = rand(gen2, UInt64)
      if u1 != u2
         println(i, " ", u1, " ", gen1, " ", rand(gen1, UInt64))
         println(i, " ", u2, " ", gen2, " ", rand(gen2, UInt64))
      end
      @test isequal(u1, u2)
   end
end

@testset "Mixmax17Skip" begin
testSkip(17)
end

@testset "Mixmax240Skip" begin
testSkip(240)
end

@testset "Mixmax17SkipTwice" begin
   N = 17
   myrange = [1535, 1536, 2000, 2046, 2047, 2048, 2049, 2050, 100000, 100001]
   gen1 = AQFED.Random.Mixmax17()
   gen2 = AQFED.Random.Mixmax17()
   for i in myrange
      n = UInt64(i)
      AQFED.Random.skip(gen1, n)
      u1 = rand(gen1, UInt64)
      for j = 1:n
         rand(gen2, UInt64)
      end
      u2 = rand(gen2, UInt64)
      if u1 != u2
         println(i, " ", u1, " ", gen1, " ", rand(gen1, UInt64))
         println(i, " ", u2, " ", gen2, " ", rand(gen2, UInt64))
      end
      @test isequal(u1, u2)
   end
end

@testset "Mixmax17Skip20" begin
   gen = AQFED.Random.Mixmax17()
   N = 17
   a0 = Vector(undef, N)
   AQFED.Random.iterateRawVec(gen)
   a0 = gen.v
   eArray = UInt64[
      843356255388265439,
      2303128525550447885,
      1403805614380726589,
      1646775209363665369,
      1656489684429545237,
      136962756483849705,
      578486145972684395,
      1186788852536058806,
      744606065061259291,
      129235498787499648,
      987271125576440333,
      2226012411473073226,
      1064839699140268466,
      1194475380559801816,
      1591779369875059359,
      1167968540715842034,
      2137709594798027207,
   ]
   a2p20 = zeros(UInt64, N)
   for i = 1:N
      for j = 1:N
         a2p20[j] = AQFED.Random.fmodmulM61(a2p20[j], eArray[i], a0[j])
      end
      AQFED.Random.iterateRawVec(gen)
      a0 = gen.v
   end
   println("a2p20 ", a2p20)
   gen = AQFED.Random.Mixmax17()
   for i = 1:(1+2^20)
      #skips 2^20 * 17 points
      AQFED.Random.iterateRawVec(gen)
      a0 = gen.v
      #end
   end
   for i = 1:N
      @test isequal(a0[i], a2p20[i])
   end
end


text:
using nauty_jll
using Parameters

const LIB_FILE = libnautyL0
const WORDSIZE = 64

@with_kw struct OptionBlk
    getcanon::Cint=false
    digraph::Cint=false
    writeautoms::Cint=false
    writemarkers::Cint=false
    defaultptn::Cint=true
    cartesian::Cint=false
    linelength::Cint=78
    outfile::Ptr{Cvoid}=C_NULL
    userrefproc::Ptr{Cvoid}=C_NULL
    userautomproc::Ptr{Cvoid}=C_NULL
    userlevelproc::Ptr{Cvoid}=C_NULL
    usernodeproc::Ptr{Cvoid}=C_NULL
    usercanonproc::Ptr{Cvoid}=C_NULL
    invarproc::Ptr{Cvoid}=C_NULL
    tc_level::Cint=100
    mininvarlevel::Cint=0
    maxinvarlevel::Cint=1
    invararg::Cint=0
    dispatch::Ptr{Cvoid}=C_NULL
    schreier::Cint=false
    extra_options::Ptr{Cvoid}=C_NULL
end

@with_kw struct StatsBlk
    grpsize1::Cdouble=0
    grpsize2::Cint=0
    numorbits::Cint=0
    numgenerators::Cint=0
    errstatus::Cint=0
    numnodes::Culong=0
    numbadleaves::Culong=0
    maxlevel::Cint=0
    tctotal::Culong=0
    canupdates::Culong=0
    invapplics::Culong=0
    invsuccesses::Culong=0
    invarsuclevel::Cint=0
end

function matrix2graph(matrix)
    n = size(matrix, 1)
    nw = WORDSIZE*cld(n, WORDSIZE)
    expanded = vcat(matrix, falses(nw-n, n))
    expanded.chunks .= expanded[end:-1:1, end:-1:1].chunks[end:-1:1]
    return expanded
end

function graph2matrix(graph)
    n = size(graph, 2)
    nw = WORDSIZE*cld(n, WORDSIZE)
    reverted = BitArray{2}(undef, size(graph)...)
    reverted.chunks = graph.chunks[end:-1:1]
    return reverted[end:-1:end+1-n, end:-1:end+1-n]
end

function canonical_form(matrix, colors; atol=0.1)
    n = size(matrix, 1)
    m = cld(n, WORDSIZE)
    lab = Vector{Cint}(sortperm(colors))
    ptn = Vector{Cint}(abs.(diff(colors[lab])) .<= atol)
    lab .= lab .- 1
    push!(ptn, 0)
    orbits = Vector{Cint}(undef, n)
    g = matrix2graph(matrix)
    cg = similar(g)
    dispatch_graph = cglobal((:dispatch_graph, LIB_FILE), Nothing)
    option_block = OptionBlk(getcanon=true, defaultptn=false, dispatch=dispatch_graph)
    ccall(
        (:densenauty, LIB_FILE),
        Cvoid,
        (Ptr{UInt64}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ref{OptionBlk}, Ref{StatsBlk}, Cint, Cint, Ptr{UInt64}),
        g.chunks, lab, ptn, orbits, option_block, StatsBlk(), m, n, cg.chunks
    )
    return lab .+ 1, graph2matrix(cg)
end


text:
module Enzyme

export Forward, Reverse
export Const, Active, Duplicated, DuplicatedNoNeed, BatchDuplicated, BatchDuplicatedNoNeed
export autodiff, jacobian, gradient, gradient!
export markType, batch_size, onehot, chunkedonehot

using LinearAlgebra

# Independent code, must be loaded before "compiler.jl"
include("pmap.jl")

using Adapt

"""
    abstract type Annotation{T}

Abstract type for [`autodiff`](@ref) function argument wrappers like
[`Const`](@ref), [`Active`](@ref) and [`Duplicated`](@ref).
"""
abstract type Annotation{T} end
Base.eltype(::Type{<:Annotation{T}}) where T = T

"""
    Const(x)

Mark a function argument `x` of [`autodiff`](@ref) as constant,
Enzyme will not auto-differentiate in respect `Const` arguments.
"""
struct Const{T} <: Annotation{T}
    val::T
end
Adapt.adapt_structure(to, x::Const) = Const(adapt(to, x.val))

# To deal with Const(Int) and prevent it to go to `Const{DataType}(T)`
Const(::Type{T}) where T = Const{Type{T}}(T)

"""
    Active(x)

Mark a function argument `x` of [`autodiff`](@ref) as active,
Enzyme will auto-differentiate in respect `Active` arguments.

!!! note

    Enzyme gradients with respect to integer values are zero.
    [`Active`](@ref) will automatically convert plain integers to floating
    point values, but cannot do so for integer values in tuples and structs.
"""
struct Active{T} <: Annotation{T}
    val::T
end
Adapt.adapt_structure(to, x::Active) = Active(adapt(to, x.val))

Active(i::Integer) = Active(float(i))

"""
    Duplicated(x, ∂f_∂x)

Mark a function argument `x` of [`autodiff`](@ref) as duplicated, Enzyme will
auto-differentiate in respect to such arguments, with `dx` acting as an
accumulator for gradients (so ``\\partial f / \\partial x`` will be *added to*)
`∂f_∂x`.
"""
struct Duplicated{T} <: Annotation{T}
    val::T
    dval::T
end
Adapt.adapt_structure(to, x::Duplicated) = Duplicated(adapt(to, x.val), adapt(to, x.dval))

"""
    DuplicatedNoNeed(x, ∂f_∂x)

Like [`Duplicated`](@ref), except also specifies that Enzyme may avoid computing
the original result and only compute the derivative values.
"""
struct DuplicatedNoNeed{T} <: Annotation{T}
    val::T
    dval::T
end
Adapt.adapt_structure(to, x::DuplicatedNoNeed) = DuplicatedNoNeed(adapt(to, x.val), adapt(to, x.dval))

"""
    BatchDuplicated(x, ∂f_∂xs)

Like [`Duplicated`](@ref), except contains several shadows to compute derivatives
for all at once. Argument `∂f_∂xs` should be a tuple of the several values of type `x`.
"""
struct BatchDuplicated{T,N} <: Annotation{T}
    val::T
    dval::NTuple{N,T}
end
Adapt.adapt_structure(to, x::BatchDuplicated) = BatchDuplicated(adapt(to, x.val), adapt(to, x.dval))

"""
    BatchDuplicatedNoNeed(x, ∂f_∂xs)

Like [`DuplicatedNoNeed`](@ref), except contains several shadows to compute derivatives
for all at once. Argument `∂f_∂xs` should be a tuple of the several values of type `x`.
"""
struct BatchDuplicatedNoNeed{T,N} <: Annotation{T}
    val::T
    dval::NTuple{N,T}
end
batch_size(::BatchDuplicated{T,N}) where {T,N} = N
batch_size(::BatchDuplicatedNoNeed{T,N}) where {T,N} = N
Adapt.adapt_structure(to, x::BatchDuplicatedNoNeed) = BatchDuplicatedNoNeed(adapt(to, x.val), adapt(to, x.dval))

"""
    abstract type Mode

Abstract type for what differentiation mode will be used.
"""
abstract type Mode end

"""
    struct Reverse <: Mode

Reverse mode differentiation
"""
struct ReverseMode <: Mode
end
const Reverse = ReverseMode()
guess_activity(::Type{T}, ::ReverseMode) where T = guess_activity(T)

"""
    struct Forward <: Mode

Forward mode differentiation
"""
struct ForwardMode <: Mode
end
const Forward = ForwardMode()
guess_activity(::Type{T}, ::ForwardMode) where T = guess_activity(T, API.DEM_ForwardMode)

import LLVM
include("api.jl")

@inline function guess_activity(::Type{T}, Mode::API.CDerivativeMode=API.DEM_ReverseModeCombined) where {T}
    return Const{T}
end
@inline function guess_activity(::Type{T}, Mode::API.CDerivativeMode=API.DEM_ReverseModeCombined) where {T<:AbstractFloat}
    if Mode == API.DEM_ForwardMode
        return DuplicatedNoNeed{T}
    else
        return Active{T}
    end
end
@inline function guess_activity(::Type{T}, Mode::API.CDerivativeMode=API.DEM_ReverseModeCombined) where {T<:Complex{<:AbstractFloat}}
    if Mode == API.DEM_ForwardMode
        return DuplicatedNoNeed{T}
    else
        return Active{T}
    end
end

@inline function guess_activity(::Type{T}, Mode::API.CDerivativeMode=API.DEM_ReverseModeCombined) where {T<:AbstractArray}
    if Mode == API.DEM_ForwardMode
        return DuplicatedNoNeed{T}
    else
        return Duplicated{T}
    end
end

include("logic.jl")
include("typeanalysis.jl")
include("typetree.jl")
include("utils.jl")
include("compiler.jl")

import .Compiler: CompilationException

# @inline annotate() = ()
# @inline annotate(arg::A, args::Vararg{Any, N}) where {A<:Annotation, N} = (arg, annotate(args...)...)
# @inline annotate(arg, args::Vararg{Any, N}) where N = (Const(arg), annotate(args...)...)

@inline function annotate(args::Vararg{Any, N}) where N
    ntuple(Val(N)) do i
        Base.@_inline_meta
        arg = @inbounds args[i]
        if arg isa Annotation
            return arg
        else
            return Const(arg)
        end
    end
end

@inline function same_or_one(args::Vararg{Any, N}) where N
    current = -1
    for arg in args
        if arg isa BatchDuplicated
            if current == -1
                current = batch_size(arg)
            else
                @assert current == batch_size(arg)
            end
        elseif arg isa BatchDuplicatedNoNeed
            if current == -1
                current = batch_size(arg)
            else
                @assert current == batch_size(arg)
            end
        end
    end

    if current == -1
        current = 1
    end

    return current
end

"""
    autodiff(::ReverseMode, f, Activity, args...)

Auto-differentiate function `f` at arguments `args` using reverse mode.

Limitations:

* `f` may only return a `Real` (of a built-in/primitive type) or `nothing`,
  not an array, struct, `BigFloat`, etc. To handle vector-valued return
  types, use a mutating `f!` that returns `nothing` and stores it's return
  value in one of the arguments, which must be wrapped in a
  [`Duplicated`](@ref).

`args` may be numbers, arrays, structs of numbers, structs of arrays and so
on. Enzyme will only differentiate in respect to arguments that are wrapped
in an [`Active`](@ref) (for immutable arguments like primitive types and
structs thereof) or [`Duplicated`](@ref) (for mutable arguments like arrays,
`Ref`s and structs thereof). Non-annotated arguments will automatically be
treated as [`Const`](@ref).

`Activity` is the Activity of the return value, it may be `Const` or `Active`.

Example:

```jldoctest
a = 4.2
b = [2.2, 3.3]; ∂f_∂b = zero(b)
c = 55; d = 9

f(a, b, c, d) = a * √(b[1]^2 + b[2]^2) + c^2 * d^2
∂f_∂a, ∂f_∂d = autodiff(Reverse, f, Active, Active(a), Duplicated(b, ∂f_∂b), c, Active(d))

# output

(3.966106403010388, 54450.0)
```

here, `autodiff` returns a tuple
``(\\partial f/\\partial a, \\partial f/\\partial d)``,
while ``\\partial f/\\partial b`` will be *added to* `∂f_∂b` (but not returned).
`c` will be treated as `Const(c)`.

!!! note

    Enzyme gradients with respect to integer values are zero.
    [`Active`](@ref) will automatically convert plain integers to floating
    point values, but cannot do so for integer values in tuples and structs.
"""
@inline function autodiff(::ReverseMode, f::F, ::Type{A}, args...) where {F, A<:Annotation}
    args′  = annotate(args...)
    tt′    = Tuple{map(Core.Typeof, args′)...}
    width = Val(same_or_one(args...))
    if A <: Active
        tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
        rt = Core.Compiler.return_type(f, tt)
        if !allocatedinline(rt)
            forward, adjoint = Enzyme.Compiler.thunk(f, #=df=#nothing, Duplicated{rt}, tt′, #=Split=# Val(API.DEM_ReverseModeGradient), width, #=ModifiedBetween=#Val(false))
            res = forward(args′...)
            @assert length(res) == 2
            tape = res[1]
            if res[2] isa Base.RefValue
                res[2][] += one(eltype(typeof(res[2])))
            else
                res[2] += one(eltype(typeof(res[2])))
            end
            return adjoint(args′..., tape)
        end
    elseif A <: Duplicated || A<: DuplicatedNoNeed || A <: BatchDuplicated || A<: BatchDuplicatedNoNeed
        throw(ErrorException("Duplicated Returns not yet handled"))
    end
    thunk = Enzyme.Compiler.thunk(f, #=df=#nothing, A, tt′, #=Split=# Val(API.DEM_ReverseModeCombined), width)
    if A <: Active
        tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
        rt = eltype(Compiler.return_type(thunk))
        args′ = (args′..., one(rt))
    end
    thunk(args′...)
end

@inline function autodiff(::ReverseMode, dupf::Duplicated{F}, ::Type{A}, args...) where {F, A<:Annotation}
    args′  = annotate(args...)
    tt′    = Tuple{map(Core.Typeof, args′)...}
    width = Val(same_or_one(args...))
    thunk = Enzyme.Compiler.thunk(#=f=#dupf.val, #=df=#dupf.dval, A, tt′, #=Split=# Val(API.DEM_ReverseModeCombined), width)
    if A <: Active
        rt = eltype(Compiler.return_type(thunk))
        args′ = (args′..., one(rt))
    end
    thunk(args′...)
end

"""
    autodiff(mode::Mode, f, args...)
    autodiff(f, mode::Mode, args...)

Like [`autodiff`](@ref) but will try to guess the activity of the return value.
"""
@inline function autodiff(mode::Mode, f::F, args...) where {F}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(f, tt)
    A     = guess_activity(rt, mode)
    autodiff(mode, f, A, args′...)
end

@inline function autodiff(mode::Mode, dupf::Duplicated{F}, args...) where {F}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(dupf.val, tt)
    A     = guess_activity(rt, mode)
    autodiff(mode, dupf, A, args′...)
end

"""
    autodiff(::ForwardMode, f, Activity, args...)

Auto-differentiate function `f` at arguments `args` using forward mode.

`args` may be numbers, arrays, structs of numbers, structs of arrays and so
on. Enzyme will only differentiate in respect to arguments that are wrapped
in a [`Duplicated`](@ref) or similar argument. Non-annotated arguments will
automatically be treated as [`Const`](@ref). Unlike reverse mode in
[`autodiff`](@ref), [`Active`](@ref) arguments are not allowed here, since
all

`Activity` is the Activity of the return value, it may be:
* `Const` if the return is not to be differentiated with respect to
* `Duplicated`, if the return is being differentiated with respect to and
  both the original value and the derivative return are desired
* `DuplicatedNoNeed`, if the return is being differentiated with respect to
  and only the derivative return is desired.

Example returning both original return and derivative:

```jldoctest
a = 4.2
b = [2.2, 3.3]; ∂f_∂b = zero(b)
c = 55; d = 9

f(x) = x*x
res, ∂f_∂x = autodiff(Forward, f, Duplicated, Duplicated(3.14, 1.0))

# output

(9.8596, 6.28)
```

Example returning just the derivative:

```jldoctest
a = 4.2
b = [2.2, 3.3]; ∂f_∂b = zero(b)
c = 55; d = 9

f(x) = x*x
∂f_∂x = autodiff(Forward, f, DuplicatedNoNeed, Duplicated(3.14, 1.0))

# output

(6.28,)
```
"""
@inline function autodiff(::ForwardMode, f::F, ::Type{A}, args...) where {F, A<:Annotation}
    args′  = annotate(args...)
    tt′    = Tuple{map(Core.Typeof, args′)...}
    width = Val(same_or_one(args...))
    if A <: Active
        throw(ErrorException("Active Returns not allowed in forward mode"))
    end
    ReturnPrimal = Val(A <: Duplicated || A <: BatchDuplicated)
    thunk = Enzyme.Compiler.thunk(f, #=df=#nothing, A, tt′, #=Mode=# Val(API.DEM_ForwardMode), width,
                                     #=ModifiedBetween=#Val(false), ReturnPrimal)
    thunk(args′...)
end

@inline function autodiff(::ForwardMode, dupf::Duplicated{F}, ::Type{A}, args...) where {F, A<:Annotation}
    args′  = annotate(args...)
    tt′    = Tuple{map(Core.Typeof, args′)...}
    width = Val(same_or_one(args...))
    if A <: Active
        throw(ErrorException("Active Returns not allowed in forward mode"))
    end
    ReturnPrimal = Val(A <: Duplicated || A <: BatchDuplicated)
    thunk = Enzyme.Compiler.thunk(#=f=#dupf.val, #=df=#dupf.dval, A, tt′, #=Mode=# Val(API.DEM_ForwardMode), width, #=ModifiedBetween=#Val(false), ReturnPrimal)
    thunk(args′...)
end

# Compat
@deprecate fwddiff(f::F, ::Type{A}, args...) where {F, A<:Annotation} autodiff(Forward, f, A, args...)
@deprecate fwddiff(f::F, args...) where {F} autodiff(Forward, f, args...)

# Compat
@inline autodiff(f::F, ::Type{A}, args...) where {F, A<:Annotation} = autodiff(Reverse, f, A, args...)
@inline autodiff(f::F, args...) where {F} = autodiff(Reverse, f, args...)

# F as first arg for `do` syntax
@inline autodiff(dupf::Duplicated{F}, mode::Mode, ::Type{A}, args...) where {F,A<:Annotation} = autodiff(mode, dupf, A, args...)
@inline autodiff(f::F, mode::Mode, ::Type{A}, args...) where {F,A<:Annotation} = autodiff(mode, f, A, args...)
@inline autodiff(dupf::Duplicated{F}, mode::Mode, args...) where {F} = autodiff(mode, dupf, args...)
@inline autodiff(f::F, mode::Mode, args...) where {F} = autodiff(mode, f, args...)

"""
    autodiff_deferred(f, Activity, args...)

Same as [`autodiff`](@ref) but uses deferred compilation to support usage in GPU
code, as well as high-order differentiation.
"""
@inline function autodiff_deferred(f::F, ::Type{A}, args...) where {F, A<:Annotation}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    width = Val(same_or_one(args...))
    if A isa UnionAll
        tt = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
        rt = Core.Compiler.return_type(f, tt)
        rt = A{rt}
    else
        @assert A isa DataType
        rt = A
    end

    if eltype(rt) == Union{}
        error("Return type inferred to be Union{}. Giving up.")
    end

    ReturnPrimal = Val(false)
    ptr   = Compiler.deferred_codegen(f, Val(tt′), Val(rt), #=dupClosure=#Val(false), Val(API.DEM_ReverseModeCombined), width, #=ModifiedBetween=#Val(false), ReturnPrimal)
    thunk = Compiler.CombinedAdjointThunk{F, rt, tt′, typeof(width), Nothing, ReturnPrimal}(f, ptr, #=df=#nothing)
    if rt <: Active
        args′ = (args′..., one(eltype(rt)))
    elseif A <: Duplicated || A<: DuplicatedNoNeed || A <: BatchDuplicated || A<: BatchDuplicatedNoNeed
        throw(ErrorException("Duplicated Returns not yet handled"))
    end
    thunk(args′...)
end

"""
    autodiff_deferred(f, args...)

Like [`autodiff_deferred`](@ref) but will try to guess the activity of the return value.
"""
@inline function autodiff_deferred(f::F, args...) where {F}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(f, tt)
    rt    = guess_activity(rt)
    autodiff_deferred(f, rt, args′...)
end

"""
    fwddiff_deferred(f, Activity, args...)

Same as `autodiff(::ForwardMode, ...)` but uses deferred compilation to support usage in GPU
code, as well as high-order differentiation.
"""
@inline function fwddiff_deferred(f::F, ::Type{A}, args...) where {F, A<:Annotation}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    width = Val(same_or_one(args...))
    if A isa UnionAll
        tt = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
        rt = Core.Compiler.return_type(f, tt)
        rt = A{rt}
    else
        @assert A isa DataType
        rt = A
    end

    if eltype(rt) == Union{}
        error("Return type inferred to be Union{}. Giving up.")
    end

    if A <: Active
        throw(ErrorException("Active Returns not allowed in forward mode"))
    end

    ReturnPrimal = Val(A <: Duplicated || A <: BatchDuplicated)
    ptr   = Compiler.deferred_codegen(f, Val(tt′), Val(rt), #=dupClosure=#Val(false), Val(API.DEM_ForwardMode), width, #=ModifiedBetween=#Val(false), ReturnPrimal)
    thunk = Compiler.ForwardModeThunk{F, rt, tt′, typeof(width), Nothing, ReturnPrimal}(f, ptr, #=df=#nothing)
    thunk(args′...)
end

"""
    fwddiff_deferred(f, args...)

Like [`fwddiff_deferred`](@ref) but will try to guess the activity of the return value.
"""
@inline function fwddiff_deferred(f::F, args...) where {F}
    args′ = annotate(args...)
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(f, tt)
    rt    = guess_activity(rt, API.DEM_ForwardMode)
    fwddiff_deferred(f, rt, args′...)
end

# White lie, should be `Core.LLVMPtr{Cvoid, 0}` but that's not supported by ccallable
Base.@ccallable function __enzyme_float(x::Ptr{Cvoid})::Cvoid
    return nothing
end

Base.@ccallable function __enzyme_double(x::Ptr{Cvoid})::Cvoid
    return nothing
end

@inline function markType(::Type{T}, ptr::Ptr{Cvoid}) where T
    markType(Base.unsafe_convert(Ptr{T}, ptr))
end

@inline function markType(data::Array{T}) where T
    GC.@preserve data markType(pointer(data))
end

# TODO(WM): We record the type of a single index here, we could give it a range
@inline function markType(data::SubArray)
    GC.@preserve data markType(pointer(data))
end

@inline function markType(data::Ptr{Float32})
    Base.llvmcall(("declare void @__enzyme_float(i8* nocapture) nounwind define void @c(i64 %q) nounwind alwaysinline { %p = inttoptr i64 %q to i8* call void @__enzyme_float(i8* %p) ret void }", "c"), Cvoid, Tuple{Ptr{Float32}}, data)
    nothing
end

@inline function markType(data::Ptr{Float64})
    Base.llvmcall(("declare void @__enzyme_double(i8* nocapture) nounwind define void @c(i64 %q) nounwind alwaysinline { %p = inttoptr i64 %q to i8* call void @__enzyme_double(i8* %p) ret void }", "c"), Cvoid, Tuple{Ptr{Float64}}, data)
    nothing
end

@inline function onehot(x, start=1, endl=length(x))
    ntuple(Val(endl-start+1)) do i
        Base.@_inline_meta
        res = similar(x)
        for idx in 1:length(x)
            @inbounds res[idx] = (i + start - 1== idx) ? 1.0 : 0.0
        end
        return res
    end
end

@inline function onehot(x::NTuple{N, T}, start=1, endl=N) where {T, N}
    ntuple(Val(endl-start+1)) do i
        Base.@_inline_meta
        ntuple(N) do idx
            Base.@_inline_meta
            return (i + start - 1 == idx) ? 1.0 : 0.0
        end
    end
end

"""
    gradient(::ReverseMode, f, x)

Compute the gradient of an array-input function `f` using reverse mode.
This will allocate and return new array with the gradient result.

Example:

```jldoctest
f(x) = x[1]*x[2]

grad = gradient(Reverse, f, [2.0, 3.0])

# output

2-element Vector{Float64}:
 3.0
 2.0
```
"""
@inline function gradient(::ReverseMode, f, x)
    dx = zero(x)
    autodiff(Reverse, f, Duplicated(x, dx))
    dx
end

"""
    gradient!(::ReverseMode, dx, f, x)

Compute the gradient of an array-input function `f` using reverse mode,
storing the derivative result in an existing array `dx`.

Example:

```jldoctest
f(x) = x[1]*x[2]

dx = [0.0, 0.0]
gradient!(Reverse, dx, f, [2.0, 3.0])

# output

2-element Vector{Float64}:
 3.0
 2.0
```
"""
@inline function gradient!(::ReverseMode, dx, f, x)
    dx .= 0
    autodiff(Reverse, f, Duplicated(x, dx))
    dx
end

"""
    gradient(::ForwardMode, f, x; shadow=onehot(x))

Compute the gradient of an array-input function `f` using forward mode. The
optional keyword argument `shadow` is a vector of one-hot vectors of type `x`
which are used to forward-propagate into the return. For performance reasons,
this should be computed once, outside the call to `gradient`, rather than
within this call.

Example:

```jldoctest
f(x) = x[1]*x[2]

grad = gradient(Forward, f, [2.0, 3.0])

# output

(3.0, 2.0)
```
"""
@inline function gradient(::ForwardMode, f, x; shadow=onehot(x))
    only(autodiff(Forward, f, BatchDuplicatedNoNeed, BatchDuplicated(x, shadow)))
end

@inline function chunkedonehot(x, ::Val{chunk}) where chunk
    sz = length(x)
    num = ((sz + chunk - 1) ÷ chunk)
    ntuple(Val(num)) do i
        Base.@_inline_meta
        onehot(x, (i-1)*chunk+1, i == num ? sz : (i*chunk) )
    end
end

@inline tupleconcat(x) = x
@inline tupleconcat(x, y) = (x..., y...)
@inline tupleconcat(x, y, z...) = (x..., tupleconcat(y, z...)...)

"""
    gradient(::ForwardMode, f, x, ::Val{chunk}; shadow=onehot(x))

Compute the gradient of an array-input function `f` using vector forward mode.
Like [`gradient`](@ref), except it uses a chunk size of `chunk` to compute
`chunk` derivatives in a single call.

Example:

```jldoctest
f(x) = x[1]*x[2]

grad = gradient(Forward, f, [2.0, 3.0], Val(2))

# output

(3.0, 2.0)
```
"""
@inline function gradient(::ForwardMode, f, x, ::Val{chunk}; shadow=chunkedonehot(x, Val(chunk))) where chunk
    tmp = ntuple(length(shadow)) do i
        autodiff(Forward, f, BatchDuplicatedNoNeed, BatchDuplicated(x, shadow[i]))[1]
    end
    tupleconcat(tmp...)
end

@inline function gradient(::ForwardMode, f, x, ::Val{1}; shadow=onehot(x))
    ntuple(length(shadow)) do i
        autodiff(Forward, f, DuplicatedNoNeed, Duplicated(x, shadow[i]))[1]
    end
end

"""
    jacobian(::ForwardMode, f, x; shadow=onehot(x))
    jacobian(::ForwardMode, f, x, ::Val{chunk}; shadow=onehot(x))

Compute the jacobian of an array-input function `f` using (potentially vector)
forward mode. This is a simple rename of the [`gradient`](@ref) function,
and all relevant arguments apply here.

Example:

```jldoctest
f(x) = [x[1]*x[2], x[2]]

grad = jacobian(Forward, f, [2.0, 3.0])

# output

2×2 Matrix{Float64}:
 3.0  2.0
 0.0  1.0
```
"""
@inline function jacobian(::ForwardMode, args...; kwargs...)
    cols = gradient(Forward, args...; kwargs...)
    reduce(hcat, cols)
end

"""
    jacobian(::ReverseMode, f, x, ::Val{num_outs}, ::Val{chunk})

Compute the jacobian of an array-input function `f` using (potentially vector)
reverse mode. The `chunk` argument denotes the chunk size to use and `num_outs`
denotes the number of outputs `f` will return in an array.

Example:

```jldoctest
f(x) = [x[1]*x[2], x[2]]

grad = jacobian(Reverse, f, [2.0, 3.0], Val(2))

# output

2×2 Matrix{Float64}:
 3.0  2.0
 0.0  1.0
```
"""
@inline function jacobian(::ReverseMode, f, x, n_outs::Val{n_out_val}, ::Val{chunk}) where {chunk, n_out_val}
    num = ((n_out_val + chunk - 1) ÷ chunk)

    tt′   = Tuple{BatchDuplicated{Core.Typeof(x), chunk}}
    tt    = Tuple{Core.Typeof(x)}
    rt = Core.Compiler.return_type(f, tt)
    primal, adjoint = Enzyme.Compiler.thunk(f, #=df=#nothing, BatchDuplicatedNoNeed{rt}, tt′, #=Split=# Val(API.DEM_ReverseModeGradient), #=width=#Val(chunk), #=ModifiedBetween=#Val(false))

    if num * chunk == n_out_val
        last_size = chunk
        primal2, adjoint2 = primal, adjoint
    else
        last_size = n_out_val - (num-1)*chunk
        tt′ = Tuple{BatchDuplicated{Core.Typeof(x), last_size}}
        primal2, adjoint2 = Enzyme.Compiler.thunk(f, #=df=#nothing, BatchDuplicatedNoNeed{rt}, tt′, #=Split=# Val(API.DEM_ReverseModeGradient), #=width=#Val(last_size), #=ModifiedBetween=#Val(false))
    end

    tmp = ntuple(num) do i
        Base.@_inline_meta
        dx = ntuple(i == num ? last_size : chunk) do idx
            Base.@_inline_meta
            zero(x)
        end
        res = (i == num ? primal2 : primal)(BatchDuplicated(x, dx))
        tape = res[1]
        j = 0
        for shadow in res[2]
            j += 1
            @inbounds shadow[(i-1)*chunk+j] += one(eltype(typeof(shadow)))
        end
        (i == num ? adjoint2 : adjoint)(BatchDuplicated(x, dx), tape)
        return dx
    end
    rows = tupleconcat(tmp...)
    mapreduce(LinearAlgebra.adjoint, vcat, rows)
end

@inline function jacobian(::ReverseMode, f, x, n_outs::Val{n_out_val}, ::Val{1} = Val(1)) where {n_out_val}
    tt′   = Tuple{Duplicated{Core.Typeof(x)}}
    tt    = Tuple{Core.Typeof(x)}
    rt = Core.Compiler.return_type(f, tt)
    primal, adjoint = Enzyme.Compiler.thunk(f, #=df=#nothing, DuplicatedNoNeed{rt}, tt′, #=Split=# Val(API.DEM_ReverseModeGradient), #=width=#Val(1), #=ModifiedBetween=#Val(false))
    rows = ntuple(n_outs) do i
        Base.@_inline_meta
        dx = zero(x)
        res = primal(Duplicated(x, dx))
        tape = res[1]
        @inbounds res[2][i] += one(eltype(typeof(res[2])))
        adjoint(Duplicated(x, dx), tape)
        return dx
    end
    mapreduce(LinearAlgebra.adjoint, vcat, rows)
end

end # module


text:
using LambdaMaker
using Test

SAMPLE_FILENAME = "Dockerfile"
PROJECT_NAME = "foobar"

@testset "_copy_file" begin
    mktempdir() do RESOURCE_DIR
        file_path = joinpath(RESOURCE_DIR, SAMPLE_FILENAME)

        @testset "dest file does not exist" begin
            LambdaMaker._copy_file(RESOURCE_DIR, SAMPLE_FILENAME)
            @test isfile(file_path)
        end

        @testset "dest file exists" begin
            @test_throws ArgumentError LambdaMaker._copy_file(RESOURCE_DIR, SAMPLE_FILENAME)
        end

        @testset "force=true" begin
            LambdaMaker._copy_file(RESOURCE_DIR, SAMPLE_FILENAME; force=true)
            @test isfile(file_path)
        end
    end
end

@testset "_render_file" begin
    mktempdir() do RESOURCE_DIR
        render_args = Dict("PKG" => PROJECT_NAME)
        LambdaMaker._render_file(SAMPLE_FILENAME, RESOURCE_DIR, render_args)

        result_file_path = joinpath(RESOURCE_DIR, SAMPLE_FILENAME)
        @test isfile(result_file_path)

        contents = read(result_file_path, String)
        @test contains(contents, "CMD [\"$(PROJECT_NAME).handle_event\"]")
    end
end

@testset "create_lambda_package" begin
    function get_project_files()
        for (_, _, files) in walkdir(joinpath(@__DIR__, "..", "template_files"))
            return files
        end
    end

    function test_file_existance()
        for project_file in get_project_files()
            @test isfile(joinpath(PROJECT_NAME, project_file))
        end
    end

    @testset "project does not exist" begin
        try
            create_lambda_package(PROJECT_NAME)
            test_file_existance()
        finally
            rm(PROJECT_NAME; force=true, recursive=true)
        end
    end

    @testset "project exists & force=false" begin
        try
            mkdir(PROJECT_NAME)
            @test_throws DirectoryExists create_lambda_package(PROJECT_NAME)
        finally
            rm(PROJECT_NAME)
        end
    end

    @testset "project exists & force=true" begin
        try
            mkdir(PROJECT_NAME)
            create_lambda_package(PROJECT_NAME; force=true)
            test_file_existance()
        finally
            rm(PROJECT_NAME; force=true, recursive=true)
        end
    end
end


text:
length("Hello, world!") # gives 13
length("Hellö, wørld!") # gives 13

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# NMF (Non-negative Matrix Factorization)
# Julia実装版

module NMF

export nmf_euc, nmf_kl, nmf_is

# EUC基準NMF処理（通常版）
# V     分解対象の行列
# k     因子分解数
# iter  イテレーション回数
function nmf_euc(V::AbstractMatrix,k::Int64=3;iter::Int64=100)
    # 変数チェック処理
    @assert countnz(V.< 0) == 0 "Input Matrix V is not non-negative"
    @assert k > 0 "Factor Num : k must be greater than 0"
    @assert iter > 0 "Iteration Num : iter must be greater than 0"

    # W,Hを指定の形で乱数作成
    W = rand(size(V)[1],k)
    H = rand(k,size(V)[2])

    # 反復乗法更新
    for i in 1:iter
        W = W .* (V * H') ./ (W * H * H')
        H = H .* (W' * V) ./ (W' * W * H)
    end
    return W,H
end

# EUC基準NMF処理（テンプレートあり版）
# V     分解対象の行列
# W     テンプレート辞書行列
# iter  イテレーション回数
function nmf_euc(V::AbstractMatrix,W::AbstractMatrix,;iter::Int64=100)
    #変数チェック処理
    @assert countnz(V.< 0) == 0 "Input Matrix V is not non-negative."
    @assert countnz(W.< 0) == 0 "Input Matrix W is not non-negative."
    @assert iter > 0 "Iteration Num : iter must be greater than 0."
    @assert size(V,1) == size(W,1) "Dimensions of V:$(size(V)) and W:$(size(W)) do not match. Align the first dimention."
    @assert size(W,2) > 0 "Factor Num : size(W,2) is greater than 0."

    # 指定の形でHを作成
    H = rand(size(W,2),size(V)[2])

    # H だけを更新
    for i in 1:iter
        H = H .* (W' * V) ./ (W' * W * H)
    end
    return H
end

# KL基準NMF処理（通常版）
# V     分解対象の行列
# k     因子分解数
# iter  イテレーション回数
function nmf_kl(V::AbstractMatrix,k::Int64=3;iter::Int64=100)
    #変数チェック処理
    @assert countnz(V.< 0) == 0 "Input Matrix V is not non-negative"
    @assert k > 0 "Factor Num : k must be greater than 0"
    @assert iter > 0 "Iteration Num : iter must be greater than 0"

    # W,Hを指定の形で乱数作成
    W = rand(size(V)[1],k)
    H = rand(k,size(V)[2])

    # 反復乗法更新
    for i in 1:iter
        W = W .* ((V ./ (W*H)) * H') ./ sum(H',1)
        H = H .* (W' * (V ./ (W*H))) ./ sum(W',2)
    end
    return W,H
end

# KL基準NMF処理（テンプレートあり版）
# V     分解対象の行列
# W     テンプレート辞書行列
# iter  イテレーション回数
function nmf_kl(V::AbstractMatrix,W::AbstractMatrix;iter::Int64=100)
    #変数チェック処理
    @assert countnz(V.< 0) == 0 "Input Matrix V is not non-negative."
    @assert countnz(W.< 0) == 0 "Input Matrix W is not non-negative."
    @assert iter > 0 "Iteration Num : iter must be greater than 0."
    @assert size(V,1) == size(W,1) "Dimensions of V:$(size(V)) and W:$(size(W)) do not match. Align the first dimention."
    @assert size(W,2) > 0 "Factor Num : size(W,2) is greater than 0."

    # 指定の形でHを作成
    H = rand(size(W,2),size(V)[2])

    # H だけを更新
    for i in 1:iter
        H = H .* (W' * (V ./ (W*H))) ./ sum(W',2)
    end
    return H
end

# IS基準NMF処理（通常版）
# V     分解対象の行列
# k     因子分解数
# iter  イテレーション回数
function nmf_is(V::AbstractMatrix,k::Int64=3;iter::Int64=100)
    #変数チェック処理
    @assert countnz(V.< 0) == 0 "Input Matrix V is not non-negative"
    @assert k > 0 "Factor Num : k must be greater than 0"
    @assert iter > 0 "Iteration Num : iter must be greater than 0"

    # W,Hを指定の形で乱数作成
    W = rand(size(V)[1],k)
    H = rand(k,size(V)[2])

    # 反復乗法更新
    for i in 1:iter
        W = W .* sqrt.(((V ./ (W*H))*(H./sum(W*H,1))') ./(sum((H./sum(W*H,1))',1)))
        H = H.*sqrt.(((W ./ sum(W*H,2))' * (V./(W*H)) ./ (sum((W ./ sum(W*H,2))',2))))
    end
    return W,H
end

# IS基準NMF処理（テンプレートあり版）
# V     分解対象の行列
# W     テンプレート辞書行列
# iter  イテレーション回数
function nmf_is(V::AbstractMatrix,W::AbstractMatrix;iter::Int64=100)
    #変数チェック処理
    @assert countnz(V.< 0) == 0 "Input Matrix V is not non-negative."
    @assert countnz(W.< 0) == 0 "Input Matrix W is not non-negative."
    @assert iter > 0 "Iteration Num : iter must be greater than 0."
    @assert size(V,1) == size(W,1) "Dimensions of V:$(size(V)) and W:$(size(W)) do not match. Align the first dimention."
    @assert size(W,2) > 0 "Factor Num : size(W,2) is greater than 0."

    # 指定の形でHを作成
    H = rand(size(W,2),size(V)[2])

    # H だけを更新
    for i in 1:iter
        H = H.*sqrt.(((W ./ sum(W*H,2))' * (V./(W*H)) ./ (sum((W ./ sum(W*H,2))',2))))
    end
    return H
end

end


text:
"""
    localizationsplot(localizations::Vector{Localizations}, color::Symbol)
    localizationsplot(localizations::Vector{Vector{Localization}}, colors::Vector{Symbol})

Plot localizations by spatial coordinates. Assumes that the y-coordinate is inverted. Defaults to 0-to-40960 on each
axis, but option may be specified using keyword arguments xlims=(x1, x2) and ylims=(y1, y2). Colors are specified in the
documentation for Plots.jl. Can tune opacity with opacity= keword argument. factor= can be tuned for display in the pane
or for saving and printing (recommended 4, default 1).
"""
plotlocalizations(locs::Vector{Localization}, color::Symbol; kargs...) =
    plotlocalizations([locs], [color]; kargs...)

function plotlocalizations(localizations::Vector{Vector{Localization}}, colors::Vector{Symbol};
    xlims=(0, 40960), ylims=(0,40960), opacity=0.75, factor=1)

    channelcount = length(localizations)
    plot(framestyle=:none, size=((1024,1024) .* factor))
    for i ∈ 1:channelcount
        coords = localizations[i] |> extractcoordinates
        scatter!(coords[1,:], coords[2,:], marker=(1 * factor, stroke(0), opacity, colors[i]))
    end
    plot!(aspect_ratio=:equal, xlims=xlims, yaxis=(ylims, :flip), legend=:none, grid=:hide, ticks=(0))
end


text:
struct SineWave
    amplitude::Float32
    phase_shift::Float32
end
SineWave() = SineWave(rand(Uniform(0.1, 5)), rand(Uniform(0, 2pi)))

(s::SineWave)(x::AbstractArray) = s.amplitude .* sin.(x .+ s.phase_shift)

function Base.show(io::IO, s::SineWave)
    print(io, "SineWave(amplitude = ", s.amplitude, ", phase shift = ", s.phase_shift, ")")
end

max_stars_repo_name: JuliaTagBot/MetaLearning

max_issues_repo_name: JuliaTagBot/MetaLearning

max_forks_repo_name: JuliaTagBot/MetaLearning


text:
export anagram_iter

function tabulate_decompositions(bag::Bag, decomp_memo=Dict{Bag,Vector{Bag}}(),
                                           vocab_dict=VOCAB_DICT)

    for subbag in AllSubBagIter(bag)

        # Check whether this subbag has any valid decompositions
        for a in AllSubBagIter(subbag)
            if haskey(vocab_dict, a)
                b = bag_complement(subbag, a)
                if haskey(decomp_memo, b) | all(b .== 0)
                    # push it onto the list of decompositions
                    if haskey(decomp_memo, subbag)
                        push!(decomp_memo[subbag], a)
                    else
                        decomp_memo[subbag] = Bag[a]
                    end
                end
            end
        end
    end

    return decomp_memo
end

function anagram_iter(text::String)
    bag = word_to_bag(text)
    decomp_cache = tabulate_decompositions(bag)
    return AnagramIter(bag, decomp_cache)
end


text:
# Many derivative free methods of different orders
#
# TODO: rework Order5 #https://pdfs.semanticscholar.org/ce50/3210d96f653a14b28da96600d5990d2abe97.pdf
# https://content.sciendo.com/view/journals/tmj/10/4/article-p103.xml 7 and 8
# order8: https://www.hindawi.com/journals/ijmms/2012/493456/ref/

##################################################
## Guard against non-robust algorithms
##
## By default, we do this by deciding if we
## should take a secant step instead of the algorithm For example, for
## Steffensen which is quadratically convergent and the Secant method
## which is only superlinear,
## the error, e_{n+1} = x_{n+1} - alpha, may be smaller after a secant
## step than a Steffensen step. (It is only once x_n is close enough
## to alpha that the method is quadratically convergent.
## The Steffensen error is
## Δn+1 = f[x,x+fx, alpha]/f[x, x+fx] * (1 - f[x, alpha]) (x-alpha)^2
##      ≈ f''/(2f') * ( 1 + f') Δn^2
## The Secant error is
## Δn+1 = f[x,x_{-1},alpha] / f[x,x_{-1}] * (x-alpha) * (x_{-1} - alpha)
##      ≈  f''/(2f')  Δn ⋅ Δn-1
## The ratio is ≈ (1 + f')(Δn / Δn-1)
## It seems reasonable, that a Steffensen step is preferred when
## the ratio satisfies -1 < (1+f') ⋅ Δn /Δn-1 < 1
## We could use f' ~ fp = (fx1-fx0)/(x1-x0); but our proxy for
## Δn/Δn-1 is problematic, as we don't know alpha, and using xn-x_{n-1}
## can be an issue when only x1 and not x0 is specified. This needs
## working around.
##
## Instead, as Steffensen is related to Newton as much as
## (f(x+fx) - fx)/fx  ≈ f'(x), we take a Steffensen step if |fx|
## is small enough. For this we use |fx| <= x/1000; which
## seems to work reasonably well over several different test cases.
@inline function do_guarded_step(M::AbstractSecant, o::AbstractUnivariateZeroState{T,S}) where {T, S}
    x, fx = o.xn1, o.fxn1
    1000 * abs(fx) >  max(oneunit(S), abs(x) * oneunit(S) /oneunit(T)) * one(T)
end

# check if we should guard against step for method M; call N if yes, P if not
function update_state_guarded(M::AbstractSecant,N::AbstractUnivariateZeroMethod, P::AbstractUnivariateZeroMethod, fs, o, options, l=NullTracks())
    if do_guarded_step(M, o)
        return update_state(N, fs, o, options, l)
    else
        update_state(P, fs, o, options, l)
    end
end

##################################################
initial_fncalls(::AbstractSecant) = 2

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

## Secant
## https://en.wikipedia.org/wiki/Secant_method
"""
    Secant()
    Order1()
    Orderφ()

The `Order1()` method is an alias for `Secant`. It specifies the
[secant method](https://en.wikipedia.org/wiki/Secant_method).
This method keeps two values in its state, `xₙ` and `xₙ₋₁`. The
updated point is the intersection point of ``x`` axis with the secant line
formed from the two points. The secant method uses ``1`` function
evaluation per step and has order `φ≈ (1+sqrt(5))/2`.

The error, `eᵢ = xᵢ - α`, satisfies
`eᵢ₊₂ = f[xᵢ₊₁,xᵢ,α] / f[xᵢ₊₁,xᵢ] * (xᵢ₊₁-α) * (xᵢ - α)`.

"""
struct Secant <: AbstractSecant end
const Order1 = Secant
const Orderφ = Secant

# init_state(M,F,x) --> call init_state(M,F,x₀,x₁,fx₀, fx₁)
function init_state(M::AbstractSecant, F::Callable_Function, x)
    x₀,x₁ = x₀x₁(x)
    fx₀, fx₁ = first(F(x₀)), first(F(x₁))
    state = init_state(M, F, x₀, x₁, fx₀, fx₁)
end

# initialize from xs, fxs
function init_state(::AbstractSecant, F, x₀, x₁, fx₀, fx₁)
    UnivariateZeroState(x₁, x₀, fx₁, fx₀)
end

function update_state(::Order1, F, o::AbstractUnivariateZeroState, options, l=NullTracks())
    xn0, xn1 = o.xn0, o.xn1
    fxn0, fxn1 = o.fxn0, o.fxn1
    Δ = fxn1 * (xn1 - xn0) / (fxn1 - fxn0)

    if isissue(Δ)
        log_message(l, "Increment `Δx` has issues. ")
        return o, true
     end

    x0,x1 = xn1, xn1 - Δ
    fx0, fx1 = fxn1, F(x1)
    incfn(l)

    @set! o.xn0 = x0
    @set! o.xn1 = x1
    @set! o.fxn0 = fx0
    @set! o.fxn1 = fx1

    return o, false

end

### Steffensen
## https://en.wikipedia.org/wiki/Steffensen's_method#Simple_description

"""
    Steffensen()
    Order2()

The quadratically converging
[Steffensen](https://en.wikipedia.org/wiki/Steffensen's_method#Simple_description)
method is used for the derivative-free `Order2()` algorithm. Unlike
the quadratically converging Newton's method, no derivative is
necessary, though like Newton's method, two function calls per step
are. Steffensen's algorithm is more sensitive than Newton's method to
poor initial guesses when `f(x)` is large, due to how `f'(x)` is
approximated. The `Order2` method replaces a Steffensen step with a secant
step when `f(x)` is large.

The error, `eᵢ - α`, satisfies
`eᵢ₊₁ = f[xᵢ, xᵢ+fᵢ, α] / f[xᵢ,xᵢ+fᵢ] ⋅ (1 - f[xᵢ,α] ⋅ eᵢ²`
"""
struct Steffensen <: AbstractSecant end
struct Order2 <: AbstractSecant end

function update_state(method::Order2, fs, o::UnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}
     update_state_guarded(method, Secant(), Steffensen(), fs, o, options, l)
end

function update_state(method::Steffensen, F, o::AbstractUnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}

    x0, x1 = o.xn0, o.xn1
    fx0, fx1 = o.fxn0, o.fxn1

    sgn = sign((fx1 - fx0) / (x1 - x0))
    x2 = x1 - sgn * fx1 / oneunit(S) * oneunit(T)

    f0 = fx1
    f1::S = F(x2)
    incfn(l, 1)

    delta = -sgn * f0 * f0 / (f1 - f0) * oneunit(T) / oneunit(S)

    if isissue(delta)
        log_message(l, "Increment `Δx` has issues. ")
        return o, true
    end

    x0, x1 = x1, x1 - delta
    fx0, fx1 = fx1, F(x1)
    incfn(l)

    @set! o.xn0 = x0
    @set! o.xn1 = x1
    @set! o.fxn0 = fx0
    @set! o.fxn1 = fx1

    return o, false
end

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

"""
    Order5()
    KumarSinghAkanksha()

Implements an order 5 algorithm from *A New Fifth Order Derivative
Free Newton-Type Method for Solving Nonlinear Equations* by Manoj
Kumar, Akhilesh Kumar Singh, and Akanksha, Appl. Math. Inf. Sci. 9,
No. 3, 1507-1513 (2015), DOI: [10.12785/amis/090346](https://doi.org/10.12785/amis/090346). Four function
calls per step are needed.  The `Order5` method replaces a Steffensen step with a secant
step when `f(x)` is large.

The error, `eᵢ = xᵢ - α`, satisfies
`eᵢ₊₁ = K₁ ⋅ K₅ ⋅ M ⋅ eᵢ⁵ + O(eᵢ⁶)`

"""
struct Order5 <: AbstractSecant end
struct KumarSinghAkanksha <: AbstractSecant end

function update_state(method::Order5, fs, o::UnivariateZeroState, options, l=NullTracks())
    update_state_guarded(method, Secant(), KumarSinghAkanksha(), fs, o, options, l)
end

function update_state(M::KumarSinghAkanksha, F, o::AbstractUnivariateZeroState{T,S},
                      options, l=NullTracks()) where {T, S}

    xn = o.xn1
    fxn = o.fxn1

    wn::T = steff_step(M, o.xn1, o.fxn1)

    fwn::S = F(wn)
    incfn(l)

    fp, issue = _fbracket(o.xn1, wn, o.fxn1, fwn)
    if issue
        log_message(l, "Issue with divided difference f[xn, wn]. ")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = wn
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fwn

        return o, true
    end

    yn::T = o.xn1 - o.fxn1 / fp
    fyn::S = F(yn)
    incfn(l)

    zn::T = xn - (fxn + fyn) / fp
    fzn::S = F(zn)
    incfn(l)

    fp, issue = _fbracket_ratio(yn, o.xn1, wn, fyn, o.fxn1, fwn)
    if issue
        log_message(l, "Issue with f[xn,yn]*f[yn,wn] / f[xn, wn]")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = yn
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fyn

        return o, true

    end

    @set! o.xn0 = o.xn1
    @set! o.fxn0 = o.fxn1
    @set! o.xn1 = zn  - fzn  / fp
    @set! o.fxn1 = F(o.xn1)
    incfn(l)

    return o, false

#    nothing
end

struct Order5Derivative <: AbstractSecant end
fn_argout(::Order5Derivative) = 2
function update_state(method::Order5Derivative, f,
                      o::AbstractUnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}

    xn, fxn = o.xn1, o.fxn1
    a::T, b::S = f(xn)
    fpxn = a/b
    incfn(l)

    if isissue(fpxn)
        return o, true
    end

    yn::T = xn - fxn / fpxn
    fyn::S, Δyn::T = f(yn)
    fpyn = fyn / Δyn
    incfn(l, 2)

    if isissue(fpyn)
        log_message(l, "Issue computing `fpyn`")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = yn
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fyn

        return o, true
    end

    zn::T = xn  - (fxn + fyn) / fpxn
    fzn::S, _ = f(zn)
    incfn(l, 2)

    xn1 = zn - fzn / fpyn
    fxn1,_ = f(xn1)
    incfn(l, 2)

    @set! o.xn0 = xn
    @set! o.xn1 = xn1
    @set! o.fxn0 = fxn
    @set! o.fxn1 = fxn1

    return o

end

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

## cf also: https://doi.org/10.1515/tmj-2017-0049
"""
    Order8()
    Thukral8()

Implements an eighth-order algorithm from *New Eighth-Order
Derivative-Free Methods for Solving Nonlinear Equations* by Rajinder
Thukral, International Journal of Mathematics and Mathematical
Sciences Volume 2012 (2012), Article ID 493456, 12 pages DOI:
[10.1155/2012/493456](https://doi.org/10.1155/2012/493456). Four
function calls per step are required.  The `Order8` method replaces a
Steffensen step with a secant step when `f(x)` is large.

The error, `eᵢ = xᵢ - α`, is expressed as `eᵢ₊₁ = K ⋅ eᵢ⁸` in
(2.25) of the paper for an explicit `K`.

"""
struct Order8 <: AbstractSecant end
struct Thukral8 <: AbstractSecant end

function update_state(method::Order8, fs, o::UnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}
    update_state_guarded(method, Secant(), Thukral8(), fs, o, options, l)
end

function update_state(M::Thukral8, F, o::AbstractUnivariateZeroState{T,S}, options, l=NullTracks()) where {T, S}

    xn = o.xn1
    fxn = o.fxn1

    wn::T = steff_step(M, xn, fxn)
    fwn::S = F(wn)
    incfn(l)

    if isissue(fwn)
        log_message(l, "issue with Steffensen step fwn")
        @set! o.xn0 = xn
        @set! o.xn1 = wn
        @set! o.fxn0 = fxn
        @set! o.fxn1 = fwn

        return o, true
    end

    fp, issue = _fbracket(xn, wn, fxn, fwn)

    if issue
        log_message(l, "issue with divided difference f[xn, wn]. ")
        return o, true
    end

    yn::T = xn - fxn / fp
    fyn::S = F(yn)
    incfn(l)

    fp, issue = _fbracket(yn, xn, fyn, fxn)
    if issue #fp
        log_message(l, "issue with divided difference f[xn, yn]. ")
        @set! o.xn0 = xn
        @set! o.xn1 = yn
        @set! o.fxn0 = fxn
        @set! o.fxn1 = fyn

        return o, true
    end

    phi = (1 + fyn / fwn)           # pick one of options
    zn =  yn - phi * fyn / fp
    fzn::S = F(zn)
    incfn(l)

    fp, issue =  _fbracket_diff(xn, yn, zn, fxn, fyn, fzn)
    if issue
        log_message(l, "issue with divided difference  f[y,z] - f[x,y] + f[x,z]. ")
        @set! o.xn0 = xn
        @set! o.xn1 = zn
        @set! o.fxn0 = fxn
        @set! o.fxn1 = fzn

        return o, true
    end

    w = 1 / (1 - fzn/fwn)

    xi = (1 - 2fyn*fyn*fyn / (fwn * fwn * fxn))

    xn1 = zn - w * xi * fzn / fp
    fxn1::S = F(xn1)
    incfn(l)

    @set! o.xn0 = xn
    @set! o.xn1 = xn1
    @set! o.fxn0 = fxn
    @set! o.fxn1 = fxn1

    return o, false

end

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

"""
    Order16()
    Thukral16()

Implements the order 16 algorithm from
*New Sixteenth-Order Derivative-Free Methods for Solving Nonlinear Equations*
by R. Thukral,
American Journal of Computational and Applied Mathematics
p-ISSN: 2165-8935;    e-ISSN: 2165-8943; 2012;  2(3): 112-118
DOI: [10.5923/j.ajcam.20120203.08](https://doi.org/10.5923/j.ajcam.20120203.08).

Five function calls per step are required. Though rapidly converging,
this method generally isn't faster (fewer function calls/steps) over
other methods when using `Float64` values, but may be useful for
solving over `BigFloat`.  The `Order16` method replaces a Steffensen step with a secant
step when `f(x)` is large.

The error, `eᵢ = xᵢ - α`, is expressed as `eᵢ₊₁ = K⋅eᵢ¹⁶` for an explicit `K`
in equation (50) of the paper.

"""
struct Order16 <: AbstractSecant end
struct Thukral16 <: AbstractSecant end

function update_state(method::Order16, fs, o::UnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}
    update_state_guarded(method, Secant(), Thukral16(), fs, o, options, l)
end

function update_state(M::Thukral16, F, o::AbstractUnivariateZeroState{T,S}, options, l=NullTracks()) where {T, S}
    xn = o.xn1
    fxn = o.fxn1

    wn::T = steff_step(M, xn, fxn)
    fwn::S = F(wn)
    incfn(l)

    fp, issue = _fbracket(xn, wn, fxn, fwn)

    if issue
        log_message(l, "issue with f[xn,wn]")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = wn
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fwn
        return o, true
    end

    yn::T = xn - fxn / fp
    fyn::S = F(yn)
    incfn(l)

    fp, issue = _fbracket_ratio(yn, xn, wn, fyn, fxn, fwn)
    if issue
        log_message(l, "issue with f[xn,yn]*f[yn,wn]/f[xn,wn]")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = yn
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fyn

        return o, true
    end

    zn = yn - fyn / fp
    fzn::S = F(zn)
    incfn(l)

    fp, issue = _fbracket_diff(xn, yn, zn, fxn, fyn, fzn)
    u2, u3, u4 = fzn/fwn, fyn/fxn, fyn/fwn

    eta = 1 / (1 + 2*u3*u4^2) / (1 - u2)
    if issue
        log_message(l, "Approximate derivative failed")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = zn
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fzn

        return o, true
    end

    an = zn - eta * fzn / fp
    fan::S = F(an)
    incfn(l)

    fp, issue = _fbracket_ratio(an, yn, zn, fan, fyn, fzn)
    if issue
        log_message(l, "Approximate derivative failed")
        @set! o.xn0 = o.xn1
        @set! o.xn1 = an
        @set! o.fxn0 = o.fxn1
        @set! o.fxn1 = fan

        return o, true
    end

    u1, u5, u6 = fzn/fxn, fan/fxn, fan/fwn

    fp1, issue = _fbracket(xn,yn, fxn, fyn)

    sigma =  1 + u1*u2 - u1*u3*u4^2 + u5 + u6 + u1^2*u4 +
    u2^2*u3 + 3*u1*u4^2*(u3^2 - u4^2)/(fp1/oneunit(fp1))

    xn1 = an - sigma * fan / fp
    fxn1::S = F(xn1)
    incfn(l)

    @set! o.xn0 = xn
    @set! o.xn1 = xn1
    @set! o.fxn0 = fxn
    @set! o.fxn1 = fxn1

    return o, false
end

##################################################
## some means of guarding against large fx when taking a secant step
## TODO: rework this
function steff_step(M::Union{Order5, Order8, Order16}, x::S, fx::T) where {S, T}

    xbar, fxbar = real(x/oneunit(x)), fx/oneunit(fx)
    thresh =  max(1, abs(xbar)) * sqrt(eps(one(xbar))) #^(1/2) # max(1, sqrt(abs(x/fx))) * 1e-6

    out = abs(fxbar) <= thresh ? fxbar  : sign(fx) * thresh
    x + out * oneunit(x)

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export circulargenomeview

"""
    circulargenomeview(;kwargs...)

A CircularGenomeView component.
CircularGenomeView renders the JBrowse 2 React Circular Genome View.
Any JB2 CGV configuration is also valid configuration for this
component.
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `aggregateTextSearchAdapters` (Array; optional): The text search adapters to use for the browser
- `assembly` (optional): The configuration for the assembly to use for the browser.. assembly has the following type: lists containing elements 'name', 'aliases', 'sequence', 'refNameAliases'.
Those elements have the following types:
  - `name` (String; required)
  - `aliases` (Array; optional)
  - `sequence` (Dict; optional)
  - `refNameAliases` (Dict; optional)
- `configuration` (optional): The theme configuration object. configuration has the following type: lists containing elements 'theme'.
Those elements have the following types:
  - `theme` (Dict; optional)
- `defaultSession` (optional): The default session to use for the browser.. defaultSession has the following type: lists containing elements 'name', 'view'.
Those elements have the following types:
  - `name` (String; required)
  - `view` (Dict; optional)
- `tracks` (Array; optional): The tracks to use for the browser.
"""
function circulargenomeview(; kwargs...)
        available_props = Symbol[:id, :aggregateTextSearchAdapters, :assembly, :configuration, :defaultSession, :tracks]
        wild_props = Symbol[]
        return Component("circulargenomeview", "CircularGenomeView", "dash_jbrowse", available_props, wild_props; kwargs...)
end


text:
# This script generates a `.hpp` file (C++ header) containing
# an LDPC code stored in compressed sparse column (CSC) format.
# See command line help for how to use it.

using SparseArrays
using LinearAlgebra

# Must be installed (use the Pkg package manager to instantiate the proviced manifest)
using ArgParse
using LDPCStorage  # custom project. Helper methods for reading and writing LDPC matrices to files

header = raw"""
// This is an automatically generated file.
// A sparse LDPC matrix (containing only zeros and ones) is saved in compressed sparse column (CSC) format.
// Since the matrix (and LDPC code) is known at compile time, there is no need to save it separately in a file.
// This significantly blows up the executable size (the memory would still have to be used when saving the matrix).
// The method seems to be reasonably fast (on a standard laptop).

#include <cstdint>
#include <array>

namespace AutogenLDPC {

"""

struct ValidPath
    val::String
    ValidPath(s::AbstractString) = isfile(s) ? new(s) : error("No such file: '$s'")
end

function parse_my_args(args)
    s = ArgParseSettings(
        "Create a C++ code-file containing an LDPC matrix in compressed sparse column (CSC) format.
        The matrix is given as an `.alist` file."
        ; allow_ambiguous_opts=false)

    @add_arg_table! s begin
        "--input_code_path"
            required = true
            arg_type = ValidPath
            help = "Path to .alist or .cscmat file containing the LDPC code.
                For .cscmat files, both qc-exponents as well as raw LDPC storage is supported."
        "--output_path"
            default = "autogen_ldpc_matrix_csc.hpp"
            arg_type = String
            help = "Path to put the automatically generated C++ file containing the LDPC code.
            If the file extension is .cscmat, a compressed sparse column storage of the full matrix
            (actual binary matrix, not QC-exponents) is written instead."
        "--debug_use_only_part"
            help = "If true, only a very small part of the provided matrix is processed."
            action = :store_true
    end

    parsed_args = parse_args(args, s)
    println("Parsed args:")
    for (key, val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end

    return abspath(parsed_args["input_code_path"].val),
            abspath(parsed_args["output_path"]),
            parsed_args["debug_use_only_part"]
end

function write_cpp_constexpr_CSC(
    H::AbstractArray{Int8, 2},
    output_cpp_path::AbstractString
    ; only_debug_mode = false
    )
    if only_debug_mode  # make matrix smaller...
        cut_at = minimum([10, size(H)...])
        H = sparse(H[1:ceil(Int, cut_at/2), 1:cut_at])
    else
        H = sparse(H)
    end

    _, _, values = findnz(H)

    all(values .== 1) || throw(ArgumentError("Expected matrix containing only zeros and ones."))

    num_nonzero = length(values)
    if log2(num_nonzero) < 16
        colptr_cpp_type = "std::uint16_t"
    elseif log2(num_nonzero) < 32
        colptr_cpp_type = "std::uint32_t"
    elseif log2(num_nonzero) < 64
        colptr_cpp_type = "std::uint64_t"
    else
        throw(ArgumentError("Input matrix not sparse? Has $num_nonzero entries..."))
    end

    if log2(size(H, 1)) < 16
        row_idx_type = "std::uint16_t"
    else
        row_idx_type = "std::uint32_t"
    end

    open(output_cpp_path, "w+") do f
        print(f, header)

        N = size(H, 2)

        println(f, """
        constexpr std::size_t M = $(size(H, 1));
        constexpr std::size_t N = $(size(H, 2));
        constexpr std::size_t num_nz = $num_nonzero;

        constexpr std::array<$colptr_cpp_type, N + 1> colptr = {""")

        for (i, idx) in enumerate(H.colptr)
            print(f, "0x$(string(idx - 1, base=16))")  # Convert index to base zero
            if i != length(H.colptr)
                print(f, ",")
            end
            if mod(i, 100) == 0
                println(f, "")  # for formatting.
            end
        end
        println(f, "\n};\n")

        println(f, "// ------------------------------------------------------- \n")
        println(f, "constexpr std::array<$row_idx_type, num_nz> row_idx = {")

        for (i, idx) in enumerate(H.rowval)
            print(f, "0x$(string(idx - 1, base=16))")  # Convert index to base zero
            if i != length(H.rowval)
                print(f, ",")
            end
            if mod(i, 100) == 0
                println(f, "")  # for formatting.
            end
        end
        println(f, "\n};\n\n")

        println(f, "} // namespace AutogenLDPC")
    end

end

function load_cscmat_standard_or_qc_exponents(code_path; mode=:AUTO)
    if mode == :AUTO
        header = LDPCStorage.read_file_header(code_path)
        if contains(header, "Quasi cyclic exponents")
            mode = :QC_EXPONENTS
        else
            mode = :UNSTRUCTURED
        end
    end

    if mode == :QC_EXPONENTS
        return load_matrix_from_qc_cscmat_file(code_path)
    elseif mode == :UNSTRUCTURED
        H = load_cscmat(code_path)
        if any(x-> 0 <= x <= 1, H)
            @warn "Failed to obtain binary LDPC matrix from CSCMAT file!"
        end
        return H
    else
        @error "unsupported mode $mode"
    end
end

"""
Converts all `.cscmat` files in the folder given by `input_folder_path` to same name `.cscmat`
files in the `output_folder_path`. Each input file is assumed to contain QC exponents.
The output files store the full binary LDPC matrix, which can be used as input to the C++ programs.

Warning: this may overwrite pre-existing files in the output folder.
"""
function convert_all_to_loadable_cscmat(
    input_folder_path::AbstractString, output_folder_path::AbstractString;
    verbose=false
    )
    isdir(input_folder_path) || error("`$input_folder_path` is not a valid folder path.")
    isdir(output_folder_path) || mkdir(output_folder_path)

    file_paths = readdir(input_folder_path; join=true)
    filter!(p -> LDPCStorage.file_extension(p) == ".cscmat", file_paths)

    for p in file_paths
        H = load_cscmat_standard_or_qc_exponents(p)
        save_to_cscmat(H, joinpath(output_folder_path, basename(p));
            allow_omit_entries_if_only_stored_ones=true,)
    end

    verbose && @info "Converted $(size(file_paths)) files."
    return nothing
end

function main(args)
    code_path, output_path, only_debug_mode = parse_my_args(args)
    if LDPCStorage.file_extension(code_path) == ".alist"
        H = load_alist(code_path)
    elseif LDPCStorage.file_extension(code_path) == ".cscmat"

        H = load_cscmat_standard_or_qc_exponents(code_path)
    else
        @error "Unsupported input file at path $code_path. Expected file extension `.cscmat` or `.alist`."
    end

    if LDPCStorage.file_extension(output_path) == ".hpp"
        write_cpp_constexpr_CSC(H, output_path; only_debug_mode)
    elseif LDPCStorage.file_extension(output_path) == ".cscmat"
        save_to_cscmat(H, output_path; allow_omit_entries_if_only_stored_ones=true,)
    else
        @error "Unsupported input file at path $output_path. Expected file extension `.hpp` or `.cscmat`."
    end

    @info "Saved output file ($(Base.stat(output_path).size) bytes) at '$(output_path)'."
end

main(ARGS)


text:
module FinSight
using Roots: find_zero
using Dates: DateTime, Day

"""
Calculates the net present value of an investment by using a discount rate and a
series of future payments (negative values) and income (positive values).

rate - the rate of discount over the length of one period.
values - arguments representing the payments and income.

From https://support.microsoft.com/en-us/office/npv-function-8672cb67-2576-4d07-b67b-ac28acf2a568
"""
function npv(rate::Real, values::AbstractVector{T}) where T <:Real
    @assert length(values) > 0
    out = 0
    for (i, val) in enumerate(values)
        out += val / (1 + rate)^i
    end
    return out
end
export npv

"""
Returns the internal rate of return for a series of cash flows represented by
the numbers in values. These cash flows do not have to be even, as they would be
for an annuity.

values - An array or a reference to cells that contain numbers for
which you want to calculate the internal rate of return.
rate_guess -  A number that you guess is close to the result of IRR

From https://support.microsoft.com/en-us/office/irr-function-64925eaa-9988-495b-b290-3ad0c163c1bc
"""
function irr(values::AbstractVector{T}, rate_guess::Real = 0.1) where T <:Real
    @assert length(values) > 0
    rate = find_zero(x -> npv(x, values),  rate_guess, )
    return rate
end
export irr

"""
Returns the net present value for a schedule of cash flows that is not
necessarily periodic. To calculate the net present value for a series of cash
flows that is periodic, use the NPV function.
"""
function xnpv(rate::Real, values::AbstractVector{T}, dates::Vector{DateTime}) where T <:Real
    @assert length(values) > 0
    @assert length(values) == length(dates)
    out = 0
    for (i, val) in enumerate(values)
        exponent = Day(dates[i] - first(dates)) / Day(365)
        out += val / (1 + rate)^exponent
    end
    return out
end
export xnpv

"""
Returns the internal rate of return for a schedule of cash flows that is not
necessarily periodic. To calculate the internal rate of return for a series of
periodic cash flows, use the IRR function.
"""
function xirr(values::AbstractVector{T}, dates::Vector{DateTime}, rate_guess::Real = 0.1) where T <:Real
    @assert length(values) > 0
    @assert length(values) == length(dates)
    rate = find_zero(x -> xnpv(x, values, dates),  rate_guess)
    return rate
end
export xirr

end

max_stars_repo_name: davidtaylor-fastlanebi/FinSight

max_issues_repo_name: davidtaylor-fastlanebi/FinSight

max_forks_repo_name: davidtaylor-fastlanebi/FinSight


text:
"""An image, taken though a particular filter band"""
type Image
    # The image height.
    H::Int

    # The image width.
    W::Int

    # An HxW matrix of pixel intensities.
    pixels::Matrix{Float32}

    # The band id (takes on values from 1 to 5).
    b::Int

    # World coordinates
    wcs::WCSTransform

    # The components of the point spread function.
    psf::Vector{PsfComponent}

    # SDSS-specific identifiers. A field is a particular region of the sky.
    # A Camcol is the output of one camera column as part of a Run.
    run_num::Int
    camcol_num::Int
    field_num::Int

    # The background noise in nanomaggies. (varies by position)
    epsilon_mat::Array{Float32, 2}

    # The expected number of photons contributed to this image
    # by a source 1 nanomaggie in brightness. (varies by row)
    iota_vec::Array{Float32, 1}

    # storing a RawPSF here isn't ideal, because it's an SDSS type
    # not a Celeste type
    raw_psf_comp::RawPSF
end

"""
Tiles of pixels that share the same set of
relevant sources (or other calculations).  It contains all the information
necessary to compute the ELBO and derivatives in this patch of sky.

Note that this cannot be of type Image for the purposes of Celeste
because the raw wcs object is a C++ pointer which julia resonably
refuses to parallelize.

Attributes:
- h_range: The h pixel locations of the tile in the original image
- w_range: The w pixel locations of the tile in the original image
- pixels: The pixel values
- remainder: the same as in the Image type.
"""
immutable ImageTile
    b::Int

    h_range::UnitRange{Int}
    w_range::UnitRange{Int}
    pixels::Matrix{Float32}

    epsilon_mat::Matrix{Float32}
    iota_vec::Vector{Float32}
end

"""
Constructs an image tile from an image.

Args:
    - img: The Image to be broken into tiles
    - hh: The tile row index (in 1:number of tile rows)
    - ww: The tile column index (in 1:number of tile columns)
    - tile_width: The width and height of a tile in pixels
"""
function ImageTile(hh::Int, ww::Int, img::Image, tile_width::Int)
    h1 = 1 + (hh - 1) * tile_width
    h2 = min(hh * tile_width, img.H)
    h_range = h1:h2

    w1 = 1 + (ww - 1) * tile_width
    w2 = min(ww * tile_width, img.W)
    w_range = w1:w2

    pixels = img.pixels[h_range, w_range]
    epsilon_mat = img.epsilon_mat[h_range, w_range]
    iota_vec = img.iota_vec[h_range]

    @assert 1 <= img.b <= B
    ImageTile(img.b,
              h_range, w_range,
              pixels, epsilon_mat, iota_vec)
end

"""An image, taken though a particular filter band"""
type TiledImage
    # The image height.
    H::Int

    # The image width.
    W::Int

    # subimages
    tiles::Matrix{ImageTile}

    # all tiles have the same height and width
    tile_width::Int

    # The band id (takes on values from 1 to 5).
    b::Int

    # World coordinates
    wcs::WCSTransform

    # The components of the point spread function.
    psf::Vector{PsfComponent}

    # SDSS-specific identifiers. A field is a particular region of the sky.
    # A Camcol is the output of one camera column as part of a Run.
    run_num::Int
    camcol_num::Int
    field_num::Int

    # storing a RawPSF here isn't ideal, because it's an SDSS type
    # not a Celeste type
    raw_psf_comp::RawPSF
end

function TiledImage(img::Image; tile_width=20)
    WW = ceil(Int, img.W / tile_width)
    HH = ceil(Int, img.H / tile_width)
    tiles = ImageTile[ImageTile(hh, ww, img, tile_width) for hh=1:HH, ww=1:WW]
    TiledImage(img.H, img.W, tiles, tile_width, img.b, img.wcs, img.psf,
               img.run_num, img.camcol_num, img.field_num,
               img.raw_psf_comp)
end

""" Returns the tile containing (or nearest to) the specified pixel"""
function get_containing_tile(pixel_crds::Vector{Float64}, img::TiledImage)
    ht0 = ceil(Int, pixel_crds[1] / img.tile_width)
    wt0 = ceil(Int, pixel_crds[2] / img.tile_width)
    ht = max(1, min(size(img.tiles, 1), ht0))
    wt = max(1, min(size(img.tiles, 2), wt0))
    img.tiles[ht, wt]
end


text:
#
#   Kepler.jl
#
#   Solves Kepler's problem for TwoBody orbits.
#

"""
Solves Kepler's Problem for `orbit` and `Δtᵢ`.
"""
function kepler(r::AbstractVector, v::AbstractVector, μ::Number, Δt::Number = period(semimajor_axis(r,v,μ), μ); tol=1e-6, max_iter=100)

    e, a, i, Ω, ω, ν = keplerian(r, v, μ)
    T = period(a, μ)
    conic_section = conic(e)

    # Guess χ₀
    if conic_section == Hyperbolic
        χ₀ = sign(Δt) * √(-a) * log(ℯ, (-2 * μ / a * Δt) / (r ⋅ v + (sign(Δt) * √(-μ * a) * (1 - norm(r) / a))))
    elseif conic_section == Parabolic
        χ₀ = √(a) * tan(ν / 2)
    else
        Δt = mod(Δt, T)
        χ₀ = √(μ) * Δt / a
    end

    # Iteratively solve for χ
    # TODO: Compare loop vs. recursion performance here.
    # There shouldn't be too large of a difference, since this tends
    # to converge with only a few iterations.
    χₙ, rₙ, ψ, C₂, C₃ = χₖ(χ₀, Δt, r, v, a, μ, tol=tol, max_iter=max_iter)

    # Convert to a Orbit
    f = 1 - χₙ^2 / norm(r) * C₂
    ḟ = √(μ) / (norm(r) * rₙ) * χₙ * (ψ * C₃ - 1)
    g = Δt - (χₙ^3 / √(μ)) * C₃
    ġ = 1 - (χₙ^2 / rₙ) * C₂

    return ((f * r) .+ (g * v), (ḟ * r) .+ (ġ * v))
end

function χₖ(χₙ, Δt, rᵢ₀, vᵢ₀, a, μ; iter=1, tol=1e-14, max_iter=100)

    r₀ = norm(rᵢ₀)
    ψ = upreferred(χₙ^2 / a)

    if ψ > tol
        C₂ = (1 - cos(√(ψ))) /  ψ
        C₃ = (√(ψ) - sin(√(ψ))) / √(ψ^3)
    elseif ψ < -tol
        C₂ = (1 - cosh(√(-ψ))) / ψ
        C₃ = (sinh(√(-ψ)) - √(-ψ)) / √((-ψ)^3)
    else
        C₂ = 1.0 / 2.0
        C₃ = 1.0 / 6.0
    end

    r = χₙ^2 * C₂ + (rᵢ₀ ⋅ vᵢ₀) * χₙ / √(μ) * (1 - ψ*C₃) + r₀ * (1 - ψ * C₂)
    χₙ₊₁ = χₙ + ((√(μ) * Δt - χₙ^3 * C₃ - (rᵢ₀ ⋅ vᵢ₀) / √(μ) * χₙ^2 * C₂ - r₀ * χₙ * (1 - ψ * C₃)) / r)

    if iter > max_iter
        @error "Failed to converge!"
        return χₙ, r, ψ, C₂, C₃
    elseif abs(χₙ₊₁ - χₙ) < oneunit(χₙ) * tol
        return  χₙ, r, ψ, C₂, C₃
    else
        return χₖ(χₙ₊₁, Δt, rᵢ₀, vᵢ₀, a, μ; iter=iter+1, tol=tol, max_iter=max_iter)
    end

end