text:
function setup_etc(backend::GPUBackend, layer::BinaryAccuracyLayer, inputs)
  etc = make_blob(backend, Float32, 1)
  return etc
end
function shutdown(backend::GPUBackend, state::BinaryAccuracyLayerState)
  destroy(state.etc)
end

function forward(backend::GPUBackend, state::BinaryAccuracyLayerState, inputs::Vector{Blob})
  pred = inputs[1]
  label = inputs[2]

  const N = length(pred)
  const x_block = div(N+CUDA.THREADS_PER_BLOCK_X-1, CUDA.THREADS_PER_BLOCK_X)

  const data_type = eltype(pred)
  if data_type == Float32
    kernel = backend.mocha.binary_accuracy_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.binary_accuracy_forward_double
  else
    error("Unsupported data type $data_type")
  end

  threshold = convert(data_type, state.layer.threshold)
  erase!(state.etc)
  CUDA.launch(kernel, (x_block,1),(CUDA.THREADS_PER_BLOCK_X,1),
      (pred.ptr.p, label.ptr.p, N, threshold, state.etc.ptr.p));

  n_wrong = Float32[0.0f0]
  copy!(n_wrong, state.etc)
  state.n_wrong += n_wrong[1]

  # accumulate accuracy
  state.n_accum += N
  state.accuracy = (state.n_accum-state.n_wrong) / state.n_accum
end


text:
function linear_interpolation(initial_state, final_state, horizon)
    n = length(initial_state)
    X = [copy(Array(initial_state)) for t = 1:horizon]
    for t = 1:horizon
        for i = 1:n
            X[t][i] = (final_state[i] - initial_state[i]) / (horizon - 1) * (t - 1) + initial_state[i]
        end
    end
    return X
end

ext: jl


text:
using MLBasedESC
using MLBasedESC: posdef, skewsym
using Flux
using LinearAlgebra
import DiffEqFlux
using DiffEqFlux: FastChain, FastDense

const NIN = 2

inmap(q,::Any=nothing) = [cos(q[1]); sin(q[1]); cos(q[2]); sin(q[2])]
function MLBasedESC.jacobian(::typeof(inmap), q,::Any=nothing)
    qbar = inmap(q)
    [
        -qbar[2] 0
         qbar[1] 0
         0 -qbar[4]
         0  qbar[3]
    ]
end

function runalltests()
    c_single_output()
    c_single_output_odd()
    c_single_output_inmap()
    c_single_output_nonneg()
    c_single_output_nonneg_odd()
    c_single_output_nonneg_inmap()
    fc_single_output()
    fc_single_output_odd()
    fc_single_output_inmap()
    fc_single_output_nonneg()
    fc_single_output_nonneg_odd()
    fc_single_output_nonneg_inmap()
    c_multi_output()
    c_multi_output_odd()
    c_multi_output_inmap()
    c_multi_output_nonneg()
    c_multi_output_nonneg_odd()
    c_multi_output_nonneg_inmap()
    fc_multi_output()
    fc_multi_output_odd()
    fc_multi_output_inmap()
    fc_multi_output_nonneg()
    fc_multi_output_nonneg_odd()
    fc_multi_output_nonneg_inmap()
    c_posdef_output()
    c_posdef_output_odd()
    c_posdef_output_inmap()
    c_posdef_output_nonneg()
    c_posdef_output_nonneg_odd()
    c_posdef_output_nonneg_inmap()
    fc_posdef_output()
    fc_posdef_output_odd()
    fc_posdef_output_inmap()
    fc_posdef_output_nonneg()
    fc_posdef_output_nonneg_odd()
    fc_posdef_output_nonneg_inmap()
    c_skewsym_output()
    c_skewsym_output_odd()
    c_skewsym_output_inmap()
    c_skewsym_output_nonneg()
    c_skewsym_output_nonneg_odd()
    c_skewsym_output_nonneg_inmap()
    fc_skewsym_output()
    fc_skewsym_output_odd()
    fc_skewsym_output_inmap()
    fc_skewsym_output_nonneg()
    fc_skewsym_output_nonneg_odd()
    fc_skewsym_output_nonneg_inmap()
end

#==============================================================================
Test 1 ===> Chain w/ single output
==============================================================================#
function c_single_output()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function c_single_output_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_single_output_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function c_single_output_nonneg()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function c_single_output_nonneg_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_single_output_nonneg_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end

#==============================================================================
Test 2 ===> FastChain w/ single output
==============================================================================#
function fc_single_output()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function fc_single_output_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_single_output_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function fc_single_output_nonneg()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function fc_single_output_nonneg_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_single_output_nonneg_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end

#==============================================================================
Test 3 ===> Chain w/ multiple output
==============================================================================#
function c_multi_output()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 10),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function c_multi_output_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 10),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_multi_output_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 10),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function c_multi_output_nonneg()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 10, square),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function c_multi_output_nonneg_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 10, square),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_multi_output_nonneg_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 10, square),
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end

#==============================================================================
Test 4 ===> FastChain w/ multiple output
==============================================================================#
function fc_multi_output()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 10),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function fc_multi_output_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 10),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_multi_output_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 10),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function fc_multi_output_nonneg()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 10, square),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end
function fc_multi_output_nonneg_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 10, square),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_multi_output_nonneg_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 10, square),
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    @assert isapprox(Ja, Jn, atol=1e-6, rtol=1e-6)
end

#==============================================================================
Test 5 ===> Chain w/ posdef output
==============================================================================#
function c_posdef_output()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9), posdef
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function c_posdef_output_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9), posdef
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_posdef_output_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9), posdef
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function c_posdef_output_nonneg()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9, square), posdef,
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function c_posdef_output_nonneg_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9, square), posdef,
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_posdef_output_nonneg_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9, square), posdef,
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end

#==============================================================================
Test 6 ===> FastChain w/ posdef output
==============================================================================#
function fc_posdef_output()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9), posdef
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function fc_posdef_output_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9), posdef
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_posdef_output_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9), posdef
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function fc_posdef_output_nonneg()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9, square), posdef,
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function fc_posdef_output_nonneg_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9, square), posdef,
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_posdef_output_nonneg_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9, square), posdef,
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end

#==============================================================================
Test 7 ===> Chain w/ skewsym output
==============================================================================#
function c_skewsym_output()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9), skewsym
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function c_skewsym_output_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9), skewsym
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_skewsym_output_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9), skewsym
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function c_skewsym_output_nonneg()
    x = rand(NIN)
    C = Chain(
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9, square), skewsym,
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function c_skewsym_output_nonneg_odd()
    x = rand(NIN)
    C = Chain(
        square,
        Dense(NIN, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9, square), skewsym,
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x), C(-x), atol=1e-6, rtol=1e-6)
end
function c_skewsym_output_nonneg_inmap()
    x = rand(NIN)
    C = Chain(
        inmap,
        Dense(4, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 9, square), skewsym,
    )
    Ja = jacobian(C, x)
    Jn = Flux.jacobian(C, x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end

#==============================================================================
Test 8 ===> FastChain w/ skewsym output
==============================================================================#
function fc_skewsym_output()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9), skewsym
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function fc_skewsym_output_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9), skewsym
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_skewsym_output_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9), skewsym
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function fc_skewsym_output_nonneg()
    x = rand(NIN)
    C = FastChain(
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9, square), skewsym,
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end
function fc_skewsym_output_nonneg_odd()
    x = rand(NIN)
    C = FastChain(
        square,
        FastDense(NIN, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9, square), skewsym,
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
    @assert isapprox(C(x,p), C(-x,p), atol=1e-6, rtol=1e-6)
end
function fc_skewsym_output_nonneg_inmap()
    x = rand(NIN)
    C = FastChain(
        inmap,
        FastDense(4, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 9, square), skewsym,
    )
    p = DiffEqFlux.initial_params(C)
    Ja = jacobian(C, x, p)
    Jn = Flux.jacobian(_1->C(_1,p), x)[1]
    for (jn, ja) in zip(eachcol(Jn), Ja)
        @assert isapprox(jn, ja[:], atol=1e-6, rtol=1e-6)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, StatsPlots, Plots, Measures
	Random.seed!(0)
 end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.21"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	b1, b2 = 0.5 , 2
	dist1, dist2, = Beta(b1,b1), Beta(b2,b2)

	n = 2000
	data1 = rand(dist1,n)
	data2 = rand(dist2,n)
end

# ╔═╡ 65067674-250d-11eb-12ab-f7aeb1fb1a16
begin
	stephist(data1, bins=15, label = "beta($b1,$b1)", c = :red, normed = true)
	p1 = stephist!(data2, bins=15, label = "beta($b2,$b2)",
			c = :blue, xlabel="x", ylabel="Density",normed = true)

	p2 = qqplot(data1, data2, c=:black, ms=1, msw =0,
			xlabel="Quantiles for beta($b1,$b1) sample",
			ylabel="Quantiles for beta($b2,$b2) sample",
			legend=false)

	plot(p1, p2, size=(800,400), margin = 5mm)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.21"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═65067674-250d-11eb-12ab-f7aeb1fb1a16
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
#=
    %FILE%
    Copyright © %YEAR% %USER% <%MAIL%>

    Distributed under terms of the %LICENSE% license.
=#

%HERE%


text:
function construct_sample_Problem(k::Union{SDEKernel,SDEKernel!}, u0, Z, alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm})
  @unpack f, g, trange, p, noise_rate_prototype = k

  prob = SDEProblem(f, g, u0, get_tspan(trange), p, noise=Z,
                    noise_rate_prototype=noise_rate_prototype)

  return prob
end

function construct_sample_Problem(k::Union{SDEKernel,SDEKernel!}, u0, Z, alg::AbstractInternalSolver, ll0=nothing)
  @unpack f, g, trange, p, noise_rate_prototype = k

  @assert u0 isa AbstractVector

  Z = compute_Z(Z,noise_rate_prototype,trange,u0)

  if ll0===nothing
    u = (1, trange[1], deepcopy(u0))
    du = (0, trange[1], zero(u0))
  else
    u = (1, trange[1], deepcopy(u0), deepcopy(ll0))
    du = (0, trange[1], zero(u0), zero(ll0))
  end
  dz = (0, trange[1], Z[1])

  return u, dz, du, Z
