text:
function handle_macro(@nospecialize(x), state) end
function handle_macro(x::EXPR, state)
    !CSTParser.ismacrocall(x) && return
    if headof(x.args[1]) === :globalrefdoc
        if length(x.args) == 4
            if isidentifier(x.args[4]) && !resolve_ref(x.args[4], state)
                if state isa Toplevel
                    push!(state.resolveonly, x)
                end
            elseif CSTParser.is_func_call(x.args[4])
                sig = (x.args[4])
                if sig isa EXPR
                    hasscope(sig) && return # We've already done this, don't repeat
                    setscope!(sig, Scope(sig))
                    mark_sig_args!(sig)
                end
                if state isa Toplevel
                    push!(state.resolveonly, x)
                end
            end
        end
    elseif CSTParser.ismacroname(x.args[1])
        state(x.args[1])
        if _points_to_Base_macro(x.args[1], Symbol("@deprecate"), state) && length(x.args) == 4
            if bindingof(x.args[3]) !== nothing
                return
            elseif CSTParser.is_func_call(x.args[3])
                # add deprecated method
                # add deprecated function binding and args in new scope
                mark_binding!(x.args[3], x)
                mark_sig_args!(x.args[3])
                s0 = state.scope # store previous scope
                state.scope = Scope(s0, x, Dict(), nothing, nothing)
                setscope!(x, state.scope) # tag new scope to generating expression
                state(x.args[3])
                state(x.args[4])
                state.scope = s0
            elseif isidentifier(x.args[3])
                mark_binding!(x.args[3], x)
            end
        elseif _points_to_Base_macro(x.args[1], Symbol("@deprecate_binding"), state) && length(x.args) == 4 && isidentifier(x.args[3]) && isidentifier(x.args[4])
            setref!(x.args[3], refof(x.args[4]))
        elseif _points_to_Base_macro(x.args[1], Symbol("@eval"), state) && length(x.args) == 3 && state isa Toplevel
            # Create scope around eval'ed expression. This ensures anybindings are
            # correctly hoisted to the top-level scope.
            setscope!(x, Scope(x))
            setparent!(scopeof(x), state.scope)
            s0 = state.scope
            state.scope = scopeof(x)
            interpret_eval(x.args[3], state)
            state.scope = s0
        elseif _points_to_Base_macro(x.args[1], Symbol("@irrational"), state) && length(x.args) == 5
            mark_binding!(x.args[3], x)
        elseif _points_to_Base_macro(x.args[1], Symbol("@enum"), state)
            for i = 3:length(x.args)
                if bindingof(x.args[i]) !== nothing
                    break
                end
                if i == 4 && headof(x.args[4]) === :block
                    for j in 1:length(x.args[4].args)
                        mark_binding!(x.args[4].args[j], x)
                    end
                    break
                end
                mark_binding!(x.args[i], x)
            end
        elseif _points_to_Base_macro(x.args[1], Symbol("@goto"), state)
            if length(x.args) == 3 && isidentifier(x.args[3])
                setref!(x.args[3], Binding(noname, nothing, nothing, EXPR[]))
            end
        elseif _points_to_Base_macro(x.args[1], Symbol("@label"), state)
            if length(x.args) == 3 && isidentifier(x.args[3])
                mark_binding!(x.args[3])
            end
        elseif _points_to_Base_macro(x.args[1], Symbol("@NamedTuple"), state) && length(x.args) > 2 && headof(x.args[3]) == :braces
            for a in x.args[3].args
                if CSTParser.isdeclaration(a) && isidentifier(a.args[1]) && !hasref(a.args[1])
                    setref!(a.args[1], Binding(noname, nothing, nothing, EXPR[]))
                end
            end
        elseif is_nospecialize(x.args[1])
            for i = 2:length(x.args)
                if bindingof(x.args[i]) !== nothing
                    break
                end
                mark_binding!(x.args[i], x)
            end
        # elseif _points_to_arbitrary_macro(x.args[1], :Turing, :model, state) && length(x) == 3 &&
        #     isassignment(x.args[3]) &&
        #     headof(x.args[3].args[2]) === CSTParser.Begin && length(x.args[3].args[2]) == 3 && headof(x.args[3].args[2].args[2]) === :block
        #     for i = 1:length(x.args[3].args[2].args[2])
        #         ex = x.args[3].args[2].args[2].args[i]
        #         if isbinarycall(ex, "~")
        #             mark_binding!(ex)
        #         end
        #     end
        # elseif _points_to_arbitrary_macro(x.args[1], :JuMP, :variable, state)
        #     if length(x.args) < 3
        #         return
        #     elseif length(x) >= 5 && ispunctuation(x[2])
        #         _mark_JuMP_binding(x[5])
        #     else
        #         _mark_JuMP_binding(x[3])
        #     end
        # elseif (_points_to_arbitrary_macro(x[1], :JuMP, :expression, state) ||
        #     _points_to_arbitrary_macro(x[1], :JuMP, :NLexpression, state) ||
        #     _points_to_arbitrary_macro(x[1], :JuMP, :constraint, state) || _points_to_arbitrary_macro(x[1], :JuMP, :NLconstraint, state)) && length(x) > 1
        #     if ispunctuation(x[2])
        #         if length(x) == 8
        #             _mark_JuMP_binding(x[5])
        #         end
        #     else
        #         if length(x) == 4
        #             _mark_JuMP_binding(x[3])
        #         end
        #     end
        end
    end
end

function _rem_ref(x::EXPR)
    if headof(x) === :ref && length(x.args) > 0
        return x.args[1]
    end
    return x
end

is_nospecialize(x) = isidentifier(x) && valofid(x) == "@nospecialize"

function _mark_JuMP_binding(arg)
    if isidentifier(arg) || headof(arg) === :ref
        mark_binding!(_rem_ref(arg))
    elseif isbinarycall(arg, "==") || isbinarycall(arg, "<=")  || isbinarycall(arg, ">=")
        if isidentifier(arg.args[1]) || headof(arg.args[1]) === :ref
            mark_binding!(_rem_ref(arg.args[1]))
        else
            mark_binding!(_rem_ref(arg.args[3]))
        end
    elseif headof(arg) === :comparision && length(arg.args) == 5
        mark_binding!(_rem_ref(arg.args[3]))
    end
end

function _points_to_Base_macro(x::EXPR, name, state)
    CSTParser.is_getfield_w_quotenode(x) && return _points_to_Base_macro(x.args[2].args[1], name, state)
    haskey(getsymbols(state)[:Base], name) || return false
    targetmacro =  maybe_lookup(getsymbols(state)[:Base][name], state)
    isidentifier(x) && Symbol(valofid(x)) == name && (ref = refof(x)) !== nothing &&
    (ref == targetmacro || (ref isa Binding && ref.val == targetmacro))
end

function _points_to_arbitrary_macro(x::EXPR, module_name, name, state)
    length(x.args) == 2 && isidentifier(x.args[2]) && valof(x.args[2]) == name && haskey(getsymbols(state), Symbol(module_name)) && haskey(getsymbols(state)[Symbol(module_name)], Symbol("@", name)) && (refof(x.args[2]) == maybe_lookup(getsymbols(state)[Symbol(module_name)][Symbol("@", name)], state) ||
    (refof(x.args[2]) isa Binding && refof(x.args[2]).val == maybe_lookup(getsymbols(state)[Symbol(module_name)][Symbol("@", name)], state)))
end

maybe_lookup(x, env::ExternalEnv) = x isa SymbolServer.VarRef ? SymbolServer._lookup(x, getsymbols(env), true) : x
maybe_lookup(x, state::State) = maybe_lookup(x, state.env)

function maybe_eventually_get_id(x::EXPR)
    if isidentifier(x)
        return x
    elseif isbracketed(x)
        return maybe_eventually_get_id(x.args[1])
    end
    return nothing
end

is_eventually_interpolated(x::EXPR) = isbracketed(x) ? is_eventually_interpolated(x.args[1]) : isunarysyntax(x) && valof(headof(x)) == "\$"
isquoted(x::EXPR) = headof(x) === :quotenode && hastrivia(x) && isoperator(x.trivia[1]) && valof(x.trivia[1]) == ":"
maybeget_quotedsymbol(x::EXPR) = isquoted(x) ? maybe_eventually_get_id(x.args[1]) : nothing

function is_loop_iterator(x::EXPR)
    CSTParser.is_range(x) &&
    ((parentof(x) isa EXPR && headof(parentof(x)) === :for) ||
    (parentof(x) isa EXPR && parentof(parentof(x)) isa EXPR && headof(parentof(parentof(x))) === :for))
end

"""
    maybe_quoted_list(x::EXPR)

Try and get a list of quoted symbols from x. Return nothing if not possible.
"""
function maybe_quoted_list(x::EXPR)
    names = EXPR[]
    if headof(x) === :vect || headof(x) === :tuple
        for i = 1:length(x.args)
            name = maybeget_quotedsymbol(x.args[i])
            if name !== nothing
                push!(names, name)
            else
                return nothing
            end
        end
        return names
    end
end

"""
interpret_eval(x::EXPR, state)

Naive attempt to interpret `x` as though it has been eval'ed. Lifts
any bindings made within the scope of `x` to the toplevel and replaces
(some) interpolated binding names with the value where possible.
"""
function interpret_eval(x::EXPR, state)
    # make sure we have bindings etc
    state(x)
    tls = retrieve_toplevel_scope(x)
    for ex in collect_expr_with_bindings(x)
        b = bindingof(ex)
        if isidentifier(b.name)
            # The name of the binding is fixed
            add_binding(ex, state, tls)
        elseif isunarysyntax(b.name) && valof(headof(b.name)) == "\$"
            # The name of the binding is variable, we need to work out what the
            # interpolated symbol points to.
            variable_name = b.name.args[1]
            resolve_ref(variable_name, state.scope, state)
            if (ref = refof(variable_name)) isa Binding
                if isassignment(ref.val) && (rhs = maybeget_quotedsymbol(ref.val.args[2])) !== nothing
                    # `name = :something`
                    toplevel_binding = Binding(rhs, b.val, nothing, [])
                    settype!(toplevel_binding, b.type)
                    infer_type(toplevel_binding, tls, state)
                    if scopehasbinding(tls, valofid(toplevel_binding.name))
                        tls.names[valofid(toplevel_binding.name)] = toplevel_binding # TODO: do we need to check whether this adds a method?
                    else
                        tls.names[valofid(toplevel_binding.name)] = toplevel_binding
                    end
                elseif is_loop_iterator(ref.val) && (names = maybe_quoted_list(rhs_of_iterator(ref.val))) !== nothing
                    # name is of a collection of quoted symbols
                    for name in names
                        toplevel_binding = Binding(name, b.val, nothing, [])
                        settype!(toplevel_binding, b.type)
                        infer_type(toplevel_binding, tls, state)
                        if scopehasbinding(tls, valofid(toplevel_binding.name))
                            tls.names[valofid(toplevel_binding.name)] = toplevel_binding # TODO: do we need to check whether this adds a method?
                        else
                            tls.names[valofid(toplevel_binding.name)] = toplevel_binding
                        end
                    end
                end
            end
        end
    end
end

function rhs_of_iterator(x::EXPR)
    if isassignment(x)
        x.args[2]
    else
        x.args[3]
    end
end

function collect_expr_with_bindings(x, bound_exprs=EXPR[])
    if hasbinding(x)
        push!(bound_exprs, x)
        # Assuming here that if an expression has a binding we don't want anything bound to chlid nodes.
    elseif x.args !== nothing && !((CSTParser.defines_function(x) && !is_eventually_interpolated(x.args[1])) || CSTParser.defines_macro(x) || headof(x) === :export)
        for a in x.args
            collect_expr_with_bindings(a, bound_exprs)
        end
    end
    return bound_exprs
end


text:
"""
    create_surfaces_neuron(filename, setup)

Create neuron surface mesh.
A neuron surface mesh is loaded or create and loaded. An ECS can be added.
"""
function create_surfaces_neuron(filename, setup::Setup)

    ecs_shape = setup.ecs_shape
    ecs_ratio = setup.ecs_ratio

    if !isfile(filename * "_elements.txt") || !isfile(filename * "_nodes.txt")
        gmesh2fem(filename)
    end
    println("Reading from neuron FE mesh from " * filename)

    # Read points
    nline = open(filename * "_nodes.txt", "r") do io
        countlines(io)
    end
    points = open(filename * "_nodes.txt", "r") do io
        points = zeros(3, nline)
        for iline = 1:nline
            points[:, iline] = parse.(Float64, split(readline(io)))
        end
        points
    end

    # Read elements
    nline = open(filename * "_elements.txt", "r") do io
        countlines(io)
    end
    elements = open(filename * "_elements.txt", "r") do io
        elements = zeros(Int, 4, nline)
        for iline = 1:nline
            elements[:, iline] = parse.(Int, split(readline(io)))
        end
        elements
    end

    # Extract boundary facets from elements (remove interior facets)
    facets = [elements[[1, 2, 3], :] elements[[1, 2, 4], :] elements[[1, 3, 4], :] elements[
        [2, 3, 4],
        :,
    ]]
    sort!(facets, dims = 1)
    facets_unq = unique(facets, dims = 2)
    ikeep = [count(f -> f == fu, eachcol(facets)) == 1 for fu ∈ eachcol(facets_unq)]
    facets = facets_unq[:, ikeep]

    volumes, centers = get_mesh_volumes(points, elements)
    centermass = centers * volumes / sum(volumes)

    facetmarkers = ones(Int, size(facets, 2))
    regions = centermass

    pmin = minimum(points, dims = 2)
    pmax = maximum(points, dims = 2)

    if ecs_shape != "no_ecs"

        # ecs_gap = ecs_ratio * min(pmax - pmin)
        ecs_gap = ecs_ratio * 10

        if ecs_shape == "box"
            points_ecs = [
                pmin(1) pmax(1) pmax(1) pmin(1) pmin(1) pmax(1) pmax(1) pmin(1)
                pmin(2) pmin(2) pmax(2) pmax(2) pmin(2) pmin(2) pmax(2) pmax(2)
                pmin(3) pmin(3) pmin(3) pmin(3) pmax(3) pmax(3) pmax(3) pmax(3)
            ]
            facets_ecs = [
                1 1 1 1 2 2 3 3 4 4 5 5
                2 3 2 6 3 7 4 8 1 5 6 7
                3 4 6 5 7 6 8 7 5 8 7 8
            ]
        elseif ecs_shape == "convex_hull"
            error("Not implemented")
        elseif ecs_shape == "tight_wrap"
            error("Not implemented")
        end

        facetmarkers_ecs = 2 * ones(Int, size(facets_ecs, 2))
        _, ind = sort(points[1, :])
        ind = ind[1]
        regions_ecs = points[:, ind] - ecs_gap / 10 * [1; 0; 0]

        npoint_out = size(points, 2)
        points = [points points_ecs]
        facets = [facets facets_ecs + npoint_out]
        append!(facetmarkers, facetmarkers_ecs)
        append!(regions, regions_ecs)
    end

    (; points, facets, facetmarkers, regions)
