text:
""" can_flatten(::Type{T}) """ # TODO document new flatten
@inline flatten(A, inds) = _flatten(can_flatten(inds), A, inds)
@inline _flatten(::True, A, inds) = flatten_indices(A, inds)
_flatten(::False, A, inds) = inds

can_flatten(x) = can_flatten(typeof(x))
can_flatten(::Type{T}) where {T} = static(false)
can_flatten(::Type{T}) where {T<:AbstractArray{<:AbstractCartesianIndex}} = static(false)
can_flatten(::Type{T}) where {T<:CartesianIndices} = static(true)
can_flatten(::Type{T}) where {T<:AbstractArray{Bool}} = static(ndims(T) > 1)
can_flatten(::Type{T}) where {T<:AbstractCartesianIndex} = static(true)

can_flatten(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}} = _can_flatten(T, static(N))
_can_flatten(::Type{T}, ::StaticInt{0}) where {T} = static(false)
@inline function _can_flatten(::Type{T}, n::StaticInt{N}) where {T,N}
    return can_flatten(_get_tuple(T, n)) | _can_flatten(T, n - static(1))
end

"""
    flatten_indices(A, inds) -> flatten_indices(axes(A), inds)

Flatten multi-dimension spanning argument into separate arguments for each dimension.
"""
flatten_indices(A, i::Tuple) = flatten_indices(lazy_axes(A), i)
@inline function flatten_indices(a::Tuple, i::Tuple{I,Vararg{Any}}) where {I}
    return (first(i), flatten_indices(tail(a), tail(i))...)
end
@inline function flatten_indices(a::Tuple, i::Tuple{I,Vararg{Any}}) where {N,I<:AbstractCartesianIndex{N}}
    _, atail = Base.IteratorsMD.split(a, Val(N))
    return (Tuple(first(i))..., flatten_indices(atail, tail(i))...)
end
@inline function flatten_indices(a::Tuple, i::Tuple{I,Vararg{Any}}) where {N,I<:CartesianIndices{N}}
    _, atail = Base.IteratorsMD.split(a, Val(N))
    return (axes(first(i))..., flatten_indices(atail, tail(i))...)
end
# we preserve CartesianIndices{0} for dropping dimensions
@inline function flatten_indices(a::Tuple, i::Tuple{I,Vararg{Any}}) where {I<:CartesianIndices{0}}
    return (first(i), flatten_indices(tail(a), tail(i))...)
end
flatten_indices(::Tuple, ::Tuple{}) = ()

## is_canonical ##
is_canonical(x) = is_canonical(typeof(x))
is_canonical(::Type{T}) where {T} = static(false)
is_canonical(::Type{Int}) = static(true)
is_canonical(::Type{<:StaticInt}) = static(true)
is_canonical(::Type{StepRange{Int,Int}}) = static(true)
is_canonical(::Type{UnitRange{Int}}) = static(true)
is_canonical(::Type{OneTo{Int}}) = static(true)
is_canonical(::Type{<:OptionallyStaticRange}) = static(true)
is_canonical(::Type{Vector{Int}}) = static(true)
is_canonical(::Type{<:Slice}) = static(true)
@inline is_canonical(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}} = _is_canon(T, static(N))
_is_canon(::Type{T}, ::StaticInt{0}) where {T} = static(true)
@inline function _is_canon(::Type{T}, n::StaticInt{N}) where {T,N}
    return is_canonical(_get_tuple(T, n)) & _is_canon(T, n - static(1))
end

"""
    canonicalize(x)

Checks if `x` is in a canonical form for indexing. If `x` is already in a canonical form
then it is returned unchanged. If `x` is not in a canonical form then it is passed to
`canonical_convert`.
"""
canonicalize(x) = _canonicalize(is_canonical(x), x)
_canonicalize(::True, x) = x
_canonicalize(::False, x) = canonical_convert(x)

canonical_convert(x::Integer) = convert(Int, x)
function canonical_convert(x::AbstractRange)
    return OptionallyStaticStepRange(static_first(x), static_step(x), static_last(x))
end
function canonicalize_convert(x::AbstractUnitRange{<:Integer})
    return OptionallyStaticUnitRange(static_first(x), static_last(x))
end

_layout(::IndexLinear, x::Tuple) = LinearIndices(x)
_layout(::IndexCartesian, x::Tuple) = CartesianIndices(x)

"""
    ArrayStyle(::Type{A})

Used to customize the meaning of indexing arguments in the context of a given array `A`.

See also: [`argdims`](@ref), [`UnsafeIndex`](@ref)
"""
abstract type ArrayStyle end

struct DefaultArrayStyle <: ArrayStyle end

ArrayStyle(A) = ArrayStyle(typeof(A))
ArrayStyle(::Type{A}) where {A} = DefaultArrayStyle()

"""
    argdims(::ArrayStyle, ::Type{T})

What is the dimensionality of the indexing argument of type `T`?
"""
argdims(x, arg) = argdims(x, typeof(arg))
argdims(x, ::Type{T}) where {T} = argdims(ArrayStyle(x), T)
argdims(s::ArrayStyle, arg) = argdims(s, typeof(arg))
# single elements initially map to 1 dimension but that dimension is subsequently dropped.
argdims(::ArrayStyle, ::Type{T}) where {T} = static(0)
argdims(::ArrayStyle, ::Type{T}) where {T<:Colon} = static(1)
argdims(::ArrayStyle, ::Type{T}) where {T<:AbstractArray} = static(ndims(T))
argdims(::ArrayStyle, ::Type{T}) where {N,T<:AbstractCartesianIndex{N}} = static(N)
argdims(::ArrayStyle, ::Type{T}) where {N,T<:AbstractArray{<:AbstractCartesianIndex{N}}} = static(N)
argdims(::ArrayStyle, ::Type{T}) where {N,T<:AbstractArray{<:Any,N}} = static(N)
argdims(::ArrayStyle, ::Type{T}) where {N,T<:LogicalIndex{<:Any,<:AbstractArray{Bool,N}}} = static(N)
_argdims(s::ArrayStyle, ::Type{I}, i::StaticInt) where {I} = argdims(s, _get_tuple(I, i))
function argdims(s::ArrayStyle, ::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}}
    return eachop(_argdims, nstatic(Val(N)), s, T)
end

_is_element_index(i) = _is_element_index(typeof(i))
_is_element_index(::Type{T}) where {T} = static(false)
_is_element_index(::Type{T}) where {T<:AbstractCartesianIndex} = static(true)
_is_element_index(::Type{T}) where {T<:Integer} = static(true)
__is_element_index(::Type{T}, i::StaticInt) where {T} = _is_element_index(_get_tuple(T, i))
function _is_element_index(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}}
    return static(all(eachop(__is_element_index, nstatic(Val(N)), T)))
end
# empty tuples refer to the single element of 0-dimensional arrays
_is_element_index(::Type{Tuple{}}) = static(true)

# are the indexing arguments provided a linear collection into a multidim collection
is_linear_indexing(A, args::Tuple{Arg}) where {Arg} = argdims(A, Arg) < 2
is_linear_indexing(A, args::Tuple{Arg,Vararg{Any}}) where {Arg} = false

"""
    to_indices(A, args::Tuple) -> to_indices(A, axes(A), args)
    to_indices(A, axes::Tuple, args::Tuple)

Maps arguments `args` to the axes of `A`. This is done by iteratively passing each
axis and argument to [`to_index`](@ref). Unique behavior based on the type of `A` may be
accomplished by overloading `to_indices(A, args)`. Unique axis-argument behavior can
be accomplished using `to_index(axis, arg)`.
"""
@propagate_inbounds to_indices(A, args::Tuple) = _to_indices(A, flatten(A, args))
@propagate_inbounds function to_indices(A, args::Tuple{LinearIndices})
    return to_indices(A, lazy_axes(A), axes(first(args)))
end
@propagate_inbounds _to_indices(A, args::Tuple) = __to_indices(is_canonical(args), A, args)
@propagate_inbounds function __to_indices(::False, A, args)
    if is_linear_indexing(A, args)
        return (to_index(eachindex(IndexLinear(), A), first(args)),)
    else
        return to_indices(A, lazy_axes(A), args)
    end
end
# don't waste time reconstructing tuple if already converted
@propagate_inbounds function __to_indices(::True, A, args)
    if is_linear_indexing(A, args)
        @boundscheck if !checkindex(Bool, eachindex(IndexLinear(), A), first(args))
            throw(BoundsError(A, args))
        end
        return args
    else
        @boundscheck if !Base.checkbounds_indices(Bool, lazy_axes(A), args)
            throw(BoundsError(A, args))
        end
        return args
    end
end

@propagate_inbounds to_indices(A, args::Tuple{}) = to_indices(A, lazy_axes(A), ())
@propagate_inbounds function to_indices(A, axs::Tuple, args::Tuple{I,Vararg{Any}},) where {I}
    return _to_indices(argdims(A, I), A, axs, args)
end
@propagate_inbounds function _to_indices(::StaticInt{0}, A, axs::Tuple, args::Tuple)
    return (to_index(first(axs), first(args)), to_indices(A, _maybe_tail(axs), _maybe_tail(args))...)
end
@propagate_inbounds function _to_indices(::StaticInt{1}, A, axs::Tuple, args::Tuple)
    return (to_index(first(axs), first(args)), to_indices(A, _maybe_tail(axs), _maybe_tail(args))...)
end
@propagate_inbounds function _to_indices(::StaticInt{N}, A, axs::Tuple, args::Tuple) where {N}
    axes_front, axes_tail = Base.IteratorsMD.split(axs, Val(N))
    return (to_index(_layout(IndexStyle(A), axes_front), first(args)), to_indices(A, axes_tail, _maybe_tail(args))...)
end
@propagate_inbounds function to_indices(A, axs::Tuple, args::Tuple{})
    @boundscheck if length(first(axs)) != 1
        error("Cannot drop dimension of size $(length(first(axs))).")
    end
    return to_indices(A, tail(axs), args)
end
@propagate_inbounds function to_indices(A, ::Tuple{}, args::Tuple{I,Vararg{Any}}) where {I}
    return (to_index(OneTo(1), first(args)), to_indices(A, (), tail(args))...)
end
to_indices(A, axs::Tuple{}, args::Tuple{}) = ()
_maybe_tail(::Tuple{}) = ()
_maybe_tail(x::Tuple) = tail(x)

"""
    to_index([::IndexStyle, ]axis, arg) -> index

Convert the argument `arg` that was originally passed to `getindex` for the dimension
corresponding to `axis` into a form for native indexing (`Int`, Vector{Int}, etc.). New
axis types with unique behavior should use an `IndexStyle` trait:

```julia
to_index(axis::MyAxisType, arg) = to_index(IndexStyle(axis), axis, arg)
to_index(::MyIndexStyle, axis, arg) = ...
```
"""
@propagate_inbounds to_index(axis, arg) = to_index(IndexStyle(axis), axis, arg)
function to_index(s, axis, arg)
    throw(ArgumentError("invalid index: IndexStyle $s does not support indices of " *
                        "type $(typeof(arg)) for instances of type $(typeof(axis))."))
end
to_index(::IndexLinear, axis, arg::Colon) = indices(axis) # Colons get converted to slices by `indices`
to_index(::IndexLinear, axis, arg::CartesianIndices{0}) = arg
to_index(::IndexLinear, axis, arg::CartesianIndices{1}) = axes(arg, 1)
@propagate_inbounds function to_index(::IndexLinear, axis, arg::AbstractCartesianIndex{1})
    return to_index(axis, first(Tuple(arg)))
end
function to_index(::IndexLinear, x, arg::AbstractCartesianIndex{N}) where {N}
    inds = Tuple(arg)
    o = offsets(x)
    s = size(x)
    return first(inds) + (offset1(x) - first(o)) + _subs2int(first(s), tail(s), tail(o), tail(inds))
end
@inline function _subs2int(stride, s::Tuple{Any,Vararg}, o::Tuple{Any,Vararg}, inds::Tuple{Any,Vararg})
    i = ((first(inds) - first(o)) * stride)
    return i + _subs2int(stride * first(s), tail(s), tail(o), tail(inds))
end
function _subs2int(stride, s::Tuple{Any}, o::Tuple{Any}, inds::Tuple{Any})
    return (first(inds) - first(o)) * stride
end
# trailing inbounds can only be 1 or 1:1
_subs2int(stride, ::Tuple{}, ::Tuple{}, ::Tuple{Any}) = static(0)

@propagate_inbounds function to_index(::IndexLinear, x, arg::Union{Array{Bool}, BitArray})
    @boundscheck checkbounds(x, arg)
    return LogicalIndex{Int}(arg)
