text:
"""
    Integrator(;name, k=1, x_start=0.0)

Outputs `y = ∫k*u dt`, corresponding to the transfer function `1/s`.

# Connectors:
- `input`
- `output`

# Parameters:
- `k`: Gain of integrator
- `x_start`: Initial value of integrator
"""
function Integrator(;name, k=1, x_start=0.0)
    @named siso = SISO()
    @unpack u, y = siso
    sts = @variables x(t)=x_start
    pars = @parameters k=k
    eqs = [
        D(x) ~ k * u
        y ~ x
    ]
    extend(ODESystem(eqs, t, sts, pars; name=name), siso)
end

"""
    Derivative(; name, k=1, T, x_start=0.0)

Outputs an approximate derivative of the input. The transfer function of this block is
```
k       k
─ - ──────────
T    2 ⎛    1⎞
    T ⋅⎜s + ─⎟
       ⎝    T⎠
```
and a state-space realization is given by `ss(-1/T, 1/T, -k/T, k/T)`
where `T` is the time constant of the filter.
A smaller `T` leads to a more ideal approximation of the derivative.

# Parameters:
- `k`: Gain
- `T`: [s] Time constants (T>0 required; T=0 is ideal derivative block)
- `x_start`: Initial value of state

# Connectors:
- `input`
- `output`
"""
function Derivative(; name, k=1, T, x_start=0.0)
    T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive"))
    @named siso = SISO()
    @unpack u, y = siso
    sts = @variables x(t)=x_start
    pars = @parameters T=T k=k
    eqs = [
        D(x) ~ (u - x) / T
        y ~ (k / T) * (u - x)
    ]
    extend(ODESystem(eqs, t, sts, pars; name=name), siso)
end

"""
    FirstOrder(; name, k=1, T, x_start=0.0)

A first-order filter with a single real pole in `s = -T` and gain `k`. The transfer function
is given by `Y(s)/U(s) = `
```
   k
───────
sT + 1
```

# Parameters:
- `k`: Gain
- `T`: [s] Time constants (T>0 required)
- `x_start`: Initial value of state

# Connectors:
- `input`
- `output`

See also [`SecondOrder`](@ref)
"""
function FirstOrder(; name, k=1, T, x_start=0.0)
    T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive"))
    @named siso = SISO()
    @unpack u, y = siso
    sts = @variables x(t)=x_start
    pars = @parameters T=T k=k
    eqs = [
        D(x) ~ (k*u - x) / T
        y ~ x
    ]
    extend(ODESystem(eqs, t, sts, pars; name=name), siso)
end

"""
    SecondOrder(; name, k=1, w, d, x_start=0.0, xd_start=0.0)

A second-order filter with gain `k`, a bandwidth of `w` rad/s and relative damping `d`. The transfer function
is given by `Y(s)/U(s) = `
```
      k*w^2
─────────────────
s² + 2d*w*s + w^2
```
Critical damping corresponds to `d=1`, which yields the fastest step response without overshoot, d < 1` results in an under-damped filter while `d > 1` results in an over-damped filter.
`d = 1/√2` corresponds to a Butterworth filter of order 2 (maximally flat frequency response).

# Parameters:
- `k`: Gain
- `w`: Angular frequency
- `d`: Damping
- `x_start`: Initial value of state (output)
- `xd_start`: Initial value of derivative of state (output)

# Connectors:
- `input`
- `output`
"""
function SecondOrder(; name, k=1, w, d, x_start=0.0, xd_start=0.0)
    @named siso = SISO()
    @unpack u, y = siso
    sts = @variables x(t)=x_start xd(t)=xd_start
    pars = @parameters k=k w=w d=d
    eqs = [
        D(x) ~ xd
        D(xd) ~ w*(w*(k*u - x) - 2*d*xd)
        y ~ x
    ]
    extend(ODESystem(eqs, t, sts, pars; name=name), siso)
end

"""
    PI(;name, k=1, T, x_start=0.0)

Textbook version of a PI-controller without actuator saturation and anti-windup measure.

# Parameters:
- `k`: Gain
- `T`: [s] Integrator time constant (T>0 required)
- `x_start`: Initial value for the integrator

# Connectors:
- `err_input`
- `ctr_output`

See also [`LimPI`](@ref)
"""
function PI(;name, k=1, T, x_start=0.0)
    T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive"))
    @named err_input = RealInput() # control error
    @named ctr_output = RealOutput() # control signal
    @named gainPI = Gain(k)
    @named addPI = Add()
    @named int = Integrator(k=1/T, x_start=x_start)
    sys = [err_input, ctr_output, gainPI, addPI, int]
    eqs = [
        connect(err_input, addPI.input1),
        connect(addPI.output, gainPI.input),
        connect(gainPI.output, ctr_output),
        connect(err_input, int.input),
        connect(int.output, addPI.input2),
    ]
    ODESystem(eqs, t, [], []; name=name, systems=sys)
end

"""
    PID(;name, k=1, Ti=false, Td=false, Nd=10, xi_start=0, xd_start=0)

Text-book version of a PID-controller without actuator saturation and anti-windup measure.

# Parameters:
- `k`: Gain
- `Ti`: [s] Integrator time constant (Ti>0 required). If set to false no integral action is used.
- `Td`: [s] Derivative time constant (Td>0 required). If set to false no derivative action is used.
- `Nd`: [s] Time constant for the derivative approximation (Nd>0 required; Nd=0 is ideal derivative).
- `x_start`: Initial value for the integrator.
- `xd_start`: Initial value for the derivative state.

# Connectors:
- `err_input`
- `ctr_output`

See also [`LimPID`](@ref)
"""
function PID(;name, k=1, Ti=false, Td=false, Nd=10, xi_start=0, xd_start=0)
    with_I = !isequal(Ti, false)
    with_D = !isequal(Td, false)
    @named err_input = RealInput() # control error
    @named ctr_output = RealOutput() # control signal
    !isequal(Ti, false) && (Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0")))
    !isequal(Td, false) && (Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0")))
    Nd > 0 || throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0"))

    @named gainPID = Gain(k)
    @named addPID = Add3()
    if with_I
        @named int = Integrator(k=1/Ti, x_start=xi_start)
    else
        @named Izero = Constant(k=0)
    end
    if with_D
        @named der = Derivative(k=1/Td, T=1/Nd, x_start=xd_start)
    else
        @named Dzero = Constant(k=0)
    end
    sys = [err_input, ctr_output, gainPID, addPID]
    if with_I
        push!(sys, int)
    else
        push!(sys, Izero)
    end
    if with_D
        push!(sys, der)
    else
        push!(sys, Dzero)
    end
    eqs = [
        connect(err_input, addPID.input1),
        connect(addPID.output, gainPID.input),
        connect(gainPID.output, ctr_output)
    ]
    if with_I
        push!(eqs, connect(err_input, int.input))
        push!(eqs, connect(int.output, addPID.input2))
    else
        push!(eqs, connect(Izero.output, addPID.input2))
    end
    if with_D
        push!(eqs, connect(err_input, der.input))
        push!(eqs, connect(der.output, addPID.input3))
    else
        push!(eqs, connect(Dzero.output, addPID.input3))
    end
    ODESystem(eqs, t, [], []; name=name, systems=sys)
end

"""
    LimPI(;name, k=1, T, u_max=1, u_min=-u_max, Ta)

Text-book version of a PI-controller with actuator saturation and anti-windup measure.

# Parameters:
- `k`: Gain
- `T`: [s] Integrator time constant (T>0 required)
- `Ta`: [s] Tracking time constant (Ta>0 required)
- `x_start`: Initial value for the integrator

# Connectors:
- `err_input`
- `ctr_output`
"""
function LimPI(;name, k=1, T, u_max=1, u_min=-u_max, Ta, x_start=0.0)
    Ta > 0 || throw(ArgumentError("Time constant `Ta` has to be strictly positive"))
    T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive"))
    u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max"))
    @named err_input = RealInput() # control error
    @named ctr_output = RealOutput() # control signal
    @named gainPI = Gain(k)
    @named addPI = Add()
    @named addTrack = Add()
    @named int = Integrator(k=1/T, x_start=x_start)
    @named limiter = Limiter(y_max=u_max, y_min=u_min)
    @named addSat = Add(k1=1, k2=-1)
    @named gainTrack = Gain(1/Ta)
    sys = [err_input, ctr_output, gainPI, addPI, int, addTrack, limiter, addSat, gainTrack]
    eqs = [
        connect(err_input, addPI.input1),
        connect(addPI.output, gainPI.input),
        connect(gainPI.output, limiter.input),
        connect(limiter.output, ctr_output),
        connect(limiter.input, addSat.input2),
        connect(limiter.output, addSat.input1),
        connect(addSat.output, gainTrack.input),
        connect(err_input, addTrack.input1),
        connect(gainTrack.output, addTrack.input2),
        connect(addTrack.output, int.input),
        connect(int.output, addPI.input2),
    ]
    ODESystem(eqs, t, [], []; name=name, systems=sys)
end

"""
    LimPID(; k, Ti=false, Td=false, wp=1, wd=1, Ni, Nd=12, u_max=Inf, u_min=-u_max, gains = false, name)

Proportional-Integral-Derivative (PID) controller with output saturation, set-point weighting and integrator anti-windup.

The equation for the control signal is roughly
```
k(ep + 1/Ti * ∫e + 1/Td * d/dt(ed))
e = u_r - u_y
ep = wp*u_r - u_y
ed = wd*u_r - u_y
```
where the transfer function for the derivative includes additional filtering, see `? Derivative` for more details.

# Parameters:
- `k`: Proportional gain
- `Ti`: [s] Integrator time constant. Set to `false` to turn off integral action.
- `Td`: [s] Derivative time constant. Set to `false` to turn off derivative action.
- `wp`: [0,1] Set-point weighting in the proportional part.
- `wd`: [0,1] Set-point weighting in the derivative part.
- `Nd`: [1/s] Derivative limit, limits the derivative gain to Nd/Td. Reasonable values are ∈ [8, 20]. A higher value gives a better approximation of an ideal derivative at the expense of higher noise amplification.
- `Ni`: `Ni*Ti` controls the time constant `Ta` of anti-windup tracking. A common (default) choice is `Ta = √(Ti*Td)` which is realized by `Ni = √(Td / Ti)`. Anti-windup can be effectively turned off by setting `Ni = Inf`.
` `gains`: If `gains = true`, `Ti` and `Td` will be interpreted as gains with a fundamental PID transfer function on parallel form `ki=Ti, kd=Td, k + ki/s + kd*s`

# Connectors:
- `reference`
- `measurement`
- `ctr_output`
"""
function LimPID(; name, k=1, Ti=false, Td=false, wp=1, wd=1,
    Ni= Ti == 0 ? Inf : √(max(Td / Ti, 1e-6)),
    Nd=10,
    u_max=Inf,
    u_min=u_max > 0 ? -u_max : -Inf,
    gains=false,
    xi_start=0.0,
    xd_start=0.0,
    )
    with_I = !isequal(Ti, false)
    with_D = !isequal(Td, false)
    with_AWM = Ni != Inf
    if gains
        Ti = k / Ti
        Td = Td / k
    end
    0 ≤ wp ≤ 1 || throw(ArgumentError("wp out of bounds, got $(wp) but expected wp ∈ [0, 1]"))
    0 ≤ wd ≤ 1 || throw(ArgumentError("wd out of bounds, got $(wd) but expected wd ∈ [0, 1]"))
    !isequal(Ti, false) && (Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0")))
    !isequal(Td, false) && (Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0")))
    u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max"))
    Nd > 0 || throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0"))

    @named reference = RealInput()
    @named measurement = RealInput()
    @named ctr_output = RealOutput() # control signal
    @named addP = Add(k1=wp, k2=-1)
    @named gainPID = Gain(k)
    @named addPID = Add3()
    @named limiter = Limiter(y_max=u_max, y_min=u_min)
    if with_I
        if with_AWM
            @named addI = Add3(k1=1, k2=-1, k3=1)
            @named addSat = Add(k1=1, k2=-1)
            @named gainTrack = Gain(1/(k * Ni))
        else
            @named addI = Add(k1=1, k2=-1)
        end
        @named int = Integrator(k=1/Ti, x_start=xi_start)
    else
        @named Izero = Constant(k=0)
    end
    if with_D
        @named der = Derivative(k=1/Td, T=1/Nd, x_start=xd_start)
        @named addD = Add(k1=wd, k2=-1)
    else
        @named Dzero = Constant(k=0)
    end

    sys = [reference, measurement, ctr_output, addP, gainPID, addPID, limiter]
    if with_I
        if with_AWM
            push!(sys, [addSat, gainTrack]...)
        end
        push!(sys, [addI, int]...)
    else
        push!(sys, Izero)
    end
    if with_D
        push!(sys, [addD, der]...)
    else
        push!(sys, Dzero)
    end

    eqs = [
        connect(reference, addP.input1),
        connect(measurement, addP.input2),
        connect(addP.output, addPID.input1),
        connect(addPID.output, gainPID.input),
        connect(gainPID.output, limiter.input),
        connect(limiter.output, ctr_output),
    ]
    if with_I
        push!(eqs, connect(reference, addI.input1))
        push!(eqs, connect(measurement, addI.input2))
        if with_AWM
            push!(eqs, connect(limiter.input, addSat.input2))
            push!(eqs, connect(limiter.output, addSat.input1))
            push!(eqs, connect(addSat.output, gainTrack.input))
            push!(eqs, connect(gainTrack.output, addI.input3))
        end
        push!(eqs, connect(addI.output, int.input))
        push!(eqs, connect(int.output, addPID.input3))
    else
        push!(eqs, connect(Izero.output, addPID.input3))
    end
    if with_D
        push!(eqs, connect(reference, addD.input1))
        push!(eqs, connect(measurement, addD.input2))
        push!(eqs, connect(addD.output, der.input))
        push!(eqs, connect(der.output, addPID.input2))
    else
        push!(eqs, connect(Dzero.output, addPID.input2))
    end

    ODESystem(eqs, t, [], []; name=name, systems=sys)
end

"""
    StateSpace(A, B, C, D=0; x_start=zeros(size(A,1)), name)

A linear, time-invariant state-space system on the form.
```
ẋ = Ax + Bu
y = Cx + Du
```
Transfer functions can also be simulated by converting them to a StateSpace form.
"""
function StateSpace(;A, B, C, D=nothing, x_start=zeros(size(A,1)), name)
    nx, nu, ny = size(A,1), size(B,2), size(C,1)
    size(A,2) == nx || error("`A` has to be a square matrix.")
    size(B,1) == nx || error("`B` has to be of dimension ($nx x $nu).")
    size(C,2) == nx || error("`C` has to be of dimension ($ny x $nx).")
    if B isa AbstractVector
        B = reshape(B, length(B), 1)
    end
    if isnothing(D) || iszero(D)
        D = zeros(ny, nu)
    else
        size(D) == (ny,nu) || error("`D` has to be of dimension ($ny x $nu).")
    end
    @named input = RealInput(nin=nu)
    @named output = RealOutput(nout=ny)
    @variables x[1:nx](t)=x_start
    # pars = @parameters A=A B=B C=C D=D # This is buggy
    eqs = [ # FIXME: if array equations work
        [Differential(t)(x[i]) ~ sum(A[i,k] * x[k] for k in 1:nx) + sum(B[i,j] * input.u[j] for j in 1:nu) for i in 1:nx]..., # cannot use D here
        [output.u[j] ~ sum(C[j,i] * x[i] for i in 1:nx) + sum(D[j,k] * input.u[k] for k in 1:nu) for j in 1:ny]...,
    ]
    compose(ODESystem(eqs, t, vcat(x...), [], name=name), [input, output])
end

StateSpace(A, B, C, D=nothing; kwargs...) = StateSpace(; A, B, C, D, kwargs...)


text:
using KernelAbstractions, Test, CUDA

if has_cuda_gpu()
    CUDA.allowscalar(false)
end

@testset "Error propagation" begin
    let event = Event(()->error(""))
        @test_throws TaskFailedException wait(event)
    end

    let event = Event(error, "")
        @test_throws CompositeException wait(MultiEvent(event))
    end

    let event = Event(error, "")
        event = Event(wait, MultiEvent(event))
        @test_throws TaskFailedException wait(event)
    end

    let event = Event(error, "")
        event = Event(()->nothing, dependencies=event)
        @test_throws TaskFailedException wait(event)
    end
end

if has_cuda_gpu()
    barrier = Base.Threads.Event()
    cpu_event = Event(wait, barrier)

    KernelAbstractions.unsafe_wait(CUDADevice(), cpu_event) # Event edge on CuDefaultStream
    gpu_event = Event(CUDADevice()) # Event on CuDefaultStream

    notify(barrier)
    wait(gpu_event)
end

ext: jl


text:
using Documenter, PackageCompilerX

makedocs(
    format = Documenter.HTML(
        # prettyurls on travis
        prettyurls = haskey(ENV, "HAS_JOSH_K_SEAL_OF_APPROVAL"),
    ),
    sitename = "PackageCompilerX",
    pages = Any[
        "Home" => "index.md",
        "Manual" => [
            "prereq.md"
            "sysimages.md"
            "apps.md"
        ],

        "Examples" => Any[
            "examples/ohmyrepl.md",
        ],
        "References" => "refs.md",
    ]
)

deploydocs(
    repo = "github.com/KristofferC/PackageCompilerX.jl.git",
)


text:
module TestFit

using Test
using Dates
using LightGBM

# test fixtures
train_matrix = rand(5000,70) # create random dataset
train_labels = rand([0, 1], 5000)
train_dataset = LightGBM.LGBM_DatasetCreateFromMat(train_matrix, "")
LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)

test_matrix = rand(2000,70) # create random dataset
test_labels = rand([0, 1], 2000)
test_dataset = LightGBM.LGBM_DatasetCreateFromMat(test_matrix, "", train_dataset)
LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)

test2_matrix = rand(1500,70) # create second random dataset
test2_labels = rand([0, 1], 1500)
test2_dataset = LightGBM.LGBM_DatasetCreateFromMat(test2_matrix, "", train_dataset)
LightGBM.LGBM_DatasetSetField(test2_dataset, "label", test2_labels)

@testset "test fit! with dataset -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(objective = "binary", num_class = 1)

    # Act
    LightGBM.fit!(estimator, train_dataset, test_dataset, verbosity = -1);
    p_binary = LightGBM.predict(estimator, train_matrix, verbosity = -1)
    binary_classes = LightGBM.predict_classes(estimator, train_matrix, verbosity = -1)

    # Assert
    # check also that binary outputs and classes are invariant in prediction output shapes, i.e. they're matrices
    @test length(size(p_binary)) == 2
    @test length(size(binary_classes)) == 2

end

@testset "test fit! with dataset without testset -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(objective = "binary", num_class = 1)

    # Act
    LightGBM.fit!(estimator, train_dataset, verbosity = -1);
    p_binary = LightGBM.predict(estimator, train_matrix, verbosity = -1)
    binary_classes = LightGBM.predict_classes(estimator, train_matrix, verbosity = -1)

    # Assert
    # check also that binary outputs and classes are invariant in prediction output shapes, i.e. they're matrices
    @test length(size(p_binary)) == 2
    @test length(size(binary_classes)) == 2

end

@testset "test train! single test set -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        objective = "binary",
        num_class = 1,
        is_training_metric = true,
        metric = ["auc"],
    )

    bst_parameters = LightGBM.stringifyparams(estimator; verbosity=-1)
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)

    # Act
    results = LightGBM.train!(estimator, 10, ["test_set_with_whatever_label"], -1, DateTime(2020, 1, 11, 11, 00));

    # Assert
    @test length(results["metrics"]["training"]["auc"]) > 0
    @test length(results["metrics"]["test_set_with_whatever_label"]["auc"]) > 0

