text:
using SparseArrays, Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

@testset "Diagonal element" begin
    k = 0
    for j in 1:10
        for i in 1:j
            k += 1
            @test MOIU.is_diagonal_vectorized_index(k) == (i == j)
        end
    end
    for side_dim in 1:10
        set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
        vec_dim = MOI.dimension(set)
        @test MOIU.side_dimension_for_vectorized_dimension(vec_dim) == side_dim
    end
end

@testset "Side dimension" begin
    for side_dim in 1:10
        set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
        vec_dim = MOI.dimension(set)
        @test MOIU.side_dimension_for_vectorized_dimension(vec_dim) == side_dim
    end
end

@testset "Constant" begin
    @test MOI.constant(MOI.EqualTo(3)) == 3
    @test MOI.constant(MOI.GreaterThan(6)) == 6
    @test MOI.constant(MOI.LessThan(2)) == 2
end

@testset "Shifts" begin
    @test MOIU.supports_shift_constant(MOI.EqualTo{Int})
    @test MOIU.shift_constant(MOI.EqualTo(3), 1) == MOI.EqualTo(4)
    @test MOIU.supports_shift_constant(MOI.GreaterThan{Int})
    @test MOIU.shift_constant(MOI.GreaterThan(6), -1) == MOI.GreaterThan(5)
    @test MOIU.supports_shift_constant(MOI.LessThan{Int})
    @test MOIU.shift_constant(MOI.LessThan(2), 2) == MOI.LessThan(4)
    @test MOIU.supports_shift_constant(MOI.Interval{Int})
    @test MOIU.shift_constant(MOI.Interval(-2, 3), 1) == MOI.Interval(-1, 4)
    @test MOIU.supports_shift_constant(MOI.ZeroOne) == false
    @test_throws MethodError MOIU.shift_constant(MOI.ZeroOne(), 1.0)
end

@testset "Dimension" begin
    @test MOI.dimension(MOI.EqualTo(3.0)) === 1
    @test MOI.dimension(MOI.Reals(8)) === 8
    @test MOI.dimension(MOI.NormInfinityCone(5)) === 5
    @test MOI.dimension(MOI.NormOneCone(5)) === 5
    @test MOI.dimension(MOI.DualExponentialCone()) === 3
    @test MOI.dimension(MOI.RelativeEntropyCone(5)) === 5
    @test MOI.dimension(MOI.NormSpectralCone(2, 3)) === 7
    @test MOI.dimension(MOI.NormNuclearCone(2, 3)) === 7
    @test MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(4)) === 10
    @test MOI.dimension(MOI.PositiveSemidefiniteConeSquare(5)) === 25
    @test MOI.dimension(MOI.RootDetConeTriangle(6)) === 22
    @test MOI.dimension(MOI.LogDetConeTriangle(6)) === 23
    @test MOI.dimension(MOI.RootDetConeSquare(4)) === 17
    @test MOI.dimension(MOI.LogDetConeSquare(4)) === 18
    @test MOI.dimension(MOI.SOS2(collect(1:6))) === 6
end

@testset "Set dot" begin
    vec = zeros(6)
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 0
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 0
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 2
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1
    vec[5] = 0
    vec[3] = 1
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1

    vec = zeros(7)
    @test MOIU.set_dot(vec, vec, MOI.RootDetConeTriangle(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.RootDetConeTriangle(3)) == 2
    vec = zeros(8)
    @test MOIU.set_dot(vec, vec, MOI.LogDetConeTriangle(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.LogDetConeTriangle(3)) == 1

    sp_vec = spzeros(6)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 0
    @test MOIU.set_dot(
        sp_vec,
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == 0
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) ==
          0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(
        sp_vec,
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == 2
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) ==
          1
    sp_vec[5] = 0
    sp_vec[3] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(
        sp_vec,
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) ==
          1

    sp_vec = spzeros(7)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.RootDetConeTriangle(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.RootDetConeTriangle(3)) == 2
    sp_vec = spzeros(8)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.LogDetConeTriangle(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.LogDetConeTriangle(3)) == 1
end

@testset "dot coefficients" begin
    vec = zeros(6)
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) ==
          vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) ==
          vec ./ 2
    vec[5] = 0
    vec[3] = 1
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) ==
          vec

    vec = zeros(7)
    @test MOIU.dot_coefficients(vec, MOI.RootDetConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.RootDetConeTriangle(3)) == vec ./ 2
    vec = zeros(8)
    @test MOIU.dot_coefficients(vec, MOI.LogDetConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.LogDetConeTriangle(3)) == vec

    sp_vec = spzeros(6)
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == sp_vec ./ 2
    sp_vec[5] = 0
    sp_vec[3] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == sp_vec

    sp_vec = spzeros(7)
    @test MOIU.dot_coefficients(sp_vec, MOI.RootDetConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.RootDetConeTriangle(3)) ==
          sp_vec ./ 2
    sp_vec = spzeros(8)
    @test MOIU.dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec
end


text:
# Not in FiniteDiff because `u` -> scalar isn't used anywhere else,
# but could be upstreamed.
mutable struct UGradientWrapper{fType,tType,P} <: Function
  f::fType
  t::tType
  p::P
end

(ff::UGradientWrapper)(uprev) = ff.f(uprev,ff.p,ff.t)

Base.@pure function determine_chunksize(u,alg::DiffEqBase.AbstractSensitivityAlgorithm)
  determine_chunksize(u,get_chunksize(alg))
end

Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
                   fx::Union{Nothing,AbstractArray{<:Number}}, alg::DiffEqBase.AbstractSensitivityAlgorithm, jac_config)
  if alg_autodiff(alg)
    if fx === nothing
      ForwardDiff.jacobian!(J, f, x)
    else
      ForwardDiff.jacobian!(J, f, fx, x, jac_config)
    end
  else
    FiniteDiff.finite_difference_jacobian!(J, f, x, jac_config)
  end
  nothing
end

function gradient!(df::AbstractArray{<:Number}, f,
                   x::Union{Number,AbstractArray{<:Number}},
                   alg::DiffEqBase.AbstractSensitivityAlgorithm, grad_config)
    if alg_autodiff(alg)
        ForwardDiff.gradient!(df, f, x, grad_config)
    else
        FiniteDiff.finite_difference_gradient!(df, f, x, grad_config)
    end
    nothing
end

"""
  jacobianvec!(Jv, f, x, v, alg, (buffer, seed)) -> nothing

``Jv <- J(f(x))v``
"""
function jacobianvec!(Jv::AbstractArray{<:Number}, f, x::AbstractArray{<:Number},
                      v, alg::DiffEqBase.AbstractSensitivityAlgorithm, config)
  if alg_autodiff(alg)
    buffer, seed = config
    TD = typeof(first(seed))
    T  = typeof(first(seed).partials)
    @. seed = TD(x, T(tuple(v)))
    f(buffer, seed)
    Jv .= ForwardDiff.partials.(buffer, 1)
  else
    buffer1, buffer2 = config
    f(buffer1,x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    f(buffer2,x)
    @. x -= ϵ*v
    @. du = (buffer2 - buffer1)/ϵ
  end
  nothing
end

function vecjacobian!(dλ, λ, p, t, S::SensitivityFunction;
                      dgrad=nothing, dy=nothing)
  _vecjacobian!(dλ, λ, p, t, S, S.sensealg.autojacvec, dgrad, dy)
  return
end

function _vecjacobian!(dλ, λ, p, t, S::SensitivityFunction, isautojacvec::Bool, dgrad, dy)
  @unpack y, sensealg = S
  prob = getprob(S)
  f = prob.f
  if isautojacvec isa Bool && !isautojacvec
    @unpack J, uf, f_cache, jac_config = S.diffcache
    if DiffEqBase.has_jac(f)
      f.jac(J,y,p,t) # Calculate the Jacobian into J
    else
      uf.t = t
      uf.p = p
      jacobian!(J, uf, y, f_cache, sensealg, jac_config)
    end
    mul!(dλ',λ',J)
    if dgrad !== nothing
      @unpack pJ, pf, paramjac_config = S.diffcache
      if DiffEqBase.has_paramjac(f)
        # Calculate the parameter Jacobian into pJ
        f.paramjac(pJ,y,prob.p,t)
      else
        pf.t = t
        pf.u = y
        jacobian!(pJ, pf, prob.p, f_cache, sensealg, paramjac_config)
      end
      mul!(dgrad',λ',pJ)
    end
    dy !== nothing && f(dy, y, p, t)
  elseif DiffEqBase.isinplace(prob)
    _vecjacobian!(dλ, λ, p, t, S, ReverseDiffVJP(), dgrad, dy)
  else
    _vecjacobian!(dλ, λ, p, t, S, ZygoteVJP(), dgrad, dy)
  end
  return
end

function _vecjacobian!(dλ, λ, p, t, S::SensitivityFunction, isautojacvec::TrackerVJP, dgrad, dy)
  @unpack y, sensealg = S
  prob = getprob(S)
  f = prob.f
  isautojacvec = get_jacvec(sensealg)
  if DiffEqBase.isinplace(prob)
    _dy, back = Tracker.forward(y, prob.p) do u, p
      out_ = map(zero, u)
      f(out_, u, p, t)
      Tracker.collect(out_)
    end
    tmp1, tmp2 = Tracker.data.(back(λ))
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(Tracker.data(_dy)))
  else
    _dy, back = Tracker.forward(y, prob.p) do u, p
      Tracker.collect(f(u, p, t))
    end
    tmp1, tmp2 = Tracker.data.(back(λ))
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(Tracker.data(_dy)))
  end
  return
end

function _vecjacobian!(dλ, λ, p, t, S::SensitivityFunction, isautojacvec::ReverseDiffVJP, dgrad, dy)
  @unpack y, sensealg = S
  prob = getprob(S)
  f = prob.f
  isautojacvec = get_jacvec(sensealg)
  tape = S.diffcache.paramjac_config

  tu, tp, tt = ReverseDiff.input_hook(tape)
  output = ReverseDiff.output_hook(tape)
  ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
  ReverseDiff.unseed!(tp)
  ReverseDiff.unseed!(tt)
  ReverseDiff.value!(tu, y)
  ReverseDiff.value!(tp, prob.p)
  ReverseDiff.value!(tt, [t])
  ReverseDiff.forward_pass!(tape)
  ReverseDiff.increment_deriv!(output, λ)
  ReverseDiff.reverse_pass!(tape)
  copyto!(vec(dλ), ReverseDiff.deriv(tu))
  dgrad !== nothing && copyto!(vec(dgrad), ReverseDiff.deriv(tp))
  dy !== nothing && copyto!(vec(dy), ReverseDiff.value(output))
  return
end

function _vecjacobian!(dλ, λ, p, t, S::SensitivityFunction, isautojacvec::ZygoteVJP, dgrad, dy)
  @unpack y, sensealg = S
  prob = getprob(S)
  f = prob.f
  isautojacvec = get_jacvec(sensealg)
  if DiffEqBase.isinplace(prob)
    _dy, back = Zygote.pullback(y, prob.p) do u, p
      out_ = Zygote.Buffer(similar(u))
      f(out_, u, p, t)
      vec(copy(out_))
    end
    tmp1,tmp2 = back(λ)
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(_dy))
  else
    _dy, back = Zygote.pullback(y, prob.p) do u, p
      vec(f(u, p, t))
    end
    tmp1,tmp2 = back(λ)
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(_dy))
  end
  return
end

function accumulate_dgdu!(dλ, y, p, t, S::SensitivityFunction)
  @unpack dg, dg_val, g, g_grad_config = S.diffcache
  if dg != nothing
    dg(dg_val,y,p,t)
  else
    g.t = t
    gradient!(dg_val, g, y, S.sensealg, g_grad_config)
  end
  dλ .+= vec(dg_val)
  return nothing
end

function build_jac_config(alg,uf,u)
  if alg_autodiff(alg)
    jac_config = ForwardDiff.JacobianConfig(uf,u,u,
                 ForwardDiff.Chunk{determine_chunksize(u,alg)}())
  else
    if diff_type(alg) != Val{:complex}
      jac_config = FiniteDiff.JacobianCache(similar(u),similar(u),
                                                 similar(u),diff_type(alg))
    else
      tmp = Complex{eltype(u)}.(u)
      du1 = Complex{eltype(u)}.(du1)
      jac_config = FiniteDiff.JacobianCache(tmp,du1,nothing,diff_type(alg))
    end
  end
  jac_config
end

function build_param_jac_config(alg,uf,u,p)
  if alg_autodiff(alg)
    jac_config = ForwardDiff.JacobianConfig(uf,u,p,
                 ForwardDiff.Chunk{determine_chunksize(p,alg)}())
  else
    if diff_type(alg) != Val{:complex}
      jac_config = FiniteDiff.JacobianCache(similar(p),similar(u),
                                                 similar(u),diff_type(alg))
    else
      tmp = Complex{eltype(p)}.(p)
      du1 = Complex{eltype(u)}.(u)
      jac_config = FiniteDiff.JacobianCache(tmp,du1,nothing,diff_type(alg))
    end
  end
  jac_config
