text:
# Regular Expression
# ==================
#
# Regular expression sequence search tools.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# String Decorators
# -----------------

macro biore_str(pat, opt...)
    if isempty(opt)
        error("symbol option is required: d(na), r(na), or a(a)")
    end
    @assert length(opt) == 1
    opt = opt[1]

    if opt ∈ ("d", "dna")
        :($(RE.Regex{DNA}(pat, :pcre)))
    elseif opt ∈ ("r", "rna")
        :($(RE.Regex{RNA}(pat, :pcre)))
    elseif opt ∈ ("a", "aa")
        :($(RE.Regex{AminoAcid}(pat, :pcre)))
    else
        error("invalid symbol option: $(opt)")
    end
end

macro prosite_str(pat)
    :($(RE.Regex{AminoAcid}(pat, :prosite)))
end

module RE

using Printf: @sprintf

import BioGenerics, BioSequences

# Syntax tree
# -----------

mutable struct SyntaxTree
    head::Symbol
    args::Vector{Any}

    function SyntaxTree(head, args)
        @assert head ∈ (
            :|,
            :*, :+, Symbol("?"), :range,
            Symbol("*?"), Symbol("+?"), Symbol("??"), Symbol("range?"),
            :set, :compset, :sym, :bits,
            :capture, :nocapture, :concat, :head, :last)
        return new(head, args)
    end
end

expr(head, args) = SyntaxTree(head, args)

function charset(s)
    set = Set{Char}()
    union!(set, s)
    union!(set, lowercase(s))
    return set
end

# list of symbols available for each symbol type
const symbols = IdDict{Any, Set{Char}}(
    BioSequences.DNA => charset("ACGTMRWSYKVHDBN"),
    BioSequences.RNA => charset("ACGUMRWSYKVHDBN"),
    BioSequences.AminoAcid     => charset("ARNDCQEGHILKMFPSTWYVOUBJZX"))

macro check(ex, err)
    esc(quote
        if !$ex
            throw($err)
        end
    end)
end

function parse(::Type{T}, pat::AbstractString) where {T}
    parens = Char[]  # stack of parens
    ex, _ = parserec(T, pat, iterate(pat), parens)
    @check isempty(parens) ArgumentError("'(' is not closed")
    return ex
end

function parserec(::Type{T}, pat, s, parens) where {T}
    args = []
    while s !== nothing
        c, state = s
        s = iterate(pat, state)
        if c == '*'
            @check !isempty(args) ArgumentError("unexpected '*'")
            arg = pop!(args)
            push!(args, expr(:*, [arg]))
        elseif c == '+'
            @check !isempty(args) ArgumentError("unexpected '+'")
            arg = pop!(args)
            push!(args, expr(:+, [arg]))
        elseif c == '?'
            @check !isempty(args) ArgumentError("unexpected '?'")
            arg = pop!(args)
            if arg.head ∈ (:*, :+, Symbol("?"), :range)
                # lazy quantifier
                push!(args, expr(Symbol(arg.head, '?'), arg.args))
            else
                # zero-or-one quantifier
                push!(args, expr(Symbol("?"), [arg]))
            end
        elseif c == '{'
            @check !isempty(args) ArgumentError("unexpected '{'")
            rng, s = parserange(pat , s)
            arg = pop!(args)
            push!(args, expr(:range, [rng, arg]))
        elseif c == '|'
            arg1 = expr(:concat, args)
            arg2, s = parserec(T, pat, s, parens)
            args = []
            push!(args, expr(:|, [arg1, arg2]))
        elseif c == '['
            setexpr, s = parseset(T, pat, s)
            push!(args, setexpr)
        elseif c == '('
            push!(parens, '(')
            head = :capture
            if peek(pat, s) == '?'
                _, state = s
                s = iterate(pat, state)
                if peek(pat, s) == ':'
                    head = :nocapture
                    _, state = s
                    s = iterate(pat, state)
                end
            end
            arg, s = parserec(T, pat, s, parens)
            push!(args, expr(head, [arg]))
        elseif c == ')'
            @check !isempty(parens) && parens[end] == '(' ArgumentError("unexpected ')'")
            pop!(parens)
            return expr(:concat, args), s
        elseif c == '^'
            push!(args, expr(:head, []))
        elseif c == '$'
            push!(args, expr(:last, []))
        elseif isspace(c)
            # skip
        elseif c ∈ symbols[T]
            push!(args, expr(:sym, [convert(T, c)]))
        else
            throw(ArgumentError("unexpected input: '$c'"))
        end
    end
    return expr(:concat, args), s
end

function parserange(pat, s)
    lo = hi = -1
    comma = false
    while s !== nothing
        c, state = s
        s = iterate(pat, state)
        if isdigit(c)
            d = c - '0'
            if comma
                if hi < 0
                    hi = 0
                end
                hi = 10hi + d
            else
                if lo < 0
                    lo = 0
                end
                lo = 10lo + d
            end
        elseif c == ','
            comma = true
        elseif c == '}'
            break
        else
            throw(ArgumentError("unexpected input: '$c'"))
        end
    end
    if comma
        if lo < 0 && hi < 0
            throw(ArgumentError("invalid range"))
        elseif lo ≥ 0 && hi < 0
            return (lo,), s
        elseif lo < 0 && hi ≥ 0
            return (0, hi), s
        else  # lo ≥ 0 && hi ≥ 0
            if lo > hi
                throw(ArgumentError("invalid range"))
            end
            return (lo, hi), s
        end
    else
        return lo, s
    end
end

function peek(pat, s)
    if s === nothing
        throw(ArgumentError("unexpected end of pattern"))
    end
    return s[1]
end

function parseset(::Type{T}, pat, s) where {T}
    if peek(pat, s) == '^'
        head = :compset
        _, state = s
        s = iterate(pat, state)
    else
        head = :set
    end
    set = T[]
    while s !== nothing
        c, state = s
        s = iterate(pat, state)
        if c ∈ symbols[T]
            push!(set, convert(T, c))
        elseif c == ']'
            break
        else
            throw(ArgumentError("unexpected input: '$c'"))
        end
    end
    return expr(head, set), s
end

function parse_prosite(pat)
    s = iterate(pat)
    args = []
    while s !== nothing
        c, state = s
        s = iterate(pat, state)
        if c == '['
            set, s = parseset_prosite(pat, s, ']')
            push!(args, set)
        elseif c == '{'
            set, s = parseset_prosite(pat, s, '}')
            push!(args, set)
        elseif c == '('
            @check !isempty(args) ArgumentError("unexpected '('")
            rng, s = parserange_prosite(pat, s)
            arg = pop!(args)
            push!(args, expr(:range, [rng, arg]))
        elseif c == '-'
            # concat
            continue
        elseif c == 'x'
            push!(args, expr(:sym, [BioSequences.AA_X]))
        elseif c == '<'
            push!(args, expr(:head, []))
        elseif c == '>'
            push!(args, expr(:last, []))
        elseif c ∈ symbols[BioSequences.AminoAcid]
            push!(args, expr(:sym, [convert(BioSequences.AminoAcid, c)]))
        else
            throw(ArgumentError("unexpected input: '$c'"))
        end
    end
    return expr(:concat, args)
end

function parserange_prosite(pat, s)
    lo = hi = -1
    comma = false
    while s !== nothing
        c, state = s
        s = iterate(pat, state)
        if isdigit(c)
            d = c - '0'
            if comma
                if hi < 0
                    hi = 0
                end
                hi = 10hi + d
            else
                if lo < 0
                    lo = 0
                end
                lo = 10lo + d
            end
        elseif c == ','
            comma = true
        elseif c == ')'
            break
        else
            throw(ArgumentError("unexpected input: '$c'"))
        end
    end
    if comma
        if lo < 0 || hi < 0
            throw(ArgumentError("invalid range"))
        end
        return (lo, hi), s
    else
        if lo < 0
            throw(ArgumentError("invalid range"))
        end
        return lo, s
    end
end

function parseset_prosite(pat, s, close)
    set = BioSequences.AminoAcid[]
    while s !== nothing
        c, state = s
        s = iterate(pat, state)
        if c ∈ symbols[BioSequences.AminoAcid]
            push!(set, convert(BioSequences.AminoAcid, c))
        elseif c == close
            if close == ']'
                return expr(:set, set), s
            elseif close == '}'
                return expr(:compset, set), s
            end
            @assert false
        else
            throw(ArgumentError("unexpected input: '$c'"))
        end
    end
end

function bits2sym(::Type{T}, bits::UInt32) where {T}
    for x in BioSequences.alphabet(T)
        if BioSequences.compatbits(x) == bits
            return x
        end
    end
    error("bits are not found")
end

mask(::Type{T}) where {T<:BioSequences.NucleicAcid} = (UInt32(1) << 4) - one(UInt32)
@assert convert(Int, BioSequences.AA_U) + 1 == 22  # check there are 22 unambiguous amino acids
mask(::Type{BioSequences.AminoAcid}) = (UInt32(1) << 22) - one(UInt32)

function desugar(::Type{T}, tree::SyntaxTree) where {T}
    head = tree.head
    args = tree.args
    if head == :+
        # e+ => ee*
        head = :concat
        args = [args[1], expr(:*, [args[1]])]
    elseif head == Symbol("+?")
        # e+? => ee*?
        head = :concat
        args = [args[1], expr(Symbol("*?"), [args[1]])]
    elseif head == Symbol("?")
        # e? => e|
        head = :|
        args = [args[1], expr(:concat, [])]
    elseif head == Symbol("??")
        # e?? => |e
        head = :|
        args = [expr(:concat, []), args[1]]
    elseif head == :sym
        head = :bits
        args = [BioSequences.compatbits(args[1])]
    elseif head == :set
        bits = UInt32(0)
        for arg in args
            bits |= BioSequences.compatbits(arg)
        end
        head = :bits
        args = [bits]
    elseif head == :compset
        bits = UInt32(0)
        for arg in args
            bits |= BioSequences.compatbits(arg)
        end
        head = :bits
        args = [~bits & mask(T)]
    elseif head == :nocapture
        head = :concat
    elseif head == :range || head == Symbol("range?")
        rng = args[1]
        pat = args[2]
        greedy = head == :range
        if isa(rng, Int)
            # e{m} => eee...e
            #         |<-m->|
            head = :concat
            args = fill(pat, rng)
        elseif isa(rng, Tuple{Int})
            # e{m,} => eee...ee*   (greedy)
            #          |<-m->|
            # e{m,} => eee...ee*?  (lazy)
            #          |<-m->|
            head = :concat
            args = fill(pat, rng[1])
            if greedy
                push!(args, expr(:*, [pat]))
            else
                push!(args, expr(Symbol("*?"), [pat]))
            end
        elseif isa(rng, Tuple{Int,Int})
            # e{m,n} => eee...eee|eee...ee|...|eee...e  (greedy)
            #           |<- n ->|              |<-m->|
            # e{m,n} => eee...e|eee...ee|...|eee...eee  (lazy)
            #           |<-m->|              |<- n ->|
            m, n = rng
            @assert m ≤ n
            if m == n
                head = :concat
                args = fill(pat, m)
            else
                head = :|
                f = (k, t) -> expr(:|, [expr(:concat, fill(pat, k)), t])
                if greedy
                    args = foldr(f, n:-1:m+1, init=expr(:concat, fill(pat, m))).args
                else
                    args = foldr(f, m:1:n-1, init=expr(:concat, fill(pat, n))).args
                end
            end
        else
            @assert false "invalid AST"
        end
    end
    i = 1
    for i in 1:lastindex(args)
        args[i] = desugar(T, args[i])
    end
    return expr(head, args)
end

function desugar(::Type{T}, atom) where {T}
    return atom
end

# Compiler
# --------

#  tag  | operation |           meaning
# ------|-----------|----------------------------
# 0b000 | match     | the pattern matches
# 0b001 | bits b    | matches b in bit-wise way
# 0b010 | jump l    | jump to l
# 0b011 | push l    | push l and go to next
# 0b100 | save i    | save string state to i
# 0b101 | head      | matches the head of string
# 0b110 | last      | matches the last of string
# 0b111 | fork l    | push next and go to l

primitive type Op 32 end

const MatchTag = UInt32(0b000) << 29
const BitsTag  = UInt32(0b001) << 29
const JumpTag  = UInt32(0b010) << 29
const PushTag  = UInt32(0b011) << 29
const SaveTag  = UInt32(0b100) << 29
const HeadTag  = UInt32(0b101) << 29
const LastTag  = UInt32(0b110) << 29
const ForkTag  = UInt32(0b111) << 29

# constructors
match() = reinterpret(Op, MatchTag)
bits(b::UInt32) = reinterpret(Op, BitsTag | b)
jump(l::Int) = reinterpret(Op, JumpTag | UInt32(l))
push(l::Int) = reinterpret(Op, PushTag | UInt32(l))
save(l::Int) = reinterpret(Op, SaveTag | UInt32(l))
head() = reinterpret(Op, HeadTag)
last() = reinterpret(Op, LastTag)
fork(l::Int) = reinterpret(Op, ForkTag | UInt32(l))

