text:
"""
REINFORCE learner (Williams 1992)

Specifically for 2 actions and reward only at the end of episode!

Note
----
The policy parameters θ as well as the feature vector ϕ are
column vectors. Here we use a (states x actions) matrix for simpler access.
also, the feature vector ϕ is a one-hot vector in our case. therefore we
do not need it to keep a vector but access the entries by the
state/action pair [s_t, a_t]
"""
mutable struct REINFORCE
    nr_states::Int
    nr_actions::Int
    α::Float64
    α_eff::Float64
    γ::Float64
    γ_k::Float64
    action_selection_prob_in_episode::Array{Any,1}
    states_in_episode::Array{Any,1}
    actions_in_episode::Array{Any,1}
    θ::Array{Float64, 2}
    target::Array{Float64, 1}
    modulatingfactor_in_episode::Array{Any,1}
    do_debug_print::Bool
end
function REINFORCE(; nr_states = 10, nr_actions = 4, α = 0.1, γ = .9,
    do_debug_print = false)
    θ = zeros(Float64, nr_states, nr_actions)

    REINFORCE(nr_states, nr_actions, α, α, γ, 1., [], [], [], θ, [0.], [], do_debug_print)
end
export REINFORCE
function do_debug_print_final(learner::Union{REINFORCE, REINFORCEModulated, REINFORCEbaseline})
    println("Final:")
    for s in 1:size(learner.θ,1)
        @show learner.θ[s, :]
    end
    println("--------------")
end
""" Update: main function """
function update!(learner::Union{REINFORCE, REINFORCEModulated, REINFORCEbaseline},
                state, action, reward, p_action, modulatingfactor)
    push!(learner.action_selection_prob_in_episode, p_action)
    push!(learner.states_in_episode, state)
    push!(learner.actions_in_episode, action)
    push!(learner.modulatingfactor_in_episode, modulatingfactor)
    if reward>0
         # Backward pass to update the policy
        for k = (length(learner.states_in_episode) : -1 : 1)

            getlearningrate!(learner, learner.modulatingfactor_in_episode[k])

            gettarget!(learner, learner.states_in_episode[k], reward)

            learner.θ[learner.states_in_episode[k], learner.actions_in_episode[k]] +=
                        learner.α_eff * learner.target[1] * (1. - learner.action_selection_prob_in_episode[k])
            learner.θ[learner.states_in_episode[k], 3-learner.actions_in_episode[k]] +=
                        learner.α_eff * learner.target[1] * (0. - (1. - learner.action_selection_prob_in_episode[k]))
            learner.γ_k *= learner.γ
        end
        # Reset
        learner.γ_k = 1.
        learner.action_selection_prob_in_episode = []
        learner.states_in_episode = []
        learner.actions_in_episode = []
        learner.modulatingfactor_in_episode = []
    end
    if learner.do_debug_print
        do_debug_print_final(learner)
    end
end
export update!
function gettarget!(learner::Union{REINFORCE, REINFORCEModulated}, state, reward)
    learner.target[1] = reward * learner.γ_k  #  note: assuming a single reward, set during the forward pass above
end


text:
# MLJ Demo from the June 2019 MLJ/sktime Sprint

# ### Load environment and seed RNG

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using Random
Random.seed!(1234);

# ## 1. Basics

# ### Constructing a task

using MLJ
models()

#-

using RDatasets
boston = dataset("MASS", "Boston");
first(boston)

#-

scitypes(boston)

#-

task = supervised(data=boston,
                  target=:MedV,
                  ignore=:Chas,
                  types=Dict(:Rad=>Continuous,:Tax=>Continuous),
                  is_probabilistic=true)

#-

shuffle!(task)

#-

models(task)

#-

task.is_probabilistic=false
models(task)

# ### Binding a task to a model

@load RidgeRegressor
model = RidgeRegressor()

#-

model.lambda = 0.1

#-

mach = machine(model, task)

# ### Evaluation

evaluate!(mach, resampling=Holdout(fraction_train=0.7),
          measure=rms)

# ### Fitting

train, test = partition(1:nrows(task), 0.7)

#-

fit!(mach, rows=train)

#-

fitted_params(mach)
report(mach)

# ### Predicting on data in the task:

yhat = predict(mach, rows=test);
yhat = predict(mach, task)[test];
yhat = predict(mach, task[test]);
yhat[1:3]
#-

y = task[test].y;

#-

rms(yhat, y)

# ### Predicting on "new" data

X = task[test].X;
yhat = predict(mach, X);
yhat[1:4]

# ## 2. Random Forests - Case Study in Ensembling and Tuning

# ### Quick look at single tree

task = load_iris() # built-in task
shuffle!(task)
@load DecisionTreeClassifier
tree = DecisionTreeClassifier()
mach = machine(tree, task)
evaluate!(mach)

#-

yhat = predict(mach, rows=1:3)

#-

[pdf(d, "virginica") for d in yhat]

#-

predict_mode(mach, rows=1:3)

# ### A random forest

forest = EnsembleModel(atom=tree)

#-

mach = machine(forest, task)

r = range(forest, :n, lower=10, upper=1000, scale=:log10)

iterator(r,5)

#-

curves = learning_curve!(mach,
                         resampling=Holdout(fraction_train=0.8),
                         nested_range=r,
                         measure=cross_entropy, n=4,
                         verbosity=0)
#-

using Plots
plot(curves.parameter_values, curves.measurements,
     xlab="number of trees", ylab="rms ")

#-

forest.n = 100; # 750 better

#-

params(forest) # all hyperparameters, as a named tuple

#-

r1 = range(tree, :n_subfeatures, lower=1, upper=4);
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0);
nested_ranges = (atom=(n_subfeatures=r1,),
                 bagging_fraction=r2)

#-

tuned_forest = TunedModel(model=forest,
                          tuning=Grid(resolution=12),
                          resampling=CV(nfolds=6),
                          nested_ranges=nested_ranges,
                          measure=cross_entropy)
#-

params(tuned_forest)

#-

mach = machine(tuned_forest, task)

#-

evaluate!(mach, resampling=Holdout(fraction_train=0.7),
          measure=[cross_entropy, misclassification_rate], verbosity=2)

#-

fitted_params(mach)

#-

best = fitted_params(mach).best_model
@show best.bagging_fraction best.atom.n_subfeatures;

#-

plot(mach)

# ## 3. Learning Networks

task = load_reduced_ames()

#-

X, y = task.X, task.y;
scitypes(X)

#-

models(task)

# ## Workflow without a pipeline (static data)

hot = OneHotEncoder()
hotM = machine(hot, X)
fit!(hotM)
Xc = transform(hotM, X)

knn = KNNRegressor()
knnM = machine(knn, Xc, y)
fit!(knnM)
yhat= predict(knnM, Xc)

#-

# ## Pipelining (dynamic data)

X = source(X)
y = source(y)

# Identical code to build network (`fit!`'s can be dropped)

hot = OneHotEncoder()
hotM = machine(hot, X)
Xc = transform(hotM, X)

knn = KNNRegressor()
knnM = machine(knn, Xc, y)
yhat = predict(knnM, Xc)

# Fit the network in one go:

fit!(yhat)

#-

hot.drop_last=true

#-

fit!(yhat)

#-

knn.K = 7

#-

fit!(yhat)

# Instead of `yhat[1:10]` we have

yhat(rows=1:10)

# Or call on new data:

yhat(task.X)

# The new data is "plugged into" the orgin node, which must be unique:

origins(yhat) == [X,]

# ## Exporting pipeline as stand-alone model

comp = @from_network SmartKNNRegressor(one_hot_encoding=hot, knn_regressor=knn) <= (X, y, yhat)

#-

params(comp)
mach = machine(comp, task)
evaluate!(mach, measure=rms)

## # More complicated example

using CategoricalArrays
x1 = rand(6);
x2 = categorical([mod(rand(Int),2) for i in 1:6]);
x3 = rand(6);
y = exp.(x1 -2x3 + 0.1*rand(6))
X = (x1=x1, x2=x2, x3=x3)

# Here's a learning network that: (i) One-hot encodes the input table `X`; (ii)
# Log transforms the continuous target `y`; (iii) Fits specified
# K-nearest neighbour and ridge regressor models to the data; (iv)
# Computes an average of individual model predictions; and (v) Inverse
# transforms (exponentiates) the blended predictions.

X = source(X)
y = source(y)

hot = machine(OneHotEncoder(), X)

# `W`, `z`, `zhat` and `yhat` are nodes in the network:

W = transform(hot, X) # one-hot encode the input
z = log(y) # transform the target

ridge = RidgeRegressor(lambda=0.1)
knn = KNNRegressor()

ridgeM = machine(ridge, W, z)
knnM = machine(knn, W, z)

# Average the predictions of the KNN and ridge models:

zhat = 0.5*predict(ridgeM, W) + 0.5*predict(knnM, W)

# Inverse the target transformation:

yhat = exp(zhat)

# A tree "splat" of the learning network terminating at `yhat`:

MLJ.tree(yhat)

#-

blend = @from_network Blen(ridge=ridge, knn=knn) <= (X, y, yhat)
blend.ridge.lambda = 0.2
mach = machine(blend, load_reduced_ames())
evaluate!(mach, measure=rms)

ext: jl


text:
# GeometryBasics' HyperRectangle seems cumbersome to index for specific points.
# And I also need to handle points from VoronoiDelaunay
struct Rectangle{T}
    LowerLeft::T
    UpperRight::T

    function Rectangle(p1::T, p2::T) where T
        left, right = minmax(getx(p1), getx(p2))
        lower, upper = minmax(gety(p1), gety(p2))

        if left >= right || lower >= upper
            throw(ArgumentError("Empty rectangle"))
        end

        lower_left = T(left, lower)
        upper_right = T(right, upper)

        new{T}(lower_left, upper_right)
    end
end

function Rectangle(p1::GeometryBasics.Point2{T}, p2::GeometryBasics.Point2{T}) where T <: Integer
    Rectangle(float.(p1), float.(p2))
end

Base.eltype(::Rectangle{T}) where T = T

left(rect::Rectangle) = getx(rect.LowerLeft)
right(rect::Rectangle) = getx(rect.UpperRight)
lower(rect::Rectangle) = gety(rect.LowerLeft)
upper(rect::Rectangle) = gety(rect.UpperRight)

upper_right(rect::Rectangle{T}) where T = rect.UpperRight
upper_left(rect::Rectangle{T}) where T = T(left(rect), upper(rect))
lower_right(rect::Rectangle{T}) where T = T(right(rect), lower(rect))
lower_left(rect::Rectangle{T}) where T = rect.LowerLeft

function area(rect::Rectangle)
    (right(rect) - left(rect)) * (upper(rect) - lower(rect))
end

"""
    corner_nearest_neighbor(points, rect)

For each corner in the rectangle `rect`, find the point(s) in `points` that are closest.
The result is a `Dict` where keys are `rect`'s corners and the values are the nearest neighbors.
"""
function corner_nearest_neighbor(points, rect)
    rect_corners = corners(rect)
    neighbors = Dict(1:4 .=> [Vector{Int64}(undef, 0)])
    corner_distances = [Inf for _ in 1:4]

    for (index, point) in enumerate(points)
        for corner_index in 1:4
            corner_dist = abs2(point, rect_corners[corner_index])

            if corner_dist ≈ corner_distances[corner_index]
                push!(neighbors[corner_index], index)
            elseif corner_dist < corner_distances[corner_index]
                corner_distances[corner_index] = corner_dist
                neighbors[corner_index] = [index]
            end
        end

    end

    return neighbors
end

function corners(rect)
    [upper_right(rect), upper_left(rect), lower_left(rect), lower_right(rect)]
end

@inline function isinside(point, rect)
    isinside_x = left(rect) <= getx(point) <= right(rect)
    isinside_y = lower(rect) <= gety(point) <= upper(rect)

    isinside_x && isinside_y
end

