text:
len = length

using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
filename = "nhept33.nc"
data = ncdata.get(filename)

specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"]),[Symbol(i)for i in data["rc"]])

rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]

rates = rates[7:len(rates)]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]

# reading in data and cleaning
joint = "\n"*(join(names(rates),"\n"))*"\n"

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

#df = DataFrame(source = Char[], target = Char[], weight = Float16[])

flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  #println(i)
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
      push!(edges,[match( r"([\d\s\.]*)(\D[\d\D]*)", j)[2],match( r"([\d\s\.]*)(\D[\d\D]*)", l)[2],i])
    end

  end
end

normalise(x) = x/norm(x)

f = open("mclsaved.out","w")

for t in 1:len(specs[1])

  links = filter(i -> flux[i][t]>0 , 1:len(flux))
  tflux = [log10(flux[i][t]) for i in links]
  weight = 1+normalise(tflux)

  source = [i[1] for i in edges]
  target=[i[2] for i in edges]
  weighted = [weight[i[3]] for i in edges]

  @rput source
  @rput target
  @rput weighted

  R"library(igraph)"

  R"el <- structure(list(V1 = source, V2 = target, weight = weighted), .Names = c('V1',
  'V2', 'weight'), class = 'data.frame', row.names = c(1:$(length(edges))
  ))"

  R"g <- graph.data.frame(el)"

  R"net = g"

  label = R"a = cluster_spinglass(g, spins=50)"

  print(t)

  for i in 1:length(label)
    R"str = paste(unlist(a[$i]),collapse='\t')"
    @rget str
    if str != ""
      write(f,str)
      write(f,",")
    end
  end
  write(f,"\n")
end
close(f)

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
@discord_object struct ActivityEmoji
    name::String
    id::Snowflake
    animated::Bool
end

ext: jl


text:
#using Pkg
#Pkg.add("DrWatson")

# Load DrWatson (scientific project manager)
using DrWatson

#using SharedArrays

# Load common code to all script in ICNBenchmarks
using ICNBenchmarks

# Activate the ICNBenchmarks project
@quickactivate "ICNBenchmarks"

# Pkg.instantiate()

using Base:
    IOError, UV_EEXIST, UV_ESRCH,
    Process

using Base.Filesystem:
    File, open, JL_O_CREAT, JL_O_RDWR, JL_O_RDONLY, JL_O_EXCL,
    samefile

using Random

# Helper function to open files in O_EXCL mode
function tryopen_exclusive(path::String, mode::Integer = 0o666)
    try
        return open(path, JL_O_RDWR | JL_O_CREAT | JL_O_EXCL, mode)
    catch ex
        (isa(ex, IOError) && ex.code == UV_EEXIST) || rethrow(ex)
    end
    return nothing
end


text:
using Base.Test
using Astrodynamics