const operand_mask = (UInt32(1) << 29) - one(UInt32)

function tag(op::Op)
    return reinterpret(UInt32, op) & ~operand_mask
end

function operand(op::Op)
    return reinterpret(UInt32, op) &  operand_mask
end

function Base.show(io::IO, op::Op)
    t = tag(op)
    x = operand(op)
    if t == MatchTag
        print(io, "match")
    elseif t == BitsTag
        print(io, "bits ", @sprintf("0x%08x", x))
    elseif t == JumpTag
        print(io, "jump ", x)
    elseif t == PushTag
        print(io, "push ", x)
    elseif t == SaveTag
        print(io, "save ", x)
    elseif t == HeadTag
        print(io, "head")
    elseif t == LastTag
        print(io, "last")
    elseif t == ForkTag
        print(io, "fork ", x)
    else
        @assert false
    end
end

function print_program(prog::Vector{Op})
    L = lastindex(prog)
    for l in 1:L
        print(lpad(l, ndigits(L)), ": ", prog[l])
        if l != lastindex(prog)
            println()
        end
    end
end

function compile(tree::SyntaxTree)
    code = Op[]
    push!(code, save(1))
    compilerec!(code, tree, 2)
    push!(code, save(2))
    push!(code, match())
    return code
end

function compilerec!(code, tree::SyntaxTree, k)
    h = tree.head
    args = tree.args
    if h == :bits
        push!(code, bits(UInt32(args[1])))
    elseif h == :concat
        for arg in args
            k = compilerec!(code, arg, k)
        end
    elseif h == :*
        push!(code, push(0))  # placeholder
        l = length(code)
        k = compilerec!(code, args[1], k)
        push!(code, jump(l))
        code[l] = push(length(code) + 1)
    elseif h == Symbol("*?")
        push!(code, fork(0))  # placeholder
        l = length(code)
        k = compilerec!(code, args[1], k)
        push!(code, jump(l))
        code[l] = fork(length(code) + 1)
    elseif h == :|
        push!(code, push(0))  # placeholder
        l = length(code)
        k = compilerec!(code, args[1], k)
        push!(code, jump(0))  # placeholder
        code[l] = push(length(code) + 1)
        l = length(code)
        k = compilerec!(code, args[2], k)
        code[l] = jump(length(code) + 1)
    elseif h == :capture
        k′ = k
        push!(code, save(2k′-1))
        k = compilerec!(code, args[1], k + 1)
        push!(code, save(2k′))
    elseif h == :head
        push!(code, head())
    elseif h == :last
        push!(code, last())
    else
        @assert false "invalid tree"
    end
    return k
end

# Virtual machine
# ---------------

"""
Regular expression for `DNA`, `RNA`, and `AminoAcid`.
"""
struct Regex{T}
    pat::String       # regular expression pattern (for printing)
    code::Vector{Op}  # compiled code
    nsaves::Int       # the number of `save` operations in `code`

    function Regex{T}(pat::AbstractString, syntax=:pcre) where T
        if syntax == :pcre
            ast = desugar(T, parse(T, pat))
        elseif syntax == :prosite
            if T != BioSequences.AminoAcid
                throw(ArgumentError("alphabet must be AminoAcid for PROSITE syntax"))
            end
            ast = desugar(BioSequences.AminoAcid, parse_prosite(pat))
        else
            throw(ArgumentError("invalid syntax: $syntax"))
        end
        code = compile(ast)
        nsaves = 0
        for op in code
            nsaves += tag(op) == SaveTag
        end
        @assert iseven(nsaves)
        return new{T}(pat, code, nsaves)
    end
end

function Base.show(io::IO, re::Regex{T}) where {T}
    if T == BioSequences.DNA
        opt = "dna"
    elseif T == BioSequences.RNA
        opt = "rna"
    elseif T == BioSequences.AminoAcid
        opt = "aa"
    else
        assert(false)
    end
    print(io, "biore\"", re.pat, "\"", opt)
end

"""
Result of matching by `Regex`.
"""
struct RegexMatch{S}
    seq::S
    captured::Vector{Int}
end

function Base.show(io::IO, m::RegexMatch)
    print(io, "RegexMatch(")
    for k in 1:div(length(m.captured), 2)
        if k > 1
            print(io, ", ", k - 1, '=')
        end
        print(io, '"', m.seq[m.captured[2k-1]:m.captured[2k]-1], '"')
    end
    print(io, ')')
end

"""
    matched(match)

Return a matched pattern.
"""
function matched(m::RegexMatch{S}) where {S}
    return m.seq[m.captured[1]:m.captured[2]-1]
end

"""
    captured(match)

Retrun a vector of captured patterns.
"""
function captured(m::RegexMatch{S}) where {S}
    return [m.captured[2k-1] != 0 && m.captured[2k] != 0 ?
            m.seq[m.captured[2k-1]:m.captured[2k]-1] : nothing
            for k in 2:div(length(m.captured), 2)]
end

function checkeltype(re::Regex{T}, seq::BioSequences.BioSequence) where {T}
    if eltype(seq) != T
        throw(ArgumentError("element type of sequence doesn't match with regex"))
    end
end

function Base.match(re::Regex{T}, seq::BioSequences.BioSequence, start::Integer=1) where {T}
    checkeltype(re, seq)

    # use the first unambiguous symbol in the regular expression to find the
    # next starting position of pattern matching; this improves the performance
    if length(re.code) ≥ 2 && tag(re.code[2]) == BitsTag && count_ones(operand(re.code[2])) == 1
        firstsym = bits2sym(T, operand(re.code[2]))
    else
        firstsym = BioSequences.gap(T)
    end

    # a thread is `(<program counter>, <sequence's iterator state>)`
    threads = Stack{Tuple{Int,Int}}()
    captured = Vector{Int}(undef, re.nsaves)
    s = start
    while true
        if firstsym != BioSequences.gap(T)
            s = findnext(firstsym, seq, s)
            if s == nothing
                break
            end
        end
        empty!(threads)
        push!(threads, (1, s))
        fill!(captured, 0)
        if runmatch!(threads, captured, re, seq)
            return RegexMatch(seq, captured)
        end
        s += 1
        if s > lastindex(seq)
            break
        end
    end
    return nothing
end

function Base.findfirst(re::Regex{T}, seq::BioSequences.BioSequence, start::Integer=1) where {T}
    checkeltype(re, seq)
    m = Base.match(re, seq, start)
    if m === nothing
        return nothing
    else
        return m.captured[1]:m.captured[2]-1
    end
end

struct RegexMatchIterator{T,S}
    re::Regex{T}
    seq::S
    overlap::Bool

    function RegexMatchIterator{T,S}(re::Regex{T}, seq::S, overlap::Bool) where {T,S}
        checkeltype(re, seq)
        return new{T,S}(re, seq, overlap)
    end
end

function Base.IteratorSize(::Type{RegexMatchIterator{T,S}}) where {T,S}
    return Base.SizeUnknown()
end

function Base.eltype(::Type{RegexMatchIterator{T,S}}) where {T,S}
    return RegexMatch{S}
end

function Base.iterate(iter::RegexMatchIterator)
    threads = Stack{Tuple{Int,Int}}()
    captured = Vector{Int}(undef, iter.re.nsaves)
    s = 1
    push!(threads, (1, s))
    fill!(captured, 0)
    state = advance!(threads, captured, s, iter.re, iter.seq, iter.overlap)
    return iterate(iter, state)
end

function Base.iterate(iter::RegexMatchIterator, state)
    item, threads, captured, s = state
    if item === nothing
        return nothing
    else
        item, advance!(threads, captured, s, iter.re, iter.seq, iter.overlap)
    end
end

function advance!(threads, captured, s, re, seq, overlap)
    while true
        if runmatch!(threads, captured, re, seq)
            if !overlap
                empty!(threads)
                s = captured[2]
                if s <= lastindex(seq)
                    push!(threads, (1, s))
                end
            end
            return RegexMatch(seq, copy(captured)), threads, captured, s
        end
        if !isempty(seq)
            s += 1
        end
        if s > lastindex(seq)
            break
        end
        push!(threads, (1, s))
        fill!(captured, 0)
    end
    return nothing, threads, captured, s
end

function Base.eachmatch(re::Regex{T}, seq::BioSequences.BioSequence, overlap::Bool = true) where {T}
    checkeltype(re, seq)
    return RegexMatchIterator{T,typeof(seq)}(re, seq, overlap)
end

function Base.occursin(re::Regex{T}, seq::BioSequences.BioSequence) where {T}
    return Base.match(re, seq) !== nothing
end

# simple stack
mutable struct Stack{T}
    top::Int
    data::Vector{T}

    function Stack{T}(sz::Int=0) where T
        return new{T}(0, Vector{T}(undef, sz))
    end
end

function Base.isempty(stack::Stack)
    return stack.top == 0
end

@inline function Base.push!(stack::Stack{T}, x::T) where {T}
    if stack.top + 1 > length(stack.data)
        push!(stack.data, x)
    else
        stack.data[stack.top+1] = x
    end
    stack.top += 1
    return stack
end

@inline function Base.pop!(stack::Stack)
    item = stack.data[stack.top]
    stack.top -= 1
    return item
end

@inline function Base.empty!(stack::Stack)
    stack.top = 0
    return stack
end

# run pattern maching using the current `threads` stack. Captured positions are
# stored in the preallocated `captured`. If match is found, this returns `true`
# immediately; otherwise returns `false`.
function runmatch!(threads::Stack{Tuple{Int,Int}},
                   captured::Vector{Int},
                   re::Regex, seq::BioSequences.BioSequence)
    while !isempty(threads)
        pc::Int, s = pop!(threads)
        while true
            op = re.code[pc]
            t = tag(op)
            if t == BitsTag
                if s > lastindex(seq)
                    break
                end
                sym = seq[s]
                s += 1
                if BioSequences.compatbits(sym) & operand(op) != 0
                    pc += 1
                else
                    break
                end
            elseif t == JumpTag
                pc = operand(op)
            elseif t == PushTag
                push!(threads, (convert(Int, operand(op)), s))
                pc += 1
            elseif t == ForkTag
                push!(threads, (pc + 1, s))
                pc = operand(op)
            elseif t == SaveTag
                captured[operand(op)] = s
                pc += 1
            elseif t == HeadTag
                if s == 1
                    pc += 1
                else
                    break
                end
            elseif t == LastTag
                if s == lastindex(seq) + 1
                    pc += 1
                else
                    break
                end
            elseif t == MatchTag
                return true
            end
        end
    end
    return false
end

end  # module RE

# exported from BioSequences
const matched = RE.matched
const captured = RE.captured


text:
# see https://mycarta.wordpress.com/2013/02/21/perceptual-rainbow-palette-the-method/
# and https://mycarta.wordpress.com/2013/03/06/perceptual-rainbow-palette-the-goodies/