function map_to_computation_rectangle(points::Vector{T}, from::Rectangle{T}, to::Rectangle) where T
    offsetx_from = left(from)
    offsety_from = lower(from)

    offsetx_to = left(to)
    offsety_to = lower(to)

    slopex = (right(to) - left(to)) / (right(from) - left(from))
    slopey = (upper(to) - lower(to)) / (upper(from) - lower(from))

    no_points = length(points)
    transformed_points = Vector{IndexablePoint2D}(undef, no_points)
    for (index, point) in enumerate(points)
        if !isinside(point, from)
            throw(error("Point is not inside rectangle"))
        end

        transformed_points[index] = IndexablePoint2D(
            offsetx_to + (getx(point) - offsetx_from) * slopex,
            offsety_to + (gety(point) - offsety_from) * slopey,
            index
        )
    end

    return transformed_points
end

function map_rectangle(points::Vector{T}, from::Rectangle{T}, to::Rectangle{S}) where T where S
    offsetx_from = left(from)
    offsety_from = lower(from)

    offsetx_to = left(to)
    offsety_to = lower(to)

    slopex = (right(to) - left(to)) / (right(from) - left(from))
    slopey = (upper(to) - lower(to)) / (upper(from) - lower(from))

    no_points = length(points)
    transformed_points = Vector{S}(undef, no_points)
    for (index, point) in enumerate(points)
        if !isinside(point, from)
            throw(error("Point is not inside rectangle"))
        end

        transformed_points[index] = S(
            offsetx_to + (getx(point) - offsetx_from) * slopex,
            offsety_to + (gety(point) - offsety_from) * slopey
        )
    end

    return transformed_points
end


text:
using Plots, Dates, Measures, Interpolations, CSV, DataFrames, LaTeXStrings, SugarKelp

# Loading data scraped from their paper
offset = 4 * 31 + 28 + 30 * 2.0
# collected from fig2 of origional paper with https://www.graphreader.com/
temp_file = CSV.read("examples/bs2012/temp.csv", DataFrame);sort!(temp_file);temp_t = temp_file.day;temp = temp_file.temp
no3_file = CSV.read("examples/bs2012/no3.csv", DataFrame);sort!(no3_file);no3_t = no3_file.day;no3 = no3_file.no3
irr_file = CSV.read("examples/bs2012/irr.csv", DataFrame);sort!(irr_file);irr_t = irr_file.day;irr = irr_file.irr
temp_t .+= offset ;no3_t .+= offset ;irr_t .+= offset

# irr .*= exp(-0.07*13) # Background attenuation is set to 0.07m^-1(?)

# Useful part of the example starts here
t_i = offset + 12
nd = 365 + 16
lat = 60.257
u = 0.15

u_arr = Interpolations.LinearInterpolation([t_i:t_i + nd;], fill(u, Int(nd + 1)))

temp_arr = Interpolations.LinearInterpolation(temp_t, temp, extrapolation_bc=Flat())
no3_arr = Interpolations.LinearInterpolation(no3_t, no3, extrapolation_bc=Flat())
irr_arr = Interpolations.LinearInterpolation(irr_t, irr, extrapolation_bc=Flat())

a_0 = 30;n_0 = 0.01;c_0 = 0.6

solution, results = SugarKelp.solve(t_i, nd, u_arr, temp_arr, irr_arr, no3_arr, lat, a_0, n_0, c_0, "../src/parameters/2013.jl", 1)

# Horrible plotting and loading/inline horriblness from their paper from here down
pyplot()
l = @layout [[a{0.25h};b{0.25h};c{0.5h}] grid(2, 1)]
plot(layout=l)

plot!(temp_t,temp,sp=1,ylabel=L"Temperature/$^o$C",xlabel="Month",legend=false,label="Temperature")
plt = twinx()
plot!([no3_t[1],no3_t[1]],[no3[1]no3[1]],sp=6,label="Temperature")
plot!(no3_t,no3,sp=6,ylabel=L"Nitrate/$\mu$ mol / L",legend=:right,label="Nitrate")

plot!(irr_t,irr,sp=2,ylabel=L"Irradiance/mol photons/m$^2$/day",legend=false,xlabel="Month")

# paper data
c_t = [10.016,100.156,130.973,164.872,192.607,224.965,266.568,291.222,343.611,409.097] .+ offset
c = [0.312,0.349,0.292,0.195,0.224,0.266,0.281,0.307,0.363,0.373]
n_t = [10.255,96.239,135.681,160.924,195.633,222.454,269.785,288.717,345.514,406.255] .+ offset
n = [0.007,0.014,0.021,0.024,0.025,0.022,0.023,0.017,0.012,0.013]
area_file = CSV.read("examples/bs2012/area.csv", DataFrame);sort!(area_file);area_t = area_file.day .+ offset;area = area_file.area
c2_file = CSV.read("examples/bs2012/carbon.csv", DataFrame);sort!(c2_file);c2_t = c2_file.day .+ offset;c2 = c2_file.carbon
n2_t = [14.8696,32.8696,57.913,75.1304,90.7826,104.8696,125.2174,145.5652,156.5217,162.7826,170.6087,178.4348,184.6957,194.087,200.3478,214.4348,237.913,248.8696,259.8261,267.6522,270.7826,273.913,278.6087,288,297.3913,308.3478,316.1739,327.1304,336.5217,352.1739,364.6957,378.7826,396.7826] .+ offset
n2 = [0.0088,0.0095,0.0109,0.0124,0.0146,0.017,0.0209,0.0249,0.027,0.0284,0.0291,0.0285,0.0281,0.0265,0.025,0.0229,0.0189,0.0175,0.0165,0.0154,0.0145,0.0136,0.0127,0.0112,0.0102,0.0095,0.0091,0.0086,0.0083,0.0083,0.0084,0.0086,0.0088]
# structural to dry weight conversion (paper plots g/g dry weight where as g/g structural weight is used in calculations), I think they are also plotting the total carbon not just the reserve
n_factor = (results.nitrogen .- N_min) .* K_N
c_factor = (results.carbon .- C_min) .* K_C
w_factor = 1 .+ n_factor .+ c_factor .+ C_min .+ N_min

plot!(results.time,results.area,sp=3,xlabel="Month", ylabel=L"Frond Area/dm$^2$",ylim=(29, 46),label="Model",legend=false)
plot!(area_t,area,label="B&S, 2012",sp=3)

plot!(results.time,(results.nitrogen .+ N_struct) ./ w_factor,sp=4, xlabel="Month", ylabel="Nitrogen reserve/gN/g dw",label="Model")
plot!(n2_t,n2,label="B&S, 2012",sp=4)
plot!(n_t,n,sp=4,seriestype=:scatter, label="Sjotun, 1993",color="#ebcb8b")

plot!(results.time,(results.carbon .+ C_struct) ./ w_factor,sp=5, xlabel="Month", ylabel="Carbon reserve/gC/g dw",label="Model",legend=false)
plot!(c2_t,c2,label="B&S, 2012",sp=5)
plot!(c_t,c,sp=5,seriestype=:scatter, label="Sjotun, 1993",color="#ebcb8b")

t_ticks = []
val_ticks = []
for day in t_i:t_i + nd
    date = Date(1981, 1, 1) + Dates.Day(day)
    if Dates.format(date, "d") == "1"
        push!(t_ticks, day)
        push!(val_ticks, Dates.format(date, "U")[1])
    end
end

plot!(xticks=(t_ticks, val_ticks))
display(display(plot!()))


text:
@testset "Encoder" begin

  o = Jib.Requests.Encoder.Enc()

  o(nothing, true, 1, 1.1, Jib.PRICE, :a, "test", (a=1, b="2"))

  m = split(String(take!(o.buf)), '\0')

  @test length(m) == 9
  @test m[1] == ""            # nothing
  @test m[2] == "1"           # Bool
  @test m[3] == "1"           # Int
  @test m[4] == "1.1"         # Float64
  @test m[5] == "1"           # Enum{Int32}
  @test m[6] == "a"           # Symbol
  @test m[7] == "test"        # String
  @test m[8] == "a=1;b=2;"    # NamedTuple
  @test m[9] == ""

  # Condition
  o(Jib.ConditionTime("o", true, "yyyymmdd"))
  @test String(take!(o.buf)) == "3\0o\x001\0yyyymmdd\0"

  # splat
  c = Jib.ComboLeg(conId=1, action="action")
  o(Jib.Requests.Encoder.splat(c, [1,3]),
    Jib.Requests.Encoder.splat(c))

  m = split(String(take!(o.buf)), '\0')
  @test m == ["1", "action", "1", "0", "action", "", "0", "0", "", "-1", ""]

  # Unsuported types
  @test_throws ErrorException o(Int32(2))
  @test_throws ErrorException o(Float32(1.))

end


