text:
using LinearAlgebra # For rmul!

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const AFFEQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}

@enum VariableType FREE NNEG QUAD PSD

struct VariableInfo
    variable_type::VariableType
    cone_index::Int
    index_in_cone::Int
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    b::Vector{Float64}

    variable_info::Vector{VariableInfo}

    # PSDCone variables
    psdc_dims::Vector{Int}
    psdc_Cvar::Vector{Vector{Int}}
    psdc_Cval::Vector{Vector{Float64}}
    psdc_Avar::Vector{Vector{Int}}
    psdc_Acon::Vector{Vector{Int}}
    psdc_Aval::Vector{Vector{Float64}}

    # QUADratic/second-order cone variables
    quad_dims::Vector{Int}
    quad_Cvar::Vector{Vector{Int}}
    quad_Cval::Vector{Vector{Float64}}
    quad_Avar::Vector{Vector{Int}}
    quad_Acon::Vector{Vector{Int}}
    quad_Aval::Vector{Vector{Float64}}

    # NonNEGatives variables
    num_nneg::Int
    nneg_info::Vector{Int} # Similar to `info` field of `MOI.Bridges.Variable.Map`.
    nneg_Cvar::Vector{Int}
    nneg_Cval::Vector{Float64}
    nneg_Avar::Vector{Int}
    nneg_Acon::Vector{Int}
    nneg_Aval::Vector{Float64}

    # FREE variables
    num_free::Int
    free_Cvar::Vector{Int}
    free_Cval::Vector{Float64}
    free_Avar::Vector{Int}
    free_Acon::Vector{Int}
    free_Aval::Vector{Float64}

    objective_sense::MOI.OptimizationSense
    objective_constant::Float64

    primal_objective_value::Float64
    dual_objective_value::Float64
    free_X::Vector{Float64}
    nneg_X::Vector{Float64}
    quad_X::Vector{Vector{Float64}}
    psdc_X::Vector{Vector{Float64}}
    y::Vector{Float64}
    free_Z::Vector{Float64}
    nneg_Z::Vector{Float64}
    quad_Z::Vector{Vector{Float64}}
    psdc_Z::Vector{Vector{Float64}}
    info::Dict{String, Any}
    runhist::Dict{String, Any}
    status::Union{Nothing, Int}
    solve_time::Float64

    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
        optimizer = new(
            Float64[], VariableInfo[],
            Int[], Vector{Int}[], Vector{Float64}[], Vector{Int}[], Vector{Int}[], Vector{Float64}[],
            Int[], Vector{Int}[], Vector{Float64}[], Vector{Int}[], Vector{Int}[], Vector{Float64}[],
            0, Int[], Int[], Float64[], Int[], Int[], Float64[],
            0, Int[], Float64[], Int[], Int[], Float64[],
            MOI.FEASIBILITY_SENSE, 0.0,
            NaN, NaN,
            Float64[], Float64[], Vector{Float64}[], Vector{Float64}[], # X
            Float64[], # y
            Float64[], Float64[], Vector{Float64}[], Vector{Float64}[], # Z
            Dict{String, Any}(), Dict{String, Any}(),
            nothing, NaN,
            false, Dict{Symbol, Any}())
        if !isempty(kwargs)
            @warn("""Passing optimizer attributes as keyword arguments to
            SDPT3.Optimizer is deprecated. Use
                MOI.set(model, MOI.RawOptimizerAttribute("key"), value)
            or
                JuMP.set_optimizer_attribute(model, "key", value)
            instead.
            """)
        end
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawOptimizerAttribute(string(key)), value)
        end
        return optimizer
    end
end

function MOI.supports(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    return param.name in ALLOWED_OPTIONS
end
function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    if !MOI.supports(optimizer, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    optimizer.options[Symbol(param.name)] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    # TODO: This gives a poor error message if the name of the parameter is invalid.
    return optimizer.options[Symbol(param.name)]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.get(::Optimizer, ::MOI.SolverName) = "SDPT3"
function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec)
    return optimizer.solve_time
end

function MOI.is_empty(optimizer::Optimizer)
    return isempty(optimizer.b) &&
        isempty(optimizer.variable_info) &&
        iszero(optimizer.psdc_dims) &&
        iszero(optimizer.quad_dims) &&
        iszero(optimizer.num_nneg) &&
        iszero(optimizer.num_free) &&
        optimizer.objective_sense == MOI.FEASIBILITY_SENSE &&
        iszero(optimizer.objective_constant)
end

function MOI.empty!(optimizer::Optimizer)
    empty!(optimizer.b)

    empty!(optimizer.variable_info)

    empty!(optimizer.psdc_dims)
    empty!(optimizer.psdc_Cvar)
    empty!(optimizer.psdc_Cval)
    empty!(optimizer.psdc_Avar)
    empty!(optimizer.psdc_Acon)
    empty!(optimizer.psdc_Aval)

    empty!(optimizer.quad_dims)
    empty!(optimizer.quad_Cvar)
    empty!(optimizer.quad_Cval)
    empty!(optimizer.quad_Avar)
    empty!(optimizer.quad_Acon)
    empty!(optimizer.quad_Aval)

    optimizer.num_nneg = 0
    empty!(optimizer.nneg_info)
    empty!(optimizer.nneg_Cvar)
    empty!(optimizer.nneg_Cval)
    empty!(optimizer.nneg_Avar)
    empty!(optimizer.nneg_Acon)
    empty!(optimizer.nneg_Aval)

    optimizer.num_free = 0
    empty!(optimizer.free_Cvar)
    empty!(optimizer.free_Cval)
    empty!(optimizer.free_Avar)
    empty!(optimizer.free_Acon)
    empty!(optimizer.free_Aval)

    optimizer.objective_sense = MOI.FEASIBILITY_SENSE
    optimizer.objective_constant = 0.0

    optimizer.primal_objective_value = NaN
    optimizer.dual_objective_value = NaN
    empty!(optimizer.free_X)
    empty!(optimizer.nneg_X)
    empty!(optimizer.quad_X)
    empty!(optimizer.psdc_X)
    empty!(optimizer.y)
    empty!(optimizer.free_Z)
    empty!(optimizer.nneg_Z)
    empty!(optimizer.quad_Z)
    empty!(optimizer.psdc_Z)
    empty!(optimizer.info)
    empty!(optimizer.runhist)
    optimizer.status = nothing
    optimizer.solve_time = NaN
end

function MOI.supports(
    optimizer::Optimizer,
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

const SupportedSets = Union{MOI.Nonnegatives, MOI.SecondOrderCone, MOI.PositiveSemidefiniteConeTriangle}
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables},
    ::Type{<:SupportedSets})
    return true
end
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{Cdouble}},
    ::Type{MOI.EqualTo{Cdouble}})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    return MOIU.default_copy_to(dest, src)
end
MOI.supports_incremental_interface(::Optimizer) = true

# Variables
function MOI.add_variable(optimizer::Optimizer)
    optimizer.num_free += 1
    push!(optimizer.variable_info, VariableInfo(FREE, 1, optimizer.num_free))
    return MOI.VariableIndex(length(optimizer.variable_info))
end
function MOI.add_variables(optimizer::Optimizer, n::Integer)
    return [MOI.add_variable(optimizer) for i in 1:n]
end
function _add_nonneg_variable(optimizer::Optimizer)
    optimizer.num_nneg += 1
    push!(optimizer.variable_info, VariableInfo(NNEG, 1, optimizer.num_nneg))
    return MOI.VariableIndex(length(optimizer.variable_info))
end
function _add_quad_variable(optimizer::Optimizer, index_in_cone)
    push!(optimizer.variable_info, VariableInfo(QUAD, length(optimizer.quad_dims), index_in_cone))
    return MOI.VariableIndex(length(optimizer.variable_info))
end
function _add_psdc_variable(optimizer::Optimizer, index_in_cone)
    push!(optimizer.variable_info, VariableInfo(PSD, length(optimizer.psdc_dims), index_in_cone))
    return MOI.VariableIndex(length(optimizer.variable_info))
end
function _add_constrained_variables(optimizer::Optimizer, set::MOI.Nonnegatives)
    push!(optimizer.nneg_info, MOI.dimension(set))
    for i in 2:MOI.dimension(set)
        push!(optimizer.nneg_info, i)
    end
    return [_add_nonneg_variable(optimizer) for i in 1:MOI.dimension(set)]
end
function _add_constrained_variables(optimizer::Optimizer, set::MOI.SecondOrderCone)
    push!(optimizer.quad_dims, MOI.dimension(set))
    push!(optimizer.quad_Cvar, Int[])
    push!(optimizer.quad_Cval, Float64[])
    push!(optimizer.quad_Avar, Int[])
    push!(optimizer.quad_Acon, Int[])
    push!(optimizer.quad_Aval, Float64[])
    return [_add_quad_variable(optimizer, i) for i in 1:MOI.dimension(set)]
end
function _add_constrained_variables(optimizer::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
    push!(optimizer.psdc_dims, MOI.side_dimension(set))
    push!(optimizer.psdc_Cvar, Int[])
    push!(optimizer.psdc_Cval, Float64[])
    push!(optimizer.psdc_Avar, Int[])
    push!(optimizer.psdc_Acon, Int[])
    push!(optimizer.psdc_Aval, Float64[])
    return [_add_psdc_variable(optimizer, i) for i in 1:MOI.dimension(set)]
end
function MOI.add_constrained_variables(optimizer::Optimizer, set::SupportedSets)
    vis = _add_constrained_variables(optimizer, set)
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(first(vis).value)
    return vis, ci
end

# Objective
function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveSense)
    return optimizer.objective_sense