cubeYFrgb = [
    0.515100000000000 0.048200000000000 0.669700000000000
    0.515898208833707 0.056090405818443 0.678507004736174
    0.516602625686953 0.063937810746495 0.687238897518437
    0.517219112558518 0.071741790110892 0.695894349148579
    0.517753531447181 0.079501919238369 0.704472030428390
    0.518211744351720 0.087217773455659 0.712970612159658
    0.518599613270914 0.094888928089498 0.721388765144175
    0.518923000203542 0.102514958466620 0.729725160183728
    0.519187767148382 0.110095439913759 0.737978468080108
    0.519399776104213 0.117629947757650 0.746147359635104
    0.519564889069815 0.125118057325028 0.754230505650506
    0.519688968043965 0.132559343942626 0.762226576928104
    0.519777875025443 0.139953382937181 0.770134244269685
    0.519837472013027 0.147299749635426 0.777952178477042
    0.519873621005496 0.154598019364096 0.785679050351962
    0.519892184001628 0.161847767449925 0.793313530696235
    0.519899023000204 0.169048569219648 0.800854290311651
    0.519900000000000 0.176200000000000 0.808300000000000
    0.519719353197973 0.183322319854640 0.815720757279215
    0.519195325357450 0.190431820903915 0.823172569049792
    0.518354785326770 0.197522126559592 0.830632970187359
    0.517224601954270 0.204586860233435 0.838079495567547
    0.515831644088289 0.211619645337211 0.845489680065986
    0.514202780577165 0.218614105282685 0.852841058558304
    0.512364880269235 0.225563863481624 0.860111165920133
    0.510344812012837 0.232462543345792 0.867277537027101
    0.508169444656310 0.239303768286956 0.874317706754838
    0.505865647047992 0.246081161716881 0.881209209978974
    0.503460288036220 0.252788347047333 0.887929581575140
    0.500980236469334 0.259418947690079 0.894456356418963
    0.498452361195669 0.265966587056883 0.900767069386075
    0.495903531063566 0.272424888559511 0.906839255352105
    0.493360614921361 0.278787475609729 0.912650449192683
    0.490850481617393 0.285047971619304 0.918178185783438
    0.488400000000000 0.291200000000000 0.923400000000000
    0.485834103257308 0.297234005081926 0.928568605531681
    0.482983048947002 0.303152030445487 0.933909879830708
    0.479881692471761 0.308965098377897 0.939361103417085
    0.476564889234263 0.314684231166368 0.944859556810815
    0.473067494637187 0.320320451098114 0.950342520531901
    0.469424364083211 0.325884780460347 0.955747275100348
    0.465670352975012 0.331388241540280 0.961011101036158
    0.461840316715269 0.336841856625128 0.966071278859334
    0.457969110706660 0.342256648002101 0.970865089089881
    0.454091590351863 0.347643637958414 0.975329812247802
    0.450242611053557 0.353013848781280 0.979402728853099
    0.446457028214420 0.358378302757910 0.983021119425778
    0.442769697237129 0.363748022175520 0.986122264485840
    0.439215473524364 0.369134029321320 0.988643444553290
    0.435829212478802 0.374547346482525 0.990521940148131
    0.432645769503121 0.379998995946347 0.991695031790367
    0.429700000000000 0.385500000000000 0.992100000000000
    0.426955627689028 0.391040776795806 0.992033430311701
    0.424341857441582 0.396601282269409 0.991837733911093
    0.421841766388028 0.402176845383000 0.991518929794611
    0.419438431658726 0.407762795098774 0.991083036958688
    0.417114930384041 0.413354460378921 0.990536074399758
    0.414854339694334 0.418947170185636 0.989884061114255
    0.412639736719970 0.424536253481110 0.989133016098613
    0.410454198591311 0.430117039227536 0.988288958349266
    0.408280802438720 0.435684856387107 0.987357906862647
    0.406102625392560 0.441235033922016 0.986345880635191
    0.403902744583195 0.446762900794454 0.985258898663332
    0.401664237140986 0.452263785966615 0.984102979943502
    0.399370180196297 0.457733018400691 0.982884143472138
    0.397003650879491 0.463165927058875 0.981608408245671
    0.394547726320931 0.468557840903360 0.980281793260537
    0.391985483650979 0.473904088896337 0.978910317513168
    0.389300000000000 0.479200000000000 0.977500000000000
    0.386485815291723 0.484459106596542 0.975760221040460
    0.383556495925275 0.489697308284900 0.973433158841523
    0.380523658744115 0.494913484852387 0.970574029230938
    0.377398920591702 0.500106516086313 0.967238048036456
    0.374193898311494 0.505275281773991 0.963480431085827
    0.370920208746951 0.510418661702730 0.959356394206802
    0.367589468741531 0.515535535659843 0.954921153227132
    0.364213295138692 0.520624783432641 0.950229923974566
    0.360803304781893 0.525685284808435 0.945337922276855
    0.357371114514594 0.530715919574536 0.940300363961749
    0.353928341180252 0.535715567518256 0.935172464857000
    0.350486601622327 0.540683108426906 0.930009440790356
    0.347057512684277 0.545617422087797 0.924866507589570
    0.343652691209560 0.550517388288241 0.919798881082390
    0.340283754041636 0.555381886815548 0.914861777096569
    0.336962318023963 0.560209797457031 0.910110411459855
    0.333700000000000 0.565000000000000 0.905600000000000
    0.330482489796562 0.569753771265260 0.901247329027488
    0.327285486357999 0.574473796931106 0.896925424471454
    0.324105620204068 0.579161070071084 0.892623861587620
    0.320939521854528 0.583816583758738 0.888332215631710
    0.317783821829136 0.588441331067614 0.884040061859447
    0.314635150647651 0.593036305071257 0.879736975526553
    0.311490138829831 0.597602498843214 0.875412531888754
    0.308345416895434 0.602140905457029 0.871056306201771
    0.305197615364218 0.606652517986248 0.866657873721327
    0.302043364755942 0.611138329504417 0.862206809703147
    0.298879295590364 0.615599333085080 0.857692689402953
    0.295702038387241 0.620036521801784 0.853105088076468
    0.292508223666332 0.624450888728074 0.848433580979416
    0.289294481947396 0.628843426937495 0.843667743367520
    0.286057443750189 0.633215129503593 0.838797150496503
    0.282793739594471 0.637566989499912 0.833811377622089
    0.279500000000000 0.641900000000000 0.828700000000000
    0.275984612289327 0.646244368674252 0.823419276244787
    0.272107656090870 0.650621840812154 0.817945310154802
    0.267943300383189 0.655020715915885 0.812297444761467
    0.263565714144843 0.659429293487622 0.806495023096205
    0.259049066354389 0.663835873029545 0.800557388190440
    0.254467525990387 0.668228754043831 0.794503883075594
    0.249895262031396 0.672596236032659 0.788353850783091
    0.245406443455974 0.676926618498208 0.782126634344354
    0.241075239242681 0.681208200942655 0.775841576790806
    0.236975818370074 0.685429282868180 0.769518021153870
    0.233182349816713 0.689578163776960 0.763175310464969
    0.229769002561157 0.693643143171175 0.756832787755527
    0.226809945581964 0.697612520553001 0.750509796056966
    0.224379347857693 0.701474595424619 0.744225678400710
    0.222551378366903 0.705217667288206 0.737999777818181
    0.221400206088152 0.708830035645940 0.731851437340803
    0.221000000000000 0.712300000000000 0.725800000000000
    0.221158864597120 0.715631854287593 0.719835177732689
    0.221618060578300 0.718845821346623 0.713926672306033
    0.222351491228271 0.721954094030321 0.708063811496819
    0.223333059831763 0.724968865191914 0.702235923081833
    0.224536669673507 0.727902327684631 0.696432334837864
    0.225936224038233 0.730766674361700 0.690642374541697
    0.227505626210672 0.733574098076348 0.684855369970121
    0.229218779475555 0.736336791681806 0.679060648899922
    0.231049587117612 0.739066948031301 0.673247539107888
    0.232971952421574 0.741776759978061 0.667405368370806
    0.234959778672170 0.744478420375315 0.661523464465462
    0.236986969154133 0.747184122076292 0.655591155168645
    0.239027427152191 0.749906057934218 0.649597768257141
    0.241055055951077 0.752656420802324 0.643532631507737
    0.243043758835520 0.755447403533838 0.637385072697221
    0.244967439090251 0.758291198981987 0.631144419602380
    0.246800000000000 0.761200000000000 0.624800000000000
    0.248620409474781 0.764181339323832 0.618327943764038
    0.250519929338712 0.767226691237550 0.611721670549682
    0.252490903747682 0.770324160916551 0.604996442805570
    0.254525676857579 0.773461853536234 0.598167522980336
    0.256616592824292 0.776627874271995 0.591250173522619
    0.258755995803710 0.779810328299233 0.584259656881055
    0.260936229951722 0.782997320793345 0.577211235504280
    0.263149639424216 0.786176956929728 0.570120171840931
    0.265388568377081 0.789337341883780 0.563001728339645
    0.267645360966205 0.792466580830900 0.555871167449058
    0.269912361347478 0.795552778946483 0.548743751617806
    0.272181913676788 0.798584041405929 0.541634743294527
    0.274446362110024 0.801548473384634 0.534559404927857
    0.276698050803075 0.804434180057995 0.527532998966433
    0.278929323911828 0.807229266601412 0.520570787858891
    0.281132525592174 0.809921838190281 0.513688034053868
    0.283300000000000 0.812500000000000 0.506900000000000
    0.285441030391755 0.814970358148540 0.500173127924337
    0.287571511027906 0.817353839911381 0.493464278901938
    0.289693702570368 0.819660032200184 0.486774343653453
    0.291809865681056 0.821898521926612 0.480104212899533
    0.293922261021885 0.824078896002326 0.473454777360828
    0.296033149254770 0.826210741338987 0.466826927757990
    0.298144791041628 0.828303644848257 0.460221554811668
    0.300259447044372 0.830367193441798 0.453639549242513
    0.302379377924918 0.832410974031271 0.447081801771176
    0.304506844345182 0.834444573528338 0.440549203118308
    0.306644106967079 0.836477578844660 0.434042644004559
    0.308793426452524 0.838519576891898 0.427563015150580
    0.310957063463431 0.840580154581716 0.421111207277021
    0.313137278661718 0.842668898825773 0.414688111104533
    0.315336332709298 0.844795396535732 0.408294617353766
    0.317556486268087 0.846969234623253 0.401931616745372
    0.319800000000000 0.849200000000000 0.395600000000000
    0.322015737312730 0.851508952431154 0.388971277701778
    0.324162879421189 0.853901005610339 0.381797811217510
    0.326259163698117 0.856361227341876 0.374201767441860
    0.328322327516257 0.858874685430085 0.366305313269494
    0.330370108248352 0.861426447679288 0.358230615595075
    0.332420243267142 0.864001581893804 0.350099841313269
    0.334490469945371 0.866585155877956 0.342035157318741
    0.336598525655780 0.869162237436064 0.334158730506156
    0.338762147771111 0.871717894372449 0.326592727770178
    0.340999073664106 0.874237194491431 0.319459316005472
    0.343327040707508 0.876705205597332 0.312880662106703
    0.345763786274057 0.879106995494473 0.306978932968536
    0.348327047736497 0.881427631987174 0.301876295485636
    0.351034562467569 0.883652182879756 0.297694916552668
    0.353904067840016 0.885765715976541 0.294556963064295
    0.356953301226579 0.887753299081848 0.292584601915185
    0.360200000000000 0.889600000000000 0.291900000000000
    0.363752939278655 0.891356573845279 0.291951234025230
    0.367692949971565 0.893082120648976 0.292100691608455
    0.371995853137442 0.894774225590923 0.292342006010974
    0.376637469834999 0.896430473850951 0.292668810494086
    0.381593621122946 0.898048450608892 0.293074738319093
    0.386840128059997 0.899625741044576 0.293553422747292
    0.392352811704864 0.901159930337836 0.294098497039985
    0.398107493116257 0.902648603668503 0.294703594458471
    0.404079993352890 0.904089346216408 0.295362348264050
    0.410246133473474 0.905479743161382 0.296068391718021
    0.416581734536721 0.906817379683258 0.296815358081685
    0.423062617601343 0.908099840961867 0.297596880616341
    0.429664603726052 0.909324712177039 0.298406592583290
    0.436363513969560 0.910489578508608 0.299238127243830
    0.443135169390580 0.911592025136403 0.300085117859262
    0.449955391047822 0.912629637240256 0.300941197690885
    0.456800000000000 0.913600000000000 0.301800000000000
    0.463933237700690 0.914547747009667 0.302782114404216
    0.471590354506306 0.915512698672784 0.303992700421243
    0.479703684829818 0.916485220383418 0.305402246162532
    0.488205563084191 0.917455677535638 0.306981239739531
    0.497028323682396 0.918414435523513 0.308700169263692
    0.506104301037399 0.919351859741111 0.310529522846465
    0.515365829562169 0.920258315582500 0.312439788599298
    0.524745243669674 0.921124168441750 0.314401454633643
    0.534174877772882 0.921939783712928 0.316385009060949
    0.543587066284760 0.922695526790103 0.318360939992666
    0.552914143618277 0.923381763067344 0.320299735540244
    0.562088444186400 0.923988857938719 0.322171883815134
    0.571042302402099 0.924507176798296 0.323947872928785
    0.579708052678340 0.924927085040144 0.325598190992647
    0.588018029428092 0.925238948058332 0.327093326118170
    0.595904567064322 0.925433131246928 0.328403766416804
    0.603300000000000 0.925500000000000 0.329500000000000
    0.610324263437393 0.925500000000000 0.330448063816230
    0.617152714711127 0.925500000000000 0.331337369435639
    0.623800821432705 0.925500000000000 0.332173466902187
    0.630284051213629 0.925500000000000 0.332961906259830
    0.636617871665404 0.925500000000000 0.333708237552528
    0.642817750399530 0.925500000000000 0.334418010824239
    0.648899155027512 0.925500000000000 0.335096776118921
    0.654877553160852 0.925500000000000 0.335750083480533
    0.660768412411053 0.925500000000000 0.336383482953033
    0.666587200389617 0.925500000000000 0.337002524580379
    0.672349384708049 0.925500000000000 0.337612758406530
    0.678070432977850 0.925500000000000 0.338219734475444
    0.683765812810524 0.925500000000000 0.338829002831079
    0.689450991817573 0.925500000000000 0.339446113517394
    0.695141437610500 0.925500000000000 0.340076616578347
    0.700852617800808 0.925500000000000 0.340726062057896
    0.706600000000000 0.925500000000000 0.341400000000000
    0.712355184459726 0.925500000000000 0.342087365110182
    0.718079256852663 0.925500000000000 0.342773389281185
    0.723771912921851 0.925500000000000 0.343458068337796
    0.729432848410332 0.925500000000000 0.344141398104801
    0.735061759061146 0.925500000000000 0.344823374406989
    0.740658340617336 0.925500000000000 0.345503993069147
    0.746222288821943 0.925500000000000 0.346183249916061
    0.751753299418008 0.925500000000000 0.346861140772519
    0.757251068148571 0.925500000000000 0.347537661463308
    0.762715290756675 0.925500000000000 0.348212807813215
    0.768145662985361 0.925500000000000 0.348886575647028
    0.773541880577669 0.925500000000000 0.349558960789533
    0.778903639276642 0.925500000000000 0.350229959065518
    0.784230634825320 0.925500000000000 0.350899566299770
    0.789522562966745 0.925500000000000 0.351567778317076
    0.794779119443958 0.925500000000000 0.352234590942224
    0.800000000000000 0.925500000000000 0.352900000000000
]