text:
@testset "constructors" begin
    @testset "Mag" begin
        # Results for doubles are only upper bounds, even for integers
        # they are quite crude
        @test Mag() == Mag(UInt64(0)) == zero(Mag) == zero(Mag()) == Mag(0.0)
        @test Mag(UInt64(1)) == one(Mag) == one(Mag()) <= Mag(1.0)
        @test π < Float64(Mag(π)) < 3.15
    end

    @testset "Arf" begin
        for T in [UInt, Int, Float64, Mag, Arf, BigFloat]
            @test Arf(zero(T)) == zero(Arf)
            @test Arf(one(T)) == one(Arf)
            @test precision(Arf(zero(T), prec = 80)) == 80
        end

        @test Arf(zero(Arf).arf) == zero(Arf)
        @test precision(Arf(zero(Arf).arf, prec = 80)) == 80

        @test precision(Arf(Arf(prec = 80))) == 80
        @test precision(Arf(BigFloat(0, precision = 80))) == 80

        @test precision(zero(Arf(prec = 80))) == 80
        @test precision(one(Arf(prec = 80))) == 80
    end

    @testset "Arb" begin
        for T in [UInt, Int, Float64, Arf, Arb, BigFloat, Rational{Int}]
            @test Arb(zero(T)) == zero(Arb)
            @test Arb(one(T)) == one(Arb)
            @test precision(Arb(zero(T), prec = 80)) == 80
        end

        @test Arb(zero(Arb).arb) == Arb("0.0") == zero(Arb)
        @test Arb(one(Arb).arb) == Arb("1.0") == one(Arb)
        @test precision(Arb(zero(Arb).arb, prec = 80)) ==
              precision(Arb("0.0", prec = 80)) ==
              80

        @test precision(Arb(Arf(prec = 80))) == 80
        @test precision(Arb(Arb(prec = 80))) == 80
        @test precision(Arb(BigFloat(0, precision = 80))) == 80

        @test precision(zero(Arb(prec = 80))) == 80
        @test precision(one(Arb(prec = 80))) == 80

        @test Arb(3) < Arb(π) < Arb(4)
        @test Arb(2) < Arb(ℯ) < Arb(3)
        @test Arb(0) < Arb(MathConstants.γ) < Arb(1)
        @test Arb(0) < Arb(MathConstants.catalan) < Arb(1)
        @test Arblib.overlaps(Arb(BigFloat(MathConstants.φ)), Arb(MathConstants.φ))
        @test precision(Arb(π, prec = 80)) == 80
        @test precision(Arb(ℯ, prec = 80)) == 80
        @test precision(Arb(MathConstants.γ, prec = 80)) == 80
        @test precision(Arb(MathConstants.catalan, prec = 80)) == 80
        @test precision(Arb(MathConstants.φ, prec = 80)) == 80
    end

    @testset "Acb" begin
        for T in [UInt, Int, Float64, Arf, Arb, BigFloat, Rational{Int}]
            @test Acb(zero(T)) == zero(Acb)
            @test Acb(one(T)) == one(Acb)
            @test precision(Acb(zero(T), prec = 80)) == 80

            @test imag(Acb(zero(T), one(T))) == one(Arb)
            @test precision(Acb(zero(T), zero(T), prec = 80)) == 80

            @test imag(Acb(Complex(zero(T), one(T)))) == one(Arb)
            @test precision(Acb(Complex(zero(T), zero(T)), prec = 80)) == 80
        end

        @test Acb(zero(Acb).acb) == Acb("0.0") == zero(Acb)
        @test Acb(one(Acb).acb) == Acb("1.0") == one(Acb)
        @test precision(Acb(zero(Acb).acb, prec = 80)) ==
              precision(Acb("0.0", prec = 80)) ==
              80

        @test imag(Acb(zero(Arb).arb, one(Arb).arb)) == imag(Acb("0.0", "1.0")) == one(Arb)
        @test precision(Acb(zero(Arb).arb, zero(Arb).arb, prec = 80)) ==
              precision(Acb("0.0", "0.0", prec = 80)) ==
              80

        @test precision(Acb(Arf(prec = 80))) == 80
        @test precision(Acb(Arb(prec = 80))) == 80
        @test precision(Acb(Acb(prec = 80))) == 80
        @test precision(Acb(BigFloat(0, precision = 80))) == 80

        @test precision(Acb(Arf(prec = 80), Arf(prec = 100))) == 100
        @test precision(Acb(Arb(prec = 80), Arb(prec = 100))) == 100
        @test precision(Acb(BigFloat(0, precision = 80), BigFloat(0, precision = 100))) ==
              100

        @test precision(Acb(Complex(Arf(prec = 80), Arf(prec = 100)))) == 100
        @test precision(Acb(Complex(Arb(prec = 80), Arb(prec = 100)))) == 100
        @test precision(Acb(Complex(
            BigFloat(0, precision = 80),
            BigFloat(0, precision = 100),
        ))) == 100

        @test precision(zero(Arb(prec = 80))) == 80
        @test precision(one(Arb(prec = 80))) == 80

        @test isequal(real(Acb(π)), Arb(π))
        @test isequal(real(Acb(ℯ)), Arb(ℯ))
        @test isequal(real(Acb(MathConstants.γ)), Arb(MathConstants.γ))
        @test precision(Acb(π, prec = 80)) == 80
        @test precision(Acb(ℯ, prec = 80)) == 80
        @test precision(Acb(MathConstants.γ, prec = 80)) == 80

        @test isequal(Acb(π, -1), Acb(Arb(π), Arb(-1)))
        @test isequal(Acb(BigFloat(1.3), ℯ), Acb(Arb(1.3), Arb(ℯ)))
    end

    @testset "Others" begin
        # TODO: Move these to other file
        @test Float64(Arf(2.0)) isa Float64
        @test Float64(Arf(2.0)) == 2.0
        @test convert(Float64, Arf(2.0)) isa Float64
        @test convert(Float64, Arf(2.0)) == 2.0
        @test Int(Arf(2.0)) isa Int
        @test Int(Arf(2.0)) == 2
        @test convert(Int, Arf(2.0)) isa Int
        @test convert(Int, Arf(2.0)) == 2
        @test Int(Arf(2.0)) == 2
        @test Arblib.get_si(Arf(0.5)) == 0
        @test Arblib.get_si(Arf(0.5); rnd = Arblib.ArbRoundFromZero) == 1

        @test Int(Arf(2)) == 2
        @test_throws InexactError Int(Arf(2.5))
        @test Float64(Mag(2.5)) ≥ 2.5
        @test Float64(Arf(2.5)) == 2.5
        @test Float64(Arb(2.5)) == 2.5
        @test ComplexF64(Acb(2.25)) == 2.25 + 0im
        @test BigFloat(Arf(2.5)) == 2.5
        @test BigFloat(Arb(2.5)) == 2.5
        @test precision(BigFloat(Arf(2.5; prec = 96))) == 96
    end

    @testset "zeros/ones" begin
        for T in [Arf, Arb, Acb]
            @test zeros(T, 2) == [zero(T), zero(T)]
            @test ones(T, 2) == [one(T), one(T)]
        end

        v = zeros(Arb, 2)
        @test length(v) == 2
        @test v[1] == v[2] == zero(Arb)
        Arblib.set!(v[2], 1)
        @test v[1] == zero(Arb) && v[2] == one(Arb)

        v = ones(Arb, 2)
        @test length(v) == 2
        @test v[1] == v[2] == one(Arb)
        Arblib.set!(v[2], 0)
        @test v[1] == one(Arb) && v[2] == zero(Arb)
    end
end


text:
using GitHub, JSON, HTTP, MbedTLS
using Compat, Dates, Test
using GitHub: Branch, name

using Base64

include("ghtype_tests.jl")
include("event_tests.jl")
include("read_only_api_tests.jl")
include("auth_tests.jl")


text:
"""
    SentenceBatch

    Keeps input/outputs of a sentence with encoded fields.
"""
struct SentenceBatch
    encodedIOs::Vector{EncodedIO}
    seqInputs::Vector{Int}
    seqSizes::Vector{Int}
    seqOutputs::Matrix{Int}
    masks::Matrix{Int}
    tagRange::Vector{UnitRange{Int}}
    unsortingIndices::Vector{Int}
    wordInputs::Vector{Int} # words
    compositeOutputs::Vector{Int}
end

"""
    miniBatch(data, v::Vocabulary; sorted=false)
    miniBatch(sentence::Vector{EncodedAnalysis}, v::Vocabulary)

    Creates SentenceBatche given encodedIO array for a sentence
"""
miniBatch(data, v::Vocabulary; sorted=false) =
map(d->miniBatch(d,v), sorted ? sort(data; by=length) : data)

function miniBatch(sentence::Vector{EncodedIO}, v::Vocabulary)
    #complete inputs
    wordInputs = map(a->a.completeWord,sentence)
    #sequence inputs
    seqInputs, seqSizes, unsortingIndices = sequenceInputBatch(sentence)
    #get correct analyses
    analyses = map(a->first(a.analyses),sentence)
    #composite output
    compositeOutputs = map(a->a.compositeTag,analyses)
    #sequence outputs
    seqOutputs, masks, tagRange = sequenceOutputBatch(analyses, v)
    SentenceBatch(sentence, seqInputs, seqSizes, seqOutputs, masks, tagRange,
                            unsortingIndices, wordInputs, compositeOutputs)
end

function sequenceInputBatch(sentence::Vector{EncodedIO})
    #sequence input
    words = map(a->a.chars,sentence)
    sortingIndices = sortperm(length.(words); rev=true)
    unsortingIndices = sortperm(sortingIndices)
    seqInputs, seqSizes = _pack_sequence(words[sortingIndices])
    return seqInputs, seqSizes, unsortingIndices
end

function sequenceOutputBatch(analyses::Vector{EncodedAnalysis}, v::Vocabulary)
    #sequence output
    lengths  = map(a->length(a.lemma)+length(a.tags)+1,analyses)
    tagRange = map(a->length(a.lemma)+1:length(a.lemma)+length(a.tags), analyses)
    maxoutlen = maximum(lengths)
    seqOutputs = fill(v.specialIndices.mask, length(analyses), maxoutlen)
    masks = ones(Int, length(analyses), maxoutlen)
    @inbounds for (k,a) in enumerate(analyses)
        seqOutputs[k,1:length(a.lemma)] .= a.lemma
        seqOutputs[k,tagRange[k]] .= a.tags
        seqOutputs[k,lengths[k]] = v.specialIndices.eow
        if a.isValid
            masks[k,lengths[k]+1:maxoutlen] .= 0
        else
            masks[k,:] .= 0
        end
    end
    return seqOutputs, masks, tagRange
end


text:
using PowerNumbers, Test, IntervalArithmetic, HypergeometricFunctions
import PowerNumbers: PowerNumber, LogNumber, realpart, apart, bpart, alpha, beta

@testset "RiemannDual -> PowerNumber" begin
    for h in (0.1,0.01), a in (2exp(0.1im),1.1)
        @test log(PowerNumber(0,a,0,1))(h) ≈ log(h*a)
        @test log(PowerNumber(a,Inf,-1,0))(h) ≈ log(a/h)
    end

    for h in (0.1,0.01), a in (2exp(0.1im),1.1)
        @test log1p(PowerNumber(-1,a,0,1))(h) ≈ log(h*a)
        @test log1p(PowerNumber(a,Inf,-1,0))(h) ≈ log(a/h)
    end

    h=0.0001
    for z in (PowerNumber(1,3exp(0.2im),0,1), PowerNumber(1,0.5exp(-1.3im),0,1),
              PowerNumber(-1,3exp(0.2im),0,1), PowerNumber(-1,0.5exp(-1.3im),0,1),
              PowerNumber(-1,1,0,1), PowerNumber(1,-1,0,1))
        @test atanh(z)(h) ≈  atanh(apart(z)*h^alpha(z)+bpart(z)*h^beta(z)) atol = 1E-4
    end

    #z = PowerNumber(1,-0.25,1)
    h = 0.0000001
    #@test speciallog(z)(h) ≈ speciallog(realpart(z)+epsilon(z)*h^alpha(z)) atol=1E-4

    z = PowerNumber(1,2,0,1)
    @test 1/(1-z(h)) ≈ (1/(1-z))(h) rtol=0.0001

    @test real(LogNumber(2im,im+1)) == LogNumber(0,1)
    @test imag(LogNumber(2im,im+1)) == LogNumber(2,1)
    @test conj(LogNumber(2im,im+1)) == LogNumber(-2im,1-im)

end

@testset "PowerNumbers arithmetic" begin
    h = 0.000000001
    @test inv(PowerNumber(1.0,2.0,-1/2,0))(h) ≈ inv(h^(-1/2) + 2) rtol = 0.001
    @test inv(PowerNumber(0.0,1.0,0,1/2))(h) ≈ inv(h^(1/2)) rtol = 0.001
    @test inv(PowerNumber(2.0,1.0,0,1/2))(h) ≈ inv(h^(1/2)+2) rtol = 0.001

    @test real(PowerNumber(2im,im+1,0,0.5)) == PowerNumber(0,1,0,0.5)
    @test imag(PowerNumber(2im,im+1,0,0.5)) == PowerNumber(2,1,0,0.5)
    @test conj(PowerNumber(2im,im+1,0,0.5)) == PowerNumber(-2im,1-im,0,0.5)

    @test (3*PowerNumber(0.01+3im,0.2im+1,0,0.4) - 4*PowerNumber(0.6-1.5im,1.5-0.3im,0,0.4))/3 == PowerNumber(-0.79+5.0im,-1.0 + 0.6im,0,0.4)
end

@testset "sin" begin
    ε = PowerNumber(1,1)
    @test sin(sqrt(ε))^2 === PowerNumber(1.0,1.0)
    @test sin(sqrt(ε))^2.0 === PowerNumber(1.0,1.0)
    @test sin(ε)/ε === PowerNumber(1.0,0.0,0,0)
    @test sin(sqrt(ε))/sqrt(ε) === PowerNumber(1.0,0.0,0.0,0.0)
end

@testset "LogNumber" begin
    ε = PowerNumber(1.0,1.0)
    z = 1-ε
    @test log1p(-z) isa LogNumber
    @test exp(LogNumber(2,3)) == exp(3)*ε^2
    HypergeometricFunctions.expm1(LogNumber(2,3))
end

@testset "Rational" begin
    ε = PowerNumber(1.0,1.0)
    @test PowerNumber(1.,0,0.,1.) + PowerNumber(-1.,0,2.,2.) == PowerNumber(1.,0.,0.,1)
    @test_broken (1 + 1/ε + 1/ε^2) / (1 + 1/ε + 1/ε^2) == 1
end

@testset "HypergeometricFunctions" begin
    a,b,c = 1.154,1.2543,1.3543345
    ε = PowerNumber(1.0,1.0)
    z = 1-ε
    _₂F₁(a,b,c,z)

    a,b,c = 1.1,1.2,1.3
    @test_broken _₂F₁(a,b,c,z)
    log1p(-z)
end

@testset "IntervalArithmetic" begin
    a = @interval(1.0)
    p = PowerNumber(a,a,0,1)
    @test (a^4+a^2-a) == (p^4+p^2-p)(0)
end