end
sense_to_sign(sense::MOI.OptimizationSense) = sense == MOI.MAX_SENSE ? -1 : 1
function MOI.set(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense != optimizer.objective_sense
        sign = sense_to_sign(sense)
        rmul!(optimizer.free_Cval, -1)
        rmul!(optimizer.nneg_Cval, -1)
        for i in eachindex(optimizer.quad_dims)
            rmul!(optimizer.quad_Cval[i], -1)
        end
        for i in eachindex(optimizer.psdc_dims)
            rmul!(optimizer.psdc_Cval[i], -1)
        end
    end
    optimizer.objective_sense = sense
end
function MOI.set(optimizer::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
                 func::MOI.ScalarAffineFunction{Float64})
    optimizer.objective_constant = MOI.constant(func)
    empty!(optimizer.free_Cvar)
    empty!(optimizer.free_Cval)
    empty!(optimizer.nneg_Cvar)
    empty!(optimizer.nneg_Cval)
    for i in eachindex(optimizer.quad_dims)
        empty!(optimizer.quad_Cvar[i])
        empty!(optimizer.quad_Cval[i])
    end
    for i in eachindex(optimizer.psdc_dims)
        empty!(optimizer.psdc_Cvar[i])
        empty!(optimizer.psdc_Cval[i])
    end
    sign = sense_to_sign(optimizer.objective_sense)
    for term in func.terms
        info = optimizer.variable_info[term.variable.value]
        if info.variable_type == FREE
            push!(optimizer.free_Cvar, info.index_in_cone)
            push!(optimizer.free_Cval, sign * term.coefficient)
        elseif info.variable_type == NNEG
            push!(optimizer.nneg_Cvar, info.index_in_cone)
            push!(optimizer.nneg_Cval, sign * term.coefficient)
        elseif info.variable_type == QUAD
            push!(optimizer.quad_Cvar[info.cone_index], info.index_in_cone)
            push!(optimizer.quad_Cval[info.cone_index], sign * term.coefficient)
        else
            @assert info.variable_type == PSD
            push!(optimizer.psdc_Cvar[info.cone_index], info.index_in_cone)
            push!(optimizer.psdc_Cval[info.cone_index], sign * term.coefficient)
        end
    end
end

# Constraints
function is_diagonal_index(k)
    # See https://jump.dev/MathOptInterface.jl/v0.9.3/apireference/#MathOptInterface.AbstractSymmetricMatrixSetTriangle
    i = div(1 + isqrt(8k - 7), 2)
    j = k - div((i - 1) * i, 2)
    return i == j
end
function MOI.add_constraint(optimizer::Optimizer, func::MOI.ScalarAffineFunction{Float64}, set::MOI.EqualTo{Float64})
    if !iszero(MOI.constant(func))
        # We use the fact that the initial function constant was zero to
        # implement getters for `MOI.ConstraintPrimal`.
        throw(MOI.ScalarFunctionConstantNotZero{
             Float64, typeof(func), typeof(set)}(MOI.constant(func)))
    end
    push!(optimizer.b, MOI.constant(set))
    con = length(optimizer.b)
    for term in func.terms
        info = optimizer.variable_info[term.variable.value]
        if info.variable_type == FREE
            push!(optimizer.free_Avar, info.index_in_cone)
            push!(optimizer.free_Acon, con)
            push!(optimizer.free_Aval, term.coefficient)
        elseif info.variable_type == NNEG
            push!(optimizer.nneg_Avar, info.index_in_cone)
            push!(optimizer.nneg_Acon, con)
            push!(optimizer.nneg_Aval, term.coefficient)
        elseif info.variable_type == QUAD
            push!(optimizer.quad_Avar[info.cone_index], info.index_in_cone)
            push!(optimizer.quad_Acon[info.cone_index], con)
            push!(optimizer.quad_Aval[info.cone_index], term.coefficient)
        else
            @assert info.variable_type == PSD
            push!(optimizer.psdc_Avar[info.cone_index], info.index_in_cone)
            push!(optimizer.psdc_Acon[info.cone_index], con)
            coef = is_diagonal_index(info.index_in_cone) ? term.coefficient : term.coefficient / √2
            push!(optimizer.psdc_Aval[info.cone_index], coef)
        end
    end
    return AFFEQ(con)
end

# TODO could do something more efficient here
#      `SparseMatrixCSC` is returned in SumOfSquares.jl test `sos_horn`
symvec(Q::SparseMatrixCSC) = symvec(Matrix(Q))
function symvec(Q::Matrix)
    n = LinearAlgebra.checksquare(Q)
    vec_dim = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
    q = Vector{eltype(Q)}(undef, vec_dim)
    k = 0
    for j in 1:n
        for i in 1:j
            k += 1
            q[k] = Q[i, j]
        end
    end
    @assert k == length(q)
    return q
end

function MOI.optimize!(optimizer::Optimizer)
    options = optimizer.options
    if optimizer.silent
        options = copy(options)
        options[:printlevel] = 0
    end

    blk1 = Any[]
    blk2 = Any[]
    m = length(optimizer.b)
    At = SparseArrays.SparseMatrixCSC{Float64,Int}[]
    # FIXME I get a strange failure with sparse vectors, need to investigate
    C = Union{Matrix{Float64}, Vector{Float64}}[]
    if !isempty(optimizer.psdc_dims)
        for (i, dim) in enumerate(optimizer.psdc_dims)
            vec_dim = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(dim))
            push!(At, sparse(optimizer.psdc_Avar[i], optimizer.psdc_Acon[i], optimizer.psdc_Aval[i], vec_dim, m))
            c = Vector(sparsevec(optimizer.psdc_Cvar[i], optimizer.psdc_Cval[i], vec_dim))
            Ci = zeros(dim, dim)
            k = 0
            for col in 1:dim
                for row in 1:(col - 1)
                    k += 1
                    Ci[row, col] = c[k] / 2
                    Ci[col, row] = c[k] / 2
                end
                k += 1
                Ci[col, col] = c[k]
            end
            push!(C, Ci)
            push!(blk1, "s")
            push!(blk2, Float64(dim))
        end
    end
    if !isempty(optimizer.quad_dims)
        for (i, dim) in enumerate(optimizer.quad_dims)
            push!(At, sparse(optimizer.quad_Avar[i], optimizer.quad_Acon[i], optimizer.quad_Aval[i], dim, m))
            push!(C, Vector(sparsevec(optimizer.quad_Cvar[i], optimizer.quad_Cval[i], dim)))
            push!(blk1, "q")
            push!(blk2, Float64(dim))
        end
    end
    if !iszero(optimizer.num_nneg)
        push!(At, sparse(optimizer.nneg_Avar, optimizer.nneg_Acon, optimizer.nneg_Aval, optimizer.num_nneg, m))
        push!(C, Vector(sparsevec(optimizer.nneg_Cvar, optimizer.nneg_Cval, optimizer.num_nneg)))
        push!(blk1, "l")
        push!(blk2, Float64(optimizer.num_nneg))
    end
    if !iszero(optimizer.num_free)
        push!(At, sparse(optimizer.free_Avar, optimizer.free_Acon, optimizer.free_Aval, optimizer.num_free, m))
        push!(C, Vector(sparsevec(optimizer.free_Cvar, optimizer.free_Cval, optimizer.num_free)))
        push!(blk1, "u")
        push!(blk2, Float64(optimizer.num_free))
    end
    blk = [blk1 blk2]

    options = optimizer.options
    if optimizer.silent
        options = copy(options)
        options[:printlevel] = 0
    end

    obj, X, y, Z, optimizer.info, optimizer.runhist = sdpt3(
        blk, At, C, optimizer.b; options...)

    optimizer.primal_objective_value, optimizer.dual_objective_value = obj
    k = 0
    optimizer.psdc_X = symvec.(X[k .+ eachindex(optimizer.psdc_dims)])
    optimizer.psdc_Z = symvec.(Z[k .+ eachindex(optimizer.psdc_dims)])
    k += length(optimizer.psdc_dims)
    optimizer.quad_X = X[k .+ eachindex(optimizer.quad_dims)]
    optimizer.quad_Z = Z[k .+ eachindex(optimizer.quad_dims)]
    k += length(optimizer.quad_dims)
    if iszero(optimizer.num_nneg)
        empty!(optimizer.nneg_X)
        empty!(optimizer.nneg_Z)
    else
        k += 1
        optimizer.nneg_X = X[k]
        optimizer.nneg_Z = Z[k]
    end
    if iszero(optimizer.num_free)
        empty!(optimizer.free_X)
        empty!(optimizer.free_Z)
    else
        k += 1
        optimizer.free_X = X[k]
        optimizer.free_Z = Z[k]
    end
    optimizer.y = y
    optimizer.status = optimizer.info["termcode"]
    optimizer.solve_time = optimizer.info["cputime"]
end

const RAW_STATUS = [
    "norm(X) or norm(Z) diverging",
    "dual   problem is suspected to be infeasible",
    "primal problem is suspected to be infeasible",
    "max(relative gap,infeasibility) < gaptol",
    "relative gap < infeasibility",
    "lack of progress in predictor or corrector",
    "X or Z not positive definite",
    "difficulty in computing predictor or corrector direction",
    "progress in relative gap or infeasibility is bad",
    "maximum number of iterations reached",
    "primal infeasibility has deteriorated too much",
    "progress in relative gap has deteriorated",
    "lack of progress in infeasibility"
]
function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    if optimizer.status === nothing
        throw(MOI.OptimizeNotCalled())
    else
        return RAW_STATUS[4 - optimizer.status]
    end
end

const TERMINATION_STATUS = [
    MOI.NUMERICAL_ERROR,
    MOI.DUAL_INFEASIBLE,
    MOI.INFEASIBLE,
    MOI.OPTIMAL,
    MOI.OTHER_ERROR, # TODO what does `relative gap < infeasibility` mean ?
    MOI.SLOW_PROGRESS,
    MOI.NUMERICAL_ERROR,
    MOI.NUMERICAL_ERROR,
    MOI.SLOW_PROGRESS,
    MOI.ITERATION_LIMIT,
    MOI.NUMERICAL_ERROR,
    MOI.NUMERICAL_ERROR,
    MOI.SLOW_PROGRESS
]
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    if optimizer.status === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    else
        return TERMINATION_STATUS[4 - optimizer.status]
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = optimizer.status
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif status == 2
        return MOI.INFEASIBILITY_CERTIFICATE
    else
        # TODO is there solution available in some case here ?
        return MOI.NO_SOLUTION
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = optimizer.status
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif status == 1
        return MOI.INFEASIBILITY_CERTIFICATE
    else
        # TODO is there solution available in some case here ?
        return MOI.NO_SOLUTION
    end
end

MOI.get(::Optimizer, ::MOI.ResultCount) = 1
function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    sign = sense_to_sign(optimizer.objective_sense)
    return sign * optimizer.primal_objective_value + optimizer.objective_constant
end
function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    sign = sense_to_sign(optimizer.objective_sense)
    return sign * optimizer.dual_objective_value + optimizer.objective_constant
end

function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(optimizer, attr)
    info = optimizer.variable_info[vi.value]
    if info.variable_type == FREE
        return optimizer.free_X[info.index_in_cone]
    elseif info.variable_type == NNEG
        return optimizer.nneg_X[info.index_in_cone]
    elseif info.variable_type == QUAD
        return optimizer.quad_X[info.cone_index][info.index_in_cone]
    else
        @assert info.variable_type == PSD
        return optimizer.psdc_X[info.cone_index][info.index_in_cone]
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    info = optimizer.variable_info[ci.value]
    if info.variable_type == FREE
        error("No constraint primal for free variables.")
    elseif info.variable_type == NNEG
        dim = optimizer.nneg_info[info.index_in_cone]
        return optimizer.nneg_X[(info.index_in_cone - 1) .+ (1:dim)]
    elseif info.variable_type == QUAD
        return optimizer.quad_X[info.cone_index]
    else
        @assert info.variable_type == PSD
        return optimizer.psdc_X[info.cone_index]
    end
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.b[ci.value]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    info = optimizer.variable_info[ci.value]
    if info.variable_type == FREE
        error("No constraint dual for free variables.")
    elseif info.variable_type == NNEG
        dim = optimizer.nneg_info[info.index_in_cone]
        return optimizer.nneg_Z[(info.index_in_cone - 1) .+ (1:dim)]
    elseif info.variable_type == QUAD
        return optimizer.quad_Z[info.cone_index]
    else
        @assert info.variable_type == PSD
        return optimizer.psdc_Z[info.cone_index]
    end
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.y[ci.value]
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
using Base.Cartesian

pad0(N::Int, t::Tuple{Vararg{Int}}) = (t..., zeros(Int, N-length(t))...)::NTuple{N,Int}
pad1(N::Int, t::Tuple{Vararg{Int}}) = (t..., ones(Int, N-length(t))...)::NTuple{N,Int}

@generated function im2col{N}(fsize::NTuple{N,Int}, xsize::NTuple{N,Int}, stride::NTuple{N,Int})
  quote
    n = 0
    c = zeros(Int, prod(fsize), 1)
    @nloops $N i j->1:fsize[j] begin
      n += 1
      @nexprs $N j->c[n] += (i_j-1)*prod(xsize[1:(j-1)])
    end

    outsize = ([max(0, div(xsize[i]-fsize[i],stride[i])+1) for i in 1:$N]...)
    ind = ones(Int, 1, prod(outsize))
    n = 0
    @nloops $N i k->1:outsize[k] begin
      n += 1
      @nexprs $N j->ind[n] += (i_j-1)*stride[j]*prod(xsize[1:(j-1)])
    end

    ind.+c
  end
end

macro ind(exp, i, fun)
  @assert exp.head == :ref
  append!(exp.args, @eval map($fun, $i))
  exp
end

function ind(exp, i, fun)
  @assert exp.head == :ref
  append!(exp.args, map(fun, i))
  exp
end

@cuda_include "cudaConv.cu"

# Direct Convolution

immutable ConvOpt{N}
  xsize::NTuple{N,Int}
  fsize::NTuple{N,Int}
  stride::NTuple{N,Int}
  padding::NTuple{N,Int}

  @compat function (::Type{ConvOpt})(fsize::Tuple{Vararg{Int}}, xsize::Tuple{Vararg{Int}},
    stride::Tuple{Vararg{Int}}, padding::Tuple{Vararg{Int}})
    NN = length(xsize)
    @assert length(fsize) <= NN
    @assert length(stride) <= NN
    @assert length(padding) <= NN
    new{NN}(xsize, pad1(NN, fsize), pad1(NN, stride), pad0(NN, padding))
  end
end

function outsize{N}(o::ConvOpt{N}, numout::Int=1)
  ([(div(2*o.padding[i]+o.xsize[i]-o.fsize[i],o.stride[i])+1)*(i > 1 ? 1 : numout) for i in 1:N]...)
end

function Base.conv2(filter::ANode, x::ANode;
  stride=(ones(Int, ndims(x))...), padding=(zeros(Int, ndims(x))...))
  opt = ConvOpt(size(filter)[2:end], size(x), stride, padding)
  ANode(:conv2, (filter, x), outsize(opt, size(filter, 1)), opt)
end

gradient_node(n::ANode{:conv2}, ::Type{Val{1}}, b) = ANode(:conv2_bp_filter, (n.input[2], b), size(n.input[1]), n.arg)
gradient_node(n::ANode{:conv2}, ::Type{Val{2}}, b) = ANode(:conv2_bp_data, (n.input[1], b), size(n.input[2]), n.arg)

@generated function do_forward!{T,N}(n::ANode{:conv2}, y::Array{T,N}, filter::Array{T}, x::Array{T,N})
  quote
    s = n.arg.stride
    p = n.arg.padding
    fill!(y, zero(T))
    @fastmath @inbounds begin
      @nloops $N i d->1:(div(size(x,d) - size(filter,d+1) + 2p[d], s[d]) + 1) begin
        @nloops $N j d->1:size(filter, d+1) begin
          xi = (@nref $N x d->clamp((i_d-1)*s[d]+j_d-p[d], 1, size(x, d)))
          @simd for j_0 in 1:size(filter, 1)
            $(ind(:(y[j_0 + (i_1-1)*size(filter, 1)]), 2:N, d->Symbol("i_",d))) += $(ind(:(filter[]), 0:N, d->Symbol("j_", d))) * xi
          end
        end
      end
    end
  end
end

@generated function do_forward!{T,N}(n::ANode{:conv2_bp_filter}, y::Array{T}, x::Array{T,N}, b::Array{T,N})
  quote
    s = n.arg.stride
    p = n.arg.padding
    fill!(y, zero(T))
    @fastmath @inbounds begin
      @nloops $N i d->1:(div(size(x,d) - size(y,d+1) + 2p[d], s[d]) + 1) begin
        @nloops $N j d->1:size(y, d+1) begin
          xi = (@nref $N x d->clamp((i_d-1)*s[d]+j_d-p[d], 1, size(x, d)))
          @simd for j_0 in 1:size(y, 1)
            bi = $(ind(:(b[j_0 + (i_1-1)*size(y, 1)]), 2:N, d->Symbol("i_",d)))
            $(ind(:(y[]), 0:N, d->Symbol("j_", d))) += bi*xi
          end
        end
      end
    end
  end
end

