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

using Core: CodeInfo, SimpleVector

const Callable = Union{Function,Type}

const Bottom = Union{}

"""
    AbstractSet{T}

Supertype for set-like types whose elements are of type `T`.
[`Set`](@ref), [`BitSet`](@ref) and other types are subtypes of this.
"""
abstract type AbstractSet{T} end

"""
    AbstractDict{K, V}

Supertype for dictionary-like types with keys of type `K` and values of type `V`.
[`Dict`](@ref), [`IdDict`](@ref) and other types are subtypes of this.
"""
abstract type AbstractDict{K,V} end

# The real @inline macro is not available until after array.jl, so this
# internal macro splices the meta Expr directly into the function body.
macro _inline_meta()
    Expr(:meta, :inline)
end
macro _noinline_meta()
    Expr(:meta, :noinline)
end

macro _gc_preserve_begin(arg1)
    Expr(:gc_preserve_begin, esc(arg1))
end

macro _gc_preserve_end(token)
    Expr(:gc_preserve_end, esc(token))
end

"""
    @nospecialize

Applied to a function argument name, hints to the compiler that the method
should not be specialized for different types of that argument,
but instead to use precisely the declared type for each argument.
This is only a hint for avoiding excess code generation.
Can be applied to an argument within a formal argument list,
or in the function body.
When applied to an argument, the macro must wrap the entire argument expression.
When used in a function body, the macro must occur in statement position and
before any code.

When used without arguments, it applies to all arguments of the parent scope.
In local scope, this means all arguments of the containing function.
In global (top-level) scope, this means all methods subsequently defined in the current module.

Specialization can reset back to the default by using [`@specialize`](@ref).

```julia
function example_function(@nospecialize x)
    ...
end

function example_function(@nospecialize(x = 1), y)
    ...
end

function example_function(x, y, z)
    @nospecialize x y
    ...
end

@nospecialize
f(y) = [x for x in y]
@specialize
```
"""
macro nospecialize(vars...)
    if nfields(vars) === 1
        # in argument position, need to fix `@nospecialize x=v` to `@nospecialize (kw x v)`
        var = getfield(vars, 1)
        if isa(var, Expr) && var.head === :(=)
            var.head = :kw
        end
    end
    return Expr(:meta, :nospecialize, vars...)
end

"""
    @specialize

Reset the specialization hint for an argument back to the default.
For details, see [`@nospecialize`](@ref).
"""
macro specialize(vars...)
    if nfields(vars) === 1
        # in argument position, need to fix `@specialize x=v` to `@specialize (kw x v)`
        var = getfield(vars, 1)
        if isa(var, Expr) && var.head === :(=)
            var.head = :kw
        end
    end
    return Expr(:meta, :specialize, vars...)
end

macro _pure_meta()
    return Expr(:meta, :pure)
end
# another version of inlining that propagates an inbounds context
macro _propagate_inbounds_meta()
    return Expr(:meta, :inline, :propagate_inbounds)
end

"""
    convert(T, x)

Convert `x` to a value of type `T`.

If `T` is an [`Integer`](@ref) type, an [`InexactError`](@ref) will be raised if `x`
is not representable by `T`, for example if `x` is not integer-valued, or is outside the
range supported by `T`.

# Examples
```jldoctest
julia> convert(Int, 3.0)
3

julia> convert(Int, 3.5)
ERROR: InexactError: Int64(3.5)
Stacktrace:
[...]
```

If `T` is a [`AbstractFloat`](@ref) or [`Rational`](@ref) type,
then it will return the closest value to `x` representable by `T`.

```jldoctest
julia> x = 1/3
0.3333333333333333

julia> convert(Float32, x)
0.33333334f0

julia> convert(Rational{Int32}, x)
1//3

julia> convert(Rational{Int64}, x)
6004799503160661//18014398509481984
```

If `T` is a collection type and `x` a collection, the result of
`convert(T, x)` may alias all or part of `x`.
```jldoctest
julia> x = Int[1, 2, 3];

julia> y = convert(Vector{Int}, x);

julia> y === x
true
```
"""
function convert end

convert(::Type{Any}, @nospecialize(x)) = x
convert(::Type{T}, x::T) where {T} = x
convert(::Type{Type}, x::Type) = x # the ssair optimizer is strongly dependent on this method existing to avoid over-specialization
                                   # in the absence of inlining-enabled
                                   # (due to fields typed as `Type`, which is generally a bad idea)

"""
    @eval [mod,] ex

Evaluate an expression with values interpolated into it using `eval`.
If two arguments are provided, the first is the module to evaluate in.
"""
macro eval(ex)
    :(Core.eval($__module__, $(Expr(:quote,ex))))
end
macro eval(mod, ex)
    :(Core.eval($(esc(mod)), $(Expr(:quote,ex))))
end

argtail(x, rest...) = rest

"""
    tail(x::Tuple)::Tuple

Return a `Tuple` consisting of all but the first component of `x`.

# Examples
```jldoctest
julia> Base.tail((1,2,3))
(2, 3)

julia> Base.tail(())
ERROR: ArgumentError: Cannot call tail on an empty tuple.
```
"""
tail(x::Tuple) = argtail(x...)
tail(::Tuple{}) = throw(ArgumentError("Cannot call tail on an empty tuple."))

tuple_type_head(T::Type) = (@_pure_meta; fieldtype(T::Type{<:Tuple}, 1))

function tuple_type_tail(T::Type)
    @_pure_meta
    if isa(T, UnionAll)
        return UnionAll(T.var, tuple_type_tail(T.body))
    elseif isa(T, Union)
        return Union{tuple_type_tail(T.a), tuple_type_tail(T.b)}
    else
        T.name === Tuple.name || throw(MethodError(tuple_type_tail, (T,)))
        if isvatuple(T) && length(T.parameters) == 1
            va = T.parameters[1]
            (isa(va, DataType) && isa(va.parameters[2], Int)) || return T
            return Tuple{Vararg{va.parameters[1], va.parameters[2]-1}}
        end
        return Tuple{argtail(T.parameters...)...}
    end
end

tuple_type_cons(::Type, ::Type{Union{}}) = Union{}
function tuple_type_cons(::Type{S}, ::Type{T}) where T<:Tuple where S
    @_pure_meta
    Tuple{S, T.parameters...}
end

function unwrap_unionall(@nospecialize(a))
    while isa(a,UnionAll)
        a = a.body
    end
    return a
end

function rewrap_unionall(@nospecialize(t), @nospecialize(u))
    if !isa(u, UnionAll)
        return t
    end
    return UnionAll(u.var, rewrap_unionall(t, u.body))
end

# replace TypeVars in all enclosing UnionAlls with fresh TypeVars
function rename_unionall(@nospecialize(u))
    if !isa(u,UnionAll)
        return u
    end
    body = rename_unionall(u.body)
    if body === u.body
        body = u
    else
        body = UnionAll(u.var, body)
    end
    var = u.var::TypeVar
    nv = TypeVar(var.name, var.lb, var.ub)
    return UnionAll(nv, body{nv})
end

const _va_typename = Vararg.body.body.name
function isvarargtype(@nospecialize(t))
    t = unwrap_unionall(t)
    return isa(t, DataType) && (t::DataType).name === _va_typename
end

function isvatuple(@nospecialize(t))
    t = unwrap_unionall(t)
    if isa(t, DataType)
        n = length(t.parameters)
        return n > 0 && isvarargtype(t.parameters[n])
    end
    return false
end

function unwrapva(@nospecialize(t))
    # NOTE: this returns a related type, but it's NOT a subtype of the original tuple
    t2 = unwrap_unionall(t)
    return isvarargtype(t2) ? rewrap_unionall(t2.parameters[1], t) : t
end

function unconstrain_vararg_length(@nospecialize(va))
    # construct a new Vararg type where its length is unconstrained,
    # but its element type still captures any dependencies the input
    # element type may have had on the input length
    T = unwrap_unionall(va).parameters[1]
    return rewrap_unionall(Vararg{T}, va)
end

typename(a) = error("typename does not apply to this type")
typename(a::DataType) = a.name
function typename(a::Union)
    ta = typename(a.a)
    tb = typename(a.b)
    ta === tb || error("typename does not apply to unions whose components have different typenames")
    return tb
end
typename(union::UnionAll) = typename(union.body)

const AtLeast1 = Tuple{Any, Vararg{Any}}

# converting to empty tuple type
convert(::Type{Tuple{}}, ::Tuple{}) = ()
convert(::Type{Tuple{}}, x::AtLeast1) = throw(MethodError(convert, (Tuple{}, x)))

# converting to tuple types with at least one element
convert(::Type{T}, x::T) where {T<:AtLeast1} = x
convert(::Type{T}, x::AtLeast1) where {T<:AtLeast1} =
    (convert(tuple_type_head(T), x[1]), convert(tuple_type_tail(T), tail(x))...)

# converting to Vararg tuple types
convert(::Type{Tuple{Vararg{V}}}, x::Tuple{Vararg{V}}) where {V} = x
convert(T::Type{Tuple{Vararg{V}}}, x::Tuple) where {V} =
    (convert(tuple_type_head(T), x[1]), convert(T, tail(x))...)

# TODO: the following definitions are equivalent (behaviorally) to the above method
# I think they may be faster / more efficient for inference,
# if we could enable them, but are they?
# TODO: These currently can't be used (#21026, #23017) since with
#     z(::Type{<:Tuple{Vararg{T}}}) where {T} = T
#   calling
#     z(Tuple{Val{T}} where T)
#   fails, even though `Type{Tuple{Val}} == Type{Tuple{Val{S}} where S}`
#   and so T should be `Val` (aka `Val{S} where S`)
#convert(_::Type{Tuple{S}}, x::Tuple{S}) where {S} = x
#convert(_::Type{Tuple{S}}, x::Tuple{Any}) where {S} = (convert(S, x[1]),)
#convert(_::Type{T}, x::T) where {S, N, T<:Tuple{S, Vararg{S, N}}} = x
#convert(_::Type{Tuple{S, Vararg{S, N}}},
#        x::Tuple{Any, Vararg{Any, N}}) where
#       {S, N} = cnvt_all(S, x...)
#convert(_::Type{Tuple{Vararg{S, N}}},
#        x::Tuple{Vararg{Any, N}}) where
#       {S, N} = cnvt_all(S, x...)
# TODO: These currently can't be used since
#   Type{NTuple} <: (Type{Tuple{Vararg{S}}} where S) is true
#   even though the value S doesn't exist
#convert(_::Type{Tuple{Vararg{S}}},
#        x::Tuple{Any, Vararg{Any}}) where
#       {S} = cnvt_all(S, x...)
#convert(_::Type{Tuple{Vararg{S}}},
#        x::Tuple{Vararg{Any}}) where
#       {S} = cnvt_all(S, x...)
#cnvt_all(T) = ()
#cnvt_all(T, x, rest...) = (convert(T, x), cnvt_all(T, rest...)...)
# TODO: These may be necessary if the above are enabled
#convert(::Type{Tuple{}}, ::Tuple{}) = ()
#convert(::Type{Tuple{Vararg{S}}} where S, ::Tuple{}) = ()