end

ext: jl


text:
# ---
# title: "`bothat`"
# id: op_bothat
# cover: assets/bothat.png
# ---

# The (black) tophat operator, also known as bottom hat, is defined as `closing(img) - img`.
# Intuitively, this filter can be used to extract small black elements and details from an
# image.

using ImageMorphology
using TestImages
using ImageBase
using ImageShow

#-

img = restrict(testimage_dip3e("fig1038")) # fingerprint

out1 = bothat(img) # default: all spatial dimensions, r=1, a box-shape SE
out2 = bothat(img; dims=(2,)) # only apply to the second dimension
out3 = bothat(img; r=5) # half-size r=5
out4 = bothat(img, strel_diamond((5, 5))) # generic SE -- this version doesn't accept `r` and `dims`

mosaic(img, out1, out2, out3, out4; nrow=1)

# ## See also

# To extract white small details, use the white [`tophat`](@ref op_tophat):

tophat(img) == bothat(complement.(img))

# For a comprehensive and more accurate documentation, please check the [`bothat`](@ref)
# reference page.

# save cover image #src
using FileIO #src
mkpath("assets") #src
img = Gray.(restrict(testimage("blobs"))) #src
save("assets/bothat.png", bothat(img; r=2)) #src


text:
module Optimizer

using ..ADCME

using LinearAlgebra
using Libdl
#=
The following code is adapted from Flux.jl: https://github.com/FluxML/Flux.jl
The original license from Flux.jl:

The Flux.jl package is licensed under the MIT "Expat" License:
Copyright (c) 2016-19: Julia Computing, INc., Mike Innes and Contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

const ϵ = 1e-8

# TODO: should use weak refs

"""
Descent(η = 0.001)

Classic gradient descent optimiser with learning rate `η`.
For each parameter `p` and its gradient `δp`, this runs `p -= η*δp`

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.

# Examples
```julia
opt = Descent()

opt = Descent(0.3)

ps = params(model)

gs = gradient(ps) do
  loss(x, y)
end

Flux.Optimise.update!(opt, ps, gs)
```
"""
mutable struct Descent
  eta::Float64
end

Descent() = Descent(0.001)

function apply!(o::Descent, x, Δ)
  Δ *= o.eta
end

"""
Momentum(η = 0.01, ρ = 0.9)

Gradient descent optimizer with learning rate `η` and momentum `ρ`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Momentum (`ρ`): Controls the acceleration of gradient descent in the
prominent direction, in effect dampening oscillations.

# Examples
```julia
opt = Momentum()

opt = Momentum(0.01, 0.99)
```
"""
mutable struct Momentum
  eta::Float64
  rho::Float64
  velocity::Union{Missing, Array}
end

Momentum(η = 0.01, ρ = 0.9) = Momentum(η, ρ, missing)

function apply!(o::Momentum, x, Δ)
  η, ρ = o.eta, o.rho
  v = coalesce(o.velocity, zero(x))
  @. v = ρ * v - η * Δ
  ismissing(o.velocity) && (o.velocity = v)
  Δ = -v
end

"""
Nesterov(η = 0.001, ρ = 0.9)

Gradient descent optimizer with learning rate `η` and Nesterov momentum `ρ`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Nesterov momentum (`ρ`): Controls the acceleration of gradient descent in the
prominent direction, in effect dampening oscillations.

# Examples
```julia
opt = Nesterov()

opt = Nesterov(0.003, 0.95)
```
"""
mutable struct Nesterov
  eta::Float64
  rho::Float64
  velocity::Union{Array, Missing}
end

Nesterov(η = 0.001, ρ = 0.9) = Nesterov(η, ρ, missing)

function apply!(o::Nesterov, x, Δ)
  η, ρ = o.eta, o.rho
  v = coalesce(o.velocity, zero(x))
  d = @. ρ^2 * v - (1+ρ) * η * Δ
  @. v = ρ*v - η*Δ
  ismissing(o.velocity) && (o.velocity = v)
  Δ = -d
end

"""
RMSProp(η = 0.001, ρ = 0.9)

Optimizer using the
[RMSProp](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
algorithm. Often a good choice for recurrent networks. Parameters other than learning rate
generally don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Momentum (`ρ`): Controls the acceleration of gradient descent in the
prominent direction, in effect dampening oscillations.

# Examples
```julia
opt = RMSProp()

opt = RMSProp(0.002, 0.95)
```
"""
mutable struct RMSProp
  eta::Float64
  rho::Float64
  acc::Union{Missing, Array}
end

RMSProp(η = 0.001, ρ = 0.9) = RMSProp(η, ρ, missing)

function apply!(o::RMSProp, x, Δ)
  η, ρ = o.eta, o.rho
  acc = coalesce(o.acc, zero(x))
  @. acc = ρ * acc + (1 - ρ) * Δ^2
  ismissing(o.acc) && (o.acc = acc)
  Δ = @. Δ * η / (√acc + ϵ)
end

"""
ADAM(η = 0.001, β::Tuple = (0.9, 0.999))

[ADAM](https://arxiv.org/abs/1412.6980v8) optimiser.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
second (β2) momentum estimate.

# Examples
```julia
opt = ADAM()

opt = ADAM(0.001, (0.9, 0.8))
```
"""
mutable struct ADAM
  eta::Float64
  beta::Tuple{Float64,Float64}
  state::Union{Missing, Tuple}
end

ADAM(η = 0.001, β = (0.9, 0.999)) = ADAM(η, β, missing)

function apply!(o::ADAM, x, Δ)
  η, β = o.eta, o.beta
  mt, vt, βp = coalesce(o.state, (zero(x), zero(x), β))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ^2
  Δ =  @. mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + ϵ) * η
  o.state = (mt, vt, βp .* β)
  return Δ
end

"""
RADAM(η = 0.001, β::Tuple = (0.9, 0.999))

[Rectified ADAM](https://arxiv.org/pdf/1908.03265v1.pdf) optimizer.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
second (β2) momentum estimate.

# Examples
```julia
opt = RADAM()

opt = RADAM(0.001, (0.9, 0.8))
```
"""
mutable struct RADAM
  eta::Float64
  beta::Tuple{Float64,Float64}
  state::Union{Tuple, Missing}
end

RADAM(η = 0.001, β = (0.9, 0.999)) = RADAM(η, β, missing)

function apply!(o::RADAM, x, Δ)
  η, β = o.eta, o.beta
  ρ∞ = 2/(1-β[2])-1
  mt, vt, βp, t = coalesce(o.state, (zero(x), zero(x), β, 1))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ^2
  ρ = ρ∞ - 2t*βp[2]/(1-βp[2])
  if ρ > 4
    r = sqrt((ρ-4)*(ρ-2)*ρ∞/((ρ∞-4)*(ρ∞-2)*ρ))
    Δ =  @. mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + ϵ) * η * r
  else
    Δ =  mt / (1 - βp[1]) * η
  end
  o.state = (mt, vt, βp .* β, t+1)
  return Δ
end

"""
AdaMax(η = 0.001, β::Tuple = (0.9, 0.999))

[AdaMax](https://arxiv.org/abs/1412.6980v9) is a variant of ADAM based on the ∞-norm.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
second (β2) momentum estimate.

# Examples
```julia
opt = AdaMax()

opt = AdaMax(0.001, (0.9, 0.995))
```
"""
mutable struct AdaMax
  eta::Float64
  beta::Tuple{Float64,Float64}
  state::Union{Tuple, Missing}
end

AdaMax(η = 0.001, β = (0.9, 0.999)) = AdaMax(η, β, missing)

function apply!(o::AdaMax, x, Δ)
  η, β = o.eta, o.beta
  mt, ut, βp = coalesce(o.state,  (zero(x), zero(x), β))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. ut = max(β[2] * ut, abs(Δ))
  Δ = @. (η/(1 - βp[1])) * mt/(ut + ϵ)
  o.state = (mt, ut, βp .* β)
  return Δ
end

"""
ADAGrad(η = 0.1)

[ADAGrad](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf) optimizer. It has
parameter specific learning rates based on how frequently it is updated.
Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.

# Examples
```julia
opt = ADAGrad()

opt = ADAGrad(0.001)
```
"""
mutable struct ADAGrad
  eta::Float64
  acc::Union{Array, Missing}
end

ADAGrad(η = 0.001) = ADAGrad(η, missing)

function apply!(o::ADAGrad, x, Δ)
  η = o.eta
  acc = coalesce(o.acc, fill!(zero(x), ϵ))
  @. acc += Δ^2
  ismissing(o.acc) && (o.acc = acc)
  Δ = @. Δ * η / (√acc + ϵ)
end

"""
ADADelta(ρ = 0.9)

[ADADelta](https://arxiv.org/abs/1212.5701) is a version of ADAGrad adapting its learning
rate based on a window of past gradient updates.
Parameters don't need tuning.

# Parameters
- Rho (`ρ`): Factor by which the gradient is decayed at each time step.

# Examples
```julia
opt = ADADelta()

opt = ADADelta(0.89)
```
"""
mutable struct ADADelta
  rho::Float64
  state::Union{Tuple, Missing}
end

ADADelta(ρ = 0.9) = ADADelta(ρ, missing)

function apply!(o::ADADelta, x, Δ)
  ρ = o.rho
  acc, Δacc = coalesce(o.state, (zero(x), zero(x)))
  @. acc = ρ * acc + (1 - ρ) * Δ^2
  Δ = @. Δ * √Δacc/ (√acc + ϵ)
  Δacc = @. ρ * Δacc + (1 - ρ) * Δ^2
  if ismissing(o.state)
    o.state = (acc, Δacc)
  end
  return Δ
end

"""
AMSGrad(η = 0.001, β::Tuple = (0.9, 0.999))

The [AMSGrad](https://openreview.net/forum?id=ryQu7f-RZ) version of the ADAM
optimiser. Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
second (β2) momentum estimate.

# Examples
```julia
opt = AMSGrad()

opt = AMSGrad(0.001, (0.89, 0.995))
```
"""
mutable struct AMSGrad
  eta::Float64
  beta::Tuple{Float64, Float64}
  state::Union{Tuple, Missing}
end

AMSGrad(η = 0.001, β = (0.9, 0.999)) = AMSGrad(η, β, missing)

function apply!(o::AMSGrad, x, Δ)
  η, β = o.eta, o.beta
  mt, vt, v̂t = coalesce(o.state, (fill!(zero(x), ϵ), fill!(zero(x), ϵ), fill!(zero(x), ϵ)))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ ^ 2
  @. v̂t = max(v̂t, vt)
  if ismissing(o.state)
    o.state = (mt, vt, v̂t)
  end
  Δ = @. η * mt / (√v̂t + ϵ)
end

"""
NADAM(η = 0.001, β::Tuple = (0.9, 0.999))

[NADAM](http://cs229.stanford.edu/proj2015/054_report.pdf) is a Nesterov variant of ADAM.
Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
second (β2) momentum estimate.

# Examples
```julia
opt = NADAM()

opt = NADAM(0.002, (0.89, 0.995))
```
"""
mutable struct NADAM
  eta::Float64
  beta::Tuple{Float64, Float64}
  state::Union{Tuple, Missing}
end

NADAM(η = 0.001, β = (0.9, 0.999)) = NADAM(η, β, missing)

function apply!(o::NADAM, x, Δ)
  η, β = o.eta, o.beta
  mt, vt, (β1p, β2p) = coalesce(o.state, (zero(x), zero(x), o.beta))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ^2
  Δ = @. (β[1] * mt / (1 - β[1] * β1p) + (1 - β[1]) * Δ / (1 - β1p)) / (√(vt * β[2] / (1 - β2p)) + ϵ) * η
  o.state = (mt, vt, (β1p * β[1], β2p * β[2]))
  return Δ
end