#0.19999999999999996, 0.10000000000000009, 1.3, 1.0 + (-1.0)ϵ^1.0 + o(ϵ^1.0)
#(HypergeometricFunctions._₂F₁general)(0.10000000000000009, 0.19999999999999996, 1.3, 1.0 + (-1.0)ϵ^1.0 + o(ϵ^1.0))
#(HypergeometricFunctions._₂F₁Inf)(0.10000000000000009, 1.1, 1.3, (-1.0)ϵ^-1.0 + (1.0)ϵ^0.0 + o(ϵ^0.0))
#(HypergeometricFunctions.AInf)(0.10000000000000009, 1.1, 1.3, (-1.0)ϵ^1.0 + (-1.0)ϵ^2.0 + o(ϵ^2.0), 1, 0.0)
#z = (-1.0)ϵ^1.0 + (-1.0)ϵ^2.0 + o(ϵ^2.0)
# z = PowerNumber(-1.,-1.,1.,2.)
# @enter (HypergeometricFunctions.BInf)(0.10000000000000009, 1.1, 1.3, z, 1, 0.0)
# @enter (HypergeometricFunctions.recInfβ₀)(0.10000000000000009, 1.1, 1.3, z, 1, 0.0)

max_forks_repo_name: idlemat/PowerNumbers


text:
# Index manipulations
#---------------------
Base.@deprecate(
    permuteind(t::TensorMap,p1::IndexTuple, p2::IndexTuple=(); copy::Bool = false),
    permute(t, p1, p2; copy = copy))

"""
    permute(tsrc::AbstractTensorMap{S}, p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int} = ())
        -> tdst::TensorMap{S,N₁,N₂}

Permute the indices of `tsrc::AbstractTensorMap{S}` such that a new tensor
`tdst::TensorMap{S,N₁,N₂}` is obtained, with indices in `p1` playing the role of the
codomain or range of the map, and indices in `p2` indicating the domain.

To permute into an existing `tdst`, see [`add!`](@ref)
"""
function permute(t::TensorMap{S},
                    p1::IndexTuple, p2::IndexTuple=();
                    copy::Bool = false) where {S}
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))

    if !copy
        # share data if possible
        if p1 === codomainind(t) && p2 === domainind(t)
            return t
        elseif has_shared_permute(t, p1, p2)
            return TensorMap(reshape(t.data, dim(cod), dim(dom)), cod, dom)
        end
    end
    # general case
    @inbounds begin
        return add!(true, t, false, similar(t, cod←dom), p1, p2)
    end
end

function permute(t::AdjointTensorMap{S}, p1::IndexTuple{N₁}, p2::IndexTuple{N₂}=();
                    copy::Bool = false) where {S,N₁,N₂}
    p1′ = map(n->adjointtensorindex(t, n), p2)
    p2′ = map(n->adjointtensorindex(t, n), p1)
    adjoint(permute(adjoint(t), p1′, p2′; copy = copy))
end

function has_shared_permute(t::TensorMap, p1, p2)
    if p1 === codomainind(t) && p2 === domainind(t)
        return true
    elseif sectortype(t) === Trivial
        stridet = i->stride(t[], i)
        sizet = i->size(t[], i)
        canfuse1, d1, s1 = TO._canfuse(sizet.(p1), stridet.(p1))
        canfuse2, d2, s2 = TO._canfuse(sizet.(p2), stridet.(p2))
        return canfuse1 && canfuse2 && s1 == 1 && (d2 == 1 || s2 == d1)
    else
        return false
    end
end

function has_shared_permute(t::AdjointTensorMap, p1, p2)
    p1′ = adjointtensorindices(t, p2)
    p2′ = adjointtensorindices(t, p1)
    return has_shared_permute(t', p1′, p2′)
end

Base.@deprecate(permuteind!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, p1, p2),
                permute!(tdst, tsrc, p1, p2))

"""
    permute(tsrc::AbstractTensorMap{S}, p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int} = ())
        -> tdst::TensorMap{S,N₁,N₂}

Permute the indices of `tsrc::AbstractTensorMap{S}` such that a new tensor
`tdst::TensorMap{S,N₁,N₂}` is obtained, with indices in `p1` playing the role of the
codomain or range of the map, and indices in `p2` indicating the domain.

To permute into an existing `tdst`, see [`add!`](@ref)
"""
@propagate_inbounds Base.permute!(tdst::AbstractTensorMap{S,N₁,N₂},
                                    tsrc::AbstractTensorMap{S},
                                    p1::IndexTuple{N₁},
                                    p2::IndexTuple{N₂}=()) where {S,N₁,N₂} =
    add!(true, tsrc, false, tdst, p1, p2)

# Braid
function braid(t::TensorMap{S}, levels::IndexTuple,
                    p1::IndexTuple, p2::IndexTuple=();
                    copy::Bool = false) where {S}
    @assert length(levels) == numind(t)
    if BraidingStyle(sectortype(S)) isa SymmetricBraiding
        return permute(t, p1, p2; copy = copy)
    end
    if !copy && p1 == codomainind(t) && p2 == domainind(t)
        return t
    end
    # general case
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
    @inbounds begin
        return add!(true, t, false, similar(t, cod←dom), p1, p2, levels)
    end
end
@propagate_inbounds braid!(tdst::AbstractTensorMap{S,N₁,N₂},
                                tsrc::AbstractTensorMap{S},
                                levels::IndexTuple,
                                p1::IndexTuple{N₁},
                                p2::IndexTuple{N₂}=()) where {S,N₁,N₂} =
    add!(true, tsrc, false, tdst, p1, p2, levels)

# Transpose
function LinearAlgebra.transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap)
    codomain(tdst) == domain(tsrc)' && domain(tdst) == codomain(tsrc)' ||
        throw(SpaceMismatch())
    levels = (codomainind(tsrc)..., domainind(tsrc)...)
    braid!(tdst, tsrc, levels, reverse(domainind(tsrc)), reverse(codomainind(tsrc)))
    if BraidingStyle(sectortype(tdst)) != Bosonic()
        for (c,b) in blocks(tdst)
            rmul!(b, twist(c))
        end
    end
    return tdst
end

LinearAlgebra.transpose(t::AbstractTensorMap) =
    transpose!(similar(t, domain(t)', codomain(t)'), t)

# Twist
twist(t::AbstractTensorMap, i::Int; inv::Bool = false) = twist!(copy(t), i; inv = inv)

function twist!(t::AbstractTensorMap, i::Int; inv::Bool = false)
    N₁ = numout(t)
    for (f1,f2) in fusiontrees(t)
        θ = i <= N₁ ? twist(f1.uncoupled[i]) : twist(f2.uncoupled[i-N₁])
        rmul!(t[f1,f2], θ)
    end
    return t
end

# Fusing and splitting


text:
#!/usr/bin/env julia

using ArgParse

function parseargs()
	s = ArgParseSettings()
	s.description=
	"""
	Program to do something amazing.
	"""
	@add_arg_table s begin
		"--reference","-r"
			help = "Reference sequence to use"
			required=true
		"--querydb","-q"
			help = "Query database of sequences"
			required=true
		"--snpsout","-s"
			help = "Name of output SNP file."
			required=true
		"--reverse"
			help = "Reverse complement sequence"
			action = :store_true
	end
	return parse_args(s)
end

function main()
	parsed_args = parseargs()
	println("Parsed args: ")
	for (arg,val) in parsed_args
		println("  $arg  ==>   $val")
	end

	println("Reference: ",parsed_args["reference"])
end

main()

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
const SDL_MIN_SINT8 = reinterpret(Int8, ~0x7F)
const SDL_MIN_SINT16 = reinterpret(Int16, ~0x7FFF)
const SDL_MIN_SINT32 = reinterpret(Int32, ~0x7FFFFFFF)
const SDL_MIN_SINT64 = reinterpret(Int64, ~0x7FFFFFFFFFFFFFFF)

const SDL_ICONV_ERROR = reinterpret(Csize_t, -1)
const SDL_ICONV_E2BIG = reinterpret(Csize_t, -2)
const SDL_ICONV_EILSEQ = reinterpret(Csize_t, -3)
const SDL_ICONV_EINVAL = reinterpret(Csize_t, -4)

const SDL_TOUCH_MOUSEID = reinterpret(UInt32, Int32(-1))


text:
module DependencyGraphTest

using Test
using ForneyLab: LinkedList, DependencyGraph, addVertex!, addEdge!, children, neighbors

mutable struct MockVertex
    id::Symbol
end

@testset "LinkedList" begin
    ll = LinkedList{typeof("test")}()
    @test isempty(ll) == true
    @test_throws Exception push!(ll, 3) # Incorrect element type
    push!(ll, "test1")
    push!(ll, "test2")
    @test isempty(ll) == false
    @test collect(ll) == ["test1"; "test2"]
end

@testset "DependencyGraph" begin
    @testset "DependencyGraph()" begin
        # should construct a dependency graph
        dg = DependencyGraph{MockVertex}()
        @test isa(dg, DependencyGraph)
        @test dg.vertices == MockVertex[]
        @test dg.neighbors == LinkedList[]

        # Construct:
        # (1)-->(2)

        # addVertex! should add a vertex
        v1 = MockVertex(:v1)
        v2 = MockVertex(:v2)
        addVertex!(dg, v1)
        addVertex!(dg, v2)
        @test dg.vertices == MockVertex[v1, v2]

        # addEdge! should construct an edge
        addEdge!(dg, v1, v2)
        @test !(v1 in neighbors(v2, dg))
        @test (v2 in neighbors(v1, dg)) # Neighbours is not a commutative relation
        @test !(v1 in neighbors(v1, dg)) # A vertex is not a neighbor of itself
        @test !(v2 in neighbors(v2, dg))

        # Extend:
        # (1)-->(2)-->(3)
        v3 = MockVertex(:v3)
        addVertex!(dg, v3)
        addEdge!(dg, v2, v3)
        @test !(v2 in neighbors(v3, dg))
        @test (v3 in neighbors(v2, dg))
        @test !(v1 in neighbors(v3, dg)) # Unconnected vertices are not neighbours
        @test !(v3 in neighbors(v1, dg))
    end

    @testset "children" begin
        # should find all children of a vertex and return a topologically sorted array

        # Construct:
        #
        #  (1)-->(4)<--(5)
        #   |     |     |
        #   v     v     v
        #  (2)   (6)   (7)
        #   |
        #   v
        #  (3)

        dg = DependencyGraph{MockVertex}()
        vertices = MockVertex[]
        for i = 1:7
            v = MockVertex(Symbol("v$i"))
            push!(vertices, v)
            addVertex!(dg, v)
        end
        addEdge!(dg, vertices[1], vertices[2])
        addEdge!(dg, vertices[2], vertices[3])
        addEdge!(dg, vertices[1], vertices[4])
        addEdge!(dg, vertices[4], vertices[6])
        addEdge!(dg, vertices[5], vertices[4])
        addEdge!(dg, vertices[5], vertices[7])

        @test children(vertices[1], dg) == MockVertex[vertices[3], vertices[2], vertices[6], vertices[4], vertices[1]]
        @test children(vertices[1], dg, breaker_sites=Set{MockVertex}([vertices[4]])) == MockVertex[vertices[3], vertices[2], vertices[1]]
        @test children(vertices[1], dg, restrict_to=Set{MockVertex}(vertices[1:4])) == MockVertex[vertices[3], vertices[2], vertices[4], vertices[1]]
    end

    @testset "Loopy dependengy graph" begin
        # should find all children of a vertex and return a topologically sorted array in a cyclic graph

        # Construct:
        #
        #  (1)-->(2)
        #   ^    /
        #    \  v
        #     (3)

        dg = DependencyGraph{MockVertex}()
        vertices = MockVertex[]
        for i = 1:3
            v = MockVertex(Symbol("v$i"))
            push!(vertices, v)
            addVertex!(dg, v)
        end
        addEdge!(dg, vertices[1], vertices[2])
        addEdge!(dg, vertices[2], vertices[3])
        addEdge!(dg, vertices[3], vertices[1])

        @test (vertices[1] in neighbors(vertices[3], dg))
        @test !(vertices[2] in neighbors(vertices[3], dg))
        @test (vertices[2] in neighbors(vertices[1], dg))
        @test !(vertices[3] in neighbors(vertices[1], dg))
        @test (vertices[3] in neighbors(vertices[2], dg))
        @test !(vertices[1] in neighbors(vertices[2], dg))

        @test_throws Exception children(vertices[1], dg) # children function errors when cycles are not allowed
        @test children(vertices[1], dg, allow_cycles=true) == MockVertex[vertices[3], vertices[2], vertices[1]]
    end