end

function sample(k::Union{SDEKernel,SDEKernel!}, u0, alg=EM(false), Z=nothing; kwargs...)
  return _sample(k, u0, alg, Z; kwargs...)
end

function _sample(k::Union{SDEKernel,SDEKernel!}, u0, alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}, Z; kwargs...)
  @unpack f, g, trange, p = k

  prob = construct_sample_Problem(k, u0, Z, alg)
  sol = solve(prob, alg, tstops = trange; kwargs...)
  return sol[end], (sol.t, sol.u, sol.W.W)
end

function _sample(k::Union{SDEKernel,SDEKernel!}, u0, alg::AbstractInternalSolver, Z; P=k, save=true, kwargs...)
  @unpack f, g, trange, p = k

  u, dz, du, Z = construct_sample_Problem(k, u0, Z, alg)
  if save
    uu = typeof(u)[]
  else
    uu = nothing
  end
  uu, uT = solve!(alg, uu, u, Z, P)
  if save
    rest = (getindex.(uu,2), getindex.(uu,3), getindex.(Z,3))
  else
    rest = (nothing, nothing, getindex.(Z,3))
  end
  return uT[end], rest
end

# Ensemble Problem
function sample(k::Union{SDEKernel,SDEKernel!}, u0, numtraj::Number, alg=EM(false), Z=nothing, ensemblealg=EnsembleThreads(),
    output_func=(sol,i) -> (sol[end],false); kwargs...)
  @unpack f, g, trange, p = k

  prob = construct_sample_Problem(k, u0, Z, alg)
  ensembleprob = EnsembleProblem(prob, output_func = output_func)

  sol = solve(ensembleprob, alg, ensemblealg, dt = get_dt(trange), trajectories=numtraj; kwargs...)
  return sol.u
end

# DefaultSampler for SDEKernel
function tangent!(du, u, dz, P::SDEKernel)
  @unpack f, g, p, noise_rate_prototype = P
  k1 = f(u[3],p,u[2])
  g1 = g(u[3],p,u[2])
  if noise_rate_prototype===nothing
    @. du[3] = k1*dz[2] + g1*dz[3]
  else
    du[3] .= k1*dz[2] + g1*dz[3]
  end
  (dz[1], dz[2], du[3])
end

function exponential_map!(u, du, P::SDEKernel)
  x = u[3]
  @. x += du[3]
  (u[1] + du[1], u[2] + du[2], x)
end

# g!(dx, x, p, t, dw) # dx .+= σ(t, x)*dW # in place

# du [index, time, space]

# DefaultSampler for SDEKernel
function tangent!(du, u, dz, P::SDEKernel!)
  @unpack f, gstep!, ws, p, noise_rate_prototype = P
  f(du[3], u[3], p, u[2])
  du[3] .*= dz[2]
  gstep!(du[3], ws, u[3], p, u[2], dz[3], noise_rate_prototype)
  (dz[1], dz[2], du[3])
end

function exponential_map!(u, du, P::SDEKernel!)
  x = u[3]
  @. x += du[3]
  (u[1] + du[1], u[2] + du[2], x)
end


text:
using ApproxFun, Compat, Base.Test

## Rectangle PDE

dx=dy=Interval()
d=dx*dy
x=Fun(identity,dx);y=Fun(identity,dy)

#dirichlet(d) is u[-1,:],u[1,:],u[:,-1],u[:,1]

G=[real(exp(-1+1.0im*y));
                        real(exp(1+1im*y));
                        real(exp(x-1im));
                        real(exp(x+1im));0.];

A=[dirichlet(d);lap(d)]
u=A\G
@test_approx_eq u(.1,.2) real(exp(0.1+0.2im))

A=[dirichlet(d);lap(d)+0.0I]
u=A\G
@test_approx_eq_eps u(.1,.2) real(exp(0.1+0.2im)) 1E-11

println("    Poisson tests")

## Poisson

f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2))  #default is [-1,1]^2
d=domain(f)
OS=S=schurfact([dirichlet(d);lap(d)],10)
u=OS\[zeros(∂(d));f]
@test_approx_eq u(.1,.2) -0.042393137972085826

d=PeriodicInterval()^2
f=ProductFun((x,y)->exp(-10(sin(x/2)^2+sin(y/2)^2)),d)
A=lap(d)+.1I
u=A\f
@test (lap(u)+.1u-f)|>coefficients|>norm < 1000000eps()

println("    Kron tests")

@static if is_apple()
    ## Kron

    dx=dy=Interval()
    d=dx*dy
    x=Fun(identity,dx);y=Fun(identity,dy)

    #dirichlet(d) is u[-1,:],u[1,:],u[:,-1],u[:,1]

    G=[real(exp(-1+1.0im*y));
                            real(exp(1+1im*y));
                            real(exp(x-1im));
                            real(exp(x+1im));0.];

    A=[dirichlet(d);lap(d)]

    S=schurfact(A,40)

    uex=A\G

    nx=ny=40
    K=kronfact(A,nx,ny)

    uex2=K\G

    @test (uex-uex2|>coefficients|>norm)<10000eps()

    # dirichlet bcs

    import ApproxFun.ChebyshevDirichlet

    S=ChebyshevDirichlet()⊗ChebyshevDirichlet();
    A=[dirichlet(S);lap(S)]
    nx=ny=20;
    KD=kronfact(A,nx,ny);

    #dirichlet(d) is u[-1,:],u[1,:],u[:,-1],u[:,1]
    x=Fun(identity);y=Fun(identity);
    G=[Fun(real(exp(-1+1.0im*y)),S[2]);
        Fun(real(exp(1+1im*y)),S[2]);
        Fun(real(exp(x-1im)),S[1]);
                            Fun(real(exp(x+1im)),S[1]);0.];

    uD=KD\G;

    @test_approx_eq uD(.1,.2) real(exp(.1+.2im))

    # fourth order
    dx=dy=Interval()
    d=dx*dy
    Dx=Derivative(dx);Dy=Derivative(dy)
    L=Dx^4⊗I+2*Dx^2⊗Dy^2+I⊗Dy^4

    K=kronfact([dirichlet(d);
         neumann(d);
         L],100,100)

    x=Fun(identity,dx);y=Fun(identity,dy)

    G=[real(exp(-1+1.0im*y));
                    real(exp(1+1im*y));
                    real(exp(x-1im));
                    real(exp(x+1im));
                    real(exp(-1+1.0im*y));
                    real(exp(1+1im*y));
                    -imag(exp(x-1im));
                    -imag(exp(x+1im))
       ]
    u=K\G
    @test_approx_eq u(.1,.2) real(exp(.1+.2im))

    # mixed

    K=kronfact([(ldirichlet(dx)+lneumann(dx))⊗I;
            (rdirichlet(dx)+rneumann(dx))⊗I;
            I⊗(ldirichlet(dy)+lneumann(dy));
            I⊗(rdirichlet(dy)+rneumann(dy));
            (ldirichlet(dx)-lneumann(dx))⊗I;
            (rdirichlet(dx)-rneumann(dx))⊗I;
            I⊗(ldirichlet(dy)-lneumann(dy));
            I⊗(rdirichlet(dy)-rneumann(dy));
             L],100,100)
    G=[2real(exp(-1+1.0im*y));
                    2real(exp(1+1im*y));
                    real(exp(x-1im))-imag(exp(x-1im));
                    real(exp(x+1im))-imag(exp(x+1im));
                    0;
                    0;
                    real(exp(x-1im))+imag(exp(x-1im));
                    real(exp(x+1im))+imag(exp(x+1im))
       ]
    u=K\G

    @test_approx_eq u(.1,.2) real(exp(.1+.2im))
end

## Test periodic x interval

d=PeriodicInterval()*Interval()
g=Fun(z->real(cos(z)),∂(d))  # boundary data
u=[dirichlet(d);lap(d)]\g

@test_approx_eq u(.1,.2) real(cos(.1+.2im))

dθ=PeriodicInterval(-2.,2.);dt=Interval(0,3.)
d=dθ*dt
Dθ=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
u=[I⊗ldirichlet(dt);Dt+Dθ]\Fun(θ->exp(-20θ^2),dθ)

@test_approx_eq u(.1,.2) exp(-20(0.1-0.2)^2)

println("   Domain Decompositon tests")

## Domain Decomposition
d=Interval(0,1)^2
n,m=20,80
A=discretize([dirichlet(d);lap(d)],n)
∂d=∂(d)
g=Fun(z->real(exp(z)),∂d)
f=[Fun([zeros(k-1);1.0],∂d) for k=1:m].'
U=A\f
@test_approx_eq dot(real(g.coefficients),U[1:ncoefficients(g)])(.1,.2) real(exp(.1+.2im))

Rectangle(a,b,c,d)=Interval(a,b)*Interval(c,d)
Γ=Rectangle(0,1,0,1)∪Rectangle(1,2,0,1)
Fun(identity,∂(Γ))|>values

## Small diffusoion

dx=Interval();dt=Interval(0,1.)
d=dx*dt
Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
x=Fun(identity,dx)
B=0.0
C=0.0
V=B+C*x
ε=0.001
f=Fun(x->exp(-20x^2),dx)
u=[timedirichlet(d);Dt-ε*Dx^2-V*Dx]\f
@test_approx_eq u(.1,.2) 0.8148207991358946
B=0.1
C=0.2
V=B+C*x
u=[timedirichlet(d);Dt-ε*Dx^2-V*Dx]\f
@test_approx_eq u(.1,.2) 0.7311625132209619

## Schrodinger

dx=Interval(0.,1.);dt=Interval(0.0,.1)
d=dx*dt

V=Fun(x->x^2,dx)

Dt=Derivative(d,[0,1]);Dx=Derivative(d,[1,0])