@generated function do_forward!{T,N}(n::ANode{:conv2_bp_data}, y::Array{T,N}, filter::Array{T}, b::Array{T,N})
  quote
    s = n.arg.stride
    p = n.arg.padding
    fill!(y, zero(T))
    @fastmath @inbounds begin
      @nloops $N i d->1:(div(size(y,d) - size(filter,d+1) + 2p[d], s[d]) + 1) begin
        @nloops $N j d->1:size(filter, d+1) begin
          @simd for j_0 in 1:size(filter, 1)
            bi = $(ind(:(b[j_0 + (i_1-1)*size(filter, 1)]), 2:N, d->Symbol("i_",d)))
            fi = $(ind(:(filter[]), 0:N, d->Symbol("j_", d)))
            (@nref $N y d->clamp((i_d-1)*s[d]+j_d-p[d], 1, size(y, d))) += bi*fi
          end
        end
      end
    end
  end
end

@cuda begin
  const conv3_Float32 = CuFunction(md, "conv3_Float32")
  const conv3_Float64 = CuFunction(md, "conv3_Float64")
  const conv3_bp_data_Float32 = CuFunction(md, "conv3_bp_data_Float32")
  const conv3_bp_data_Float64 = CuFunction(md, "conv3_bp_data_Float64")
  const conv3_bp_filter_Float32 = CuFunction(md, "conv3_bp_filter_Float32")
  const conv3_bp_filter_Float64 = CuFunction(md, "conv3_bp_filter_Float64")

  @generated function do_forward!{T}(n::ANode{:conv2}, y::CudaArray{T,4}, filter::CudaArray{T,4}, x::CudaArray{T,4})
     quote
       o = outsize(n.arg, size(filter,1))
       bs = (8,8,8)
       gs = ((o[1]+bs[1]-1)÷bs[1], (o[2]+bs[2]-1)÷bs[2], o[4]*((o[3]+bs[3]-1)÷bs[3]))
       opt = CudaArray(UInt8, sizeof(n.arg))
       copy!(opt, reinterpret(UInt8, [n.arg]))
       shmem = (bs[1]*n.arg.stride[1]+n.arg.fsize[1])*(bs[2]*n.arg.stride[2]+n.arg.fsize[2])*(bs[3]*n.arg.stride[3]+n.arg.fsize[3])
       cudacall($(Symbol("conv3_", T)), gs, bs, (Ptr, Int, Ptr{T}, Ptr{T}, Ptr{T}), opt.ptr, size(filter, 1), y.ptr, filter.ptr, x.ptr; shmem_bytes=sizeof(T)*shmem)
       free(opt)
     end
  end

  @generated function do_forward!{T}(n::ANode{:conv2_bp_data}, dx::CudaArray{T,4}, f::CudaArray{T,4}, b::CudaArray{T,4})
     quote
       o = size(dx)
       bs = (8,8,8)
       gs = ((o[1]+bs[1]-1)÷bs[1], (o[2]+bs[2]-1)÷bs[2], o[4]*((o[3]+bs[3]-1)÷bs[3]))
       opt = CudaArray(UInt8, sizeof(n.arg))
       copy!(opt, reinterpret(UInt8, [n.arg]))
       cudacall($(Symbol("conv3_bp_data_", T)), gs, bs, (Ptr, Int, Ptr{T}, Ptr{T}, Ptr{T}), opt.ptr, size(f, 1), dx.ptr, f.ptr, b.ptr)
       free(opt)
     end
  end

  @generated function do_forward!{T}(n::ANode{:conv2_bp_filter}, df::CudaArray{T,4}, x::CudaArray{T,4}, b::CudaArray{T,4})
     quote
       o = size(df)
       bs = (8,8,8)
       gs = ((o[1]*o[2]+bs[1]-1)÷bs[1], (o[3]+bs[2]-1)÷bs[2], (o[4]+bs[3]-1)÷bs[3])
       opt = CudaArray(UInt8, sizeof(n.arg))
       copy!(opt, reinterpret(UInt8, [n.arg]))
       shmem = (bs[1]*n.arg.stride[1]+n.arg.fsize[1])*(bs[2]*n.arg.stride[2]+n.arg.fsize[2])*(bs[3]*n.arg.stride[3]+n.arg.fsize[3])
       cudacall($(Symbol("conv3_bp_filter_", T)), gs, bs, (Ptr, Int, Ptr{T}, Ptr{T}, Ptr{T}), opt.ptr, size(df, 1), df.ptr, x.ptr, b.ptr; shmem_bytes=sizeof(T)*shmem)
       free(opt)
     end
  end
end

# Pooling

function pool(x::ANode, sz::Tuple{Vararg{Int}}; stride=sz, padding=(0,))
  opt = ConvOpt(sz, size(x), stride, padding)
  ANode(:pool, (x,), outsize(opt), opt)
end

gradient_node(n::ANode{:pool}, ::Type{Val{1}}, b) = ANode(:pool_bp, (n.input[1], n, b), size(n.input[1]), n.arg)

@generated function do_forward!{T,N}(n::ANode{:pool}, y::Array{T,N}, x::Array{T,N})
  quote
    s = n.arg.stride
    p = n.arg.padding
    sz = n.arg.fsize
    @fastmath @inbounds begin
      @nloops $N i d->1:(div(size(x,d) - sz[d] + 2p[d], s[d]) + 1) begin
        a = @nref $N x d->clamp((i_d-1)*s[d]+1, 1, size(x, d))
        @nloops $N j d->1:sz[d] begin
          b = @nref $N x d->clamp((i_d-1)*s[d]+j_d-p[d], 1, size(x, d))
          if b > a
            a = b
          end
        end
        (@nref $N y d->i_d) = a
      end
    end
  end
end

@generated function do_forward!{T,N}(n::ANode{:pool_bp}, y::Array{T,N}, x::Array{T,N}, z::Array{T,N}, b::Array{T,N})
  quote
    s = n.arg.stride
    p = n.arg.padding
    sz = n.arg.fsize
    fill!(y, zero(T))
    @fastmath @inbounds begin
      @nloops $N i d->1:(div(size(x,d) - sz[d] + 2p[d], s[d]) + 1) begin
        zi = (@nref $N z d->i_d)
        @nloops $N j d->1:sz[d] begin
          if (@nref $N x d->clamp(((i_d-1)*s[d]+j_d), 1, size(x, d))) == zi
            (@nref $N y d->clamp((i_d-1)*s[d]+j_d, 1, size(y, d))) += (@nref $N b d->i_d)
          end
        end
      end
    end
  end
end

@cuda begin
  const pool_Float32 = CuFunction(md, "pool_Float32")
  const pool_Float64 = CuFunction(md, "pool_Float64")
  const pool_bp_Float32 = CuFunction(md, "pool_bp_Float32")
  const pool_bp_Float64 = CuFunction(md, "pool_bp_Float64")

  @generated function do_forward!{T}(n::ANode{:pool}, y::CudaArray{T,4}, x::CudaArray{T,4})
     quote
       o = size(y)
       bs = (8,8,8)
       gs = ((o[1]+bs[1]-1)÷bs[1], (o[2]+bs[2]-1)÷bs[2], o[4]*((o[3]+bs[3]-1)÷bs[3]))
       opt = CudaArray(UInt8, sizeof(n.arg))
       copy!(opt, reinterpret(UInt8, [n.arg]))
       shmem = (bs[1]*n.arg.stride[1]+n.arg.fsize[1])*(bs[2]*n.arg.stride[2]+n.arg.fsize[2])*(bs[3]*n.arg.stride[3]+n.arg.fsize[3])
       cudacall($(Symbol("pool_", T)), gs, bs, (Ptr, Ptr{T}, Ptr{T}), opt.ptr, y.ptr, x.ptr; shmem_bytes=sizeof(T)*shmem)
       free(opt)
     end
  end

  @generated function do_forward!{T}(n::ANode{:pool_bp}, dx::CudaArray{T,4}, x::CudaArray{T,4}, y::CudaArray{T,4}, b::CudaArray{T,4})
     quote
       fill!(dx, zero(eltype(dx)))
       o = size(dx)
       bs = (8,8,8)
       gs = ((o[1]+bs[1]-1)÷bs[1], (o[2]+bs[2]-1)÷bs[2], o[4]*((o[3]+bs[3]-1)÷bs[3]))
       opt = CudaArray(UInt8, sizeof(n.arg))
       copy!(opt, reinterpret(UInt8, [n.arg]))
       cudacall($(Symbol("pool_bp_", T)), gs, bs, (Ptr, Ptr{T}, Ptr{T}, Ptr{T}, Ptr{T}), opt.ptr, dx.ptr, x.ptr, y.ptr, b.ptr)
       free(opt)
     end
  end
end

ext: jl


text:
using ProjectEuler
using Test
cd("src")
@testset "ProjectEuler.jl" begin
    @test solution1() == 233168
    @test solution2() == 4613732
    @test solution3() == 6857
    @test solution4() == 906609
    @test solution5() == 232792560
    @test solution6() == 25164150
    @test solution7() == 104743
    @test solution8() == 23514624000
    @test solution9() == 31875000
    @test solution10() == 142913828922
    @test solution11() == 70600674
    @test solution12() == 76576500
    @test solution13() == 5537376230
    @test solution14() == 837799
    @test solution15() == 137846528820
    @test solution16() == 1366
    @test solution17() == 21124
    @test solution18() == 1074
    @test solution19() == 171
    @test solution20() == 648
    @test solution21() == 31626
    @test solution22() == 871198282
    @test solution23() == 4179871
    @test solution24() == 2783915460
    @test solution25() == 4782
    @test solution26() == 983
    @test solution27() == -59231
    @test solution28() == 669171001
    @test solution29() == 9183
    @test solution30() == 443839
    @test solution31() == 73682
    @test solution32() == 45228
    @test solution33() == 100
    @test solution34() == 40730
    @test solution35() == 55
    @test solution36() == 872187
    @test solution37() == 748317
end


text:
using AtomicLevels
import AtomicLevels: @hi_str
using WignerSymbols
using Test

include("parity.jl")
include("orbitals.jl")
include("configurations.jl")
include("excited_configurations.jl")
include("terms.jl")
include("jj_terms.jl")
include("couple_terms.jl")
include("csfs.jl")
include("levels.jl")
include("jj2lsj.jl")


text:
using NNlib

# Activation functions

@cufunc σ(x) = ifelse(x < -80, zero(x), one(x) / (one(x) + exp(-x)))

@cufunc function logσ(x)
  max_v = max(zero(x), -x)
  z = exp(-max_v) + exp(-x-max_v)
  -(max_v + log(z))
end

@cufunc elu(x, α = one(x)) =
  ifelse(x ≥ 0, x/1, α * (exp(x) - one(x)))

@cufunc swish(x) = x * σ(x)

@cufunc function gelu(x)
  λ = oftype(x/1, √(2/π))
  α = oftype(x/1, 0.044715)
  h = oftype(x/1, 0.5)
  h * x * (one(x) + tanh(λ * (x + α * x^3)))
end

@cufunc function selu(x)
  λ = oftype(x/1, 1.0507009873554804934193349852946)
  α = oftype(x/1, 1.6732632423543772848170429916717)
  λ * ifelse(x > 0, x/1, α * (exp(x) - 1))
end

@cufunc softplus(x) = ifelse(x > 0, x + log1p(exp(-x)), log1p(exp(x)))

# Batched matrix multiplication

const batched_gemm_args = [
    (:(CuArray{T, 3}), 'N'),
    (:(NNlib.BatchedTranspose{T, <:CuArray{T, 3}}), 'T'),
    (:(NNlib.BatchedAdjoint{T, <:CuArray{T, 3}}), 'C')
]

for (TA, transA) in batched_gemm_args, (TB, transB) in batched_gemm_args
    @eval function NNlib.batched_mul!(C::CuArray{T, 3}, A::$TA, B::$TB) where {T<:CUBLAS.CublasFloat}
        CuArrays.CUBLAS.gemm_strided_batched!($transA, $transB, one(T), NNlib._unbatch(A), NNlib._unbatch(B), zero(T), C)
        C
    end
end


text:
# BlackScholesModel

function initial_state(ctx::CompilerContext{B}) where {B<:BlackScholesModel}
    get_pricing_date(ctx)
end

# An EuropeanOption is recognized as a Call
# when the `Both` contract has a `Unit{Stock} leg.
function is_european_call_option(c::EuropeanOption)
    both_contract = c.c.c1::Both
    return (isa(both_contract.c1, Unit{Stock}) || isa(both_contract.c2, Unit{Stock})) && (isa(both_contract.c1, Give) || isa(both_contract.c2, Give))
end

# An EuropeanOption is recognized as a Put
# when the `Both` contract has a Give(Unit{Stock}) leg.
function is_european_put_option(c::EuropeanOption)
    both_contract = c.c.c1::Both

    if !(isa(both_contract.c1, Give) || isa(both_contract.c2, Give))
        return false
    end

    give_contract = isa(both_contract.c1, Give) ? both_contract.c1 : both_contract.c2
    @assert isa(give_contract, Give)
    return isa(give_contract.c, Unit{Stock})
end

is_european_call_option(c::Contract) = false
is_european_put_option(c::Contract) = false

struct EuropeanOptionLegs
    call_or_put::Symbol
    stock::Unit{Stock}
    strike::Contract
end

function EuropeanOptionLegs(c::EuropeanOption)
    both_contract = c.c.c1::Both

    if is_european_call_option(c)
        call_or_put = :call
        if isa(both_contract.c1, Unit{Stock})
            return EuropeanOptionLegs(call_or_put, both_contract.c1, both_contract.c2.c)
        else
            return EuropeanOptionLegs(call_or_put, both_contract.c2, both_contract.c1.c)
        end
    elseif is_european_put_option(c)
        call_or_put = :put
        if isa(both_contract.c1, Give)
            return EuropeanOptionLegs(call_or_put, both_contract.c1.c, both_contract.c2)
        else
            return EuropeanOptionLegs(call_or_put, both_contract.c2.c, both_contract.c1)
        end
    end

    error("Contract was not recognized as an European Option: $c.")