end

end # module


text:
abstract type PrecipitationBC end

"""
    OutflowPrecipitation() :: PrecipitationBC

Free flux out of the domain.
"""
struct OutflowPrecipitation <: PrecipitationBC end
function atmos_precipitation_boundary_state!(
    nf,
    bc_precipitation::OutflowPrecipitation,
    atmos,
    _...,
) end
function atmos_precipitation_normal_boundary_flux_second_order!(
    nf,
    bc_precipitation::OutflowPrecipitation,
    atmos,
    _...,
) end


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

struct NLSolversFG!{F<:GradientFunction}
    gradfunc::F
end

NLSolversBase.only_fg!(f::GradientFunction) = NLSolversBase.only_fg!(NLSolversFG!(f))

function (gf::NLSolversFG!)(val_f::Real, grad_f::Any, v::Any)
    return gf.gradfunc(!, grad_f, v)
end

function (gf::NLSolversFG!)(val_f::Nothing, grad_f::Any, v::Any)
    gf.gradfunc(!, grad_f, v)
    return Nothing
end

function _bat_findmode_impl_optim(target::AnySampleable, algorithm::AbstractModeEstimator)
    density_notrafo = convert(AbstractDensity, target)
    trafoalg = bat_default(bat_transform, Val(:algorithm), algorithm.trafo, density_notrafo)
    shaped_density, trafo = bat_transform(algorithm.trafo, density_notrafo, trafoalg)
    shape = varshape(shaped_density)
    density = unshaped(shaped_density)

    rng = bat_determ_rng()
    x_init = collect(unshaped(bat_initval(rng, density, apply_trafo_to_init(trafo, algorithm.init)).result))

    f = -logdensityof(density)
    r_optim = Optim.MaximizationWrapper(_run_optim(f, x_init, algorithm))
    mode_trafo_unshaped = Optim.minimizer(r_optim.res)
    mode_trafo = shape(mode_trafo_unshaped)
    mode_notrafo = inv(trafo)(mode_trafo)
    (result = mode_notrafo, result_trafo = mode_trafo, trafo = trafo, info = r_optim)
end

# Wrapper for type stability of optimize result (why does this work?):
function _optim_optimize(f, x0::AbstractArray, method::Optim.AbstractOptimizer, options = Optim.Options(); kwargs...)
    Optim.optimize(f, x0, method, options; kwargs...)
end

"""
    struct MaxDensityNelderMead <: AbstractModeEstimator

Estimates the mode of a probability density using Nelder-Mead optimization
(currently via [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl),
subject to change).

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MaxDensityNelderMead{
    TR<:AbstractDensityTransformTarget,
    IA<:InitvalAlgorithm
} <: AbstractModeEstimator
    trafo::TR = NoDensityTransform()
    init::IA = InitFromTarget()
end
export MaxDensityNelderMead

function _run_optim(f::Function, x_init::AbstractArray{<:Real}, algorithm::MaxDensityNelderMead)
    _optim_optimize(f, x_init, Optim.NelderMead())
end

bat_findmode_impl(target::AnySampleable, algorithm::MaxDensityNelderMead) = _bat_findmode_impl_optim(target, algorithm)

"""
    struct MaxDensityLBFGS <: AbstractModeEstimator

Estimates the mode of a probability density using LBFGS optimization
(currently via [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl),
subject to change).

The gradient of the target density is computed via auto-differentiation.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MaxDensityLBFGS{
    TR<:AbstractDensityTransformTarget,
    IA<:InitvalAlgorithm
} <: AbstractModeEstimator
    trafo::TR = PriorToGaussian()
    init::IA = InitFromTarget()
end
export MaxDensityLBFGS

bat_findmode_impl(target::AnySampleable, algorithm::MaxDensityLBFGS) = _bat_findmode_impl_optim(target, algorithm)

function _run_optim(f::Function, x_init::AbstractArray{<:Real}, algorithm::MaxDensityLBFGS)
    fg = valgradof(f)
    _optim_optimize(Optim.only_fg!(fg), x_init, Optim.LBFGS())
end


text:
mutable struct pose_t <: LCMType
    utime::Int64
    pos::SVector{3, Float64}
    vel::SVector{3, Float64}
    orientation::SVector{4, Float64}
    rotation_rate::SVector{3, Float64}
    accel::SVector{3, Float64}
end

@lcmtypesetup(pose_t)


text:
# Test layers and data/model movements on and off the GPU
# Add tests for layers and their gradients on the GPU
# Most of the forward passes should be fine being applied
# to bitstype objects, but this gives higher coverage for our use-cases
# Check that getting the gradients does not throw

# generic movement tests
@testset "Basic GPU Movement" begin
  @test gradient(x -> sum(gpu(x)), rand(3,3)) isa Tuple
  @test gradient(x -> sum(cpu(x)), gpu(rand(3,3))) isa Tuple
end

# TODO: These layers get into scalar indexing
# `AlphaDropout` throws a compilation error on GPUs,
# whereas, the rest are scalar indexing issues.
const BROKEN_LAYERS = Union{DepthwiseConv,
                            AlphaDropout,
                            InstanceNorm,
                            GroupNorm}

function gpu_gradtest(name::String, layers::Vector, x_cpu=nothing, args...; test_cpu=true)
  isnothing(x_cpu) && error("Missing input to test the layers against.")
  @testset "$name GPU grad tests" begin
    for layer in layers
      @testset "$layer GPU grad test" begin

        # compute output and grad of parameters
        l_cpu = layer(args...)
        ps_cpu = Flux.params(l_cpu)
        y_cpu, back_cpu = pullback(() -> sum(l_cpu(x_cpu)), ps_cpu)
        gs_cpu = back_cpu(1f0)

        x_gpu = gpu(x_cpu)
        l_gpu = l_cpu |> gpu
        ps_gpu = Flux.params(l_gpu)

        if l_gpu isa BROKEN_LAYERS
          @test_broken gradient(() -> sum(l_gpu(x_gpu)), ps_gpu) isa Flux.Zygote.Grads
        else
          y_gpu, back_gpu = pullback(() -> sum(l_gpu(x_gpu)), ps_gpu)
          gs_gpu = back_gpu(1f0) # TODO many layers error out when backprop int 1, should fix

          # compute grad of input
          xg_cpu = gradient(x -> sum(l_cpu(x)), x_cpu)[1]
          xg_gpu = gradient(x -> sum(l_gpu(x)), x_gpu)[1]

          # test
          if test_cpu
            @test y_gpu ≈ y_cpu   rtol=1e-4 atol=1e-4
            @test Array(xg_gpu) ≈ xg_cpu   rtol=1e-4 atol=1e-4
          end
          @test gs_gpu isa Flux.Zygote.Grads
          for (p_cpu, p_gpu) in zip(ps_cpu, ps_gpu)
            @test gs_gpu[p_gpu] isa Flux.CUDA.CuArray
            if test_cpu
              @test Array(gs_gpu[p_gpu]) ≈ gs_cpu[p_cpu]   rtol=1e-4 atol=1e-4
            end
          end
        end
      end
    end
  end
end

# Just to give testset in gradtest meaningful labels
ConvNoBias(args...) = Conv(args...; bias=false)
ConvTransposeNoBias(args...) = ConvTranspose(args...; bias=false)
CrossCorNoBias(args...) = CrossCor(args...; bias=false)
DepthwiseConvNoBias(args...) = DepthwiseConv(args...;bias=false)
r = rand(Float32, 28, 28, 1, 1)
conv_layers = [Conv, ConvNoBias, ConvTranspose, ConvTransposeNoBias, CrossCor, CrossCorNoBias, DepthwiseConv, DepthwiseConvNoBias]
gpu_gradtest("Conv", conv_layers, r, (2,2), 1=>3)

pooling_layers = [MaxPool, MeanPool]
gpu_gradtest("Pooling", pooling_layers, r, (2,2))

adaptive_pooling_layers = [AdaptiveMaxPool, AdaptiveMeanPool]
gpu_gradtest("AdaptivePooling", adaptive_pooling_layers, r, (7,7))

dropout_layers = [Dropout, AlphaDropout]
gpu_gradtest("Dropout", dropout_layers, r, 0.5f0; test_cpu=false) # dropout is not deterministic

layer_norm = [LayerNorm]
gpu_gradtest("LayerNorm 1", layer_norm, rand(Float32, 28,28,3,4), 1, test_cpu=false) #TODO fix errors
gpu_gradtest("LayerNorm 2", layer_norm, rand(Float32, 5,4), 5)

batch_norm = [BatchNorm]
gpu_gradtest("BatchNorm 1", batch_norm, rand(Float32, 28,28,3,4), 3, test_cpu=false) #TODO fix errors
gpu_gradtest("BatchNorm 2", batch_norm, rand(Float32, 5,4), 5)

instancenorm = [InstanceNorm]
gpu_gradtest("InstanceNorm", instancenorm, r, 1)

groupnorm = [GroupNorm]
gpu_gradtest("GroupNorm", groupnorm, rand(Float32, 28,28,3,1), 3, 1)

@testset "function layers" begin
  x = rand(3,3)
  gpu_gradtest(x -> sum(Flux.normalise(x; dims=1)), x)
  gpu_gradtest(x -> sum(Flux.normalise(x; dims=2)), x)
  gpu_gradtest(x -> sum(Flux.normalise(x)), x)
end

@testset "Zeros mapped for $cl" for cl in (Conv, ConvTranspose, CrossCor, DepthwiseConv)
  l = cl((2,2), 1=>3, bias = false) |> gpu
  ip = zeros(Float32, 28,28,1,1) |> gpu
  if l isa BROKEN_LAYERS
    @test_broken sum(l(ip)) ≈ 0.f0
    @test_broken gradient(() -> sum(l(ip)), Flux.params(l)) isa Flux.Zygote.Grads
  else
    @test sum(l(ip)) ≈ 0.f0
    gs = gradient(() -> sum(l(ip)), Flux.params(l))
    @test l.bias ∉ gs.params
  end
end

@testset "Dense with Zeros bias" begin
  l = Dense(ones(Float32, 4,3), Flux.Zeros()) |> gpu
  ip = zeros(Float32, 3, 7) |> gpu

  @test sum(l(ip)) ≈ 0.f0
  gs = gradient(() -> sum(l(ip)), Flux.params(l))
  @test l.b ∉ gs.params
end


text:
########################################################################################
# Algorithm
########################################################################################
abstract type GaussianODEFilter <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm end
abstract type AbstractEK <: GaussianODEFilter end

"""
    EK0(; prior=:ibm, order=1, diffusionmodel=:dynamic, smooth=true)

Gaussian ODE filtering with zeroth order extended Kalman filter.

Currently, only the integrated Brownian motion prior `:ibm` is supported.
For the diffusionmodel, chose one of
`[:dynamic, :dynamicMV, :fixed, :fixedMV, :fixedMAP]`.

See also: [`EK1`](@ref)

# References:
- M. Schober, S. Särkkä, and P. Hennig: **A Probabilistic Model for the Numerical Solution
  of Initial Value Problems**
- F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig: **Probabilistic Solutions To Ordinary
  Differential Equations As Non-Linear Bayesian Filtering: A New Perspective**
"""
Base.@kwdef struct EK0 <: AbstractEK
    prior::Symbol = :ibm
    order::Int = 1
    diffusionmodel::Symbol = :dynamic
    smooth::Bool = true
end

"""
    EK1(; prior=:ibm, order=1, diffusionmodel=:dynamic, smooth=true)

Gaussian ODE filtering with first order extended Kalman filter

Currently, only the integrated Brownian motion prior `:ibm` is supported.
For the diffusionmodel, chose one of
`[:dynamic, :dynamicMV, :fixed, :fixedMV, :fixedMAP]`.

See also: [`EK0`](@ref)

# References:
- F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig: **Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective**
"""
Base.@kwdef struct EK1 <: AbstractEK
    prior::Symbol = :ibm
    order::Int = 1
    diffusionmodel::Symbol = :dynamic
    smooth::Bool = true
end


text:
# https://adventofcode.com/2020/day/6