"""
    oftype(x, y)

Convert `y` to the type of `x` (`convert(typeof(x), y)`).

# Examples
```jldoctest
julia> x = 4;

julia> y = 3.;

julia> oftype(x, y)
3

julia> oftype(y, x)
4.0
```
"""
oftype(x, y) = convert(typeof(x), y)

unsigned(x::Int) = reinterpret(UInt, x)
signed(x::UInt) = reinterpret(Int, x)

"""
    cconvert(T,x)

Convert `x` to a value to be passed to C code as type `T`, typically by calling `convert(T, x)`.

In cases where `x` cannot be safely converted to `T`, unlike [`convert`](@ref), `cconvert` may
return an object of a type different from `T`, which however is suitable for
[`unsafe_convert`](@ref) to handle. The result of this function should be kept valid (for the GC)
until the result of [`unsafe_convert`](@ref) is not needed anymore.
This can be used to allocate memory that will be accessed by the `ccall`.
If multiple objects need to be allocated, a tuple of the objects can be used as return value.

Neither `convert` nor `cconvert` should take a Julia object and turn it into a `Ptr`.
"""
function cconvert end

cconvert(T::Type, x) = convert(T, x) # do the conversion eagerly in most cases
cconvert(::Type{<:Ptr}, x) = x # but defer the conversion to Ptr to unsafe_convert
unsafe_convert(::Type{T}, x::T) where {T} = x # unsafe_convert (like convert) defaults to assuming the convert occurred
unsafe_convert(::Type{T}, x::T) where {T<:Ptr} = x  # to resolve ambiguity with the next method
unsafe_convert(::Type{P}, x::Ptr) where {P<:Ptr} = convert(P, x)

"""
    reinterpret(type, A)

Change the type-interpretation of a block of memory.
For arrays, this constructs a view of the array with the same binary data as the given
array, but with the specified element type.
For example,
`reinterpret(Float32, UInt32(7))` interprets the 4 bytes corresponding to `UInt32(7)` as a
[`Float32`](@ref).

# Examples
```jldoctest
julia> reinterpret(Float32, UInt32(7))
1.0f-44

julia> reinterpret(Float32, UInt32[1 2 3 4 5])
1×5 reinterpret(Float32, ::Array{UInt32,2}):
 1.4013e-45  2.8026e-45  4.2039e-45  5.60519e-45  7.00649e-45
```
"""
reinterpret(::Type{T}, x) where {T} = bitcast(T, x)
reinterpret(::Type{Unsigned}, x::Float16) = reinterpret(UInt16,x)
reinterpret(::Type{Signed}, x::Float16) = reinterpret(Int16,x)

"""
    sizeof(T::DataType)
    sizeof(obj)

Size, in bytes, of the canonical binary representation of the given `DataType` `T`, if any.
Size, in bytes, of object `obj` if it is not `DataType`.

# Examples
```jldoctest
julia> sizeof(Float32)
4

julia> sizeof(ComplexF64)
16

julia> sizeof(1.0)
8

julia> sizeof([1.0:10.0;])
80
```

If `DataType` `T` does not have a specific size, an error is thrown.

```jldoctest
julia> sizeof(AbstractArray)
ERROR: Abstract type AbstractArray does not have a definite size.
Stacktrace:
[...]
```
"""
sizeof(x) = Core.sizeof(x)

# simple Array{Any} operations needed for bootstrap
@eval setindex!(A::Array{Any}, @nospecialize(x), i::Int) = arrayset($(Expr(:boundscheck)), A, x, i)

"""
    precompile(f, args::Tuple{Vararg{Any}})

Compile the given function `f` for the argument tuple (of types) `args`, but do not execute it.
"""
function precompile(@nospecialize(f), args::Tuple)
    ccall(:jl_compile_hint, Int32, (Any,), Tuple{Core.Typeof(f), args...}) != 0
end

function precompile(argt::Type)
    ccall(:jl_compile_hint, Int32, (Any,), argt) != 0
end

"""
    esc(e)

Only valid in the context of an [`Expr`](@ref) returned from a macro. Prevents the macro hygiene
pass from turning embedded variables into gensym variables. See the [Macros](@ref man-macros)
section of the Metaprogramming chapter of the manual for more details and examples.
"""
esc(@nospecialize(e)) = Expr(:escape, e)

"""
    @boundscheck(blk)

Annotates the expression `blk` as a bounds checking block, allowing it to be elided by [`@inbounds`](@ref).

!!! note
    The function in which `@boundscheck` is written must be inlined into
    its caller in order for `@inbounds` to have effect.

# Examples
```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
julia> @inline function g(A, i)
           @boundscheck checkbounds(A, i)
           return "accessing (\$A)[\$i]"
       end;

julia> f1() = return g(1:2, -1);

julia> f2() = @inbounds return g(1:2, -1);

julia> f1()
ERROR: BoundsError: attempt to access 2-element UnitRange{Int64} at index [-1]
Stacktrace:
 [1] throw_boundserror(::UnitRange{Int64}, ::Tuple{Int64}) at ./abstractarray.jl:455
 [2] checkbounds at ./abstractarray.jl:420 [inlined]
 [3] g at ./none:2 [inlined]
 [4] f1() at ./none:1
 [5] top-level scope

julia> f2()
"accessing (1:2)[-1]"
```

!!! warning

    The `@boundscheck` annotation allows you, as a library writer, to opt-in to
    allowing *other code* to remove your bounds checks with [`@inbounds`](@ref).
    As noted there, the caller must verify—using information they can access—that
    their accesses are valid before using `@inbounds`. For indexing into your
    [`AbstractArray`](@ref) subclasses, for example, this involves checking the
    indices against its [`size`](@ref). Therefore, `@boundscheck` annotations
    should only be added to a [`getindex`](@ref) or [`setindex!`](@ref)
    implementation after you are certain its behavior is correct.
"""
macro boundscheck(blk)
    return Expr(:if, Expr(:boundscheck), esc(blk))
end

"""
    @inbounds(blk)

Eliminates array bounds checking within expressions.

In the example below the in-range check for referencing
element `i` of array `A` is skipped to improve performance.

```julia
function sum(A::AbstractArray)
    r = zero(eltype(A))
    for i = 1:length(A)
        @inbounds r += A[i]
    end
    return r
end
```

!!! warning

    Using `@inbounds` may return incorrect results/crashes/corruption
    for out-of-bounds indices. The user is responsible for checking it manually.
    Only use `@inbounds` when it is certain from the information locally available
    that all accesses are in bounds.
"""
macro inbounds(blk)
    return Expr(:block,
        Expr(:inbounds, true),
        Expr(:local, Expr(:(=), :val, esc(blk))),
        Expr(:inbounds, :pop),
        :val)
end

"""
    @label name

Labels a statement with the symbolic label `name`. The label marks the end-point
of an unconditional jump with [`@goto name`](@ref).
"""
macro label(name::Symbol)
    return esc(Expr(:symboliclabel, name))
end

"""
    @goto name

`@goto name` unconditionally jumps to the statement at the location [`@label name`](@ref).

`@label` and `@goto` cannot create jumps to different top-level statements. Attempts cause an
error. To still use `@goto`, enclose the `@label` and `@goto` in a block.
"""
macro goto(name::Symbol)
    return esc(Expr(:symbolicgoto, name))
end

# SimpleVector

function getindex(v::SimpleVector, i::Int)
    @boundscheck if !(1 <= i <= length(v))
        throw(BoundsError(v,i))
    end
    t = @_gc_preserve_begin v
    x = unsafe_load(convert(Ptr{Ptr{Cvoid}},pointer_from_objref(v)) + i*sizeof(Ptr))
    x == C_NULL && throw(UndefRefError())
    o = unsafe_pointer_to_objref(x)
    @_gc_preserve_end t
    return o
end

function length(v::SimpleVector)
    t = @_gc_preserve_begin v
    l = unsafe_load(convert(Ptr{Int},pointer_from_objref(v)))
    @_gc_preserve_end t
    return l
end
firstindex(v::SimpleVector) = 1
lastindex(v::SimpleVector) = length(v)
iterate(v::SimpleVector, i=1) = (length(v) < i ? nothing : (v[i], i + 1))
eltype(::Type{SimpleVector}) = Any
keys(v::SimpleVector) = OneTo(length(v))
isempty(v::SimpleVector) = (length(v) == 0)
axes(v::SimpleVector) = (OneTo(length(v)),)
axes(v::SimpleVector, d::Integer) = d <= 1 ? axes(v)[d] : OneTo(1)

function ==(v1::SimpleVector, v2::SimpleVector)
    length(v1)==length(v2) || return false
    for i = 1:length(v1)
        v1[i] == v2[i] || return false
    end
    return true
end

map(f, v::SimpleVector) = Any[ f(v[i]) for i = 1:length(v) ]

getindex(v::SimpleVector, I::AbstractArray) = Core.svec(Any[ v[i] for i in I ]...)

"""
    isassigned(array, i) -> Bool

Test whether the given array has a value associated with index `i`. Return `false`
if the index is out of bounds, or has an undefined reference.

# Examples
```jldoctest
julia> isassigned(rand(3, 3), 5)
true

julia> isassigned(rand(3, 3), 3 * 3 + 1)
false

julia> mutable struct Foo end

julia> v = similar(rand(3), Foo)
3-element Array{Foo,1}:
 #undef
 #undef
 #undef

julia> isassigned(v, 1)
false
```
"""
function isassigned end

function isassigned(v::SimpleVector, i::Int)
    @boundscheck 1 <= i <= length(v) || return false
    t = @_gc_preserve_begin v
    x = unsafe_load(convert(Ptr{Ptr{Cvoid}},pointer_from_objref(v)) + i*sizeof(Ptr))
    @_gc_preserve_end t
    return x != C_NULL
end

"""
    Colon()

Colons (:) are used to signify indexing entire objects or dimensions at once.

Very few operations are defined on Colons directly; instead they are converted
by [`to_indices`](@ref) to an internal vector type (`Base.Slice`) to represent the
collection of indices they span before being used.

The singleton instance of `Colon` is also a function used to construct ranges;
see [`:`](@ref).
"""
struct Colon <: Function
end
const (:) = Colon()