end

@testset "test train! multiple test set multiple metrics -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        objective = "binary",
        num_class = 1,
        is_training_metric = true,
        metric = ["auc", "l2"],
    )

    bst_parameters = LightGBM.stringifyparams(estimator; verbosity=-1)
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test2_dataset)

    # Act
    results = LightGBM.train!(estimator, 10, ["test_1", "test_2"], -1, DateTime(2020, 1, 11, 11, 00));

    # Assert
    @test length(results["metrics"]["training"]["auc"]) > 0
    @test length(results["metrics"]["training"]["l2"]) > 0
    @test length(results["metrics"]["test_1"]["auc"]) > 0
    @test length(results["metrics"]["test_1"]["l2"]) > 0
    @test length(results["metrics"]["test_2"]["auc"]) > 0
    @test length(results["metrics"]["test_2"]["l2"]) > 0

end

@testset "test train! multiple test set multiple metrics -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        objective = "binary",
        num_class = 1,
        is_training_metric = true,
        metric = ["auc", "l2"],
    )

    bst_parameters = LightGBM.stringifyparams(estimator; verbosity=-1)
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test2_dataset)

    # Act
    results = LightGBM.train!(estimator, 10, ["test_1", "test_2"], -1, DateTime(2020, 1, 11, 11, 00));

    # Assert
    @test length(results["metrics"]["training"]["auc"]) > 0
    @test length(results["metrics"]["training"]["l2"]) > 0
    @test length(results["metrics"]["test_1"]["auc"]) > 0
    @test length(results["metrics"]["test_1"]["l2"]) > 0
    @test length(results["metrics"]["test_2"]["auc"]) > 0
    @test length(results["metrics"]["test_2"]["l2"]) > 0

end

@testset "test store_scores! adds results" begin
    # Arrange
    results_fixture = Dict("metrics" => Dict{String,Dict{String,Vector{Float64}}}())

    # Act
    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_metric",0.5)
    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_metric",0.55)

    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_extra_metric",0.6)
    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_extra_metric",0.65)

    LightGBM.store_scores!(results_fixture, "whatever_second_dataset","some_metric",0.7)
    LightGBM.store_scores!(results_fixture, "whatever_second_dataset","some_metric",0.75)

    # Assert
    @test results_fixture["metrics"]["whatever_dataset"]["some_metric"] == [0.5, 0.55]
    @test results_fixture["metrics"]["whatever_dataset"]["some_extra_metric"] == [0.6, 0.65]
    @test results_fixture["metrics"]["whatever_second_dataset"]["some_metric"] == [0.7, 0.75]

end

@testset "test eval_metrics! early stop disabled" begin
    # Arrange
    results_fixture = Dict(
        "best_iter" => 0,
        "metrics" => Dict{String,Dict{String,Vector{Float64}}}(),
    )
    estimator = LightGBM.LGBMClassification(
        num_iterations = 10, objective = "binary", num_class = 1,
        is_training_metric = false, metric = ["auc"],
        early_stopping_round = 0 # default value, but stating explicitly to test!
    )

    bst_parameters = LightGBM.stringifyparams(estimator; verbosity=-1)
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)

    bigger_is_better = Dict{String,Float64}("auc" => 1.0)
    best_scores = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))
    best_iterations = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))

    # Act and assert each round returns false (despite iter 1 have best score always)
    for iter in 1:10
        LightGBM.LGBM_BoosterUpdateOneIter(estimator.booster)
        output = LightGBM.eval_metrics!(
            results_fixture, estimator, ["test_bla"], iter, -1,
            bigger_is_better, best_scores, best_iterations, ["auc"]
        )

        @test output == false

    end
end

#=
This test assumes:
* early_stopping_round set to 5
* total iterations = 10
* best iter is ALWAYS round 1 of metric score 1 (auc)
Criteria: early_stopping should kick in on round 6
=#
@testset "test eval_metrics! stops correctly" begin
    # Arrange
    results_fixture = Dict(
        "best_iter" => 0,
        "metrics" => Dict{String,Dict{String,Vector{Float64}}}(),
    )
    estimator = LightGBM.LGBMClassification(
        num_iterations = 10, objective = "binary", num_class = 1,
        is_training_metric = false, metric = ["auc"],
        early_stopping_round = 5
    )

    bst_parameters = LightGBM.stringifyparams(estimator; verbosity=-1)
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)

    bigger_is_better = Dict{String,Float64}("auc" => 1.0)
    best_scores = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))
    best_iterations = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))

    # Act and assert each round returns expected output
    for iter in 1:10
        LightGBM.LGBM_BoosterUpdateOneIter(estimator.booster)
        output = LightGBM.eval_metrics!(
            results_fixture, estimator, ["test_bla"], iter, -1,
            bigger_is_better, best_scores, best_iterations, ["auc"]
        )

        # reset scores to round 1 being best
        best_scores = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))
        best_iterations = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))

        if iter < 6
            @test output == false
        else
            @test output == true
        end
    end
end

@testset "test truncate_model!" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(num_iterations = 100)
    verbosity = "verbose=-1"

    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, verbosity)

    # update learning 100 times
    for _ in [1:100;]
        finished = LightGBM.LGBM_BoosterUpdateOneIter(estimator.booster)
    end
    @test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) == 100

    # Act
    LightGBM.truncate_model!(estimator, 87)

    # Assert
    @test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) == 87
end

@testset "test fit when truncate_booster" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        num_class = 1,
        num_iterations = 500,
        early_stopping_round = 5,
        metric = ["auc"],
        objective = "binary",
    )
    verbosity = "verbose=-1"

    # Act
    output = LightGBM.fit!(estimator, train_dataset, test_dataset; truncate_booster=true, verbosity=-1)

    # Assert
    eval_metrics_run_count = length(output["metrics"]["test_1"]["auc"])
    @test eval_metrics_run_count >= 5 # at least this will be more than or equal to early_stopping_round
    @test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) == output["best_iter"] # test this matches as expected

end

@testset "test fit when NOT truncate_booster" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        num_class = 1,
        num_iterations = 200,
        early_stopping_round = 5,
        metric = ["auc"],
        objective = "binary",
    )
    verbosity = "verbose=-1"

    # Act
    output = LightGBM.fit!(estimator, train_dataset, test_dataset; truncate_booster=false, verbosity=-1)

    # Assert
    eval_metrics_run_count = length(output["metrics"]["test_1"]["auc"])
    @test eval_metrics_run_count >= 5 # at least this will be more than or equal to early_stopping_round
    @test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) >= 5 # given no truncation, this must be greater than early_stopping_round
    @test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) > output["best_iter"] # unless this runs to the very last num_iterations this should pass...should buy the lottery if this fails

end

@testset "test fit when truncate_booster inactive when no early_stopping" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        num_class = 1,
        num_iterations = 100,
        early_stopping_round = 0,
        metric = ["auc"],
        objective = "binary",
    )
    verbosity = "verbose=-1"

    # Act
    output = LightGBM.fit!(estimator, train_dataset, test_dataset; truncate_booster=true, verbosity=-1)

    # Assert
    eval_metrics_run_count = length(output["metrics"]["test_1"]["auc"])
    @test eval_metrics_run_count == 100
    @test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) == 100
    @test output["best_iter"] == 0

end

end # module


text:
module Bench_test

using SpikeSorting, JLD

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

a=load(string(dirname(Base.source_path()),"/data/spikes2.jld"))
time_stamps=a["time_stamps"]
fv=a["fv"]

detect=DetectNeg()
cluster=ClusterOSort(100,50)
align=AlignMin()
feature=FeatureTime()
reduce=ReductionNone()
thres=ThresholdMeanN(2.0)
s1=create_multi(detect,cluster,align,feature,reduce,thres,1)

ss=SpikeSorting.benchmark(fv,time_stamps,s1[1],180.0,20000)