ϵ=1.
u0=Fun(x->exp(-100*(x-.5)^2)*exp(-1./(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
L=ϵ*Dt+(.5im*ϵ^2*Dx^2)
ny=200;u=pdesolve([timedirichlet(d);L],u0,ny)
@test_approx_eq_eps u(.2,.1) (0.2937741918470843 + 0.22130344715160255im )  0.000001

## Periodic

d=PeriodicInterval()^2
f=Fun((θ,ϕ)->exp(-10(sin(θ/2)^2+sin(ϕ/2)^2)),d)
A=lap(d)+.1I
u=A\f
@test_approx_eq u(.1,.2) u(.2,.1)

d=PeriodicInterval()*Interval()
g=Fun(z->real(cos(z)),∂(d))  # boundary data
u=[dirichlet(d);lap(d)]\g
@test_approx_eq u(.1,.2) real(cos(.1+.2im))

dθ=PeriodicInterval(-2.,2.);dt=Interval(0,3.)
d=dθ*dt
Dθ=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
u=[I⊗ldirichlet(dt);Dt+Dθ]\Fun(θ->exp(-20θ^2),dθ)

# Check bug in cache
CO=cache(ldirichlet(dt))
ApproxFun.resizedata!(CO,:,2)
ApproxFun.resizedata!(CO,:,4)
@test_approx_eq CO*Fun(exp,dt) 1.0

dθ=PeriodicInterval(-2.,2.);dt=Interval(0,3.)
d=dt*dθ
Dt=Derivative(d,[1,0]);Dθ=Derivative(d,[0,1])
A=[ldirichlet(dt)⊗I;Dt+Dθ]
f=Fun(θ->exp(-20θ^2),dθ)
ut=A\f

@test_approx_eq u(.1,.2) ut(.2,.1)

# Beam

dθ=PeriodicInterval(0.0,1.0);dt=Interval(0,0.03)
d=dθ*dt
Dθ=Derivative(d,[1,0]);Dt=Derivative(d,[0,1]);

B=[I⊗ldirichlet(dt),I⊗lneumann(dt)]
u=pdesolve([B;Dt^2+Dθ^4],Fun(θ->exp(-200(θ-.5).^2),dθ),200)

@test_approx_eq_eps u(.1,.01) -0.2479768394633227  1E-8 #empirical

## Rectangle PDEs

# Screened Poisson

d=Interval()^2
u=[neumann(d);lap(d)-100.0I]\ones(∂(d))
@test_approx_eq u(.1,.9) 0.03679861429138079

# PiecewisePDE

a=Fun([1,0.5,1],[-1.,0.,0.5,1.])
s=space(a)
dt=Interval(0,2.)
Dx=Derivative(s);Dt=Derivative(dt)
Bx=[ldirichlet(s);continuity(s,0)]

# test resize bug
CO=cache(Bx[2])
ApproxFun.resizedata!(CO,:,2)
ApproxFun.resizedata!(CO,:,4)
@test_approx_eq CO.data*collect(1:4) [3.,-1.]

u=pdesolve([I⊗ldirichlet(dt);Bx⊗I;I⊗Dt+(a*Dx)⊗I],Any[Fun(x->exp(-20(x+0.5)^2),s)],200)
@test_approx_eq_eps u(-.1,.2) exp(-20(-.2-.1+0.5)^2) 0.00001

## Test error

dx=Interval();dt=Interval(0,2.)
d=dx*dt
Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
x=Fun(identity,dx)
u=[I⊗ldirichlet(dt);Dt+x*Dx]\Fun(x->exp(-20x^2),dx)

@test_approx_eq u(0.1,0.2) 0.8745340845783758  # empirical

dθ=PeriodicInterval();dt=Interval(0,10.)
d=dθ*dt
ε=.01
Dθ=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])

# Parentheses are a hack to get rank 2 PDE
u=[I⊗ldirichlet(dt);Dt-ε*Dθ^2-Dθ]\Fun(θ->exp(-20θ^2),dθ)

@test_approx_eq_eps u(0.1,0.2) 0.1967278179230314 1000eps()


text:
using Printf
using GLMakie
using Oceananigans
using Oceananigans.Units

include("cyclic_interpolate_utils.jl")

# 2.8125 degree resolution
Nx = 128
Ny = 60
Nz = 18
reference_density = 1035

#####
##### Load forcing files roughly from CORE2 paper
##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html)
#####

east_west_stress_path = "off_TAUXvar1.bin"
north_south_stress_path = "off_TAUY.bin"
sea_surface_temperature_path="sst25_128x60x12.bin"

thirty_days = 30days
Nmonths = 12
bytes = sizeof(Float32) * Nx * Ny

τˣ = - reshape(bswap.(reinterpret(Float32, read(east_west_stress_path, Nmonths * bytes))), (Nx, Ny, Nmonths)) ./ reference_density
τʸ = - reshape(bswap.(reinterpret(Float32, read(north_south_stress_path, Nmonths * bytes))), (Nx, Ny, Nmonths)) ./ reference_density
T★ = reshape(bswap.(reinterpret(Float32, read(sea_surface_temperature_path, Nmonths * bytes))), (Nx, Ny, Nmonths))

times = 0:1days:24*30days
discrete_times = 0:30days:11*30days

#=
i = 80
j = 12 # Southern ocean?
Nt = length(times)
τ_ij = zeros(Nt)

for (m, t) in enumerate(times)
    n₁ = current_time_index(t)
    n₂ = next_time_index(t)
    τ_ij[m] = cyclic_interpolate(τˣ[i, j, n₁], τˣ[i, j, n₂], t)
end

fig = Figure(resolution = (1200, 600))

ax = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Target SST at (i, j) = ($i, $j)")
scatter!(ax, discrete_times ./ day, τˣ[i, j, :], markersize=20, marker=:utriangle, color=:red)
lines!(ax, times ./ day, τ_ij)

display(fig)
=#

#####
##### Analyze forcing
#####

fig = Figure(resolution = (1200, 800))

t = Node(0.0)

T★_continuous(t) = cyclic_interpolate(T★, t)
τˣ_continuous(t) = cyclic_interpolate(τˣ, t)
τʸ_continuous(t) = cyclic_interpolate(τʸ, t)

T★_t = @lift T★_continuous($t)
τˣ_t = @lift τˣ_continuous($t)
τʸ_t = @lift τʸ_continuous($t)

max_τ = maximum(abs, τˣ)
max_T = maximum(T★)
min_T = minimum(T★)

fig = Figure(resolution = (1200, 600))

ax_τˣ = Axis(fig[1, 1], title="East-west wind stress (m² s⁻²)")
hm_τˣ = heatmap!(ax_τˣ, τˣ_t, colorrange=(-max_τ, max_τ), colormap=:balance)
cb_τˣ = Colorbar(fig[1, 2], hm_τˣ)

ax_τʸ = Axis(fig[2, 1], title="North-south wind stress (m² s⁻²)")
hm_τʸ = heatmap!(ax_τʸ, τʸ_t, colorrange=(-max_τ, max_τ), colormap=:balance)
cb_τʸ = Colorbar(fig[2, 2], hm_τʸ)

ax_T = Axis(fig[3, 1], title="Sea surface temperature (ᵒC)")
hm_T = heatmap!(ax_T, T★_t, colorrange=(min_T, max_T), colormap=:thermal)
cb_T = Colorbar(fig[3, 2], hm_T)

title_str = @lift "Earth day = " * prettytime($t)
ax_t = fig[0, :] = Label(fig, title_str)

display(fig)

record(fig, "annual_cycle_forcing.mp4", times, framerate=8) do tn
    t[] = tn
end


text:
using Benchmark

tests = ["benchmarks.jl", "compare.jl", "consistency.jl"]

@printf "Running tests:\n"

for t in tests
    include(t)
    @printf " * %s\n" t
end


text:
function options_ASB07_decomposed()
    𝑂spec = Vector{OptionSpec}()

    push!(𝑂spec, OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
        domain_check=(v  ->  v > 0.), info="time step"))

    push!(𝑂spec, OptionSpec(:max_order, 10, domain=Int,
        info="maximum allowed order of zonotopes"))

    push!(𝑂spec, OptionSpec(:order_discretization, 2, domain=Int,
        info="order of Taylor approximation in discretization"))

    push!(𝑂spec, OptionSpec(:set_operations_discretization, "zonotope",
        domain=String, domain_check=(v  ->  v ∈ ["zonotope", "lazy"]),
        info="type of set operations applied during discretization"))

    push!(𝑂spec, OptionSpec(:block_options_init, Zonotope, domain=Any,
        info="option for the decomposition of the initial states"))

    push!(𝑂spec, OptionSpec(:partition, [Int[]],
        domain=AbstractVector{<:AbstractVector{Int}}, domain_check=ispartition,
        info="block partition; a block is represented by a vector containing " *
             "its indices"))

    push!(𝑂spec, OptionSpec(:vars, Int[], domain=AbstractVector{Int},
        domain_check=(v  ->  length(v) > 0 && all(e -> e > 0, v)),
        info="variables of interest; default: all variables"))

    return 𝑂spec
end


text:
module UInt12ArraysBase

export UInt12Array, UInt12Vector, UInt12Matrix

import ..UInt12Arrays: default_eltype
import ..UInt12s: UInt12
import ..UInt24s: UInt24

abstract type AbstractUInt12Array{T, N} <: AbstractArray{T,N} end

"""
    UInt12Array{T, B, N}(data::B [, size::NTuple{N,Int} ])
    UInt12Array{T, B, N}(undef    , size::NTuple{N,Int}  )

UInt12Array represents an array of densely packed 12-bit integers.

UInt12AArray has the following parameters:
`T` - Element type of the UInt12Array. Typically this might be either UInt12 or
    UInt16. Default is UInt12Arrays.default_eltype.
`B` - Base type for the UIn12Array's internal vector storage.
`N` - Number of dimensions of the array.
"""
mutable struct UInt12Array{T, B <: AbstractVector, N} <: AbstractUInt12Array{T,N}
    data::B
    size::NTuple{N,Int}
end

const UInt12Vector{T,B} = UInt12Array{T, B, 1}
const UInt12Matrix{T,B} = UInt12Array{T, B, 2}

# Constructors

# Copy constructor, still a view on the underlying data
UInt12Array(A::UInt12Array{T, B, N}) where {T,B,N} =
    UInt12Array{T, B, N}(A.data, A.size)