mu = 398600.448073446
rv = [8.59072560e+02, -4.13720368e+03, 5.29556871e+03, 7.37289205e+00, 2.08223573e+00, 4.39999794e-01]
rvmat = [rv';rv';rv']
ele = [6.78085866e+03,   1.30546073e-03,   9.00610999e-01, 3.46237544e+00,   6.85258061e-01,   8.13258968e-01]
elemat = [ele';ele';ele']
@test_approx_eq ele rvtokepler(rv, mu)
@test_approx_eq elemat rvtokepler(rvmat, mu)
@test_approx_eq rv keplertorv(ele, mu)
@test_approx_eq rvmat keplertorv(elemat, mu)


text:
"""
    ExprParsers = EP

A module for parsing common `Base.Expr` patterns.

The main interface encompass just three concepts, which seamlessly interact with oneanother
- macro `@exprparser`: easily create definitions for highly flexible and nestable parsers
- function `parse_expr`: compares a parser with a value, and returns a parsed value
- function `to_expr`: transforms parsed values back to `Base.Expr`
"""
module ExprParsers
export EP, @passert, parse_expr, to_expr

"""
    EP = ExprParsers

Short alias for the `ExprParsers` package.
"""
const EP = ExprParsers

using Compat

include("Utils.jl")
using .Utils
include("exceptions.jl")
include("syntax.jl")

# Main Interface
# ==============

"""
    parse_expr(parser, value)

Match parser against a value, will throw ParseError if the parser does not match.

Defaults to comparing with `==`, if matches, will return the value.
Parsers will be called instead and return their parsed value.

# Examples
```jldoctest
julia> using ExprParsers

julia> parse_expr(:a, :a)
:a
julia> parse_expr([1,2,3], [1,2,3])
3-element Array{Int64,1}:
 1
 2
 3
julia> parse_expr([1,2,3], [1,2,35])
ERROR: ParseError: Using default `==` comparison, but parser `3` ≠ value `35`.
julia> parse_expr(:(a = 4), :(a = 5))
ERROR: ParseError: Using default `==` comparison, but parser `4` ≠ value `5`.
julia> parse_expr([1,2,3,4], [1,2])
ERROR: ParseError: length(parser) == length(values) = false
  length(parser) = 4
  parser = [1, 2, 3, 4]
  length(values) = 2
  values = [1, 2]
```
"""
function parse_expr(parser, value)
  @passert parser == value "Using default `==` comparison, but parser `$(repr(parser))` ≠ value `$(repr(value))`."
  value
end
# vectors are treated especially
function parse_expr(parser::Union{Vector, Tuple, Iterator}, values::Union{Vector, Tuple})
  if !isa(Base.IteratorSize(parser), Base.IsInfinite)
    @passert length(parser) == length(values)
  end
  map(parser, values) do p, v
    parse_expr(p, v)
  end
end
# also Base.Expr are treated especially
function parse_expr(parser::Base.Expr, values::Base.Expr)
  @passert parser.head == values.head
  Base.Expr(parser.head, parse_expr(parser.args, values.args)...)
end

"""
    to_expr(parsed)

Converts parsed information back to Expr.

Defaults to returning same value, however if something knows about how it can be translated back,
just overload the function.
"""
to_expr(a) = a
# reverses also Iterator, as parse_expr on Iterator will return Vector
to_expr(a::Union{Vector, Tuple}) = map(to_expr, a)
to_expr(a::Base.Expr) = Base.Expr(a.head, to_expr(a.args)...)

"""
    EP.ExprParser

All parsers in the `ExprParsers` package inherit from this type.
"""
abstract type ExprParser end

# we support broadcasting, a ExprParser is always like a single value
Base.length(::ExprParser) = 1
Base.iterate(parser::ExprParser) = parser, nothing
Base.iterate(::ExprParser, ::Nothing) = nothing

"""
    EP.ExprParserWithParsed

This Type is mainly for internal usage. Please use [`@exprparser`](@ref) instead for the public interface.

It is a subtype of [`ExprParser`](@ref) which indicates that this parser actually constructs a [`ExprParsed`](@ref)
object when calling `parse_expr(parser, expr)`.

The resulting `ExprParsed` object is a struct with identical fields like the parser,
where then the parsed values will be stored.

`ExprParsed(parser::ExprParserWithParsed)` will return the corresponding `ExprParsed` type.
"""
abstract type ExprParserWithParsed <: ExprParser end

"""
    EP.ExprParsed(ParserType::Type{<:ExprParserWithParsed}) -> Type{<:ExprParsed}

Maps Parser Type to respective Parsed Type, and is also abstract super type of all Parsed types.

Example
```
EP.ExprParsed(EP.Assignment) == EP.Assignment_Parsed
```
"""
abstract type ExprParsed end

Base.convert(::Base.Type{Base.Expr}, parsed::ExprParsed) = to_expr(parsed)

# Parsers
# =======

include("expr_parsers_core.jl")
include("expr_parsers_with_parsed.jl")
include("expr_parsers_meta.jl")

end # module


text:
using Colors
mutable struct HeatWaveEffect{T<:AbstractFloat} <: ConfigurableEffect
    config::EffectConfig{T}
    colormap::Array{UInt8, 2}
    colorgrad::Array{HSL{T}}
    num_colors::T
    spec::Array{T}
    maxspec::T
    function HeatWaveEffect{T}(leddata::LEDArray, config::EffectConfig{T}, ana::AudioAnalysis) where {T<:AbstractFloat}
        init_spec, maxspec = binFFT(ana, 25)
        num_colors = 100*config.scaling
        init_colorgrad = linspace(config.primary_color, config.secondary_color, num_colors)
        init_colormap = cat(1, [hsl_to_rgb(c.h/360, c.s, c.l) for c in init_colorgrad]'...)
        return new(config, init_colormap, init_colorgrad, num_colors, init_spec, maxspec)
    end
end
HeatWaveEffect(leddata::LEDArray, config::EffectConfig{T}, ana::AudioAnalysis) where {T<:AbstractFloat} = HeatWaveEffect{T}(leddata, config, ana)
function update_vars!(leddata::LEDArray, ana::AudioAnalysis, effect::HeatWaveEffect)
    effect.num_colors = effect.config.scaling * 100
    effect.colorgrad = linspace(effect.config.primary_color, effect.config.secondary_color, effect.num_colors)
    effect.colormap = cat(1, [hsl_to_rgb(c.h/360, c.s, c.l) for c in effect.colorgrad]'...)
    effect.spec, effect.maxspec = binFFT(ana, 25)
end
function update!(leddata::LEDArray, ana::AudioAnalysis, effect::HeatWaveEffect)
    update_vars!(leddata, ana, effect)
    color_idx = @bounded(round(Int, mean(effect.spec[1:5])*size(effect.colormap,1)), 1, size(effect.colormap,1) )
    for channel in leddata.channels
        for i in 1:size(channel, 1)
            channel[i, :] = effect.colormap[color_idx, :]
        end
        update!(channel)
    end
end


text:
# ---
# title: 845. Longest Mountain in Array
# id: problem845
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers
# link: <https://leetcode.com/problems/longest-mountain-in-array/description/>
# hidden: true
# ---
#
# You may recall that an array `arr` is a **mountain array** if and only if:
#
#   * `arr.length >= 3`
#   * There exists some index `i` ( **0-indexed** ) with `0 < i < arr.length - 1` such that:
#     * `arr[0] < arr[1] < ... < arr[i - 1] < arr[i]`
#     * `arr[i] > arr[i + 1] > ... > arr[arr.length - 1]`
#
# Given an integer array `arr`, return _the length of the longest subarray,
# which is a mountain_. Return `0` if there is no mountain subarray.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [2,1,4,7,3,2,5]
#     Output: 5
#     Explanation: The largest mountain is [1,4,7,3,2] which has length 5.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [2,2,2]
#     Output: 0
#     Explanation: There is no mountain.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 104`
#   * `0 <= arr[i] <= 104`
#
#
#
# **Follow up:**
#
#   * Can you solve it using only one pass?
#   * Can you solve it in `O(1)` space?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "131.palindrome-partitioning.jl" begin
    @test partition_ps("aab") == [["a", "a", "b"], ["aa", "b"]]
    @test partition_ps("a") == [["a"]]
end


text:
# basic parameter and prediction tests

using DelimitedFiles
using JLD2
using TensorCast
using Flux
using RoME

using Gadfly
Gadfly.set_default_plot_size(35cm,20cm)

include(joinpath((@__DIR__), "LoadPyNNTxt.jl"))

allModels = []
for i in 0:99
# /home/dehann/data/racecar/results/conductor/models/retrained_network_weights0
  # push!(allModels, loadPose2OdoNNModelIntoFlux(ENV["HOME"]*"/data/racecar/results/conductor/models/retrained_network_weights$i") )
  push!(allModels, loadPose2OdoNNModelIntoFlux(ENV["HOME"]*"/data/racecar/results/conductor/sim_models/sim_network_weights$i") )
end

allModels[1]

# load sim data

# thestring = joinpath(ENV["HOME"],"Documents","learning","sim5data.jld2") |> string
@load "/home/dehann/Documents/learning/sim5data.jld2"

sim5Inputs

# check consistency of all input data
for i in 2:4, j in 1:50
  @assert sim5Inputs[1][j] - sim5Inputs[i][j] |> norm < 1e-10 "input data inconsistent"
end

## check on GPU

# allModelsGPU = deepcopy(allModels)

# loadPose2OdoNNModelIntoFluxGPU
# sim1_gpu = Float32.(sim5Inputs[1][1]) |> gpu
# fnc = allModelsGPU[1]
# @btime fnc(sim1_gpu)

## compare to CPU

# sim1_cpu = Float32.(sim5Inputs[1][1])
# fnc = allModels[1]
# @btime fnc($sim1_cpu)

# turns out much faster while the scalar slow down on GPU is not fixed
##

modelnum = 4
pred5 = (allModels[modelnum]).(sim5Inputs[modelnum])

@cast out5new[row,col] := pred5[1,row][col]
pln = Gadfly.plot(x=out5new[:,1],y=out5new[:,2],Geom.path)

@cast out5py[row,col] := sim5Outs[$modelnum][row][1,col]
plp = Gadfly.plot(x=out5py[:,1],y=out5py[:,2],Geom.path)

hstack(pln, plp)

# sim5Outs[1]

#


text:
"""
    _Model{V <: Variable{<:AbstractDomain},C <: Constraint{<:Function},O <: Objective{<:Function}}
A struct to model a problem as a set of variables, domains, constraints, and objectives.
```
struct _Model{V <: Variable{<:AbstractDomain},C <: Constraint{<:Function},O <: Objective{<:Function}}
    variables::Dictionary{Int,V}
    constraints::Dictionary{Int,C}
    objectives::Dictionary{Int,O}

    # counter to add new variables: vars, cons, objs
    max_vars::Ref{Int}
    max_cons::Ref{Int}
    max_objs::Ref{Int}

    # Bool to indicate if the _Model instance has been specialized (relatively to types)
    specialized::Ref{Bool}

    # Symbol to indicate the kind of model for specialized methods such as pretty printing
    kind::Symbol
end
```
"""
struct _Model{V <: Variable{<:AbstractDomain},C <: Constraint{<:Function},O <: Objective{<:Function}}# <: MOI.ModelLike
    variables::Dictionary{Int,V}
    constraints::Dictionary{Int,C}
    objectives::Dictionary{Int,O}

    # counter to add new variables: vars, cons, objs
    max_vars::Ref{Int} # TODO: UInt ?
    max_cons::Ref{Int}
    max_objs::Ref{Int}

    # Sense (min = 1, max = -1)
    sense::Ref{Int}

    # Bool to indicate if the _Model instance has been specialized (relatively to types)
    specialized::Ref{Bool}

    # Symbol to indicate the kind of model for specialized methods such as pretty printing
    kind::Symbol

    # Best known bound
    best_bound::Union{Nothing,Float64}

    # Time of construction (seconds) since epoch
    time_stamp::Float64
end

"""
    model()
Construct a _Model, empty by default. It is recommended to add the constraints, variables, and objectives from an empty _Model. The following keyword arguments are available,
- `vars=Dictionary{Int,Variable}()`: collection of variables
- `cons=Dictionary{Int,Constraint}()`: collection of cosntraints
- `objs=Dictionary{Int,Objective}()`: collection of objectives
- `kind=:generic`: the kind of problem modeled (useful for specialized methods such as pretty printing)
"""
function model(;
    vars=Dictionary{Int,Variable}(),
    cons=Dictionary{Int,Constraint}(),
    objs=Dictionary{Int,Objective}(),
    kind=:generic,
    best_bound=nothing,
)

    max_vars = Ref(zero(Int))
    max_cons = Ref(zero(Int))
    max_objs = Ref(zero(Int))
    sense = Ref(one(Int)) # minimization

    specialized = Ref(false)

    _Model(vars, cons, objs, max_vars, max_cons, max_objs, sense, specialized, kind, best_bound, time())
end

"""
    _max_vars(m::M) where M <: Union{Model, AbstractSolver}
Access the maximum variable id that has been attributed to `m`.
"""
_max_vars(m::_Model) = m.max_vars.x

"""
    _max_cons(m::M) where M <: Union{Model, AbstractSolver}
Access the maximum constraint id that has been attributed to `m`.
"""
_max_cons(m::_Model) = m.max_cons.x

"""
    _max_objs(m::M) where M <: Union{Model, AbstractSolver}
Access the maximum objective id that has been attributed to `m`.
"""
_max_objs(m::_Model) = m.max_objs.x

_best_bound(m::_Model) = m.best_bound

"""
    _inc_vars!(m::M) where M <: Union{Model, AbstractSolver}
Increment the maximum variable id that has been attributed to `m`.
"""
_inc_vars!(m::_Model) = m.max_vars.x += 1

"""
    _inc_vars!(m::M) where M <: Union{Model, AbstractSolver}
Increment the maximum constraint id that has been attributed to `m`.
"""
_inc_cons!(m::_Model) = m.max_cons.x += 1

"""
    _inc_vars!(m::M) where M <: Union{Model, AbstractSolver}
Increment the maximum objective id that has been attributed to `m`.
"""
_inc_objs!(m::_Model) = m.max_objs.x += 1

"""
    get_variables(m::M) where M <: Union{Model, AbstractSolver}
Access the variables of `m`.
"""
get_variables(m::_Model) = m.variables

"""
    get_constraints(m::M) where M <: Union{Model, AbstractSolver}
Access the constraints of `m`.
"""
get_constraints(m::_Model) = m.constraints

"""
    get_objectives(m::M) where M <: Union{Model, AbstractSolver}
Access the objectives of `m`.
"""
get_objectives(m::_Model) = m.objectives

"""
    get_kind(m::M) where M <: Union{Model, AbstractSolver}
Access the kind of `m`, such as `:sudoku` or `:generic` (default).
"""
get_kind(m::_Model) = m.kind

"""
    get_variable(m::M, x) where M <: Union{Model, AbstractSolver}
Access the variable `x`.
"""
get_variable(m::_Model, x) = get_variables(m)[x]

"""
    get_constraint(m::M, c) where M <: Union{Model, AbstractSolver}
Access the constraint `c`.
"""
get_constraint(m::_Model, c) = get_constraints(m)[c]

"""
    get_objective(m::M, o) where M <: Union{Model, AbstractSolver}
Access the objective `o`.
"""
get_objective(m::_Model, o) = get_objectives(m)[o]

"""
    get_domain(m::M, x) where M <: Union{Model, AbstractSolver}
Access the domain of variable `x`.
"""
get_domain(m::_Model, x) = get_domain(get_variable(m, x))

"""
    get_name(m::M, x) where M <: Union{Model, AbstractSolver}
Access the name of variable `x`.
"""
get_name(::_Model, x) = "x$x"

"""
    get_cons_from_var(m::M, x) where M <: Union{Model, AbstractSolver}
Access the constraints restricting variable `x`.
"""
get_cons_from_var(m::_Model, x) = _get_constraints(get_variable(m, x))

"""
    get_vars_from_cons(m::M, c) where M <: Union{Model, AbstractSolver}
Access the variables restricted by constraint `c`.
"""
get_vars_from_cons(m::_Model, c) = _get_vars(get_constraint(m, c))

"""
    get_time_stamp(m::M) where M <: Union{Model, AbstractSolver}
Access the time (since epoch) when the model was created. This time stamp is for internal performance measurement.
"""
get_time_stamp(m::_Model) = m.time_stamp

"""
    length_var(m::M, x) where M <: Union{Model, AbstractSolver}
Return the domain length of variable `x`.
"""
length_var(m::_Model, x) = length(get_variable(m, x))

"""
    length_cons(m::M, c) where M <: Union{Model, AbstractSolver}
Return the length of constraint `c`.
"""
length_cons(m::_Model, c) = _length(get_constraint(m, c))

"""
    length_objs(m::M) where M <: Union{Model, AbstractSolver}
Return the number of objectives in `m`.
"""
length_objs(m::_Model) = length(get_objectives(m))

"""
    length_vars(m::M) where M <: Union{Model, AbstractSolver}
Return the number of variables in `m`.
"""
length_vars(m::_Model) = length(get_variables(m))

"""
    length_cons(m::M) where M <: Union{Model, AbstractSolver}
Return the number of constraints in `m`.
"""
length_cons(m::_Model) = length(get_constraints(m))

"""
    draw(m::M, x) where M <: Union{Model, AbstractSolver}
Draw a random value of `x` domain.
"""
draw(m::_Model, x) = rand(get_variable(m, x))
draw(m::_Model, x, n) = rand(get_variable(m, x), n)

"""
    constriction(m::M, x) where M <: Union{Model, AbstractSolver}
Return the constriction of variable `x`.
"""
constriction(m::_Model, x) = _constriction(get_variable(m, x))

"""
    delete_value(m::M, x, val) where M <: Union{Model, AbstractSolver}
Delete `val` from `x` domain.
"""
delete_value!(m::_Model, x, val) = delete!(get_variable(m, x), val)

"""
    delete_var_from_cons(m::M, c, x) where M <: Union{Model, AbstractSolver}
Delete `x` from the constraint `c` list of restricted variables.
"""
delete_var_from_cons!(m::_Model, c, x) = _delete!(get_constraint(m, c), x)

"""
    add_value!(m::M, x, val) where M <: Union{Model, AbstractSolver}
Add `val` to `x` domain.
"""
add_value!(m::_Model, x, val) = add!(get_variable(m, x), val)

"""
    add_var_to_cons!(m::M, c, x) where M <: Union{Model, AbstractSolver}
Add `x` to the constraint `c` list of restricted variables.
"""
add_var_to_cons!(m::_Model, c, x) = _add!(get_constraint(m, c), x)

"""    mts = - get_time_stamp(model)
return TimeStamps(mts, mts, mts, mts, mts, mts, mts)
end

    add!(m::M, x) where M <: Union{Model, AbstractSolver}
    add!(m::M, c) where M <: Union{Model, AbstractSolver}
    add!(m::M, o) where M <: Union{Model, AbstractSolver}
Add a variable `x`, a constraint `c`, or an objective `o` to `m`.
"""
function add!(m::_Model, x::Variable)
    _inc_vars!(m)
    insert!(get_variables(m), _max_vars(m), x)
end

function add!(m::_Model, c::Constraint)
    _inc_cons!(m)
    insert!(get_constraints(m), _max_cons(m), c)
    foreach(x -> _add_to_constraint!(m.variables[x], _max_cons(m)), c.vars)
end

function add!(m::_Model, o::Objective)
    _inc_objs!(m)
    insert!(get_objectives(m), _max_objs(m), o)
end

"""
    variable!(m::M, d) where M <: Union{Model, AbstractSolver}
Add a variable with domain `d` to `m`.
"""
function variable!(m::_Model, d=domain())
    add!(m, variable(d))
    return _max_vars(m)
end

"""
    constraint!(m::M, func, vars) where M <: Union{Model, AbstractSolver}
Add a constraint with an error function `func` defined over variables `vars`.
"""
function constraint!(m::_Model, func, vars::V) where {V <: AbstractVector{<:Number}}
    add!(m, constraint(func, vars))
    return _max_cons(m)
end

"""
    objective!(m::M, func) where M <: Union{Model, AbstractSolver}
Add an objective evaluated by `func`.
"""
function objective!(m::_Model, func)
    add!(m, objective(func, "o" * string(_max_objs(m) + 1)))
end

"""
    describe(m::M) where M <: Union{Model, AbstractSolver}
Describe the model.
"""
function describe(m::_Model) # TODO: rewrite _describe
    objectives = ""
    if Dictionaries.length(m.objectives) == 0
        objectives = "Constraint Satisfaction Program (CSP)"
    else
        objectives = "Constraint Optimization Program (COP) with Objective(s)\n"
        objectives *=
            mapreduce(o -> "\t\t" * o.name * "\n", *, get_objectives(m); init="")[1:end - 1]
    end
    variables = mapreduce(
        x -> "\t\tx$(x[1]): " * string(get_domain(x[2])) * "\n",
        *, pairs(m.variables); init=""
    )[1:end - 1]
    constraints = mapreduce(c -> "\t\tc$(c[1]): " * string(c[2].vars) * "\n", *, pairs(m.constraints); init="")[1:end - 1]

    str =
    """
    _Model description
        $objectives
        Variables: $(length(m.variables))
    $variables
        Constraints: $(length(m.constraints))
    $constraints
    """

    return str
end

"""
    is_sat(m::M) where M <: Union{Model, AbstractSolver}
Return `true` if `m` is a satisfaction model.
"""
is_sat(m::_Model) = length_objs(m) == 0

"""
    is_specialized(m::M) where M <: Union{Model, AbstractSolver}
Return `true` if the model is already specialized.
"""
function is_specialized(m::_Model)
    return m.specialized.x
end

"""
    specialize(m::M) where M <: Union{Model, AbstractSolver}
Specialize the structure of a model to avoid dynamic type attribution at runtime.
"""
function specialize(m::_Model)
    vars_types = Set{Type}()
    cons_types = Set{Type}()
    objs_types = Set{Type}()

    foreach(x -> push!(vars_types, typeof(x.domain)), get_variables(m))
    foreach(c -> push!(cons_types, typeof(c.f)), get_constraints(m))
    foreach(o -> push!(objs_types, typeof(o.f)), get_objectives(m))

    vars_union = _to_union(vars_types)
    cons_union = _to_union(cons_types)
    objs_union = _to_union(objs_types)

    vars = similar(get_variables(m), Variable{vars_union})
    cons = similar(get_constraints(m), Constraint{cons_union})
    objs = similar(get_objectives(m), Objective{objs_union})

    foreach(x -> vars[x] = Variable(vars_union, get_variable(m, x)), keys(vars))
    foreach(c -> cons[c] = Constraint(cons_union, get_constraint(m, c)), keys(cons))
    foreach(o -> objs[o] = Objective(objs_union, get_objective(m, o)), keys(objs))

    max_vars = Ref(_max_vars(m))
    max_cons = Ref(_max_cons(m))
    max_objs = Ref(_max_objs(m))

    specialized = Ref(true)

    _Model(vars, cons, objs, max_vars, max_cons, max_objs, m.sense, specialized,
        get_kind(m), _best_bound(m), get_time_stamp(m))
end

"""
    _is_empty(m::Model)

DOCSTRING
"""
function _is_empty(m::_Model)
    return length_objs(m) + length_vars(m) == 0
end

"""
    _set_domain!(m::Model, x, values)

DOCSTRING

# Arguments:
- `m`: DESCRIPTION
- `x`: DESCRIPTION
- `values`: DESCRIPTION
"""
function _set_domain!(m::_Model, x, values)
    d = domain(values)
    m.variables[x] = Variable(d, get_cons_from_var(m, x))
end

function _set_domain!(m::_Model, x, a::Tuple, b::Tuple)
    d = domain(a, b)
    m.variables[x] = Variable(d, get_cons_from_var(m, x))
end

# function _set_domain!(m::_Model, x,r::R) where {R <: AbstractRange}
#     d = domain(r)
#     m.variables[x] = Variable(d, get_cons_from_var(m, x))
# end

function _set_domain!(m::_Model, x, d::D) where {D <: AbstractDomain}
    m.variables[x] = Variable(d, get_cons_from_var(m, x))
end

"""
    domain_size(m::Model, x) = begin

DOCSTRING
"""
domain_size(m::_Model, x) = domain_size(get_variable(m, x))

"""
    max_domains_size(m::Model, vars) = begin

DOCSTRING
"""
max_domains_size(m::_Model, vars) = maximum(map(x -> domain_size(m, x), vars))

"""
    empty!(m::Model)

DOCSTRING
"""
function Base.empty!(m::_Model)
    empty!(m.variables)
    empty!(m.constraints)
    empty!(m.objectives)
    m.max_vars[] = 0
    m.max_cons[] = 0
    m.max_objs[] = 0
    m.specialized[] = false
end

# Non modificating cost and objective computations

draw(m::_Model) = map(rand, get_variables(m))

compute_cost(c::Constraint, values, X) = apply(c, map(x -> values[x], c.vars), X)
function compute_costs(m, values, X)
    return sum(c -> compute_cost(c, values, X), get_constraints(m); init = 0.0)
end
function compute_costs(m, values, cons, X)
    return  sum(c -> compute_cost(c, values, X), view(get_constraints(m), cons); init = 0.0)
end

compute_objective(m, values; objective = 1) = apply(get_objective(m, objective), values)

function update_domain!(m, x, d)
    if isempty(get_variable(m,x))
        old_d = get_variable(m, x).domain
        _set_domain!(m, x, d.domain)
    else
        old_d = get_variable(m, x).domain
        are_continuous = typeof(d) <: ContinuousDomain && typeof(old_d) <: ContinuousDomain
        new_d = if are_continuous
            intersect_domains(old_d, d)
        else
            intersect_domains(convert(RangeDomain,old_d), convert(RangeDomain, d))
        end
        _set_domain!(m, x, new_d)
    end
end

sense(m::_Model) = m.sense[]
sense!(m::_Model, ::Val{:min}) = m.sense[] = 1
sense!(m::_Model, ::Val{:max}) = m.sense[] = -1


text:
function read_file(datafile, datarow=2, rename=true, delim=';')
    df = CSV.read(datafile.filename, delim=delim, datarow=datarow, copycols=true)
    rename && rename!(df, datafile)

    return df
end

function header(df, delim)
    buffer = join(names(df), delim)

    return buffer * (@static Sys.iswindows() ? "\r\n" : '\n')
end

function write_file(datafile::DataFile, df, delim)
    ncols = length(eachcol(df))
    info = comment_value(datafile)
    units = join(datafile.units, delim)
    new_line = join(repeat([info], ncols), delim)

    write_file(df, new_line, datafile.savename, delim)
end

function write_file(df::AbstractDataFrame, extra, filename, delim)
    buffer = IOBuffer()
    nl = @static Sys.iswindows() ? "\r\n" : '\n'
    df |> CSV.write(buffer, delim=delim, writeheader=false, newline=nl)

    h = header(df, delim)
    file = String(take!(buffer))
    if extra isa Tuple
        content = h * join(extra, nl) * nl
    else
        content = h * extra * nl
    end

    open(filename, "w") do f
        write(f, content, file)
    end
    return nothing
end


text:
"""
Abstract type to represent Electroencephalography (EEG) data.

The following types inherit from the EEG type and can be used to process your data:

- `GeneralEEG`: Used to store data without assumption of any experimental paradigm.
- `SSR`: Used to store data acquired with a steady state response experiment paradigm.
- `TR`: Used to store data acquired with a transient response experiment paradigm.

# Examples
```julia
data = # load your EEG data using for example read_EEG()

samplingrate(data)  # Returns the sampling rate
channelnames(data)  # Returns the channel names
```
"""
abstract type EEG <: NeuroimagingMeasurement end

import Base.show
function Base.show(io::IO, a::EEG)
    time_length = round.(size(a.data, 1) / samplingrate(Float64, a) / 60)
    println(
        io,
        "EEG measurement of $time_length mins with $(size(a.data,2)) channels sampled at $(a.samplingrate)",
    )
end

"""
Type for storing general EEG data without assumption of any experimental paradigm.

# Examples
```julia
s = read_EEG(filename)
s = rereference(s, "Cz")
s = remove_channel!(s, "Cz")
```
"""
mutable struct GeneralEEG <: EEG
    data::Array
    sensors::Array{Sensor}
    triggers::Dict
    system_codes::Dict
    samplingrate::typeof(1.0u"Hz")
    reference_channel::Array{AbstractString,1}
    file_path::AbstractString
    file_name::AbstractString
    processing::Dict
    header::Dict
end
GeneralEEG(args...; kwargs...) = read_EEG(args...; kwargs...)

"""
    samplingrate(t::Type, s::EEG)
    samplingrate(s::EEG)

Return the sampling rate of an EEG type in Hz as the requested type.
If no type is provided, the sampling rate is returned as a floating point number.

# Examples
```julia
s = read_EEG(filename)
samplingrate(s)
```
"""
samplingrate(s::EEG) = s.samplingrate |> u"Hz"
samplingrate(t, s::EEG) = convert(t, s.samplingrate |> u"Hz" |> ustrip)

"""
    times(t::Type, s::EEG)
    times(s::EEG)

Return the times associated with each sample in seconds.

# Examples
```julia
s = read_EEG(filename)
times(s)
```
"""
times(s::EEG) = collect(1:size(data(s), 1)) ./ samplingrate(Float64, s) .* u"s"
times(t, s::EEG) = convert.(t, times(s) |> ustrip)

"""
    channelnames(s::EEG)

Return the names of sensors in EEG measurement.

# Examples
```julia
s = read_EEG(filename)
channelnames(s)
```
"""
channelnames(s::EEG) = labels(s.sensors)

"""
    channelnames(s::EEG, i::Int, l::AbstractString)
    channelnames(s::EEG, l::AbstractVector{AbstractString})

Change the names of `i`th sensors in an EEG measurement `s` to `l`.
Or change the name of all sensors by pass a vector of strings.

# Examples
```julia
s = read_EEG(filename)
channelnames(s, 1, "Fp1")
```
"""
function channelnames(s::EEG, i::Int, l::S) where {S<:AbstractString}

    s.sensors[i].label = l
    return s
end
function channelnames(s::EEG, l::AbstractVector{S}) where {S<:AbstractString}
    @assert length(l) == length(channelnames(s))
    for li = 1:length(l)
        s = channelnames(s, li, l[li])
    end
    return s
end

"""
    data(s::EEG)
    data(s::EEG, channel::AbstractString)
    data(s::EEG, channels::AbstractVector{AbstractString})

Return the data stored in the type object.
This may be useful for integration with custom processing.

# Examples
```julia
s = read_EEG(filename)
data(s)
data(s, "Cz")
data(s, ["Cz", "P2"])
```
"""
data(s::EEG) = s.data
data(s::EEG, channel::AbstractString) = data(s, [channel])
function data(s::EEG, channels::AbstractVector{S}) where {S<:AbstractString}
    return data(s, _channel_indices(s, channels))
end
function data(s::EEG, channels::AbstractVector{S}) where {S<:Int}
    return deepcopy(s.data)[:, channels]
end

"""
    sensors(s::EEG)
    channelnames(s::EEG, l::AbstractVector{AbstractString})

Returns the sensors for an EEG recording.

# Examples
```julia
s = read_EEG(filename)
sensors(s)
```
"""
sensors(s::EEG) = s.sensors
electrodes(s::EEG) = s.sensors

#######################################
#
# EEG type operations
#
#######################################

import Base.hcat
"""
    hcat(a::EEG, b::EEG)

Concatenate two EEG measurements together, effectively creating a single long measurement.

# Examples
```julia
hcat(a, b)
```
"""
function hcat(a::EEG, b::EEG)

    if channelnames(a) != channelnames(b)
        throw(
            ArgumentError(
                string("Channels do not match $(channelnames(a)) != $(channelnames(b))"),
            ),
        )
    end

    if haskey(a.processing, "epochs")
        @warn("Epochs have already been extracted and will no longer be valid")
    end
    if haskey(a.processing, "statistics")
        @warn("Statistics have already been calculated and will no longer be valid")
    end

    @debug(
        "Appending two EEGs with $(size(a.data, 2)) .& $(size(b.data, 2)) channels and lengths $(size(a.data, 1)) $(size(b.data, 1))"
    )

    join_triggers(a, b)
    a.data = [a.data; b.data]

    return a
end

"""
Append the trigger information of one EEG type to another.
Places the trigger information at the end of first file

#### Example

```julia
join_triggers(a, b)
```
"""
function join_triggers(a, b; offset = size(a.data, 1))

    a.triggers["Index"] = [a.triggers["Index"]; (b.triggers["Index"] .+ offset)]
    a.triggers["Code"] = [a.triggers["Code"]; b.triggers["Code"]]
    a.triggers["Duration"] = [a.triggers["Duration"]'; b.triggers["Duration"]']'

    a
end

#######################################
#
# Change reference channels
#
#######################################

"""
    rereference(a::EEG, refChan::Union{AbstractString, Array{AbstractString}}; kwargs...)

Reference data to specified channel(s).

#### Example

```julia
a = rereference(a, "Cz")
# or
a = rereference(a, ["P9", "P10"])
```

"""
function rereference(
    a::EEG,
    refChan::Union{S,Array{S}};
    kwargs...,
) where {S<:AbstractString}

    a.data = rereference(a.data, refChan, channelnames(a))

    a.reference_channel = [refChan]

    return a
end

#######################################
#
# Manipulate channels
#
#######################################

"""
    add_channel(a::EEG, data::Vector, chanLabel::AbstractString)

Add a channel to the EEG type with specified channel names.

# Examples
```julia
s = read_EEG(filename)
new_channel = mean(s.data, 2)
s = add_channel(s, new_channel, "MeanChannelData")
```
"""
function add_channel(a::EEG, data::Vector, chanLabel::AbstractString; kwargs...)

    @info("Adding channel $chanLabel")

    a.data = hcat(a.data, data)
    push!(a.sensors, Electrode(chanLabel, Talairach(NaN, NaN, NaN), Dict()))

    return a
end

"""
    remove_channel!(a::EEG, channelname::AbstractString)
    remove_channel!(a::EEG, channelnames::Array{AbstractString})
    remove_channel!(a::EEG, channelidx::Int)
    remove_channel!(a::EEG, channelidxs::Array{Int})

Remove channel(s) from EEG as specifed by `channelname` or `channelidx`.

# Examples
```julia
a = read_EEG(filename)
remove_channel!(a, ["TP8", "Cz"])
```
"""
function remove_channel!(a::EEG, channel_name::S; kwargs...) where {S<:AbstractString}
    @debug("Removing channel $(channel_name)")
    remove_channel!(a, [channel_name])
end

function remove_channel!(
    a::EEG,
    channel_names::Array{S};
    kwargs...,
) where {S<:AbstractString}
    @debug("Removing channels $(join(channel_names, " "))")
    remove_channel!(a, _channel_indices(a, channel_names))
end

remove_channel!(a::EEG, channel_names::Int; kwargs...) =
    remove_channel!(a, [channel_names]; kwargs...)

function remove_channel!(a::EEG, channel_idx::Array{Int}; kwargs...)

    channel_idx = channel_idx[channel_idx.!=0]

    @debug("Removing channel(s) $channel_idx")
    if any(channel_idx .== 0)
        @warn("Failed to remove a channel")
    end

    keep_idx = [1:size(a.data)[end];]
    for c in sort(channel_idx, rev = true)
        try
            splice!(keep_idx, c)
        catch
            # Nothing
        end
    end

    if haskey(a.processing, "epochs")
        if size(a.processing["epochs"], 3) == size(a.data, 2)
            @debug("Removing channel(s) from epoch data")
            a.processing["epochs"] = a.processing["epochs"][:, :, keep_idx]
        end
    end
    if haskey(a.processing, "sweeps")
        if size(a.processing["sweeps"], 3) == size(a.data, 2)
            @debug("Removing channel(s) from sweep data")
            a.processing["sweeps"] = a.processing["sweeps"][:, :, keep_idx]
        end
    end

    a.data = a.data[:, keep_idx]
    a.sensors = a.sensors[keep_idx]

    return a
end

"""
    keep_channel!(a::EEG, channelname::AbstractString)
    keep_channel!(a::EEG, channelnames::Array{AbstractString})
    keep_channel!(a::EEG, channelidxs::Array{Int})

Remove all channels except those requested from EEG.

# Examples
```julia
a = read_EEG(filename)
keep_channel!(a, ["P8", "Cz"])
```
"""
function keep_channel!(a::EEG, channel_name::AbstractString; kwargs...)
    keep_channel!(a, [channel_name]; kwargs...)
end

function keep_channel!(
    a::EEG,
    channel_names::AbstractVector{S};
    kwargs...,
) where {S<:AbstractString}
    @info("Keeping channel(s) $(join(channel_names, " "))")
    keep_channel!(a, _channel_indices(a, channel_names))
end

function keep_channel!(a::EEG, channel_idx::AbstractVector{Int}; kwargs...)

    remove_channels = [1:size(a.data, 2);]

    channel_idx = sort(channel_idx, rev = true)
    for c in channel_idx
        splice!(remove_channels, c)
    end

    remove_channel!(a, remove_channels; kwargs...)
end

#######################################
#
# Trim channels
#
#######################################

"""
    trim_channel(a::EEG, stop::Int; start::Int=1)

Trim EEG recording by removing data after `stop` specifed samples
and optionally before `start` samples.

# Examples
```julia
s = trim_channel(s, 8192*300, start=8192)
```
"""
function trim_channel(a::EEG, stop::Int; start::Int = 1, kwargs...)

    @info("Trimming $(size(a.data)[end]) channels between $start and $stop")

    a.data = a.data[start:stop, :]

    to_keep = findall((a.triggers["Index"] .>= start) .& (a.triggers["Index"] .<= stop))
    a.triggers["Index"] = a.triggers["Index"][to_keep]
    a.triggers["Duration"] = a.triggers["Duration"][to_keep]
    a.triggers["Code"] = a.triggers["Code"][to_keep]
    a.triggers["Index"] .-= (start - 1)

    to_keep =
        findall((a.system_codes["Index"] .>= start) .& (a.system_codes["Index"] .<= stop))
    a.system_codes["Index"] = a.system_codes["Index"][to_keep]
    a.system_codes["Duration"] = a.system_codes["Duration"][to_keep]
    a.system_codes["Code"] = a.system_codes["Code"][to_keep]
    a.system_codes["Index"] .-= (start - 1)

    return a
end

#######################################
#
# Merge channels
#
#######################################

"""
    merge_channels(a::EEG, merge_Chans::Array{S}, new_name::S) where {S<:AbstractString}
    merge_channels(a::EEG, merge_Chans::S, new_name::S) where {S<:AbstractString}

Average `EEG` channels listed in `merge_Chans` and label the averaged channel as `new_name`.

# Examples
```julia
s = merge_channels(s, ["P6", "P8"], "P68")
```
"""
function merge_channels(
    a::EEG,
    merge_Chans::Array{S},
    new_name::S;
    kwargs...,
) where {S<:AbstractString}

    @debug("Number of original channels: $(length(channelnames(a)))")

    keep_idxs = _channel_indices(a, merge_Chans)

    if sum(keep_idxs .== 0) > 0
        @warn(
            "Could not merge as these channels don't exist: $(join(vec(merge_Chans[keep_idxs .== 0]), " "))"
        )
        keep_idxs = keep_idxs[keep_idxs.>0]
    end

    @info("Merging channels $(join(vec(channelnames(a)[keep_idxs,:]), " "))")
    @debug("Merging channels $keep_idxs")

    a = add_channel(
        a,
        vec(Statistics.mean(a.data[:, keep_idxs], dims = 2)),
        new_name;
        kwargs...,
    )
end

function merge_channels(
    a::EEG,
    merge_Chans::S,
    new_name::S;
    kwargs...,
) where {S<:AbstractString}
    a = merge_channels(a, [merge_Chans], new_name; kwargs...)
end

"""
    read_EEG(fname::AbstractString)
    read_EEG(args...)

Read a file or IO stream and store the data in an `GeneralEEG` type.

# Arguments

- `fname`: Name of the file to be read
- `min_epoch_length`: Minimum epoch length in samples. Shorter epochs will be removed (0)
- `max_epoch_length`: Maximum epoch length in samples. Longer epochs will be removed (0 = all)
- `valid_triggers`: Triggers that are considered valid, others are removed ([1,2])
- `stimulation_amplitude`: Amplitude of stimulation (NaN)
- `remove_first`: Number of epochs to be removed from start of recording (0)
- `max_epochs`: Maximum number of epochs to retain (0 = all)

# Supported file formats

- BIOSEMI (.bdf)
"""
function read_EEG(
    fname::AbstractString;
    valid_triggers::Array{Int} = [1, 2],
    min_epoch_length::Int = 0,
    max_epoch_length::Int = 0,
    remove_first::Int = 0,
    max_epochs::Int = 0,
    kwargs...,
)

    @info("Importing EEG from file: $fname")
    file_path, file_name, ext = fileparts(fname)

    #
    # Read file data
    #

    # Import raw data
    if ext == "bdf"
        data, triggers, system_codes, samplingrate, reference_channel, header =
            import_biosemi(fname; kwargs...)
    else
        warn("File type $ext is unknown")
    end

    # Create electrodes
    elecs = Electrode[]
    for e in header["chanLabels"]
        push!(elecs, Electrode(e, Talairach(NaN * u"m", NaN * u"m", NaN * u"m"), Dict()))
    end

    # Create EEG type
    if unit(samplingrate) == unit(1.0u"Hz")
        #nothing
    else
        samplingrate = samplingrate * 1.0u"Hz"
    end

    a = GeneralEEG(
        data,
        elecs,
        triggers,
        system_codes,
        samplingrate,
        [reference_channel],
        file_path,
        file_name,
        Dict(),
        header,
    )

    #
    # Clean up
    #

    # Remove status channel information
    remove_channel!(a, "Status")

    # Clean epoch index
    a.triggers = clean_triggers(
        a.triggers,
        valid_triggers,
        min_epoch_length,
        max_epoch_length,
        remove_first,
        max_epochs,
    )

    # Try and match sensor names to known locations
    locs = read_elp(joinpath(datadep"BioSemi64Locations", "biosemi64.elp"))
    new_sens, idx = match_sensors(locs, labels(sensors(a)))
    if length(new_sens) == length(sensors(a))
        @debug("Sucsessfully matches location of all sensors, using new locations")
        a.sensors = new_sens
    end

    return a
end

function trigger_channel(a::EEG; kwargs...)

    create_channel(a.triggers, a.data, samplingrate(Float64, a))
end

function system_code_channel(a::EEG; kwargs...)

    create_channel(a.system_codes, a.data, samplingrate(Float64, a))
end

"""
    epoch_rejection(a::EEG; retain_percentage::Number = 0.95, kwargs...)

Reject epochs such that `retain_percentage` is retained.
"""
function epoch_rejection(a::EEG; retain_percentage::Number = 0.95, kwargs...)

    a.processing["epochs"] =
        epoch_rejection(a.processing["epochs"], retain_percentage; kwargs...)

    return a
end


text:
# Very basic setup, useful for testing
using DFTK

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]

model = model_LDA(lattice, atoms)
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])

scfres = self_consistent_field(basis, tol=1e-8)


text:
module BearSSL

const bearssl = "/Users/jacobquinn/BearSSL/build/libbearssl.so"

include("trustanchors.jl")

br_x509_minimal_context_size = 3168
br_ssl_engine_context_size = 3576
br_ssl_client_context_size = 3680

struct br_x509_minimal_context
    data::NTuple{br_x509_minimal_context_size, UInt8}
    br_x509_minimal_context() = new(ntuple(x->0x00, br_x509_minimal_context_size))
end
Base.show(io::IO, b::br_x509_minimal_context) = print(io, "BearSSL.br_x509_minimal_context()")

struct br_ssl_engine_context
    data::NTuple{br_ssl_engine_context_size, UInt8}
    br_ssl_engine_context() = new(ntuple(x->0x00, br_ssl_engine_context_size))
end
Base.show(io::IO, b::br_ssl_engine_context) = print(io, "BearSSL.br_ssl_engine_context()")

struct br_ssl_client_context
    eng::br_ssl_engine_context
    data::NTuple{br_ssl_client_context_size - br_ssl_engine_context_size, UInt8}
    function br_ssl_client_context()
        eng = br_ssl_engine_context()
        data = ntuple(x->0x00, br_ssl_client_context_size - br_ssl_engine_context_size)
        client = new(eng, data)
    end
end
Base.show(io::IO, b::br_ssl_client_context) = print(io, "BearSSL.br_ssl_client_context()")

struct br_sslio_context
    eng::br_ssl_engine_context
    low_read::Ptr{Void}
    low_write::Ptr{Void}
    br_sslio_context(eng) = new(eng, C_NULL, C_NULL)
end

struct SSLConfig
    x509::br_x509_minimal_context
end

struct SSLContext <: IO
    client::br_ssl_client_context
    bio::TCPSocket
end

function sock_read_func(ctx::Ptr{Void}, buf, len)
    while true
        rlen = ccall(:read, Cssize_t, (Cint, Ptr{UInt8}, Csize_t), unsafe_load(convert(Ptr{Cint}, ctx)), buf, len)
        if rlen <= 0
            rlen < 0 && errno == EINTR && continue
            return Cint(-1)
        end
        return Cint(rlen)
    end
end
const sock_read = cfunction(sock_read_func, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t))