cmap = []
for i = 1:256
    push!(cmap, RGB(cubeYFrgb[i, 1], cubeYFrgb[i, 2], cubeYFrgb[i, 3]))
end

function cubeYF()
    return cmap
end

ext: jl

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
# remove deps.jl if it exists, in case build.jl fails
isfile("deps.jl") && rm("deps.jl")

@static if is_unix() begin
    const LIBPARDISONAMES = [
    "libpardiso500-GNU461-X86-64",
    "libpardiso500-GNU472-X86-64",
    "libpardiso500-GNU481-X86-64",
    "libpardiso"
]
end
end

@static if is_windows() begin
    const LIBPARDISONAMES = ["libpardiso500-WIN-X86-64.dll", "libpardiso"]
end
end

const PATH_PREFIXES = [
    dirname(@__FILE__),
    ""
]

# print to stderr, since that is where Pkg prints its messages
eprintln(x...) = println(STDERR, x...)

function find_paradisolib()
    found_lib = false
    for prefix in PATH_PREFIXES
        for libname in LIBPARDISONAMES
            try
                path = joinpath(prefix, libname)
                Libdl.dlopen(path, Libdl.RTLD_GLOBAL)
                global PARDISO_LIB_FOUND = true
                eprintln("found libpardiso at $path, using it")
                return path, true
            end
        end
    end
    eprintln("did not find libpardiso, assuming PARDISO 5.0 is not installed")
    return "", false
end

function find_mklparadiso()
    if haskey(ENV, "MKLROOT")
        eprintln("found MKLROOT key, using it")
        return ENV["MKLROOT"], true
    end
    eprintln("did not find MKLROOT key, assuming MKL is not installed")
    return "", false
end

pardisopath, found_pardisolib = find_paradisolib()
mklroot, found_mklpardiso = find_mklparadiso()

if !(found_mklpardiso || found_pardisolib)
    warn("no Pardiso library managed to load")
end

open("deps.jl", "w") do f
    print(f,
"""
const MKL_PARDISO_LIB_FOUND = $found_mklpardiso
const PARDISO_LIB_FOUND = $found_pardisolib
const MKLROOT = R"$mklroot"
const PARDISO_PATH = R"$pardisopath"
"""
)

end


text:
@doc raw"""
The Newton form for an interpolating polynomial represents ``f(x)`` at
points ``\tau``` by a sum.

``\sum_{i=0}^n f[x_0,\ldots,x_i]\prod_{j=0}^{i-1}(x-x_j)

Here, the `a` are the divided differences and ``\tau`` are the `x[j]`.
"""
struct NewtonForm{T}
    τ::Vector
    a::Vector{T}
end

"""
Construct a NewtonForm. From Conte and deBoor 1980, Eqn. 4.8.
"""
function NewtonForm(τ::Vector, f::Vector{T}) where {T <: Real}
    k = len(τ)
    a = zeros(T)
    dd = DividedDifference{T}()
    for j = 1:k
        a[j] = divided_difference(1, j, τ, f, dd)
    end
    NewtonForm{T}(τ, a)
end

function (nf::NewtonForm{T})(x) where {T}
    k = len(nf.a)
    total = zero(T)
    for j = k:-1:2
        total += a[j]
        total *= (x - τ[j - 1])
    end
    total += a[1]
    total
end

"""
Construct a polynomial that interpolates a function, using polynomials with zeroes.
This should match the constructor above. It comes from deBoor equation 4.2.
"""
function single_polynomial(τ::Vector{X}, f::Vector{T}) where {T <: Real, X<: Real}
    a = zeros(T, length(f))
    for i = 1:length(f)
        total = one(T)
        for l1 = 1:(i - 1)
            total *= 1 / (τ[l1] - τ[i])
        end
        for l2 = (i + 2):length(f)
            total *= 1 / (τ[l2] - τ[i])
        end
        a[i] = f[i] * total
    end
    NewtonForm{T}(τ, a)
end


text:
"""
    recoverableProblem(C, c, X, α, dg)

Solves recoverable problem ``REC(c) = min_{x \\in X} min_{y \\in X^{\\alpha}_x} (Cx + cy)``.

Check section 4 _Solving the problems by MIP formulations_ of [publication](https://arxiv.org/abs/1811.06719)
for more information about this algorithm.

# Arguments
- `C`: is a vector of nonnegative first stage costs.
- `c`: is a vector of a nonnegative nominal second stage costs.
- `X`: is a set of feasible solutions represented as a list functions, each of which
    accepts a list of JuMP variables as an argument and returns a JuMP linear constraint.
- `α`: fixed number belonging to ``[0, 1]``
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
function recoverableProblem(C, c, X, α, dg)
    @assert size(C) == size(c)

    n = size(C, 1)

    model = Model(solver=CplexSolver(CPXPARAM_ScreenOutput = getProperty("recoverableProblem.cplexLogging")))
    if ndims(C) == 1
        @variable(model, x[1:n], Bin)
        @variable(model, y[1:n], Bin)
        if hasEqualCardinalityProperty(dg)
            @variable(model, z[1:n] >= 0)
        else
            @variable(model, z[1:n], Bin)
        end
    elseif ndims(C) == 2
        m = size(C, 2)
        @variable(model, x[1:n, 1:m], Bin)
        @variable(model, y[1:n, 1:m], Bin)
        if hasEqualCardinalityProperty(dg)
            @variable(model, z[1:n, 1:m] >= 0)
        else
            @variable(model, z[1:n, 1:m], Bin)
        end
    else
        throw(ArgumentError("Unsupported number of dimensions ($(ndims(C)))"))
    end

    @objective(model, Min, vecdot(C, x) + vecdot(c, y))

    if hasEqualCardinalityProperty(dg)
        l = ceil(n * (1 - α))
        @constraint(model, sum(z) >= l)
    else
        @constraint(model, sum(x-z) <= α * sum(x))
        @constraint(model, z .>= x + y - 1)
    end

    @constraint(model, z .<= x)
    @constraint(model, z .<= y)

    # Constraints defining set of feasible solutions
    for constraint in X
        JuMP.addconstraint(model, constraint(x))
    end

    for constraint in X
        JuMP.addconstraint(model, constraint(y))
    end

    status = solve(model)
    (getvalue(x), getvalue(y), getobjectivevalue(model))
end


text:
using XGBoost, MCMCChains, StanSample

"""
    rstar(chains::Chains; subset = 0.8, niter = 1_000, eta = 0.5, XGBoostParams)
    rstar(chains::Chains, iterations::Int; subset = 0.8, niter = 1_000, eta = 0.5, XGBoostParams)
    rstar(x::AbstractMatrix, y::AbstractVector, nchains::Int, iterations::Int; subset = 0.8, niter = 1_000, eta = 0.5, XGBoostParams)

Compute the R* statistic for convergence diagnostic of MCMC. This implementation is an adaption of Algorithm 1 & 2, described in [Lambert & Vehtari]. Note that the correctness of the statistic depends on the convergence of the classifier used internally in the statistic. You can track if the training of the classifier converged by inspection of the printed RMSE values from the XGBoost backend. To adjust the number of iterations used to train the classifier set `niter` accordingly.

# Usage

```julia-repl
using XGBoost
# You need to load XGBoost before using MCMCChains.Rstar

...

chn = ...

# Compute R⋆ using defaults settings for the  gradient boosting classifier used to compute the statistic.
# This is the recomended use.
R = rstar(chn)

# Compute 100 samples of the R⋆ statistic using sampling from according to the prediction probabilities.
# This approach can be slow and results in a less accurate estimation of the R* statistic.
# See discussion in Section 3.1.3 in the paper.
Rs = rstar(chn, 100)

# estimate Rstar
R = mean(Rs)

# visualize distribution
histogram(Rs)
```

## References:
[Lambert & Vehtari] Ben Lambert and Aki Vehtari. "R∗: A robust MCMC convergence diagnostic with uncertainty using gradient-boostined machines." Arxiv 2020.
"""
function rstar(x::AbstractMatrix, y::AbstractVector, nchains::Int, iterations::Int; subset = 0.8, niter = 1_000, eta = 0.5, xgboostparams...)

    N = length(y)

    # randomly sub-select training and testing set
    Ntrain = round(Int, N*subset)
    Ntest = N - Ntrain

    ids = Random.randperm(N)
    train_ids = view(ids, 1:Ntrain)
    test_ids = view(ids, (Ntrain+1):N)

    @assert length(test_ids) == Ntest

    # use predicted probabilities?
    mode = iterations > 1 ? "multi:softprob" : "multi:softmax"

    @info "Training classifier"
    # train classifier using XGBoost
    classif = XGBoost.xgboost(x[train_ids,:], niter; label = y[train_ids],
                      objective = mode, num_class = nchains,
                      xgboostparams...)

    @info "Computing R* statistics"
    Rstats = ones(iterations) * Inf
    for i in 1:iterations
        # predict labels for "test" data
        p = XGBoost.predict(classif, x[test_ids,:])

        pred = if length(p) == Ntest*nchains
            probs = reshape(p, Ntest, nchains)
            map(s -> rand(Categorical(s / sum(s))), eachrow(probs))
        else
            p
        end

        # compute statistic
        a = mean(pred .== y[test_ids])
        Rstats[i] = nchains*a
    end

    return Rstats
end

function rstar(chn::Chains, iterations; kwargs...)
    nchains = size(chn, 3)
    @assert nchains > 1

    # collect data
    x = mapreduce(c -> Array(chn[:,:,c]), vcat, chains(chn))
    y = mapreduce(c -> ones(Int, length(chn))*c, vcat, chains(chn)) .- 1

    return rstar(x, y, nchains, iterations, kwargs...)
end

rstar(chn::Chains; kwargs...) = first(rstar(chn, 1; kwargs...))


text:
"""
    create_replicator_eqs(game::Dict{String,<:AbstractArray}, s::Array{SymPy.Sym,1})

Creates differential equations for replicator

**Input parameters**
* `game::Dict{String,<:AbstractArray}` - dictionary of players and their payoff matrices
* `s::Array{SymPy.Sym,1}` - collection of symbolic actions probabilities
"""
function create_replicator_eqs(game::Dict{String,<:AbstractArray},
               s::Array{SymPy.Sym,1})
    n = size(game["player1"])[1]
    id_m = Matrix(1I, n, n)
    res = []
    for i in 1:n
        append!(res, sympy.simplify(s[i]*(get_payoff(game, [map(x->sympify(x), id_m[i, :]), s])["player1"] - get_payoff(game, [s, s])["player1"])))
    end
    return res
end


text:
export interval

"""
    interval(period::Union{Int, Nothing} = nothing)

Creation operator for the `TimerObservable`. `interval` returns an Observable that emits an infinite sequence of ascending integers,
with a constant interval of time of your choosing between those emissions.
The first emission is not sent immediately, but only after the first period has passed.

# Arguments
- `interval`: the interval size in milliseconds

# Examples

```
using Rocket

source = interval(50)

subscription = subscribe!(source, logger())
sleep(0.215)
unsubscribe!(subscription)
sleep(0.215)
subscription = subscribe!(source, logger())
sleep(0.185)
unsubscribe!(subscription)

close(source)
;

# output

[LogActor] Data: 0
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 8
[LogActor] Data: 9
[LogActor] Data: 10
```

See also: [`timer`](@ref), [`TimerObservable`](@ref), [`Subscribable`](@ref)
"""
interval(period::Int) = timer(period, period)


text:
"""
    IntervalConstantSolver(; θ = 0.5, timestep)

BTPDE solver specialized on intervalwise constant `ScalarGradient`s, e.g [`PGSE`](@ref),
[`DoublePGSE`](@ref). Raises an error for other gradients.
"""
Base.@kwdef struct IntervalConstantSolver{T}
    θ::T = 0.5
    timestep::T
end


text:
using Markdown
# using HorizonSideRobots
#########################################################
md"""
    #Базовые функции, связанные с передвижением робота