"""
    UInt12Vector(data::Vector{UInt8})

    The default constructor has an element type of UInt16 and is backed by a
    byte (UInt8) vector that consists of packed 12-bit integers. Every three
    bytes in the vector represents two 12-bit integers.

    The following packing is assumed:
    The 8 lower bits of the first integer are contained in the first byte.
    The 4 higher bits of the first integer are first 4 bits of the second byte.
    The 4 lower bits of the second integer are the last 4 bits of the second byte.
    The 8 higher bits of the second integer are contained in the third byte.

# Example
```jldoctest
julia> A = UInt12Array(UInt8[0xba, 0xdc, 0xfe])
2-element UInt12Array{UInt16,Array{UInt8,1},1}:
 0x0cba
 0x0fed

julia> A[1]
0x0cba

julia> copy(A)
2-element Array{UInt16,1}:
 0x0cba
 0x0fed

julia>
```
"""
UInt12Vector(data::B) where B = UInt12Vector{default_eltype}(data)
UInt12Vector(data::B, length::Int) where B = UInt12Vector{default_eltype}(data, length)
UInt12Matrix(data::B, row::Int, col::Int) where B =
    UInt12Matrix{default_eltype}(data, row, col)
UInt12Array(data::B, row::Int, col::Int, depth::Int) where B =
    UInt12Array{default_eltype}(data, row, col, depth)

UInt12Array(data::B) where B <: AbstractVector = UInt12Vector(data)

"""
    UInt12Vector{T}(data::Vector{UInt8})

    The single parameter constructor allows for specification of an element type. For example
    you may want the element type to actually be a UInt12 rather than the default UInt16.
"""
UInt12Vector{T}(data::B) where {T,B <: AbstractVector} =
    UInt12Vector{T}( reinterpret(UInt8, data) )
UInt12Vector{T}(data::B) where {T,B <: AbstractVector{UInt8}} =
    UInt12Vector{T,B}(data, size(data) .*2 .÷ 3)
UInt12Vector{T}(data::B) where {T,B <: AbstractArray{UInt8}} =
    UInt12Vector{T,B}(data, size(data) .*2 .÷ 3)
UInt12Vector{T}(data::B, length::Int) where {T,B} =
    UInt12Vector{T,B}(data, (length,))
UInt12Matrix{T}(data::B, row::Int, col::Int) where {T,B} =
    UInt12Matrix{T,B}(data, (row,col))
UInt12Array{T}(data::B, row::Int, col::Int, depth::Int) where {T,B} =
    UInt12Array{T,B,3}(data, (row, col, depth))

# Undefined initialization

## Element type not defined, default to UInt16
UInt12Vector(::UndefInitializer, length::Int) = UInt12Vector{default_eltype}(undef, length)
UInt12Matrix(::UndefInitializer, row::Int, col::Int) = UInt12Matrix{default_eltype}(undef, row, col)

UInt12Array(::UndefInitializer, row::Int, col::Int, depth::Int) =
    UInt12Array{UInt16}(undef, row, col, depth)

UInt12Array(::UndefInitializer, d::Vararg{Int, N}) where N =
    UInt12Array{default_eltype}(undef, d...)

## Element type defined
UInt12Vector{T}(::UndefInitializer, length::Int) where T = UInt12Array{T}(undef, length)
UInt12Matrix{T}(::UndefInitializer, row::Int, col::Int) where T = UInt12Array{T}(undef, row, col)

UInt12Array{T}(::UndefInitializer, length::Int) where T =
    UInt12Array{T,Vector{UInt8},1}( Vector{UInt8}(undef, ceil(Int, length*3/2)), (length,) )

UInt12Array{T}(::UndefInitializer, rows::Int, cols::Int) where T =
    UInt12Array{T,Vector{UInt8},2}( Vector{UInt8}(undef, ceil(Int, rows*cols*3/2)), (rows, cols,) )

UInt12Array{T}(::UndefInitializer, rows::Int, cols::Int, depth::Int) where T =
    UInt12Array{T,Vector{UInt8},3}( Vector{UInt8}(undef, ceil(Int, rows * cols * depth * 3/2)), (rows, cols, depth) )

UInt12Array{T}(::UndefInitializer, d::Vararg{Int, N}) where {T,N} =
    UInt12Array{T,Vector{UInt8},N}( Vector{UInt8}(undef, ceil(Int, prod(d)*3/2) ), d)

function UInt12Array{T, B, N}(::UndefInitializer, size::NTuple{N,Int}) where {T,B <: AbstractVector,N}
    E = eltype(B)
    szE = sizeof(E)
    UInt12Array{T, B, N}(B(undef, div.( size .* 3 .÷ szE, 2, RoundUp) ), size)
end

# Generic methods for all parameters
Base.size(data::UInt12Array) = data.size
Base.IndexStyle(::Type{<: UInt12Array}) = IndexLinear()

# Methods for B <: AbstractVector{UInt8}, this assumes a straightforward, dense packing
function Base.getindex(A::UInt12Array{T,<: AbstractVector{UInt8},N}, i::Int) where {T,N}
    @boundscheck checkbounds(A, i)
    first_byte_index = ( (i - 1) ÷ 2 ) *3 + 1
    ord = mod1(i,2)
    if ord == 1
        out = A.data[ first_byte_index ] | ( UInt16(A.data[ first_byte_index + 1] & 0x0f) << 8 )
    else
        out = (A.data[ first_byte_index + 1 ] & 0xf0) >> 4 | ( UInt16(A.data[ first_byte_index + 2]) << 4 )
    end
    convert(T, out)
end
function Base.setindex!(A::UInt12Array{T,<: AbstractVector{UInt8},N}, v, i::Int) where {T,N}
    @boundscheck checkbounds(A, i)
    if v > 0xfff
        throw( InexactError(:UInt12, UInt12, v))
    end
    first_byte_index = ( (i - 1) ÷ 2 ) * 3 + 1
    ord = mod1(i,2)
    p = A.data[ first_byte_index + 1 ]
    if ord == 1
        A.data[ first_byte_index ] = 0xff & v
        A.data[ first_byte_index + 1 ] = ( ( 0xf00 & v ) >> 8 ) | (p & 0xf0)
    else
        A.data[ first_byte_index + 1] = ( ( 0x0f & v) << 4 ) | (p & 0x00f)
        A.data[ first_byte_index + 2] = ( 0xff0 & v ) >> 4
    end
end
function Base.resize!(A::UInt12Array{T,B,N}, nl::Integer) where {T,B <: AbstractVector{UInt8},N}
    resize!(A.data, ceil(Int, nl*3/2) )
    A.size = (nl,)
end

# Methods for B <: AbstractVector{UInt24}
function Base.getindex(A::UInt12Array{T,<: AbstractVector{UInt24},N}, i::Int) where {T,N}
    #@boundscheck checkbounds(A, i)
    uint24_index = ( (i - 1) ÷ 2 + 1)
    ord = mod1(i,2)
    if ord == 1
        out = A.data[ uint24_index ] & 0xfff
    else
        out = A.data[ uint24_index ] >>> 12
    end
    convert(T, out)
end
function Base.setindex!(A::UInt12Array{T,<: AbstractVector{UInt24},N}, v, i::Int) where {T,N}
    @boundscheck checkbounds(A, i)
    if v > 0xfff
        throw( InexactError(:UInt12, UInt12, v))
    end
    uint24_index = ( (i - 1) ÷ 2 + 1)
    ord = mod1(i,2)
    if ord == 1
        A.data[ uint24_index ] = (0xfff & v) | (0xfff000 & A.data[ uint24_index ])
    else
        A.data[ uint24_index ] = (0x000fff & v) << 12 | (0x000fff & A.data[ uint24_index ])
    end
end
function Base.resize!(A::UInt12Array{T,B,N}, nl::Integer) where {T,B <: AbstractVector{UInt24},N}
    resize!(A.data, div(nl, 2, RoundUp) )
    A.size = (nl,)
end
function Base.parent(A::UInt12Array{T,B,N}) where {T,B,N}
    A.data
end
"""
    map_idex_to_byte(idx, byte_idx = 1)

    Map array index to an underlying byte. byte_idx ∈ 1:3
"""
map_idx_to_byte(idx, byte_idx = 1) = (idx - 1) ÷ 2 * 3 + byte_idx

function Base.parentindices(A::UInt12Array{T,B,N}, indices::Base.OneTo = (eachindex(A),)) where {T, B <: AbstractVector{UInt8}, N}
    # Only calculate the last index since we know the first index must be `1`
    last_indices = last.(indices)
    # If odd, map the 2nd byte. If even, map to the 3rd byte.
    last_byte_idx = map_idx_to_byte.(last_indices, iseven.(last_indices) .+ 2)
    Base.OneTo.(last_byte_idx)
end
function Base.parentindices(A::UInt12Array{T,B,N}, indices) where {T, B <: AbstractVector{UInt8}, N}
    first_byte_idx = map_idx_to_byte.(first.(indices))
    last_indices = last.(indices)
    # If odd, map the 2nd byte. If even, map to the 3rd byte.
    last_byte_idx = map_idx_to_byte.(last_indices, iseven.(last_indices) .+ 2)
    ntuple(i->first_byte_idx[i]:last_byte_idx[i], length(indices))
end
function Base.convert(::Type{UInt12Array{T}}, S::SubArray{T,N,UInt12Array{T,B,N},I,L}) where {T,B,N,I,L}
    p_idx = parentindices(S)
    @assert all(isodd.(first.(p_idx))) "The initial parent indices of each dimension must be odd."
    idx = parentindices(parent(S), p_idx)
    UInt12Vector{T}(view(parent(parent(S)), idx...))
end
Base.convert(::Type{UInt12Array}, S::SubArray{T,N,UInt12Array{T,B,N},I,L}) where {T,B,N,I,L} = convert(UInt12Array{T}, S)

end # module UInt12ArraysBase


text:
mutable struct Dish
    count::Int64
    cluster_params::cluster_parameters
end

mutable struct Table
    count::Int64
    dish_index::Int64
end

mutable struct Restaurant
    tables::Vector{Table}
    points::AbstractArray{Float64,2}
    labels::AbstractArray{Int64,2}
end

mutable struct Crf_model
    dishes::Vector{Dish}
    restaurants::Vector{Restaurant}
    prior::distribution_hyper_params
