text:
using HDF5Arrays, HDF5, Test

@testset "HDF5Arrays.jl" begin

@testset "Constructors" begin
    mktemp() do path, io
        path *= ".h5"
        h5write(path, "group/x", rand(Int, 200, 50))
        f = h5open(path, "r+")
        g = f["group"]
        dset = g["x"]
        x = read(dset)

        @test HDF5Array(f, "group/x") == x
        @test HDF5Array(path, "group/x") == x
        @test HDF5Array(g, "x") == x
        @test HDF5Array(dset) == x

        @test_throws ErrorException get_chunk(HDF5Array(dset))

        g["c", "chunk", (5, 6, 7)] = rand(15, 8, 10)
        @test get_chunk(HDF5Array(g, "c")) == get_chunk(g["c"]) == (5, 6, 7)

        # Implementation detail
        @test HDF5Array(g, "c") isa HDF5Array{Float64, 3, HDF5Dataset, NTuple{3, Int}}

        y = HDF5Array{Float64, 2}(path, "group/y", undef, 100, 5)
        @test g["y"] isa HDF5Dataset
        @test eltype(g["y"]) == Float64
        @test size(g["y"]) == (100, 5)

        z = HDF5Array{Int32, 3}(g, "z", undef, 3, 5, 7, chunk=(3, 4, 5))
        @test g["z"] isa HDF5Dataset
        @test eltype(g["z"]) == Int32
        @test size(g["z"]) == (3, 5, 7)
        @test get_chunk(g["z"]) == (3, 4, 5)

        a = HDF5Array{Float32, 2}(g, "a", undef, 20, 7, chunk=(7, 7), max_dims=(-1, 14))
        @test g["a"] isa HDF5Dataset
        @test eltype(g["a"]) == Float32
        @test size(g["a"]) == (20, 7)
        @test get_chunk(g["a"]) == (7, 7)

        @test_throws ErrorException set_dims!(g["a"], (20, 15))

        set_dims!(g["a"], (1000, 14))
        @test size(a) == size(g["a"]) == (1000, 14)

        b = HDF5Array{Int, 4}(g, "b", undef, 6, 6, 6, 6, chunk=(3, 3, 3, 3), compress=5)
        @test g["b"] isa HDF5Dataset
        @test eltype(g["b"]) == Int
        @test size(g["b"]) == (6, 6, 6, 6)
        @test get_chunk(g["b"]) == (3, 3, 3, 3)
    end
end

@testset "Read / Write" begin
    mktemp() do path, io
        path *= ".h5"
        h5write(path, "group/x", rand(10, 15, 20))
        dset = h5open(path, "r+")["group/x"]
        x = read(dset)

        arr = HDF5Array(dset)

        @test all(arr .== x)

        @test arr[:, 2, 10:end] == x[:, 2, 10:end]

        arr[5, 5, 5] = 12345.6
        @test arr[5, 5, 5] == 12345.6
        @test arr[CartesianIndex(5, 5, 5)] == 12345.6
        @test dset[5, 5, 5] == 12345.6

        @inferred arr[5, 5, 5]

        arr[:, 1, 1] = 1:10
        @test arr[:, 1, 1] == 1:10
        @test all(arr[CartesianIndices((1:10, 1:1, 1:1))] .== 1:10)
        @test dset[:, 1, 1] == 1:10

        @test_broken @inferred arr[:, 1, 1]

        yarr = HDF5Array{Float64, 2}(path, "group/y", undef, 30, 5)
        ydset = h5open(path, "r+")["group/y"]
        ydset[1, 5] = 5432.1

        @test yarr[1, 5] == 5432.1

        yarr[:, 2] = 1:30
        @test yarr[:, 2] == 1:30
    end
end

@testset "Iteration Interface" begin
    mktemp() do path, io
        path *= ".h5"
        h5open(path, "w")["x", "chunk", (3, 2)] = Array(reshape(1:60, 15, 4))

        arr = HDF5Array(path, "x")

        x, s = iterate(arr)
        @test firstindex(arr) == CartesianIndex(1, 1)
        @test x == first(x) == 1

        x, s = iterate(arr, s) # 2
        x, s = iterate(arr, s) # 3
        x, s = iterate(arr, s) # next col => x = 16
        @test x == 16

        x, s = iterate(arr, s) # 17
        x, s = iterate(arr, s) # 18
        x, s = iterate(arr, s) # next chunk => x = 4
        @test x == 4

        c = collect(y for y in arr)
        @test length(c) == length(arr)
        @test lastindex(arr) == CartesianIndex(15, 4)
        @test last(c) == last(arr)
    end
end

# @testset "Chained" begin
#     path1, path2, path3 = tempname() * ".h5", tempname() * ".h5", tempname() * ".h5"
#     h5write(path1, "group/x", rand(50))
#     h5write(path2, "group/x", rand(100))
#     h5write(path3, "group/x", rand(75))
#
#     x = vcat(
#         h5read(path1, "group/x"),
#         h5read(path2, "group/x"),
#         h5read(path3, "group/x")
#     )
#
#     arr = HDF5Array((path1, path2, path3), "group/x")
#     @test arr == x
#
#     arr[26:75] = 1:50
#     @test h5read(path1, "group/x")[26:50] == 1:25
#     @test h5read(path2, "group/x")[1:25] == 26:50
# end

end


text:
using Test
using MPI

# TODO: enable CUDA tests once OpenMPI has full support
ArrayType = Array

MPI.Init()

comm = MPI.COMM_WORLD

const rank = MPI.Comm_rank(comm)
const N = MPI.Comm_size(comm)

buf = ArrayType{Int}(fill(Int(rank),N))
received = ArrayType{Int}(fill(-1,N))

# Test Get using fence
win = MPI.Win_create(buf, comm)
MPI.Win_fence(0, win)
MPI.Get!(received, (rank+1)%N, win)
MPI.Win_fence(0, win)

@test received == ArrayType{Int}(fill(Int(rank+1)%N, size(received)))

# Test locked window
if rank != 0
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 0, 0, win)
  received[1] = rank
  MPI.Put!(view(received,1:1), 0, rank, win)
  MPI.Win_unlock(0, win)
else
  buf[1] = 0
end

MPI.Win_fence(0, win)

if rank == 0
  @test buf == collect(0:N-1)
end

MPI.Barrier(comm)

if rank == 1
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 1, 0, win)
  fill!(buf,3)
  MPI.Win_unlock(1,win)
end

MPI.Barrier(comm)

if rank == 0
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 0, 0, win)
  fill!(buf, 2)
  MPI.Win_unlock(0,win)
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 1, 0, win)
  result = similar(buf)
  MPI.Get_accumulate!(buf, result, 1, 0, MPI.SUM, win)
  MPI.Win_unlock(1,win)
  @test all(result .== 3)
end

MPI.Barrier(comm)

if rank == 1
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 1, 0, win)
  @test all(buf .== 5)
  fill!(buf,-2)
  MPI.Win_unlock(1,win)
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 0, 0, win)
  MPI.Accumulate!(buf, 0, 0, MPI.SUM, win)
  MPI.Win_unlock(0,win)
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 1, 0, win)
  fill!(buf,1)
  MPI.Win_unlock(1,win)
end

MPI.Barrier(comm)

if rank == 0
  MPI.Win_lock(MPI.LOCK_EXCLUSIVE, 0, 0, win)
  @test all(buf .== 0)
  MPI.Win_unlock(0,win)
end

MPI.Barrier(comm)
MPI.free(win)
MPI.Barrier(comm)

# Try a dynamic window
win = MPI.Win_create_dynamic(comm)
MPI.Win_attach!(win, buf)

address_buf = Vector{Cptrdiff_t}(undef, 1)
address_win = MPI.Win_create(address_buf, comm)

# Address for dynamic window is global instead of relative to the start of the buffer, so we need to share it
MPI.Win_lock(MPI.LOCK_EXCLUSIVE, rank, 0, address_win)
address_buf[1] = MPI.Get_address(buf)
MPI.Win_unlock(rank, address_win)

MPI.Barrier(comm)

# Test atomic operation
if rank == 0
  received = Ref(0)
  to_send = Ref(0)
  for r in 0:N-1
    address = Ref{Cptrdiff_t}(0)
    MPI.Win_lock(MPI.LOCK_EXCLUSIVE, r, 0, address_win)
    MPI.Get!(address, r, address_win)
    MPI.Win_flush(r, address_win)
    to_send[] = r+5
    MPI.Win_lock(MPI.LOCK_EXCLUSIVE, r, 0, win)
    MPI.Fetch_and_op!(to_send, received, r, address[]+r*sizeof(Int), MPI.REPLACE, win)
    MPI.Win_flush(r, win)
    @test received[] == r
    MPI.Win_unlock(r, win)
    MPI.Win_unlock(r, address_win)
  end
end

MPI.Barrier(comm)

@test buf[rank+1] == rank+5

MPI.Barrier(comm)
MPI.Win_detach!(win, buf)
MPI.free(win)
MPI.free(address_win)

MPI.Finalize()
@test MPI.Finalized()

ext: jl


text:
using Base: @get!

@nograd readline

@adjoint copy(x::AbstractArray) = copy(x), ȳ -> (ȳ,)

grad_mut(d::AbstractDict) = Dict()

# TODO perhaps look up mutable gradients in `forward`
function accum(a::AbstractDict, b::AbstractDict)
  @assert a === b
  return a
end

@adjoint function getindex(d::AbstractDict, k)
  d[k], function (Δ)
    grad = grad_mut(__context__, d)
    grad[k] = accum(get(grad, k, nothing), Δ)
    return (grad, nothing)
  end
end

@adjoint! function setindex!(d::AbstractDict, v, k)
  setindex!(d, v, k), function (Δ)
    (nothing, get(grad_mut(__context__, d), k, nothing), nothing)
  end
end


text:
### CONCRETE TYPE: DIRECT PROX EVALUATION
# prox! is computed using a Cholesky factorization of A'A + I/(lambda*gamma)
# or AA' + I/(lambda*gamma), according to which matrix is smaller.
# The factorization is cached and recomputed whenever gamma changes

using LinearAlgebra
using SparseArrays
using SuiteSparse