detected_spikes=sum([length(ss[1][i]) for i=1:length(ss[1])])
tpfp=ss[3]+sum(ss[4])
tpfn=ss[3]+sum(ss[5])

    @test detected_spikes == tpfp
    @test sum(ss[2]) == tpfn

ss2=SpikeSorting.benchmark(fv,time_stamps,s1[1],60.0,20000)

detected_spikes2=sum([length(ss2[1][i]) for i=1:length(ss2[1])])
tpfp2=ss2[3]+sum(ss2[4])
tpfn2=ss2[3]+sum(ss2[5])

    @test detected_spikes2 == tpfp2
    @test sum(ss2[2]) == tpfn2
    @test detected_spikes < less_than(detected_spikes2)

end


text:
using UUIDs,Base.Threads,Observables
d = Dict()

function g(p,uuid)
    sleep(1)
    for i in 1:30
      if i % 1 == 0
        s = 0.0
        for j in 1:100_000_000
            s += rand()
        end
        p[]=T(uuid,(i,s))
       end
   end
end

function ad()
    o = Observable(T(0,0))
    on(o) do p
        d[p.a] = p.b
    end
    g(o,uuid)
end

uuid = UUIDs.uuid4()
@async ad(uuid)

function aad()
    uuid = UUIDs.uuid4()
    t = @async ad()
    (t, uuid)
end

#-- async using channels
using UUIDs,Base.Threads,Observables

struct T
    a
    b
end

progress = Dict()
data = Dict()

observer = Observable(T(0,0))
on(observer) do p
    progress[p.a] = p.b
end

inqueue = Channel{T}(32);
outqueue = Channel{T}(32);

function calc()
    global observer
    while true
        job = take!(inqueue)
        r = 0.0
        for i in 1:job.b
            r += rand()
            if i % 100 == 0
                observer[]= T(job.a,i)
            end
        end
        put!( outqueue, T( job.a, r ))
    end
end

function submitjob()
    uuid = UUIDs.uuid4()
    put!( inqueue, T(uuid, rand(1:1_000_000)))
    return uuid
end

function makejobs(n)
    for i in 1:n
        put!( inqueue, T(UUIDs.uuid4(), rand(1:1_000_000)))
    end
end

# i in 1:4 # start 4 tasks to process requests in parallel
#    errormonitor(@async calc())
# end

function clearjobs(n)
     while true
        res = take!( outqueue )
        data[res.a] = res
        print( "from clearjobs ")
        println( res )
    end
end

n = 4
for i in 1:n # start 4 tasks to process requests in parallel
    errormonitor(@async calc())
end
errormonitor(@async clearjobs(n))


text:
using Base: RefValue

# Interfaces

accum() = nothing
accum(x) = x

accum(x, y) =
  x === nothing ? y :
  y === nothing ? x :
  x + y

accum(x, y, zs...) = accum(accum(x, y), zs...)

accum(x::Tuple, y::Tuple) = accum.(x, y)
accum(x::AbstractArray, y::AbstractArray) = accum.(x, y)

@generated function accum(x::NamedTuple, y::NamedTuple)
  grad(x) = x in fieldnames(y) ? :(y.$x) : :nothing
  Expr(:tuple, [:($f=accum(x.$f, $(grad(f)))) for f in fieldnames(x)]...)
end

function accum(x::RefValue, y::RefValue)
  @assert x === y
  return x
end

# Core functions

@nograd Core.apply_type, Core.typeof, nfields, fieldtype, Core.TypeVar, Core.UnionAll,
  (==), (===), (<=), (>=), (<), (>), isempty, supertype, Base.typename,
  eps, Meta.parse, Base.eval, sleep, isassigned

@adjoint deepcopy(x) = deepcopy(x), ȳ -> (ȳ,)

@adjoint (::Type{V})(x...) where V<:Val = V(x...), _ -> nothing

@adjoint ifelse(cond::Bool, t, f) =
  ifelse(cond, t, f),
  Δ -> cond ? (nothing, Δ, zero(Δ)) : (nothing, zero(Δ), Δ)

@adjoint Base.typeassert(x, T) = Base.typeassert(x, T), Δ -> (Δ, nothing)

@generated function accum_param(cx::Context, x, Δ)
  isbitstype(x) && return :(Δ)
  quote
    if haskey(cache(cx), x)
      cache(cx)[x] = accum(cache(cx)[x],Δ)
      return
    else
      return Δ
    end
  end
end

function accum_global(cx::Context, ref, x̄)
  (x̄ === nothing || isconst(ref.mod, ref.name)) && return
  gs = cache(cx)
  gs[ref] = accum(get(gs, ref, nothing), x̄)
  return
end

unwrap(x) = x

@adjoint unwrap(x) = unwrap(x), x̄ -> (accum_param(__context__, x, x̄),)

unwrap(ref, x) = x

@adjoint unwrap(ref, x) = unwrap(x), function (x̄)
  accum_global(__context__, ref, x̄)
  (accum_param(__context__, x, x̄),)
end

function global_set(ref, val)
  ccall(:jl_set_global, Cvoid, (Any, Any, Any),
        ref.mod, ref.name, val)
end

@adjoint! function global_set(ref, x)
  global_set(ref, x), function (x̄)
    gs = cache(__context__)
    x̄ = accum(get(gs, ref, nothing), x̄)
    gs[ref] = nothing
    return (nothing, x̄)
  end
end

# Tuples

using Base: tail

@adjoint tuple(xs...) = xs, identity

@adjoint function literal_getindex(xs::NTuple{N,Any}, ::Val{i}) where {N,i}
  val = xs[i]
  function back(Δ)
    accum_param(__context__, val, Δ) === nothing && return
    return ntuple(j -> i == j ? Δ : nothing, Val(N)), nothing
  end
  val, back
end

@adjoint function getindex(xs::NTuple{N,Any}, i::Integer) where N
  val = xs[i]
  function back(Δ)
    accum_param(__context__, val, Δ) === nothing && return
    return ntuple(j -> i == j ? Δ : nothing, Val(N)), nothing
  end
  return val, back
end

@adjoint getindex(xs::NTuple{N,Any}, r::AbstractUnitRange) where N =
  (xs[r], Δ -> (ntuple(j -> j in r ? Δ[findfirst(isequal(j), r)] : nothing, Val(N)), nothing))

function _pullback(cx::Context, ::typeof(literal_indexed_iterate), xs::Tuple, ::Val{i}) where i
  y, b = _pullback(cx, literal_getindex, xs, Val(i))
  back(::Nothing) = nothing
  back(ȳ) = b(ȳ[1])
  (y, i+1), back
end

function _pullback(cx::Context, ::typeof(literal_indexed_iterate), xs::Tuple, ::Val{i}, st) where i
  y, b = _pullback(cx, literal_getindex, xs, Val(i))
  back(::Nothing) = nothing
  back(ȳ) = (b(ȳ[1])..., nothing)
  (y, i+1), back
end

# Needed for iteration lowering
@adjoint Core.getfield(xs::NTuple{N,Any}, i::Integer) where N =
  (xs[i], Δ -> (ntuple(j -> i == j ? Δ : nothing, Val(N)), nothing))

@adjoint Core.getfield(xs::NamedTuple{K,<:NTuple{N,Any}}, i::Integer) where {K,N} =
  (xs[i], Δ -> (NamedTuple{K}(ntuple(j -> i == j ? Δ : nothing, Val(N))), nothing))

@adjoint function Base.first(xs::Tuple)
  drest = map(_->nothing, tail(xs))
  first(xs), Δ -> ((Δ, drest...),)
end

@adjoint Base.tail(xs::Tuple) = tail(xs), x̄s -> ((nothing, x̄s...),)

_empty(x) = length(x)
_empty(x::Union{Tuple,NamedTuple}) = map(_->nothing, x)

_unapply(t::Integer, xs) = xs[1:t], xs[t+1:end]
_unapply(t, xs) = first(xs), tail(xs)
_unapply(t::Tuple{}, xs) = (), xs

function _unapply(t::Tuple, xs)
  t1, xs1 = _unapply(first(t), xs)
  t2, xs2 = _unapply(tail(t), xs1)
  (t1, t2...), xs2
end

function _unapply(t::NamedTuple{K}, xs) where K
  t, rst = _unapply(Tuple(t), xs)
  NamedTuple{K}(t), rst
end

unapply(t, xs) = _unapply(t, xs)[1]

@adjoint! function Core._apply(f, args...)
  y, back = Core._apply(_pullback, (__context__, f), args...)
  st = map(_empty, args)
  y, function (Δ)
    Δ = back(Δ)
    Δ === nothing ? nothing :
      (first(Δ), unapply(st, Base.tail(Δ))...)
  end
end

if VERSION >= v"1.4.0-DEV.304"
  @adjoint! function Core._apply_iterate(::typeof(iterate), f, args...)
    y, back = Core._apply(_pullback, (__context__, f), args...)
    st = map(_empty, args)
    y, function (Δ)
      Δ = back(Δ)
      Δ === nothing ? nothing :
        (nothing, first(Δ), unapply(st, Base.tail(Δ))...)
    end
  end
end

# Structs

deref!(x) = x

function deref!(x::Ref)
  d = x[]
  x[] = nothing
  return d
end

@generated nt_nothing(x) = Expr(:tuple, [:($f=nothing) for f in fieldnames(x)]...)

@generated pair(::Val{k}, v) where k = :($k = v,)

@adjoint function literal_getproperty(x, ::Val{f}) where f
  val = getproperty(x, f)
  function back(Δ)
    accum_param(__context__, val, Δ) === nothing && return
    if isimmutable(x)
      ((;nt_nothing(x)...,pair(Val(f), Δ)...), nothing)
    else
      dx = grad_mut(__context__, x)
      dx[] = (;dx[]...,pair(Val(f),accum(getfield(dx[], f), Δ))...)
      return (dx,nothing)
    end
  end
  unwrap(val), back
end

_pullback(cx::Context, ::typeof(getproperty), x, f::Symbol) =
  _pullback(cx, literal_getproperty, x, Val(f))

_pullback(cx::Context, ::typeof(getfield), x, f::Symbol) =
  _pullback(cx, literal_getproperty, x, Val(f))

_pullback(cx::Context, ::typeof(literal_getindex), x::NamedTuple, ::Val{f}) where f =
  _pullback(cx, literal_getproperty, x, Val(f))

_pullback(cx::Context, ::typeof(literal_getproperty), x::Tuple, ::Val{f}) where f =
  _pullback(cx, literal_getindex, x, Val(f))

grad_mut(x) = Ref{Any}(nt_nothing(x))

function grad_mut(cx::Context, x)
  ch = cache(cx)
  if haskey(ch, x)
    ch[x]
  else
    ch[x] = grad_mut(x)
  end
end

@adjoint! function setfield!(x, f, val)
  y = setfield!(x, f, val)
  g = grad_mut(__context__, x)
  y, function (_)
    Δ = getfield(g[], f)
    g[] = (;g[]...,pair(Val(f),nothing)...)
    (nothing, nothing, Δ)
  end
end

struct Jnew{T,G,splat}
  g::G
end

Jnew{T}(g) where T = Jnew{T,typeof(g)}(g)

@adjoint! function __new__(T, args...)
  x = __new__(T, args...)
  g = !T.mutable || fieldcount(T) == 0 ? nothing : grad_mut(__context__, x)
  x, Jnew{T,typeof(g),false}(g)
end

@adjoint! function __splatnew__(T, args)
  x = __splatnew__(T, args)
  g = !T.mutable || fieldcount(T) == 0 ? nothing : grad_mut(__context__, x)
  x, Jnew{T,typeof(g),true}(g)
end

# TODO captured mutables + multiple calls to `back`
@generated function (back::Jnew{T,G,false})(Δ::Union{NamedTuple,Nothing,RefValue}) where {T,G}
  !T.mutable && Δ == Nothing && return :nothing
  Δ = G == Nothing ? :Δ :
      Δ <: RefValue ? :(back.g[]) :
      :(accum(back.g[], Δ))
  quote
    x̄ = $Δ
    $(G == Nothing || :(back.g[] = nt_nothing($Δ)))
    (nothing, $(map(f -> :(x̄.$f), fieldnames(T))...))
  end
end

@generated function (back::Jnew{T,G,true})(Δ::Union{NamedTuple,Nothing,RefValue}) where {T,G}
  !T.mutable && Δ == Nothing && return :nothing
  Δ = G == Nothing ? :Δ : :(back.g)
  quote
    x̄ = $Δ
    $(G == Nothing || :($Δ = nt_nothing($Δ)))
    (nothing, ($(map(f -> :(x̄.$f), fieldnames(T))...),))
  end
end

(back::Jnew{T})(Δ) where T = error("Need an adjoint for constructor $T. Gradient is of type $(typeof(Δ))")


text:
#
#      Hn - The manifold of the n-dimensional sphere
#  Point is a Point on the n-dimensional sphere.
#
import LinearAlgebra: norm, dot
import Base: exp, log, show
export Hyperbolic, HnPoint, HnTVector, getValue
export distance, dot, exp, log, manifoldDimension, norm, parallelTransport
export typeofMPoint, typeofTVector, MinkowskiDot
export validateMPoint, validateTVector, zeroTVector
#
# Type definitions
#

@doc doc"""
    Hyperbolic <: Manifold
The manifold $\mathbb H^n$ is the set

```math
\mathbb H^n = \Bigl\{x\in\mathbb R^{n+1}
\ \Big|\ \langle x,x \rangle_{\mathrm{M}}= -x_{n+1}^2
+ \displaystyle\sum_{k=1}^n x_k^2 = -1, x_{n+1} > 0\Bigr\},
```

where $\langle\cdot,\cdot\rangle_{\mathrm{M}}$ denotes the [`MinkowskiDot`](@ref)
is Minkowski inner product, and this inner product in the embedded space yields
the Riemannian metric when restricted to the tangent bundle $T\mathbb H^n$.

This manifold is a matrix manifold (see [`IsMatrixM`](@ref)) and embedded (see
[`IsEmbeddedM`](@ref)).

# Abbreviation

`Hn`

# Constructor

    Hyperbolic(n)

generates the `n`-dimensional hyperbolic manifold embedded in $\mathbb R^{n+1}$.
"""
struct Hyperbolic <: Manifold
  name::String
  dimension::Int
  abbreviation::String
  Hyperbolic(dimension::Int) = new("$dimension-Hyperbolic Space",dimension,"Hn($(dimension-1))")