function part1(input)
    groups = map(x -> Set(replace(x, r"[^a-z]" => "")), split(input, "\n\n"))
    return sum(map(length, groups))
end

example = "
abc

a
b
c

ab
ac

a
a
a
a

b
"

@assert part1(example) == 11

function part2(input)
    groups = split(input, "\n\n")
    groups = map(s -> filter(ch -> ch != "", split(s, '\n')), groups)
    common_items = map(g -> reduce(intersect, g), groups)
    return sum(map(length, common_items))
end

@assert part2(example) == 6

test = read("data/day-06.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 6530
@assert result2 == 3323

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
#=
Visual illustration of the central limit theorem in 3d

@author : Spencer Lyon <spencer.lyon@nyu.edu>

References
----------
Based off the original python file clt3d.py
=#
using PyPlot
using Distributions
using KernelDensity

beta_dist = Beta(2.0, 2.0)

function gen_x_draws(k)
    bdraws = rand(beta_dist, (3, k))

    # == Transform rows, so each represents a different distribution == #
    bdraws[1, :] -= 0.5
    bdraws[2, :] += 0.6
    bdraws[3, :] -= 1.1

    # == Set X[i] = bdraws[j, i], where j is a random draw from {1, 2, 3} == #
    js = rand(1:3, k)
    X = Array(Float64, k)
    for i=1:k
        X[i]=  bdraws[js[i], i]
    end

    # == Rescale, so that the random variable is zero mean == #
    m, sigma = mean(X), std(X)
    return (X .- m) ./ sigma
end

nmax = 5
reps = 100000
ns = 1:nmax

# == Form a matrix Z such that each column is reps independent draws of X == #
Z = Array(Float64, reps, nmax)
for i=ns
    Z[:, i] = gen_x_draws(reps)
end

# == Take cumulative sum across columns
S = cumsum(Z, 2)

# == Multiply j-th column by sqrt j == #
Y = S .* (1. ./ sqrt(ns))'

# == Plot == #
fig = figure()
ax = fig[:gca](projection="3d")

a, b = -3, 3
gs = 100
xs = linspace(a, b, gs)

# == Build verts == #
greys = linspace(0.3, 0.7, nmax)
verts = {}
for n=ns
    kde_est = kde(Y[:, n])
    _xs, ys = kde_est.x, kde_est.density
    ax[:plot3D](_xs, ys, n, zdir="x", color=string(greys[n]),
                fillstyle="bottom")
end

ax[:set_xlim3d](1, nmax)
ax[:set_xticks](ns)
ax[:set_xlabel]("n")
ax[:set_yticks]((-3, 0, 3))
ax[:set_ylim3d](a, b)
ax[:set_zlim3d](0, 0.4)
ax[:set_zticks]((0.2, 0.4))


text:
function get_generator_data(data::OPFData; use_gpu=false)
    ngen = length(data.generators)

    if use_gpu
        pgmin = CuArray{Float64}(undef, ngen)
        pgmax = CuArray{Float64}(undef, ngen)
        qgmin = CuArray{Float64}(undef, ngen)
        qgmax = CuArray{Float64}(undef, ngen)
        c2 = CuArray{Float64}(undef, ngen)
        c1 = CuArray{Float64}(undef, ngen)
        c0 = CuArray{Float64}(undef, ngen)
    else
        pgmin = Array{Float64}(undef, ngen)
        pgmax = Array{Float64}(undef, ngen)
        qgmin = Array{Float64}(undef, ngen)
        qgmax = Array{Float64}(undef, ngen)
        c2 = Array{Float64}(undef, ngen)
        c1 = Array{Float64}(undef, ngen)
        c0 = Array{Float64}(undef, ngen)
    end

    Pmin = Float64[data.generators[g].Pmin for g in 1:ngen]
    Pmax = Float64[data.generators[g].Pmax for g in 1:ngen]
    Qmin = Float64[data.generators[g].Qmin for g in 1:ngen]
    Qmax = Float64[data.generators[g].Qmax for g in 1:ngen]
    coeff0 = Float64[data.generators[g].coeff[3] for g in 1:ngen]
    coeff1 = Float64[data.generators[g].coeff[2] for g in 1:ngen]
    coeff2 = Float64[data.generators[g].coeff[1] for g in 1:ngen]
    copyto!(pgmin, Pmin)
    copyto!(pgmax, Pmax)
    copyto!(qgmin, Qmin)
    copyto!(qgmax, Qmax)
    copyto!(c0, coeff0)
    copyto!(c1, coeff1)
    copyto!(c2, coeff2)

    return pgmin,pgmax,qgmin,qgmax,c2,c1,c0
end

function get_bus_data(data::OPFData; use_gpu=false)
    nbus = length(data.buses)

    FrIdx = [l for b=1:nbus for l in data.FromLines[b]]
    ToIdx = [l for b=1:nbus for l in data.ToLines[b]]
    GenIdx = [g for b=1:nbus for g in data.BusGenerators[b]]
    FrStart = accumulate(+, vcat([1], [length(data.FromLines[b]) for b=1:nbus]))
    ToStart = accumulate(+, vcat([1], [length(data.ToLines[b]) for b=1:nbus]))
    GenStart = accumulate(+, vcat([1], [length(data.BusGenerators[b]) for b=1:nbus]))

    Pd = Float64[data.buses[i].Pd for i=1:nbus]
    Qd = Float64[data.buses[i].Qd for i=1:nbus]

    if use_gpu
        cuFrIdx = CuArray{Int}(undef, length(FrIdx))
        cuToIdx = CuArray{Int}(undef, length(ToIdx))
        cuGenIdx = CuArray{Int}(undef, length(GenIdx))
        cuFrStart = CuArray{Int}(undef, length(FrStart))
        cuToStart = CuArray{Int}(undef, length(ToStart))
        cuGenStart = CuArray{Int}(undef, length(GenStart))
        cuPd = CuArray{Float64}(undef, nbus)
        cuQd = CuArray{Float64}(undef, nbus)

        copyto!(cuFrIdx, FrIdx)
        copyto!(cuToIdx, ToIdx)
        copyto!(cuGenIdx, GenIdx)
        copyto!(cuFrStart, FrStart)
        copyto!(cuToStart, ToStart)
        copyto!(cuGenStart, GenStart)
        copyto!(cuPd, Pd)
        copyto!(cuQd, Qd)

        return cuFrStart,cuFrIdx,cuToStart,cuToIdx,cuGenStart,cuGenIdx,cuPd,cuQd
    else
        return FrStart,FrIdx,ToStart,ToIdx,GenStart,GenIdx,Pd,Qd
    end
end

function get_branch_data(data::OPFData; use_gpu=false)
    buses = data.buses
    lines = data.lines
    BusIdx = data.BusIdx
    nline = length(data.lines)
    ybus = Ybus{Array{Float64}}(computeAdmitances(data.lines, data.buses, data.baseMVA; VI=Array{Int}, VD=Array{Float64})...)
    frBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].from]].Vmin^2, buses[BusIdx[lines[l].from]].Vmax^2) ]
    toBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].to]].Vmin^2, buses[BusIdx[lines[l].to]].Vmax^2) ]

    if use_gpu
        cuYshR = CuArray{Float64}(undef, length(ybus.YshR))
        cuYshI = CuArray{Float64}(undef, length(ybus.YshI))
        cuYffR = CuArray{Float64}(undef, nline)
        cuYffI = CuArray{Float64}(undef, nline)
        cuYftR = CuArray{Float64}(undef, nline)
        cuYftI = CuArray{Float64}(undef, nline)
        cuYttR = CuArray{Float64}(undef, nline)
        cuYttI = CuArray{Float64}(undef, nline)
        cuYtfR = CuArray{Float64}(undef, nline)
        cuYtfI = CuArray{Float64}(undef, nline)
        cuFrBound = CuArray{Float64}(undef, 2*nline)
        cuToBound = CuArray{Float64}(undef, 2*nline)
        copyto!(cuYshR, ybus.YshR)
        copyto!(cuYshI, ybus.YshI)
        copyto!(cuYffR, ybus.YffR)
        copyto!(cuYffI, ybus.YffI)
        copyto!(cuYftR, ybus.YftR)
        copyto!(cuYftI, ybus.YftI)
        copyto!(cuYttR, ybus.YttR)
        copyto!(cuYttI, ybus.YttI)
        copyto!(cuYtfR, ybus.YtfR)
        copyto!(cuYtfI, ybus.YtfI)
        copyto!(cuFrBound, frBound)
        copyto!(cuToBound, toBound)

        return cuYshR, cuYshI, cuYffR, cuYffI, cuYftR, cuYftI,
               cuYttR, cuYttI, cuYtfR, cuYtfI, cuFrBound, cuToBound
    else
        return ybus.YshR, ybus.YshI, ybus.YffR, ybus.YffI, ybus.YftR, ybus.YftI,
               ybus.YttR, ybus.YttI, ybus.YtfR, ybus.YtfI, frBound, toBound
    end
end

function init_solution!(env::AdmmEnv, sol::SolutionOneLevel, ybus::Ybus, rho_pq, rho_va)
    data, model = env.data, env.model

    lines = data.lines
    buses = data.buses
    BusIdx = data.BusIdx
    ngen = length(data.generators)
    nline = length(data.lines)

    YffR = ybus.YffR; YffI = ybus.YffI
    YttR = ybus.YttR; YttI = ybus.YttI
    YftR = ybus.YftR; YftI = ybus.YftI
    YtfR = ybus.YtfR; YtfI = ybus.YtfI

    for g=1:ngen
        pg_idx = model.gen_mod.gen_start + 2*(g-1)
        CUDA.@allowscalar sol.v_curr[pg_idx] = 0.5*(data.generators[g].Pmin + data.generators[g].Pmax)
        CUDA.@allowscalar sol.v_curr[pg_idx+1] = 0.5*(data.generators[g].Qmin + data.generators[g].Qmax)
    end

    sol.rho .= rho_pq
    fill!(sol.u_curr, 0.0)
    fill!(sol.v_curr, 0.0)

    for l=1:nline
        wij0 = (buses[BusIdx[lines[l].from]].Vmax^2 + buses[BusIdx[lines[l].from]].Vmin^2) / 2
        wji0 = (buses[BusIdx[lines[l].to]].Vmax^2 + buses[BusIdx[lines[l].to]].Vmin^2) / 2
        wR0 = sqrt(wij0 * wji0)

        pij_idx = model.line_start + 8*(l-1)
        CUDA.@allowscalar sol.u_curr[pij_idx] = YffR[l] * wij0 + YftR[l] * wR0
        CUDA.@allowscalar sol.u_curr[pij_idx+1] = -YffI[l] * wij0 - YftI[l] * wR0
        CUDA.@allowscalar sol.u_curr[pij_idx+2] = YttR[l] * wji0 + YtfR[l] * wR0
        CUDA.@allowscalar sol.u_curr[pij_idx+3] = -YttI[l] * wji0 - YtfI[l] * wR0
        #=
        u_curr[pij_idx+4] = wij0
        u_curr[pij_idx+5] = wji0
        u_curr[pij_idx+6] = 0.0
        u_curr[pij_idx+7] = 0.0
        =#
        # wRIij[2*(l-1)+1] = wR0
        # wRIij[2*l] = 0.0

        CUDA.@allowscalar sol.v_curr[pij_idx+4] = wij0
        CUDA.@allowscalar sol.v_curr[pij_idx+5] = wji0
        CUDA.@allowscalar sol.v_curr[pij_idx+6] = 0.0
        CUDA.@allowscalar sol.v_curr[pij_idx+7] = 0.0

        sol.rho[pij_idx+4:pij_idx+7] .= rho_va
    end

    sol.l_curr .= 0
    return
end

function copy_data_kernel(n::Int, dest::CuDeviceArray{Float64,1}, src::CuDeviceArray{Float64,1})
    tx = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if tx <= n
        dest[tx] = src[tx]
    end
    return
end

function update_multiplier_kernel(n::Int, l_curr::CuDeviceArray{Float64,1},
                                  u_curr::CuDeviceArray{Float64,1},
                                  v_curr::CuDeviceArray{Float64,1},
                                  rho::CuDeviceArray{Float64,1})
    tx = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if tx <= n
        l_curr[tx] += rho[tx] * (u_curr[tx] - v_curr[tx])
    end
    return