"""
#########################################################
"""
Делает попытку шагнуть в сторону {side} и возвращает {false} при неудачной попытке, иначе {true}
"""
function carefull_move!(r::Robot, side::HorizonSide)
    if !isborder(r, side)
        move!(r, side)
        return true
    else
        return false
    end
end

"""
Двигает робота до перегородки в сторону {side} и возвращает кол-во шагов {steps}
"""
function move_till_border!(r, side::HorizonSide)
    steps = 0
    while isborder(r, side) == 0
        move!(r, side)
        steps += 1
    end
    return steps
end

"""
Двигает робота на {x} клеток в сторону {side}
"""
function move_for!(r, side::HorizonSide, x::Int)
    for _ in 1:x
        if !isborder(r, side)
            move!(r, side)
        else
            return false
        end
    end
    return true
end

"""
Перемещает робота в угол ({side_1}, {side_2}) и возвращает проделанный при этом путь {path}
"""
function move_to_corner!(r, side_1::HorizonSide, side_2::HorizonSide,)
    path = []
    while !isborder(r, side_1) || !isborder(r, side_2)
        push!(path, move_till_border!(r, side_1))
        push!(path, move_till_border!(r, side_2))
    end
    return path
end

"""
Перемещает робота из угла в стартовую позицию по {back_path}
"""
function move_back_from_corner!(r, side_1, side_2, back_path)
    side_1 = opposite_side(side_1)
    side_2 = opposite_side(side_2)
    for i in 1:2:size(back_path)[1]
        move_for!(r, side_2, back_path[i])
        move_for!(r, side_1, back_path[i+1])
    end
end

"""
Двигает робота на клетку (x, y)
"""
function move_to!(r::Robot, x, y)
    move_till_border!(r, West)
    move_till_border!(r, Sud)
    for _ in 1:y-1
        move!(r, Nord)
    end
    for _ in 1:x-1
        move!(r, Ost)
    end
end

"""Возвращает противоположную {side} сторону."""
opposite_side(side::HorizonSide) = HorizonSide((Int(side) + 2) % 4)

"""Возвращает следующую сторону от {side} по часовой стрелке"""
next_side(side::HorizonSide) = opposite_side(HorizonSide((Int(side) + 1) % 4))

"""Возвращает следующую сторону от {side} против часовой стрелке"""
previous_side(side::HorizonSide) = HorizonSide((Int(side) + 1) % 4)

"""Возвращает сторону, в которой нет перегородки"""
function get_no_border_side(r)
    for side in [Nord, Ost, Sud, West]
        if !isborder(r, side)
            return side
        end
    end
    return false
end

"""Возвращает координаты клетки сверху/снизу/слева/справа от переданной"""
function side_to_i_j(side::HorizonSide, cell)
    i, j = cell
    if side == Nord
        return (i - 1, j)
    elseif side == Sud
        return (i + 1, j)
    elseif side == West
        return (i, j - 1)
    elseif side == Ost
        return (i, j + 1)
    end
end

"""Поиск в глубину, робот проходится по всем клеткам и маркирует переданные в {cells}"""
function mark_some_cells_with_dfs!(r, visited, cell, side, cells)
    is_not_at_start = side != false

    if is_not_at_start
        cheker = isborder(r, side) == false
    else
        cheker = true
    end

    if (cell in visited) == false && cheker && length(cells) > 0
        if is_not_at_start
            move!(r, side)
        end
        push!(visited, cell)

        if cell in cells
            putmarker!(r)
            println("!Marker at ", cell)
            filter!(el -> el ≠ cell, cells)
        end

        for tempSide in [Nord, Ost, Sud, West]
            if mark_some_cells_with_dfs!(r, visited, side_to_i_j(tempSide, cell), tempSide, cells)
                move!(r, opposite_side(tempSide))
            end
        end

        return true

    else
        return false
    end
end

"""Проверяет {x}-овую клетку в сторону {side} от места старта робота на наличие выхода в стороне {border_side}"""
function find_or_move!(r, side, x, border_side)
    move_for!(r, side, x)
    if isborder(r, border_side) == false
        return true
    end
    move_for!(r, opposite_side(side), x)
    return false
end

"""Перемещает робота к выходу со стороны {border_side}"""
function escape_find(r, border_side)
    count = 1
    sides = (next_side(border_side), opposite_side(next_side(border_side)))
    while true                                  # Робот двигается поочередно влево-вправо (с каждой итерацией проверяя на одну клетку с каждой стороны больше), пока не наткнется на выход
        for side in sides
            if find_or_move!(r, side, count, border_side)
                return (count, side)
            end
        end
        count += 1
    end
    return false
end

"""Обходит любую одинарную прямоугольную перегородку в стороне {side}"""
function try_move!(r::Robot, side::HorizonSide)
    if isborder(r, side)
        n = 0;
        ort_side = next_side(side)
        while isborder(r, side)
            if !isborder(r, ort_side)
                move!(r, ort_side)
                n+=1
            else
                move_for!(r, opposite_side(ort_side), n)
                return false
            end
        end

        move!(r, side)
        while isborder(r, opposite_side(ort_side))
            move!(r, side)
        end
        move_for!(r, opposite_side(ort_side), n)
        return true
    else
        move!(r, side)
        return true
    end
end

"""
Обходит любую одинарную прямоугольную перегородку в стороне {side}
Возвращает:
        {-1} - преград нет
         {0} - неудача
         {1} - перегородка линейна
         {2} - перегородка прямоугольна
"""
function try_move_extra!(r::Robot, side::HorizonSide)
    result = 0
    if isborder(r, side)
        n = 0;
        ort_side = next_side(side)
        while isborder(r, side)
            if !isborder(r, ort_side)
                move!(r, ort_side)
                n+=1
            else
                move_for!(r, opposite_side(ort_side), n)
                return result
            end
        end

        move!(r, side)
        while isborder(r, opposite_side(ort_side))
            result = 1
            move!(r, side)
        end
        move_for!(r, opposite_side(ort_side), n)
        result += 1
    else
        move!(r, side)
        result = -1
    end
    return result
end

"""Посещает все доступные клетки и возвращет {ans}, имеющий значение соответствующее {[i_min, i_max, j_min, j_max]}"""
function visit_all_cells!(r, side=Nord, cell=(0, 0), visited=Set(), ans=[0, 0, 0, 0,])

    if (cell in visited) == false && !isborder(r, side)
        move!(r, side)
        push!(visited, cell)

        ans[1] = min(cell[1], ans[1])
        ans[2] = max(cell[1], ans[2])
        ans[3] = min(cell[2], ans[3])
        ans[4] = max(cell[2], ans[4])

        for tempSide in [Nord, Ost, Sud, West]
            tempAns = visit_all_cells!(r, tempSide, side_to_i_j(tempSide, cell), visited, ans)
            if tempAns != false
                move!(r, opposite_side(tempSide))
                ans = tempAns
            end
        end
        return ans
    else
        return false
    end
end

"""Возвращает размер поля {(M, N)} и текущую позицию робота {r} {(i, j)}"""
function get_field_size_and_position!(r)
    side = get_no_border_side(r)
    if side == false
        return (1, 1), (1, 1)
    end
    ans = visit_all_cells!(r, side, side_to_i_j(side, (0,0)))
    move!(r, opposite_side(side))
    M = ans[2]-ans[1]+1
    N = ans[4]-ans[3]+1
    return (M, N), (M-ans[2], N-ans[4])
end

#########################################################
md"""
    #Базовые функции, связанные с базовыми типами данных, и сами типы и связанные структуры
"""
#########################################################

import HorizonSideRobots
"""Переопределение стандартных функций робота для стандартного робота {::Robot}"""
move!(robot::Robot, side::HorizonSide) = HorizonSideRobots.move!(robot, side)
isborder(robot::Robot,  side::HorizonSide) = HorizonSideRobots.isborder(robot, side)
putmarker!(robot::Robot) = HorizonSideRobots.putmarker!(robot)
ismarker(robot::Robot) = HorizonSideRobots.ismarker(robot)
temperature(robot::Robot) = HorizonSideRobots.temperature(robot)
show(robot::Robot) = HorizonSideRobots.show(robot)
show!(robot::Robot) = HorizonSideRobots.show!(robot)

#---------------------------------#

"""Структура координат {i, j}, где i - номер строки поля, а {j} - номер столбца поля (для робота)"""
mutable struct Coordinates
    i::Integer
    j::Integer

    Coordinates() = new(0,0)
    Coordinates(i::Integer,j::Integer) = new(i,j)
    Coordinates((i, j)::Tuple{Int64, Int64}) = new(i, j)
end

"""Изменяет координаты соответственно перемещению робота"""
function move!(coordinates::Coordinates, side::HorizonSide)
    if side==Nord
        coordinates.i -= 1
    elseif side==Sud
        coordinates.i += 1
    elseif side==Ost
        coordinates.j += 1
    else
        coordinates.j -= 1
    end
end

"""Обновляет значения координат соответственно на {i} и {j}"""
function update!(coordinates::Coordinates, (i, j)::Tuple{Int64, Int64})
    coordinates.i = i
    coordinates.j = j
end

#---------------------------------#

"""Базовый абстрактный робот"""
abstract type AbstractRobot end
"""Определение стандартных функций робота для робота {::AbstractRobot}"""
move!(robot::AbstractRobot, side::HorizonSide) = HorizonSideRobots.move!(get(robot), side)
isborder(robot::AbstractRobot,  side::HorizonSide) = HorizonSideRobots.isborder(get(robot), side)
putmarker!(robot::AbstractRobot) = HorizonSideRobots.putmarker!(get(robot))
ismarker(robot::AbstractRobot) = HorizonSideRobots.ismarker(get(robot))
temperature(robot::AbstractRobot) = HorizonSideRobots.temperature(get(robot))
show(robot::AbstractRobot) = HorizonSideRobots.show(get(robot))
show!(robot::AbstractRobot) = HorizonSideRobots.show!(get(robot))
"""Возвращает {nothing} для робота типа {AbstractRobot}"""
get(robot::AbstractRobot) = nothing

"""Двигает робота змейкой по всему полю"""
function move_as_snake!(r::AbstractRobot, side=Ost)
    do_movements!(r, side) #этаче, должно быть {do_movements} (которая еще и переопределяется для разных типов в зависимости от задач получается(?))
    while isborder(r, Nord)
        move!(r, Nord)
        side = opposite_side(side)
        do_movements!(r, side)
    end
end

#---------------------------------#

"""Абстрактный тип робота, обходящего перегородки"""
abstract type AbstractBorderRobot <: AbstractRobot end

"""Переопределенная {move_till_border!} для {AbstractBorderRobot}:
                        + Использует {try_move!}, что позволяет обходить одиночные перегородки по пути
"""
function move_till_border!(r::AbstractBorderRobot, side::HorizonSide)
    counter = 0

    while try_move!(get(r), side)
        counter += 1
    end

    return n
end

#---------------------------------#

"""Маркирующий робот от {AbstractBorderRobot}"""
struct PutMarkersBorderRobot <: AbstractBorderRobot
    robot::Robot
end

"""Возвращает ссылку на робота {::Robot} для робота типа {PutMarkersBorderRobot}"""
get(r::PutMarkersBorderRobot) = r.robot

"""
Переопределенная {try_move!}:
                        + Ставит маркеры
"""
function try_move!(r::PutMarkersBorderRobot, side::HorizonSide)
    try_move!(get(r), side)
    putmarker!(r)
end

"""Конструктор {PutMarkersBorderRobot}"""
function PutMarkersBorderRobot(r::Robot)
    putmarker!(r)
    new(r)
end

"""Маркирует все поле используя {move_as_snake!}"""
function mark_field(r::Robot)
    # утв: робот в ЮЗ углу
    r = PutMarkersBorderRobot(r)
    move_as_snake!(r)
end

#---------------------------------#

"""Тип робота, использующего глобальные координаты"""
struct SuperCheatRobot <: AbstractRobot
    robot::Robot
    position::Coordinates
    field_size::Tuple{Int64, Int64}

    SuperCheatRobot(r, (i, j), (M, N)) = new(r, Coordinates(i, j), (M, N))
    function SuperCheatRobot(r)
        (M, N), (i, j) = get_field_size_and_position!(r)
        return SuperCheatRobot(r, (i, j), (M, N))
    end
end

"""Возвращает ссылку на робота {::Robot} для робота типа {SuperCheatRobot}"""
get(robot::SuperCheatRobot) = robot.robot

"""Обновляет данные о текущей позиции {r.position} робота {r} (происходит вызов {get_field_size_and_position}!)"""
function update!(r::SuperCheatRobot)
    update!(r.position, get_field_size_and_position!(r.robot)[2])
end

"""Устанавливает данные о текущей позиции {r.position} робота {r} равными {(i, j)}"""
function setPositionTo!(r::SuperCheatRobot, (i, j))
    update!(r.position, (i,j))
end

"""
Переопределенная {move!(::Robot)}:
                        + Изменяет текущие координаты робота {r} согласно перемещению