function sock_write_func(ctx::Ptr{Void}, buf, len)
    while true
        wlen = ccall(:write, Cssize_t, (Cint, Ptr{UInt8}, Csize_t), unsafe_load(convert(Ptr{Cint}, ctx)), buf, len)
    end
end
const sock_write = cfunction(sock_write_func, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t))

function SSLContext()
    tcp = connect("www.google.com", 80)
    fd = Base._fd(tcp).fd
    client = BearSSL.br_ssl_client_context()
    x509 =   BearSSL.br_x509_minimal_context()
    # iobuf_in = zeros(UInt8, 16709)
    # iobuf_out = zeros(UInt8, 16469)
    iobuf = zeros(UInt8, 16709+16469);
    ioc = BearSSL.br_sslio_context(client.eng)

    ccall((:br_ssl_client_init_full, BearSSL.bearssl), Void,
        (Ptr{BearSSL.br_ssl_client_context}, Ptr{BearSSL.br_x509_minimal_context}, Ptr{BearSSL.br_x509_trust_anchor}, Csize_t),
        pointer_from_objref(client), pointer_from_objref(x509), pointer_from_objref(BearSSL.TAs), BearSSL.TAs_NUM)

    ccall((:br_ssl_engine_set_buffer, BearSSL.bearssl), Void,
        (Ptr{BearSSL.br_ssl_engine_context}, Ptr{Void}, Csize_t, Cint),
        pointer_from_objref(client.eng), iobuf, sizeof(iobuf), 1)

    ccall((:br_ssl_client_reset, BearSSL.bearssl), Cint,
        (Ptr{BearSSL.br_ssl_client_context}, Ptr{UInt8}, Cint),
        pointer_from_objref(client), C_NULL, 0)

    ccall((:br_sslio_init, BearSSL.bearssl), Void,
        (Ptr{BearSSL.br_sslio_context}, Ptr{BearSSL.br_ssl_engine_context}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
        pointer_from_objref(ioc), pointer_from_objref(client.eng), BearSSL.sock_read, Ref(fd), BearSSL.sock_write, Ref(fd))

    ccall((:br_sslio_write_all, BearSSL.bearssl), Cint,
        (Ptr{BearSSL.br_sslio_context}, Ptr{UInt8}, Csize_t),
        pointer_from_objref(ioc), "GET", 4)
    ccall((:br_sslio_write_all, BearSSL.bearssl), Cint,
        (Ptr{BearSSL.br_sslio_context}, Ptr{UInt8}, Csize_t),
        pointer_from_objref(ioc), "/", 1)
    ccall((:br_sslio_write_all, BearSSL.bearssl), Cint,
        (Ptr{BearSSL.br_sslio_context}, Ptr{UInt8}, Csize_t),
        pointer_from_objref(ioc), " HTTP/1.0\4\nHost: ", 17)
    ccall((:br_sslio_write_all, BearSSL.bearssl), Cint,
        (Ptr{BearSSL.br_sslio_context}, Ptr{UInt8}, Csize_t),
        pointer_from_objref(ioc), "domorig.io", 10)
    ccall((:br_sslio_write_all, BearSSL.bearssl), Cint,
        (Ptr{BearSSL.br_sslio_context}, Ptr{UInt8}, Csize_t),
        pointer_from_objref(ioc), "\r\n\r\n", 4)

    ccall((:br_sslio_flush, BearSSL.bearssl), Void,
        (Ptr{BearSSL.br_sslio_context},),
        pointer_from_objref(ioc))

    while true
        tmp = zeros(UInt8, 512);
        rlen = ccall((:br_sslio_read, BearSSL.bearssl), Cint,
            (Ptr{BearSSL.br_sslio_context}, Ptr{Void}, Csize_t),
            pointer_from_objref(ioc), tmp, sizeof(tmp))
        rlen < 0 && break
        println(unsafe_string(pointer(tmp), rlen))
    end
end

# ccall((:br_ssl_engine_last_error, BearSSL.bearssl), Cint,
#     (Ptr{BearSSL.br_ssl_engine_context},),
#     pointer_from_objref(client.eng))

# ccall((:br_ssl_engine_set_buffer, BearSSL.bearssl), Void,
#     (Ptr{BearSSL.br_ssl_engine_context},Ptr{Void},Cint,Csize_t),
#     pointer_from_objref(client.eng), C_NULL, 0, 0)

# ccall((:br_ssl_engine_init_rand, BearSSL.bearssl), Cint,
#     (Ptr{BearSSL.br_ssl_engine_context},),
#     pointer_from_objref(client.eng))

# unsafe_load(convert(Ptr{Cint}, pointer_from_objref(client.eng)))

# ccall((:br_ssl_engine_current_state, BearSSL.bearssl), Cuint,
#     (Ptr{BearSSL.br_ssl_engine_context},),
#     pointer_from_objref(client.eng))
# TLS interface methods
# SSLConfig(verify::Bool)
# SSLContext()
# TLS.setup!(stream, opts.tlsconfig::TLS.SSLConfig)
# TLS.hostname!(stream, hostname)
# TLS.associate!(stream, socket)
# TLS.handshake!(stream)

# IO interface methods
Base.isopen(io::SSLContext) = isopen(io.bio)

Base.close(io::SSLContext) = close(io.bio)

function Base.write(io::SSLContext, bytes::Vector{UInt8})

end

function Base.readavailable(io::SSLContext)

end

end # module


text:
export SymmorphicIrrepComponent
export get_irrep_components
export get_irrep_iterator

"""
    SymmorphicIrrepComponent{S1<:SymmetryOrEmbedding, S2<:SymmetryOrEmbedding}

Irrep component of a symmorphic symmetry (embedding).

# Fields
* `normal::IrrepComponent{S1}`
* `rest::IrrepComponent{S2}`
"""
struct SymmorphicIrrepComponent{S1<:SymmetryOrEmbedding, S2<:SymmetryOrEmbedding}<:AbstractSymmetryIrrepComponent
    normal::IrrepComponent{S1}
    rest::IrrepComponent{S2}

    @doc """
        SymmorphicIrrepComponent(normal, rest)

    Construct a `SymmorphicIrrepComponent` of `normal` and `rest`.

    # Arguments
    * `normal::IrrepComponent{S1}`
    * `rest::IrrepComponent{S2}`
    )
    """
    function SymmorphicIrrepComponent(
        normal::IrrepComponent{S1},
        rest::IrrepComponent{S2},
    ) where {S1<:AbstractSymmetry, S2<:AbstractSymmetry}
        if !iscompatible(normal.symmetry, normal.irrep_index, rest.symmetry)
            throw(ArgumentError(
                "symmetry $(symmetry_name(rest.symmetry)) is not compatible with "*
                "symmetry $(symmetry_name(normal.symmetry)) at irrep $(normal.irrep_index)"
            ))
        end
        return new{S1, S2}(normal, rest)
    end

    """
        SymmorphicIrrepComponent(normal, rest)

    # Arguments
    * `normal::IrrepComponent{SymmetryEmbedding{S1}}`
    * `rest::IrrepComponent{SymmetryEmbedding{S2}}`
    """
    function SymmorphicIrrepComponent(
        normal::IrrepComponent{SymmetryEmbedding{S1}},
        rest::IrrepComponent{SymmetryEmbedding{S2}},
    ) where {S1<:AbstractSymmetry, S2<:AbstractSymmetry}
        if !iscompatible(normal.symmetry, normal.irrep_index, rest.symmetry)
            throw(ArgumentError(
                "symmetry $(symmetry_name(rest.symmetry)) is not compatible with "*
                "symmetry $(symmetry_name(normal.symmetry)) at irrep $(normal.irrep_index)"
            ))
        end
        return new{SymmetryEmbedding{S1}, SymmetryEmbedding{S2}}(normal, rest)
    end
end

"""
    group_order(arg::SymmorphicIrrepComponent)

Group order of the underlying symmorphic group of `arg`.
"""
group_order(arg::SymmorphicIrrepComponent) = group_order(arg.normal) * group_order(arg.rest)

"""
    get_irrep_components(sym::SymmorphicSymmetry)

Return an iterator which iterates over the irrep component of the symmorphic symmetry `sym`.
"""
function get_irrep_components(sym::SymmorphicSymmetry)
    return (
        SymmorphicIrrepComponent(normal_sic, rest_sic)
        for normal_sic in get_irrep_components(sym.normal)
        for rest_sic in get_irrep_components(little_symmetry(normal_sic, sym.rest))
    )
end

"""
    get_irrep_components(sym::SymmorphicSymmetryEmbedding)

Return an iterator which iterates over the irrep component of the symmorphic symmetry
embedding `sym`.
"""
function get_irrep_components(sym::SymmorphicSymmetryEmbedding)
    return (
        SymmorphicIrrepComponent(normal_sic, rest_sic)
        for normal_sic in get_irrep_components(sym.normal)
        for rest_sic in get_irrep_components(little_symmetry(normal_sic, sym.rest))
    )
end

"""
    get_irrep_iterator(ssic::SymmorphicIrrepComponent)

Return an iterator which iterates over the elements, together with their corresponding irrep
coefficient of the given irrep component.

# Return
* [(r⋅n, cᵣ⋅cₙ) for all r for all n]
"""
function get_irrep_iterator(ssic::SymmorphicIrrepComponent)
    normal_elements, normal_irrep_components = let
        nsym = ssic.normal.symmetry
        nii  = ssic.normal.irrep_index
        nic  = ssic.normal.irrep_component
        elements(nsym), [m[nic, nic] for m in irrep(nsym, nii)]
    end
    rest_elements, rest_irrep_components = let
        rsym = ssic.rest.symmetry
        rii  = ssic.rest.irrep_index
        ric  = ssic.rest.irrep_component
        elements(rsym), [m[ric, ric] for m in irrep(rsym, rii)]
    end
    return (
        (relem * nelem, rphase * nphase)
        for (relem, rphase) in zip(rest_elements, rest_irrep_components)
        for (nelem, nphase) in zip(normal_elements, normal_irrep_components)
    )
end


text:
module NEOs

# __precompile__(false)

export propagate, observer_position, delay_doppler, ismonostatic,
    mas2rad, t2c_rotation_iau_76_80,
    process_radar_data_jpl, RadarDataJPL, julian2etsecs, etsecs2julian,
    RNp1BP_pN_A_J23E_J2S_ng_eph!, RNp1BP_pN_A_J23E_J2S_ng_eph_threads!,
    utcepoch, delay, delay_sigma, delay_units, doppler, doppler_sigma,
    doppler_units, freq, rcvr, xmit, bouncepoint, valsecchi_circle,
    radec, radec_table, nrms, chi2, newtonls, newtonls_6v, diffcorr,
    newtonls_Q, readfwf, readmp, w8sveres17, bopik, yarkp2adot, pv2kep,
    x0_JPL_s197, x0_JPL_s199

using Distributed
using TaylorIntegration
using Printf, DelimitedFiles, Test, LinearAlgebra
using Dates: DateTime, julian2datetime, datetime2julian
import PlanetaryEphemeris
using PlanetaryEphemeris: daysec, su, ea, α_p_sun, δ_p_sun,
    t2c_jpl_de430, pole_rotation, au, c_au_per_day, R_sun,
    c_cm_per_sec, c_au_per_sec, yr, RE, TaylorInterpolant, Rx, Ry, Rz,
    semimajoraxis, eccentricity, inclination, longascnode, argperi,
    timeperipass, nbodyind
using JLD
using EarthOrientation, SPICE
using Dates
using Quadmath
using Healpix: ang2pixRing, Resolution
using LazyArtifacts
using DataFrames
import SatelliteToolbox
using SatelliteToolbox: nutation_fk5, J2000toGMST, rECEFtoECI,
    get_ΔAT, JD_J2000, EOPData_IAU1980, rECItoECI, DCM,
    TOD, GCRF, ITRF, rECItoECI, PEF, satsv, EOPData_IAU2000A
import RemoteFiles
using StaticArrays: SArray, @SVector
using InteractiveUtils

# integration parameters
const order = 30
const abstol = 1.0E-30

# vector of GM's (DE430 values)
const μ_DE430 = PlanetaryEphemeris.μ
const μ_B16_DE430 = μ_DE430[12:27] # DE430 GM's of 16 most massive asteroids
const μ_ast343_DE430 = μ_DE430[12:end] # DE430 GM's of 343 main belt asteroids included in DE430 integration

# standard value of nominal mean angular velocity of Earth (rad/sec), ESAA 2014 Sec 7.4.3.3 p. 296
const ω = 7.2921151467e-5 # 7.292115e-5 rad/sec
# The relationship of the angular velocity of the earth omega with LOD is (https://www.iers.org/IERS/EN/Science/EarthRotation/UT1LOD.html)
# where `omega` is in units of rad/sec, and `lod` is in units of milliseconds
omega(lod) = (1e-12)*(72921151.467064 - 0.843994809lod)

const A_sun = 1.06e8 # Solar corona parameter A [cm^-3] (ESAA 2014, Table 8.5 p. 329)
const a_sun = 4.89e5 # Solar corona parameter a [cm^-3] (ESAA 2014, Table 8.5 p. 329)
const b_sun = 3.91e5 # Solar corona parameter b [cm^-3] (ESAA 2014, Table 8.5 p. 329)

const S0_sun = 63.15E6 # Sun radiated power intensity at photosphere surface, Watt/meter^2
# const m2_s3_to_au2_day3 = 1e-6daysec^3/au^2 # conversion factor from m^2/sec^3 to au^2/day^3

# vector of J2*R^2 values
const Λ2 = zeros(11)
Λ2[ea] = 1.9679542578489185e-12
# vector of J3*R^3 values
const Λ3 = zeros(11)
Λ3[ea] = -1.962633335678878e-19
const clightkms = 2.99792458E5 # speed of light, km/sec

# isless overloads workaround to make rECEFtoECI work with TaylorSeries
import Base.isless
isless(a::Union{Taylor1,TaylorN}, b::Union{Taylor1,TaylorN}) = isless(constant_term(a), constant_term(b))
isless(a::Real, b::Union{Taylor1,TaylorN}) = isless(a, constant_term(b))
isless(a::Union{Taylor1,TaylorN}, b::Real) = isless(constant_term(a), b)

import ReferenceFrameRotations: orthonormalize
# This method extends orthonormalize to handle Taylor1 and TaylorN
# The ReferenceFrameRotations.jl package is licensed under the MIT "Expat" License
# Copyright (c) 2014-2019: Ronan Arraes Jardim Chagas.
function orthonormalize(dcm::SArray{Tuple{3,3},T,2,9} where {T<:Union{Taylor1,TaylorN}})
    e₁ = dcm[:,1]
    e₂ = dcm[:,2]
    e₃ = dcm[:,3]

    en₁  = e₁/norm(e₁)
    enj₂ =   e₂ - (en₁⋅e₂)*en₁
    en₂  = enj₂ / norm(enj₂)
    enj₃ =   e₃ - (en₁⋅e₃)*en₁
    enj₃ = enj₃ - (en₂⋅enj₃)*en₂
    en₃  = enj₃ / norm(enj₃)

    SArray(  hcat(hcat(en₁, en₂), en₃)  )
end

function yarkp2adot(A2, a, e, μ_S)
    return 2A2/(sqrt(a)*(1-e^2)*sqrt(μ_S))
end

function pv2kep(xas, μ_S, jd=JD_J2000)
    ec0 = eccentricity(xas..., μ_S, 0.0)
    a0 = semimajoraxis(xas..., μ_S, 0.0)
    qr0 = a0*(1-ec0)
    tp0 = timeperipass(jd, xas..., μ_S, 0.0)
    Ω0 = rad2deg(longascnode(xas...))
    ω0 = rad2deg(argperi(xas..., μ_S, 0.0))
    i0 = rad2deg(inclination(xas...))
    return (e=ec0, q=qr0, tp=tp0, Ω=Ω0, ω=ω0, i=i0, a=a0)
end

include("process_radar_data_jpl.jl")
include("topocentric.jl")
include("delay_doppler.jl")
include("radec.jl")
include("asteroid_dynamical_models.jl")
include("initial_conditions.jl")
include("integration_methods.jl")
include("propagation.jl")
include("b_plane.jl")
include("least_squares.jl")

end


text:
using StreamDatasets
using Test

@testset "StreamDatasets.jl" begin
    # Write your tests here.
    @test true = true
end


text:
# ---
# title: 1252. Cells with Odd Values in a Matrix
# id: problem1252
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/cells-with-odd-values-in-a-matrix/description/>
# hidden: true
# ---
#
# Given `n` and `m` which are the dimensions of a matrix initialized by zeros
# and given an array `indices` where `indices[i] = [ri, ci]`. For each pair of
# `[ri, ci]` you have to increment all cells in row `ri` and column `ci` by 1.
#
# Return _the number of cells with odd values_ in the matrix after applying the
# increment to all `indices`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/10/30/e1.png)
#
#
#
#     Input: n = 2, m = 3, indices = [[0,1],[1,1]]
#     Output: 6
#     Explanation: Initial matrix = [[0,0,0],[0,0,0]].
#     After applying first increment it becomes [[1,2,1],[0,1,0]].
#     The final matrix will be [[1,3,1],[1,3,1]] which contains 6 odd numbers.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/10/30/e2.png)
#
#
#
#     Input: n = 2, m = 2, indices = [[1,1],[0,0]]
#     Output: 0
#     Explanation: Final matrix = [[2,2],[2,2]]. There is no odd number in the final matrix.
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 50`
#   * `1 <= m <= 50`
#   * `1 <= indices.length <= 100`
#   * `0 <= indices[i][0] < n`
#   * `0 <= indices[i][1] < m`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@inline function ExaTronKernels.dgpstep(n::Int,x::ROCDeviceArray{Float64,1},xl::ROCDeviceArray{Float64,1},
                         xu::ROCDeviceArray{Float64,1},alpha,w::ROCDeviceArray{Float64,1},
                         s::ROCDeviceArray{Float64,1})
    tx = workitemIdx().x

    if tx <= n
        @inbounds begin
            # It might be better to process this using just a single thread,
            # rather than diverging between multiple threads.

            if x[tx] + alpha*w[tx] < xl[tx]
                s[tx] = xl[tx] - x[tx]
            elseif x[tx] + alpha*w[tx] > xu[tx]
                s[tx] = xu[tx] - x[tx]
            else
                s[tx] = alpha*w[tx]
            end
        end
    end
    AMDGPU.sync_workgroup()

    return
end


text:
# Copyright 2016, 2017, 2018 Martin Holters
# See accompanying license file.

export SimpleSolver, HomotopySolver, CachingSolver
import Base.copy!

struct ParametricNonLinEq{F_eval<:Function,F_setp<:Function,F_calcjp<:Function,Scratch}
    func::F_eval
    set_p::F_setp
    calc_Jp::F_calcjp
    res::Vector{Float64}
    Jp::Matrix{Float64}
    J::Matrix{Float64}
    scratch::Scratch
    function ParametricNonLinEq(func::F_eval, set_p::F_setp, calc_Jp::F_calcjp,
            scratch::Scratch, nn::Integer, np::Integer
        ) where {F_eval<:Function,F_setp<:Function,F_calcjp<:Function,Scratch}
        res = zeros(nn)
        Jp = zeros(nn, np)
        J = zeros(nn, nn)
        return new{F_eval,F_setp,F_calcjp,Scratch}(func, set_p, calc_Jp, res, Jp, J, scratch)
    end
end
ParametricNonLinEq(func::Function, nn::Integer, np::Integer) =
    ParametricNonLinEq(func, default_set_p, default_calc_Jp,
                       (zeros(np), zeros(nn, np)), nn, np)

default_set_p(scratch, p) = (copyto!(scratch[1], p); nothing)
default_calc_Jp(scratch, Jp) = (copyto!(Jp, scratch[2]); nothing)

nn(nleq::ParametricNonLinEq) = length(nleq.res)
np(nleq::ParametricNonLinEq) = size(nleq.Jp, 2)

set_p!(nleq::ParametricNonLinEq, p) = nleq.set_p(nleq.scratch, p)
calc_Jp!(nleq::ParametricNonLinEq) = nleq.calc_Jp(nleq.scratch, nleq.Jp)
evaluate!(nleq::ParametricNonLinEq, z) =
    nleq.func(nleq.res, nleq.J, nleq.scratch, z)

struct LinearSolver
    factors::Matrix{Float64}
    ipiv::Vector{Int}
    function LinearSolver(n::Int)
        new(zeros(n, n), zeros(Int, n))
    end
end

function setlhs!(solver::LinearSolver, A::Matrix{Float64})
    m, n = size(solver.factors)
    if (m, n) ≠ size(A)
        throw(DimensionMismatch("matrix has size $(size(A)), but must have size $(size(solver.factors))"))
    end
    copyto!(solver.factors, A)

    # based on Julia's generic_lufact!, but storing inverses on the diagonal;
    # faster than calling out to dgetrf for sizes up to about 60×60
    factors = solver.factors
    minmn = min(m,n)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            amax = 0.0
            for i = k:m
                absi = abs(factors[i,k])
                if absi > amax
                    kp = i
                    amax = absi
                end
            end
            solver.ipiv[k] = kp
            if factors[kp,k] != 0.0
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = factors[k,i]
                        factors[k,i] = factors[kp,i]
                        factors[kp,i] = tmp
                    end
                end
                # Scale first column
                fkkinv = factors[k,k] = inv(factors[k,k])
                for i = k+1:m
                    factors[i,k] *= fkkinv
                end
            else
                return false
            end
            # Update the rest
            for j = k+1:n
                for i = k+1:m
                    factors[i,j] -= factors[i,k]*factors[k,j]
                end
            end
        end
    end
    return true
end

function solve!(solver::LinearSolver, x::Vector{Float64}, b::Vector{Float64})
    n = size(solver.factors, 2)
    if n ≠ length(x)
        throw(DimensionMismatch("x has length $(length(x)), but needs $n"))
    end
    if x !== b
        if n ≠ length(b)
            throw(DimensionMismatch("b has length $(length(b)), but needs $n"))
        end
        copyto!(x, b)
    end

    # native Julia implementation seems to be faster than dgetrs up to about
    # n=45 (and not slower up to about n=70)
    @inbounds for i in 1:n
        x[i], x[solver.ipiv[i]] = x[solver.ipiv[i]], x[i]
    end
    # taken from Julia's naivesub!(::UnitLowerTriangular, ...)
    @inbounds for j in 1:n
        xj = x[j]
        for i in j+1:n
            x[i] -= solver.factors[i,j] * xj
        end
    end
    # based on Julia's naivesub!(::UpperTriangular, ...), but with factors[j,j]
    # holding inverses
    @inbounds for j in n:-1:1
        xj = x[j] = solver.factors[j,j] * x[j]
        for i in 1:j-1
            x[i] -= solver.factors[i,j] * xj
        end
    end

    return nothing
end

function copy!(dest::LinearSolver, src::LinearSolver)
    copyto!(dest.factors, src.factors)
    copyto!(dest.ipiv, src.ipiv)
end

"""
    SimpleSolver

The `SimpleSolver` is the simplest available solver. It uses Newton iteration
which features fast local convergence, but makes no guarantees about global
convergence. The initial solution of the iteration is obtained by extrapolating
the last solution found (or another solution provided externally) using the
available Jacobians. Due to the missing global convergence, the `SimpleSolver`
is rarely useful as such.
"""
mutable struct SimpleSolver{NLEQ<:ParametricNonLinEq}
    nleq::NLEQ
    z::Vector{Float64}
    linsolver::LinearSolver
    last_z::Vector{Float64}
    last_p::Vector{Float64}
    last_Jp::Matrix{Float64}
    last_linsolver::LinearSolver
    iters::Int
    resmaxabs::Float64
    tol::Float64
    tmp_nn::Vector{Float64}
    tmp_np::Vector{Float64}
    function SimpleSolver(nleq::NLEQ, initial_p::Vector{Float64},
            initial_z::Vector{Float64}) where {NLEQ<:ParametricNonLinEq}
        z = zeros(nn(nleq))
        linsolver = LinearSolver(nn(nleq))
        last_z = zeros(nn(nleq))
        last_p = zeros(np(nleq))
        last_Jp = zeros(nn(nleq), np(nleq))
        last_linsolver = LinearSolver(nn(nleq))
        tmp_nn = zeros(nn(nleq))
        tmp_np = zeros(np(nleq))
        solver = new{NLEQ}(nleq, z, linsolver, last_z, last_p, last_Jp,
            last_linsolver, 0, 0.0, 1e-10, tmp_nn, tmp_np)
        set_extrapolation_origin(solver, initial_p, initial_z)
        return solver
    end
end

set_resabstol!(solver::SimpleSolver, tol) = solver.tol = tol

function set_extrapolation_origin(solver::SimpleSolver, p, z)
    set_p!(solver.nleq, p)
    evaluate!(solver.nleq, z)
    setlhs!(solver.linsolver, solver.nleq.J)
    calc_Jp!(solver.nleq)
    set_extrapolation_origin(solver, p, z, solver.nleq.Jp, solver.linsolver)
end

function set_extrapolation_origin(solver::SimpleSolver, p, z, Jp, linsolver)
    copy!(solver.last_linsolver, linsolver)
    copyto!(solver.last_Jp, Jp)
    copyto!(solver.last_p, p)
    copyto!(solver.last_z, z)
end

get_extrapolation_origin(solver::SimpleSolver) = solver.last_p, solver.last_z

get_extrapolation_jacobian(solver::SimpleSolver) =
    -solver.nleq.J \ solver.nleq.Jp

hasconverged(solver::SimpleSolver) = solver.resmaxabs < solver.tol

needediterations(solver::SimpleSolver) = solver.iters

function solve(solver::SimpleSolver, p::AbstractVector{Float64}, maxiter=500)
    set_p!(solver.nleq, p)
    #solver.z = solver.last_z - solver.last_J\(solver.last_Jp * (p-solver.last_p))
    copyto!(solver.tmp_np, p)
    BLAS.axpy!(-1.0, solver.last_p, solver.tmp_np)
    BLAS.gemv!('N', 1.,solver.last_Jp, solver.tmp_np, 0., solver.tmp_nn)
    solve!(solver.last_linsolver, solver.tmp_nn, solver.tmp_nn)
    copyto!(solver.z, solver.last_z)
    BLAS.axpy!(-1.0, solver.tmp_nn, solver.z)

    for solver.iters=1:maxiter
        evaluate!(solver.nleq, solver.z)
        solver.resmaxabs = isempty(solver.nleq.res) ? 0.0 : maximum(abs, solver.nleq.res)
        if !isfinite(solver.resmaxabs) || !all(isfinite, solver.nleq.J)
            return solver.z
        end
        if !setlhs!(solver.linsolver, solver.nleq.J) # J was singular
            return solver.z
        end
        hasconverged(solver) && break
        #solver.z -= solver.nleq.J\solver.nleq.res
        solve!(solver.linsolver, solver.tmp_nn, solver.nleq.res)
        BLAS.axpy!(-1.0, solver.tmp_nn, solver.z)
    end
    if hasconverged(solver)
        calc_Jp!(solver.nleq)
        set_extrapolation_origin(solver, p, solver.z, solver.nleq.Jp, solver.linsolver)
    end
    return solver.z
end

"""
    HomotopySolver{BaseSolver}

The `HomotopySolver` extends an existing solver (provided as the type parameter)
by applying homotopy to (at least theoretically) ensure global convergence. It
can be combined with the `SimpleSolver` as `HomotopySolver{SimpleSolver}` to
obtain a useful Newton homtopy solver with generally good convergence
properties.
"""
mutable struct HomotopySolver{BaseSolver}
    basesolver::BaseSolver
    start_p::Vector{Float64}
    pa::Vector{Float64}
    iters::Int
    HomotopySolver(basesolver::BaseSolver, np::Integer) where {BaseSolver} =
        new{BaseSolver}(basesolver, zeros(np), zeros(np), 0)
    function HomotopySolver{BaseSolver}(nleq::ParametricNonLinEq,
            initial_p::Vector{Float64}, initial_z::Vector{Float64}
        ) where {BaseSolver}
        basesolver = BaseSolver(nleq, initial_p, initial_z)
        return HomotopySolver(basesolver, np(nleq))
    end
end

set_resabstol!(solver::HomotopySolver, tol) =
    set_resabstol!(solver.basesolver, tol)

set_extrapolation_origin(solver::HomotopySolver, p, z) =
    set_extrapolation_origin(solver.basesolver, p, z)

function solve(solver::HomotopySolver, p)
    z = solve(solver.basesolver, p)
    solver.iters = needediterations(solver.basesolver)
    if !hasconverged(solver)
        a = 0.5
        best_a = 0.0
        copyto!(solver.start_p, get_extrapolation_origin(solver.basesolver)[1])
        while best_a < 1
            # copyto!(solver.pa, (1-a) * solver.start_p + a * p)
            copyto!(solver.pa, solver.start_p)
            rmul!(solver.pa, 1-a)
            axpy!(a, p, solver.pa)
            z = solve(solver.basesolver, solver.pa)
            solver.iters += needediterations(solver.basesolver)
            if hasconverged(solver)
                best_a = a
                a = 1.0
            else
                new_a = (a + best_a) / 2
                if !(best_a < new_a < a)
                    # no floating point value inbetween best_a and a
                    break
                end
                a = new_a
            end
        end
    end
    return z
end

hasconverged(solver::HomotopySolver) = hasconverged(solver.basesolver)
needediterations(solver::HomotopySolver) = solver.iters

get_extrapolation_jacobian(solver::HomotopySolver) =
    get_extrapolation_jacobian(solver.basesolver)

"""
    CachingSolver{BaseSolver}

The `CachingSolver` extends an existing solver (provided as the type parameter)
by storing found solutions in a k-d tree to use as initial solutions in the
future. Whenever the underlying solver needs more than a preset number of
iterations (defaults to five), the solution will be stored. Storing new
solutions is a relatively expensive operation, so until the stored solutions
suffice to ensure convergence in few iterations throughout, use of a
`CachingSolver` may actually slow things down.

See [M. Holters, U. Zölzer, "A k-d Tree Based Solution Cache for the Non-linear
Equation of Circuit Simulations"](http://www.eurasip.org/Proceedings/Eusipco/Eusipco2016/papers/1570255150.pdf)
for a more detailed discussion.
"""
mutable struct CachingSolver{BaseSolver}
    basesolver::BaseSolver
    ps_tree::KDTree{Vector{Float64}, Matrix{Float64}}
    zs::Matrix{Float64}
    num_ps::Int
    new_count::Int
    new_count_limit::Int
    alts::Alts{Float64}
    function CachingSolver(basesolver::BaseSolver, initial_p::Vector{Float64},
            initial_z::Vector{Float64}, nn::Integer) where {BaseSolver}
         ps_tree = KDTree(hcat(initial_p))
         zs = reshape(copy(initial_z), nn, 1)
         alts = Alts(initial_p)
         return new{BaseSolver}(basesolver, ps_tree, zs, 1, 0, 2, alts)
    end
    function CachingSolver{BaseSolver}(nleq::ParametricNonLinEq,
            initial_p::Vector{Float64}, initial_z::Vector{Float64}) where {BaseSolver}
        basesolver = BaseSolver(nleq, initial_p, initial_z)
        return CachingSolver(basesolver, initial_p, initial_z, nn(nleq))
    end
end

set_resabstol!(solver::CachingSolver, tol) =
    set_resabstol!(solver.basesolver, tol)

hasconverged(solver::CachingSolver) = hasconverged(solver.basesolver)
needediterations(solver::CachingSolver) = needediterations(solver.basesolver)

function solve(solver::CachingSolver, p)
    origin_p = get_extrapolation_origin(solver.basesolver)[1]
    best_diff = 0.0
    for i in eachindex(origin_p)
        best_diff += abs2(p[i] - origin_p[i])
    end
    idx = 0
    for i in (solver.num_ps-solver.new_count+1):solver.num_ps
        diff = 0.
        for j in 1:size(solver.ps_tree.ps, 1)
            diff += abs2(solver.ps_tree.ps[j,i] - p[j])
        end
        if diff < best_diff
            best_diff = diff
            idx = i
        end
    end

    init!(solver.alts, best_diff, idx)
    idx = indnearest(solver.ps_tree, p, solver.alts)

    if idx ≠ 0
        set_extrapolation_origin(solver.basesolver,
                                 solver.ps_tree.ps[:,idx], solver.zs[:,idx])
    end

    z = solve(solver.basesolver, p)
    if needediterations(solver.basesolver) > 5 && hasconverged(solver.basesolver)
        solver.num_ps += 1
        if solver.num_ps > size(solver.ps_tree.ps, 2)
            solver.ps_tree.ps =
                copyto!(zeros(size(solver.ps_tree.ps, 1), 2solver.num_ps),
                      solver.ps_tree.ps)
            solver.zs =
                copyto!(zeros(size(solver.zs, 1), 2solver.num_ps), solver.zs)
        end
        solver.ps_tree.ps[:,solver.num_ps] = p
        solver.zs[:,solver.num_ps] = z
        solver.new_count += 1
    end
    if solver.new_count > 0
        solver.new_count_limit -= 1
    end
    if solver.new_count > solver.new_count_limit
        solver.ps_tree = KDTree(solver.ps_tree.ps, solver.num_ps)
        solver.new_count = 0
        solver.new_count_limit = 2size(solver.ps_tree.ps, 2)
    end
    return z
end

get_extrapolation_origin(solver::CachingSolver) =
    get_extrapolation_origin(solver.basesolver)

set_extrapolation_origin(solver::CachingSolver, p, z) =
    set_extrapolation_origin(solver.basesolver, p, z)

get_extrapolation_jacobian(solver::CachingSolver) =
    get_extrapolation_jacobian(solver.basesolver)

function set_resabs2tol!(solver, tol)
    Base.depwarn(string("set_resabs2tol!(solver, tol) is deprecated, use set_resabstol!(solver, sqrt(tol)) instead."),
                 :set_resabs2tol!)
    set_resabstol!(solver, sqrt(tol))
end

function linearize(solver, p::AbstractVector{Float64})
    z = solve(solver, p)
    set_extrapolation_origin(solver, p, z)
    if !hasconverged(solver)
        throw(ArgumentError("Cannot linearize because no solution found at p=$p"))
    end
    return z, get_extrapolation_jacobian(solver)
end


text:
@testset "dense LinearAlgebra" begin
    @testset "dot" begin
        @testset "Vector{$T}" for T in (Float64, ComplexF64)
            test_frule(dot, randn(T, 3), randn(T, 3))
            test_rrule(dot, randn(T, 3), randn(T, 3))
        end
        @testset "Array{$T, 3}" for T in (Float64, ComplexF64)
            test_frule(dot, randn(T, 3, 4, 5), randn(T, 3, 4, 5))
            test_rrule(dot, randn(T, 3, 4, 5), randn(T, 3, 4, 5))
        end
        @testset "mismatched shapes" begin
            # forward
            test_frule(dot, randn(3, 5), randn(5, 3))
            test_frule(dot, randn(15), randn(5, 3))
            # reverse
            test_rrule(dot, randn(3, 5), randn(5, 3))
            test_rrule(dot, randn(15), randn(5, 3))
        end
        @testset "3-arg dot, Array{$T}" for T in (Float64, ComplexF64)
            test_frule(dot, randn(T, 3), randn(T, 3, 4), randn(T, 4))
            test_rrule(dot, randn(T, 3), randn(T, 3, 4), randn(T, 4))
        end
        permuteddimsarray(A) = PermutedDimsArray(A, (2, 1))
        @testset "3-arg dot, $F{$T}" for T in (Float32, ComplexF32),
            F in (adjoint, permuteddimsarray)

            A = F(rand(T, 4, 3)) ⊢ F(rand(T, 4, 3))
            test_frule(dot, rand(T, 3), A, rand(T, 4); rtol=1f-3)
            test_rrule(dot, rand(T, 3), A, rand(T, 4); rtol=1f-3)
        end
        @testset "different types" begin
            test_rrule(dot, rand(2), rand(2, 2), rand(ComplexF64, 2))
            test_rrule(dot, rand(2), Diagonal(rand(2)), rand(ComplexF64, 2))

            # Inference failure due to https://github.com/JuliaDiff/ChainRulesCore.jl/issues/407
            test_rrule(dot, Diagonal(rand(2)), rand(2, 2); check_inferred=false)
        end
    end

    @testset "cross" begin
        test_frule(cross, randn(3), randn(3))
        test_frule(cross, randn(ComplexF64, 3), randn(ComplexF64, 3))
        test_rrule(cross, randn(3), randn(3))
        # No complex support for rrule(cross,...

        # mix types
        test_rrule(cross, rand(3), rand(Float32, 3); rtol=1.0e-7, atol=1.0e-7)
    end
    @testset "pinv" begin
        @testset "$T" for T in (Float64, ComplexF64)
            test_scalar(pinv, randn(T))
            @test frule((ZeroTangent(), randn(T)), pinv, zero(T))[2] ≈ zero(T)
            @test rrule(pinv, zero(T))[2](randn(T))[2] ≈ zero(T)
        end
        @testset "Vector{$T}" for T in (Float64, ComplexF64)
            test_frule(pinv, randn(T, 3), 0.0)
            test_frule(pinv, randn(T, 3), 0.0)

            # Checking types. TODO do we still need this?
            x = randn(T, 3)
            ẋ = randn(T, 3)
            Δy = copyto!(similar(pinv(x)), randn(T, 3))
            @test frule((ZeroTangent(), ẋ), pinv, x)[2] isa typeof(pinv(x))
            @test rrule(pinv, x)[2](Δy)[2] isa typeof(x)
        end

        @testset "$F{Vector{$T}}" for T in (Float64, ComplexF64), F in (Transpose, Adjoint)
            test_frule(pinv, F(randn(T, 3)))
            check_inferred = VERSION ≥ v"1.5"
            test_rrule(pinv, F(randn(T, 3)); check_inferred=check_inferred)

            # Check types.
            # TODO: Do we need this still?
            x, ẋ, x̄ = F(randn(T, 3)), F(randn(T, 3)), F(randn(T, 3))
            y = pinv(x)
            Δy = copyto!(similar(y), randn(T, 3))

            y_fwd, ∂y_fwd = frule((ZeroTangent(), ẋ), pinv, x)
            @test y_fwd isa typeof(y)
            @test ∂y_fwd isa typeof(y)

            y_rev, back = rrule(pinv, x)
            @test y_rev isa typeof(y)
            @test back(Δy)[2] isa typeof(x)
        end
        @testset "Matrix{$T} with size ($m,$n)" for T in (Float64, ComplexF64),
            m in 1:3,
            n in 1:3

            test_frule(pinv, randn(T, m, n))
            test_rrule(pinv, randn(T, m, n))
        end
    end
    @testset "$f" for f in (det, logdet)
        @testset "$f(::$T)" for T in (Float64, ComplexF64)
            b = (f === logdet && T <: Real) ? abs(randn(T)) : randn(T)
            test_scalar(f, b)
        end
        @testset "$f(::Matrix{$T})" for T in (Float64, ComplexF64)
            B = generate_well_conditioned_matrix(T, 4)
            if f === logdet && float(T) <: Float32
                test_frule(f, B; atol=1e-5, rtol=1e-5)
                test_rrule(f, B; atol=1e-5, rtol=1e-5)
            else
                test_frule(f, B)
                test_rrule(f, B)
            end
        end
    end
    @testset "logabsdet(::Matrix{$T})" for T in (Float64, ComplexF64)
        B = randn(T, 4, 4)
        test_frule(logabsdet, B)
        test_rrule(logabsdet, B)
        # test for opposite sign of determinant
        test_frule(logabsdet, -B)
        test_rrule(logabsdet, -B)
    end
    @testset "tr" begin
        test_frule(tr, randn(4, 4))
        test_rrule(tr, randn(4, 4))
    end
    @testset "sylvester" begin
        @testset "T=$T, m=$m, n=$n" for T in (Float64, ComplexF64), m in (2, 3), n in (1, 3)
            A = randn(T, m, m)
            B = randn(T, n, n)
            C = randn(T, m, n)
            test_frule(sylvester, A, B, C)
            test_rrule(sylvester, A, B, C)
        end
    end
    @testset "lyap" begin
        n = 3
        @testset "Float64" for T in (Float64, ComplexF64)
            A = randn(T, n, n)
            C = randn(T, n, n)
            test_frule(lyap, A, C)
            test_rrule(lyap, A, C)
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Libglvnd"
version = v"1.3.0"

# Collection of sources required to build Libglvnd
sources = [
    ArchiveSource("https://gitlab.freedesktop.org/glvnd/libglvnd/uploads/d164b4e6bed74290f4d60e9a5b9bc31c/libglvnd-$(version).tar.gz",
                  "0f43bd9f6c20e6a75ff8cd57736bd78071da7a68e078ed39c81bcc710af30dc7"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libglvnd-*/
export CPPFLAGS="-I${prefix}/include"
./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
# The license is embedded in the README file
install_license README.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if p isa Union{Linux,FreeBSD}]

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libEGL", :libEGL),
    LibraryProduct("libGLdispatch", :libGLdispatch),
    LibraryProduct("libGLESv1_CM", :libGLESv1_CM),
    LibraryProduct("libGLESv2", :libGLESv2),
    LibraryProduct("libGL", :libGL),
    LibraryProduct("libGLX", :libGLX),
    LibraryProduct("libOpenGL", :libOpenGL),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Xorg_libX11_jll"),
    Dependency("Xorg_libXext_jll"),
    BuildDependency("Xorg_glproto_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
nchild = 4
keylist = Int[i for i=1:20]

# Chunk up the keylist into an acceptable number of keys per process.
max_keys_per_process = iceil(length(keylist)/nchild)
partial_proc = max_keys_per_process * nchild - length(keylist) # num of unfull processors
full_proc = nchild - partial_proc # number of full processors

# Feed max_keys_per_process keys to each of the full_proc, then feed
# (max_keys_per_process - 1) to each of the partial_proc
chunk_keylist = Array(Any, nchild)
start_key = 1
end_key = start_key + (max_keys_per_process - 1)
for i=1:full_proc
    println("Choosing from $start_key to $end_key")
    chunk_keylist[i] = keylist[start_key:end_key]
    start_key += max_keys_per_process
    end_key += max_keys_per_process
end

# Inch back to prep for filling the partially filled processes
end_key -= 1
# Now fill the partial_proc
for i=(full_proc+1):nchild
    println("Choosing from $start_key to $end_key")
    chunk_keylist[i] = keylist[start_key:end_key]
    start_key += (max_keys_per_process - 1)
    end_key += (max_keys_per_process - 1)
end

println(chunk_keylist)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
typealias SCIP_Vartype Int8

SCIP_VARTYPE_BINARY     = int8( 0)
SCIP_VARTYPE_INTEGER    = int8( 1)
SCIP_VARTYPE_IMPLINT    = int8( 2)
SCIP_VARTYPE_CONTINUOUS = int8( 3)

SCIP_VARTYPE = Dict{SCIP_Vartype, String}()
SCIP_VARTYPE[SCIP_VARTYPE_BINARY]     = "binary variable: x in {0,1}"
SCIP_VARTYPE[SCIP_VARTYPE_INTEGER]    = "integer variable: x in {lb, ..., ub}"
SCIP_VARTYPE[SCIP_VARTYPE_IMPLINT]    = "implicit integer variable: continuous variable, that is always integral"
SCIP_VARTYPE[SCIP_VARTYPE_CONTINUOUS] = "continuous variable: x in [lb,ub]"


text:
using Test

@testset "DeepcopyModules.jl" begin
  include("DeepcopyModules.jl")
end


text:
using Base.Test
using FactCheck
# test linear system solvers

function test_linear_solver(solver::abstract_linear_system_solver, sym::Symbol)
    start_advanced_timer()
    n = 40;
    mat = speye(n)
    mat = mat + mat'
    mat[1,1] = -100;
    mat[2,2] = -100;

    initialize!(solver, mat, sym);

    test_linear_solver_inertia(solver, sym, mat, 1, 38, 2)

    rhs = 1.0*rand(n, 2);
    A = rand(n,n)
    mat[:,:] = A + A'

    test_linear_solver_accuracy(solver, sym, mat, rhs, 38, 2)

    finalize!(solver)
    pause_advanced_timer()
end

function test_linear_solver_inertia(solver::abstract_linear_system_solver, sym::Symbol, mat::SparseMatrixCSC{Float64,Int64}, correct_inertia::Int64, n::Int64, m::Int64)

   @test solver._SparseMatrix == mat

   inertia = ls_factor!(solver, n, m)
   @fact inertia --> correct_inertia
end

function test_linear_solver_accuracy(solver::abstract_linear_system_solver, sym::Symbol, mat::SparseMatrixCSC{Float64,Int64}, rhs::AbstractArray, n::Int64, m::Int64)
  ls_factor!(solver, n, m)

  sol = ls_solve(solver,rhs);

  #sol = ls_solve(solver,rhs);
  err = norm(mat*sol - rhs,1)/norm(sol,1)
  @fact err --> less_than(1e-6)
end

facts("Linear solvers") do

  # test mumps solver
  context("MUMPS") do
    ls_solver_mumps = linear_solver_MUMPS();
    test_linear_solver( ls_solver_mumps, :symmetric )
    #MPI.Finalize()
  end

  # test julia solver
  context("julia") do
    # test julia lu factor
    ls_solver_julia = linear_solver_JULIA();
    test_linear_solver( ls_solver_julia , :unsymmetric )

    # test julia cholesky factor
    ls_solver_julia = linear_solver_JULIA();
    test_linear_solver( ls_solver_julia, :definite )
  end

end

# test matlab ldl solver
#context("matlab") do
#  ls_solver_matlab = linear_solver_MATLAB();
#  test_linear_solver( ls_solver_matlab, :symmetric )
#end

max_forks_repo_name: ohinder/Everest


text:
for (T, msg) in zip([:PhysicalDeviceFeatures, :PhysicalDeviceVulkan11Features, :PhysicalDeviceVulkan12Features, :PhysicalDeviceVulkan13Features], ["", " Vulkan 1.1", " Vulkan 1.2", " Vulkan 1.3"])
    hasnext = T ≠ :PhysicalDeviceFeatures
    fdecl = hasnext ? :($T(features::AbstractArray; next=C_NULL)) : :($T(features::AbstractArray))
    call = hasnext ? :($T(args...; next)) : :($T(args...))
    @eval $fdecl = begin
        names = filter(!in((:s_type, :next)), fieldnames($T))
        diff = setdiff(features, names)
        if length(diff) ≠ 0
            error(string("Invalid physical device", $msg, " features: ", join(diff, ", ")))
        end
        args = map(in(features), names)
        $call
    end
end

"""
Return a `PhysicalDeviceFeatures` object with the provided `features` set to true.

```jldoctest
julia> PhysicalDeviceFeatures()
PhysicalDeviceFeatures()

julia> PhysicalDeviceFeatures(:wide_lines, :sparse_binding)
PhysicalDeviceFeatures(wide_lines, sparse_binding)
```
"""
PhysicalDeviceFeatures(features::Symbol...) = PhysicalDeviceFeatures(collect(features))

"""
Return a `PhysicalDeviceVulkan11Features` object with the provided `features` set to true.

```jldoctest
julia> PhysicalDeviceVulkan11Features(; next = C_NULL)
PhysicalDeviceVulkan11Features(next=Ptr{Nothing} @0x0000000000000000)

julia> PhysicalDeviceVulkan11Features(:multiview, :variable_pointers, next = C_NULL)
PhysicalDeviceVulkan11Features(next=Ptr{Nothing} @0x0000000000000000, multiview, variable_pointers)
```
"""
PhysicalDeviceVulkan11Features(features::Symbol...; next = C_NULL) = PhysicalDeviceVulkan11Features(collect(features); next)

"""
Return a `PhysicalDeviceVulkan12Features` object with the provided `features` set to true.

```jldoctest
julia> PhysicalDeviceVulkan12Features(; next = C_NULL)
PhysicalDeviceVulkan12Features(next=Ptr{Nothing} @0x0000000000000000)

julia> PhysicalDeviceVulkan12Features(:draw_indirect_count, :descriptor_binding_variable_descriptor_count)
PhysicalDeviceVulkan12Features(next=Ptr{Nothing} @0x0000000000000000, draw_indirect_count, descriptor_binding_variable_descriptor_count)
```
"""
PhysicalDeviceVulkan12Features(features::Symbol...; next = C_NULL) = PhysicalDeviceVulkan12Features(collect(features); next)

"""
Return a `PhysicalDeviceVulkan13Features` object with the provided `features` set to true.

```jldoctest
julia> PhysicalDeviceVulkan13Features(; next = C_NULL)
PhysicalDeviceVulkan13Features(next=Ptr{Nothing} @0x0000000000000000)

julia> PhysicalDeviceVulkan13Features(:dynamic_rendering, :descriptor_binding_variable_descriptor_count)
PhysicalDeviceVulkan13Features(next=Ptr{Nothing} @0x0000000000000000, dynamic_rendering)
```
"""
PhysicalDeviceVulkan13Features(features::Symbol...; next = C_NULL) = PhysicalDeviceVulkan13Features(collect(features); next)

"""
Find a queue index (starting at 0) from `physical_device` which matches the provided `queue_capabilities`.

```jldoctest
julia> find_queue_family(physical_device, QUEUE_COMPUTE_BIT & QUEUE_GRAPHICS_BIT)
0
```
"""
function find_queue_family(physical_device::PhysicalDevice, queue_capabilities::QueueFlag)
    qf_props = get_physical_device_queue_family_properties(physical_device)
    index = findfirst(x -> queue_capabilities in x.queue_flags, qf_props)
    if isnothing(index)
        error("No queue with the desired capabilities could be found.")
    end
    index - 1
end

device(device::Device)::Device = device
device(handle::Handle)::Device = device(parent(handle))

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
# Autogenerated wrapper script for embree_jll for x86_64-apple-darwin
export libembree3

using oneTBB_jll
using GLFW_jll
using JpegTurbo_jll
using libpng_jll
JLLWrappers.@generate_wrapper_header("embree")
JLLWrappers.@declare_library_product(libembree3, "@rpath/libembree3.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header(oneTBB_jll, GLFW_jll, JpegTurbo_jll, libpng_jll)
    JLLWrappers.@init_library_product(
        libembree3,
        "lib/libembree3.3.13.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()