end

function is_european_option(c::EuropeanOption)
    return is_european_call_option(c) || is_european_put_option(c)
end

is_european_option(c::Contract) = false
check_is_european_option(c::EuropeanOption) = @assert is_european_option(c) "Contract is not an European Option: $c."

# WhenAt -> Either{Both, Worthless}
# Call: Both -> Unit{Stock}, Give(c)
# Put:  Both -> c, Give(Unit{Stock})
function lower!(ctx::CompilerContext{B}, c::EuropeanOption, at::Date) :: OptimizingIR.ImmutableValue where {B<:BlackScholesModel}

    check_is_european_option(c)
    assert_at_initial_state(ctx, at)

    maturity = c.maturity
    legs = EuropeanOptionLegs(c)
    underlying_stock_observable = legs.stock.o

    yf_to_maturity = yearfraction(ctx, maturity)
    arg_years_to_maturity = lower!(ctx, Konst(yf_to_maturity), at)
    arg_stock = lower!(clone_context_with_model(ctx, ctx.model.static_model), WhenAt(maturity, legs.stock), at)
    arg_strike = lower!(clone_context_with_model(ctx, ctx.model.static_model), WhenAt(maturity, legs.strike), at)
    arg_vol = lower!(ctx, Volatility(underlying_stock_observable), at)

    op_bs = legs.call_or_put == :call ? OP_BLACK_FORMULA_CALL : OP_BLACK_FORMULA_PUT
    return add_instruction!(ctx, op_bs, arg_stock, arg_strike, arg_years_to_maturity, arg_vol)
end

function yearfraction(ctx::CompilerContext{B}, maturity::Date) where {B<:BlackScholesModel}
    # TODO : maybe look at the convention of the riskfree curve of the functional currency...
    if get_functional_currency(ctx.model) == Currencies.BRL
        return InterestRates.yearfractionvalue(InterestRates.BDays252(CAL_BRL), get_pricing_date(ctx), maturity)
    else
        error("Functional currency $(get_functional_currency(ctx.model)) not supported by `Compiler.yearfraction`.")
    end
end

function lower!(ctx::CompilerContext{B}, v::Volatility, at::Date) :: OptimizingIR.ImmutableValue where {B<:BlackScholesModel}
    assert_at_initial_state(ctx, at)
    return add_input_riskfactor!(ctx, v)
end


text:
#= Polyline encoder and decoder

This module performes an encoding and decoding for
GPS coordinates into a polyline using the algorithm
detailed in:
https://developers.google.com/maps/documentation/utilities/polylinealgorithm

Example:
    julia> enc = encodePolyline([[38.5 -120.2]; [40.7 -120.95]; [43.252 -126.453]])
    "_p~iF~ps|U_ulLnnqC_mqNvxq`@"

    julia> dec = decodePolyline("_p~iF~ps|U_ulLnnqC_mqNvxq`@")
    [[38.5 -120.2]; [40.7 -120.95]; [43.252 -126.453]]

=#

__precompile__()

module Polyline

using HTTP

# coordinate is a type to represent a GPS data point
struct coordinate{T}
    Lat::T
    Lng::T
end

include("encode.jl")
include("decode.jl")
include("loadGPX.jl")
include("googleMaps.jl")

export decodePolyline, encodePolyline, getMapImage, mapsURL

end # module


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

# StreamInfo.jl: type and method definitions for steam information object

import Base.unsafe_convert
import Base.show

export StreamInfo
export name, type, channel_count, nominal_srate, channel_format, source_id
export version, created_at, session_id, uid, hostname
export desc, XML

"""
    StreamInfo

The StreamInfo type stores the declaration of a data stream.

Represents the following information:
 a) stream data format (#channels, channel format)
 b) core information (stream name, content type, sampling rate)
 c) optional meta-data about the stream content (channel labels,
    measurement units, etc.)

Whenever a program wants to provide a new stream on the lab network it will typically first
create a StreamInfo to describe its properties and then construct a StreamOutlet with it to
create the stream on the network. Recipients who discover the outlet can query the
StreamInfo; it is also written to disk when recording the stream (playing a similar role as
a file header).
"""
mutable struct StreamInfo{T}
  handle::lib.lsl_streaminfo
  format::Type{T}

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

end

# Define conversion to pointer
unsafe_convert(::Type{lib.lsl_streaminfo}, info::StreamInfo) = info.handle

# Define constructor
"""
    StreamInfo(name = "untitled";
               type = "",
               channel_count = 1,
               nominal_srate = LSL_IRREGULAR_RATE,
               channel_format = Float32,
               source_id = "")

Construct a new stream information object.

Core stream information is specified here. Any remaining meta-data can be added later.

# Keyword arguments
- `name::String`: Name of the stream. Describes the device (or product series)
                  that this stream makes available (for use by programs, experimenters or
                  data analysts). Cannot be empty.
- `type::String`: Content type of the stream. Please see https://github.com/sccn/xdf/wiki/Meta-Data
                  (or web search for: XDF meta-data) for pre-defined content-type names, but
                  you can also make up your own. The content type is the preferred way to
                  find streams (as opposed to searching by name).
- `channel_count::Integer`:   Number of channels per sample. This stays constant for the
                              lifetime of the stream.
- `nominal_srate::Number`:    The sampling rate (in Hz) as advertised by the datasource, if
                              regular (otherwise set to LSL_IRREGULAR_RATE).
- `channel_format::DataType`: Format/type of each channel.If your channels have different
                              formats, consider supplying multiple streams or use the
                              largest type that can hold them all (such as Float64).
 - `source_id::String`:       Unique identifier of the source or device, if available
                              (such as the serial number). Allows recipients to recover from
                              failure even after the serving app or device crashes. May in
                              some cases also be constructed from device settings.
"""
function StreamInfo(; name::String = "untitled",
                      type::String = "",
                      channel_count::Integer = 1,
                      nominal_srate::Number = IRREGULAR_RATE,
                      channel_format::DataType = Float32,
                      source_id::String = "")

  # Get LSL format, create and test handle
  format = _lsl_channel_format(channel_format)
  handle = lib.lsl_create_streaminfo(name, type, channel_count, nominal_srate, format, source_id)
  if handle == C_NULL
    error("liblsl library returned NULL pointer during StreamInfo creation")
  end

  return StreamInfo(handle, channel_format)
end

# Destroy a StreamInfo handle
function _destroy(info::StreamInfo)
  if info.handle != C_NULL
    lib.lsl_destroy_streaminfo(info)
    info.handle = C_NULL
  end
  return nothing
end

# function Base.show(io::IO, info::StreamInfo)
#   print(io, "Stream ", source_id(info), "\tname: ", name(info), "\ttype: ", type(info))
# end

#
# Core information retrieval (assigned at construction)
#

"""
    name(info::StreamInfo)

Return name of the stream `info`.

This is a human-readable name. For streams offered by device modules, it refers to the type
of device or product series that is generating the data of the stream. If the source is an
application, the name may be a more generic or specific identifier. Multiple streams with
the same name can coexist, though potentially at the cost of ambiguity (for the recording
app or experimenter).
"""
name(info::StreamInfo) = unsafe_string(lib.lsl_get_name(info))

"""
    type(info::StreamInfo)

Return the content type of the stream `info`.

The content type is a short string such as "EEG", "Gaze" which describes the content carried
by the channel (if known). If a stream contains mixed content this value need not be
assigned but may instead be stored in the description of channel types. To be useful to
applications and automated processing systems using the recommended content types is
preferred. Content types usually follow those pre-defined in https://github.com/sccn/xdf/wiki/Meta-Data
(or web search for: XDF meta-data).
"""
type(info::StreamInfo) = unsafe_string(lib.lsl_get_type(info))

"""
    channel_count(info::StreamInfo)

Return mumber of channels of the stream `info`.

A stream has at least one channel; the channel count stays constant for all samples.
"""
channel_count(info::StreamInfo) = lib.lsl_get_channel_count(info)

"""
    nominal_srate(info::StreamInfo)

Return sampling rate of the stream `info`, according to the source (in Hz).

If a stream is irregularly sampled, this should be set to IRREGULAR_RATE.

Note that no data will be lost even if this sampling rate is incorrect or if a device has
temporary hiccups, since all samples will be recorded anyway (except for those dropped by
the device itself). However, when the recording is imported into an application, a good
importer may correct such errors more accurately if the advertised sampling rate was close
to the specs of the device.
"""
nominal_srate(info::StreamInfo) = lib.lsl_get_nominal_srate(info)

"""
    channel_format(info::StreamInfo)

Return (Julia) channel format of the stream `info`.

All channels in a stream have the same format. However, a device might offer multiple
time-synched streams each with its own format.

Note that the format is returned from the parametric type of the StreamInfo object, it is
not queried from the library.
"""
channel_format(info::StreamInfo{T}) where T = T

"""
    source_id(info::StreamInfo)

Return unique identifier of the stream's source, if available.

The unique source (or device) identifier is an optional piece of information that, if
available, allows that endpoints (such as the recording program) can re-acquire a stream
automatically once it is back online.
"""
source_id(info::StreamInfo) = unsafe_string(lib.lsl_get_source_id(info))

#
# Hosting information (assigned when bound to an outlet/inlet)
#

"""
    version(info::SteramInfo)

Return protocol version used to deliver the stream.
"""
function version(info::StreamInfo)
  major, minor = divrem(lib.lsl_get_version(info), 100)
  return VersionNumber(major, minor)
end

"""
    created_at(info::StreamInfo)

Return creation time stamp of the stream.

This is the time stamp when the stream was first created (as determined via local_clock()
on the providing machine).
"""
created_at(info::StreamInfo) = lib.lsl_get_created_at(info)

"""
    uid(info::StreamInfo)

Return unique ID of the stream outlet (once assigned).

This is a unique identifier of the stream outlet, and is guaranteed to be different across
multiple instantiations of the same outlet (e.g., after a re-start).
"""
uid(info::StreamInfo) = unsafe_string(lib.lsl_get_uid(info))

"""
    session_id(info::StreamInfo)

Return session ID for the given stream.

The session id is an optional human-assigned identifier of the recording session. While it
is rarely used, it can be used to prevent concurrent recording activitites in the same
sub-network (e.g., in multiple experiment areas) from seeing each other's streams (assigned
via a configuration file by the experimenter, see Network Connectivity on the LSL wiki).
"""
session_id(info::StreamInfo) = unsafe_string(lib.lsl_get_session_id(info))

"""
    hostname(info::StreamInfo)

Return hostname of the providing machine (once bound to an outlet). Modification is not
permitted.
"""
hostname(info::StreamInfo) = unsafe_string(lib.lsl_get_hostname(info))

#
# Data description (can be modified)
#

"""
    desc(info::StreamInfo)

Get extended description of the stream as an XMLElement.

It is highly recommended that at least the channel labels are described here. See code
examples on the LSL wiki. Other information, such as amplifier settings, measurement units
if deviating from defaults, setup information, subject information, etc., can be specified
here, as well. Meta-data recommendations follow the XDF file format project
(github.com/sccn/xdf/wiki/Meta-Data or web search for: XDF meta-data).

Important: if you use a stream content type for which meta-data recommendations exist, please
try to lay out your meta-data in agreement with these recommendations for compatibility with
other applications.
"""
desc(info::StreamInfo) = LSLXMLElement(lib.lsl_get_desc(info))

"""
    XML(info::StreamInfo)

Return XML document (from LightXML) containing the entire StreamInfo.

This yields an XML document whose top-level element is `<info>`. The info element contains
one element for each field of the StreamInfo instance, including:

- the core elements `<name>`, `<type>`, `<channel_count`, `<nominal_srate>`,
  `<channel_format>`, `<source_id>`
- the misc elements `<version>`, `<created_at>`, `<uid>`, `<session_id>`, `<v4address>`,
  `<v4data_port>`, `<v4service_port>`, `<v6address>`, `<v6data_port>`, `<v6service_port>`
- the extended description element `<desc>` with user-defined sub-elements.
"""
function XML(info::StreamInfo)
    xml_string_ptr = lib.lsl_get_xml(info)
    if xml_string_ptr == C_NULL
      error("liblsl returned NULL pointer when requesting XML stream information")
    end
    xmldoc = LightXML.parse_string(unsafe_string(xml_string_ptr))
    lib.lsl_destroy_string(xml_string_ptr)
    return xmldoc
end

# /// Number of bytes occupied by a channel (0 for string-typed channels).
# extern LIBLSL_C_API int32_t lsl_get_channel_bytes(lsl_streaminfo info);

# /// Number of bytes occupied by a sample (0 for string-typed channels).
# extern LIBLSL_C_API int32_t lsl_get_sample_bytes(lsl_streaminfo info);

# /**
#  * Tries to match the stream info XML element @p info against an
#  * XPath query.
#  *
#  * Example query strings:
#  * @code
#  * channel_count>5 and type='EEG'
#  * type='TestStream' or contains(name,'Brain')
#  * name='ExampleStream'
#  * @endcode
#  */
# extern LIBLSL_C_API int lsl_stream_info_matches_query(lsl_streaminfo info, const char* query);

# /// Create a streaminfo object from an XML representation
# extern LIBLSL_C_API lsl_streaminfo lsl_streaminfo_from_xml(const char *xml);


text:
module AdvectionEquationMod