end
@doc doc"""
    HnPoint <: MPoint

A point $x$ on the manifold $\mathbb H^n$ represented by a vector
$x\in\mathbb R^{n+1}$ with Minkowski inner product

$\langle x,x\rangle_{\mathrm{M}} = -x_{n+1}^2 + \sum_{k=1}^n x_k^2 = -1$.
"""
struct HnPoint{T<:AbstractFloat} <: MPoint
  value::Vector{T}
  HnPoint{T}( value::Vector{T} ) where {T<:AbstractFloat} = new(value)
  HnPoint{T}( value::T) where {T<:AbstractFloat} = new([value])
end
HnPoint(value::Vector{T}) where {T<:AbstractFloat} = HnPoint{T}(value)
HnPoint(value::T) where {T <: AbstractFloat} = HnPoint{T}(value)

getValue(x::HnPoint) = length(x.value)==1 ? x.value[1] : x.value

@doc doc"""
    HnTVector <: TVector

A tangent vector $\xi \in T_x\mathbb H^n$ to a [`HnPoint`](@ref) $x$ on the
$n$-dimensional [`Hyperbolic`](@ref) space $\mathbb H^n$. To be precise
$\xi\in\mathbb R^{n+1}$ is hyperbocally orthogonal to $x\in\mathbb R^{n+1}$,
i.e. orthogonal with respect to the Minkowski inner product

$\langle \xi, x \rangle_{\mathrm{M}} = -\xi_{n+1}x_{n+1} + \sum_{k=1}^n \xi_k x_k = 0$
"""
struct HnTVector{T <: AbstractFloat}  <: TVector
    value::Vector{T}
    HnTVector{T}(value::Vector{T})  where {T <: AbstractFloat}  = new(value)
    HnTVector{T}(value::T) where {T <: AbstractFloat}  = new([value])
end
HnTVector(value::T) where {T <: AbstractFloat} = HnTVector{T}(value)
HnTVector(value::Vector{T})  where {T <: AbstractFloat}  = HnTVector{T}(value)

getValue(ξ::HnTVector) = length(ξ.value)==1 ? ξ.value[1] : ξ.value

# Traits
# ---
# (a) Hn is a MatrixManifold
@traitimpl IsMatrixM{Hyperbolic}
@traitimpl IsMatrixP{HnPoint}
@traitimpl IsMatrixTV{HnTVector}
# (b) Hn is a MatrixManifold
@traitimpl IsEmbeddedM{Hyperbolic}
@traitimpl IsEmbeddedP{HnPoint}
@traitimpl IsEmbeddedV{HnTVector}

# Functions
# ---
@doc doc"""
    distance(M,x,y)

compute the Riemannian distance on the [`Hyperbolic`](@ref) space $\mathbb H^n$
embedded in $\mathbb R^{n+1}$ can be computed as

```math
d_{\mathbb H^n}(x,y)
= \operatorname{acosh} \bigl(-\langle x,y\rangle_{\mathrm{M}}\bigr),
```

where $\langle x,y\rangle_{\mathrm{M}} = -x_{n+1}y_{n+1} +
\displaystyle\sum_{k=1}^n x_ky_k$ denotes the [`MinkowskiDot`](@ref) Minkowski
inner product on $\mathbb R^{n+1}$.
"""
distance(M::Hyperbolic,x::HnPoint{T},y::HnPoint{T}) where {T <: AbstractFloat} = acosh(  max(1.,-MinkowskiDot(getValue(x), getValue(y)))  )

@doc doc"""
    dot(M,x,ξ,ν)

compute the Riemannian inner product for two [`HnTVector`](@ref)s `ξ` and `ν`
from $T_x\mathcal M$ of the [`Hyperbolic`](@ref) space $\mathbb H^n$ given by
$\langle \xi, \nu \rangle_{\mathrm{M}}$ the [`MinkowskiDot`](@ref) Minkowski
inner product on $\mathbb R^{n+1}$.
"""
dot(M::Hyperbolic, x::HnPoint{T}, ξ::HnTVector{T}, ν::HnTVector{T}) where {T <: AbstractFloat} = MinkowskiDot( getValue(ξ), getValue(ν) )

@doc doc"""
    exp(M,x,ξ,[t=1.0])

computes the exponential map on the [`Hyperbolic`](@ref) space $\mathbb H^n$ with
respect to the [`HnPoint`](@ref) `x` and the [`HnTVector`](@ref) `ξ`, which can
be shortened with `t` to `tξ`. The formula reads

$\exp_x\xi = \cosh(\sqrt{\langle\xi,\xi\rangle_{\mathrm{M}}})x + \operatorname{sinh}(\sqrt{\langle\xi,\xi\rangle_{\mathrm{M}}})\frac{\xi}{\sqrt{\langle\xi,\xi\rangle_{\mathrm{M}}}}.$
"""
function exp(M::Hyperbolic,x::HnPoint{T},ξ::HnTVector{T},t::Float64=1.0)  where {T <: AbstractFloat}
  len = sqrt(  max( MinkowskiDot( getValue(ξ), getValue(ξ) ), 0.)  )
  if len < eps(Float64)
  	return x
  else
  	return HnPoint( cosh(t*len)*getValue(x) + sinh(t*len)/len * getValue(ξ) )
  end
end
@doc doc"""
    log(M,x,y)

computes the logarithmic map on the [`Hyperbolic`](@ref) space $\mathbb H^n$,
i.e., the [`HnTVector`](@ref) whose corresponding
[`geodesic`](@ref) starting from [`HnPoint`](@ref) `x` reaches the
[`HnPoint`](@ref) `y` after time 1 on the [`Hyperbolic`](@ref) space
$\mathbb H^n$.
The formula reads for $x\neq y$

```math
\log_x y = d_{\mathbb H^n}(x,y)\frac{y-\langle x,y\rangle_{\mathrm{M}} x}{\lVert y-\langle x,y\rangle_{\mathrm{M}} x \rVert_2}
```
and is zero otherwise.
"""
function log(M::Hyperbolic,x::HnPoint{T},y::HnPoint{T}) where {T <: AbstractFloat}
  scp = MinkowskiDot( getValue(x), getValue(y) )
  ξvalue = getValue(y) + scp*getValue(x)
  ξnorm = sqrt(max( scp.^2-1, 0.))
  if (ξnorm > eps(Float64))
    return HnTVector( acosh(max(1.,-scp))/ξnorm * ξvalue )
  else
    return zeroTVector(M,x)
  end
end
@doc doc"""
    manifoldDimension(x)

returns the dimension of the [`Hyperbolic`](@ref) space $\mathbb H^n$, the
[`HnPoint`](@ref) `x`, itself embedded in $\mathbb R^{n+1}$, belongs to.
"""
manifoldDimension(x::HnPoint)::Integer = length( getValue(x) )-1
@doc doc"""
    manifoldDimension(M)

returns the dimension of the [`Hyperbolic`](@ref) space $\mathbb H^n$.
"""
manifoldDimension(M::Hyperbolic)::Integer = M.dimension
@doc doc"""
    norm(M,x,ξ)
Computes the norm of the [`HnTVector`](@ref) `ξ` in the tangent space
$T_x\mathcal M$ at [`HnPoint`](@ref) `x` of the
[`Hyperbolic`](@ref) space $\mathbb H^n$.
"""
norm(M::Hyperbolic, x::HnPoint{T}, ξ::HnTVector{T})  where {T <: AbstractFloat}= sqrt(max( dot(M,x,ξ,ξ), 0.) )
@doc doc"""
    parallelTransport(M,x,y,ξ)

Compute the paralllel transport of the [`HnTVector`](@ref) `ξ` from
the tangent space $T_x\mathcal M$ at [`HnPoint`](@ref) `x` to
$T_y\mathcal M$ at [`HnPoint`](@ref) `y` on the
[`Hyperbolic`](@ref) space $\mathbb H^n$ along the unique [`geodesic`](@ref)
$g(\cdot;x,y)$.
The formula reads

$P_{x\to y}(\xi) = \xi - \frac{\langle \log_xy,\xi\rangle_x}
{d^2_{\mathbb H^n}(x,y)}\bigl(\log_xy + \log_yx \bigr).$
"""
function parallelTransport(M::Hyperbolic, x::HnPoint{T}, y::HnPoint{T}, ξ::HnTVector{T})  where {T <: AbstractFloat}
  ν = log(M,x,y)
  νL = norm(M,x,ν)
	return ξ - ( νL > 0 ? dot(M,x,ν,ξ)*(ν + log(M,y,x))/νL^2 : zeroTVector(M,x) )
end

typeofTVector(::Type{HnPoint{T}}) where T = HnTVector{T}
typeofMPoint(::Type{HnTVector{T}}) where T = HnPoint{T}

@doc doc"""
    typicalDistance(M)
returns the typical distance on the [`Hyperbolic`](@ref) space `M`: $\sqrt{n}$.
"""
typicalDistance(M::Hyperbolic) = sqrt(M.dimension);
@doc doc"""
    validateMPoint(M,x)

validate, that the [`HnPoint`](@ref) `x` is a valid point on the
[`Hyperbolic`](@ref) space `M`, i.e. that the dimension of $x\in\mathbb H^n$ is
correct and that its [`MinkowskiDot`](@ref) inner product is $\langle x,x\rangle_{\mathrm{M}} = -1$.
"""
function validateMPoint(M::Hyperbolic, x::HnPoint)
    if length(getValue(x)) ≠ M.dimension+1
    throw( ErrorException(
      "The Point $x is not on the $(M.name), since the vector dimension ($(length(getValue(x)))) is not $(M.dimension+1)."
    ))
  end
  if abs(MinkowskiDot(getValue(x),getValue(x))+1) >= 10^(-15)
    throw( ErrorException(
      "The Point $x is not on the $(M.name) since its minkowski inner product <x,x>_MN is $(MinkowskiDot(getValue(x),getValue(x))) and not -1"
    ))
  end
  return true
end
@doc doc"""
    validateTVector(M,x,ξ)

check that the [`HnTVector`](@ref) `ξ` is a valid tangent vector in the tangent
space of the [`HnPoint`](@ref) `x` on the [`Hyperbolic`](@ref) space `M`, i.e. `x`
is a valid point on `M`, the vectors within `ξ` and `x` agree in length and the
Minkowski inner product [`MinkowskiDot`](@ref)`(x,ξ) `is zero.
"""
function validateTVector(M::Hyperbolic, x::HnPoint, ξ::HnTVector)
    if !validateMPoint(M,x)
        return false
     end
    if length(getValue(x)) != length(getValue(ξ))
        throw( ErrorException(
            "The lengths of x ($(length(getValue(x)))) and ξ ($(length(getValue(ξ)))) do not agree, so ξ can not be a tangent vector of x."
        ))
    end
    if abs( MinkowskiDot(getValue(x),getValue(ξ)) ) > 10^(-15)
        throw( ErrorException(
            "The tangent vector should be (hyperbolically) orthogonal to its base, but the (Minkowski) inner product yields $(  MinkowskiDot(getValue(x),getValue(ξ)) )."
        ))
    end
    return true
end
@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`HnPoint`](@ref) $x\in\mathbb H^n$ on the [`Hyperbolic`](@ref) space `M`.
"""
zeroTVector(M::Hyperbolic, x::HnPoint{T}) where {T <: AbstractFloat} = HnTVector(  zero( getValue(x) )  );
# Display
# ---
show(io::IO, M::Hyperbolic) = print(io, "The $(M.name).")
show(io::IO, p::HnPoint) = print(io, "Hn($( getValue(p) ))")
show(io::IO, ξ::HnTVector) = print(io, "HnT($( getValue(ξ) ))")
# Helper
#
@doc doc"""
    MinkowskiDot(a,b)
computes the Minkowski inner product of two Vectors `a` and `b` of same length
`n+1`, i.e.

$\langle a,b\rangle_{\mathrm{M}} = -a_{n+1}b_{n+1} +
\displaystyle\sum_{k=1}^n a_kb_k.$
"""
MinkowskiDot(a::Vector,b::Vector) = -a[end]*b[end] + sum( a[1:end-1].*b[1:end-1] )


text:
function permutahedrontest(lib::Polyhedra.Library)
    A = [1 1 1; 1 0 0; 0 1 0; 0 0 1; -1 0 0; 0 -1 0; 0 0 -1]
    b = [6, 3, 3, 3, -1, -1, -1]
    linset = BitSet([1])
    V = [2 3 1; 1 3 2; 3 1 2; 3 2 1; 2 1 3; 1 2 3]
    ine = hrep(A, b, linset)
    poly = polyhedron(ine, lib)
    @test !isempty(poly)
    inequality_fulltest(poly, A, b, linset)
    generator_fulltest(poly, V, Matrix{Int}(undef, 0, 3))

    isa(lib, Polyhedra.Library) && return
    # x1___x4____________1
    #      |         |
    #      V         V
    # x2___x5___x6_______2
    #           |
    #           V
    # x3_________________3
    Alift = [-1  0  0  1  0  0;
              0 -1  0  1  0  0;
             -1 -1  0  1  1  0;
              1  1  0 -1 -1  0;
              0  0 -1  0  0  1;
              0  0  0  0 -1  1;
              0  0 -1  0 -1  1;
              0  0  1  0  1 -1;
              0  0  0 -1  0  0;
              0  0  0  0  0 -1;
              0  0  0 -1  0 -1;
              0  0  0  1  0  1]
    blift = [0; 0; 0; 0; 0; 0; -3; 3; -1; -1; -(1+2); (1+2)]
    linsetlift = BitSet()
    inelift = hrep(Alift, blift, linsetlift)
    polylift2 = polyhedron(inelift, lib)
    #polylift1 = project(polylift2, 1:4)
    polylift1 = eliminate(polylift2, BitSet([5, 6]))
    polylift0 = eliminate(polylift1)
    inequality_fulltest(polylift0, A, b, linset)
    generator_fulltest(polylift0, V)

    # removeredundantinequalities!(polylift0)
    # inelift0d = Polyhedra.HRepresentation{Int}(Polyhedra.Representation(inelift0))
    # inelift0df = Polyhedra.HRepresentation{Int}(Polyhedra.Representation(inelift0f))
    # @test inelift0d.linset == BitSet([1])
    # @test length(inelift0d.b) == 7
    # @test inelift0d.b[1] / sign(inelift0d.b[1]) == 6
    # @test vec(Array{Int}(inelift0d.A[1,:] / sign(inelift0d.b[1]))) == [1; 1; 1] # Array{Int} cast and vec are for julia 0.4
    # @test inelift0df.linset == BitSet([1])
    # @test length(inelift0df.b) == 7
    # @test inelift0df.b[1] / sign(inelift0df.b[1]) == 6
    # @test vec(Array{Int}(inelift0df.A[1,:] / sign(inelift0df.b[1]))) == [1; 1; 1] # Array{Int} cast and vec are for julia 0.4
    # polylift = CDDPolyhedra(inelift0)
    # polyliftf = CDDPolyhedra(inelift0f)
    # extunlift = Polyhedra.VRepresentation{Int}(Polyhedra.Representation(copygenerators(polylift)))
    # extunliftf = Polyhedra.VRepresentation{Int}(round(Polyhedra.Representation(copygenerators(polyliftf))))
    # generator_fulltest(extunlift, V, Array(Int, 0, 3))
    # generator_fulltest(extunliftf, V, Array(Int, 0, 3))
end


text:
#######################################################################
# Exceptions, Errors, Error processing
#######################################################################

abstract type PkgDownloaderException <: Exception end

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

max_stars_repo_name: julbinb/julia-pkgsrc-downloader

max_issues_repo_name: julbinb/julia-pkgsrc-downloader

max_forks_repo_name: julbinb/julia-pkgsrc-downloader


text:
to_func_name(x::Symbol) = string(x) |> lowercase |> Symbol

"""

# type recipe are really simple and just overload the argument conversion pipeline.
convert_arguments(P::Type{<:AbstractPlot}, x::MyType) = convert_arguments(P, rand(10, 10))
# only apply this for a certain plot type:
convert_arguments(P::Type{<:Scatter}, x::MyType) = convert_arguments(P, rand(10, 10))
# optionally define the plotting type, when someone just says plot(x::MyType)
# - will fall back to whatever the standard plot type for what you return from argument_convert is!
plottype(::MyType) = Surface

# full recipes:

# (x, y, z) && themes are optional
@recipe(MyPlot, x, y, z) do scene
    Theme(
        plot_color => :red
    )
end
# this macro defines the following:

const Myplot{ArgTypes} = Combined{myplot, ArgTypes}
# the type parameter contains the function instead of e.g. a symbol. This way
# the mapping from Myplot to myplot is safer and simpler.
# Downside, we always need a function. Would there be rejection against that?

# all the signatures to make it nice to use:
myplot(args...; kw_args...) = ...
myplot!(scene, args...; kw_args...) = ...
myplot(kw_args::Dict, args...) = ...
myplot!(scene, kw_args::Dict, args...) = ...
#etc (not 100% settled what signatures there will be)

# the next one is optional, but it will allow to say plot_object[:x] to fetch the first argument
# from the call `plot_object = myplot(rand(10), rand(10), rand(10))`.
# If you leave out the (x, y, z)
# in the recipe macro, it will default to plot_object[1] (always works), and plot_object[:arg1]
argument_names(::Type{<: MyPlot}) = (:x, :y, :z)

# this function will insert the theme into any scene that plots Myplot.
function default_theme(scene, ::Myplot)
    Theme(
        plot_color => :red
    )
end
#-----------------------------
# Implement the recipe
function plot!(plot::MyPlot)
    # normal plotting code, building up on any previously defined recipes
    # or atomic plotting operations, and adds it to the combined `plot`:
    lines!(plot, rand(10), color = plot[:plot_color])
    plot!(plot, plot[:x], plot[:y])
    plot
end

# Add specialization
const MyVolume = MyPlot{Tuple{ArgTypes <: AbstractArray{<: AbstractFloat, 3}}}
argument_names(::Type{<: MyVolume}) = (:volume,) # again, optional
function plot!(plot::MyVolume)
    # plot a volume with a colormap going from fully transparent to plot_color
    volume!(plot, plot[:volume], colormap = :transparent => plot[:plot_color])
    plot
end
"""
macro recipe(theme_func, Tsym::Symbol, args::Symbol...)
    funcname = to_func_name(Tsym)
    funcname! = esc(Symbol("$(funcname)!"))
    T = esc(Tsym)
    funcname = esc(funcname)
    expr = quote
        $funcname() = not_implemented_for($funcname)
        if !@isdefined($(Tsym))# make this work with interactive usage
            const $(T){$(esc(:ArgType))} = Combined{$funcname, $(esc(:ArgType))}
        end
        Base.show(io::IO, ::Type{<: $T}) = print(io, $(string(Tsym)), "{...}")
        $funcname(args...; attributes...) = plot($T, args...; attributes...)
        $funcname!(scene::SceneLike, args...; attributes...) = plot!(scene, $T, args...; attributes...)
        $funcname(attributes::Attributes, args...) = plot($T, attributes, args...)
        $funcname!(scene::SceneLike, attributes::Attributes, args...) = plot!(scene, $T, attributes, args...)
        $funcname!(attributes::Attributes, args...) = plot!(current_scene(), $T, attributes, args...)
        $funcname!(args...; kw_args...) = plot!(current_scene(), $T, Attributes(kw_args), args...)

        Base.@__doc__($funcname)
        AbstractPlotting.default_theme(scene, ::Type{<: $T}) = $(esc(theme_func))(scene)

        export $T, $funcname, $funcname!
    end
    if !isempty(args)
        # quoted_args = map(QuoteNode, args)
        push!(expr.args, :($(esc(:argument_names))(::Type{<: $T}, len::Integer) = $args))
    end
    expr
end

macro atomic(theme_func, Tsym::Symbol)
    funcname = to_func_name(Tsym)
    funcname! = esc(Symbol("$(funcname)!"))
    T = esc(Tsym)
    funcname = esc(funcname)
    quote
        $funcname() = not_implemented_for($funcname)
        Base.@__doc__($funcname)
        const $T{$(esc(:ArgType))} = Atomic{$funcname, $(esc(:ArgType))}
        Base.show(io::IO, ::Type{<: $T}) = print(io, $(string(Tsym)), "{...}")
        $funcname(args...; attributes...) = plot($T, args...; attributes...)
        $funcname!(scene::SceneLike, args...; attributes...) = plot!(scene, $T, args...; attributes...)
        $funcname(attributes::Attributes, args...) = plot($T, attributes, args...)
        $funcname!(scene::SceneLike, attributes::Attributes, args...) = plot!(scene, $T, attributes, args...)

        $funcname!(attributes::Attributes, args...) = plot!(current_scene(), $T, attributes, args...)
        $funcname!(args...; kw_args...) = plot!(current_scene(), $T, Attributes(kw_args), args...)

        AbstractPlotting.default_theme(scene, ::Type{<: $T}) = $theme_func(scene)

        export $T, $funcname, $funcname!
    end
end


text:
module CipherModes

import Base
importall BlockCiphers
importall Padding
importall Iteration
importall Rand

function make_output_array(padding, extra, key, in)
    out = zeros(Uint8, extra + padded_size(block_size(key), in))
    out[(1 + extra):(length(in) + extra)] = in
    pad!(padding, block_size(key), out, length(in) + extra)
    out
end

# --- ecb mode
# warning -- exposed for amusement purposes only, DO NOT USE
function ecb_encrypt!(padding, key, in :: Array{Uint8})
    for block = slicerator(in, block_size(key))
        encrypt!(key, block)
    end
    in
end
ecb_encrypt(padding, key, in :: Array{Uint8}) = ecb_encrypt!(padding, key, make_output_array(padding, 0, key, in))

function ecb_decrypt!(padding, key, in :: Array{Uint8})
    for block = slicerator(in, block_size(key))
        decrypt!(key, block)
    end
    unpad!(padding, in)
    in
end
ecb_decrypt(padding, key, in :: Array{Uint8}) = ecb_decrypt!(padding, key, copy(in))

function xor!(a, b)
    for i=1:length(a)
        a[i] $= b[i]
    end
    a
end

# --- cbc mode
function cbc_encrypt!(padding, key, in :: Array{Uint8})
    in[1:block_size(key)] = cryptographic_rand(block_size(key))

    previous = in[1:block_size(key)]
    for in_block = offset(1, slicerator(in, block_size(key)))
        xor!(in_block, previous)
        encrypt!(key, in_block)
        previous[1:block_size(key)] = in_block[1:block_size(key)]
    end
    in
end
cbc_encrypt(padding, key, in) = cbc_encrypt!(padding, key, make_output_array(padding, block_size(key), key, in))

function cbc_decrypt!(padding, key, in :: Array{Uint8})
    iv = in[1:block_size(key)]
    for i=1:block_size(key)
        shift!(in)
    end

    tmp = zeros(Uint8, block_size(key))
    previous = iv
    for block = slicerator(in, block_size(key))
        tmp[1:block_size(key)] = block
        decrypt!(key, block, block)
        xor!(block, previous)
        previous[1:block_size(key)] = tmp
    end
    unpad!(padding, in)
    in
end
cbc_decrypt(padding, key, in) = cbc_decrypt!(padding, key, copy(in))

# --- pcb mode
function pcbc_encrypt!(padding, key, in :: Array{Uint8})
    in[1:block_size(key)] = cryptographic_rand(block_size(key))

    unencrypted = zeros(Uint8, block_size(key))
    tmp = in[1:block_size(key)]
    for block = offset(1, slicerator(in, block_size(key)))
        unencrypted[1:block_size(key)] = block
        xor!(block, tmp)
        encrypt!(key, block)
        tmp[1:block_size(key)] = unencrypted $ block
    end
    in
end
pcbc_encrypt(padding, key, in) = pcbc_encrypt!(padding, key, make_output_array(padding, block_size(key), key, in))

function pcbc_decrypt!(padding, key, in :: Array{Uint8})
    iv = in[1:block_size(key)]
    for i=1:block_size(key)
        shift!(in)
    end

    encrypted = zeros(Uint8, block_size(key))
    tmp = copy(iv)
    for block = slicerator(in, block_size(key))
        encrypted[1:block_size(key)] = block
        decrypt!(key, block)
        xor!(block, tmp)
        tmp[1:block_size(key)] = block $ encrypted
    end
    unpad!(padding, in)
    in
end
pcbc_decrypt(padding, key, in) = pcbc_decrypt!(padding, key, copy(in))

# --- ctr mode
function ctr_encrypt64!(key, nonce :: Int64, in :: Array{Uint8})
    seed = zeros(Uint8, block_size(key))
    seed[5:8] = [uint8(nonce << (offset * 8)) for offset=0:3]

    tmp = zeros(Uint8, block_size(key))
    for ctr = 0:int32(ceil(length(in) / block_size(key)))
        tmp[1:8] = seed
        tmp[1:4] = [uint8(ctr << (offset * 8)) for offset=0:3]

        encrypt!(key, tmp)
        xor!(slice(in, ((ctr * block_size(key)) + 1) : min(length(in), ((ctr + 1) * block_size(key)))), tmp)
    end
    in
end

function ctr_encrypt128!(key, nonce :: Int64, in :: Array{Uint8})
    seed = zeros(Uint8, block_size(key))
    seed[9:16] = [uint8(nonce << (offset * 8)) for offset=0:7]

    tmp = zeros(Uint8, block_size(key))
    for ctr = 0:int64(ceil(length(in) / block_size(key)))
        tmp[1:16] = seed
        tmp[1:8] = [uint8(ctr << (offset * 8)) for offset=0:7]

        encrypt!(key, tmp)
        xor!(slice(in, ((ctr * block_size(key)) + 1) : min(length(in), ((ctr + 1) * block_size(key)))), tmp)
    end
    in
end
function ctr_encrypt!(key, nonce :: Int64, in :: Array{Uint8})
    if(block_size(key) == 8)
        ctr_encrypt64!(key, nonce, in)
    else
        ctr_encrypt128!(key, nonce, in)
    end
end
ctr_encrypt(key, nonce :: Int64, in :: Array{Uint8}) = ctr_encrypt!(key, nonce, copy(in))

ctr_decrypt!(key, nonce :: Int64, in :: Array{Uint8}) = ctr_encrypt!(key, nonce, in)
ctr_decrypt(key, nonce :: Int64, in :: Array{Uint8}) = ctr_encrypt(key, nonce, in)

# --- cfb mode
function cfb_encrypt!(key, in :: Array{Uint8})
    in[1:block_size(key)] = cryptographic_rand(block_size(key))

    tmp = in[1:block_size(key)]
    for idx = 1:int64(floor(length(in) / block_size(key)))
        encrypt!(key, tmp)
        xor!(slice(in, ((idx * block_size(key)) + 1) : min(length(in), ((idx + 1) * block_size(key)))), tmp)
        tmp = slice(in, ((idx * block_size(key)) + 1) : min(length(in), ((idx + 1) * block_size(key))))
    end
    in
end
function cfb_encrypt(key, in)
    out = zeros(Uint8, length(in) + block_size(key))
    cfb_encrypt!(key, out)
end

function cfb_decrypt!(key, in :: Array{Uint8})
    iv = in[1:block_size(key)]
    for i=1:block_size(key)
        shift!(in)
    end

    tmp = iv
    for idx = 0:int64(floor(length(in) / block_size(key)))
        encrypt!(key, tmp)
        tmp = copy(slice(in, ((idx * block_size(key)) + 1) : min(length(in), ((idx + 1) * block_size(key)))))
        xor!(slice(in, ((idx * block_size(key)) + 1) : min(length(in), ((idx + 1) * block_size(key)))), tmp)
    end
    in
end
cfb_decrypt(key, in) = cfb_decrypt!(key, copy(in))

# --- ofb mode
function ofb_encrypt!(key, in :: Array{Uint8})
    in[1:block_size(key)] = cryptographic_rand(block_size(key))

    tmp = in[1:block_size(key)]
    for idx = 1:int64(floor(length(in) / block_size(key)))
        encrypt!(key, tmp)
        xor!(slice(in, ((idx * block_size(key)) + 1) : min(length(in), ((idx + 1) * block_size(key)))), tmp)
    end
    in
end
function ofb_encrypt(key, in)
    out = zeros(Uint8, length(in) + block_size(key))
    ofb_encrypt!(key, out)
end

function ofb_decrypt!(key, in :: Array{Uint8})
    iv = in[1:block_size(key)]
    for i=1:block_size(key)
        shift!(in)
    end

    tmp = iv
    for idx = 0:int64(floor(length(in) / block_size(key)))
        encrypt!(key, tmp)
        xor!(slice(in, ((idx * block_size(key)) + 1) : min(length(in), ((idx + 1) * block_size(key)))), tmp)
    end
    in
end
ofb_decrypt(key, in) = ofb_decrypt!(key, copy(in))

export ecb_encrypt!, ecb_encrypt, ecb_decrypt!, ecb_decrypt,
       cbc_encrypt, cbc_decrypt,
       pcbc_encrypt, pcbc_decrypt,
       ctr_encrypt, ctr_decrypt,
       cfb_encrypt, cfb_decrypt,
       ofb_encrypt, ofb_decrypt

end # module CipherModes


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

"""
    ConjArray(array)