end

function build_grad_config(alg,tf,du1,t)
  if alg_autodiff(alg)
    grad_config = ForwardDiff.GradientConfig(tf,du1,
                                    ForwardDiff.Chunk{determine_chunksize(du1,alg)}())
  else
    grad_config = FiniteDiff.GradientCache(du1,t,diff_type(alg))
  end
  grad_config
end

ext: jl


text:
# Replace Flux's LSTMCell with BasicLSTMCell from TensorFlow 1.
# It's implementation is slightly different.
# https://github.com/tensorflow/tensorflow/blob/v2.4.1/tensorflow/python/keras/layers/legacy_rnn/rnn_cell_impl.py#L786-L787
# https://github.com/FluxML/Flux.jl/blob/master/src/layers/recurrent.jl#L141-L142
# The replacement is required because the pre-trained models use BasicLSTMCell.
struct BasicLSTMCell{A,V,S}
    Wi::A
    Wh::A
    b::V
    state0::S
end

function BasicLSTMCell(in::Integer, out::Integer;
                       init = Flux.glorot_uniform,
                       initb = zeros,
                       init_state = zeros)
    cell = BasicLSTMCell(init(out * 4, in), init(out * 4, out), initb(out * 4), (init_state(out,1), init_state(out,1)))
    cell.b[gate(out, 2)] .= 1
    return cell
end

function (m::BasicLSTMCell)((h, c), x) where {A,V,T}
    b, o = m.b, size(h, 1)
    g = m.Wi*x .+ m.Wh*h .+ b
    input = σ.(gate(g, o, 1))
    cell = tanh.(gate(g, o, 2))
    forget = σ.(gate(g, o, 3))
    output = σ.(gate(g, o, 4))
    c = forget .* c .+ input .* cell
    h′ = output .* tanh.(c)
    sz = size(x)
    return (h′, c), reshape(h′, :, sz[2:end]...)
end

Flux.@functor BasicLSTMCell

Base.show(io::IO, l::BasicLSTMCell) =
    print(io, "BasicLSTMCell(", size(l.Wi, 2), ", ", size(l.Wi, 1)÷4, ")")

"""
    BasicLSTM(in::Integer, out::Integer)
Implements the BasicLSTMCell from TensorFlow 1
https://github.com/tensorflow/tensorflow/blob/v2.4.1/tensorflow/python/keras/layers/legacy_rnn/rnn_cell_impl.py#L648-L814
"""
BasicLSTM(a...; ka...) = Flux.Recur(BasicLSTMCell(a...; ka...))
Flux.Recur(m::BasicLSTMCell) = Flux.Recur(m, m.state0)


text:
"""
    agents_data_per_step(properties::Array{Symbol}, aggregators::Array)

Collect data from a `property` of agents (a `fieldname`) and apply `aggregators` function to them.

If a fieldname of agents returns an array, this will use the `mean` of the array on which to apply aggregators.

If a fieldname is `:agent`, it applies the aggregator function to the list of agents.

If a fieldname is `:model`, it applies the aggregator function to the model object.
"""
function agents_data_per_step(properties::AbstractArray{Symbol}, aggregators::AbstractArray, model::AbstractModel; step=1)
  output = Array{Any}(undef, length(properties) * length(aggregators) + 1)
  output[1] = step
  agentslen = nagents(model)
  counter = 2
  for fn in properties
    if fn == :pos  && typeof(model.agents[1].pos) <: Tuple
      temparray = [coord_to_vertex(model.agents[i], model) for i in 1:agentslen]
    elseif fn == :agent
      temparray = model.agents
    elseif fn == :model
      temparray = model
    elseif typeof(getproperty(model.agents[1], fn)) <: AbstractArray
      temparray = [mean(getproperty(model.agents[i], fn)) for i in 1:agentslen]
    else
      temparray = [getproperty(model.agents[i], fn) for i in 1:agentslen]
    end
    for agg in aggregators
      output[counter] = agg(temparray)
      counter += 1
    end
  end
  colnames = hcat(["step"], [join([string(i[1]), split(string(i[2]), ".")[end]], "_") for i in product(properties, aggregators)])
  return output, colnames
end

"""
    agents_data_per_step(propagg::Dict, model::AbstractModel; step=1)

Collect data from keys of `propagg` and apply the function in the value of each key to aggregate those data.
"""
function agents_data_per_step(propagg::Dict, model::AbstractModel; step=1)
  ncols = 1
  colnames = ["step"]
  for (k,v) in propagg
    ncols += length(v)
    for vv in v
      push!(colnames, join([string(k), vv], "_"))
    end
  end
  output = Array{Any}(undef, ncols)
  output[1] = step
  agentslen = nagents(model)
  counter = 2
  for (fn, aggs) in propagg
    if fn == :pos && typeof(model.agents[1].pos) <: Tuple
      temparray = [coord_to_vertex(model.agents[i], model) for i in 1:agentslen]
    elseif fn == :agent
      temparray = model.agents
    elseif fn == :model
      temparray = model
    elseif typeof(getproperty(model.agents[1], fn)) <: AbstractArray
      temparray = [mean(getproperty(model.agents[i], fn)) for i in 1:agentslen]
    else
      temparray = [getproperty(model.agents[i], fn) for i in 1:agentslen]
    end
    for agg in aggs
      output[counter] = agg(temparray)
      counter += 1
    end
  end
  return output, colnames
end

"""
    agents_data_complete(properties::Array{Symbol}, model::AbstractModel)

Collect data from a `property` of agents (a `fieldname`) into a dataframe.

If a fieldname of agents returns an array, this will use the `mean` of the array. If you want to record positions of the agents, use the `pos` field.
"""
function agents_data_complete(properties::Array{Symbol}, model::AbstractModel; step=1)
  # colnames = [join([string(i[1]), split(string(i[2]), ".")[end]], "_") for i in product(properties, aggregators)]
  dd = DataFrame()
  agentslen = nagents(model)
  for fn in properties
    if fn == :pos  && typeof(model.agents[1].pos) <: Tuple
      temparray = [coord_to_vertex(model.agents[i], model) for i in 1:agentslen]
    elseif typeof(getproperty(model.agents[1], fn)) <: AbstractArray
      temparray = [mean(getproperty(model.agents[i], fn)) for i in 1:agentslen]
    else
      temparray = [getproperty(model.agents[i], fn) for i in 1:agentslen]
    end
    begin
      dd[!, :id] = [i.id for i in model.agents]
    end
    fieldname = Symbol(join([string(fn), step], "_"))
    begin
      dd[!, fieldname] = temparray
    end
  end
  return dd
end

"""
    data_collector(properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer)  where T<: Integer

Used in the `step!` function.
"""
function data_collector(properties::AbstractArray{Symbol}, aggregators::AbstractArray, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer) where T<: Integer
  d, colnames = agents_data_per_step(properties, aggregators, model, step=step)
  dict = Dict(Symbol(colnames[i]) => d[i] for i in 1:length(d))
  df = DataFrame(dict)
  return df
end

"""
    data_collector(properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer, df::DataFrame) T<:Integer

Used in the `step!` function.
"""
function data_collector(properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer, df::DataFrame) where T<:Integer
  d, colnames = agents_data_per_step(properties, aggregators, model, step=step)
  dict = Dict(Symbol(colnames[i]) => d[i] for i in 1:length(d))
  push!(df, dict)
  return df
end

"""
    data_collector(propagg::Dict, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer) where T<: Integer

Used in the `step!` function.
"""
function data_collector(propagg::Dict, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer) where T<: Integer
  d, colnames = agents_data_per_step(propagg, model, step=step)
  dict = Dict(Symbol(colnames[i]) => d[i] for i in 1:length(d))
  df = DataFrame(dict)
  return df
end

"""
    data_collector(propagg::Dict, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer, df::DataFrame) where T<:Integer

Used in the `step!` function.
"""
function data_collector(propagg::Dict, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer, df::DataFrame) where T<:Integer
  d, colnames = agents_data_per_step(propagg, model, step=step)
  dict = Dict(Symbol(colnames[i]) => d[i] for i in 1:length(d))
  push!(df, dict)
  return df
end

"""
    data_collector(properties::Array{Symbol}, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer) where T<:Integer

Used in the `step!` function.
"""
function data_collector(properties::Array{Symbol}, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer) where T<:Integer
  df = agents_data_complete(properties, model, step=step)
  return df
end

"""
    data_collector(properties::Array{Symbol}, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer, df::DataFrame)  where T<:Integer
Used in the `step!` function.
"""
function data_collector(properties::Array{Symbol}, steps_to_collect_data::AbstractArray{T}, model::AbstractModel, step::Integer, df::DataFrame) where T<:Integer
  d = agents_data_complete(properties, model, step=step)
  df = join(df, d, on=:id, kind=:outer)
  return df
end

"""
    combine_columns(data::DataFrame, column_names::Array{Symbol}, aggregator::AbstractVector)

Combines columns of the data that contain the same type of info from different steps of the model into one column using an aggregator, e.g. mean. You should either supply all column names that contain the same type of data, or one name (as a string) that precedes a number in different columns, e.g. "pos_"{some number}.
"""
function combine_columns!(data::DataFrame, column_names::Array{Symbol}, aggregators::AbstractVector)
  for ag in aggregators
    d = by(data, :step, column_names => x-> (ag([getproperty(x, i) for i in column_names])))
    colname = Symbol(string(column_names[1])[1:end-1] * string(ag))
    data[!, colname] = d[!, names(d)[end]]
  end
  return data
end

"""
    combine_columns!(data::DataFrame, column_base_name::String, aggregators::AbstractVector)

Combines columns of the data that contain the same type of info from different steps of the model into one column using an aggregator, e.g. mean. You should either supply all column names that contain the same type of data, or one name (as a string) that precedes a number in different columns, e.g. "pos_"{some number}.
"""
function combine_columns!(data::DataFrame, column_base_name::String, aggregators::AbstractVector)
  column_names = vcat([column_base_name], [column_base_name*"_"*string(i) for i in 1:size(data, 2)])
  datanames = [string(i) for i in names(data)]
  final_names = Array{Symbol}(undef, 0)
  for cn in column_names
    if cn in datanames
      push!(final_names, Symbol(cn))
    end
  end
  combine_columns!(data, final_names, aggregators)
end

"""
    write_to_file(;df::DataFrame, filename::AbstractString)

Writes a DataFrame to file.
"""
function write_to_file(;df::DataFrame, filename::AbstractString)
  CSV.write(filename, df, append=false, delim="\t", writeheader=true)
end


text:
using EconPDEs, Distributions

struct WachterModel{T<: Distribution}
    # consumption process parameters
    μ::Float64
    σ::Float64

    # idiosyncratic
    λbar::Float64
    κλ::Float64
    νλ::Float64
    ZDistribution::T

    # utility parameters
    ρ::Float64
    γ::Float64
    ψ::Float64
    ϕ::Float64
end

# distirbution of jump comes from Ian Martin higher order cumulants paper
function WachterModel(;μ = 0.025, σ = 0.02, λbar = 0.0355, κλ = 0.08, νλ = 0.067, ZDistribution = Normal(-0.4, 0.25), ρ = 0.012, γ = 3.0, ψ = 1.1, ϕ = 2.6)
    WachterModel(μ, σ, λbar, κλ, νλ, ZDistribution, ρ, γ, ψ, ϕ)
end

function initialize_stategrid(m::WachterModel; n = 30)
    μ = m.μ ; σ = m.σ ; λbar = m.λbar ; κλ = m.κλ ; νλ = m.νλ ; ZDistribution = m.ZDistribution ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ
    OrderedDict(:λ => range(0.0, stop = 0.1, length = n))
end

function initialize_y(m::WachterModel, stategrid::OrderedDict)
    OrderedDict(:p => ones(length(stategrid[:λ])))
end

function (m::WachterModel)(state::NamedTuple, y::NamedTuple)
    μ = m.μ ; σ = m.σ ; λbar = m.λbar ; κλ = m.κλ ; νλ = m.νλ ; ZDistribution = m.ZDistribution ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ ; ϕ = m.ϕ
    λ = state.λ
    p, pλ, pλλ = y.p, y.pλ, y.pλλ

    # Drift and volatility of λ, p
    μλ = κλ * (λbar - λ)
    σλ = νλ * sqrt(λ)
    μp = pλ / p * μλ + 0.5 * pλλ / p * σλ^2
    σp_Zλ = pλ / p * σλ

    # Market price of risk
    κ_Zc = γ * σ
    κ_Zλ = (γ * ψ - 1) / (ψ - 1) * σp_Zλ
    η =  λ * (mgf(ZDistribution, 1) - 1 + mgf(ZDistribution, -γ) - 1 - (mgf(ZDistribution, 1 - γ) - 1))

    # Interest rate r
    r = ρ + μ / ψ - (1 + 1 / ψ) / 2 * γ * σ^2 - λ * (mgf(ZDistribution, -γ) - 1) + (1 / ψ - γ) / (1 - γ) * λ * (mgf(ZDistribution, 1 - γ) - 1) - (γ * ψ - 1) / (2 * (ψ - 1)) * σp_Zλ^2

    # Market Pricing
    pt = p * (1 / p  + μ + μp + λ * (mgf(ZDistribution, 1) - 1) - r - κ_Zc * σ - κ_Zλ * σp_Zλ - η)

    return (pt,), (μλ,), (p = p, r = r, κ_Zc = κ_Zc, κ_Zλ = κ_Zλ, η = η)