using PDESolver  # setupf LOAD_PATH to find PDESolver components
using SolverCommon
using ArrayViews
import ArrayViews.view
using ODLCommonTools
using SummationByParts
using PdePumiInterface
#using NonlinearSolvers
using Jacobian
using LinearSolvers
using MPI
using Utils
#using PETSc

import ODLCommonTools.sview
using Input

export AdvectionData, AdvectionData_, run_advection #getMass, array3DTo1D, array1DTo3D

# include("advection.jl")
# include("getMass.jl")

"""
  Direct subtype of [`AbstractSolutionData`](@ref), inheriting `Tsol` and
  `Tres` as static parameter
"""
abstract type AbstractAdvectionData{Tsol, Tres} <: AbstractSolutionData{Tsol, Tres} end

"""
  Subtype of [`AbstractAdvectionData`](@ref), inheriting its static parameters
  and adding `Tdim`.
"""
abstract type AdvectionData{Tsol, Tres, Tdim} <: AbstractAdvectionData{Tsol, Tres} end

include("types.jl")
include("functionals.jl")
include(joinpath(Pkg.dir("PDESolver"), "src/solver/debug.jl"))  # debug macro
include("advection.jl")
include("common_funcs.jl")
include("bc.jl")
include("bc_solvers.jl")
include("ic.jl")
include("GLS.jl")
include("GLS2.jl")
include("boundary_functional.jl")
include("functional_deriv.jl")
include("source.jl")
include("flux.jl")
include("check_options.jl")
include("eqn_deepcopy.jl")
include("startup_func.jl")  # function to invoke the solver

# register this physics module
"""
  This physics is named `Advection`
"""
global const PhysicsName = "Advection"
register_physics(PhysicsName, AdvectionEquationMod, createObjects, checkOptions)

end # end module


text:
immutable NoncentralBeta{T<:Real} <: ContinuousUnivariateDistribution
    α::T
    β::T
    λ::T

    function (::Type{NoncentralBeta{T}}){T}(α::T, β::T, λ::T)
        @check_args(NoncentralBeta, α > zero(α) && β > zero(β))
        @check_args(NoncentralBeta, λ >= zero(λ))
        new{T}(α, β, λ)
    end
end

NoncentralBeta{T<:Real}(α::T, β::T, λ::T) = NoncentralBeta{T}(α, β, λ)
NoncentralBeta(α::Real, β::Real, λ::Real) = NoncentralBeta(promote(α, β, λ)...)
NoncentralBeta(α::Integer, β::Integer, λ::Integer) = NoncentralBeta(Float64(α), Float64(β), Float64(λ))

@distr_support NoncentralBeta 0.0 1.0

### Parameters

params(d::NoncentralBeta) = (d.α, d.β, d.λ)
@inline partype{T<:Real}(d::NoncentralBeta{T}) = T

### Evaluation & Sampling

# TODO: add mean and var

@_delegate_statsfuns NoncentralBeta nbeta α β λ

function rand(d::NoncentralBeta)
    a = rand(NoncentralChisq(2d.α, d.β))
    b = rand(Chisq(2d.β))
    a / (a + b)
end


text:
function leaves2roottrunc(tensor::Array{Float64},dimtree,ranks)
    # Source: Algorithm 2 of Hierarchical SVD of Tensors by Lars Grasedyck
    N  = ndims(tensor);
    sz = size(tensor);
    Ut = Matrix{Float64};
    Cp = unfold(tensor,[1],"left");
    for t = 1:N
        At = Cp;
        F  = svd!(At);
        Ut = F.U[:,1:ranks[t]];
        Cp = Ut'*Cp;
        Cp = reshape(Cp,(ranks[1:t]..., sz[t+1:end]...));
        if t<N
            Cp = unfold(Cp,[t+1],"left");
        end
    end
    Cl = Cp;
    for l = p-1:-1:1
        for t = 1:dimtree[l]

        end
    end

end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using Forecasts
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "NaiveForecast" begin
    data = [1, 2, 3]
    fc = naivef(data)
    resultA = collect(take(fc, 1))
    @test resultA == [3]
    resultB = collect(take(fc, 5))
    @test resultB == [3, 3, 3, 3, 3]
end

@testset "SeasonalNaiveForecast" begin
    data = [1., 2., 3., 4., 1., 2., 3., 4.]
    fc = snaivef(data, 4)
    result = collect(take(fc, 1))
    @test result == [1.]
    result = collect(take(fc, 4))
    @test result == [1., 2., 3., 4.]
end

@testset "MeanForecast" begin
    data = [1., 2., 3., 4., 1., 2., 3., 4.]
    fc = meanf(data)
    result = collect(take(fc, 1))
    @test result == [2.5]
    result = collect(take(fc, 3))
    @test result == [2.5, 2.5, 2.5]
end

@testset "RwDriftForecast" begin
    data = collect(Float64, 1:10)
    fc = rwf(data)
    result = collect(take(fc, 3))
    @test result == [11., 12., 13.]
end

@testset "AR" begin
end

@testset "forecast" begin
    data = collect(Float64, 1:10)
    fc = rwf(data)
    result = forecast(fc, 3)
    @test result == [11., 12., 13.]
end


text:
using Distributed
using Printf

@show procs()
@show rmprocs(procs()[2:end])
@show addprocs()
@show procs()
println()

function calc_julia_distributed()
    total = @distributed (+) for i in 1:10^5
        tmp2 = 0

        for j in 1:10^4
            tmp = if j % 4 == 0
                     i + j
                  elseif j % 4 == 1
                     i - j
                  elseif j % 4 == 2
                     i * j
                  else
                     i ÷ j
                  end
            tmp2 += tmp
        end
        tmp2
    end

    return total
end

@time total = calc_julia_distributed()
@time total = calc_julia_distributed()
@time total = calc_julia_distributed()

@printf("total = %d\n", total)


text:
using DitherPunk
using Test

@testset "DitherPunk.jl" begin
    @testset "Bayer matrices" begin
        include("test_bayer.jl")
    end
    @testset "Gradient image" begin
        include("test_gradient.jl")
    end
    @testset "Color image" begin
        @testset "Custom palette" begin
            include("test_color.jl")
        end
        @testset "Fixed palette" begin
            include("test_fixed_color.jl")
        end
    end
    @testset "Utilities" begin
        include("test_utils.jl")
    end
end


text:
################################################################################
# PANOC iterator (with L-BFGS directions)

mutable struct PANOCIterator{I <: Integer, R <: Real, D, CS, FS, AS, CQ, FQ, AQ, G, HH} <: ProximalAlgorithm{I,D}
    x::D
    fs::FS
    As::AS
    fq::FQ
    Aq::AQ
    g::G
    gamma::R
    maxit::I
    tol::R
    adaptive::Bool
    verbose::I
    verbose_freq::I
    alpha::R
    beta::R
    tau::R
    y::D # gradient step
    xbar::D # proximal-gradient step
    H::HH # inverse Jacobian approximation
    FPR_x::D
    Aqx::CQ
    Asx::CS
    Aqxnew::CQ
    Asxnew::CS
    Aqd::CQ
    Asd::CS
    Aqfb::CQ
    Asfb::CS
    gradfq_Aqx::CQ
    gradfs_Asx::CS
    Aqxbar::CQ
    Asxbar::CS
    gradfq_Aqxbar::CQ
    gradfs_Asxbar::CS
    fs_Asx::R
    fq_Aqx::R
    f_Ax::R
    At_gradf_Ax::D
    Aqt_gradfq_Aqx::D
    Ast_gradfs_Asx::D
    g_xbar::R
    FBE_x::R
    FBE_x_prev::R
    x_prev::D
    FPR_x_prev::D
    xnew::D
    xnewbar::D
    FPR_xnew::D
    d::D
    normFPR_x::R
end

################################################################################
# Constructor

function PANOCIterator(x0::D;
                       fs::FS=Zero(), As::AS=Identity(size(x0)),
                       fq::FQ=Zero(), Aq::AQ=Identity(size(x0)),
                       g::G=Zero(),
                       gamma::R=-1.0, maxit::I=10000, tol::R=1e-4, adaptive::Bool=false,
                       memory::I=10, verbose::I=1, verbose_freq::I=100,
                       alpha::R=0.95, beta::R=0.5) where {I, R, D, FS, AS, FQ, AQ, G}
    x       = copy(x0)
    xbar    = zero(x0)
    x_prev  = zero(x0)
    xnew    = zero(x0)
    xnewbar = zero(x0)
    y = zero(x0)
    FPR_x = zero(x0)
    FPR_x_prev = zero(x0)
    FPR_xnew = zero(x0)
    Aqx = Aq*x
    Asx = As*x
    Aqxnew = zero(Aqx)
    Asxnew = zero(Asx)
    Aqd = zero(Aqx)
    Asd = zero(Asx)
    Aqfb = zero(Aqx)
    Asfb = zero(Asx)
    gradfq_Aqx = zero(Aqx)
    gradfs_Asx = zero(Asx)
    Aqxbar = zero(Aqx)
    Asxbar = zero(Asx)
    gradfq_Aqxbar = zero(Aqx)
    gradfs_Asxbar = zero(Asx)
    At_gradf_Ax = zero(x0)
    Aqt_gradfq_Aqx = zero(x0)
    Ast_gradfs_Asx = zero(x0)
    d = zero(x0)
    H = LBFGS(x, memory)
    CQ = typeof(Aqx)
    CS = typeof(Asx)
    HH = typeof(H)
    PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}(
               x, fs, As,
               fq, Aq, g,
               gamma, maxit, tol,
               adaptive, verbose,
               verbose_freq, alpha, beta,
               one(R), y, xbar,
               H, FPR_x,
               Aqx, Asx, Aqxnew, Asxnew, Aqd, Asd, Aqfb, Asfb, gradfq_Aqx, gradfs_Asx,
               Aqxbar, Asxbar, gradfq_Aqxbar, gradfs_Asxbar,
               zero(R), zero(R), zero(R),
               At_gradf_Ax, Aqt_gradfq_Aqx, Ast_gradfs_Asx,
               zero(R), zero(R),
               zero(R), x_prev, FPR_x_prev,
               xnew, xnewbar, FPR_xnew,
               d, zero(R))
end

################################################################################
# Utility methods

maxit(sol::PANOCIterator{I}) where {I} = sol.maxit

converged(sol::PANOCIterator{I,R,D}, it::I)  where {I,R,D}= it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol

verbose(sol::PANOCIterator) = sol.verbose > 0
verbose(sol::PANOCIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0))

function display(sol::PANOCIterator)
    @printf("%6s | %10s | %10s | %10s | %10s |\n ", "it", "gamma", "fpr", "tau", "FBE")
    @printf("------|------------|------------|------------|------------|\n")
end

function display(sol::PANOCIterator, it)
    @printf("%6d | %7.4e | %7.4e | %7.4e | %7.4e | \n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma, sol.tau, sol.FBE_x)
end

function Base.show(io::IO, sol::PANOCIterator)
    println(io, "PANOC" )
    println(io, "fpr        : $(maximum(abs,sol.FPR_x))")
    println(io, "gamma      : $(sol.gamma)")
    println(io, "tau        : $(sol.tau)")
    print(  io, "FBE        : $(sol.FBE_x)")
end

################################################################################
# Initialization

function initialize!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}) where {I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}

    # reset L-BFGS operator (would be nice to have this option)
    # TODO add function reset!(::LBFGS) in AbstractOperators
    sol.H.currmem, sol.H.curridx = 0, 0
    sol.H.H = 1.0

    # compute first forward-backward step here
    mul!(sol.Aqx, sol.Aq, sol.x)
    sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqx)
    mul!(sol.Asx, sol.As, sol.x)
    sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx)
    sol.At_gradf_Ax .= sol.As'*sol.gradfs_Asx .+ sol.Aq'*sol.gradfq_Aqx
    sol.f_Ax = sol.fs_Asx + sol.fq_Aqx

    if sol.gamma <= 0.0 # estimate L in this case, and set gamma = 1/L
        # this part should be as follows:
        # 1) if adaptive = false and only fq is present then L is "accurate"
        # 2) otherwise L is "inaccurate" and set adaptive = true
        # TODO: implement case 1), now 2) is always performed
        xeps = sol.x .+ sqrt(eps())
        Aqxeps = sol.Aq*xeps
        gradfq_Aqxeps, = gradient(sol.fq, Aqxeps)
        Asxeps = sol.As*xeps
        gradfs_Asxeps, = gradient(sol.fs, Asxeps)
        At_gradf_Axeps = sol.As'*gradfs_Asxeps .+ sol.Aq'*gradfq_Aqxeps
        L = norm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*length(xeps)))
        sol.adaptive = true
        # in both cases set gamma = 1/L
        sol.gamma = sol.alpha/L
    end

    sol.y .= sol.x .-sol.gamma .* sol.At_gradf_Ax
    sol.g_xbar = prox!(sol.xbar, sol.g, sol.y, sol.gamma)
    sol.FPR_x .= sol.x .- sol.xbar

    sol.normFPR_x = norm(sol.FPR_x)
    sol.FBE_x = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*sol.normFPR_x^2 + sol.g_xbar

    return sol.xbar

end

################################################################################
# Iteration