A lazy-view wrapper of an `AbstractArray`, taking the elementwise complex conjugate. This
type is usually constructed (and unwrapped) via the [`conj`](@ref) function (or related
[`ctranspose`](@ref)), but currently this is the default behavior for `RowVector` only. For
other arrays, the `ConjArray` constructor can be used directly.

# Examples

```jldoctest
julia> [1+im, 1-im]'
1×2 RowVector{Complex{Int64},ConjArray{Complex{Int64},1,Array{Complex{Int64},1}}}:
 1-1im  1+1im

julia> ConjArray([1+im 0; 0 1-im])
2×2 ConjArray{Complex{Int64},2,Array{Complex{Int64},2}}:
 1-1im  0+0im
 0+0im  1+1im
```
"""
immutable ConjArray{T, N, A <: AbstractArray} <: AbstractArray{T, N}
    parent::A
end

@inline ConjArray{T,N}(a::AbstractArray{T,N}) = ConjArray{conj_type(T), N, typeof(a)}(a)

ConjVector{T, V <: AbstractVector} = ConjArray{T, 1, V}
@inline ConjVector{T}(v::AbstractVector{T}) = ConjArray{conj_type(T), 1, typeof(v)}(v)

ConjMatrix{T, M <: AbstractMatrix} = ConjArray{T, 2, M}
@inline ConjMatrix{T}(m::AbstractMatrix{T}) = ConjArray{conj_type(T), 2, typeof(m)}(m)

# This type can cause the element type to change under conjugation - e.g. an array of complex arrays.
@inline conj_type(x) = conj_type(typeof(x))
@inline conj_type{T}(::Type{T}) = promote_op(conj, T)

@inline parent(c::ConjArray) = c.parent
@inline parent_type(c::ConjArray) = parent_type(typeof(c))
@inline parent_type{T,N,A}(::Type{ConjArray{T,N,A}}) = A

@inline size(a::ConjArray) = size(a.parent)
IndexStyle{CA <: ConjArray}(::CA) = IndexStyle(parent_type(CA))
IndexStyle{CA <: ConjArray}(::Type{CA}) = IndexStyle(parent_type(CA))

@propagate_inbounds getindex{T,N}(a::ConjArray{T,N}, i::Int) = conj(getindex(a.parent, i))
@propagate_inbounds getindex{T,N}(a::ConjArray{T,N}, i::Vararg{Int,N}) = conj(getindex(a.parent, i...))
@propagate_inbounds setindex!{T,N}(a::ConjArray{T,N}, v, i::Int) = setindex!(a.parent, conj(v), i)
@propagate_inbounds setindex!{T,N}(a::ConjArray{T,N}, v, i::Vararg{Int,N}) = setindex!(a.parent, conj(v), i...)

@inline similar{T,N}(a::ConjArray, ::Type{T}, dims::Dims{N}) = similar(parent(a), T, dims)

# Currently, this is default behavior for RowVector only
@inline conj(a::ConjArray) = parent(a)

# Helper functions, currently used by RowVector
@inline _conj(a::AbstractArray) = ConjArray(a)
@inline _conj{T<:Real}(a::AbstractArray{T}) = a
@inline _conj(a::ConjArray) = parent(a)
@inline _conj{T<:Real}(a::ConjArray{T}) = parent(a)


text:
using BinDeps
@BinDeps.setup

visa = library_dependency("visa", aliases = ["visa64","VISA","/Library/Frameworks/VISA.framework/VISA"])

@BinDeps.install Dict(:visa=>:libvisa)


text:
## To do:
# * exceptions: no file, missing fam/bim, out-of-bounds indexing, not a bed file, old bed file
# * working directory

include("../src/BEDMatrices.jl")
using BEDMatrices
using Base.Test
include("testtools.jl")

const examplepath = "./data/"
const bedfile = "example.bed"
const rawfile = "example.raw"

const singleindices = [(13, 1), (13, 2), (13, 3), (13, 4),
                       (50, 998), (50, 999), (50, 1000),
                       (47, 3), (48, 3), (49, 3), (50, 3)]
const rowarrays = [Vector{Int}(0), collect(1:50), [1, 11, 12], [1, 2], [48, 7, 49, 50]]
const colarrays = [Vector{Int}(0), collect(1:1000), [491, 13, 73], [1], [33, 31]]
const rowranges = [1:0, 1:2, 1:5, 1:50, 2:50, 3:50, 4:50]
const colranges = [1:0, 1:2, 1:101, 990:1000]
const linearranges = [1:75, 2:75, 50:51, 49:101, 40_000:50_000]
const rowlogicals = [fill(true, 50), fill(false, 50), map(x -> x > 25, 1:50)]
const collogicals = [fill(true, 1000), fill(false, 1000), map(x -> x > 777, 1:1000)]

const examplerows, examplecols, exampledata = readRAW(examplepath*rawfile)
const bed = BEDMatrix(examplepath*bedfile)
const simbeds = [simulatedBEDMatrix(4, 5), simulatedBEDMatrix(5, 5, Int, 255), simulatedBEDMatrix(6, 5, Float64),
                 simulatedBEDMatrix(7, 5, Int8, -128), simulatedBEDMatrix(1, 5, Int8, -128, false), simulatedBEDMatrix(2, 1, Int8, -128, false),
                 simulatedBEDMatrix(3, 1, UInt, 0b11, false), simulatedBEDMatrix(4, 1, Int8, -128, false),
                 simulatedBEDMatrix(13, 11, Int8, 127, true, true), simulatedBEDMatrix(13, 11, Float32, -128.0, false, true)]

@testset "Basic Functional Testing" begin
    @testset "Reading BED file directly into matrix" begin
        @test BEDintomatrix(examplepath*bedfile) == exampledata
    end

    @testset "BEDMatrix" begin
        @test bed[:, :] == exampledata

        @test path(bed) == abspath(examplepath*bedfile)

        @test rownames(bed) == examplerows
        @test colnames(bed) == examplecols
    end

    @testset "simulated BEDMatrix tests" begin
        # Mostly to check for possible mod 4 or type issues
        for (simdata, simbed) in simbeds
            @test simbed[:, :] == simdata
            for j in eachindex(simdata)
                @test simbed[j] === simdata[j]
            end
        end
    end
end

@testset "Numeric Indexing" begin
    @testset "single indices" begin
        for (r, c) in singleindices
            @test bed[r, c] == exampledata[r, c]
            @test bed[50*(c-1) + r] == exampledata[50*(c-1) + r]
        end
    end
    @testset "row ranges" begin
        for c in 1:1000
            for rrange in rowranges
                @test bed[rrange, c] == exampledata[rrange, c]
            end
            @test bed[4:end, c] == exampledata[4:end, c]
            @test bed[:, c] == exampledata[:, c]
        end
    end
    @testset "column ranges" begin
        for crange in colranges
            for r in 1:50
                @test bed[r, crange] == exampledata[r, crange]
            end
        end
        @test bed[17, 2:end] == exampledata[17, 2:end]
        @test bed[19, :] == exampledata[19, :]
    end
    @testset "linear ranges" begin
        for lrange in linearranges
            @test bed[lrange] == exampledata[lrange]
        end
    end
    @testset "column arrays" begin
        for cvec in colarrays
            @test bed[11, cvec] == exampledata[11, cvec]
        end
    end
    @testset "row arrays" begin
        for rvec in rowarrays
            @test bed[rvec, 11] == exampledata[rvec, 11]
        end
    end
    @testset "row arrays with column ranges" begin
        for crange in colranges
            for rvec in rowarrays
                @test bed[rvec, crange] == exampledata[rvec, crange]
            end
        end
    end
end

@testset "Name Indexing" begin
    @testset "single indices" begin
        for (r, c) in singleindices
            rname = examplerows[r]
            cname = examplecols[c]
            @test bed[rname, c] == exampledata[r, c]
            @test bed[r, cname] == exampledata[r, c]
            @test bed[rname, cname] == exampledata[r, c]
            @test bed[rname, :] == convert(Vector{eltype(bed)}, exampledata[r, :])
            @test bed[:, cname] == convert(Vector{eltype(bed)}, exampledata[:, c])
        end
    end

    @testset "arrays" begin
        for cvec in colarrays
            cnames = [examplecols[c] for c in cvec]
            @test bed[11, cnames] == convert(Vector{eltype(bed)}, exampledata[11, cvec])
            @test bed[:, cnames] == convert(Matrix{eltype(bed)}, exampledata[:, cvec])
            @test bed[2:end, cnames] == convert(Matrix{eltype(bed)}, exampledata[2:end, cvec])

            for rvec in rowarrays
                rnames = [examplerows[r] for r in rvec]
                @test bed[rnames, cnames] == exampledata[rvec, cvec]
            end
        end
        for rvec in rowarrays
            rnames = [examplerows[r] for r in rvec]
            @test bed[rnames, 11] == convert(Vector{eltype(bed)}, exampledata[rvec, 11])
            @test bed[rnames, :] == convert(Matrix{eltype(bed)}, exampledata[rvec, :])
            @test bed[rnames, 2:end] == convert(Matrix{eltype(bed)}, exampledata[rvec, 2:end])
        end
    end
end

@testset "Logical Indexing" begin
    for (r, c) in singleindices
        for clogic in collogicals
            @test bed[r, clogic] == exampledata[r, clogic]
        end
        for rlogic in rowlogicals
            @test bed[rlogic, c] == exampledata[rlogic, c]
        end
    end
    for clogic in collogicals
        for rlogic in rowlogicals
            @test bed[rlogic, clogic] == exampledata[rlogic, clogic]
        end
    end
end

# @testset "tocontiguous" begin
#     for rrange in rowranges
#         rows = collect(rrange)
#         if isempty(rows)
#             @test BEDMatrices.tocontiguous(rows) == Vector{Int}(0)
#         else
#             @test BEDMatrices.tocontiguous(rows)[1] == rrange
#         end
#     end
#
#     for rows in rowarrays[2:end]
#         @test mapreduce(collect, vcat, BEDMatrices.tocontiguous(rows)) == rows
#     end
#
#     for rlogic in rowlogicals
#         if any(rlogic)
#             @test [x in mapreduce(collect, vcat, BEDMatrices.tocontiguous(rlogic)) for x in 1:50] == rlogic
#         else
#             @test BEDMatrices.tocontiguous(rlogic) == []
#         end
#     end
# end

@testset "Column Tools" begin
    @testset "typealiases" begin
        m = zeros(10, 20)

        v = view(bed, :, 1)
        w = view(bed, 1:25, 1)
        x = view(bed, [1, 2, 5], 1)

        u = view(m, :, 1)
        t = view(m, 1:4, 1)
        s = view(m, [1, 2, 5], 1)

        @test isa(v, BEDMatrices.Column) == true
        @test isa(v, BEDMatrices.ColumnUnitRange) == true
        @test isa(v, BEDMatrices.BEDColumn) == true
        @test isa(v, BEDMatrices.BEDColumnUnitRange) == true
        @test isa(v, BEDMatrices.BEDSubColumn) == true

        @test isa(w, BEDMatrices.Column) == false
        @test isa(w, BEDMatrices.ColumnUnitRange) == true
        @test isa(w, BEDMatrices.BEDColumn) == false
        @test isa(w, BEDMatrices.BEDColumnUnitRange) == true
        @test isa(w, BEDMatrices.BEDSubColumn) == true

        @test isa(x, BEDMatrices.Column) == false
        @test isa(x, BEDMatrices.ColumnUnitRange) == false
        @test isa(x, BEDMatrices.BEDColumn) == false
        @test isa(x, BEDMatrices.BEDColumnUnitRange) == false
        @test isa(x, BEDMatrices.BEDSubColumn) == true

        @test isa(u, BEDMatrices.Column) == true
        @test isa(u, BEDMatrices.ColumnUnitRange) == true
        @test isa(u, BEDMatrices.BEDColumn) == false
        @test isa(u, BEDMatrices.BEDColumnUnitRange) == false
        @test isa(u, BEDMatrices.BEDSubColumn) == false

        @test isa(t, BEDMatrices.Column) == false
        @test isa(t, BEDMatrices.ColumnUnitRange) == true
        @test isa(t, BEDMatrices.BEDColumn) == false
        @test isa(t, BEDMatrices.BEDColumnUnitRange) == false
        @test isa(t, BEDMatrices.BEDSubColumn) == false

        @test isa(s, BEDMatrices.Column) == false
        @test isa(s, BEDMatrices.ColumnUnitRange) == false
        @test isa(s, BEDMatrices.BEDColumn) == false
        @test isa(s, BEDMatrices.BEDColumnUnitRange) == false
        @test isa(s, BEDMatrices.BEDSubColumn) == false
    end

    @testset "hasNAs" begin
        # extra test to confirm that simulatedBEDMatrix is working correctly
        nomissingdata, nomissingbed = simulatedBEDMatrix(50, 1000, Int, -128, false)
        @test any(e -> e == -128, nomissingdata) == false
        @test hasNAs(nomissingbed) == false

        @test hasNAs(bed) == true
        for (simdata, simbed) in simbeds
            @test hasNAs(simbed) == any(e -> e === NArep(simbed), simdata)
        end

        for col in 1:20
            @test hasNAs(view(bed, :, col)) == any(e -> e === BEDMatrices.NA_byte, exampledata[:, col])
            @test hasNAs(view(bed, 1:25, col)) == any(e -> e === BEDMatrices.NA_byte, exampledata[1:25, col])
        end

        @test hasNAs(view(bed, 1:6, 9)) == any(e -> e === BEDMatrices.NA_byte, exampledata[1:6, 9])
        @test hasNAs(view(bed, 1:7, 9)) == any(e -> e === BEDMatrices.NA_byte, exampledata[1:7, 9])
        @test hasNAs(view(bed, 8:50, 9)) == any(e -> e === BEDMatrices.NA_byte, exampledata[8:end, 9])
        @test hasNAs(view(bed, 7:50, 9)) == any(e -> e === BEDMatrices.NA_byte, exampledata[7:end, 9])

        @test hasNAs(view(bed, 1:6, 7:10)) == any(e -> e === BEDMatrices.NA_byte, exampledata[1:6, 7:10])
        @test hasNAs(view(bed, 7:10, 1:20)) == any(e -> e === BEDMatrices.NA_byte, exampledata[7:10, 1:20])
    end

    @testset "countNAs" begin
        @test countNAs(bed) == count(e -> e === BEDMatrices.NA_byte, exampledata)

        for col in 1:25
            @test countNAs(view(bed, :, col)) == count(e -> e === BEDMatrices.NA_byte, exampledata[:, col])
            @test countNAs(view(bed, 1:25, col)) == count(e -> e === BEDMatrices.NA_byte, exampledata[1:25, col])
            @test countNAs(view(bed, [2, 3, 5, 7, 11, 13], col)) == count(e -> e === BEDMatrices.NA_byte, exampledata[[2, 3, 5, 7, 11, 13], col])
        end

        for (simdata, simbed) in simbeds
            @test countNAs(simbed) == count(e -> e === NArep(simbed), simdata)
            @test countNAs(view(simbed, :, 1)) == count(e -> e === NArep(simbed), simdata[:, 1])
        end
    end

    @testset "sum" begin
        for col in 1:20
            @test sum(view(bed, :, col)) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col])
            @test sum(view(bed, 1:25, col)) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[1:25, col])
            @test sum(view(bed, [2, 3, 5, 7, 11], col)) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[[2, 3, 5, 7, 11], col])

            @test sum(x -> (x-1)^10, view(bed, :, col)) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : (e - 1)^10, exampledata[:, col])
            @test sum(x -> x === NArep(bed) ? 37 : (x-1)^2, view(bed, :, col); skipna=false) == sum(e -> e === BEDMatrices.NA_byte ? 37 : (e - 1)^2, exampledata[:, col])
        end
    end

    @testset "column_sum(func, ...)" begin
        for col in 1:20
            @test column_sum(abs2, bed, col) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : abs2(e), exampledata[:, col])
            @test column_sum(x -> (x-1)^10, bed, col, 1:25) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : (e-1)^10, exampledata[1:25, col])
        end
    end

    @testset "BEDdot" begin
        v = collect(1:50)
        for col1 in 20:30
            @test BEDdot(view(bed, :, col1), v) == dot(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col1]), v)
            @test BEDdot(v, view(bed, :, col1)) == dot(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col1]), v)
            for col2 in 31:36
                @test BEDdot(view(bed, :, col1), view(bed, :, col2)) == dot(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col1]),
                                                                            map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col2]))
            end
        end
    end

    @testset "column_dot" begin
        v = collect(1:50)
        for col1 in 20:30
            @test column_dot(bed, col1, v) == dot(map(e -> e === NArep(bed) ? zero(e) : e, bed[:, col1]), v)
            for col2 in 26:30
                @test column_dot(bed, col1, col2) == dot(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col1]),
                                                         map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col2]))
            end
        end
    end

    @testset "column_dist" begin
        for col in 21:40
            dist = BEDMatrices.column_dist(bed, col)
            @test dist[1] == count(e -> e == 0, exampledata[:, col])
            @test dist[2] == count(e -> e == 1, exampledata[:, col])
            @test dist[3] == count(e -> e == 2, exampledata[:, col])
            @test dist[4] == count(e -> e === BEDMatrices.NA_byte, exampledata[:, col])
        end
    end

    @testset "column_dist_dot" begin
        v = collect(1:50)
        for col in 21:40
            distdot = BEDMatrices.column_dist_dot(bed, col, v)
            @test distdot[1] == count(e -> e == 0, exampledata[:, col])
            @test distdot[2] == count(e -> e == 1, exampledata[:, col])
            @test distdot[3] == count(e -> e == 2, exampledata[:, col])
            @test distdot[4] == count(e -> e === BEDMatrices.NA_byte, exampledata[:, col])
            @test distdot[5] == dot(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col]), v)
            @test distdot[6] == dot(map(e -> e === BEDMatrices.NA_byte ? one(e) : zero(e), exampledata[:, col]), v)
            @test distdot[7] == dot(map(e -> e === BEDMatrices.NA_byte ? one(e) : zero(e), exampledata[:, col]), v.^2)
        end
    end

    @testset "column_NAsup_dot" begin
        v = collect(1:50)
        for col in 1:20
            @test column_NAsup_dot(bed, col, v) == dot(map(e -> e === BEDMatrices.NA_byte ? one(e) : zero(e), exampledata[:, col]), v)
        end
    end

    @testset "sumabs2" begin
        for col in 1:20
            @test sumabs2(view(bed, :, col)) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : abs2(e), exampledata[:, col])
            @test sumabs2(view(bed, 25:50, col)) == sum(e -> e === BEDMatrices.NA_byte ? zero(e) : abs2(e), exampledata[25:50, col])
        end
    end

    @testset "column_norm" begin
        for col in 20:40
            @test column_norm(bed, col) == norm(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col]))
            @test column_norm(bed, col, :, 3) ≈ norm(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col]), 3)
            @test column_norm(bed, col, :, 4) ≈ norm(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[:, col]), 4)
            @test column_norm(bed, col, 1:30, 4) ≈ norm(map(e -> e === BEDMatrices.NA_byte ? zero(e) : e, exampledata[1:30, col]), 4)
        end
    end
end


text:
function MarkAllCellsNew(r)
    side = Ost
    counter = 0
    B = Array{Int, 1}(undef, 1000)
    B, counter = movetothecorner(r, West, B)

    markrow(r, side)
    putmarker!(r)
    while isborder(r, Nord) == false
        move!(r, Nord)
        side = inverse(side)
        markrow(r, side)
        putmarker!(r)
    end

    movements(r, Sud)
    movements(r, West)
    returnfromthecorner(r, inverse(West), B, counter)
end

function movements(r, side)
    while isborder(r, side) == false
        move!(r, side)
    end
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function markrow(r, side)
    putmarker!(r)
    while trymove!(r, side) == false
        putmarker!(r)
        trymove!(r, side)
        putmarker!(r)
    end
end

function movetothecorner(r, side, B)
    counter = 1
    while ((!isborder(r, side)) || (!isborder(r, HorizonSide((Int(side) + 1) % 4))))
        if !isborder(r, side)
            move!(r, side)
            B[counter] = 1
            counter += 1
        end
        if !isborder(r, HorizonSide((Int(side) + 1) % 4))
            move!(r, HorizonSide((Int(side) + 1) % 4))
            B[counter] = 2
            counter += 1
        end
    end
    return B, counter
end

function returnfromthecorner(r, side, B, counter)
    for i in 1:(counter - 1)
        if B[counter - i] == 1
            move!(r, side)
        end
        if B[counter - i] == 2
            move!(r, HorizonSide((Int(side) + 1) % 4))
        end
    end
end

function trymove!(r::Robot, side::HorizonSide)
    x = 0
    z = false
    if isborder(r,side)
        while isborder(r, side)
            if isborder(r,HorizonSide((Int(side) + 1) % 4))
                z = true
                break
            end
            move!(r,HorizonSide((Int(side) + 1) % 4))
            x = x + 1
        end
        if !z
            move!(r,side)
            while isborder(r,HorizonSide((Int(side) + 3) % 4))
                move!(r,side)
            end
        end

        while x != 0
            move!(r,HorizonSide((Int(side) + 3) % 4))
            x = x - 1
        end
        return z
    else
        move!(r,side)
        return false
    end
end


text:
using LinearAlgebra

# List Operations
x = [ 1.0, 0 ];   y = [ 1.0, -1.0 ];   z = [ 0, 1.0];

list = [ x, y, z ]              # lists 3 arrays w/elements
list[2]                         # lists elements of 2nd array
list[3]                         # lists elements of 3rd array
list = ( x, y, z )              # list of 3 arrays w/elements

# Creating vectors with Julia functions
zeros(3)                        # creates zero vector of dim 3

i = 3; n = 4;                   # sets values
ei = zeros(n);                  # creates zero vector of dim n
ei[i] = 1;                      # sets 3rd value to 1.0
ei                              # produces [ 0.0, 0.0, 1.0, 0.0]

unit_vector(i,n) = [zeros(i-1); 1 ; zeros(n-i)]
unit_vector(3,4)                # produces [ 0.0, 0.0, 1.0, 0.0]

ones(2)                         # creates one vector of dim 2
rand(2)                         # creates rand vector of dim 2

# Add or subtract vectors
[ 0, 7, 3] + [ 1, 2, 0]
[ 1, 9 ] - [ 1, 1]

# Scalar vector operations
x = [ 0, 2, -1 ];
2.2 * x
x / 3
[ 1.1, -3.7, 0.3 ] .- 1.4       # '.' applies op to each element
0.7 .+ [1,-1]                   # '.' applies op to each element

p_1= [ 22.15, 89.32, 56.77 ];
p_2 = [ 23.05, 87.32, 57.13 ];
r = (p_2 - p_1) ./ p_1          # '.' applies op to each element

w = [1,2,2]; z = [1,2,3];
w == z                          # returns false
w .== z                         # returns true, true, false

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
##### Get prognostic variable list

prognostic_vars(::TotalEnergyModel) = (Energy(),)
prognostic_vars(::θModel) = (ρθ_liq_ice(),)

prognostic_vars(::DryModel) = ()
prognostic_vars(::EquilMoist) = (TotalMoisture(),)
prognostic_vars(::NonEquilMoist) =
    (TotalMoisture(), LiquidMoisture(), IceMoisture())

prognostic_vars(::NoPrecipitation) = ()
prognostic_vars(::RainModel) = (Rain(),)
prognostic_vars(::RainSnowModel) = (Rain(), Snow())

prognostic_vars(::NoTracers) = ()
prognostic_vars(::NTracers{N}) where {N} = (Tracers{N}(),)

prognostic_vars(atmos::AtmosModel) = prognostic_vars(atmos.physics)

prognostic_vars(m::AtmosPhysics) = (
    Mass(),
    Momentum(),
    prognostic_vars(energy_model(m))...,
    prognostic_vars(moisture_model(m))...,
    prognostic_vars(precipitation_model(m))...,
    prognostic_vars(tracer_model(m))...,
    prognostic_vars(turbconv_model(m))...,
)

get_prog_state(state, ::Mass) = (state, :ρ)
get_prog_state(state, ::Momentum) = (state, :ρu)
get_prog_state(state, ::Energy) = (state.energy, :ρe)
get_prog_state(state, ::ρθ_liq_ice) = (state.energy, :ρθ_liq_ice)
get_prog_state(state, ::TotalMoisture) = (state.moisture, :ρq_tot)
get_prog_state(state, ::LiquidMoisture) = (state.moisture, :ρq_liq)
get_prog_state(state, ::IceMoisture) = (state.moisture, :ρq_ice)
get_prog_state(state, ::Rain) = (state.precipitation, :ρq_rai)
get_prog_state(state, ::Snow) = (state.precipitation, :ρq_sno)
get_prog_state(state, ::Tracers{N}) where {N} = (state.tracers, :ρχ)

get_specific_state(state, ::Mass) = (state, :ρ)
get_specific_state(state, ::Momentum) = (state, :u)
get_specific_state(state, ::Energy) = (state.energy, :e)
get_specific_state(state, ::ρθ_liq_ice) = (state.energy, :θ_liq_ice)
get_specific_state(state, ::TotalMoisture) = (state.moisture, :q_tot)
get_specific_state(state, ::LiquidMoisture) = (state.moisture, :q_liq)
get_specific_state(state, ::IceMoisture) = (state.moisture, :q_ice)
get_specific_state(state, ::Rain) = (state.precipitation, :q_rai)
get_specific_state(state, ::Snow) = (state.precipitation, :q_sno)
get_specific_state(state, ::Tracers{N}) where {N} = (state.tracers, :χ)

prognostic_vars(m::AtmosLinearModel) = (Mass(), Momentum(), Energy())


text:
@testset "delta" begin
    A = rand(10, 5)
    B = rand(20, 5)
    d = KernelFunctions.Delta()
    @test pairwise(d, A; dims=1) == Matrix(I, size(A, 1), size(A, 1))
    @test pairwise(d, A, B; dims=1) == zeros(size(A, 1), size(B, 1))
    @test d(1, 2) == 0
    @test d(1, 1) == 1
    @test_throws DimensionMismatch d(rand(3), rand(4))
end


text:
using CSV
using Query
using DataFrames

casesf = reformat(CSV.read("time_series_covid19_confirmed_global.csv"), :Date, :Cases)
deathf = reformat(CSV.read("time_series_covid19_deaths_global.csv"), :Date, :Deaths)
recovf = reformat(CSV.read("time_series_covid19_recovered_global.csv"), :Date, :Recoveries)

chinaf = joindates(select_country.([casesf, deathf, recovf], "China")...)
hubeif = joindates(select_province.([casesf, deathf, recovf], "China", "Hubei")...)

parsedates!(hubeif)

# cdf = @from i in casesf begin
#     @join j in deathf on i.Date equals j.Date #and i.Province equals j.Province
#     @where i.Province == j.Province  && i.Country == j.Country
#     @select {i.Date,i.Cases,j.Deaths,j.Province, j.Country}
#     @collect DataFrame
# end

caseplot(hubeif)

df = @from i in joindates(casesf, deathf, recovf) begin
  @where i.Cases != 0
  @select i
  @collect DataFrame
end

CSV.write("covid.csv", df)

italyf = joindates(select_country.([casesf, deathf, recovf], "Italy")...)
francef = joindates(select_mainland.([casesf, deathf, recovf], "France")...)
ukf = joindates(select_mainland.([casesf, deathf, recovf], "United Kingdom")...)
euf = vcat(italyf, francef, ukf)

CSV.write("eucovid.csv", euf)


text:
export proxLLR!, normLLR

"""
    proxLLR!(x::Vector{T}, λ::Float64=1e-6; kargs...) where T