end

function create_new_dish(prior)
    #pts = point[:,:]
    suff = create_sufficient_statistics(prior,[])
    post = calc_posterior(prior,suff)
    dist = sample_distribution(prior)
    cluster_params = cluster_parameters(prior,dist,suff,post)
    dish = Dish(1,cluster_params)
    return dish
end

function resample_dish!(pts_vec,dish)
    if length(pts_vec) == 0
        suff = create_sufficient_statistics(dish.cluster_params.hyperparams,[])
        post = calc_posterior(dish.cluster_params.hyperparams,suff)
        dist = sample_distribution(dish.cluster_params.hyperparams)
        cluster_params = cluster_parameters(dish.cluster_params.hyperparams,dist,suff,post)
        dish.cluster_params = cluster_params
    else
        all_pts = reduce(hcat,pts_vec)
        suff = create_sufficient_statistics(dish.cluster_params.hyperparams,dish.cluster_params.hyperparams,all_pts)
        post = calc_posterior(dish.cluster_params.hyperparams,suff)
        dist = sample_distribution(post)
        cluster_params = cluster_parameters(dish.cluster_params.hyperparams,dist,suff,post)
        dish.cluster_params = cluster_params
    end
end

function init_crf_model(data,prior)
    restaurants = []
    dishes = []
    for i=1:length(data)
        v = data[i]
        table = Table(size(v,2),1)
        restaurant = Restaurant([table],v,ones(size(v,2),1))
        push!(restaurants,restaurant)
    end
    first_dish = create_new_dish(prior)
    first_dish.count = length(data) #Amount of tables at init
    model = Crf_model([first_dish],restaurants,prior)
    return model
end

function get_tables_probability_vec(tables,α)
    counts = [x.count for x in tables]
    push!(counts,α)
    weights = [x/sum(counts) for x in counts]
    return weights
end

function sample_dish(dishes, pts, γ,new_dish,cur_dish = -1)
    counts = [x.count for x in dishes]
    if cur_dish > 0
        counts[cur_dish] -= 1
    end
    push!(counts,γ)
    weights = [x/sum(counts) for x in counts]
    dists = [x.cluster_params.distribution for x in dishes]
    push!(dists, new_dish.cluster_params.distribution)
    parr = zeros(length(size(pts,2)), length(dists))
    for (k,v) in enumerate(dists)
        log_likelihood!(reshape((@view parr[:,k]),:,1), pts,v)
    end
    parr = sum(parr, dims = 1)[:,:] + reshape(log.(weights),1,:)
    labels = [0][:,:]
    sample_log_cat_array!(labels,parr)
    return labels[1,1]
end

function sample_table(dists,weights,pts)
    parr = zeros(length(size(pts,2)), length(dists))
    for (k,v) in enumerate(dists)
        log_likelihood!(reshape((@view parr[:,k]),:,1), pts,v)
    end
    parr = sum(parr, dims = 1) + reshape(log.(weights),1,:)

    labels = [0][:,:]
    sample_log_cat_array!(labels,parr)
    return labels[1,1]
end

function crf_hdp_fit(data,α,γ,prior,iters)
    model = init_crf_model(data,prior)
    #Each iteration is a single pass on all points
    for iter=1:iters
        println("Iter " * string(iter) * "  Dish count: " * string(length(model.dishes)))
        #Resample dishes distribution
        cur_new_dish = create_new_dish(prior)
        for (dish_index,dish) in enumerate(model.dishes)
            pts = []
            for restaurant in model.restaurants
                for (table_index,table) in enumerate(restaurant.tables)
                    if table.dish_index == dish_index
                        push!(pts, restaurant.points[:,restaurant.labels[:] .== table_index])
                    end
                end
            end
            resample_dish!(pts,dish)
        end
        #Samples points
        for restaurant in model.restaurants
            # println([x.count for x in restaurant.tables])
            for i=1:size(restaurant.points,2)
                point = restaurant.points[:,i][:,:]
                prev_table = restaurant.labels[i,1]
                restaurant.tables[prev_table].count -= 1
                weights = get_tables_probability_vec(restaurant.tables,α)
                new_table_dish = sample_dish(model.dishes, point,γ,cur_new_dish)
                dists = [model.dishes[x.dish_index].cluster_params.distribution for x in restaurant.tables]
                if new_table_dish > length(model.dishes)
                    push!(dists, cur_new_dish.cluster_params.distribution)
                else
                    push!(dists,model.dishes[new_table_dish].cluster_params.distribution)
                end
                try
                    log.(weights)
                catch e
                    println("error print:")
                    println([x.count for x in restaurant.tables])
                    rethrow([e])
                end
                new_table = sample_table(dists,weights,point)
                restaurant.labels[i,1] = new_table
                if new_table > length(restaurant.tables)
                    push!(restaurant.tables,Table(1,new_table_dish))
                    if new_table_dish > length(model.dishes)
                        cur_new_dish.count = 1
                        push!(model.dishes, cur_new_dish)
                        cur_new_dish = create_new_dish(prior)
                    else
                        model.dishes[new_table_dish].count +=1
                    end
                else
                    restaurant.tables[new_table].count += 1
                end
                if restaurant.tables[prev_table].count == 0
                    model.dishes[restaurant.tables[prev_table].dish_index].count -= 1
                end
            end
        end
        #Samples tables
        for restaurant in model.restaurants
            for (index,table) in enumerate(restaurant.tables)
                if table.count == 0
                    continue
                end
                point = restaurant.points[:,(restaurant.labels .== index)[:]]
                prev_dish = table.dish_index
                new_table_dish = sample_dish(model.dishes, point,γ,cur_new_dish,prev_dish)
                model.dishes[prev_dish].count -=1
                if new_table_dish > length(model.dishes)
                    cur_new_dish.count = 1
                    push!(model.dishes, cur_new_dish)
                    cur_new_dish = create_new_dish(prior)
                else
                    model.dishes[new_table_dish].count +=1
                end
                # println("old dish:" * string(prev_dish) * " new dish:" * string(new_table_dish))
                table.dish_index = new_table_dish
            end
        end
    end
    return model
end

function get_dish_labels_from_model(model::Crf_model)
    labels_dict = Dict()
    for (k,v) in enumerate(model.restaurants)
        labels = [v.tables[v.labels[i,1]].dish_index for i=1:size(v.points,2)]
        labels_dict[k] = labels
    end
    return labels_dict
end


text:
mutable struct SolverOptions
    backtrack           :: Bool
    max_bt_steps        :: Int64
    backtrack_sorting   :: Bool
    keep_logs           :: Bool
end

function get_default_options()
    backtrack           = true
    max_bt_steps        = typemax(Int64)
    backtrack_sorting   = true
    keep_logs           = false

    return SolverOptions(backtrack, max_bt_steps, backtrack_sorting, keep_logs)
end

function combine_options(options)
    defaults = get_default_options()
    options_dict = Dict{Symbol,Any}()
    for kv in options
        if !in(kv[1], fieldnames(SolverOptions))
            @warn "Option "*string(kv[1])*" is not available"
        else
            options_dict[kv[1]] = kv[2]
        end
    end

    for fname in fieldnames(SolverOptions)
        if haskey(options_dict, fname)
            setfield!(defaults, fname, options_dict[fname])
        end
    end
    return defaults
end


text:
function add_rules(from, to, n)
    for α in [0, 1, 2, 3]
        from_rot = rotr90(from, α)
        rules[n][from_rot] = to
        rules[n][flipdim(from_rot, 1)] = to
        rules[n][flipdim(from_rot, 2)] = to
    end
end
parse_array(v, n) = reshape(BitArray([x == '.' ? false : true for x in v]), (n, n))'

d = split.(replace.(readlines("input.txt"), "/", ""), " => ")
rules = Dict(2 => Dict(), 3 => Dict())
for l in d
    n = Int(sqrt(length(l[1])))
    from = parse_array(l[1], n)
    to = parse_array(l[2], n + 1)
    add_rules(from, to, n)
end

function iter(grid)
    s = size(grid)[1]
    t = s % 2 == 0 ? 2 : 3
    squares = Int(s / t)
    new_grid = falses(squares * (t+1), squares * (t+1))
    for x in 1:squares
        for y in 1:squares
            update = rules[t][grid[t*(x-1)+1:t*x, t*(y-1)+1:t*y]]
            new_grid[(t+1)*(x-1)+1:(t+1)*x, (t+1)*(y-1)+1:(t+1)*y] = update
        end
    end
    return new_grid
end

grid = parse_array(".#...####", 3)
for i in 1:18
    grid = iter(grid)
    if i in [5, 18]
        println(sum(grid))
    end
end


text:
## =========================== ## ## =========================== ##
## 	   Package of collective motion models			 ##
##	   Martin Zumaya Hernandez 						             ##
##	   Release Version 0.1.0 						             ##
##     29 / 05 / 2018									         ##
## =========================== ## ## =========================== ##

module CollectiveDynamics

## SIMPLE VICSEK MODEL 2D
module SVM2D
using LinearAlgebra
export Box, Flock, set_output_data_structure_vsk, evolve_system
include("simple_vicsek_model_2D.jl")
end

## SIMPLE VICSEK MODEL 3D
# module SVM3D
# using Quaternions
# export Box, Flock, set_output_data_structure_vsk, evolve_system
# include("simple_vicsek_model_3D.jl")
# end

## BEHAVIOURAL RULES MODEL
# module BehaviouralRules
# using Quaternions
# export set_output_data_structure, evolve_system
# include("behavioural_rules_model.jl")
# end

## SHORT AND LONG RANGE INTERACTIONS MODEL
# module ShortLongRange
# using Distributions, Quaternions
# export set_output_data_structure, evolve_metric_system, evolve_metric_system_2D, evolve_topological_system, evolve_topological_system_2D
# include("short_long_range_model.jl")
# end

## INTERIAL SPIN MODEL MODEL
# module InertialSpin
# using Quaternions
# export InertialFlock, InertialExtFlock, SharedInertialExtFlock, InertialParameters, set_output_data_structure, set_output_data_structure_lr, evolve_system, evolve_extended_system, evolve_extended_system_sh
# include("inertial_spin_model.jl")
# end