end

m = WachterModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0)

#========================================================================================

Solve for levered equity claim

========================================================================================#

function pde_levered(m, state, y, r, κ_Zc, κ_Zλ)
    μ = m.μ ; σ = m.σ ; λbar = m.λbar ; κλ = m.κλ ; νλ = m.νλ ; ZDistribution = m.ZDistribution ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ ; ϕ = m.ϕ
    λ = state.λ

    pe, peλ, peλλ = y.pe, y.peλ, y.peλλ
    μλ = κλ * (λbar - λ)
    σλ = νλ * sqrt(λ)
    μpe = peλ / pe * μλ + 0.5 * peλλ / pe * σλ^2
    σpe_Zλ = peλ / pe * σλ
    η =  λ * (mgf(ZDistribution, ϕ) - 1 + mgf(ZDistribution, -γ) - 1 - (mgf(ZDistribution, ϕ - γ) - 1))
    pet = pe * (1 / pe  + μ + μpe + λ * (mgf(ZDistribution, ϕ) - 1) - r(λ) - κ_Zc(λ) * ϕ * σ - κ_Zλ(λ) * σpe_Zλ - η)
    return (pet,), (μλ,)
end
#using Interpolations
#r = interpolate((state[:λ],), result[:r], Gridded(Linear()))
#κ_Zc = interpolate((state[:λ],), result[:κ_Zc], Gridded(Linear()))
#κ_Zλ = interpolate((state[:λ],), result[:κ_Zλ], Gridded(Linear()))
#y2, _, distance =  pdesolve((state, y) -> pde_levered(m, state, y, r, κ_Zc, κ_Zλ), stategrid, OrderedDict(:pe => y[:p]))


text:
struct KFTimeUpdate{X,P} <: AbstractTimeUpdate{X,P}
    state::X
    covariance::P
end

struct KFTUIntermediate{T}
    x_apri::Vector{T}
    p_apri::Matrix{T}
    fp::Matrix{T}
end

KFTUIntermediate(T::Type, num_x::Number) =
    KFTUIntermediate(
        Vector{T}(undef, num_x),
        Matrix{T}(undef, num_x, num_x),
        Matrix{T}(undef, num_x, num_x)
    )

KFTUIntermediate(num_x::Number) = KFTUIntermediate(Float64, num_x)

struct KFMeasurementUpdate{X,P,R,S,K} <: AbstractMeasurementUpdate{X,P}
    state::X
    covariance::P
    innovation::R
    innovation_covariance::S
    kalman_gain::K
end

struct KFMUIntermediate{T,K<:Union{<:AbstractVector{T},<:AbstractMatrix{T}}}
    innovation::Vector{T}
    innovation_covariance::Matrix{T}
    kalman_gain::K
    pht::K
    s_chol::Matrix{T}
    x_posterior::Vector{T}
    p_posterior::Matrix{T}
end

function KFMUIntermediate(T::Type, num_x::Number, num_y::Number)
    return KFMUIntermediate(
        Vector{T}(undef, num_y),
        Matrix{T}(undef, num_y, num_y),
        Matrix{T}(undef, num_x, num_y),
        Matrix{T}(undef, num_x, num_y),
        Matrix{T}(undef, num_y, num_y),
        Vector{T}(undef, num_x),
        Matrix{T}(undef, num_x, num_x)
    )
end

KFMUIntermediate(num_x::Number, num_y::Number) = KFMUIntermediate(Float64, num_x, num_y)

calc_apriori_state(x, F) = F * x
calc_apriori_covariance(P, F, Q) = F * P * F' .+ Q

function time_update(x, P, F::Union{Number, AbstractMatrix}, Q)
    x_apri = calc_apriori_state(x, F)
    P_apri = calc_apriori_covariance(P, F, Q)
    KFTimeUpdate(x_apri, P_apri)
end

function time_update!(tu::KFTUIntermediate, x, P, F::AbstractMatrix, Q)
    x_apri = calc_apriori_state!(tu.x_apri, x, F)
    P_apri = calc_apriori_covariance!(tu.p_apri, tu.fp, P, F, Q)
    KFTimeUpdate(x_apri, P_apri)
end

function measurement_update(x, P, y, H::Union{Number, AbstractVector, AbstractMatrix}, R; consider = nothing)
    ỹ = calc_innovation(H, x, y)
    PHᵀ = calc_P_xy(P, H)
    S = calc_innovation_covariance(H, P, R)
    K = calc_kalman_gain(PHᵀ, S, consider)
    x_post = calc_posterior_state(x, K, ỹ, consider)
    P_post = calc_posterior_covariance(P, PHᵀ, K, consider)
    KFMeasurementUpdate(x_post, P_post, ỹ, S, K)
end

@inline calc_P_xy(P, H) = P * H'
@inline calc_innovation(H, x, y) = y .- H * x
@inline calc_innovation_covariance(H, P, R) = H * P * H' .+ R
@inline calc_kalman_gain(PHᵀ, S, consider::Nothing) = PHᵀ / S
@inline calc_posterior_state(x, K, ỹ, consider::Nothing) = x .+ K * ỹ
@inline calc_posterior_covariance(P, PHᵀ, K, consider::Nothing) = P .- PHᵀ * K' # (I - K * H) * P ?

function measurement_update!(mu::KFMUIntermediate, x, P, y, H::AbstractMatrix, R)
    ỹ = calc_innovation!(mu.innovation, H, x, y)
    PHᵀ = calc_P_xy!(mu.pht, P, H)
    S = calc_innovation_covariance!(mu.innovation_covariance, H, PHᵀ, R)
    K = calc_kalman_gain!(mu.s_chol, mu.kalman_gain, PHᵀ, S)
    x_post = calc_posterior_state!(mu.x_posterior, x, K, ỹ)
    P_post = calc_posterior_covariance!(mu.p_posterior, P, PHᵀ, K)
    KFMeasurementUpdate(x_post, P_post, ỹ, S, K)
end