end

function primal_residual_kernel(n::Int, rp::CuDeviceArray{Float64,1},
                                u_curr::CuDeviceArray{Float64,1},
                                v_curr::CuDeviceArray{Float64,1})
    tx = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if tx <= n
        rp[tx] = u_curr[tx] - v_curr[tx]
    end

    return
end

function dual_residual_kernel(n::Int, rd::CuDeviceArray{Float64,1},
                              v_prev::CuDeviceArray{Float64,1},
                              v_curr::CuDeviceArray{Float64,1},
                              rho::CuDeviceArray{Float64,1})
    tx = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if tx <= n
        rd[tx] = -rho[tx] * (v_curr[tx] - v_prev[tx])
    end

    return
end

function check_linelimit_violation(data::OPFData, u)
    lines = data.lines
    nline = length(data.lines)
    line_start = 2*length(data.generators) + 1

    rateA_nviols = 0
    rateA_maxviol = 0.0
    rateC_nviols = 0
    rateC_maxviol = 0.0

    for l=1:nline
        pij_idx = line_start + 8*(l-1)
        ij_val = 0.0
        ji_val = 0.0
        CUDA.@allowscalar ij_val = u[pij_idx]^2 + u[pij_idx+1]^2
        CUDA.@allowscalar ji_val = u[pij_idx+2]^2 + u[pij_idx+3]^2

        limit = (lines[l].rateA / data.baseMVA)^2
        if limit > 0 && limit < 1e10
            if ij_val > limit || ji_val > limit
                rateA_nviols += 1
                rateA_maxviol = max(rateA_maxviol, max(ij_val - limit, ji_val - limit))
            end
        end

        limit = (lines[l].rateC / data.baseMVA)^2
        if limit > 0 && limit < 1e10
            if ij_val > limit || ji_val > limit
                rateC_nviols += 1
                rateC_maxviol = max(rateC_maxviol, max(ij_val - limit, ji_val - limit))
            end
        end
    end
    rateA_maxviol = sqrt(rateA_maxviol)
    rateC_maxviol = sqrt(rateC_maxviol)

    return rateA_nviols, rateA_maxviol, rateC_nviols, rateC_maxviol
end

"""
    admm_restart!

This function restarts the ADMM with a given `env::AdmmEnv` containing solutions and all the other parameters.
"""
admm_restart!(env::AdmmEnv; options...) = admm_solve!(env, env.solution; options...)

function admm_solve!(env::AdmmEnv, sol::SolutionOneLevel; iterlim=800, scale=1e-4)
    if env.use_gpu
        CUDA.device!(env.gpu_no)
    end

    data, par, mod = env.data, env.params, env.model

    shift_lines = 0
    shmem_size = sizeof(Float64)*(14*mod.n+3*mod.n^2) + sizeof(Int)*(4*mod.n)

    nblk_gen = div(mod.gen_mod.ngen, 32, RoundUp)
    nblk_br = mod.nline
    nblk_bus = div(mod.nbus, 32, RoundUp)

    it = 0
    time_gen = time_br = time_bus = 0.0

    @time begin
    while it < iterlim
        it += 1

        if !env.use_gpu
            sol.u_prev .= sol.u_curr
            sol.v_prev .= sol.v_curr
            sol.l_prev .= sol.l_curr

            tcpu = generator_kernel(mod.gen_mod, data.baseMVA, sol.u_curr, sol.v_curr, sol.l_curr, sol.rho)
            time_gen += tcpu.time

            if env.use_polar
                tcpu = @timed auglag_it, tron_it = polar_kernel_cpu(mod.n, mod.nline, mod.line_start, scale,
                                                                    sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                    shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                    mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            else
                tcpu = @timed auglag_it, tron_it = auglag_kernel_cpu(mod.n, mod.nline, it, par.max_auglag, mod.line_start, par.mu_max,
                                                                     sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                     shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                     mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            end
            time_br += tcpu.time

            tcpu = @timed bus_kernel_cpu(data.baseMVA, mod.nbus, mod.gen_mod.gen_start, mod.line_start,
                                         mod.FrStart, mod.FrIdx, mod.ToStart, mod.ToIdx, mod.GenStart,
                                         mod.GenIdx, mod.Pd, mod.Qd, sol.u_curr, sol.v_curr, sol.l_curr, sol.rho, mod.YshR, mod.YshI)
            time_bus += tcpu.time

            sol.l_curr .+= sol.rho .* (sol.u_curr .- sol.v_curr)
            sol.rd .= -sol.rho .* (sol.v_curr .- sol.v_prev)
            sol.rp .= sol.u_curr .- sol.v_curr
            #sol.rp_old .= sol.u_prev .- sol.v_prev

            primres = norm(sol.rp)
            dualres = norm(sol.rd)

            eps_pri = sqrt(length(sol.l_curr))*par.ABSTOL + par.RELTOL*max(norm(sol.u_curr), norm(-sol.v_curr))
            eps_dual = sqrt(length(sol.u_curr))*par.ABSTOL + par.RELTOL*norm(sol.l_curr)

            (par.verbose > 0) && @printf("[CPU] %10d  %.6e  %.6e  %.6e  %.6e  %6.2f  %6.2f\n",
                    it, primres, dualres, eps_pri, eps_dual, auglag_it, tron_it)

            if primres <= eps_pri && dualres <= eps_dual
                break
            end
        else
            @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) copy_data_kernel(mod.nvar, sol.u_prev, sol.u_curr)
            @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) copy_data_kernel(mod.nvar, sol.v_prev, sol.v_curr)
            @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) copy_data_kernel(mod.nvar, sol.l_prev, sol.l_curr)
            CUDA.synchronize()

            tgpu = generator_kernel(mod.gen_mod, data.baseMVA, sol.u_curr, sol.v_curr, sol.l_curr, sol.rho)

            time_gen += tgpu.time
            if env.use_polar
                tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_br shmem=shmem_size polar_kernel(mod.n, mod.nline, mod.line_start, scale,
                                                                                sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                                shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                                mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            else
                tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_br shmem=shmem_size auglag_kernel(mod.n, it, par.max_auglag, mod.line_start, scale, par.mu_max,
                                                                                sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                                shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                                mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            end
            time_br += tgpu.time
            tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_bus bus_kernel(data.baseMVA, mod.nbus, mod.gen_mod.gen_start, mod.line_start,
                                                                           mod.FrStart, mod.FrIdx, mod.ToStart, mod.ToIdx, mod.GenStart,
                                                                           mod.GenIdx, mod.Pd, mod.Qd, sol.u_curr, sol.v_curr, sol.l_curr,
                                                                           sol.rho, mod.YshR, mod.YshI)
            time_bus += tgpu.time

            @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) update_multiplier_kernel(mod.nvar, sol.l_curr, sol.u_curr, sol.v_curr, sol.rho)
            @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) primal_residual_kernel(mod.nvar, sol.rp, sol.u_curr, sol.v_curr)
            @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) dual_residual_kernel(mod.nvar, sol.rd, sol.v_prev, sol.v_curr, sol.rho)
            CUDA.synchronize()

            gpu_primres = CUDA.norm(sol.rp)
            gpu_dualres = CUDA.norm(sol.rd)

            gpu_eps_pri = sqrt(length(sol.l_curr))*par.ABSTOL + par.RELTOL*max(CUDA.norm(sol.u_curr), CUDA.norm(sol.v_curr))
            gpu_eps_dual = sqrt(length(sol.u_curr))*par.ABSTOL + par.RELTOL*CUDA.norm(sol.l_curr)

            (par.verbose > 0) && @printf("[GPU] %10d  %.6e  %.6e  %.6e  %.6e\n", it, gpu_primres, gpu_dualres, gpu_eps_pri, gpu_eps_dual)

            if gpu_primres <= gpu_eps_pri && gpu_dualres <= gpu_eps_dual
                break
            end
        end
    end
    end

    u_curr = zeros(mod.nvar)
    copyto!(u_curr, sol.u_curr)
    objval = sum(data.generators[g].coeff[data.generators[g].n-2]*(data.baseMVA*u_curr[mod.gen_mod.gen_start+2*(g-1)])^2 +
                 data.generators[g].coeff[data.generators[g].n-1]*(data.baseMVA*u_curr[mod.gen_mod.gen_start+2*(g-1)]) +
                 data.generators[g].coeff[data.generators[g].n]
                 for g in 1:mod.gen_mod.ngen)::Float64
    sol.objval = objval

    if it < iterlim
        sol.status = HAS_CONVERGED
    else
        sol.status = MAXIMUM_ITERATIONS
    end

    if par.verbose > 0
        rateA_nviols, rateA_maxviol, rateC_nviols, rateC_maxviol = check_linelimit_violation(data, u_curr)
        @printf(" ** Line limit violations\n")
        @printf("RateA number of violations = %d (%d)\n", rateA_nviols, mod.nline)
        @printf("RateA maximum violation    = %.2f\n", rateA_maxviol)
        @printf("RateC number of violations = %d (%d)\n", rateC_nviols, mod.nline)
        @printf("RateC maximum violation    = %.2f\n", rateC_maxviol)

        @printf(" ** Time\n")
        @printf("Generator = %.2f\n", time_gen)
        @printf("Branch    = %.2f\n", time_br)
        @printf("Bus       = %.2f\n", time_bus)
        @printf("Total     = %.2f\n", time_gen + time_br + time_bus)

        @printf("Objective value = %.6e\n", objval)
    end
    return
end

function admm_rect_gpu(case::String; iterlim=800, rho_pq=400.0, rho_va=40000.0, scale=1e-4,
                       use_gpu=false, use_polar=true, gpu_no=0, verbose=1)
    env = AdmmEnv(case, use_gpu, rho_pq, rho_va; use_polar=use_polar, gpu_no=gpu_no, verbose=verbose,)
    admm_restart!(env, iterlim=iterlim, scale=scale)
    return env
end