"""
    Val(c)

Return `Val{c}()`, which contains no run-time data. Types like this can be used to
pass the information between functions through the value `c`, which must be an `isbits`
value. The intent of this construct is to be able to dispatch on constants directly (at
compile time) without having to test the value of the constant at run time.

# Examples
```jldoctest
julia> f(::Val{true}) = "Good"
f (generic function with 1 method)

julia> f(::Val{false}) = "Bad"
f (generic function with 2 methods)

julia> f(Val(true))
"Good"
```
"""
struct Val{x}
end

Val(x) = (@_pure_meta; Val{x}())

"""
    invokelatest(f, args...; kwargs...)

Calls `f(args...; kwargs...)`, but guarantees that the most recent method of `f`
will be executed.   This is useful in specialized circumstances,
e.g. long-running event loops or callback functions that may
call obsolete versions of a function `f`.
(The drawback is that `invokelatest` is somewhat slower than calling
`f` directly, and the type of the result cannot be inferred by the compiler.)
"""
function invokelatest(@nospecialize(f), @nospecialize args...; kwargs...)
    if isempty(kwargs)
        return Core._apply_latest(f, args)
    end
    # We use a closure (`inner`) to handle kwargs.
    inner() = f(args...; kwargs...)
    Core._apply_latest(inner)
end

# TODO: possibly make this an intrinsic
inferencebarrier(@nospecialize(x)) = Ref{Any}(x)[]

"""
    isempty(collection) -> Bool

Determine whether a collection is empty (has no elements).

# Examples
```jldoctest
julia> isempty([])
true

julia> isempty([1 2 3])
false
```
"""
function isempty(itr)
    d = isdone(itr)
    d !== missing && return d
    iterate(itr) === nothing
end

"""
    values(iterator)

For an iterator or collection that has keys and values, return an iterator
over the values.
This function simply returns its argument by default, since the elements
of a general iterator are normally considered its "values".

# Examples
```jldoctest
julia> d = Dict("a"=>1, "b"=>2);

julia> values(d)
Base.ValueIterator for a Dict{String,Int64} with 2 entries. Values:
  2
  1

julia> values([2])
1-element Array{Int64,1}:
 2
```
"""
values(itr) = itr

"""
    Missing

A type with no fields whose singleton instance [`missing`](@ref) is used
to represent missing values.
"""
struct Missing end

"""
    missing

The singleton instance of type [`Missing`](@ref) representing a missing value.
"""
const missing = Missing()

"""
    ismissing(x)

Indicate whether `x` is [`missing`](@ref).
"""
ismissing(::Any) = false
ismissing(::Missing) = true

function popfirst! end
function peek end

"""
    @__LINE__ -> Int

Expand to the line number of the location of the macrocall.
Return `0` if the line number could not be determined.
"""
macro __LINE__()
    return __source__.line
end

# Just for bootstrapping purposes below
macro __FILE_SYMBOL__()
    return Expr(:quote, __source__.file)
end

# Iteration
"""
    isdone(itr, state...) -> Union{Bool, Missing}

This function provides a fast-path hint for iterator completion.
This is useful for mutable iterators that want to avoid having elements
consumed, if they are not going to be exposed to the user (e.g. to check
for done-ness in `isempty` or `zip`). Mutable iterators that want to
opt into this feature should define an isdone method that returns
true/false depending on whether the iterator is done or not. Stateless
iterators need not implement this function. If the result is `missing`,
callers may go ahead and compute `iterate(x, state...) === nothing` to
compute a definite answer.
"""
isdone(itr, state...) = missing

"""
    iterate(iter [, state]) -> Union{Nothing, Tuple{Any, Any}}

Advance the iterator to obtain the next element. If no elements
remain, `nothing` should be returned. Otherwise, a 2-tuple of the
next element and the new iteration state should be returned.
"""
function iterate end

"""
    isiterable(T) -> Bool

Test if type `T` is an iterable collection type or not,
that is whether it has an `iterate` method or not.
"""
function isiterable(T)::Bool
    return hasmethod(iterate, Tuple{T})
end

max_stars_repo_name: ajozefiak/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using .GraphSearch: TreeSearchNode
import Base: push!, length, keys, pop!, peek, get, delete!, haskey, values, empty!

mutable struct MultipleOpenLists{S,A}
    key_to_nodehandle_map::Dict{S, Tuple{TreeSearchNode{S,A}, Float64, Float64, Int}}
    weights::Vector{Float64}
    pqs::Vector{MutableBinaryMinHeap{Tuple{Float64, Int, S}}}
    stats::Tuple{Float64,Float64,Float64,Float64,Float64,Int}
    counter::Int
    function MultipleOpenLists{S,A}(weights::AbstractArray{T}) where {S,A,T<:Real}
        mol = new{S,A}()
        mol.key_to_nodehandle_map = Dict{S, Tuple{TreeSearchNode{S,A}, Float64, Float64, Int}}()
        mol.weights =  convert(Vector{Float64}, weights)
        if !(1.0 in mol.weights)
            mol.weights = vcat(1.0, mol.weights)
        end
        mol.pqs = MutableBinaryMinHeap{Tuple{Float64, Int, S}}[MutableBinaryMinHeap{Tuple{Float64, Int, S}}() for w in mol.weights]
        mol.stats = (0.0, 0.0, 0.0, 0.0, 0.0, 0)
        mol.counter = 0
        return mol
    end
end

length(mol::MultipleOpenLists) = length(mol.key_to_nodehandle_map)
keys(mol::MultipleOpenLists) = keys(mol.key_to_nodehandle_map)
values(mol::MultipleOpenLists) = values(mol.key_to_nodehandle_map)
haskey(mol::MultipleOpenLists{S,A}, key::S) where {S,A} = key in keys(mol.key_to_nodehandle_map)

function push!(mol::MultipleOpenLists{S,A}, key::S, node::TreeSearchNode{S,A}, g::Float64, h::Float64)::Nothing where {S,A}
    if haskey(mol, key)
        error("already present")
    end
    mol.counter -= 1  # depth first tie breaking
    handle::Int = 0
    @inbounds for i in 1:length(mol.weights)
        w = mol.weights[i]
        f_w = g + w * h
        handle = push!(mol.pqs[i], (f_w, mol.counter, key))
    end
    entry = (node, g, h, handle)
    mol.key_to_nodehandle_map[key] = entry
    mol.stats = update_mean_std_corr(mol.stats..., g, h)
    return nothing
end

@inline function idxof(weight::Float64, weights::AbstractArray{Float64})::Int
    findfirst(isequal(weight), weights)
end

function pop!(mol::MultipleOpenLists{S,A}, weight::Float64)::Tuple{TreeSearchNode{S,A}, Float64, Float64} where {S,A}
    @inbounds pq = mol.pqs[idxof(weight, mol.weights)]
    (_, _, key), handle = top_with_handle(pq)
    pop!(pq)
    @inbounds for i in 1:length(mol.weights)
        w = mol.weights[i]
        w ≠ weight && delete!(mol.pqs[i], handle)
    end
    node, g, h, handle = mol.key_to_nodehandle_map[key]
    delete!(mol.key_to_nodehandle_map, key)
    mol.stats = update_mean_std_corr(mol.stats..., g, h, remove=true)
    return node, g, h
end

function peek(mol::MultipleOpenLists{S,A}, w::Float64)::Tuple{TreeSearchNode{S,A}, Float64, Float64} where {S,A}
    @inbounds pq = mol.pqs[idxof(w, mol.weights)]
    (_, _, key), handle = top_with_handle(pq)
    node, g, h, handle = mol.key_to_nodehandle_map[key]
    return node, g, h
end

function get(mol::MultipleOpenLists{S,A}, key::S)::Tuple{TreeSearchNode{S,A}, Float64, Float64} where {S,A}
    return mol.key_to_nodehandle_map[key]
end

function delete!(mol::MultipleOpenLists{S,A}, key::S)::Nothing where {S,A}
    _, g, h, handle = mol.key_to_nodehandle_map[key]
    @inbounds for i in 1:length(mol.weights)
        delete!(mol.pqs[i], handle)
    end
    delete!(mol.key_to_nodehandle_map, key)
    mol.stats = update_mean_std_corr(mol.stats..., g, h, remove=true)
    return nothing
end

function empty!(mol::MultipleOpenLists{S, A})::Nothing where {S, A}
    empty!(mol.key_to_nodehandle_map)
    @inbounds for i in 1:length(mol.weights)
        mol.pqs[i] = MutableBinaryMinHeap{Tuple{Float64, Int, S}}()  # TODO: is there an equivalent of empty! ?
    end
    mol.stats = (0.0, 0.0, 0.0, 0.0, 0.0, 0)
    mol.counter = 0
    return nothing
end

ext: jl


