text:
Base.lock(h::Hist1D) = lock(h.hlock)
Base.unlock(h::Hist1D) = unlock(h.hlock)

"""
    sample(h::Hist1D, n::Int=1)

Sample a histogram's with weights equal to bin count, `n` times.
The sampled values are the bins' lower edges.
"""
function sample(h::Hist1D; n::Int=1)
    StatsBase.sample(binedges(h)[1:end-1], Weights(bincounts(h)), n)
end

@inline function _edge_binindex(r::AbstractRange, x::Real)
    return floor(Int, (x - first(r)) * inv(step(r))) + 1
    # # 20% faster and assigns -Inf, Inf, NaN to typemin(Int64)
    # return Base.unsafe_trunc(Int, round((x - first(r)) * inv(step(r)), RoundDown)) + 1
end

@inline function _edge_binindex(r::AbstractRange{<:Integer}, x::Integer)
    return (x - first(r))÷step(r) + 1
end

@inline function _edge_binindex(v::AbstractVector, x::Real)
    return searchsortedlast(v, x)
end

"""
    bincounts(h::Hist1D)

Get the bin counts of a histogram.
"""
@inline bincounts(h::Hist1D) = h.hist.weights

"""
    binedges(h::Hist1D)

Get the bin edges of a histogram.
"""
@inline binedges(h::Hist1D) = h.hist.edges[1]

"""
    bincenters(h::Hist1D)

Get the bin centers of a histogram.
"""
function bincenters(h::Hist1D)
    StatsBase.midpoints(binedges(h))
end

"""
    binerrors(f::T, h::Hist1D) where T<:Function = f.(h.sumw2)
    binerrors(h::Hist1D) = binerrors(sqrt, h)

Calculate the bin errors from `sumw2` with a Gaussian default.
"""
binerrors(f::T, h::Hist1D) where T<:Function = f.(h.sumw2)
binerrors(h::Hist1D) = binerrors(sqrt, h)

"""
    nbins(h::Hist1D)

Get the number of bins of a histogram.
"""
function nbins(h::Hist1D)
    length(bincounts(h))
end

"""
    integral(h::Hist1D)

Get the integral a histogram.
"""
function integral(h::Hist1D)
    sum(bincounts(h))
end

"""
    empty!(h::Hist1D)

Resets a histogram's bin counts and `sumw2`.
"""
function Base.empty!(h::Hist1D{T,E}) where {T,E}
    h.hist.weights .= zero(T)
    h.sumw2 .= 0.0
    return h
end

"""
    push!(h::Hist1D, val::Real, wgt::Real=1)
    atomic_push!(h::Hist1D, val::Real, wgt::Real=1)

Adding one value at a time into histogram.
`sumw2` (sum of weights^2) accumulates `wgt^2` with a default weight of 1.
`atomic_push!` is a slower version of `push!` that is thread-safe.

N.B. To append multiple values at once, use broadcasting via
`push!.(h, [-3.0, -2.9, -2.8])` or `push!.(h, [-3.0, -2.9, -2.8], 2.0)`
"""
@inline function atomic_push!(h::Hist1D{T,E}, val::Real, wgt::Real=1) where {T,E}
    lock(h)
    push!(h, val, wgt)
    unlock(h)
    return nothing
end

@inline function Base.push!(h::Hist1D{T,E}, val::Real, wgt::Real=1) where {T,E}
    r = binedges(h)
    L = nbins(h)
    binidx = _edge_binindex(r, val)
    if h.overflow
        binidx = clamp(binidx, 1, L)
        @inbounds h.hist.weights[binidx] += wgt
        @inbounds h.sumw2[binidx] += wgt^2
    else
        if unsigned(binidx - 1) < L
            @inbounds h.hist.weights[binidx] += wgt
            @inbounds h.sumw2[binidx] += wgt^2
        end
    end
    return nothing
end

Base.broadcastable(h::Hist1D) = Ref(h)

"""
    Hist1D(elT::Type{T}=Float64; binedges, overflow) where {T}

Initialize an empty histogram with bin content typed as `T` and bin edges.
To be used with [`push!`](@ref). Default overflow behavior (`false`)
will exclude values that are outside of `binedges`.
"""
function Hist1D(elT::Type{T}=Float64; bins, overflow=_default_overflow) where {T}
    counts = zeros(elT, length(bins) - 1)
    return Hist1D(Histogram(bins, counts); overflow=overflow)
end

"""
    Hist1D(array, edges::AbstractRange; overflow)
    Hist1D(array, edges::AbstractVector; overflow)

Create a `Hist1D` with given bin `edges` and vlaues from
array. Weight for each value is assumed to be 1.
"""
function Hist1D(A::AbstractVector, r::AbstractRange; overflow=_default_overflow)
    h = Hist1D(Int; bins=r, overflow=overflow)
    push!.(h, A)
    return h
end
function Hist1D(A::AbstractVector, edges::AbstractVector; overflow=_default_overflow)
    if _is_uniform_bins(edges)
        r = range(first(edges), last(edges), length=length(edges))
        return Hist1D(A, r; overflow=overflow)
    else
        h = Hist1D(Int; bins=edges, overflow=overflow)
        push!.(h, A)
        return h
    end
end

"""
    Hist1D(array, wgts::AbstractWeights, edges::AbstractRange; overflow)
    Hist1D(array, wgts::AbstractWeights, edges::AbstractVector; overflow)

Create a `Hist1D` with given bin `edges` and vlaues from
array. `wgts` should have the same `size` as `array`.
"""
function Hist1D(A, wgts::AbstractWeights, r::AbstractRange; overflow=_default_overflow)
    @boundscheck @assert size(A) == size(wgts)
    h = Hist1D(eltype(wgts); bins=r, overflow=overflow)
    push!.(h, A, wgts)
    return h
end
function Hist1D(A, wgts::AbstractWeights, edges::AbstractVector, overflow=_default_overflow)
    @inbounds if _is_uniform_bins(edges)
        r = range(first(edges), last(edges), length=length(edges))
        return Hist1D(A, wgts, r; overflow=overflow)
    else
        h = Hist1D(eltype(wgts); bins=edges, overflow=overflow)
        push!.(h, A, wgts)
        return h
    end
end

"""
    Hist1D(A::AbstractVector{T}; nbins::Integer=_sturges(A), overflow) where T
    Hist1D(A::AbstractVector{T}, wgts::AbstractWeights; nbins::Integer=_sturges(A), overflow) where T

Automatically determine number of bins based on `Sturges` algo.
"""
function Hist1D(A::AbstractVector{T}; nbins::Integer=_sturges(A), overflow=_default_overflow) where {T}
    F = float(T)
    lo, hi = minimum(A), maximum(A)
    r = StatsBase.histrange(F(lo), F(hi), nbins)
    return Hist1D(A, r; overflow=overflow)
end

function Hist1D(
    A::AbstractVector{T},
    wgts::AbstractWeights;
    nbins::Integer=_sturges(A),
    overflow=_default_overflow,
) where {T}
    F = float(T)
    lo, hi = minimum(A), maximum(A)
    r = StatsBase.histrange(F(lo), F(hi), nbins)
    return Hist1D(A, wgts, r; overflow=overflow)
end

"""
    Statistics.mean(h::Hist1D)
    Statistics.std(h::Hist1D)
    Statistics.median(h::Hist1D)
    Statistics.quantile(h::Hist1D, p)

Compute statistical quantities based on the bin centers weighted
by the bin counts.
"""
Statistics.mean(h::Hist1D) = Statistics.mean(bincenters(h), Weights(bincounts(h)))
Statistics.std(h::Hist1D) = sqrt(Statistics.var(bincenters(h), Weights(bincounts(h))))
Statistics.median(h::Hist1D) = Statistics.median(bincenters(h), Weights(bincounts(h)))
Statistics.quantile(h::Hist1D, p) = Statistics.quantile(bincenters(h), Weights(bincounts(h)), p)

"""
    function lookup(h::Hist1D, x)

For given x-axis value `x`, find the corresponding bin and return the bin content.
If a value is out of the histogram range, return `missing`.
"""
function lookup(h::Hist1D, x)
    r = binedges(h)
    !(first(r) <= x <= last(r)) && return missing
    return bincounts(h)[_edge_binindex(r, x)]
end

"""
    normalize(h::Hist1D)

Create a normalized histogram via division by `integral(h)`.
"""
function normalize(h::Hist1D)
    return h*(1/integral(h))
end

"""
    cumulative(h::Hist1D; forward=true)

Create a cumulative histogram. If `forward`, start
summing from the left.
"""
function cumulative(h::Hist1D; forward=true)
    # https://root.cern.ch/doc/master/TH1_8cxx_source.html#l02608
    f = forward ? identity : reverse
    h = deepcopy(h)
    h.hist.weights .= f(cumsum(f(h.hist.weights)))
    h.sumw2 .= f(cumsum(f(h.sumw2)))
    return h
end

"""
    rebin(h::Hist1D, n::Int=1)
    rebin(n::Int) = h::Hist1D -> rebin(h, n)

Merges `n` consecutive bins into one.
The returned histogram will have `nbins(h)/n` bins.
"""
function rebin(h::Hist1D, n::Int=1)
    @assert nbins(h) % n == 0
    p = x->Iterators.partition(x, n)
    counts = sum.(p(bincounts(h)))
    sumw2 = sum.(p(h.sumw2))
    edges = first.(p(binedges(h)))
    if _is_uniform_bins(edges)
        edges = range(first(edges), last(edges), length=length(edges))
    end
    return Hist1D(Histogram(edges, counts), sumw2; overflow=h.overflow)
end
rebin(n::Int) = h::Hist1D -> rebin(h, n)

"""
    restrict(h::Hist1D, low=-Inf, high=Inf)
    restrict(low=-Inf, high=Inf) = h::Hist1D -> restrict(h, low, high)

Returns a new histogram with a restricted x-axis.
`restrict(h, 0, 3)` (or `h |> restrict(0, 3)`)
will return a slice of `h` where the bin centers are in `[0, 3]` (inclusive).
"""
function restrict(h::Hist1D, low=-Inf, high=Inf)
    sel = low .<= bincenters(h) .<= high
    @assert sum(sel) > 0 "No bin centers contained in [$(low), $(high)]"
    edgesel = push!(copy(sel), false)

    # include the right edge of the rightmost selected bin
    lastidx = findlast(edgesel)
    if lastidx != nothing
        edgesel[lastidx+1] = 1
    end

    c = bincounts(h)[sel]
    edges = binedges(h)[edgesel]
    sumw2 = h.sumw2[sel]
    if _is_uniform_bins(edges)
        edges = range(first(edges), last(edges), length=length(edges))
    end
    Hist1D(Histogram(edges, c), sumw2; overflow=h.overflow)
end
restrict(low=-Inf, high=Inf) = h::Hist1D->restrict(h, low, high)

function _svg(h::Hist1D)
    paddingx, paddingy = 0.05, 0.10
    framewidth, frameheight = 250, 200
    fill = "#ffffff00" # 0 alpha
    strokecolor, strokewidth = "black", 1
    strokewidth = 1
    _c, _e = copy(bincounts(h)), copy(binedges(h))

    hasneg = minimum(_c) < 0
    zerolineyfrac = 0.0
    if hasneg
        zerolineyfrac = -minimum(_c)/(maximum(_c)-minimum(_c))
        _c .-= minimum(_c)
    end
    ys = frameheight * ((2*paddingy-1)/maximum(_c) .* _c .+ (1-paddingy))
    xs = framewidth * (
        (1 - 2 * paddingx) / (maximum(_e) - minimum(_e))
        .* (_e .- minimum(_e))
        .+ paddingx
    )
    points = [(paddingx*framewidth, (1-paddingy)*frameheight)] # bottom left
    for i in 1:(length(xs)-1)
        push!(points, (xs[i],ys[i])) # left bin edge
        push!(points, (xs[i+1],ys[i])) # right bin edge
    end
    push!(points, ((1-paddingx)*framewidth,(1-paddingy)*frameheight))
    if !hasneg
        push!(points, points[1]) # close path
    end
    pathstr = join(["$(x),$(y)" for (x,y) in points],",")
    zeroliney = hasneg ? frameheight*(1-zerolineyfrac) : frameheight*(1-paddingy)
    return """
    <svg width="$(framewidth)" height="$(frameheight)" version="1.1" xmlns="http://www.w3.org/2000/svg">
        <polyline points="$(pathstr)" stroke="$(strokecolor)" fill="$(fill)" stroke-width="$(strokewidth)"/>
        <polyline points="$(framewidth*paddingx),$(zeroliney),$(framewidth*(1-paddingx)),$(zeroliney)" stroke="black" stroke-width="1"/>
        <text x="$(framewidth*paddingx)" y="$(frameheight*(1-0.5*paddingy))" dominant-baseline="middle" text-anchor="start" fill="black">$(minimum(_e))</text>
        <text x="$(framewidth*(1-paddingx))" y="$(frameheight*(1-0.5*paddingy))" dominant-baseline="middle" text-anchor="end" fill="black">$(maximum(_e))</text>
    </svg>
    """
end

function Base.show(io::IO, h::Hist1D)
    compact = get(io, :compact, false)
    if compact
        print(io, "Hist1D{$(eltype(bincounts(h)))}, ")
        print(io, "edges=$(repr(binedges(h), context=:limit => true)), ")
        print(io, "integral=$(integral(h))")
    else
        if (nbins(h) < 50) && all(bincounts(h) .>= 0)
            _e = binedges(h)
            _h = Histogram(float.(_e), bincounts(h))
            show(io, UnicodePlots.histogram(_h; width=30, xlabel=""))
        end
        println(io)
        println(io, "edges: ", binedges(h))
        println(io, "bin counts: ", bincounts(h))
        print(io, "total count: ", integral(h))
    end
end

function Base.show(io::IO, m::MIME"text/html", h::Hist1D)
    println(io, """

$(_svg(h))

edges: $(repr(binedges(h), context=:limit => true))

bin counts: $(repr(bincounts(h), context=:limit => true))

maximum count: $(maximum(bincounts(h)))

total count: $(integral(h))

    """)
end


text:
@testset "Random primitives" begin

@testset "rand(::Type{<:Integer}) statements" begin