## NETWORK BASED MODEL
# module NetworkModel
# using Distributions
# export set_output_data_structure, set_directed_random_network, evolve_system_2D
# include("network_model.jl")
# end

end


text:
"""
    interpolation_search(array::Array{Int64,1}, target::Int64)

The Interpolation search algorithm (ISA) finds a target-position in a sorted array by using
a numerical procedure. The sorting procedure uses a linear fitting for finding the target
position of the remaining search space in the array in more detail. The array's target
position is calculated by the straight slope between the lowest and largest boundary of the
remaining array and the lowest array position itself during each optimization cycle. If the
target-position cannot be found, the array-space will be shrink for the lower or higher
boundary region based on a comparison. For more information see: [https://en.wikipedia.org/wiki/Interpolation_search](https://en.wikipedia.org/wiki/Interpolation_search)

# Arguments
- `array::Array{Int64,1}`: Sorted array of integers
- `target::Int64`: Target-value to find the position

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [10,11, 12, 14, 16, 18, 19, 20, 21, 22, 23, 24, 33, 35, 42, 47]
julia> target = 12
julia> ClassicAlgorithmsCollections.interpolation_search(arr, target)
3
```
"""
function interpolation_search(array::Array{Int64,1}, target::Int64)
    # Find indexs of two low and high array-boundary
    low = 1
    high = length(array)

    # Since array is sorted, an element present in array must be in range defined by boundary
    while low <= high && target >= array[low] && target <= array[high]

        # Check if array contains only one element
        if low == high
            if array[low] == target
                return low
            end
            return nothing
        end

        # Linear fitting of the new position pos and transforming to an integer value.
        pos =
            low + floor(
                Int64,
                ((float(high - low) / (array[high] - array[low])) * (target - array[low])),
            )

        #  Check if target is now found
        if array[pos] == target
            return pos
        end

        # If x is larger, x is in upper part
        if array[pos] < target
            low = pos + 1
            # Else x has to be smaller, x is in lower part
        else
            high = pos - 1
        end
    end
    # If target is not in the list
    return nothing
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# machine code generation

function machine(cap::VersionNumber, triple::String)
    InitializeNVPTXTarget()
    InitializeNVPTXTargetInfo()
    t = Target(triple)

    InitializeNVPTXTargetMC()
    cpu = "sm_$(cap.major)$(cap.minor)"
    if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
        # in the case of CUDA 9, we use sync intrinsics from PTX ISA 6.0+
        feat = "+ptx60"
    else
        feat = ""
    end
    tm = TargetMachine(t, triple, cpu, feat)
    asm_verbosity!(tm, true)

    return tm
end

# final preparations for the module to be compiled to PTX
# these passes should not be run when e.g. compiling to write to disk.
function prepare_execution!(job::CompilerJob, mod::LLVM.Module)
    let pm = ModulePassManager()
        global current_job
        current_job = job

        global_optimizer!(pm)

        add!(pm, ModulePass("ResolveCPUReferences", resolve_cpu_references!))

        global_dce!(pm)
        strip_dead_prototypes!(pm)

        run!(pm, mod)
        dispose(pm)
    end

    return
end

# some Julia code contains references to objects in the CPU run-time,
# without actually using the contents or functionality of those objects.
#
# prime example are type tags, which reference the address of the allocated type.
# since those references are ephemeral, we can't eagerly resolve and emit them in the IR,
# but at the same time the GPU can't resolve them at run-time.
#
# this pass performs that resolution at link time.
function resolve_cpu_references!(mod::LLVM.Module)
    job = current_job::CompilerJob
    changed = false

    for f in functions(mod)
        fn = LLVM.name(f)
        if isdeclaration(f) && intrinsic_id(f) == 0 && startswith(fn, "jl_")
            # eagerly resolve the address of the binding
            address = ccall(:jl_cglobal, Any, (Any, Any), fn, UInt)
            dereferenced = unsafe_load(address)
            dereferenced = LLVM.ConstantInt(dereferenced, JuliaContext())

            function replace_bindings!(value)
                changed = false
                for use in uses(value)
                    val = user(use)
                    if isa(val, LLVM.ConstantExpr)
                        # recurse
                        changed |= replace_bindings!(val)
                    elseif isa(val, LLVM.LoadInst)
                        # resolve
                        replace_uses!(val, dereferenced)
                        unsafe_delete!(LLVM.parent(val), val)
                        # FIXME: iterator invalidation?
                        changed = true
                    end
                end
                changed
            end

            changed |= replace_bindings!(f)
        end
    end

    return changed
end