text:
function density(obj::MultivariateRealDistribution, arg0::Vector{jdouble})
    return jcall(obj, "density", jdouble, (Vector{jdouble},), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_dimension(obj::AbstractMultivariateRealDistribution)
    return jcall(obj, "getDimension", jint, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function reseed_random_generator(obj::AbstractMultivariateRealDistribution, arg0::jlong)
    return jcall(obj, "reseedRandomGenerator", void, (jlong,), arg0)
end

function sample(obj::AbstractMultivariateRealDistribution)
    return jcall(obj, "sample", Vector{jdouble}, ())
end

function sample(obj::AbstractMultivariateRealDistribution, arg0::jint)
    return jcall(obj, "sample", Vector{Vector{jdouble}}, (jint,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using PyPlot

xkcd()
fig = figure(figsize=(10,5))
ax = axes()

x = linspace(1,10,101)
y = sin(3x + cos(5x))
p = plot(x,y)

ax[:set_xlim]([0.0,6])
annotate("A little\nsketchy",xy=[0.98,.001],
arrowprops=["arrowstyle"=>"->"],xytext=[1.3,-0.3])

xticks([])
yticks([])
xlabel("TIME")
ylabel("PRODUCTIVITY")
title("An xkcd-style plot in Julia")

ax[:spines]["top"][:set_color]("none")
ax[:spines]["right"][:set_color]("none")

max_stars_repo_path: Module 3/Chapter07/xkcd-graph.jl

max_issues_repo_path: Module 3/Chapter07/xkcd-graph.jl

max_forks_repo_path: Module 3/Chapter07/xkcd-graph.jl


text:
#-----------------------------------------------------------------------
#   MULTIPLE EMISSION SET VECTOR GENERATOR
#-----------------------------------------------------------------------
"""
    multivec( X::RealMtx , k::Int )

This function takes a ``m×n`` matrix `X`, and for `k```<n``,
constructs a vector of ``m×```k` matrices constructed from all choices
of `k` columns from the ``m×n`` matrix `X`.

"""
function multivec( X::RealMtx , k::Int )
    tpfl = typeof(X[1,1])
    d    = size(X)[1]
    np   = size(X)[2]
    Z    = [zeros(tpfl,d) for _ in 1:np]

    for i=1:np
        Z[i] = X[:,i]
    end

    w   = collect(combinations(Z,k))

    l   = length(w)
    z   = [zeros(tpfl,d,k) for _ in 1:l]
    for i=1:l, j=1:k
        z[i][:,j] = w[i][j]
    end

    return z
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   MULTIPLE EMISSION SET LOCATOR
#-----------------------------------------------------------------------
"""
    mlocator( X::RealMtx )

The function `mlocator` computes a single intersection point ``X_c``
from an ``m×n`` matrix `X` formed from `n>4` emission points, which form
the column vectors of `X`.

"""
function mlocator( X::RealMtx )
    #   Calculates location for more than five emission points
    tpfl = typeof(X[1,1])

    l = size(X)

    if l[2] == 5
        return locator5RTC21(X)
    elseif l[2] > 5
        W   = multivec(X,5)
        k   = length(W)

        w   = [zeros(tpfl,4) for _ =1:2*k ]

        for a=1:k
            w[a] = locator5RTC21(W[a])
        end

        #   Re-sort location points according to smallest Minkowski norm
        δW = zeros(tpfl,k)
        for a=1:k
            for i=1:4
                δW[a] += abs(ηdot(X[:,i]-w[a],X[:,i]-w[a]))
            end
        end
        iW = sortperm(δW)

        return w[iW][1]
    else
        return zeros(tpfl,4)
    end
end     #---------------------------------------------------------------


text:
AbsoluteRectangle(mini::Vec{N,T}, maxi::Vec{N,T}) where {N,T} = HyperRectangle{N,T}(mini, maxi-mini)

AABB(a) = AABB{Float32}(a)
function (B::Type{AABB{T}})(a::Pyramid) where T
    w,h = a.width/T(2), a.length
    m = Vec{3,T}(a.middle)
    B(m-Vec{3,T}(w,w,0), m+Vec{3,T}(w, w, h))
end
(B::Type{AABB{T}})(a::Cube) where {T} = B(origin(a), widths(a))
(B::Type{AABB{T}})(a::AbstractMesh) where {T} = B(vertices(a))
(B::Type{AABB{T}})(a::NativeMesh) where {T} = B(gpu_data(a.data[:vertices]))

function (B::Type{AABB{T}})(
        positions, scale, rotation,
        primitive::AABB{T}
    ) where T

    ti = TransformationIterator(positions, scale, rotation)
    B(ti, primitive)
end
function (B::Type{AABB{T}})(instances::Instances) where T
    ti = TransformationIterator(instances)
    B(ti, B(instances.primitive))
end

function transform(translation, scale, rotation, points)
    _max = Vec3f0(typemin(Float32))
    _min = Vec3f0(typemax(Float32))
    for p in points
        x = scale.*Vec3f0(p)
        rv = rotation*Vec4f0(x[1], x[2], x[3], 1f0)
        x = Vec3f0(rv[1], rv[2], rv[3])
        x = Vec3f0(translation)+x
        _min = min.(_min, x)
        _max = max.(_max, x)
    end
    AABB{Float32}(_min, _max-_min)
end

function (B::Type{AABB{T}})(
      ti::TransformationIterator, primitive::AABB{T}
    ) where T
    v_state = iterate(ti)
    v_state === nothing && return primitive

    tsr, state = v_state
    points = decompose(Point3f0, primitive)::Vector{Point3f0}
    bb = transform(tsr[1], tsr[2], tsr[3], points)
    v_state = iterate(ti, state)
    while v_state !== nothing
        tsr, state = v_state
        translatet_bb = transform(tsr[1], tsr[2], tsr[3], points)
        bb = union(bb, translatet_bb)
        v_state = iterate(ti, state)
    end
    bb
end


text:
# Creation macros #
# ========================== #
"""
    @scenario(def)

Define a scenario type compatible with StochasticPrograms using the syntax
```julia
@scenario name = begin
    ...structdef...

    [@zero begin
        ...
        return zero(scenario)
    end]

    [@expectation begin
        ...
        return expected(scenarios)
     end]
end
```
The generated type is referenced through `name` and a default constructor is always generated. This constructor accepts the keyword `probability` to set the probability of the scenario occuring. Otherwise, any internal variables and specialized constructors are defined in the @scenario block as they would be in any Julia struct.

If possible, a `zero` method and an `expected` method will be generated for the defined type. Otherwise, or if the default implementation is not desired, these can be user provided through [`@zero`](@ref) and [`@expectation`](@ref).

The defined scenario type will be available on all Julia processes.

## Examples

The following defines a simple scenario ``ξ`` with a single value.

```jldoctest
@scenario ExampleScenario = begin
    ξ::Float64
end

ExampleScenario(1.0, probability = 0.5)

# output

ExampleScenario with probability 0.5
  ξ: 1.0

```

See also: [`@zero`](@ref), [`@expectation`](@ref), [`@sampler`](@ref)
"""
macro scenario(arg)
    @capture(arg, scenarioname_Symbol = scenariodef_) || error("Invalid syntax. Expected: scenarioname = begin scenariodef end")
    vars = Vector{Symbol}()
    vartypes = Vector{Union{Expr,Symbol}}()
    vardefs = Vector{Expr}()
    zerodefs = Vector{Expr}()
    expectdefs = Vector{Expr}()
    def = postwalk(prettify(scenariodef)) do x
        if @capture(x, var_Symbol::vartype_)
            push!(vars, var)
            push!(vartypes, vartype)
            push!(vardefs, x)
            return x
        end
        if @capture(x, @zero zerodef_)
            push!(zerodefs, zerodef)
            return @q begin end
        end
        if @capture(x, @expectation expectdef_)
            push!(expectdefs, expectdef)
            return @q begin end
        end
        return x
    end
    # Handle zero definition
    if length(zerodefs) > 1
        error("Only provide one zero implementation")
    end
    provided_zerodef, zerodef = if length(zerodefs) == 1
        zerodef = zerodefs[1]
        true, zerodef
    else
        false, Expr(:block)
    end
    # Handle expectation definition
    if length(expectdefs) > 1
        error("Only provide one expectation implementation")
    end
    provided_expectdef, expectdef = if length(expectdefs) == 1
        expectdef = expectdefs[1]
        expectdef = postwalk(prettify(expectdef)) do x
            @capture(x, constructor_Symbol(args__)) && constructor == scenarioname && return :(StochasticPrograms.ExpectedScenario($x))
            return x
        end
        true, expectdef
    else
        false, Expr(:block)
    end
    # Prepare automatic expectation definition
    combine = Expr(:call)
    push!(combine.args, scenarioname)
    push!(combine.args, :($(Expr(:parameters, :(probability = 1.0)))))
    for var in vars
        push!(combine.args, :(probability(s1)*s1.$var + probability(s2)*s2.$var))
    end
    # Prepare automatic scenariotext definition
    textdef = Expr(:block)
    for var in vars
        key = Meta.quot(var)
        push!(textdef.args, :(print(io, "\n  $($key): $(scenario.$var)")))
    end
    # Define scenario type
    code = @q begin
        if StochasticPrograms.supports_expected([$(vartypes...)], $provided_expectdef) || $provided_expectdef
            if StochasticPrograms.supports_zero([$(vartypes...)], $provided_zerodef) || $provided_zerodef
                struct $scenarioname <: AbstractScenario
                    probability::Probability
                    $def
                    function (::Type{$scenarioname})($(vardefs...); probability::AbstractFloat = 1.0)
                        return new(Probability(probability), $(vars...))
                    end
                end
                if $provided_zerodef
                    function Base.zero(::Type{$scenarioname})
                        $zerodef
                    end
                else
                    function Base.zero(::Type{$scenarioname})
                        return $scenarioname(zero.([$(vartypes...)])...; probability = 1.0)
                    end
                end
                if $provided_expectdef
                    function StochasticPrograms.expected(scenarios::Vector{$scenarioname})
                        isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero($scenarioname))
                        $expectdef
                    end
                else
                    function StochasticPrograms.expected(scenarios::Vector{$scenarioname})
                        isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero($scenarioname))
                        return StochasticPrograms.ExpectedScenario(reduce(scenarios) do s1, s2
                                                                   $combine
                                                                   end)
                    end
                end
                function StochasticPrograms.scenariotext(io::IO, scenario::$scenarioname)
                    $textdef
                    return io
                end
            else
                @warn "The scenario type $(string($(Meta.quot(scenarioname)))) was not defined. A user-provided implementation \n\n    function zero(::Type{{$(string($(Meta.quot(scenarioname))))})\n        ...\n    end\n\nis required."
            end
        else
            @warn "The scenario type $(string($(Meta.quot(scenarioname)))) was not defined. A user-provided implementation \n\n    function expected(scenarios::Vector{$(string($(Meta.quot(scenarioname))))})\n        ...\n    end\n\nis required."
        end
    end
    code = prettify(code)
    code = @q begin
        if (@__MODULE__) == $(esc(:Main))
            @everywhere begin
                $code
            end
        else
            $(esc(code))
        end
    end
    return prettify(code)
end
"""
    @zero(def)

Define the additive zero scenario inside a @scenario block using the syntax:
```julia
@zero begin
    ...
    return zero_scenario
end
```

## Examples

The following defines a zero scenario for the example scenario defined in [`@scenario`](@ref)

```julia
@zero begin
    return ExampleScenario(0.0)
end
```

See also [`@scenario`](@ref)
"""
macro zero(def) @warn "@zero should be used inside a @scenario block." end
"""
    @expectation(def)

Define how to form the expected scenario inside a @scenario block. The scenario collection is accessed through the reserved keyword `scenarios`.

```julia
@zero begin
    ...
    return zero_scenario
end
```

## Examples

The following defines expectation for the example scenario defined in [`@scenario`](@ref)

```julia
@expectation begin
    return ExampleScenario(sum([probability(s)*s.ξ for s in scenarios]))
end
```

See also [`@scenario`](@ref)
"""
macro expectation(def) @warn "@expectation should be used inside a @scenario block." end
"""
    @container_scenario([i=..., j=..., ...], expr, probability = 1.0)

Wraps JuMP's `@container` macro to create [`Scenario`](@ref) instances with `DenseAxisArray` or `SparseAxisArray` as underlying data. See `@container` for syntax.

"""
macro container_scenario(args...)
    args, kw_args, probability, requested_container = _extract_kw_args(args)
    @assert length(args) == 2
    @assert isempty(kw_args)
    var, value = args
    code = Containers.parse_container(error, var, esc(value), requested_container)
    return :(Scenario($code, probability = $(esc(probability))))
end
function _extract_kw_args(args)
    kw_args = filter(x -> isexpr(x, :(=)) && x.args[1] != :container && x.args[1] != :probability, collect(args))
    flat_args = filter(x->!isexpr(x, :(=)), collect(args))
    requested_container = :Auto
    probability = 1.0
    for kw in args
        if isexpr(kw, :(=)) && kw.args[1] == :container
            requested_container = kw.args[2]
        end
        if isexpr(kw, :(=)) && kw.args[1] == :probability
            probability = kw.args[2]
        end
    end
    return flat_args, kw_args, probability, requested_container
end
"""
    @sampler(def)

Define a sampler type for some `scenariotype` compatible with StochasticPrograms using the syntax
```julia
@sampler samplername = begin
    ...internals...

    @sample scenariotype begin
        ...
        return scenario
    end
end
```
Any internal state required by the sampler, as well as any specialized constructor, are defined in the @sampler block as they would be in any Julia struct. Define the sample operation inside the [`@sample`](@ref) block and specify the `scenariotype` that the sampler returns. The defined sampler will be available on all Julia processes.

## Examples

The following defines a simple dummy sampler, with some internal weight value, for the scenario defined in [`@scenario`](@ref), and samples one scenario.

```jldoctest; setup = :(@scenario ExampleScenario = ξ::Float64), filter = r".*"
@sampler ExampleSampler = begin
    w::Float64

    ExampleSampler(w::AbstractFloat) = new(w)

    @sample ExampleScenario begin
        @parameters w
        return ExampleScenario(w*randn(), probability = rand())
    end
end
s = ExampleSampler(2.0)
s()

# output

ExampleScenario with probability 0.29
  ξ: 1.48

```

See also: [`@sample`](@ref), [`@scenario`](@ref)
"""
macro sampler(arg)
    @capture(prettify(arg), samplername_Symbol = samplerdef_) || error("Invalid syntax. Expected: sampler = begin samplerdef end")
    scenarioname = :undefined
    sampledefs = Vector{Expr}()
    def = postwalk(prettify(samplerdef)) do x
        if @capture(x, @parameters args__)
            code = Expr(:block)
            for param in args
                push!(code.args, :($param = sampler.$param))
            end
            return code
        elseif @capture(x, @sample sname_ sampledef_)
            scenarioname = sname
            push!(sampledefs, sampledef)
            return @q begin end
        else
            return x
        end
    end
    # Handle sample definition
    length(sampledefs) == 1 || error("Must provide exactly one @sample operation.")
    sampledef = sampledefs[1]
    # Define sampler type
    code = @q begin
        struct $samplername <: AbstractSampler{$scenarioname}
            $def
        end
        function (sampler::$samplername)()
            $sampledef
        end
    end
    code = prettify(code)
    code = @q begin
        if (@__MODULE__) == $(esc(:Main))
            @everywhere begin
                $code
            end
        else
            $(esc(code))
        end
    end
    return prettify(code)
end
"""
    @sample(def)

Define the sample operaton inside a @sampler block, using the syntax
```julia
@sample begin
    ...
    return sampled_scenario
end
```
The sampler object is referenced through the reserved keyword `sampler`, from which any internals can be accessed.
"""
macro sample(def) @warn "@sample should be used inside a @sampler block." end
"""
    @first_stage(def)

Add a first stage model generation recipe to `stochasticprogram` using the syntax
```julia
@first_stage stochasticprogram::StochasticProgram = begin
    ...
end [defer]
```
where JuMP syntax is used inside the block to define the first stage model. During definition, the first stage model is referenced through the reserved keyword `model`.

## Examples

The following defines the first stage model given by:

```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```

```julia
@first_stage sp = begin
    @decision(model, x₁ >= 40)
    @decision(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end
```

See also: [`@second_stage`](@ref)
"""
macro first_stage(arg)
    @capture(arg, sp_Symbol = def_) || error("Invalid syntax. Expected stage, multistage = begin JuMPdef end")
    return esc(:(@stage 1 $arg))
end
"""
    @second_stage(def)

Add a second stage model generation recipe to `stochasticprogram` using the syntax
```julia
@second_stage stochasticprogram::StochasticProgram = begin
    @known var1 var2 ...
    ...
end
```
where JuMP syntax is used inside the block to define the second stage model. During definition, the second stage model is referenced through the reserved keyword `model`.

## Examples

The following defines the second stage model given by:
```math
  minimize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ`` and ``x₁, x₂`` are first stage variables. Two scenarios are added so that two second stage models are generated.

```julia
@second_stage sp = begin
    @known x₁ x₂
    @uncertain q₁ q₂ d₁ d₂
    @variable(model, 0 <= y₁ <= d₁)
    @variable(model, 0 <= y₂ <= d₂)
    @objective(model, Min, q₁*y₁ + q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end
```

See also: [`@first_stage`](@ref)
"""
macro second_stage(arg)
    return esc(:(@stage 2 $arg))
end
"""
    @parameters(def)

Define the problem parameters in a @stage block
```julia
@parameters param1, param2, ...
```
possibly with default values. Any defined parameter without a default value must be supplied as a keyword argument to [`instantiate`](@ref) when creating models.

## Examples

```julia
@parameters d

@parameters begin
    Crops = [:wheat, :corn, :beets]
    Cost = Dict(:wheat=>150, :corn=>230, :beets=>260)
    Budget = 500
end
```

See also [`@decision`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro parameters(def) @warn "@parameters should be used inside a @stage block." end
"""
    @decision(model, expr, args..., kw_args...)

Add a decision variable to `model` described by the expression `expr`. If used inside a [`@stage`](@ref) block, the created variable can be used in subsequent stage blocks. See `@variable` for syntax details.

## Examples

```julia
@decision(model, x >= 40)
```

See also [`@parameters`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro decision(def...) @warn "@decision should be used inside a @stage block." end
macro _decision(args...)
    args = [args...]
    # Check context of decision definition
    known = if args[1] isa Symbol && args[1] == :known
        known = true
    else
        args[1] == :unknown || error("Incorrect usage of @_decision.")
        known = false
    end
    # Remove context from args
    deleteat!(args, 1)
    # Move any set definitions
    set = NoSpecifiedConstraint()
    for (i, arg) in enumerate(args)
        set, arg, found = extract_set(arg)
        if found
            if arg == :()
                deleteat!(args, i)
            else
                args[i] = arg
            end
            break
        end
    end
    # Return expression based on context
    if known
        return esc(:(@variable $((args)...) set = StochasticPrograms.KnownSet()))
    else
        return esc(:(@variable $((args)...) set = StochasticPrograms.DecisionSet(constraint = $set)))
    end
end

"""
    @known(def)

Annotate each decision taken in the previous stage. Any [`@decision`](@ref) included in a [`@stochastic_model`](@ref) definition will implicitly add `@known` annotations to subsequent stages.

## Examples

```julia
@known x₁, x₂
```

See also [`@decision`](@ref), [`@parameters`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro known(args...)
    code = Expr(:block)
    for var in args
        varkey = Meta.quot(var)
        push!(code.args, :($(esc(var)) = $(esc(:model))[$varkey]))
    end
    return code
end
"""
    @uncertain(def)

In a @stage block, annotate each uncertain variable using the syntax
```julia
@uncertain var1, var2, ...
```
or using JuMP's container syntax
```julia
@uncertain ξ[i=..., j=..., ...]
```
This assumes that the [`Scenario`] type is used. Alternatively, user-defined scenarios can be specified by annotating the type. Also, inside a @stochastic_model block, user-defined scenarios can be created during the @uncertain annotation, following [`@scenario`](@ref).

## Examples

```julia
@uncertain q₁ q₂ d₁ d₂

@scenario Simple = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end
@uncertain ξ::SimpleScenario

@stochastic_model begin
    ...
    @uncertain ξ::SimpleScenario = begin
        q₁::Float64
        q₂::Float64
        d₁::Float64
        d₂::Float64
    end
    ...
end
```

See also [`@scenario`](@ref), [`@parameters`](@ref), [`@decision`](@ref), [`@stage`](@ref)
"""
macro uncertain(def) @warn "@uncertain should be used inside a @stage block." end
"""
    @stage(def)

Add a stage model generation recipe to `stochasticprogram` using the syntax
```julia
@stage stage stochasticprogram::StochasticProgram = begin
    @parameters param1 param2 ...
    @decision(model, var) ...
    @uncertain ξ
    ... JuMPdef ...
    ...
end
```
where JuMP syntax is used inside the block to define the stage model. During definition, the second stage model is referenced through the reserved keyword `model`.

## Examples

The following defines the first stage model given by:
```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```
and the second-stage model given by:
```math
  maximize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ`` and ``x₁, x₂`` are first stage variables. Two scenarios are added so that two second stage models are generated.

```jldoctest
ξ₁ = Scenario(q₁ = 24.0, q₂ = 28.0, d₁ = 500.0, d₂ = 100.0, probability = 0.4)
ξ₂ = Scenario(q₁ = 28.0, q₂ = 32.0, d₁ = 300.0, d₂ = 300.0, probability = 0.6)

sp = StochasticProgram([ξ₁, ξ₂])

@stage 1 sp = begin
    @decision(model, x₁ >= 40)
    @decision(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end

@stage 2 sp = begin
    @uncertain q₁ q₂ d₁ d₂
    @variable(model, 0 <= y₁ <= d₁)
    @variable(model, 0 <= y₂ <= d₂)
    @objective(model, Max, q₁*y₁ + q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end

# output

Stochastic program with:
 * 2 decision variables
 * 2 scenarios of type Scenario
Solver is default solver

```

See also: [`@parameters`](@ref), [`@decision`](@ref), [`@uncertain`](@ref)
"""
macro stage(stage, args)
    @capture(args, sp_Symbol = def_) || error("Invalid syntax. Expected stage, multistage = begin JuMPdef end")
    # Decision definitions might require parameter calculations,
    # so we first need to extract and save any such lines
    vardefs = Expr(:block)
    decisiondefs = postwalk(def) do x
        if  @capture(x, @variable(m_Symbol, variabledef__)) ||
            @capture(x, @known(m_Symbol, variabledef__)) ||
            @capture(x, @constraint(m_Symbol, constdef__)) ||
            @capture(x, @expression(m_Symbol, expressiondef__)) ||
            @capture(x, @objective(m_Symbol, objdef__)) ||
            @capture(x, @parameters args__) ||
            @capture(x, @uncertain args__)
            # Skip any line related to the JuMP model, stochastics, or unhandled @parameter lines
            return Expr(:block)
        elseif @capture(x, @decision args__)
            # Handle @decision
            return @q begin
                StochasticPrograms.@_decision known $((args)...)
            end
        else
            # Anything else could be required for decision variable construction, and is therefore saved
            return x
        end
    end
    # Next, handle @decision annotations in main definition
    def = postwalk(def) do x
        if @capture(x, @decision args__)
            return @q begin
                StochasticPrograms.@_decision unknown $((args)...)
            end
        end
        return x
    end
    # Handle parameters
    def = postwalk(def) do x
        if @capture(x, @parameters arg_)
            code = Expr(:block)
            for paramdef in block(prettify(arg)).args
                if @capture(paramdef, key_Symbol = val_) || @capture(paramdef, key_Symbol)
                    push!(code.args, :($key = stage.$key))
                else
                    error("Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            # Extracted parameters might be required for decision variable construction
            pushfirst!(decisiondefs.args, code)
            return code
        elseif @capture(x, @parameters args__)
            code = Expr(:block)
            for param in args
                push!(code.args, :($param = stage.$param))
            end
            # Extracted paremeters might be required for decision variable construction
            pushfirst!(decisiondefs.args, code)
            return code
         elseif @capture(x, @uncertain var_Symbol::t_Symbol)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            return :($var::$t = scenario)
        elseif @capture(x, @uncertain vars__ from scenvar_Symbol::t_Symbol)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                $scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty($scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain vars__ from t_Symbol)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty(scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain var_Symbol::t_Symbol = scenariodef_)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            return @q begin
                $var::$t = scenario
            end
        elseif @capture(x, @uncertain vars__ from scenvar_Symbol::t_Symbol = scenariodef_)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                $scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty($scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain vars__ from t_Symbol = scenariodef_)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty(scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain scenvar_Symbol[args__])
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            idx = Expr(:tuple)
            for idxdef in args
                if @capture(idxdef, p_Symbol = set_) || @capture(idxdef, p_Symbol in set_)
                    push!(idx.args, p)
                end
            end
            code = @q begin
                typeof(scenario) <: Scenario || error("@uncertain declarations of type `@uncertain ξ[i = ..., j = ..., ...]` only support scenarios of type `Scenario`. Consider declaring a custom scenario type.")
                $scenvar = try
                    $scenvar = Containers.@container([$(args...)], scenario.data[$idx...])
                catch err
                    error("Given scenario $scenario does not match @uncertain declaration.")
                end
            end
            return code
        elseif @capture(x, @uncertain args__)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                typeof(scenario) <: Scenario || error("@uncertain declarations of type `@uncertain var1, var2, ...` only support scenarios of type `Scenario`. Consider declaring a custom scenario type.")
            end
            for var in args
                varkey = Meta.quot(var)
                code = @q begin
                    $code
                    $var = try
                        $var = scenario.data[$varkey]
                    catch err
                        error("Given scenario $scenario does not match @uncertain declaration.")
                    end
                end
            end
            return code
        else
            return x
        end
    end
    # The first stage is treated different from the subsequent stages
    generatordefs = if stage == 1
        @q begin
            $(esc(sp)).generator[:stage_1] = ($(esc(:model))::JuMP.Model, $(esc(:stage))) -> begin
                $(esc(def))
	            return $(esc(:model))
            end
        end
    else
        @q begin
            $(esc(sp)).generator[Symbol(:stage_,$stage)] = ($(esc(:model))::JuMP.Model, $(esc(:stage)), $(esc(:scenario))::AbstractScenario) -> begin
                $(esc(def))
	            return $(esc(:model))
            end
        end
    end
    # Create definition code
    code = @q begin
        isa($(esc(sp)), StochasticProgram) || error("Given object is not a stochastic program.")
        $(esc(sp)).generator[Symbol(:stage_,$stage,:_decisions)] = ($(esc(:model))::JuMP.Model, $(esc(:stage))) -> begin
            $(esc(decisiondefs))
	        return $(esc(:model))
        end
        # Stage model generation code
        $generatordefs
        $(esc(sp))
    end
    # Return code
    return code
end
# ========================== #


text:
@recipe function f(fCCF::CCF)

    # plot configuration
    legend:= false

    titlefont:= font(8, "Courier")
    xtickfont:= font(3, "Courier")
    ytickfont:= font(3, "Courier")
    yguidefont:= font(5, "Courier")
    xguidefont:= font(5, "Courier")

    N = fCCF.N
    if split(fCCF.type,'_')[1] == "pacf"
        type = (fCCF.type) == "pacf_step" ? "Stepwise PACF" : "Real PACF"
        type = (fCCF.type) == "pacf_stepwise-real" ? "Stepwise and Real PACF   |s|r|" : type
        auto_prefix = ""
    else
        type = (fCCF.type == "cor") ? "Correlation" : "Covariance"
        auto_prefix = fCCF.auto ? "Auto-" : "Cross-"
    end

    yguide := auto_prefix*type
    xguide := "Lag"

    # Center ticks
    lr = length(fCCF.ccf)
    gap = div(lr,11)
    mp = div(lr+1,2)
    rs = mp:gap:lr
    ls = mp-gap:-gap:1

    if fCCF.auto
        xticks := 1:gap:lr
    else
        xt = vcat(reverse(collect(ls)),collect(rs))
        xticks := (xt, xt .- xt[div(end+1,2)])
    end

    @series begin
        seriestype := :sticks
        ymirror := false
        guide_position := :left
        linewidth := 100/(fCCF.lag+1)

        collect(1:size(fCCF.ccf)[1]+1) .- ((ndims(fCCF.ccf) == 2) ? .1 : 0),
        ndims(fCCF.ccf) == 2 ? fCCF.ccf[:,1] : fCCF.ccf
    end

    if ndims(fCCF.ccf) == 2
        @series begin
            seriestype := :sticks
            ymirror := false
            guide_position := :left
            linewidth := 100/(fCCF.lag+1)

            collect(1:size(fCCF.ccf)[1]+1) .+ .1,
            fCCF.ccf[:,2]
        end
    end

    if (fCCF.type == "cor") | (split(fCCF.type,'_')[1] == "pacf")
        lg = fCCF.lag
        a1 = fCCF.alpha[1]; c1 = fCCF.ci[1]
        a2 = fCCF.alpha[2]; c2 = fCCF.ci[2]
        dc = c2-c1

        @series begin
            seriestype := :hline
            seriescolor := :black
            linealpha := 0.5
            linestyle := :dash
            [-c1,c1]
        end

        @series begin
            seriestype := :hline
            seriescolor := :black
            linealpha := 0.5
            linestyle := :dot
            ax = fCCF.auto ? lg : 2*lg + 1
            annotations :=
                [(ax,c1+dc/8,string("CI %",a1*100),font(3, "Courier")),
                 (ax,c2+dc/8,string("CI %",a2*100),font(3, "Courier"))]
            [-c2,c2]
        end
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export button

"""
    button(;kwargs...)
    button(children::Any;kwargs...)
    button(children_maker::Function;kwargs...)

A Button component.
Carbon Button
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): Specify the content of your Button
- `id` (String; required): The id of the button
- `style` (Dict; optional): jsx Style
- `className` (String; optional): Specify an optional className to be added to your Button
- `disabled` (Bool; optional): Specify whether the Button should be disabled, or not
- `dangerDescription` (String; optional): Specify the message read by screen readers for the danger button variant
- `href` (String; optional): Optionally specify an href for your Button to become an `<a>` element
- `isSelected` (Bool; optional): Specify whether the Button is currently selected
- `kind` (a value equal to: 'primary', 'secondary', 'danger', 'ghost', 'danger--primary', 'danger--ghost', 'danger--tertiary', 'tertiary'; required): Specify the kind of Button you want to create
- `role` (String; optional): Optional prop to specify the role of the Button
- `size` (a value equal to: 'default', 'field', 'small', 'sm', 'lg', 'xl'; optional): Specify the size of the button, from a list of available sizes.
For `default` buttons, this prop can remain unspecified.
- `tooltipAlignment` (a value equal to: 'start', 'center', 'end'; optional): Specify the alignment of the tooltip to the icon-only button.
Can be one of: start, center, or end.
- `tooltipPosition` (a value equal to: 'top', 'right', 'bottom', 'left'; optional): Specify the direction of the tooltip for icon-only buttons.
Can be either top, right, bottom, or left.
- `type` (a value equal to: 'button', 'reset', 'submit'; optional): Optional prop to specify the type of the Button
- `n_clicks` (Real; optional): Number of times this button has been clicked
"""
function button(; kwargs...)
        available_props = Symbol[:children, :id, :style, :className, :disabled, :dangerDescription, :href, :isSelected, :kind, :role, :size, :tooltipAlignment, :tooltipPosition, :type, :n_clicks]
        wild_props = Symbol[]
        return Component("button", "Button", "dash_carbon_components", available_props, wild_props; kwargs...)
end

button(children::Any; kwargs...) = button(;kwargs..., children = children)
button(children_maker::Function; kwargs...) = button(children_maker(); kwargs...)

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
function isprerelease(version::VersionNumber)::Bool
    prerelease = version.prerelease
    isempty(prerelease) && return false
    for x in prerelease
        if !isempty(strip(x))
            return true
        end
    end
    return false
end

function compare_versions(master_version::VersionNumber,
                          head_version::VersionNumber;
                          allow_unchanged_prerelease::Bool)::Nothing
    if head_version > master_version
        @info("Version number has increased")
    elseif head_version == master_version
        if isprerelease(head_version) && isprerelease(master_version) && allow_unchanged_prerelease
            @info("Version number did not change, but it is a prerelease so this is allowed")
        else
            throw(ErrorException("Version number is unchanged, which is not allowed"))
        end
    else
        throw(ErrorException("Version number decreased, which is not allowed"))
    end
    return nothing
end


text:
# This file was generated, do not modify it. # hide
# hideall
c = """
col1,col2,col3,col4,col5,col6,col7,col8
,1,1.0,1,one,2019-01-01,2019-01-01T00:00:00,true
,2,2.0,2,two,2019-01-02,2019-01-02T00:00:00,false
,3,3.0,3.14,three,2019-01-03,2019-01-03T00:00:00,true
"""
fpath, = mktemp()
write(fpath, c);


text:
function tron_dscal(n,da,dx,incx)

    # Scales a vector by a constant.
    # It uses unrolled loops for increment equal to one.
    # Jack Dongarra, LINPACK, 3/11/78.
    # Modified 3/93 to return if incx <= 0.
    # Modified 12/3/93, array(1) declarations changed to array(*)

    if n <= 0 || incx <= 0
        return
    end

    if incx != 1

        # Code for increment not equal to 1

        nincx = n*incx
        for i=1:incx:nincx
            dx[i] = da*dx[i]
        end

        return

    else

        # Code for increment equal to 1
        # Clean-up loops

        m = mod(n,5)
        if m != 0
            for i=1:m
                dx[i] = da*dx[i]
            end
            if n < 5
                return
            end
        end

        mp1 = m + 1
        for i=mp1:5:n
            dx[i] = da*dx[i]
            dx[i + 1] = da*dx[i + 1]
            dx[i + 2] = da*dx[i + 2]
            dx[i + 3] = da*dx[i + 3]
            dx[i + 4] = da*dx[i + 4]
        end

        return
    end
end

if isequal(BLAS_LIBRARY, :OpenBlas)
    dscal(n,da,dx,incx) = BLAS.scal!(n, da, dx, incx)
    dscal(n,da,dx::CuArray,incx) = CUBLAS.scal!(n, da, dx)
else
    dscal(n,da,dx,incx) = tron_dscal(n, da, dx, incx)
end

@inline function dscal(n::Int,da::Float64,dx::CuDeviceArray{Float64,1},incx::Int)
    tx = threadIdx().x
    ty = threadIdx().y

    # Ignore incx for now.
    if tx <= n && ty == 1
        @inbounds dx[tx] = da*dx[tx]
    end
    CUDA.sync_threads()

    return
end


text:
type LTIModel <: DischargeModel
    M::Int
    β::Vector{Float64}
    H::Matrix{Float64}
    Q::Vector{Float64}
end

LTIModel(M,N) = LTIModel(M,zeros(M),zeros(M,N),zeros(N))

function StatsBase.fit(M::LTIModel,H::Matrix{Float64},Q::Vector{Float64})
    LTIModel(M.M,H'\Q,H,Q)
end

function StatsBase.fit!(M::LTIModel,H::Matrix{Float64},Q::Vector{Float64})
    M.β = H'\Q
    M.H = H
    M.Q = Q
    M
end

StatsBase.coef(M::LTIModel) = M.β
StatsBase.nobs(M::LTIModel) = length(M.Q)
StatsBase.model_response(M::LTIModel) = M.Q
StatsBase.predict(M::LTIModel) = M.H'M.β
StatsBase.predict(M::LTIModel,H) = H'M.β
StatsBase.residuals(M::LTIModel) = M.Q-predict(M)
StatsBase.residuals(M::LTIModel,H,Q) = Q-predict(M,H)


text:
# Single Models ---------------------------------------------------------------------------------

function start_simple(model::Union{Sem, SemForwardDiff, SemFiniteDiff}; kwargs...)
    return start_simple(
        model.observed,
        model.imply,
        model.diff,
        model.loss.functions...,
        kwargs...)
end

function start_simple(observed, imply::Union{RAM, RAMSymbolic}, diff, args...; kwargs...)
    return start_simple(imply.ram_matrices; kwargs...)
end

# Ensemble Models --------------------------------------------------------------------------------
function start_simple(model::SemEnsemble; kwargs...)

    start_vals = []

    for sem in model.sems
        push!(start_vals, start_simple(sem; kwargs...))
    end

    has_start_val = [.!iszero.(start_val) for start_val in start_vals]

    start_val = similar(start_vals[1])
    start_val .= 0.0

    for (j, indices) in enumerate(has_start_val)
        start_val[indices] .= start_vals[j][indices]
    end

    return start_val

end

function start_simple(
    ram_matrices::RAMMatrices;
    start_loadings = 0.5,
    start_regressions = 0.0,
    start_variances_observed = 1,
    start_variances_latent = 0.05,
    start_covariances_observed = 0.0,
    start_covariances_latent = 0.0,
    start_covariances_obs_lat = 0.0,
    start_means = 0.0,
    kwargs...)

    A_ind, S_ind, F_ind, M_ind, parameters =
        ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.F_ind, ram_matrices.M_ind, ram_matrices.parameters

    n_par = length(parameters)
    start_val = zeros(n_par)
    n_var, n_nod = ram_matrices.size_F

    C_indices = CartesianIndices((n_nod, n_nod))

    for i in 1:n_par
        if length(S_ind[i]) != 0
            # use the first occurence of the parameter to determine starting value
            c_ind = C_indices[S_ind[i][1]]
            if c_ind[1] == c_ind[2]
                if c_ind[1] ∈ F_ind
                    start_val[i] = start_variances_observed
                else
                    start_val[i] = start_variances_latent
                end
            else
                o1 = c_ind[1] ∈ F_ind
                o2 = c_ind[2] ∈ F_ind
                if o1 & o2
                    start_val[i] = start_covariances_observed
                elseif !o1 & !o2
                    start_val[i] = start_covariances_latent
                else
                    start_val[i] = start_covariances_obs_lat
                end
            end
        elseif length(A_ind[i]) != 0
            c_ind = C_indices[A_ind[i][1]]
            if (c_ind[1] ∈ F_ind) & !(c_ind[2] ∈ F_ind)
                start_val[i] = start_loadings
            else
                start_val[i] = start_regressions
            end
        elseif !isnothing(M_ind) && (length(M_ind[i]) != 0)
            start_val[i] = start_means
        end
    end
    return start_val
end


text:
"""
    contact_point
Find approximate contact point. Current available methods (cp) are `threshold' (i.e. the contact point is identified by fixing a force threshold) or `hertz' (i.e. the contact point is identified by fitting a hertz model to the curve).
"""
function contact_point(self::RheoTimeData, interface::Interface, cp::String, param::NamedTuple; savelog = true,  sec = :)

    contactmodels = Dict("threshold" => contact_threshold,
                         "hertz" => contact_hertz);

    data = self[interface];

    cp_index = contactmodels[cp](data[:f], data[:d], param);

    if cp_index != nothing
        t = self.t[cp_index:end] .- self.t[cp_index];
        d = data[:d][cp_index:end] .- data[:d][cp_index];
        f = data[:f][cp_index:end] .- data[:f][cp_index];
    else
        t = self.t[:]
        d = data[:d][:]
        f = data[:f][:]
    end

    ϵ,σ = interface.to_ϵσ(d, f);

    log = self.log == nothing ? nothing : [self.log; RheoLogItem( (type=:process, funct=:contact_point, params=(cp = cp, param = param ), keywords=(sec = sec, ) ),
    (comment="Contact point detection and data shift",) ) ]

 return RheoTimeData(σ, ϵ, t, log)

end

"""
    contact_threshold(f::Array{Float64,1}, d::Array{Float64,1}; _param::NamedTuple)

Applies force threshold to find approximate contact point.
"""
function contact_threshold(force::Array{RheoFloat}, disp::Array{RheoFloat}, param::NamedTuple)

    # get index of contact
    cp_index = findfirst(force .>= param.threshold);

    # return
    return cp_index

end

"""
    contact_hertz(f::Array{Float64,1}, δ::Array{Float64,1}; _param::Float64 = NaN64)

    Apply spherical Hertz contact model to find approximate contact point.
"""
function contact_hertz(f::Array{RheoFloat}, δ::Array{RheoFloat}, param::NamedTuple)

    # regularise displacement
    δ = δ .- minimum(δ);

    # radius of sphere
    R = param.R;

    # poisson ratio - assume incompressible
    ν = 0.5;

    # section of the curve to be used to identify the first approximation
    s = param.s

    # get approximate δ0 and YM values
    (init_δ₀, init_YM, tilt_approx, offset_approx) =  hertz_approx(f, δ, R, ν, s);

    # least squares fit
    opt = Opt(:LN_SBPLX, 4); # 4 parameters to fit

    bounds_param = 0.25 # +/- limits allowance from initial contact point guess
    lower_bounds!(opt, [init_δ₀*(1 - bounds_param), 0.0, -5.0, -5.0]);
    upper_bounds!(opt, [init_δ₀*(1 + bounds_param), 1e10, 5.0, 5.0]);

    xtol_rel!(opt,1e-5);

    min_objective!(opt, (params, grad) -> hertz_cost(f, δ, R, ν, params[1], params[2], params[3], params[4]));

    params_init = [init_δ₀, init_YM, tilt_approx, offset_approx];

    (minf, minx, ret) = optimize(opt, params_init);

    # return index of contact point
    δ₀_index = argmin(abs.(δ .- minx[1]));

    return (δ₀_index)

end

###########################################
# Convenience functions for contact_hertz #
###########################################
"""
Generate approximate Hertz fit a, b, E, cp to use as initial conditions in proper fit.
"""

function hertz_approx(f::Array{RheoFloat,1}, δ::Array{RheoFloat,1}, R::RheoFloat, ν::RheoFloat, s::RheoFloat)

    # approximate contact point, assume linear near end of approach section
    # and trace line down to axis
    deriv = RHEOS.derivCD(RHEOS.derivCD(f, δ),δ);
    value_max, index_max = findmax(deriv);

    # #coeff = (1-index_max/length(f))
    last_section = round(Int64, s*length(f)):length(f)

    avg_gradient = sum(RHEOS.derivBD(f[last_section], δ[last_section]))/length(last_section);

    c = f[end] - avg_gradient*δ[end];

    approx_δ₀ = -c/avg_gradient;

    # transform to find approx YM
    δᵦ = (heaviside(δ .- approx_δ₀).*δ).^(3/2);

    prefactor = (4/3)*sqrt(R)/(1 - ν^2);

    avg_gradient_transformed = sum(RHEOS.derivBD(f[last_section], δᵦ[last_section]))/length(last_section);

    approx_YM = avg_gradient_transformed/prefactor;

    # estimate offset and tilt
    first_section = 1:round(Int64, 0.5*length(f));

    df_dh = RHEOS.derivBD(f, δ);

    tilt_approx = sum(df_dh[first_section])/length(first_section);

    offset_approx = sum(f[first_section])/length(first_section);

    if approx_δ₀<0
        approx_δ₀ = 0
    end

    if approx_YM<0
        approx_YM = 50
    end

    # return
    approx_δ₀, approx_YM, tilt_approx, offset_approx;
end

"""
Heaviside function
"""
function heaviside(x)
   0.5*(sign.(x) .+ 1);
end

"""
Hertz spherical model
"""
function hertz_model(δ::Array{RheoFloat,1}, R::RheoFloat, ν::RheoFloat, δ₀::RheoFloat, E::RheoFloat, tilt::RheoFloat, offset::RheoFloat)
    # define prefactor for Hertz sphere model
    prefactor = (4/3)*sqrt(R)/(1 - ν^2);

    δ_shift = δ - δ₀;

    # force
    F = tilt*δ + offset + prefactor*E*((δ_shift.*heaviside(δ_shift)).^(3/2));
end

"""
Hertz cost function
"""
function hertz_cost(f::Array{RheoFloat,1}, δ::Array{RheoFloat,1}, R::RheoFloat, ν::RheoFloat, δ₀::RheoFloat, E::RheoFloat, tilt::RheoFloat, offset::RheoFloat)
    # difference
    Δ = f - hertz_model(δ, R, ν, δ₀, E, tilt, offset);

    # regularisation parameter for offset and tilt
    λ = 0.1;

    # sum of squares
    cost = sum(Δ.^2); # Maybe not needed + λ*(tilt + offset)^2
end


text:
module testMarkowitz

using Base.Test
using DataFrames
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running markowitz tests\n")

df = convert(DataFrame, rand(500, 4))
mus = AssetMgmt.mean(df, 1)
covMatr = AssetMgmt.cov(df)
wgts = AssetMgmt.gmv(mus, covMatr) # should be close to
											  # [0.25 0.25 0.25 0.25]

###############
## real data ##
###############

include("/home/chris/research/julia/EconDatasets/src/EconDatasets.jl")
logRet = EconDatasets.dataset("SP500")

## transform to discrete non-percentage returns
discRet = exp(logRet/100).-1

mus = mean(discRet, 1)
covMatr = cov(discRet)

wgts = AssetMgmt.gmv(mus, covMatr)

end


text:
using SparseArrays
using MatrixNetworks

function hypergraph_largest_component(H::SparseMatrixCSC)
    m,n = size(H)
    A = [spzeros(m,m) H; sparse(H') spzeros(n,n)]
    lcc, pcc = largest_component(A)
    p_nodes = pcc[m+1:end]
    p_edges = pcc[1:m]
    H = H[p_edges,p_nodes]
    return H, p_edges, p_nodes
end

# These are some useful functions for converting back and forth between
# different ways to store hypergraphs.
function NeighborList(He2n::SparseArrays.SparseMatrixCSC{Float64,Int64},Hn2e::SparseArrays.SparseMatrixCSC{Float64,Int64})
    """
    NeighborList: given hypergraph (edge x node) indicence and (node x edge)
    incidence, return a list of neighbors for each node.
    """
    m,n = size(He2n)
    rp = He2n.rowval
    ci = He2n.colptr
    Neighbs = Vector{Vector{Int64}}()
    for i = 1:n
        # Get list of edges that node i is in
        Edges = rp[ci[i]:ci[i+1]-1]

        # build a set of neighboring nodes
        NeighbSet = Vector{Int64}()
        for j = 1:length(Edges)
            ej = Edges[j]               # For each node ej adjacent to I...
            nodes = getnodes(Hn2e,ej)   # ...get nodes in edge ej
            append!(NeighbSet,nodes)
        end
        push!(Neighbs, sort(setdiff(unique(NeighbSet),i)))
    end

    return Neighbs
end

function NeighborList(node2edge::Vector{Vector{Int64}},edge2node::Vector{Vector{Int64}})
    """
    NeighborList: given edge to node list and node to edge list, compute
    neighbors
    """
    n = length(node2edge)
    m = length(edge2node)
    Neighbs = Vector{Vector{Int64}}()
    for i = 1:n
        # Get list of edges that node i is in
        Edges = node2edge[i]

        # build a set of neighboring nodes
        NeighbSet = Vector{Int64}()
        for j = 1:length(Edges)
            ej = Edges[j]               # For each node ej adjacent to I...
            nodes = edge2node[ej]       # ...get nodes in edge ej
            append!(NeighbSet,nodes)
        end
        push!(Neighbs, sort(setdiff(unique(NeighbSet),i)))
    end

    return Neighbs
end

function getedges(He2n::SparseArrays.SparseMatrixCSC{Float64,Int64},I::Int64)
    """
    Given the edge-by-node incidence matrix He2n, return the set of hyperedges
    that node I is in.
    """
    first = He2n.colptr[I]
    last = He2n.colptr[I+1]-1
    edges = He2n.rowval[first:last]

    return edges
end

function getnodes(Hn2e::SparseArrays.SparseMatrixCSC{Float64,Int64},J::Int64)
    """
    Given the node-by-edge incidence matrix Hn2e, return the set of nodes
    contained in hyperedge J
    """
    first = Hn2e.colptr[J]
    last = Hn2e.colptr[J+1]-1
    nodes = Hn2e.rowval[first:last]
    mult = Hn2e.nzval[first:last]
    edge = Vector{Int64}()
    # need to adjust for multiplicities
    for t = 1:length(nodes)
        node = nodes[t]
        for k = 1:mult[t]
            push!(edge,node)
        end
    end

    return edge
end

function incidence2elist(H::SparseArrays.SparseMatrixCSC{Float64,Int64},nodelist::Bool=false)
    """
    This converts a hypergraph in incidence matrix form to hyperedge list form.
    Incidence matrix form:  H[e,u] = 1  iff node u is in hyperedge e
    Edgelist form: Hyperedges[j] = array of nodes in hyperedge j

    nodelist == true: means you actually want to get a map from node IDs to hyperedge ids
    """
    if ~nodelist
        # unless you want the node2edge map, transpose first
        H = SparseArrays.sparse(H')
    end
    rp = H.rowval
    ci = H.colptr
    nz = H.nzval
    Hyperedges = Vector{Vector{Int64}}()
    n,m = size(H)

    for i = 1:m
        startedge = ci[i]
        endedge = ci[i+1]-1
        nodes = rp[startedge:endedge]
        mult = nz[startedge:endedge]
        edge = Vector{Int64}()
        # need to adjust for multiplicities
        for t = 1:length(nodes)
            node = nodes[t]
            for k = 1:mult[t]
                push!(edge,node)
            end
        end

        push!(Hyperedges,edge)
    end
    return Hyperedges
end

function elist2incidence(Hyperedges::Vector{Vector{Int64}}, N::Int64)
    """
    Converts a hyperedge list into a binary incidence matrix for the hypergraph.

    This is the exact inverse of incidence2elist
    """
    U = Vector{Int64}()
    E = Vector{Int64}()
    M = length(Hyperedges)
    for enum = 1:length(Hyperedges)
        e = Hyperedges[enum]
        for node in e
            push!(U,node)
            push!(E,enum)
        end
    end

    H = SparseArrays.sparse(E,U,ones(length(U)),M,N)
    return H
end

function RandomSplittingFunctions(Edges,probs)
    ## Generate random splitting functions
    EdgesW = Vector{Vector{Float64}}()
    for edge in Edges
        k = length(edge)
        w = random_scb_function(k,probs[1],probs[2])
        while ~check_cb_submodular(w,0)
            w = random_scb_function(k,probs[1],probs[2])
        end
        push!(EdgesW,w)
    end
    return EdgesW
end

function CliqueExpansion(Edges::Vector{Vector{Int64}},n::Int64,weighted::Bool=true,binary::Bool=false)
    """
    Weighted clique expansion where a hyperedge e is expanded to a
    weighted clique with each edge having weight 1/(|e| - 1)
    """
    I = Int64[]
    J = Int64[]
    V = Float64[]
    for edge in Edges
        k = length(edge)
        for i = 1:k-1
            ei = edge[i]
            for j = i+1:k
                ej = edge[j]
                push!(I,ei)
                push!(J,ej)
                if weighted
                    push!(V, 1/k)
                else
                    push!(V, 1)
                end
            end
        end
    end
    A = SparseArrays.sparse(I,J,V,n,n)
    for i = 1:n
        A[i, i] = 0.0
    end
    SparseArrays.dropzeros!(A)
    A = SparseArrays.sparse(A+A')
    if binary
        I, J, V = SparseArrays.findnz(A)
        A = SparseArrays.sparse(I, J, 1, n, n)
    end
    return A
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
module Objects
  using CG;
  export houseOfSantaClaus;
  function houseOfSantaClaus()
    return Object(Vec4f(-1, -1, 0, 1),
    Vec4f(1, -1, 0, 1),
    Vec4f(-1, 1, 0, 1),
    Vec4f(0, 2, 0, 1),
    Vec4f(1, 1, 0, 1),
    Vec4f(1, 1, 0, 1),
    Vec4f(-1, -1, 0, 1),
    Vec4f(-1, 1, 0, 1),
    Vec4f(1, 1, 0, 1),
    Vec4f(1, -1, 0, 1));
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
##############################################################################
##
## Gensys solver adapted from phactsolver.m
##
##############################################################################

function gensys(Γ0, Γ1, c, Ψ, Π; clean = true, continuous = true, check_existence = true, check_uniqueness = true)

    if clean
        @sprintf "Converting to Reduced Form"
        redundant = (maxabs(Γ0, 2) .== 0) & (maxabs(Ψ, 2) .== 0)
        base = nullspace(Γ1[redundant, :])
        try
            Γ0 = lufact!(At_mul_B(base, Γ0 * base))
            Γ1 = Γ0 \ At_mul_B(base, Γ1 * base)
            Ψ = Γ0 \ At_mul_B(base, Ψ)
            Π = Γ0 \ At_mul_B(base, Π)
            c = Γ0 \ At_mul_B(base, c)
        catch
            error("Wrong Form. Try running Gensys")
        end
    else
        Γ1 = Γ0 \ Γ1
    end
    n = size(Γ1, 1)

    # Schur Decomposition
    Γ1 = schurfact!(Γ1)
    if continuous
        select = real(Γ1[:values]) .< 0
    else
        select = abs(Γ1[:values]) .< 1
    end
    ordschur!(Γ1, select)
    n1 = sum(select)
    Γ1vectors = Γ1[:vectors]

    # Compute G1
    G1 = real(
        A_mul_Bt(Γ1vectors * Γ1[:Schur] * diagm(vcat(ones(n1), zeros(n - n1))), Γ1vectors))

    # Compute impact
    u2 = Γ1vectors[:, (n1 + 1):n]
    etawt = svdfact!(At_mul_B(u2, Π))
    ueta, deta, veta  = etawt[:U], etawt[:S], etawt[:V]
    impact = real(-Π * veta * (diagm(deta) \ ueta') * At_mul_B(u2, Ψ) + Ψ)

    # check existence
    if check_existence
        temp = svdfact!(At_mul_B(u2, Ψ))
        uz, dz, vz =  temp[:U], temp[:S], temp[:V]
        existence = vecnorm(uz - ueta * At_mul_B(ueta, uz), 2) < (sqrt(eps()) * 10 * n)
    end

    # check uniqueness
    if check_uniqueness
        u1 = Γ1vectors[:, 1:n1]
        temp = svdfact!(At_mul_B(u1, Π))
        dont, deta1, veta1 = temp[:U], temp[:S], temp[:V]
        uniqueness = vecnorm(veta1 - veta * At_mul_B(veta, veta1), 2) < (sqrt(eps()) * 10  *  n)
    end

    if clean
        G1 = base * A_mul_Bt(G1, base)
        impact = base * impact
    end

    return G1, impact, existence, uniqueness
end