mutable struct LeastSquaresDirect{R <: Real, RC <: RealOrComplex{R}, M <: AbstractMatrix{RC}, V <: AbstractVector{RC}, F <: Factorization} <: LeastSquares
  A::M # m-by-n matrix
  b::V
  Atb::V
  lambda::R
  gamma::R
  shape::Symbol
  S::M
  res::Vector{RC} # m-sized buffer
  q::Vector{RC} # n-sized buffer
  fact::F
  function LeastSquaresDirect{R, RC, M, V, F}(A::M, b::V, lambda::R) where {R <: Real, RC <: RealOrComplex{R}, M <: AbstractMatrix{RC}, V <: AbstractVector{RC}, F <: Factorization}
    if size(A, 1) != length(b)
      error("A and b have incompatible dimensions")
    end
    if lambda <= 0
      error("lambda must be positive")
    end
    m, n = size(A)
    if m >= n
      S = A'*A
      shape = :Tall
    else
      S = A*A'
      shape = :Fat
    end
    new(A, b, A'*b, lambda, -1, shape, S, zeros(RC, m), zeros(RC, n))
  end
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, RC <: Union{R, Complex{R}}, M <: DenseMatrix{RC}, V <: AbstractVector{RC}}
  LeastSquaresDirect{R, RC, M, V, Cholesky{RC, M}}(A, b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, RC <: Union{R, Complex{R}}, I <: Integer, M <: SparseMatrixCSC{RC, I}, V <: AbstractVector{RC}}
  LeastSquaresDirect{R, RC, M, V, SuiteSparse.CHOLMOD.Factor{RC}}(A, b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, RC <: Union{R, Complex{R}}, M <: AbstractMatrix{RC}, V <: AbstractVector{RC}}
  warn("Could not infer type of Factorization for $M in LeastSquaresDirect, this type will be type-unstable")
  LeastSquaresDirect{R, RC, M, V, Factorization}(A, b, lambda)
end

function (f::LeastSquaresDirect{R, RC, M, V, F})(x::AbstractVector{RC}) where {R, RC, M, V, F}
  mul!(f.res, f.A, x)
  f.res .-= f.b
  return (f.lambda/2)*norm(f.res, 2)^2
end

function prox!(y::AbstractVector{D}, f::LeastSquaresDirect{R, RC, M, V, F}, x::AbstractVector{D}, gamma::R=one(R)) where {R, RC, M, V, F, D <: RealOrComplex{R}}
  # if gamma different from f.gamma then call factor_step!
  if gamma != f.gamma
    factor_step!(f, gamma)
  end
  solve_step!(y, f, x, gamma)
  mul!(f.res, f.A, y)
  f.res .-= f.b
  return (f.lambda/2)*norm(f.res, 2)^2
end

function factor_step!(f::LeastSquaresDirect{R, RC, M, V, F}, gamma::R) where {R, RC, M <: DenseMatrix, V, F}
  lamgam = f.lambda*gamma
  f.fact = cholesky(f.S + I/lamgam)
  f.gamma = gamma
end

function factor_step!(f::LeastSquaresDirect{R, RC, M, V, F}, gamma::R) where {R, RC, M <: SparseMatrixCSC, V, F}
  lamgam = f.lambda*gamma
  f.fact = ldlt(f.S; shift = 1.0/lamgam)
  f.gamma = gamma
end

function solve_step!(y::AbstractVector{D}, f::LeastSquaresDirect{R, RC, M, V, F}, x::AbstractVector{D}, gamma::R) where {R, RC, M, V, F <: Cholesky{RC, M}, D <: RealOrComplex{R}}
  lamgam = f.lambda*gamma
  f.q .= f.Atb .+ x./lamgam
  # two cases: (1) tall A, (2) fat A
  if f.shape == :Tall
    # y .= f.fact\f.q
    y .= f.q
    LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, y)
    LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, y)
  else # f.shape == :Fat
    # y .= lamgam*(f.q - (f.A'*(f.fact\(f.A*f.q))))
    mul!(f.res, f.A, f.q)
    LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, f.res)
    LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, f.res)
    mul!(y, adjoint(f.A), f.res)
    y .-= f.q
    y .*= -lamgam
  end
end

function solve_step!(y::AbstractVector{D}, f::LeastSquaresDirect{R, RC, M, V, F}, x::AbstractVector{D}, gamma::R) where {R, RC, M, V, F <: SuiteSparse.CHOLMOD.Factor{RC}, D <: RealOrComplex{R}}
  lamgam = f.lambda*gamma
  f.q .= f.Atb .+ x./lamgam
  # two cases: (1) tall A, (2) fat A
  if f.shape == :Tall
    y .= f.fact\f.q
  else # f.shape == :Fat
    # y .= lamgam*(f.q - (f.A'*(f.fact\(f.A*f.q))))
    mul!(f.res, f.A, f.q)
    f.res .= f.fact\f.res
    mul!(y, adjoint(f.A), f.res)
    y .-= f.q
    y .*= -lamgam
  end
end

function gradient!(y::AbstractVector{D}, f::LeastSquaresDirect{R, RC, M, V, F}, x::AbstractVector{D}) where {R, RC, M, V, F, D <: Union{R, Complex{R}}}
  mul!(f.res, f.A, x)
  f.res .-= f.b
  mul!(y, adjoint(f.A), f.res)
  y .*= f.lambda
  fy = (f.lambda/2)*dot(f.res, f.res)
end