function foo(n::Int)
    a = rand(Int8)
    b = rand(Int16, n)
    c = rand(Int32, (n, n))
end

genfoo = genify(foo, Int)

# Check types
trace = simulate(genfoo, (10,))
@test trace[:a] isa Int8
@test trace[:b] isa Vector{Int16}
@test trace[:c] isa Matrix{Int32}

# Test array sizes
@test size(trace[:a]) == ()
@test size(trace[:b]) == (10,)
@test size(trace[:c]) == (10, 10)

# Test sample range
@test all(typemin(Int8) .<= trace[:a] .<= typemax(Int8))
@test all(typemin(Int16) .<= trace[:b] .<= typemax(Int16))
@test all(typemin(Int32) .<= trace[:c] .<= typemax(Int32))

# Test generate with constraints
trace, weight = generate(genfoo, (10,), choicemap(:a => 0))
@test weight ≈ -8*log(2)
trace, weight = generate(genfoo, (10,), choicemap(:b => zeros(Int16, 10)))
@test weight ≈ -16*log(2) * 10
trace, weight = generate(genfoo, (10,), choicemap(:c => zeros(Int16, 10, 10)))
@test weight ≈ -32*log(2) * 10 * 10

end

@testset "rand(::Type{<:AbstractFloat}) statements" begin

function foo(n::Int)
    a = rand(Float16)
    b = rand(Float32, n)
    c = rand(Float64, n, n)
end

genfoo = genify(foo, Int)

# Check types
trace = simulate(genfoo, (10,))
@test trace[:a] isa Float16
@test trace[:b] isa Vector{Float32}
@test trace[:c] isa Matrix{Float64}

# Test array sizes
@test size(trace[:a]) == ()
@test size(trace[:b]) == (10,)
@test size(trace[:c]) == (10, 10)

# Test ranges
@test all(0 .<= trace[:a] .<= 1)
@test all(0 .<= trace[:b] .<= 1)
@test all(0 .<= trace[:c] .<= 1)

# Test generate with constraints
trace, weight = generate(genfoo, (10,), choicemap(:a => 0.5))
@test weight == 0
trace, weight = generate(genfoo, (10,), choicemap(:b => [-1; fill(0, 9)]))
@test weight == -Inf
trace, weight = generate(genfoo, (10,), choicemap(:c => rand(Float64, 10, 10)))
@test weight == 0

end

@testset "rand(::Indexable|Setlike) statements" begin

function foo(dims...)
    a = rand([:h, :e, :l, :l, :o], dims...)
    b = rand("world", dims...)
    c = rand(1:5, dims...)
    d = rand((9, 6, 42), dims...)
    e = rand(Set("hello"), dims...)
    f = rand(Dict(zip("world", "hello")), dims...)
end

genfoo = genify(foo)

# Check types
trace = simulate(genfoo, ())
@test trace[:a] isa Symbol
@test trace[:b] isa Char
@test trace[:c] isa Int
@test trace[:d] isa Int
@test trace[:e] isa Char
@test trace[:f] isa Pair{Char,Char}

# Test array sizes
trace = simulate(genfoo, (10, 10))
@test size(trace[:a]) == (10, 10)
@test size(trace[:b]) == (10, 10)
@test size(trace[:c]) == (10, 10)
@test size(trace[:d]) == (10, 10)
@test size(trace[:e]) == (10, 10)
@test size(trace[:f]) == (10, 10)

# Test ranges
@test all(x in [:h, :e, :l, :l, :o] for x in trace[:a])
@test all(x in "world" for x in trace[:b])
@test all(x in 1:5 for x in trace[:c])
@test all(x in (9, 6, 42) for x in trace[:d])
@test all(x in Set("hello") for x in trace[:e])
@test all(x in Dict(zip("world", "hello")) for x in trace[:f])

# Test generate with constraints
trace, weight = generate(genfoo, (), choicemap(:a => :l))
@test weight ≈ log(2/5)
trace, weight = generate(genfoo, (), choicemap(:b => 'w'))
@test weight ≈ log(1/5)
trace, weight = generate(genfoo, (), choicemap(:c => 3))
@test weight ≈ log(1/5)
trace, weight = generate(genfoo, (), choicemap(:d => 42))
@test weight ≈ log(1/3)
trace, weight = generate(genfoo, (), choicemap(:e => 'l'))
@test weight ≈ log(1/4)
trace, weight = generate(genfoo, (), choicemap(:f => ('o' => 'e')))
@test weight ≈ log(1/5)

end

@testset "rand(::Distributions.Distribution) statements" begin

function foo(dims...)
    a = rand(Dirichlet(5, 2.0))
    b = rand(Categorical(a))
    c = rand(Wishart(b+2, [1.0 0.0; 0.0 1.0]))
    d = rand(Normal(0, 1), dims...)
end

genfoo = genify(foo)

# Check types
trace = simulate(genfoo, (10,))
@test trace[:a] isa Vector{Float64}
@test trace[:b] isa Int
@test trace[:c] isa Matrix{Float64}
@test trace[:d] isa Vector{Float64}

# Test array sizes
@test size(trace[:a]) == (5,)
@test size(trace[:b]) == ()
@test size(trace[:c]) == (2, 2)
@test size(trace[:d]) == (10,)

# Test ranges
@test all(0 .<= trace[:a] .<= 1)
@test all(1 .<= trace[:b] .<= 5)
@test insupport(Wishart, trace[:c])

# Test generate with constraints
trace, weight = generate(genfoo, (), choicemap(:a => ones(5)/5))
@test weight ≈ Distributions.logpdf(Dirichlet(5, 2), ones(5)/5)
trace, weight = generate(genfoo, (), choicemap(:a => ones(5)/5, :b => 5))
@test weight ≈ Distributions.logpdf(Dirichlet(5, 2), ones(5)/5) + log(1/5)
trace, weight = generate(genfoo, (), choicemap(:c => [1 0; 0 1]))
@test weight ≈ Distributions.logpdf(Wishart(trace[:b]+2, [1 0; 0 1]), [1 0; 0 1])
trace, weight = generate(genfoo, (10, 10), choicemap(:d => zeros(10, 10)))
@test weight ≈ Gen.logpdf(normal, 0, 0, 1) * 100

end

@testset "randn() and randexp() statements" begin

function foo(dims...)
    a = randn()
    b = randexp()
    c = randn(Float64, dims)
    d = randexp(dims...)
end

genfoo = genify(foo)

# Check types
trace = simulate(genfoo, (10,))
@test trace[:a] isa Float64
@test trace[:b] isa Float64
@test trace[:c] isa Vector{Float64}
@test trace[:d] isa Vector{Float64}

# Test array sizes
@test size(trace[:a]) == ()
@test size(trace[:b]) == ()
@test size(trace[:c]) == (10,)
@test size(trace[:d]) == (10,)

# Test ranges
@test all(trace[:b] .>= 0)
@test all(trace[:d] .>= 0)

# Test generate with constraints
trace, weight = generate(genfoo, (), choicemap(:a => 1, :b => 1))
@test weight ≈ Gen.logpdf(normal, 0, 1, 1) + Gen.logpdf(exponential, 1, 1)
trace, weight = generate(genfoo, (10,), choicemap(:c => ones(10), :d => ones(10)))
@test weight ≈ 10 * (Gen.logpdf(normal, 0, 1, 1) + Gen.logpdf(exponential, 1, 1))

end

@testset "randperm(), randcycle(), and shuffle() statements" begin

function foo(n::Int)
    for i in 1:100
        a = randperm(n)
        b = randcycle(n)
    end
    c = shuffle([:a, :b, :c, :d, :e])
    d = shuffle(['h', 'e', 'l', 'l', 'o', 'o'])
end

genfoo = genify(foo, Int)

# Check types
trace = simulate(genfoo, (5,))
@test all(trace[:a => i] isa Vector{Int} for i in 1:100)
@test all(trace[:b => i] isa Vector{Int} for i in 1:100)
@test trace[:c] isa Vector{Symbol}
@test trace[:d] isa Vector{Char}

# Test array sizes
@test all(length(trace[:a => i]) == 5 for i in 1:100)
@test all(length(trace[:b => i]) == 5 for i in 1:100)
@test length(trace[:c]) == 5
@test length(trace[:d]) == 6

# Test ranges
@test all(isperm(trace[:a => i]) for i in 1:100)
@test all(Genify.iscycle(trace[:b => i]) for i in 1:100)
@test Set(trace[:c]) == Set([:a, :b, :c, :d, :e])
@test sort(trace[:d]) == ['e', 'h', 'l', 'l', 'o', 'o']

# Test generate with constraints
choices = choicemap((:a => 1, randperm(5)), (:b => 1, randcycle(5)))
trace, weight = generate(genfoo, (5,), choices)
@test weight ≈ log(1/factorial(5)) + log(1/factorial(4))
trace, weight = generate(genfoo, (5,), choicemap((:c,  [:e, :b, :d, :c, :a])))
@test weight == log(1/factorial(5))
trace, weight = generate(genfoo, (5,), choicemap((:d,  ['h', 'l', 'l', 'o', 'e', 'o'])))
@test weight == log((factorial(2) * factorial(2)) / factorial(6))

trace, weight = generate(genfoo, (5,), choicemap((:a => 1,  ones(Int, 5))))
@test weight == -Inf
trace, weight = generate(genfoo, (5,), choicemap((:b => 1,  [1, 2, 3, 4, 5])))
@test weight == -Inf
trace, weight = generate(genfoo, (5,), choicemap((:c,  [:e, :b, :d, :c, :f])))
@test weight == -Inf
trace, weight = generate(genfoo, (5,), choicemap((:d,  ['h', 'l', 'l', 'e', 'o'])))
@test weight == -Inf

end

@testset "sample() statements" begin

function foo()
    a = sample([:h, :e, :l, :l, :o], weights([1, 1, 5, 1, 1]))
    b = sample(collect("world"), (10, 10))
end

genfoo = genify(foo)

# Check types
trace = simulate(genfoo, ())
@test trace[:a] isa Symbol
@test trace[:b] isa Matrix{Char}

# Test array sizes
trace = simulate(genfoo, ())
@test size(trace[:b]) == (10, 10)

# Test ranges
@test trace[:a] in [:h, :e, :l, :l, :o]
@test all(x in "world" for x in trace[:b])

# Test generate with constraints
trace, weight = generate(genfoo, (), choicemap(:a => :l))
@test weight ≈ log(6/9)
trace, weight = generate(genfoo, (), choicemap(:b => fill('w', (10, 10))))
@test weight ≈ log(1/5) * 100

end

@testset "Primitives with RNGs" begin

function foo(n::Int)
    rng = MersenneTwister(0)
    a = rand(rng, Int8)
    b = rand(rng, Float32, n)
    c = rand(rng, "world")
    d = rand(rng, Wishart(5, [1.0 0.0; 0.0 1.0]))
    e = randexp(rng)
    f = randn(rng)
    g = randperm(rng, n)
    h = randcycle(rng, n)
end

genfoo = genify(foo, Int)

# Check types
trace = simulate(genfoo, (10,))
@test trace[:a] isa Int8
@test trace[:b] isa Vector{Float32}
@test trace[:c] isa Char
@test trace[:d] isa Matrix{Float64}
@test trace[:e] isa Float64
@test trace[:f] isa Float64
@test trace[:g] isa Vector{Int}
@test trace[:h] isa Vector{Int}

end

@testset "In-place primitives" begin

function foo(A::Vector{Float64}, B::Vector{Int})
    a = rand!(A)
    b = randn!(A)
    c = randexp!(A)
    d = randperm!(B)
    e = randcycle!(B)
    f = shuffle!(B)
    g = sample!(B, A)
end

genfoo = genify(foo, Vector{Float64}, Vector{Int})

# Check that arrays were modified
A, B = Vector{Float64}(undef, 10), Vector{Int}(undef, 10)
trace = simulate(genfoo, (A, B))
@test A == trace[:g]
@test B == trace[:f]

# Check types
@test trace[:a] isa Vector{Float64}
@test trace[:b] isa Vector{Float64}
@test trace[:c] isa Vector{Float64}
@test trace[:d] isa Vector{Int}
@test trace[:e] isa Vector{Int}
@test trace[:f] isa Vector{Int}
@test trace[:g] isa Vector{Int}
@test trace[] isa Vector{Float64}

# Test ranges
@test all(0 .<= trace[:a] .<= 1)
@test all(-Inf .<= trace[:b] .<= Inf)
@test all(0 .<= trace[:c] .<= Inf)
@test isperm(trace[:d])
@test Genify.iscycle(trace[:e])
@test isperm(trace[:f])
@test all(1 .<= trace[:g] .<= 10)

end

end

ext: jl


text:
# Invertible CNN layer from Putzky & Welling (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test

# Input
nx = 28
ny = 28
n_in = 8
n_hidden = 8
batchsize = 2

# Input images
X = randn(Float32, nx, ny, n_in, batchsize)
X0 = randn(Float32, nx, ny, n_in, batchsize)
dX = X - X0

# Invertible layers
L = CouplingLayerIRIM(nx, ny, n_in, n_hidden, batchsize)
L01 = CouplingLayerIRIM(nx, ny, n_in, n_hidden, batchsize)
L02 = CouplingLayerIRIM(nx, ny, n_in, n_hidden, batchsize)

###################################################################################################
# Test invertibility

X_ = L.inverse(L.forward(X))
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1e-2)

X_ = L.forward(L.inverse(X))
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1e-2)

###################################################################################################
# Gradient tests

# Loss Function
function loss(L, X, Y)
    Y_ = L.forward(X)
    ΔY = Y_ - Y

    f = .5f0*norm(ΔY)^2
    ΔX = L.backward(ΔY, Y_)[1]

    # Pass back gradients w.r.t. input X and from the residual block and 1x1 conv. layer
    return f, ΔX, L.C.v1.grad, L.C.v2.grad, L.C.v3.grad, L.RB.W1.grad, L.RB.W2.grad, L.RB.W3.grad
end