"""
    InvDecay(γ = 0.001)
Apply inverse time decay to an optimiser, so that the effective step size at
iteration `n` is `eta / (1 + γ * n)` where `eta` is the initial step size.
The wrapped optimiser's step size is not modified.
# Examples
```julia
decay = InvDecay()
d = apply!(decay, x0, g)
```
"""
mutable struct InvDecay
  gamma::Float64
  state::Union{Int64, Missing}
end

InvDecay(γ = 0.001) = InvDecay(γ, missing)

function apply!(o::InvDecay, x, Δ)
  γ = o.gamma
  n = coalesce(o.state, 1)
  Δ *= 1 / (1 + γ * n)
  o.state = n + 1
  return Δ
end

"""
    ExpDecay(η = 0.001, decay = 0.1, decay_step = 1000, clip = 1e-4)
Discount the learning rate `η` by the factor `decay` every `decay_step` steps till
a minimum of `clip`.
# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- `decay`: Factor by which the learning rate is discounted.
- `decay_step`: Schedule decay operations by setting the number of steps between
                two decay operations.
- `clip`: Minimum value of learning rate.
# Examples
To apply exponential decay to an optimiser:
```julia
decay = ExpDecay()
d = apply!(decay, x0, g)
```
"""
mutable struct ExpDecay
  eta::Float64
  decay::Float64
  step::Int64
  clip::Float64
  current::Union{Int64, Missing}
end

ExpDecay(opt = 0.001, decay = 0.1, decay_step = 1000, clip = 1e-4) = ExpDecay(opt, decay, decay_step, clip, missing)

function apply!(o::ExpDecay, x, Δ)
  η, s, decay = o.eta, o.step, o.decay
  n = o.current = coalesce(o.current, 0) + 1
  if o.current%s == 0 && o.current÷s == 1
    η = max(η * decay, o.clip)
    o.eta = η
  end
  Δ *= η
  return Δ
end

#=
The following code is adapted from Optim.jl: https://github.com/JuliaNLSolvers/Optim.jl
Original license:

Optim.jl is licensed under the MIT License:
Copyright (c) 2012: John Myles White, Tim Holy, and other contributors. Copyright (c) 2016: Patrick Kofod Mogensen, John Myles White, Tim Holy, and other contributors. Copyright (c) 2017: Patrick Kofod Mogensen, Asbjørn Nilsen Riseth, John Myles White, Tim Holy, and other contributors.
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
© 2020 GitHub, Inc.
=#

function twoloop!(s,
  gr,
  rho,
  dx_history,
  dg_history,
  m::Integer,
  pseudo_iteration::Integer,
  alpha,
  q,
  scaleinvH0::Bool,
  precon)
  # Count number of parameters
  n = length(s)

  # Determine lower and upper bounds for loops
  lower = pseudo_iteration - m
  upper = pseudo_iteration - 1

  # Copy gr into q for backward pass
  copyto!(q, gr)

  # Backward pass
  for index in upper:-1:lower
    if index < 1
      continue
    end
    i   = mod1(index, m)
    dgi = dg_history[i]
    dxi = dx_history[i]
    @inbounds alpha[i] = rho[i] * real(dot(dxi, q))
    @inbounds q .-= alpha[i] .* dgi
  end

  # Copy q into s for forward pass
  # apply preconditioner if precon != nothing
  # (Note: preconditioner update was done outside of this function)
  if scaleinvH0 == true && pseudo_iteration > 1
    # Use the initial scaling guess from
    # Nocedal & Wright (2nd ed), Equation (7.20)

    #=
    pseudo_iteration > 1 prevents this scaling from happening
    at the first iteration, but also at the first step after
    a reset due to invH being non-positive definite (pseudo_iteration = 1).
    TODO: Maybe we can still use the scaling as long as iteration > 1?
    =#
    i = mod1(upper, m)
    dxi = dx_history[i]
    dgi = dg_history[i]
    scaling = real(dot(dxi, dgi)) / sum(abs2, dgi)
    @. s = scaling*q
  else
    if isnothing(precon)
      copyto!(s, q)
    else
      ldiv!(s, precon, q)
    end
  end
  # Forward pass
  for index in lower:1:upper
    if index < 1
      continue
    end
    i = mod1(index, m)
    dgi = dg_history[i]
    dxi = dx_history[i]
    @inbounds beta = rho[i] * real(dot(dgi, s))
    @inbounds s .+= dxi .* (alpha[i] - beta)
  end

  # # Negate search direction
  # rmul!(s, eltype(s)(-1))

  return
end

mutable struct LBFGS
  m::Int
  P
  precondprep!::Function
  scaleinvH0::Bool
  pseudo_iteration::Int64
  rho::Union{Missing, Vector{<:Real}}
  dx_history
  dg_history
  twoloop_alpha
  twoloop_q
  dx::Union{Missing, Array{<:Real}}
  dg::Union{Missing, Array{<:Real}}
  g_previous::Union{Missing, Array{<:Real}}
  x_previous::Union{Missing, Array{<:Real}}
  s::Union{Missing, Array{<:Real}} # current search direction
end

"""
LBFGS(; m::Integer = 10,
P=nothing,
precondprep = (P, x) -> nothing,
scaleinvH0::Bool = true && (typeof(P) <: Nothing) )

Returns a LBFGS optimizer. `P` is the preconditioner for the LBFGS algorithm. When `P` is not nothing,
`precondprep` modifies `P` in place. When `P` is nothing, a scalar matrix is used as the preconditioner.
"""
function LBFGS(; m::Integer = 10,
  P=nothing,
  precondprep = (P, x) -> nothing,
  scaleinvH0::Bool = true && (typeof(P) <: Nothing) )
  LBFGS(Int(m), P, precondprep, scaleinvH0, 0,
  missing, missing, missing, missing, missing, missing, missing, missing, missing, missing)
end

function apply!(o::LBFGS, x, Δ)
  if !ismissing(o.rho)
    update_h!(o, x, Δ)
  end

  # the first time some fields in LBFGS have access to the data
  if ismissing(o.rho)
    T = eltype(x)
    o.rho = Vector{T}(undef, o.m)
    o.dx_history = [similar(x) for i = 1:o.m]
    o.dg_history = [similar(x) for i = 1:o.m]
    o.twoloop_alpha = Vector{T}(undef, o.m)
    o.twoloop_q = similar(x)
    o.dx = similar(x)
    o.dg = similar(x)
    o.g_previous = similar(x)
    o.s = similar(x)
    o.x_previous = copy(x)
  end

  # Increment the number of steps we've had to perform
  o.pseudo_iteration += 1

  # update the preconditioner
  o.precondprep!(o.P, x)

  # Determine the L-BFGS search direction # FIXME just pass state and method?
  twoloop!(o.s, Δ, o.rho, o.dx_history, o.dg_history,
  o.m, o.pseudo_iteration,
  o.twoloop_alpha, o.twoloop_q, o.scaleinvH0, o.P)

  copyto!(o.g_previous, Δ)

  return o.s
end

function update_h!(o::LBFGS, x, d)
  n = length(x)
  # Measure the change in the gradient
  o.dg .= d .- o.g_previous
  o.dx, o.x_previous = x - o.x_previous, copy(x)

  # Update the L-BFGS history of positions and gradients
  rho_iteration = one(eltype(o.dx)) / real(dot(o.dx, o.dg))
  if isinf(rho_iteration)
    # TODO: Introduce a formal error? There was a warning here previously
    error(DomainError("rho is Inf"))
  end
  idx = mod1(o.pseudo_iteration, o.m)
  @inbounds o.dx_history[idx] .= o.dx
  @inbounds o.dg_history[idx] .= o.dg
  @inbounds o.rho[idx] = rho_iteration
end

mutable struct AndersonAcceleration
    mem::Int64
    atype::Int64
    initilized::Bool
end

"""
    AndersonAcceleration(; mem::Int64 = 50, atype::Int64 = 1)

Creates an Anderson accelerator. To apply the Anderson Accelerator, do

```julia
y = apply!(aa, x, Δ)
```

Here `aa` is an AndersonAcceleration instant, `x` is the current state, and `Δ` is the (modified) gradient.

`x` will be modified in place, but is also returned for convenience.
"""
function AndersonAcceleration(; mem::Int64 = 50, atype::Int64 = 1)
  if Sys.iswindows()
    if !isfile(joinpath(ADCME.LIBDIR, "liblapack.lib"))
      PWD = pwd()
      download("https://github.com/kailaix/ADCME.jl/releases/download/blas_lapack/aa_win.zip", joinpath(ADCME.LIBDIR, "aa_win.zip"))
      cd(ADCME.LIBDIR)
      run(`cmd /c unzip aa_win.zip`)
      files = readdir("aa_win")
      for f in files
        mv(joinpath("aa_win", f), f, force=true)
      end
      cd(PWD)
    end
  end
  lib = Sys.iswindows() ? "AA.dll" : "libAA.$dlext"
  pth = joinpath(@__DIR__, "..", "deps", "CustomOps", "AndersonAcceleration", "build", lib)
  if !isfile(pth)
    PWD = pwd()
    cd(joinpath(@__DIR__, "..", "deps", "CustomOps", "AndersonAcceleration"))
    compile()
    cd(PWD)
  end
  AndersonAcceleration(mem, atype, false)
end

function apply!(o::AndersonAcceleration, x, Δ)
  n = length(x)
  x_prev = x
  x -= Δ

  lib = Sys.iswindows() ? "AA.dll" : "libAA.$dlext"
  pth = joinpath(@__DIR__, "..", "deps", "CustomOps", "AndersonAcceleration", "build", lib)
  if !o.initilized
    @eval ccall((:init_aa, $pth), Cvoid, (Cint, Cint, Cint), $n, $o.mem, $o.atype)
  end
  @eval ccall((:apply_aa, $pth), Cvoid, (Ref{Cdouble}, Ref{Cdouble}), $x, $x_prev)
  return x
end

function Base.:finalizer(f, o::AndersonAcceleration)
  lib = Sys.iswindows() ? "AA.dll" : "libAA.$dlext"
  pth = joinpath(@__DIR__, "..", "deps", "CustomOps", "AndersonAcceleration", "build", lib)
  @eval ccall((:finalize_aa, $pth), Cvoid, ())
end

end


text:
export merged

import Base: show

# -------------------- #
# Merge observable     #
# -------------------- #

struct MergeObservable{D, S} <: Subscribable{D}
    sources :: S
end

function on_subscribe!(observable::MergeObservable{D}, actor) where D
    merge_main    = __make_merge_main_actor(D, actor, length(observable.sources))

    subscriptions = map(enumerate(observable.sources)) do (index, source)
        return subscribe!(source, __make_merge_child_actor_factory(index, merge_main))
    end

    return MergeSubscription(subscriptions)
end

# -------------------- #
# Merge main actor     #
# -------------------- #

struct MergeMainActor{D, A} <: Actor{D}
    actor             :: A
    completion_status :: BitArray{1}
end

__make_merge_main_actor(::Type{D}, actor::A, length::Int) where { D, A } = MergeMainActor{D, A}(actor, falses(length))

on_next!(actor::MergeMainActor{D}, data::L) where { D, L <: D } = next!(actor.actor, data)
on_error!(actor::MergeMainActor, err)                           = error!(actor.actor, err)
on_complete!(actor::MergeMainActor)                             = begin
    if all(actor.completion_status)
        complete!(actor.actor)
    end
end

# -------------------- #
# Merge child actor    #
# -------------------- #

struct MergeChildActor{D, I, A} <: Actor{D}
    main :: A
end

on_next!(actor::MergeChildActor{D}, data::D) where D     = next!(actor.main, data)
on_error!(actor::MergeChildActor, err)                   = error!(actor.main, err)
on_complete!(actor::MergeChildActor{D, I}) where { D, I } = begin
    actor.main.completion_status[I] = true
    complete!(actor.main)
end

struct MergeChildActorFactory{I, A} <: AbstractActorFactory
    main :: A
end

__make_merge_child_actor_factory(index::Int, main::A) where A = MergeChildActorFactory{index, A}(main)

create_actor(::Type{L}, factory::MergeChildActorFactory{I, A}) where { L, I, A } = MergeChildActor{L, I, A}(factory.main)

# -------------------- #
# Merge subscription   #
# -------------------- #

struct MergeSubscription{S} <: Teardown
    subscriptions :: S
end

as_teardown(::Type{<:MergeSubscription}) = UnsubscribableTeardownLogic()

on_unsubscribe!(subscription::MergeSubscription) = foreach(s -> unsubscribe!(s), subscription.subscriptions)

"""
    merged(sources::T) where { T <: Tuple }

Creation operator for the `MergeObservable` with a given `sources` collected in a tuple.
`merge` subscribes to each given input Observable (as arguments), and simply forwards (without doing any transformation) all the values from all the input
Observables to the output Observable. The output Observable only completes once all input Observables have completed.
Any error delivered by an input Observable will be immediately emitted on the output Observable.

# Examples

```jldoctest
using Rocket

observable = merged((from(1:4), of(2.0), from("Hello")))

subscribe!(observable, logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 4
[LogActor] Data: 2.0
[LogActor] Data: H
[LogActor] Data: e
[LogActor] Data: l
[LogActor] Data: l
[LogActor] Data: o
[LogActor] Completed
```

```
using Rocket

subject = Subject(Int, scheduler = AsyncScheduler())

observable = merged((subject, of(2.0), from("Hello")))

actor = sync(logger())

subscribe!(observable, actor)

setTimeout(200) do
    next!(subject, 1)
    complete!(subject)
end

wait(actor)
;

# output
[LogActor] Data: 2.0
[LogActor] Data: H
[LogActor] Data: e
[LogActor] Data: l
[LogActor] Data: l
[LogActor] Data: o
[LogActor] Data: 1
[LogActor] Completed
```

See also: [`Subscribable`](@ref)
"""
merged(sources::T) where { T <: Tuple }         = MergeObservable{ Union{ subscribable_extract_type.(sources)... }, T }(sources)
merged(sources::T) where { T <: AbstractArray } = error("Rocket.merge takes a tuple of sources as an argument, not an array")