# TODO: This needs revised to use AdmmEnv.
function admm_rect_gpu_mpi(
    case::String;
    iterlim=800, rho_pq=400.0, rho_va=40000.0, scale=1e-4, use_gpu=false, use_polar=true, gpu_no=1,
    comm::MPI.Comm=MPI.COMM_WORLD,
)
    data = opf_loaddata(case)

    ngen = length(data.generators)
    nline = length(data.lines)
    nbus = length(data.buses)
    nvar = 2*ngen + 8*nline
    if use_gpu
        @assert MPI.has_cuda()
    end

    # MPI settings
    root = 0
    is_master = MPI.Comm_rank(comm) == root
    n_processes = MPI.Comm_size(comm)
    nlines_local = div(nline, n_processes, RoundUp)
    nlines_padded = n_processes * nlines_local

    nvar_padded = 2*ngen + 8 * nlines_padded

    baseMVA = data.baseMVA
    n = (use_polar == true) ? 4 : 10
    mu_max = 1e8
    rho_max = 1e6
    rho_min_pq = 5.0
    rho_min_w = 5.0
    eps_rp = 1e-4
    eps_rp_min = 1e-5
    rt_inc = 2.0
    rt_dec = 2.0
    eta = 0.99
    Kf = 100
    Kf_mean = 10

    if use_gpu
        CUDA.device!(MPI.Comm_rank(comm) % CUDA.ndevices())
    end

    ybus = Ybus{Array{Float64}}(computeAdmitances(data.lines, data.buses, data.baseMVA; VI=Array{Int}, VD=Array{Float64})...)

    pgmin, pgmax, qgmin, qgmax, c2, c1, c0 = get_generator_data(data)
    YshR, YshI, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, FrBound, ToBound = get_branch_data(data)
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx, Pd, Qd = get_bus_data(data)

    cu_pgmin, cu_pgmax, cu_qgmin, cu_qgmax, cu_c2, cu_c1, cu_c0 = get_generator_data(data; use_gpu=use_gpu)
    cuYshR, cuYshI, cuYffR, cuYffI, cuYftR, cuYftI, cuYttR, cuYttI, cuYtfR, cuYtfI, cuFrBound, cuToBound = get_branch_data(data; use_gpu=use_gpu)
    cu_FrStart, cu_FrIdx, cu_ToStart, cu_ToIdx, cu_GenStart, cu_GenIdx, cu_Pd, cu_Qd = get_bus_data(data; use_gpu=use_gpu)

    gen_start = 1
    line_start = 2*ngen + 1

    # Allocations
    u_curr = zeros(nvar_padded)
    v_curr = zeros(nvar_padded)
    l_curr = zeros(nvar_padded)
    u_prev = zeros(nvar_padded)
    v_prev = zeros(nvar_padded)
    l_prev = zeros(nvar_padded)
    rho = zeros(nvar_padded)
    rd = zeros(nvar_padded)
    rp = zeros(nvar_padded)
    rp_old = zeros(nvar_padded)
    rp_k0 = zeros(nvar_padded)
    param = zeros(31, nlines_padded)
    wRIij = zeros(2*nline)

    init_values(data, ybus, gen_start, line_start,
                rho_pq, rho_va, u_curr, v_curr, l_curr, rho, wRIij)

    if use_gpu
        cu_u_curr = CuArray{Float64}(undef, nvar_padded)
        cu_v_curr = CuArray{Float64}(undef, nvar_padded)
        cu_l_curr = CuArray{Float64}(undef, nvar_padded)
        cu_u_prev = CuArray{Float64}(undef, nvar_padded)
        cu_v_prev = CuArray{Float64}(undef, nvar_padded)
        cu_l_prev = CuArray{Float64}(undef, nvar_padded)
        cu_rho = CuArray{Float64}(undef, nvar_padded)
        cu_rd = CuArray{Float64}(undef, nvar_padded)
        cu_rp = CuArray{Float64}(undef, nvar_padded)
        cu_rp_old = CuArray{Float64}(undef, nvar_padded)
        cu_rp_k0 = CuArray{Float64}(undef, nvar_padded)
        cuParam = CuArray{Float64}(undef, (31, nlines_padded))
        cuWRIij = CuArray{Float64}(undef, 2*nline)

        copyto!(cu_u_curr, u_curr)
        copyto!(cu_v_curr, v_curr)
        copyto!(cu_l_curr, l_curr)
        copyto!(cu_rho, rho)
        copyto!(cuParam, param)
        copyto!(cuWRIij, wRIij)
    end

    # MPI: Global info
    if use_gpu
        u_lines_root = @view cu_u_curr[line_start:end]
        l_lines_root = @view cu_l_curr[line_start:end]
        v_lines_root = @view cu_v_curr[line_start:end]
        rho_lines_root = cu_rho[line_start:end]
    else
        u_lines_root = @view u_curr[line_start:end]
        l_lines_root = @view l_curr[line_start:end]
        v_lines_root = @view v_curr[line_start:end]
        rho_lines_root = @view rho[line_start:end]
    end
    # MPI: Local info
    # We need only to transfer info about lines
    if use_gpu
        u_local = CUDA.zeros(Float64, 8 * nlines_local)
        v_local = CUDA.zeros(Float64, 8 * nlines_local)
        l_local = CUDA.zeros(Float64, 8 * nlines_local)
        rho_local = CUDA.zeros(Float64, 8 * nlines_local)
    else
        u_local = zeros(8 * nlines_local)
        v_local = zeros(8 * nlines_local)
        l_local = zeros(8 * nlines_local)
        rho_local = zeros(8 * nlines_local)
    end

    MPI.Scatter!(u_lines_root, u_local, root, comm)
    MPI.Scatter!(rho_lines_root, rho_local, root, comm)

    max_auglag = 50

    nblk_gen = div(ngen, 32, RoundUp)
    nblk_br = nline
    nblk_br_local = nlines_local
    nblk_bus = div(nbus, 32, RoundUp)

    ABSTOL = 1e-6
    RELTOL = 1e-5

    it = 0
    time_gen = time_br = time_bus = 0
    time_br_scatter = time_br_gather = 0

    h_u_curr = zeros(nvar)
    h_param = zeros(31, nline)
    h_wRIij = zeros(2*nline)

    shift_lines = MPI.Comm_rank(comm) * nlines_local
    # GPU settings
    shmem_size = sizeof(Float64)*(14*n+3*n^2) + sizeof(Int)*(4*n)

    while it < iterlim
        it += 1

        # CPU code
        if !use_gpu

            if is_master
                u_prev .= u_curr
                v_prev .= v_curr
                l_prev .= l_curr
                tcpu = @timed generator_kernel_cpu(baseMVA, ngen, gen_start, u_curr, v_curr, l_curr, rho,
                                                pgmin, pgmax, qgmin, qgmax, c2, c1)
                time_gen += tcpu.time
            end

            # MPI routines to be implemented:
            #  - Broadcast cu_v_curr and cu_l_curr to GPUs.
            #  - Collect cu_u_curr.
            #  - div(nblk_br / number of GPUs, RoundUp)
            # scatter / gather

            tcpu_mpi = @timed begin
                MPI.Scatter!(v_lines_root, v_local, root, comm)
                MPI.Scatter!(l_lines_root, l_local, root, comm)
            end
            time_br_scatter += tcpu_mpi.time

            nlines_actual = min(nlines_local, nline - shift_lines)
            tcpu = @timed auglag_it, tron_it = polar_kernel_cpu(n, nlines_actual, 1, scale,
                                                                u_local, v_local, l_local, rho_local,
                                                                shift_lines, param, YffR, YffI, YftR, YftI,
                                                                YttR, YttI, YtfR, YtfI, FrBound, ToBound)

            time_br += tcpu.time
            tcpu_mpi = @timed MPI.Gather!(u_local, u_lines_root, root, comm)
            time_br_gather += tcpu_mpi.time

            if is_master
                tcpu = @timed bus_kernel_cpu(baseMVA, nbus, gen_start, line_start,
                                            FrStart, FrIdx, ToStart, ToIdx, GenStart,
                                            GenIdx, Pd, Qd, u_curr, v_curr, l_curr, rho, YshR, YshI)
                time_bus += tcpu.time

                l_curr .+= rho .* (u_curr .- v_curr)
                rd .= -rho .* (v_curr .- v_prev)
                rp .= u_curr .- v_curr
                #rp_old .= u_prev .- v_prev

                primres = norm(rp)
                dualres = norm(rd)

                eps_pri = sqrt(length(l_curr))*ABSTOL + RELTOL*max(norm(u_curr), norm(-v_curr))
                eps_dual = sqrt(length(u_curr))*ABSTOL + RELTOL*norm(l_curr)

                @printf("[CPU] %10d  %.6e  %.6e  %.6e  %.6e  %6.2f  %6.2f\n",
                        it, primres, dualres, eps_pri, eps_dual, auglag_it, tron_it)
            end
        # GPU code
        else
            if is_master
                @cuda threads=64 blocks=(div(nvar-1, 64)+1) copy_data_kernel(nvar, cu_u_prev, cu_u_curr)
                @cuda threads=64 blocks=(div(nvar-1, 64)+1) copy_data_kernel(nvar, cu_v_prev, cu_v_curr)
                @cuda threads=64 blocks=(div(nvar-1, 64)+1) copy_data_kernel(nvar, cu_l_prev, cu_l_curr)
                CUDA.synchronize()

                tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_gen generator_kernel(baseMVA, ngen, gen_start,
                                                                                cu_u_curr, cu_v_curr, cu_l_curr, cu_rho,
                                                                                cu_pgmin, cu_pgmax, cu_qgmin, cu_qgmax, cu_c2, cu_c1)

                time_gen += tgpu.time
            end

            #  - Broadcast cu_v_curr and cu_l_curr to GPUs.
            tgpu_mpi = @timed begin
                MPI.Scatter!(v_lines_root, v_local, root, comm)
                MPI.Scatter!(l_lines_root, l_local, root, comm)
            end
            time_br_scatter += tgpu_mpi.time

            #  - div(nblk_br / number of GPUs, RoundUp)
            nblk_br_actual = min(nblk_br_local, nline - shift_lines)
            tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_br_actual shmem=shmem_size polar_kernel(n, nline, 1, scale,
                u_local, v_local, l_local, rho_local,
                shift_lines, cuParam, cuYffR, cuYffI, cuYftR, cuYftI,
                cuYttR, cuYttI, cuYtfR, cuYtfI, cuFrBound, cuToBound
            )
            time_br += tgpu.time

            #  - Collect cu_u_curr.
            tgpu_mpi = @timed MPI.Gather!(u_local, u_lines_root, root, comm)
            time_br_gather += tgpu_mpi.time

            if is_master
                tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_bus bus_kernel(baseMVA, nbus, gen_start, line_start,
                                                                            cu_FrStart, cu_FrIdx, cu_ToStart, cu_ToIdx, cu_GenStart,
                                                                            cu_GenIdx, cu_Pd, cu_Qd, cu_u_curr, cu_v_curr, cu_l_curr,
                                                                            cu_rho, cuYshR, cuYshI)
                time_bus += tgpu.time
                @cuda threads=64 blocks=(div(nvar-1, 64)+1) update_multiplier_kernel(nvar, cu_l_curr, cu_u_curr, cu_v_curr, cu_rho)
                @cuda threads=64 blocks=(div(nvar-1, 64)+1) primal_residual_kernel(nvar, cu_rp, cu_u_curr, cu_v_curr)
                @cuda threads=64 blocks=(div(nvar-1, 64)+1) dual_residual_kernel(nvar, cu_rd, cu_v_prev, cu_v_curr, cu_rho)
                CUDA.synchronize()

                gpu_primres = CUDA.norm(cu_rp)
                gpu_dualres = CUDA.norm(cu_rd)

                gpu_eps_pri = sqrt(length(l_curr))*ABSTOL + RELTOL*max(CUDA.norm(cu_u_curr), CUDA.norm(cu_v_curr))
                gpu_eps_dual = sqrt(length(u_curr))*ABSTOL + RELTOL*CUDA.norm(cu_l_curr)

                @printf("[GPU] %10d  %.6e  %.6e  %.6e  %.6e\n", it, gpu_primres, gpu_dualres, gpu_eps_pri, gpu_eps_dual)
            end
        end
    end

    if use_gpu
        copyto!(u_curr, cu_u_curr)
    end

    rank = MPI.Comm_rank(comm)
    @printf(" ** Time\n")
    @printf("[%d] Generator = %.2f\n", rank, time_gen)
    @printf("[%d] Branch    = %.2f\n", rank, time_br)
    @printf("[%d] Bus       = %.2f\n", rank, time_bus)
    @printf("[%d] G+Br+Bus  = %.2f\n", rank, time_gen + time_br + time_bus)
    @printf("[%d] Scatter   = %.2f\n", rank, time_br_scatter)
    @printf("[%d] Gather    = %.2f\n", rank, time_br_gather)
    @printf("[%d] MPI(S+G)  = %.2f\n", rank, time_br_scatter + time_br_gather)

    if is_master
        objval = sum(data.generators[g].coeff[data.generators[g].n-2]*(baseMVA*u_curr[gen_start+2*(g-1)])^2 +
                    data.generators[g].coeff[data.generators[g].n-1]*(baseMVA*u_curr[gen_start+2*(g-1)]) +
                    data.generators[g].coeff[data.generators[g].n]
                    for g in 1:ngen)
        @printf("Objective value = %.6e\n", objval)
    end

    return
end


text:
export ReLU

"""
$(TYPEDEF)

Represents a ReLU operation.

`p(x)` is shorthand for [`relu(x)`](@ref) when `p` is an instance of
`ReLU`.
"""
struct ReLU <: Layer
    tightening_algorithm::Union{TighteningAlgorithm, Nothing}
end

ReLU() = ReLU(nothing)

Base.hash(a::ReLU, h::UInt) = hash(:ReLU, h)

function Base.show(io::IO, p::ReLU)
    print(io, "ReLU()")
end

(p::ReLU)(x::Array{<:Real}) = relu(x)
(p::ReLU)(x::Array{<:JuMPLinearType}) = (Memento.info(MIPVerify.LOGGER, "Applying $p ..."); relu(x, nta = p.tightening_algorithm))