function mcgen(job::CompilerJob, mod::LLVM.Module, f::LLVM.Function)
    tm = machine(job.cap, triple(mod))

    InitializeNVPTXAsmPrinter()
    return String(emit(tm, mod, LLVM.API.LLVMAssemblyFile))
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -32.35079088597315 <= q <= 33.04707617725232)

                     add_NL_constraint(m, :(swish(0.7516814527569986 + 0.07211062597860263*swish(0.800956367833014 + 0.1575422515944389*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.28099609060437203*swish(-0.8241959750660866 + 0.48752067922981546*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.19401527858092393*swish(-0.014178711859067938 + -0.6496362145686132*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.43991506890783594*swish(-0.38559964579408224 + -0.6882459841770032*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + swish(0.5442086655639207 + -0.859033717055464*swish(0.800956367833014 + 0.1575422515944389*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.6994107397450926*swish(-0.8241959750660866 + 0.48752067922981546*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.21447765318835765*swish(-0.014178711859067938 + -0.6496362145686132*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9247182813449877*swish(-0.38559964579408224 + -0.6882459841770032*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + swish(-0.33457802196328057 + 0.8511359109817849*swish(0.800956367833014 + 0.1575422515944389*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9056512987027685*swish(-0.8241959750660866 + 0.48752067922981546*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.901983268277478*swish(-0.014178711859067938 + -0.6496362145686132*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.02854306019224717*swish(-0.38559964579408224 + -0.6882459841770032*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + swish(0.85223750045835 + -0.7468516096711739*swish(0.800956367833014 + 0.1575422515944389*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.1269236474608202*swish(-0.8241959750660866 + 0.48752067922981546*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.5407492362268407*swish(-0.014178711859067938 + -0.6496362145686132*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.8670343803826994*swish(-0.38559964579408224 + -0.6882459841770032*swish(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*swish(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*swish(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*swish(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) - $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:
module Fixtures

include("fixture.jl")
include("registered_fixtures.jl")
include("fact_check_support.jl")
include("matchers.jl")
include("file_fixtures.jl")
include("mock/mock.jl")
include("mock/patch.jl")

end


text:
_generic_lp_key(solver) = Symbol(string(typeof(solver)))

function solve_linear(solver, reward::Dict{T, Float64}) where T
  # Generic implementation that should work for all solvers following the same conventions:
  #  - `has_lp_formulation` returns `true` for the solver
  #  - `supports_solve_budgeted_linear` returns `true` for the solver
  #  - `get_lp_formulation` returns a LP-like formulation (may even be a MIP)
  # Lazy formulations are allowed.
  @assert has_lp_formulation(solver)
  @assert supports_solve_budgeted_linear(solver)

  m, obj, vars = get_lp_formulation(solver, reward)
  @objective(m, Max, obj)

  # Maybe solve_budgeted_linear has already been called on this model.
  key = _generic_lp_key(solver)
  if key in keys(m.ext)
    budget_constraint = m.ext[key][:budget_constraint]
    set_normalized_rhs(budget_constraint, 0)
  end

  set_silent(m)
  optimize!(m)

  if termination_status(m) != MOI.OPTIMAL
    return T[]
  end
  return T[arm for arm in keys(reward) if value(vars[arm]) > 0.5]
end

function solve_budgeted_linear(solver,
                               reward::Dict{T, Float64},
                               weight::Dict{T, S},
                               budget::Int) where {T, S<:Number} # Handle both Int and Float64
  # Generic implementation that should work for all solvers following the same conventions:
  #  - `has_lp_formulation` returns `true` for the solver
  #  - `supports_solve_budgeted_linear` returns `true` for the solver
  #  - `get_lp_formulation` returns a LP-like formulation (may even be a MIP)
  #    whose variables are alse used for the budget constraint
  # Lazy formulations are allowed.
  @assert has_lp_formulation(solver)
  @assert supports_solve_budgeted_linear(solver)

  m, obj, vars = get_lp_formulation(solver, reward)
  @objective(m, Max, obj)

  # Add the budget constraint (or change the existing constraint).
  key = _generic_lp_key(solver)
  if key in keys(m.ext)
    budget_constraint = m.ext[key][:budget_constraint]
    set_normalized_rhs(budget_constraint, budget)
  else
    budget_constraint = @constraint(m, sum(weight[arm] * vars[arm] for arm in keys(reward)) >= budget)
    m.ext[key] = Dict(:budget_constraint => budget_constraint)
  end

  set_silent(m)
  optimize!(m)

  if termination_status(m) != MOI.OPTIMAL
    return T[]
  end
  return T[arm for arm in keys(reward) if value(vars[arm]) > 0.5]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""RuntimeClass defines a class of container runtime supported in the cluster. The RuntimeClass is used to determine which container runtime is used to run all containers in a pod. RuntimeClasses are (currently) manually defined by a user or cluster provisioner, and referenced in the PodSpec. The Kubelet is responsible for resolving the RuntimeClassName reference before running the pod.  For more details, see https://git.k8s.io/enhancements/keps/sig-node/runtime-class.md

    IoK8sApiNodeV1beta1RuntimeClass(;
        apiVersion=nothing,
        handler=nothing,
        kind=nothing,
        metadata=nothing,
        overhead=nothing,
        scheduling=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - handler::String : Handler specifies the underlying runtime and configuration that the CRI implementation will use to handle pods of this class. The possible values are specific to the node & CRI configuration.  It is assumed that all handlers are available on every node, and handlers of the same name are equivalent on every node. For example, a handler called \"runc\" might specify that the runc OCI runtime (using native Linux containers) will be used to run the containers in a pod. The Handler must conform to the DNS Label (RFC 1123) requirements, and is immutable.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
    - overhead::IoK8sApiNodeV1beta1Overhead : Overhead represents the resource overhead associated with running a pod for a given RuntimeClass. For more details, see https://git.k8s.io/enhancements/keps/sig-node/20190226-pod-overhead.md This field is alpha-level as of Kubernetes v1.15, and is only honored by servers that enable the PodOverhead feature.
    - scheduling::IoK8sApiNodeV1beta1Scheduling : Scheduling holds the scheduling constraints to ensure that pods running with this RuntimeClass are scheduled to nodes that support it. If scheduling is nil, this RuntimeClass is assumed to be supported by all nodes.
"""
mutable struct IoK8sApiNodeV1beta1RuntimeClass <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    handler::Any # spec type: Union{ Nothing, String } # spec name: handler
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    overhead::Any # spec type: Union{ Nothing, IoK8sApiNodeV1beta1Overhead } # spec name: overhead
    scheduling::Any # spec type: Union{ Nothing, IoK8sApiNodeV1beta1Scheduling } # spec name: scheduling

    function IoK8sApiNodeV1beta1RuntimeClass(;apiVersion=nothing, handler=nothing, kind=nothing, metadata=nothing, overhead=nothing, scheduling=nothing)
        o = new()
        validate_property(IoK8sApiNodeV1beta1RuntimeClass, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiNodeV1beta1RuntimeClass, Symbol("handler"), handler)
        setfield!(o, Symbol("handler"), handler)
        validate_property(IoK8sApiNodeV1beta1RuntimeClass, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiNodeV1beta1RuntimeClass, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiNodeV1beta1RuntimeClass, Symbol("overhead"), overhead)
        setfield!(o, Symbol("overhead"), overhead)
        validate_property(IoK8sApiNodeV1beta1RuntimeClass, Symbol("scheduling"), scheduling)
        setfield!(o, Symbol("scheduling"), scheduling)
        o
    end
end # type IoK8sApiNodeV1beta1RuntimeClass

const _property_map_IoK8sApiNodeV1beta1RuntimeClass = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("handler")=>Symbol("handler"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("overhead")=>Symbol("overhead"), Symbol("scheduling")=>Symbol("scheduling"))
const _property_types_IoK8sApiNodeV1beta1RuntimeClass = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("handler")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("overhead")=>"IoK8sApiNodeV1beta1Overhead", Symbol("scheduling")=>"IoK8sApiNodeV1beta1Scheduling")
Base.propertynames(::Type{ IoK8sApiNodeV1beta1RuntimeClass }) = collect(keys(_property_map_IoK8sApiNodeV1beta1RuntimeClass))
Swagger.property_type(::Type{ IoK8sApiNodeV1beta1RuntimeClass }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNodeV1beta1RuntimeClass[name]))}
Swagger.field_name(::Type{ IoK8sApiNodeV1beta1RuntimeClass }, property_name::Symbol) =  _property_map_IoK8sApiNodeV1beta1RuntimeClass[property_name]

function check_required(o::IoK8sApiNodeV1beta1RuntimeClass)
    (getproperty(o, Symbol("handler")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiNodeV1beta1RuntimeClass }, name::Symbol, val)
end


text:
using Test
using PDBTools

include("./select.jl")
include("./properties.jl")


text:
# Tabular Function file

using LinearAlgebra

struct LinearFunction{TP, TB} <: AbstractFuncApprox where {TP<:Real,TB<:AbstractFuncApprox}
    θ::Array{TP, 2}
    ϕ::TB
    feats::Array{TP, 1}
    LinearFunction(input_dim::Int, output_dim::Int) = new{Float64, IdentityBasis}(zeros(Float64,input_dim, output_dim), IdentityBasis(input_dim), zeros(Float64, 1))
    LinearFunction(::Type{T}, input_dim::Int, output_dim::Int) where {T<:Real} = new{T, IdentityBasis}(zeros(T,input_dim, output_dim), IdentityBasis(input_dim), zeros(T, 1))
    LinearFunction(::Type{T}, basis::AbstractFuncApprox, output_dim::Int) where {T<:Real} = new{T, typeof(basis)}(zeros(T,get_num_outputs(basis), output_dim), clone(basis), zeros(T, get_num_outputs(basis)))
    LinearFunction(θ::Array{T,2}, basis::AbstractFuncApprox, feats::Array{T,1}) where {T<:Real} = new{T, typeof(basis)}(θ, basis, feats)
end

function call(fun::LinearFunction{TP,IdentityBasis}, x::Int)::Array{TP,1} where {TP<:Real}
    return fun.θ[x, :]
end

function call!(out::Array{TP, 1}, fun::LinearFunction{TP,IdentityBasis}, x::Int) where {TP<:Real}
    out .= fun.θ[x, :]
end

function call(fun::LinearFunction{TP,IdentityBasis}, x::Int, y::Int)::TP where{TP<:Real}
    return fun.θ[x, y]
end

function call(fun::LinearFunction{TP,IdentityBasis}, x::Array{TP,1})::Array{TP,1} where {TP<:Real}
    return fun.θ'*x
end

function call!(out::Array{TP, 1}, fun::LinearFunction{TP,IdentityBasis}, x::Array{TP,1}) where {TP<:Real}
    # out .= fun.θ'*x
    for i in 1:size(fun.θ)[2]
        out[i] .= dot(view(fun.θ, :, i), x)
    end
end

function call(fun::LinearFunction{TP}, x::Array{TP,1})::Array{TP,1} where {TP<:Real}
    call!(fun.feats, fun.ϕ, x)
    return fun.θ'*fun.feats
end

function call!(out::Array{TP, 1}, fun::LinearFunction{TP}, x::Array{TP,1}) where {TP<:Real}
    call!(fun.feats, fun.ϕ, x) #
    # out .= fun.θ'*fun.feats
    for i in 1:size(fun.θ)[2]
        out[i] = dot(view(fun.θ, :, i), fun.feats)
    end
end

function call(fun::LinearFunction{TP,IdentityBasis}, x::Array{TP,1}, y::Int)::TP where {TP<:Real}
    return dot(view(fun.θ,:, y),x)
end

function call(fun::LinearFunction{TP}, x::Array{TP,1}, y::Int) where {TP<:Real}
    call!(fun.feats, fun.ϕ, x)
    return dot(view(fun.θ,:, y), fun.feats)
end

function gradient!(grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x::Int) where {TP<:Real}
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    grad[x, :] .= 1.
end

function gradient!(grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x::Int, y::Int) where {TP<:Real}
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    grad[x, y] = 1.
end

function gradient!(grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x::Array{TP, 1}) where {TP<:Real}
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    grad .= reshape(repeat(x, size(fun.θ)[end]), size(fun.θ))
end

function gradient!(grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x::Array{TP, 1}, y::Int) where {TP<:Real}
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    grad[:, y] .= x
end

function gradient!(grad::Array{TP}, fun::LinearFunction{TP}, x::Array{TP, 1}) where {TP<:Real}
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    call!(fun.feats, fun.ϕ, x)
    grad .= reshape(repeat(fun.feats, size(fun.θ)[end]), size(fun.θ))
end

function gradient!(grad::Array{TP}, fun::LinearFunction{TP}, x::Array{TP, 1}, y::Int) where {TP<:Real}
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    call!(fun.feats, fun.ϕ, x)
    grad[:, y] .= fun.feats
end

function gradient(fun::LinearFunction{TP}, x) where {TP<:Real}
    grad::Array{TP, 2} = zeros(TP, size(fun.θ))
    gradient!(grad, fun, x)
    return grad
end

function gradient(fun::LinearFunction{TP}, x, y::Int) where {TP<:Real}
    grad::Array{TP, 2} = zeros(TP, size(fun.θ))
    gradient!(grad, fun, x, y)
    return grad
end

function call_gradient!(out::Array{TP, 1}, grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x) where {TP<:Real}
    call!(out, fun, x)
    gradient!(grad, fun, x)
end

function call_gradient!(out::Array{TP, 1}, grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x::Array{TP, 1}) where {TP<:Real}
    call!(out, fun, x)
    gradient!(grad, fun, x)
end

function call_gradient!(grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x, y::Int)::TP where {TP<:Real}
    out::TP = call(fun, x, y)
    gradient!(grad, fun, x, y)
    return out
end

function call_gradient!(grad::Array{TP}, fun::LinearFunction{TP, IdentityBasis}, x::Array{TP, 1}, y::Int)::TP where {TP<:Real}
    out::TP = call(fun, x, y)
    gradient!(grad, fun, x, y)
    return out
end

function call_gradient!(out::Array{TP, 1}, grad::Array{TP}, fun::LinearFunction{TP}, x::Array{TP,1}) where {TP<:Real}
    call!(out, fun, x)
    g = reshape(grad, size(fun.θ))
    fill!(g, 0.)
    g .= reshape(repeat(fun.feats, size(fun.θ)[end]), size(fun.θ))
end

function call_gradient!(grad::Array{TP}, fun::LinearFunction{TP}, x::Array{TP, 1}, y::Int) where {TP<:Real}
    out = call(fun, x, y)
    grad = reshape(grad, size(fun.θ))
    fill!(grad, 0.)
    grad[:, y] .= fun.feats
    return out
end

function get_params(fun::LinearFunction{TP}) where {TP<:Real}
    return fun.θ
end

function copy_params!(params::Array{TP}, fun::LinearFunction{TP}) where {TP<:Real}
    p = reshape(params, size(fun.θ))
    copyto!(p, fun.θ)
end

function copy_params(fun::LinearFunction{TP}) where {TP<:Real}
    params = deepcopy(fun.θ)
    return params
end

function set_params!(fun::LinearFunction{TP}, θ) where {TP<:Real}
    fun.θ .= reshape(θ,size(fun.θ))
end

function add_to_params!(fun::LinearFunction{TP}, θ) where {TP<:Real}
    fun.θ .+= reshape(θ, size(fun.θ))
end

function get_num_params(fun::LinearFunction)::Int64
    return length(fun.θ)
end

function clone(fun::LinearFunction{TP, TB})::LinearFunction{TP, TB} where {TP<:Real,TB<:AbstractFuncApprox}
    fun2 = LinearFunction(deepcopy(fun.θ), clone(fun.ϕ), zeros(TP, size(fun.feats)))
    fun2.feats .= deepcopy(fun.feats)
    return fun2
end

function get_num_outputs(fun::LinearFunction{TP})::Int where {TP<:Real}
    return size(fun.θ)[end]
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
@inline function anydifferent(xs::AbstractArray)
    length(xs) < 2 && return false
    el = xs[1]
    @inbounds for i in 2:length(xs)
        xs[i] == el || return true
    end
    false
end

@inline function scanl(f::Function, x::T, n::Int) where T
    result = Vector{T}(undef, n+1)
    result[1] = x
    @inbounds for i in 1:n
        result[i+1] = f(result[i])
    end
    result
end


text:
include("generator_utils.jl")

emoji_data       = get_emoji_data()
emoji_name_table = reverse(emoji_data)

write_dict_to_file("emoji_name_table.jl", "emoji_name_table", emoji_name_table)


text:
############################################################################################
"""
$(TYPEDEF)

Contains useful information for post-sampling analysis. Also allows to continue sampling with given sampler.

# Fields
$(TYPEDFIELDS)
"""
struct TraceInfo{
    A<:TemperingMethod,
    B<:Union{DataTune, Vector{<:DataTune}},
    D<:BaytesCore.SampleDefault,
    S<:SamplingInfo
    }
    "Tuning container for temperature tempering"
    tempertune::A
    "Tuning container for data tempering"
    datatune::B
    "Default information used for sample function"
    default::D
    "Information about trace used for postprocessing."
    sampling::S
    "Progress Log while sampling."
    progress::ProgressMeter.Progress
    function TraceInfo(
        tempertune::A,
        datatune::B,
        default::D,
        sampling::S,
        progress::ProgressMeter.Progress,
    ) where {
        A<:TemperingMethod,
        B<:Union{DataTune, Vector{<:DataTune}},
        D<:BaytesCore.SampleDefault,
        S<:SamplingInfo
    }
        ## Return info
        return new{A,B,D,S}(tempertune, datatune, default, sampling, progress)
    end
end

############################################################################################
"""
$(TYPEDEF)

Contains sampling chain and diagnostics for given algorithms.

# Fields
$(TYPEDFIELDS)
"""
struct Trace{C<:TraceInfo, A<:NamedTuple,B}
    "Model samples ~ out vector for corresponding chain, inner vector for iteration"
    val::Vector{Vector{A}}
    "Algorithm diagnostics ~ out vector for corresponding chain, inner vector for iteration"
    diagnostics::Vector{B}
    "Information about trace used for postprocessing."
    info::C
    function Trace(
        val::Vector{Vector{A}},
        diagnostics::Vector{B},
        info::C,
    ) where {A,B,C<:TraceInfo}
        ## Return trace
        return new{C,A,B}(val, diagnostics, info)
    end
end

function Trace(
    _rng::Random.AbstractRNG,
    algorithmᵛ::A,
    model::ModelWrapper,
    data::D,
    info::TraceInfo,
) where {A,D}
    @unpack iterations, Nchains = info.sampling
    ## Create Model Parameter buffer
    val = [Vector{typeof(model.val)}(undef, iterations) for _ in Base.OneTo(Nchains)]
    ## Create Diagnostics buffer for each algorithm used
    diagtypes = infer(_rng, AbstractDiagnostics, algorithmᵛ, model, data)
    diagnostics = diagnosticsbuffer(diagtypes, iterations, Nchains, algorithmᵛ)
    ## Return trace
    return Trace(val, diagnostics, info)
end

############################################################################################
#3 Proposal step
"""
$(SIGNATURES)
Propose new parameter for each model in `modelᵛ` given `data` with algorithms `algorithmᵛ`. A separate model for each chain is provided to avoid pointer issues.
If `args` is a single `SMC` algorithm, chains will not be separately allocated but instead used within the algorithm.
Note that smc still works as intended if used alongside other mcmc sampler in `args`.

# Examples
```julia
```

"""
function propose!(
    _rng::Random.AbstractRNG,
    trace::Trace{<:TraceInfo{<:BaytesCore.IterationTempering}},
    algorithmᵛ::AbstractVector,
    modelᵛ::Vector{M},
    data::D,
) where {M<:ModelWrapper,D}
    @unpack default, tempertune, datatune, sampling, progress = trace.info
    @unpack iterations, Nchains, Nalgorithms, captured = sampling
    @unpack log = default.report
    ## Propagate through data
    Base.Threads.@threads for Nchain in Base.OneTo(Nchains)
        ## Compute initial temperature
        temperature = BaytesCore.initial(tempertune)
        for iter in Base.OneTo(iterations)
            ## Update data iteration and current representation
            BaytesCore.update!(_rng, datatune[Nchain])
            dataₜ = BaytesCore.adjust(datatune[Nchain], data)
            ## Propose next step
            for Nalgorithm in Base.OneTo(Nalgorithms)
                _, trace.diagnostics[Nchain][Nalgorithm][iter] = propose!(
                    _rng,
                    algorithmᵛ[Nchain][Nalgorithm],
                    modelᵛ[Nchain],
                    dataₜ,
                    temperature,
                    captured
                )
                update!(progress, log, trace.diagnostics[Nchain][Nalgorithm][iter])
            end
            trace.val[Nchain][iter] = modelᵛ[Nchain].val
            ## Compute new temperature
            temperature = update!(tempertune, trace, algorithmᵛ, iter)
        end
    end
    return nothing
end

function propose!(
    _rng::Random.AbstractRNG,
    trace::Trace{<:TraceInfo{<:BaytesCore.JointTempering}},
    algorithmᵛ::AbstractVector,
    modelᵛ::Vector{M},
    data::D,
) where {M<:ModelWrapper,D}
    @unpack default, tempertune, datatune, sampling, progress = trace.info
    @unpack iterations, Nchains, Nalgorithms, captured = sampling
    @unpack log = default.report
    ## Compute initial temperature
    temperature = BaytesCore.initial(tempertune)
    ## Propagate through data
    for iter in Base.OneTo(iterations)
        Base.Threads.@threads for Nchain in Base.OneTo(Nchains)
            ## Update data iteration and current representation
            BaytesCore.update!(_rng, datatune[Nchain])
            dataₜ = BaytesCore.adjust(datatune[Nchain], data)
            ## Propose next step
            for Nalgorithm in Base.OneTo(Nalgorithms)
                _, trace.diagnostics[Nchain][Nalgorithm][iter] = propose!(
                    _rng,
                    algorithmᵛ[Nchain][Nalgorithm],
                    modelᵛ[Nchain],
                    dataₜ,
                    temperature,
                    captured
                )
                update!(progress, log, trace.diagnostics[Nchain][Nalgorithm][iter])
            end
            trace.val[Nchain][iter] = modelᵛ[Nchain].val
        end
        ## Compute new temperature
        temperature = update!(tempertune, trace, algorithmᵛ, iter)
    end
    return nothing
end

function propose!(
    _rng::Random.AbstractRNG,
    trace::Trace,
    algorithmᵛ::SMC,
    modelᵛ::M,
    data::D,
) where {T<:Trace,M<:ModelWrapper,D}
    @unpack default, tempertune, datatune, sampling, progress = trace.info
    @unpack iterations, Nchains, Nalgorithms, captured = sampling
    @unpack log = default.report
    ## Compute initial temperature
    temperature = BaytesCore.initial(tempertune)
    ## Propagate through data
    for iter in Base.OneTo(iterations)
        ## Update data iteration and current representation
        BaytesCore.update!(_rng, datatune)
        dataₜ = BaytesCore.adjust(datatune, data)
        ## Propagate through data
        _, trace.diagnostics[iter] = propose!(_rng, algorithmᵛ, modelᵛ, dataₜ, temperature, captured)
        for Nchain in eachindex(algorithmᵛ.particles.model)
            trace.val[Nchain][iter] = algorithmᵛ.particles.model[Nchain].val
            update!(progress, log, trace.diagnostics[iter])
        end
        temperature = update!(tempertune, trace, algorithmᵛ, iter)
    end
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Update tempering tune based on output diagnostics of last kernels in each chain.

# Examples
```julia
```

"""
function update!(tempertune::IterationTempering, trace::Trace, algorithm, iter::Integer)
    return BaytesCore.update!(tempertune, iter)
end

function update!(tempertune::JointTempering, trace::Trace, algorithm::AbstractVector, iter::Integer)
    # Update log weights in weights buffer
    for Nchain in eachindex(tempertune.weights)
        tempertune.weights[Nchain] = BaytesCore.untemper(trace.diagnostics[Nchain][end][iter].base)
    end
    # Update and return temperature - tempertune.weights will be normalized weights in original scale.
    update!(tempertune, tempertune.weights)
    return tempertune.val.current
end
function update!(tempertune::JointTempering, trace::Trace, algorithm::SMC, iter::Integer)
    # Update weights
    scaling = 1.0/trace.diagnostics[iter].base.temperature
    for Nchain in eachindex(tempertune.weights)
        tempertune.weights[Nchain] = scaling * trace.diagnostics[iter].ℓweights[Nchain]
    end
    # Update and return temperature
    update!(tempertune, tempertune.weights)
    return tempertune.val.current
end

############################################################################################
"""
$(SIGNATURES)
Safe `trace`, `model` and `algorithm` to current working directory.

# Examples
```julia
```

"""
function safetrace(trace::Trace, model::ModelWrapper, algorithm, chains::Integer, iterations::Integer, burnin::Integer)
    JLD2.jldsave(
        join((
            Base.nameof(typeof(model.id)),
            "_",
            Base.nameof(typeof(algorithm)),
            "_",
            Dates.today(),
            "_H",
            Dates.hour(Dates.now()),
            "M",
            Dates.minute(Dates.now()),
            "_Nchains",
            chains,
            "_Iter",
            iterations,
            "_Burnin",
            burnin,
            ".jld2",
        ));
        trace=trace,
        model=model,
        algorithm=algorithm,
    )
    return nothing
end

############################################################################################
#export
export TraceInfo, Trace, propose!