end
@propagate_inbounds function to_index(::IndexLinear, x, arg::AbstractArray{<:AbstractCartesianIndex})
    @boundscheck _multi_check_index(axes(x), arg) || throw(BoundsError(x, arg))
    return arg
end
@propagate_inbounds function to_index(::IndexLinear, x, arg::LogicalIndex)
    @boundscheck checkbounds(Bool, x, arg) || throw(BoundsError(x, arg))
    return arg
end
@propagate_inbounds function to_index(::IndexLinear, x, arg::Integer)
    @boundscheck checkindex(Bool, indices(x), arg) || throw(BoundsError(x, arg))
    return canonicalize(arg)
end
@propagate_inbounds function to_index(::IndexLinear, axis, arg::AbstractArray{Bool})
    @boundscheck checkbounds(axis, arg)
    return LogicalIndex(arg)
end
@propagate_inbounds function to_index(::IndexLinear, x, arg::AbstractArray{<:Integer})
    @boundscheck checkindex(Bool, x, arg) || throw(BoundsError(x, arg))
    return arg
end
@propagate_inbounds function to_index(::IndexLinear, x, arg::AbstractRange{Integer})
    @boundscheck checkindex(Bool, indices(axis), arg) || throw(BoundsError(axis, arg))
    return static_first(arg):static_step(arg):static_last(arg)
end

## IndexCartesian ##
to_index(::IndexCartesian, x, arg::Colon) = CartesianIndices(x)
to_index(::IndexCartesian, x, arg::CartesianIndices{0}) = arg
to_index(::IndexCartesian, x, arg::AbstractCartesianIndex) = arg
function to_index(::IndexCartesian, x, arg)
    @boundscheck _multi_check_index(axes(x), arg) || throw(BoundsError(x, arg))
    return arg
end
@propagate_inbounds function to_index(::IndexCartesian, x, arg::AbstractArray{<:AbstractCartesianIndex})
    @boundscheck _multi_check_index(axes(x), arg) || throw(BoundsError(x, arg))
    return arg
end
@propagate_inbounds function to_index(::IndexCartesian, x, arg::AbstractArray{Bool})
    @boundscheck checkbounds(x, arg)
    return LogicalIndex(arg)
end

function _multi_check_index(axs::Tuple, arg::AbstractArray{T}) where {T<:AbstractCartesianIndex}
    b = true
    for i in arg
        b &= Base.checkbounds_indices(Bool, axs, (i,))
    end
    return b
end

@propagate_inbounds function to_index(::IndexCartesian, x, arg::Union{Array{Bool}, BitArray})
    @boundscheck checkbounds(x, arg)
    return LogicalIndex{Int}(arg)
end
to_index(::IndexCartesian, x, i::Integer) = NDIndex(_int2subs(offsets(x), size(x), i - offset1(x)))
@inline function _int2subs(o::Tuple{Any,Vararg{Any}}, s::Tuple{Any,Vararg{Any}}, i)
    len = first(s)
    inext = div(i, len)
    return (canonicalize(i - len * inext + first(o)), _int2subs(tail(o), tail(s), inext)...)
end
_int2subs(o::Tuple{Any}, s::Tuple{Any}, i) = canonicalize(i + first(o))

"""
    unsafe_reconstruct(A, data; kwargs...)

Reconstruct `A` given the values in `data`. New methods using `unsafe_reconstruct`
should only dispatch on `A`.
"""
function unsafe_reconstruct(axis::OneTo, data; kwargs...)
    if axis === data
        return axis
    else
        return OneTo(data)
    end
end
function unsafe_reconstruct(axis::UnitRange, data; kwargs...)
    if axis === data
        return axis
    else
        return UnitRange(first(data), last(data))
    end
end
function unsafe_reconstruct(axis::OptionallyStaticUnitRange, data; kwargs...)
    if axis === data
        return axis
    else
        return OptionallyStaticUnitRange(static_first(data), static_last(data))
    end
end
function unsafe_reconstruct(A::AbstractUnitRange, data; kwargs...)
    return static_first(data):static_last(data)
end

"""
    to_axes(A, inds)

Construct new axes given the corresponding `inds` constructed after
`to_indices(A, args) -> inds`. This method iterates through each pair of axes and
indices calling [`to_axis`](@ref).
"""
@inline function to_axes(A, inds::Tuple)
    if ndims(A) === 1
        return (to_axis(axes(A, 1), first(inds)),)
    elseif is_linear_indexing(A, inds)
        return (to_axis(eachindex(IndexLinear(), A), first(inds)),)
    else
        return to_axes(A, axes(A), inds)
    end
end
to_axes(A, a::Tuple, i::Tuple{I,Vararg{Any}},) where {I} = _to_axes(argdims(A, I), A, a, i)
# drop this dimension
_to_axes(::StaticInt{0}, A, axs::Tuple, inds::Tuple) = to_axes(A, tail(inds), tail(inds))
function _to_axes(::StaticInt{1}, A, axs::Tuple, inds::Tuple)
    return (to_axis(first(axs), first(inds)), to_axes(A, tail(axs), tail(inds))...)
end
@propagate_inbounds function _to_axes(::StaticInt{N}, A, axs::Tuple, inds::Tuple) where {N}
    axes_front, axes_tail = Base.IteratorsMD.split(axs, Val(N))
    return (
        to_axis(_layout(IndexStyle(A), axes_front), first(inds)),
        to_axes(A, axes_tail, tail(inds))...
    )
end
to_axes(A, ::Tuple{Ax,Vararg{Any}}, ::Tuple{}) where {Ax} = ()
to_axes(A, ::Tuple{}, ::Tuple{}) = ()

"""
    to_axis(old_axis, index) -> new_axis

Construct an `new_axis` for a newly constructed array that corresponds to the
previously executed `to_index(old_axis, arg) -> index`. `to_axis` assumes that
`index` has already been confirmed to be in bounds. The underlying indices of
`new_axis` begins at one and extends the length of `index` (i.e., one-based indexing).
"""
@inline function to_axis(axis, inds)
    if !can_change_size(axis) &&
       (known_length(inds) !== nothing && known_length(axis) === known_length(inds))
        return axis
    else
        return to_axis(IndexStyle(axis), axis, inds)
    end
end

# don't need to check size b/c slice means it's the entire axis
@inline function to_axis(axis, inds::Slice)
    if can_change_size(axis)
        return copy(axis)
    else
        return axis
    end
end
to_axis(S::IndexLinear, axis, inds) = StaticInt(1):static_length(inds)

################
### getindex ###
################
"""
    ArrayInterface.getindex(A, args...)

Retrieve the value(s) stored at the given key or index within a collection. Creating
another instance of `ArrayInterface.getindex` should only be done by overloading `A`.
Changing indexing based on a given argument from `args` should be done through,
[`to_index`](@ref), or [`to_axis`](@ref).
"""
@propagate_inbounds getindex(A, args...) = unsafe_get_index(A, to_indices(A, args))
@propagate_inbounds function getindex(A; kwargs...)
    return unsafe_get_index(A, to_indices(A, order_named_inds(dimnames(A), kwargs.data)))
end
@propagate_inbounds getindex(x::Tuple, i::Int) = getfield(x, i)
@propagate_inbounds getindex(x::Tuple, ::StaticInt{i}) where {i} = getfield(x, i)

## unsafe_get_index ##
unsafe_get_index(A, inds::Tuple) = _unsafe_get_index(_is_element_index(inds), A, inds)
_unsafe_get_index(::False, A, inds::Tuple) = unsafe_get_collection(A, inds)
_unsafe_get_index(::True, A, inds::Tuple) = __unsafe_get_index(A, inds)
__unsafe_get_index(A, inds::Tuple{}) = unsafe_get_element(A, ())
__unsafe_get_index(A, inds::Tuple{Any}) = unsafe_get_element(A, first(inds))
__unsafe_get_index(A, inds::Tuple{Any,Vararg{Any}}) = unsafe_get_element(A, NDIndex(inds))

"""
    unsafe_get_element(A::AbstractArray{T}, inds::Tuple) -> T

Returns an element of `A` at the indices `inds`. This method assumes all `inds`
have been checked for being in bounds. Any new array type using `ArrayInterface.getindex`
must define `unsafe_get_element(::NewArrayType, inds)`.
"""
unsafe_get_element(a::A, inds) where {A} = _unsafe_get_element(has_parent(A), a, inds)
_unsafe_get_element(::True, a, inds) = unsafe_get_element(parent(a), inds)
_unsafe_get_element(::False, a, inds) = @inbounds(parent(a)[inds])
_unsafe_get_element(::False, a::AbstractArray2, i) = unsafe_get_element_error(a, i)

## Array ##
unsafe_get_element(A::Array, ::Tuple{}) = Base.arrayref(false, A, 1)
unsafe_get_element(A::Array, i::Integer) = Base.arrayref(false, A, Int(i))
unsafe_get_element(A::Array, i::NDIndex) = unsafe_get_element(A, to_index(A, i))

## LinearIndices ##
unsafe_get_element(A::LinearIndices, i::Integer) = Int(i)
unsafe_get_element(A::LinearIndices, i::NDIndex) = unsafe_get_element(A, to_index(A, i))

unsafe_get_element(A::CartesianIndices, i::NDIndex) = CartesianIndex(i)
unsafe_get_element(A::CartesianIndices, i::Integer) = unsafe_get_element(A, to_index(A, i))

unsafe_get_element(A::ReshapedArray, i::Integer) = unsafe_get_element(parent(A), i)
function unsafe_get_element(A::ReshapedArray, i::NDIndex)
    return unsafe_get_element(parent(A), to_index(IndexLinear(), A, i))
end

unsafe_get_element(A::SubArray, i) = @inbounds(A[i])
function unsafe_get_element_error(@nospecialize(A), @nospecialize(i))
    throw(MethodError(unsafe_get_element, (A, i)))
end

# This is based on Base._unsafe_getindex from https://github.com/JuliaLang/julia/blob/c5ede45829bf8eb09f2145bfd6f089459d77b2b1/base/multidimensional.jl#L755.
"""
    unsafe_get_collection(A, inds)

Returns a collection of `A` given `inds`. `inds` is assumed to have been bounds-checked.
"""
function unsafe_get_collection(A, inds)
    axs = to_axes(A, inds)
    dest = similar(A, axs)
    if map(Base.unsafe_length, axes(dest)) == map(Base.unsafe_length, axs)
        _unsafe_get_index!(dest, A, inds...) # usually a generated function, don't allow it to impact inference result
    else
        Base.throw_checksize_error(dest, axs)
    end
    return dest
end

function _generate_unsafe_get_index!_body(N::Int)
    quote
        Base.@_inline_meta
        D = eachindex(dest)
        Dy = iterate(D)
        @inbounds Base.Cartesian.@nloops $N j d -> I[d] begin
            # This condition is never hit, but at the moment
            # the optimizer is not clever enough to split the union without it
            Dy === nothing && return dest
            (idx, state) = Dy
            dest[idx] = unsafe_get_element(src, NDIndex(Base.Cartesian.@ntuple($N, j)))
            Dy = iterate(D, state)
        end
        return dest
    end
end
@generated function _unsafe_get_index!(dest, src, I::Vararg{Any,N}) where {N}
    return _generate_unsafe_get_index!_body(N)
end

_ints2range(x::Integer) = x:x
_ints2range(x::AbstractRange) = x
@inline function unsafe_get_collection(A::CartesianIndices{N}, inds) where {N}
    if (length(inds) === 1 && N > 1) || stride_preserving_index(typeof(inds)) === False()
        return Base._getindex(IndexStyle(A), A, inds...)
    else
        return CartesianIndices(to_axes(A, _ints2range.(inds)))
    end
end
@inline function unsafe_get_collection(A::LinearIndices{N}, inds) where {N}
    if is_linear_indexing(A, inds)
        return @inbounds(eachindex(A)[first(inds)])
    elseif stride_preserving_index(typeof(inds)) === True()
        return LinearIndices(to_axes(A, _ints2range.(inds)))
    else
        return Base._getindex(IndexStyle(A), A, inds...)
    end
end

#################
### setindex! ###
#################
"""
    ArrayInterface.setindex!(A, args...)

Store the given values at the given key or index within a collection.
"""
@propagate_inbounds function setindex!(A, val, args...)
    if can_setindex(A)
        return unsafe_set_index!(A, val, to_indices(A, args))
    else
        error("Instance of type $(typeof(A)) are not mutable and cannot change elements after construction.")
    end
end
@propagate_inbounds function setindex!(A, val; kwargs...)
    return unsafe_set_index!(A, val, to_indices(A, order_named_inds(dimnames(A), kwargs.data)))
end