function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it::I) where {I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}

    if sol.adaptive
        for it_gam = 1:100 # TODO: replace/complement with lower bound on gamma
            mul!(sol.Aqxbar, sol.Aq, sol.xbar)
            fq_Aqxbar = gradient!(sol.gradfq_Aqxbar, sol.fq, sol.Aqxbar)
            mul!(sol.Asxbar, sol.As, sol.xbar)
            fs_Asxbar = gradient!(sol.gradfs_Asxbar, sol.fs, sol.Asxbar)
            f_Axbar = fs_Asxbar + fq_Aqxbar

            uppbnd = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) +
                 0.5/sol.gamma*sol.normFPR_x^2
            if f_Axbar > uppbnd + 1e-6*abs(sol.f_Ax)
                sol.gamma = 0.5*sol.gamma
                sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax
                sol.g_xbar = prox!(sol.xbar, sol.g, sol.y, sol.gamma)
                sol.FPR_x .= sol.x .- sol.xbar
                sol.normFPR_x = norm(sol.FPR_x)
            else
                sol.FBE_x = uppbnd + sol.g_xbar
                break
            end
        end
    end

    if it > 1
        update!(sol.H, sol.x, sol.x_prev, sol.FPR_x, sol.FPR_x_prev)
    end
    mul!(sol.d, sol.H, ( x -> .-x ).(sol.FPR_x)) # TODO: not nice

    sol.FPR_x_prev, sol.FPR_x = sol.FPR_x, sol.FPR_x_prev
    sol.x_prev .= sol.x

    sigma = 0.5*sol.beta/sol.gamma*(1.0-sol.alpha)
    maxit_tau = 10

    # tau = 1
    sol.tau = one(R)

    mul!( sol.Aqd, sol.Aq, sol.d)
    mul!( sol.Asd, sol.As, sol.d)

    # xnew = x + tau*d
    sol.xnew .= sol.x .+ sol.tau .* sol.d
    # Aq*xnew = Aq*x + tau*Aq*d
    sol.Aqxnew .= sol.Aqx .+ sol.tau .* sol.Aqd
    # As*xnew = As*x + tau*As*d
    sol.Asxnew .= sol.Asx .+ sol.tau .* sol.Asd

    # calculate new FBE in xnew
    sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqxnew)
    sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asxnew)

    mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx)
    mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx)

    sol.At_gradf_Ax .= sol.Aqt_gradfq_Aqx .+ sol.Ast_gradfs_Asx
    sol.f_Ax = sol.fs_Asx + sol.fq_Aqx

    # gradient step
    sol.y .= sol.xnew .- sol.gamma .* sol.At_gradf_Ax
    # prox step
    sol.g_xbar = prox!(sol.xnewbar, sol.g, sol.y, sol.gamma)

    sol.FPR_xnew .= sol.xnew .- sol.xnewbar
    norm_FPRxnew = norm(sol.FPR_xnew)

    FBE_xnew = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_xnew)) +
                   0.5/sol.gamma*norm_FPRxnew^2 + sol.g_xbar

    if FBE_xnew > sol.FBE_x - sigma*sol.normFPR_x^2
        # start using convex combination of FB direction and d

        mul!(sol.Aqfb, sol.Aq, sol.FPR_x_prev)
        mul!(sol.Asfb, sol.As, sol.FPR_x_prev)

        for it_tau = 1:maxit_tau # TODO: replace/complement with lower bound on tau

            sol.tau *= 0.5

            # xnew = x + tau*d
            sol.xnew .= sol.x .+ sol.tau .* sol.d
            # xnew = x + tau*d - (1-tau)*fb
            sol.xnew .= sol.xnew .+ (sol.tau-1.0) .* sol.FPR_x_prev

            # Aq*xnew = Aq*x + tau*Aq*d
            sol.Aqxnew .= sol.Aqx .+ sol.tau .* sol.Aqd
            # Aq*xnew = Aq*x + tau*Aq*d - (1-tau)*Aq*fb
            sol.Aqxnew .= sol.Aqxnew .+ (sol.tau-1.0) .* sol.Aqfb

            # As*xnew = As*x + tau*As*d
            sol.Asxnew .= sol.Asx .+ sol.tau .* sol.Asd
            # As*xnew = As*x + tau*As*d - (1-tau)*As*fb
            sol.Asxnew .= sol.Asxnew .+ (sol.tau-1.0) .* sol.Asfb

            # calculate new FBE in xnew
            sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqxnew)
            sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asxnew)

            mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx)
            mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx)

            sol.At_gradf_Ax .= sol.Aqt_gradfq_Aqx .+ sol.Ast_gradfs_Asx
            sol.f_Ax = sol.fs_Asx + sol.fq_Aqx

            # gradient step
            sol.y .= sol.xnew .- sol.gamma .* sol.At_gradf_Ax
            # prox step
            sol.g_xbar = prox!(sol.xnewbar, sol.g, sol.y, sol.gamma)

            sol.FPR_xnew .= sol.xnew .- sol.xnewbar
            norm_FPRxnew = norm(sol.FPR_xnew)

            FBE_xnew = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_xnew)) +
                   0.5/sol.gamma*norm_FPRxnew^2 + sol.g_xbar

            if FBE_xnew <= sol.FBE_x - sigma*sol.normFPR_x^2
                break
            end

        end

    end

    sol.normFPR_x = norm_FPRxnew
    sol.FPR_x, sol.FPR_xnew = sol.FPR_xnew, sol.FPR_x
    sol.FBE_x = FBE_xnew
    sol.x, sol.xbar, sol.xnew, sol.xnewbar = sol.xnew, sol.xnewbar, sol.x, sol.xbar
    sol.Aqx, sol.Aqxnew = sol.Aqxnew, sol.Aqx
    sol.Asx, sol.Asxnew = sol.Asxnew, sol.Asx

    return sol.xbar

end

################################################################################
# Solver interface(s)

function PANOC(x0; kwargs...)
    sol = PANOCIterator(x0; kwargs...)
    it, point = run!(sol)
    x0 .= point
    return (it, point, sol)
end


text:
@testset "ZpnGModules" begin

  @testset "Minimal Submodules" begin

    F,a=Nemo.FiniteField(3,1,"a")
    R=ResidueRing(FlintZZ,9)

    V=DiagonalGroup([3,3,9,9])

    l=[1,1,3,0,2,1,3,3,1,1,1,1,0,0,0,1]
    l1=[1,1,1,0,2,1,1,1,0,0,1,1,0,0,0,1]
    A=MatrixSpace(R,4,4)(l)
    A1=MatrixSpace(F,4,4)(l1)

    M=ZpnGModule(V,[A])
    M1=FqGModule([A1])

    ls=minimal_submodules(M)
    ls1=minimal_submodules(M1)

    @test length(ls)==length(ls1)
    for x in ls
      @test Hecke.issubmodule(M,x)
    end
  end

  @testset "Dual Module" begin

    R=ResidueRing(FlintZZ,9)
    V=DiagonalGroup([3,3,9,9])
    V.issnf=true
    V.snf=[3,3,9,9]
    l=[1,1,3,0,2,1,3,3,1,1,1,1,0,0,0,1]
    A=MatrixSpace(R,4,4)(l)
    M=ZpnGModule(V,[A])
    N= Hecke.dual_module(M)
    ls=submodules(N)
    v=fmpz[3,3,1,1]
    for x in ls
      @test Hecke.issubmodule(M,Hecke._dualize(x,V,v))
    end

  end

  @testset "submodules with given structure" begin

    R=ResidueRing(FlintZZ,8)
    V=DiagonalGroup([2,4,8,8])
    V.issnf=true
    V.snf=[2,4,8,8]
    l=[1,2,4,0,1,1,0,2,1,1,1,1,0,2,0,1]
    l1=[1,0,0,0,0,3,4,2,1,0,0,1,0,0,1,0]
    l2=[1,0,0,4,1,1,0,0,0,2,1,0,1,1,1,1]
    A=MatrixSpace(R,4,4)(l)
    B=MatrixSpace(R,4,4)(l1)
    C=MatrixSpace(R,4,4)(l2)
    M=ZpnGModule(V,[A,B,C])
    ls=submodules(M,typesub=[2,3])
    y=subgroups(V,quotype=[4,8])

    mp1=Hecke.GrpAbFinGenMap(V,V,lift(A))
    mp2=Hecke.GrpAbFinGenMap(V,V,lift(B))
    mp3=Hecke.GrpAbFinGenMap(V,V,lift(C))
    act=[mp1,mp2,mp3]

    i=0
    for el in y
      if Hecke.is_stable(act,el[2])
        i+=1
      end
    end
    @test i==length(ls)

    ls=submodules(M,typesub=[3])
    y=subgroups(V,quotype=[8])
    i=0
    for el in y
      if Hecke.is_stable(act,el[2])
        i+=1
      end
    end
    @test i==length(ls)

  end

  @testset "submodules" begin

    R=ResidueRing(FlintZZ,4)
    V=DiagonalGroup([2,2,4])
    V.issnf=true
    V.snf=[2,2,4]
    A=MatrixSpace(R,3,3)(1)
    M=ZpnGModule(V,[A])
    ls=submodules(M)
    lsub=subgroups(V)
    @test length(collect(ls))==length(collect(lsub))

  end

end


text:
# Let's try parameter estimation for an ODE, here the famous Lorentz attractor.
# We base our implementation on the code from Paulo Marques available at:
#   https://github.com/pjpmarques/Julia-Modeling-the-World/blob/master/Lorenz%20Attractor.ipynb

function lorentz_ode!(dr::AbstractVector{Float64},
                      params::AbstractVector{Float64},
                      r::AbstractVector{Float64})
    @assert length(dr) == length(r) == 3

    # Get individual params for this ODE
    @inbounds sigma, rho, beta = params

    # Extract the coordinates from the r vector
    @inbounds x, y, z = r

    # The Lorenz equations, put the derivative into dr
    @inbounds dr[1] = sigma*(y - x)       # dx/dt
    @inbounds dr[2] = x*(rho - z) - y     # dy/dt
    @inbounds dr[3] = x*y - beta*z        # dz/dt

    return dr
end

# Define time vector and interval grid.
t = range(0.0, step=0.001, stop=100.0)
# But in order to compare to [Xiang2015] paper we use their choice instead:
t_Xiang2015 = range(0.0, step=0.01, length=301)

# Initial position in space
r0 = [0.1, 0.0, 0.0]

# Constants sigma, rho and beta. In a real ODE problem these would not be known
# and would be estimated.
sigma = 10.0
rho   = 28.0
beta  = 8.0/3.0
real_params = [sigma, rho, beta]

# "Play" an ODE from a starting point and into the future given a sequence of time steps.
# Put the restults into `states`
function calc_states!(states::AbstractMatrix{Float64},
    params::AbstractVector{Float64}, odefunc!::Function,
    startx::AbstractVector{Float64}, times::AbstractVector{Float64})

    numsamples = length(times)
    @assert size(states) == (length(startx), numsamples)

    states[:, 1] = startx
    tprev = times[1]
    deriv = similar(startx)
    for i in 2:numsamples
        tnow = times[i]
        stateprev = view(states, :, i-1)
        odefunc!(deriv, params, stateprev)
        @inbounds states[:, i] .= stateprev .+ deriv .* (tnow - tprev)
        tprev = tnow
    end

    return states
end

calc_states(params::AbstractVector{Float64}, odefunc!::Function,
    startx::AbstractVector{Float64}, times::AbstractVector{Float64}) =
    calc_states!(Matrix{Float64}(undef, length(startx), length(times)),
                 params, odefunc!, startx, times)

# RSS = Residual Sum of Squares, columnwise
rss(actual::AbstractMatrix{Float64}, estimated::AbstractMatrix{Float64}) =
    @inbounds sum(i -> abs2(actual[i] - estimated[i]), eachindex(actual, estimated))

# Calculate the actual/original state vectors that we will use for parameter
# estimation:
origstates = calc_states(real_params, lorentz_ode!, r0, t)
origstates_Xiang2015 = calc_states(real_params, lorentz_ode!, r0, t_Xiang2015)

function subsample(origstates::AbstractMatrix{Float64},
                   times::AbstractVector{Float64};
                   lenratio = 0.25)
    N = length(times)
    @assert size(origstates, 2) == N
    stopidx = round(Int, lenratio*N)
    indexes = 1:stopidx
    return origstates[:, indexes], times[indexes]
end

# The [Xiang2015] paper, https://www.hindawi.com/journals/ddns/2015/740721/,
# used these param bounds:
Xiang2015Bounds = [(9., 11.), (20., 30.), (2., 3.)]

# Now we can optimize using BlackBoxOptim
using BlackBoxOptim

# store temporary states for fitness calculation
const tmpstates_pool = Dict{NTuple{2, Int}, Vector{Matrix{Float64}}}()

# get the RSS between the origstates trajectory and ODE solution for given params
function lorentz_fitness(params::AbstractVector{Float64},
                         origstates::AbstractMatrix{Float64},
                         times::AbstractVector{Float64})
    # get the state matrix from the pool of the proper size
    statesize = size(origstates)
    states_pool = get!(() -> Vector{Matrix{Float64}}(), tmpstates_pool, statesize)
    states = !isempty(states_pool) ? pop!(states_pool) : Matrix{Float64}(undef, statesize...)
    # solve ODE for given params
    calc_states!(states, params, lorentz_ode!, r0, times)
    push!(states_pool, states) # return states matrix back to the pool
    return rss(origstates, states)
end

# But optimizing all states in each optimization step is too much so lets
# sample a small subset and use for first opt iteration.
origstates1, times1 = subsample(origstates, t; lenratio=0.04); # Sample only first 4%
origstates1_Xiang2015, times1_Xiang2015 = subsample(origstates_Xiang2015, t_Xiang2015; lenratio=1.00);

res1 = bboptimize(params -> lorentz_fitness(params, origstates1, times1);
    SearchRange = Xiang2015Bounds, MaxSteps = 8e3)