"""
function move!(robot::SuperCheatRobot, side::HorizonSide)
    move!(robot.robot, side)
    move!(robot.position, side)
end

"""Возвращает текущую позицию робота {r.position} в виде {(i, j)}"""
position(r) = (r.position.i, r.position.j)


text:
# https://gist.github.com/silgon/0ba43e00e0749cdf4f8d244e67cd9d6a
# https://riptutorial.com/julia-lang/example/19443/parsing-json
# https://github.com/JuliaIO/JSON.jl

import Pkg;
Pkg.add("HTTP")
Pkg.add("JSON")
Pkg.add("JSON")
Pkg.add("JSON3")
Pkg.add("LazyJSON")

import HTTP
import JSON
import JSON;
import JSON3;
import LazyJSON;

using JSON
using HTTP
using JSON;
using JSON3;
using LazyJSON;

url = "https://raw.githubusercontent.com/xamarin/AndroidX/20200915-mono.cecil-fix/config.json"
response = HTTP.get(url)

response_body_as_string = response_body_as_stringing(response.body)
println("response_body_as_string = ", response_body_as_string)

json_object_JSON = JSON.parse(response_body_as_string)
println("json_object_JSON = ", json_object_JSON)

json_object_JSON3 = JSON3.read(response_body_as_string)
println("json_object_JSON3 = ", json_object_JSON)

json_object_LazyJSON = LazyJSON.value(response_body_as_string)
println("json_object_LazyJSON = ", json_object_LazyJSON)


text:
#! /usr/bin/env julia

using DifferentialEquations

b, c = 0.25, 5.0
pend2 = @ode_def Pendulum begin
  dθ = ω
  dω = (-b * ω) - (c * sin(θ))
end

prob = ODEProblem(pend, y0, (0.0, 10.0))
sol = solve(prob)   # solve on [0,10]
t, y = sol.t, hcat(sol.u...)'

using Winston
plot(t, y[:, 1], t, y[:, 2])
title("2D Differential Equation")
savefig("figures/Pendulum_solution.png")

max_stars_repo_name: pierre-haessig/slides

max_issues_repo_name: pierre-haessig/slides

max_forks_repo_name: pierre-haessig/slides


text:
function calc_nchk(n::Int,k::Int)
    #=
    Calculate n choose k
    =#
    if n == 0 || k == 0
		return 1
	end

    #@myassert(n>=k)
    accum::BigInt = 1
    for i in 1:k
        accum = accum * (n-k+i) ÷ i
    end
    return accum
end

function build_nchk_table(N)
	binom_coeff = Array{BigInt,2}(undef,N+1,N+1)
	for i in 0:N
    		for j in i:N
        		binom_coeff[j+1,i+1] = calc_nchk(j,i)
    		end
    	end
	return binom_coeff
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
#!/usr/bin/env julia

println("Hello world!")

max_stars_repo_name: smenjas/programming-languages-compared

max_issues_repo_name: smenjas/programming-languages-compared

max_forks_repo_name: smenjas/programming-languages-compared


text:
include("PivotFactorization.jl")
include("HierarchicalVector.jl")
include("HierarchicalMatrix.jl")
include("FormattedOperations.jl")
include("hierarchicalsolve.jl")


text:
######################################################
### This test implements example 13.5 in Greene (2012)
######################################################

# this is the data from Example C.1 in Greene 2012
income_data = [20.5, 31.5, 47.7, 26.2, 44, 8.28, 30.8, 17.2, 19.9, 9.96,
               55.8, 25.2, 29, 85.5, 15.1, 28.5, 21.4, 17.7, 6.42, 84.9]

function g_gamma(θ)
    P, λ = θ[1], θ[2]
    hcat(income_data - P / λ,
         income_data.^2 - P*(P+1) / (λ^2),
         log(income_data) + log(λ) - digamma(P),
         1./income_data - λ/(P - 1))
end

two_step = GMM.gmm(g_gamma, [2.40, 0.08], [0.0, 0.0], [Inf, Inf], eye(4);
                   mgr=GMM.TwoStepGMM())


text:
let
    counter = 0
    global clearcounter() = counter = 0
    global counter!() = (counter += 1; (counter > 100) && throw("Excessive Recursion Error"))
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
     module Simulate_FreeBodyRotation

Simulate body with unconstrained rotational motion in 3D (using quaternions; = index-1 model).
"""
module Simulate_FreeBodyRotation

include(joinpath("models", "FreeBodyRotation.jl"))
import .FreeBodyRotation
import ModiaMath

model  = FreeBodyRotation.Model()
result = ModiaMath.simulate!(model, stopTime=5.0, log=true)

ModiaMath.plot(result, [("Q[1]", "Q[2]", "Q[3]", "Q[4]"),
                        ("w[1]", "w[2]", "w[3]"),
                        ("der(w[1])", "der(w[2])", "der(w[3])")], heading="simulationWithoutMacro/Simulate_FreeBodyRotation.jl")

end


text:
# =====================================
# Projectile
#
# system type: affine continuous system
# state dimension: 4
# =====================================
using MathematicalSystems, SparseArrays

function projectile_model()
    # system matrix
    A = sparse([1, 3], [2, 4], [0.5, 0.7], 4, 4)

    # affine term
    b = sparsevec([4], [-9.81], 4)

    # affine continuous system
    S = @system(x' = Ax + b)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
@testset "minimal.atan2_test" begin

    @test atan(emptyinterval(), emptyinterval()) === emptyinterval()

    @test atan(emptyinterval(), entireinterval()) === emptyinterval()

    @test atan(entireinterval(), emptyinterval()) === emptyinterval()

    @test atan(interval(0.0, 0.0), interval(0.0, 0.0)) === emptyinterval()

    @test atan(entireinterval(), entireinterval()) === Interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)

    @test atan(interval(0.0, 0.0), interval(-Inf, 0.0)) === Interval(0x1.921FB54442D18p1, 0x1.921FB54442D19p1)

    @test atan(interval(0.0, 0.0), interval(0.0, Inf)) === Interval(0.0, 0.0)

    @test atan(interval(0.0, Inf), interval(0.0, 0.0)) === Interval(0x1.921FB54442D18p0, 0x1.921FB54442D19p0)

    @test atan(interval(-Inf, 0.0), interval(0.0, 0.0)) === Interval(-0x1.921FB54442D19p0, -0x1.921FB54442D18p0)

    @test atan(interval(-0x1p-1022, 0.0), interval(-0x1p-1022, -0x1p-1022)) === Interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)

    @test atan(interval(1.0, 1.0), interval(-1.0, -1.0)) === Interval(0x1.2D97C7F3321D2p1, 0x1.2D97C7F3321D3p1)

    @test atan(interval(1.0, 1.0), interval(1.0, 1.0)) === Interval(0x1.921FB54442D18p-1, 0x1.921FB54442D19p-1)

    @test atan(interval(-1.0, -1.0), interval(1.0, 1.0)) === Interval(-0x1.921FB54442D19p-1, -0x1.921FB54442D18p-1)

    @test atan(interval(-1.0, -1.0), interval(-1.0, -1.0)) === Interval(-0x1.2D97C7F3321D3p1, -0x1.2D97C7F3321D2p1)

    @test atan(interval(-0x1p-1022, 0x1p-1022), interval(-0x1p-1022, -0x1p-1022)) === Interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)

    @test atan(interval(-0x1p-1022, 0x1p-1022), interval(0x1p-1022, 0x1p-1022)) === Interval(-0x1.921FB54442D19p-1, +0x1.921FB54442D19p-1)

    @test atan(interval(-0x1p-1022, -0x1p-1022), interval(-0x1p-1022, 0x1p-1022)) === Interval(-0x1.2D97C7F3321D3p1, -0x1.921FB54442D18p-1)

    @test atan(interval(0x1p-1022, 0x1p-1022), interval(-0x1p-1022, 0x1p-1022)) === Interval(0x1.921FB54442D18p-1, 0x1.2D97C7F3321D3p1)

    @test atan(interval(-2.0, 2.0), interval(-3.0, -1.0)) === Interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)

    @test atan(interval(0.0, 2.0), interval(-3.0, -1.0)) === Interval(0x1.0468A8ACE4DF6p1, 0x1.921FB54442D19p1)

    @test atan(interval(1.0, 3.0), interval(-3.0, -1.0)) === Interval(0x1.E47DF3D0DD4Dp0, 0x1.68F095FDF593Dp1)

    @test atan(interval(1.0, 3.0), interval(-2.0, 0.0)) === Interval(0x1.921FB54442D18p0, 0x1.56C6E7397F5AFp1)

    @test atan(interval(1.0, 3.0), interval(-2.0, 2.0)) === Interval(0x1.DAC670561BB4Fp-2, 0x1.56C6E7397F5AFp1)

    @test atan(interval(1.0, 3.0), interval(0.0, 2.0)) === Interval(0x1.DAC670561BB4Fp-2, 0x1.921FB54442D19p0)

    @test atan(interval(1.0, 3.0), interval(1.0, 3.0)) === Interval(0x1.4978FA3269EE1p-2, 0x1.3FC176B7A856p0)

    @test atan(interval(0.0, 2.0), interval(1.0, 3.0)) === Interval(0x0p0, 0x1.1B6E192EBBE45p0)

    @test atan(interval(-2.0, 2.0), interval(1.0, 3.0)) === Interval(-0x1.1B6E192EBBE45p0, +0x1.1B6E192EBBE45p0)

    @test atan(interval(-2.0, 0.0), interval(1.0, 3.0)) === Interval(-0x1.1B6E192EBBE45p0, 0x0p0)

    @test atan(interval(-3.0, -1.0), interval(1.0, 3.0)) === Interval(-0x1.3FC176B7A856p0, -0x1.4978FA3269EE1p-2)

    @test atan(interval(-3.0, -1.0), interval(0.0, 2.0)) === Interval(-0x1.921FB54442D19p0, -0x1.DAC670561BB4Fp-2)

    @test atan(interval(-3.0, -1.0), interval(-2.0, 2.0)) === Interval(-0x1.56C6E7397F5AFp1, -0x1.DAC670561BB4Fp-2)

    @test atan(interval(-3.0, -1.0), interval(-2.0, 0.0)) === Interval(-0x1.56C6E7397F5AFp1, -0x1.921FB54442D18p0)

    @test atan(interval(-3.0, -1.0), interval(-3.0, -1.0)) === Interval(-0x1.68F095FDF593Dp1, -0x1.E47DF3D0DD4Dp0)

    @test atan(interval(-2.0, 0.0), interval(-3.0, -1.0)) === Interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)

    @test atan(interval(-5.0, 0.0), interval(-5.0, 0.0)) === Interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)

    @test atan(interval(0.0, 5.0), interval(-5.0, 0.0)) === Interval(0x1.921FB54442D18p0, 0x1.921FB54442D19p1)

    @test atan(interval(0.0, 5.0), interval(0.0, 5.0)) === Interval(0x0p0, 0x1.921FB54442D19p0)

    @test atan(interval(-5.0, 0.0), interval(0.0, 5.0)) === Interval(-0x1.921FB54442D19p0, 0x0p0)

end


text:
# check speed of new and old sparse filter

using Revise
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pandas
using Seaborn
using Dates
using Logging

using Profile
using ProfileView

b3, something = LarSurf.get_boundary3([30,30,30])
LarSurf.sparse_filter!(b3, 1, 1, 0)
LarSurf.sparse_filter_old!(b3, 1, 1, 0)
println("start")
@time LarSurf.sparse_filter!(b3, 1, 1, 0)
@time LarSurf.sparse_filter!(b3, 1, 1, 0)
@time LarSurf.sparse_filter!(b3, 1, 1, 0)
@time LarSurf.sparse_filter!(b3, 1, 1, 0)
@time LarSurf.sparse_filter!(b3, 1, 1, 0)

println("asdf")
# b3 = LarSurf.get_boundary3([30,30,30])
@time LarSurf.sparse_filter_old!(b3, 1, 1, 0)
# b3 = LarSurf.get_boundary3([30,30,30])
@time LarSurf.sparse_filter_old!(b3, 1, 1, 0)
@time LarSurf.sparse_filter_old!(b3, 1, 1, 0)
@time LarSurf.sparse_filter_old!(b3, 1, 1, 0)
@time LarSurf.sparse_filter_old!(b3, 1, 1, 0)


text:
function derivative(f, ftype::Symbol, dtype::Symbol)
  if ftype == :scalar
    return x::Number -> finite_difference(f, float(x), dtype)
  elseif ftype == :vector
    return x::Vector -> finite_difference(f, float(x), dtype)
  else
    error("ftype must :scalar or :vector")
  end
end
derivative(f, x::Union{T, Vector{T}}, dtype::Symbol = :central) where {T <: Number} = finite_difference(f, float(x), dtype)
derivative(f, dtype::Symbol = :central) = derivative(f, :scalar, dtype)

gradient(f, x::Union{T, Vector{T}}, dtype::Symbol = :central) where {T <: Number} = finite_difference(f, float(x), dtype)
gradient(f, dtype::Symbol = :central) = derivative(f, :vector, dtype)