unsafe_set_index!(A, v, inds::Tuple) = _unsafe_set_index!(_is_element_index(inds), A, v, inds)
_unsafe_set_index!(::False, A, v, inds::Tuple) = unsafe_set_collection!(A, v, inds)
_unsafe_set_index!(::True, A, v, inds::Tuple) = __unsafe_set_index!(A, v, inds)
__unsafe_set_index!(A, v, inds::Tuple{}) = unsafe_set_element!(A, v, ())
function __unsafe_set_index!(A, v, inds::Tuple{Any})
    return unsafe_set_element!(A, v, to_index(A, first(inds)))
end
function __unsafe_set_index!(A, v, inds::Tuple{Any,Vararg{Any}})
    return unsafe_set_element!(A, v, to_index(A, NDIndex(inds)))
end

"""
    unsafe_set_element!(A, val, inds::Tuple)

Sets an element of `A` to `val` at indices `inds`. This method assumes all `inds`
have been checked for being in bounds. Any new array type using `ArrayInterface.setindex!`
must define `unsafe_set_element!(::NewArrayType, val, inds)`.
"""
unsafe_set_element!(a, val, inds) = _unsafe_set_element!(has_parent(a), a, val, inds)
_unsafe_set_element!(::True, a, val, inds) = unsafe_set_element!(parent(a), val, inds)
_unsafe_set_element!(::False, a, val, inds) = @inbounds(parent(a)[inds] = val)

function _unsafe_set_element!(::False, a::AbstractArray2, val, inds)
    unsafe_set_element_error(a, val, inds)
end
unsafe_set_element_error(A, v, i) = throw(MethodError(unsafe_set_element!, (A, v, i)))

function unsafe_set_element!(A::Array{T}, val, ::Tuple{}) where {T}
    Base.arrayset(false, A, convert(T, val)::T, 1)
end
function unsafe_set_element!(A::Array{T}, val, i::Integer) where {T}
    return Base.arrayset(false, A, convert(T, val)::T, Int(i))
end

# This is based on Base._unsafe_setindex!.
"""
    unsafe_set_collection!(A, val, inds)

Sets `inds` of `A` to `val`. `inds` is assumed to have been bounds-checked.
"""
@inline unsafe_set_collection!(A, v, i) = _unsafe_setindex!(A, v, i...)

function _generate_unsafe_setindex!_body(N::Int)
    quote
        x′ = Base.unalias(A, x)
        Base.Cartesian.@nexprs $N d -> (I_d = Base.unalias(A, I[d]))
        idxlens = Base.Cartesian.@ncall $N Base.index_lengths I
        Base.Cartesian.@ncall $N Base.setindex_shape_check x′ (d -> idxlens[d])
        Xy = iterate(x′)
        @inbounds Base.Cartesian.@nloops $N i d->I_d begin
            # This is never reached, but serves as an assumption for
            # the optimizer that it does not need to emit error paths
            Xy === nothing && break
            (val, state) = Xy
            unsafe_set_element!(A, val, NDIndex(Base.Cartesian.@ntuple($N, i)))
            Xy = iterate(x′, state)
        end
        A
    end
end

@generated function _unsafe_setindex!(A, x, I::Vararg{Any,N}) where {N}
    return _generate_unsafe_setindex!_body(N)
end


text:
using DSGE
using JLD2

# Initialize a new instance of an An Schorfheide model
m = AnSchorfheide()

# Data can be constructed automatically, i.e. calling estimate(m),
# but for ease of example we will load in data
@load "../../test/reference/smc.jld2" data

# SMC can take advantage of parallel computing resources, should
# they be available
m <= Setting(:use_parallel_workers, false)

# SMC related settings
m <= Setting(:n_particles, 800)
m <= Setting(:n_Φ, 50) # number of stages in SMC
m <= Setting(:λ, 2.0) # the curvature of the tempering schedule
m <= Setting(:adaptive_tempering_target_smc, false) # Setting a fixed tempering schedule

# No need to find the mode for SMC
# If set to true, the mode will be calculated during estimation
m <= Setting(:reoptimize, false)
m <= Setting(:sampling_method, :SMC)

println("estimation beginning")
# If verbose=:low instead, the mean and sd of the model's
# parameters will not be displayed at each stage of SMC
estimate(m, data, verbose = :low)
println("estimation finished")
println("sample saved in $(rawpath(m, "estimate", "smcsave.h5")))")
println("particle cloud saved in $(rawpath(m, "estimate", "smc_cloud.jld2")))")

# Compute and save a table of moments
moment_tables(m)
println("tex table of moments saved in $(tablespath(m, "estimate"))")

ext: jl


text:
using Base: Float64
using OrdinaryDiffEq
using Trixi
using Plots

# Vgl. elixir_euler_free_stream_curved
###############################################################################
CFL = 0.5          # 2
tspan = (0.0, 2.0)
N = 3
c = 8
mu = 0.0

# Source terms and initial conditions for EOC tests
function initial_condition_convergence_test_viscous(x, t, equations::Trixi.AbstractEquations)
    # μ = equations.mu
    # Pr = equations.Pr
    # κ = equations.gamma
    # k = 2 # k = c_p * μ / Pr
    # k = 114 * μ / Pr
    k = 0
    c = 4
    ω = 0.5
    ini = sin(k * (x[1] + x[2]) - ω * t) + c

    rho = ini
    rho_v1 = ini
    rho_v2 = ini
    rho_e = ini^2

    return SVector(rho, rho_v1, rho_v2, rho_e)
end

@inline function source_terms_convergence_test_viscous(u, x, t, equations::Trixi.AbstractEquations)
    # Same settings as in `initial_condition`
    μ = equations.mu
    Pr = equations.Pr
    κ = equations.gamma
    # k = 2 # k = c_p * μ / Pr
    k = 114 * μ / Pr
    # k = 0
    c = 4
    ω = 0.5

    x, y = x
    co = cos(k * (x + y) - ω * t)
    si = sin(2 * k * (x + y) - 2 * ω * t)
    tmp1 = 2 * k - ω
    tmp2 = - ω + 2 * k * c * κ - k * κ + 3 * ω - 2 * k * c
    tmp3 = k * κ - k
    tmp4 = 4 * k * c * κ + 2 * k - 2 * k * κ - 2 * ω * c
    tmp5 = 2 * k * κ - ω

    du1 = co * tmp1
    du2 = co * tmp2 + si * tmp3
    du3 = co * tmp2 + si * tmp3
    du4 = co * tmp4 + si * tmp5 + sin(k * (x + y) - ω * t) * (2 * k^2 * μ * κ / Pr)

    return SVector(du1, du2, du3, du4)
end

function boundary_condition_convergence_test_viscous(u_inner, orientation, direction, x, t,
                                          surface_flux_function, equations::Trixi.AbstractEquations)
    u_boundary = initial_condition_convergence_test_viscous(x, t, equations)

    # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end

    return flux
end

function solution(x, t)
    # μ = equations.mu
    # Pr = equations.Pr
    # κ = equations.gamma
    # k = 2 # k = c_p * μ / Pr
    # k = 114 * μ / Pr
    k = 0
    c = 4
    ω = 0.5
    ini = sin(k * (x[1] + x[2]) - ω * t) + c

    rho = ini
    rho_v1 = ini
    rho_v2 = ini
    rho_e = ini^2

    return SVector(rho, rho_v1, rho_v2, rho_e)
end

function wrap_array(u_ode::AbstractVector, mesh::Union{TreeMesh{2},CurvedMesh{2},UnstructuredQuadMesh}, equations, dg::DG, cache)
    @boundscheck begin
        @assert length(u_ode) == nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)
    end
    unsafe_wrap(Array{eltype(u_ode),ndims(mesh) + 2}, pointer(u_ode),
                (nvariables(equations), nnodes(dg), nnodes(dg), nelements(dg, cache)))
end

function WR2_initial_condition_convergence_test(x, t, equations::Trixi.AbstractEquations)
    c = 2
    A = 0.1
    L = 2
    f = 1 / L
    ω = 2 * pi * f
    ini = c + A * sin(ω * (x[1] + x[2] - t))

    rho = ini
    rho_v1 = ini
    rho_v2 = ini
    rho_e = ini^2

    return SVector(rho, rho_v1, rho_v2, rho_e)
end

@inline function WR2_source_terms_convergence_test(u, x, t, equations::Trixi.AbstractEquations)
    # Same settings as in `initial_condition`
    c = 2
    A = 0.1
    L = 2
    f = 1 / L
    ω = 2 * pi * f
    γ = equations.gamma

    x1, x2 = x
    si, co = sincos((x1 + x2 - t) * ω)
    tmp1 = co * A * ω
    tmp2 = si * A
    tmp3 = γ - 1
    tmp4 = (2 * c - 1) * tmp3
    tmp5 = (2 * tmp2 * γ - 2 * tmp2 + tmp4 + 1) * tmp1
    tmp6 = tmp2 + c

    du1 = tmp1
    du2 = tmp5
    du3 = tmp5
    du4 = 2 * ((tmp6 - 1) * tmp3 + tmp6 * γ) * tmp1

    return SVector(du1, du2, du3, du4)
end

function WR2_boundary_condition_convergence_test(u_inner, orientation, direction, x, t,
    surface_flux_function,
    equations::Trixi.AbstractEquations)
u_boundary = WR2_initial_condition_convergence_test(x, t, equations)

# Calculate boundary flux
if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
else # u_boundary is "left" of boundary, u_inner is "right" of boundary
flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
end

return flux
end

# semidiscretization of the compressible Euler equations
equations = CompressibleEulerEquations2D(1.4, viscous = true, mu = mu)

# initial_condition = initial_condition_convergence_test_viscous
# source_terms = source_terms_convergence_test_viscous
# boundary_conditions = boundary_condition_convergence_test_viscous

initial_condition = WR2_initial_condition_convergence_test
source_terms = WR2_source_terms_convergence_test
boundary_conditions = WR2_boundary_condition_convergence_test

surface_flux = FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive))
basis = LobattoLegendreBasis(N)
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

# Mapping
function mapping1zu1(xi_, eta_)

    x = xi_
    y = eta_

    return SVector(x, y)
end

cells_per_dimension = (c, c)

# mesh = CurvedMesh(cells_per_dimension, mappingCos, periodicity = true)
# mesh = CurvedMesh(cells_per_dimension, (-1.0, -1.0), (1.0, 1.0))
# mesh = CurvedMesh(cells_per_dimension, (-pi/2, -pi/2), (pi/2, pi/2))
mesh = CurvedMesh(cells_per_dimension, (0.0, 0.0), (2.0, 2.0))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                      # extra_analysis_errors=(:conservation_error,),
                                      # extra_analysis_integrals=(entropy, energy_total,
                                      # energy_kinetic, energy_internal)
                                      )

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        #alive_callback,
                        #save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

###############################################################################
#

# solver = DGSEM(basis, surface_flux, volume_integral)
# initial_condition_sol = initial_condition_convergence_test_viscous
# semi_sol = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sol, solver, boundary_conditions=boundary_conditions)
# ode_sol = semidiscretize(semi_sol, (0.1, 0.1))

# u_sol = wrap_array(ode_sol.u0, mesh, equations, solver, semi_sol.cache)
# u_dg = wrap_array(sol[2], mesh, equations, solver, semi.cache)
# u_diff = abs.(u_sol-u_dg)

# u_diff_max1 = maximum(abs.(u_sol[1,:,:,:]-u_dg[1,:,:,:]))
# u_diff_max2 = maximum(abs.(u_sol[2,:,:,:]-u_dg[2,:,:,:]))
# u_diff_max3 = maximum(abs.(u_sol[3,:,:,:]-u_dg[3,:,:,:]))
# u_diff_max4 = maximum(abs.(u_sol[4,:,:,:]-u_dg[4,:,:,:]))

# maximum(abs.(ode_sol.u0-sol.u[2]))
maximum(abs.(sol.u[1] - sol.u[2]))

max_stars_repo_path: src/WR2/MS5 EOC Test.jl

max_issues_repo_path: src/WR2/MS5 EOC Test.jl

max_forks_repo_path: src/WR2/MS5 EOC Test.jl


text:
# TODO packages in "development" mode symlinked under depot/development/PackageName?

