text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Random123: Philox4x, Threefry4x
using Random123: random123_r, gen_seed

"""
    bat_rng()

*Experimental feature, not part of stable public API.*

Return a new BAT-compatible random number generator, with a random seed drawn
from the system entropy pool.
"""
function bat_rng end
export bat_rng

bat_rng() = Philox4x()::Philox4x{UInt64,10}

bat_determ_rng() = Philox4x((0, 0))::Philox4x{UInt64,10}

struct RNGPartition{R<:AbstractRNG,S,C,I}
    seed::S
    partctrsbase::C
    depth::Int
    idxs::I
end

function RNGPartition(parent_rng::AbstractRNG, partidxs::AbstractUnitRange{<:Integer})
    seed = rngpart_getseed(parent_rng)
    partcounters, depth = rngpart_getpartctrs(parent_rng)

    child_rngpart_partctrsbase = _rngpart_inc_partctrs(partcounters, depth, 1)

    parent_mod_partcounters = _rngpart_inc_partctrs(child_rngpart_partctrsbase, depth, length(partidxs))
    rngpart_setpartctrs!(parent_rng, parent_mod_partcounters, depth)

    R = typeof(parent_rng)
    S = typeof(seed)
    C = typeof(child_rngpart_partctrsbase)
    I = typeof(partidxs)
    RNGPartition{R,S,C,I}(seed, child_rngpart_partctrsbase, depth, partidxs)
end

Base.eachindex(rngpart::RNGPartition) = rngpart.idxs
Base.length(rngpart::RNGPartition) = length(eachindex(rngpart))
Base.size(rngpart::RNGPartition) = (length(rngpart),)

function set_rng!(rng::R, rngpart::RNGPartition{R}, i::Integer) where R <: AbstractRNG
    idxs = eachindex(rngpart)
    Base.checkindex(Bool, idxs, i) || throw(ArgumentError("Index $i not in partition indices $idxs of $rngpart"))

    j = i - minimum(idxs)
    mod_partcounters = _rngpart_inc_partctrs(rngpart.partctrsbase, rngpart.depth, j)
    mod_depth = rngpart.depth + 1

    Random.seed!(rng, rngpart.seed)
    rngpart_setpartctrs!(rng, mod_partcounters, mod_depth)

    rng
end

Random.AbstractRNG(rngpart::RNGPartition{R}, i::Integer) where R =
    set_rng!(rngpart_createrng(R), rngpart, i)

rngpart_createrng(::Type{T}) where {T <: Philox4x} = T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

rngpart_getseed(rng::Philox4x) = (rng.key1, rng.key2)

# function set_rng!(target::Philox4x, src::Philox4x)
#     target.ctr1 = src.ctr1
#     target.ctr2 = src.ctr2
#     target.ctr3 = src.ctr3
#     target.ctr4 = src.ctr4
#     target.key1 = src.key1
#     target.key2 = src.key2
#     target.p = src.p
#     target.x1 = src.x1
#     target.x2 = src.x2
#     target.x3 = src.x3
#     target.x4 = src.x4
#     target
# end

rngpart_createrng(::Type{T}) where {T <: Threefry4x} = T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

rngpart_getseed(rng::Threefry4x) = (rng.key1, rng.key2, rng.key3, rng.key4)

# function set_rng!(target::Threefry4x, src::Threefry4x)
#     target.x1 = src.x1
#     target.x2 = src.x2
#     target.x3 = src.x3
#     target.x4 = src.x4
#     target.key1 = src.key1
#     target.key2 = src.key2
#     target.key3 = src.key3
#     target.key4 = src.key4
#     target.ctr1 = src.ctr1
#     target.ctr2 = src.ctr2
#     target.ctr3 = src.ctr3
#     target.ctr4 = src.ctr4
#     target.p = src.p
# end

function rngpart_getpartctrs(rng::Union{Philox4x,Threefry4x})
    partctrinfo = _rngpart_split_uints(rng.ctr4, rng.ctr3, rng.ctr2)
    depth = _rngpart_getdepth(partctrinfo)
    partctrs = map(_rngpart_getpartctr, partctrinfo)
    (partctrs = partctrs, depth = depth)
end

function rngpart_setpartctrs!(rng::Union{Philox4x,Threefry4x}, partctrs::NTuple{6,UInt32}, depth::Integer)
    tagged_partctrs = _rngpart_settopbit(partctrs, depth)
    merged_partctrs = _rngpart_merge_uints(tagged_partctrs...)::NTuple{3,UInt64}

    rng.ctr4 = merged_partctrs[1]
    rng.ctr3 = merged_partctrs[2]
    rng.ctr2 = merged_partctrs[3]
    rng.ctr1 = 0

    random123_r(rng)

    rng
end

@inline _rngpart_split_uints() = ()
@inline _rngpart_split_uints(x::UInt64, xs::UInt64...) =
    (UInt32(x >> 32), UInt32(x << 32 >> 32), _rngpart_split_uints(xs...)...)

@inline _rngpart_merge_uints() = ()
@inline _rngpart_merge_uints(x1::UInt32, x2::UInt32, xs::UInt32...) =
    (UInt64(x1) << 32 | UInt64(x2), _rngpart_merge_uints(xs...)...)

_rngpart_topbit_mask(::Type{T}) where {T <: Unsigned} = ((typemax(T) >> 1) + one(T))

_rngpart_lowbits_mask(::Type{T}) where {T <: Unsigned} = typemax(T) >> 1

_rngpart_haspartctrtag(x::T) where {T<:Unsigned} = (x & _rngpart_topbit_mask(T)) > 0

_rngpart_getpartctr(x::T) where {T<:Unsigned} = (x & _rngpart_lowbits_mask(T))

function _rngpart_getdepth(partctrinfo::NTuple{N,T}) where {N,T<:Unsigned}
    cycle::Int = 1
    for i in eachindex(partctrinfo)
        x = partctrinfo[i]
        if _rngpart_haspartctrtag(x)
            cycle = i
        else
            x == 0 || throw(ArgumentError("Inconsistent partition counter information"))
        end
    end
    return cycle
end

function _rngpart_inc_partctrs(partctrs::NTuple{N,T}, depth::Integer, x::Integer) where {N,T<:Unsigned}
    1 <= depth <= length(partctrs) || throw(ArgumentError("Partition depth out of allowed range"))
    m = ntuple(i -> i == depth ? T(x) : zero(T), Val(length(partctrs)))
    partctrs .+ m
end

function _rngpart_settopbit(partctrs::NTuple{N,T}, depth::Integer) where {N,T<:Unsigned}
    1 <= depth <= length(partctrs) || throw(ArgumentError("Partition depth out of allowed range"))
    any(_rngpart_haspartctrtag, partctrs) && throw(ArgumentError("Partition counter(s) out of allowed range"))
    m = ntuple(i -> i <= depth ? _rngpart_topbit_mask(T) : zero(T), Val(length(partctrs)))
    partctrs .| m
end


text:
include("C:/Users/cgamboa/Documents/MIPSolver2017_08.jl/src/branch_and_bound.jl")
include("C:/Users/cgamboa/Documents/MIPTests.jl/miptests.jl")

function runtest(BNB::Function, solver::MathProgBase.AbstractMathProgSolver = JuMP.UnsetSolver())
@testset "Testes" begin
solver = GurobiSolver(OutputFlag=0)

      test1(BNB,solver)
      test2(BNB,solver)
      test3(BNB,solver)
      test3_2(BNB,solver)
      test3_3(BNB,solver)
      testCaminho(BNB,solver)
      testInfeasibleKnapsack(BNB,solver)
      testInfeasibleUC(BNB,solver)
      testRobustCCUC(BNB,solver)
      testSudoku(BNB,solver)
      testSudoku4x4(BNB,solver)
      testUnboundedKnapsack(BNB,solver)
      test_Minimal_UC(BNB,solver)
      test_P1_Brito(BNB,solver)
      test_PL_Infeasible_Brito(BNB,solver)
      test_PL_Infeasible_Raphael(BNB,solver)
      test_PL_Simples_Brito(BNB,solver)
      test_PL_Simples_Raphael(BNB,solver)
      test_PL_Unbounded_Brito(BNB,solver)
      test_feature_selection_grande(BNB,solver)
      test_feature_selection_medio(BNB,solver)
      test_feature_selection_pequeno_inviavel(BNB,solver)
      test_feature_selection_pequeno_viavel(BNB,solver)
      teste_PL_andrew_inviavel(BNB,solver)
      teste_PL_andrew_unbounded(BNB,solver)
      teste_PL_andrew_viavel(BNB,solver)
    end
end

runtest(BNB,GurobiSolver(OutputFlag=0))

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

abstract type AbstractCartesianIndex{N} end # This is a hacky forward declaration for CartesianIndex
const ViewIndex = Union{Real, AbstractArray}
const ScalarIndex = Real

"""
    SubArray{T,N,P,I,L} <: AbstractArray{T,N}

`N`-dimensional view into a parent array (of type `P`) with an element type `T`, restricted by a tuple of indices (of type `I`). `L` is true for types that support fast linear indexing, and `false` otherwise.

Construct `SubArray`s using the [`view`](@ref) function.
"""
struct SubArray{T,N,P,I,L} <: AbstractArray{T,N}
    parent::P
    indices::I
    offset1::Int       # for linear indexing and pointer, only valid when L==true
    stride1::Int       # used only for linear indexing
    function SubArray{T,N,P,I,L}(parent, indices, offset1, stride1) where {T,N,P,I,L}
        @_inline_meta
        check_parent_index_match(parent, indices)
        new(parent, indices, offset1, stride1)
    end
end
# Compute the linear indexability of the indices, and combine it with the linear indexing of the parent
function SubArray(parent::AbstractArray, indices::Tuple)
    @_inline_meta
    SubArray(IndexStyle(viewindexing(indices), IndexStyle(parent)), parent, ensure_indexable(indices), index_dimsum(indices...))
end
function SubArray(::IndexCartesian, parent::P, indices::I, ::NTuple{N,Any}) where {P,I,N}
    @_inline_meta
    SubArray{eltype(P), N, P, I, false}(parent, indices, 0, 0)
end
function SubArray(::IndexLinear, parent::P, indices::I, ::NTuple{N,Any}) where {P,I,N}
    @_inline_meta
    # Compute the stride and offset
    stride1 = compute_stride1(parent, indices)
    SubArray{eltype(P), N, P, I, true}(parent, indices, compute_offset1(parent, stride1, indices), stride1)
end

check_parent_index_match(parent, indices) = check_parent_index_match(parent, index_ndims(indices...))
check_parent_index_match(parent::AbstractArray{T,N}, ::NTuple{N, Bool}) where {T,N} = nothing
check_parent_index_match(parent, ::NTuple{N, Bool}) where {N} =
    throw(ArgumentError("number of indices ($N) must match the parent dimensionality ($(ndims(parent)))"))

# This computes the linear indexing compatibility for a given tuple of indices
viewindexing(I::Tuple{}) = IndexLinear()
# Leading scalar indices simply increase the stride
viewindexing(I::Tuple{ScalarIndex, Vararg{Any}}) = (@_inline_meta; viewindexing(tail(I)))
# Slices may begin a section which may be followed by any number of Slices
viewindexing(I::Tuple{Slice, Slice, Vararg{Any}}) = (@_inline_meta; viewindexing(tail(I)))
# A UnitRange can follow Slices, but only if all other indices are scalar
viewindexing(I::Tuple{Slice, AbstractUnitRange, Vararg{ScalarIndex}}) = IndexLinear()
viewindexing(I::Tuple{Slice, Slice, Vararg{ScalarIndex}}) = IndexLinear() # disambiguate
# In general, ranges are only fast if all other indices are scalar
viewindexing(I::Tuple{AbstractRange, Vararg{ScalarIndex}}) = IndexLinear()
# All other index combinations are slow
viewindexing(I::Tuple{Vararg{Any}}) = IndexCartesian()
# Of course, all other array types are slow
viewindexing(I::Tuple{AbstractArray, Vararg{Any}}) = IndexCartesian()

# Simple utilities
size(V::SubArray) = (@_inline_meta; map(n->Int(unsafe_length(n)), axes(V)))

similar(V::SubArray, T::Type, dims::Dims) = similar(V.parent, T, dims)

sizeof(V::SubArray) = length(V) * sizeof(eltype(V))
sizeof(V::SubArray{<:Any,<:Any,<:Array}) = length(V) * elsize(V.parent)

copy(V::SubArray) = V.parent[V.indices...]

parent(V::SubArray) = V.parent
parentindices(V::SubArray) = V.indices

"""
    parentindices(A)

Return the indices in the [`parent`](@ref) which correspond to the array view `A`.

# Examples
```jldoctest
julia> A = [1 2; 3 4];

julia> V = view(A, 1, :)
2-element view(::Matrix{Int64}, 1, :) with eltype Int64:
 1
 2

julia> parentindices(V)
(1, Base.Slice(Base.OneTo(2)))
```
"""
parentindices(a::AbstractArray) = map(OneTo, size(a))

## Aliasing detection
dataids(A::SubArray) = (dataids(A.parent)..., _splatmap(dataids, A.indices)...)
_splatmap(f, ::Tuple{}) = ()
_splatmap(f, t::Tuple) = (f(t[1])..., _splatmap(f, tail(t))...)
unaliascopy(A::SubArray) = typeof(A)(unaliascopy(A.parent), map(unaliascopy, A.indices), A.offset1, A.stride1)

# When the parent is an Array we can trim the size down a bit. In the future this
# could possibly be extended to any mutable array.
function unaliascopy(V::SubArray{T,N,A,I,LD}) where {T,N,A<:Array,I<:Tuple{Vararg{Union{Real,AbstractRange,Array}}},LD}
    dest = Array{T}(undef, index_lengths(V.indices...))
    copyto!(dest, V)
    SubArray{T,N,A,I,LD}(dest, map(_trimmedindex, V.indices), 0, Int(LD))
end
# Transform indices to be "dense"
_trimmedindex(i::Real) = oftype(i, 1)
_trimmedindex(i::AbstractUnitRange) = oftype(i, OneTo(length(i)))
_trimmedindex(i::AbstractArray) = oftype(i, reshape(eachindex(IndexLinear(), i), axes(i)))

## SubArray creation
# We always assume that the dimensionality of the parent matches the number of
# indices that end up getting passed to it, so we store the parent as a
# ReshapedArray view if necessary. The trouble is that arrays of `CartesianIndex`
# can make the number of effective indices not equal to length(I).
_maybe_reshape_parent(A::AbstractArray, ::NTuple{1, Bool}) = reshape(A, Val(1))
_maybe_reshape_parent(A::AbstractArray{<:Any,1}, ::NTuple{1, Bool}) = reshape(A, Val(1))
_maybe_reshape_parent(A::AbstractArray{<:Any,N}, ::NTuple{N, Bool}) where {N} = A
_maybe_reshape_parent(A::AbstractArray, ::NTuple{N, Bool}) where {N} = reshape(A, Val(N))
"""
    view(A, inds...)

Like [`getindex`](@ref), but returns a lightweight array that lazily references
(or is effectively a _view_ into) the parent array `A` at the given index or indices
`inds` instead of eagerly extracting elements or constructing a copied subset.
Calling [`getindex`](@ref) or [`setindex!`](@ref) on the returned value
(often a [`SubArray`](@ref)) computes the indices to access or modify the
parent array on the fly.  The behavior is undefined if the shape of the parent array is
changed after `view` is called because there is no bound check for the parent array; e.g.,
it may cause a segmentation fault.

Some immutable parent arrays (like ranges) may choose to simply
recompute a new array in some circumstances instead of returning
a `SubArray` if doing so is efficient and provides compatible semantics.

# Examples
```jldoctest
julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> b = view(A, :, 1)
2-element view(::Matrix{Int64}, :, 1) with eltype Int64:
 1
 3

julia> fill!(b, 0)
2-element view(::Matrix{Int64}, :, 1) with eltype Int64:
 0
 0

julia> A # Note A has changed even though we modified b
2×2 Matrix{Int64}:
 0  2
 0  4

julia> view(2:5, 2:3) # returns a range as type is immutable
3:4
```
"""
function view(A::AbstractArray, I::Vararg{Any,N}) where {N}
    @_inline_meta
    J = map(i->unalias(A,i), to_indices(A, I))
    @boundscheck checkbounds(A, J...)
    unsafe_view(_maybe_reshape_parent(A, index_ndims(J...)), J...)
end

# Ranges implement getindex to return recomputed ranges; use that for views, too (when possible)
function view(r1::OneTo, r2::OneTo)
    @_propagate_inbounds_meta
    getindex(r1, r2)
end
function view(r1::AbstractUnitRange, r2::AbstractUnitRange{<:Integer})
    @_propagate_inbounds_meta
    getindex(r1, r2)
end
function view(r1::AbstractUnitRange, r2::StepRange{<:Integer})
    @_propagate_inbounds_meta
    getindex(r1, r2)
end
function view(r1::StepRange, r2::AbstractRange{<:Integer})
    @_propagate_inbounds_meta
    getindex(r1, r2)
end
function view(r1::StepRangeLen, r2::OrdinalRange{<:Integer})
    @_propagate_inbounds_meta
    getindex(r1, r2)
end
function view(r1::LinRange, r2::OrdinalRange{<:Integer})
    @_propagate_inbounds_meta
    getindex(r1, r2)
end

function unsafe_view(A::AbstractArray, I::Vararg{ViewIndex,N}) where {N}
    @_inline_meta
    SubArray(A, I)
end
# When we take the view of a view, it's often possible to "reindex" the parent
# view's indices such that we can "pop" the parent view and keep just one layer
# of indirection. But we can't always do this because arrays of `CartesianIndex`
# might span multiple parent indices, making the reindex calculation very hard.
# So we use _maybe_reindex to figure out if there are any arrays of
# `CartesianIndex`, and if so, we punt and keep two layers of indirection.
unsafe_view(V::SubArray, I::Vararg{ViewIndex,N}) where {N} =
    (@_inline_meta; _maybe_reindex(V, I))
_maybe_reindex(V, I) = (@_inline_meta; _maybe_reindex(V, I, I))
_maybe_reindex(V, I, ::Tuple{AbstractArray{<:AbstractCartesianIndex}, Vararg{Any}}) =
    (@_inline_meta; SubArray(V, I))
# But allow arrays of CartesianIndex{1}; they behave just like arrays of Ints
_maybe_reindex(V, I, A::Tuple{AbstractArray{<:AbstractCartesianIndex{1}}, Vararg{Any}}) =
    (@_inline_meta; _maybe_reindex(V, I, tail(A)))
_maybe_reindex(V, I, A::Tuple{Any, Vararg{Any}}) = (@_inline_meta; _maybe_reindex(V, I, tail(A)))
function _maybe_reindex(V, I, ::Tuple{})
    @_inline_meta
    @inbounds idxs = to_indices(V.parent, reindex(V.indices, I))
    SubArray(V.parent, idxs)
end

## Re-indexing is the heart of a view, transforming A[i, j][x, y] to A[i[x], j[y]]
#
# Recursively look through the heads of the parent- and sub-indices, considering
# the following cases:
# * Parent index is array  -> re-index that with one or more sub-indices (one per dimension)
# * Parent index is Colon  -> just use the sub-index as provided
# * Parent index is scalar -> that dimension was dropped, so skip the sub-index and use the index as is

AbstractZeroDimArray{T} = AbstractArray{T, 0}

reindex(::Tuple{}, ::Tuple{}) = ()

# Skip dropped scalars, so simply peel them off the parent indices and continue
reindex(idxs::Tuple{ScalarIndex, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1], reindex(tail(idxs), subidxs)...))

# Slices simply pass their subindices straight through
reindex(idxs::Tuple{Slice, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (subidxs[1], reindex(tail(idxs), tail(subidxs))...))

# Re-index into parent vectors with one subindex
reindex(idxs::Tuple{AbstractVector, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1]], reindex(tail(idxs), tail(subidxs))...))

# Parent matrices are re-indexed with two sub-indices
reindex(idxs::Tuple{AbstractMatrix, Vararg{Any}}, subidxs::Tuple{Any, Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1], subidxs[2]], reindex(tail(idxs), tail(tail(subidxs)))...))

# In general, we index N-dimensional parent arrays with N indices
@generated function reindex(idxs::Tuple{AbstractArray{T,N}, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) where {T,N}
    if length(subidxs.parameters) >= N
        subs = [:(subidxs[$d]) for d in 1:N]
        tail = [:(subidxs[$d]) for d in N+1:length(subidxs.parameters)]
        :(@_propagate_inbounds_meta; (idxs[1][$(subs...)], reindex(tail(idxs), ($(tail...),))...))
    else
        :(throw(ArgumentError("cannot re-index SubArray with fewer indices than dimensions\nThis should not occur; please submit a bug report.")))
    end
end

# In general, we simply re-index the parent indices by the provided ones
SlowSubArray{T,N,P,I} = SubArray{T,N,P,I,false}
function getindex(V::SubArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[reindex(V.indices, I)...]
    r
end

# But SubArrays with fast linear indexing pre-compute a stride and offset
FastSubArray{T,N,P,I} = SubArray{T,N,P,I,true}
function getindex(V::FastSubArray, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + V.stride1*i]
    r
end
# We can avoid a multiplication if the first parent index is a Colon or AbstractUnitRange,
# or if all the indices are scalars, i.e. the view is for a single value only
FastContiguousSubArray{T,N,P,I<:Union{Tuple{Union{Slice, AbstractUnitRange}, Vararg{Any}},
                                      Tuple{Vararg{ScalarIndex}}}} = SubArray{T,N,P,I,true}
function getindex(V::FastContiguousSubArray, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + i]
    r
end
# For vector views with linear indexing, we disambiguate to favor the stride/offset
# computation as that'll generally be faster than (or just as fast as) re-indexing into a range.
function getindex(V::FastSubArray{<:Any, 1}, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + V.stride1*i]
    r
end
function getindex(V::FastContiguousSubArray{<:Any, 1}, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + i]
    r
end

# Indexed assignment follows the same pattern as `getindex` above
function setindex!(V::SubArray{T,N}, x, I::Vararg{Int,N}) where {T,N}
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds V.parent[reindex(V.indices, I)...] = x
    V
end
function setindex!(V::FastSubArray, x, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + V.stride1*i] = x
    V
end
function setindex!(V::FastContiguousSubArray, x, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + i] = x
    V
end
function setindex!(V::FastSubArray{<:Any, 1}, x, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + V.stride1*i] = x
    V
end
function setindex!(V::FastContiguousSubArray{<:Any, 1}, x, i::Int)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + i] = x
    V
end

IndexStyle(::Type{<:FastSubArray}) = IndexLinear()
IndexStyle(::Type{<:SubArray}) = IndexCartesian()

# Strides are the distance in memory between adjacent elements in a given dimension
# which we determine from the strides of the parent
strides(V::SubArray) = substrides(strides(V.parent), V.indices)

substrides(strds::Tuple{}, ::Tuple{}) = ()
substrides(strds::NTuple{N,Int}, I::Tuple{ScalarIndex, Vararg{Any}}) where N = (substrides(tail(strds), tail(I))...,)
substrides(strds::NTuple{N,Int}, I::Tuple{Slice, Vararg{Any}}) where N = (first(strds), substrides(tail(strds), tail(I))...)
substrides(strds::NTuple{N,Int}, I::Tuple{AbstractRange, Vararg{Any}}) where N = (first(strds)*step(I[1]), substrides(tail(strds), tail(I))...)
substrides(strds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubArrays with indices of type $(typeof(I[1]))"))

stride(V::SubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]

compute_stride1(parent::AbstractArray, I::NTuple{N,Any}) where {N} =
    (@_inline_meta; compute_stride1(1, fill_to_length(axes(parent), OneTo(1), Val(N)), I))
compute_stride1(s, inds, I::Tuple{}) = s
compute_stride1(s, inds, I::Tuple{Vararg{ScalarIndex}}) = s
compute_stride1(s, inds, I::Tuple{ScalarIndex, Vararg{Any}}) =
    (@_inline_meta; compute_stride1(s*unsafe_length(inds[1]), tail(inds), tail(I)))
compute_stride1(s, inds, I::Tuple{AbstractRange, Vararg{Any}}) = s*step(I[1])
compute_stride1(s, inds, I::Tuple{Slice, Vararg{Any}}) = s
compute_stride1(s, inds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("invalid strided index type $(typeof(I[1]))"))

elsize(::Type{<:SubArray{<:Any,<:Any,P}}) where {P} = elsize(P)

iscontiguous(A::SubArray) = iscontiguous(typeof(A))
iscontiguous(::Type{<:SubArray}) = false
iscontiguous(::Type{<:FastContiguousSubArray}) = true

first_index(V::FastSubArray) = V.offset1 + V.stride1 # cached for fast linear SubArrays
function first_index(V::SubArray)
    P, I = parent(V), V.indices
    s1 = compute_stride1(P, I)
    s1 + compute_offset1(P, s1, I)
end

# Computing the first index simply steps through the indices, accumulating the
# sum of index each multiplied by the parent's stride.
# The running sum is `f`; the cumulative stride product is `s`.
# If the parent is a vector, then we offset the parent's own indices with parameters of I
compute_offset1(parent::AbstractVector, stride1::Integer, I::Tuple{AbstractRange}) =
    (@_inline_meta; first(I[1]) - stride1*first(axes1(I[1])))
# If the result is one-dimensional and it's a Colon, then linear
# indexing uses the indices along the given dimension.
# If the result is one-dimensional and it's a range, then linear
# indexing might be offset if the index itself is offset
# Otherwise linear indexing always starts with 1.
compute_offset1(parent, stride1::Integer, I::Tuple) =
    (@_inline_meta; compute_offset1(parent, stride1, find_extended_dims(1, I...), find_extended_inds(I...), I))
compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{Slice}, I::Tuple) =
    (@_inline_meta; compute_linindex(parent, I) - stride1*first(axes(parent, dims[1])))  # index-preserving case
compute_offset1(parent, stride1::Integer, dims, inds::Tuple{AbstractRange}, I::Tuple) =
    (@_inline_meta; compute_linindex(parent, I) - stride1*first(axes1(inds[1]))) # potentially index-offsetting case
compute_offset1(parent, stride1::Integer, dims, inds, I::Tuple) =
    (@_inline_meta; compute_linindex(parent, I) - stride1)  # linear indexing starts with 1
function compute_linindex(parent, I::NTuple{N,Any}) where N
    @_inline_meta
    IP = fill_to_length(axes(parent), OneTo(1), Val(N))
    compute_linindex(1, 1, IP, I)
end
function compute_linindex(f, s, IP::Tuple, I::Tuple{ScalarIndex, Vararg{Any}})
    @_inline_meta
    Δi = I[1]-first(IP[1])
    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
end
function compute_linindex(f, s, IP::Tuple, I::Tuple{Any, Vararg{Any}})
    @_inline_meta
    Δi = first(I[1])-first(IP[1])
    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
end
compute_linindex(f, s, IP::Tuple, I::Tuple{}) = f

find_extended_dims(dim, ::ScalarIndex, I...) = (@_inline_meta; find_extended_dims(dim + 1, I...))
find_extended_dims(dim, i1, I...) = (@_inline_meta; (dim, find_extended_dims(dim + 1, I...)...))
find_extended_dims(dim) = ()
find_extended_inds(::ScalarIndex, I...) = (@_inline_meta; find_extended_inds(I...))
find_extended_inds(i1, I...) = (@_inline_meta; (i1, find_extended_inds(I...)...))
find_extended_inds() = ()

function unsafe_convert(::Type{Ptr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{RangeIndex}}}) where {T,N,P}
    return unsafe_convert(Ptr{T}, V.parent) + _memory_offset(V.parent, map(first, V.indices)...)
end

pointer(V::FastSubArray, i::Int) = pointer(V.parent, V.offset1 + V.stride1*i)
pointer(V::FastContiguousSubArray, i::Int) = pointer(V.parent, V.offset1 + i)

function pointer(V::SubArray{<:Any,<:Any,<:Array,<:Tuple{Vararg{RangeIndex}}}, is::AbstractCartesianIndex{N}) where {N}
    index = first_index(V)
    strds = strides(V)
    for d = 1:N
        index += (is[d]-1)*strds[d]
    end
    return pointer(V.parent, index)
end

# indices are taken from the range/vector
# Since bounds-checking is performance-critical and uses
# indices, it's worth optimizing these implementations thoroughly
axes(S::SubArray) = (@_inline_meta; _indices_sub(S.indices...))
_indices_sub(::Real, I...) = (@_inline_meta; _indices_sub(I...))
_indices_sub() = ()
function _indices_sub(i1::AbstractArray, I...)
    @_inline_meta
    (unsafe_indices(i1)..., _indices_sub(I...)...)
end


text:
using PyPlot

function binarySearch(; f, left, right, niter=50, eps=1e-6)
    mid = left
    for i = 1 : niter
        mid = left + (right - left) / 2

        if abs((right - left) / mid) < eps
            return mid
        elseif f(mid) < 0
            left = mid
        elseif f(mid) > 0
            right = mid
        else
            return mid
        end
    end
    return mid
end

function plotFunction(f, x₀)
    fig = PyPlot.figure(figsize=(10, 10))
    x = collect(500:1:1000)
    y = f.(x)
    grid("on")
    vlines(x₀, -.16, .125)
    plot(x, y)
    show()
end

function findPressure(; left, right, niter, m₁, m₂, V, eps)
    ## We want to find "any" zero of this function
    function f(ρ₂)
        ρ₁ = m₁ * ρ₂ / (V * ρ₂ - m₂)
        P = 8.314 * 298 / 0.029 * ρ₁
        ρ₀ = 616.18
        P₀ = 1e5
        return (ρ₂ - ρ₀) / ρ₂ - 0.2105 * log10((35e6 + P) / (35e6 + P₀))
    end

    x₀ = binarySearch(; f, left, right, eps)
    plotFunction(f, x₀)

    ρ₁ = m₁ * x₀ / (V * x₀ - m₂)
    println("ρ₁", ρ₁)
    P = 8.314 * 298 / 0.029 * ρ₁
    return P
end

let
    args = (
            left = 500,
            right = 1000,
            m₁ = 0.02,
            m₂ = 0.1,
            V = 0.01
            eps = 1e-6
           )
    println(findPressure(; args...))
end

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
using NMR
using Base.Test

ϵ = 0.001

# Phenyl boronic acid in CDCl3 (300 MHz)
dataset_path = joinpath(@__DIR__, "data", "PhB(OH)2", "1")

# write your own tests here
@testset "Read Bruker" begin
    s = Spectrum(dataset_path, 1)
    @test s.default_proc == 1
    @test length(s) == length(s[1])  == s["SI"] == 32768
    @test s["TD"] == 18486
    @test typeof(s[:]) == typeof(s[1][:]) <: Array{T,1} where T <: Number
    @test limits(s)[1] ≈ -0.4997 atol=ϵ
    @test limits(s)[2] ≈ 13.4996 atol=ϵ
    @test all( [1.2, 3.5, 7.9, 12.9] .∈ Ref(s) )
    @test all( [14.2, 30.9, -2.0, -0.7] .∉ Ref(s) )
end

@testset "Unit conversions" begin
    s = Spectrum(dataset_path, 1)
    sf = s["SF"]
    bf = s["BF1"]
    @test NMR.sr(sf, bf) ≈ 0.0
    @test -2 < ppmtoindex(s, 6.4502) - 16502 < 2
    let rng = ppmtoindex(s, limits(s))
        @test 1 ∈ rng
        @test s["SI"] - 1 ∈ rng
    end
end


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

@doc raw"""JSONSchemaProps is a JSON-Schema following Specification Draft 4 (http://json-schema.org/).

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps(;
        ref=nothing,
        schema=nothing,
        additionalItems=nothing,
        additionalProperties=nothing,
        allOf=nothing,
        anyOf=nothing,
        default=nothing,
        definitions=nothing,
        dependencies=nothing,
        description=nothing,
        enum=nothing,
        example=nothing,
        exclusiveMaximum=nothing,
        exclusiveMinimum=nothing,
        externalDocs=nothing,
        format=nothing,
        id=nothing,
        items=nothing,
        maxItems=nothing,
        maxLength=nothing,
        maxProperties=nothing,
        maximum=nothing,
        minItems=nothing,
        minLength=nothing,
        minProperties=nothing,
        minimum=nothing,
        multipleOf=nothing,
        not=nothing,
        nullable=nothing,
        oneOf=nothing,
        pattern=nothing,
        patternProperties=nothing,
        properties=nothing,
        required=nothing,
        title=nothing,
        type=nothing,
        uniqueItems=nothing,
        x_kubernetes_embedded_resource=nothing,
        x_kubernetes_int_or_string=nothing,
        x_kubernetes_list_map_keys=nothing,
        x_kubernetes_list_type=nothing,
        x_kubernetes_map_type=nothing,
        x_kubernetes_preserve_unknown_fields=nothing,
    )

    - ref::String
    - schema::String
    - additionalItems::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool
    - additionalProperties::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool
    - allOf::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}
    - anyOf::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}
    - default::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON : default is a default value for undefined object fields. Defaulting is a beta feature under the CustomResourceDefaulting feature gate. CustomResourceDefinitions with defaults must be created using the v1 (or newer) CustomResourceDefinition API.
    - definitions::Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}
    - dependencies::Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrStringArray}
    - description::String
    - enum::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON}
    - example::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON
    - exclusiveMaximum::Bool
    - exclusiveMinimum::Bool
    - externalDocs::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation
    - format::String : format is an OpenAPI v3 format string. Unknown formats are ignored. The following formats are validated:  - bsonobjectid: a bson object ID, i.e. a 24 characters hex string - uri: an URI as parsed by Golang net/url.ParseRequestURI - email: an email address as parsed by Golang net/mail.ParseAddress - hostname: a valid representation for an Internet host name, as defined by RFC 1034, section 3.1 [RFC1034]. - ipv4: an IPv4 IP as parsed by Golang net.ParseIP - ipv6: an IPv6 IP as parsed by Golang net.ParseIP - cidr: a CIDR as parsed by Golang net.ParseCIDR - mac: a MAC address as parsed by Golang net.ParseMAC - uuid: an UUID that allows uppercase defined by the regex (?i)^[0-9a-f]{8}-?[0-9a-f]{4}-?[0-9a-f]{4}-?[0-9a-f]{4}-?[0-9a-f]{12}$ - uuid3: an UUID3 that allows uppercase defined by the regex (?i)^[0-9a-f]{8}-?[0-9a-f]{4}-?3[0-9a-f]{3}-?[0-9a-f]{4}-?[0-9a-f]{12}$ - uuid4: an UUID4 that allows uppercase defined by the regex (?i)^[0-9a-f]{8}-?[0-9a-f]{4}-?4[0-9a-f]{3}-?[89ab][0-9a-f]{3}-?[0-9a-f]{12}$ - uuid5: an UUID5 that allows uppercase defined by the regex (?i)^[0-9a-f]{8}-?[0-9a-f]{4}-?5[0-9a-f]{3}-?[89ab][0-9a-f]{3}-?[0-9a-f]{12}$ - isbn: an ISBN10 or ISBN13 number string like \"0321751043\" or \"978-0321751041\" - isbn10: an ISBN10 number string like \"0321751043\" - isbn13: an ISBN13 number string like \"978-0321751041\" - creditcard: a credit card number defined by the regex ^(?:4[0-9]{12}(?:[0-9]{3})?|5[1-5][0-9]{14}|6(?:011|5[0-9][0-9])[0-9]{12}|3[47][0-9]{13}|3(?:0[0-5]|[68][0-9])[0-9]{11}|(?:2131|1800|35\\d{3})\\d{11})$ with any non digit characters mixed in - ssn: a U.S. social security number following the regex ^\\d{3}[- ]?\\d{2}[- ]?\\d{4}$ - hexcolor: an hexadecimal color code like \"#FFFFFF: following the regex ^#?([0-9a-fA-F]{3}|[0-9a-fA-F]{6})$ - rgbcolor: an RGB color code like rgb like \"rgb(255,255,2559\" - byte: base64 encoded binary data - password: any kind of string - date: a date string like \"2006-01-02\" as defined by full-date in RFC3339 - duration: a duration string like \"22 ns\" as parsed by Golang time.ParseDuration or compatible with Scala duration format - datetime: a date time string like \"2014-12-15T19:30:20.000Z\" as defined by date-time in RFC3339.
    - id::String
    - items::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray
    - maxItems::Int64
    - maxLength::Int64
    - maxProperties::Int64
    - maximum::Float64
    - minItems::Int64
    - minLength::Int64
    - minProperties::Int64
    - minimum::Float64
    - multipleOf::Float64
    - not::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps
    - nullable::Bool
    - oneOf::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}
    - pattern::String
    - patternProperties::Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}
    - properties::Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}
    - required::Vector{String}
    - title::String
    - type::String
    - uniqueItems::Bool
    - x_kubernetes_embedded_resource::Bool : x-kubernetes-embedded-resource defines that the value is an embedded Kubernetes runtime.Object, with TypeMeta and ObjectMeta. The type must be object. It is allowed to further restrict the embedded object. kind, apiVersion and metadata are validated automatically. x-kubernetes-preserve-unknown-fields is allowed to be true, but does not have to be if the object is fully specified (up to kind, apiVersion, metadata).
    - x_kubernetes_int_or_string::Bool : x-kubernetes-int-or-string specifies that this value is either an integer or a string. If this is true, an empty type is allowed and type as child of anyOf is permitted if following one of the following patterns:  1) anyOf:    - type: integer    - type: string 2) allOf:    - anyOf:      - type: integer      - type: string    - ... zero or more
    - x_kubernetes_list_map_keys::Vector{String} : x-kubernetes-list-map-keys annotates an array with the x-kubernetes-list-type `map` by specifying the keys used as the index of the map.  This tag MUST only be used on lists that have the \"x-kubernetes-list-type\" extension set to \"map\". Also, the values specified for this attribute must be a scalar typed field of the child structure (no nesting is supported).
    - x_kubernetes_list_type::String : x-kubernetes-list-type annotates an array to further describe its topology. This extension must only be used on lists and may have 3 possible values:  1) `atomic`: the list is treated as a single entity, like a scalar.      Atomic lists will be entirely replaced when updated. This extension      may be used on any type of list (struct, scalar, ...). 2) `set`:      Sets are lists that must not have multiple items with the same value. Each      value must be a scalar, an object with x-kubernetes-map-type `atomic` or an      array with x-kubernetes-list-type `atomic`. 3) `map`:      These lists are like maps in that their elements have a non-index key      used to identify them. Order is preserved upon merge. The map tag      must only be used on a list with elements of type object. Defaults to atomic for arrays.
    - x_kubernetes_map_type::String : x-kubernetes-map-type annotates an object to further describe its topology. This extension must only be used when type is object and may have 2 possible values:  1) `granular`:      These maps are actual maps (key-value pairs) and each fields are independent      from each other (they can each be manipulated by separate actors). This is      the default behaviour for all maps. 2) `atomic`: the list is treated as a single entity, like a scalar.      Atomic maps will be entirely replaced when updated.
    - x_kubernetes_preserve_unknown_fields::Bool : x-kubernetes-preserve-unknown-fields stops the API server decoding step from pruning fields which are not specified in the validation schema. This affects fields recursively, but switches back to normal pruning behaviour if nested properties or additionalProperties are specified in the schema. This can either be true or undefined. False is forbidden.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps <: SwaggerModel
    ref::Any # spec type: Union{ Nothing, String } # spec name: \$ref
    schema::Any # spec type: Union{ Nothing, String } # spec name: \$schema
    additionalItems::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool } # spec name: additionalItems
    additionalProperties::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool } # spec name: additionalProperties
    allOf::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps} } # spec name: allOf
    anyOf::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps} } # spec name: anyOf
    default::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON } # spec name: default
    definitions::Any # spec type: Union{ Nothing, Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps} } # spec name: definitions
    dependencies::Any # spec type: Union{ Nothing, Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrStringArray} } # spec name: dependencies
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    enum::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON} } # spec name: enum
    example::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON } # spec name: example
    exclusiveMaximum::Any # spec type: Union{ Nothing, Bool } # spec name: exclusiveMaximum
    exclusiveMinimum::Any # spec type: Union{ Nothing, Bool } # spec name: exclusiveMinimum
    externalDocs::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation } # spec name: externalDocs
    format::Any # spec type: Union{ Nothing, String } # spec name: format
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    items::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray } # spec name: items
    maxItems::Any # spec type: Union{ Nothing, Int64 } # spec name: maxItems
    maxLength::Any # spec type: Union{ Nothing, Int64 } # spec name: maxLength
    maxProperties::Any # spec type: Union{ Nothing, Int64 } # spec name: maxProperties
    maximum::Any # spec type: Union{ Nothing, Float64 } # spec name: maximum
    minItems::Any # spec type: Union{ Nothing, Int64 } # spec name: minItems
    minLength::Any # spec type: Union{ Nothing, Int64 } # spec name: minLength
    minProperties::Any # spec type: Union{ Nothing, Int64 } # spec name: minProperties
    minimum::Any # spec type: Union{ Nothing, Float64 } # spec name: minimum
    multipleOf::Any # spec type: Union{ Nothing, Float64 } # spec name: multipleOf
    not::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps } # spec name: not
    nullable::Any # spec type: Union{ Nothing, Bool } # spec name: nullable
    oneOf::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps} } # spec name: oneOf
    pattern::Any # spec type: Union{ Nothing, String } # spec name: pattern
    patternProperties::Any # spec type: Union{ Nothing, Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps} } # spec name: patternProperties
    properties::Any # spec type: Union{ Nothing, Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps} } # spec name: properties
    required::Any # spec type: Union{ Nothing, Vector{String} } # spec name: required
    title::Any # spec type: Union{ Nothing, String } # spec name: title
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    uniqueItems::Any # spec type: Union{ Nothing, Bool } # spec name: uniqueItems
    x_kubernetes_embedded_resource::Any # spec type: Union{ Nothing, Bool } # spec name: x-kubernetes-embedded-resource
    x_kubernetes_int_or_string::Any # spec type: Union{ Nothing, Bool } # spec name: x-kubernetes-int-or-string
    x_kubernetes_list_map_keys::Any # spec type: Union{ Nothing, Vector{String} } # spec name: x-kubernetes-list-map-keys
    x_kubernetes_list_type::Any # spec type: Union{ Nothing, String } # spec name: x-kubernetes-list-type
    x_kubernetes_map_type::Any # spec type: Union{ Nothing, String } # spec name: x-kubernetes-map-type
    x_kubernetes_preserve_unknown_fields::Any # spec type: Union{ Nothing, Bool } # spec name: x-kubernetes-preserve-unknown-fields

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps(;ref=nothing, schema=nothing, additionalItems=nothing, additionalProperties=nothing, allOf=nothing, anyOf=nothing, default=nothing, definitions=nothing, dependencies=nothing, description=nothing, enum=nothing, example=nothing, exclusiveMaximum=nothing, exclusiveMinimum=nothing, externalDocs=nothing, format=nothing, id=nothing, items=nothing, maxItems=nothing, maxLength=nothing, maxProperties=nothing, maximum=nothing, minItems=nothing, minLength=nothing, minProperties=nothing, minimum=nothing, multipleOf=nothing, not=nothing, nullable=nothing, oneOf=nothing, pattern=nothing, patternProperties=nothing, properties=nothing, required=nothing, title=nothing, type=nothing, uniqueItems=nothing, x_kubernetes_embedded_resource=nothing, x_kubernetes_int_or_string=nothing, x_kubernetes_list_map_keys=nothing, x_kubernetes_list_type=nothing, x_kubernetes_map_type=nothing, x_kubernetes_preserve_unknown_fields=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("\$ref"), ref)
        setfield!(o, Symbol("ref"), ref)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("\$schema"), schema)
        setfield!(o, Symbol("schema"), schema)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("additionalItems"), additionalItems)
        setfield!(o, Symbol("additionalItems"), additionalItems)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("additionalProperties"), additionalProperties)
        setfield!(o, Symbol("additionalProperties"), additionalProperties)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("allOf"), allOf)
        setfield!(o, Symbol("allOf"), allOf)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("anyOf"), anyOf)
        setfield!(o, Symbol("anyOf"), anyOf)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("default"), default)
        setfield!(o, Symbol("default"), default)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("definitions"), definitions)
        setfield!(o, Symbol("definitions"), definitions)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("dependencies"), dependencies)
        setfield!(o, Symbol("dependencies"), dependencies)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("enum"), enum)
        setfield!(o, Symbol("enum"), enum)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("example"), example)
        setfield!(o, Symbol("example"), example)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("exclusiveMaximum"), exclusiveMaximum)
        setfield!(o, Symbol("exclusiveMaximum"), exclusiveMaximum)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("exclusiveMinimum"), exclusiveMinimum)
        setfield!(o, Symbol("exclusiveMinimum"), exclusiveMinimum)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("externalDocs"), externalDocs)
        setfield!(o, Symbol("externalDocs"), externalDocs)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("format"), format)
        setfield!(o, Symbol("format"), format)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("maxItems"), maxItems)
        setfield!(o, Symbol("maxItems"), maxItems)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("maxLength"), maxLength)
        setfield!(o, Symbol("maxLength"), maxLength)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("maxProperties"), maxProperties)
        setfield!(o, Symbol("maxProperties"), maxProperties)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("maximum"), maximum)
        setfield!(o, Symbol("maximum"), maximum)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("minItems"), minItems)
        setfield!(o, Symbol("minItems"), minItems)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("minLength"), minLength)
        setfield!(o, Symbol("minLength"), minLength)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("minProperties"), minProperties)
        setfield!(o, Symbol("minProperties"), minProperties)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("minimum"), minimum)
        setfield!(o, Symbol("minimum"), minimum)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("multipleOf"), multipleOf)
        setfield!(o, Symbol("multipleOf"), multipleOf)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("not"), not)
        setfield!(o, Symbol("not"), not)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("nullable"), nullable)
        setfield!(o, Symbol("nullable"), nullable)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("oneOf"), oneOf)
        setfield!(o, Symbol("oneOf"), oneOf)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("pattern"), pattern)
        setfield!(o, Symbol("pattern"), pattern)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("patternProperties"), patternProperties)
        setfield!(o, Symbol("patternProperties"), patternProperties)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("required"), required)
        setfield!(o, Symbol("required"), required)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("title"), title)
        setfield!(o, Symbol("title"), title)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("uniqueItems"), uniqueItems)
        setfield!(o, Symbol("uniqueItems"), uniqueItems)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("x-kubernetes-embedded-resource"), x_kubernetes_embedded_resource)
        setfield!(o, Symbol("x_kubernetes_embedded_resource"), x_kubernetes_embedded_resource)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("x-kubernetes-int-or-string"), x_kubernetes_int_or_string)
        setfield!(o, Symbol("x_kubernetes_int_or_string"), x_kubernetes_int_or_string)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("x-kubernetes-list-map-keys"), x_kubernetes_list_map_keys)
        setfield!(o, Symbol("x_kubernetes_list_map_keys"), x_kubernetes_list_map_keys)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("x-kubernetes-list-type"), x_kubernetes_list_type)
        setfield!(o, Symbol("x_kubernetes_list_type"), x_kubernetes_list_type)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("x-kubernetes-map-type"), x_kubernetes_map_type)
        setfield!(o, Symbol("x_kubernetes_map_type"), x_kubernetes_map_type)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps, Symbol("x-kubernetes-preserve-unknown-fields"), x_kubernetes_preserve_unknown_fields)
        setfield!(o, Symbol("x_kubernetes_preserve_unknown_fields"), x_kubernetes_preserve_unknown_fields)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps = Dict{Symbol,Symbol}(Symbol("\$ref")=>Symbol("ref"), Symbol("\$schema")=>Symbol("schema"), Symbol("additionalItems")=>Symbol("additionalItems"), Symbol("additionalProperties")=>Symbol("additionalProperties"), Symbol("allOf")=>Symbol("allOf"), Symbol("anyOf")=>Symbol("anyOf"), Symbol("default")=>Symbol("default"), Symbol("definitions")=>Symbol("definitions"), Symbol("dependencies")=>Symbol("dependencies"), Symbol("description")=>Symbol("description"), Symbol("enum")=>Symbol("enum"), Symbol("example")=>Symbol("example"), Symbol("exclusiveMaximum")=>Symbol("exclusiveMaximum"), Symbol("exclusiveMinimum")=>Symbol("exclusiveMinimum"), Symbol("externalDocs")=>Symbol("externalDocs"), Symbol("format")=>Symbol("format"), Symbol("id")=>Symbol("id"), Symbol("items")=>Symbol("items"), Symbol("maxItems")=>Symbol("maxItems"), Symbol("maxLength")=>Symbol("maxLength"), Symbol("maxProperties")=>Symbol("maxProperties"), Symbol("maximum")=>Symbol("maximum"), Symbol("minItems")=>Symbol("minItems"), Symbol("minLength")=>Symbol("minLength"), Symbol("minProperties")=>Symbol("minProperties"), Symbol("minimum")=>Symbol("minimum"), Symbol("multipleOf")=>Symbol("multipleOf"), Symbol("not")=>Symbol("not"), Symbol("nullable")=>Symbol("nullable"), Symbol("oneOf")=>Symbol("oneOf"), Symbol("pattern")=>Symbol("pattern"), Symbol("patternProperties")=>Symbol("patternProperties"), Symbol("properties")=>Symbol("properties"), Symbol("required")=>Symbol("required"), Symbol("title")=>Symbol("title"), Symbol("type")=>Symbol("type"), Symbol("uniqueItems")=>Symbol("uniqueItems"), Symbol("x-kubernetes-embedded-resource")=>Symbol("x_kubernetes_embedded_resource"), Symbol("x-kubernetes-int-or-string")=>Symbol("x_kubernetes_int_or_string"), Symbol("x-kubernetes-list-map-keys")=>Symbol("x_kubernetes_list_map_keys"), Symbol("x-kubernetes-list-type")=>Symbol("x_kubernetes_list_type"), Symbol("x-kubernetes-map-type")=>Symbol("x_kubernetes_map_type"), Symbol("x-kubernetes-preserve-unknown-fields")=>Symbol("x_kubernetes_preserve_unknown_fields"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps = Dict{Symbol,String}(Symbol("\$ref")=>"String", Symbol("\$schema")=>"String", Symbol("additionalItems")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool", Symbol("additionalProperties")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool", Symbol("allOf")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}", Symbol("anyOf")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}", Symbol("default")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON", Symbol("definitions")=>"Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}", Symbol("dependencies")=>"Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrStringArray}", Symbol("description")=>"String", Symbol("enum")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON}", Symbol("example")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON", Symbol("exclusiveMaximum")=>"Bool", Symbol("exclusiveMinimum")=>"Bool", Symbol("externalDocs")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation", Symbol("format")=>"String", Symbol("id")=>"String", Symbol("items")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray", Symbol("maxItems")=>"Int64", Symbol("maxLength")=>"Int64", Symbol("maxProperties")=>"Int64", Symbol("maximum")=>"Float64", Symbol("minItems")=>"Int64", Symbol("minLength")=>"Int64", Symbol("minProperties")=>"Int64", Symbol("minimum")=>"Float64", Symbol("multipleOf")=>"Float64", Symbol("not")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps", Symbol("nullable")=>"Bool", Symbol("oneOf")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}", Symbol("pattern")=>"String", Symbol("patternProperties")=>"Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}", Symbol("properties")=>"Dict{String, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps}", Symbol("required")=>"Vector{String}", Symbol("title")=>"String", Symbol("type")=>"String", Symbol("uniqueItems")=>"Bool", Symbol("x-kubernetes-embedded-resource")=>"Bool", Symbol("x-kubernetes-int-or-string")=>"Bool", Symbol("x-kubernetes-list-map-keys")=>"Vector{String}", Symbol("x-kubernetes-list-type")=>"String", Symbol("x-kubernetes-map-type")=>"String", Symbol("x-kubernetes-preserve-unknown-fields")=>"Bool")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps)
    true
end

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


text:
include("../../ApproximateCOSMO.jl/src/COSMO.jl")
using Main.COSMO
include("../src/sedumi.jl"); include("../utils/plot_utils.jl")
using SparseArrays, MAT

A, b, c, K = load_sedumi_file(filename)
x1, value1, status1 = solve_sedumi_jump(A, b, c, K, COSMO.Optimizer)
x2, value2, status2 = solve_sedumi_jump_dual(A, b, c, K, COSMO.Optimizer, max_iter = 500000, check_termination = 1000)


text:
include("gtf_parser.jl")

# open a gtf stream, mode should be either r or w
function gtf_open(filename::AbstractString, mode::AbstractString="r")
    return opengene_open(filename, mode)
end

gtf_close(stream) = close(stream)

"""
Write a Gtf object into a file
Gtf data field is tab-separated
"""
function gtf_write(filename::AbstractString, obj::Gtf)
    stream = gtf_open(filename, "w")
    if stream == false
        return false
    end
    gtf_write_header(stream, obj.header)
    gtf_write_data(stream, obj.data)
    close(stream)
end

"""
Load a GTF file as a DataFrame
"""
function gtf_read(filename::AbstractString; loaddata = true)
    stream = gtf_open(filename)
    if stream == false
        return false
    end
    hd = gtf_read_header(stream)
    hd_rows = length(hd)

    # reopen the file
    # and skip the header to read contents
    stream = gtf_open(filename)
    skipped = 0
    while skipped < hd_rows
        readline(stream)
        skipped += 1
    end

    # if loaddata is set true, then only return the gtf object
    # else also return the stream
    if loaddata
        data = gtf_read_data(stream)
        return return Gtf(hd, data)
    else
        data = GtfData()
        return Gtf(hd, data), stream
    end
end

# read a row from a gtf stream, and just return it
function gtf_read_row(stream)
    if eof(stream)
        return false
    end
    line = readline(stream)
    gtfitem = make_gtf_row(line)
    return gtfitem
end

# read a row from a gtf stream, and write it into the gtf struct
function gtf_read_row!(gtf::Gtf, stream)
    if eof(stream)
        return false
    end
    line = readline(stream)
    gtfitem = make_gtf_row(line)
    if gtfitem != false
        push!(gtf.data, gtfitem)
    end
    return gtfitem
end

function gtf_read_data(stream)
    data = GtfData()
    while true
        if eof(stream)
            return data
        end
        line = readline(stream)
        gtfitem = make_gtf_row(line)
        if gtfitem != false
            push!(data, gtfitem)
        end
    end
    return data
end

function make_gtf_row(line)
    line = rstrip(line, '\n')
    items = split(line, "\t")
    if length(items) < 9
        return false
    end
    seqname = items[1]
    source = items[2]
    feature = items[3]
    start_pos = parse(Int64,items[4])
    end_pos = parse(Int64,items[5])
    score = items[6]
    strand = items[7]
    frame = items[8]
    attributes = gtf_parse_attributes(ASCIIString(items[9]))
    gtfitem = GtfItem(seqname, source, feature, start_pos, end_pos, score, strand, frame, attributes)
    return gtfitem
end

function gtf_read_header(stream)
    header = GtfHeader()
    while true
        if eof(stream)
            return header
        end
        line = readline(stream)
        if startswith(line, "##")
            line = lstrip(rstrip(line, '\n'), '#')
            items = split(line, ":")
            if length(items) < 2
                continue
            end
            key = strip(items[1], ' ')
            value = strip(items[2], ' ')
            header[key] = value
        else
            break
        end
    end
    return header
end

# write a gtf header to a file
function gtf_write_header(stream, header::GtfHeader)
    try
        for (k, v) in header
            print(stream, "##", k, ": ", v, "\n")
        end
    catch e
        warn(e)
        return false
    end
end

# write gtf data field into stream
function gtf_write_data(stream, data::GtfData)
    for item in data
        print(stream, item.seqname, "\t")
        print(stream, item.source, "\t")
        print(stream, item.feature, "\t")
        print(stream, item.start_pos, "\t")
        print(stream, item.end_pos, "\t")
        print(stream, item.score, "\t")
        print(stream, item.strand, "\t")
        print(stream, item.frame, "\t")
        gtf_print_attributes(item.attributes, stream)
        print(stream, "\n")
    end
end

function gtf_print_attributes(attributes, stream)
    i = 0
    for (key, value) in attributes
        vars = split(value, "/")
        for v in 1:length(vars)
            var = vars[v]
            print(stream, "$key $var;")
            if v < length(vars)
                print(stream, "  ")
            end
        end
        i += 1
        if i < length(attributes)
            print(stream, "  ")
        end
    end
end


text:
# Model
include_model("particle")

# Dimensions
nq = 3                    # configuration dimension
nu = 3                    # control dimension
nc = 1                    # number of contact points
nf = 2                    # number of parameters for friction cone
nb = nc * nf
ns = 1

# Parameters
μ = 0.5      # coefficient of friction
mass = 1.0   # mass
g = 9.81

n = 2 * nq
m = nu + nc + nb + 3 + ns
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:3)
idx_η = (1:0)#nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + 3 .+ (1:ns)

# Methods
function M_func(model::Particle, q)
    @SMatrix [model.mass 0.0 0.0;
              0.0 model.mass 0.0;
              0.0 0.0 model.mass]
end

G_func(model::Particle, q) = @SVector [0.0, 0.0, model.mass * model.g]

function ϕ_func(::Particle, q)
    @SVector[q[3]]
end

B_func(::Particle, q) = @SMatrix [1.0 0.0 0.0;
                                  0.0 1.0 0.0;
                                  0.0 0.0 1.0]

N_func(::Particle, q) = @SMatrix [0.0 0.0 1.0]

function P_func(model::Particle, q)
   return @SMatrix [1.0 0.0 0.0;
                    0.0 1.0 0.0]
end

# function P_func(model::Particle, q)
#    return @SMatrix [1.0 0.0 0.0;
#                     0.0 1.0 0.0;
#                     -1.0 0.0 0.0;
#                     0.0 -1.0 0.0]
# end

function friction_cone(model::Particle, u)
	λ = u[model.idx_λ]
	b = u[model.idx_b]
	return @SVector [model.μ * λ[1] - norm(b' * b)]
end

function maximum_dissipation(model::Particle, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]

	λ = u[model.idx_λ]
	b = u[model.idx_b]
	ψ = u[model.idx_ψ]

	return [P_func(model, q3) * (q3 - q2) / h - ψ[2:3];
		model.μ * λ[1] * ψ[2:3] + ψ[1] * b]
end

function fd(model::Particle{Discrete, FixedTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

    [q2⁺ - q2⁻;
    ((1.0 / h) * (M_func(model, q1) * (SVector{3}(q2⁺) - SVector{3}(q1))
    - M_func(model, q2⁺) * (SVector{3}(q3) - SVector{3}(q2⁺)))
    + h * (transpose(B_func(model, q3)) * SVector{3}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{1}(λ)
    + transpose(P_func(model, q3)) * SVector{2}(b)
    - G_func(model, q2⁺)))]
end

model = Particle{Discrete, FixedTime}(n, m, d,
				 mass, g, μ,
				 nq,
				 nu,
				 nc,
				 nf,
				 nb,
				 ns,
				 idx_u,
				 idx_λ,
				 idx_b,
				 idx_ψ,
				 idx_η,
				 idx_s)

# Horizon
T = 11

# Time step
tf = 1.0
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 0.0
_ul = zeros(model.m)
_ul[model.idx_b] .= -Inf
_ul[model.idx_ψ] .= -Inf
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.0, 1.0]
v1 = [1.0, 1.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q2]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj = PenaltyObjective(1.0e3, model.m)

# Constraints
include_constraints("contact_soc")
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
x0 = [0.01 * rand(model.n) for t = 1:T] #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-5, c_tol = 1.0e-5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model,
    state_to_configuration(x̄),
    Δt = h)

open(vis)


text:
@testset "Robots: point foot quadruped" begin
    # TODO: add tests
    q0 = nominal_configuration(point_foot_quadruped)
    # visualizer
    vis = RoboDojo.Visualizer();
    @test visualize!(vis, box, [point_foot_quadruped], Δt=0.1);
end


text:
#######################################
#
# BDF file
#
#######################################

"""
Import Biosemi files
"""
function import_biosemi(fname::Union{AbstractString, IO}; kwargs...)

    @info("Importing BIOSEMI data file")

    # Read raw data using BDF.jl
    data, triggers, trigger_channel, system_code_channel = readBDF(identity(fname), transposeData=true; kwargs...)
    header = readBDFHeader(identity(fname))

    # Check the sample rate
    sample_rate = header["sampRate"]
    if sum(diff(sample_rate)) != 0
        @warn("Sampling rate varies across channels")
        sample_rate = NaN
    else
        sample_rate = sample_rate[1]
    end

    reference_channel = "Raw"

    # Tidy the trigger channel to standard names
    triggers = Dict("Code"     => triggers["code"],
                "Index"    => triggers["idx"],
                "Duration" => triggers["dur"])

    # Tidy channel names if required
    if any(header["chanLabels"] .== "B16")  # Cz is usually present
        @debug("  Converting names from BIOSEMI to 10-20")
        header["chanLabels"] = channelNames_biosemi_1020(header["chanLabels"])
    end

    system_codes = create_events(system_code_channel, sample_rate)

    return  data, triggers, system_codes, sample_rate, reference_channel, header
end

# Create events from channel
############################

function create_events(channel::Array{Int16,1}, fs::Number)

    startPoints = vcat(1, findall(diff(channel) .!= 0).+1)
    stopPoints = vcat(findall(diff(channel) .!= 0), length(channel))
    trigDurs = (stopPoints - startPoints)/fs

    evt = channel[startPoints]
    evtTab = Dict("Code" => evt,
                  "Index" => startPoints,
                  "Duration" => trigDurs)

end

# Create channel from events
############################

function create_channel(t::Dict, data::Array, fs::Number; kwargs...)

    create_channel(t, maximum(size(data)), fs; kwargs...)
end

function create_channel(t::Dict, l::Int, fs::Number; code::AbstractString="Code", index::AbstractString="Index",
                        duration::AbstractString="Duration")

    @debug("Creating trigger channel from data. Length: $l Triggers: $(length(t[index])) Fs: $fs")

    channel = Array{Int16}(l)

    # Initialise array to 252 code
    for i = 1:l ; channel[i] = 252; end

    for i = 1:length(t[index])-1
        channel[t[index][i] : t[index][i] + round.(Int, t[duration][i] * fs)] = t[code][i]
    end

    return channel
end

# Change biosemi labels to 1020
#######################################

function channelNames_biosemi_1020(original::S) where S <: AbstractString

    biosemi_1020 = ["A01" "Fp1"
                    "A1"  "Fp1"
                    "A05" "F3"
                    "A5"  "F3"
                    "A09" "FC5"
                    "A9"  "FC5"
                    "A13" "C3"
                    "A17" "CP5"
                    "A21" "P3"
                    "A25" "PO7"
                    "A29" "Oz"
                    "B01" "Fpz"
                    "B1"  "Fpz"
                    "B05" "AFz"
                    "B5"  "AFz"
                    "B9"  "F6"
                    "B13" "FC4"
                    "B17" "C2"
                    "B21" "TP8"
                    "B25" "P2"
                    "B29" "P10"
                    "A02" "AF7"
                    "A2"  "AF7"
                    "A06" "F5"
                    "A6"  "F5"
                    "A10" "FC3"
                    "A14" "C5"
                    "A18" "CP3"
                    "A22" "P5"
                    "A26" "PO3"
                    "A30" "POz"
                    "B02" "Fp2"
                    "B2"  "Fp2"
                    "B06" "Fz"
                    "B6"  "Fz"
                    "B10" "F8"
                    "B14" "FC2"
                    "B18" "C4"
                    "B22" "CP6"
                    "B26" "P4"
                    "B30" "PO8"
                    "A03" "AF3"
                    "A3"  "AF3"
                    "A07" "F7"
                    "A7"  "F7"
                    "A11" "FC1"
                    "A15" "T7"
                    "A19" "CP1"
                    "A23" "P7"
                    "A27" "O1"
                    "A31" "Pz"
                    "B03" "AF8"
                    "B3"  "AF8"
                    "B07" "F2"
                    "B7"  "F2"
                    "B11" "FT8"
                    "B15" "FCz"
                    "B19" "C6"
                    "B23" "CP4"
                    "B27" "P6"
                    "B31" "PO4"
                    "A04" "F1"
                    "A4"  "F1"
                    "A08" "FT7"
                    "A8"  "FT7"
                    "A12" "C1"
                    "A16" "TP7"
                    "A20" "P1"
                    "A24" "P9"
                    "A28" "Iz"
                    "A32" "CPz"
                    "B04" "AF4"
                    "B4"  "AF4"
                    "B08" "F4"
                    "B8"  "F4"
                    "B12" "FC6"
                    "B16" "Cz"
                    "B20" "T8"
                    "B24" "CP2"
                    "B28" "P8"
                    "B32" "O2"
                    "Status" "Status"]

    idx = something(findfirst(isequal(original), biosemi_1020), 0)

    if idx == 0
        error("Channel $original is unknown")
    end

    converted = biosemi_1020[:, 2][idx]
end

function channelNames_biosemi_1020(original::Array{S}) where S <: AbstractString

    converted = Array{AbstractString}(size(original))

    @info("Fixing channel names of $(length(original)) channels")

    for i = 1:length(original)
        converted[i] = channelNames_biosemi_1020(original[i])
    end

    return converted
end


text:
export RectilinearGrid
export xcoords, ycoords, zcoords, xmin, ymin, zmin, xmax, ymax, zmax

struct RectilinearGrid{X,Y,Z,T}
    x::Vec{X,T}
    y::Vec{Y,T}
    z::Vec{Z,T}

    function RectilinearGrid{X,Y,Z,T}(x, y, z) where {X,Y,Z,T}
        if X < 2
            error("Must have at least 2 X-divisions.")
        end
        if Y < 2
            error("Must have at least 2 Y-divisions.")
        end
        if Z !== 0
            if Z < 2
                error("Must have at least 2 Z-divisions.")
            end
            if any(i->z[i+1] - z[i] < 0, 1:Z-1)
                error("Divisions must be monotonically increasing.")
            end
        end
        if any(i->x[i+1] - x[i] < 0, 1:X-1) ||
           any(i->y[i+1] - y[i] < 0, 1:Y-1)
            error("Divisions must be monotonically increasing.")
        end
        return new{X,Y,Z,T}(x, y, z)
    end
end

# constructors
function RectilinearGrid(x::Vec{X,T}, y::Vec{Y,T}) where {X,Y,T}
    return RectilinearGrid{X,Y,0,T}(x, y, Vec{0,T}())
end

function RectilinearGrid(x::Vec{X,T}, y::Vec{Y,T}, z::Vec{Z,T}) where {X,Y,Z,T}
    return RectilinearGrid{X,Y,Z,T}(x, y, z)
end

function RectilinearGrid(x, y)
    X = length(x)
    Y = length(y)
    sx = Vec{X}(x)
    sy = Vec{Y}(y)
    return RectilinearGrid(sx, sy)
end

function RectilinearGrid(x, y, z)
    X = length(x)
    Y = length(y)
    Z = length(z)
    sx = Vec{X}(x)
    sy = Vec{Y}(y)
    sz = Vec{Z}(z)
    return RectilinearGrid(sx, sy, sz)
end

Base.issubset(g1::RectilinearGrid, g2::RectilinearGrid) = g1.x ⊆ g2.x &&
                                                          g1.y ⊆ g2.y &&
                                                          g1.z ⊆ g2.z
xcoords(rg::RectilinearGrid) = rg.x
ycoords(rg::RectilinearGrid) = rg.y
zcoords(rg::RectilinearGrid) = rg.z
xmin(rg::RectilinearGrid) = rg.x[begin]
ymin(rg::RectilinearGrid) = rg.y[begin]
zmin(rg::RectilinearGrid) = rg.z[begin]
xmax(rg::RectilinearGrid) = rg.x[end]
ymax(rg::RectilinearGrid) = rg.y[end]
zmax(rg::RectilinearGrid) = rg.z[end]

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Pkg
Pkg.activate(".")
using Revise
using WeightedNetwork

# generat random (undirected) network:
g1 = WeightedNetwork.rand_network(6, 0.2)

# it's type ist SimpleWeightedGraph.
g1.weights

# plot graph:
plot_network(g1)

# generate random (directed) network:
g2 = WeightedNetwork.rand_directed_network(6, 0.2)

# plot graph:
plot_network(g2)

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
using MakieGallery, AbstractPlotting, GLMakie, Test

empty!(MakieGallery.plotting_backends)
push!(MakieGallery.plotting_backends, "Makie")
database = MakieGallery.load_database([
                    "tutorials.jl",
                    "attributes.jl",
                    "short_tests.jl"
                    ])

test_record_path = joinpath(@__DIR__, "test_recordings")
recordings = MakieGallery.record_examples(test_record_path)
@assert length(recordings) == length(database)
@info "Precompile script has completed execution."


text:
function write(file::AbstractString, jws::JSONWorksheet; kwargs...)
    open(file, "w") do io
        write(io, jws; kwargs...)
    end
end
function write(io::IO, jws::JSONWorksheet; pretty = true, drop_null = false)
    if pretty
        JSON3.pretty(io, Tables.rows(jws))
    else
        JSON3.write(io, Tables.rows(jws))
    end
    # drop null array such as [null, null, ....]
    if drop_null
        replace!(io, r"(\"[\w]*\":null,)|(,?\"[\w]*\":null)" => "")
    end
    return io
end

function write(path::String, jwb::JSONWorkbook; kwargs...)
    f = splitext(basename(xlsxpath(jwb)))[1]
    for s in sheetnames(jwb)
        write(joinpath(path, "$(f)_$(s).json"), jwb[s]; kwargs...)
    end
end

function write_xlsx(file::String, jwb::JSONWorkbook; delim = ";", anchor_cell = "A1")
    XLSX.openxlsx(file, mode="w") do xf

        for (i, s) in enumerate(sheetnames(jwb))
            jws = jwb[i]
            if i == 1
                sheet = xf[1]
                XLSX.rename!(sheet, s)
            else
                sheet = XLSX.addsheet!(xf, s)
            end

            colnames = pointer_to_colname.(jws.pointer)
            columns = []
            for p in jws.pointer
                data = get.(jws.data, Ref(p), missing)
                if eltype(data) <: Array
                    data = join.(data, delim)
                end
                push!(columns, data)
            end

            XLSX.writetable!(sheet, columns, colnames, anchor_cell=XLSX.CellRef(anchor_cell))
        end
    end
end


text:
# In M, a n×n distance matrix, find s such that:
# ∀i ∈ 1:n, ∃j ∈ 1:n, i≠j; M[i, j] ≤ s.
function sampledensity(M::AbstractMatrix{T}) where T
    s = typemin(T)
    for i in 1:size(M, 1)
        nearest = typemax(T)
        for j in 1:size(M, 1)
            i == j && continue
            dist = M[j, i]
            nearest = dist < nearest ? dist : nearest
        end
        s = s > nearest ? s : nearest
    end
    s
end

sampledensity(g::AbstractGraph) = maximum(weights(g))

function edgetoorder(G::AbstractGraph)
    w = adjacency_matrix(G)
    for (i, e) in enumerate(edges(G))
        w[src(e), dst(e)] = i
        w[dst(e), src(e)] = i
    end
    Symmetric(w)
end
function edgetoorder(G::AbstractSimpleWeightedGraph)
    es = collect(edges(G))
    sort!(es, by = weight)
    I = Vector{Int}(length(es))
    J = similar(I)
    i = 1
    for e in es
        I[i] = src(e)
        J[i] = dst(e)
        i += 1
    end
    n = nv(G)
    Symmetric(sparse(I, J, 1:length(es), n, n), :U)
end

function spannergraph(pts::AbstractMatrix{T},
                      metric = Euclidean(), r = nothing
                      ) ::SimpleWeightedGraph{Int, T} where T
    s    = r == nothing ? 2 * sampledensity(pairwise(metric, pts)) : T(r)
    n, m = size(pts)
    _pts = unique(reinterpret(SVector{n, T}, pts, (m,)))
    spannergraph(_pts, metric, s)
end
function spannergraph(pts::AbstractVector{SVector{D, T}},
                      metric = Euclidean(), r = nothing
                      ) ::SimpleWeightedGraph{Int, T} where {D, T}
    n = length(pts)
    s = r == nothing ?
        2 * T(sampledensity(pairwise(metric, reinterpret(T, pts, (D, n))))) :
        T(r)
    btree  = BallTree(pts, metric, reorder = true)
    neighs = inrange(btree, pts, s, true)
    is = Int[]; js = Int[]; ds = T[]
    for (i, ns) in enumerate(neighs)
        for j in ns
            i < j || continue
            append!(is, [i, j])
            append!(js, [j, i])
            d = evaluate(metric, pts[i], pts[j])
            append!(ds, [d, d])
        end
    end
    n = length(pts)
    SimpleWeightedGraph(sparse(is, js, ds, n, n))
end

function persistence(pts::AbstractVecOrMat;
                     show_trace = false,
                     use_hash   = false,
                     skip_small = true,
                     resolution = nothing,
                     metric     = Euclidean())

    graph = spannergraph(pts, metric, resolution)
    res   = persistence(graph,
                        eltype(weights(graph)),
                        show_trace = show_trace,
                        skip_small = skip_small,
                        use_hash   = use_hash)
    add_data(res, pts)
end

function persistence(graph::AbstractGraph,
                     T::DataType = eltype(weights(graph));
                     skip_small ::Bool = true,
                     show_trace ::Bool = false,
                     use_hash   ::Bool = false) # <- T mora bit eltype(weights(G))
                                                        #    preveri če mora bit parameter

    is_directed(graph) && throw(ArgumentError("Directed graphs not supported!"))

    show_trace && print("# ===================== #\n",
                        "| Intrinsic persistence |\n",
                        "# ===================== #\n\n",
                        "Preprocessing...")

    s = sampledensity(graph)
    n = nv(graph)

    # Constants.
    order     = edgetoorder(graph)
    dijkstra  = parallel_multisource_dijkstra_shortest_paths(graph)
    distmx    = dijkstra.dists
    parents   = dijkstra.parents
    nrows     = ne(graph)
    triangles = equilaterals(distmx, tol = s)

    # Result.
    generators = IntrinsicGenerator{T}[]

    # State.
    reductionmx      = fill(Int[], nrows)
    st               = start(triangles)
    ngenerators_left = nrows - n + length(connected_components(graph))

    # Debug variables.
    r = zero(T)
    currindex = 0
    nskipped  = 0
    nreduced  = 0

    # Buffers.
    column = Int[]; sizehint!(column, n)
    buff   = Int[]; sizehint!(buff, n)
    paths  = Int[], Int[], Int[]

    if use_hash
        hashes = Set{UInt}()
    end

    if show_trace
        @time ntriangles = iterlength(triangles)
        print("\r")
        println("* Sample density:       $(round(s, 5))")
        println("* Spanner density:      $(density(graph))")
        println("* Connected components: $(length(connected_components(graph)))")
        println("* Number of triangles:  $(ntriangles)")
        println("* Number of generators: $(ngenerators_left)")
        println()
    end

    maxcollen = 0
    while ngenerators_left > 0 && !done(triangles, st)
        currindex += 1
        Δ, st = next(triangles, st)
        i, j, k = Δ.vertices
        r = Δ.r

        if show_trace && currindex % (ntriangles ÷ 100) == 0
            showprogress(currindex, r, ntriangles, nreduced,
                         length(generators), ngenerators_left, nrows)
        end

        empty!(column); empty!(paths[1]); empty!(paths[2]); empty!(paths[3])
        collectpath!(column, paths[1], i, j, parents, order)
        collectpath!(column, paths[2], j, k, parents, order)
        collectpath!(column, paths[3], k, i, parents, order)
        sort!(column)

        # Skipping.
        #r < 2s && length(column) > 3           && continue
        use_hash && !add_hash!(hashes, column) && continue
        !alluniquesorted(column)               && continue
        !islocallyminimal(paths, distmx, s)    && continue

        # Matrix reduction.
        low = getlow(column)
        stop = false
        while !stop && low ≠ 0 && !isempty(reductionmx[low])
            column, buff = sumcols!(buff, column, reductionmx[low])
            low = getlow(column)
            maxcollen = max(length(column), maxcollen)

            if use_hash && !add_hash!(hashes, column)
                stop = true
            end
        end
        stop && continue

        # Output.
        if low ≠ 0
            reductionmx[low] = column
            column = Int[]

            if (!skip_small && sum(length.(paths)) > 3) || r > s
                path = Int[]
                for i in 1:3
                    append!(path, paths[i])
                end
                lengths = similar(path, T)
                for p_i in eachindex(path)
                    p_j = p_i < length(path) ? p_i + 1 : 1
                    lengths[p_i] = distmx[path[p_i], path[p_j]]
                end

                push!(generators, IntrinsicGenerator(Δ, path, lengths))
            end
            ngenerators_left -= 1
            nskipped += 1
        else
            nreduced += 1
        end
    end

    if show_trace
        showprogress(currindex, r, ntriangles, nreduced,
                     length(generators), ngenerators_left, nrows)
        if use_hash
            println("no. hashes: $(length(hashes))")
        end
        println()
        println()
        @show maxcollen
    end
    ngenerators_left > 0 && warn("Not all generators found!")

    IntrinsicPersistenceResult(generators, graph)
end

function collectpath!(∂_j, side, v1, v2, parents, order)
    # sem not gre lahko čekiranje če je prejšnji trikotnik vsebovan v poti
    index = v1
    while parents[v2, index] ≠ 0
        push!(side, index)
        nxtindex = parents[v2, index]
        push!(∂_j, order[index, nxtindex])
        index = nxtindex
    end
end

function add_hash!(hashes, column)
    h = hash(column)
    if h in hashes
        false
    else
        push!(hashes, h)
        true
    end
end

function alluniquesorted(vec)
    for i in 1:length(vec)-1
        vec[i] == vec[i+1] && return false
    end
    true
end

# Check for each pair of points on opposing sides of the triangle whether the
# path through triangle is the shortest path.
function islocallyminimal(sides, distmx, s)
    @inbounds for (i,j,k) in ((1,2,3), (2,3,1), (3,1,2))
        a = sides[i]; b = sides[j]; c = sides[k]
        common = b[1]
        opp1   = a[1]
        opp2   = c[1]

        for i in 2:length(a), j in 2:length(b)
            direct = distmx[a[i], b[j]]
            around = min(distmx[a[i], common] + distmx[common, b[j]],
                         distmx[a[i], opp1] + distmx[opp1, opp2] + distmx[opp2, b[j]])
            if direct + s < around
                return false
            end
        end
    end
    true
end

function getlow(col::AbstractVector{T}) where T
    isempty(col) ? zero(T) : last(col)
end

@inline function sumcols!(buff, col1, col2)
    i = j = 1
    n = length(col1)
    m = length(col2)
    @inbounds while i ≤ n && j ≤ m
        if col1[i] < col2[j]
            push!(buff, col1[i])
            i += 1
        elseif col1[i] > col2[j]
            push!(buff, col2[j])
            j += 1
        else
            i += 1
            j += 1
        end
    end
    @inbounds while j ≤ m
        push!(buff, col2[j])
        j += 1
    end
    @inbounds while i ≤ n
        push!(buff, col1[i])
        i += 1
    end
    empty!(col1)
    buff, col1
end

function iterlength(itr)
    len = 0
    for _ in itr
        len += 1
    end
    len
end

function showprogress(j, r, ntriangles, nreduced, ngenerators, ngenerators_left, nrows)
    print("\rTriangles visited:",
          lpad(round(Int, 100 * j/ntriangles), 3), "%, ℓ ≈ ",
          rpad(round(3r, 5), 7, 0),
          ", reduced:",
          lpad(round(Int, 100 * nreduced/j), 3), "%",
          ". Generators found:",
          lpad(ngenerators, 2),
          ", left: ",
          rpad(ngenerators_left, floor(Int, log10(nrows) + 1)))
    flush(STDOUT::IO)
end


text:
"""
    struct GadgetCylinder{T} <: AbstractGadgetGeometry
        pos_start::Vector{T}
        pos_end::Vector{T}
        radius::T
    end

Defines a cylinder by two end points `pos_start` and `pos_start` as well as a `radius`.
To be used for [`read_particles_in_geometry`](@ref)
"""
struct GadgetCylinder{T} <: AbstractGadgetGeometry
    pos_start::Vector{T}
    pos_end::Vector{T}
    radius::T
end

"""
    function get_geometry_center(cylinder::GadgetCylinder)

Returns the center of `cylinder`.
"""
get_geometry_center(cylinder::GadgetCylinder) = 1 // 2 .* cylinder.pos_start .+ cylinder.pos_end

"""
    get_geometry_box_corners(cylinder::GadgetCylinder)

Returns a tuple with the lower left and upper right corner of a box which contains the `cylinder`.
"""
function get_geometry_box_corners(cylinder::GadgetCylinder)
    min.(cylinder.pos_start, cylinder.pos_end) .- cylinder.radius,
    max.(cylinder.pos_start, cylinder.pos_end) .+ cylinder.radius
end

"""
    get_geometry_mask(cylinder::GadgetCylinder, pos::Matrix{T}) where T

Returns the indices of all particles contained in the `cylinder`.
"""
function get_geometry_mask(cylinder::GadgetCylinder, pos::Matrix{T}) where T

    # allocate mask array
    mask = Vector{Bool}(undef, size(pos,2))

    # check every position array if it's contained in the cylinder
    # and update the mask array
    @inbounds for i = 1:size(pos,2)
        mask[i] = check_in_cylinder(pos[:,i],
                                    cylinder.pos_start,
                                    cylinder.pos_end,
                                    cylinder.radius)
    end

    return findall(mask)
end


text:
type MergedIter
    it1   # iterator 1
    it2   # iterator 2
    lt    # custom `<` operator to compare elements of 2 iterators
    MergeIter(it1, it2) = new(it1, it2, lt)
end

type MergedState
    hd1   # current head of 1st iterator
    hd2   # current head of 2nd iterator
    s1    # current state of 1st iterator
    s2    # current state of 2nd iterator
end

function Base.start(it::MergedIter)
    s1 = start(it.it1)
    s2 = start(it.it2)
    # doing 1 step ahead and caching next element of each iter
    hd1, s1 = safenext(it.it1, s1)
    hd2, s2 = safenext(it.it2, s2)
    return MergedState(hd1, hd2, s1, s2)
end

"""
Like `next`, but may be called on iterator that is already done. In this case
new item is `nothing` and new state is the same as previous.
"""
function safenext(it, s)
    if done(it, s)
        return nothing, s
    else
        return next(it, s)
    end
end

function Base.next(it::MergedIter, s::MergedState)
    lt = it.lt
    if done(it, s) throw(ArgumentError("Iterator is done")) end
    it1 = it.it1; it2 = it.it2; hd1 = s.hd1; hd2 = s.hd2; s1 = s.s1; s2 = s.s2
    while hd1 == hd2 && hd2 != nothing
        # move second iterator till heads are different or iterator is done
        hd2, s2 = safenext(it2, s2)
    end
    if complete(it1, s1, hd1)
        # 1st iterator is done, but we have one save value in hd1
        if hd1 == nothing || lt(hd2, hd1) # hd1 > hd2
            next_hd2, next_s2 = safenext(it2, s2)
            return hd2, MergedState(hd1, next_hd2, s1, next_s2)
        else  # hd1 cached and hd1 < hd2
            return hd1, MergedState(nothing, hd2, s1, s2)
        end
    elseif complete(it2, s2, hd2)
        # 2st iterator is done, but we have one save value in hd2
        if hd2 == nothing || lt(hd1, hd2) # hd2 > hd1
            next_hd1, next_s1 = safenext(it1, s1)
            return hd1, MergedState(next_hd1, hd2, next_s1, s2)
        else  # hd2 cached and hd2 < hd1
            return hd2, MergedState(hd1, nothing, s1, s2)
        end
    elseif lt(hd1, hd2) # hd1 < hd2
        next_hd1, next_s1 = safenext(it1, s1)
        return hd1, MergedState(next_hd1, hd2, next_s1, s2)
    else
        next_hd2, next_s2 = safenext(it2, s2)
        return hd2, MergedState(hd1, next_hd2, s1, next_s2)
    end
end

"""
Check if iterator is done and no head doesn't contain a cached value.
"""
function complete(it, s, hd)
    return done(it, s) && hd == nothing
end

function Base.done(it::MergedIter, s::MergedState)
    return (done(it.it1, s.s1) && done(it.it2, s.s2) &&
            s.hd1 == nothing && s.hd2 == nothing)
end

"""
Merge 2 sorted iterators. Let `hd1` be a head of 1st iterator and `hd2` -
a head of 2nd iterator. Then merged iterator works as follows:

 * if hd1 < hd2, hd1 is emitted
 * if hd1 > hd2, hd2 is emitted
 * if hd1 == hd2, hd1 is emitted and hd2 is discarded (so no duplicates
                  is produced)
 * if one of iterators is done, tail of another is taken

Comparison is done using `lt` option that defaults to `isless` function.
"""
function mergesorted(it1, it2; lt = isless)
    return MergedIter(it1, it2, lt)
end


text:
export SuffixArray

struct SuffixArray
    data::Vector
    index::Vector{Int}
end

function SuffixArray(str::String)
    data = map(Int32, collect(str))
    k = maximum(data)
    index = sais(data, Int(k)+1)
    SuffixArray(data, index)
end

Base.length(sa::SuffixArray) = length(sa.index)

function prefixsearch(sa::SuffixArray, query::Vector{Int32})
    ii = 1
    ij = length(sa)
    while true
        ik = (ii+ij) ÷ 2
        xi = sa.index[ik]

        qi = 1
        while sa.data[xi+qi-1] == query[qi]
            if qi > length(query)
                i = ik
                while getlcp() == lcp
                    i--
                end
                i = ik
                while true
                    i++
                end

                for i = ik:-1:1
                    getlcp(sa, i)
                end
            end
            qi == length(query) && break
            xi+qi-1 == length(sa) && break
            qi += 1
        end
        lcp = qi - 1
        ii == ij && return xi:xi+lcp-1

        if sa.data[xi+qi-1] < query[qi]
            ii = ik + 1
        else
            ij = ik - 1
        end
    end
end

function getlcp(sa::SuffixArray, xi::Int, query::Vector)
    qi = 1
    while sa.data[xi+qi-1] == query[qi]
        qi == length(query) && return qi
        xi+qi-1 == length(sa) && return qi
        qi += 1
    end
    qi - 1
end

"""
    lcparray

Kasai's algorithm for linear-time construction of LCP array from Suffix Array
"""
function lcparray(sa::Vector{Int}, data::Vector{T}) where T<:Integer
    n = length(sa)
    lcps = similar(sa)
    rank = similar(sa)
    for i = 1:n
        rank[sa[i]] = i
    end

    lcp = 0
    for i = 1:n
        rank[i] > 1 || continue
        j = sa[rank[i]-1]
        while i+lcp <= n && j+lcp <= n && data[i+lcp] == data[j+lcp]
            lcp += 1
        end
        lcps[rank[i]] = lcp
        lcp > 0 && (lcp -= 1)
    end
    lcps[1] = 0
    lcps
end


text:
module Tests

using LRUCache
using Base.Test

function test_order(lru, keys, vals)
    if length(lru) != 0
        f = first(lru.q)
        @test f.k == keys[1]
        @test f.v == vals[1]
        n = f.next
        i = 2
        while !is(n, f)
            @test n.k == keys[i]
            @test n.v == vals[i]
            i += 1
            n = n.next
        end
    end
end

const CACHE = LRU{Int, Int}(20)
# Test insertion ordering
kvs = 1:10
for i in reverse(kvs)
    CACHE[i] = i
end
test_order(CACHE, 1:10, 1:10)

# Test reinsertion ordering
kvs = [1, 3, 5, 7, 9, 2, 4, 6, 8, 10]
for i in reverse(kvs)
    CACHE[i] = i
end
test_order(CACHE, kvs, kvs)

# Test least recently used items are evicted
resize!(CACHE, 5)
test_order(CACHE, kvs[1:5], kvs[1:5])

resize!(CACHE, 10)
test_order(CACHE, kvs[1:5], kvs[1:5])

kvs = 1:11
for i in reverse(kvs)
    CACHE[i] = i
end
test_order(CACHE, kvs[1:end-1], kvs[1:end-1])

# Test lookups, and that lookups reorder
kvs = 1:10
for i in reverse(kvs)
    CACHE[i] = 2*i
end
for i in kvs
    @test CACHE[i] == 2*i
end
test_order(CACHE, reverse(kvs), reverse(2*kvs))

# Test empty!
empty!(CACHE)
@test length(CACHE) == 0

# Test get! with default parameter
@test get!(CACHE, 1, 2)  == 2    # not in cache
@test get!(CACHE, 1, 2)  == 2    # in cache

# Test get! with function call
val = get!(CACHE, 2) do
          3
      end
@test val == 3
get!(CACHE, 2) do
    error("this shouldn't have been called!")
end

# Test @get! with begin block
val = @get! CACHE 3 begin
    4
end
@test val == 4
@get! CACHE 3 begin
    error("this shouldn't have been called!")
end

# Test Abstract typed cache. All we're checking for here is that the container
# is able to hold abstract types without issue. Insertion order is already
# tested above.
const CACHE2 = LRU{String, Integer}(5)
CACHE2["test"] = 4
CACHE2[utf8("test2")] = BigInt(5)
@test CACHE2["test"] == 4
@test CACHE2["test2"] == 5

end


text:
# Autogenerated wrapper script for CSPICE_jll for i686-w64-mingw32
export brief, chronos, ckbrief, commnt, dskbrief, dskexp, frmdiff, inspekt, libcspice, libcsupport, mkdsk, mkspk, msopck, spacit, spkdiff, spkmerge, tobin, toxfr, version

JLLWrappers.@generate_wrapper_header("CSPICE")
JLLWrappers.@declare_executable_product(brief)
JLLWrappers.@declare_executable_product(chronos)
JLLWrappers.@declare_executable_product(ckbrief)
JLLWrappers.@declare_executable_product(commnt)
JLLWrappers.@declare_executable_product(dskbrief)
JLLWrappers.@declare_executable_product(dskexp)
JLLWrappers.@declare_executable_product(frmdiff)
JLLWrappers.@declare_executable_product(inspekt)
JLLWrappers.@declare_library_product(libcspice, "libcspice.dll")
JLLWrappers.@declare_library_product(libcsupport, "libcsupport.dll")
JLLWrappers.@declare_executable_product(mkdsk)
JLLWrappers.@declare_executable_product(mkspk)
JLLWrappers.@declare_executable_product(msopck)
JLLWrappers.@declare_executable_product(spacit)
JLLWrappers.@declare_executable_product(spkdiff)
JLLWrappers.@declare_executable_product(spkmerge)
JLLWrappers.@declare_executable_product(tobin)
JLLWrappers.@declare_executable_product(toxfr)
JLLWrappers.@declare_executable_product(version)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        brief,
        "bin\\brief.exe",
    )

    JLLWrappers.@init_executable_product(
        chronos,
        "bin\\chronos.exe",
    )

    JLLWrappers.@init_executable_product(
        ckbrief,
        "bin\\ckbrief.exe",
    )

    JLLWrappers.@init_executable_product(
        commnt,
        "bin\\commnt.exe",
    )

    JLLWrappers.@init_executable_product(
        dskbrief,
        "bin\\dskbrief.exe",
    )

    JLLWrappers.@init_executable_product(
        dskexp,
        "bin\\dskexp.exe",
    )

    JLLWrappers.@init_executable_product(
        frmdiff,
        "bin\\frmdiff.exe",
    )

    JLLWrappers.@init_executable_product(
        inspekt,
        "bin\\inspekt.exe",
    )

    JLLWrappers.@init_library_product(
        libcspice,
        "bin\\libcspice.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsupport,
        "bin\\libcsupport.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mkdsk,
        "bin\\mkdsk.exe",
    )

    JLLWrappers.@init_executable_product(
        mkspk,
        "bin\\mkspk.exe",
    )

    JLLWrappers.@init_executable_product(
        msopck,
        "bin\\msopck.exe",
    )

    JLLWrappers.@init_executable_product(
        spacit,
        "bin\\spacit.exe",
    )

    JLLWrappers.@init_executable_product(
        spkdiff,
        "bin\\spkdiff.exe",
    )

    JLLWrappers.@init_executable_product(
        spkmerge,
        "bin\\spkmerge.exe",
    )

    JLLWrappers.@init_executable_product(
        tobin,
        "bin\\tobin.exe",
    )

    JLLWrappers.@init_executable_product(
        toxfr,
        "bin\\toxfr.exe",
    )

    JLLWrappers.@init_executable_product(
        version,
        "bin\\version.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module PhylogeneticTrees

    import JuMP, Gurobi, MathProgBase
    import LightGraphs
    import CSV, IterTools, DataFrames
    import Printf

    include("tree.jl")
    include("popdata.jl")
    include("model.jl")
    include("model-nodebased.jl")
    include("model-common.jl")
    include("symmetries.jl")
    include("display.jl")
end