Base.show(io::IO, ::MergeObservable{D}) where D = print(io, "MergeObservable($D)")
Base.show(io::IO, ::MergeSubscription)          = print(io, "MergeSubscription()")


text:
import CumulantsFeatures: makeblocksize, computeblock

@testset "axiliary for cum2mat" begin
    Random.seed!(42)
    t = rand(SymmetricTensor{Float64, 3}, 9, 4)
    @test makeblocksize(t, (1,1)) == (4,4)
    @test makeblocksize(t, (1,3)) == (4,1)
    @test makeblocksize(t, (3,1)) == (1,4)
    t1 = rand(SymmetricTensor{Float64, 3}, 9, 3)
    @test makeblocksize(t1, (1,1)) == (3,3)
    dims = (fill(t.bln, 1)...,)
    @test computeblock(t, (1,1), dims) ≈ [21.8503 16.2408 20.5519 18.4874; 16.2408 25.1741 18.1145 19.4899; 20.5519 18.1145 31.5073 22.8366; 18.4874 19.4899 22.8366 28.156] atol = 0.01
end

@testset "cum2mat tests on cumulants" begin
    x = rand(20,13)
    c = cumulants(x,5,2)
    s = cum2mat(c[3])
    X = unfold(Array(c[3]), 1)
    M = X*transpose(X)
    @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)

    s = cum2mat(c[4])
    X = unfold(Array(c[4]), 1)
    M = X*transpose(X)
    @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)

    s = cum2mat(c[5])
    X = unfold(Array(c[5]), 1)
    M = X*transpose(X)
    @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
end

@testset "cum2mat tests on random tensor one core" begin
    @testset "order 3" begin
        Random.seed!(42)
        t = rand(SymmetricTensor{Float64, 3}, 12, 4)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
    @testset "order 4" begin
        t = rand(SymmetricTensor{Float64, 4}, 12, 3)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
    @testset "order 5" begin
        t = rand(SymmetricTensor{Float64, 5}, 10, 2)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
    @testset "order 5 not squared" begin
        t = rand(SymmetricTensor{Float64, 5}, 11, 3)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
    @testset "order 6" begin
        t = rand(SymmetricTensor{Float64, 6}, 9, 3)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
    @testset "order 6 not squared" begin
        t = rand(SymmetricTensor{Float64, 6}, 7, 3)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
    @testset "larger data set" begin
        t = rand(SymmetricTensor{Float64, 5}, 26, 4)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-8)
    end
end


text:
using JavaCall
DBtableTypeorString = Union{DBtableType, AbstractString}

function jaccard(A::DBtableTypeorString, ADeg::DBtableTypeorString, Rname::AbstractString; filterRowCol::ValidQueryTypes="")

# Original Java docstring:
# /**
# * From input <b>unweighted, undirected</b> adjacency table Aorig,
# * put the Jaccard coefficients in the upper triangle of Rfinal.
# * @param Aorig Unweighted, undirected adjacency table.
# * @param ADeg Degree table name.
# * @param Rfinal Should not previously exist. Writes the Jaccard table into Rfinal,
# *               using a couple combiner-like iterators.
# * @param filterRowCol Filter applied to rows and columns of Aorig
# *                     (must apply to both rows and cols because A is undirected Adjacency table).
# * @param Aauthorizations Authorizations for scanning Atable. Null means use default: Authorizations.EMPTY
# * @param RNewVisibility Visibility label for new entries created in Rtable. Null means no visibility label.
# * @return number of partial products sent to Rtable during the Jaccard coefficient calculation
# */
# public long Jaccard(String Aorig, String ADeg, String Rfinal,
# String filterRowCol, Authorizations Aauthorizations, String RNewVisibility)

    JAuthorizations = @jimport "org.apache.accumulo.core.security.Authorizations"
    emptyauth = JAuthorizations(  (), )
    newVisibility = "" # TODO default or null

    if isa(A, DBtablePair)
        Aname = A.name1
    elseif isa(A, DBtable)
        Aname = A.name
    else
        Aname = A
    end

    if isa(ADeg, DBtablePair)
        ADegname = ADeg.name1
    elseif isa(ADeg, DBtable)
        ADegname = ADeg.name
    else
        ADegname = ADeg
    end

    DB = A.DB

    return jcall(DB.Graphulo,"Jaccard", jlong,
        (JString,JString,JString,JString,JAuthorizations,JString,),
        Aname,ADegname,Rname,toDBstring(filterRowCol),emptyauth, newVisibility)

end


text:
# Physical grid tests
testdx(g) = abs(sum(g.x[2:end].-g.x[1:end-1]) - (g.nx-1)*g.dx) < 1e-15*(g.nx-1)
testdy(g) = abs(sum(g.y[2:end].-g.y[1:end-1]) - (g.ny-1)*g.dy) < 1e-15*(g.nx-1)

testx(g) = g.x[end]-g.x[1] == g.Lx-g.dx
testy(g) = g.y[end]-g.y[1] == g.Ly-g.dy

# Test proper arrangement of fft wavenumbers
testk(g) = sum(g.k[2:g.nkr-1] .+ flipdim(g.k[g.nkr+1:end], 1)) == 0.0
testkr(g) = sum(g.k[1:g.nkr] .- g.kr) == 0.0
testl(g) = sum(g.l[:, 2:Int(g.ny/2)] .+ flipdim(g.l[:, Int(g.ny/2+2):end], 2)) == 0.0

testX(g) = isapprox(sum(g.X[end, :].-g.X[1, :]), (g.Lx-g.dx)*g.ny, rtol = 1e-15*g.ny)
testY(g) = isapprox(sum(g.Y[:, end].-g.Y[:, 1]), (g.Ly-g.dy)*g.nx, rtol = 1e-15*g.nx)
testK(g) = isapprox(sum(g.K[2:g.nkr-1, :] .+ flipdim(g.K[g.nkr+1:end, :], 1)), 0, rtol=1e-14)
testL(g) = isapprox(sum(g.L[:, 2:Int(g.ny/2)] .+ flipdim(g.L[:, Int(g.ny/2+2):end], 2)), 0, rtol=1e-14)
testKr(g) = isapprox(sum(g.K[1:g.nkr, :] .- g.Kr), 0, rtol=1e-14)
testLr(g) = isapprox(sum(g.L[1:g.nkr, :] .- g.Lr), 0, rtol=1e-14)

g2 = CuTwoDGrid(32, 2π, 24, 4π)

@test testdx(g2)
@test testdy(g2)
@test testx(g2)
@test testk(g2)
@test testkr(g2)

@test testy(g2)
@test testl(g2)

@test testX(g2)
@test testY(g2)
@test testK(g2)
@test testL(g2)
@test testKr(g2)
@test testLr(g2)


text:
# See https://github.com/ecmwf/cfgrib.jl for more information about CfGRIB
# You can install the HEAD version directly from github by hitting ] to enter package mode and then
#.   add https://github.com/ecmwf/cfgrib.jl/
# followed by backspace
using CfGRIB
using AxisArrays
using Printf
using Dates
using HTTP
using GZip

"""
readMRMS(fname)

Extracts MRMS data in grib format from a specified file and then extracts specific
data elements into CSV formatted files in a directory call csv (one file per location
with non-zero precipitation)

Saves data with detectable precipiation to folder ./csv,
based on truncated lat-lons.
Example: ./csv/48_70.csv
"""
function readMRMS(fname)
    fh = GZip.open(fname, "r")
    data = read(fh)
    io = open("data.grib2", "w")
    write(io, data)
    close(io)

    dataset = CfGRIB.DataSet("data.grib2")

    # Get datetime from dataset, but we could also parse the fname instead
    datetime = replace(string(unix2datetime(dataset.variables["time"].data)), "T" => " ")

    lats = dataset.variables["latitude"].data
    lons = dataset.variables["longitude"].data
    data = dataset.variables["paramId0"].data[:, :]

    files = Dict()

    # MRMS is 3500 x 7000
    for i = 1:3500, j = 1:7000
        lon = lons[j] - 360

        # Conterminous US data only
        if 25 <= lats[i] <= 50 && -125 <= lon <= -67
            # Exclude non-detectable precipitation to save space
            if data[j, i] > 0
                fn = "$(trunc(Int, lats[i]))_$(trunc(Int, -lon))"
                if !haskey(files, fn)
                    files[fn] = open("csv/$fn.csv", "w");
                    write(files[fn], "precipitation,date,lat,lon\n")
                end
                write(files[fn], "$(data[j, i]),$(datetime),$(@sprintf("%.4f", lats[i])),$(@sprintf("%.4f", lon))\r\n")
            end
        end
    end

    for file in values(files)
        close(file)
    end

    # not yet implemented
    # writeToPostgres()
end

"""
processMRMS()

Downloads data from the Iowa State archive of MRMS data.
Calls readMRMS() to extract and process the data.
Reruns after a 15 minute delay.
"""
function processMRMS()
    for date in Dates.Date(2021, 1, 1):Day(1):today()
        year = Dates.year(date)
        m2 = lpad(Dates.month(date), 2, "0")
        d2 = lpad(Dates.day(date),   2, "0")
        for hour in 0:23
            h2 = lpad(hour, 2, "0")
            fname = "MultiSensor_QPE_01H_Pass2_00.00_$year$m2$d2-$(h2)0000.grib2.gz"
            if !isfile("mrms/$fname")
                url = "https://mtarchive.geol.iastate.edu/$year/$m2/$d2/mrms/ncep/MultiSensor_QPE_01H_Pass2/$fname"
                try
                    download(url, "mrms/$fname")
                catch
                    rm("mrms/$fname") # file gets saved even if 404 error, so it must be removed
                    println("Couldn't download $url")
                    continue
                end

                println("Downloaded $fname")
                readMRMS("mrms/$fname")
            end
        end
    end
end

while true
    processMRMS()
    println("Pausing for 15 minutes")
    sleep(60 * 15)
end

max_stars_repo_name: gunnarnewell/weather-server

max_issues_repo_name: gunnarnewell/weather-server

max_forks_repo_name: gunnarnewell/weather-server


text:
###########################
### DYNAMIC PROGRAMMING ###
###########################

### STEP 0
α = 0.3
δ = 0.1
β = 0.9

kupper = 2
klower = 0.001
n = 10
kgrid = collect(range(klower, stop=kupper, length=n))

### STEP 1
tolerance = 0.001
imax = 1000
vnew = zeros(length(kgrid)) # initial value function guess
v = vnew .+ 2*tolerance # initialize v in a way that ensures loop will start
global cartesianindex = Array{CartesianIndex{2}, n} #vector of cartesion coordinates
i = 1