function create_package(name::String, registry::String; path::String = joinpath(registry_depot(registry), "development", name), license::Union{Void, String} = nothing, repo::Union{Void, String} = nothing)
    #depot = registry_depot(registry)
    if isdir(path)
        error("Package path exists")
    else
        mkpath(path)
    end

    # Make README.md, LICENSE.md, dependencies.toml, compatibility.toml, src/, test/
    open(joinpath(path, "README.md"), "w") do io
        write(io, "# $name\n")
    end

    if license === nothing
        error("Must specify license (\"MIT\" or \"proprietary\")")
    end
    if license == "MIT" # TODO we need the name of the copyright holder
        open(joinpath(path, "LICENSE.md"), "w") do io
            write(io, """
            The $name.jl package is licensed under the MIT "Expat" License:

            > Copyright (c) 2018.
            >
            > Permission is hereby granted, free of charge, to any person obtaining a copy
            > of this software and associated documentation files (the "Software"), to deal
            > in the Software without restriction, including without limitation the rights
            > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
            > copies of the Software, and to permit persons to whom the Software is
            > furnished to do so, subject to the following conditions:
            >
            > The above copyright notice and this permission notice shall be included in all
            > copies or substantial portions of the Software.
            >
            > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
            > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
            > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
            > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
            > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
            > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
            > SOFTWARE.
            >
            """)
        end
    elseif license == "proprietary" # TODO we need the name of the copyright holder
        open(joinpath(path, "LICENSE.md"), "w") do io
            write(io, "Copyright (c) 2018.\n\nAll rights reserved.\n")
        end
    end

    open(joinpath(path, "dependencies.toml"), "w") do io
    end

    open(joinpath(path, "compatibility.toml"), "w") do io
    end

    mkdir(joinpath(path, "src"))
    open(joinpath(path, "src", "$name.jl"), "w") do io
        write(io, """
            module $name

            # package code goes here

            end # module
            """)
    end

    mkdir(joinpath(path, "test"))
    open(joinpath(path, "test", "runtests.jl"), "w") do io
        write(io, """
            using Pkg3
            using $name
            using Base.Test

            @test 1 == 0
            """)
    end

    register_package(name, registry; path = path)
end

function register_package(pkgname::String, registry::String; path::String = "", uuid::String = string(Base.Random.uuid4()), shard::Bool = (registry == "Uncurated"))
    # TODO
end


text:
using LSH
using Distances
using DataFrames

function constructLSHTable(params,data::Array{Float64,2})
    T = LSHTable(0.6,LSH.createUHashesAM04(params[:d],4.0,params[:k],
                        params[:L],0.6,params[:m]),
    Vector{Float64}[ vec(data[i,:]) for i=1:size(data,1) ]; progress=true)
end

function constructLSHTable(params,data::Vector{Vector{Float64}})
    T = LSHTable(0.6,LSH.createUHashesAM04(params[:d],4.0,params[:k],
                        params[:L],0.6,params[:m]),data; progress=true)
end

function runQueries(T,data,query)
    if typeof(query) == Array{Float64,2}
        query = Vector{Float64}[ vec(query[i,:]) for i=1:size(query,1) ]
    end
    if typeof(data) == Array{Float64,2}
        data = Vector{Float64}[ vec(data[i,:]) for i=1:size(data,1) ]
    end

    # JIT
    r3 = Dict{Any,Any}()
    for (i,q) in enumerate(query)
        local found_points
        time = @elapsed (found_points = T[q])
        nns = length(found_points)
        push!(r3,i,(time,nns))
    end

    r2 = Dict{Any,Any}()
    for (i,q) in enumerate(query)
        local found_points
        time = @elapsed (found_points = T[q])
        nns = length(found_points)
        push!(r2,i,(time,nns))
    end

    exact = Dict{Any,Any}()
    for (i,q) in enumerate(query)
        count = 0
        time = @elapsed for p in data
            if euclidean(p,q) < 0.6
                count += 1
            end
        end
        push!(exact,i,(time,count))
    end

    (r2,exact)
end

function to_dataframes(a,b,c)
    A = DataFrame(pointno=[i for i = 1:length(a)],
            time=[a[i-1][1] for i = 1:length(a)],
            no=[a[i-1][2] for i = 1:length(a)],
            kind=["E2LSH" for i =1:length(a)])
    B = DataFrame(pointno=[i for i = 1:length(b)],
            time=[b[i][1] for i = 1:length(b)],
            no=[b[i][2] for i = 1:length(b)],
            kind=["LSH.jl" for i =1:length(b)])
    C = DataFrame(pointno=[i for i = 1:length(c)],
            time=[c[i][1] for i = 1:length(c)],
            no=[c[i][2] for i = 1:length(c)],
            kind=["Exact" for i =1:length(c)])
    vcat(A,B,C)
end

function compare_to_E2LSH(data, query, jlremote=false)
    if typeof(query) == Array{Float64,2}
        query = Vector{Float64}[ vec(query[i,:]) for i=1:size(query,1) ]
    end
    if typeof(data) == Array{Float64,2}
        data = Vector{Float64}[ vec(data[i,:]) for i=1:size(data,1) ]
    end
    params, e2lshdata = fetch(@spawnat 2 run_e2lsh_on_data(data,query))
    T = constructLSHTable(params,data)
    r2, exact = runQueries(T,data,query)
    e2lshdata, r2, exact
end

function compare_to_E2LSH_mnist(data, query, jlremote=false)
    params, e2lshdata = remotecall_fetch(2,run_mnist_e2lsh_on_data,data,query)
    data = normalized([ trainfeatures(i) for i in data ])
    query = normalized([ testfeatures(i) for i in query ])
    T = constructLSHTable(params,data)
    r2, exact = runQueries(T,data,query)
    e2lshdata, r2, exact
end


text:
@doc raw"""
    PolyaGammaNegativeMultinomial(y::BitVector, c::AbstractVector{<:Real}, p::AbstractVector{<:Real})

A multivariate distribution, used as hierachical prior as:
```math
    p(\boldsymbol{\omega}, \boldsymbol{n}) = \operatorname{NM}(\boldsymbol{n}|1, \boldsymbol{p})\prod_{i=1}^K\operatorname{PG}(\omega|y_i + n_i, c).
```

Random samples as well as statistics from the distribution will returned as a `NamedTuple` : `(;ω, n)`.

This structured distributions is needed for the [`CategoricalLikelihood`](https://juliagaussianprocesses.github.io/GPLikelihoods.jl/stable/api#GPLikelihoods.CategoricalLikelihood) with a `LogisticSoftMaxLink`.
"""
struct PolyaGammaNegativeMultinomial{Ty,Tc,Tp} <: AbstractNTDist
    y::Ty # Intermediate first parameter for PG(y + n, c)
    c::Tc # Second parameter for PG
    p::Tp # Negative Multinomial parameters
end

function Base.:(==)(p::PolyaGammaNegativeMultinomial, q::PolyaGammaNegativeMultinomial)
    return (p.y == q.y) && (p.c == q.c) && (p.p == q.p)
end

NegativeMultinomial(d::PolyaGammaNegativeMultinomial) = NegativeMultinomial(1, d.p)

Distributions.length(::PolyaGammaNegativeMultinomial) = 2 * length(d.p)

function ntrand(rng::AbstractRNG, d::PolyaGammaNegativeMultinomial)
    n = rand(rng, NegativeMultinomial(d))
    ω = rand.(rng, PolyaGamma.(n + d.y, d.c))
    return (; ω, n)
end

function MeasureBase.logdensity_def(d::PolyaGammaNegativeMultinomial, x::NamedTuple)
    logpdf_n = logpdf(NegativeMultinomial(d), x.n)
    logpdf_ω = sum(1:length(x)) do i
        return logpdf(PolyaGamma(d.y[i] + x.n[i], d.c[i]), x.ω[i])
    end
    return logpdf_ω + logpdf_n
end

function tvmean(ds::AbstractVector{<:PolyaGammaNegativeMultinomial})
    n = nestedview(mapreduce(d -> mean(NegativeMultinomial(d)), hcat, ds))
    ω = nestedview(
        mapreduce(hcat, ds, n) do d, n
            mean.(PolyaGamma.(d.y + n, d.c))
        end,
    )
    return TupleVector(; ω, n)
end

function ntmean(d::PolyaGammaNegativeMultinomial)
    n = mean(NegativeMultinomial(d))
    return (; ω=mean.(PolyaGamma.(d.y + n, d.c)), n)
end

function Distributions.kldivergence(
    q::PolyaGammaNegativeMultinomial, p::PolyaGammaNegativeMultinomial
)
    # TODO: Optimize this
    (all(==(0), p.c) && all(p.y .== q.y)) ||
        error("no kl divergence available for this prior. qΩ = $q, pΩ = $p")
    b = q.y + mean(NegativeMultinomial(q))
    return sum(kldivergence.(PolyaGamma.(b, q.c), PolyaGamma.(b, 0))) +
           kldivergence(NegativeMultinomial(q), NegativeMultinomial(p))
end


text:
## LLE pressure solver
function x0_LLE_pressure(model::EoSModel,T,x)
    xx = 1 .-x

    pure = split_model(model)

    eachx = eachcol(Diagonal(ones(eltype(x),length(x))))

    sat = sat_pure.(pure,T)

    P_sat = [tup[1] for tup in sat]
    V_l_sat = [tup[2] for tup in sat]

    V0_l = sum(x.*V_l_sat)
    V0_ll = sum(xx.*V_l_sat)

    prepend!(xx,log10.([V0_l,V0_ll]))
    return xx
end

function LLE_pressure(model::EoSModel, T, x; v0 =nothing)
    TYPE = promote_type(eltype(T),eltype(x))