proximal map for LLR regularization using singular-value-thresholding

# Arguments
* `x::Vector{T}`                - Vector to apply proximal map to
* `λ::Float64`                  - regularization parameter
* `shape::Tuple{Int}=[]`        - dimensions of the image
* `blockSize::Tuple{Int}=[2;2]` - size of patches to perform singluar value thresholding on
* `randshift::Bool=true`        - randomly shifts the patches to ensure translation invariance
"""
function proxLLR!(x::Vector{T}, λ; shape::NTuple{N,TI}=error(),
   blockSize::NTuple{N,TI}=ntuple(_-> 2, N), randshift::Bool=true) where {T, N,TI <: Integer}

  x = reshape(x, tuple(shape..., length(x) ÷ prod(shape)))

  block_idx = CartesianIndices(blockSize)
  K = size(x)[end]

  if randshift
    # Random.seed!(1234)
    shift_idx = (Tuple(rand(block_idx))..., 0)
    xs = circshift(x, shift_idx)
  else
    xs = x
  end

  ext = mod.(shape,blockSize)
  pad = mod.(blockSize .- ext, blockSize)
  if any(pad .!= 0)
    xp = zeros(T, (shape .+ pad)..., K)
    xp[CartesianIndices(x)] .= xs
  else
    xp = xs
  end

  xᴸᴸᴿ = Array{T}(undef, prod(blockSize), K)
  for i ∈ CartesianIndices(StepRange.(0, blockSize, shape .- 1))
    @views xᴸᴸᴿ .= reshape(xp[i .+ block_idx,:], :, K)
    # threshold singular values
    SVDec = svd!(xᴸᴸᴿ)
    proxL1!(SVDec.S,λ)
    xp[i .+ block_idx,:] .= reshape(SVDec.U * Diagonal(SVDec.S) * SVDec.Vt, blockSize..., :)
  end

  if any(pad .!= 0)
    xs .= xp[CartesianIndices(xs)]
  end

  if randshift
    x .= circshift(xs, -1 .* shift_idx)
  end

  x = vec(x)
  return x
end

"""
    normLLR(x::Vector{T}, λ::Float64; kargs...) where T