@static if isdefined(Compat.LinearAlgebra, :gradient)
function Compat.LinearAlgebra.gradient(f, x::Union{T, Vector{T}}, dtype::Symbol = :central) where T <: Number
    Base.depwarn("The finite difference methods from Calculus.jl no longer extend " *
                 "Base.gradient and should be called as Calculus.gradient instead. " *
                 "This usage is deprecated.", :gradient)
    Calculus.gradient(f,x,dtype)
end

function Compat.LinearAlgebra.gradient(f, dtype::Symbol = :central)
    Base.depwarn("The finite difference methods from Calculus.jl no longer extend " *
                 "Base.gradient and should be called as Calculus.gradient instead. " *
                 "This usage is deprecated.", :gradient)
    Calculus.gradient(f,dtype)
end
end

function jacobian(f, x::Vector{T}, dtype::Symbol) where T <: Number
    finite_difference_jacobian(f, x, dtype)
end
function jacobian(f, dtype::Symbol)
    g(x::Vector) = finite_difference_jacobian(f, x, dtype)
    return g
end
jacobian(f) = jacobian(f, :central)

function second_derivative(f, g, ftype::Symbol, dtype::Symbol)
  if ftype == :scalar
    return x::Number -> finite_difference_hessian(f, g, x, dtype)
  elseif ftype == :vector
    return x::Vector -> finite_difference_hessian(f, g, x, dtype)
  else
    error("ftype must :scalar or :vector")
  end
end
function second_derivative(f, g, x::Union{T, Vector{T}}, dtype::Symbol) where T <: Number
  finite_difference_hessian(f, g, x, dtype)
end
function hessian(f, g, x::Union{T, Vector{T}}, dtype::Symbol) where T <: Number
  finite_difference_hessian(f, g, x, dtype)
end
function second_derivative(f, g, x::Union{T, Vector{T}}) where T <: Number
  finite_difference_hessian(f, g, x, :central)
end
function hessian(f, g, x::Union{T, Vector{T}}) where T <: Number
  finite_difference_hessian(f, g, x, :central)
end
function second_derivative(f, x::Number, dtype::Symbol)
  finite_difference_hessian(f, derivative(f), x, dtype)
end
function hessian(f, x::Number, dtype::Symbol)
  finite_difference_hessian(f, derivative(f), x, dtype)
end
function second_derivative(f, x::Vector{T}, dtype::Symbol) where T <: Number
  finite_difference_hessian(f, gradient(f), x, dtype)
end
function hessian(f, x::Vector{T}, dtype::Symbol) where T <: Number
  finite_difference_hessian(f, gradient(f), x, dtype)
end
function second_derivative(f, x::Number)
  finite_difference_hessian(f, derivative(f), x, :central)
end
function hessian(f, x::Number)
  finite_difference_hessian(f, derivative(f), x, :central)
end
function second_derivative(f, x::Vector{T}) where T <: Number
  finite_difference_hessian(f, gradient(f), x, :central)
end
function hessian(f, x::Vector{T}) where T <: Number
  finite_difference_hessian(f, gradient(f), x, :central)
end
second_derivative(f, g, dtype::Symbol) = second_derivative(f, g, :scalar, dtype)
second_derivative(f, g) = second_derivative(f, g, :scalar, :central)
second_derivative(f) = second_derivative(f, derivative(f), :scalar, :central)
hessian(f, g, dtype::Symbol) = second_derivative(f, g, :vector, dtype)
hessian(f, g) = second_derivative(f, g, :vector, :central)
hessian(f) = second_derivative(f, gradient(f), :vector, :central)


text:
import WAV

export Signal
export getdata
export getframerate
export analyse
export synthesize
export performdft
export performfft
export countseconds
export countchunks
export getchunk
export getfftchunk

"""
    Represents an audio signal.

fields: data and framerate.
"""
struct Signal
    data::Array{Complex,1}
    framerate::Float64
end

"""
    Signal(filepath [, channel])

Read a wave file and instantiate a Signal with the given `filepath` and
`channel`. Stereo audio files should contain 2 channels.
"""
Signal(filepath::String; channel::Int = 1) = begin
    data, framerate = WAV.wavread(filepath)
    Signal(Complex.(data[:, channel]), framerate)
end

getdata(signal::Signal) = signal.data
getframerate(signal::Signal) = signal.framerate

"""
    countseconds(signal)

Return the period of an audio signal in seconds with the given `signal`.
"""
function countseconds(signal::Signal)
    Int(length(getdata(signal)) ÷ getframerate(signal))
end

"""
    countchunks(signal, chunkspersecond)

Return the number of data chunks that an audio signal contains with the given
`signal` and chunks per second `chunkspersecond`.
"""
function countchunks(signal::Signal, chunkspersecond::Int)
    countseconds(signal) * chunkspersecond
end

"""
    analyse(signal)

Convert an audio signal from time domain to frequency domain with the given
`signal`.
"""
function analyse(signal::Signal)
    data = getdata(signal)
    framerate = getframerate(signal)
    N = length(data)
    f(ω) = begin
        scale = 1 / sqrt(N)
        estimate = 0
        for t = 0:N-1
            estimate += data[t+1] * exp(-im * 2π * ω / N * t)
        end
        scale * estimate
    end
    Ω = Array{Complex,1}(undef, N)
    for ω = 0:N-1
        Ω[ω+1] = f(ω)
    end
    Signal(Ω, framerate)
end

"""
    synthesize(signal)

Convert an audio signal from frequency domain to time domain with the given
`signal`.
"""
function synthesize(signal::Signal)
    data = getdata(signal)
    framerate = getframerate(signal)
    N = length(data)
    f(t) = begin
        scale = 1 / sqrt(N)
        estimate = 0
        for ω = 0:N-1
            estimate += data[ω+1] * exp(im * 2π * ω / N * t)
        end
        scale * estimate
    end
    T = Array{Complex,1}(undef, N)
    for t = 0:N-1
        T[t+1] = f(t)
    end
    Signal(T, framerate)
end

"""
    performdft(signal)

Perform a Discrete Fourier Transform (DFT) on an audio signal with the given
`signal`.
"""
function performdft(signal::Signal)
    analyse(signal)
end

"""
    performfft(signal)

Perform Fast Fourier Transform (FFT) on an audio signal with the given `signal`.
"""
function performfft(signal::Signal)
    data = getdata(signal)
    framerate = getframerate(signal)
    N = length(data)
    if N % 2 > 0
        println("The length of a data chunk must be a power of 2.")
        return analyse(signal)
    elseif N ≤ 2
        return analyse(signal)
    else
        e = getdata(performfft(Signal(view(data, 2:2:N), framerate)))
        o = getdata(performfft(Signal(view(data, 1:2:N-1), framerate)))
        r = convert(Array{Int64}, floor.(range(0, stop=N-1, length=N)))
        v = exp.(-im .* 2pi .* r ./ N)
        i = Integer(N ÷ 2)
        return Signal([e .+ v[1:i] .* o; e .+ v[i+1:end] .* o], framerate)
    end
end

"""
    getchunk(signal, i, [, chunkspersecond [, offset]])

Return a chunk of an audio signal with the given `signal`, the chunk ordinal
number `i` and the number of chunks per second `chunkspersecond`. The optional
arguments `offset` shifts the sampling window in the time dimension. Supply
`offset` for collecting multiple consecutive chunks and then taking the average.
"""
function getchunk(signal::Signal, i::Int; chunkspersecond::Int = 30,
                  offset::Int = 0)
    data = getdata(signal)
    framerate = getframerate(signal)
    samplesperchunk = Int(framerate ÷ chunkspersecond)
    window = 2^(Int(ceil(log2(samplesperchunk))))
    origin = (i - 1) * samplesperchunk + 1
    start = origin + offset
    finish = origin + window + offset
    totalsamples = size(data, 1)
    if finish ≤ totalsamples
        return Signal(data[start:finish-1], framerate)
    else
        return Signal(data[total_samples-window+1:end], framerate)
    end
end

"""
    getfftchunk(signal, i, chunkspersecond , samples)

Get a chunck of an audio signal under FFT by averaging over a number of
consecutive samples such that two consecutive samples differ only in a single
frame in the time dimension, with the given `signal`, the chunk cardinal number
`i`, the number of chunks per second `chunkspersecond` and the number of
consecutive `samples` to use before taking the average.
"""
function getfftchunk(signal::Signal, i::Int, chunkspersecond::Int, samples::Int)
    framerate = getframerate(signal)
    array = []
    for j in 1:samples
        sₜ = getdata(getchunk(signal, i, chunkspersecond = chunkspersecond,
                              offset = j))
        push!(array, getdata(performfft(Signal(sₜ, framerate))))
    end
    sum = array[1]
    for j in 2:samples
        sum = sum .+ array[j]
    end
    average = sum ./ samples
    Signal(average, framerate)
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
@testset "Mismatches" begin
    function test_mismatches(a, b)
        count = 0
        for (x, y) in zip(a, b)
            count += x != y
        end
        @test mismatches(a, b) === mismatches(b, a) === count
    end

    for len in 1:32, _ in 1:10
        a = random_dna_kmer(len)
        b = random_dna_kmer(len)
        test_mismatches(DNAKmer(a), DNAKmer(b))

        a = random_rna_kmer(len)
        b = random_rna_kmer(len)
        test_mismatches(RNAKmer(a), RNAKmer(b))
    end
end


text:
"""
    space_dimension(S::AbstractSet)

Return the space dimension of the given set.

### Input

- `S` -- set

### Output

An integer.

### Notes

The `space_dimension` refers to the ambient dimension of the set. For example,
a line in an ``n``-dimensional space has `space_dimension` equal to ``n``, even
if the dimension of its affine hull is ``1``.

See also [`dimension(::AbstractSet)`](@ref) for the dimension of the set.
"""
function space_dimension(::AbstractSet) end

"""
    dimension(S::AbstractSet)

Return the dimension of the given set.

### Input

- `S` -- set

### Output

An integer.

### Notes

The actual notion of `dimension` adopted is left to the package implementing
this function.

For instance, in many applications `dimension` refers to the dimension of the
affine hull of the set. For example, a line in an ``n``-dimensional space has
`dimension` equal to ``1``, even if the dimension of the space is ``n``.
However, in other domains, e.g. semialgebraic sets, `dimension` is not the
dimension of its affine hull.

See also [`space_dimension(::AbstractSet)`](@ref) for the ambient dimension of
the set.
"""
function dimension(::AbstractSet) end


text:
push!(LOAD_PATH, "src")

using SCIPCInterface

println("Loaded SCIP version: $(SCIPversion())")
scip = SCIPcreate()
SCIPincludeDefaultPlugins(scip)
SCIPcreateProbBasic(scip, "Test problem")
x1 = SCIPcreateVarBasic(scip, "x1", 0.0, 1.0, 4.0, SCIPCInterface.SCIP_VARTYPE_BINARY)
SCIPaddVar(scip, x1)
x2 = SCIPcreateVarBasic(scip, "x1", 0.0, 1.0, 3.0, SCIPCInterface.SCIP_VARTYPE_BINARY)
SCIPaddVar(scip, x2)

cons1 = SCIPcreateConsBasicLinear(scip, "cons1", [x1, x2], [1.0,2.0], 0.0, 1.0)
@assert SCIPgetRhsLinear(scip, cons1) == 1.0

println("done")


text:
struct SymbolicIntervalHeur{F} <: AbstractSymbolicIntervalBounds{F}
    sym::SymbolicInterval{F}
    lbs::Vector{Vector{Float64}}
    ubs::Vector{Vector{Float64}}
    var_los::Matrix{Float64} # symbolic lower bounds of fresh vars
    var_his::Matrix{Float64} # symbolic upper bounds of fresh vars
    max_vars::Int64
    importance::Vector{Float64}
end

function init_symbolic_interval_heur(net::NetworkNegPosIdx, input_set; max_vars=10)
    sym = init_deep_poly_symbolic_interval(input_set)
    lbs = [fill(-Inf, n_nodes(l)) for l in net.layers]
    ubs = [fill( Inf, n_nodes(l)) for l in net.layers]

    # symbolic bounds on variables have number of inputs plus constant
    var_his = zeros(max_vars, dim(input_set) + 1)
    var_los = zeros(max_vars, dim(input_set) + 1)

    # importance holds scores for each of the input variables
    importance = zeros(dim(input_set))
    return SymbolicIntervalHeur(sym, lbs, ubs, var_los, var_his, max_vars, importance)
end

function init_symbolic_interval_heur(s::SymbolicIntervalHeur, input_set::AbstractHyperrectangle; max_vars=10)
    sym = init_deep_poly_symbolic_interval(input_set)
    lbs = [copy(lb) for lb in s.lbs]
    ubs = [copy(ub) for ub in s.ubs]
    var_his = zeros(max_vars, dim(input_set) + 1)
    var_los = zeros(max_vars, dim(input_set) + 1)
    importance = zeros(dim(input_set))
    return SymbolicIntervalHeur(sym, lbs, ubs, var_los, var_his, max_vars, importance)