# Gradient test w.r.t. input X0
Y = L.forward(X)
f0, ΔX = loss(L, X0, Y)[1:2]
h = 0.1f0
maxiter = 5
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test invertible layer\n")
for j=1:maxiter
    f = loss(L, X0 + h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test w.r.t. weights of residual block
Y = L.forward(X)
Lini = deepcopy(L02)
dW1 = L.RB.W1.data - L02.RB.W1.data
dW2 = L.RB.W2.data - L02.RB.W2.data
dW3 = L.RB.W3.data - L02.RB.W3.data

f0, ΔX, Δv1, Δv2, Δv3, ΔW1, ΔW2, ΔW3 = loss(L02, X, Y)
h = 0.1f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test invertible layer\n")
for j=1:maxiter
    L02.RB.W1.data = Lini.RB.W1.data + h*dW1
    L02.RB.W2.data = Lini.RB.W2.data + h*dW2
    L02.RB.W3.data = Lini.RB.W3.data + h*dW3
    f = loss(L02, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW1, ΔW1) - h*dot(dW2, ΔW2) - h*dot(dW3, ΔW3))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test w.r.t. 1x1 conv weights
Y = L.forward(X)
Lini = deepcopy(L01)
dv1 = L.C.v1.data - L01.C.v1.data
dv2 = L.C.v2.data - L01.C.v2.data
dv3 = L.C.v3.data - L01.C.v3.data

f0, ΔX, Δv1, Δv2, Δv3, ΔW1, ΔW2, ΔW3 = loss(L01, X, Y)
h = 0.1f0
maxiter = 4
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)

print("\nGradient test invertible layer\n")
for j=1:maxiter
    L01.C.v1.data = Lini.C.v1.data + h*dv1
    L01.C.v2.data = Lini.C.v2.data + h*dv2
    L01.C.v3.data = Lini.C.v3.data + h*dv3
    f = loss(L01, X, Y)[1]
    err5[j] = abs(f - f0)
    err6[j] = abs(f - f0 - h*dot(dv1, Δv1) - h*dot(dv2, Δv2) - h*dot(dv3, Δv3))
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
## 19.jl
# 各行の1コラム目の文字列の出現頻度を求め，出現頻度の高い順に並べる

f = open("./dat/hightemp.txt")
d = Dict()
for line in readlines(f)
  col1 = split(chomp(line), "\t")[1]
  if haskey(d,col1)
    d[col1] += 1
  else
    d[col1] = 1
  end
end

for (key,value) in reverse(sort(collect(d), by = x-> x[2]))
  @printf "%s : %d\n" key value
end
close(f)


text:
function create_agent(
    ::Val{:DDPG};
    env,
    seed=37,
    kwargs...
)
    rng = StableRNG(seed)

    ns = length(state(env))
    na = length(action_space(env))

    init = glorot_uniform(rng)

    create_actor() = Chain(
        Dense(ns, 30, relu; init = init),
        Dense(30, 30, relu; init = init),
        Dense(30, na; init = init),
        softmax
    )

    create_critic() = Chain(
        Dense(ns + na, 30, relu; init = init),
        Dense(30, 30, relu; init = init),
        Dense(30, 1; init = init),
    )

    return Agent(
        policy = DDPGPolicy(
            behavior_actor = NeuralNetworkApproximator(
                model = create_actor(),
                optimizer = ADAM(),
            ),
            behavior_critic = NeuralNetworkApproximator(
                model = create_critic(),
                optimizer = ADAM(),
            ),
            target_actor = NeuralNetworkApproximator(
                model = create_actor(),
                optimizer = ADAM(),
            ),
            target_critic = NeuralNetworkApproximator(
                model = create_critic(),
                optimizer = ADAM(),
            ),
            γ = 0.99f0,
            ρ = 0.995f0,
            na = na,
            batch_size = 64,
            start_steps = 1000,
            start_policy = RandomPolicy(Space([-1.0..1.0 for i in 1:na]); rng = rng),
            update_after = 1000,
            update_freq = 1,
            act_limit = 1.0,
            act_noise = 0.01,
            rng = rng,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 10000,
            state = Vector{Float32} => (ns,),
            action = Float32 => (na,),
        ),
    )
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
@testset "reported issues" begin
    @testset "barplot nl in text (#102)" begin
        p = barplot(
            ["I need a\nbreak", "I don't but could", :Hello, :Again],
            [30, 40, 20, 10],
        )
        test_ref("references/issues/barplot_nl.txt", @show_col(p))
    end

    @testset "plot hang (tan → ∞) (#129)" begin
        p = lineplot([cos, sin, tan], -π / 2, 2π)
        test_ref("references/issues/cos_sin_tan.txt", @show_col(p))
    end

    @testset "integer edges (#139)" begin
        p_int = histogram(fit(Histogram, rand(10), 0:3))
        p_float = histogram(fit(Histogram, rand(10), 0.0:1.0:3.0))
        @test @print_col(p_int) == @print_col(p_float)
    end

    @testset "isolated colorbar (#169)" begin
        p = Plot([0], [0], colorbar = true, colormap = :cividis, width = 0)
        test_ref("references/issues/isolated_colorbar.txt", @show_col(p))
    end

    @testset "steps (#211)" begin
        x = 1:5
        y = [1.0, 2.0, 3.0, 2.0, 1.0]
        p = scatterplot(x, y, marker = :circle)
        x = [1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0]
        y = [1.0, 2.0, 2.0, 3.0, 3.0, 2.0, 2.0, 1.0, 1.0]
        p = lineplot!(p, x, y)
        test_ref("references/issues/steps.txt", @show_col(p))
    end
end