res2 = bboptimize(params -> lorentz_fitness(params, origstates1_Xiang2015, times1_Xiang2015);
    SearchRange = Xiang2015Bounds, MaxSteps = 11e3) # They allow 12k fitness evals for 3-param estimation

# But lets also try with less tight bounds
LooserBounds = [(0., 22.), (0., 60.), (1., 6.)]
res3 = bboptimize(params -> lorentz_fitness(params, origstates1_Xiang2015, times1_Xiang2015);
    SearchRange = LooserBounds, MaxSteps = 11e3) # They allow 12k fitness evals for 3-param estimation

@info "Results on the long time sequence from Paulo Marques:"
estfitness = lorentz_fitness(best_candidate(res1), origstates, t)
@show estfitness best_candidate(res1) best_fitness(res1)
origfitness = lorentz_fitness(real_params, origstates, t)
@show origfitness real_params

@info "Results on the short time sequence used in [Xiang2015] paper:"
estfitness = lorentz_fitness(best_candidate(res2), origstates_Xiang2015, t_Xiang2015)
@show estfitness best_candidate(res2) best_fitness(res2)
origfitness = lorentz_fitness(real_params, origstates_Xiang2015, t_Xiang2015)
@show origfitness real_params

@info "Results on the short time sequence used in [Xiang2015] paper, but with looser bounds:"
estfitness = lorentz_fitness(best_candidate(res3), origstates_Xiang2015, t_Xiang2015)
@show estfitness best_candidate(res3) best_fitness(res3)
origfitness = lorentz_fitness(real_params, origstates_Xiang2015, t_Xiang2015)
@show origfitness real_params


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "cdplayer_model.jl")
include(@relpath "cdplayer_specifications.jl")

S = cdplayer_model()
X0, options = cdplayer_specification()
options = Dict(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# partition with 1D blocks
partition = [i:i for i in 1:120]

# reachability algorithms
algorithm_reach = BFFPSV18(:δ => 1e-3, :vars => [1], :partition => partition,
                           :assume_sparse => true)
algorithm_check = BFFPSV18(:δ => 1e-3, :vars => [1, 2], :partition => partition,
                           :assume_sparse => true)

# check property (does not verify for the given initial states)
options[:mode] = "check"
solution = solve(problem, options; op=algorithm_check)
# @assert solution.satisfied

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options; op=algorithm_reach)

# project flowpipe to time and x₁
solution.options[:plot_vars] = [0, 1]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("cdplayer")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
coast(region=[110 140 20 35],
      proj=(name=:Albers, center=[125 20], parallels=[25 45]),
      frame=:ag,
      resolution=:low,
      area=250,
      land=:green,
      shore=:thinnest,
      fmt=:png, savefig="1")

coast(region=[-88 -70 18 24], proj=(name=:eqdc, center=[-79 21], parallels=[19 23]),
      frame=:ag, res=:intermediate, borders=(type=1,pen=("thick","red")), land=:green,
      shore=:thinnest, fmt=:png, savefig="2")

coast(region=[-130 -70 24 52], proj=(name=:lambertConic, center=[-100 35], parallels=[33 45]),
      frame=:ag, res=:low, borders=((type=1, pen=("thick","red")), (type=2, pen=("thinner",))),
      area=500, land=:tan, water=:blue, shore=(:thinnest,:white), fmt=:png, savefig="3")

coast(region=(-180,-20,0,90), proj=:poly, xaxis=(annot=30,grid=10), yaxis=(annot=10,grid=10),
      res=:crude, area=1000, land=:lightgray, shore=:thinnest, figsize=10, fmt=:png, savefig="4")

coast(region="0/-40/60/-10+r", proj=(name=:laea, center=[30,-30]), frame=:ag, res=:low,
      area=500, land=(pattern=10,dpi=300), shore=:thinnest, figsize=10, fmt=:png, savefig="5")
gmt("destroy")

coast(region=:g, proj=(name=:laea, center=[280,30]), frame=:g, res=:crude, area=1000,
      land=:navy, figsize=8, fmt=:png, savefig="6")

coast(region=(-30,30,60,72), proj=(name=:Stereographic, center=[0,90], paralles=60),
      frame=:a10g, res=:low, area=250, land=:royalblue, water=:seashell,
      figscale="1:30000000", fmt=:png, savefig="7")

coast(region="-25/59/70/72+r", proj=(name=:stereographic, center=(10,90)), frame=:a20g, res=:low,
      area=250, land=:darkbrown, shore=:thinnest, water=:lightgray, figsize=11, fmt=:png, savefig="8")

coast(region="100/-42/160/-8r", proj=(name=:stereographic, center=(130,-30)), frame=:ag, res=:low,
      area=500, land=:green, ocean=:lightblue, shore=:thinnest, figsize=10, fmt=:png, savefig="1")

coast(region=:g, proj="G4/52/230/90/60/180/60/60", xaxis=(annot=2,grid=2), yaxis=(annot=1,grid=1),
      rivers=:all, res=:intermediate, land=:lightbrown, ocean=:lightblue, shore=:thinnest, figsize=10,
      par=(:MAP_ANNOT_MIN_SPACING,0.65), fmt=:png, savefig="9")

coast(region=:g, proj=(name=:ortho, center=(-75,41)), frame=:g, res=:crude, area=5000,
      land=:pink, ocean=:thistle, figsize=10, fmt=:png, savefig="10")

coast(region=:g, proj=(name=:azimuthalEquidistant, center=(-100,40)), frame=:g,
      res=:crude, area=10000, land=:lightgray, shore=:thinnest, figsize=10, fmt=:png, savefig="11")

coast(region=:g, proj=(name=:Gnomonic, center=(-120,35), horizon=60),
      frame=(annot=30, grid=15), res=:crude, area=10000, land=:tan, ocean=:cyan,
      shore=:thinnest, figsize=10, fmt=:png, savefig="12")

coast(region=(0,360,-70,70), proj=:Mercator, xaxis=(annot=60,ticks=15), yaxis=(annot=30,ticks=15),
      res=:crude, area=:5000, land=:red, scale=0.03, par=(:MAP_FRAME_TYPE,"fancy+"), fmt=:png, savefig="13")

coast(region="20/30/50/45r", proj=(name=:tmerc, center=35), frame=:ag, res=:low,
      area=250, land=:lightbrown, ocean=:seashell, shore=:thinnest, scale=0.45, fmt=:png, savefig="14")

coast(region=(0,360,-80,80), proj=(name=:tmerc, center=[330 -45]),
      frame=(annot=30, grid=:auto, axes=:WSne), res=:crude, area=2000, land=:black,
      water=:lightblue, figsize=9, fmt=:png, savefig="15")

coast(region="270/20/305/25+r", proj=(name=:omercp, center=[280 25.5], parallels=[22 69]),
      frame=:ag, res=:i, area=250, shore=:thinnest, land=:burlywood, water=:azure,
      rose="jTR+w1+f2+l+o0.4", figsize=12, par=(FONT_TITLE=8, MAP_TITLE_OFFSET=0.12), fmt=:png, savefig="16")

coast(region="7:30/38:30/10:30/41:30r", proj=(name=:Cassini, center=[8.75 40]),
      frame=:afg, map_scale="jBR+c40+w100+f+o0.4/0.5", land=:springgreen,
      res=:high, water=:azure, shore=:thinnest, rivers=(:all,:thinner), figsize=6,
      par=(:FONT_LABEL,12), fmt=:png, savefig="17")

coast(region=:g, proj=:equidistCylindrical, frame=(annot=60, ticks=30, grid=30),
      res=:crude, area=5000, land=:tan4, water=:lightcyan, figsize=12, fmt=:png, savefig="18")

coast(region=(-145,215,-90,90), proj=(name=:cylindricalEqualArea, center=(35,30)),
      frame=(annot=45, grid=45), res=:crude, area=10000, water=:dodgerblue,
      shore=:thinnest, figsize=12, fmt=:png, savefig="19")

coast(region=(-90,270,-80,90), proj=:Miller, xaxis=(annot=45,grid=45),
      yaxis=(annot=30,grid=30), res=:crude, area=10000, land=:khaki, water=:azure,
      shore=:thinnest, scale="1:400000000", fmt=:png, savefig="20")

coast(region=(-180,180,-60,80), proj=(name=:cylindricalStereographic, center=(0,45)),
      xaxis=(annot=60,ticks=30, grid=30), yaxis=(annot=30,grid=30), res=:crude,
      area=5000, shore=:black, land=:seashell4, ocean=:antiquewhite1, figsize=12, fmt=:png, savefig="21")

coast(region=:g, proj=:Hammer, frame=:g, res=:crude, area=10000, land=:black,
	  ocean=:cornsilk, figsize=12, fmt=:png, savefig="22")

coast(region=:d, proj=:Mollweide, frame=:g, res=:crude, area=10000, land=:tomato1,
	  water=:skyblue, figsize=12, fmt=:png, savefig="23")

coast(region=:d, proj=:Winkel, frame=:g, res=:crude, area=10000, land=:burlywood4,
	  water=:wheat1, figsize=12, fmt=:png, savefig="24")

coast(region=:d, proj=:Robinson, frame=:g, res=:crude, area=10000, land=:goldenrod,
	  water=:snow2, figsize=12, fmt=:png, savefig="25")

coast(region=:d, proj=:EckertIV, frame=:g, res=:crude, area=10000, land=:ivory,
	  water=:bisque3, shore=:thinnest, figsize=12, fmt=:png, savefig="26")

coast(region=:d, proj=:EckertVI, frame=:g, res=:crude, area=10000, land=:ivory,
	  water=:bisque3, shore=:thinnest, figsize=12, fmt=:png, savefig="27")

coast(region=:d, proj=:Sinusoidal, xaxis=(grid=30,), yaxis=(grid=15,), res=:crude,
      area=10000, land=:coral4, water=:azure3, figsize=12, fmt=:png, savefig="28")

coast(region=(200,340,-90,90), proj=:Sinusoidal, frame=:g, res=:crude, area=10000,
      land=:darkred, water=:azure, scale=0.03333)
coast!(region=(-20,60,-90,90), frame=:g, res=:crude, area=10000, land=:darkgreen,
       water=:azure, xoff=4.666)
coast!(region=(60,200,-90,90), frame=:g, res=:crude, area=10000, land=:darkblue,
       water=:azure, xoff=2.6664, fmt=:png, savefig="29")

coast(region=:g, proj=:VanderGrinten, xaxis=(grid=30,), yaxis=(grid=15,),res=:crude,
      land=:lightgray, water=:cornsilk, area=10000, shore=:thinnest, figsize=10, fmt=:png, savefig="30")


text:
counts = Dict{AbstractString, UInt64}()

for line in eachline(stdin)
  for m in eachmatch(r"[a-z\']+", lowercase(line))
    counts[m.match] = get(counts, m.match, 0) + 1
  end
end
for (word, count) in sort(collect(counts))
  println("$word $count")
end

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


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

@traitdef IsGrid{D}
@traitimpl IsGrid{D} <- isgrid(D)

"""
    isgrid(D)

Tells whether or not a domain or data of type `D` behaves as a grid
(or array) of elements. This information can be used to optimize
algorithms (e.g. slices, bounding boxes).
"""
isgrid(D) = false


text:
using DimensionalData
using PyCall
using IntervalSets
using HDF5

function downloadlfp(S::AbstractSession, probeid::Int)
    @assert(any(getprobeids(S) .== probeid), "Probe $probeid does not belong to session $(getid(S))")
    @assert(subset(getprobes(S), :id=>ByRow(==(probeid)))[!, :has_lfp_data][1], @error "Probe $probeid does not have LFP data")
    _ = S.pyObject.get_lfp(probeid)
    return nothing
end

function getlfppath(session::AbstractSession, probeid)
    path = joinpath(datadir, "Ecephys", "session_"*string(getid(session)), "probe_"*string(probeid)*"_lfp.nwb");
end

function getlfptimes(session::AbstractSession, probeid)
    path = getlfppath(session, probeid)
    if !isfile(path)
        @error "Probe lfp data file does not exist. This may be becuase you have not downloaded the probe data. Try using `getlfp`"
    end

    f = h5open(path)
    times = f["acquisition"][splitext(basename(path))[1]][splitext(basename(path))[1]*"_data"]["timestamps"][:]
    close(f)
    return times
end

function getlfpchannels(session::AbstractSession, probeid)
    if !isfile(getlfppath(session, probeid))
        downloadlfp(session, probeid)
    end
    f = h5open(getlfppath(session, probeid))
    channels = f["general"]["extracellular_ephys"]["electrodes"]["id"][:]
    close(f)
    return channels
end

"""
Get the lfp data for a probe,  providing *indices* for channels and times. See function below for indexing by channel ids and time values/intervals
"""
function _getlfp(session::AbstractSession, probeid::Int; channelidxs=1:length(getlfpchannels(session, probeid)), timeidxs=1:length(getlfptimes(session, probeid)))
    @assert(any(getprobeids(session) .== probeid), "Probe $probeid does not belong to session $(getid(session))")
    @assert(subset(getprobes(session), :id=>ByRow(==(probeid)))[!, :has_lfp_data][1], @error "Probe $probeid does not have LFP data")
    path = getlfppath(session, probeid)
    if !isfile(path)
        downloadlfp(session, probeid)
    end
    if !isfile(path)
        @error "LFP data did not download correctly"
    end

    f = h5open(path)

    timedata = getlfptimes(session, probeid)
    timedata = timedata[timeidxs]
    dopermute = true
    channelids = getlfpchannels(session, probeid)
    channelids = channelids[channelidxs]
    if (channelidxs isa Union{Int64, AbstractRange{Int64}}) & (timeidxs isa Union{Int64, AbstractRange{Int64}}) # Can use HDF5 slicing
        lfp = f["acquisition"][splitext(basename(path))[1]][splitext(basename(path))[1]*"_data"]["data"][channelidxs, timeidxs]
    elseif timeidxs isa Union{Int64, AbstractRange{Int64}}
        lfp = [f["acquisition"][splitext(basename(path))[1]][splitext(basename(path))[1]*"_data"]["data"][i, timeidxs] for i ∈ channelidxs]
        lfp = hcat(lfp...)
        dopermute = false
    else
        lfp = read(f["acquisition"][splitext(basename(path))[1]][splitext(basename(path))[1]*"_data"]["data"])
        lfp = lfp[channelidxs, timeidxs]
    end
    if lfp isa Vector
       lfp = reshape(lfp, 1, length(lfp))
    end
    if channelids isa Number
        channelids = [channelids]
    end
    if dopermute
        lfp = permutedims(lfp, reverse(1:ndims(lfp)))
    end
    X = DimArray(lfp, (Ti(timedata),  Dim{:channel}(channelids)))
    close(f)
    return X