returns the value of the LLR-regularization term.
Arguments are the same is in `proxLLR!`
"""
function normLLR(x::Vector{T}, λ::Float64; shape::NTuple{N,TI}, L=1, blockSize::NTuple{N,TI}=ntuple(_-> 2, N), randshift::Bool=true, kargs...) where {N, T, TI <: Integer}

  Nvoxel = prod(shape)
  K = floor(Int,length(x)/(Nvoxel*L))
  normᴸᴸᴿ = 0.
  for i = 1:L
    normᴸᴸᴿ +=  blockNuclearNorm(x[(i-1)*Nvoxel*K+1:i*Nvoxel*K], shape; blockSize=blockSize, randshift=randshift, kargs...)
  end

  return λ*normᴸᴸᴿ
end

function blockNuclearNorm(x::Vector{T}, shape::NTuple{N,TI}; blockSize::NTuple{N,TI}=ntuple(_-> 2, N),
      randshift::Bool=true, kargs...) where {N, T, TI <: Integer}
    x = reshape( x, tuple( shape...,floor(Int64, length(x)/prod(shape)) ) )

    Wy = blockSize[1]
    Wz = blockSize[2]

    if randshift
      srand(1234)
      shift_idx = [rand(1:Wy) rand(1:Wz) 0]
      x = circshift(x, shift_idx)
    end

    ny, nz, K = size(x)

    # reshape into patches
    L = floor(Int,ny*nz/Wy/Wz) # number of patches, assumes that image dimensions are divisble by the blocksizes

    xᴸᴸᴿ = zeros(T,Wy*Wz,L,K)
    for i=1:K
      xᴸᴸᴿ[:,:,i] = im2colDistinct(x[:,:,i], (Wy,Wz))
    end
    xᴸᴸᴿ = permutedims(xᴸᴸᴿ,[1 3 2])

    # L1-norm of singular values
    normᴸᴸᴿ = 0.
    for i = 1:L
      SVDec = svd(xᴸᴸᴿ[:,:,i])
      normᴸᴸᴿ += norm(SVDec.S,1)
    end

    return normᴸᴸᴿ
end


text:
# # Basic Usage

using Convex
using LinearAlgebra
if VERSION < v"1.2.0-DEV.0"
    (I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))
end

using SCS
## passing in verbose=0 to hide output from SCS
solver = SCSSolver(verbose=0)

# ### Linear program
#
# $$
# \begin{array}{ll}
#   \mbox{maximize} & c^T x \\
#   \mbox{subject to} & A x \leq b\\
#   & x \geq 1 \\
#   & x \leq 10 \\
#   & x_2 \leq 5 \\
#   & x_1 + x_4 - x_2 \leq 10 \\
# \end{array}
# $$
#

x = Variable(4)
c = [1; 2; 3; 4]
A = I(4)
b = [10; 10; 10; 10]
p = minimize(dot(c, x)) # or c' * x
p.constraints += A * x <= b
p.constraints += [x >= 1; x <= 10; x[2] <= 5; x[1] + x[4] - x[2] <= 10]
solve!(p, solver)

println(round(p.optval, digits=2))
println(round.(x.value, digits=2))
println(evaluate(x[1] + x[4] - x[2]))

# ### Matrix Variables and promotions
#
# $$
# \begin{array}{ll}
#   \mbox{minimize} & \| X \|_F + y \\
#   \mbox{subject to} & 2 X \leq 1\\
#   & X' + y \geq 1 \\
#   & X \geq 0 \\
#   & y \geq 0 \\
# \end{array}
# $$
#

X = Variable(2, 2)
y = Variable()
## X is a 2 x 2 variable, and y is scalar. X' + y promotes y to a 2 x 2 variable before adding them
p = minimize(norm(X) + y, 2 * X <= 1, X' + y >= 1, X >= 0, y >= 0)
solve!(p, solver)
println(round.(X.value, digits=2))
println(y.value)
p.optval

# ### Norm, exponential and geometric mean
#
# $$
# \begin{array}{ll}
#   \mbox{satisfy} & \| x \|_2 \leq 100 \\
#   & e^{x_1} \leq 5 \\
#   & x_2 \geq 7 \\
#   & \sqrt{x_3 x_4} \geq x_2
# \end{array}
# $$
#

x = Variable(4)
p = satisfy(norm(x) <= 100, exp(x[1]) <= 5, x[2] >= 7, geomean(x[3], x[4]) >= x[2])
solve!(p, SCSSolver(verbose=0))
println(p.status)
x.value

# ### SDP cone and Eigenvalues

y = Semidefinite(2)
p = maximize(lambdamin(y), tr(y)<=6)
solve!(p, SCSSolver(verbose=0))
p.optval

#-

x = Variable()
y = Variable((2, 2))
## SDP constraints
p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1)
solve!(p, solver)
y.value

# ### Mixed integer program
#
# $$
# \begin{array}{ll}
#   \mbox{minimize} & \sum_{i=1}^n x_i \\
#     \mbox{subject to} & x \in \mathbb{Z}^n \\
#   & x \geq 0.5 \\
# \end{array}
# $$
#

using GLPKMathProgInterface
x = Variable(4, :Int)
p = minimize(sum(x), x >= 0.5)
solve!(p, GLPKSolverMIP())
x.value

#-


text:
using NNlib: conv

@generated sub2(::Type{Val{N}}) where N = :(Val{$(N-2)})

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

"""
    Conv(size, in=>out)
    Conv(size, in=>out, relu)

Standard convolutional layer. `size` should be a tuple like `(2, 2)`.
`in` and `out` specify the number of input and output channels respectively.

Data should be stored in WHCN order. In other words, a 100×100 RGB image would
be a `100×100×3` array, and a batch of 50 would be a `100×100×3×50` array.

Takes the keyword arguments `pad`, `stride` and `dilation`.
"""
struct Conv{N,F,A,V}
  σ::F
  weight::A
  bias::V
  stride::NTuple{N,Int}
  pad::NTuple{N,Int}
  dilation::NTuple{N,Int}
end

Conv(w::AbstractArray{T,N}, b::AbstractVector{T}, σ = identity;
     stride = 1, pad = 0, dilation = 1) where {T,N} =
  Conv(σ, w, b, expand.(sub2(Val{N}), (stride, pad, dilation))...)

Conv(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity; init = initn,
     stride = 1, pad = 0, dilation = 1) where N =
  Conv(param(init(k..., ch...)), param(zeros(ch[2])), σ,
       stride = stride, pad = pad, dilation = dilation)

Flux.treelike(Conv)

function (c::Conv)(x)
  # TODO: breaks gpu broadcast :(
  # ndims(x) == ndims(c.weight)-1 && return squeezebatch(c(reshape(x, size(x)..., 1)))
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  σ.(conv(x, c.weight, stride = c.stride, pad = c.pad, dilation = c.dilation) .+ b)
end

function Base.show(io::IO, l::Conv)
  print(io, "Conv(", size(l.weight)[1:ndims(l.weight)-2])
  print(io, ", ", size(l.weight, ndims(l.weight)-1), "=>", size(l.weight, ndims(l.weight)))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end


text:
"""
The span of a dictionary is the set of all possible expansions in that
dictionary, with coefficient eltype determined uniquely by the Dictionary.

The span of a dictionary is a function space, mapping `S` to `T`. Here, `S` is
the domain type of the dictionary. The `T` is the codomain type.
"""
struct Span{S,T} <: FunctionSpace{S,T}
    dictionary  ::  Dictionary{S}
end

Span(dict::Dictionary{S}) where S = Span{S,span_codomaintype(dict)}(dict)

span_codomaintype(dict::Dictionary) =
    span_codomaintype(coefficienttype(dict), codomaintype(dict))
span_codomaintype(dict::Dictionary, coefficients) =
    span_codomaintype(eltype(coefficients), codomaintype(dict))
# - When the types are the same, that type is the result
span_codomaintype(::Type{T}, ::Type{T}) where {T <: Number} = T
# - the coefficient types are complex and the set itself is real
span_codomaintype(::Type{Complex{T}}, ::Type{T}) where {T <: Number} = Complex{T}
# - or vice-versa
span_codomaintype(::Type{T}, ::Type{Complex{T}}) where {T <: Number} = Complex{T}
# Default fallback
span_codomaintype(::Type{A}, ::Type{Z}) where {Z,A} = typeof(zero(A) * zero(Z))

coefficienttype(span::Span) = coefficienttype(dictionary(span))

# Convenient shorthand
coefficienttype = coefficienttype

dictionary(s::Span) = s.dictionary

similar(s::Span, ::Type{T}, dims) where {T} = Span(similar(dictionary(s), T, dims))

promote_domaintype(span::Span, ::Type{T}) where {T} = similar(span, T, size(span))

random_expansion(span::Span) = Expansion(dictionary(span), rand(dictionary(span)))

zero(span::Span) = Expansion(dictionary(span), zeros(dictionary(span)))

tensorproduct(s1::Span, s2::Span) = Span(tensorproduct(dictionary(s1), dictionary(s2)))

size(span::Span) = size(dictionary(span))
length(span::Span) = length(dictionary(span))

show(io::IO, mime::MIME"text/plain", s::Span) = composite_show(io, mime, s)
Display.displaystencil(s::Span) = ["Span(", dictionary(s), ")"]
show(io::IO, s::Span) = print(io, "Space : $(domaintype(s)) → $(codomaintype(s)) (span of $(repr(dictionary(s))))")