end

function convert2symbolic_interval_fv(s::SymbolicIntervalHeur)
    return SymbolicIntervalFV(s.sym, s.lbs, s.ubs, s.var_los, s.var_his, s.max_vars)
end

# TODO: only calculate upper bound, don't use bounds which calculates lo and up
# TODO: merge ρ for SymbolicIntervalFV and SymbolicIntervalHeur
function LazySets.ρ(d::AbstractArray{T,1} where T, sym::A where A<:SymbolicIntervalHeur)
    n_sym = size(sym.sym.Low, 2) - 1
    n_in = dim(domain(sym))
    current_n_vars = n_sym - n_in

    # substitute before interval_map as symbolic intervals are exact for affine
    # transformations and after substitution, dimension of symlo, symhi is smaller
    subs_sym_lo, subs_sym_hi = substitute_variables(sym.sym.Low, sym.sym.Up, sym.var_los, sym.var_his, n_in, current_n_vars)
    tmp_sym = SymbolicInterval(subs_sym_lo, subs_sym_hi, domain(sym))
    d_Low, d_Up = interval_map(reshape(d, 1, length(d)), tmp_sym.Low, tmp_sym.Up)
    sym_prime = SymbolicInterval(d_Low, d_Up, domain(sym))
    lo, up = bounds(sym_prime)
    # the result should only be one-dimensional
    return up[1]
end

# TODO: factor out into some utils file, don't want to have SymbolicIntervalFV function here
function maximizer(s::Union{SymbolicIntervalHeur, SymbolicIntervalFV})
    n_sym = size(s.sym.Low, 2) - 1
    n_in = dim(domain(s))
    current_n_vars = n_sym - n_in

    subs_sym_lo, subs_sym_hi = substitute_variables(s.sym.Low, s.sym.Up, s.var_los, s.var_his, n_in, current_n_vars)

    return maximizer(subs_sym_hi, low(domain(s)), high(domain(s)))
end

######### Splitting

function split_symbolic_interval_bounds(s::SymbolicIntervalHeur{<:Hyperrectangle}, index::Int64)
    lbs, ubs = low(domain(s)), high(domain(s))
    split_val = 0.5 * (lbs[index] + ubs[index])

    high1 = copy(ubs)
    high1[index] = split_val
    low2 = copy(lbs)
    low2[index] = split_val

    domain1 = Hyperrectangle(low=lbs, high=high1)
    domain2 = Hyperrectangle(low=low2, high=ubs)

    s1 = init_symbolic_interval_heur(s, domain1, max_vars=s.max_vars)
    s2 = init_symbolic_interval_heur(s, domain2, max_vars=s.max_vars)

    return [s1, s2]
end

"""
Splitting based on intermediate importance scores (i.e. coefficients in all crossing ReLUs)
"""
function split_important_interval(s::SymbolicIntervalHeur)
    radius = high(domain(s)) - low(domain(s))
    # if there are no more crossing ReLUs, importance will be zero vector -> if we
    # multiply it with radius, all inputs are equally important, and argmax always
    # returns the first index -> infinitely loop splitting
    most_important_dim = sum(s.importance) == 0. ? argmax(radius) : argmax(s.importance .* radius)
    return split_symbolic_interval_bounds(s, most_important_dim)
end


text:
precompile(read_abinit_density, (AbstractString,))
precompile(readWAVECAR, ())
precompile(readDOSCAR, ())

max_stars_repo_name: brainandforce/Xtal

max_issues_repo_name: brainandforce/Xtal

max_forks_repo_name: brainandforce/Xtal


text:
"""
    format_mangal_coordinates(d::Dict{T,Any}) where {T <: AbstractString}

Returns a set of coordinates in a `GeoInterface` object, which can be a `Point`
or a `Polygon`.
"""
function format_mangal_coordinates(d::Dict{T,Any}) where {T <: AbstractString}
    point_type = d["geom"]["type"] == "Point" ? Point : Polygon
    if point_type == Polygon
        coords = [float.(x) for x in first(d["geom"]["coordinates"])]
        point_coordinates = point_type(coords)
    else
        coords = float.(d["geom"]["coordinates"])
        point_coordinates = point_type(coords...)
    end
    return point_coordinates
end

function format_mangal_response(::Type{MangalDataset}, d::Dict{T,Any}) where {T <: AbstractString}
    obj_id = d["id"]
    obj_public = d["public"]
    obj_name = d["name"]
    obj_created = DateTime(d["created_at"][1:19])
    obj_updated = DateTime(d["updated_at"][1:19])
    obj_reference = isnothing(d["ref_id"]) ? missing : reference(d["ref_id"])
    obj_user = d["user_id"]
    obj_description = d["description"]

    return MangalDataset(obj_id, obj_public, obj_name, obj_created,
        obj_updated, obj_reference, obj_user, obj_description
        )
end

function format_mangal_response(::Type{MangalNetwork}, d::Dict{T,Any}) where {T <: AbstractString}
    obj_id = d["id"]
    obj_public = d["public"]
    obj_name = d["name"]
    obj_date = isnothing(d["date"]) ? missing : DateTime(d["date"][1:19])
    obj_position = isnothing(d["geom"]) ? missing : format_mangal_coordinates(d)
    obj_created = DateTime(d["created_at"][1:19])
    obj_updated = DateTime(d["updated_at"][1:19])
    obj_user = isnothing(d["user_id"]) ? missing : d["user_id"]
    obj_description = d["description"]
    obj_complete = d["all_interactions"]
    obj_dataset = isnothing(d["dataset_id"]) ? missing : dataset(d["dataset_id"])

    return MangalNetwork(obj_id, obj_public, obj_name, obj_date, obj_position,
        obj_created, obj_updated, obj_user, obj_description,
        obj_complete, obj_dataset)

end

function format_mangal_response(::Type{MangalNode}, d::Dict{T,Any}) where {T <: AbstractString}

    obj_id = d["id"]
    obj_name = d["original_name"]
    obj_created = DateTime(d["created_at"][1:19])
    obj_updated = DateTime(d["updated_at"][1:19])
    obj_taxon = isnothing(d["taxonomy"]) ? missing : Mangal.format_mangal_response(MangalReferenceTaxon, d["taxonomy"])

    return MangalNode(obj_id, obj_name, obj_created, obj_updated, obj_taxon)

end

function format_mangal_response(::Type{MangalReferenceTaxon}, d::Dict{T,Any}) where {T <: AbstractString}

    obj_id = d["id"]
    obj_name = d["name"]
    obj_bold = isnothing(d["bold"]) ? missing : d["bold"]
    obj_tsn = isnothing(d["tsn"]) ? missing : d["tsn"]
    obj_ncbi = isnothing(d["ncbi"]) ? missing : d["ncbi"]
    obj_eol = isnothing(d["eol"]) ? missing : d["eol"]
    obj_gbif = isnothing(d["gbif"]) ? missing : d["gbif"]
    obj_created = DateTime(d["created_at"][1:19])
    obj_updated = DateTime(d["updated_at"][1:19])

    return MangalReferenceTaxon(obj_id, obj_name, obj_bold, obj_tsn,
        obj_ncbi, obj_eol, obj_gbif, obj_created, obj_updated)

end

function format_mangal_response(::Type{MangalInteraction}, d::Dict{T,Any}) where {T <: AbstractString}
    obj_id = d["id"]
    obj_network = network(d["network_id"])
    obj_from = node(d["node_from"])
    obj_to = node(d["node_to"])
    obj_date = isnothing(d["date"]) ? missing : DateTime(d["date"][1:19])
    obj_position = isnothing(d["geom"]) ? missing : format_mangal_coordinates(d)
    obj_directed = d["direction"] == "directed"
    obj_interaction = Symbol(d["type"])
    obj_method = isnothing(d["method"]) ? missing : d["method"]
    obj_strength = isnothing(d["value"]) ? missing : d["value"]
    obj_created = DateTime(d["created_at"][1:19])
    obj_updated = DateTime(d["updated_at"][1:19])
    obj_attribute = isnothing(d["attr_id"]) ? missing : format_mangal_response(MangalAttribute, d["attribute"])

    return MangalInteraction(obj_id, obj_network, obj_from, obj_to, obj_date, obj_position,
        obj_directed, obj_interaction, obj_method, obj_strength, obj_created, obj_updated,
        obj_attribute)

end

function format_mangal_response(::Type{MangalReference}, d::Dict{T,Any}) where {T <: AbstractString}
    obj_id = d["id"]
    obj_year = d["year"] == "NA" ? missing : parse(Int64, d["year"][1:4])
    obj_doi = isnothing(d["doi"]) ? missing : d["doi"]
    obj_jstor = isnothing(d["jstor"]) ? missing : d["jstor"]
    obj_pmid = isnothing(d["pmid"]) ? missing : d["pmid"]
    obj_bibtex = isnothing(d["bibtex"]) ? missing : d["bibtex"]
    obj_paper = isnothing(d["paper_url"]) ? missing : d["paper_url"]
    obj_data = isnothing(d["data_url"]) ? missing : d["data_url"]

    return MangalReference(obj_id, obj_year, obj_doi, obj_jstor, obj_pmid,
        obj_bibtex, obj_paper, obj_data)
end

function format_mangal_response(::Type{MangalAttribute}, d::Dict{T,Any}) where {T <: AbstractString}
    obj_id = d["id"]
    obj_name = d["name"]
    obj_description = d["description"]
    obj_unit = isnothing(d["unit"]) ? missing : d["unit"]

    return MangalAttribute(obj_id, obj_name, obj_description, obj_unit)
end


text:
function create_system_data(; with_forecasts=false, time_series_in_memory=false)
    data = IS.SystemData(; time_series_in_memory=time_series_in_memory)

    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(data, component)

    if with_forecasts
        file = joinpath(FORECASTS_DIR, "ComponentsAsColumnsNoTime.json")
        IS.add_forecasts!(IS.InfrastructureSystemsType, data, file)

        forecasts = get_all_forecasts(data)
        @assert length(forecasts) > 0
    end

    return data
end

function get_all_forecasts(data)
    return collect(IS.iterate_forecasts(data))
end

function create_time_series_data()
    dates = collect(Dates.DateTime("1/1/2020 00:00:00", "d/m/y H:M:S") : Dates.Hour(1):Dates.DateTime("1/1/2020 23:00:00", "d/m/y H:M:S"))
    data = collect(1:24)
    component_name = "gen"
    ta = TimeSeries.TimeArray(dates, data, [component_name])
    return IS.TimeSeriesData(ta)
end


text:
using CSV, DataFrames, Statistics

results_dir = ARGS[1]
results = CSV.read("$results_dir/results.csv")

by(results, [:sampler, :Nd], :tau_ess=>mean, :A_ess=>mean, :k_ess=>mean, Symbol("v[1]_ess")=>mean, Symbol("v[2]_ess")=>mean, Symbol("v[3]_ess")=>mean, :tree_depth=>mean, :epsilon=>mean) |> println


text:
import DocSeeker: baseURL, finddocsURL, readmepath

@testset "finddocs" begin
    @testset "finddocsURL" begin
        @test baseURL(finddocsURL("base")) == "https://docs.julialang.org"
    end

    @testset "readmepath" begin
        @test readmepath("DocSeeker") == abspath(joinpath(@__DIR__, "..", "README.md"))
    end
end


text:
# Abstract type for environments
abstract type Fym end


text:
using MultimodeCavity
const mm = MultimodeCavity

Nparticles = 2
Nlevels = 5
E0 = 0.1
P = mm.potentials.BoxPotential

index = 3
Nmodes = 10
delta = -1
eta = 1
kappa = 2
U0 = 1

s = 0.4

particles_bosons = mm.quantum.Bosons(Nparticles, Nlevels, E0, P)
particles_fermions = mm.quantum.Fermions(Nparticles, Nlevels, E0, P)
mode1 = mm.quantum.CavityMode(index, Nmodes, delta, eta, kappa, U0)
mode2 = mm.quantum.CavityMode(index, Nmodes, delta, eta, kappa, U0)
system = mm.quantum.MultimodeSystem(s, particles_bosons, [mode1, mode2])
system = mm.quantum.MultimodeSystem(s, particles_fermions, [mode1, mode2])

b = mm.quantum.basis(particles_bosons)
b = mm.quantum.basis(particles_fermions)
b = mm.quantum.basis(mode1)
b = mm.quantum.basis(system)

An = mm.quantum.An(1, system)
An = mm.quantum.An(2, system)

Bn = mm.quantum.Bn(1, system)
Bn = mm.quantum.Bn(2, system)

H = mm.quantum.Hamiltonian(particles_fermions)
H = mm.quantum.Hamiltonian(particles_bosons)
H = mm.quantum.Hamiltonian(mode1)
H = mm.quantum.Hamiltonian(system)

J = mm.quantum.Jump_operators(system)