while maximum(abs.(v - vnew)) > tolerance && i <= imax
   ### STEP 2a - calculate frist by columns in the first row and then jump to the other rwo
    # pre-allocate memory
    v = vnew
    c = zeros(n, n)
    for i in 1:n
        for j in 1:n
            c[i,j] = kgrid[i]^α + (1-δ)*kgrid[i] - kgrid[j]
            if c[i,j] < 0
                c[i,j] = 0
            end
        end
    ### STEP 2b
    (vnew, cartesianindex) = findmax(log.(c) .+ β*v', dims = 2)
    end
    i += 1 #adding counter
end

vnew

scatter(kgrid, vnew, title = "v(k)")

# Policy function
kprimeindex = getindex.(cartesianindex, 2) #take the second index of cartesiona index
kprime = kgrid[kprimeindex]
scatter(kgrid, kprime, title = "k'(k)")

max_stars_repo_name: sebaschoo/Dynamic-Programming

max_issues_repo_name: sebaschoo/Dynamic-Programming

max_forks_repo_name: sebaschoo/Dynamic-Programming


text:
## Allows to iterate over kernels
Base.length(::Kernel) = 1
Base.iterate(k::Kernel) = (k, nothing)
Base.iterate(k::Kernel, ::Any) = nothing

printshifted(io::IO, o, shift::Int) = print(io, o)

# Fallback implementation of evaluate for `SimpleKernel`s.
(k::SimpleKernel)(x, y) = kappa(k, evaluate(metric(k), x, y))


text:
using Arrow, DataFrames
using PyPlot

include("plot_styles.jl")

basedir = "/work/sschaub/JuliaForHEP/final_plotting/"
output_dir = joinpath(basedir, "lola+all_1")
recorded_measures = DataFrame(Arrow.Table(joinpath(output_dir, "recorded_measures.arrow")))

fig, ax = subplots(figsize=(12, 9))
ax.plot([recorded_measures.train_loss recorded_measures.test_loss recorded_measures.validation_loss])
ax.set_title("Loss Values")
ax.set_xlabel("Epoch")
ax.set_ylabel("Loss")
ax.legend(["Training Loss", "Testing Loss", "Validation Loss"])
annotate_cms(ax)

fig.savefig(joinpath(output_dir, "losses.pdf"))


text:
sessions = Dict()

function start_http_and_websocket_servers(ip_addr, http_port, http_router, websocket_port, websocket_router)
    websocket_task = @task HTTP.WebSockets.listen(ip_addr, UInt16(websocket_port)) do ws
        while !eof(ws)
            local data
            process = true
            try
                data = readavailable(ws)
            catch e
                process = false
            end
            process && report_ws_errors(ws, revise_function(websocket_router))(ws, data)
        end
    end

    try
        schedule(websocket_task)
        HTTP.serve(report_http_errors(revise_function(http_router)), ip_addr, http_port; reuse_limit=0,)
    finally
        schedule(websocket_task, "exit"; error=true)
    end
end

# Allow live editing of routes and request-handling code
# NOTE: Passed function must be named (not anonymous) for live revisions to work
function revise_function(func)
    (args...; kw...) -> begin
        try
            Revise.revise()
        catch x
            printstyled(x, '\n'; color=:red)
            for level in stacktrace(catch_backtrace())
                println(level)
            end
            return HTTP.Response(500, "There was a revision error")
        end

        Base.invokelatest(func, args...; kw...)
    end
end

# Used by route_pages function
const SCHEMES = Dict{String, Val}("http" => Val{:http}(), "https" => Val{:https}())
const EMPTYVAL = Val{()}()

#Match URL's to handler functions
#Attempts to get to private pages are rerouted to the login page if not currently logged in
function route_http_request(request::HTTP.Request; public=Dict(), private=Dict(), login=missing, default=missing, public_path=Dict())
    sessionid = get_sessionid(request)
    m = Val(Symbol(request.method))
    uri = HTTP.URI(request.target)
    s = get(SCHEMES, uri.scheme, EMPTYVAL)
    h = Val(Symbol(uri.host))
    p = uri.path
    segments = split(uri.path, '/'; keepempty=false)
    # Check if URL is private (requires a login)
    if haskey(private, uri.path)
        if isloggedin(sessionid)
            return private[uri.path](request, sessionid)
        else
            printstyled("Falling back to login\n", color=:red)
            if !ismissing(login)
                public[login](request, sessionid)
            else
                println("404 $(uri.path)")
                HTTP.Response(404, "No such page")
            end
        end
    # Handle public URL
    elseif haskey(public, uri.path)
        public[uri.path](request, sessionid)
    # Handle public URLs with certain prefixes
    else
        for (prefix, func) in pairs(public_path)
            if startswith(uri.path, prefix)
                return func(request, uri.path[length(prefix):end], sessionid)
            end
        end
        if !ismissing(default)
            default(request, sessionid)
        else
            println("404 $(uri.path)")
            HTTP.Response(404, "No such page")
        end
    end
end

function route_websocket_request(ws, data; public=Dict(), private=Dict(), default=missing)
    local msg
    is_json = true
    try
        msg = JSON3.read(data)
    catch e
        is_json = false
    end

    if is_json
        _type = msg[:type]
        if haskey(private, _type) private[_type](ws, msg)
        elseif haskey(public, _type) public[_type](ws, msg)
        elseif !ismissing(default) default(ws, msg)
        end
    end
end

# Output a response including the sessionid cookie
function response(sessionid, code, headers; body="")
    HTTP.Response(code, vcat(["Set-Cookie" => String(HTTP.Cookie("sessionid", string(sessionid)))], headers); body=body)
end

# Determine the current valid sessionid or empty string
function get_sessionid(request::HTTP.Request)
    for cookie in HTTP.Cookies.readcookies(request.headers, "sessionid")
        haskey(sessions, cookie.value) && return cookie.value
    end
    return ""
end

# Determine if user is logged in
isloggedin(sessionid) = haskey(sessions, sessionid)

# Show generic error in browser but details in console
function report_ws_errors(ws, func)
    (args...; kw...) -> begin
        try
            func(args...; kw...)
        catch x
            printstyled(x, '\n'; color=:red)
            for level in stacktrace(catch_backtrace())
                println(level)
            end
            write(ws, Dict(:type=>"error", :msg=>"There was an error"))
        end
    end
end

# Show generic error in browser but details in console
function report_http_errors(func)
    (args...; kw...) -> begin
        try
            func(args...; kw...)
        catch x
            printstyled(x, '\n'; color=:red)
            for level in stacktrace(catch_backtrace())
                println(level)
            end
            HTTP.Response(500, "There was an error")
        end
    end
end


text:
__precompile__()

"""
# Bootstrap package

Statistical bootstrapping for Julia
"""
module Bootstrap

using Statistics
using Random
using StatsBase
using Distributions
using DataFrames
using StatsModels

import Base: iterate, eltype, length, size
import StatsBase: confint, stderror

export
    Model,
    nvar,
    nrun,
    draw!,
    rademacher,
    mammen,
    iquantile,
    bootstrap,
    BasicSampling,
    AntitheticSampling,
    BalancedSampling,
    ExactSampling,
    MaximumEntropySampling,
    ResidualSampling,
    WildSampling,
    BootstrapSample,
    NonParametricBootstrapSample,
    ParametricBootstrapSample,
    original,
    straps,
    data,
    model,
    statistic,
    sampling,
    noise,
    bias,
    stderror,
    confint,
    BasicConfInt,
    PercentileConfInt,
    NormalConfInt,
    BCaConfInt,
    StudentConfInt,
    level

include("stats.jl")
include("nobs.jl")
include("draw.jl")
include("bootsampling.jl")
include("get.jl")
include("show.jl")
include("confint.jl")
include("datasets/Datasets.jl")
include("deprecates.jl")

end


text:
"""
    with_lengthscale(kernel::Kernel, lengthscale::Real)

Construct a transformed kernel with `lengthscale`.

# Examples

```jldoctest
julia> kernel = with_lengthscale(SqExponentialKernel(), 2.5);

julia> x = rand(2);

julia> y = rand(2);

julia> kernel(x, y) ≈ (SqExponentialKernel() ∘ ScaleTransform(0.4))(x, y)
true
```
"""
function with_lengthscale(kernel::Kernel, lengthscale::Real)
    return kernel ∘ ScaleTransform(inv(lengthscale))
end

"""
    with_lengthscale(kernel::Kernel, lengthscales::AbstractVector{<:Real})

Construct a transformed "ARD" kernel with different `lengthscales` for each dimension.

# Examples

```jldoctest
julia> kernel = with_lengthscale(SqExponentialKernel(), [0.5, 2.5]);

julia> x = rand(2);

julia> y = rand(2);

julia> kernel(x, y) ≈ (SqExponentialKernel() ∘ ARDTransform([2, 0.4]))(x, y)
true
```
"""
function with_lengthscale(kernel::Kernel, lengthscales::AbstractVector{<:Real})
    return kernel ∘ ARDTransform(map(inv, lengthscales))
end


text:
"""
Non-macro adaptation from `BenchmarkTools.@benchmarkable`.
"""
function benchmarkable(ex, params = [])
  core, setup, teardown, _ = benchmarkable_parts([:(Ref($ex)[])])

  # extract any variable bindings shared between the core and setup expressions
  setup_vars = isa(setup, Expr) ? collectvars(setup) : []
  core_vars = isa(core, Expr) ? collectvars(core) : []
  out_vars = filter(in(setup_vars), core_vars)

  # generate the benchmark definition
  bench = generate_benchmark_definition(
    @__MODULE__,
    out_vars,
    setup_vars,
    core,
    setup,
    teardown,
    Parameters(; params...),
  )
end

function add_benchmark!(group::BenchmarkGroup, ex::Expr)
  group[string(ex)] = benchmarkable(ex)
end

function add_benchmark!(group::BenchmarkGroup, exs)
  add_benchmark!.(Ref(group), exs)
  group
end

function display_trial(trial::Trial)
  trialallocs = allocs(trial)
  trialmin = minimum(trial)

  println(
    "\e[33;1;1m",
    prettytime(time(trial)),
    "\e[m (", trialallocs , " allocation",
    trialallocs == 1 ? "" : "s", ": ",
    prettymemory(memory(trial)), ")",
  )
end

function display_group(group::BenchmarkGroup)
    nchars_left = maximum(length, keys(group.data))
    foreach(group.data) do (str, val)
      if val isa Trial
        print(' '^2, rpad(str, nchars_left), " :  ")
        display_trial(val)
      else
        println("\e[34;1;1m", str, "\e[m")
        display_group(val)
      end
    end
end


text:
"""
Mask out data above a threshold `t`
"""
function ThresholdMask(q::Quantity,t)
    Mask(times(q),Array{Bool}(quantity(q).<t))
end

@quantity MaskInt Int

MaskInt(m::Mask) = MaskInt(times(m),Array{Int}(quantity(m)))

"""
Tidally average a mask and threshold

This provides a mask for excluding tides
from computations.

- `dm`: the mask that you want to average
- `m` : the mask that segments the tides
- `t` : the percentage of good data points required
"""
function average_mask(dm::Mask,m::Mask,t)
    Tm = tidalaverage(MaskInt(dm),m)
    Mask(times(Tm),Array{Bool}(quantity(Tm).>t))
end


text:
using Test
using Elemental

for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    A = Elemental.DistMatrix(T)
    Elemental.gaussian!(A, 10)
    @test convert(Elemental.DistMatrix{T}, ones(T, 5)) == ones(5,1)
    B = Elemental.DistMatrix(T)
    Array(copy!(B, T[2 1; 1 2])) == T[2 1; 1 2]
end


text:
function decompose_RX_on_5qubits()

    println(">>>>> RX Gate on fourth qubit using U3Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 5,
    "maximum_depth" => 3,
    "elementary_gates" => ["H_1xCNot_2_3xI_4xI_5", "RX_2", "CU3_3_4", "Identity"],
    "target_gate" => QCOpt.kron_two_qubit_gate(5, QCOpt.CRXGate(π/4), "q3", "q4"),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "CU3_θ_discretization" => [0, π/4],
    "CU3_ϕ_discretization" => [0, -π/2],
    "CU3_λ_discretization" => [0, π/2],
    "RX_discretization" => [π/2],
    )

end


text:
"""
    ImbensWagerOptRD(B::Float64, solver, num_buckets=2000, variance=Homoskedastic()) <: SharpRD

A `SharpRD` estimator.

`B` is a bound on the second derivative of the mean of
the outcome variable, conditional on the running variable above and below the cutoff.

The `solver` object is a `MathOptInterface` compatible optimizer such as `Mosek.Optimizer` that will be
used to solve the numerical optimization problem described in Imbens & Wager (2019).

The number of bins for the discretization of the running variable that is required to solve
the numerical optimization problem is `num_buckets`.
"""
Base.@kwdef struct ImbensWagerOptRD <: SharpRD
    B::Float64
    solver::Any
    num_buckets = 2000
    variance = Homoskedastic()
end

Base.@kwdef struct FittedOptRD
    tau_est::AbstractFloat
    se_est::AbstractFloat
    maxbias::AbstractFloat
    ci::Array{Float64}
    weights::Array{Float64}
    B::Float64
    coeftable::Any
end

function estimate_var(data::RDData)
    fitted_lm = fit(LinearModel, @formula(Ys ~ Ws * ZsC), data)
    yhat = predict(fitted_lm)
    σ² = mean((data.Ys - yhat) .^ 2) * length(data.Zs) / (length(data.Zs) - 4)
    return σ²
end

function bias_adjusted_gaussian_ci(se; maxbias = 0.0, level = 0.95)
    rel_bias = maxbias / se
    zz = fzero(
        z -> cdf(Normal(), rel_bias - z) + cdf(Normal(), -rel_bias - z) + level - 1,
        0,
        rel_bias - quantile(Normal(), (1 - level) / 2.1),
    )
    zz * se
end

function Statistics.std(::Homoskedastic, γ::Vector, sig2::Float64)
    sqrt(sum(γ .^ 2) * sig2)
end

"""
    fit(method::ImbensWagerOptRD, data::RDData; level=0.95)

Find the min-max optimal RD estimator over a convex class of conditional mean
functions with a bounded second derivative, based on the numerical optimization
procedure described in Imbens & Wager (2019).

Returns a FittedOptRD object with fields `tau_est`, `se_est`, `maxbias`, `ci`,
`weights` the RD estimator weights for each observation, `B` the bound on
the second derivative of the conditional mean functions used in estimation, and
`coeftable` providing an object of type CoefTable for the estimator.
"""
function fit(method::ImbensWagerOptRD, data::RDData; level = 0.95)
    B = method.B
    ZsD = DiscretizedRunningVariable(data.ZsR, method.num_buckets)
    X = ZsD.Zs
    h = ZsD.h
    d = length(X)
    n = ZsD.weights
    ixc = argmin(abs.(ZsD.Zs .- ZsD.cutoff))
    c = ZsD.cutoff
    W = ZsD.Ws
    sig2 = estimate_var(data)
    σ² = zeros(d) .+ sig2
    model = Model(method.solver)

    @variable(model, f[1:d])
    @variable(model, l1 >= 0)
    @variable(model, l2)
    @variable(model, l3)
    @variable(model, l4)
    @variable(model, l5)
    @variable(model, s)

    G = @expression(
        model,
        2 * B .* f .+ l2 .* (1 .- W) +
        l3 .* W +
        l4 .* (X .- c) +
        l5 .* (W .- 0.5) .* (X .- c)
    )

    @constraint(model, f[ixc] == 0)
    @constraint(model, (f[ixc+1] - f[ixc]) / h[ixc] == 0)

    for i = 3:d
        @constraint(model, (f[i] - 2 * f[i-1] + f[i-2]) / (h[i-1] * h[i-2]) - l1 <= 0)
        @constraint(model, (f[i] - 2 * f[i-1] + f[i-2]) / (h[i-1] * h[i-2]) + l1 >= 0)
    end

    qobj = @expression(model, [1 / 2 .* sqrt.(n) .* G ./ sqrt.(σ²); l1])
    @constraint(model, [0.5; s; qobj] in RotatedSecondOrderCone())

    @objective(model, Min, s - l2 + l3)
    optimize!(model)
    γ_xx = -value.(G) ./ (2 .* σ²)
    γ = γ_xx[ZsD.binmap]

    τ = sum(γ .* data.Ys)
    maxbias = value(l1)
    se = Statistics.std(method.variance, γ, sig2)
    ci = bias_adjusted_gaussian_ci(se; maxbias = maxbias, level = level)
    ci = [τ - ci, τ + ci]
    levstr = isinteger(level * 100) ? string(Integer(level * 100)) : string(level * 100)
    res = [τ se maxbias first(ci) last(ci)]
    colnms = ["τ̂"; "se"; "max bias"; "Lower $levstr%"; "Upper $levstr%"]
    rownms = ["Sharp RD estimand"]
    coeftbl = CoefTable(res, colnms, rownms, 0, 0)

    FittedOptRD(
        tau_est = sum(γ .* data.Ys),
        se_est = se,
        maxbias = maxbias,
        ci = ci,
        weights = γ,
        B = B,
        coeftable = coeftbl,
    )
end

function Base.show(io::IO, rdd_fit::RegressionDiscontinuity.FittedOptRD)
    println(io, "Imbens-Wager (2019) optimized regression discontinuity design")
    println(io, "Max Second Derivative Bound: ", string(rdd_fit.B))
    Base.show(io, rdd_fit.coeftable)
end


text:
# A test for solvers
# Maintained by @zsunberg

mutable struct TestSimulator
    rng::AbstractRNG
    max_steps::Int
end

function simulate(sim::TestSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_distribution::Any)

    s = rand(sim.rng, initial_distribution)
    b = initialize_belief(updater, initial_distribution)

    disc = 1.0
    r_total = 0.0

    step = 1

    while !isterminal(pomdp, s) && step <= sim.max_steps # TODO also check for terminal observation
        a = action(policy, b)

        (sp, o, r) = gen(DDNOut(:sp, :o, :r), pomdp, s, a, sim.rng)

        r_total += disc*r

        b = update(updater, b, a, o)

        disc *= discount(pomdp)
        s = sp
        step += 1
    end

    return r_total
end

function simulate(sim::TestSimulator, mdp::MDP, policy::Policy, s)

    disc = 1.0
    r_total = 0.0

    step = 1

    while !isterminal(mdp, s) && step <= sim.max_steps # TODO also check for terminal observation
        a = action(policy, s)

        (sp, r) = gen(DDNOut(:sp, :r), mdp, s, a, sim.rng)

        r_total += disc*r

        disc *= discount(mdp)
        s = sp
        step += 1
    end

    return r_total
end

"""
    test_solver(solver::Solver, problem::POMDP)
    test_solver(solver::Solver, problem::MDP)

Use the solver to solve the specified problem, then run a simulation.

This is designed to illustrate how solvers are expected to function. All solvers should be able to complete this standard test with the simple models in the POMDPModels package.

To run this with a solver called YourSolver, run
```
using POMDPToolbox
using POMDPModels

solver = YourSolver(# initialize with parameters #)
test_solver(solver, BabyPOMDP())

```
"""
function test_solver(solver::Solver, problem::POMDP; max_steps=10, updater=nothing)

    policy = solve(solver, problem)
    if updater==nothing
        updater = POMDPs.updater(policy)
    end

    sim = TestSimulator(MersenneTwister(1), max_steps)

    simulate(sim, problem, policy, updater, initialstate_distribution(problem))
end

function test_solver(solver::Solver, problem::MDP; max_steps=10)

    policy = solve(solver, problem)

    sim = TestSimulator(MersenneTwister(1), max_steps)

    simulate(sim, problem, policy, initialstate(problem, MersenneTwister(0)))
end


text:
function _helperEkmanTransport(
    τ_x :: Float64,
    τ_y :: Float64,
    ϵ   :: Float64,
    f   :: Float64,
)
    s2 = ϵ^2 + f^2
    return (ϵ * τ_x + f * τ_y) / s, (ϵ * τ_y - f * τ_x) / s
end

function calEkmanFlow!(
    τ_x :: AbstractArray{Float64, 2},
    τ_y :: AbstractArray{Float64, 2},
    v_x :: AbstractArray{Float64, 2},
    v_y :: AbstractArray{Float64, 2},
    f   :: AbstractArray{Float64, 2},
    ϵ   :: AbstractArray{Float64, 2},
    h   :: AbstractFloat,
)
    for i = 1:occ.Nx, j = 1:occ.Ny
        M_x[i, j], M_y[i, j] = _helperEkmanTransport(
            τ_x[i, j],
            τ_y[i, j],
            ϵ[i, j],
            f[i, j],
        )
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using MLDataUtils
using StatsBase

# Implement functions needed for DataSubset to work
MLDataUtils.getobs(s::ImageSource, idx) = s[idx]
StatsBase.nobs(s::ImageSource) = length(s)


text:
using CUDA

# `DI/DJ ~ [-1, 0, 1]`, number of threads should be `(nx÷3) * (ny÷3)`.
@i @inline function i_one_step_kernel1!(Δt, hx, hy, u!, w, wold, φ!, φ0, ψ!, ψ0, σ, τ, c::AbstractMatrix{T}, vi::Val{DI}, vj::Val{DJ}) where {T,DI,DJ}
    # update u!
    @routine @invcheckoff begin
        @zeros Float64 Δt2 Δt_hx Δt_hy Δt_hx2 Δt_hy2 Δt_2
        Δt2 += Δt ^ 2
        Δt_hx += Δt / hx
        Δt_hy += Δt / hy
        Δt_hx2 += Δt_hx ^ 2
        Δt_hy2 += Δt_hy ^ 2
        Δt_2 += Δt/2

        i ← ((blockIdx().x-1) * blockDim().x + threadIdx().x)*3 + DI - 1
        j ← ((blockIdx().y-1) * blockDim().y + threadIdx().y)*3 + DJ - 1
    end

    @invcheckoff @inbounds if i<size(u!, 1) && j<size(u!,2)
        @routine begin
            @zeros T σpτ σpτΔt_2 cΔt_hx2 cΔt_hy2 dwx dwy dφx dψy στ anc1 anc2 anc3 uij σij τij φ0ipj φ0imj ψ0ijp ψ0ijm wipj wimj wijp wijm woldij wij cij
            σij += σ[i,j]
            τij += τ[i,j]
            wipj += w[i+1,j]
            wimj += w[i-1,j]
            wijp += w[i,j+1]
            wijm += w[i,j-1]
            wij += w[i,j]
            woldij += wold[i,j]
            φ0ipj += φ0[i+1,j]
            φ0imj += φ0[i-1,j]
            ψ0ijp += ψ0[i,j+1]
            ψ0ijm += ψ0[i,j-1]
            cij += c[i,j]
            σpτ += σij + τij
            σpτΔt_2 += σpτ * Δt_2
            cΔt_hx2 += Δt_hx2 * cij
            cΔt_hy2 += Δt_hy2 * cij
            dwx += wipj + wimj
            dwy += wijp + wijm
            dφx += φ0ipj - φ0imj
            dψy += ψ0ijp - ψ0ijm
            στ += σij * τij

            anc1 += 2
            anc1 -= στ * Δt2
            anc1 -= 2 * cΔt_hx2
            anc1 -= 2 * cΔt_hy2
            anc2 += Δt_2 * Δt_hx
            anc3 += Δt_2 * Δt_hy
            uij += anc1 * wij
            uij += cΔt_hx2  *  dwx
            uij += cΔt_hy2  *  dwy
            uij += anc2 * dφx
            uij += anc3 * dψy
            uij -= woldij
            uij += σpτΔt_2 * woldij
            σpτΔt_2 += 1
        end
        u![i,j] += uij / σpτΔt_2
        ~@routine
    end
    ~@routine
end

@i @inline function i_one_step_kernel2!(Δt, hx, hy, u!, w, wold, φ!, φ0, ψ!, ψ0, σ, τ, c::AbstractMatrix{T}, vi::Val{DI}, vj::Val{DJ}) where {T,DI,DJ}
    @routine @invcheckoff begin
        i ← ((blockIdx().x-1) * blockDim().x + threadIdx().x)*3 + DI - 1
        j ← ((blockIdx().y-1) * blockDim().y + threadIdx().y)*3 + DJ - 1
        @zeros Float64 Δt_hx Δt_hy
        Δt_hx += Δt / hx
        Δt_hy += Δt / hy
    end
    @invcheckoff @inbounds if i<size(u!, 1) && j<size(u!,2)
        @routine begin
            @zeros T σmτ σmτ_2 dux duy cσmτ_2 σΔt τΔt anc1 anc2 σij τij uipj uimj uijp uijm φ0ij ψ0ij cij
            σij += σ[i,j]
            τij += τ[i,j]
            uipj += u![i+1,j]
            uimj += u![i-1,j]
            uijp += u![i,j+1]
            uijm += u![i,j-1]
            φ0ij += φ0[i,j]
            ψ0ij += ψ0[i,j]
            cij += c[i,j]
            σmτ += σij - τij
            σmτ_2 += σmτ / 2
            dux += uipj - uimj
            duy += uijp - uijm
            cσmτ_2 += cij * σmτ_2
            σΔt += Δt * σij
            τΔt += Δt * τij
            σΔt -= 1
            τΔt -= 1
            anc1 += Δt_hx * cσmτ_2
            anc2 += Δt_hy * cσmτ_2
        end
        φ![i,j] -= σΔt * φ0ij
        φ![i,j] -=  anc1 * dux
        ψ![i,j] -= τΔt * ψ0ij
        ψ![i,j] += anc2 * duy
        ~@routine
    end
    ~@routine
end

let
    exprs = [:((threads, blocks) ← cudiv(ceil(Int,param.NX/3), ceil(Int,param.NY/3)))]
    for KF in [:(i_one_step_kernel1!), :(i_one_step_kernel2!)]
        for (DI, DJ) in Base.Iterators.product((0,1,2), (0,1,2))
            push!(exprs, :(@cuda threads=threads blocks=blocks $KF(
                    param.DELTAT, param.DELTAX, param.DELTAY, u, w, wold,
                    φ, φ0, ψ, ψ0, param.Σx, param.Σy, c,
                    Val($DI), Val($DJ))))
        end
    end
    ex = Expr(:block, exprs..., :((threads, blocks) → cudiv(ceil(Int,param.NX/3), ceil(Int,param.NY/3))))
    @eval @i function i_one_step_parallel!(param::AcousticPropagatorParams, u, w, wold, φ, φ0, ψ, ψ0, c::AbstractMatrix{T}) where T
        $ex
    end
end

@i function bennett_step!(dest::T, src::T, param::AcousticPropagatorParams, srcloc, srcv, c) where T<:SeismicState{<:CuArray}
    @routine begin
        d2 ← zero(param.DELTAT)
        d2 += param.DELTAT^2
    end
    dest.upre += src.u
    dest.step[] += src.step[] + 1
    @safe CUDA.synchronize()
    i_one_step_parallel!(param, dest.u, src.u, src.upre,
        dest.φ, src.φ, dest.ψ, src.ψ, c)
    dest.u[SafeIndex(srcloc)] += srcv[dest.step[]] * d2
    ~@routine
end

function treeverse_step(s::CuSeismicState, param, src, srcv, c::CuMatrix)
    unext, u, φ, ψ = zero(s.u), copy(s.u), copy(s.φ), copy(s.ψ)
    one_step!(param, unext, u, s.upre, φ, ψ, param.Σx, param.Σy, c)
    s2 = SeismicState(u, unext, φ, ψ, Ref(s.step[]+1))
    s2.u[SafeIndex(src)] += srcv[s2.step[]]*param.DELTAT^2
    return s2
end

function treeverse_grad(x::CuSeismicState, g::CuSeismicState, param, src, srcv, gsrcv, c::CuMatrix, gc::CuMatrix)
    println("gradient: $(x.step[]+1) -> $(x.step[])")
    #CUDA.memory_status()
    y = treeverse_step(x, param, src, srcv, c)  # this function is not inplace!
    gt = SeismicState([GVar(getfield(y, field), getfield(g, field)) for field in fieldnames(SeismicState)[1:end-1]]..., Ref(y.step[]))
    x = GVar(x)
    srcv = GVar(srcv, gsrcv)
    c = GVar(c, gc)
    (~bennett_step!)(gt, x, param, src, srcv, c)
    gv = grad(srcv)
    gc2 = grad(c)
    gs = grad(x)
    (gs, gv, gc2)
end


text:
using CmdStan, DataFrames, Test, Statistics

tempdir = pwd()

nwarmup, nsamples, nchains = 1000, 1000, 4

stan_code = "
data {
    int<lower = 1> J; // number of cases
    real<lower = 0> onset_day_reported[J]; //day of illness onset
}

parameters {
    real<lower = 0> incubation[J];

    real mu;
    real<lower = 0> sigma;
}

model {
    mu ~ normal(0, 5);
    sigma ~ cauchy(0, 5);

    incubation ~ lognormal(mu,sigma);

    for (j in 1:J) {
        target += normal_lpdf(incubation[j] | onset_day_reported[j]+0.5, 0.5);
    }

}

generated quantities {
    real incubation_mean = exp(mu+sigma^2/2.0);
    real incubation_sd = sqrt((exp(sigma^2)-1.0)*exp(2.0*mu+sigma^2));
}
"

df_inc = DataFrame(Dict("n" => [1, 5, 7, 9, 17, 18, 16, 6, 8, 10, 7, 1, 2, 4, 4, 1],
    "day" => [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 25]))
days_inc = vcat(fill.(df_inc[!,:day],df_inc[!,:n])...)

stan_data = Dict("J" => length(days_inc), "onset_day_reported" => days_inc)

# stan_init = Dict("mu" => 2.5, "sigma" => 0.2, "incubation" => days_inc .+ 0.5)
stan_init = Dict("incubation[$i]"=>v+.5 for (i,v) in enumerate(days_inc))
temp_init = Dict("mu" => 2.5, "sigma" => 1.0)
merge!(stan_init, temp_init)

stanmodel = Stanmodel(
    name = "init_array",
    model = stan_code,
    nchains = nchains,
    #init = [stan_init],
    num_warmup = nwarmup,
    num_samples = nsamples);

rc, a3d, cnames = stan(stanmodel, stan_data,
        init = stan_init, summary = true);

if rc == 0
    sdf  = read_summary(stanmodel)
    @test sdf[sdf.parameters .== :mu, :mean][1] ≈ 2.5 rtol=0.1
end


text:
X = Interval(0,1)
ll = LazyBox(Dict(1=>Interval(0,1)))
ll[3]
@test domaineq(ll[3],Interval(0,1))
@test ndims(ll) == 2
ll[5]
@test ndims(ll) == 3
ll[5]
@test ndims(ll) == 3

@test length(convert(Vector{Interval},ll)) == ndims(ll)
@test length(convert(Vector{Interval},ll,[1,2])) == 2

l1 = LazyBox(Float64)
l1[1]
l1[2] = Interval(0.0,10.0)
rand(l1)
l2 = LazyBox(Int64)
l2[1]
l2[2] = Interval(2,5)
rand(l2)

@test 0 <= mid(l1)[2] <= 10
@test 2 <= mid(l2)[2] <= 5

# Splitting
l1split = mid_split(l1)
@test length(l1split) == 4
# partial_split_box(l2,Dict(2=>3)) Need to fix splitting intervals
l1split = partial_split_box(l1,Dict(1=>0.5))
@test length(l1split) == 2


text:
module RestrictProlong

using Base: tail

export restrict, restrict!, prolong

const DimsLike = Union{Dims,AbstractVector{Int}}

### restrict, for reducing the image size by 2-fold

"""
    Ar = restrict(A[, dims])

Perform two-fold reduction in size along the dimensions listed in
`dims`, or all coordinates if `dims` is omitted.  It anti-aliases A,
so is more accurate than a naive summation over
2×2×... blocks. `restrict` is normalized so as to approximately
preserve the mean value of `A`.

Thought of as an operator, `restrict` is equal to the transpose of `prolong`.
"""
restrict(A::AbstractArray, dims::DimsLike=1:ndims(A)) =  _restrict(A, dims...)

restrict(A::AbstractArray, dim::Integer) = _restrict(A, dim)

@inline _restrict(A, d1, d2, d...) = _restrict(_restrict(A, d1), d2, d...)
_restrict(A) = A

"""
    Ap = prolong(A, [sz::Tuple])

Perform two-fold expansion in size. If `sz` is omitted, all
dimensions are expanded. If `sz` is specified, it can be one of the
following:

- `sz[d] = size(A, d)`: no expansion is performed
- `sz[d] = 2*size(A, d) - 1`: expansion to an odd-size, adding (by interpolation) points on the half-grid between all points along dimension `d`
- `sz[d] = 2*size(A, d)`: expansion to an even-size, estimating the points on the 1/4 and 3/4 grid (including one point beyond either edge of the current span)

Any other choice results in an error. The default choice is to perform
odd-expansion. `prolong` is normalized so as to approximately preserve
the sum of `A`.

Thought of as an operator, `prolong` is equal to the transpose of `restrict`.
"""
prolong(A::AbstractArray) = prolong(A, map(n->2*n-1, _size(A)))

prolong(A::AbstractVector, sz::Integer) = prolong(A, sz, 1)

prolong(A::AbstractArray{T, N}, sz::Dims{N}) where {T, N} = _prolong(A, sz, ntuple(identity, Val(N)))

@inline _prolong(A, sz, dims) = _prolong(prolong(A, sz[1], dims[1]), tail(sz), tail(dims))
_prolong(A, ::Tuple{}, ::Tuple{}) = A

### Inner routines

function _restrict(A::AbstractArray{T, N}, dim::Integer) where {T, N}
    indsA = axes(A)
    sz = dim <= N ? length(indsA[dim]) : 1
    if sz <= 2
        return copy(A)
    end
    newinds = ntuple(i -> i==dim ? oftype(indsA[i], Base.OneTo(restrict_size(sz))) : indsA[i], Val(N))
    A1 = first(A)
    out = similar(A, typeof(A1/4+A1/2), newinds)
    restrict!(out, A, dim)
    out
end

# This properly anti-aliases. The only "oddity" is that edges tend
# towards zero under repeated restriction.

function restrict!(out, A, dim)
    checkdims(out, A, dim)
    indsA = axes(A)
    Rpre = CartesianIndices(indsA[1:dim-1])
    Rpost = CartesianIndices(indsA[dim+1:end])
    _restrict!(out, axes(out, dim), A, Rpre, indsA[dim], Rpost)
end

# Normalized so that out has roughly the same mean as A
@noinline function _restrict!(out::AbstractArray{T}, indout, A, Rpre::CartesianIndices, indA, Rpost::CartesianIndices) where T
    l = length(indA)
    fill!(out, zero(T))
    if isodd(l)
        half = convert(eltype(T), 0.5)
        quarter = convert(eltype(T), 0.25)
        rngA = first(indA):2:last(indA)-1
        rngout = first(indout):last(indout)-1
        for Ipost in Rpost
            for (i,j) in zip(rngA, rngout)
                @inbounds for Ipre in Rpre
                    tmp = quarter*A[Ipre, i+1, Ipost]
                    out[Ipre, j, Ipost]   += half*A[Ipre, i, Ipost] + tmp
                    out[Ipre, j+1, Ipost] += tmp
                end
            end
            @inbounds for Ipre in Rpre
                out[Ipre, last(indout), Ipost] += half*A[Ipre, last(indA), Ipost]
            end
        end
    else
        threeeighths = convert(eltype(T), 0.375)
        oneeighth = convert(eltype(T), 0.125)
        half = convert(eltype(T), 0.5)
        rngA = first(indA)+1:2:last(indA)-1
        rngout = first(indout):last(indout)-1
        for Ipost in Rpost
            @inbounds for Ipre in Rpre
                out[Ipre, first(indout), Ipost] = half*A[Ipre, first(indA), Ipost]
            end
            @inbounds for (i,j) in zip(rngA, rngout)
                for Ipre in Rpre
                    tmpA1, tmpA2 = A[Ipre, i, Ipost], A[Ipre, i+1, Ipost]
                    out[Ipre, j, Ipost]   += threeeighths*tmpA1 + oneeighth*tmpA2
                    out[Ipre, j+1, Ipost] += oneeighth*tmpA1 + threeeighths*tmpA2
                end
            end
            @inbounds for Ipre in Rpre
                out[Ipre, last(indout), Ipost] += half*A[Ipre, last(indA), Ipost]
            end
        end
    end
    out
end

"""
    prolong(A, sz::Integer, dim::Integer)

Perform expansion to size `sz` along dimension `dim`. `sz` must be one
of the valid choices for prolongation size.
"""
function prolong(A::AbstractArray{T, N}, sz::Integer, dim::Integer) where {T, N}
    if dim > N
        sz == 1 || throw(DimensionMismatch("cannot prolong $N-dimensional array to size $sz along dimension $dim (size must be 1)"))
        return copy(A)
    end
    indsA = axes(A)
    l = dim <= ndims(A) ? length(indsA[dim]) : 1
    (sz == l) | (sz == 2*l-1) | (sz == 2*l) || throw(DimensionMismatch("along dimension $dim, sz must be one of $l, $(2*l-1), or $(2*l), got $sz"))
    newinds = ntuple(i -> i==dim ? oftype(indsA[i], Base.OneTo(sz)) : indsA[i], Val(N))
    A1 = first(A)
    out = similar(A, typeof(A1/2+A1/2), newinds)
    Rpre = CartesianIndices(indsA[1:dim-1])
    Rpost = CartesianIndices(indsA[dim+1:end])
    _prolong!(out, newinds[dim], A, Rpre, indsA[dim], Rpost)
end

@noinline function _prolong!(out::AbstractArray{T}, indout, A, Rpre::CartesianIndices, indA, Rpost::CartesianIndices) where T
    l = length(indout)
    fill!(out, zero(T))
    if isodd(l)
        half = convert(eltype(T), 0.5)
        quarter = convert(eltype(T), 0.25)
        rngout = first(indout):2:last(indout)-1
        for Ipost in Rpost
            i = first(indA)
            for j in rngout
                @inbounds for Ipre in Rpre
                    tmpA1, tmpA2 = A[Ipre, i, Ipost], A[Ipre, i+1, Ipost]
                    out[Ipre, j, Ipost]   = half*tmpA1
                    out[Ipre, j+1, Ipost] = quarter*tmpA1+quarter*tmpA2
                end
                i += 1
            end
            @inbounds for Ipre in Rpre
                out[Ipre, last(indout), Ipost] = half*A[Ipre, last(indA), Ipost]
            end
        end
    else
        threeeighths = convert(eltype(T), 0.375)
        oneeighth = convert(eltype(T), 0.125)
        half = convert(eltype(T), 0.5)
        rngout = first(indout)+1:2:last(indout)-1
        rngA = first(indA):last(indA)-1
        for Ipost in Rpost
            @inbounds for Ipre in Rpre
                out[Ipre, first(indout), Ipost] = half*A[Ipre, first(indA), Ipost]
            end
            @inbounds for (i,j) in zip(rngA, rngout)
                for Ipre in Rpre
                    tmpA1, tmpA2 = A[Ipre, i, Ipost], A[Ipre, i+1, Ipost]
                    out[Ipre, j, Ipost]   = threeeighths*tmpA1 + oneeighth*tmpA2
                    out[Ipre, j+1, Ipost] = oneeighth*tmpA1 + threeeighths*tmpA2
                end
            end
            @inbounds for Ipre in Rpre
                out[Ipre, last(indout), Ipost] = half*A[Ipre, last(indA), Ipost]
            end
        end
    end
    out
end

restrict_size(len::Integer) = (len+1)>>1

function checkdims(out, A, dim)
    1 <= dim <= ndims(A) || throw(DimensionMismatch("dim must be between 1 and $(ndims(A)), got $dim"))
    ndims(out) == ndims(A) || throw(DimensionMismatch("A and out must have the same dimensions, got $(ndims(A)) and $(ndims(out))"))
    indsA = axes(A)
    indsout = axes(out)
    for i = 1:ndims(A)
        if i != dim
            indsout[i] == indsA[i] || throw(DimensionMismatch("A and out must have the same axces for all non-reduced dimensions, got $(indsA[i]) and $(indsout[i]) for dimension $i"))
        end
    end
    indAr = indsA[dim]
    indout = indsout[dim]
    length(indout) == restrict_size(length(indAr)) || throw(DimensionMismatch("A and out disagree along the restricted dimension, have axes $indAr and $indout but latter must be of length $(length(indout))"))
    nothing
end

_size(A::AbstractArray) = map(length, axes(A))

end # module


text:
using CUDA.CUDNN:
    cudnnSoftmaxForward,
    cudnnSoftmaxForward!,
    cudnnSoftmaxBackward,
    cudnnSoftmaxAlgorithm_t,
        CUDNN_SOFTMAX_FAST,     # 0, /* straightforward implementation */
        CUDNN_SOFTMAX_ACCURATE, # 1, /* subtract max from every point to avoid overflow */
        CUDNN_SOFTMAX_LOG,      # 2
    cudnnSoftmaxMode_t,
        CUDNN_SOFTMAX_MODE_INSTANCE, # 0, /* compute the softmax over all C, H, W for each N */
        CUDNN_SOFTMAX_MODE_CHANNEL   # 1  /* compute the softmax over all C for each H, W, N */

@testset "cudnn/softmax" begin
    ax,ay = randn(Float32,10,10),randn(Float32,10,10)
    cx,cy = CuArray.((ax,ay))

    function softmaxtest(;
        alpha=1,
        beta=0,
        mode=CUDNN_SOFTMAX_MODE_INSTANCE,
        algo=CUDNN_SOFTMAX_FAST
    )
        d = mode === CUDNN_SOFTMAX_MODE_INSTANCE ? 1 : 2
        x = ax .- maximum(ax, dims=d)
        y = x .- log.(sum(exp.(x), dims=d))
        if algo !== CUDNN_SOFTMAX_LOG; y = exp.(y); end
        add1(x)=reshape(x, (size(x)..., 1))
        if mode === CUDNN_SOFTMAX_MODE_CHANNEL
            y,cx1,cy1 = add1.((y,cx,cy))
        else
            cx1,cy1 = cx,cy
        end
        y0 = alpha * y
        y1 = y0 .+ beta * ay
        @test y0 ≈ cudnnSoftmaxForward(cx1; algo, mode, alpha) |> Array
        @test y1 ≈ cudnnSoftmaxForward!(copy(cy1), cx1; algo, mode, alpha, beta) |> Array
    end

    softmaxtest()
    softmaxtest(alpha=2)
    softmaxtest(beta=2)
    softmaxtest(mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    softmaxtest(mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    softmaxtest(algo=CUDNN_SOFTMAX_FAST)
    softmaxtest(algo=CUDNN_SOFTMAX_ACCURATE)
    softmaxtest(algo=CUDNN_SOFTMAX_LOG)
end


text:
module TBBSim

using Random, Distributions, DataFrames, Parameters, CSV

export Market, MarketParams
export runsims, res_agg, printsimresults

include("tbb_types.jl")
include("simulation.jl")

end # module


text:
using StringEncodings
using UrlDownload
using Unitful
using StaticArrays
import Unitful: Length, Temperature, Quantity, Units
import ZipFile

function load_glass_db(directory::String)
    glass_catalog = Dict{String,Dict}()
    for path in readdir(directory)
        if uppercase(splitext(path)[2]) == ".AGF"
            parse_glass_file!(glass_catalog, joinpath(directory, path))
        end
    end
    return glass_catalog
end

function generate_cat_jl(cat, jlpath)
    catalogs = []
    io = open(jlpath, "w")
    idnum = 1
    glassnames = []
    for nameandcatalog in cat
        catalog_name, catalog = nameandcatalog
        push!(catalogs, catalog_name)
        eval_string = ["module $catalog_name"]
        push!(eval_string, "using ..GlassCat: Glass, GlassID, AGF")
        push!(eval_string, "using StaticArrays: SVector")
        for x in catalog
            let N = 0
                glass_name, glass_info = x
                push!(glassnames, "$catalog_name.$glass_name")
                vals = []
                for fn in string.(fieldnames(Glass))
                    if fn == "ID"
                        push!(vals, "GlassID(AGF, $idnum)")
                    elseif fn in ["D₀", "D₁", "D₂", "E₀", "E₁", "λₜₖ"]
                        push!(vals, repr(get(glass_info, fn, 0.0)))
                    elseif fn == "temp"
                        push!(vals, repr(get(glass_info, fn, 20.0)))
                    elseif fn == "transmission"
                        v = get(glass_info, "transmission", nothing)
                        if isnothing(v)
                            push!(vals, repr(nothing))
                        else
                            str = join(["($(join(a, ", ")))" for a in v], ", ")
                            push!(vals, "[$str]")
                        end
                    elseif fn == "transmissionN"
                        continue
                    else
                        push!(vals, repr(get(glass_info, fn, NaN)))
                    end
                end
                raw_name = glass_info["raw_name"] == glass_name ? "" : " ($(glass_info["raw_name"]))"
                doc_string = "\"\"\"    $catalog_name.$glass_name$raw_name\n"
                doc_string *= "```\n$(rpad("ID:", 25))AGF:$idnum\n"
                doc_string *= "$(rpad("RI @ 587nm:", 25))$(get(glass_info, "Nd", 0.0))\n"
                doc_string *= "$(rpad("Abbe Number:", 25))$(get(glass_info, "Vd", 0.0))\n"
                doc_string *= "$(rpad("ΔPgF:", 25))$(get(glass_info, "ΔPgF", 0.0))\n"
                doc_string *= "$(rpad("TCE (÷1e-6):", 25))$(get(glass_info, "TCE", 0.0))\n"
                doc_string *= "$(rpad("Density:", 25))$(get(glass_info, "p", 0.0))g/m³\n"
                doc_string *= "$(rpad("Valid wavelengths:", 25))$(get(glass_info, "λmin", 0.0))μm to $(get(glass_info, "λmax", 0.0))μm\n"
                doc_string *= "$(rpad("Reference Temp:", 25))$(get(glass_info, "temp", 20.0))°C\n"
                doc_string *= "```\n\"\"\""
                push!(eval_string, doc_string)
                push!(eval_string, "const $glass_name = Glass($(join(vals, ", "))) \n export $glass_name")
            end
            idnum += 1
        end
        push!(eval_string, "end #module \n export $catalog_name \n") # module
        eval_string = join(eval_string, "\n")
        write(io, eval_string * "\n")
    end
    write(io, "const AGF_GLASS_NAMES = [$(join(repr.(glassnames), ", "))]\n")
    write(io, "const AGF_GLASSES = [$(join(glassnames, ", "))]\n")
    close(io)
end

"""This function will download glass catalogs from publicly available sources and extract them to glassdirectory.  You can execute this in the Julia repl or run it as a Julia script from the command line. If you do the latter then cd to the GlassCat directory and enter this at the command line: julia "src/DownloadGlasses.jl" "path to the directory where you want your glass files stored". Schott, Sumita, and NHG have publicly available glass catalogs but you will have to manually download them and extract them to the GlassFiles directory."""
function downloadcatalogs(glassdirectory::String)
    catalogs = ("https://www.nikon.com/products/optical-glass/assets/pdf/nikon_zemax_data.zip",
   # "http://hbnhg.com/down/data/nhgagp.zip", #this zip file has invalid UTF-8 characters in one of the files contained within it which causes ZipFile to crash. There are two files in the zip, one of which has wacky characters in its filename (perhaps Chinese) so you will have to manually extract just the one with the non-wacky characters to the GlassFiles directory.
    "https://www.oharacorp.com/xls/OHARA_201130_CATALOG.zip",
    "https://hoyaoptics.com/wp-content/uploads/2019/10/HOYA20170401.zip",
    "https://www.schott.com/d/advanced_optics/6959f9a4-0e4f-4ef2-a302-2347468a82f5/1.31/schott-optical-glass-overview-zemax-format.zip")
    # can't download directly have to click on box https://www.sumita-opt.co.jp/en/download/
    # can't download directly have to click on selection https://www.schott.com/advanced_optics/english/download/index.html

    getzip(url)  = urldownload(url,compress = :zip, parser = identity)

    function writeglassfile(url,filename::String)
        try
            zipcat = getzip(url)
            # filename = replace(zipcat.name, r"""[ ,.:;?!()&-]""" => "_")
            # filename = replace(filename, "_agf" => ".agf")
            # filename = replace(filename, "_AGF" => ".agf")
            catname = joinpath(glassdirectory,filename)

            @info "reading $filename from web"
            temp = ZipFile.read(zipcat,String)

            @info "writing $filename to $catname"
            write(joinpath(glassdirectory,filename),temp)
        catch err
            @info "Couldn't download $url"
        end
    end

    #write the glass files with standard names so the examples in OpticSim.jl will work
    writeglassfile(catalogs[1],"NIKON.agf")
    writeglassfile(catalogs[2],"OHARA.agf")
    writeglassfile(catalogs[3],"HOYA.agf")
    writeglassfile(catalogs[4],"SCHOTT.agf")
end

function string_list_to_float_list(x)
    npts = length(x)
    if (npts == 0) || ((npts == 1) && (strip(x[1]) == "-"))
        return (repeat([-1.0], 10))
    end
    res = []
    for a in x
        if (strip(a) == "-")
            push!(res, -1.0)
        else
            try
                push!(res, parse(Float64, a))
            catch
                push!(res, NaN)
            end
        end
    end
    return (res)
end

function parse_glass_file!(glass_catalog, filename::String)
    if !isfile(filename)
        throw(error("AGF file doesn't exist"))
    end
    catalog_name = splitext(basename(filename))[1]
    # remove invalid characters
    catalog_name = replace(catalog_name, r"""[ ,.:;?!()&-]""" => "_")
    try
        # cant have module names which are just numbers so add a _ to the start
        parse(Int, catalog_name[1])
        catalog_name = "_" * catalog_name
    catch
        ()
    end
    glass_catalog[catalog_name] = Dict{String,Any}()
    # check whether the file is UTF8 or UTF16 encoded
    if !isvalid(readuntil(filename, " "))
        fo = open(filename, enc"UTF-16LE", "r")
    else
        fo = open(filename, "r")
    end
    # read the file
    glass_name = ""
    let transmission_data = nothing
        for line in readlines(fo)
            if strip(line) == "" || length(strip(line)) == 0 || startswith(line, "CC ") || startswith(line, "GC ")
                continue
            end
            if startswith(line, "NM ")
                transmission_data = Vector{SVector{3,Float64}}(undef, 0)
                nm = split(line)
                glass_name = nm[2]
                original_glass_name = glass_name
                # remove invalid characters
                glass_name = replace(glass_name, "*" => "_STAR")
                glass_name = replace(glass_name, r"""[ ,.:;?!()&-]""" => "_")
                try
                    # cant have module names which are just numbers so add a _ to the start
                    parse(Int, glass_name[1])
                    glass_name = "_" * glass_name
                catch
                    ()
                end
                glass_catalog[catalog_name][glass_name] = Dict{String,Any}()
                glass_catalog[catalog_name][glass_name]["raw_name"] = original_glass_name
                glass_catalog[catalog_name][glass_name]["dispform"] = Int(parse(Float64, nm[3]))
                glass_catalog[catalog_name][glass_name]["Nd"] = parse(Float64, nm[5])
                glass_catalog[catalog_name][glass_name]["Vd"] = parse(Float64, nm[6])
                if length(nm) < 7
                    glass_catalog[catalog_name][glass_name]["exclude_sub"] = 0
                else
                    glass_catalog[catalog_name][glass_name]["exclude_sub"] = Int(parse(Float64, nm[7]))
                end
                if length(nm) < 8
                    glass_catalog[catalog_name][glass_name]["status"] = 0
                else
                    glass_catalog[catalog_name][glass_name]["status"] = Int(parse(Float64, nm[8]))
                end
                if length(nm) < 9 || "-" ∈ nm
                    glass_catalog[catalog_name][glass_name]["meltfreq"] = 0
                else
                    glass_catalog[catalog_name][glass_name]["meltfreq"] = Int(parse(Float64, nm[9]))
                end
            elseif startswith(line, "ED ")
                ed = split(line)
                glass_catalog[catalog_name][glass_name]["TCE"] = parse(Float64, ed[2])
                glass_catalog[catalog_name][glass_name]["p"] = parse(Float64, ed[4])
                glass_catalog[catalog_name][glass_name]["ΔPgF"] = parse(Float64, ed[5])
                if (length(ed) < 6)
                    glass_catalog[catalog_name][glass_name]["ignore_thermal_exp"] = 0
                else
                    glass_catalog[catalog_name][glass_name]["ignore_thermal_exp"] = Int(parse(Float64, ed[6]))
                end
            elseif startswith(line, "CD ")
                cd = parse.(Float64, split(line)[2:end])
                glass_catalog[catalog_name][glass_name]["C1"] = get(cd, 1, NaN)
                glass_catalog[catalog_name][glass_name]["C2"] = get(cd, 2, NaN)
                glass_catalog[catalog_name][glass_name]["C3"] = get(cd, 3, NaN)
                glass_catalog[catalog_name][glass_name]["C4"] = get(cd, 4, NaN)
                glass_catalog[catalog_name][glass_name]["C5"] = get(cd, 5, NaN)
                glass_catalog[catalog_name][glass_name]["C6"] = get(cd, 6, NaN)
                glass_catalog[catalog_name][glass_name]["C7"] = get(cd, 7, NaN)
                glass_catalog[catalog_name][glass_name]["C8"] = get(cd, 8, NaN)
                glass_catalog[catalog_name][glass_name]["C9"] = get(cd, 9, NaN)
                glass_catalog[catalog_name][glass_name]["C10"] = get(cd, 10, NaN)
            elseif startswith(line, "TD ")
                td = parse.(Float64, split(line)[2:end])
                glass_catalog[catalog_name][glass_name]["D₀"] = get(td, 1, 0.0)
                glass_catalog[catalog_name][glass_name]["D₁"] = get(td, 2, 0.0)
                glass_catalog[catalog_name][glass_name]["D₂"] = get(td, 3, 0.0)
                glass_catalog[catalog_name][glass_name]["E₀"] = get(td, 4, 0.0)
                glass_catalog[catalog_name][glass_name]["E₁"] = get(td, 5, 0.0)
                glass_catalog[catalog_name][glass_name]["λₜₖ"] = get(td, 6, 0.0)
                glass_catalog[catalog_name][glass_name]["temp"] = get(td, 7, 20.0)
            elseif startswith(line, "OD ")
                od = string_list_to_float_list(split(line)[2:end])
                glass_catalog[catalog_name][glass_name]["relcost"] = get(od, 1, -1)
                glass_catalog[catalog_name][glass_name]["CR"] = get(od, 2, -1)
                glass_catalog[catalog_name][glass_name]["FR"] = get(od, 3, -1)
                glass_catalog[catalog_name][glass_name]["SR"] = get(od, 4, -1)
                glass_catalog[catalog_name][glass_name]["AR"] = get(od, 5, -1)
                glass_catalog[catalog_name][glass_name]["PR"] = get(od, 6, -1)
            elseif startswith(line, "LD ")
                ld = parse.(Float64, split(line)[2:end])
                glass_catalog[catalog_name][glass_name]["λmin"] = ld[1]
                glass_catalog[catalog_name][glass_name]["λmax"] = ld[2]
            elseif startswith(line, "IT ")
                it_row = parse.(Float64, split(line)[2:end])
                if length(it_row) == 3 && it_row[1] != 0.0
                    entry = SVector{3,Float64}(it_row[1], it_row[2], it_row[3])
                    push!(transmission_data, entry)
                end
                glass_catalog[catalog_name][glass_name]["transmission"] = transmission_data
            end
        end
    end
end