function prox_naive(f::LeastSquaresDirect, x::AbstractVector, gamma=1.0)
  lamgam = f.lambda*gamma
  y = (f.A'*f.A + I/lamgam)\(f.Atb + x/lamgam)
  fy = (f.lambda/2)*norm(f.A*y-f.b)^2
  return y, fy
end


text:
# taken from Peaks.jl (https://github.com/halleysfifthinc/Peaks.jl)

export maxima,
       minima

for (funcname, comp, notcomp) in ((:maxima, :<, :>),
                         (:minima, :>, :<))

    @eval begin
        function ($funcname)(x::AbstractVector,
                        w::Integer=1,
                        strictbounds::Bool=true)
            w > 0 || throw(ArgumentError("window cannot be negative"))
            xlen = length(x)
            idxs = Int[]

            # There can't be more than one peak every `w` elements
            maxN = strictbounds ? cld(xlen,2*w) : fld(xlen,2*w)
            N = 0

            ww = [ collect(-w:-1); collect(1:w) ]

            sizehint!(idxs,maxN)

            firsti = firstindex(x)
            lasti = lastindex(x)
            i = firsti
            maxI = lasti - w

            if !strictbounds
                while i <= firsti + w
                    xi = x[i]
                    peak = true
                    if ismissing(xi)
                        i = something(findnext(!ismissing, x, i+1), lasti+1)
                    elseif isnan(xi)
                        i = something(findnext(!isnan, x, i+1), lasti+1)
                    else
                        for j in max(i-w, firsti):(i+w)
                            i === j && continue
                            xj = x[j]
                            if coalesce(($comp)(xi, xj), strictbounds) || (strictbounds && (ismissing(xj) || isnan(xj)))
                                peak &= false
                                break # No reason to continue checking if the rest of the elements qualify
                            elseif xi === xj
                                k = findnext(y -> xi !== y, x, j+1)
                                if isnothing(k) # x is constant till the end, not a peak
                                    if !strictbounds
                                        push!(idxs,i)
                                        N += 1
                                        i = lasti+1
                                        peak &= false
                                        break
                                    else
                                        i = lasti+1
                                        peak &= false
                                        break
                                    end
                                else
                                    xk = x[k]
                                    if coalesce(($comp)(xi, xk), strictbounds) || (strictbounds && (ismissing(xk) || isnan(xk))) # x moves towards a peak
                                        peak &= false
                                        break
                                    else # Push new peak here to shift the right number of elements
                                        push!(idxs,i)
                                        N += 1
                                        i = max(k,i+w)
                                        peak &= false
                                        break
                                    end
                                end
                            end
                        end

                        if peak
                            push!(idxs,i)
                            N += 1
                            i += w # There can't be another peak for at least `w` more elements
                        end
                        i += 1
                    end
                end
            end

            i = max(i, firsti + w, something(findfirst(y -> y !== x[firsti], x), firsti))

            while i <= maxI
                xi = x[i]
                peak = true
                if ismissing(xi)
                    i = something(findnext(!ismissing, x, i+1), lasti+1)
                elseif isnan(xi)
                    i = something(findnext(!isnan, x, i+1), lasti+1)
                else
                    for j in ww # For all elements within the window
                        xj = x[i+j]
                        if coalesce(($comp)(xi, xj), strictbounds) || (strictbounds && (ismissing(xj) || isnan(xj)))
                            peak &= false
                            break # No reason to continue checking if the rest of the elements qualify
                        elseif xi === xj
                            xi1 = x[i-1]
                            if coalesce(($notcomp)(xi, xi1), !strictbounds) || (!strictbounds && isnan(xi1))
                                k = findnext(y -> xi !== y, x, i+j+1)
                                if isnothing(k)
                                    if !strictbounds # x is constant till the end, not a peak
                                        push!(idxs,i)
                                        N += 1
                                        i = lasti+1
                                        peak &= false
                                        break
                                    else
                                        i = lasti+1
                                        peak &= false
                                        break
                                    end
                                else
                                    xk = x[k]
                                    if coalesce(($comp)(xi, xk), strictbounds) || (strictbounds && (ismissing(xk) || isnan(xk)))
                                        # x moves towards a peak or ends with a missing or NaN
                                        peak &= false
                                        break
                                    elseif j === 1 # Push first element of plateau as peak here to shift the correct number of elements
                                        push!(idxs,i)
                                        N += 1
                                        i = max(k,i+w)
                                        peak &= false
                                        break
                                    end
                                end
                            else
                                peak &= false
                                break
                            end
                        end
                    end

                    if peak
                        push!(idxs,i)
                        N += 1
                        i += w # There can't be another peak for at least `w` more elements
                    end
                    i += 1
                end
            end

            if !strictbounds
                while i <= lasti
                    xi = x[i]
                    peak = true
                    if ismissing(xi)
                        i = something(findnext(!ismissing, x, i+1), lasti+1)
                    elseif isnan(xi)
                        i = something(findnext(!isnan, x, i+1), lasti+1)
                    else
                        for j in (i-w):min((i+w),lasti)
                            i === j && continue
                            xj = x[j]
                            if coalesce(($comp)(xi, xj), strictbounds) || (strictbounds && (ismissing(xj) || isnan(xj)))
                                peak &= false
                                break # No reason to continue checking if the rest of the elements qualify
                            elseif xi === xj
                                xi1 = x[i-1]
                                if coalesce(($notcomp)(xi, xi1), !strictbounds) || (!strictbounds && isnan(xi1))
                                    k = findnext(y -> xi !== y, x, j+1)
                                    if isnothing(k)
                                        if !strictbounds # x is constant till the end, not a peak
                                            push!(idxs,i)
                                            N += 1
                                            i = lasti+1
                                            peak &= false
                                            break
                                        else
                                            i = lasti+1
                                            peak &= false
                                            break
                                        end
                                    else
                                        xk = x[k]
                                        if coalesce(($comp)(xi, xk), strictbounds) || (strictbounds && (ismissing(xk) || isnan(xk))) # x moves towards a peak
                                            peak &= false
                                            break
                                        else # Push new peak here to shift the right number of elements
                                            push!(idxs,i)
                                            N += 1
                                            i = max(k,i+w)
                                            peak &= false
                                            break
                                        end
                                    end
                                else
                                    peak &= false
                                    break
                                end
                            end
                        end

                        if peak
                            push!(idxs,i)
                            N += 1
                            i += w # There can't be another peak for at least `w` more elements
                        end
                        i += 1
                    end
                end
            end

            resize!(idxs,N)
            sizehint!(idxs,N)

            return idxs
        end
    end
end

@doc """
    maxima(x[, w=1, strictbounds=true])

Find the indices of the local maxima of `x` where each maxima is either the maximum of
`x[-w:w]` or the first index of a plateau.

If `strictbounds` is `true`, all elements of `x[-w:w]` must exist and may not be `missing`
or `NaN`. If `strictbounds` is `false`, elements of `x[-w:w]` may not exist (eg peaks may
be less than `w` indices from either end of `x`), or may be `missing` or `NaN`. `missing`
or `NaN` must not be maxima.
"""
maxima

@doc """
    minima(x[, w=1, strictbounds=false])

Find the indices of the local minima of `x` where each minima is either the minimum of
`x[-w:w]` or the first index of a plateau.

If `strictbounds` is `true`, all elements of `x[-w:w]` must exist and may not be `missing`
or `NaN`. If `strictbounds` is `false`, elements of `x[-w:w]` may not exist (eg peaks may
be less than `w` indices from either end of `x`), or may be `missing` or `NaN`. `missing`
or `NaN` must not be minima.
"""
minima


text:
function DropTip(args...; phy::Phylo=tree,
    tip,
    trim_internal = true, subtree = false,
    root_edge::Int = 0,
    rooted = true, interactive = false)

    if isRooted(tree) == false
        error("isRooted:wrongtpe\n", "# Rooted phylogenetic tree required!")
    end
    Ntip = length(phy["tip_label"])
    if interactive == true
        xy = locator()
        nToDrop = length(xy["x"])
        tip = convert(Int, nToDrop)
        lastPP = get("last_plot.phylo", envir = PlotPhyloEnv)
        for i = 1:nToDrop
            d = sqrt(xy["x"][i]-lastPP["xx"]^2+(xy["y"][i]-lastPP["yy"])^2)
            tip[i] = indmin(d)
        end
    else
        if typeof(tip)==String
            tip = findin(phy["tip_label"],tip)
        end
    end

    out_of_range = []
    for i = 1:length(tip)
        push!(out_of_range, tip[i] > Ntip)
    end

    if any(out_of_range)
        warn("some tip numbers were larger than the number of tips: they were ignored")
        tip = tip[out_of_range .!=true]
    end
    if length(tip)==0
        return phy
    end
    if length(tip) == Ntip
        warn("drop all tips of the tree: returning NULL")
        return NULL
    end

    wbl = length(phy["edge_length"]) != 0
    if isRooted(tree) == false and subtree == false
        # phy = root(phy, (1:Ntip)[-tip][1])  # need to write the root function
        root_edge = 0
    end
    # phy = reorder(phy)
    NEWROOT = ROOT = Ntip + 1
    Nnode = phy["Nnode"]
    Nedge = size(phy["edge"])[1]
    if subtree == true
        trim_internal = true
        tr = reorder(phy, "postorder")
        N = ccall((:node_depth, "libc"), Array(Float64, 1),
            (Int, Int, Int, Int, Int, Int, Int),
            Ntip, Nnode, phy["edge"][:,1],phy["edge"][:,2],Nedge,
            convert(Array{Int64,1},zeros(Ntip+Nnode)),1)[6]
    end
    edge1 = phy["edge"][:,1]
    edge2 = phy["edge"][:,2]
    keep = trues(Nedge)
    mathidx = []
    for i = 1:length(edge2)
        if tip[i] == edge2[i]
            push!(mathidx, i)
        end
    end
    keep[mathidx] = false
    if trim_internal == true
        ints = []
        for i = 1:length(edge2)
            push!(ints, edge2[i] > Ntip)
        end
# need to add other options for subtree etc.
    end
end


text:
const GAS_CONSTANT_SI = 8.31446261815324 # Exact since May20, 2019; https://en.wikipedia.org/wiki/Gas_constant
const PA_IN_BAR = 1e5               # Сколько паскалей в баре
const BOLTZMANN_SI = 1.380649e-23 # Exact since May 20, 2019; https://en.wikipedia.org/wiki/Boltzmann_constant

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
export CUFFTError

struct CUFFTError <: Exception
    code::cufftResult
    msg::AbstractString
end
Base.show(io::IO, err::CUFFTError) = print(io, "CUFFTError(code $(err.code), $(err.msg))")

function CUFFTError(code::cufftResult)
    msg = status_message(code)
    return CUFFTError(code, msg)
end

function status_message(status)
    if status == CUFFT_STATUS_SUCCESS
        return "the operation completed successfully"
    elseif status == CUFFT_STATUS_INVALID_PLAN
        return "cuFFT was passed an invalid plan handle"
    elseif status == CUFFT_STATUS_ALLOC_FAILED
        return "cuFFT failed to allocate GPU or CPU memory"
    elseif status == CUFFT_STATUS_INVALID_TYPE
        return "cuFFT invalid type " # No longer used
    elseif status == CUFFT_STATUS_INVALID_VALUE
        return "User specified an invalid pointer or parameter"
    elseif status == CUFFT_STATUS_INTERNAL_ERROR
        return "Driver or internal cuFFT library error"
    elseif status == CUFFT_STATUS_EXEC_FAILED
        return "Failed to execute an FFT on the GPU"
    elseif status == CUFFT_STATUS_SETUP_FAILED
        return "The cuFFT library failed to initialize"
    elseif status == CUFFT_STATUS_INVALID_SIZE
        return "User specified an invalid transform size"
    elseif status == CUFFT_STATUS_UNALIGNED_DATA
        return "cuFFT unaligned data" # No longer used
    elseif status == CUFFT_STATUS_INCOMPLETE_PARAMETER_LIST
        return "Missing parameters in call"
    elseif status == CUFFT_STATUS_INVALID_DEVICE
        return "Execution of a plan was on different GPU than plan creation"
    elseif status == CUFFT_STATUS_PARSE_ERROR
        return "Internal plan database error"
    elseif status == CUFFT_STATUS_NO_WORKSPACE
        return "No workspace has been provided prior to plan execution"
    elseif status == CUFFT_STATUS_NOT_IMPLEMENTED
        return "Function does not implement functionality for parameters given."
    elseif status == CUFFT_STATUS_LICENSE_ERROR
        return "cuFFT license error" # Used in previous versions.
    elseif status == CUFFT_STATUS_NOT_SUPPORTED
        return "Operation is not supported for parameters given."
    else
        return "unknown status"
    end
end

macro check(fft_func)
    quote
        local err::cufftResult
        err = $(esc(fft_func::Expr))
        if err != CUFFT_STATUS_SUCCESS
            throw(CUFFTError(err))
        end
        err
    end
end


text:
export Model, loglik, elbo, predict

"""
    struct Model

# Fields
- `encoder`: Encoder.
- `decoder`: Decoder.
"""
struct Model
    encoder
    decoder
end

@Flux.functor Model

function (model::Model)(xc::AA, yc::AA, xt::AA; num_samples::Integer=1, kws...)
    size(xc, 1) == 0 && (xc = yc = nothing)  # Handle empty set case.
    xz, pz = code(model.encoder, xc, yc, xt; kws...)
    z = sample(pz, num_samples=num_samples)
    _, d = code(model.decoder, xz, z, xt)
    return d
end

"""
    loglik(
        model::Model,
        epoch::Integer,
        xc::AA,
        yc::AA,
        xt::AA,
        yt::AA;
        num_samples::Integer,
        batch_size::Integer=1024,
        importance_weighted::Bool=false,
        fixed_σ::Float32=1f-2,
        fixed_σ_epochs::Integer=0,
        kws...
    )

Log-expected-likelihood loss. This is a biased estimate of the log-likelihood.

# Arguments
- `model::Model`: Model.
- `epoch::Integer`: Current epoch.
- `xc::AA`: Locations of context set of shape `(n, dims, batch)`.
- `yc::AA`: Observed values of context set of shape `(n, channels, batch)`.
- `xt::AA`: Locations of target set of shape `(m, dims, batch)`.
- `yt::AA`: Observed values of target set of shape `(m, channels, batch)`.

# Keywords
- `num_samples::Integer`: Number of samples.
- `batch_size::Integer=1024`: Batch size to use in sampling.
- `importance_weighted::Bool=false`: Do an importance-weighted estimate.
- `fixed_σ::Float32=1f-2`: Hold the observation noise fixed to this value initially.
- `fixed_σ_epochs::Integer=0`: Number of iterations to hold the observation noise fixed for.
- `kws...`: Further keywords to pass on.

# Returns
- `Tuple{Real, Integer}`: Average negative log-expected likelihood and the "size" of the
    loss.
"""
function loglik(
    model::Model,
    epoch::Integer,
    xc::AA,
    yc::AA,
    xt::AA,
    yt::AA;
    num_samples::Integer,
    batch_size::Integer=1024,
    importance_weighted::Bool=false,
    fixed_σ::Float32=1f-2,
    fixed_σ_epochs::Integer=0,
    kws...
)
    n_target = size(xt, 1)

    # Determine batches.
    num_batches, batch_size_last = divrem(num_samples, batch_size)
    batches = Int[batch_size for _ = 1:num_batches]
    batch_size_last > 0 && push!(batches, batch_size_last)

    # Initialise variable that accumulates the log-pdfs.
    logpdfs = nothing

    # Concatenate inputs for IW estimate.
    if importance_weighted
        x_all = cat(xc, xt, dims=1)
        y_all = cat(yc, yt, dims=1)
    end

    # Handle empty set case.
    size(xc, 1) == 0 && (xc = yc = nothing)

    # Perform encoding.
    xz, pz, h = code_track(model.encoder, xc, yc, xt; kws...)

    # Construct posterior over latent variable for IW estimate.
    if importance_weighted
        qz = recode_stochastic(model.encoder, pz, x_all, y_all, h; kws...)
    end

    # Compute the loss in a batched way.
    for batch in batches
        if importance_weighted
            # Sample from posterior.
            z = sample(qz, num_samples=batch)

            # Do an importance weighted estimate.
            weights = sum(logpdf(pz, z), dims=(1, 2)) .- sum(logpdf(qz, z), dims=(1, 2))
        else
            # Sample from the prior.
            z = sample(pz, num_samples=batch)

            # Do a regular Monte Carlo estimate.
            weights = 0
        end

        # Perform decoding.
        _, d = code(model.decoder, xz, z, xt)

        # Fix the noise for the early epochs to force the model to fit.
        if epoch <= fixed_σ_epochs
            d = Gaussian(mean(d), [fixed_σ] |> gpu)
        end

        # Add a diagonal to correlated covariances to help the Cholesky decompositions.
        d = _regularise_multivariate_gaussian(d, epoch)

        # Perform Monte Carlo estimate.
        batch_logpdfs = weights .+ sum(logpdf(d, yt), dims=(1, 2))

        # Accumulate sum.
        logpdfs = isnothing(logpdfs) ? batch_logpdfs : cat(logpdfs, batch_logpdfs, dims=4)
        logpdfs = logsumexp(logpdfs, dims=4)
    end

    # Turn log-sum-exp into a log-mean-exp.
    logpdfs = logpdfs .- Float32(log(num_samples))

    # Return average over batches and the "size" of the loss.
    return -mean(logpdfs), size(xt, 1)
end

_regularise_multivariate_gaussian(d, epoch) = d
function _regularise_multivariate_gaussian(d::MultivariateGaussian, epoch; σ²=nothing)
    if isnothing(σ²)
        if epoch == 1
            σ² = 1f-1 * maximum(Tracker.data(std(d)))^2
        end
    end
    if !isnothing(σ²)
        n = size(mean(d), 1)
        ridge = Matrix(σ² * I, n, n) |> gpu
        return MultivariateGaussian(mean(d), var(d) .+ ridge)
    else
        return d
    end
end

"""
    elbo(
        model::Model,
        epoch::Integer,
        xc::AA,
        yc::AA,
        xt::AA,
        yt::AA;
        num_samples::Integer,
        fixed_σ::Float32=1f-2,
        fixed_σ_epochs::Integer=0,
        kws...
    )

Neural process ELBO-style loss. Subsumes the context set into the target set.

# Arguments
- `model::Model`: Model.
- `epoch::Integer`: Current epoch.
- `xc::AA`: Locations of context set of shape `(n, dims, batch)`.
- `yc::AA`: Observed values of context set of shape `(n, channels, batch)`.
- `xt::AA`: Locations of target set of shape `(m, dims, batch)`.
- `yt::AA`: Observed values of target set of shape `(m, channels, batch)`.

# Keywords
- `num_samples::Integer`: Number of samples.
- `fixed_σ::Float32=1f-2`: Hold the observation noise fixed to this value initially.
- `fixed_σ_epochs::Integer=0`: Number of iterations to hold the observation noise fixed for.
- `kws...`: Further keywords to pass on.

# Returns
- `Tuple{Real, Integer}`: Average negative NP loss and the "size" of the loss.
"""
function elbo(
    model::Model,
    epoch::Integer,
    xc::AA,
    yc::AA,
    xt::AA,
    yt::AA;
    num_samples::Integer,
    fixed_σ::Float32=1f-2,
    fixed_σ_epochs::Integer=0,
    kws...
)
    # We subsume the context set into the target set for this ELBO.
    x_all = cat(xc, xt, dims=1)
    y_all = cat(yc, yt, dims=1)

    # Handle empty set case.
    size(xc, 1) == 0 && (xc = yc = nothing)

    # Perform deterministic and latent encoding.
    xz, pz, h = code_track(model.encoder, xc, yc, x_all; kws...)

    # Construct posterior over latent variable.
    qz = recode_stochastic(model.encoder, pz, x_all, y_all, h; kws...)

    # Sample latent variable and perform decoding.
    z = sample(qz, num_samples=num_samples)
    _, d = code(model.decoder, xz, z, x_all)

    # Fix the noise for the early epochs to force the model to fit.
    if epoch <= fixed_σ_epochs
        d = Gaussian(mean(d), [fixed_σ] |> gpu)
    end

    # Estimate ELBO from samples.
    elbos = mean(sum(logpdf(d, y_all), dims=(1, 2)), dims=4) .- sum(kl(qz, pz), dims=(1, 2))

    # Return average over batches and the "size" of the loss.
    return -mean(elbos), size(x_all, 1)
end

"""
    predict(model::Model, xc::AV, yc::AV, xt::AV; num_samples::Integer=10)

# Arguments
- `model::Model`: Model.
- `xc::AV`: Locations of observed values of shape `(n)`.
- `yc::AV`: Observed values of shape `(n)`.
- `xt::AV`: Locations of target values of shape `(m)`.

# Keywords
- `num_samples::Integer=10`: Number of posterior samples.
- `epoch::Integer=100`: Current epoch.
- `kw...`: Further keywords to pass on.

# Returns
- `Tuple`:  Tuple containing means, lower and upper 95% central credible bounds, and
    `num_samples` posterior samples.
"""
function predict(
    model::Model,
    xc::AV,
    yc::AV,
    xt::AV;
    num_samples::Integer=10,
    epoch::Integer=100,
    kws...
)
    # Run model.
    d = model(
        expand_gpu.((xc, yc, xt))...;
        # Use at least 20 samples to estimate uncertainty.
        num_samples=max(num_samples, 20),
        kws...
    )
    μ = mean(d)[:, 1, 1, :] |> cpu
    σ = std(d)[:, 1, 1, :] |> cpu

    if size(μ, 2) >= num_samples
        samples = μ[:, 1:num_samples]
    elseif d isa MultivariateGaussian
        # Regularise, because the covariance can be unstable.
        d = _regularise_multivariate_gaussian(d, epoch; σ²=1f-4)
        samples = sample(d, num_samples=num_samples)[:, 1, 1, :] |> cpu
    else
        # There are no samples.
        samples = nothing
    end

    # Estimate functional uncertainty. Do not use the correction because there may only
    # be one sample. We Gaussianise to make it appear smoother.
    μ_σ = std(μ, dims=2, corrected=false)
    ε = 2 .* mean(sqrt.(μ_σ.^2 .+ σ.^2), dims=2) # Add variances to make total error.

    # Compute bounds.
    μ = mean(μ, dims=2)
    lowers = μ .- ε
    uppers = μ .+ ε

    # Remove singleton dimensions.
    μ = dropdims(μ, dims=2)
    lowers = dropdims(lowers, dims=2)
    uppers = dropdims(uppers, dims=2)

    return μ, lowers, uppers, samples
end


text:
using AbstractMover
using Test

@testset "AbstractMover.jl" begin
    # Write your tests here.
end


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 89efa14e-b00c-11ec-09fe-15316ef94166
1+1

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═89efa14e-b00c-11ec-09fe-15316ef94166
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: seoyeonc/julia

max_issues_repo_name: seoyeonc/julia

max_forks_repo_name: seoyeonc/julia


text:
using DataStructures

const kernel_cache = Dict{Symbol,Nullable{Function}}()
const translate_queue = Tuple[]

"""
Function to be implemented by the various targets:
translate_kernel(target, kernel_id, signature)
"""
function translate_kernel
end

"""
Dispatch point for kernels.
"""
function kernel_call(kernel_id::Symbol, target, f::Function, args)
    # Keep me as fast as possible! Currently the overhead
    # is about 2-3x that of regular dispatch.
    nf = get(kernel_cache, kernel_id, Nullable{Function}())
    if isnull(nf)
        # Got a new kernel!
        signature = (typeof(f), map(typeof, args))
        nf = translate_cached(target, kernel_id, signature)
    end
    nf.value(f,args)
end

function translate_cached(target, kernel_id, signature)
    info("Compiling new kernel $kernel_id $signature for target $target")
    kernel_func = translate_kernel(target, kernel_id, signature)

    # Cache it:
    nf = Nullable{Function}(kernel_func)
    kernel_cache[kernel_id] = nf
    nf
end

"""
Main entry point (kernel-per-kernel translation).
"""
@generated function translated(f::Function, target, args...)
    # We use the @generated mechanism to associate a kernel ID
    # to each new call signature. The mechanism is robust to
    # methods being redefined, ie. it will correctly generate
    # a new kernel ID in that case. EDIT: NOT! need to wait
    # for fix to #265.
    # This is much faster (2-3x) than looking up the method in
    # Julia's dispatch table and then looking up the kernel ID
    # for that method at runtime.
    kernel_id = gensym("kernel")
    quote
        ASTIR.kernel_call($(QuoteNode(kernel_id)), target, f, args)
    end
end

"""
Main entry point with batch translation.

Batch translation is very advantageous when the translation involves a command-
line tool with high startup overhead, eg. compilers.

This entry point (ab)uses type inference to get a peek at kernels that *may*
be about to be executed, *before* they're actually executed. It builds a list
of candidates that can be batch-translated the first time the code is run.

Note: this technique may be in fact be proscribed, see:
https://github.com/JuliaLang/julia/issues/18568
"""
@generated function batch_translated{target}(f::Function, t::Type{target}, args...)
    kernel_id = gensym("kernel")
    signature = (f, args)
    push!(translate_queue, (target, kernel_id, signature))
    # This function is typically called only once per new type
    # signature, but may occasionally be called a handful of times,
    # which may lead to some redundant translations being triggered
    # in batch mode. The overhead is usually small enough to be ignored.
    quote
        ASTIR.translate_all()
        ASTIR.kernel_call($(QuoteNode(kernel_id)), t, f, args)
    end
end

function translate_all()
    if isempty(translate_queue)
        return # fast trivial path
    end

    # Group queued kernels by target:
    by_target = DefaultDict(DataType,Vector,[])
    while !isempty(translate_queue)
        target, kernel_id, signature = pop!(translate_queue)
        push!(by_target[target], (kernel_id, signature))
    end

    # Translate each target group:
    for (target, kernels) in by_target
        funcs = translate_batch(target, kernels)
        for (kernel_func, (kernel_id, signature)) in zip(funcs, kernels)
            # Cache it:
            nf = Nullable{Function}(kernel_func)
            kernel_cache[kernel_id] = nf
        end
    end
end

"""
Default batch translate implementation.
Targets can override this for better efficiency when translating many kernels at once.
"""
function translate_batch(target::Any, kernels)
    [translate_kernel(target, args...) for args in kernels]
end


text:
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# StreamOutlet.jl: type and method definitions for steam outlet

export StreamOutlet

export push_sample, push_chunk, have_consumers, wait_for_consumers

mutable struct StreamOutlet{T}
  handle::lib.lsl_outlet
  info::StreamInfo{T}

  function StreamOutlet(handle, info::StreamInfo{T}) where T
    outlet = new{T}(handle, info)
    finalizer(_destroy, outlet)
  end

end

# Destroy a StreamInfo handle
function _destroy(outlet::StreamOutlet)
  if outlet.handle != C_NULL
    lib.lsl_destroy_outlet(outlet)
    outlet.handle = C_NULL
  end
  return nothing
end

# Define conversion to pointer
unsafe_convert(::Type{lib.lsl_outlet}, outlet::StreamOutlet) = outlet.handle

# Define constructor
"""
    StreamOutlet(info; chunk_size = 0, max_buffered = 360)

Establish a new stream outlet. This makes the stream discoverable.

# arguments
- `info::StreamInfo`: The stream information to use for creating this stream. Stays constant
                      over the lifetime of the outlet.

# Keyword arguments
- `chunk_size::Integer`: Desired chunk granularity (in samples) for transmission. If
                         specified as 0, each push operation yields one chunk. Stream
                         recipients can have this setting bypassed.
- `max_buffered::Integer`: Maximum amount of data to buffer (in seconds if there is a
                           nominal sampling rate, otherwise x100 in samples). A good default
                           is 360, which corresponds to 6 minutes of data. Note that, for
                           high-bandwidth data you will almost certainly want to use a lower
                           value here to avoid running out of RAM.
"""
function StreamOutlet(info::StreamInfo{T}; chunk_size = 0, max_buffered = 360) where T

  # Create and test handle
  handle = lib.lsl_create_outlet(info, chunk_size, max_buffered)
  if handle == C_NULL
    error("liblsl library returned NULL pointer during StreamOutlet creation")
  end

  # Create an internal information object
  info = lib.lsl_get_info(handle)

  return StreamOutlet(handle, StreamInfo(info, T))
end

#
# Push sample
#

# Type mappings
const _lsl_typestring_map = Dict{DataType, String}(Float32 => "f",
                                                   Float64 => "d",
                                                   Clong => "l",
                                                   Int32 => "i",
                                                   Int16 => "s",
                                                   Cchar => "c",
                                                   String => "str",
                                                   Cvoid => "v")

_lsl_push_sample_tp(o, d::Vector{Float32}, ts, pt)  = lib.lsl_push_sample_ftp(o, d, ts, pt)
_lsl_push_sample_tp(o, d::Vector{Float64}, ts, pt)  = lib.lsl_push_sample_dtp(o, d, ts, pt)
@static if !Sys.iswindows() && Sys.WORD_SIZE == 64
  _lsl_push_sample_tp(o, d::Vector{Clong}, ts, pt)    = lib.lsl_push_sample_ltp(o, d, ts, pt)
end
_lsl_push_sample_tp(o, d::Vector{Int32}, ts, pt)    = lib.lsl_push_sample_itp(o, d, ts, pt)
_lsl_push_sample_tp(o, d::Vector{Int16}, ts, pt)    = lib.lsl_push_sample_stp(o, d, ts, pt)
_lsl_push_sample_tp(o, d::Vector{Cchar}, ts, pt)    = lib.lsl_push_sample_ctp(o, d, ts, pt)
_lsl_push_sample_tp(o, d::Vector{Cvoid}, ts, pt)    = lib.lsl_push_sample_vtp(o, d, ts, pt)
_lsl_push_sample_tp(o, d::Vector{String}, ts, pt)   = lib.lsl_push_sample_strtp(o, d, ts, pt)

 """
    push_sample(outlet::StreamOutlet{T},
                data::Vector{T};
                timestamp = 0.0,
                pushthrough = true)

Push a sample into the outlet. Each entry in the vector corresponds to one channel.

# Keyword arguments
-`timestamp::Number`: Capture time of the sample, in agreement with `local_clock()`, if
                      ommitted, the current time is used.
- `passthrough::Bool`: Whether to push the sample through to the receivers instead of
                       buffering it with subsequent samples. Note that the chunk_size, if
                       specified at outlet construction, takes precedence over the
                       pushthrough flag.
"""
function push_sample(outlet::StreamOutlet{T},
                     data::Vector{T};
                     timestamp = 0.0,
                     pushthrough = true) where T

  length(data) == channel_count(outlet.info) || error("data length ≂̸ channel count")
  handle_error(_lsl_push_sample_tp(outlet, data, timestamp, pushthrough))
end

#
# Push chunk
#

# Type mappings
_lsl_push_chunk_tp(o, d::Matrix{Float32}, ts::Number, pt) = lib.lsl_push_chunk_ftp(o, d, length(d), ts, pt)
_lsl_push_chunk_tp(o, d::Matrix{Float64}, ts::Number, pt) = lib.lsl_push_chunk_dtp(o, d, length(d), ts, pt)
@static if !Sys.iswindows() && Sys.WORD_SIZE == 64
  _lsl_push_chunk_tp(o, d::Matrix{Clong}, ts::Number, pt)   = lib.lsl_push_chunk_ltp(o, d, length(d), ts, pt)
end
_lsl_push_chunk_tp(o, d::Matrix{Int32}, ts::Number, pt)   = lib.lsl_push_chunk_itp(o, d, length(d), ts, pt)
_lsl_push_chunk_tp(o, d::Matrix{Int16}, ts::Number, pt)   = lib.lsl_push_chunk_stp(o, d, length(d), ts, pt)
_lsl_push_chunk_tp(o, d::Matrix{Cchar}, ts::Number, pt)   = lib.lsl_push_chunk_ctp(o, d, length(d), ts, pt)
_lsl_push_chunk_tp(o, d::Matrix{Cvoid}, ts::Number, pt)   = lib.lsl_push_chunk_vtp(o, d, length(d), ts, pt)
#_lsl_push_chunk_tp(o, d::Vector{String}, ts::Number, pt)  = lsl_push_chunk_strtp(o, d, length(d), ts, pt)

_lsl_push_chunk_tnp(o, d::Matrix{Float32}, ts::Vector, pt) = lib.lsl_push_chunk_ftnp(o, d, length(d), ts, pt)
_lsl_push_chunk_tnp(o, d::Matrix{Float64}, ts::Vector, pt) = lib.lsl_push_chunk_dtnp(o, d, length(d), ts, pt)
@static if !Sys.iswindows() && Sys.WORD_SIZE == 64
  _lsl_push_chunk_tnp(o, d::Matrix{Clong}, ts::Vector, pt)   = lib.lsl_push_chunk_ltnp(o, d, length(d), ts, pt)
end
_lsl_push_chunk_tnp(o, d::Matrix{Int32}, ts::Vector, pt)   = lib.lsl_push_chunk_itnp(o, d, length(d), ts, pt)
_lsl_push_chunk_tnp(o, d::Matrix{Int16}, ts::Vector, pt)   = lib.lsl_push_chunk_stnp(o, d, length(d), ts, pt)
_lsl_push_chunk_tnp(o, d::Matrix{Cchar}, ts::Vector, pt)   = lib.lsl_push_chunk_ctnp(o, d, length(d), ts, pt)
_lsl_push_chunk_tnp(o, d::Matrix{Cvoid}, ts::Vector, pt)   = lib.lsl_push_chunk_vtnp(o, d, length(d), ts, pt)
#_lsl_push_chunk_tnp(o, d::Vector{String}, ts::Vector, pt)  = lsl_push_chunk_strtnp(o, d, length(d), ts, pt)

"""
push_chunk(outlet::StreamOutlet{T},
           data::Matrix{T};
           timestamp = 0.0,
           pushthrough = true)

Push a chunk of samples into the outlet.

Each sample consists of the column of the matrix `data` such that size(data) = (M,N) where
M is the channel count of the outlet, and N is the number of samples in the chunk.

# Keyword arguments
-`timestamp::Number`: Capture time of the sample, in agreement with `local_clock()`, if
                  ommitted, the current time is used.
- `passthrough::Bool`: Whether to push the sample through to the receivers instead of
                   buffering it with subsequent samples. Note that the chunk_size, if
                   specified at outlet construction, takes precedence over the
                   pushthrough flag.
"""
function push_chunk(outlet::StreamOutlet{T},
                    data::Matrix{T};
                    timestamp = 0.0,
                    pushthrough = true) where T

  size(data,1) == channel_count(outlet.info) || error("data length ≂̸ channel count")
  handle_error(_lsl_push_chunk_tp(outlet, data, timestamp, pushthrough))
end

"""
push_chunk(outlet::StreamOutlet{T},
           data::Matrix{T},
           timestamp::Vector{Float64};
           pushthrough = true)

Push a chunk of samples into the outlet with individual timestamps.

Each sample consists of the column of the matrix `data` such that size(data) = (M,N) where
M is the channel count of the outlet, and N is the number of samples in the chunk.

# Arguments
-`timestamps`: Capture time of the sample, in agreement with `local_clock()`, if ommitted,
               the current time is used.

# Keyword arguments
-`passthrough::Bool`: Whether to push the sample through to the receivers instead of
                      buffering it with subsequent samples. Note that the chunk_size, if
                      specified at outlet construction, takes precedence over the
                      pushthrough flag.
"""
function push_chunk(outlet::StreamOutlet{T},
                    data::Matrix{T},
                    timestamp::Vector;
                    pushthrough = true) where T

  size(data,2) == length(timestamp) || error("number of timestamps != chunk_count")
  size(data,1) == channel_count(outlet.info) || error("data length ≂̸ channel count")
  handle_error(_lsl_push_chunk_tnp(outlet, data, timestamp, pushthrough))
end

# """
# push_chunk(outlet::StreamOutlet{T},
#            data::Vector{Vector{T}};
#            timestamp = 0.0,
#            pushthrough = true)

# Push a vector of samples into the outlet.

# Each entry in the vector corresponds sample, and each entry in the underlying vector
# corresponds to one channel.

# Note that this function must allocate a new vector to concatenate the individual samples, it
# is preferable to specify the chunks as a matrix of N columns, where there are N samples in
# the chunk.

# # Keyword arguments
# -`timestamp::Number`: Capture time of the most recent sample, in agreement with
#                       local_clock(); if omitted, the current time is used. The time stamps
#                       of other samples are automatically derived according to the sampling
#                       rate of the stream.
# - `passthrough::Bool`: Whether to push the sample through to the receivers instead of
#                        buffering it with subsequent samples. Note that the chunk_size, if
#                        specified at outlet construction, takes precedence over the
#                        pushthrough flag.
# """
# function push_chunk(outlet::StreamOutlet{T},
#                     data::Vector{Vector{T}};
#                     timestamp = 0.0,
#                     pushthrough = true) where T

#   # Check all samples are the same, and the correct length
#   if !all(length(d) == channel_count(outlet.info) for d in data)
#     error("at least one data length ≂̸ channel count")
#   end
#   dmat = hcat(data...)
#   push_chunk(outlet, dmat, timestamp = timestamp, passthrough = passthrough)
# end

#
# Utility functions
#

"""
    have_consumers(outlet::StreamOutlet)

Return true if consumers are currently registered, false otherwise.
"""
have_consumers(outlet::StreamOutlet) = lib.lsl_have_consumers(outlet) > 0 ? true : false

"""
    wait_for_consumers(outlet::SteamOutlet)

Wait until some consumer shows up (without wasting resources).

Returns true if succesful, false if timeout ocurred.
"""
function wait_for_consumers(outlet::StreamOutlet, timeout)
  status = lib.lsl_wait_for_consumers(outlet, timeout)
  return status == 1 ? true : false
end

"""
    info(outlet::StreamOutlet)

Return a StreamInfo record provided by the outlet.

This is what was used to create the stream (and also has the Additional Network Information
fields assigned).
"""
function info(outlet::StreamOutlet{T}) where T
  handle = lib.lsl_get_info(outlet)
  return StreamInfo(handle, T)
end


text:
using DataFrames, CSV, Plots

dir = "../OCTHaGOn_benchmarks/gams/global/"
filenames = readdir(dir)
GAMS_DIR = dir

pd = DataFrame()
n_max = 100
pd = CSV.read(dir * "problem_stats.csv", DataFrame)
pd = pd[pd.all_bounded .>= 0.5, :]
pd = pd[pd.n_vars .<= n_max, :]
pd = pd[1:80, :]

globalresults_idxs = [i for i = 1:80 if !isinf(pd.optimum[i])]

# Plotting the actual optimality gaps (NOTE: perhaps make a log plot?)
ylims = (-0.15, 0.35)
plt = plot(ylims = ylims, xticks = 0:10:80, yticks = minimum(ylims):0.05:maximum(ylims), plot_title	= "Optimality gaps of global benchmarks")
xlabel!("Problem number")
ylabel!("Optimality gap")
for i = 1:length(globalresults_idxs)
    color = "blue"
    if !pd.oct_feas[i]
        color = "red"
    end
    error = 0
    if abs(pd.optimum[i]) <= 1
        error = pd.oct_cost[i] - pd.optimum[i]
    else
        error = (pd.oct_cost[i] - pd.optimum[i]) / pd.optimum[i]
    end
    if error >= maximum(ylims)
        plt = quiver!([i], [maximum(ylims)*7/8], quiver=([0], [maximum(ylims)*1/8]), arrow = true, color = color, linewidth = 1)
    elseif error <= minimum(ylims)
        plt = quiver!([i], [-maximum(ylims)*7/8], quiver=([0], [-maximum(ylims)*1/8]), arrow = true, color = color, linewidth = 1)
    else
        plt = scatter!([i], [error], color = color, legend = false)
    end
end
display(plt)

# # Checking that problems are imported properly...
# valid_filenames = []
# for filename in pd.name
#     filename = filename * ".gms"
#     @info "Trying " * filename * "."
#     model = JuMP.Model()
#     # Parsing GAMS Files
#     lexed = GAMSFiles.lex(GAMS_DIR * filename)
#     gams = []
#     try
#         gams = GAMSFiles.parsegams(GAMS_DIR * filename)
#     catch e
#         if e isa KeyError
#             @warn filename * " failed due to KeyError."
#         elseif e isa LoadError
#             @warn filename * " failed due to an InvalidCharError. "
#         else
#             @warn filename * " failed due to an unknown Error."
#         end
#         continue
#     end
#     GAMSFiles.parseconsts!(gams)

#     vars = GAMSFiles.getvars(gams["variables"])
#     sets = Dict{String, Any}()
#     if haskey(gams, "sets")
#         sets = gams["sets"]
#     end
#     preexprs, bodyexprs = Expr[], Expr[]
#     if haskey(gams, "parameters") && haskey(gams, "assignments")
#         GAMSFiles.parseassignments!(preexprs, gams["assignments"], gams["parameters"], sets)
#     end

#     # Getting variables
#     vardict, constdict = generate_variables!(model, gams) # Actual JuMP variables
#     if length(vardict) <= n_max
#         @info filename * " added with $(length(vardict)) variables."
#         push!(valid_filenames, filename)
#     end
# end


text:
"""
Macro to extract row `i` from matrix `M`, for compatibility across Julia releasesq
"""
macro row(M, i)
    if VERSION < @v_str("0.5-")
        return :( squeeze($M[$i,:], 1) )
    else
        return :( $M[$i,:] )
    end
end


text:
struct DemuxLogger <: AbstractLogger
    loggers::Vector{AbstractLogger}
end

"""
    DemuxLogger(loggers...; include_current=true)

Sets the global_logger to demux,
so that messages are sent to all the loggers.
If `include_current_global` is true, then messages are also sent to the global logger
(or rather to what was the global logger when this was constructed).
Normally this would be the ConsoleLogger in the REPL etc.
"""
function DemuxLogger(loggers::Vararg{AbstractLogger}; include_current_global=true)
    loggers = Vector{AbstractLogger}(collect(loggers))
    if include_current_global
        push!(loggers, global_logger())
    end
    DemuxLogger(loggers)
end

function handle_message(demux::DemuxLogger, args...; kwargs...)
    for logger in demux.loggers
        if comp_handle_message_check(logger, args...; kwargs...)
            handle_message(logger, args...; kwargs...)
        end
    end
end

function shouldlog(demux::DemuxLogger, args...)
    any(comp_shouldlog(logger, args...) for logger in demux.loggers)
end

function min_enabled_level(demux::DemuxLogger)
    minimum(min_enabled_level(logger) for logger in demux.loggers)
end

function catch_exceptions(demux::DemuxLogger)
    any(catch_exceptions(logger) for logger in demux.loggers)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.453305806680525 <= q <= 37.480198521751326)

                     add_NL_constraint(m, :(swish(-0.5060529352779994 + 0.34605798298901425*swish(-0.49627455608242954 + -0.46375151869975495*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.8504066842806357*swish(-0.2321501407718678 + 0.5600091754803458*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.462211279142271*swish(-0.9807289687619161 + -0.5968595287514287*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.0428411186394575*swish(0.7899732054652109 + -0.5363355140855166*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + swish(0.02992678336456356 + -0.8781954986928167*swish(-0.49627455608242954 + -0.46375151869975495*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.84297522225665*swish(-0.2321501407718678 + 0.5600091754803458*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.44937020244474857*swish(-0.9807289687619161 + -0.5968595287514287*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.012040065223148932*swish(0.7899732054652109 + -0.5363355140855166*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + swish(-0.25761034404046734 + -0.19911102915318102*swish(-0.49627455608242954 + -0.46375151869975495*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.7457041629864394*swish(-0.2321501407718678 + 0.5600091754803458*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.7638178262546895*swish(-0.9807289687619161 + -0.5968595287514287*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.07889233901952375*swish(0.7899732054652109 + -0.5363355140855166*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + swish(0.1153458350428731 + 0.9297971624965538*swish(-0.49627455608242954 + -0.46375151869975495*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.7017891371439648*swish(-0.2321501407718678 + 0.5600091754803458*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.8093911002605614*swish(-0.9807289687619161 + -0.5968595287514287*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.43763294749909276*swish(0.7899732054652109 + -0.5363355140855166*swish(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*swish(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*swish(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*swish(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("DataFS.jl")
include("ConstantsFS.jl")
include("StateFS.jl")
include("TunersFS.jl")

include("update_p_feature_select.jl")
include("update_omega_feature_select.jl")
include("update_W_feature_select.jl")
include("update_r_feature_select.jl")
include("update_lam_feature_select.jl")

include("update_feature_select.jl")
include("fit_feature_select.jl")


text:
module CDLasso

import Base.LinAlg.BLAS
import Base.LinAlg.SVD

export lasso!, LassoOptions

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

type LassoOptions
  max_inner_iter::Int64
  max_outer_iter::Int64
  xtol::Float64
  gradtol::Float64
  zerothr::Float64
  printEvery::Int64
end

function LassoOptions(;max_inner_iter::Integer     = 200,
                      max_outer_iter::Integer      = 200,
                      xtol::Float64                = 1e-6,
                      gradtol::Float64             = 1e-6,
                      zerothr::Float64             = 1e-4,
                      printEvery::Integer          = -1
                      )
  LassoOptions(
    max_inner_iter,
    max_outer_iter,
    xtol,
    gradtol,
    zerothr,
    printEvery
    )
end

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

include("lasso.jl")
include("glasso.jl")

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

function _penalty_loadings!{T<:AbstractFloat}(
    γ::StridedVector{T},
    λ::StridedVector{T},
    X::StridedMatrix{T},
    res::StridedVector{T}
    )
  n, p = size(X)
  for j=1:p
    s = zero(T)
    m = zero(T)
    for i=1:n
      s += X[i, j]^2. * res[i]^2.
      m += X[i, j] * res[i]
    end
    γ[j] = sqrt(s / n - (m/n)^2.) * λ[j]
  end
end

# numIterations that estimate Gamma matrix
function feasible_lasso!{T<:AbstractFloat}(
    x::StridedVector{T},
    X::StridedMatrix{T},
    Y::StridedVector{T},
    λ::Vector{T};
    options::LassoOptions   =   LassoOptions(),
    refit::Bool             =   true,
    numIterations::Int      =   1,
    xtol::T                 =   1e-4
    )

  n, p = size(X)
  A = X'X
  scale!(A, 1./n)
  b = X'Y
  scale!(b, -1./n)

  xold = copy(x)
  γ = copy(λ)
  res = copy(Y)

  BLAS.gemv!('N', -1., X, x, 1., res)
  _penalty_loadings!(γ, λ, X, res)

  for iter=1:numIterations
    _, activeset = lasso!(x, A, b, γ)

    # check if done
    s = zero(T)
    for j=1:p
      s = max(s, abs(x[j] - xold[j]))
    end
    if s < xtol
      break
    end

    if refit
      indNZ = activeset.indexes[1:activeset.numActive]
      #       xv = sub(x, indNZ)
      #       Av = sub(A, indNZ, indNZ)
      #       bv = sub(b, indNZ)
      #       xv[:] = Av \ bv
      #       scale!(xv, -1)
      x[indNZ] = X[:, indNZ] \ Y[:]
    end

    # update γ
    if iter < numIterations
      # re-estimate gamma unless last iteration
      copy!(res, Y)
      BLAS.gemv!('N', -1., X, x, 1., res)
      _penalty_loadings!(γ, λ, X, res)

      copy!(xold, x)
    end
    #
  end
  x
end

end

max_stars_repo_name: mlakolar/CDLasso

max_issues_repo_name: mlakolar/CDLasso

max_forks_repo_name: mlakolar/CDLasso


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module EA

using ExtractMacro
using ..Interface
using ..Common
using ..DFloats

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphEA, GraphEANormalDiscretized, GraphEANormal

import ..Interface: energy, delta_energy, neighbors, allΔE,
                    update_cache!, update_cache_residual!

import ..DFloats: MAXDIGITS
sentinel(::Type{ET}) where {ET} = typemin(ET)

discr(::Type{ET}, x::Real) where {ET} = convert(ET, round(x, digits=MAXDIGITS))
discr(::Type{DFloat64}, x::Real) = x
discr(::Type{ET}, x::Integer) where {ET<:Integer} = convert(ET, x)

function gen_EA(L::Integer, D::Integer)
    L ≥ 2 || throw(ArgumentError("L must be ≥ 2, given: $L"))
    D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
    N = L^D
    A = [Int[] for x = 1:N]
    dims = ntuple(d->L, D)

    for cind in CartesianIndices(dims)
        x = LinearIndices(dims)[cind.I...]
        for d in 1:D
            I1 = ntuple(k->(k==d ? mod1(cind.I[k] + 1, L) : cind.I[k]), D)
            y = LinearIndices(dims)[I1...]
            push!(A[x], y)
            push!(A[y], x)
        end
    end
    map!(sort!, A, A)
    tA = NTuple{2D,Int}[tuple(Ax...) for Ax in A]
    return tA
end

function gen_J(f, ET::Type, N::Integer, A::Vector{NTuple{twoD,Int}}) where {twoD}
    @assert all(issorted, A)
    m = sentinel(ET)
    J = Vector{ET}[zeros(ET,twoD) for i = 1:N]
    map(Jx->fill!(Jx, m), J)
    for x = 1:N
        Jx = J[x]
        Ax = A[x]
        for k = 1:length(Ax)
            y = Ax[k]
            if x < y
                Jxy = f()
                @assert Jx[k] == m
                Jx[k] = Jxy
                J[y][findfirst(==(m), J[y])] = Jxy
            # else # this check fails for L=2
            #     l = findfirst(A[y], x)
            #     Jxy = J[y][l]
            #     @assert Jxy ≠ m
            #     @assert J[x][k] == Jxy
            end
        end
    end
    @assert all(Jx->all(Jx .≠ m), J)
    tJ = NTuple{twoD,ET}[tuple(Jx...) for Jx in J]
    return tJ
end

function gen_AJ(fname::AbstractString)
    D = 2
    twoD = 2D
    f = open(fname)
    local N, L, A, J, tJ
    try
        @assert startswith(strip(readline(f)), "type:")
        ls = split(readline(f))
        @assert length(ls) == 2
        @assert ls[1] == "size:"
        L = parse(Int, ls[2])
        @assert startswith(strip(readline(f)), "name:")
        A = EA.gen_EA(L, D)
        N = length(A)

        m = EA.sentinel(Float64)
        J = Vector{Float64}[zeros(twoD) for i = 1:N]
        map(Jx->fill!(Jx, m), J)

        for l in eachline(f)
            ls = split(l)
            @assert length(ls) == 3
            x, y, Jxy = parse(Int, ls[1]), parse(Int, ls[2]), parse(Float64, ls[3])

            Jx = J[x]
            Ax = A[x]
            k = findfirst(Ax, y)
            @assert k ≠ 0
            @assert Jx[k] == m
            Jx[k] = Jxy

            Jy = J[y]
            Ay = A[y]
            k = findfirst(Ay, x)
            @assert k ≠ 0
            @assert Jy[k] == m
            Jy[k] = Jxy
        end
        @assert all(Jx->all(Jx .≠ m), J)
    finally
        close(f)
    end
    tJ = NTuple{twoD,Float64}[tuple(Jx...) for Jx in J]

    return L, D, A, tJ
end

function get_vLEV(LEV, ET::Type)
    LEV isa Tuple{Real,Vararg{Real}} || throw(ArgumentError("invalid level spec, expected a Tuple of Reals, given: $LEV"))
    length(unique(LEV)) == length(LEV) || throw(ArgumentError("repeated levels in LEV: $LEV"))

    vLEV = Array{ET}(undef, length(LEV))
    try
        for i = 1:length(LEV)
            vLEV[i] = LEV[i]
        end
    catch
        throw(ArgumentError("incompatible energy type and level spec, conversion failed, given: ET=$ET LEV=$LEV"))
    end
    m = sentinel(ET)
    m ∈ vLEV && throw(ArgumentError("illegal level value: $m"))
    any(x->discr(ET, x) ≠ x, vLEV) && throw(ArgumentError("up to $MAXDIGITS decimal digits supported in levels, given: $LEV"))
    return vLEV
end

mutable struct GraphEA{ET,LEV,twoD} <: DiscrGraph{ET}
    N::Int
    #L::Int
    A::Vector{NTuple{twoD,Int}}
    J::Vector{NTuple{twoD,ET}}
    uA::Vector{Vector{Int}}
    cache::LocalFields{ET}
    function GraphEA{ET,LEV,twoD}(A::Vector{NTuple{twoD,Int}}, J::Vector{NTuple{twoD,ET}}) where {ET,LEV,twoD}
        isa(twoD, Integer) || throw(ArgumentError("twoD must be integer, given a: $(typeof(twoD))"))
        iseven(twoD) || throw(ArgumentError("twoD must be even, given: $twoD"))
        D = twoD ÷ 2
        D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))

        N = length(A)
        N ≥ 2 || throw(ArgumentError("invalid A, expected length ≥ 2, found: $N"))
        all(a->length(a) == twoD, A) || throw(ArgumentError("invalid A inner length, expected $twoD, given: $(unique(map(length,A)))"))

        isL2 = (A[1][1] == A[1][2])
        all(a->(issorted(a) && unique(a) == collect(a[1:(1+isL2):end])), A) || throw(ArgumentError("invalid A, does not look like an EA graph"))

        uA = [collect(a[1:(1+isL2):end]) for a in A]

        vLEV = get_vLEV(LEV, ET)
        all(Jx->all(Jxy->Jxy ∈ vLEV, Jx), J) || throw(ArgumentError("the given J is incompatible with levels $LEV"))
        length(J) == N || throw(ArgumentError("incompatible lengths of A and J: $(length(A)), $(length(J))"))
        all(j->length(j) == twoD, J) || throw(ArgumentError("invalid J inner length, expected $twoD, given: $(unique(map(length,J)))"))

        cache = LocalFields{ET}(N)

        return new{ET,LEV,twoD}(N, A, J, uA, cache)
    end
end

"""
    GraphEA(L::Integer, D::Integer, LEV = (-1,1)) <: DiscrGraph

An Edwards-Anderson `DiscrGraph`: spins are arranged on a square lattice of size `L`
in `D` dimensions (i.e. there are \$L^D\$ total spins), with periodic boundary
conditions.

The interactions are extracted at random from `LEV`, which must be a `Tuple` of `Real`s.
No external fields.
"""
function GraphEA(L::Integer, D::Integer, LEV::Tuple{ET,Vararg{ET}}) where {ET<:Real}
    A = gen_EA(L, D)
    vLEV = get_vLEV(LEV, ET)
    N = length(A)
    J = gen_J(ET, N, A) do
        rand(vLEV)
    end
    return GraphEA{ET,LEV,2D}(A, J)
end
GraphEA(L::Integer, D::Integer, LEV::Tuple{Real,Vararg{Real}}) = GraphEA(L, D, promote(LEV...))
GraphEA(L::Integer, D::Integer) = GraphEA(L, D, (-1,1))

GraphEA(L::Integer, D::Integer, LEV::Tuple{Float64,Vararg{Float64}}) = GraphEA(L, D, map(DFloat64, LEV))

function energy(X::GraphEA{ET}, C::Config) where {ET}
    @assert X.N == C.N
    @extract C : s
    @extract X : A J cache
    @extract cache : lfields lfields_last
    n = zero(ET)
    for x = 1:length(A)
        Jx = J[x]
        σx = 2 * s[x] - 1
        Ax = A[x]
        lf = zero(ET)
        for k = 1:length(Ax)
            y = Ax[k]
            σy = 2 * s[y] - 1
            Jxy = Jx[k]

            lf -= Jxy * σx * σy
        end
        n += lf
        lfields[x] = discr(ET, 2lf)
    end
    #@assert n % 2 == 0
    #n ÷= 2
    n /= 2
    cache.move_last = 0
    fill!(lfields_last, zero(ET))
    return discr(ET, n)
end

function update_cache!(X::GraphEA{ET}, C::Config, move::Int) where {ET}
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : A uA J cache

    @extract cache : lfields lfields_last move_last
    if move_last == move
        @inbounds begin
            Ux = uA[move]
            for y in Ux
                lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
            end
            lfields[move] = -lfields[move]
            lfields_last[move] = -lfields_last[move]
        end
        return
    end

    @inbounds begin
        Ux = uA[move]
        for y in Ux
            lfields_last[y] = lfields[y]
        end
        Jx = J[move]
        sx = s[move]
        Ax = A[move]
        for k = 1:length(Ax)
            y = Ax[k]
            σxy = 1 - 2 * (sx ⊻ s[y])
            Jxy = Jx[k]
            lfields[y] = discr(ET, lfields[y] - 4 * σxy * Jxy)
        end
        lfm = lfields[move]
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache.move_last = move

    return
end

function delta_energy(X::GraphEA{ET}, C::Config, move::Int) where {ET}
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    #@extract C : s
    #@extract X : A J
    @extract X : cache
    @extract cache : lfields

    @inbounds Δ = -lfields[move]
    return Δ

    # @inbounds begin
    #     Δ = zero(ET)
    #     Jx = J[move]
    #     σx = 2 * s[move] - 1
    #     Ax = A[move]
    #     for k = 1:length(Ax)
    #         y = Ax[k]
    #         σy = 2 * s[y] - 1
    #         Jxy = Jx[k]
    #         Δ += 2 * Jxy * σx * σy
    #     end
    # end
    # return discr(ET, Δ)
end

neighbors(X::GraphEA, i::Int) = return X.uA[i]
@generated allΔE(::Type{GraphEA{Int,(-1,1),twoD}}) where {twoD} = Expr(:tuple, ntuple(d1->(4 * (d1 - 1)), (twoD÷2)+1)...)

@generated function allΔE(::Type{GraphEA{ET,LEV,twoD}}) where {ET,LEV,twoD}
    list = Set{ET}()
    L = length(LEV)
    es = Set{ET}(zero(ET))
    for n = 1:twoD
        newes = Set{ET}()
        for n in es, l in LEV
            push!(newes, discr(ET, n + l))
            push!(newes, discr(ET, n - l))
        end
        es = newes
    end
    deltas = sort!(unique((x->2 * abs(x)).(collect(es))))
    return Expr(:tuple, deltas...)
end

mutable struct GraphEANormalDiscretized{ET,LEV,twoD} <: DoubleGraph{DiscrGraph{ET},Float64}
    N::Int
    X0::GraphEA{ET,LEV,twoD}
    A::Vector{NTuple{twoD,Int}}
    rJ::Vector{NTuple{twoD,Float64}}
    uA::Vector{Vector{Int}}
    cache::LocalFields{Float64}
    function GraphEANormalDiscretized{ET,LEV,twoD}(L::Integer) where {ET,LEV,twoD}
        isa(twoD, Integer) || throw(ArgumentError("twoD must be integer, given a: $(typeof(twoD))"))
        iseven(twoD) || throw(ArgumentError("twoD must be even, given: $twoD"))
        D = twoD ÷ 2
        D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
        A = gen_EA(L, D)
        N = length(A)
        cJ = gen_J(Float64, N, A) do
            randn()
        end

        dJ = Array{NTuple{twoD,ET}}(undef, N)
        rJ = Array{NTuple{twoD,Float64}}(undef, N)
        for (x, cJx) in enumerate(cJ)
            dJ[x], rJ[x] = discretize(cJx, LEV)
        end

        X0 = GraphEA{ET,LEV,twoD}(A, dJ)
        cache = LocalFields{Float64}(N)
        return new{ET,LEV,twoD}(N, X0, A, rJ, X0.uA, cache)
    end
end

"""
    GraphEANormalDiscretized(L::Integer, D::Integer, LEV) <: DoubleGraph{DiscrGraph,Float64}

An Edwards-Anderson `DoubleGraph`: spins are arranged on a square lattice of size `L`
in `D` dimensions (i.e. there are \$L^D\$ total spins), with periodic boundary
conditions.

The interactions are extracted at random from a normal distribution
with unit variance, and are then discretized using the values in `LEV`,
which must be a `Tuple` of `Real`s. No external fields.

Same as [`GraphEANormal`](@ref), but works differently when used with [`rrrMC`](@ref).
"""
GraphEANormalDiscretized(L::Integer, D::Integer, LEV::Tuple{ET,Vararg{ET}}) where {ET<:Real} = GraphEANormalDiscretized{ET,LEV,2D}(L)
GraphEANormalDiscretized(L::Integer, D::Integer, LEV::Tuple{Real,Vararg{Real}}) = GraphEANormalDiscretized(L, D, promote(LEV...))

GraphEANormalDiscretized(L::Integer, D::Integer, LEV::Tuple{Float64,Vararg{Float64}}) = GraphEANormalDiscretized{DFloat64,map(DFloat64,LEV),2D}(L)

function energy(X::GraphEANormalDiscretized, C::Config)
    @assert X.N == C.N
    @extract X : X0 A J=rJ cache
    @extract C : N s
    @extract cache : lfields lfields_last

    E0 = energy(X0, C)

    E1 = 0.0
    for x = 1:length(A)
        Jx = J[x]
        σx = 2 * s[x] - 1
        Ax = A[x]
        lf = 0.0
        for k = 1:length(Ax)
            y = Ax[k]
            σy = 2 * s[y] - 1
            Jxy = Jx[k]

            lf -= Jxy * σx * σy
        end
        E1 += lf
        lfields[x] = 2lf
    end
    E1 /= 2
    cache.move_last = 0
    fill!(lfields_last, 0.0)

    return convert(Float64, E0 + E1)
end

function update_cache!(X::GraphEANormalDiscretized{ET}, C::Config, move::Int) where {ET}
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N

    if X.cache.move_last ≠ X.X0.cache.move_last
        update_cache!(X.X0, C, move)
        update_cache_residual!(X, C, move)
        return
    end

    @extract C : N s
    @extract X : A uA rJ X0 cache
    @extract X0 : J0=J cache0=cache

    @extract cache : lfields lfields_last move_last
    @extract cache0 : lfields0=lfields lfields_last0=lfields_last move_last0=move_last
    @assert move_last0 == move_last
    if move_last == move
        @inbounds begin
            Ux = uA[move]
            for y in Ux
                lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
                lfields0[y], lfields_last0[y] = lfields_last0[y], lfields0[y]
            end
            lfields[move] = -lfields[move]
            lfields_last[move] = -lfields_last[move]
            lfields0[move] = -lfields0[move]
            lfields_last0[move] = -lfields_last0[move]
        end
        return
    end

    @inbounds begin
        Ux = uA[move]
        for y in Ux
            lfields_last[y] = lfields[y]
            lfields_last0[y] = lfields0[y]
        end
        Jx0 = J0[move]
        Jx = rJ[move]
        sx = s[move]
        Ax = A[move]
        for k = 1:length(Ax)
            y = Ax[k]
            σxy = 1 - 2 * (sx ⊻ s[y])

            Jxy0 = Jx0[k]
            lfields0[y] = discr(ET, lfields0[y] - 4 * σxy * Jxy0)

            Jxy = Jx[k]
            lfields[y] -= 4 * σxy * Jxy
        end
        lfm0 = lfields0[move]
        lfields_last0[move] = lfm0
        lfields0[move] = -lfm0

        lfm = lfields[move]
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache0.move_last = move
    cache.move_last = move

    return
end

function update_cache_residual!(X::GraphEANormalDiscretized, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : A uA J=rJ cache

    @extract cache : lfields lfields_last move_last
    if move_last == move
        @inbounds begin
            Ux = uA[move]
            for y in Ux
                lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
            end
            lfields[move] = -lfields[move]
            lfields_last[move] = -lfields_last[move]
        end
        return
    end

    @inbounds begin
        Ux = uA[move]
        for y in Ux
            lfields_last[y] = lfields[y]
        end
        Jx = J[move]
        sx = s[move]
        Ax = A[move]
        for k = 1:length(Ax)
            y = Ax[k]
            σxy = 1 - 2 * (sx ⊻ s[y])
            Jxy = Jx[k]
            lfields[y] -= 4 * σxy * Jxy
        end
        lfm = lfields[move]
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache.move_last = move

    return
end

function delta_energy_residual(X::GraphEANormalDiscretized, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    #@extract C : s
    @extract X : cache
    @extract cache : lfields

    @inbounds Δ = -lfields[move]
    return Δ

    # @inbounds begin
    #     Δ = 0.0
    #     Jx = J[move]
    #     σx = 2 * s[move] - 1
    #     Ax = A[move]
    #     for k = 1:length(Ax)
    #         y = Ax[k]
    #         σy = 2 * s[y] - 1
    #         Jxy = Jx[k]
    #         Δ += 2 * Jxy * σx * σy
    #     end
    # end
    # return Δ
end

function delta_energy(X::GraphEANormalDiscretized, C::Config, move::Int)
    ΔE0 = delta_energy(X.X0, C, move)
    ΔE1 = delta_energy_residual(X, C, move)
    return convert(Float64, ΔE0 + ΔE1)
end

neighbors(X::GraphEANormalDiscretized, i::Int) = return X.uA[i]

## GraphEANormal

mutable struct GraphEANormal{twoD} <: SimpleGraph{Float64}
    N::Int
    A::Vector{NTuple{twoD,Int}}
    J::Vector{NTuple{twoD,Float64}}
    uA::Vector{Vector{Int}}
    cache::LocalFields{Float64}
    function GraphEANormal{twoD}(L::Integer, A, J) where {twoD}
        isa(twoD, Integer) || throw(ArgumentError("twoD must be integer, given a: $(typeof(twoD))"))
        iseven(twoD) || throw(ArgumentError("twoD must be even, given: $twoD"))
        D = twoD ÷ 2
        D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
        N = length(A)
        # TODO: check A and J

        uA = [unique(a) for a in A] # needed for the case L=2

        cache = LocalFields{Float64}(N)
        return new{twoD}(N, A, J, uA, cache)
    end
end

"""
    GraphEANormal(L::Integer, D::Integer) <: SimpleGraph{Float64}

An Edwards-Anderson `SimpleGraph`: spins are arranged on a square lattice of size `L`
in `D` dimensions (i.e. there are \$L^D\$ total spins), with periodic boundary
conditions.

Same as [`GraphEA`](@ref), but the interactions are extracted from a normal distribution
with unit variance.
"""
function GraphEANormal(L::Integer, D::Integer; genJf=randn)
    D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
    A = gen_EA(L, D)
    N = length(A)
    J = gen_J(Float64, N, A) do
        genJf()
    end

    return GraphEANormal{2D}(L, A, J)
end

function GraphEANormal(fname::AbstractString)
    L, D, A, J = EA.gen_AJ(fname)
    N = length(A)
    @assert N == L^D
    return GraphEANormal{2D}(L, A, J)
end

function energy(X::GraphEANormal, C::Config)
    @assert X.N == C.N
    @extract X : A J cache
    @extract C : N s
    @extract cache : lfields lfields_last

    E1 = 0.0
    for x = 1:length(A)
        Jx = J[x]
        σx = 2 * s[x] - 1
        Ax = A[x]
        lf = 0.0
        for k = 1:length(Ax)
            y = Ax[k]
            σy = 2 * s[y] - 1
            Jxy = Jx[k]

            lf -= Jxy * σx * σy
        end
        E1 += lf
        lfields[x] = 2lf
    end
    E1 /= 2
    cache.move_last = 0
    fill!(lfields_last, 0.0)

    return E1
end

function update_cache!(X::GraphEANormal, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : A uA J cache

    @extract cache : lfields lfields_last move_last
    if move_last == move
        @inbounds begin
            Ux = uA[move]
            for y in Ux
                lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
            end
            lfields[move] = -lfields[move]
            lfields_last[move] = -lfields_last[move]
        end
        return
    end

    @inbounds begin
        Ux = uA[move]
        for y in Ux
            lfields_last[y] = lfields[y]
        end
        Jx = J[move]
        sx = s[move]
        Ax = A[move]
        for k = 1:length(Ax)
            y = Ax[k]
            σxy = 1 - 2 * (sx ⊻ s[y])
            Jxy = Jx[k]
            lfields[y] -= 4 * σxy * Jxy
        end
        lfm = lfields[move]
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache.move_last = move

    return
end

function delta_energy(X::GraphEANormal, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    #@extract C : s
    @extract X : cache
    @extract cache : lfields

    @inbounds Δ = -lfields[move]
    return Δ

    # @inbounds begin
    #     Δ = 0.0
    #     Jx = J[move]
    #     σx = 2 * s[move] - 1
    #     Ax = A[move]
    #     for k = 1:length(Ax)
    #         y = Ax[k]
    #         σy = 2 * s[y] - 1
    #         Jxy = Jx[k]
    #         Δ += 2 * Jxy * σx * σy
    #     end
    # end
    # return Δ
end

neighbors(X::GraphEANormal, i::Int) = return X.uA[i]

end


text:
abstract type DCLine <: Branch end

struct HVDCLine <: Branch
    name::String
    available::Bool
    connectionpoints::From_To_Bus
    activepowerlimits_from::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #MW
    activepowerlimits_to::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #MW
    reactivepowerlimits_from::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #MVar
    reactivepowerlimits_to::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #MVar
    loss::NamedTuple{(:l0, :l1),Tuple{Float64,Float64}}
end

HVDCLine(; name ="init",
        available = true,
        connectionpoints = (from = Bus(), to = Bus()),
        activepowerlimits_from = (min=0.0, max=0.0),
        activepowerlimits_to = (min=0.0, max=0.0),
        reactivepowerlimits_from = (min=0.0, max=0.0),
        reactivepowerlimits_to = (min=0.0, max=0.0),
        loss = (l0=0.0, l1=0.0)
    ) = HVDCLine(name, available, connectionpoints, activepowerlimits_from, activepowerlimits_to, reactivepowerlimits_from, reactivepowerlimits_to,loss )

"""
As implemented in Milano's Book Page 397
"""
struct VSCDCLine <: DCLine
    name::String
    available::Bool
    connectionpoints::From_To_Bus
    rectifier_taplimits::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #pu
    rectifier_xrc::Float64
    rectifier_firingangle::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #radians
    inverter_taplimits::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #pu
    inverter_xrc::Float64
    inverter_firingangle::NamedTuple{(:min, :max),Tuple{Float64,Float64}} #radians
end

VSCDCLine(; name ="init",
        available = true,
        connectionpoints = (from = Bus(), to = Bus()),
        rectifier_taplimits = (min=0.0, max=0.0),
        rectifier_xrc = 0.0,
        rectifier_firingangle = (min=0.0, max=0.0),
        inverter_taplimits = (min=0.0, max=0.0),
        inverter_xrc = 0.0,
        inverter_firingangle = (min=0.0, max=0.0),
    ) = VSCDCLine(name, available, connectionpoints,rectifier_taplimits, rectifier_xrc, rectifier_firingangle, inverter_taplimits, inverter_xrc, inverter_firingangle)


text:
module Nettle
using Nettle_jll, Libdl

include("hash_common.jl")
include("hash.jl")
include("hmac.jl")
include("cipher.jl")

# SnoopCompile acceleration
include("precompile.jl")
_precompile_()

end # module