end

function isinvalidtime(session, probeids=getprobeids(session), times=NaN)
    if isempty(session.pyObject.get_invalid_times()) # No invalid times in this session!
        return false
    end
    intervals = [session.pyObject.get_invalid_times().start_time.values, session.pyObject.get_invalid_times().stop_time.values]
    intervals = [[intervals[1][i], intervals[2][i]] for i ∈ 1:length(intervals[1])]
    tags = session.pyObject.get_invalid_times().tags.values
    tags = vcat([tryparse.(Int64, i) for i ∈ tags]...)
    badprobes = tags[.!isnothing.(tags)]
    if times isa Interval
        isininterval = [any(i .∈ (times,)) for i ∈ intervals]
    else
        isininterval = [any((times .> i[1]) .& (times .< i[2])) for i ∈ intervals]
    end
    return any(probeids .∈ (badprobes,)) & any(isininterval)
end

"""
This is the one you should be using. Get lfp data by channel id and time intervals or vector. Also, throw error if you try to access an invalid time interval.
"""
function getlfp(session::AbstractSession, probeid::Int; channels=getlfpchannels(session, probeid), times=ClosedInterval(extrema(getlfptimes(session, probeid))...), inbrain=false)
    if isinvalidtime(session, probeid, times)
        @error "Requested LFP data contains an invalid time..."
    end

    if inbrain isa Symbol || inbrain isa Real || inbrain
        depths = getchanneldepths(session, probeid, channels)
        if inbrain isa Real # A depth cutoff
            channels = channels[depths .> inbrain]
        elseif inbrain isa Symbol # A mode
        else # Just cutoff at the surface
            channels = channels[depths .> 0]
        end
    end

    channelidxs = getlfpchannels(session, probeid)
    channelidxs = indexin(channels, channelidxs)
    channelidxs = filter(!isnothing, channelidxs)
    @assert length(channels) == length(channelidxs)

    timeidxs = getlfptimes(session, probeid)
    if !(times isa Interval) && length(times) == 2
        times = ClosedInterval(times...)
    end
    if times isa Interval
        timeidxs = findall(timeidxs .∈ (times,))
    else
        timeidxs = indexin(times, timeidxs)
        timeidxs = filter(!isnothing, timeidxs)
        @assert length(timeidxs) == length(times)
    end

    # See if we can convert to unitranges for faster HDF5 reading via slicing
    if collect(UnitRange(extrema(timeidxs)...)) == timeidxs
        timeidxs = UnitRange(extrema(timeidxs)...)
    end
    if collect(UnitRange(extrema(channelidxs)...)) == channelidxs
        channelidxs = UnitRange(extrema(channelidxs)...)
    end

    @info "Accessing LFP data"
    _getlfp(session, probeid; channelidxs, timeidxs)
end

"""
If you want to downsample the LFP data, its quicker to use this function and then perform slicing afterwards (since getlfp() has to check all of the time coordinates you supply, which can be slow).
"""
function getdownsampledlfp(session, probeid; downsample=100, timerange=ClosedInterval(extrema(getlfptimes(session, probeid))...), channels=getlfpchannels(session, probeid))
    if !(timerange isa Interval) && length(timerange) == 2
        timerange = ClosedInterval(timerange...)
    end
    timevals = getlfptimes(session, probeid)
    tidxs = timevals .∈ (timerange,)
    times = findfirst(tidxs):downsample:findlast(tidxs)
    _getlfp(session, probeid; timeidxs = times)[:, At(channels)]
end

"""
Now we can overload `getlfp()` to index by structure
"""
function getlfp(session::AbstractSession, probeid::Int, structures::Union{Vector{String}, String}; kwargs...)
    if structures isa String
        structures = [structures]
    end
    channels = subset(getchannels(session, probeid), :ecephys_structure_acronym=>ByRow(∈(structures)), skipmissing=true)
    channels = channels.id ∩ getlfpchannels(session, probeid)
    getlfp(session, probeid; channels, kwargs...)
end

function getlfp(session, probeids::Vector{Int}, args...; kwargs...)
    LFP = [getlfp(session, probeid, args...; kwargs...) for probeid ∈ probeids]
end

function getchannels(data::AbstractDimArray)
    dims(data, :channel).val
end

"""
At the moment this is just a proxy: distance along the probe to the cortical surface
"""
function getchanneldepths(session, probeid, channels)
    cdf = getchannels(session, probeid)
    return _getchanneldepths(cdf, channels)
end
# function getchanneldepths(session, channels)
#     cdf = getchannels(session) # Slightly slower than the above
#     #cdf = cdf[indexin(channels, cdf.id), :]
#     cdfs = groupby(cdf, :probe_id)
#     depths = vcat([_getchanneldepths(c, c.id) for c ∈ cdfs]...)
#     depths = depths[indexin(channels, vcat(cdfs...).id)]
# end
function _getchanneldepths(cdf, channels)
    surfaceposition = minimum(subset(cdf, :ecephys_structure_acronym=>ByRow(ismissing)).probe_vertical_position)
    # Assume the first `missing` channel corresponds to the surfaceprobe_vertical_position
    idxs = indexin(channels, cdf.id)[:]
    alldepths = surfaceposition .- cdf.probe_vertical_position # in μm
    depths = fill(NaN, size(idxs))
    depths[.!isnothing.(idxs)] = alldepths[idxs[.!isnothing.(idxs)]]
    return depths
end

getdim(X::AbstractDimArray, dim) = dims(X, dim).val
gettimes(X::AbstractDimArray) = getdim(X, Ti)

function sortbydepth(session, probeid, LFP::AbstractDimArray)
    depths = getchanneldepths(session, probeid, getchannels(LFP))
    indices = Array{Any}([1:size(LFP, i) for i in 1:length(size(LFP))])
    indices[findfirst(isa.(dims(LFP), Dim{:channel}))] = sortperm(depths)
    return LFP[indices...]
end


text:
"""
    BipartiteView{T<:Real} <: Graphs.SimpleGraphs.AbstractSimpleGraph{Int}

Represents a bipartite view of a hypergraph `h`.
Note this is a view - changes to the original hypergraph will be automatically reflected in the view.

**Constructors**

BipartiteView(::Hypergraph)

The bipartite view of a hypergraph is suitable for processing with the Graphs.jl package.
Several Graphs methods are provided for the compability.

"""
struct BipartiteView{T<:Real} <: Graphs.SimpleGraphs.AbstractSimpleGraph{Int}
    h::Hypergraph{T}
end

"""
  Return the number of vertices in a bipartite view `b` of a hypergraph.
"""
Graphs.nv(b::BipartiteView) = length(b.h.v2he)+length(b.h.he2v)

Graphs.vertices(b::BipartiteView) = Base.OneTo(Graphs.nv(b))

"""
  Return the number of edges in a bipartite view `b` of a hypergraph.
"""
Graphs.ne(b::BipartiteView) = sum(length.(b.h.v2he))

function Graphs.all_neighbors(b::BipartiteView, v::Integer)
    n1 = length(b.h.v2he)
    if v <= n1
      n1 .+ keys(b.h.v2he[v])
    else
      collect(keys(b.h.he2v[v-n1]))
    end
end

function Graphs.has_edge(b::BipartiteView, s, d)
    n1 = length(b.h.v2he)
    if s <= n1
        d > n1 && has_key(b.v2he[s], d - n1)
    else
        d <= n1 && has_key(b.h.he2v[s-n1], d)
    end
end

Graphs.has_vertex(b::BipartiteView, v::Integer) = 1 <= v <= Graphs.nv(b)

Graphs.outneighbors(b::BipartiteView, v::Integer) = Graphs.all_neighbors(b::BipartiteView, v)

Graphs.inneighbors(b::BipartiteView, v::Integer) = Graphs.all_neighbors(b::BipartiteView, v)

"""
    Graphs.SimpleGraph(b::BipartiteView)

Creates a `Graphs.SimpleGraph` representation of a `BipartiteView` b.

This creates a copy of the date. Note that the weights information is not stored
in the created `SimpleGraph`.
"""
function Graphs.SimpleGraph(b::BipartiteView)
    g = SimpleGraph(nv(b))
    for v in keys(b.h.v2he)
        for he in keys(b.h.v2he[v])
            add_edge!(g, v, length(b.h.v2he) + he)
        end
    end
    g
end

Graphs.is_directed(b::BipartiteView{T}) where T = false

Graphs.is_directed(::Type{BipartiteView{T}}) where T = false

Base.eltype(::BipartiteView{T}) where T = Int

"""
    shortest_path(b::BipartiteView,source::Int, target::Int)

Finds a single shortest path in a graph `b` between vertices
`source` and `target`.
Note that if several paths of the same length exist, only one
will be returned.

"""
function shortest_path(b::BipartiteView,source::Int, target::Int)
    checkbounds(b.h.v2he, source)
    checkbounds(b.h.v2he, target)
    dj = dijkstra_shortest_paths(b, source)
    enumerate_paths(dj)[target][1:2:end]
end

"""
    Graphs.SimpleGraphs.fadj(b::BipartiteView)

Generates an adjency list for this view of a hypergraph.
"""
function Graphs.SimpleGraphs.fadj(b::BipartiteView)
    res = Vector{Vector{Int}}(undef, Graphs.nv(b))

    h_nv = length(b.h.v2he)
    for i in 1:h_nv
       res[i] = h_nv .+ sort!(collect(keys(b.h.v2he[i])))
    end
    for i in 1:length(b.h.he2v)
        res[i+h_nv] = sort!(collect(keys(b.h.he2v[i])))
    end
    res
end

Graphs.SimpleGraphs.fadj(b::BipartiteView, v::Integer) = Graphs.all_neighbors(b,v)
Graphs.edges(b::BipartiteView) = Graphs.SimpleGraphs.SimpleEdgeIter(b)

Graphs.edgetype(b::BipartiteView{T}) where T = Graphs.SimpleGraphs.SimpleEdge{Int}

Graphs.zero(t::BipartiteView{T}) where T = BipartiteView(Hypergraph{T}(0,0))
Graphs.zero(::Type{BipartiteView{T}}) where T = BipartiteView(Hypergraph{T}(0,0))


text:
using ArcadeLearningEnvironment
using IICGP
using Test

ROM_NAMES = ["boxing", "assault", "freeway"]

@testset "Reset ALE" begin
    for rom_name in ROM_NAMES
        for _ in 1:10
            seed = 0
            game = Game(rom_name, seed)
            s0 = get_state(game, true, false)
            for t in 1:100
                action = 0
                act(game.ale, action)
            end
            reset!(game)
            s1 = get_state(game, true, false)
            @test s0 == s1
            #display(IICGP.implot(s0[1]))
            #display(IICGP.implot(s1[1]))
            close!(game)
        end
    end
end


text:
#!/usr/bin/env julia

# uses too many queries for testset 3

function solve(N, Q)
    function ask(query)
        Q < 0 && exit()
        println(join(query, ' '))
        res = parse(Int, readline())
        Q -= 1
        res == -1 && exit()
        res
    end

    # bootstrap so we have 3 to begin with
    query = [1, 2, 3]
    mid = ask(query)
    popat!(query, mid)
    S = [query[1], mid, query[2]]

    function bisect(val, left, right)
        left == right && return insert!(S, left, val)

        if right - left == 1
            mid = ask((S[left], val, S[right]))
            mid == S[left] && return insert!(S, left, val)
            mid == S[right] && return insert!(S, right + 1, val)
            return insert!(S, right, val)
        end

        mid = (left + right) ÷ 2
        #@info val left S[left] right S[right] mid S[mid]

        lefttomid = ask((S[left], val, S[mid]))
        lefttomid == S[left] && return insert!(S, left, val)
        lefttomid == val && return bisect(val, left, mid)
        midtoright = ask((S[mid], val, S[right]))
        midtoright == S[right]  && return insert!(S, right + 1, val)
        midtoright == mid && insert!(S, mid, val)
        bisect(val, mid + 1, right)
    end

    for val in 4:N
        bisect(val, firstindex(S), lastindex(S))
        #@info S
    end

    @assert ask(S) == 1
end

function main()
    T, N, Q = map(x -> parse(Int, x), split(readline()))
    for case in 1:T
        solve(N, Q)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
@everywhere begin

function checkwid(x...)
    @assert myid() == 2
    return 1
end

end

@testset "Scheduler options" begin
    @testset "single worker" begin
        a = delayed(checkwid)(1)
        b = delayed(checkwid)(2)
        c = delayed(checkwid)(a,b)

        @test collect(Context(), c; single=2) == 1
    end
end