#     lb_v = lb_volume(model,x)
    ts = T_scales(model,x)
    pmix = p_scale(model,x)
    if v0 === nothing
        v0 = x0_LLE_pressure(model,T,x)
    end
    len = length(v0[1:end-1])
    #xcache = zeros(eltype(x0),len)
    Fcache = zeros(eltype(v0[1:end-1]),len)
    f! = (F,z) -> Obj_bubble_pressure(model, F, T, exp10(z[1]), exp10(z[2]), x,z[3:end],ts,pmix)
    r  =Solvers.nlsolve(f!,v0[1:end-1],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    v_l = exp10(sol[1])
    v_ll = exp10(sol[2])
    xx = FractionVector(sol[3:end])
    P_sat = pressure(model,v_l,T,x)
    return (P_sat, v_l, v_ll, xx)
end


text:
using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

num_trials = 100::Int64
perc = 1.
grid_size = 10::Int64
time_depth = 8::Int64
num_pairs = 50::Int64
asynchronus_weight = 100::Int64

##### MEMORY RATE #####
G = GridNetwork(grid_size, grid_size)

# Extend graph in time without memory

raw_data = []
for j in 1:num_trials
    # Extend graph in time without memory
    T = QuNet.TemporalGraph(G, time_depth, memory_prob=0.0)
    # Make a network with some probability of memory
    T_mem = QuNet.TemporalGraph(G, time_depth, memory_prob=1.0, memory_costs = unit_costvector())

    # Get i random userpairs with asynchronus src and dst nodes.
    mem_user_pairs = make_user_pairs(T, num_pairs, src_layer=1, dst_layer=-1)
    user_pairs = make_user_pairs(T, num_pairs, src_layer=-1, dst_layer=-1)

    # Add async nodes
    QuNet.add_async_nodes!(T_mem, user_pairs, ϵ=asynchronus_weight)
    QuNet.add_async_nodes!(T, user_pairs, ϵ=asynchronus_weight)

    # Get pathset data
    # NOTE: the line below used to have mem_user_pairs as an argument.
    pathset_mem, dum1, dum2 = QuNet.greedy_multi_path!(T_mem, QuNet.purify, mem_user_pairs, 4)
    pathset, dum1, dum2 = QuNet.greedy_multi_path!(T, QuNet.purify, user_pairs, 4)

    max_depth_mem = QuNet.max_timedepth(pathset_mem, T)
    max_depth = QuNet.max_timedepth(pathset, T)

    # Get the bandwidth of this quantity
    push!(raw_data, max_depth / max_depth_mem)
end

println("Mean raw_data")
println(mean(raw_data))

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule librtlsdr_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("librtlsdr")
JLLWrappers.@generate_main_file("librtlsdr", UUID("d1c10d9e-0af6-511b-a44d-62fc8e59b084"))
end  # module librtlsdr_jll


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 75fe8fd0-2563-11ec-30ef-d15cdc3d4e63
using PlutoUI

# ╔═╡ d4b5618f-f1a7-4483-b738-fb94e1e2bbbe
html"""

&nbsp&nbsp💻&nbspLógica de Programação por Franco Naghetini é licenciado sob CC BY 4.0

"""

# ╔═╡ deee1735-37ee-4c97-a4df-8fcb61e95d6a
PlutoUI.TableOfContents(aside=true, title="Sumário",
						indent=true, depth=2)

# ╔═╡ 488dcacd-109d-41f2-b904-3d17193e6190
md"""
![ufmg-logo](https://logodownload.org/wp-content/uploads/2015/02/ufmg-logo-2.png)
"""

# ╔═╡ 6ca3e113-02db-4cef-ad9e-3941ac7d7a6d
md"""
# 💻 Lógica de Programação

O ato de **programar** pode ser entendido como uma forma de se comunicar com as máquinas a partir de um conjunto de instruções não ambíguas com o intuito de se realizar uma determinada tarefa. Essa comunicação com as máquinas é realizada por meio das **linguagens de programação**, como Python, R, C e Julia.

Nesse contexto, a **lógica de programação** é de suma importância, uma vez que ela define a estrutura das instruções dadas às máquinas, ou seja, define a forma que você se comunica com as máquinas. Se quer aprender mais sobre a importância da lógica de programação assista [este vídeo](https://www.youtube.com/watch?v=l7lIPXij85I) do canal Programação Dinâmica.

Neste módulo, estudaremos alguns conceitos básicos de lógica de programação. Para isso, utilizaremos [Julia](https://julialang.org/), uma linguagem de programação moderna e de propósito geral que possui diversos recursos favoráveis ao desenvolvimento de rotinas geoestatísticas.

> A linguagem é *simples* como Python e *rápida* como C 🚀

A seguir, serão apresentados alguns conceitos básicos de lógica de programação, além de uma introdução sobre os recursos interativos da linguagem Julia.
"""

# ╔═╡ 14db5525-7c0a-433b-a23c-088db728f46b
md"""
>##### 📚 Sobre
>- Você pode exportar este notebook como PDF ou HTML estático. Para isso, clique no ícone 🔺🔴, localizado no canto superior direito da pagina. Entretanto, ambos os formatos não são compatíveis com os recursos interativos do notebook.
>- Caso deseje executar alguma célula do notebook, clique no ícone ▶, localizado no canto inferior direito da célula.
>- Algumas células encontram-se ocultadas (e.g. células que geram os plots). Você pode clicar no ícone 👁, localizado no canto superior esquerdo da célula, para ocultá-la ou exibí-la.
>- A explicação das células que geram os plots está fora do escopo deste notebook. Entretanto, a sintaxe é bem intuitiva e pode ser facilmente compreendida!
>- Você pode ainda clicar no ícone `...`, no canto superior direito de uma célula, para excluí-la do notebook.
>- Algumas células deste notebook encontram-se encapsuladas pela expressão `md"..."` (e.g. esta célula). Essas são células de texto chamadas de *markdown*. Caso deseje aprender um pouco mais sobre a linguagem *markdown*, clique [aqui](https://docs.pipz.com/central-de-ajuda/learning-center/guia-basico-de-markdown#open).
>- No Pluto, todos os pacotes devem ser importados/baixados na primeira célula do notebook. Clique no ícone 👁 para exibir essa célula ou consulte a seção *Pacotes utilizados* deste notebook para saber mais informações sobre os pacotes.
>- Utilize a macro ` @which` para verificar a qual pacote uma determinada função pertence.
>- Você pode utilizar este notebook da forma que quiser, basta referenciar [este link](https://github.com/fnaghetini/intro-to-geostats). Consulte a [licença]  (https://github.com/fnaghetini/intro-to-geostats/blob/main/LICENSE) para saber mais detalhes.
>- Para mais informações acesse o [README](https://github.com/fnaghetini/intro-to-geostats/blob/main/README.md) do projeto 🚀
"""

# ╔═╡ c111da75-d294-4def-93fb-56953c3585ad
md"""

## 1. Variáveis

As **variáveis** podem ser entendidas como partes da memória (ou "caixas") onde se armazena **valores** de diferentes **tipos** para serem posteriormente processados. Cada variável é identificada/rotulada por um **nome**.

As variáveis podem ser **numéricas**...

"""

# ╔═╡ 7ceb9c27-0310-4f97-8349-91286c1f9235
idade = 24

# ╔═╡ 218bb8c3-729b-43c7-9a53-d3b579cf2d21
altura = 1.58

# ╔═╡ 3b3e6955-1573-4620-8109-d80b71cf0708
md" As variáveis podem armazenar **cadeias de caracteres**... "

# ╔═╡ f2e6adcd-9095-496f-a8e3-88069537015c
nome = "Camila"

# ╔═╡ 870cbb76-b0c3-4d14-be22-86e0d34ebe58
frase = "Olá, mundo! 🌎"

# ╔═╡ 55a2c622-0f9f-45db-b534-55873c0759d4
md"""
> ⚠ As cadeias de caracteres (**strings**) devem ser encapsuladas por áspas duplas. Emojis também são considerados cadeias de caracteres! Clique [aqui](https://getemoji.com/) para copiar emojis e colar nas células deste notebook.
"""

# ╔═╡ 449c91c3-7f03-4ffd-95e7-02cdb58323fd
md"""
Um recurso muito utilizado e presente em diversas linguagens de programação é a **interpolação de strings**. Esse procedimento utiliza o valor de uma variável dentro de uma string, de modo que a string interpolada atua como um "template". Em Julia, usamos o símbolo `$` para interpolar strings.

A célula abaixo mostra um exemplo de interpolação de string, utilizando as variáveis já criadas `nome`, `idade` e `altura`...
"""

# ╔═╡ 8d463f50-d694-41eb-bbad-c40f93471852
"Olá! Meu nome é $nome, tenho $idade anos e tenho $altura m de altura."

# ╔═╡ 4d31b989-7a7d-43ca-a1f8-831db58f1c8e
md"""

Podemos utilizar *símbolos matemáticos* como nomes das variáveis. Para inserir o símbolo $\alpha$, por exemplo, digita-se `\alpha` + `TAB`. No caso do símbolo $\gamma$, digita-se `\gamma` + `TAB`.

Cique [aqui](https://docs.julialang.org/en/v1/manual/unicode-input/) para acessar a lista completa de símbolos disponíveis.

"""

# ╔═╡ b0d01238-1213-45c3-9bef-a443afc6b87c
β = 2.0

# ╔═╡ 2cce4c5e-7890-4c6c-bd7c-eecc8437bef3
ϵ = 0.5

# ╔═╡ 1b454935-c858-4817-928d-da116011d7b0
β + ϵ

# ╔═╡ 96ced3bd-69e7-490c-b8be-b06a30049a31
π

# ╔═╡ 3de2d16f-8849-4462-9ef3-0cc1830285c6
md"""
> ⚠ Perceba que o símbolo $π$ já possui um valor associado na linguagem Julia.
"""

# ╔═╡ d28119c3-6800-430d-9ffe-d6f25e8e5c2a
md"""
## 2. Funções

As **funções**, no mundo da computação, são similares às funções matemáticas. De forma geral, uma função recebe valores de **entrada**, processa-os e retorna valores de **saída**.

Em Julia, existem, essencialmente, três maneiras distintas de se construir funções:

"""

# ╔═╡ bb495116-8ff9-4a66-9b66-6f6fe990c5d9
f(χ) = 5χ + 25

# ╔═╡ a63d4a04-0ceb-4e4f-90ba-30dc798e1b48
function areatriangulo(b,h)
	A = (b * h) / 2
	return A
end

# ╔═╡ 9cf1d023-1421-4adb-af35-64aa6a8ddb53
md"""
> ⚠ O comando `return` foi utilizado aqui apenas por fins didáticos. Entretanto, poderíamos omiti-lo.
"""

# ╔═╡ 5b168a45-b5a2-4415-80f4-e2347c21730a
diagonal = l -> l * √2

# ╔═╡ 1a1ac9a0-b8f0-46d1-9d89-21b6bc0e6aec
f(2) + areatriangulo(10,4) / diagonal(1)

# ╔═╡ 784426f8-ac63-4555-adb2-cbfb9b37fdda
md"""
#### 🖊 Exercício

Crie uma função `recurso` que retorne a tonelagem em metal contido de um depósito de volume `V`, densidade média `ρ` e teor médio `T` (em %).

"""

# ╔═╡ ffc5e222-8ae8-4314-a616-e4588d7808fe
recurso(V, ρ, T) = missing

# ╔═╡ 5dc76f3b-badf-4951-9eb4-f61ca66dfb98
md"""
## 3. Coleções

As **coleções** são conjuntos de elementos, normalmente do mesmo tipo, mas não necessariamente. São exemplos *vetores*, *matrizes* e *tensores*.

Os **vetores** podem ser escritos como...

"""

# ╔═╡ d349ac4d-44bb-4505-a6c4-3096d8c48ba0
vetor1 = [1,2,3,4,5,6,7,8,9,10]

# ╔═╡ 9e181c4a-a759-4f92-8444-6a23387a20b9
vetor_2 = collect(10:20)

# ╔═╡ d39f503a-158e-4774-851e-d786014b8f2b
md" As **matrizes** podem ser escritas como... "

# ╔═╡ b5387324-df28-422d-8e74-4fefa8054e7f
A = [1 2
	 3 4]

# ╔═╡ 570f62db-edd7-429d-a084-2686f4b6923c
B = [1 2 3; 4 5 6; 7 8 9]

# ╔═╡ 7e8a4e33-0e9c-45f5-b8bc-e2adece0806d
md" Os **tensores** podem ser escritos como... "

# ╔═╡ 204bcedc-4123-4845-bcac-493ce3544582
tensor_de_zeros = zeros(5,5,2)

# ╔═╡ 92b0cb5f-41a1-420f-8b81-3b967db8138f
tensor_de_uns = ones(2,2,3)

# ╔═╡ 7471120d-6e00-46be-be11-649c8cdbb32c
md"""
Uma operação muito utilizada é o **fatiamento**, que visa obter apenas certo(s) elemento(s) de uma coleção. O fatiamento de vetores, por exemplo, é realizado a partir da seguinte sintaxe:

```julia
vetor[índice inicial:índice final]
```

Observe o vetor `v₁`...

"""

# ╔═╡ b31b883d-5a69-4e04-b8b5-cedca81dcf9a
v₁ = [10,20,30,40,50,60,70,80,90,100]

# ╔═╡ c7f94dcd-64a2-4e97-b38d-363dc6f9f757
md" O primeiro elemento pode ser obtido a partir do fatiamento de `v₁`..."

# ╔═╡ 2a25a884-44cf-4460-ab23-02e728201f90
v₁[1]

# ╔═╡ e127bdf8-2a53-4534-a628-1463a5735f41
md"""

> ⚠ Assim como R e diferentemente de Python ou JavaScript, a numeração do índice se inicia em **1** na linguagem Julia.

Os três primeiros elementos de `v₁` podem ser obtidos como...

"""

# ╔═╡ 0806feb3-33eb-402c-92d2-af0b153834ba
v₁[1:3]

# ╔═╡ 1a5f467b-7244-4587-b01c-61fd22bfbb9b
v₁[begin:3]

# ╔═╡ 4b1ec558-5f94-417a-aa3c-d13a0039a6a0
md" Os dois últimos elementos de `v₁` podem ser obtidos como..."

# ╔═╡ 6bee210a-11cc-40c4-bfb0-d14a6b2495b4
v₁[9:10]

# ╔═╡ 10cc97dd-a390-4d4a-9644-20b654ff66bb
v₁[end-1:end]

# ╔═╡ 04738501-16fa-43e1-b640-47682c155bcb
md"""

#### 🖊 Exercício

Fatie o vetor `v₂`, de modo que apenas os elementos da terceira posição em diante sejam retornados. A fatia resultante deve ser armazenada na variável `v₃`.

"""

# ╔═╡ f9317c0c-2389-43a8-8b91-f8b2f278e258
begin
	v₂ = collect(10:15)
	v₃ = missing
end

# ╔═╡ 8253149c-e1b0-4e3f-9f5b-909621dba86e
md"""
## 4. Condicionais

As **estruturas condicionais** são utilizadas em situações em que se deseja executar algum trecho de código apenas quando uma **condição** é satisfeita:

```julia
if condição
	comando_1
else
	comando_2
end
```

O exemplo abaixo mostra uma função `égranitoide` que utiliza estruturas condicionais. **Se** o conteúdo em quartzo for menor que 20%, a função retorna a frase "Não é granitoide 😐". **Caso contrário**, a frase "É granitoide!" é retornada.

"""

# ╔═╡ 4a778055-7ade-4275-aa87-95c121b31232
function égranitoide(quartzo)
	if quartzo < 20
   		rochaignea = "Não é granitoide 😐"
	else
   		rochaignea = "É granitoide!"
	end

	return rochaignea
end

# ╔═╡ ca5453f5-1768-4f59-b09d-e40bb4981622
égranitoide(30)

# ╔═╡ 217945a3-919b-4c5f-afb0-6eaadf4e9022
égranitoide(5)

# ╔═╡ 1f0abfbe-49cd-4529-8f4d-65f8286b6f0e
md"""

#### 🖊 Exercício

Agora, crie uma função `tiporocha` que recebe o nome de uma rocha e retorna o seu tipo:

* "gabro" → "ígnea"
* "gnaisse" → "metamórfica"
* "ritmito" → "sedimentar"

"""

# ╔═╡ 6fcc1433-c348-4825-8503-782304ed508b
function tiporocha(rocha)
	missing
end

# ╔═╡ 3cfbddfb-be59-4019-965f-13bdbac2295a
md"""

## 5. Laços de repetição

Os **laços de repetição** são utilizados quando se deseja repetir determinado trecho do algoritmo múltiplas vezes. O número de repetições pode ser indeterminado, mas necessariamente finito (*Dauricio, 2015*).

Uma das repetições mais utilizadas é o `for`, que apresenta a seguinte sintaxe:

```julia
for elemento in vetor
	comando
end
```

O comando `for` também pode ser utilizado para criar vetores/listas, a partir das chamadas **list comprehensions**. A sintaxe é definida da seguinte maneira:

```julia
[comando for elemento in vetor]
```

O exemplo abaixo mostra a utilização do recurso *list comprehension* para a criação de um vetor constituído por potências de 2:

"""

# ╔═╡ abc428c6-b762-47bc-b5ac-7abedd1a7021
[2 ^ num for num in 1:10]

# ╔═╡ bc8fa2b9-0803-40e7-af91-2196a11aae44
md"""

#### 🖊 Exercício

Escreva a função `raiz` que retorna a raíz quadrada de todos os elementos do vetor `v₄`.

"""

# ╔═╡ 729eaff4-7aa1-4db1-a04f-57216175b029
raiz(v₄) = missing

# ╔═╡ 5eb8ff96-2adf-465d-bf8c-bd4e64fa6342
md"""
Podemos ainda utilizar a sintaxe *list comprehension* para criarmos uma tabuada em apenas uma linha!
"""

# ╔═╡ 0d420048-c4f9-4149-b25b-149ef77f3264
[linha * coluna for linha in 1:10, coluna in 1:10]

# ╔═╡ f8f272e0-17ca-47e5-8dd0-577eb6322c90
md"""

## 6. Interatividade

O ambiente [Pluto](https://github.com/fonsp/Pluto.jl) e o pacote [PlutoUI](https://github.com/fonsp/PlutoUI.jl) proporciona uma série de recursos interativos a partir da macro `@bind`. Esses recursos permitem que o valor de uma variável seja alterado de acordo com a interação entre o usuário e o notebook.

O recurso **Slider** funciona como uma barra de deslizamento para variáveis que assumem valores contínuos. A sua sintaxe é dada por:

```julia
@bind nome_var Slider(inicio:passo:fim)
```

O *slider* abaixo mostra o teor de Cu (%) em uma amostra...

"""

# ╔═╡ e137bf5e-d795-49ac-b641-b5f7364cfac9
@bind CU Slider(0.05:0.01:3.00, default=1.0, show_value=true)

# ╔═╡ fae5e284-66e7-4599-9905-f0d5e64387ea
md" A amostra apresenta $CU % de Cu"

# ╔═╡ 9ee680a7-2866-4f1b-84af-0fdb737246ad
md"""
Um outro elemento interativo é o **TextField**, que se refere a uma caixa de texto. Sua sintaxe é definida como:

```julia
@bind nome_var TextField()
```

"""

# ╔═╡ 6c7e7d3e-d544-4b59-9504-e28fe2cc3619
md"""

#### 🖊 Exercício

Com seus conhecimentos geológicos, tente identificar qual rocha é exibida na Figura 01. Escreva o nome da rocha (apenas com letras minúsculas) na caixa de texto abaixo.

"""

# ╔═╡ 2fd0708e-b225-4f8c-8f4f-b46305e06366
md"""
![rocha](http://lh3.googleusercontent.com/-4dxg2gAuqgk/VeqttO0BqSI/AAAAAAAAKq0/j_G5leYIfIg/s720/01465%252520IMG_20150809_104742%252520mylonite.jpg)

_**Figura 01:** Identifique a rocha acima. Extraído de [sandatlas.org](https://www.sandatlas.org/mylonite/)._
"""

# ╔═╡ 98a32d2d-7173-406d-92cf-dfa847d47c49
md" Rocha: $(@bind rocha TextField())"

# ╔═╡ df1f77ac-bbf5-4d42-91bf-e634ebee4bb9
md""" Podemos também utilizar o recurso interativo **Select** que, por sua vez, atua como uma lista suspensa. Esse comando possui a seguinte sintaxe:

```julia
@bind nome_var Select(["elem_1", "elem_2", ..., "elem_n"])
```

"""

# ╔═╡ ba7ce05c-cae6-4eec-9760-e4e253986061
md"""

#### 🖊 Exercício

Na lista suspensa seguir, selecione a opção que corresponda à correta subclasse do silicato exibido na Figura 02.

"""

# ╔═╡ b844b981-2922-4056-811d-8c61b904d996
md"""
![silicato](https://geology.com/minerals/photos/titanite-crystal.jpg)

_**Figura 02:** Identifique a subclasse do silicato acima. Extraído de: [geology.com](https://geology.com/minerals/titanite.shtml)._
"""

# ╔═╡ 87832d6c-525d-4060-a1e8-d06a611b5dfe
md""" Subclasse: $(@bind silicato Select(["Tectossilicato","Filossilicato",
										  "Inossilicato","Ciclossilicato",
										  "Sorossilicato","Nesossilicato"]))
"""

# ╔═╡ a93dd751-a6b3-454d-9f70-b864fdbbd968
md"""
## Referências

*Dauricio, J. S. **Algoritmos e Lógica de Programação**. Londrina: Editora e Distribuidora Educacional S.A., 2015.*
"""

# ╔═╡ 8f866ee4-ce44-4d42-be56-7764168c1c71
md"""
## Recursos adicionais

Abaixo, são listados alguns recursos complementares a este notebook:

> [Documentação da Linguagem Julia](https://docs.julialang.org/en/v1/)

> [YouTube - Julia Language](https://www.youtube.com/c/TheJuliaLanguage)

> [YouTube - Programação Dinâmica](https://www.youtube.com/playlist?list=PL5TJqBvpXQv4TO4Y_JZExBCBzNXDLYUCj)
"""

# ╔═╡ ce463e41-8cb2-46cb-9656-1503997a883e
md"""
## Pacotes utilizados

Os seguintes pacotes foram utilizados neste notebook:

|                    Pacote                     |      Descrição      |
|:---------------------------------------------:|:-------------------:|
|[PlutoUI](https://github.com/fonsp/PlutoUI.jl) | Widgets interativos |

"""

# ╔═╡ bb899ab7-b3d8-493f-aaa6-85a4710a6690
begin
	hint(text) = Markdown.MD(Markdown.Admonition("hint", "Dica", [text]))

	almost(text) = Markdown.MD(Markdown.Admonition("warning", "Quase lá!", [text]))

	still_missing(text=md"Troque `missing` pela sua resposta.") = Markdown.MD(Markdown.Admonition("warning", "Aqui vamos nós!", [text]))

	keep_working(text=md"A resposta não está correta 😔") = Markdown.MD(Markdown.Admonition("danger", "Não desanime, você está quase lá!", [text]))

	yays = [md"Fantástico!", md"Ótimo!", md"Yay ❤", md"Legal! 🎉", md"Muito bem!", md"Bom trabalho!", md"Você conseguiu a resposta certa!", md"Vamos seguir para próxima seção."]

	correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Certa resposta!", [text]))

	not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Opa!", [md"Tenha certeza que definiu uma variável chamada **$(Markdown.Code(string(variable_name)))**"]))
end;

# ╔═╡ a4f3b825-6806-4c64-bbe0-6d7383015d68
begin
	s1 = false
	_rec = recurso(12500000,2.7,5)
	if ismissing(_rec)
		still_missing()
	elseif _rec ≈ (12500000 * 2.7 * 5) / 100
		s1 = true
		correct()
	elseif _rec isa Number
		almost(md"A fórmula não está certa...")
	else
		keep_working()
	end
end

# ╔═╡ 3ccac0bc-9f78-4ece-988b-ba832811e538
hint(md"Alguém me contou que a fórmula é $\frac{VρT}{100}$...")

# ╔═╡ a84b90c1-fa85-499a-9fa4-ac1f33ed8469
begin
	s2 = false
	if ismissing(v₃)
		still_missing()
	elseif v₃ == v₂[3:end]
		s2 = true
		correct()
	else
		keep_working()
	end
end

# ╔═╡ e2c42c9e-abbb-4a33-844d-ef0f0c9b4939
hint(md"Utilize `end` como índice final.")

# ╔═╡ 33b81b24-ca33-4dd3-a85e-0371ca8623e8
begin
	s3 = false
	_rcktype = tiporocha.(["gabro","gnaisse","ritmito"])
	if all(ismissing.(_rcktype))
		still_missing()
	elseif all(_rcktype .== ["ígnea","metamórfica","sedimentar"])
		s3 = true
		correct()
	elseif _rcktype ⊆ ["ígnea","metamórfica","sedimentar"]
		almost(md"A resposta não está 100% correta...")
	else
		keep_working()
	end
end

# ╔═╡ e462ca14-b38a-46b0-8faa-2261a2cfc150
hint(md"Basta escrever uma sequência de `if rocha == \"gabro\" return \"ígnea\" end`")

# ╔═╡ 7d12857d-d767-49e5-b7de-11e696f1d84a
begin
	s4 = false
	_sqrt = raiz([1 4; 9 16])
	if ismissing(_sqrt)
		still_missing()
	elseif _sqrt == [1 2; 3 4]
		s4 = true
		correct()
	elseif _sqrt == [1, 2, 3, 4]
		almost(md"Tente usar a dica!")
	else
		keep_working()
	end
end

# ╔═╡ a9f386e1-8c6e-40e7-ba67-b3792c50f87e
hint(md"A notação list comprehension `[f(x) for v in v₄]` pode ser bem útil!")

# ╔═╡ f1080b76-7564-4419-884f-70a49bfe9d57
begin
	s5 = false
	if rocha == ""
		almost(md"Identifique a rocha!")
	elseif rocha == "milonito"
		s5 = true
		correct()
	elseif rocha == "protomilonito" || rocha == "ultramilonito"
		almost(md"""A porcentagem de matriz "triturada" é de 50% a 90%.""")
	else
		keep_working()
	end
end

# ╔═╡ 844171ee-dfc3-4985-a6a1-7aa7986d1ac4
hint(md"Essa rocha é um típico produto de recristalização dinâmica!")

# ╔═╡ 8b0df84d-a112-49ec-9f7b-c6e5efd45c48
begin
	s6 = false
	if silicato == "Nesossilicato"
		s6 = true
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 55298aee-618b-4704-98ff-e81fdbddd9a7
hint(md"A fórmula química deste silicato é `CaTiO(SiO₄)`.")

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─75fe8fd0-2563-11ec-30ef-d15cdc3d4e63
# ╟─d4b5618f-f1a7-4483-b738-fb94e1e2bbbe
# ╟─deee1735-37ee-4c97-a4df-8fcb61e95d6a
# ╟─488dcacd-109d-41f2-b904-3d17193e6190
# ╟─6ca3e113-02db-4cef-ad9e-3941ac7d7a6d
# ╟─14db5525-7c0a-433b-a23c-088db728f46b
# ╟─c111da75-d294-4def-93fb-56953c3585ad
# ╠═7ceb9c27-0310-4f97-8349-91286c1f9235
# ╠═218bb8c3-729b-43c7-9a53-d3b579cf2d21
# ╟─3b3e6955-1573-4620-8109-d80b71cf0708
# ╠═f2e6adcd-9095-496f-a8e3-88069537015c
# ╠═870cbb76-b0c3-4d14-be22-86e0d34ebe58
# ╟─55a2c622-0f9f-45db-b534-55873c0759d4
# ╟─449c91c3-7f03-4ffd-95e7-02cdb58323fd
# ╠═8d463f50-d694-41eb-bbad-c40f93471852
# ╟─4d31b989-7a7d-43ca-a1f8-831db58f1c8e
# ╠═b0d01238-1213-45c3-9bef-a443afc6b87c
# ╠═2cce4c5e-7890-4c6c-bd7c-eecc8437bef3
# ╠═1b454935-c858-4817-928d-da116011d7b0
# ╠═96ced3bd-69e7-490c-b8be-b06a30049a31
# ╟─3de2d16f-8849-4462-9ef3-0cc1830285c6
# ╟─d28119c3-6800-430d-9ffe-d6f25e8e5c2a
# ╠═bb495116-8ff9-4a66-9b66-6f6fe990c5d9
# ╠═a63d4a04-0ceb-4e4f-90ba-30dc798e1b48
# ╟─9cf1d023-1421-4adb-af35-64aa6a8ddb53
# ╠═5b168a45-b5a2-4415-80f4-e2347c21730a
# ╠═1a1ac9a0-b8f0-46d1-9d89-21b6bc0e6aec
# ╟─784426f8-ac63-4555-adb2-cbfb9b37fdda
# ╠═ffc5e222-8ae8-4314-a616-e4588d7808fe
# ╟─a4f3b825-6806-4c64-bbe0-6d7383015d68
# ╟─3ccac0bc-9f78-4ece-988b-ba832811e538
# ╟─5dc76f3b-badf-4951-9eb4-f61ca66dfb98
# ╠═d349ac4d-44bb-4505-a6c4-3096d8c48ba0
# ╠═9e181c4a-a759-4f92-8444-6a23387a20b9
# ╟─d39f503a-158e-4774-851e-d786014b8f2b
# ╠═b5387324-df28-422d-8e74-4fefa8054e7f
# ╠═570f62db-edd7-429d-a084-2686f4b6923c
# ╟─7e8a4e33-0e9c-45f5-b8bc-e2adece0806d
# ╠═204bcedc-4123-4845-bcac-493ce3544582
# ╠═92b0cb5f-41a1-420f-8b81-3b967db8138f
# ╟─7471120d-6e00-46be-be11-649c8cdbb32c
# ╠═b31b883d-5a69-4e04-b8b5-cedca81dcf9a
# ╟─c7f94dcd-64a2-4e97-b38d-363dc6f9f757
# ╠═2a25a884-44cf-4460-ab23-02e728201f90
# ╟─e127bdf8-2a53-4534-a628-1463a5735f41
# ╠═0806feb3-33eb-402c-92d2-af0b153834ba
# ╠═1a5f467b-7244-4587-b01c-61fd22bfbb9b
# ╟─4b1ec558-5f94-417a-aa3c-d13a0039a6a0
# ╠═6bee210a-11cc-40c4-bfb0-d14a6b2495b4
# ╠═10cc97dd-a390-4d4a-9644-20b654ff66bb
# ╟─04738501-16fa-43e1-b640-47682c155bcb
# ╠═f9317c0c-2389-43a8-8b91-f8b2f278e258
# ╟─a84b90c1-fa85-499a-9fa4-ac1f33ed8469
# ╟─e2c42c9e-abbb-4a33-844d-ef0f0c9b4939
# ╟─8253149c-e1b0-4e3f-9f5b-909621dba86e
# ╠═4a778055-7ade-4275-aa87-95c121b31232
# ╠═ca5453f5-1768-4f59-b09d-e40bb4981622
# ╠═217945a3-919b-4c5f-afb0-6eaadf4e9022
# ╟─1f0abfbe-49cd-4529-8f4d-65f8286b6f0e
# ╠═6fcc1433-c348-4825-8503-782304ed508b
# ╟─33b81b24-ca33-4dd3-a85e-0371ca8623e8
# ╟─e462ca14-b38a-46b0-8faa-2261a2cfc150
# ╟─3cfbddfb-be59-4019-965f-13bdbac2295a
# ╠═abc428c6-b762-47bc-b5ac-7abedd1a7021
# ╟─bc8fa2b9-0803-40e7-af91-2196a11aae44
# ╠═729eaff4-7aa1-4db1-a04f-57216175b029
# ╟─7d12857d-d767-49e5-b7de-11e696f1d84a
# ╟─a9f386e1-8c6e-40e7-ba67-b3792c50f87e
# ╟─5eb8ff96-2adf-465d-bf8c-bd4e64fa6342
# ╠═0d420048-c4f9-4149-b25b-149ef77f3264
# ╟─f8f272e0-17ca-47e5-8dd0-577eb6322c90
# ╟─e137bf5e-d795-49ac-b641-b5f7364cfac9
# ╟─fae5e284-66e7-4599-9905-f0d5e64387ea
# ╟─9ee680a7-2866-4f1b-84af-0fdb737246ad
# ╟─6c7e7d3e-d544-4b59-9504-e28fe2cc3619
# ╟─2fd0708e-b225-4f8c-8f4f-b46305e06366
# ╟─98a32d2d-7173-406d-92cf-dfa847d47c49
# ╟─f1080b76-7564-4419-884f-70a49bfe9d57
# ╟─844171ee-dfc3-4985-a6a1-7aa7986d1ac4
# ╟─df1f77ac-bbf5-4d42-91bf-e634ebee4bb9
# ╟─ba7ce05c-cae6-4eec-9760-e4e253986061
# ╟─b844b981-2922-4056-811d-8c61b904d996
# ╟─87832d6c-525d-4060-a1e8-d06a611b5dfe
# ╟─8b0df84d-a112-49ec-9f7b-c6e5efd45c48
# ╟─55298aee-618b-4704-98ff-e81fdbddd9a7
# ╟─a93dd751-a6b3-454d-9f70-b864fdbbd968
# ╟─8f866ee4-ce44-4d42-be56-7764168c1c71
# ╟─ce463e41-8cb2-46cb-9656-1503997a883e
# ╟─bb899ab7-b3d8-493f-aaa6-85a4710a6690
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: fnaghetini/intro-to-geostats

max_issues_repo_name: fnaghetini/intro-to-geostats

max_forks_repo_name: fnaghetini/intro-to-geostats


text:
module BHTsne

#using ArgParse
#using StrPack
using PyCall
using Statistics
using LinearAlgebra
# @pyimport struct # error, replace with
pystruct = pyimport("struct")

const BH_TSNE_path = joinpath(dirname(@__FILE__), "cpp/bh_tsne")

export bh_tsne

function __init__()
    # if OS_NAME == :Windows
    #     error("Sorry, BHTsne.jl doesn't support Windows now")
    # end
    if !isfile(BH_TSNE_path)
        cd(dirname(BH_TSNE_path)) do
            run(`g++ sptree.cpp tsne.cpp -o bh_tsne -O2`)
        end
    end
end

function read_unpack(fmt, f)
    # struct.unpack(fmt, read(f, struct.calcsize(fmt)))
    pystruct[:unpack](fmt, read(f, pystruct[:calcsize](fmt)))
end

function bh_tsne(samples;no_dims=2, initial_dims=50, perplexity=50,
                 theta=0.5, randseed=-1, verbose=false)

    samples = broadcast(-, samples, mean(samples))
    cov_x = samples' * samples
    eig_val,eig_vec = eigen(cov_x)

    eig_vec = eig_vec[:,sortperm(eig_val,rev=true)]

    num_eigvec = size(eig_vec, 2)
    if initial_dims > num_eigvec
        initial_dims = num_eigvec
    end

    eig_vec = eig_vec[:,1:initial_dims]
    samples = samples * eig_vec

    sample_count, sample_dim = size(samples)

    mktempdir() do temp_dir

        open(joinpath(temp_dir, "data.dat"),"w") do data_file
            write(data_file, pystruct[:pack]("iiddi", sample_count,sample_dim, theta, perplexity,no_dims))
            nrow,ncol = size(samples)
            fmt = repeat("d",ncol)
            for i = 1:nrow
                # data_packed = mapreduce(x->struct.pack("d",x), *, samples[i,:])
                data_packed = mapreduce(x->pystruct[:pack]("d",x), *, samples[i,:])
                write(data_file, data_packed)
            end
            if randseed != -1
                # write(data_file, struct.pack("i",randseed))
                write(data_file, pystruct[:pack]("i",randseed))
            end
        end

        cd(temp_dir) do
            if verbose
                try
                    run(pipeline(`$BH_TSNE_path`, stdout=STDERR))
                catch excp
                    warn(excp)
                    error("ERROR: Call to bh_tsne exited with a non-zero return code exit status,
                          please refer to the bh_tsne output for further details")
                end
            else
                try
                    run(pipeline(`$BH_TSNE_path`, stdout=DevNull))
                catch excp
                    warn(excp)
                    error("ERROR: Call to bh_tsne exited with a non-zero return code exit status,
                          please enable verbose mode and refer to the bh_tsne output for further details")
                end
            end
        end

        open(joinpath(temp_dir, "result.dat"),"r") do output_file
            result_samples, result_dims = read_unpack("ii", output_file)
            results = [read_unpack(repeat("d", result_dims), output_file) for _ in 1:result_samples]
            results = [(read_unpack("i", output_file),e) for e in results]
            sort!(results)
            ret = Array{Float64,2}(result_dims, result_samples)
            for i in 1:result_samples
                ret[:,i] = collect(results[i][2])
            end
            return ret'
        end

    end

end

#=
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "-d","--no_dims"
            arg_type = Int64
        "-p","--perplexity"
            arg_type = Float64
            default = 50
        "-t","--theta"
            arg_type = Float64
            default = 0.5
        "-r","--randseed"
            arg_type = Int64
            default = -1
        "-n","--initial_dims"
            arg_type = Int64
            default = 50
        "-v","--verbose"
            action = :store_true
        "-i","--input"
            arg_type = ASCIIString
            default = STDIN
        "-o","--output"
            arg_type = ASCIIString
            default = STDOUT
    end

    parse_args(s)
end
=#

end


text:
#contains geomtry related functions
"""point in polygon
"""
function pointxy_inpoly_boundingbox(x,y,poly)
    poly_min=minimum(seg,dims=1)
    poly_max=maximum(seg,dims=1)
    ind=(x .> poly_max[1]) .|| (x .< poly_min[1]) .|| (y.>poly_max[2]) .|| (y.<poly_min[2])
end

"""compute the distance of xy from center and return boolean array for each radius limit
    radius is an array of radius which must be sorted

    return
"""
function point_in_circle!(x,y;x_cent::Number=0,y_cent::Number=0,radius1=0,radius2=0,radius3=0,
       dist2=nothing,sel_rad1=nothing,sel_rad2=nothing,sel_rad3=nothing)

    nobj=size(x,1)
    if(dist2==nothing)
        dist2=zeros(nobj)
    end

    @inbounds for ii in 1:nobj
        dist2[ii] = (x[ii]-x_cent)^2 + (y[ii]-y_cent)^2
    end

    if(radius1==0)
        return
    else
        @. sel_rad1 = dist2 < (radius1*radius1)
    end

    if(radius2==0)
        return
    else
        @. sel_rad2 = (!sel_rad1) & (dist2 < (radius2*radius2))
    end

    if(radius3==0)
        return
    else
        @. sel_rad3 = (!(sel_rad2 || sel_rad1)) & (dist2 < (radius3*radius3))
    end

end

"""
Estimates the winding number for point with respect to a given polygon
"""
function winding_count(pin,poly;dist=nothing)
    #Number of segments
    nseg=size(poly,1)-1

    #check if the point intersect y-axis
    ipos_y=@. poly[:,2]>pin[2]
    ineg_y=@. poly[:,2]<pin[2]

    is_up = @. ineg_y[begin:end-1] & ipos_y[begin+1:end]
    is_down = @. ipos_y[begin:end-1] & ineg_y[begin+1:end]

    #get is left
    a2=poly[begin:end-1,2] .- poly[begin+1:end,2]
    b2=poly[begin+1:end,1] .- poly[begin:end-1,1]
    c2=(poly[begin:end-1,1] .* poly[begin+1:end,2]) .- (poly[begin+1:end,1] .*poly[begin:end-1,2])

    if(dist==nothing)
        dist=zeros(nseg)
    end

    dist .=(a2 * pin[1]) + (b2 * pin[2]) + c2;
    is_left= @. dist<0
    is_right=@. dist>0

    wn = size(findall(is_up .& is_left),1) - size(findall(is_down .& is_right),1)
    return wn
    #@show wn

    #@show is_up,is_down

    if(false)
        plot(poly[:,1],poly[:,2],ls=:solid,color=:red)
        for ii in 1:nseg
            if(is_left[ii])
                plot!(poly[ii:ii+1,1],poly[ii:ii+1,2],ls=:solid,color=:blue)
            elseif(is_right[ii])
                plot!(poly[ii:ii+1,1],poly[ii:ii+1,2],ls=:solid,color=:green)
            end
        end
        plot!([pin[1]],[pin[2]],seriestype=:scatter,color=:black)
        plot!(poly[:,1],poly[:,2],seriestype=:scatter,color=:black) |> display
    end

end

"""distance of point p2 from line segment defined by p1,p2"""
function line_distance(p1,p2,p3)
    dist=(p1[2]-p2[2])*p3[1]
    dist +=  (p2[1]-p1[1])* p3[2]
    dist += ((p1[1]*p2[2]) - (p2[1]*p1[2]))
    return dist
end

""" Determines if a point is inside the polygon or not
Argument: pin: point in x,y plane of type array and length 2
 poly is array of polygon vertices with first and last vertex being same to draw all the edges

"""
function InPolyWinding_loop(pin,poly)
    #Number of segments
    nseg=size(poly,1)-1

    wn=0
    for ii in 1:nseg
        #println(ii)
        if((poly[ii,2]<pin[2]) & (poly[ii+1,2]>pin[2])) #is up
            dist=line_distance(view(poly,ii,:),view(poly,ii+1,:),pin)
            if(dist<0) # is left?
                wn +=1
            end
        elseif((poly[ii,2]>pin[2]) & (poly[ii+1,2]<pin[2])) # is down
            dist=line_distance(view(poly,ii,:),view(poly,ii+1,:),pin)
            if(dist>0) # is right?
                wn -=1
            end
        end
    end

    return wn!=0
end

"""
return the points in side the polygon and point inside the padded circle
The padded circle is defined with center at the mean of polygon vertices and
radius as the sum of maximum vertex length and a given padding (pad_circle)
The function must only set things to true and not set things to false
as it will be using arrays preallocated which might have been true for other reasons
Except when something is inpoly we want to et in circle to false
"""
function PointsInPoly_and_PaddedCircle!(x,y,poly,pad_circle;max_rad=nothing,in_poly=nothing,in_circle=nothing)
    nvert=size(poly,1)
    nobj=size(x,1)

    #memory allocation if not pre-allocated
    if(in_poly==nothing)
        in_poly=zeros(Bool,nobj)
    end
    if(in_circle==nothing)
        in_circle=zeros(Bool,nobj)
    end

    circ_cent=[mean(poly[:,1]),mean(poly[:,2])]

    #determine the maximum distance of polygon vertex from the center
    if(max_rad==nothing)
        max_rad2=maximum_radius2_polygon(poly,circ_cent)
        max_rad=sqrt(max_rad2)
    else
        max_rad2=max_rad^2
    end

    pad_circ_rad2= (max_rad+pad_circle)^2

    #calculate the distance of each point from the circle center
    @inbounds for ii in 1:nobj
        if(in_poly[ii])
            continue
        end
        #calculate the distance
        dist2 = (x[ii]-circ_cent[1])^2 + (y[ii]-circ_cent[2])^2

        if(dist2>pad_circ_rad2) #outside the padded circle
            continue
        elseif(dist2>max_rad2)#Determine if in padded circle
            in_circle[ii] = true
            continue
        elseif(InPolyWinding_loop([x[ii],y[ii]],poly)) #checks if truly inside the polygon
            in_poly[ii] = true
            #If something is inside the polygon then we must set in circle to false
            in_circle[ii] = false
        else
            in_circle[ii]=true
        end
    end

    #To make plots the result for testing and debugging
    if(false)
        plot(x[in_poly],y[in_poly],seriestype=:scatter,color=:red)
        plot!(x[in_circle],y[in_circle],seriestype=:scatter,color=:green)
        out_all= .! (in_circle .|| in_poly)
        #@show sum(out_all),sum(in_circle),sum(in_poly),size(in_poly)
        plot!(x[out_all],y[out_all],seriestype=:scatter,color=:grey)

        plot!([circ_cent[1]-sqrt(pad_circ_rad2),circ_cent[1]+sqrt(pad_circ_rad2)],
        [circ_cent[2],circ_cent[2]],s=:solid,color=:red,lw=2)

        plot!(poly[:,1],poly[:,2],ls=:solid,color=:black,lw=2) |> display

    end
end

"""Estimates the maximum radius of the polygon
"""
function maximum_radius2_polygon(poly,cent)
    max_rad2=0
    @inbounds for ii in 1:size(poly,1)-1
        this_rad2= (poly[ii,1]-cent[1])^2 +  (poly[ii,2]-cent[2])^2
        if(this_rad2>max_rad2)
            max_rad2=this_rad2
        end
    end
    return max_rad2
end

""" checks if two polygon intersects
First check for all the vertices in poly1 and check if any in poly2
Then rpeat to all vertices is in poly2 for poly1
If none of them are inside then polygon do not intersects
This function will not work for axis aligned polygon where they fully or
partially overlaps and the vetices of one is one the edges of another.
"""
function PolygonsIntersect(poly1,poly2)
    for ii in 1:size(poly1,1)
        if(InPolyWinding_loop(poly1[ii,:],poly2))
            return true
        end
    end

    for ii in 1:size(poly2,1)
        if(InPolyWinding_loop(poly2[ii,:],poly1))
            return true
        end
    end
    return false
end

"""Check if the bounding box of two input polygon is non-intersecting
just touhing is considered intersecting to be conservative

example:
poly1=zeros(5,2)
poly2=zeros(5,2)

poly1[:,1] .= [0,0,1,1,0]
poly1[:,2] .= [0,1,1,0,0]

dx=0
dy=0.8
poly2[:,1] .= [0,0,1,1,0] .+ dx
poly2[:,2] .= [0,1,1,0,0] .+ dy

PolygonsIntersect(poly1,poly2)
PolyIntersect_boundingbox(poly1,poly2)

plot(Shape(poly1[:,1],poly1[:,2]),color=:green)
plot!(Shape(poly2[:,1],poly2[:,2]),color=:red) |> display

"""
function PolyIntersect_boundingbox(poly1,poly2)
    poly1_min=minimum(poly1,dims=1)
    poly1_max=maximum(poly1,dims=1)

    poly2_min=minimum(poly2,dims=1)
    poly2_max=maximum(poly2,dims=1)

    if( ((poly1_max[1]<poly2_min[1]) || (poly2_max[1]<poly1_min[1])) ||
        ((poly1_max[2]<poly2_min[2]) || (poly2_max[2]<poly1_min[2]))
        )
        return false
    else
        return true
    end
end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
# TODO: fit(D::Distribution, x::Vector{Float}, w::Vector{Float}) --> weighted maximum likelihood estimation...

# Fitting function, add more methods here
function fit!(hmm, o; method=:baum_welch, max_iter=100, tol=1e-5, scaling=true)
	if method == :baum_welch
		return baum_welch!(hmm, o, max_iter, tol, scaling)
	end
end

function baum_welch!(hmm, o, max_iter, tol, scaling)
	# Fit hmm parameters given set of observation sequences
	# Setting tol to NaN will prevent early stopping, resulting in 'max_iter' iterations
	#n_seq = length(sequences)

    # convergence history of the fit, log-liklihood
    ch = (Float64)[]

    for k = 1:max_iter
    	push!(ch,0.0)

		# E step
		log_p_obs,alpha,beta,x,g = calc_stats(hmm,o,scaling)
		ch[end] += log_p_obs #(log_p_obs + log(length(o)))

		# M step (re-estimation)
		re_estimate!(hmm,o,x,g)

		if length(ch)>1 && (ch[end]-ch[end-1] < tol)
			println("Baum-Welch converged, stopping early")
			break
		end
	end

	return ch
end

function re_estimate!(hmm,o,x,g)
	# Update state transition matrix
    for i = 1:hmm.n
        denom = sum(g[1:end-1,i]) # Expected number of transitions from state 'i'
        for j = 1:hmm.n
            # Numerator is the expected number of transitions from state 'i' to 'j'
            hmm.A[i,j] = sum(x[:,i,j]) / denom
        end
    end

	# Update emission probability distributions
	for i = 1:hmm.n
		# weight each observation by probability of being in state i
		hmm.B[i] = fit_mle(typeof(hmm.B[i]), o, g[:,i])
	end

	# Re-estimate hmm.p (initial state probabilities)
	hmm.p = vec(g[1,:])
end

function calc_stats(hmm,o,scaling)
	## Single E-M iteration in the Baum-Welch procedure
	n_obs = length(o)

	# Calculate forward/backward probabilities
	alpha,beta,log_p_obs = forward_backward(hmm,o;scaling=scaling)

	# x[t,i,j] = probability of being in state 'i' at 't' and then in state 'j' at 't+1'
	x = zeros(n_obs-1, hmm.n, hmm.n)
	for t = 1:(n_obs-1)
		for i = 1:hmm.n
			for j = 1:hmm.n
				x[t,i,j] = alpha[t,i] * hmm.A[i,j] * pdf(hmm.B[j],o[t+1]) * beta[t+1,j]
			end
		end
		x[t,:,:] ./= sum(x[t,:,:]) # normalize to achieve probabilities
	end

	# g[t,i] = probability of being in state 'i' at step 't' given all observations
	g = alpha .* beta
	g ./= sum(g,2)   # normalize across states

	return log_p_obs, alpha, beta, x, g
end


text:
# Testcase from example given in Mocking.jl's README

import Compat: is_unix, read, unsafe_string

# Note: Function only works in UNIX environments.
function randdev(n::Integer)
    @mock open("/dev/urandom") do fp
        reverse(read(fp, n))
    end
end

n = 10
if is_unix()
    result = randdev(n)  # Reading /dev/urandom only works on UNIX environments
    @test eltype(result) == UInt8
    @test length(result) == n
end

# Produces a string with sequential UInt8 values from 1:n
data = unsafe_string(pointer(convert(Array{UInt8}, 1:n)))

# Generate a alternative method of `open` which call we wish to mock
patch = @patch open(fn::Function, f::AbstractString) = fn(IOBuffer(data))

# Apply the patch which will modify the behaviour for our test
apply(patch) do
    @test randdev(n) == convert(Array{UInt8}, n:-1:1)
end

if is_unix()
    # Outside of the scope of the patched environment `@mock` is essentially a no-op
    @test randdev(n) != convert(Array{UInt8}, n:-1:1)
end


text:
function getergodic(transmatrix::Array{Float64,2})
    D,V = eig(transmatrix')
    eigindex = indmin(abs.(D-1.0))
    ergodic = vec(V[:,eigindex]/sum(V[:,eigindex]))
    if any(abs.(imag(ergodic)).>1e-16)
        error("Imaginary ergodic distribution")
    else
        return real(ergodic)
    end
end

function createtransmatrix_sym(numregimes::Int)
    out = Array{SymPy.Sym}(numregimes)
    for j = 1:numregimes
        out[j] = Sym("TRANSITIONPROBABILITY_RX_RP$(j)")
    end
    return out
end

function TransitionMatrix(trans)
    # Sanity check the regime transition matrix
    numregimes = size(trans,1)
    if size(trans,2)!=size(trans,1)
        error("Regime transition matrix is not square.")
    end
    if any((abs.(sum(trans,2))-1.0).>1e-12)
        error("Rows of the regime transition matrix do not sum to 1.")
    end
    # Create a sympy version of the regime trans matrix
    syms = createtransmatrix_sym(numregimes)

    # Get the ergodic distribution of the markovian regime process
    ergodic = getergodic(trans)
    return TransitionMatrix(trans,syms,ergodic,Sym("TRANSITIONPROBABILITY"))
end


text:
using Polyhedra
using BenchmarkTools
using StaticArrays

# Intersection/Vector
# BenchmarkTools.Trial:
#   memory estimate:  133.73 KiB
#   allocs estimate:  2424
#   --------------
#   minimum time:     454.533 μs (0.00% GC)
#   median time:      461.111 μs (0.00% GC)
#   mean time:        488.875 μs (4.38% GC)
#   maximum time:     4.360 ms (87.57% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
# Intersection/SVector # FIXME why does it need more allocs ? It should be the opposite
# BenchmarkTools.Trial:
#   memory estimate:  572.52 KiB
#   allocs estimate:  21094
#   --------------
#   minimum time:     499.253 μs (0.00% GC)
#   median time:      612.736 μs (0.00% GC)
#   mean time:        866.356 μs (29.38% GC)
#   maximum time:     12.673 ms (93.11% GC)
#   --------------
#   samples:          5749
#   evals/sample:     1
# MixedMatHRep
# BenchmarkTools.Trial:
#   memory estimate:  148.64 KiB
#   allocs estimate:  2775
#   --------------
#   minimum time:     479.105 μs (0.00% GC)
#   median time:      486.236 μs (0.00% GC)
#   mean time:        508.813 μs (3.73% GC)
#   maximum time:     3.555 ms (83.56% GC)
#   --------------
#   samples:          9809
#   evals/sample:     1

let
    h = hrep([HalfSpace([ 1.,  1], 1),
              HalfSpace([ 1., -1], 0),
              HalfSpace([-1.,  0], 0)])
    display(@benchmark(doubledescription($h)))
end

let
    h = hrep([HalfSpace((@SVector [ 1,  1]), 1),
              HalfSpace((@SVector [ 1, -1]), 0),
              HalfSpace((@SVector [-1,  0]), 0)])
    display(@benchmark(doubledescription($h)))
end

let
    h = hrep([ 1  1
               1 -1
              -1  0],
             [1., 0, 0])
    display(@benchmark(doubledescription($h)))
end