function calc_P_xy!(PHᵀ, P, H)
    mul!(PHᵀ, P, H')
    PHᵀ
end

function calc_apriori_state!(x_apri, x, F)
    x_apri .= @~ F * x
end

function calc_apriori_covariance!(P_apri, FP, P, F, Q)
    FP .= @~ F * P
    P_apri .= @~ FP * F' + Q
end

function calc_innovation!(ỹ, H, x, y)
    ỹ .= @~ -1 * H * x + y # Order is important to trigger BLAS
end

function calc_innovation_covariance!(S, H, PHᵀ, R)
    S .= @~ H * PHᵀ + R
end

function calc_kalman_gain!(S_chol, K, PHᵀ, S)
    S_chol .= S
    K .= PHᵀ
    rdiv!(K, cholesky!(Hermitian(S_chol)))
end

function calc_posterior_state!(x_posterior, x, K, ỹ)
    x_posterior .= @~ K * ỹ + x
end

function calc_posterior_covariance!(P_posterior, P, PHᵀ, K)
    P_posterior .= @~ -1 * PHᵀ * K' + P # Order is important to trigger BLAS
end


text:
using StaticLinearMaps
using StaticArrays
using Base.Test

const ℝ{L} = SVector{L,Float64}

# write your own tests here
include("testmain.jl")


text:
function formatdf(a::Array{ASCIIString,1})
    i = 0
    s = "    "
    for c in a
        s *= @sprintf "%6s" c
        i += 1
        i %= 8
        if i == 0
            s *= "\n    "
        end
    end
    return s
end

cnum = 20
d = [COLORS[rand(1:3)] for i in 1:cnum]
while d == dutchsort(d)
    d = [COLORS[rand(1:3)] for i in 1:cnum]
end

println("The original list is:")
println(formatdf(d))

print("Sorting with dutchsort, ")
@time e = dutchsort(d)
println(formatdf(e))

print("Sorting conventionally, ")
@time e = sort(d, by=x->findfirst(COLORS, x))
println(formatdf(e))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct VectorPriorityQueue{K,V}
    keys::Vector{K}
    values::Vector{V}
end

VectorPriorityQueue{K,V}() where {K,V} = VectorPriorityQueue{K,V}(K[], V[])

function DataStructures.enqueue!(pq::VectorPriorityQueue{K,V}, k::K, v::V) where {K,V}
    left = searchsortedfirst(pq.values, v)
    insert!(pq.keys, left, k)
    insert!(pq.values, left, v)
    return nothing
end

function DataStructures.dequeue!(pq::VectorPriorityQueue)
    k = popfirst!(pq.keys)
    popfirst!(pq.values)
    return k
end

Base.first(pq::VectorPriorityQueue) = first(pq.keys)
Base.peek(pq::VectorPriorityQueue) = first(pq)
Base.length(pq::VectorPriorityQueue) = length(pq.keys)
Base.keys(pq::VectorPriorityQueue) = pq.keys
Base.haskey(pq::VectorPriorityQueue{K}, k::K) where {K} = k in keys(pq)
Base.values(pq::VectorPriorityQueue) = pq.values
Base.isempty(pq::VectorPriorityQueue) = isempty(pq.keys)
Base.pairs(pq::VectorPriorityQueue) = (k => v for (k, v) in zip(pq.keys, pq.values))

function Base.deleteat!(pq::VectorPriorityQueue, args...)
    deleteat!(pq.keys, args...)
    deleteat!(pq.values, args...)
    return nothing
end


text:
module ParamsIO
#
# This is a failed experiment in reading and writing the model parameters in JSON - one
# of several time-consuming failures. Kept around because json I/O would be lovely
# if it could be made to work.
#
# THIS DOESN'T WORK !!!1
# json3 turns the structures below the top-level into Dicts
# and I don't know why.
#
# See scripts/serialisation_experiments.jl
#

    using ScottishTaxBenefitModel
    using .STBParameters
    using .Definitions

    using JSON3
    using JSON
    using StructTypes
    using TimeSeries

    const T = Float64 # FIXME change this as needed

    export load, to_file, from_file, toJSON, fromJSON

    StructTypes.StructType(::Type{TimeSeries.TimeArray}) = StructTypes.ArrayType()

    StructTypes.StructType(::Type{TimeSeries.TimeArray{Int64,2,Dates.Date,Array{Int64,2}}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{IncomeTaxSys{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{AgeLimits}) = StructTypes.Struct()

    StructTypes.StructType(::Type{NationalInsuranceSys{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{PersonalAllowances{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{Premia{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{WorkingTaxCredit{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{ChildTaxCredit{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{HoursLimits}) = StructTypes.Struct()

    StructTypes.StructType(::Type{IncomeRules{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{MinimumWage{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{SavingsCredit{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{HousingBenefits{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{LegacyMeansTestedBenefitSystem{T}}) = StructTypes.Struct()

    StructTypes.StructType(::Type{HousingRestrictions{T}}) = StructTypes.Struct()

    ## .. and so on

    StructTypes.StructType(::Type{TaxBenefitSystem{T}}) = StructTypes.Struct()

    function to_file( filename :: String, t :: TaxBenefitSystem )
        io = open( filename, "w")
        JSON3.pretty( io, JSON3.write( t ), 4 )
        close( io )
    end

    function from_file( filename :: String )::TaxBenefitSystem
        io = open( filename, "r")
        t = JSON3.read( io, TaxBenefitSystem{T} )
        close( io )
        return t
    end

    function toJSON( t :: TaxBenefitSystem{T}) :: String
        JSON3.write( t )
    end

    function fromJSON( s :: String ) :: TaxBenefitSystem
        t = TaxBenefitSystem{T}()# JSON3.read( s, TaxBenefitSystem{T} )
        #sp :: Dict = JSON.parse( s )
        #t.it = JSON3.read( sp["it"], IncomeTaxSys{T} )
    end
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
@testset "operations" begin
    mock_db = DB(joinpath(@__DIR__, "mock_db"))

    @testset "getindex" begin
        @test mock_db['n', "'hood"].word == "'hood"
        @test mock_db["'hood", 'n'].word == "'hood"
    end

    @testset "synsets" begin
        ss = synsets(mock_db, mock_db['n', "'hood"])

        @test length(ss) == 1
        @test ss[1].gloss == "(slang) a neighborhood"
    end

    @testset "expanded_hypernyms" begin
        ss = synsets(mock_db, mock_db['n', "'hood"])[1]
        @test length(expanded_hypernyms(mock_db, ss)) == 8
    end

    @testset "antonyms" begin
        ss = synsets(mock_db, mock_db['n', "'hood"])[1]
        @test length(antonyms(mock_db, ss)) == 0
    end

    @testset "multiple hypernyms" begin
        ss = synsets(mock_db, mock_db['n', "lactate"])[1]
        @test length(hypernyms(mock_db, ss)) == 2
    end

    @testset "similar tos" begin
        aquatic = synsets(mock_db, mock_db['a', "aquatic"])[1]
        marine = synsets(mock_db, mock_db['a', "marine"])[1]
        @test WordNet.relation(mock_db, aquatic, WordNet.SIMILAR_TO) == [marine]
        @test WordNet.relation(mock_db, marine, WordNet.SIMILAR_TO) == [aquatic]
    end
end


text:
@Base.kwdef struct Color
    color::Union{Symbol, Int} = :normal
    bold::Bool = false
    underline::Bool = false
    blink::Bool = false
    reverse::Bool = false
    hidden::Bool = false
end

function Color(s::AbstractString)
    return Color(split(s, ":"))
end

Color(c::Color) = c

function Color(v::Vector)
    bold = false
    underline = false
    blink = false
    reverse = false
    hidden = false
    color = :normal
    for item in v
        if item == "bold"
            bold = true
        elseif item == "underline"
            underline = true
        elseif item == "blink"
            blink = true
        elseif item == "reverse"
            reverse = true
        elseif item == "hidden"
            hidden = true
        elseif !isempty(item)
            val = tryparse(Int, item)
            if val === nothing
                color = Symbol(item)
            else
                color = val
            end
        end
    end

    return Color(color, bold, underline, blink, reverse, hidden)
end

function coal(a, b, c)
    a === nothing ? b === nothing ? c : b : a
end


text:
@testset "Attributes" begin
    @testset "is_set_by_optimize" begin
        @test MOI.is_set_by_optimize(MOI.TerminationStatus())
        @test !MOI.is_set_by_optimize(MOI.ConstraintSet())
        @test !MOI.is_set_by_optimize(MOI.ObjectiveSense())
        @test MOI.is_set_by_optimize(MOI.CallbackNodeStatus(1))
    end
    @testset "is_copyable" begin
        @test !MOI.is_copyable(MOI.TerminationStatus())
        @test !MOI.is_copyable(MOI.ConstraintSet())
        @test MOI.is_copyable(MOI.ObjectiveSense())
    end
    @testset "supports" begin
        model = DummyModel()
        @test_throws ArgumentError MOI.supports(model, MOI.TerminationStatus())
        @test_throws ArgumentError begin
            MOI.supports(model, MOI.ConstraintSet(),
                         MOI.ConstraintIndex{MOI.SingleVariable,
                                             MOI.EqualTo{Float64}})
        end
        @test MOI.supports(model, MOI.ObjectiveSense())
    end
    @testset "set vector" begin
        attr = MOI.VariablePrimalStart()
        err = DimensionMismatch("Number of indices (1) does not match the " *
                                "number of values (2) set to `$attr`.")
        model = DummyModel()
        x = MOI.VariableIndex(1)
        @test_throws err MOI.set(model, MOI.VariablePrimalStart(), [x],
                                 ones(2))
    end
end

@testset "test_integration_compute_conflict" begin
    optimizer = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.Model{Float64}(),
        MOI.Bridges.full_bridge_optimizer(optimizer, Float64),
    )
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(0.0))
    c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    MOI.optimize!(model)
    @test MOI.get(optimizer, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    MOI.set(optimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
    MOI.set(
        optimizer,
        MOI.ConstraintConflictStatus(),
        MOI.get(
            optimizer,
            MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.LessThan{Float64}}(),
        )[1],
        MOI.NOT_IN_CONFLICT,
    )
    MOI.set(
        optimizer,
        MOI.ConstraintConflictStatus(),
        MOI.get(
            optimizer,
            MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.GreaterThan{Float64}}(),
        )[1],
        MOI.IN_CONFLICT,
    )
    MOI.compute_conflict!(model)
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.NOT_IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.IN_CONFLICT
end

MOI.Utilities.@model(
    OnlyScalarConstraints,
    (),
    (MOI.GreaterThan, MOI.LessThan),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)

@testset "test_integration_compute_conflict" begin
    optimizer = MOI.Utilities.MockOptimizer(OnlyScalarConstraints{Float64}())
    model = MOI.Bridges.full_bridge_optimizer(optimizer, Float64)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.Interval(0.0, 1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(optimizer, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    MOI.set(optimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
    MOI.compute_conflict!(model)
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test_throws ArgumentError MOI.get(model, MOI.ConstraintConflictStatus(), c)
end


text:
module Weave
using Compat

"""
`list_out_formats()`

List supported output formats
"""
function list_out_formats()
  for format = keys(formats)
      println(string(format,": ",  formats[format].description))
  end
end

"""
`tangle(source ; out_path=:doc, informat="noweb")`

Tangle source code from input document to .jl file.

* `informat`: `"noweb"` of `"markdown"`
* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document, `:pwd`: Julia working directory,  `"somepath"`, directory name as a string e.g `"/home/mpastell/weaveout"`
or filename as string e.g. ~/outpath/outfile.jl.
"""
function tangle(source ; out_path=:doc, informat=:auto)
    doc = read_doc(source, informat)
    doc.cwd = get_cwd(doc, out_path)

    outname = get_outname(out_path, doc, ext = "jl")

    open(outname, "w") do io
    for chunk in doc.chunks
      if typeof(chunk) == CodeChunk
          options = merge(rcParams[:chunk_defaults], chunk.options)
          if options[:tangle]
            write(io, chunk.content*"\n")
          end
      end
    end
  end
  doc.cwd == pwd()  && (outname = basename(outname))
  info("Writing to file $outname")
end

"""
`function weave(source ; doctype = :auto, plotlib="Gadfly",
        informat=:auto, out_path=:doc, fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache=:off)`

Weave an input document to output file.

* `doctype`: :auto = set based on file extension or specify one of the supported formats.
  See `list_out_formats()`
* `plotlib`: `"PyPlot"`, `"Gadfly"` or `nothing`
* `informat`: :auto = set based on file extension or set to  `"noweb"`, `"markdown"` or  `script`
* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document, `:pwd`: Julia working directory, `"somepath"`: output directory as a String e.g `"/home/mpastell/weaveout"` or filename as string e.g. ~/outpath/outfile.tex.
* `fig_path`: where figures will be generated, relative to out_path
* `fig_ext`: Extension for saved figures e.g. `".pdf"`, `".png"`. Default setting depends on `doctype`.
* `cache_path`: where of cached output will be saved.
* `cache`: controls caching of code: `:off` = no caching, `:all` = cache everything,
  `:user` = cache based on chunk options, `:refresh`, run all code chunks and save new cache.

**Note:** Run Weave from terminal and not using IJulia, Juno or ESS, they tend to mess with capturing output.
"""
function weave(source ; doctype = :auto, plotlib=:auto,
        informat=:auto, out_path=:doc, fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache=:off)

    doc = read_doc(source, informat)
    doc = run(doc, doctype = doctype, plotlib=plotlib,
            out_path=out_path,
            fig_path = fig_path, fig_ext = fig_ext, cache_path = cache_path, cache=cache)
    formatted = format(doc)

    formatted = join(formatted, "\n")

    outname = get_outname(out_path, doc)

    open(outname, "w") do io
        write(io, formatted)
    end

    #Convert using pandoc
    if doc.doctype == "md2html"
        outname = get_outname(out_path, doc, ext = "html")
        pandoc2html(formatted, doc, outname)
    elseif doc.doctype == "md2pdf"
        outname = get_outname(out_path, doc, ext = "pdf")
        pandoc2pdf(formatted, doc, outname)
    end

    doc.cwd == pwd() && (outname = basename(outname))

    info("Report weaved to $outname")
end

function weave(doc::AbstractString, doctype::AbstractString)
    weave(doc, doctype=doctype)
end

#Hooks to run before and after chunks, this is form IJulia,
#but note that Weave hooks take the chunk as input
const preexecute_hooks = Function[]
push_preexecute_hook(f::Function) = push!(preexecute_hooks, f)
pop_preexecute_hook(f::Function) = splice!(preexecute_hooks, findfirst(pretexecute_hooks, f))

const postexecute_hooks = Function[]
push_postexecute_hook(f::Function) = push!(postexecute_hooks, f)
pop_postexecute_hook(f::Function) = splice!(postexecute_hooks, findfirst(postexecute_hooks, f))

export weave, list_out_formats, tangle,
        set_chunk_defaults, get_chunk_defaults, restore_chunk_defaults

include("config.jl")
include("chunks.jl")
include("display_methods.jl")
include("readers.jl")
include("run.jl")
include("cache.jl")
include("formatters.jl")
include("pandoc.jl")
end


text:
module CTCopy

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

##################Dependencies

using  DataFrames, Distributions, StatsBase, GLM, CategoricalArrays,
  Dates, NLopt, ForwardDiff, Formatting, LinearAlgebra

######################Methods####################
export CTLM, #Regression methods
  CT2SLS,
  project!,
  getResid!,
  pullModelMatrix,
  getCoeff!,
  getHomoskedΣ!,
  getModWhiteΣ!,
  getWhiteΣ!,
  getNeweyWest!,
  getModWhiteΣSlow,
  getHomoskedΣSlow,
  getWhiteΣSlow,
  getNeweyWestSlow,
  dropNullsFromDF,
  CTExpr,
  CTSym,
  CTQR,
  get1stStage,
  getTerm,
  getR,
  getR²,
  getNeweyWestFunc,
  getClustered!,

  texTable, #IO Mthods
  writeTables2File,
  array2String,
  num2Str,
  vec2String,

  skewnessStat, #CTstat methods
  kurtosisStat,
  JBStat,
  UnivariateTest,
  autoρ,
  ljungBox,
  test,
  pAutoρ,
  ARMA,
  ARMAFillYHat!,
  ARMALogLike,
  ARMAEstimate,
  generalizeARMA!,
  narrowARMA!,
  lagDF!,

  CTΦInv,#CTNumerical
  CTΦ,

  VAR, #CTVAR
  VAREstimate,
  VARPrediction,
  varΣ,
  propagateImpulse

##################Custom types
###NOTE: 0.7 Compat Hacks, delte on upgrade
#const Nothing = Void
#const LinearAlgebra = Base.LinAlg #NOTE:Compat Hack
#const AbstractRange = Range #NOTE compat hack

###CTReg Types
const CTExpr = Union{Symbol,Expr,Nothing}
const CTSym = Union{Symbol,Nothing}
const CTData = Union{Float64,Int,Date,Symbol, CategoricalValue, Missing, Nothing}#,
abstract type CTModel end

###VAR Types
const Structure = Union{Nothing, Vector{Float64}}
const CoefficientErrors = Union{Nothing, Vector{Float64}}
const RegressionSS = Union{Nothing, Float64}
const RegressionN = Union{Nothing, Int}

################Constants

#shared cosntants
const DEFAULT_STAR_LEGEND = "*p<0.1, **p<0.05, ***p<0.01"
const DEBUG_CTMOD = false
const DEFAULT_DECIMALS = 4

#Stat constants
const MAX_TIME = 30.
const F_TOL_ABS = 10. ^ -12.
const F_TOL_REL = 10. ^ -10.
const OPT_BOUND = 2.

#numerical constants
const sqrt2Inv = 1.0/(2.0^.5)
const sqrt2PiInv = 1.0/(2.0*π)^0.5
#  MFloat64, MInt, MBool, MSymbol}

###################Files############

include("CTReg.jl")
include("CTIO.jl")
include("CTStat.jl")
include("CTVAR.jl")
include("CTNumerical.jl")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#MDDatasets: Measure binary signals
#-------------------------------------------------------------------------------
#=NOTE:
 - Many functions make use of "cross" functions.
 - Naming assumes x-values are "time".
=#

#==Helper functions
===============================================================================#

function _buildckvector(ctx::String, tstart::Real, tstop::Real, tsample::Real)
	if isinf(tstart)
		msg = "$ctx: Must specify finite clock start time."
		throw(ArgumentError(msg))
	end

	return collect(tstart:tsample:tstop)
end

#==Main functions
===============================================================================#

#measdelay: Measure delay between crossing events of two signals:
#-------------------------------------------------------------------------------
function measdelay(dref::DataF1, dmain::DataF1; nmax::Integer=0,
	tstart_ref::Real=-Inf, tstart_main::Real=-Inf,
	xing_ref::CrossType=CrossType(), xing_main::CrossType=CrossType())

	xref = xcross(dref, nmax=nmax, xstart=tstart_ref, allow=xing_ref)
	xmain = xcross(dmain, nmax=nmax, xstart=tstart_main, allow=xing_main)
	npts = min(length(xref), length(xmain))
	delay = xmain.y[1:npts] - xref.y[1:npts]
	x = xref.x[1:npts]
	return DataF1(x, delay)
end
function measdelay(::DS{:event}, dref::DataF1, dmain::DataF1, args...; kwargs...)
	d = measdelay(dref, dmain, args...;kwargs...)
	return DataF1(collect(1:length(d.x)), d.y)
end

#measck2q: Measure clock-to-Q delay
#-------------------------------------------------------------------------------

#=_measck2q: Core algorithm to measure clock-to-Q delay
Inputs
   delaymin: Minimum circuit delay used to align clock & q edges
=#
function _measck2q(xingck::Vector, xingq::Vector, delaymin::Real)
	xq = copy(xingq) .- delaymin
	qlen = length(xq) #Maximum # of q-events
	x = copy(xingq) #Allocate space for delay starts
	Δ = copy(xingq) #Allocate space for delays
	cklen = length(xingck)
	npts = 0
	stop = false

	if qlen < 1 || cklen < 2 #Need to know if q is between 2 ck events
		xt = eltype(x)
		return DataF1(Vector{xt}(), Vector{xt}())
	end

	iq = 1
	ick = 1
	xqi = xq[iq]
	xcki = xingck[ick]
	xcki1 = xingck[ick+1]
	while xcki > xqi #Find first q event after first ck event.
		iq += 1
		xqi = xq[iq]
	end

	while iq <= qlen
		xqi = xq[iq]
		#Find clock triggering q event:
		while xcki1 <= xqi
			ick +=1
			if ick < cklen
				xcki = xingck[ick]
				xcki1 = xingck[ick+1]
			else #Not sure if this xqi corresponds to xcki
				stop = true
				break
			end
		end
		if stop; break; end

		#Compute delay (re-insert removed minimum delay):
		npts += 1
		x[npts] = xcki
		Δ[npts] = xqi - xcki + delaymin

		#Consider next q transition:
		iq += 1
	end

	return DataF1(x[1:npts], Δ[1:npts])
end

#=Measure clock-to-Q delay with non-ideal clock.
Inputs
   delaymin: Minimum circuit delay used to align clock & q edges
             Needed when delay is larger than time between ck events.
=#
function measck2q(ck::DataF1, q::DataF1; delaymin::Real=0,
	tstart_ck::Real=-Inf, tstart_q::Real=-Inf,
	xing_ck::CrossType=CrossType(), xing_q::CrossType=CrossType())

	xingck = xcross(ck, xstart=tstart_ck, allow=xing_ck)
	xingq = xcross(q, xstart=tstart_q, allow=xing_q)
	return _measck2q(xingck.x, xingq.x, delaymin)
end
#Measure clock-to-Q delay an ideal sampling clock (tsample).
function measck2q(q::DataF1, tsample::Real; delaymin::Real=0,
	tstart_ck::Real=-Inf, tstart_q::Real=-Inf,
	xing_q::CrossType=CrossType())

	xingck = _buildckvector("measck2q", tstart_ck, (q.x[end]+tsample), tsample)
	xingq = xcross(q, xstart=tstart_q, allow=xing_q)
	return _measck2q(xingck, xingq.x, delaymin)
end

#Measure rise/fall times of a signal
#-------------------------------------------------------------------------------
function measedgewidth(sig::DataF1, tstart::Real, thresh1::Real, thresh2::Real, cross::CrossType)
	x1 = xveccross(sig-thresh1, 0, tstart, cross)
	tstart = x1[1]
	x2 = xveccross(sig-thresh2, 0, tstart, cross)
	_inf = convert(eltype(x1), Inf)

	n1 = length(x1)
	n2 = length(x2)
	y = similar(x1)

	i2 = 1
	x2i = x2[i2]
	for i in 1:n1
		x1i = x1[i]
		while x2i < x1i
			if i2 <= n2
				x2i = x2[i2]
				i2 += 1
			else
				x2i = _inf
			end
		end

		y[i] = x2i - x1i
	end
	#Clamp rise events that failed to cross thresh2:
	for i in 1:(n1-1)
		Δ = x1[i+1] - x1[i]
		if Δ < y[i]
			y[i] = Δ
		end
	end
	#Remove last point if infinite:
	if isinf(y[end])
		n1 -= 1
		resize!(x1, n1)
		resize!(y, n1)
	end
	return DataF1(x1, y)
end

function _measrise(sig::DataF1, tstart::Real, lthresh::Real, hthresh::Real)
	return measedgewidth(sig, tstart, lthresh, hthresh, CrossType(:rise))
end

_measrise(sig::DataF1, tstart::Real, lthresh::Real, hthresh) =
	throw("measrise: Must provide a real value for hthresh")
_measrise(sig::DataF1, tstart::Real, lthresh, hthresh::Real) =
	throw("measrise: Must provide a real value for lthresh")
measrise(sig::DataF1; tstart::Real=-Inf, lthresh=nothing, hthresh=nothing) =
	_measrise(sig, tstart, lthresh, hthresh)

function _measfall(sig::DataF1, tstart::Real, lthresh::Real, hthresh::Real)
	return measedgewidth(sig, tstart, hthresh, lthresh, CrossType(:fall))
end

_measfall(sig::DataF1, tstart::Real, lthresh::Real, hthresh) =
	throw("measrise: Must provide a real value for hthresh")
_measfall(sig::DataF1, tstart::Real, lthresh, hthresh::Real) =
	throw("measrise: Must provide a real value for lthresh")
measfall(sig::DataF1; tstart::Real=-Inf, lthresh=nothing, hthresh=nothing) =
	_measfall(sig, tstart, lthresh, hthresh)

#measskew: Measure skew statistics of a signal
#-------------------------------------------------------------------------------
function _getskewstats(Δr::DataMD, Δf::DataMD)
	μΔr = mean(Δr)
	μΔf = mean(Δf)
	μΔ = (μΔr + μΔf) / 2
	Δmax = max(maximum(Δr), maximum(Δf))
	Δmin = min(minimum(Δr), minimum(Δf))

	#Use ASCII symbols to avoid issues with UTF8:
	return Dict{Symbol, Any}(
		:mean_delrise => μΔr,
		:min_delrise => minimum(Δr),
		:max_delrise => maximum(Δr),
		:mean_delfall => μΔf,
		:min_delfall => minimum(Δf),
		:max_delfall => maximum(Δf),
		:mean_del => μΔ,
		:mean_skew => μΔf-μΔr,
		:max_skew => Δmax - Δmin,
		:std_delrise => std(Δr),
		:std_delfall => std(Δf),
	)
end

#Measure delay skew between a signal and its reference
#returns various statistics.
function measskew(ref::DataMD, sig::DataMD;
	tstart_ref=-Inf, tstart_sig=-Inf)
	xrise = CrossType(:rise)
	xfall = CrossType(:fall)

	Δr = measdelay(ref, sig, tstart_ref=tstart_ref, tstart_main=tstart_sig,
		xing_ref=xrise, xing_main=xrise)
	Δf = measdelay(ref, sig, tstart_ref=tstart_ref, tstart_main=tstart_sig,
		xing_ref=xfall, xing_main=xfall)
	return _getskewstats(Δr, Δf)
end

#Last line


text:
module MetoceanTools
using PyPlot
using Optim
using StatsBase

export tsParm, tsBurst, waveSpec

type dirnSpec
 isGoingTo::Bool      # tells whether direction is going to (true) or coming from (false)
 isClockwise::Bool    # tells whether direction rotate c/w with increase value (true) or roates cc/w with increasing value (false)
 refDirn::Float64   # Gives some reference direction to which a directional parameter refers
 function dirnSpec()
   return new()
 end
end

"""
 Time series type for basic parameterised metocean data
"""
type tsParm
 sDate::Array{DateTime,1}       # Each row is another time corresponding to rows in tsData
 Data::Array{Float64,2}         # Columns are different parameters, rows are times
 iName::Dict                    # Lookup column index by Name iName["Hs"] returns column index with Hs
 Name::Array{AbstractString,1}  # Names of parameters in each column
 Meta::Dict                     # Arbitrary storage of meta-data
 Units::Array{AbstractString,1} # Units of parameters
 isDirn::Array{Bool,1}          # for Each column gives whether this parameter is directional (e.g. [0-360])
 isVector::Array{Bool,1}    	# for Each column gives whether this parameter is a complex vector
 dirnSpec::Array{dirnSpec,1}    # for Each column, If isDirn==true or isVector==true then gives direction specification
 function tsParm()
   return new()
 end
end
### Add function to set the iName automatically if empty when setting Name array

"""
 Time series type for burst record data  -  like surface elevations from a wave buoy (or daily market data)
 Unlike parameter data which can store multiple parameters in one record, tsBurst can store only one parameter
"""
type tsBurst
 sDate::Array{DateTime,1}
 sTime::Array{Float64,2}
 Data::Array{Float64,2}
 function tsBurst()
   return new()
 end
end

"""
 Basic wave spectrum records data which can contain 1D or 2D spectra
"""
type waveSpec
 f::Array{Float64,1}
 df::Array{Float64,1}
 th::Array{Float64,1}
 dth::Array{Float64,1}
 Sf::Array{Float64,1}
 Sth::Array{Float64,1}
 S::Array{Float64,2}
 function waveSpec()
   t=new()
   t.f=makeWaveSpecFreq()
   t.df=makeWaveSpecDelFreq(t.f)
   return t
 end
end

# dlp includes:
#  dlpRead for reading in *.dlp files
include("dlp.jl")

#  tsPlot includes plotting routines for plotting tsParm data series
#    tsQuickPlot - plots time series of all data in tsParm data series as stacked plot
include("tsPlot.jl")

# Given an array will make a count of unique values
include("countUniqueValues.jl")

#   findGaps for finding gaps in data and duration of good data
include("findGaps.jl")

#   showSampling for displaying and getting sampling within time series data
include("showSampling.jl")

# Provides basic zero-crossing analysis data types and methods
include("tsZeroCross.jl")

# General functions for matlab conversion
include("matlabHelpers.jl")

# General functions for Julia library
include("Helpers.jl")

# wave Spectrum functions
include("waveSpec.jl")

end # module

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
export AngularMeasurements, AngularSpectra, parent, angles, angleidx, closest_angle, between_indices
import Base.parent
# export maximumamp, minimumamp, normalizeamp, normalizeamp!
# export normalizeampangle, normalizeampangle!
# export normalizeampfreq, normalizeampfreq!
# export relativeto, relativeto!, smooth, interpolate, mirror, mirror!

struct AngularMeasurements{T,P<:AbstractSampleArray{T},A<:AbstractVector{<:Real}} <: AbstractSampleArray{T}
    parent::P
    angles::A

    function AngularMeasurements(parent, angles)
        nchannels(parent) == length(angles) || throw(ArgumentError("number of channels ($(nchannels(a))) != number of angles ($(length(angles)))"))
        s = deepcopy(parent) # because we sort it below
        a = normdeg.(Base.vect(angles...)) # normalize degrees to (-180°, 180°]
        length(a) == length(unique(a)) || throw(ArgumentError("non-unique angles!"))
        idxs = sortperm(a)
        new{eltype(parent), typeof(parent), typeof(angles)}(s[:, idxs], collect(a[idxs]))
    end
end

const AngularSpectra{T,P<:AbstractSpectrumArray{T},A<:AbstractVector{<:Real}} = AngularMeasurements{T,P,A}

function Interpolations.interpolate(am::AngularMeasurements, template::T) where {U, T <: AbstractSampleArray{U}}
    AngularMeasurements(Interpolations.interpolate(parent(am), template), angles(am))
end

Base.parent(am::AngularMeasurements) = am.parent
angles(am::AngularMeasurements) = am.angles

SampleArrays.rate(am::AngularMeasurements) = rate(parent(am))
SampleArrays.data(am::AngularMeasurements) = data(parent(am))
SampleArrays.data_no0(as::AngularSpectra) = data_no0(parent(as))

SampleArrays.domain(am::AngularMeasurements) = domain(parent(am))
SampleArrays.domain_no0(as::AngularSpectra) = domain_no0(parent(as))

SampleArrays.names(am::AngularMeasurements) = names(parent(am))
SampleArrays.names!(am::AngularMeasurements) = names!(parent(am))
SampleArrays.nchannels(am::AngularMeasurements) = nchannels(parent(am))
SampleArrays.nframes(am::AngularMeasurements) = nframes(parent(am))

SampleArrays.ntimeframes(as::AngularSpectra) = ntimeframes(parent(as))

function Base.show(io::IO, t::MIME"text/plain", am::AngularMeasurements)
    as = angles(am)
    println(io, "AngularMeasurements: $(minimum(as))°..$(maximum(as))°:")
    show(io, t, parent(am))
end

function angleidx(angles::AbstractVector{<:Real}, α; func=(==))
    # find first angle index (and the actual angle), for which func is true
    for (i, a) in enumerate(angles)
        func(a, α) && return (idx=i, α=a)
    end
    error("angle ($(α)) not found!")
end

angleidx(am::AngularMeasurements, α) = angleidx(angles(am), α)

function closest_angle(angles::AbstractVector{<:Real}, α)
    closeness = abs.(angles .- α)
    _, idx = findmin(closeness) # index of angle spectra closest to the angle α
    (idx=idx, α=angles[idx])
end

closest_angle(am::AngularMeasurements, α) = closest_angle(angles(am), α)

function between_indices(am::AngularMeasurements, α)
    α = normdeg(α)
    nc = nchannels(am)
    if nc == 1
        return 1, 1
    end
    as = [am.angles..., am.angles[1]]
    for i in 1:nc
        j = i % nc + 1
        d1, d2 = deltadeg(am.angles[i], α), deltadeg(α, am.angles[j])
        if d1 >= 0 && d2 >= 0
            if d1 ≈ 0
                return i, i
            elseif d2 ≈ 0
                return j, j
            end
            return i, j
        end
    end
end

function Base.getindex(am::AngularMeasurements, I::R, J::S) where {R <: SampleArrays.FrameIndex, S <: SampleArrays.ChannelIndex}
    idxs = SampleArrays.tochannelidx(parent(am), J)
    p = parent(am)[I, idxs]
    AngularMeasurements(p, angles(am)[idxs])
end
Base.getindex(am::AngularMeasurements, I::R) where {R <: SampleArrays.FrameIndex} = am[I, :]

Base.getindex(am::AngularMeasurements, α::AbstractFloat) = am[:, angleidx(am, α).idx]
function Base.getindex(am::AngularMeasurements, α::AbstractFloat, β::AbstractFloat)
    α <= β || error("α > β!")
    am[:, last(between_indices(am, α)):first(between_indices(am, β))]
end

function Base.similar(am::AngularMeasurements, t::Type, dims::Dims)
    AngularMeasurements(similar(parent(am), t, dims), angles(am))
end

# maximumamp(ds::DirectivitySpectra) = maximum(abs.(ds.spectra))
# minimumamp(ds::DirectivitySpectra) = minimum(abs.(ds.spectra))

# normalizeamp(ds::DirectivitySpectra) = normalizeamp!(deepcopy(ds))

# function normalizeamp!(ds::DirectivitySpectra)
#     maxamp = maximumamp(ds)
#     maxamp == 0 && return ds
#     ds.spectra ./= maxamp
#     return ds
# end

# normalizeampangle(ds::DirectivitySpectra) = normalizeampfreq!(deepcopy(ds))

# function normalizeampangle!(ds::DirectivitySpectra)
#     for α in 1:length(ds)
#         maxamp = maximum(abs.(ds.spectra[:, α]))
#         maxamp == 0 && continue
#         ds.spectra[:, α] ./= maxamp
#     end
#     return ds
# end

# normalizeampfreq(ds::DirectivitySpectra) = normalizeampfreq!(deepcopy(ds))

# function normalizeampfreq!(ds::DirectivitySpectra)
#     for fidx in 1:nfreqs(ds)
#         maxamp = maximum(abs.(data(spectra(ds))[fidx, :]))
#         maxamp == 0 && continue
#         data(spectra(ds))[fidx, :] ./= maxamp
#     end
#     return ds
# end

function Interpolations.interpolate(am::AngularMeasurements, α)
    i, j = between_indices(am, α)
    if i == j
        return am.parent[:, i]
    end
    β, γ = angles(am)[[i, j]]
    b, g = Float32(deltadeg(β, α)), Float32(deltadeg(α, γ)) # TODO better way to keep Float32 of parent(am)!
    @assert b >= 0 && g >= 0
    s = b + g
    @assert s > 0
    return mean([2g/s, 2b/s]' .* parent(am)[:, [i, j]])
end

# relativeto(ds::DirectivitySpectra, to) = relativeto!(deepcopy(ds), to)

# function relativeto!(ds::DirectivitySpectra, to::AbstractSpectrumArray)
#     @assert nchannels(to) == 1
#     to = deepcopy(to.data)
#     any(to .== 0.0) && @warn("relativeto!: zero in denominator!")
#     ds.spectra.data ./= to
#     return ds
# end

# relativeto!(ds::DirectivitySpectra, α::Real) = relativeto!(ds, ds[Float64(α)])

# """
#     mirror(ds::DirectivitySpectra)

# Mirror all spectra, i.e., append copies of spectra with negative angles (if not already in the collection).
# """
# function mirror(ds::DirectivitySpectra)
#     newspectra = []
#     newangles = deepcopy(ds.angles)
#     for i in 1:nchannels(ds)
#         α = normdeg(-newangles[i])
#         if α ∉ newangles
#             push!(newangles, α)
#             push!(newspectra, ds[i])
#         end
#     end
#     newspectra = hcat(ds.spectra, hcat(newspectra...))
#     idxs = sortperm(newangles)
#     newspectra .= newspectra[:, idxs]
#     newangles .= newangles[idxs]
#     return DirectivitySpectra(newspectra, newangles)
# end

# smooth(ds::DirectivitySpectra{T}, octave::Real = 1//3) where {T <: RFFTSpectrumArray} =
#     DirectivitySpectra(smooth(ds.spectra, octave), angles(ds))

# mean(ds::DirectivitySpectra) = mean(ds.spectra)


text:
struct PlanetaryConstants{T<:AbstractFloat} <: ConstantsCollection
    Ω::T  # Rotation rate of the planet [rad/s].
    f::T  # Nominal value for the Coriolis frequency [rad/s].
    g::T  # Standard acceleration due to gravity [m/s²].
end

function Earth()
    Ω = 7.2921150e-5
    f = 1e-4
    g = 9.80665
    PlanetaryConstants(Ω, f, g)
end

function EarthStationary()
    Ω = 0
    f = 0
    g = 9.80665
    PlanetaryConstants(Ω, f, g)
end


text:
function onedLinForm( Ff, test, wg )
#= """
  onedLinForm - routine to compute \int{ f*test }

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    F = onedLinForm( Ff, test, w_g )
  ```

  Arguments:
  - `Ff`: Function values at the Gauss points
  - `test`: Matrix of test functions evaluated at the Gauss points
            (dim: n_gauss, n_dof)
  - `w_g`: Column vector of Gauss weights

  Output arguments:
  - `F`: the element's contribution to a right-hand-side
""" =#

  F = test'*(wg.*Ff)

  return F

end


text:
#---------------------------------------------------------
# # [ImageGeom](@id ImageGeom)
#---------------------------------------------------------

# This page explains
# [`ImageGeom`](https://github.com/JeffFessler/MIRT.jl/blob/master/src/fbp/image_geom.jl)
# in
# [MIRT.jl](http://github.com/JeffFessler/MIRT.jl).

# ### Setup

# Packages needed here.

using MIRT
using MIRTjim: jim, prompt
using Plots; default(markerstrokecolor=:auto)

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# When performing tomographic image reconstruction,
# one must specify the geometry of the grid of image pixels.
# (In contrast, in image denoising and image deblurring problems,
# one works with the given discrete image
# and no physical coordinates are needed.)

# The key parameters of a grid of image pixels are
# * the size (dimensions) of the grid, e.g., `128 × 128`,
# * the spacing of the pixels, e.g., `1mm × 1mm`,
# * the offset of the pixels relative to the origin, e.g., `(0,0)`

# The data type `ImageGeom` describes such a geometry,
# for arbitrary dimensions (2D, 3D, etc.).

# There are several ways to construct this structure.
# The default is a `128 × 128` grid
# with pixel size ``\Delta_X = \Delta_Y = 1`` (unitless) and zero offset:

ig = ImageGeom()

# Here is a 3D example with non-cubic voxel size:

ig = ImageGeom( (512,512,128), (1,1,2), (0,0,0) )

# To avoid remembering the order of the arguments,
# named keyword pairs are also supported:

ig = ImageGeom( dims=(512,512,128), deltas=(1,1,2), offsets=(0,0,0) )

# ### Units

# The pixel dimensions `deltas` can (and should!) be values with units.

# Here is an example for a video (2D+time) with 12 frames per second:
using UnitfulRecipes
using Unitful: mm, s

ig = ImageGeom( dims=(640,480,1000), deltas=(1mm,1mm,(1//12)s) )

# ### Properties

# An ImageGeom object has *many* useful properties;
# there are too many to remember, so there is built-in help:

ig.help

# This small 2D example illustrates how some properties are used:

ig = ImageGeom( dims=(12,10), deltas=(1mm,1mm), offsets=(0.5,0.5) )

showgrid = (ig) -> # x,y grid locations of pixel centers
    scatter(ig.xg, ig.yg, label="", xlabel="x", ylabel="y",
        xlims = maximum(abs, ig.x) * 1.2 .* (-1,1),
        xticks=[ig.x[1], zero(eltype(ig.x)), ig.x[end]],
        ylims = maximum(abs, ig.y) * 1.2 .* (-1,1),
        yticks=[ig.y[1], zero(eltype(ig.y)), ig.y[end]],
        aspect_ratio=1, title="offsets $(ig.offsets)")
showgrid(ig)

#-
prompt();

# ### Offsets (unitless translation of grid)

# The default `offsets` are zeros,
# corresponding to symmetric sampling around origin:

ig = ImageGeom( dims=(12,10), deltas=(1mm,1mm) )
p = showgrid(ig)

#-
prompt();

# That default offset is natural for tomography
# when considering finite pixel size:

square = (x,y,Δ) -> plot!(p, label="", color=:black,
    x .+ Δ[1] * ([0,1,1,0,0] .- 0.5),
    y .+ Δ[2] * ([0,0,1,1,0] .- 0.5),
)
showgrid(ig)
square2 = (x,y) -> square(x, y, ig.deltas)
square2.(ig.xg,ig.yg)
plot!(p)

#-
prompt();

# In that default geometry, the center `(0,0)` of the image
# is at a corner of the middle 4 pixels (for even image sizes).
# That default is typical for tomographic imaging (e.g., CT, PET, SPECT).
# One must be careful when using operations like `imrotate` or `fft`.

# ### Mask

# In tomographic image reconstruction, patients are usually more "round"
# than "square" so often we only want to estimate the pixels inside some
# support `mask`: a Bool array indicating which pixels are to be estimated.
# (The rest are constrained to be zero.)
# The `ImageGeom` struct has an entry to store this `mask`.
# The default is `Trues(dims)` which is a "lazy" Bool `AbstractArray`
# from the `FillArrays` package that is conceptually similar to `trues(dims)`
# but requires `O(1)` storage.  So there is essentially no memory penalty
# to storing this entry in the `ImageGeom` for users who do not want
# to think about a `mask`.
# For users who do want a `mask`, fortunately Julia uses a special `BitArray`
# type to store Bool arrays, so the storage is 8× less than using bytes
# in most other languages.

# Often we use a "circle inscribed in the square" as a generic support mask,
# and one of the built-in properties can generate such a circular mask:

mask = ig.circ()
ig = ImageGeom(ig.dims, ig.deltas, ig.offsets, mask)
ig.plot(jim)

# Note that `jim` displays the axes with the units naturally;
# see [MIRTjim.jl](http://github.com/JeffFessler/MIRTjim.jl).

# There is an older interface `image_geom` similar to the function
# of the same name in
# [the Matlab version of MIRT](https://github.com/JeffFessler/mirt)
# provided for backward compatibility,
# but using `ImageGeom` is recommended for Julia work.

# ### AxisArrays

# There is a natural connection between `ImageGeom` and `AxisArrays`.
# Note the automatic labeling by
# [MIRTjim.jim](https://github.com/JeffFessler/MIRTjim.jl).

using AxisArrays
using Unitful: mm
ig = ImageGeom( dims=(60,48), deltas=(1.5mm,1mm) )
za = AxisArray( ig.circ() * 10/mm ; x=ig.x, y=ig.y)
jim(za, "AxisArray example")


text:
using KUparser, KUnet, HDF5, JLD, ArgParse, Compat, CUDArt
VERSION < v"0.4-" && eval(Expr(:using,:Dates))
require("fscore.jl")

type Fselect allfeats; bestfeats; bestscore; idxqueue; scores; Fselect()=new(); end

function main()
    args = parse_commandline()
    (data, ndeps) = initdata(args)
    pt = eval(parse(args["arctype"]))

    # We want a single copy of the state variables
    fs = Fselect()
    fs.allfeats = eval(parse("Flist."*args["allfeats"]))
    fs.idxqueue = randperm(length(fs.allfeats))
    fs.bestfeats = eval(parse("Flist."*args["feats"]))
    fs.bestscore = getcache(args["cache"], fs.bestfeats)
    if fs.bestscore < 0
        fs.bestscore = @fetchfrom workers()[1] fscore(pt, data, ndeps, fs.bestfeats, args)
        updatecache(args["cache"], fs.bestfeats, fs.bestscore)
    end

    function getnextidx()
        # See if we can improve fs.bestfeats with single feature flip steps from cache
        while true
            fs.scores = map(fs.allfeats) do f
                getcache(args["cache"], fs.bestfeats, f)
            end
            (smax,imax) = findmax(fs.scores)
            smax <= fs.bestscore && break
            fs.bestscore = smax
            fs.bestfeats = flip!(fs.bestfeats, fs.allfeats[imax])
            fs.idxqueue = randperm(length(fs.allfeats))
        end

        @show fs; flush(STDOUT)

        # Find the next neighbor to compute
        while (!isempty(fs.idxqueue) && (fs.scores[fs.idxqueue[1]] >= 0))
            shift!(fs.idxqueue)
        end

        return (minimum(fs.scores) >= 0 ? nothing : # terminate if all scores known
                isempty(fs.idxqueue) ? 0 :          # wait for remaining computations if queue empty
                shift!(fs.idxqueue))                # compute the next element in queue
    end

    # Feeder tasks based on multi.jl:pmap implementation:
    while !(isempty(fs.idxqueue) && (getnextidx()==nothing))
        @show fs; flush(STDOUT)
        @sync for wpid in workers()
            @async begin
                idx = getnextidx()
                # This got messed up because gc() is a leaky bucket
                # while idx != nothing
                if idx == nothing
                    # do nothing, time to terminate
                elseif idx == 0
                    sleep(10)
                else
                    @show ("$wpid gets fs.allfeats[$idx]=$(fs.allfeats[idx])"); flush(STDOUT)
                    @show feats = flip(fs.bestfeats, fs.allfeats[idx])
                    try
                        score = remotecall_fetch(wpid, fscore, pt, data, ndeps, feats, args)
                        @show ("$wpid gets score[$(fs.allfeats[idx])]=$score"); flush(STDOUT)
                        if isa(score, Number)
                            updatecache(args["cache"], feats, score)
                        else
                            push!(fs.idxqueue, idx)
                            @show ("Got $score from $wpid"); flush(STDOUT); sleep(10)
                        end
                    catch ex
                        push!(fs.idxqueue, idx)
                        @show ("Caught $ex from $wpid"); flush(STDOUT); sleep(10)
                    end
                end  # if idx == 0
                # idx = getnextidx()
                # end
            end  # @async begin
        end  # @sync for wpid in workers()
        @show fs; flush(STDOUT)
        @date Main.restartmachines(); flush(STDOUT)
        require("KUparser")
        require("fscore.jl")
    end  # while
    info("$(fs.bestscore)\t$(join(sort(fs.bestfeats), ' '))")
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "datafiles"
        help = "JLD files for input: first one will be used for training second for scoring"
        nargs = '+'
        required = true
        "--cache"
        help = "Cache file (plain text) for saving and reading results"
        default = ""
        "--nogpu"
        help = "Do not use gpu"
        action = :store_true
        # We use oparser on the training set to get training data.
        # We train the net for n epochs with fixed training data, no parsing.
        # Finally we evaluate by parsing the dev corpus using greedy or beam parser specified below.
        "--parser"
        help = "Parsing algorithm to use: g(reedy)parser, b(eam)parser"
        default = "gparser"
        "--arctype"
        help = "Move set to use: ArcEager{R1,13}, ArcHybrid{R1,13}"
        default = "ArcEager13"
        "--feats"
        help = "Features to start with from KUparser.Flist"
        default = "acl11eager"
        "--allfeats"
        help = "Complete feature set to use from KUparser.Flist"
        default = "acl11eager"
        "--hidden"
        help = "One or more hidden layer sizes"
        nargs = '+'
        arg_type = Int
        default = [20000]
        "--epochs"
        help = "Number of epochs to train"
        arg_type = Int
        default = 10
        # Don't need this, we only parse at the end:
        # "--pepochs"
        # help = "Number of epochs between parsing"
        # arg_type = Int
        # default = 1
        # We will just use whatever julia was started with
        # "--ncpu"
        # help = "Number of workers for multithreaded parsing"
        # arg_type = Int
        # default = 6
        "--pbatch"
        help = "Minibatch size for parsing"
        arg_type = Int
        default = 500
        "--tbatch"
        help = "Minibatch size for training"
        arg_type = Int
        default = 128
        "--nbeam"
        help = "Beam size for bparser"
        arg_type = Int
        default = 10
        "--adagrad"
        help = "If nonzero apply adagrad using arg as epsilon"
        arg_type = Float32
        nargs = '+'
        default = [1f-8]
        "--dropout"
        help = "Dropout probability"
        arg_type = Float32
        nargs = '+'
        # No dropout for 10-epoch training
        # default = [0.2f0, 0.7f0]
        "--l1reg"
        help = "L1 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--l2reg"
        help = "L2 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--learningRate"
        help = "Learning rate"
        arg_type = Float32
        nargs = '+'
        default = [2f-2]
        "--maxnorm"
        help = "If nonzero upper limit on weight matrix row norms"
        arg_type = Float32
        nargs = '+'
        "--momentum"
        help = "Momentum"
        arg_type = Float32
        nargs = '+'
        "--nesterov"
        help = "Apply nesterov's accelerated gradient"
        arg_type = Float32
        nargs = '+'
    end
    args = parse_args(s)
    for (arg,val) in args
        print("$arg:$val ")
    end
    println("")
    return args
end

function getcache(cachefile::String, feats::Fvec, f::Feature="")
    isempty(f) || (feats = flip(feats, f))
    fkey = join(sort(feats), ' ')
    fval = -1.0
    if isfile(cachefile)
        open(cachefile) do fp
            for l in eachline(fp)
                (score, fstr) = split(chomp(l), '\t')
                fstr == fkey && (fval = float(score); break)
            end
        end
    end
    return fval
end

function updatecache(cachefile::String, feats::Fvec, score::Real)
    open(cachefile, "a") do fp
        fkey = join(sort(feats), ' ')
        write(fp, "$score\t$fkey\n")
    end
end

function updatecache(cachefile::String, feats::Fvec, toflip::Fvec, scores::AbstractArray)
    @assert length(toflip) == length(scores)
    open(cachefile, "a") do fp
        for i=1:length(toflip)
            fkey = join(sort(flip(feats,toflip[i])), ' ')
            fval = scores[i]
            write(fp, "$fval\t$fkey\n")
        end
    end
end

function initdata(args::Dict)
    data = Corpus[]
    deprel = nothing
    ndeps = 0
    for f in args["datafiles"]
        @date d=load(f)
        push!(data, d["corpus"])
        if deprel == nothing
            deprel = d["deprel"]
            ndeps = length(deprel)
        else
            @assert deprel == d["deprel"]
        end
    end
    @assert length(data)==2 "fselect needs a training and a development set"
    (data, ndeps)
end

main()


text:
#Abstract type of a position in the space
abstract type Position end


text:
using UnsafeArrays
import Base: showarg, eltype
const DSYEVR_ = (BLAS.@blasfunc(dsyevr_),Base.liblapack_name)

# ----------------------------------------------------
# Zero cone
# ----------------------------------------------------
"""
    ZeroSet(dim)

Creates the zero set ``\\{ 0 \\}^{dim}`` of dimension `dim`. If `x` ∈ `ZeroSet` then all entries of x are zero.
"""
struct ZeroSet{T} <: AbstractConvexCone{T}
	dim::Int
	function ZeroSet{T}(dim::Int) where {T}
		dim >= 0 ? new(dim) : throw(DomainError(dim, "dimension must be nonnegative"))
	end
end
ZeroSet(dim) = ZeroSet{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, ::ZeroSet{T}) where{T}
	x .= zero(T)
	return nothing
end

function in_dual(x::AbstractVector{T}, ::ZeroSet{T}, tol::T) where{T}
	true
end

function in_pol_recc(x::AbstractVector{T}, ::ZeroSet{T}, tol::T) where{T}
	!any( x-> (abs(x) > tol), x)
end

function scale!(::ZeroSet{T}, ::SplitView{T}) where{T}
	return nothing
end

function rectify_scaling!(E,work,set::ZeroSet{T}) where{T}
	return false
end

# ----------------------------------------------------
# Nonnegative orthant
# ----------------------------------------------------
"""
    Nonnegatives(dim)

Creates the nonnegative orthant ``\\{ x \\in \\mathbb{R}^{dim} : x \\ge 0 \\}``  of dimension `dim`.
"""
struct Nonnegatives{T} <: AbstractConvexCone{T}
	dim::Int
	function Nonnegatives{T}(dim::Int) where {T}
		dim >= 0 ? new(dim) : throw(DomainError(dim, "dimension must be nonnegative"))
	end
end
Nonnegatives(dim) = Nonnegatives{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, C::Nonnegatives{T}) where{T}
	x .= max.(x, zero(T))
	return nothing
end

function in_dual(x::AbstractVector{T}, ::Nonnegatives{T}, tol::T) where{T}
	!any( x-> (x < -tol), x)
end

function in_pol_recc(x::AbstractVector{T}, ::Nonnegatives{T}, tol::T) where{T}
	!any( x-> (x > tol), x)
end

function scale!(cone::Nonnegatives{T}, ::AbstractVector{T}) where{T}
	return nothing
end

function rectify_scaling!(E, work, set::Nonnegatives{T}) where{T}
	return false
end

# ----------------------------------------------------
# Second Order Cone
# ----------------------------------------------------
"""
    SecondOrderCone(dim)

Creates the second-order cone (or Lorenz cone) ``\\{ (t,x) \\in \\mathrm{R}^{dim} : || x ||_2  \\leq t \\}``.
"""
struct SecondOrderCone{T} <: AbstractConvexCone{T}
	dim::Int
	function SecondOrderCone{T}(dim::Int) where {T}
		dim >= 0 ? new(dim) : throw(DomainError(dim, "dimension must be nonnegative"))
	end
end
SecondOrderCone(dim) = SecondOrderCone{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, ::SecondOrderCone{T}) where{T}
	t = x[1]
	xt = view(x, 2:length(x))
	norm_x = norm(xt, 2)
	if norm_x <= t
		nothing
	elseif norm_x <= -t
		x[:] .= zero(T)
	else
		x[1] = (norm_x + t) / 2
		#x(2:end) assigned via view
		@. xt = (norm_x + t) / (2 * norm_x) * xt
	end
	return nothing
end

function in_dual(x::AbstractVector{T}, ::SecondOrderCone{T}, tol::T) where{T}
	@views norm(x[2:end]) <= (tol + x[1]) #self dual
end

function in_pol_recc(x::AbstractVector{T}, ::SecondOrderCone, tol::T) where{T}
	@views norm(x[2:end]) <= (tol - x[1]) #self dual
end

function scale!(cone::SecondOrderCone{T}, ::AbstractVector{T}) where{T}
	return nothing
end

function rectify_scaling!(E, work, set::SecondOrderCone{T}) where{T}
	return rectify_scalar_scaling!(E, work)
end

# ----------------------------------------------------
# Positive Semidefinite Cone
# ----------------------------------------------------

#a type to maintain internal workspace data for the BLAS syevr function

mutable struct PsdBlasWorkspace{T}
    m::Base.RefValue{Int64}
    w::Vector{T}
    Z::Matrix{T}
    isuppz::Vector{BLAS.BlasInt}
    work::Vector{T}
    lwork::BLAS.BlasInt
    iwork::Vector{BLAS.BlasInt}
    liwork::BLAS.BlasInt
    info::Base.RefValue{Int64}

    function PsdBlasWorkspace{T}(n::Integer) where{T}

        BlasInt = BLAS.BlasInt

        #workspace data for BLAS
        m      = Ref{BlasInt}()
        w      = Vector{T}(undef,n)
        Z      = Matrix{T}(undef,n,n)
        isuppz = Vector{BlasInt}(undef, 2*n)
        work   = Vector{T}(undef, 1)
        lwork  = BlasInt(-1)
        iwork  = Vector{BlasInt}(undef, 1)
        liwork = BlasInt(-1)
        info   = Ref{BlasInt}()

        new(m,w,Z,isuppz,work,lwork,iwork,liwork,info)
    end
end

function _syevr!(A::AbstractMatrix{T}, ws::PsdBlasWorkspace) where{T <: Float64}

    #Float64 only support for now since we call dsyevr_ directly
    n       = size(A,1)
    ldz     = n
    lda     = stride(A,2)

        ccall(DSYEVR_, Cvoid,
        (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BLAS.BlasInt},
        Ptr{T}, Ref{BLAS.BlasInt}, Ref{T}, Ref{T},
        Ref{BLAS.BlasInt}, Ref{BLAS.BlasInt}, Ref{T}, Ptr{BLAS.BlasInt},
        Ptr{T}, Ptr{T}, Ref{BLAS.BlasInt}, Ptr{BLAS.BlasInt},
        Ptr{T}, Ref{BLAS.BlasInt}, Ptr{BLAS.BlasInt}, Ref{BLAS.BlasInt},
        Ptr{BLAS.BlasInt}),
        'V', 'A', 'U', n,
        A, max(1,lda), 0.0, 0.0,
        0, 0, -1.0,
        ws.m, ws.w, ws.Z, ldz, ws.isuppz,
        ws.work, ws.lwork, ws.iwork, ws.liwork,
        ws.info)
        LAPACK.chklapackerror(ws.info[])
end

function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}

    #computes the upper triangular part of the projection of X onto the PSD cone

     #allocate additional workspace arrays if the ws
     #work and iwork have not yet been sized
     if ws.lwork == -1
         _syevr!(X,ws)
         ws.lwork = BLAS.BlasInt(real(ws.work[1]))
         resize!(ws.work, ws.lwork)
         ws.liwork = ws.iwork[1]
         resize!(ws.iwork, ws.liwork)
     end

	 # below LAPACK function does the following: w,Z  = eigen!(Symmetric(X))
	 	_syevr!(X, ws)
		# compute upper triangle of: X .= Z*Diagonal(max.(w, 0.0))*Z'
		rank_k_update!(X, ws)
end

function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where{T}
  n = size(X, 1)
  X .= 0
  nnz_λ = 0
  for j = 1:length(ws.w)
    λ = ws.w[j]
    if λ > 0
      nnz_λ += 1
      @inbounds for i = 1:n
        ws.Z[i, j] = ws.Z[i, j] * sqrt(λ)
      end
    end
  end

  if nnz_λ > 0
    V = uview(ws.Z, :, (n - nnz_λ + 1):n)
    BLAS.syrk!('U', 'N', 1.0, V, 1.0, X)
  end
  return nothing
end

"""
    PsdCone(dim)

Creates the cone of symmetric positive semidefinite matrices ``\\mathcal{S}_+^{dim}``. The entries of the matrix `X` are stored column-by-column in the vector `x` of dimension `dim`.
Accordingly  ``X \\in \\mathbb{S}_+ \\Rightarrow x \\in \\mathcal{S}_+^{dim}``, where ``X = \\text{mat}(x)``.
"""
struct PsdCone{T} <: AbstractConvexCone{T}
	dim::Int
	sqrt_dim::Int
    work::PsdBlasWorkspace{T}
	function PsdCone{T}(dim::Int) where{T}
		dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
		iroot = isqrt(dim)
		iroot^2 == dim || throw(DomainError(dim, "dimension must be a square"))
		new(dim, iroot,PsdBlasWorkspace{T}(iroot))
	end
end
PsdCone(dim) = PsdCone{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, cone::PsdCone{T}) where{T}
	n = cone.sqrt_dim

    # handle 1D case
    if length(x) == 1
        x = max.(x, zero(T))
    else
        # symmetrized square view of x
        X    = reshape(x, n, n)
        symmetrize_upper!(X)
        _project!(X,cone.work)

        #fill in the lower triangular part
        for j=1:n, i=1:(j-1)
            X[j,i] = X[i,j]
        end
    end
    return nothing
end

function in_dual(x::AbstractVector{T}, cone::PsdCone{T}, tol::T) where{T}
	n = cone.sqrt_dim
	X = reshape(x, n, n)
  return is_pos_sem_def(X, tol)
end

function in_pol_recc(x::AbstractVector{T}, cone::PsdCone{T}, tol::T) where{T}
	n = cone.sqrt_dim
	X = reshape(x, n, n)
	return is_neg_sem_def(X, tol)
end

function scale!(cone::PsdCone{T}, ::AbstractVector{T}) where{T}
	return nothing
end

function rectify_scaling!(E, work, set::PsdCone{T}) where{T}
	return rectify_scalar_scaling!(E, work)
end

# ----------------------------------------------------
# Positive Semidefinite Cone (Triangle)
# ----------------------------------------------------

# Psd cone given by upper-triangular entries of matrix
"""
    PsdConeTriangle(dim)

Creates the cone of symmetric positive semidefinite matrices. The entries of the upper-triangular part of matrix `X` are stored in the vector `x` of dimension `dim`.
A ``r \\times r`` matrix has ``r(r+1)/2`` upper triangular elements and results in a vector of ``\\mathrm{dim} = r(r+1)/2``.

### Examples
The matrix
```math
\\begin{bmatrix} x_1 & x_2 & x_4\\\\ x_2 & x_3 & x_5\\\\ x_4 & x_5 & x_6 \\end{bmatrix}
```
is transformed to the vector ``[x_1, x_2, x_3, x_4, x_5, x_6]^\\top `` with corresponding constraint  `PsdConeTriangle(6)`.

"""
struct PsdConeTriangle{T} <: AbstractConvexCone{T}

    dim::Int #dimension of vector
    sqrt_dim::Int # side length of matrix
    X::Array{T,2}
    work::PsdBlasWorkspace{T}

    function PsdConeTriangle{T}(dim::Int) where{T}
        dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
        side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);

        new(dim, side_dimension, zeros(side_dimension, side_dimension),PsdBlasWorkspace{T}(side_dimension))
    end
end
PsdConeTriangle(dim) = PsdConeTriangle{DefaultFloat}(dim)

function project!(x::AbstractArray, cone::PsdConeTriangle{T}) where{T}
    # handle 1D case
    if length(x) == 1
        x = max.(x,zero(T))
    else
        populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
        _project!(cone.X,cone.work)
        extract_upper_triangle!(cone.X, x, sqrt(2) )
    end
    return nothing
end

function in_dual(x::AbstractVector{T}, cone::PsdConeTriangle{T}, tol::T) where{T}
    n = cone.sqrt_dim
    populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
    return is_pos_sem_def(cone.X, tol)
end

function in_pol_recc(x::AbstractVector{T}, cone::PsdConeTriangle{T}, tol::T) where{T}
    n = cone.sqrt_dim
    populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
    Xs = Symmetric(cone.X)
    return is_neg_sem_def(cone.X, tol)
end

function scale!(cone::PsdConeTriangle{T}, ::AbstractVector{T}) where{T}
    return nothing
end

function rectify_scaling!(E, work, set::PsdConeTriangle{T}) where{T}
    return rectify_scalar_scaling!(E,work)
end

function populate_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
 	k = 0
  	for j in 1:size(A, 2)
     	for i in 1:j
        	k += 1
        	if i != j
        		A[i, j] = scaling_factor * x[k]
        	else
        		A[i, j] = x[k]
        	end
      	end
  	end
  	nothing
end

function extract_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
	k = 0
  	for j in 1:size(A, 2)
     	for i in 1:j
        	k += 1
        	if i != j
        		x[k] = scaling_factor * A[i, j]
        	else
        		x[k] = A[i, j]
        	end
      	end
  	end
	nothing
end

# ----------------------------------------------------
# Box
# ----------------------------------------------------
"""
    Box(l, u)

Creates a box or intervall with lower boundary vector ``l \\in  \\mathbb{R}^m \\cup \\{-\\infty\\}^m`` and upper boundary vector``u \\in \\mathbb{R}^m\\cup \\{+\\infty\\}^m``.
"""
struct Box{T} <: AbstractConvexSet{T}
	dim::Int
	l::Vector{T}
	u::Vector{T}
	function Box{T}(dim::Int) where{T}
		dim >= 0 || throw(DomainError(dim, "dimension must be nonnegative"))
		l = fill!(Vector{T}(undef, dim), -Inf)
		u = fill!(Vector{T}(undef, dim), +Inf)
		new(dim, l, u)
	end
	function Box{T}(l::Vector{T}, u::Vector{T}) where{T}
		length(l) == length(u) || throw(DimensionMismatch("bounds must be same length"))
        _box_check_bounds(l,u)
        #enforce consistent bounds
		new(length(l), l, u)
	end
end
Box(dim) = Box{DefaultFloat}(dim)
Box(l, u) = Box{DefaultFloat}(l, u)

function _box_check_bounds(l,u)
    for i in eachindex(l)
        l[i] > u[i] && error("Box set: inconsistent lower/upper bounds specified at index i = ", i, ": l[i] = ",l[i],", u[i] = ",u[i])
    end
end

function project!(x::AbstractVector{T}, box::Box{T}) where{T}
	@. x = clip(x, box.l, box.u)
	return nothing
end

function support_function(x::AbstractVector{T}, B::Box{T}, tol::T) where{T}
    s = 0.
    for i in eachindex(x)
        s+= ( abs(x[i] > tol) && x[i] > 0) ? x[i]*B.u[i] : x[i]*B.l[i]
    end
    return s
end

function in_pol_recc(x::AbstractVector{T}, B::Box{T}, tol::T) where{T}
    !any(XU -> (XU[2] == Inf && XU[1] > tol), zip(x,B.u)) && !any(XL -> (XL[2] == -Inf && XL[1] < -tol), zip(x,B.l))
end

function scale!(box::Box{T}, e::AbstractVector{T}) where{T}
	@. box.l = box.l * e
	@. box.u = box.u * e
	return nothing
end

function rectify_scaling!(E, work, box::Box{T}) where{T}
	return false #no correction needed since we can scale componentwise
end

# ----------------------------------------------------
# Composite Set
# ----------------------------------------------------

#struct definition is provided in projections.jl, since it
#must be available to SplitVector, which in turn must be
#available for most of the methods here.

CompositeConvexSet(args...) = CompositeConvexSet{DefaultFloat}(args...)

function project!(x::SplitVector{T}, C::CompositeConvexSet{T}) where{T}
	@assert x.split_by === C
	foreach(xC -> project!(xC[1], xC[2]), zip(x.views, C.sets))
	return nothing
end

function support_function(x::SplitVector{T}, C::CompositeConvexSet{T}, tol::T) where{T}
	sum(xC -> support_function(xC[1], xC[2], tol),zip(x.views, C.sets))
end

function in_pol_recc(x::SplitVector{T}, C::CompositeConvexSet{T}, tol::T) where{T}
	all(xC -> in_pol_recc(xC[1], xC[2], tol), zip(x.views, C.sets))
end

function scale!(C::CompositeConvexSet{T}, e::SplitVector{T}) where{T}
	@assert e.split_by === C
	for i = eachindex(C.sets)
		scale!(C.sets[i], e.views[i])
	end
end

function rectify_scaling!(E::SplitVector{T},
	work::SplitVector{T},
	C::CompositeConvexSet{T}) where {T}
	@assert E.split_by === C
	@assert work.split_by === C
	any_changed = false
	for i = eachindex(C.sets)
		any_changed |= rectify_scaling!(E.views[i], work.views[i], C.sets[i])
	end
	return any_changed
end

#-------------------------
# general AbstractConvexCone operations
#-------------------------

# sup_{z in K_tilde_b = {-K} x {b} } <z,δy> = { <y,b> ,if y in Ktilde_polar
#                                                 +∞   ,else}

function support_function(y::SplitView{T}, cone::AbstractConvexCone{T}, tol::T) where{T}
	in_dual(-y, cone, tol) ? 0. : Inf;
end

#-------------------------
# generic set operations
#-------------------------
# function Base.showarg(io::IO, C::AbstractConvexSet{T}, toplevel) where{T}
#    print(io, typeof(C), " in dimension '", A.dim, "'")
# end

eltype(::AbstractConvexSet{T}) where{T} = T
num_subsets(C::AbstractConvexSet{T}) where{T}  = 1
num_subsets(C::CompositeConvexSet{T}) where{T} = length(C.sets)

function get_subset(C::AbstractConvexSet, idx::Int)
	idx == 1 || throw(DimensionMismatch("Input only has 1 subset (itself)"))
	return C
end
get_subset(C::CompositeConvexSet, idx::Int) = C.sets[idx]

function rectify_scalar_scaling!(E, work)
	tmp = mean(E)
	work .= tmp ./ E
	return true
end

# computes the row indices of A,b for each convex set
function get_set_indices(sets::Array{COSMO.AbstractConvexSet, 1})
	sidx = 0
	indices = Array{UnitRange{Int64}, 1}(undef, length(sets))
	for i = eachindex(sets)
		indices[i] = (sidx + 1) : (sidx + sets[i].dim)
		sidx += sets[i].dim
	end
	return indices
end


text:
# Autogenerated wrapper script for Raylib_jll for x86_64-linux-gnu
export libraylib

using Libglvnd_jll
using Xorg_libXrandr_jll
using Xorg_libX11_jll
using Xorg_libXrender_jll
using Xorg_libXi_jll
using Xorg_libXext_jll
using Xorg_libXcursor_jll
using Xorg_libXdamage_jll
using Xorg_libXfixes_jll
using Xorg_libXcomposite_jll
using Xorg_libXinerama_jll
using GLU_jll
using Mesa_jll
using alsa_jll
JLLWrappers.@generate_wrapper_header("Raylib")
JLLWrappers.@declare_library_product(libraylib, "libraylib.so.351")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, Xorg_libXrandr_jll, Xorg_libX11_jll, Xorg_libXrender_jll, Xorg_libXi_jll, Xorg_libXext_jll, Xorg_libXcursor_jll, Xorg_libXdamage_jll, Xorg_libXfixes_jll, Xorg_libXcomposite_jll, Xorg_libXinerama_jll, GLU_jll, Mesa_jll, alsa_jll)
    JLLWrappers.@init_library_product(
        libraylib,
        "lib/libraylib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const GENERAL_FOLDER = "tierra_logs/tierra_logs" # geb_logs"

const SNAPSHOTS_FOLDER = GENERAL_FOLDER * "/snapshots/"
const LOGGER_FOLDER = GENERAL_FOLDER * "/logs/"

const CALCULATED_FOLDER = GENERAL_FOLDER * "/calculated/"

const LOG_FREQUENCY = 1000 # 1e7

const NUM_THREADS = 4

const EPS = 1e-16

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark