text:
## Activation functions
#
# Some of activation functions have its wrapper function for GPU in NNlibCUDA.jl.
# https://github.com/JuliaGPU/CuArrays.jl/issues/614

ACTIVATIONS = [
    :σ, :hardσ, :hardtanh, :relu,
    :leakyrelu, :relu6, :rrelu, :elu, :gelu, :swish, :selu,
    :celu, :softplus, :softsign, :logσ, :logcosh,
    :mish, :tanhshrink, :softshrink, :trelu, :lisht,
    :tanh_fast, :sigmoid_fast,
    ]

for f in ACTIVATIONS
    @eval export $(f)
end

# Aliases
export sigmoid, hardsigmoid, logsigmoid, thresholdrelu

# of type float (to allow for integer inputs)
oftf(x, y) = oftype(float(x), y)

"""
    σ(x) = 1 / (1 + exp(-x))

Classic [sigmoid](https://en.wikipedia.org/wiki/Sigmoid_function) activation
function.
Unicode `σ` can be entered as `\\sigma` then tab, in many editors.
The ascii name `sigmoid` is also exported.

See also [`sigmoid_fast`](@ref).

```
julia> lineplot(sigmoid, -5, 5, height=7)
          ┌────────────────────────────────────────┐
        1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⣀⡠⠤⠖⠒⠒⠋⠉⠉⠉⠉⠉⠉│ σ(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⡠⠖⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⣀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⡏⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡔⠋⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠊⠁⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⠤⠤⠤⠒⠊⠉⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> sigmoid === σ
true
```
"""
function σ(x)
    t = exp(-abs(x))
    ifelse(x ≥ 0, inv(1 + t), t / (1 + t))
end

const sigmoid = σ

"""
    hardσ(x) = max(0, min(1, (x + 3) / 6))

Piecewise linear approximation of [`sigmoid`](@ref).

```
julia> lineplot(hardsigmoid, -5, 5, height=7)
          ┌────────────────────────────────────────┐
        1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢀⡠⠖⠋⠉⠉⠉⠉⠉⠉⠉⠉│ hardσ(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⣀⡤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⡠⠔⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⡗⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡤⠖⠋⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⠤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> lineplot(sigmoid, -5, 5, height=7)
          ┌────────────────────────────────────────┐
        1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⣀⡠⠤⠖⠒⠒⠋⠉⠉⠉⠉⠉⠉│ σ(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⡠⠖⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⣀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⡏⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡔⠋⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠊⠁⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⠤⠤⠤⠒⠊⠉⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
hardσ(x) = clamp((x + 3) / 6, 0, 1)

# https://pytorch.org/docs/stable/generated/torch.nn.Hardsigmoid.html

const hardsigmoid = hardσ

"""
    logσ(x)

Return `log(σ(x))` which is computed in a numerically stable way.

```
julia> lineplot(logsigmoid, -5, 5, height=7)
           ┌────────────────────────────────────────┐
         0 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡧⠤⠔⠒⠒⠒⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│ logσ(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠖⠊⠉⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⢀⡤⠖⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⣀⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⡤⠖⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -6 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
logσ(x) = -softplus(-x)

const logsigmoid = logσ

"""
    hardtanh(x) = max(-1, min(1, x))

Segment-wise linear approximation of `tanh`, much cheaper to compute.
See ["Large Scale Machine Learning"](https://ronan.collobert.com/pub/matos/2004_phdthesis_lip6.pdf).

See also [`tanh_fast`](@ref).
```
julia> lineplot(hardtanh, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⣀⠔⠋⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│ hardtanh(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⣀⡤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⢀⡤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡤⡷⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠖⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠖⠋⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⠔⠋⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> lineplot(tanh, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⣀⡠⠤⠤⠒⠒⠒⠊⠉⠉⠉│ tanh(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⢀⡠⠔⠊⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⢀⡤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡤⡷⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡤⠖⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡠⠔⠊⠁⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⣀⣀⣀⡠⠤⠤⠤⠖⠒⠊⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
hardtanh(x) = clamp(x, oftype(x, -1), oftype(x, 1))  # clamp(x, -1, 1) is type-stable, but would promote Int32, for which we have tests

"""
    relu(x) = max(0, x)

[Rectified Linear Unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
activation function.

```
julia> lineplot(relu, -2, 2, height=7)
          ┌────────────────────────────────────────┐
        2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠋│ relu(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠊⠁⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡤⠊⠁⠀⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢀⡤⠖⠁⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⡠⠖⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⡠⠖⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣇⠔⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
relu(x) = ifelse(x<0, zero(x), x)  # faster than max(zero(x), x), still preserves NaN

"""
    leakyrelu(x, a=0.01) = max(a*x, x)

Leaky [Rectified Linear Unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
activation function.
You can also specify the coefficient explicitly, e.g. `leakyrelu(x, 0.01)`.

```julia
julia> lineplot(x -> leakyrelu(x, 0.5), -2, 2, height=7)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉│ #42(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊⠉⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⣀⡤⠖⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⣀⠤⠖⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⣤⣤⡤⡧⠶⠭⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⣀⠤⠤⠒⠒⠋⠉⠁⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⣀⣀⠤⠤⠒⠒⠊⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> leakyrelu(-10f0, 1//5)
-2.0f0

julia> leakyrelu(-10f0, 1//20)
-0.5f0
```
"""
leakyrelu(x, a=oftf(x, 0.01)) = ifelse(x>0, float(x), oftf(x, a*x))  # max(a*x, x) is 3x slower

"""
    relu6(x) = min(max(0, x), 6)

[Rectified Linear Unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
activation function capped at 6.
See ["Convolutional Deep Belief Networks"](https://www.cs.toronto.edu/~kriz/conv-cifar10-aug2010.pdf) from CIFAR-10.

```
julia> lineplot(relu6, -10, 10, height=7)
          ┌────────────────────────────────────────┐
        6 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠎⠉⠉⠉⠉⠉⠉⠉⠉│ relu6(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⡤⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⡠⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⠖⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⡔⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⡧⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-10⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀10
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
relu6(x) = clamp(x, oftype(x, 0), oftype(x, 6))  # clamp promotes, but clamp(x, 0, 6) would promote x::Int32

"""
    rrelu(x, lo=1/8, hi=1/3) = max(a*x, x)
    # where `a` is randomly sampled from uniform distribution `U(lo, hi)`

Randomized Leaky Rectified Linear Unit activation function.
See ["Empirical Evaluation of Rectified Activations"](https://arxiv.org/abs/1505.00853)
You can also specify the bound explicitly, e.g. `rrelu(x, 0.0, 1.0)`.

```julia
julia> lineplot(rrelu, -20, 10, height=7)
            ┌────────────────────────────────────────┐
         10 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠖⠋│ rrelu(x)
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⢀⡠⠖⠋⠁⠀⠀⠀│
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⢀⡠⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀│
   f(x)     │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡤⠤⣤⣤⢤⣤⣤⠤⠤⠤⢼⠮⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
            │⣰⢀⣆⡄⣄⡄⡠⡰⠦⠷⡜⢢⠷⠳⠢⠊⠉⠉⠀⠀⠁⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠃⠉⠙⠘⠃⠈⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -10 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            └────────────────────────────────────────┘
            ⠀-20⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀10
            ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> extrema(rrelu.(fill(-10f0, 1000)))
(-3.3316886f0, -1.2548422f0)
```
"""
function rrelu(x::T, l=1//8, u=1//3) where T<:Number
    a = (u - l) * rand(float(T)) + l
    return leakyrelu(x, a)
end

"""
    elu(x, α=1) = x > 0 ? x : α * (exp(x) - 1)

Exponential Linear Unit activation function.
See ["Fast and Accurate Deep Network Learning by Exponential Linear Units"](https://arxiv.org/abs/1511.07289).
You can also specify the coefficient explicitly, e.g. `elu(x, 1)`.

```
julia> lineplot(elu, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉│ elu(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊⠉⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⣀⡤⠖⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⣀⠤⠖⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡤⡧⠶⠭⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⣀⠤⠔⠒⠋⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⠤⠤⠤⠤⠔⠒⠒⠒⠊⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> elu(-10f0)
-0.9999546f0

julia> elu(-10f0, 2)
-1.9999092f0
```
"""
elu(x, α=1) = ifelse(x ≥ 0, float(x), α * (exp(x) - 1))

deriv_elu(Ω, α=1) = ifelse(Ω ≥ 0, one(Ω), Ω + α)

"""
    gelu(x) = 0.5x * (1 + tanh(√(2/π) * (x + 0.044715x^3)))

Activation function from ["Gaussian Error Linear Units"](https://arxiv.org/abs/1606.08415).

```
julia> lineplot(gelu, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊│ gelu(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠁⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⣀⡠⠤⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⣤⣤⣤⣤⣤⣤⣤⣤⡤⠤⠤⠤⠤⠤⠤⠤⣤⣤⣤⡤⡧⠶⠶⠭⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠉⠉⠉⠉⠉⠉⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
function gelu(x)
    α = oftf(x, 0.044715)
    λ = oftf(x, gelu_λ)
    x/2 * (1 + tanh(λ * (x + α * x^3)))
end

const gelu_λ = √(2 / π)

"""
    swish(x) = x * σ(x)

Self-gated activation function.
See ["Swish: a Self-Gated Activation Function"](https://arxiv.org/abs/1710.05941).

```
julia> lineplot(swish, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤│ swish(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠖⠋⠁⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠖⠋⠁⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⢀⣀⡤⠔⠊⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⣤⣤⡤⡧⠴⠶⠯⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠉⠑⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠉⠉⠉⠉⠁⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
swish(x) = x * σ(x)

"""
    lisht(x) = x * tanh(x)

Activation function from
["LiSHT: Non-Parametric Linearly Scaled Hyperbolic Tangent ..."](https://arxiv.org/abs/1901.05894)

```
julia> lineplot(lisht, -2, 2, height=7)
          ┌────────────────────────────────────────┐
        2 │⠢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔│ lisht(x)
          │⠀⠈⠑⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡤⠊⠁⠀│
          │⠀⠀⠀⠀⠈⠣⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠁⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠑⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠁⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠢⡄⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢀⠔⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⢄⡀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢀⡠⠖⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⠦⣄⣀⣀⣇⣀⣀⠤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
lisht(x) = x * tanh(x)

"""
    selu(x) = λ * (x ≥ 0 ? x : α * (exp(x) - 1))

    λ ≈ 1.05070...
    α ≈ 1.67326...

Scaled exponential linear units.
See ["Self-Normalizing Neural Networks"](https://arxiv.org/abs/1706.02515).

```
julia> lineplot(selu, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ selu(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡠⠤⠔⠒│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⣀⣀⠤⠔⠒⠋⠉⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⢀⣀⡤⠤⠒⠊⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⣉⡩⠭⠛⡏⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
           │⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⡤⠤⠔⠒⠊⠉⠁⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -2 │⠒⠒⠉⠉⠉⠉⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> selu(-10f0)
-1.7580194f0
```
"""
function selu(x)
    λ = oftf(x, selu_λ)
    α = oftf(x, selu_α)
    λ * ifelse(x > 0, x, α * (exp(x) - 1))
end

const selu_λ = 1.0507009873554804934193349852946
const selu_α = 1.6732632423543772848170429916717

function deriv_selu(Ω)
    λ = oftf(Ω, selu_λ)
    α = oftf(Ω, selu_α)
    ifelse(Ω > 0, λ, Ω + α * λ)
end

"""
    celu(x, α=1) = x ≥ 0 ? x : α * (exp(x/α) - 1)

Activation function from ["Continuously Differentiable Exponential Linear Units"](https://arxiv.org/abs/1704.07483).

```
julia> lineplot(celu, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉│ celu(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊⠉⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⣀⡤⠖⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⣀⠤⠖⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡤⡧⠶⠭⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⣀⠤⠔⠒⠋⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⠤⠤⠤⠤⠔⠒⠒⠒⠊⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> celu(-10f0)
-0.9999546f0
```
"""
celu(x, α=1) = ifelse(x ≥ 0, float(x), α * (exp(x/α) - 1))

"""
    trelu(x, theta=1) = x > theta ? x : 0

Threshold gated rectified linear activation function.
See ["Zero-bias autoencoders and the benefits of co-adapting features"](https://arxiv.org/abs/1402.3337)

```
julia> lineplot(trelu, -2, 4, height=7)
          ┌────────────────────────────────────────┐
        4 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠖⠋│ trelu(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠖⠋⠁⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠴⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⣠⠤⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⡏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣇⣀⣀⣀⣀⣀⣀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀4
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
trelu(x, theta=1) = ifelse(x <= theta, zero(x), x)

const thresholdrelu = trelu

"""
    softsign(x) = x / (1 + |x|)

See ["Quadratic Polynomials Learn Better Image Features"](http://www.iro.umontreal.ca/~lisa/publications2/index.php/attachments/single/205) (2009).

```
julia> lineplot(softsign, -5, 5, height=7)
           ┌────────────────────────────────────────┐
         1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⠤⠤⠤⠤⠤│ softsign(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⣀⡤⠖⠒⠋⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⡔⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡯⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠁⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⠤⠤⠒⠋⠁⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⠒⠒⠒⠒⠒⠊⠉⠉⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> lineplot!(ans, tanh)
           ┌────────────────────────────────────────┐
         1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⡤⠖⠊⠉⠉⠉⣉⣉⣉⣉⣉⠭⠭⠭⠭⠭│ softsign(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⡔⣃⡤⠖⠒⠋⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀│ tanh(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣧⡞⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡯⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡴⠃⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⠤⠤⠒⢋⠕⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -1 │⣒⣒⣒⣒⣒⣊⣉⣉⣉⣉⣁⣀⣀⡠⠤⠒⠁⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> softsign(1f0)
0.5f0

julia> softsign(100f0)
0.990099f0
```
"""
softsign(x) = x / (1 + abs(x))

"""
    softplus(x) = log(exp(x) + 1)

See ["Deep Sparse Rectifier Neural Networks"](http://proceedings.mlr.press/v15/glorot11a/glorot11a.pdf), JMLR 2011.

```
julia> lineplot(softplus, -3, 3, height=7)
          ┌────────────────────────────────────────┐
        4 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ softplus(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠔⠊⠁⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠔⠊⠁⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⣀⡠⠤⠒⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⡧⠤⠒⠊⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⡠⠤⠤⠤⠤⠔⠒⠒⠚⠉⠉⠁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-3⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀3
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> lineplot!(ans, relu)
          ┌────────────────────────────────────────┐
        4 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ softplus(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣠│ relu(x)
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣠⡴⠞⠋⠁│
   f(x)   │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣤⡴⠞⠋⠁⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⣀⡠⢤⡲⠝⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⡧⠤⠒⠊⣉⠥⠚⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⣀⣀⣀⣀⣀⣀⣀⣠⣤⣤⣤⣤⣔⣒⣒⣚⣉⣉⣁⣀⣇⠴⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-3⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀3
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> softplus(16f0)
16.0f0
```
"""
softplus(x) = log1p(exp(-abs(x))) + max(x, 0)

"""
    logcosh(x)

Return `log(cosh(x))` which is computed in a numerically stable way.

```
julia> lineplot(logcosh, -5, 5, height=7)
          ┌────────────────────────────────────────┐
        5 │⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ logcosh(x)
          │⠉⠢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠋│
          │⠀⠀⠀⠑⠢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠁⠀⠀│
   f(x)   │⠀⠀⠀⠀⠀⠀⠑⠦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠊⠁⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠦⡀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢀⡤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⠦⡀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢀⡤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        0 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠢⢄⣀⣀⣇⣀⡠⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
          ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
logcosh(x) = x + softplus(-2x) - oftf(x, log2)

const log2 = log(2)

"""
    mish(x) = x * tanh(softplus(x))

Activation function from ["Mish: A Self Regularized Non-Monotonic Neural Activation Function"](https://arxiv.org/abs/1908.08681).

```
julia> lineplot(mish, -5, 5, height=7)
           ┌────────────────────────────────────────┐
         5 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠖⠋│ mish(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡤⠒⠁⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠔⠋⠁⠀⠀⠀⠀⠀⠀│
   f(x)    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⢀⡠⠖⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⡤⠖⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣧⣔⣊⣁⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀│
        -1 │⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-5⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀5
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x
```
"""
mish(x) = x * tanh(softplus(x))

"""
    tanhshrink(x) = x - tanh(x)

See ["Tanhshrink Activation Function"](https://www.gabormelli.com/RKB/Tanhshrink_Activation_Function).

```
julia> lineplot(tanhshrink, -3, 3, height=7)
           ┌────────────────────────────────────────┐
         3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ tanhshrink(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡠⠤⠖⠊│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢀⣀⡠⠤⠒⠊⠉⠁⠀⠀⠀⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⣤⡤⠤⠤⠤⠤⠤⠤⡷⠶⠶⠶⠶⠶⠮⠭⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⣀⡠⠴⠒⠊⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⡠⠴⠒⠊⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-3⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀3
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> tanhshrink.((-10f0, 10f0))
(-9.0f0, 9.0f0)
```
"""
tanhshrink(x) = x - tanh(x)

"""
    softshrink(x, λ=0.5) =
        (x ≥ λ ? x - λ : (-λ ≥ x ? x + λ : 0))

See ["Softshrink Activation Function"](https://www.gabormelli.com/RKB/Softshrink_Activation_Function).

```
julia> lineplot(softshrink, -2, 2, height=7)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀│ softshrink(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠉⠁│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⣀⡤⠤⠒⠋⠁⠀⠀⠀⠀⠀⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⣤⡤⠤⠤⠤⠤⠤⠤⡧⠤⠤⠤⠤⠶⠮⠭⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⠀⠀⠀⠀⠀⢀⣀⠤⠖⠒⠉⠁⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⣀⠤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -2 │⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> lineplot!(ans, tanhshrink)
           ┌────────────────────────────────────────┐
         2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀│ softshrink(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⣉⡡│ tanhshrink(x)
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⣀⡤⠤⣒⣋⠥⠤⠒⠊⠉⠁⠀│
   f(x)    │⠤⠤⠤⠤⠤⠤⠤⠤⠤⣤⣤⣤⣤⡤⠤⠤⠤⠤⠤⠤⡷⠶⠶⠶⠶⠶⠾⠿⠯⠭⠭⠤⠤⠤⠤⠤⠤⠤⠤⠤│
           │⠀⢀⣀⡠⠤⠖⢒⣋⠭⠗⠒⠉⠁⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠊⣉⠤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
        -2 │⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           ⠀-2⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀2
           ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀x

julia> softshrink.((-10f0, 10f0))
(-9.5f0, 9.5f0)
```
"""
softshrink(x, λ=oftf(x, 0.5)) = min(max(0, x - λ), x + λ)

# Provide an informative error message if activation functions are called with an array
for f in ACTIVATIONS
    @eval $(f)(x::AbstractArray, args...) =
      error("Use broadcasting (`", $(string(f)), ".(x)`) to apply activation functions to arrays.")
end

## Faster, less accurate, versions of some.

"""
    tanh_fast(x)

This is a faster but slighly less accurate version of `tanh`.

Where Julia's `tanh` function has an error under 2 eps, this
may be wrong by 5 eps, a reduction by less than one decimal digit.

For `x::Float32` this is usually about 10 times faster,
with a smaller speedup for `x::Float64`.
For any other number types, it just calls `tanh`.

See also [`sigmoid_fast`](@ref).

```
julia> tanh(0.5f0)
0.46211717f0

julia> tanh_fast(0.5f0)
0.46211714f0

julia> hard_tanh(0.5f0)
0.5f0
```
"""
@inline function tanh_fast(x::Float32)
    x2 = abs2(x)
    n = evalpoly(x2, (1.0f0, 0.1346604f0, 0.0035974074f0, 2.2332108f-5, 1.587199f-8))
    d = evalpoly(x2, (1.0f0, 0.4679937f0, 0.026262015f0, 0.0003453992f0, 8.7767893f-7))
    ifelse(x2 < 66f0, x * (n / d), sign(x))
end

@inline function tanh_fast(x::Float64)
    exp2x = @fastmath exp(x + x)
    y = (exp2x - 1) / (exp2x + 1)
    # That has large errors near zero; using `expm1` would more accurate, but about as slow as `tanh`.
    # Instead, we switch to a polynomial, which is very accurate within its range:
    x2 = x * x
    ypoly = x * evalpoly(x2, (1.0, -0.33333333333324583, 0.13333333325511604, -0.05396823125794372, 0.02186660872609521, -0.008697141630499953))
    ifelse(x2 > 900.0, sign(y), ifelse(x2 < 0.017, oftype(y, ypoly), y))
end

# These approximations are very badly behaved for Float16; none are fast.
# They are also a bit slower with ForwardDiff.Dual numbers, let's use Base:
tanh_fast(x::Real) = Base.tanh(x)

"""
    sigmoid_fast(x)

This is a faster, and very slightly less accurate, version of `sigmoid`.
For `x::Float32, perhaps 3 times faster, and maximum errors 2 eps instead of 1.

See also [`tanh_fast`](@ref).

```
julia> sigmoid(0.2f0)
0.54983395f0

julia> sigmoid_fast(0.2f0)
0.54983395f0

julia> hardσ(0.2f0)
0.53333336f0
```
"""
@inline function sigmoid_fast(x::Real)
    t = @fastmath exp(-abs(x))
    y = ifelse(x ≥ 0, inv(1 + t), t / (1 + t))
    ifelse(x > 40, one(y), ifelse(x < -80, zero(y), y))
end
# For x::Float32, this is not as quick as the rational tanh_fast(x) above,
# but that polynomial has poor relative accuracy for negative x.

sigmoid_fast(x::Float16) = sigmoid(x)  # sigmoid_fast is extremely badly behaved at large x

## Define rrules for some activation functions, along with the
## broadcasted rrule activation functions.
## TODO: add to the lists below all activations.

## This is a performance hack specifically for Zygote, because it doesn't handle fused
## broadcasts well; but it generally should be good (or at least harmless) for any AD, as
## it saves ADing the broadcasting machinery.
## Related Issue https://github.com/JuliaDiff/ChainRulesCore.jl/issues/271

UNARY_ACTS = [ # f, df
    (:relu,         :(x > 0)),
    (:hardtanh,     :(-1 < x < 1)),
    (:selu,         :(deriv_selu(Ω))),
    (:σ,            :(conj(Ω * (1 - Ω)))),
    (:elu,          :(deriv_elu(Ω))),
    (:softplus,     :(σ(x))),

    (:tanh_fast,    :(conj(1 - Ω^2))),
    (:sigmoid_fast, :(conj(Ω * (1 - Ω)))),
    ]

for (f, df) in UNARY_ACTS
    @eval @scalar_rule($f(x), $df)

    pullback = Symbol(:broadcasted_, f, :_pullback)
    @eval function rrule(::typeof(broadcasted),
                         ::typeof($f), x::Numeric)
        Ω = $f.(x)
        function $pullback(Δ)
            x_thunk = InplaceableThunk(
                dx -> @.(dx += Δ * $df),
                @thunk @.(Δ * $df)
            )
            NoTangent(), NoTangent(), x_thunk
        end
        return Ω, $pullback
    end
end

BINARY_ACTS = [ # f, df1, df2
    (:elu, :(deriv_elu(Ω, x2)), :(NoTangent())), # TODO use real deriv instead of DNE
    ]

for (f, df1, df2) in BINARY_ACTS
    @eval @scalar_rule($f(x1, x2), ($df1, $df2))

    pullback = Symbol(:broadcasted_, f, :_pullback)
    @eval function rrule(::typeof(broadcasted),
                         ::typeof($f),
                         x1::Numeric, x2::Numeric)
        Ω = $f.(x1, x2)
        function $pullback(Δ)
            NoTangent(), NoTangent(), @.(Δ * $df1), @.(Δ * $df2)
        end
        return Ω, $pullback
    end
end


text:
fs2()

@testset "ignore/fs2" begin
    gotd()
    s = """
        @def ignore = ["foo.md", "path/foo.md", "dir/", "path/dir/"]
        """
    write(joinpath(td, "config.md"), s);
    F.process_config()
    @test F.globvar("ignore") == ["foo.md", "path/foo.md", "dir/", "path/dir/"]

    write(joinpath(td, "foo.md"), "anything")
    mkpath(joinpath(td, "path"))
    write(joinpath(td, "path", "foo.md"), "anything")
    mkpath(joinpath(td, "dir"))
    write(joinpath(td, "dir", "foo1.md"), "anything")
    mkpath(joinpath(td, "path", "dir"))
    write(joinpath(td, "path", "dir", "foo2.md"), "anything")
    write(joinpath(td, "index.md"), "standard things")
    watched = F.fd_setup()
    @test length(watched.md) == 1
    @test first(watched.md).first.second == "index.md"
end

fs1()

@testset "ignore/fs1" begin
    gotd()
    s = """
        @def ignore = ["foo.md", "path/foo.md", "dir/", "path/dir/"]
        """
    mkpath(joinpath(td, "src"))
    write(joinpath(td, "src", "config.md"), s);
    F.process_config()
    @test F.globvar("ignore") == ["foo.md", "path/foo.md", "dir/", "path/dir/"]

    write(joinpath(td, "src", "foo.md"), "anything")
    mkpath(joinpath(td, "src", "path"))
    write(joinpath(td, "src", "path", "foo.md"), "anything")
    mkpath(joinpath(td, "src", "dir"))
    write(joinpath(td, "src", "dir", "foo1.md"), "anything")
    mkpath(joinpath(td, "src", "path", "dir"))
    write(joinpath(td, "src", "path", "dir", "foo2.md"), "anything")
    write(joinpath(td, "src", "index.md"), "standard things")

    mkpath(joinpath(td, "src", "pages"))
    mkpath(joinpath(td, "src", "_css"))
    mkpath(joinpath(td, "src", "_html_parts"))

    watched = F.fd_setup()
    @test length(watched.md) == 1
    @test first(watched.md).first.second == "index.md"
end

ext: jl


text:
using Compat.Test
using Compat.LinearAlgebra
using CoordinateTransformations
using ForwardDiff: Dual, partials
using StaticArrays

@testset "CoordinateTransformations" begin

    include("core.jl")
    include("coordinatesystems.jl")
    include("affine.jl")
    include("perspective.jl")

end


text:
function sbb(borrow::UInt8, a::ArbDigit, b::ArbDigit)
    difference = SignedDoubleArbDigit(a) - SignedDoubleArbDigit(b) - SignedDoubleArbDigit(borrow)
    lowBits(difference), UInt8(difference < 0)
end

function sub2!(a::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit})
    len = min(length(a), length(b))
    a_lo, a_hi = @views a[begin:len], a[len+1:end]
    b_lo, b_hi = @views b[begin:len], b[len+1:end]

    borrow = 0x0
    for i in eachindex(a_lo)
        a_lo[i], borrow = sbb(borrow, a_lo[i], b_lo[i])
    end

    if !iszero(borrow)
        for i in eachindex(a_hi)
            a_hi[i], borrow = sbb(borrow, a_hi[i], zero(ArbDigit))
            iszero(borrow) && break
        end
    end

    (!iszero(borrow) || any(!iszero, b_hi)) && throw(ArgumentError("Cannot subtract b from a because b is larger than a."))
end

function sub!(a::ArbUInt, b::ArbUInt)
    sub2!(a.data, b.data)
    normalize!(a)
end
sub!(a::ArbUInt, b::T) where T <: Unsigned = sub!(a, ArbUInt(b))

function Base.:(-)(a::ArbUInt, b::ArbUInt)
    a_copy = deepcopy(a)
    sub2!(a_copy.data, b.data)
    a_copy
end
Base.:(-)(a::ArbUInt, b::T) where T <: Unsigned = sub!(deepcopy(a), b)
Base.:(-)(a::T, b::ArbUInt) where T <: Unsigned = sub!(deepcopy(b), a)


text:
# The following is based on doi 10.1109/22.275248
"
    abcd2s(ABCD(ω);Z01=50,Z02=50)

    returns a function of frequency given a function, ABCD(ω). The result
    possible depends on frequency and (possibly complex) source (Z01) and load
    (Z02) impedances - source and load impedances should be provided as
    functions of frequency.
"
function abcd2s(abcd_matrix;Z01=w->50, Z02=w->50)
    if size(abcd_matrix(0)) != (2,2)
        error("abcd_matrix must by a 2x2 matrix.")
    end
    A(w) = abcd_matrix(w)[1,1]
    B(w) = abcd_matrix(w)[1,2]
    C(w) = abcd_matrix(w)[2,1]
    D(w) = abcd_matrix(w)[2,2]
    R01(w) = real(Z01(w)); R02(w) = real(Z02(w));

    denominator(w) = A(w)*Z02(w) + B(w) + C(w)*Z01(w)*Z02(w) + D(w)*Z01(w)

    S11(w) = A(w)*Z02(w)+B(w)-C(w)*conj(Z01(w))*Z02(w)-D(w)*conj(Z01(w))

    S12(w) = 2*(A(w)*D(w)-B(w)*C(w))*sqrt(R01(w)*R02(w))
    S21(w) = 2*sqrt(R01(w)*R02(w))
    S22(w) = (-A(w)*conj(Z02(w))+B(w)-C(w)*Z01(w)*conj(Z02(w))+D(w)*Z01(w))
    return w -> (denominator(w))^-1*[S11(w) S12(w); S21(w) S22(w)]
end


text:
# This file is a part of SearchModels.jl

module SearchModels

include("utils.jl")
include("search.jl")

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BallNeighborhood(radius, metric=Euclidean())

A ball neighborhood with `radius` and `metric`.
"""
struct BallNeighborhood{T,M} <: AbstractNeighborhood
  radius::T
  metric::M

  function BallNeighborhood{T,M}(radius, metric) where {T,M}
    @assert radius > 0 "radius must be positive"
    new(radius, metric)
  end
end

BallNeighborhood(radius::T, metric::M=Euclidean()) where {T,M} =
  BallNeighborhood{T,M}(radius, metric)

"""
    radius(ball)

Return the radius of the `ball`.
"""
radius(ball::BallNeighborhood) = ball.radius

"""
    metric(ball)

Return the metric of the `ball`.
"""
metric(ball::BallNeighborhood) = ball.metric

isneighbor(ball::BallNeighborhood, xₒ::AbstractVector, x::AbstractVector) =
  evaluate(ball.metric, xₒ, x) ≤ ball.radius

# ------------
# IO methods
# ------------
function Base.show(io::IO, ball::BallNeighborhood)
  r = ball.radius
  d = ball.metric
  print(io, "BallNeighborhood($r, $d)")
end

function Base.show(io::IO, ::MIME"text/plain", ball::BallNeighborhood)
  println(io, "BallNeighborhood")
  println(io, "  radius: ", ball.radius)
  print(  io, "  metric: ", ball.metric)
end


text:
struct BoundaryPaddedMesh
    mesh::Any
    gridsize::Any
    nfmside::Any
    numghostlayers::Any
    bottomghostcoords::Any
    rightghostcoords::Any
    topghostcoords::Any
    leftghostcoords::Any
end

function uniform_grid_coordinates(x0, widths, npts)
    xrange = range(x0[1], stop = x0[1] + widths[1], length = npts[1])
    yrange = range(x0[2], stop = x0[2] + widths[2], length = npts[2])

    xcoords = repeat(xrange, inner = npts[2])
    ycoords = repeat(yrange, outer = npts[1])

    return vcat(xcoords', ycoords')
end

function bottom_ghost_coordinates(x0, xwidth, npts, dy, numlayers)

    return uniform_grid_coordinates(
        x0 - [0.0, numlayers * dy],
        [xwidth, (numlayers - 1) * dy],
        [npts, numlayers],
    )
end

function right_ghost_coordinates(x0, widths, npts, dx, numlayers)

    return uniform_grid_coordinates(
        x0 + [widths[1] + dx, 0.0],
        [(numlayers - 1) * dx, widths[2]],
        [numlayers, npts],
    )
end

function top_ghost_coordinates(x0, widths, npts, dy, numlayers)

    return uniform_grid_coordinates(
        x0 + [0.0, widths[2] + dy],
        [widths[1], (numlayers - 1) * dy],
        [npts, numlayers],
    )
end

function left_ghost_coordinates(x0, ywidth, npts, dx, numlayers)

    return uniform_grid_coordinates(
        x0 - [numlayers * dx, 0.0],
        [(numlayers - 1) * dx, ywidth],
        [numlayers, npts],
    )
end

function BoundaryPaddedMesh(mesh, numghostlayers)
    x0 = reference_corner(mesh)
    nfmside = nodes_per_mesh_side(mesh)
    meshwidths = widths(mesh)

    gridsize = meshwidths ./ (nfmside .- 1)

    bottomghostcoords =
        bottom_ghost_coordinates(x0, meshwidths[1], nfmside[1], gridsize[2], numghostlayers)
    rightghostcoords =
        right_ghost_coordinates(x0, meshwidths, nfmside[2], gridsize[1], numghostlayers)
    topghostcoords =
        top_ghost_coordinates(x0, meshwidths, nfmside[1], gridsize[2], numghostlayers)
    leftghostcoords =
        left_ghost_coordinates(x0, meshwidths[2], nfmside[2], gridsize[1], numghostlayers)

    return BoundaryPaddedMesh(
        mesh,
        gridsize,
        nfmside,
        numghostlayers,
        bottomghostcoords,
        rightghostcoords,
        topghostcoords,
        leftghostcoords,
    )
end

struct BoundaryPaddedLevelset
    interiordist::Any
    bottomghostdist::Any
    rightghostdist::Any
    topghostdist::Any
    leftghostdist::Any
    nfmside::Any
    gridsize::Any
    numghostlayers::Any
end

function BoundaryPaddedLevelset(
    paddedmesh::BoundaryPaddedMesh,
    refseedpoints,
    spatialseedpoints,
    seedcellids,
    levelset,
    levelsetcoeffs,
    cutmesh,
    tol,
    boundingradius,
)

    bottomghostdist = distance_to_zero_levelset(
        paddedmesh.bottomghostcoords,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        cutmesh,
        tol,
        boundingradius,
    )
    rightghostdist = distance_to_zero_levelset(
        paddedmesh.rightghostcoords,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        cutmesh,
        tol,
        boundingradius,
    )
    topghostdist = distance_to_zero_levelset(
        paddedmesh.topghostcoords,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        cutmesh,
        tol,
        boundingradius,
    )
    leftghostdist = distance_to_zero_levelset(
        paddedmesh.leftghostcoords,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        cutmesh,
        tol,
        boundingradius,
    )

    return BoundaryPaddedLevelset(
        levelsetcoeffs,
        bottomghostdist,
        rightghostdist,
        topghostdist,
        leftghostdist,
        paddedmesh.nfmside,
        paddedmesh.gridsize,
        paddedmesh.numghostlayers,
    )
end

function first_order_horizontal_backward_difference(paddedlevelset)
    interior = paddedlevelset.interiordist
    left = paddedlevelset.leftghostdist
    cols, rows = paddedlevelset.nfmside
    dx = paddedlevelset.gridsize[1]

    npts = length(interior)
    derivative = zeros(npts)

    derivative[1:rows] = (interior[1:rows] - left) / dx

    for col = 2:cols
        prevstart = (col - 2) * rows + 1
        prevstop = prevstart + rows - 1

        start = (col - 1) * rows + 1
        stop = start + rows - 1

        derivative[start:stop] = (interior[start:stop] - interior[prevstart:prevstop]) / dx
    end

    return derivative
end

function first_order_horizontal_forward_difference(paddedlevelset)
    interior = paddedlevelset.interiordist
    right = paddedlevelset.rightghostdist
    cols, rows = paddedlevelset.nfmside
    dx = paddedlevelset.gridsize[1]

    npts = length(interior)
    derivative = zeros(npts)

    for col = 1:cols-1
        start = (col - 1) * rows + 1
        stop = start + rows - 1
        nextstart = stop + 1
        nextstop = nextstart + rows - 1

        derivative[start:stop] = (interior[nextstart:nextstop] - interior[start:stop]) / dx
    end

    start = (cols - 1) * rows + 1
    stop = start + rows - 1
    derivative[start:stop] = (right - interior[start:stop]) / dx

    return derivative
end

function first_order_vertical_backward_difference(paddedlevelset)
    interior = paddedlevelset.interiordist
    bottom = paddedlevelset.bottomghostdist
    cols, rows = paddedlevelset.nfmside
    dy = paddedlevelset.gridsize[2]
    numghostlayers = paddedlevelset.numghostlayers

    npts = length(interior)
    derivative = zeros(npts)

    interioridx = range(1, step = rows, length = cols)
    ghostidx = range(numghostlayers, step = numghostlayers, length = cols)

    derivative[interioridx] = (interior[interioridx] - bottom[ghostidx]) / dy

    for row = 2:rows
        interioridx = range(row, step = rows, length = cols)
        previnterioridx = range(row - 1, step = rows, length = cols)

        derivative[interioridx] = (interior[interioridx] - interior[previnterioridx]) / dy
    end

    return derivative
end

function first_order_vertical_forward_difference(paddedlevelset)
    interior = paddedlevelset.interiordist
    top = paddedlevelset.topghostdist
    cols, rows = paddedlevelset.nfmside
    dy = paddedlevelset.gridsize[2]
    numghostlayers = paddedlevelset.numghostlayers

    npts = length(interior)
    derivative = zeros(npts)

    for row = 1:rows-1
        interioridx = range(row, step = rows, length = cols)
        nextinterioridx = range(row + 1, step = rows, length = cols)

        derivative[interioridx] = (interior[nextinterioridx] - interior[interioridx]) / dy
    end

    interioridx = range(rows, step = rows, length = cols)
    ghostidx = range(1, step = numghostlayers, length = cols)

    derivative[interioridx] = (top[ghostidx] - interior[interioridx]) / dy

    return derivative
end

function first_order_nabla(u1, u2, v1, v2)
    return sqrt.(
        (max.(u1, 0.0)) .^ 2 +
        (min.(u2, 0.0)) .^ 2 +
        (max.(v1, 0.0)) .^ 2 +
        (min.(v2, 0.0)) .^ 2,
    )
end

function step_first_order_levelset(paddedlevelset, speed, dt)
    Dmx = first_order_horizontal_backward_difference(paddedlevelset)
    Dpx = first_order_horizontal_forward_difference(paddedlevelset)
    Dmy = first_order_vertical_backward_difference(paddedlevelset)
    Dpy = first_order_vertical_forward_difference(paddedlevelset)

    delplus = first_order_nabla(Dmx, Dpx, Dmy, Dpy)
    delminus = first_order_nabla(Dpx, Dmx, Dpy, Dmy)

    return paddedlevelset.interiordist -
           dt * (max.(speed, 0.0) .* delplus + min.(speed, 0.0) .* delminus)
end


text:
include("intersect/linesegment-linesegment.jl")
include("intersect/linesegment-quadraticsegment.jl")
include("intersect/linesegment-rectilinear_grid.jl")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
struct DoFHandler{N,T}
    DoF::Vector{NTuple{N,T}}
    nDofs::T
end

Base.getindex(a::DoFHandler, I::Int64) = a.DoF[I]
Base.getindex(a::DoFHandler, I::Int32) = a.DoF[I]
Base.getindex(a::DoFHandler, I::Int8) = a.DoF[I]

function DoFs_Thermal(EL2NOD, nnodel)
    DoF = [ntuple(j->EL2NOD[j,i], nnodel) for i in axes(EL2NOD,2)]
    nDofs = maximum(view(EL2NOD,1:nnodel,:))
    DoFHandler(DoF,nDofs)
end

function DoFs_Pressure(e2nP)
    DoF = [ntuple(j->e2nP[j,i],3) for i in axes(e2nP,2)]
    nDofs = maximum(e2nP)
    DoFHandler(DoF, nDofs)
end

function DoFs_Velocity(EL2NOD)
    dummy = Matrix{Int64}(undef, size(EL2NOD, 2), 12)
    @views dummy[:, 1:2:end] .= (@. 2*(EL2NOD-1)+1)'
    @views dummy[:, 2:2:end] .= (@. 2*(EL2NOD))'

    DoF = [ntuple(j->dummy[i,j], 12) for i in axes(EL2NOD,2)]
    nDofs = maximum(EL2NOD)*2
    DoFHandler(DoF,nDofs)
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
using Revise
using ActiveSetMethods
using PowerModels, JuMP, GLPK, Ipopt
# using CPLEX

PowerModels.silence()

# choose an internal LP solver
lp_solver = GLPK.Optimizer
# lp_solver = optimizer_with_attributes(
#     CPLEX.Optimizer,
#     "CPX_PARAM_SCRIND" => 0,
#     "CPX_PARAM_THREADS" => 1,
# )

include("acwr.jl")

build_acp(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), ACPPowerModel, PowerModels.build_opf)
build_acr(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), ACRPowerModel, PowerModels.build_opf)
build_iv(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), IVRPowerModel, PowerModels.build_opf_iv)
build_dcp(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), DCPPowerModel, PowerModels.build_opf_iv)

function run_opf(data_file::String, max_iter::Int = 100, algorithm = "Line Search")
    pm = build_acp(data_file)
    # pm2 = build_acp(data_file)
    # JuMP.@objective(pm2.model, Min, 0)
    init_vars(pm)
    # init_vars_from_ipopt(pm, pm2)
    result = optimize_model!(pm, optimizer = optimizer_with_attributes(
        ActiveSetMethods.Optimizer,
        "external_optimizer" => lp_solver,
        "max_iter" => max_iter,
        "algorithm" => algorithm,
    ))
    return pm, result
end

include("init_opf.jl")

max_stars_repo_name: exanauts/ActiveSetMethods

max_issues_repo_name: exanauts/ActiveSetMethods

max_forks_repo_name: exanauts/ActiveSetMethods


text:
const Param = Pair{AbstractString,AbstractString}
const Params = Vector{Param}
const Entry = Pair{AbstractString,Params}
const Entries = Vector{Entry}

@testset "parse_header_value" begin
    headervalue = "value1;q=0.5;r=b,value2;q=0"
    parsed = Spirit.parse_header_value(headervalue)

    @test isa(parsed, Entries)
    @test length(parsed) == 2

    entry1 = parsed[1]
    entry2 = parsed[2]

    # @test first(entry1) == "value1"
    # @test length(last(entry1)) == 2
    # @test last(entry1)[1] == Param("q" => "0.5")
    # @test last(entry1)[2] == Param("r" => "b")
    @test parsed == Entries([
        "value1" => [
            "q" => "0.5",
            "r" => "b"
        ],
        "value2" => [
            "q" => "0"
        ]
    ])

    headervalue = " value1 ; q = 0.5 ; r = b ,  value2 ; q =  0"
    parsed = Spirit.parse_header_value(headervalue)

    @test isa(parsed, Entries)
    @test length(parsed) == 2

    entry1 = parsed[1]
    entry2 = parsed[2]

    # @test first(entry1) == "value1"
    # @test length(last(entry1)) == 2
    # @test last(entry1)[1] == Param("q" => "0.5")
    # @test last(entry1)[2] == Param("r" => "b")
    @test parsed == Entries([
        "value1" => [
            "q" => "0.5",
            "r" => "b"
        ],
        "value2" => [
            "q" => "0"
        ]
    ])

    # no parameters
    headervalue = "value1;q=0.5;r=b,value2;q=0"
    @test_throws HTTPError Spirit.parse_header_value(headervalue; parameterized=false)

    headervalue = "value1,value2"
    parsed = Spirit.parse_header_value(headervalue)
    @test parsed == Entries([
        "value1" => [
        ],
        "value2" => [
        ]
    ])

    headervalue = "value1"
    parsed = Spirit.parse_header_value(headervalue)
    @test parsed == Entries([
        "value1" => [
        ]
    ])

    headervalue = ""
    @test_throws HTTPError Spirit.parse_header_value(headervalue; allowempty=false)

    headervalue = "bird,"
    @test_throws HTTPError Spirit.parse_header_value(headervalue; allowempty=false)

    headervalue = "bird,,salmon"
    @test_throws HTTPError Spirit.parse_header_value(headervalue; allowempty=false)

    headervalue = "bird,,salmon"
    parsed = Spirit.parse_header_value(headervalue; allowempty=true)
    @test parsed == Entries([
        "bird" => [
        ],
        "" => [
        ],
        "salmon" => [
        ]
    ])

    headervalue = "bird,,salmon"
    parsed = Spirit.parse_header_value(headervalue; allowempty=true, filterempty=true)
    @test parsed == Entries([
        "bird" => [
        ],
        "salmon" => [
        ]
    ])

    headervalue = "198"
    parsed = Spirit.parse_header_value(headervalue; permittedchars=Set(UInt8('0'):UInt8('9')))
    @test parsed == Entries([
        "198" => [
        ],
    ])
end


text:
mutable struct CRP <: Distributions.DiscreteMultivariateDistribution
    N::Int
    α::Float64
end

CRP(N, α, G₀) = CRP(N, α)
Base.show(io::IO, d::CRP) = print(io, "CRP(", d.N, ", ", d.α, ")")

function Distributions.rand(rng::AbstractRNG, d::CRP)
    nk = Int[]
    D = ChineseRestaurantProcess(DirichletProcess(d.α), nk)
    z = rand(rng, D, d.N)
end

function Distributions.logpdf(d::CRP, z::AbstractVector{<:Int})
    nk, K = histogram(d, z)
    D = ChineseRestaurantProcess(DirichletProcess(d.α), nk)
    return mapreduce(zi -> logpdf(D, zi), +, z; init=0.0)
end

# `Product`s can be treated as an array of iid variables
function (c::GibbsConditional{V, L})(θ) where {
    V<:VarName, L<:LogLikelihood{<:CRP}}

    N = c.base.dist.N
    Ωs = [1:n for n = 1:N]
    conditionals = similar(Ωs, DiscreteNonParametric)

    for index in eachindex(Ωs, conditionals)
        sub_vn = DynamicPPL.VarName(c.vn, (DynamicPPL.getindexing(c.vn)..., (index,)))
        θs_on_support = fixvalues(θ, sub_vn => Ωs[index])
        logtable = map(θs_on_support) do θ′
            c.base(θ′) + reduce(+,
                                (_likelihood(vn, β, c, θ′) for (vn, β) in c.blanket if vn == sub_vn),
                                init=0.0)
        end
        conditionals[index] = DiscreteNonParametric(Ωs[index], softmax!(vec(logtable)))
    end

    return Product(conditionals)
end

function _likelihood(vn, β, c, θ)
    try
        return β(θ)
    catch
        G₀ = c.base.args[3](θ)
        return logpdf(G₀, θ[vn])
    end
end

function update_histogram!(histogram, bin)
    if bin > length(histogram)
        push!(histogram, 1)
    else
        histogram[bin] += 1
    end

    return histogram
end

function histogram(d, z)
    nk = Int[]
    K = 0

    for n in 1:d.N
        nk = update_histogram!(nk, z[n])
        K = max(K, z[n])
    end

    return nk, K
end


text:
import Base.Matrix
import Base.+
import Base.*

export SU2

"""
    Represents the group SU(2), {[α -̅β; β ̅α] | α, β ∈ ℂ, |α|² + |β|² = 1}.

fields: α and β.
"""
struct SU2 <: Group
    α::U1
    β::U1
    function SU2(α::U1, β::U1)
        new(α, β)
    end
    function SU2(a::Array{U1,2})
        @assert(size(a) == (2, 2), "The size of Matrix must be 2×2.")
        new(a[1], a[2])
    end
    function SU2(a::Array{Complex,2})
        @assert(size(a) == (2, 2), "The size of Matrix must be 2×2.")
        new(U1(a[1]), U1(a[2]))
    end
    function SU2(a::Array{Complex{Float64},2})
        @assert(size(a) == (2, 2), "The size of Matrix must be 2×2.")
        new(U1(a[1]), U1(a[2]))
    end
end

"""
    Matrix(a::SU2)

Construct a Matrix with the given number.
"""
function Matrix(a::SU2)
    convert(Array{U1,2}, [a.α -conj(a.β); a.β conj(a.α)])
end

"""
    Base.show(io::IO, a::SU2)

Print a string representation of the given number.
"""
function Base.show(io::IO, a::SU2)
    print(io, Matrix(a))
end

"""
    +(a::SU2, b::SU2)

Add two numbers.
"""
function +(a::SU2, b::SU2)
    convert(Array{Complex,2}, Matrix(a) + Matrix(b))
end

"""
    +(a::SU2, b::Array{U1,2})

Add two numbers.
"""
function +(a::SU2, b::Array{U1,2})
    convert(Array{Complex,2}, Matrix(a) + b)
end

"""
    +(a::Array{U1,2}, b::SU2)

Add two numbers.
"""
function +(a::Array{U1,2}, b::SU2)
    convert(Array{Complex,2}, a + Matrix(b))
end

"""
    +(a::SU2, b::Array{Complex,2})

Add two numbers.
"""
function +(a::SU2, b::Array{Complex,2})
    Matrix(a) + b
end

"""
    +(a::Array{Complex,2}, b::SU2)

Add two numbers.
"""
function +(a::Array{Complex,2}, b::SU2)
    a + Matrix(b)
end

"""
    *(a::SU2, b::SU2)

Multiply two numbers.
"""
function *(a::SU2, b::SU2)
    SU2(Matrix(a) * Matrix(b))
end

"""
    *(a::SU2, b::Array{U1,2})

Multiply two numbers.
"""
function *(a::SU2, b::Array{U1,2})
    a * SU2(b)
end

"""
    *(a::Array{U1,2}, b::SU2)

Multiply two numbers.
"""
function *(a::Array{U1,2}, b::SU2)
    SU2(a) * b
end

"""
    *(a::SU2, b::Array{Complex,2})

Multiply two numbers.
"""
function *(a::SU2, b::Array{Complex,2})
    Matrix(a) * b
end

"""
    *(a::Array{Complex,2}, b::SU2)

Multiply two numbers.
"""
function *(a::Array{Complex,2}, b::SU2)
    s * Matrix(b)
end

"""
    *(a::SU2, b::U1)

Multiply two numbers.
"""
function *(a::SU2, b::U1)
    SU2(Matrix(a) .* b)
end

"""
    *(a::U1, b::SU2)

Multiply two numbers.
"""
function *(a::U1, b::SU2)
    SU2(a .* Matrix(b))
end


text:
const GLContext = Any
# struct GLContext{T} <: ShaderAbstractions.AbstractContext
#     context::T
#     opengl_version::VersionNumber
#     glsl_version::VersionNumber
#     # Use a unique id, since we can't track this via pointer identity
#     # (OpenGL may reuse the same pointers)
#     unique_id::UInt64
# end
#
# let counter = Threads.Atomic{UInt64}(0)
#     global unique_context_counter
#     function unique_context_counter()
#         # dont start at zero, so we can keep zero special
#         counter[] = counter[] + 1
#         return counter[]
#     end
# end
#
# function GLContext(::Nothing)
#     return GLContext(nothing, v"0.0.0", v"0.0.0", UInt64(0))
# end


text:
include("../src/BlueMesh7.jl")
using .BlueMesh7: generate_positions, initialize_mesh, start, plotgrid

function run()
    positions = generate_positions()

    node_roles = rand(0:1, length(positions))
    mesh = initialize_mesh(positions, node_roles)

    println("running the simulation for 1 minute (60000 steps, 1 step = 1ms)...")
    produced, received = start(mesh, minutes = 1)
    println("packets produced: $produced, received: $received, $(round(100 * received / produced; digits=2))%")

    plotgrid(mesh)
end

run()


text:
export AbstractHomotopy

"""
    AbstractHomotopy

An abstract type representing a homotopy ``H(x,t)``.

The following homotopies are available:

* [`AffineChartHomotopy`](@ref)
* [`AffineSubspaceHomotopy`](@ref)
* [`ModelKitHomotopy`](@ref)
* [`ParameterHomotopy`](@ref)
* [`StraightLineHomotopy`](@ref)
"""
abstract type AbstractHomotopy end

Base.size(H::AbstractHomotopy, i::Integer) = size(H)[i]
function set_solution!(x::AbstractVector, H::AbstractHomotopy, y::AbstractVector, t)
    x .= y
end
get_solution(H::AbstractHomotopy, x::AbstractVector, t) = copy(x)

start_parameters!(H::AbstractHomotopy, p) = H
target_parameters!(H::AbstractHomotopy, p) = H
# internal only
include("homotopies/toric_homotopy.jl")

# public, these should be linked on the top
include("homotopies/affine_chart_homotopy.jl")
include("homotopies/affine_subspace_homotopy.jl")
include("homotopies/model_kit_homotopy.jl")
include("homotopies/parameter_homotopy.jl")
include("homotopies/coefficient_homotopy.jl")
include("homotopies/straight_line_homotopy.jl")


text:
filename() = "day03_sample.txt"

# Read file. Repeat horizontally certain number of times.
lines(repeat) = [r^repeat for r in readlines(filename())]

# Part1
# Walk down the grid according to "right 3, 1 down" rule.
# Count how many trees (`#` sign) that I encounter.
# 32 times because 323 * 3 / 31 = 31.26
part1() = count([val[(row-1)*3+1] == '#' for (row,val) in enumerate(lines(32))])

# Part2
# Customized rule with `right` and `down` values.
function count_trees(lines, right, down)
    col = 1
    cnt = 0
    for row in 1:down:length(lines)
        row > 1 || continue
        col += right
        if lines[row][col] == '#'
            cnt += 1
        end
    end
    return cnt
end

function part2()
    grid = lines(73)
    counts = [count_trees(grid, x[1], x[2]) for x in ((1,1),(3,1),(5,1),(7,1),(1,2))]
    return prod(counts)
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    for f in (detect_rdata, detect_rdata_single, detectwav, detect_bedgraph,
              detecttiff, detect_noometiff, detect_ometiff, detectavi,
              detecthdf5, detect_stlascii, detect_stlbinary, detect_gadget2)
        @assert precompile(f, (IOStream,))
    end

    for F in (String, IOStream, Formatted)
        @assert precompile(query, (F,))
        @assert precompile(load, (F,))
        @assert precompile(save, (F,Nothing,))
        @assert precompile(loadstreaming, (F,))
        @assert precompile(savestreaming, (F,))
    end
    @assert precompile(action, (Symbol,Vector{Union{PkgId, Module}},Symbol,String))
    @assert precompile(action, (Symbol,Vector{Union{PkgId, Module}},Symbol,IOStream))
    @assert precompile(action, (Symbol,Vector{Union{PkgId, Module}},Formatted))
    @assert precompile(loadstreaming, (Function, Any))
    @assert precompile(savestreaming, (Function, Any))
    @assert precompile(skipmagic, (IOStream,Vector{Vector{UInt8}},))
end


text:
doc"""
  generate_makefile(io, defs)

  generates a makefile which builds the make file corresponding to the posit definitions.
"""
function generate_makefile(io, posit_defs)

  #a list of general posit file names.
  pfilenames = [:posit_conv, :posit_extended, :posit_basics, :posit_err, :posit_full]

  #a list of posit classnames
  pclassnames = []
  for posit_size in sort(collect(keys(posit_defs)))
    for exp in posit_defs[posit_size]
      push!(pclassnames, Symbol("P$(posit_size)e$(exp)"))
    end
  end

  obj_file = (s::Union{Symbol}) -> ("obj/$s.o")
  cpp_file = (s::Union{Symbol}) -> ("$s.cpp")

  obj_def = (s::Symbol) -> """
  $(obj_file(s)): $(cpp_file(s))
  \tgcc -c -Wall -Werror -fpic $(cpp_file(s)) -o $(obj_file(s))
  """

  all_objects = join(map(obj_file, pfilenames), " ") * " " * join(map(obj_file, pclassnames), " ")

  c_objects = join(map(obj_def, pfilenames),"\n")

  cpp_objects = join(map(obj_def, pclassnames), "\n")

  link_library = """
  build: $(all_objects)
  \tgcc -shared -o libfastposit.so $all_objects
  """

  write(io, """
  $(c_objects)

  $(cpp_objects)

  $(link_library)

  test: test.cpp
  \tgcc -Wall -o test test.cpp -lstdc++ -lfastposit
  \tchmod +x test
  \t./test

  .PHONY: install
  install:
  \tcp ./libfastposit.so /usr/lib
  \tchmod +r /usr/lib/libfastposit.so

  .PHONY: clean
  clean:
  \t-rm *.so
  \t-rm obj/*.o
  """)
end


text:
# RanDAGen/src/graphbits.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

abstract type AbstractGraph end

mutable struct Vertex{T<:Any}
    ID::Int64
    data::T

    neighbors::Vector{Int64}

    function Vertex{T}(
        ID::Int64=rand(Int64, 1),
        data::T=undef) where {T<:Any}
        return new(ID, data)
    end
end

mutable struct Edge{T<:Any}
    src::Int64
    data::T
    dest::Int64

    function Edge{T}(
        src::Int64=-1,
        data::T=undef,
        dest::Int64=nothing) where {T<:Any}
        return new(src, data, dest)
    end
end

mutable struct Graph{T<:Any} <: AbstractGraph

    vertices::Dict{Int64,Vertex{T}}

    edges::Dict{(Int6,Int64),Edge{T}}

    function Graph{T}(
        vertices=Vector{Vertex{T}},
        edges=Vector{Edge{T}}) where {T<:Any}
        return new(vertices, edges)
    end
end

nv(graph::Graph) = length(graph.vertices)
ne(graph::Graph) = length(graph.vertices)

vertices(graph::Graph) = [v for v in values(graph.vertices)]
edges(graph::Graph) = [v for v in values(graph.edges)]

has_vertex(graph::Grap, v::Int64) = haskey(graph.vertices, v)
has_edge(graph::Graph, e::(Int64,Int64)) = haskey(graph.edges, e)

function add_edge!(graph::Graph, edge::Edge) begin
    if has_vertex(graph, edge.src) && has_vertex(graph, edge.dest)

        # Add neighbors
        union!(graph.vertices[edge.src].neighbors, edge.dest)
        union!(graph.vertices[edge.dest].neighbors, edge.src)
        graph.edges[(edge.src,edge.dest)] = edge

        return true
    end
    return false
end

function delete_edge!(graph::Graph, edge::Tuple{Int64,Int64}) begin
    delete!(graph.edges, (edge.first,edge.second))
    delete!(graph.edges, (edge.second,edge.first))
end

function add_node!(graph::Graph, vertex::Vertex) begin
    if !has_vertex(graph, vertex.ID)
        graph.vertices[vertex.ID] = vertex

        for neighbor in vertex.neighbors
            edge = Edge{Nothing}(1,nothing,2)
            add_edge!(graph, edge)
        end

        return true
    end
    return false
end

function add_node!(graph::Graph, vertex::Vertex, edges::Array{Edge}) begin
    if !has_vertex(graph, vertex.ID)
        graph.vertices[vertex.ID] = vertex

        for edge in edges
            add_edge!(graph, edge)
        end

        return true
    end
    return false
end

function delete_node!(graph::Graph, ID::Int64) begin
    if !has_vertex(graph, ID)
        neighborhood = graph.vertices[ID].neighbors

        for neighbor in neighborhood
            delete_edge!(graph.edges, (ID,neighbor.ID))
            filter!(x->ID,graph.vertices[neighbor].neighbors)
        end

        delete!(graph.vertices, ID)

        return true
    end
    return false
end


text:
module Menu

    using ArgParse

    function main_menu()
        s = ArgParseSettings()

        @add_arg_table! s begin
            "--start", "--year-start"
                arg_type = Int
                help = "Year start"
                required = true

            "--end", "--year-end"
                arg_type = Int
                help = "Year end"
                required = true

            "--width", "--area-width"
                arg_type = Float32
                help = "Area width"
                required = true

            "--height", "--area-height"
                arg_type = Float32
                help = "Area height"
                required = true

            "--timeout"
                arg_type = Int
                default = 300  # 5 minutes
                help = "Timeout"

        end

        args = parse_args(s)

        # check the arguments correctness
        @assert args["start"] <= args["end"]

        return args
    end

end


text:
function ModifiedLambdaMethod()
    return ModifiedLambdaMethod(())
end


text: trace_norm(A::Operator) = sum(d -> √(abs(d)) , eigvals(A'*A))


text:
@testset "Save particles data into file" begin
    n_particles = 100000
    xmin = 1.0::Float64
    xmax = 4π + 1.0
    Lx = xmax - xmin
    nx = 64

    mesh = OneDGrid(xmax, xmin, nx)

    pg = ParticleGroup{1,2}(n_particles)

    params = (k=[[0.5]], α=[0.01], σ=[[0.1, 2.0]], μ=[[0.0, 0.0]])

    df = CosSumGaussian{1,2}(params...)

    n_particles = pg.n_particles

    sampler = ParticleSampler{1,2}(:sobol, true, n_particles)

    sample!(pg, sampler, df, mesh)

    GEMPIC.save("test", 1, pg)

    @test isfile("test-000001.jld2")

    σ, μ = 0.02, 0.0
    kx, α = 1.004355, 0.001
    xmin, xmax = 0, 2π / kx
    domain = [xmin, xmax, xmax - xmin]
    nx = 64
    n_particles = 1000
    mesh = OneDGrid(xmin, xmax, nx)
    spline_degree = 3

    df = CosSumGaussianSpin([[kx]], [α], [[σ]], [[μ]])

    spg = ParticleGroup{1,1}(n_particles; n_spin=3)
    sampler = ParticleSampler{1,1}(:sobol, n_particles)

    sample!(spg, sampler, df, mesh)

    GEMPIC.save("test", 2, spg)

    @test isfile("test-000002.jld2")
end


text:
using StochasticDiffEq, Test
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -9.81
end

function g(du,u,p,t)
  nothing
end

function condtion(u,t,integrator) # Event when event_f(u,p,t,k) == 0
  u[1]
end

affect! = nothing
function affect_neg!(integrator)
  integrator.u[2] = -integrator.u[2]
end

# Continuous callback
callback = ContinuousCallback(condtion,affect!,affect_neg!)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = SDEProblem(f,g,u0,tspan)

sol = solve(prob,SRIW1(),callback=callback,adaptive=false,dt=3/4)

@test minimum(sol[1,:]) > -1e-12 && minimum(sol[1,:]) < 1e-12

sol = solve(prob,SRIW1(),callback=callback,save_everystep=false)
t = sol.t[end÷2] # this is the callback time point
sol = solve(prob,SRIW1(),callback=callback,saveat=t)
@test count(x->x==t, sol.t) == 2
sol = solve(prob,SRIW1(),callback=callback,saveat=t-eps(t))
@test count(x->x==t, sol.t) == 2

function g(du,u,p,t)
  du[2] = .125*u[2]
end

prob = SDEProblem(f,g,u0,tspan)

sol = solve(prob,SRIW1(),callback=callback)

sol = solve(prob,EM(),callback=callback,dt=1/4)

# Discrete callback
tstop = [5.;8.]
condition_dc = (u,t,integrator) -> t in tstop
affect!_dc = (integrator) -> integrator.u .= 1.0
save_positions = (true,true)
callback_dc = DiscreteCallback(condition_dc, affect!_dc, save_positions=save_positions)
sol = solve(prob, SRIW1(), callback=callback_dc, tstops=tstop, saveat=tstop)
@test count(x->x==tstop[1], sol.t) == 2
@test count(x->x==tstop[2], sol.t) == 2
sol = solve(prob, SRIW1(), callback=callback_dc, tstops=tstop, saveat=prevfloat.(tstop))
@test count(x->x==tstop[1], sol.t) == 2
@test count(x->x==tstop[2], sol.t) == 2


text:
# Initial Idea, far from usability
# DO NOT use it
struct LineBuffer <: AbstractString
    data::Vector{UInt8}
end

function Base.iterate(s::LineBuffer, i::Int = 1)
    i > length(s.data) && return nothing
    return (Char(s.data[i]), i+1)
end

Base.lastindex(s::LineBuffer) = length(s.data)

Base.getindex(s::LineBuffer, i::Int) = s.data[i]

Base.sizeof(s::LineBuffer) = sizeof(s.data)

Base.length(s::LineBuffer) = length(s.data)

Base.ncodeunits(s::LineBuffer) = length(s.data)

Base.codeunit(::LineBuffer) = UInt8
Base.codeunit(s::LineBuffer, i::Integer) = s.data[i]

Base.isvalid(s::LineBuffer, i::Int) = checkbounds(Bool, s, i)

# minimum type defination for parsing TimeType data

struct DT{N} <: AbstractString
    data::Vector{UInt8}
end

strlength(s::DT{N}) where N = N
function allocateDT(::Type{DT{N}}, len) where N
    DT{N}(fill(0x00, N*len))
end


text:
#On beaglebone, run:
# include("LabConnections/src/LabConnections.jl")
# using LabConnections.BeagleBone
# run_server()

using LabConnections.Computer

stream = BeagleBoneStream(ip"192.168.7.2")
led1 = SysLED(1)
led2 = SysLED(2)
led3 = SysLED(3)
led4 = SysLED(4)

# Send info to steam about which inputs/otputs to initialize
# and adds a ref in motor and led to stream
init_devices!(stream, led1, led2, led3, led4)
ledon = true
for i = 1:100
    put!(led1, ledon)
    put!(led2, !ledon)
    put!(led3, ledon)
    put!(led4, !ledon)
    send(stream) #Sends all the outputs to the stream in one atomic call
    #read(stream)
    get(led1)
    get(led2)
    get(led3)
    #sleep(0.1)
    v1,v2,v3 = read(stream) #Sends request to read, reads all inputs for which get! was called
    v1 == !v2 == v3 == ledon ? nothing : println("ledon is $ledon, but read v1, v2, v3 = $v1, $v2, $v3")
    ledon = !ledon
end
for i = 1:40
    send(led1, ledon)
    #sleep(0.03)
    v1 = read(led1)
    send(led2, ledon)
    #sleep(0.03)
    v2 = read(led2)
    send(led3, ledon)
    #sleep(0.03)
    v3 = read(led3)
    send(led4, ledon)
    #sleep(0.03)
    v4 = read(led4)
    v1 == v2 == v3 == v4 == ledon ? nothing : println("ledon is $ledon, but read v1, v2, v3, v4 = $v1, $v2, $v3, $v4")
    ledon = !ledon
end
put!(led1, false)
put!(led2, false)
put!(led3, false)
put!(led4, false)
send(stream)
close(stream) #Tells BeagleBone to stop listening and close outputs


text:
function f()
  println("Hello")
end


text:
requiredTypes = [ "AbstractLevel" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include(  reqType * ".jl" )
    end  # if !isdefined( Symbol( ...
end

export SubJob

type SubJob <: AbstractLevel

    # Name of the SubJob ex: X, A, B
    Name::String

    # Maximum and minimum stay in this SubJob
    MinStay::Int
    MaxStay::Int

    # A list of allowed PS. Just to maintain the same structure (no next level needed)
    NextLevel::Array{String}

    SubJob() = new("", 0, 0, Array{String}())
    SubJob(nm::String;
           min::Int = 0,
           max::Int = 0,
           next::Array{String} = Array{String}()) = new(nm, min, max, next)

end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
struct VPIntegratorParameters{T}
    dt::T          # time step
    nₜ::Int        # number of time steps
    nₛ::Int        # number of saved time steps
    nₕ::Int        # number of basis functions
    nₚ::Int        # number of particles

    t::Vector{T}

    function VPIntegratorParameters(dt::T, nₜ::Int, nₛ::Int, nₕ::Int, nₚ::Int) where {T}
        t = collect(range(0, stop=dt*nₜ, length=nₛ))
        new{T}(dt,nₜ,nₛ,nₕ,nₚ,t)
    end
end

mutable struct VPIntegratorCache{T}
    x::Vector{T}
    v::Vector{T}
    a::Vector{T}
    w::Vector{T}

    ρ::Vector{T}
    ϕ::Vector{T}
    rhs::Vector{T}

    X::Matrix{T}
    V::Matrix{T}
    A::Matrix{T}
    Φ::Matrix{T}

    W::Vector{T}
    K::Vector{T}
    M::Vector{T}
end

VPIntegratorCache(IP::VPIntegratorParameters{T}) where {T} = VPIntegratorCache(
                                                                zeros(T,IP.nₚ), # x
                                                                zeros(T,IP.nₚ), # v
                                                                zeros(T,IP.nₚ), # a
                                                                zeros(T,IP.nₚ), # w
                                                                zeros(T,IP.nₕ), # ρ
                                                                zeros(T,IP.nₕ), # ϕ
                                                                zeros(T,IP.nₕ), # rhs
                                                                zeros(T,IP.nₚ,IP.nₛ), # X
                                                                zeros(T,IP.nₚ,IP.nₛ), # V
                                                                zeros(T,IP.nₚ,IP.nₛ), # A
                                                                zeros(T,IP.nₕ,IP.nₛ), # Φ
                                                                zeros(T,IP.nₛ), # W
                                                                zeros(T,IP.nₛ), # K
                                                                zeros(T,IP.nₛ)  # M
                                                            )

function save_timestep!(IC::VPIntegratorCache, poisson::PoissonSolver, ts, χ)
    IC.X[:,ts] .= IC.x
    IC.V[:,ts] .= IC.v
    IC.A[:,ts] .= IC.a
    IC.Φ[:,ts] .= poisson.ϕ

    IC.W[ts] = dot(poisson.ϕ, poisson.S, poisson.ϕ) / 2 * χ^2
    IC.K[ts] = dot(IC.w .* IC.v, IC.v) / 2
    IC.M[ts] = dot(IC.w, IC.v)
end

function solve_potential!(IC::VPIntegratorCache, poisson::PoissonSolver, ts, χ, Φₑₓₜ, given_phi)
    if given_phi
        IC.ϕ .= Φₑₓₜ[:,ts]
    else
        solve!(poisson, IC.x, IC.w)
        poisson.ϕ ./= χ^2
        IC.ϕ .= poisson.ϕ
    end
end

function integrate_vp!(P::ParticleList{T},
                       poisson::PoissonSolver{T},
                       parameters::NamedTuple,
                       IP::VPIntegratorParameters{T},
                       IC::VPIntegratorCache{T} = VPIntegratorCache(IP);
                       given_phi = false,
                       Φₑₓₜ::Array{T} = zeros(T, IP.nₕ, IP.nₚ),
                       save = true) where {T}

    nsave = div(IP.nₜ, IP.nₛ-1)

    if given_phi
        @assert IP.nₛ == IP.nₜ + 1
    end

    # simulation parameter
    χ = parameters.χ

    # initial conditions
    IC.x .= P.x
    IC.v .= P.v
    IC.w .= P.w

    # save initial conditions
    if save
        solve_potential!(IC, poisson, 1, χ, Φₑₓₜ, given_phi)
        save_timestep!(IC, poisson, 1, χ)
    end

    tₛ = 1
    for t in 1:IP.nₜ
        # half an advection step
        IC.x .+= 0.5 .* IP.dt .* IC.v .* χ

        # solve for potential
        solve_potential!(IC, poisson, t+1, χ, Φₑₓₜ, given_phi)

        # evaluate electric field
        eval_field!(IC.a, poisson, IC.x)

        # acceleration step
        IC.v .+= IP.dt .* IC.a .* χ

        # half an advection step
        IC.x .+= 0.5 .* IP.dt .* IC.v .* χ

        if save && t % nsave == 0
            solve_potential!(IC, poisson, t+1, χ, Φₑₓₜ, given_phi)
            save_timestep!(IC, poisson, tₛ+1, χ)
            tₛ += 1
        end
    end

    # P.x .= IC.x
    # P.v .= IC.v
end


text:
using Oceananigans.Grids: total_length, topology

using OffsetArrays: OffsetArray

#####
##### Creating offset arrays for field data by dispatching on architecture.
#####

"""
Return a range of indices for a field located at `Center` centers
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(loc, topo, N, H=0) = 1 - H : N + H

"""
Return a range of indices for a field located at cell `Face`s
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(::Type{Face}, ::Type{Bounded}, N, H=0) = 1 - H : N + H + 1

"""
Return a range of indices for a field along a 'reduced' dimension.
"""
offset_indices(::Type{Nothing}, topo, N, H=0) = 1 : 1

function offset_data(underlying_data::AbstractArray{FT, 3}, loc, topo, N, H) where FT
    ii = offset_indices(loc[1], topo[1], N[1], H[1])
    jj = offset_indices(loc[2], topo[2], N[2], H[2])
    kk = offset_indices(loc[3], topo[3], N[3], H[3])

    return OffsetArray(underlying_data, ii, jj, kk)
end

function offset_data(underlying_data::AbstractArray{FT, 2}, loc, topo, N, H) where FT
    ii = offset_indices(loc[1], topo[1], N[1], H[1])
    jj = offset_indices(loc[2], topo[2], N[2], H[2])

    return OffsetArray(underlying_data, ii, jj)
end

"""
    offset_data(underlying_data, grid::AbstractGrid, loc)

Returns an `OffsetArray` that maps to `underlying_data` in memory,
with offset indices appropriate for the `data` of a field on
a `grid` of `size(grid)` and located at `loc`.
"""
offset_data(underlying_data, grid::AbstractGrid, loc) =
    offset_data(underlying_data, loc, topology(grid), size(grid), halo_size(grid))

"""
    new_data([FT=Float64], arch, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT` on `arch`itecture,
with indices corresponding to a field on a `grid` of `size(grid)` and located at `loc`.
"""
function new_data(FT, arch, grid, loc)
    underlying_data = zeros(FT, arch,
                            total_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                            total_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                            total_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    return offset_data(underlying_data, grid, loc)
end

# Default to type of Grid
new_data(arch, grid, loc) = new_data(eltype(grid), arch, grid, loc)


text:
"""
    combinedims(array_of_arrays)

Combine the dimensions of a nested array structure into a new, flat array.

This is the inverse operation of `splitdims` / `splitdimsview`.

See also `combinedimsview`, the lazy version of this function.

# Example

```julia
julia> combinedims([[1, 2], [3, 4]])
2×2 Array{Int64,2}:
 1  3
 2  4
```
"""
@inline combinedims(a::AbstractArray{<:AbstractArray{<:Any, N}, M}) where {N,M} = combinedims(a, ntuple(i -> N + i, M))
@inline combinedims(a::AbstractArray{<:AbstractArray{<:Any, N}, 1}, outerdim::Int) where {N} = _combinedims(a, Val((outerdim,)))
@inline combinedims(a::AbstractArray{<:AbstractArray{<:Any, N}, M}, outerdims::NTuple{M, Int}) where {N, M} = _combinedims(a, Val(outerdims))

function _combinedims(a::AbstractArray, ::Val{outer_dims}) where {outer_dims}
    outeraxes = axes(a)
    inneraxes = _inneraxes(a)
    ndims_total = length(outeraxes) + length(inneraxes)
    newaxes = _combine_tuples(ndims_total, outer_dims, outeraxes, inneraxes)

    T = inner_eltype(a)
    out = similar(a, T, newaxes)
    for j in CartesianIndices(outeraxes)
        I = slice_inds(j, Val(outer_dims), Val(ndims_total))
        out[I..., :] = a[j]  # trailing ':' is required for zero-dimensional a[j]
    end
    return out
end

@inline function _combine_tuples(n, dims, t_then, t_else)
    i_out = Ref(0)
    i_in = Ref(0)
    ntuple(n) do i
        i ∈ dims ? t_then[i_out[] += 1] : t_else[i_in[] += 1]
    end
end

_inneraxes(a) = axes(first(a)) # Can specialize this for static arrays, for example

inner_eltype(a) = Any
inner_eltype(a::AbstractArray{<:AbstractArray{T}}) where {T} = T

## Lazy version

struct CombineDimsArray{T, N, Nout, A} <: AbstractArray{T, N}
    parent::A
    outer_dims::NTuple{Nout, Int}
end

Base.parent(a::CombineDimsArray) = a.parent

@inline dims_outer(a::CombineDimsArray) = a.outer_dims
@inline function dims_inner(a::CombineDimsArray{T, N, Nout}) where {T, N, Nout}
    # below is a type-stable version of:
    # filter(∉(dims_outer(a)), ntuple(identity, N))
    i = Ref(0)
    ntuple(N - Nout) do _
        i[] = findnext(∉(dims_outer(a)), 1:N, i[] + 1)
    end
end

Base.size(a::CombineDimsArray) = _combine_tuples(ndims(a), dims_outer(a), size(parent(a)), size(first(parent(a))))
axes(a::CombineDimsArray) = _combine_tuples(ndims(a), dims_outer(a), axes(parent(a)), axes(first(parent(a))))
Base.IndexStyle(::CombineDimsArray) = Base.IndexCartesian()
@propagate_inbounds function Base.getindex(a::CombineDimsArray{T, N}, i::Vararg{Int, N}) where {T, N}
    outer_inds = getindices(i, dims_outer(a))
    inner_inds = getindices(i, dims_inner(a))
    return parent(a)[outer_inds...][inner_inds...]
end

"""
    combinedimsview(array_of_arrays)

Lazily create a flat array view of a nested array structure.

This is the inverse operation of `splitdims` / `splitdimsview`.

See also `combinedims`, the eager version of this function.
"""
combinedimsview(a::AbstractArray{<:AbstractArray{T, N}, M}) where {T, N, M} = combinedimsview(a, ntuple(i -> N + i, Val(M)))
combinedimsview(a::AbstractArray{<:AbstractArray{T, N}, 1}, outerdim::Int) where {T, N} = combinedimsview(a, (outerdim,))
combinedimsview(a::AbstractArray{<:AbstractArray{T, N}, M}, outerdims::NTuple{M, Int}) where {T, N, M} = CombineDimsArray{inner_eltype(a), N + M, M, typeof(a)}(a, outerdims)


text:
#####################################################
# Surya
# NREL
# November 2021
# NAERM Threats
# Testing of EGRET2SIIP module
#######################################################
# Loading the required packages
#######################################################
# Load EGRET2SIIP module here
#####################################################################################
# Read EGRET System JSON
#####################################################################################
import JSON
Location_1 = dirname(pwd())*"/test/DAY_AHEAD_Model_2020-07-01_2020-07-14.json";
Location_2 = dirname(pwd())*"/test/REAL_TIME_Model_2020-07-01_2020-07-14.json";
EGRET_json_DA = JSON.parsefile(Location_1);
EGRET_json_RT = JSON.parsefile(Location_2);
#####################################################################################
# Convert EGRET JSON to CSV in a format SIIP Tabular Data understands
#####################################################################################
location, base_MVA,rt_flag = EGRET2SIIP.parse_EGRET_JSON(EGRET_json_DA,EGRET_json_RT=EGRET_json_RT);

if (rt_flag)
    sys_DA,sys_RT = EGRET2SIIP.parse_tabular_data(location,base_MVA,rt_flag,ts_pointers_file ="CSV");
else
    sys_DA = EGRET2SIIP.parse_tabular_data(location,base_MVA,rt_flag,ts_pointers_file="CSV");
end

#           (or)
sys_DA, sys_RT = EGRET2SIIP.EGRET_TO_PSY(EGRET_json_DA,EGRET_json_RT=EGRET_json_RT,ts_pointers_file="CSV");


text:
@doc raw"""
```
DeterministicCustom(p::AbstractArray, q::AbstractArray, integration_factor::Real = 10, smoothness::Real = 20, β::Real = 1.0, integration_scheme::AbstractIntegrator = Horn())
```
Same as [`Deterministic`](@ref) exept impliments a custom gradient descent and
line search algorithm.
# Example
The following example demonstrates the use of the `DeterministicCustom` method.
```julia
using ShapeFromShading, Makie, Images

#generate synthetic image and initial solution
img = generate_surface(SynthSphere(38))
pin, qin = synthetic_gradient(SynthSphere(38))

#calculate the heightmap
deterministicCustom = DeterministicCustom(p=pin,q=qin, integration_factor = 10, smoothness = 50, β = 1.0, integration_scheme = Horn())
Z,p,q = deterministicCustom(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(abs.(Z))

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)
```
# Reference
[1] A. Crouzil, X. Descombes and J. Durou, "A multiresolution approach for shape from shading coupling deterministic and stochastic optimization," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 11, pp. 1416-1421, Nov. 2003. [doi: 10.1109/TPAMI.2003.1240116](https://doi.org/10.1109/TPAMI.2003.1240116 )
"""
function (algorithm::DeterministicCustom)(img::AbstractArray)
    # Setup initial variables
    E = Array{Float64}(img)
    Eₘ = maximum(E)
    p = copy(algorithm.p)
    q = copy(algorithm.q)
    N, M = size(E)
    λᵢ = algorithm.integration_factor
    λₛ = algorithm.smoothness
    β = algorithm.β
    R = zeros(Float64,axes(E))
    Z = zeros(axes(E))
    Δϵ = zeros(Float64,(size(E)...,2))
    val = sqrt(2 * length(E))*2 + 10
    prevd = 0.0001

    # Loop untill norm of Δϵ is less then desired value
    while val > β*sqrt(2 * length(E))

        # Calculate the energy
        ϵ = 0.0
        for i = 1:(M-1)
            for j = 1:(N-1)
                ϵ += (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i])^2 + λᵢ*((p[j,i+1]-p[j,i])-(q[j+1,i]-q[j,i]))^2 + λₛ*((p[j,i+1]-p[j,i])^2 + (p[j+1,i]-p[j,i])^2 + (q[j,i+1]-q[j,i])^2 + (q[j+1,i]-q[j,i])^2)
                if i>1 && j>1
                    Δϵ[j,i,1] = -2 * Eₘ * p[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) - 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((p[j,i+1]-p[j,i]) + (p[j+1,i]-p[j,i])) + 2*λᵢ*((p[j,i]-p[j-1,i])-(q[j-1,i+1]-q[j-1,i])) + 2*λₛ*((p[j,i]-p[j,i-1]) + (p[j,i]-p[j-1,i]))
                    Δϵ[j,i,2] = -2 * Eₘ * q[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) + 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((q[j,i+1]-q[j,i]) + (q[j+1,i]-q[j,i])) - 2*λᵢ*((p[j+1,i-1]-p[j,i-1])-(q[j,i]-q[j,i-1])) + 2*λₛ*((q[j,i]-q[j,i-1]) + (q[j,i]-q[j-1,i]))
                elseif i>1
                    Δϵ[j,i,1] = -2 * Eₘ * p[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) - 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((p[j,i+1]-p[j,i]) + (p[j+1,i]-p[j,i])) + 2*λₛ*(p[j,i]-p[j,i-1])
                    Δϵ[j,i,2] = -2 * Eₘ * q[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) + 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((q[j,i+1]-q[j,i]) + (q[j+1,i]-q[j,i])) - 2*λᵢ*((p[j+1,i-1]-p[j,i-1])-(q[j,i]-q[j,i-1])) + 2*λₛ*(q[j,i]-q[j,i-1])
                elseif j>1
                    Δϵ[j,i,1] = -2 * Eₘ * p[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) - 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((p[j,i+1]-p[j,i]) + (p[j+1,i]-p[j,i])) + 2*λᵢ*((p[j,i]-p[j-1,i])-(q[j-1,i+1]-q[j-1,i])) +2*λₛ*(p[j,i]-p[j-1,i])
                    Δϵ[j,i,2] = -2 * Eₘ * q[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) + 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((q[j,i+1]-q[j,i]) + (q[j+1,i]-q[j,i])) + 2*λᵢ*(q[j,i]-q[j-1,i])
                else
                    Δϵ[j,i,1] = -2 * Eₘ * p[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) - 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((p[j,i+1]-p[j,i]) + (p[j+1,i]-p[j,i]))
                    Δϵ[j,i,2] = -2 * Eₘ * q[j,i] * (Eₘ/sqrt(1+p[j,i]^2+q[j,i]^2)-E[j,i]) * (1/(1+p[j,i]^2+q[j,i]^2)^(3/2)) + 2*λᵢ*((p[j+1,i]-p[j,i])-(q[j,i+1]-q[j,i])) - 2*λₛ*((q[j,i+1]-q[j,i]) + (q[j+1,i]-q[j,i]))
                end
            end
        end

        # Update norm value
        val = norm(Δϵ)

        # Line search to find minimizor of ϵ along the direction of gradient
        dmin = 1
        ϵmin = typemax(Float64)
        for d = prevd/10:prevd/10:prevd*10
            # Setup new gradient field
            ϵtemp = 0.0
            ptemp = copy(p)
            qtemp = copy(q)
            for i = 1:(M-1)
                for j = 1:(N-1)
                    ptemp[j,i] = p[j,i] - d*Δϵ[j,i,1]
                    qtemp[j,i] = q[j,i] - d*Δϵ[j,i,2]
                end
            end

            # Calculate new energy
            for i = 1:(M-1)
                for j = 1:(N-1)
                    ϵtemp += (Eₘ/sqrt(1+(ptemp[j,i])^2+qtemp[j,i]^2)-E[j,i])^2 + λᵢ*((ptemp[j,i+1]-ptemp[j,i])-(qtemp[j+1,i]-qtemp[j,i]))^2 + λₛ*((ptemp[j,i+1]-ptemp[j,i])^2 + (ptemp[j+1,i]-ptemp[j,i])^2 + (qtemp[j,i+1]-qtemp[j,i])^2 + (qtemp[j+1,i]-qtemp[j,i])^2)
                end
            end

            # Check if new energy is less then the minimum found
            if ϵtemp < ϵmin
                ϵmin = ϵtemp
                dmin = d
            end
        end

        # Set gradients to minimize energy
        for i = 1:(M-1)
            for j = 1:(N-1)
                p[j,i] = p[j,i] - dmin*Δϵ[j,i,1]
                q[j,i] = q[j,i] - dmin*Δϵ[j,i,2]
            end
        end
        prevd = dmin
        @show maximum(Δϵ)
        @show norm(Δϵ), ϵ, dmin
    end
    Z = algorithm.integration_scheme(p,q)
    return Z,p,q
end

# Energy function
function f(E, pq, λᵢ, λₛ)
    Eₘ = maximum(E)
    val = 0.0
    N, M = size(E)
    for j = 1:(M-1)
        for i = 1:(N-1)
            val += (Eₘ/sqrt(1+pq[j,i,1]^2+pq[j,i,2]^2)-E[j,i])^2 + λᵢ*((pq[j,i+1,1]-pq[j,i,1])-(pq[j+1,i,2]-pq[j,i,2]))^2 + λₛ*((pq[j,i+1,1]-pq[j,i,1])^2 + (pq[j+1,i,1]-pq[j,i,1])^2 + (pq[j,i+1,2]-pq[j,i,2])^2 + (pq[j+1,i,2]-pq[j,i,2])^2)
        end
    end
    return val
end

@doc raw"""
```
Deterministic(p::AbstractArray, q::AbstractArray, integration_factor::Real = 10, smoothness::Real = 20, β::Real = 1.0, integration_scheme::AbstractIntegrator = Horn())
```
Impliments a deterministic method for reconstructing the surface using a gradient
descent method. This version makes use of a [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl )
solver.
# Output
Returns a `Deterministic` functor which can be called on an image to reconstruct
the surface normals from the image.
# Details
Under the assumption that the object is illuminated from directly above ([0,0,1])
and that the surface is Lambertian, the eikonal equation, shown below, can be used
to describe the surface in an image.
```math
E_{i,j}=\frac{E_{max}}{\sqrt{p_{i,j}^2+q_{i,j}^2+1}}
```
As the surface must be integrable, the integrability condition; ``p_y=q_x`` must
hold on the surface. Using this a discreate itegrability function can be derived
in the form:
```math
p_{i,j+1}-p_{i,j}=q_{i+1,j}-q_{i,j}
```
A discreate functional for the energy can now be derived from these in the form:
```math
\begin{aligned}
\epsilon&=\sum_{i,j\in\Omega}\left[\frac{E_{max}}{\sqrt{p_{i,j}^2+q_{i,j}^2+1}}-E_{i,j}\right]^2\\
&+\lambda_i\sum_{i,j\in\Omega}[(p_{i,j+1}-p_{i,j})-(q_{i+1,j}-q_{i,j})]^2
\end{aligned}
```
Finally a smoothness term can be introduced as below. This both aids in the quality
of the reconstruction and reduces the number of minima in the functional helping
to ensure a good solution is found.
```math
\begin{aligned}
\epsilon&=\sum_{i,j\in\Omega}\left[\frac{E_{max}}{\sqrt{p_{i,j}^2+q_{i,j}^2+1}}-E_{i,j}\right]^2\\
&+\lambda_i\sum_{i,j\in\Omega}[(p_{i,j+1}-p_{i,j})-(q_{i+1,j}-q_{i,j})]^2\\
&+\lambda_i\sum_{i,j\in\Omega}[(p_{i,j+1}-p_{i,j})^2+(p_{i+1,j}-p_{i,j})^2+(q_{i+1,j}-q_{i,j})^2+(q_{i,j+1}-q_{i,j})^2]
\end{aligned}
```
The problem is then minimized using gradient decsent method to attampt to find
a solution minimizing the energy function. Gradient descent is terminated when
``|\nabla\epsilon|<\beta\sqrt{2|\Omega|}``.
# Parameters
## `p`
An `AbstractArray` which contains initial solution for p, the x gradient of the
surface. Better inital solution will produce better results.
## `q`
An `AbstractArray` which contains initial solution for q, the y gradient of the
surface. Better inital solution will produce better results.
## `integration_factor`
A `Real` which controls the impact of the integration term on the energy function.
Can be used to tune the final solution. Defults to 10.
## `smoothness`
A `Real` which controls the impact of the smoothness term on the energy function.
Can be used to tune the final solution. Defults to 20.
## `β`
A `Real` which controls the value of ``|\nabla\epsilon|`` which will terminate
the algorithm. If algorithm terminates on first iteration value will need reducing
and if algorithm fails to converge the value will need increasing. Defults to 1.0.
## `integration_scheme`
A `AbstractIntegrator` which will be used to integrate the calculated gradients.
Defults to `Horn()`.
# Example
The following example demonstrates the use of the `Deterministic` method.
```julia
using ShapeFromShading, Makie, Images

#generate synthetic image and initial solution
img = generate_surface(SynthSphere(38))
pin, qin = synthetic_gradient(SynthSphere(38))

#calculate the heightmap
deterministic = Deterministic(p=pin,q=qin, integration_factor = 10, smoothness = 50, β = 1.0, integration_scheme = Horn())
Z,p,q = deterministic(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(abs.(Z))

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)
```
# Reference
[1] A. Crouzil, X. Descombes and J. Durou, "A multiresolution approach for shape from shading coupling deterministic and stochastic optimization," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 11, pp. 1416-1421, Nov. 2003. [doi: 10.1109/TPAMI.2003.1240116](https://doi.org/10.1109/TPAMI.2003.1240116 )
"""
function (algorithm::Deterministic)(img::AbstractArray)
    # Setup initial variables
    E = Float64.(img)
    x₀ = zeros(Float64,first(size(img)),last(size(img)),2)
    x₀[:,:,1] .= algorithm.p
    x₀[:,:,2] .= algorithm.q
    λᵢ = algorithm.integration_factor
    λₛ = algorithm.smoothness
    β = algorithm.β

    # Calcualte derivatives
    od = OnceDifferentiable(pq -> f(E, pq, λᵢ, λₛ), x₀; autodiff = :forward);
    a = zeros(size(x₀))

    # Run gradient descent
    res = optimize(od, x₀, GradientDescent(), Optim.Options(show_trace = true, g_tol = β*sqrt(2*length(E))))

    # Get minimised gradients
    out = Optim.minimizer(res)
    p = out[:,:,1]
    q = out[:,:,2]
    Z = algorithm.integration_scheme(p,q)
    return Z,p,q
end


text:
module FisherPlot

using LaTeXStrings
using CairoMakie
using Makie
using Base: @kwdef

function ellipseparameterization(a::Float64, b::Float64, θ::Float64)
    t = LinRange(0,2π, 200)
    x = Array(a .* sin.(θ) .* cos.(t) + b .* cos.(θ) .* sin.(t))
    y = Array(a .* cos.(θ) .* cos.(t) - b .* sin.(θ) .* sin.(t))
    return x, y
end

function gaussian(μ::Float64, σ::Float64, x)
    return 1/(sqrt(2π*σ^2))*exp(-0.5*(x-μ)^2/σ^2)
end

function ellipseparameters(covmatrix::Matrix{Float64}, i::Int64, j::Int64)
    σi = sqrt(covmatrix[i,i])
    σj = sqrt(covmatrix[j,j])
    σij = covmatrix[i,j]
    θ = (atan(2σij,(σi^2-σj^2)))/2
    a = sqrt((σi^2+σj^2)/2+sqrt(((σi^2-σj^2)^2)/4+σij^2))
    if i == j
        b = 0.
    else
        b = sqrt((σi^2+σj^2)/2-sqrt(((σi^2-σj^2)^2)/4+σij^2))
    end
    return σi, σj, a, b, θ
end

function preparecanvas(LaTeX_array, limits, ticks, probes, colors, PlotPars::Dict)
    matrix_dim = length(LaTeX_array)
    #TODO: add the textsize to PlotPars
    figure = Makie.Figure(textsize = 40, font = PlotPars["font"])
    for i in 1:matrix_dim
        for j in 1:i
            if i == j
                ax = Axis(figure[i,i],
                    width = PlotPars["sidesquare"], height = PlotPars["sidesquare"],
                    xticklabelsize = PlotPars["dimticklabel"],
                    yticklabelsize = PlotPars["dimticklabel"], yaxisposition = (:right),
                    xlabel = L"%$((LaTeX_array)[i])", xlabelsize = PlotPars["parslabelsize"],
                    ylabel = L"P/P_\mathrm{max}",ylabelsize = PlotPars["PPmaxlabelsize"], yticks = [0,1],
                    xticklabelrotation = PlotPars["xticklabelrotation"],
                    xticks = ([ticks[i,1], 0.5*(ticks[i,1]+ticks[i,2]), ticks[i,2]],
                        [string(myi) for myi in round.([ticks[i,1], 0.5*(ticks[i,1]+ticks[i,2]), ticks[i,2]], sigdigits = 3)]))
                Makie.ylims!(ax, (-0.0,1.05))
                Makie.xlims!(ax, (limits[i,1],limits[i,2]))
                Makie.hideydecorations!(ax, ticks = false, ticklabels = false, label = false)
                if i != matrix_dim
                    ax.alignmode = Mixed(right = MakieLayout.Protrusion(0), bottom = MakieLayout.Protrusion(0), top= MakieLayout.Protrusion(0))
                    hidexdecorations!(ax, ticks = true, ticklabels = true,  label = true)
                else
                    hidexdecorations!(ax, ticks = false, ticklabels = false,  label = false)
                end
            else
                ax = Axis(figure[i,j], width = PlotPars["sidesquare"], height = PlotPars["sidesquare"],
                    xticklabelsize = PlotPars["dimticklabel"], yticklabelsize = PlotPars["dimticklabel"],
                    ylabel = L"%$(LaTeX_array[i])", xlabel = L"%$((LaTeX_array)[j])",
                    ylabelsize = PlotPars["parslabelsize"], xlabelsize = PlotPars["parslabelsize"], xticklabelrotation = PlotPars["xticklabelrotation"],
                    yticks = ([ticks[i,1], 0.5*(ticks[i,1]+ticks[i,2]), ticks[i,2]],
                        [string(myi) for myi in round.([ticks[i,1], 0.5*(ticks[i,1]+ticks[i,2]), ticks[i,2]], sigdigits = 3)]),
                    xticks = ([ticks[j,1], 0.5*(ticks[j,1]+ticks[j,2]), ticks[j,2]],
                        [string(myi) for myi in round.([ticks[j,1], 0.5*(ticks[j,1]+ticks[j,2]), ticks[j,2]], sigdigits = 3)]),
                yticklabelpad=8)
                Makie.ylims!(ax, (limits[i,1],limits[i,2]))
                Makie.xlims!(ax, (limits[j,1],limits[j,2]))
                if i == matrix_dim
                    hidexdecorations!(ax, ticks = false, ticklabels = false,  label = false)
                else
                    hidexdecorations!(ax, ticks = true, ticklabels = true,  label = true)
                end
                if j == 1
                    hideydecorations!(ax, ticks = false, ticklabels = false,  label = false)
                    Legend(figure[1,matrix_dim],
                    [PolyElement(color = color, strokecolor = color, strokewidth = 1) for color in colors],
                    probes,
                    tellheight = false, tellwidth = false, rowgap = 10,
                    halign = :right, valign = :top, framecolor = :black, labelsize =55, patchsize = (70, 40), framevisible = true)
                else
                    hideydecorations!(ax, ticks = true, ticklabels = true,  label = true)
                    ax.alignmode = Mixed(right = MakieLayout.Protrusion(0), bottom = MakieLayout.Protrusion(0), top = MakieLayout.Protrusion(0))
                end
            end
        end
    end
    Makie.resize!(figure.scene, figure.layout.layoutobservables.reportedsize[]...)
    return figure
end

function drawgaussian!(canvas, σ, i, central, color)
    ax = canvas[i,i]
    x = Array(LinRange(-4σ+central,4σ+central, 200))
    Makie.lines!(ax, x, gaussian.(central, σ, x)./gaussian.(central, σ, central), color = color, linewidth = 4)
    Makie.band!(ax, x, 0, gaussian.(central, σ, x)./gaussian.(central, σ, central) , color=(color, 0.2))
    x = Array(LinRange(-σ+central,σ+central, 200))
    Makie.band!(ax, x, 0, gaussian.(central, σ, x)./gaussian.(central, σ, central) , color=(color, 0.4))
end

function drawellipse!(canvas, i, j, x, y, central_values, color)
    ax = canvas[i,j]

    Makie.lines!(ax, x .+ central_values[j], y .+ central_values[i], color = color, linewidth = 4)
    Makie.lines!(ax, 3x .+ central_values[j], 3y .+ central_values[i], color = color, linewidth = 4)

    Makie.band!(ax, x .+ central_values[j], 0, y .+ central_values[i], color=(color, 0.4))
    Makie.band!(ax, 3x .+ central_values[j], 0, 3y .+ central_values[i], color=(color, 0.2))
end

function paintcorrmatrix!(canvas, central_values, corr_matrix, color)
    for i in 1:length(central_values)
        for j in 1:i
            if i == j
                drawgaussian!(canvas, sqrt(corr_matrix[i,i]), i, central_values[i], color)
            else
                σi, σj, a, b, θ = ellipseparameters(corr_matrix, i,j)
                x,y = ellipseparameterization(a, b, θ)
                drawellipse!(canvas, i, j, x, y, central_values, color)
            end
        end
    end
end

function save(filename, obj)
    Makie.save(filename, obj)
end

end # module

max_stars_repo_name: marcobonici/FisherPlot

max_issues_repo_name: marcobonici/FisherPlot

max_forks_repo_name: marcobonici/FisherPlot


text:
using SellmeierFit
using Downloads
using SimpleConstants
using Test

include("read.jl")
include("model.jl")
include("fit.jl")


text:
@testset "Searcher: build, utils" begin
    local cfg
    for config in CONFIG_FUNCTIONS
        cfg = mktemp() do path, io  # write and parse config file on-the-fly
                   write(io, eval(config)())
                   flush(io)
                   parse_configuration(path)
               end
    end
    # Load data and check primary id
    dbdata = cfg.data_loader()

    # Build embedders
    embedders = [Garamond.build_embedder(dbdata,
                                         embdr_config;
                                         vectors_eltype=cfg.vectors_eltype,
                                         id_key=cfg.id_key)
                 for embdr_config in cfg.embedder_configs]

    # Build searchers
    searchers = [Garamond.build_searcher(dbdata,
                                         embedders,
                                         srcher_config;
                                         id_key=cfg.id_key)
                 for srcher_config in cfg.searcher_configs]
    @test searchers isa Vector{<:AbstractSearcher{eval(cfg.vectors_eltype)}}

    for srcher in searchers
        @test srcher.config.id == id(srcher)
        @test srcher.config.description == description(srcher)
    end

    idx_searcher = 1
    disable!(searchers[idx_searcher])
    @test !isenabled(searchers[idx_searcher])
    enable!(searchers[idx_searcher])
    @test isenabled(searchers[idx_searcher])
end


text:
x = rand(5,5);
println(x);


text:
# test EuclidDistance

using IncrementalInference
using Test

##

@testset "test EuclidDistance on 1 dim" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0], Prior(Normal()))

addVariable!(fg, :x1, ContinuousScalar)

eud = EuclidDistance(Normal(10,1))
addFactor!(fg, [:x0;:x1], eud)

##

tree, _, = solveTree!(fg)

##

@test isapprox(getPPE(fg, :x0).suggested[1], 0, atol=1)

pts = getBelief(fg, :x1) |> getPoints
N = size(pts, 2)

@test 0.3*N < sum( 5 .< pts )
@test 0.3*N < sum( pts .< -5 )
@test sum( -5 .< pts .< 5 ) < 0.1*N

##

end

@testset "test EuclidDistance on 2 dim" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousEuclid{2})
addFactor!(fg, [:x0], Prior(MvNormal(zeros(2),diagm([1;1.0]))))

addVariable!(fg, :x1, ContinuousEuclid{2})

eud = EuclidDistance(Normal(10,1))
addFactor!(fg, [:x0;:x1], eud)

tree, _, = solveTree!(fg)

@test isapprox(getPPE(fg, :x0).suggested[1], 0, atol=1)
@test isapprox(getPPE(fg, :x0).suggested[1], 0, atol=1)

pts = getBelief(fg, :x1) |> getPoints
N = size(pts, 2)

pts .^= 2
@test 0.5*N < sum( 7 .< sqrt.(sum(pts, dims=1)) .< 13 )

##

end


text:
"""
simulateCorrelatedData(ρMat, n; σₙ, σ, μ)

Generate a matrix of random values from a normal distribution,
where columns are correlated according to the `ρMat` correlation matrix.

# Arguments
- `ρMat::Matrix{Real}`: coefficients of the correlation matrix
- `n::Int64`: number of samples, default `n` is `100`
- `σₙ::Real`: noise standard deviation, default `σₙ` is `0`
- `σ::Real`: simulated data standard deviation, default `σ` is `1`
- `μ::Real`: simulated data mean, default `μ` is `0`

# Examples
```julia
simulateCorrelatedData([1.0   0.5  0.25;
                        0.5   1.0  0.5 ;
                        0.25  0.5  1.0])
```
"""

function simulateCorrelatedData(ρMat::Matrix{Float64}, n::Int64=100;
                                σₙ::Real=0, σ::Real = 1, μ::Real= 0)

    # Get number of variables
    p = size(ρMat)[2]

    # Whitening transformation of the random samples matrix
    X = standardize(ZScoreTransform, rand(Normal(0, 1), n, p), dims = 1)
    X = X / LinearAlgebra.cholesky(cov(X)).U

    # Apply correlation matrix
    X = X * LinearAlgebra.cholesky(ρMat).U
    X = X *  σ .+ μ .+ rand(Normal(0, σₙ), n)

    return X

end