text:
function game_of_life(board, steps)
    return board
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(cot(x)^4, x) == :(x+-1//3*cot(x)^3+cot(x))
@test integrate(x^-4*(1+x^2)^-1, x) == :(x^-1+-1//3*x^-3+arctan(x))
@test integrate(x^-1//2*(x+x^2), x) == :(2//3*x^3//2+2//5*x^5//2)
@test integrate(cos(x), x) == :(sin(x))
@test integrate(x*exp(x^2), x) == :((1/2)*exp(x^2))
@test integrate(sec(x)^2*tan(x), x) == :((1/2)*sec(x)^2)
@test integrate(x*(1+x^2)^(1/2), x) == :(1//3*(1+x^2)^3//2)
@test integrate(exp(x)*sin(x), x) == :((1/2)*exp(x)*sin(x)+-1//2*cos(x)*exp(x))
@test integrate(csc(x)^2*sin(x)^-2*cos(x), x) == :(-1//3*csc(x)^3)
@test integrate(sin(exp(x)), x) == :(sinintegral(exp(x)))
@test integrate(y^-1*sin(y), y) == :(sinintegral(y))
@test integrate(exp(x)+sin(x), x) == :(-1*cos(x)+exp(x))
@test integrate(2*x^2*exp(x^2)+exp(x^2), x) == :(x*exp(x^2))
@test integrate((x+exp(x))^2, x) == :((1/2)*exp(2x)+-2*exp(x)+1//3*x^3+2*x*exp(x))
@test integrate(x^2+2*exp(x)+exp(2x), x) == :((1/2)*exp(2x)+2*exp(x)+1//3*x^3)
@test integrate(cos(x)*sin(x), x) == :((1/2)*sin(x)^2)
@test integrate(x*exp(x^2), x) == :((1/2)*exp(x^2))
@test integrate(x*(1+x^2)^(1/2), x) == :(1//3*(1+x^2)^3//2)
@test integrate((1+exp(x))^-1*exp(x), x) == :(log(1+exp(x)))
@test integrate(x^3//2, x) == :(2//5*x^5//2)
@test integrate(cos(3+2x), x) == :((1/2)*sin(3+2x))
@test integrate(2*y*z*exp(2x), x) == :(y*z*exp(2x))
@test integrate(cos(exp(x))^2*exp(x)*sin(exp(x)), x) == :(-1//3*cos(exp(x))^3)
@test integrate(x*(1+x)^(1/2), x) == :(-2//3*(1+x)^3//2+2//5*(1+x)^5//2)
@test integrate((-1+x^4)^-1, x) == :(-1//2*arctan(x)+-1//2*arctanh(x))
@test integrate((2+3*exp(2x))^-1*exp(x), x) == :(1//6*6^(1/2)*arctan((1/2)*6^(1/2)*exp(x)))
@test integrate((A+B*exp(4x))^-1*exp(2x), x) == :((1/2)*A^-1//2*B^-1//2*arctan(A^-1//2*B^(1/2)*exp(2x)))
@test integrate((1+exp(x))^-1*exp(1+x), x) == :(ℯ*log(1+exp(x)))
@test integrate(10^x*exp(x), x) == :((10ℯ)^x*(1+log(10))^-1)
@test integrate(x^3*sin(x^2), x) == :((1/2)*sin(x^2)+-1//2*x^2*cos(x^2))
@test integrate(x^7*(1+x^12)^-1, x) == :(-1//12*log(1+x^4)+1//24*log(1+x^8+-1*x^4)+-1//12*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2*x^4)))
@test integrate(x^(3a)*sin(x^(2a)), x) == :(-1//4*im*a^-1*x^(1+3a)*(im*x^(2a))^(-1//2*a^-1*(1+3a))*SpecialFunctions.gamma(3//2+(1/2)*a^-1,im*x^(2a))+1//4*im*a^-1*x^(1+3a)*(-1*im*x^(2a))^(-1//2*a^-1*(1+3a))*SpecialFunctions.gamma(3//2+(1/2)*a^-1,-1*im*x^(2a)))
@test integrate(cos(x^(1/2)), x) == :(2*cos(x^(1/2))+2*x^(1/2)*sin(x^(1/2)))
@test integrate(x*(1+x)^(1/2), x) == :(-2//3*(1+x)^3//2+2//5*(1+x)^5//2)
@test integrate((x^(1/2)+x^1//3)^-1, x) == :(-6*log(1+x^1//6)+-3*x^1//3+2*x^(1/2)+6*x^1//6)
@test integrate(((3+2x)^-1*(1+x))^(1/2), x) == :((1/2)*(1+x)^(1/2)*(3+2x)^(1/2)+-1//4*2^(1/2)*arcsinh(2^(1/2)*(1+x)^(1/2)))
@test integrate(x^4*(1+-1*x^2)^-5//2, x) == :(-1*x*(1+-1*x^2)^-1//2+1//3*x^3*(1+-1*x^2)^-3//2+arcsin(x))
@test integrate(x^(1/2)*(1+x)^5//2, x) == :(-5//64*arcsinh(x^(1/2))+1//4*x^3//2*(1+x)^5//2+5//24*x^3//2*(1+x)^3//2+5//32*x^3//2*(1+x)^(1/2)+5//64*x^(1/2)*(1+x)^(1/2))
@test integrate(x^4*(1+-1*x^2)^-5//2, x) == :(-1*x*(1+-1*x^2)^-1//2+1//3*x^3*(1+-1*x^2)^-3//2+arcsin(x))
@test integrate((1+-1*y^2)^-1*(A^2+B^2+-1*B^2*y^2)^(1/2), y) == :(A*arctanh(A*y*(A^2+B^2+-1*B^2*y^2)^-1//2)+B*arctan(B*y*(A^2+B^2+-1*B^2*y^2)^-1//2))
@test integrate(sin(x)^2, x) == :((1/2)*x+-1//2*cos(x)*sin(x))
@test integrate((1+cos(x))^-1, x) == :((1+cos(x))^-1*sin(x))
@test integrate(x*exp(x), x) == :(-1*exp(x)+x*exp(x))
@test integrate(x*(1+x)^-2*exp(x), x) == :((1+x)^-1*exp(x))
@test integrate((1+2*x^2)*exp(x^2), x) == :(x*exp(x^2))
@test integrate(exp(x^2), x) == :((1/2)*pi^(1/2)*erfi(x))
@test integrate(x^-1*exp(x), x) == :(expintegralei(x))
@test integrate(x*(1+x^3)^-1, x) == :(-1//3*log(1+x)+1//6*log(1+x^2+-1x)+-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2x)))
@test integrate((A^4+x^2*(B^2+-1*A^2)+-1*A^2*B^2)^-1, x) == :(A^-1*(A^2+-1*B^2)^-1*arctanh(x*A^-1))
@test integrate(x*log(x), x) == :(-1//4*x^2+(1/2)*x^2*log(x))
@test integrate(x^2*arcsin(x), x) == :(-1//9*(1+-1*x^2)^3//2+1//3*(1+-1*x^2)^(1/2)+1//3*x^3*arcsin(x))
@test integrate((1+x^2+2x)^-1, x) == :(-1*(1+x)^-1)
@test integrate((1+log(x))^-2*log(x), x) == :(x*(1+log(x))^-1)
@test integrate(x^-1*(1+log(x)^2)^-1, x) == :(arctan(log(x)))
@test integrate(log(x)^-1, x) == :(logintegral(x))
@test integrate(x*(cos(x)+sin(x)), x) == :(x*sin(x)+-1*x*cos(x)+cos(x)+sin(x))
@test integrate((x+exp(x))*exp(-1x), x) == :(x+-1*exp(-1x)+-1*x*exp(-1x))
@test integrate(x*(1+exp(x))^2, x) == :((1/2)*x^2+-2*exp(x)+-1//4*exp(2x)+(1/2)*x*exp(2x)+2*x*exp(x))
@test integrate(x*cos(x), x) == :(x*sin(x)+cos(x))
@test integrate(cos(x^(1/2)), x) == :(2*cos(x^(1/2))+2*x^(1/2)*sin(x^(1/2)))
@test integrate(x*cos(x), x) == :(x*sin(x)+cos(x))
@test integrate(x*log(x)^2, x) == :(1//4*x^2+(1/2)*x^2*log(x)^2+-1//2*x^2*log(x))
@test integrate((1+sin(x)^3)*cos(x), x) == :(1//4*sin(x)^4+sin(x))
@test integrate(x^-1*(1+log(x)^2)^-1, x) == :(arctan(log(x)))
@test integrate((1+arcsin(x)^2)^-1*(1+-1*x^2)^-1//2, x) == :(arctan(arcsin(x)))
@test integrate((cos(x)+sin(x))^-1*sin(x), x) == :((1/2)*x+-1//2*log(cos(x)+sin(x)))
@test integrate(-1*(1+-1*y^2)^-1*(A^2+B^2*(1+-1*y^2))^(1/2), y) == :(-1*A*arctanh(A*y*(A^2+B^2+-1*B^2*y^2)^-1//2)+-1*B*arctan(B*y*(A^2+B^2+-1*B^2*y^2)^-1//2))
@test integrate(B^-1*(1+-1*B^-2*sin(z)^2*(A^2+B^2))^-1*cos(z)^2*(-1*A^2+-1*B^2), z) == :(-1*A*arctanh(A*B^-1*tan(z))+-1*B*z)
@test integrate(B^-1*(1+w^2)^-2*(1+-1*B^-2*w^2*(1+w^2)^-1*(A^2+B^2))^-1*(-1*A^2+-1*B^2), w) == :(-1*A*arctanh(A*w*B^-1)+-1*B*arctan(w))
@test integrate(-1*B*(1+w^2)^-1*(B^2+-1*A^2*w^2)^-1*(A^2+B^2), w) == :(-1*A*arctanh(A*w*B^-1)+-1*B*arctan(w))
@test integrate(x^4*(1+-1*x^2)^-5//2, x) == :(-1*x*(1+-1*x^2)^-1//2+1//3*x^3*(1+-1*x^2)^-3//2+arcsin(x))
@test integrate(cos(y)^-4*sin(y)^4, y) == :(y+-1*tan(y)+1//3*tan(y)^3)
@test integrate(z^4*(1+z^2)^-1, z) == :(-1z+1//3*z^3+arctan(z))
@test integrate((1+2*x^2)*exp(x^2), x) == :(x*exp(x^2))
@test integrate((1+x^2)^-2*(1+x^3+2*x^6+4*x^2+5*x^4)*exp(x^2), x) == :(x*exp(x^2)+(2+2*x^2)^-1*exp(x^2))
@test integrate(exp(-1)*exp(-1x), x) == :(-1*exp(-1+-1x))
@test integrate((x+x^-1)*log(x), x) == :((1/2)*log(x)^2+-1//4*x^2+(1/2)*x^2*log(x))
@test integrate(x*(1+x^4)^-1, x) == :((1/2)*arctan(x^2))
@test integrate(x^5*(1+x^4)^-1, x) == :((1/2)*x^2+-1//2*arctan(x^2))
@test integrate((1+tan(x)^2)^-1, x) == :((1/2)*x+(1/2)*cos(x)*sin(x))
@test integrate(x^4*(1+-1*x^2)^-5//2, x) == :(-1*x*(1+-1*x^2)^-1//2+1//3*x^3*(1+-1*x^2)^-3//2+arcsin(x))
@test integrate(-1*x^2*(1+-1*x^2)^-3//2, x) == :(-1*x*(1+-1*x^2)^-1//2+arcsin(x))
@test integrate(exp(x)*sin(x), x) == :((1/2)*exp(x)*sin(x)+-1//2*cos(x)*exp(x))
@test integrate(x^-1, x) == :(log(x))
@test integrate((1+sec(t)^2+3*tan(t))^-1*sec(2t), t) == :(-1*(2+2*tan(t))^-1+-1//4*log(cos(t)+sin(t))+-1//12*log(-1*sin(t)+cos(t))+1//3*log(2*cos(t)+sin(t)))
@test integrate(sec(x)^-2, x) == :((1/2)*x+(1/2)*cos(x)*sin(x))
@test integrate(x^-1//2*(1+x^2), x) == :(2*x^(1/2)+2//5*x^5//2)
@test integrate(x*(5+x^2+2x)^-1//2, x) == :((5+x^2+2x)^(1/2)+-1*arcsinh(1/2+(1/2)*x))
@test integrate(sin(x)^2*cos(x), x) == :(1//3*sin(x)^3)
@test integrate((1+exp(x))^-1*exp(x), x) == :(log(1+exp(x)))
@test integrate((1+exp(x))^-1*exp(2x), x) == :(-1*log(1+exp(x))+exp(x))
@test integrate((1+-1*cos(x))^-1, x) == :(-1*(1+-1*cos(x))^-1*sin(x))
@test integrate(sec(x)^2*tan(x), x) == :((1/2)*sec(x)^2)
@test integrate(x*log(x), x) == :(-1//4*x^2+(1/2)*x^2*log(x))
@test integrate(cos(x)*sin(x), x) == :((1/2)*sin(x)^2)
@test integrate((-1*x^2+2x)^-1//2*(1+x), x) == :(-1*(-1*x^2+2x)^(1/2)+-2*arcsin(1+-1x))
@test integrate(2*(2+3*exp(2x))^-1*exp(x), x) == :(1//3*6^(1/2)*arctan((1/2)*6^(1/2)*exp(x)))
@test integrate(x^4*(1+-1*x^2)^-5//2, x) == :(-1*x*(1+-1*x^2)^-1//2+1//3*x^3*(1+-1*x^2)^-3//2+arcsin(x))
@test integrate((1+exp(4x))^-1*exp(6x), x) == :((1/2)*exp(2x)+-1//2*arctan(exp(2x)))
@test integrate(log(2+3*x^2), x) == :(-2x+x*log(2+3*x^2)+2//3*6^(1/2)*arctan((1/2)*x*6^(1/2)))
@test integrate(r^-1*(-1*a^2+2*H*r^2)^-1//2, x) == :(x*r^-1*(-1*a^2+2*H*r^2)^-1//2)
@test integrate(r^-1*(-1*a^2+-1*e^2+2*H*r^2)^-1//2, x) == :(x*r^-1*(-1*a^2+-1*e^2+2*H*r^2)^-1//2)
@test integrate(r^-1*(-1*a^2+-2*K*r^4+2*H*r^2)^-1//2, x) == :(x*r^-1*(-1*a^2+-2*K*r^4+2*H*r^2)^-1//2)
@test integrate(r^-1*(-1*a^2+-1*e^2+-2*K*r^4+2*H*r^2)^-1//2, x) == :(x*r^-1*(-1*a^2+-1*e^2+-2*K*r^4+2*H*r^2)^-1//2)
@test integrate(r^-1*(-1*a^2+-2*K*r+2*H*r^2)^-1//2, x) == :(x*r^-1*(-1*a^2+-2*r*(K+-1*H*r))^-1//2)
@test integrate(r*(-1*a^2+2*ℯ*r^2)^-1//2, x) == :(r*x*(-1*a^2+2*ℯ*r^2)^-1//2)
@test integrate(r*(-1*a^2+-1*e^2+2*ℯ*r^2)^-1//2, x) == :(r*x*(-1*a^2+-1*e^2+2*ℯ*r^2)^-1//2)
@test integrate(r*(-1*a^2+-2*K*r^4+2*ℯ*r^2)^-1//2, x) == :(r*x*(-1*a^2+-2*K*r^4+2*ℯ*r^2)^-1//2)
@test integrate(r*(-1*a^2+-1*e^2+-2*K*r^4+2*ℯ*r^2)^-1//2, x) == :(r*x*(-1*a^2+-1*e^2+-2*K*r^4+2*ℯ*r^2)^-1//2)

max_stars_repo_path: test/0 Independent test suites/Moses Problems.jl


text:
# Benders Decomposition
# 1. BendersOptim.milp: for MILP with Sub and Ray Problems
# 2. BendersOptim.lshaped: for stochastic programming
# 0. Linear programming in matrix
# Version: 5.1
# Author: Edward J. Xu, edxu96@outlook.com
# Date: April 5th, 2019
module Benders
    using JuMP
    using GLPKMathProgInterface
    using PrettyTables
    # 1.BendersOptim.milp: for MILP with Sub and Ray Problems
    include("Benders_milp.jl")
    # 2. BendersOptim.lshaped: for stochastic programming
    include("Benders_lshaped.jl")
    # 0. Linear programming
    include("lp.jl")
end

max_stars_repo_name: eitt/BendersDecomposition


text:
using Parameters: @with_kw
using MLDataUtils
using Flux
using Flux: OneHotArray
using Zygote
using CUDA
import Flux.Optimise: update!, ADAM
using Flux: logitcrossentropy, logitbinarycrossentropy
using Transformers.HuggingFace

"""
@with_kw struct HyperParams
    batchsize::Int=8
    eos_token::Int=50527
    pad_token::Int=0
    embed_size::Int=768
    classification_type::String="Binary"
    target_class=1
    num_classes::Int=2
    cached::Bool=true
    lr::Float64=1e-6
    epochs::Int=20
    log_interval::Int=100
end

Structure for `HyperParams` to select hyperparamters for training the discriminator model.

"""
@with_kw struct HyperParams
    batchsize::Int=8
    eos_token::Int=50527
    pad_token::Int=0
    embed_size::Int=768
    classification_type::String="Binary"
    target_class=1
    num_classes::Int=2
    cached::Bool=true
    lr::Float64=1e-6
    epochs::Int=20
    log_interval::Int=100
end

"""
    train!(discrim, data_loader; args=args)

Train the discriminator using the provided `data_loader` training data and arguments `args` provided.
"""
function train!(discrim, data_loader; args=nothing)
    if args==nothing
        @warn "No args provided, using the default arguments from HyperParams structure."
        args = HyperParams()
    end

    global device
    opt = ADAM(args.lr)
    ps = params(discrim)

    function loss_fn(x, y)
        if args.classification_type=="Binary"
            return logitbinarycrossentropy(x, y)
        else
            return logitcrossentropy(x, y)
        end
    end

    print("Training starting.........\n\n")
    for epoch in 1:args.epochs
        avg_loss = 0.0
        step =0
        for (data_x, data_y) in data_loader
            if args.cached == true
                data_x, data_y = (data_x, data_y) |> device
                mask = nothing
            else
               data_x, data_y, mask = data_preprocess(data_x, data_y, args.classification_type, args.num_classes; args=args) |> device
            end

            loss, back = Zygote.pullback(ps) do
                        logits = discrim(data_x, mask; args = args)
                        loss_fn(logits, data_y)
                      end
            grads = back(1f0)
            update!(opt, ps, grads)
            avg_loss += loss
            step+=1

            if step%args.log_interval==0
                print("Epoch #", epoch, " : Step #", step, " : NLL Loss: ", avg_loss/step, "\n")
            end
        end
        print("Loss Epoch #", epoch, " : ", step, " NLL Loss: ", avg_loss/step, "\n")
        print("\n")
    end
end

"""
    test!(discrim, data_loader; args=nothing)

Test the discriminator on test data provided using `data_loader`, based on Accuracy and NLL Loss.

"""
function test!(discrim, data_loader; args=nothing)
    global device
    function loss_fn(x, y)
        if args.classification_type=="Binary"
            return logitbinarycrossentropy(x, y)
        else
            return logitcrossentropy(x, y)
        end
    end

    avg_loss = 0
    if args.classification_type=="Binary"
        preds = Array{Float32}(undef, 1, 0)
        ys = Array{Float32}(undef, 1, 0)
    else
        preds = Array{Float32}(undef, args.num_classes, 0)
        ys = Array{Float32}(undef, args.num_classes, 0)
    end
    steps=0
    print("Testing.......\n\n")
    for (data_x, data_y) in data_loader
        if args.cached == true
            data_x, data_y = (data_x, data_y) |> device
            mask = nothing
        else
           data_x, data_y, mask = data_preprocess(data_x, data_y, args.classification_type, args.num_classes; args=args) |> device
        end

        logits = discrim(data_x, mask; args = args)
        loss = loss_fn(logits, data_y)
        if args.classification_type !="Binary"
            logits = softmax(logits)
        else
            logits = sigmoid.(logits)
        end
        logits, data_y = (logits, data_y) |> cpu
        preds= cat(preds, logits, dims=2)
        ys = cat(ys, data_y, dims=2)
        avg_loss += loss
        steps+=1
    end
    loss = avg_loss/steps
    if args.classification_type=="Binary"
        acc = binary_accuracy(preds, ys)
    else
        acc = categorical_accuracy(preds, ys)
    end
    print("Test Results: NLL Loss: ", loss, ", Accuracy: ", acc)
end

"""
    train_discriminator(text, labels, batchsize::Int=8, classification_type::String="Binary", num_classes::Int=2; model="gpt2", cached::Bool=true, discrim=nothing, tokenizer=nothing, truncate=true, max_length=256, train_size::Float64=0.9, lr::Float64=1e-5, epochs::Int=10, args=nothing)

Function to train discriminator for provided `text` and target `labels`, based on set of function paramters provided. Returns `discrim` discriminator after training.

Here the `cached=true` allows cacheing of contexualized embeddings (forward pass) in GPT2 model, as the model itself is non-trainable. This reduces the time of training effectively as the forward pass through GPT2 model is to be done only once.

# Example

Consider a Multiclass classification problem with class size of 5, it can trained on `text` and `labels` vectors using:

```julia
train_discriminator(text, labels, 16, "Multiclass", 5)
```

"""
function train_discriminator(text, labels, batchsize::Int=8, classification_type::String="Binary", num_classes::Int=2; model="gpt2", cached::Bool=true, discrim=nothing, tokenizer=nothing, truncate=true, max_length=256, train_size::Float64=0.9, lr::Float64=1e-5, epochs::Int=10, args=nothing)
    global device
    # Right now only support gpt2
    if tokenizer == nothing && model=="gpt2"
        model = hgf"gpt2:lmheadmodel" |> device
        config = load_config("gpt2")
        tokenizer = load_pretrained_tokenizer(GPT2)
    end

    if args == nothing
        args = HyperParams(batchsize=batchsize, embed_size=config.n_embd, classification_type=classification_type, num_classes=num_classes, cached=cached, lr=lr, epochs=epochs)
    end

    if discrim == nothing
        if args.classification_type=="Binary"
            class_size=1
        else
            class_size=num_classes
        end
        discrim = get_discriminator(model, class_size=class_size, embed_size=args.embed_size)
    end
    if cached == true
        print("Creating Cache of data and Data Loader...")
        text, labels = shuffleobs((text, labels))
        (train_x, train_y), (test_x, test_y) = splitobs((text, labels); at=train_size)

        train_loader = load_cached_data(discrim, train_x, train_y, tokenizer; truncate=truncate, max_length=max_length, shuffle=true,
                    drop_last=true, classification_type=classification_type, num_classes=num_classes, args=args)

        test_loader = load_cached_data(discrim, test_x, test_y, tokenizer; truncate=truncate, max_length=max_length, shuffle=false,
                    drop_last=true, classification_type=classification_type, num_classes=num_classes, args=args)

        train!(discrim, train_loader; args=args)
        test!(discrim, test_loader; args=args)
    else
        print("Creating Data Loader...")
        text, labels = shuffleobs((text, labels))
        (train_x, train_y), (test_x, test_y) = splitobs((text, labels); at=test_size)

        train_loader = load_data(train_x, train_y, tokenizer; truncate=truncate, max_length=max_length, shuffle=true,
                    drop_last=true, args=args)

        test_loader = load_data(test_x, test_y, tokenizer; truncate=truncate, max_length=max_length, shuffle=false,
                    drop_last=true, args=args)

        train!(discrim, train_loader; args=args)
        test!(discrim, test_loader; args=args)
    end
    return discrim
end


text:
using RandomProcesses

d = CRP(100, 0.01)
x = rand(d)
pdf(d, x)

d = CRP(100, 1.0)
x = rand(d)
pdf(d, x)
shuffle!(x)
pdf(d, x)

d = CRP(100, 10000.0)
x = rand(d)
pdf(d, x)


text:
n4f4ton0f8(x::UInt16) = reinterpret(N0f8, (x % UInt8) * 0x11)
n3f5ton0f8(x::UInt16) = reinterpret(N0f8, (x * 0x083A + 0x80) >> 0x8 % UInt8)
n2f6ton0f8(x::UInt16) = reinterpret(N0f8, (x * 0x040C + 0x80) >> 0x8 % UInt8)

function rationalized_scale(::Type{T}) where {T}
    a = zeros(T, 256)
    for d = 1:255
        for n = 1:d
            i = round(Int, n / d * 255) + 1
            if a[i] == zero(T)
                a[i] = n / d
            end
        end
    end
    return a
end

const n0f8tof32 = rationalized_scale(Float32)
const n0f8tof64 = rationalized_scale(Float64)

hue_u8(hue::Real) = unsafe_trunc(UInt8, round((hue * 0.5f0) % 180))
hue_u8(hue::Integer) = (hue >> 0x1) % UInt8

lstar_u8(l::Real) = unsafe_trunc(UInt8, round(l * 2.55f0))

abstar_u8(ab::Real) = unsafe_trunc(UInt8, round(ab + 128))

ustar_u8(u::Real) = unsafe_trunc(UInt8, round((u + 134) * 255 / 354))

vstar_u8(v::Real) = unsafe_trunc(UInt8, round((v + 140) * 255 / 262))

function ColorTypes._convert(::Type{Cout}, ::Type{C1}, ::Type{C2}, c) where
    {Cout<:HSV, C1<:HSV, C2<:HSV24}

    Cout(comp1(c), comp2(c), comp3(c))
end

function ColorTypes._convert(::Type{Cout}, ::Type{C1}, ::Type{C2}, c) where
    {Cout<:AbstractRGB, C1<:AbstractRGB, C2<:Union{HSV24, HSL24, HLS24}}
    convert(Cout, floattype(C2)(comp1(c), comp2(c), comp3(c)))
end

reinterpret(::Type{UInt16}, c::C) where {C <: AbstractRGB16} = c.color
reinterpret(::Type{C}, x::UInt16) where {C <: AbstractRGB16} = C(x, Val(:no))


text:
struct PhysicalRelation
    relpath::Vector{String}
    columns::JSON3.Array
end
_slots(relpath::Vector{String}) =
    Int[i for (i,v) in enumerate(relpath) if !startswith(v, ':')]

# reverse constructor from TuplesIterator (for convenience at REPL)
function PhysicalRelation(relpath, tuples::TuplesIterator)
    # TODO: there must be a better way to do this than to write+read?
    columns = JSON3.read(JSON3.write([
        [ tup[slot] for tup in tuples ]
        for slot in _slots(relpath)
    ]))
    return PhysicalRelation(relpath, columns)
end

_num_tuples(r::PhysicalRelation) = max(length(r.columns[1]), 1)

function tuples(r::PhysicalRelation)
    len = _num_tuples(r)
    tuples = (
        row_getter(i)
        for row_getter in (RAI._make_getrow(r.relpath, r.columns),)
        for i in 1:_num_tuples(r)
    )
    return TuplesIterator(len, tuples)
end
function tuples(rs::Vector{PhysicalRelation})
    len = sum((_num_tuples(r) for r in rs), init=0)
    tuples = (
        row_getter(i)
        for r in sort(rs, by=r->r.relpath)
        for row_getter in (RAI._make_getrow(r.relpath, r.columns),)
        for i in 1:_num_tuples(r)
    )
    return TuplesIterator(len, tuples)
end

function Base.show(io::IO, r::PhysicalRelation)
    print(io, "$PhysicalRelation($(r.relpath), ")
    # recur_io = IOContext(io, :SHOWN_SET => t)
    limit = get(io, :limit, false)::Bool
    if limit  # meaning we're inside another context like an Array,
        # Also set :compact, since the Tuples vector can be quite big
        if !haskey(io, :compact)
            io = IOContext(io, :compact => true)
        end
    end
    show(io, tuples(r))
    print(io, ")")
end


text:
"""
    scatterplot(x, y; kwargs...)

Description
============

Draws the given points on a new canvas.

The first (optional) vector `x` should contain the horizontal
positions for all the points. The second vector `y` should then
contain the corresponding vertical positions respectively. This
means that the two vectors must be of the same length and
ordering.

Usage
======

    scatterplot([x], y; name = "", marker = :pixel, title = "", xlabel = "", ylabel = "", labels = true, border = :solid, margin = 3, padding = 1, color = :auto, xlim = (0, 0), ylim = (0, 0), canvas = BrailleCanvas, grid = true)

Arguments
==========

- **`x`** : Optional. The horizontal position for each point.
  If omitted, the axes of `y` will be used as `x`.

- **`y`** : The vertical position for each point.

- **`name`** : Annotation of the current drawing to be displayed on the right

- **`marker`** : Choose a marker from $(keys(MARKERS)), a `Char`, a unit length `String` or a vector of these

$DOC_PLOT_PARAMS

- **`height`** : Number of character rows that should be used for plotting.

- **`xlim`** : Plotting range for the x axis.
  `(0, 0)` stands for automatic.

- **`ylim`** : Plotting range for the y axis.
  `(0, 0)` stands for automatic.

- **`canvas`** : The type of canvas that should be used for drawing.

- **`grid`** : If `true`, draws grid-lines at the origin.

Returns
========

A plot object of type `Plot{T<:Canvas}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
=========

```julia-repl
julia> scatterplot(randn(50), randn(50), title = "My Scatterplot")
                   My Scatterplot
      ┌────────────────────────────────────────┐
    3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡆⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⡇⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠡⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠁⠀⡇⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠐⠄⠀⠀⠀⠀⡀⠀⠂⡇⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠤⠤⠤⠤⠤⠤⠴⠤⠤⠤⠬⠤⠬⠤⠤⠤⡧⢥⡤⠤⠤⡤⠤⠤⠤⠬⠤⠤⠤⠤⠤⠤⠤⠤⡤⠤⠤⠤⠤⠄│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⡗⠐⠨⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⡇⠀⠐⠀⠨⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠄⠀⠀⠀⠀⠀⡇⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠁⠂⠰⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      └────────────────────────────────────────┘
      -2                                       3
```

See also
=========

[`Plot`](@ref), [`lineplot`](@ref), [`stairs`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function scatterplot(
    x::AbstractVector, y::AbstractVector;
    canvas::Type = BrailleCanvas, color::Union{UserColorType,AbstractVector} = :auto,
    marker::Union{MarkerType,AbstractVector} = :pixel, name = "", kw...
)
    new_plot = Plot(x, y, canvas; kw...)
    scatterplot!(new_plot, x, y; color = color, name = name, marker = marker)
end

scatterplot(y::AbstractVector; kw...) = scatterplot(axes(y, 1), y; kw...)

function scatterplot!(
    plot::Plot{<:Canvas}, x::AbstractVector, y::AbstractVector;
    color::Union{UserColorType,AbstractVector} = :auto,
    marker::Union{MarkerType,AbstractVector} = :pixel, name = ""
)
    color = (color == :auto) ? next_color!(plot) : color
    name == "" || label!(plot, :r, string(name), color isa AbstractVector ? color[1] : color)
    if marker ∈ (:pixel, :auto)
        points!(plot, x, y, color)
    else
        for (xi, yi, mi, ci) in zip(x, y, iterable(marker), iterable(color))
            annotate!(plot, xi, yi, char_marker(mi); color=ci)
        end
    end
    plot
end

scatterplot!(plot::Plot{<:Canvas}, y::AbstractVector; kw...) = scatterplot!(
  plot, axes(y, 1), y; kw...
)


text:
@static if VERSION < v"0.7-"
    using Base.Test
    using Compat: range
else
    using Test
end

using WeightedNorms
@static if VERSION >= v"1.0-"
    if !(@isdefined WN)
        const WN = WeightedNorms
    end
else
    if !isdefined(:WN)
        const WN = WeightedNorms
    end
end

const testGridLength = 100
const tolerance = 1e-10
testgridArchetype = range(0, stop=1, length=testGridLength)
testGrids = [testgridArchetype, collect(testgridArchetype)]

@testset "Grid type=$(typeof(testgrid))" for testgrid in testGrids
    zerofunc(t) = zero(t)
    constfunc(t) = one(t)
    linfunc(t) = t
    sqrfunc(t) = (t^2) / 2
    cubefunc(t) = (t^3) / 6
    signchangingfunc(t) = t - 1/2
    testFunctions = [zerofunc, constfunc, linfunc, sqrfunc, cubefunc, signchangingfunc, sin, exp]

    @testset "Equality, norm=$func" for func in WN.exportedNorms
        @testset "Equality for Equally Spaced Grids, test function=$testfunc" for testfunc in testFunctions
            testdata = map(testfunc, testgrid)
            f1 = @inferred func(testdata, testgrid)
            f2 = @inferred func(testdata, testgrid[2] - testgrid[1])
            @test isapprox(f1, f2, atol = tolerance)
        end
    end

    @testset "Known Exact Norms" for func in WN.exportedNorms
        @testset "Norm of Zero is Zero" begin
            # Norm of zero is zero
            testdata = map(zerofunc, testgrid)
            f1 = @inferred func(testdata, testgrid)
            f2 = @inferred func(testdata, testgrid[2] - testgrid[1])
            @test isapprox(f1, 0.0)
            @test isapprox(f2, 0.0)
        end
        #   Norm is non-negative
        @testset "Positivity, func=$testfunc"for testfunc in testFunctions
            testdata = map(testfunc, testgrid)
            f1 = @inferred func(testdata, testgrid)
            f2 = @inferred func(testdata, testgrid[2] - testgrid[1])
            @test f1 >= 0.0
            @test f2 >= 0.0
        end
    end

    # Other known values
    knownValues = (
                   (constfunc, WN.norm_l2, 1.0),
                   (constfunc, WN.normsq_w21_only, 0.0),
                   (constfunc, WN.normsq_w22_only, 0.0),
                   (constfunc, WN.norm_w21_only, 0.0),
                   (constfunc, WN.norm_w22_only, 0.0),
                   (linfunc, WN.norm_w22_only, 0.0),
                   (linfunc, WN.normsq_w22_only, 0.0),
                   (linfunc, WN.norm_w21_only, 1.0),
                   #(sqrfunc, WN.norm_w22_only, 1.0)
                   )
    @testset "Known Values, testfunc=$testfunc, normfunc=$normfunc" for (testfunc, normfunc, expected) in knownValues
        testdata = map(testfunc, testgrid)
        f1 = @inferred normfunc(testdata, testgrid)
        f2 = @inferred normfunc(testdata, testgrid[2] - testgrid[1])
        @test (abs(f1 - expected) < tolerance)
        @test (abs(f2 - expected) < tolerance)
    end #for
    testOnesMat = ones(testGridLength, testGridLength)
    testZerosMat = zeros(testGridLength, testGridLength)
    @testset "Matrix Norm, func=$func" for func in [WN.norm_l2, WN.norm_l1]
        v1 = @inferred func(testOnesMat, testgrid, testgrid)
        @test isapprox(v1, 1.0)

        v2 = @inferred func(testZerosMat, testgrid, testgrid)
        @test isapprox(v2, 0.0)
    end

end

@testset "Convergence for Grid Input" begin
    ConvergenceNRange = [2^n for n in 3:6]
    L1ErrorValues = zeros(length(ConvergenceNRange))
    L2ErrorValues = zeros(length(ConvergenceNRange))

    @testset "Constant Function C=$C (so norm is domain size), LX=$LX, LT=$LT" for C in -1:1, LX in 3:5, LT in 1:4
        for (i,N) in enumerate(ConvergenceNRange)
            xTestGrid = collect(range(-1, stop=(-1+LX), length=N+1))
            tTestGrid = collect(range(0, stop=(0+LT), length=N))
            testfunc = C*ones(N+1, N)

            testfuncL1NormExact = abs(C)*LX*LT
            L1ErrorValues[i] = abs(testfuncL1NormExact - WN.norm_l1(testfunc, xTestGrid, tTestGrid))

            testfuncL2NormExact = sqrt(abs(C)*LX*LT)
            L2ErrorValues[i] = abs(testfuncL2NormExact - WN.norm_l2(testfunc, xTestGrid, tTestGrid))
        end

        @test all(e <= tolerance for e in L1ErrorValues)
        @test all(e <= tolerance for e in L2ErrorValues)
    end

    fill!(L1ErrorValues, 0)
    fill!(L2ErrorValues, 0)

    # Nonlinear function with changing sign.
    for (i,N) in enumerate(ConvergenceNRange)
        xTestGrid = collect(range(-1, stop=0.5, length=N+1))
        tTestGrid = collect(range(0, stop=1, length=N))
        testfunc1 = [x * (4*t - 1) for x in xTestGrid, t in tTestGrid]

        testfunc1L1NormExact = 25//32
        L1ErrorValues[i] = abs(testfunc1L1NormExact - WN.norm_l1(testfunc1, xTestGrid, tTestGrid))

        testfunc1L2NormExact = 7//8
        L2ErrorValues[i] = abs(testfunc1L2NormExact - WN.norm_l2(testfunc1, xTestGrid, tTestGrid)^2)
    end
    @test all(d <=0 for d in diff(L1ErrorValues))
    @test all(d <=0 for d in diff(L2ErrorValues))

end

@testset "Convergence for Range Input" begin
    ConvergenceNRange = [2^n for n in 3:6]
    L1ErrorValues = zeros(length(ConvergenceNRange))
    L2ErrorValues = zeros(length(ConvergenceNRange))

    @testset "Constant Function C=$C (so norm is domain size), LX=$LX, LT=$LT" for C in -1:1, LX in 3:5, LT in 1:4
        for (i,N) in enumerate(ConvergenceNRange)
            xTestGrid = range(-1, stop=(-1+LX), length=N+1)
            tTestGrid = range(0, stop=(0+LT), length=N)
            testfunc = C*ones(N+1, N)

            testfuncL1NormExact = abs(C)*LX*LT
            L1ErrorValues[i] = abs(testfuncL1NormExact - WN.norm_l1(testfunc, xTestGrid, tTestGrid))

            testfuncL2NormExact = sqrt(abs(C)*LX*LT)
            L2ErrorValues[i] = abs(testfuncL2NormExact - WN.norm_l2(testfunc, xTestGrid, tTestGrid))
        end

        @test all(e <= tolerance for e in L1ErrorValues)
        @test all(e <= tolerance for e in L2ErrorValues)
    end

    fill!(L1ErrorValues, 0)
    fill!(L2ErrorValues, 0)

    # Nonlinear function with changing sign.
    for (i,N) in enumerate(ConvergenceNRange)
        xTestGrid = range(-1, stop=0.5, length=N+1)
        tTestGrid = range(0, stop=1, length=N)
        testfunc1 = [x * (4*t - 1) for x in xTestGrid, t in tTestGrid]

        testfunc1L1NormExact = 25//32
        L1ErrorValues[i] = abs(testfunc1L1NormExact - WN.norm_l1(testfunc1, xTestGrid, tTestGrid))

        testfunc1L2NormExact = 7//8
        L2ErrorValues[i] = abs(testfunc1L2NormExact - WN.norm_l2(testfunc1, xTestGrid, tTestGrid)^2)
    end
    @test all(d <=0 for d in diff(L1ErrorValues))
    @test all(d <=0 for d in diff(L2ErrorValues))
end

@testset "Edge case regression tests" begin
    twoValGrid = [0.0, 1.0]
    @test isnan(WN.norm_w22_only(twoValGrid, [-1.0, 1.0]))
    @test isnan(WN.norm_w22_only(twoValGrid, 1e-1))

    threeValGrid1 = [0.0, 1.0, 2.0]
    @test isapprox(WN.norm_w22_only(threeValGrid1, [1.0, 2.0, 3.0]), 0.0, atol=tolerance)
    @test isapprox(WN.norm_w22_only(threeValGrid1, 1), 0.0, atol=tolerance)

    threeValGrid2 = [1.0, 0.0, 1.0]
    @test WN.norm_w22_only(threeValGrid2, [1.0, 2.0, 3.0]) >= 0
    @test WN.norm_w22_only(threeValGrid2, 1) >= 0
    @test WN.norm_w22_only(-threeValGrid2, [1.0, 2.0, 3.0]) >= 0
    @test WN.norm_w22_only(-threeValGrid2, 1) >= 0

    @test WN.norm_w22_only(-threeValGrid2, 1) == WN.norm_w22_only(threeValGrid2, 1)
    @test WN.norm_w22_only(-threeValGrid2, [1.0, 2.0, 3.0]) == WN.norm_w22_only(threeValGrid2, [1.0, 2.0, 3.0])
end


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

function simplex_transform(img::AbstractArray, nvertices::Integer)
  # binary images are trivial
  nvertices == 1 && return [img]

  ncoords = nvertices - 1

  # simplex construction
  vertices = [eye(ncoords) ones(ncoords)*(1-sqrt(ncoords+1))/2]
  center = sum(vertices, 2) / nvertices
  vertices = vertices .- center

  # map 0 to (0,0,...,0)
  vertices = [zeros(ncoords) vertices]
  idx = map(Int, img + 1)

  result = Array{Array}(ncoords)
  for i=1:ncoords
    coords = similar(img)
    coords[:] = vertices[i,idx[:]]
    result[i] = coords
  end

  result
end


text:
#' ---
#' title: Getting Started with JuMP
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' This tutorial is aimed at providing a quick introduction to writing JuMP code. It assumes familiar with basic optimization and
#' the notion of an [AML](https://en.wikipedia.org/wiki/Algebraic_modeling_language).

#' ## What is JuMP?
#' JuMP ("Julia for Mathematical Programming") is an open-source modeling language that is embedded in Julia. It allows users to
#' users formulate various classes of optimization problems (linear, mixed-integer, quadratic, conic quadratic, semidefinite,
#' and nonlinear) with easy-to-read code. These problems can then be solved using state-of-the-art open-source and commercial solvers.
#' JuMP also makes advanced optimization techniques easily accessible from a high-level language.

#' ## Installing JuMP
#' JuMP is a package for Julia. From Julia, JuMP is installed by using the built-in package manager.
#+ tangle = false

import Pkg
Pkg.add("JuMP")

#' ## A Complete Example
#' Let's try to solve the following linear programming problem by using JuMP and GLPK (a linear and mixed integer programming
#' solver). We will first look at the complete code to solve the problem and then go through it step by step.

#' $$
#' \begin{align*}
#' & \min & 12x + 20y \\
#' & \;\;\text{s.t.} & 6x + 8y \geq 100 \\
#' & & 7x + 12y \geq 120 \\
#' & & x \geq 0 \\
#' & & y \geq 0 \\
#' \end{align*}
#' $$

using JuMP
using GLPK

model = Model(with_optimizer(GLPK.Optimizer))
@variable(model, x >= 0)
@variable(model, y >= 0)
@constraint(model, 6x + 8y >= 100)
@constraint(model, 7x + 12y >= 120)
@objective(model, Min, 12x + 20y)

optimize!(model)

@show value(x);
@show value(y);
@show objective_value(model);

#' ## Step by Step JuMP Code
#' Once JuMP is installed, to use JuMP in your programs, we just need to write-

using JuMP

#' We also need to include a Julia package which provides an appropriate solver. We want to use GLPK.Optimizer here which is
#' provided by the GLPK.jl package.

using GLPK

#' A model object is a container for variables, constraints, solver options, etc. Models are created with the Model() function.
#' The `with_optimizer` syntax is used to specify the optimizer to be used which is GLPK in this case.

model = Model(with_optimizer(GLPK.Optimizer));

#' A variable is modelled using `@variable(name of the model object, variable name and bound, variable type)`. The bound can be a
#' lower bound, an upper bound or both. If no variable type is defined, then it is treated as real.

@variable(model, x >= 0)
@variable(model, y >= 0);

#' A constraint is modelled using `@constraint(name of the model object, constraint)`.

@constraint(model, 6x + 8y >= 100)
@constraint(model, 7x + 12y >= 120);

#' The objective is set in a similar manner using `@objective(name of the model object, Min/Max, function to be optimized)`

@objective(model, Min, 12x + 20y);

#' To solve the optimization problem, we call the optimize function.

optimize!(model)

#' Let's now check the value of objective and variables.

@show value(x);
@show value(y);
@show objective_value(model);


text:
# Test that the way we specify boundary conditions works as expected
using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state

@testset "Boundary condition functions" begin
    ClimateMachine.init()

    FT = Float64

    function init_soil_water!(land, state, aux, localgeo, time)
        myfloat = eltype(state)
        state.soil.water.ϑ_l = myfloat(land.soil.water.initialϑ_l(aux))
        state.soil.water.θ_i = myfloat(land.soil.water.initialθ_i(aux))
    end
    wpf = WaterParamFunctions(FT; Ksat = 1e-7, S_s = 1e-3)
    soil_param_functions = SoilParamFunctions(FT; porosity = 0.75, water = wpf)
    bottom_flux_amplitude = FT(-3.0)
    f = FT(pi * 2.0 / 300.0)
    # nota bene: the flux is -K∇h
    bottom_flux =
        (aux, t) -> bottom_flux_amplitude * sin(f * t) * aux.soil.water.K
    surface_state = (aux, t) -> eltype(aux)(0.2)
    lateral_state = (aux, t) -> eltype(aux)(0.0)
    ϑ_l0 = (aux) -> eltype(aux)(0.2)

    bc = LandDomainBC(
        bottom_bc = LandComponentBC(soil_water = Neumann(bottom_flux)),
        surface_bc = LandComponentBC(soil_water = Dirichlet(surface_state)),
        minx_bc = LandComponentBC(soil_water = Neumann(lateral_state)),
        maxx_bc = LandComponentBC(soil_water = Neumann(lateral_state)),
        miny_bc = LandComponentBC(soil_water = Neumann(lateral_state)),
        maxy_bc = LandComponentBC(soil_water = Neumann(lateral_state)),
    )
    soil_water_model = SoilWaterModel(FT; initialϑ_l = ϑ_l0)
    soil_heat_model = PrescribedTemperatureModel()

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = ()
    m = LandModel(
        param_set,
        m_soil;
        boundary_conditions = bc,
        source = sources,
        init_state_prognostic = init_soil_water!,
    )

    N_poly = 5
    xres = FT(0.2)
    yres = FT(0.2)
    zres = FT(0.01)
    # Specify the domain boundaries.
    zmax = FT(0)
    zmin = FT(-1)
    xmax = FT(1)
    ymax = FT(1)

    driver_config = ClimateMachine.MultiColumnLandModel(
        "LandModel",
        (N_poly, N_poly),
        (xres, yres, zres),
        xmax,
        ymax,
        zmax,
        param_set,
        m;
        zmin = zmin,
    )

    t0 = FT(0)
    timeend = FT(300)
    dt = FT(0.5)

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )
    n_outputs = 30
    every_x_simulation_time = ceil(Int, timeend / n_outputs)

    state_types = (Auxiliary(), GradientFlux())
    dons_arr =
        Dict[dict_of_nodal_states(solver_config, state_types; interp = false)]
    time_data = FT[0] # store time data

    # We specify a function which evaluates `every_x_simulation_time` and returns
    # the state vector, appending the variables we are interested in into
    # `all_data`.

    callback = GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
        dons = dict_of_nodal_states(solver_config, state_types; interp = false)
        push!(dons_arr, dons)
        push!(time_data, gettime(solver_config.solver))
        nothing
    end

    # # Run the integration
    ClimateMachine.invoke!(solver_config; user_callbacks = (callback,))
    computed_bottom_∇h =
        [dons_arr[k]["soil.water.K∇h[3]"][1] for k in 2:n_outputs] ./ [dons_arr[k]["soil.water.K"][1] for k in 2:n_outputs]

    t = time_data[2:n_outputs]
    # we need a -1 out in front here because the flux BC is on -K∇h
    prescribed_bottom_∇h = t -> FT(-1) * FT(-3.0 * sin(pi * 2.0 * t / 300.0))

    MSE = mean((prescribed_bottom_∇h.(t) .- computed_bottom_∇h) .^ 2.0)
    computed_y1_∇h = maximum(abs.(
        [dons_arr[k]["soil.water.K∇h[2]"][1] for k in 2:n_outputs] ./ [dons_arr[k]["soil.water.K"][1] for k in 2:n_outputs],
    ))
    computed_x1_∇h = maximum(abs.(
        [dons_arr[k]["soil.water.K∇h[1]"][1] for k in 2:n_outputs] ./ [dons_arr[k]["soil.water.K"][1] for k in 2:n_outputs],
    ))
    @test MSE < 1e-4
    @test computed_x1_∇h < 1e-10
    @test computed_y1_∇h < 1e-10
end


text:
function estConvRate(residuals)
    #set up and solve least squares problem to estimate slope of convergence line
    n=size(residuals,1)
    A=zeros(n,2)
    A[:,1]=ones(n)
    A[:,2]=collect(0:n-1)

    #solve least squares Ax=residuals
    x=\(A,log.(residuals)/log(10))
    return x[2]
end

function linquad(M,D,K)

    #A=[-D -K; eye(size(M,1)) zeros(size(M))]
    #B=[M zeros(size(M)); zeros(size(M)) eye(size(M,1))]

    A=[-D -K; -K zeros(size(M))]
    B=[M zeros(size(M)); zeros(size(M)) -K]

    #A=[-K zeros(size(M)); zeros(size(M)) M]
    #B=[D M; M zeros(size(M))]

    return (A,B)
end

#get left and right eigenvectors of generalized, nonsymmetric eigenvalue problem
function nseig(A,B)
    if(typeof(A)!=Array{ComplexF64,2})
        A=convert(Array{ComplexF64,2},A)
    end

    if(typeof(B)!=Array{ComplexF64,2})
        B=convert(Array{ComplexF64,2},B)
    end

    alpha,beta,vl,vr=LinearAlgebra.LAPACK.ggev!('V','V',A,B)

    return vl,alpha./beta,vr
end

#trapezoidal quadrature rule
function trapezoidal(nc)
    points=zeros(nc)
    weights=zeros(nc)

    a=-1.0
    b=1.0
    dx=(b-a)/(nc)
    for i in 1:nc
        points[i]=a+dx*(i-1)
        weights[i]=(b-a)/(nc+1)
    end

    return (points,weights)
end

function getRealSpurious(lest,resnorms,emid,ra,rb; inEps=1e-2)
#get distances of estimated eigenvalues from middle of contour:
	#have to rescale the approximate eigenvalues to adjust for the fact that we are using an ellipse and not a circle
	#more complicated contours will require a different approach
	lestdif=abs.(real(lest.-emid)/abs.(ra)+im*imag(lest.-emid)/abs.(rb))
	m0=size(lest,1)
	ninside=0

    p=sortperm(resnorms)

	s=Array{Int64,1}() #indices of eigenvalues inside contour
	spurious=Array{Int64,1}() #indices of spurious eigenvalues inside contour
	other=Array{Int64,1}() #stuff outside contour

	for k in 1:m0
	    if(lestdif[k]<=1.0 && resnorms[k]<inEps)
	        ninside=ninside+1
	        push!(s,k)
	    elseif(lestdif[k]<=1.0 && !(resnorms[k]<inEps))
	        push!(spurious,k)
	    elseif (!(lestdif[k]<=1.0))
	        push!(other,k)
	    end
	end

    return s,spurious,other
end

#find maximum residual inside contour
function getMaxResInside(lest,x,res,emid,ra,rb;inEps=1e-2, getIndices=false, getSpurious=false, set_nin=false)
#get distances of estimated eigenvalues from middle of contour:
	#have to rescale the approximate eigenvalues to adjust for the fact that we are using an ellipse and not a circle
	#more complicated contours will require a different approach
	lestdif=abs.(real(lest.-emid)/abs.(ra)+im*imag(lest.-emid)/abs.(rb))
	(n,m0)=size(x)
	ninside=0

	resnorms=zeros(m0)
    for i in 1:m0
        resnorms[i]=norm(res[:,i])/norm(x[:,i])
    end

    p=sortperm(resnorms)

	s=Array{Int64,1}() #indices of eigenvalues inside contour
	spurious=Array{Int64,1}() #indices of spurious eigenvalues inside contour

	for k in 1:m0
	    if(lestdif[k]<=1.0 && resnorms[k]<inEps)
	        ninside=ninside+1
	        push!(s,k)
	    end
	    if(lestdif[k]<=1.0 && !(resnorms[k]<inEps))
	        push!(spurious,k)
	    end
	end

    maxres=0.0
	if (set_nin!=false)
	    maxres=resnorms[set_nin]
    else
        if(ninside==0)
            maxres=maximum(resnorms)
        else
            for k in 1:m0
                if(lestdif[k]<=1.0 && resnorms[k]<inEps) #using 1.0 as the radius because we rescaled everything when forming lestdif
                   tmp=resnorms[k]#norm(res[:,k])/norm(x[:,k])
                   if(tmp>maxres)
                        maxres=tmp
                   end
                end
            end
        end
    end

    #println("50th residual=$(resnorms[p[50]])")

    if(getIndices & !(getSpurious))
        return maxres,ninside,s
    elseif (getIndices & (getSpurious))
        return maxres,ninside,s,spurious
    else
        return maxres,ninside
    end
end

function getInsideIndex(lest,emid,ra,rb)

    s=Array{Int64,1}()
    m0=size(lest,1)

    lestdif=abs.(real(lest.-emid)/abs.(ra)+im*imag(lest.-emid)/abs.(rb))

    for k in 1:m0
        if(lestdif[k]<=1.0) #using 1.0 as the radius because we rescaled everything when forming lestdif
           push!(s,k)
        end
    end

    return s
end

function getLock(lest,x,res,emid,ra,rb,eps)
#get distances of estimated eigenvalues from middle of contour:
	#have to rescale the approximate eigenvalues to adjust for the fact that we are using an ellipse and not a circle
	#more complicated contours will require a different approach
	lestdif=abs.(real(lest.-emid)/abs.(ra)+im*imag(lest.-emid)/abs.(rb))
	(n,m0)=size(x)
	ninside=0
	for k in 1:m0
	    if(lestdif[k]<=1.0)
	        ninside=ninside+1
	    end
	end

    resnorms=zeros(m0)
    for i in 1:m0
        resnorms[i]=norm(res[:,i])/norm(x[:,i])
    end

    p=sortperm(resnorms)

    locked=Array{Int64,1}()

    maxres=0.0
    for k in 1:m0
        if(lestdif[k]<=1.0 && resnorms[k]<=eps) #using 1.0 as the radius because we rescaled everything when forming lestdif
            push!(locked,k)
        end
    end

    #println("50th residual=$(resnorms[p[50]])")

    return locked
end

function getNotLock(lest,x,res,emid,ra,rb,eps)

    resnorms=zeros(m0)
    for i in 1:m0
        resnorms[i]=norm(res[:,i])/norm(x[:,i])
    end

    notlocked=Array{Int64,1}()

    maxres=0.0
    for k in 1:m0
        if(resnorms[k]>eps) #using 1.0 as the radius because we rescaled everything when forming lestdif
            push!(notlocked,k)
        end
    end

    #println("50th residual=$(resnorms[p[50]])")

    return notlocked
end

#solve reduced quadratic nonlinear problem with linearization
#select the m0 solutions that are closest to the center of the contour
function Quadrrsolve(Q,M,D,K,emid,ra,rb)
        (n,m0)=size(Q)
        #rayleigh ritz
        Mq=Matrix(Q'*(M*Q))
        Dq=Matrix(Q'*(D*Q))
        Kq=Matrix(Q'*(K*Q))

        #solve projected problem
        (Aq,Bq) = linquad(Mq,Dq,Kq)
        lest,xest=eigen(Matrix(Aq),Matrix(Bq))

        #get eigenvectors from linearization
        xall=Q*(xest[m0+1:2*m0,:])

	    #get distances of estimated eigenvalues from middle of contour:
	    #have to rescale the approximate eigenvalues to adjust for the fact that we are using an ellipse and not a circle
	    #more complicated contours will require a different approach
	    lestdif=abs.(real(lest.-emid)/abs.(ra)+im*imag(lest.-emid)/abs.(rb))

        #sort eigenpairs based on those distances:
        p=sortperm(lestdif)

        #take only the eigenpairs closest to the contour center:
        x=xall[:,p[1:m0]]
        lest=lest[p[1:m0]]
        return (lest,x)
end

#Beyn's method for nonlinear problems
#TinvApp = function for applying T^-1 to matrices
#Tapp=T(z)
function beynLoop(Tf,V,ncmax,emid,ra,rb,eps)
    (n,m0)=size(V)

    l=zeros(ComplexF64,m0)
    x=zeros(ComplexF64,n,m0)

    for nc in 1:ncmax
        (l,x)=beyn(Tf,V,nc,emid,ra,rb)
        R=zeros(ComplexF64,n,m0)
        for i in 1:m0
            R[:,i]=Tf(l[i])*x[:,i]
        end
        res,ninside=getMaxResInside(l,x,R,emid,ra,rb)
        println("---Beyn cp $nc  $res")
        if(res<=eps)
            return (l,x)
        end
    end

    return (l,x)
end

function beyn(Tf,V::Array{ComplexF64,2},nc,emid,ra,rb; verbose=false)
    (n,m0)=size(V)

    (gk,wk)=trapezoidal(nc) #contour integral quadrature rule
    offset=pi/nc #offset angle of first quadrature point; makes sure it isn't on real axis for hermitian problems

    A0=zero(V)
    A1=zero(V)
    for k in 1:nc
	#for k in 1:nc
	    if(verbose)
	        println("Beyn system $k of $nc")
	    end
	    theta=gk[k]*pi+pi+offset
	    z=emid+((ra+rb)/2)*exp(im*theta)+((ra-rb)/2)*exp(-1.0*im*theta)

        Qk=zeros(ComplexF64,n,m0)
        Tfz=Tf(z)

        Qk=\(Tfz,V)

	    A0=A0+wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*Qk
	    A1=A1+wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*z*Qk
    end

    (v0,s0,w0) = svd(A0)

    Av=v0'*A1*w0*diagm(0 => 1 ./ s0)

    l,xv=eigen(Av)

    return (l,v0*xv)
end

function ibeyn(Tf,V::Array{ComplexF64,2},alpha,isMaxit,nc,emid,ra,rb; verbose=false)
    (n,m0)=size(V)

    (gk,wk)=trapezoidal(nc) #contour integral quadrature rule
    offset=pi/nc #offset angle of first quadrature point; makes sure it isn't on real axis for hermitian problems

    A0=zeros(V)
    A1=zeros(V)
    for k in 1:nc
	    if(verbose)
	        println("Beyn system $k of $nc")
	    end
	    theta=gk[k]*pi+pi+offset
	    z=emid+((ra+rb)/2)*exp(im*theta)+((ra-rb)/2)*exp(-1.0*im*theta)

        Qk=zeros(ComplexF64,n,m0)
        Tfz=Tf(z)

        #if(typeof(Tfz)==SparseMatrixCSC{Complex{Float64},Int64})
        #    Qk=\(Tfz,V)
        #else
        #    Qk=\(Tfz,V)
        #end

        for i in 1:m0
            #Qk[:,i]=bicgstabl(Tfz,V[:,i],1,max_mv_products=isMaxit,tol=alpha,initial_zero=true,log=false)
            #Qk[:,i]=gmres(Tfz,V[:,i],restart=isMaxit,tol=alpha,initially_zero=true,maxiter=isMaxit,log=false)
            Qk[:,i]=minres(Tfz,V[:,i],tol=alpha,initially_zero=true,maxiter=isMaxit,log=false)
        end

	    A0=A0+wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*Qk
	    A1=A1+wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*z*Qk
    end

    (v0,s0,w0) = svd(A0)

    Av=v0'*A1*w0*diagm(1 ./ s0)

    (l,xv)=eigen(Av)

    return (l,v0*xv)
end

function beynscan(emin,emax,ndiv,Tf,nc,n,m0; verbose=false)

    eigs=Array{ComplexF64,1}()
    resvals=Array{ComplexF64,1}()

    edges=linspace(emin,emax,ndiv+1)
    x0=rand(ComplexF64,n,m0)

    for i in 1:ndiv
        if(verbose)
            println("Div $i of $ndiv")
        end
        emid=(edges[i]+edges[i+1])/2.0
        ra=(edges[i]-edges[i+1])/2.0
        rb=ra

        (lest,xest)=beyn(Tf,x0,nc,emid,ra,rb)

        res=zeros(x0)
        for j in 1:m0
            res[:,j]=Tf(lest[j])*xest[:,j]
        end
        (maxres,ninside,insideIndices)=getMaxResInside(lest,xest,res,emid,ra,rb;getIndices=true)

        for j in 1:ninside
            push!(eigs,lest[insideIndices[j]])
            tmp=norm(res[:,insideIndices[j]])/norm(xest[:,insideIndices[j]])
            push!(resvals,tmp)
            if(verbose)
                println("   $(lest[insideIndices[j]])   $tmp")
            end
        end
    end

    return eigs,resvals
end

function trapezoidalCont(nc)
    offset=pi/nc
    a=offset
    b=2*pi+offset

    dx=(b-a)/nc

    weights=zeros(nc)
    points=zeros(nc)

    weights[1]=dx
    points[1]=a
    for i in 2:nc
        points[i]=points[i-1]+dx
        weights[i]=dx
    end

    return (points,weights)
end

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
using Revise
using PyPlot
using PyCall
using Statistics
using LinearAlgebra
using Infiltrator
using Bem2d

"""
    dislocationinabox()

Comparing half-space and dislocaiton in a box solutions
"""
function dislocationinabox()
    close("all")
    mu = 30e9
    nu = 0.25
    alpha = 1e-7 # scalar preconditioner
    npts = 50
    offset = 1
    xgrid, ygrid = obsgrid(-30e3+offset, -20e3+offset, 30e3-offset, -1-offset, npts)

    # Element geometries and data structures for half space approximation
    elshalfspace = Elements(Int(1e5))
    nfault = 1
    x1, y1, x2, y2 = discretizedline(-10e3, -10e3, 0, 0, nfault) # 45 degree dipping fault
    addelsez!(elshalfspace, x1, y1, x2, y2, "fault")
    nsurf = 100
    x1, y1, x2, y2 = discretizedline(-100e3, 0, 100e3, 0, nsurf) # Free surface
    addelsez!(elshalfspace, x1, y1, x2, y2, "surf")
    idxhalfspace = getidxdict(elshalfspace)

    # Element geometries and data structures for the box cases
    elsbox = Elements(Int(1e5))
    nfault = 1
    x1, y1, x2, y2 = discretizedline(-10e3, -10e3, 0, 0, nfault) # 45 degree dipping fault
    addelsez!(elsbox, x1, y1, x2, y2, "fault")
    nside = 40
    x1, y1, x2, y2 = discretizedline(-30000, -20000, 30000, -20000, nside) # Bottom
    addelsez!(elsbox, x1, y1, x2, y2, "B")
    x1, y1, x2, y2 = discretizedline(30000, -20000, 30000, 0, nside) # Right hand side
    addelsez!(elsbox, x1, y1, x2, y2, "R")
    x1, y1, x2, y2 = discretizedline(30000, 0, -30000, 0, nside) # Top
    addelsez!(elsbox, x1, y1, x2, y2, "T")
    x1, y1, x2, y2 = discretizedline(-30000, 0, -30000, -20000, nside) # Left hand side
    addelsez!(elsbox, x1, y1, x2, y2, "L")
    idxbox = getidxdict(elsbox)

    #
    # Halfspace BEM problem
    #
    _, H_surf_fault = PUTC(slip2dispstress, elshalfspace,
                           idxhalfspace["surf"], idxhalfspace["fault"], mu, nu)
    _, H_surf_surf = PUTC(slip2dispstress, elshalfspace,
                          idxhalfspace["surf"], idxhalfspace["surf"], mu, nu)
    faultslip = sqrt(2) / 2 * ones(2 * nfault)
    Ueff = inv(H_surf_surf) * (H_surf_fault * faultslip)

    # Forward evaluation
    Ufault, Sfault = constdispstress(slip2dispstress, xgrid, ygrid, elshalfspace,
                                     idxhalfspace["fault"], ones(nfault), ones(nfault),
                                     mu, nu)
    Usurf, Ssurf = constdispstress(slip2dispstress, xgrid, ygrid, elshalfspace,
                                   idxhalfspace["surf"], Ueff[1:2:end], Ueff[2:2:end],
                                   mu, nu)
    Utotal = @. Ufault - Usurf
    Stotal = @. Sfault - Ssurf
    plotfields(elshalfspace, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               Utotal, Stotal, "Total")

    # Alternative BEM solution
    T_fault_all, _ = PUTC(slip2dispstress, elshalfspace, idxhalfspace["fault"],
                          collect(1:1:elshalfspace.endidx), mu, nu)
    _, H_surf_all = PUTC(slip2dispstress, elshalfspace, idxhalfspace["surf"],
                         collect(1:1:elshalfspace.endidx), mu, nu)
    bcs = zeros(2 * elshalfspace.endidx)
    bcs[1:2*nfault] .= 0.5
    Ueffalt = inv([T_fault_all; alpha .* H_surf_all]) * bcs

    # Alternative forward evaluation
    Ualt, Salt = constdispstress(slip2dispstress, xgrid, ygrid, elshalfspace,
        collect(1:1:elshalfspace.endidx), Ueffalt[1:2:end], Ueffalt[2:2:end], mu, nu)
    plotfields(elshalfspace, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               Ualt, Salt, "Alternative")

    # Plot difference
    Udiff = @. Ualt - Utotal
    Sdiff = @. Salt - Stotal
    # plotfields(elshalfspace, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
    #            Udiff, Sdiff, "Alternative - Total")

    # Okada comparison if on Linux.
    # I can only PyCall to see okada_wrapper there :(
    if Sys.islinux()
        ow = pyimport("okada_wrapper") # from okada_wrapper import dc3dwrapper
        Uokada = zeros(length(xgrid), 2)
        Sokada = zeros(length(xgrid), 3)
        for i in 1:length(xgrid)
            _, u, gradu = ow.dc3dwrapper((2.0/3.0),
                                         [0.0, xgrid[i]+5e3, ygrid[i]],
                                         5e3,
                                         45,
                                         [-1000000, 1000000],
                                         [-5e3*sqrt(2), 5e3*sqrt(2)],
                                         [0.0, 1.0, 0.0])
            strain = @. 0.5 * (gradu' + gradu)
            stress = mu*LinearAlgebra.I(3)*tr(strain) + 2.0*mu*strain
            Uokada[i, 1] = u[2]
            Uokada[i, 2] = u[3]
            Sokada[i, 1] = stress[2, 2]
            Sokada[i, 2] = stress[3, 3]
            Sokada[i, 3] = stress[2, 3]
        end
        plotfields(elshalfspace, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
                   Uokada, Sokada, "Okada")

        # Let's look at some Residuals
        Utotalres = @. Uokada - Utotal
        Stotalres = @. Sokada - Stotal
        Ualtres = @. Uokada - Ualt
        Saltres = @. Sokada - Salt
        plotfields(elshalfspace, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               Utotalres, Stotalres, "Okada - Total")
        plotfields(elshalfspace, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               Ualtres, Saltres, "Okada - Alternative")
    end

    #
    # Box BEM problem
    #
    T_faultB_faultBRTL, H_pfaultB_faultBRTL = PUTC(slip2dispstress, elsbox,
                                                   [idxbox["fault"]; idxbox["B"]],
                                                   collect(1:1:elsbox.endidx),
                                                   mu, nu)
    T_RTL_faultBRTL, H_RTL_faultBRTL = PUTC(slip2dispstress, elsbox,
                                            [idxbox["R"]; idxbox["T"]; idxbox["L"]],
                                            collect(1:1:elsbox.endidx),
                                            mu, nu)
    bcsbox = zeros(2*elsbox.endidx)
    bcsbox[1:2] .= 0.5
    Ueffbox = inv([T_faultB_faultBRTL ; alpha .* H_RTL_faultBRTL]) * bcsbox
    Ubox, Sbox = constdispstress(slip2dispstress, xgrid, ygrid, elsbox,
        collect(1:1:elsbox.endidx), Ueffbox[1:2:end], Ueffbox[2:2:end], mu, nu)
    plotfields(elsbox, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               Ubox, Sbox, "Box")

end
dislocationinabox()


text:
using CUDA.CUTENSOR
using CUDA
using LinearAlgebra

# using host memory with CUTENSOR doesn't work on Windows
can_pin = !Sys.iswindows()

eltypes = ((Float16, Float16),
            #(Float16, Float32),
            (Float32, Float32),
            #(Float32, Float64),
            (Float64, Float64),
            #(ComplexF16, ComplexF16), (ComplexF16, ComplexF32),
            (ComplexF32, ComplexF32), #(ComplexF32, ComplexF64),
            (ComplexF64, ComplexF64))
@testset for N=2:5
    @testset for (eltyA, eltyC) in eltypes
        # setup
        eltyD = eltyC
        dmax  = 2^div(18,N)
        dims  = rand(2:dmax, N)
        p     = randperm(N)
        indsA = collect(('a':'z')[1:N])
        indsC = indsA[p]
        dimsA = dims
        dimsC = dims[p]
        A     = rand(eltyA, dimsA...)
        dA    = CuArray(A)
        C     = rand(eltyC, dimsC...)
        dC    = CuArray(C)
        can_pin && Mem.pin(A)
        can_pin && Mem.pin(C)

        # simple case
        opA   = CUTENSOR.CUTENSOR_OP_IDENTITY
        opC   = CUTENSOR.CUTENSOR_OP_IDENTITY
        dD    = similar(dC, eltyD)
        opAC  = CUTENSOR.CUTENSOR_OP_ADD
        dD    = CUTENSOR.elementwiseBinary!(1, dA, indsA, opA, 1, dC, indsC, opC,
                                            dD, indsC, opAC)
        D = collect(dD)
        @test D ≈ permutedims(A, p) .+ C
        if can_pin
            Dsimple = similar(C)
             Mem.pin(Dsimple)
            Dsimple = CUDA.@sync CUTENSOR.elementwiseBinary!(1, A, indsA, opA, 1, C, indsC, opC,
                                                             Dsimple, indsC, opAC)
            @test Dsimple ≈ permutedims(A, p) .+ C
        end

        # using integers as indices
        dD = CUTENSOR.elementwiseBinary!(1, dA, 1:N, opA, 1, dC, p, opC, dD, p, opAC)
        D = collect(dD)
        @test D ≈ permutedims(A, p) .+ C
        if can_pin
            Dint = zeros(eltyC, dimsC...)
            Mem.pin(Dint)
            Dint = CUDA.@sync CUTENSOR.elementwiseBinary!(1, A, 1:N, opA, 1, C, p, opC, Dint, p, opAC)
            @test Dint ≈ permutedims(A, p) .+ C
        end

        # multiplication as binary operator
        opAC = CUTENSOR.CUTENSOR_OP_MUL
        dD = CUTENSOR.elementwiseBinary!(1, dA, indsA, opA, 1, dC, indsC, opC,
                                            dD, indsC, opAC)
        D = collect(dD)
        @test D ≈ permutedims(A, p) .* C
        if can_pin
            Dmult = zeros(eltyC, dimsC...)
            Mem.pin(Dmult)
            Dmult = CUDA.@sync CUTENSOR.elementwiseBinary!(1, A, indsA, opA, 1, C, indsC, opC,
                                                           Dmult, indsC, opAC)
            @test Dmult ≈ permutedims(A, p) .* C
        end

        # with non-trivial coefficients and conjugation
        opA = eltyA <: Complex ? CUTENSOR.CUTENSOR_OP_CONJ :
                                CUTENSOR.CUTENSOR_OP_IDENTITY
        opC = CUTENSOR.CUTENSOR_OP_IDENTITY
        opAC = CUTENSOR.CUTENSOR_OP_ADD
        α = rand(eltyD)
        γ = rand(eltyD)
        dD = CUTENSOR.elementwiseBinary!(α, dA, indsA, opA, γ, dC, indsC, opC,
                                            dD, indsC, opAC)
        D = collect(dD)
        @test D ≈ α .* conj.(permutedims(A, p)) .+ γ .* C
        if can_pin
            Dnontrivial = similar(C)
            Mem.pin(Dnontrivial)
            Dnontrivial = CUDA.@sync CUTENSOR.elementwiseBinary!(α, A, indsA, opA, γ, C, indsC, opC,
                                                                 Dnontrivial, indsC, opAC)
            @test Dnontrivial ≈ α .* conj.(permutedims(A, p)) .+ γ .* C
        end

        # test in-place, and more complicated unary and binary operations
        opA = eltyA <: Complex ? CUTENSOR.CUTENSOR_OP_IDENTITY :
                                CUTENSOR.CUTENSOR_OP_SQRT
        # because we use rand, entries of A will be positive when elty is real
        opC = eltyC <: Complex ? CUTENSOR.CUTENSOR_OP_CONJ :
                                CUTENSOR.CUTENSOR_OP_IDENTITY
        opAC = eltyD <: Complex ? CUTENSOR.CUTENSOR_OP_ADD :
                                CUTENSOR.CUTENSOR_OP_MAX
        α = rand(eltyD)
        γ = rand(eltyD)
        dD = CUTENSOR.elementwiseBinary!(α, dA, indsA, opA, γ, dC, indsC, opC,
                                            dC, indsC, opAC)
        D = collect(dC)
        if eltyD <: Complex
            if eltyA <: Complex
                @test D ≈ α .* permutedims(A, p) .+ γ .* conj.(C)
            else
                @test D ≈ α .* sqrt.(convert.(eltyD, permutedims(A, p))) .+
                            γ .* conj.(C)
            end
        else
            @test D ≈ max.(α .* sqrt.(convert.(eltyD, permutedims(A, p))), γ .* C)
        end

        # using CuTensor type
        dA = CuArray(A)
        dC = CuArray(C)
        ctA = CuTensor(dA, indsA)
        ctC = CuTensor(dC, indsC)
        ctD = ctA + ctC
        hD = collect(ctD.data)
        @test hD ≈ permutedims(A, p) .+ C
        ctD = ctA - ctC
        hD = collect(ctD.data)
        @test hD ≈ permutedims(A, p) .- C

        α = rand(eltyD)
        ctD = LinearAlgebra.axpy!(α, ctA, ctC)
        hD = collect(ctD.data)
        @test hD ≈ α.*permutedims(A, p) .+ C

        γ = rand(eltyD)
        ctD = LinearAlgebra.axpby!(α, ctA, γ, ctC)
        hD = collect(ctD.data)
        @test hD ≈ α.*permutedims(A, p) .+ γ.*C
    end
end


text:
import Statistics

function SimulateInvestment(systemSettings::SystemSettings, MO::MarketOperator, allGenCos::Vector{GenCo}, systemLoad::Load,allPowerFlows::Vector{PowerFlow}, technologyPool::Vector{Technology},systemState::SystemState,systemHistory::SystemHistory)
    totalHorizon = systemSettings.totalHorizon
    milestoneYear = systemSettings.milestoneYear
    weightFactors = systemSettings.simulationSettings.representativeData[:Weights]
    year2day = systemSettings.simulationSettings.year2day
    day2hour = systemSettings.simulationSettings.day2hour
    investmentSequence = systemSettings.simulationSettings.investmentSequence

    for year in 1:totalHorizon
        ########################################################################
        #NOTE:Technology decommission and record the amount each agent has decommissioned
        decomissionedCapacity = zeros(1,length(allGenCos))
        for genCoCount in 1:length(allGenCos)
            technologyToDecomission = allGenCos[genCoCount].Decomission(year, allPowerFlows)
            for technology in technologyToDecomission
                decomissionedCapacity[genCoCount] += technology.designProperties.installedCapacity
            end
        end

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

        println("This is year ", year)

        if year % milestoneYear == 0
            push!(systemHistory.systemDecisionSequence, Array{Int64,1}(undef,0))
            ##########################################################################
            #NOTE:start investment process here
            roundCount = 0
            addNew = true
            while addNew == true
                roundCount += 1
                sequence = randperm(length(allGenCos))
                push!(systemHistory.systemDecisionSequence[end],sequence)
                addNew = false
                println("This sequence: ", sequence)
                for i in sequence
                    thisGenCo = allGenCos[i]
                    investmentIndicator,bestTechnology = thisGenCo.InvestmentDecision(allGenCos,MO,technologyPool,systemLoad,year,systemSettings,systemHistory)

                    addNew = addNew || investmentIndicator
                    if investmentIndicator == true
                        #NOTE:build essential connection for new investment decision
                        thisGenCo.ExecuteInvestmentDecision(bestTechnology,allPowerFlows,MO,year)
                        push!(systemHistory.InvestmentHistory,bestTechnology.type)
                        println(bestTechnology.type)
                    end
                    if investmentIndicator != true
                        println(investmentIndicator)
                        push!(systemHistory.InvestmentHistory,:NoNew)
                    end
                end
            end
        end

        ##########################################################################
        #NOTE: system state/environment changes each year
        yearlyTechnologyMix = Array{Technology,1}(undef,0)
        for genCo in allGenCos
            for technology in genCo.technologyPortfolio.technologies
                push!(yearlyTechnologyMix, technology)
            end
        end
        push!(systemHistory.systemCapacityMix, deepcopy(yearlyTechnologyMix))
        ##########################################################################

        ##########################################################################
        #NOTE: real world spot market simulation
        yearlyPriceStream = Array{Float64,1}(undef,0)
        for day in 1:year2day
            dailyPriceStream = Array{Float64,1}(undef,0)
            for hour in 1:day2hour
                systemState.TimeSeries(day,hour,systemLoad,yearlyTechnologyMix)
                thisMarketPrice = MO.ClearMarket(yearlyTechnologyMix, allPowerFlows, systemLoad)
                push!(dailyPriceStream, deepcopy(thisMarketPrice))
                MO.Payment(allPowerFlows,day,hour,:real,systemSettings)
            end
            dailyAveragePrice = Statistics.mean(dailyPriceStream)
            push!(yearlyPriceStream, deepcopy(dailyAveragePrice))
        end

        yearlyAveragePrice = (yearlyPriceStream' * weightFactors[:,2])/sum(weightFactors[:,2])
        append!(systemHistory.systemElectricityPrice , yearlyAveragePrice)
        ##########################################################################

        ##########################################################################
        #NOTE: record history of this year
        for genCo in allGenCos
            push!(genCo.history.investmentHistory, TechnologyPortfolio(genCo.history.yearlyInvestment))
            push!(genCo.history.technologyMixHistory, deepcopy(genCo.technologyPortfolio))
            genCo.history.yearlyInvestment = Array{Technology}(undef,0)
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

global __register__ = Set{String}()

"""
    register_abaqus_keyword(keyword::String)

Add ABAQUS keyword `s` to register. That is, after registration every time
keyword show up in `.inp` file a new section is started
"""
function register_abaqus_keyword(keyword::String)
    push!(__register__, keyword)
    return Type{Val{Symbol(keyword)}}
end

"""
    is_abaqus_keyword_registered(keyword::String)

Return true/false is ABAQUS keyword registered.
"""
function is_abaqus_keyword_registered(keyword::String)
    return keyword in __register__
end