text:
#####
##### PushVector
#####

mutable struct PushVector{T} <: AbstractVector{T}
    data::Vector{T}
    len::Int
end

PushVector{T}() where {T} = PushVector{T}(T[], 0)
Base.size(v::PushVector) = (v.len,)
maxlength(v::PushVector) = length(v.data)

Base.IndexStyle(::PushVector) = Base.IndexLinear()
Base.@propagate_inbounds function Base.getindex(A::PushVector, i::Int)
    @boundscheck checkbounds(A, i)
    return @inbounds(A.data[i])
end

Base.@propagate_inbounds function Base.setindex!(A::PushVector, v, i::Int)
    @boundscheck checkbounds(A, i)
    return @inbounds(A.data[i] = v)
end

function Base.push!(A::PushVector, v)
    len = length(A)
    if len < maxlength(A)
        @inbounds(A.data[len + 1] = v)
    else
        push!(A.data, v)
    end
    A.len += 1
    return A
end

Base.empty!(A::PushVector) = (A.len = 0)

#####
##### Node
#####

@enum Direction::UInt8 Left Right
Base.reverse(x::Direction) = Base.bitcast(Direction, xor(UInt8(x), one(UInt8)))

const SEED = UInt(0x167ef00a9)
mutable struct Node{T}
    key::T
    priority::UInt
    nchildren::Int
    left::Union{Nothing,Node{T}}
    right::Union{Nothing,Node{T}}

    function Node{T}(key::T) where {T}
        priority = hash(key, SEED)
        return new{T}(key, priority, 1, nothing, nothing)
    end
end

Node(key::T) where {T} = Node{T}(key)

function Base.show(io::IO, node::Node{T}) where {T}
    hasleft = lchild(node) !== nothing
    hasright = rchild(node) !== nothing
    print(
        io,
        "Node{$T}($(node.key), $(node.priority), $(node.nchildren), $hasleft, $hasright)",
    )
    return nothing
end

mutable struct Treap{T}
    root::Union{Nothing,Node{T}}
    path::PushVector{Tuple{Node{T},Direction}}
    nodebuffer::Vector{Node{T}}
end
Treap{T}() where {T} = Treap{T}(nothing, PushVector{Tuple{Node{T},Direction}}(), Node{T}[])

lchild(n::Node) = n.left
rchild(n::Node) = n.right
hasleft(n::Node) = lchild(n) !== nothing
hasright(n::Node) = rchild(n) !== nothing

@inline setleft!(n::Node, v::Node) = (n.left = v)
@inline setleft!(n::Node) = (n.left = nothing)
@inline setright!(n::Node, v::Node) = (n.right = v)
@inline setright!(n::Node) = (n.right = nothing)

function reassign!(node::Node{T}, key::T) where {T}
    node.key = key
    node.priority = hash(key, SEED)
    node.nchildren = 1
    node.left = nothing
    node.right = nothing
    return nothing
end

# Compute number of children.
function update!(node::Node)
    subtree = 1
    left = lchild(node)
    left === nothing || (subtree += left.nchildren)
    right = rchild(node)
    right === nothing || (subtree += right.nchildren)
    node.nchildren = subtree
    return subtree
end

function rotate_right!(root::Node)
    newroot = lchild(root)
    newroot === nothing && return root
    right = rchild(newroot)
    right === nothing ? setleft!(root) : setleft!(root, right)
    setright!(newroot, root)
    update!(root)
    update!(newroot)
    return newroot
end

function rotate_left!(root::Node)
    newroot = rchild(root)
    newroot === nothing && return root
    left = lchild(newroot)
    left === nothing ? setright!(root) : setright!(root, left)
    setleft!(newroot, root)
    update!(root)
    update!(newroot)
    return newroot
end

function search(root::Node, key)
    root.key == key && return root
    next = (key < root.key) ? lchild(root) : rchild(root)
    return next === nothing ? nothing : search(next, key)
end
Base.in(key, root::Node) = (search(root, key) !== nothing)
Base.haskey(root::Node, key) = (search(root, key) !== nothing)

function nchildren(root::Node, key)
    count = 0
    while true
        if root.key == key
            right = rchild(root)
            if right !== nothing
                count += right.nchildren
            end
            return count
        elseif key < root.key
            count += 1
            right = rchild(root)
            if right !== nothing
                count += right.nchildren
            end
            next = lchild(root)
        else
            next = rchild(root)
        end
        next === nothing && error()
        root = next
    end
end

#####
##### insertion
#####

# Node only method
function Base.insert!(root::Node, node)
    if node.key < root.key
        _left = lchild(root)
        left = setleft!(root, _left === nothing ? node : insert!(_left, node))
        if left.priority > root.priority
            root = rotate_right!(root)
        end
    else
        _right = rchild(root)
        right = setright!(root, _right === nothing ? node : insert!(_right, node))
        if right.priority > root.priority
            root = rotate_left!(root)
        end
    end
    update!(root)
    return root
end

# Using treap auxiliary data structures
function insertdown!(treap::Treap, node)
    _current = root(treap)
    _current === nothing && return nothing
    current = _current
    path = treap.path
    empty!(path)
    while true
        if node.key < current.key
            push!(path, (current, Left))
            next = lchild(current)
        else
            push!(path, (current, Right))
            next = rchild(current)
        end
        next === nothing && return nothing
        current = next
    end
end

function insertup!(treap::Treap, node::Node)
    path = treap.path
    i = length(path)
    @inbounds while i > 0
        root, direction = path[i]
        if direction == Left
            left = setleft!(root, node)
            left.priority > root.priority ? rotate_right!(root) : break
        else
            right = setright!(root, node)
            right.priority > root.priority ? rotate_left!(root) : break
        end
        update!(node)
        i -= 1
    end

    # If the path is empty, then we percolated updates all the way to the root and
    # need to update the "treap's" root.
    # Otherwise, we need to bookkeep the children numbers on the path we visited.
    if iszero(i)
        treap.root = node
    else
        @inbounds for j in Base.OneTo(i)
            root, _ = path[j]
            root.nchildren += 1
        end
    end
    return treap
end

#####
##### deletion
#####

macro _split(sym, f)
    return esc(quote
        if $sym !== nothing
            _tmp = delete!($sym, key)
            _tmp === nothing ? $f(root) : $f(root, _tmp)
        end
    end)
end

function Base.delete!(root::Node{T}, key::T) where {T}
    left = lchild(root)
    right = rchild(root)
    # Branch left if possible
    if (key < root.key)
        @_split left setleft!
        # Branch right if possible
    elseif (key > root.key)
        @_split right setright!
        # Key is at node, try removal of left or right.
    elseif left === nothing
        root = right
    elseif right === nothing
        root = left
        # Node has two children, rotate and move down.
    elseif left.priority < right.priority
        old = root
        root = rotate_left!(root)
        _tmp = delete!(old, key)
        _tmp === nothing ? setleft!(root) : setleft!(root, _tmp)
    else
        old = root
        root = rotate_right!(root)
        _tmp = delete!(old, key)
        _tmp === nothing ? setright!(root) : setright!(root, _tmp)
    end
    root === nothing || update!(root)
    return root
end

function deletefind!(treap::Treap, current::Node, key)
    path = treap.path
    empty!(path)
    while true
        current.key == key && return current
        if key < current.key
            push!(path, (current, Left))
            next = lchild(current)
        else
            push!(path, (current, Right))
            next = rchild(current)
        end
        next === nothing && return nothing
        current = next
    end
end

function handleroot!(treap::Treap, node)
    path = treap.path
    left = lchild(node)
    right = rchild(node)

    if left === nothing
        if right === nothing
            # Delete node - return true to indicate that upper level function
            # should simply return.
            treap.root = nothing
            push!(treap.nodebuffer, node)
            return true
        else
            newroot = rotate_left!(node)
            push!(path, (newroot, Left))
            treap.root = newroot
        end
    elseif right === nothing
        newroot = rotate_right!(node)
        push!(path, (newroot, Right))
        treap.root = newroot
    elseif left.priority < right.priority
        newroot = rotate_left!(node)
        push!(path, (right, Left))
        treap.root = newroot
    else
        newroot = rotate_right!(node)
        push!(path, (left, Right))
        treap.root = newroot
    end
    return false
end

function deletedown!(treap::Treap, root, node)
    updateroot = (node == root)
    if updateroot && handleroot!(treap, node)
        return nothing
    end

    # Shift the node down until it falls off the bottom of the tree.
    path = treap.path
    while true
        parent, direction = @inbounds(path[end])
        left = lchild(node)
        right = rchild(node)

        # Easy path, both children are nothing, we just need to
        if left === nothing
            if right === nothing
                direction === Left ? setleft!(parent) : setright!(parent)
                push!(treap.nodebuffer, node)
                break
            else
                rotate_left!(node)
                direction === Left ? setleft!(parent, right) : setright!(parent, right)
                push!(path, (right, Left))
            end
        elseif right === nothing
            rotate_right!(node)
            direction === Left ? setleft!(parent, left) : setright!(parent, left)
            push!(path, (left, Right))
        elseif left.priority < right.priority
            # Right will become the successor
            rotate_left!(node)
            direction === Left ? setleft!(parent, right) : setright!(parent, right)
            push!(path, (right, Left))
        else
            rotate_right!(node)
            direction === Left ? setleft!(parent, left) : setright!(parent, left)
            push!(path, (left, Right))
        end
    end

    for (n, _) in path
        n.nchildren -= 1
    end
    return nothing
end

#####
##### misc
#####

function traverse(f::F, root::Node) where {F}
    left = lchild(root)
    if left !== nothing
        traverse(f, left)
    end
    f(root)
    right = rchild(root)
    if right !== nothing
        traverse(f, right)
    end
    return nothing
end

function list(root::Node{T}) where {T}
    v = T[]
    traverse(x -> push!(v, x.key), root)
    return v
end

#####
##### Treap
#####

function Base.length(treap::Treap)
    base = root(treap)
    return base === nothing ? 0 : base.nchildren
end
root(treap::Treap) = treap.root

function nchildren(treap::Treap, key)
    base = root(treap)
    base === nothing && return -1
    return nchildren(base, key)
end

function Base.in(key0, treap::Treap{T}) where {T}
    node = root(treap)
    node === nothing && return false
    return in(convert(T, key0), node)
end

function unsafe_push!(treap::Treap{T}, key0) where {T}
    buffer = treap.nodebuffer
    key = convert(T, key0)
    if !isempty(buffer)
        node = pop!(buffer)
        reassign!(node, key)
    else
        node = Node(key)
    end

    base = root(treap)
    if base === nothing
        treap.root = node
        return treap
    end

    insertdown!(treap, node)
    insertup!(treap, node)
    return treap
end

function Base.push!(treap::Treap{T}, key0) where {T}
    key = convert(T, key0)
    in(key, treap) || unsafe_push!(treap, key)
    return treap
end

function Base.delete!(treap::Treap{T}, key0) where {T}
    key = convert(T, key0)
    # If the root doesn't exist, then obviously we can't delete the key.
    base = root(treap)
    base === nothing && return false

    # If `deletefind!` doesn't find the node - then we've done a bit of extra work by
    # keeping track of our downward path.
    #
    # But, if we DO find the node, then we can save redoing the work we just did, which
    # is potentially much more useful.
    node = deletefind!(treap, base, key)
    node === nothing && return false
    deletedown!(treap, base, node)
    return true
end

#####
##### Treap verification
#####

function istree(treap::Treap)
    base = root(treap)
    return base === nothing ? true : istree(base)
end
function istree(node::Node{T}) where {T}
    v = T[]
    traverse(node) do n
        push!(v, n.key)
    end
    return issorted(v)
end

function isheap(treap::Treap)
    base = root(treap)
    return base === nothing ? true : isheap(base)
end
function isheap(node::Node)
    passed = true
    left = lchild(node)
    if left !== nothing
        passed &= isheap(left)
        passed &= (node.priority > left.priority)
    end
    right = rchild(node)
    if right !== nothing
        passed &= isheap(right)
        passed &= (node.priority > right.priority)
    end

    return passed
end

#####
##### AbstractTrees
#####

function AbstractTrees.children(node::Node)
    if hasleft(node)
        if hasright(node)
            return (lchild(node), rchild(node))
        else
            return (lchild(node),)
        end
    else
        if hasright(node)
            return (rchild(node),)
        else
            return ()
        end
    end
end

AbstractTrees.printnode(io::IO, node::Node) = print(io, "($(node.key), $(node.priority))")


text:
"""
    Revise.track(Base)
    Revise.track(Core.Compiler)
    Revise.track(stdlib)

Track updates to the code in Julia's `base` directory, `base/compiler`, or one of its
standard libraries.
"""
function track(mod::Module; modified_files=revision_queue)
    id = PkgId(mod)
    modname = nameof(mod)
    return _track(id, modname; modified_files=modified_files)
end

function _track(id, modname; modified_files=revision_queue)
    inpath(path, dirs) = all(dir->occursin(dir, path), dirs)
    haskey(pkgdatas, id) && return nothing  # already tracked
    isbase = modname == :Base
    isstdlib = !isbase && modname ∈ stdlib_names
    if isbase || isstdlib
        # Test whether we know where to find the files
        if isbase
            srcdir = fixpath(joinpath(juliadir, "base"))
            dirs = ["base"]
        else
            stdlibv = joinpath("stdlib", "v$(VERSION.major).$(VERSION.minor)", String(modname))
            srcdir = fixpath(joinpath(juliadir, stdlibv))
            if !isdir(srcdir)
                srcdir = fixpath(joinpath(juliadir, "stdlib", String(modname)))
            end
            if !isdir(srcdir)
                # This can happen for Pkg, since it's developed out-of-tree
                srcdir = joinpath(juliadir, "usr", "share", "julia", stdlibv)
            end
            dirs = ["stdlib", String(modname)]
        end
        if !isdir(srcdir)
            @error "unable to find path containing source for $modname, tracking is not possible"
        end
        # Determine when the basesrccache was built
        mtcache = mtime(basesrccache)
        # Initialize expression-tracking for files, and
        # note any modified since Base was built
        if !haskey(pkgdatas, id)
            pkgdatas[id] = PkgData(id, srcdir)
        end
        pkgdata = pkgdatas[id]
        for (submod, filename) in Iterators.drop(Base._included_files, 1)  # stepping through sysimg.jl rebuilds Base, omit it
            ffilename = fixpath(filename)
            inpath(ffilename, dirs) || continue
            keypath = ffilename[1:last(findfirst(dirs[end], ffilename))]
            rpath = relpath(ffilename, keypath)
            fullpath = joinpath(basedir(pkgdata), rpath)
            if fullpath != filename
                cache_file_key[fullpath] = filename
                src_file_key[filename] = fullpath
            end
            push!(pkgdata, rpath=>FileInfo(submod, basesrccache))
            if mtime(ffilename) > mtcache
                with_logger(_debug_logger) do
                    @debug "Recipe for Base/StdLib" _group="Watching" filename=filename mtime=mtime(filename) mtimeref=mtcache
                end
                push!(modified_files, (pkgdata, rpath))
            end
        end
        # Add files to CodeTracking pkgfiles
        CodeTracking._pkgfiles[id] = pkgdata.info
        # Add the files to the watch list
        init_watching(pkgdata, srcfiles(pkgdata))
    elseif modname == :Compiler
        compilerdir = normpath(joinpath(juliadir, "base", "compiler"))
        if !haskey(pkgdatas, id)
            pkgdatas[id] = PkgData(id, compilerdir)
        end
        track_subdir_from_git(id, compilerdir; modified_files=modified_files)
    else
        error("no Revise.track recipe for module ", modname)
    end
    return nothing
end

# Fix paths to files that define Julia (base and stdlibs)
function fixpath(filename::AbstractString; badpath=basebuilddir, goodpath=juliadir)
    startswith(filename, badpath) || return normpath(filename)
    filec = filename
    relfilename = relpath(filename, badpath)
    relfilename0 = relfilename
    for strippath in (#joinpath("usr", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)"),
                      joinpath("usr", "share", "julia"),)
        if startswith(relfilename, strippath)
            relfilename = relpath(relfilename, strippath)
            if occursin("stdlib", relfilename0) && !occursin("stdlib", relfilename)
                relfilename = joinpath("stdlib", relfilename)
            end
        end
    end
    ffilename = normpath(joinpath(goodpath, relfilename))
    if (isfile(filename) & !isfile(ffilename))
        ffilename = normpath(filename)
    end
    return ffilename
end
_fixpath(lnn; kwargs...) = LineNumberNode(lnn.line, Symbol(fixpath(String(lnn.file); kwargs...)))
fixpath(lnn::LineNumberNode; kwargs...) = _fixpath(lnn; kwargs...)
fixpath(lnn::Core.LineInfoNode; kwargs...) = _fixpath(lnn; kwargs...)

# For tracking subdirectories of Julia itself (base/compiler, stdlibs)
function track_subdir_from_git(id::PkgId, subdir::AbstractString; commit=Base.GIT_VERSION_INFO.commit, modified_files=revision_queue)
    if !haskey(pkgdatas, id)
        pkgdatas[id] = PkgData(id)
    end
    pkgdata = pkgdatas[id]
    # diff against files at the same commit used to build Julia
    repo, repo_path = git_repo(subdir)
    if repo == nothing
        throw(GitRepoException(subdir))
    end
    prefix = relpath(subdir, repo_path)   # git-relative path of this subdir
    tree = git_tree(repo, commit)
    files = Iterators.filter(file->startswith(file, prefix) && endswith(file, ".jl"), keys(tree))
    ccall((:giterr_clear, :libgit2), Cvoid, ())  # necessary to avoid errors like "the global/xdg file 'attributes' doesn't exist: No such file or directory"
    for file in files
        fullpath = joinpath(repo_path, file)
        rpath = relpath(fullpath, pkgdata)  # this might undo the above, except for Core.Compiler
        local src
        try
            src = git_source(file, tree)
        catch err
            if err isa KeyError
                @warn "skipping $file, not found in repo"
                continue
            end
            rethrow(err)
        end
        fmod = get(juliaf2m, fullpath, Core.Compiler)  # Core.Compiler is not cached
        fmod === Core.Compiler && (endswith(fullpath, "compiler.jl") || endswith(fullpath, "tfuncs.jl")) && continue  # defines the module, skip
        if src != read(fullpath, String)
            push!(modified_files, (pkgdata, rpath))
        end
        fi = FileInfo(fmod)
        if parse_source!(fi.modexsigs, src, file, fmod) === nothing
            @warn "failed to parse Git source text for $file"
        else
            instantiate_sigs!(fi.modexsigs)
        end
        push!(pkgdata, rpath=>fi)
    end
    if !isempty(pkgdata.fileinfos)
        pkgdatas[id] = pkgdata
        CodeTracking._pkgfiles[id] = pkgdata.info
        init_watching(pkgdata, srcfiles(pkgdata))
    end
    return nothing
end

# For tracking Julia's own stdlibs
const stdlib_names = Set([
    :Base64, :CRC32c, :Dates, :DelimitedFiles, :Distributed,
    :FileWatching, :Future, :InteractiveUtils, :Libdl,
    :LibGit2, :LinearAlgebra, :Logging, :Markdown, :Mmap,
    :OldPkg, :Pkg, :Printf, :Profile, :Random, :REPL,
    :Serialization, :SHA, :SharedArrays, :Sockets, :SparseArrays,
    :Statistics, :SuiteSparse, :Test, :Unicode, :UUIDs])

# This replacement is needed because the path written during compilation differs from
# the git source path
const stdlib_rep = joinpath("usr", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)") => "stdlib"

const juliaf2m = Dict(normpath(replace(file, stdlib_rep))=>mod
    for (mod,file) in Base._included_files)

ext: jl


text:
function parabolicinterp(f::Function, xl::Number, xm::Number, xu::Number, ea::Number, maxIter::Int64)
  # 1D optimisation using parabolic interpolation

  xo = NaN; # Current estimate of the optimum
  fxm = f(xm);
  fxu = f(xu);
  fxl = f(xl);

  for i = 1:maxIter
    xoold = xo; # xoxo

    # Solution for xo using simultaneous equations
    fracNumerator = ((xm - xl)^2 * (fxm-fxu)) - ((xm - xu)^2 * (fxm - fxl));
    fracDenominator = ((xm - xl) * (fxm - fxu)) - ((xm - xu) * (fxm - fxl));
    xo = xm - 0.5 * (fracNumerator / fracDenominator);

    if xo < xm
      # The optimum is to the left of xm so make the upper estimate xm and make
      # the middle estimate xo.
      xu = xm;
      fxu = fxm;
      xm = xo;
      fxm = f(xm);
    elseif xo > xm
      # The optimum is to the right of xm so make the lower estimate xm and
      # the middle estimate xo.
      xl = xm;
      fxl = fxm;
      xm = xo;
      fxm = f(xm);
    else
      break;
    end

    # Relative error checks
    if abs((xo - xoold) / xo) * 100 < ea
      break;
    end
  end
  return xo;
end

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_name: alexjohnj/numerical-methods


text:
module Unums

using Compat

import Base: print, show, showcompact, convert, zero, one
import Base: +, -, *, /, ^, abs, ==, <, <=, sqrt, signbit, in, isfinite, isnan

export Unum, Unum22, Ubound22, Unum34, Ubound34, Bnum, Bbound, isexact, esmax, fsmax, emax, fmax, limits, print_bits

include("bnum.jl")
include("unum.jl")

end # module


text:
#
# Types supporting instantiated models and their components
#

# Supertype for variables and parameters in component instances
@class ComponentInstanceData{NT <: NamedTuple} <: MimiClass begin
    nt::NT
    comp_paths::Vector{ComponentPath}   # records the origin of each datum
end

nt(obj::AbstractComponentInstanceData) = getfield(obj, :nt)
types(obj::AbstractComponentInstanceData) = typeof(nt(obj)).parameters[2].parameters
Base.names(obj::AbstractComponentInstanceData)  = keys(nt(obj))
Base.values(obj::AbstractComponentInstanceData) = values(nt(obj))

# Centralizes the shared functionality from the two component data subtypes.
function _datum_instance(subtype::Type{<: AbstractComponentInstanceData},
                         names, types, values, paths)
    # @info "_datum_instance: names=$names, types=$types"
    NT = NamedTuple{Tuple(names), Tuple{types...}}
    return subtype(NT(values), Vector{ComponentPath}(paths))
end

@class ComponentInstanceParameters <: ComponentInstanceData begin
    function ComponentInstanceParameters(nt::NT, paths::Vector{ComponentPath}) where {NT <: NamedTuple}
        return new{NT}(nt, paths)
    end

    function ComponentInstanceParameters(names::Vector{Symbol},
                                         types::Vector{DataType},
                                         values::Vector{Any},
                                         paths)
        return _datum_instance(ComponentInstanceParameters, names, types, values, paths)
    end
end

@class ComponentInstanceVariables <: ComponentInstanceData begin
    function ComponentInstanceVariables(nt::NT, paths::Vector{ComponentPath}) where {NT <: NamedTuple}
        return new{NT}(nt, paths)
    end

    function ComponentInstanceVariables(names::Vector{Symbol},
                                        types::Vector{DataType},
                                        values::Vector{Any},
                                        paths)
        return _datum_instance(ComponentInstanceVariables, names, types, values, paths)
    end
end

# A container class that wraps the dimension dictionary when passed to run_timestep()
# and init(), so we can safely implement Base.getproperty(), allowing `d.regions` etc.
# All values in the dictionary are vectors of Ints, except the `:time` value, which is a
# vector of AbstractTimesteps, so that `d.time` returns values that can be used for indexing
# into timestep arrays.
struct DimValueDict <: MimiStruct
    dict::Dict{Symbol, Union{Vector{Int}, Vector{AbstractTimestep}}}

    function DimValueDict(dim_dict::AbstractDict, clock::Clock)
        d = Dict([name => (name == :time ? timesteps(clock) : collect(values(dim))) for (name, dim) in dim_dict])
        new(d)
    end
end

# Special case support for Dicts so we can use dot notation on dimension.
# The run_timestep() and init() funcs pass a DimValueDict of dimensions by name
# as the "d" parameter.
Base.getproperty(obj::DimValueDict, property::Symbol) = getfield(obj, :dict)[property]

# Superclass for both LeafComponentInstance and CompositeComponentInstance.
# This allows the former to be type-parameterized and the latter to not be.
@class mutable ComponentInstance <: MimiClass begin
    comp_name::Symbol
    comp_id::ComponentId
    comp_path::ComponentPath
    first::Union{Nothing, Int}
    last::Union{Nothing, Int}

    function ComponentInstance(self::AbstractComponentInstance,
                               comp_def::AbstractComponentDef,
                               time_bounds::Tuple{Int,Int},
                               name::Symbol=nameof(comp_def))
        self.comp_name = name
        self.comp_id = comp_id = comp_def.comp_id
        self.comp_path = comp_def.comp_path

        # If first or last is `nothing`, substitute first or last time period
        self.first = @or(comp_def.first, time_bounds[1])
        self.last  = @or(comp_def.last,  time_bounds[2])
    end

    function ComponentInstance(comp_def::AbstractComponentDef,
                               time_bounds::Tuple{Int,Int},
                               name::Symbol=nameof(comp_def))
        self = new()
        return ComponentInstance(self, comp_def, time_bounds, name)
    end
end

@class mutable LeafComponentInstance{TV <: ComponentInstanceVariables,
                                     TP <: ComponentInstanceParameters} <: ComponentInstance begin
    variables::TV                   # TBD: write functions to extract these from type instead of storing?
    parameters::TP
    init::Union{Nothing, Function}
    run_timestep::Union{Nothing, Function}

    function LeafComponentInstance(self::AbstractComponentInstance,
                                   comp_def::AbstractComponentDef,
                                   vars::TV, pars::TP,
                                   time_bounds::Tuple{Int,Int},
                                   name::Symbol=nameof(comp_def)) where
                                {TV <: ComponentInstanceVariables,
                                 TP <: ComponentInstanceParameters}

        # superclass initializer
        ComponentInstance(self, comp_def, time_bounds, name)

        self.variables = vars
        self.parameters = pars

        comp_module = compmodule(self)

        # The try/catch allows components with no run_timestep function (as in some of our test cases)
        # CompositeComponentInstances use a standard method that just loops over inner components.
        # TBD: use FunctionWrapper here?
        function get_func(name)

            #
            # TBD: since this class is no longer a superclass of CompositeComponentInstance
            # this test should be unnecessary. Double-check this though...
            #
            if is_composite(self)
                return nothing
            end

            func_name = Symbol("$(name)_$(nameof(comp_module))_$(self.comp_id.comp_name)")
            try
                getfield(comp_module, func_name)
            catch err
                # @info "Eval of $func_name in module $comp_module failed"
                nothing
            end
        end

        self.init         = get_func("init")
        self.run_timestep = get_func("run_timestep")

        return self
    end

    # Create an empty instance with the given type parameters
    function LeafComponentInstance{TV, TP}() where {TV <: ComponentInstanceVariables, TP <: ComponentInstanceParameters}
        return new{TV, TP}()
    end
end

function LeafComponentInstance(comp_def::AbstractComponentDef, vars::TV, pars::TP,
                               time_bounds::Tuple{Int,Int},
                               name::Symbol=nameof(comp_def)) where
        {TV <: ComponentInstanceVariables, TP <: ComponentInstanceParameters}

    self = LeafComponentInstance{TV, TP}()
    return LeafComponentInstance(self, comp_def, vars, pars, time_bounds, name)
end

# These can be called on CompositeComponentInstances and ModelInstances
compdef(obj::AbstractComponentInstance) = compdef(comp_id(obj))
Base.pathof(obj::AbstractComponentInstance) = obj.comp_path
first_period(obj::AbstractComponentInstance) = obj.first
last_period(obj::AbstractComponentInstance)  = obj.last

"""
Return the ComponentInstanceParameters/Variables exported by the given list of
component instances.
"""
function _comp_instance_vars_pars(comp_def::AbstractCompositeComponentDef,
								  comps::Vector{<: AbstractComponentInstance})
    vdict = Dict([:types => [], :names => [], :values => [], :paths => []])
    pdict = Dict([:types => [], :names => [], :values => [], :paths => []])

    comps_dict = Dict([comp.comp_name => comp for comp in comps])

    vars = ComponentInstanceVariables(Vector{Symbol}(vdict[:names]), Vector{DataType}(vdict[:types]),
                                      Vector{Any}(vdict[:values]), Vector{ComponentPath}(vdict[:paths]))

    pars = ComponentInstanceParameters(Vector{Symbol}(pdict[:names]), Vector{DataType}(pdict[:types]),
                                       Vector{Any}(pdict[:values]), Vector{ComponentPath}(pdict[:paths]))
    return vars, pars
end

@class mutable CompositeComponentInstance <: ComponentInstance begin
    comps_dict::OrderedDict{Symbol, AbstractComponentInstance}
    variables::NamedTuple
    parameters::NamedTuple

    function CompositeComponentInstance(self::AbstractCompositeComponentInstance,
                                        comps::Vector{<: AbstractComponentInstance},
                                        comp_def::AbstractCompositeComponentDef,
                                        time_bounds::Tuple{Int,Int},
                                        variables::NamedTuple,
                                        parameters::NamedTuple,
                                        name::Symbol=nameof(comp_def))

        comps_dict = OrderedDict{Symbol, AbstractComponentInstance}()
        for ci in comps
            comps_dict[ci.comp_name] = ci
        end

        ComponentInstance(self, comp_def, time_bounds, name)
        CompositeComponentInstance(self, comps_dict, variables, parameters)
        return self
    end

    # TBD: Construct vars and params from sub-components
    function CompositeComponentInstance(comps::Vector{<: AbstractComponentInstance},
                                        comp_def::AbstractCompositeComponentDef,
                                        time_bounds::Tuple{Int,Int},
                                        variables::NamedTuple,
                                        parameters::NamedTuple,
                                        name::Symbol=nameof(comp_def))
        CompositeComponentInstance(new(), comps, comp_def, time_bounds, variables, parameters, name)
    end
end

# These methods can be called on ModelInstances as well
components(obj::AbstractCompositeComponentInstance) = values(obj.comps_dict)
has_comp(obj::AbstractCompositeComponentInstance, name::Symbol) = haskey(obj.comps_dict, name)
compinstance(obj::AbstractCompositeComponentInstance, name::Symbol) = obj.comps_dict[name]

is_leaf(ci::LeafComponentInstance) = true
is_leaf(ci::AbstractCompositeComponentInstance) = false
is_composite(ci::AbstractComponentInstance) = !is_leaf(ci)

# ModelInstance holds the built model that is ready to be run
@class ModelInstance <: CompositeComponentInstance begin
    md::ModelDef
end


text:
"""
    Bgen(path; sample_path=nothing, delay_parsing=false)

Read in the Bgen file information: header, list of samples.
Variants and genotypes are read separately.

- `path`: path to the ".bgen" file.
- `sample_path`: path to  ".sample" file, if applicable.
- `idx_path`: path to ".bgi" file, defaults to `path * ".bgi`.
"""
function Bgen(path::AbstractString;
    sample_path = nothing,
    idx_path = isfile(path * ".bgi") ? path * ".bgi" : nothing
    )
    io = open(path)
    fsize = filesize(path)
    # read header
    header = Header(io)
    # read samples
    if sample_path !== nothing
        samples = get_samples(sample_path, header.n_samples)
        # disregard sample names in the header if .sample file is provided
        if header.has_sample_ids
            header_sample_length = read(io, UInt32)
            read(io, header_sample_length)
        end
    elseif header.has_sample_ids
        samples = get_samples(io, header.n_samples)
    else
        samples = get_samples(header.n_samples)
    end

    offset = header.offset + 4 # location of the first variant_ids

    if idx_path !== nothing
        if isfile(idx_path)
            idx = Index(idx_path)
        else
            @error "$idx_path is not a file"
        end
    else
        idx = nothing
    end

    Bgen(io, fsize, header, samples, idx)
end

@inline io(b::Bgen) = b.io
Base.close(b::Bgen) = close(b.io)
@inline fsize(b::Bgen)::Int = b.fsize
@inline samples(b::Bgen) = b.samples
@inline n_samples(b::Bgen)::Int = b.header.n_samples
@inline n_variants(b::Bgen)::Int = b.header.n_variants
const compression_modes = ["None", "Zlib", "Zstd"]
@inline function compression(b::Bgen)
    compression_modes[b.header.compression + 1]
end

"""
    iterator(b::Bgen; offsets=nothing, from_bgen_start=nothing)
Retrieve a variant iterator for `b`.

- If `offsets` is provided, or `.bgen.bgi` is provided and
`from_bgen_start` is `false`, it returns a `VariantIteratorFromOffsets`,
iterating over the list of offsets.
- Otherwise, it returns a `VariantIteratorFromStart`, iterating from the start
of bgen file to the end of it sequentially.
"""
function iterator(b::Bgen; offsets=nothing, from_bgen_start=false)
    if offsets === nothing
        if b.idx === nothing || from_bgen_start
            return VariantIteratorFromStart(b)
        else
            return VariantIteratorFromOffsets(b, BGEN.offsets(b.idx))
        end
    else
        return VariantIteratorFromOffsets(b, offsets)
    end
end

"""
    offset_first_variant(x)
returns the offset of the first variant
"""
@inline function offset_first_variant(x::Bgen)
    return x.header.offset + 4
end

"""
    parse_variants(b::Bgen; offsets=offsets)
Parse variants of the file.
"""
function parse_variants(b::Bgen; offsets=nothing, from_bgen_start=false)
    collect(iterator(b; offsets=offsets, from_bgen_start=from_bgen_start))
end
function parse_variants(v::VariantIterator)
    collect(v)
end

"""
    rsids(vi)
    rsids(b; offsets=nothing, from_bgen_start=false)
Get rsid list of all variants.

Arguments:
- `vi`: a collection of `Variant`s
- `bgen`: `Bgen` object
- `offsets`: offset of each variant to be returned
"""
function rsids(b::Bgen; offsets=nothing, from_bgen_start=false)
    if b.idx !== nothing && offsets === nothing
        return rsids(b.idx)
    else
        vi = iterator(b; offsets=offsets, from_bgen_start=from_bgen_start)
        rsids(vi)
    end
end
function rsids(vi::VariantIterator)
    collect(v.rsid for v in vi)
end

"""
    chroms(vi)
    chroms(bgen; offsets=nothing)
Get chromosome list of all variants.

Arguments:
- `vi`: a collection of `Variant`s
- `bgen`: `Bgen` object
- `offsets`: offset of each variant to be returned
"""
function chroms(b::Bgen; vi=nothing, offsets=nothing, from_bgen_start=false)
    if b.idx !== nothing && offsets === nothing
        return chroms(b.idx)
    else
        vi = iterator(b; offsets=offsets, from_bgen_start=from_bgen_start)
        chroms(vi)
    end
end
function chroms(vi::VariantIterator)
    collect(v.chroms for v in vi)
end

"""
    positions(vi)
    positions(bgen; offsets=nothing)
Get base pair positions of all variants.

Arguments:
- `vi`: a collection of `Variant`s
- `bgen`: `Bgen` object
- `offsets`: offset of each variant to be returned
"""
function positions(b::Bgen; offsets=nothing,
    from_bgen_start=false)::Vector{Int}
    if b.idx !== nothing && offsets === nothing
        return positions(b.idx)
    else
        vi = iterator(b; offsets=offsets, from_bgen_start=from_bgen_start)
        positions(vi)
    end
end
function positions(vi::VariantIterator)
    collect(v.pos for v in vi)
end


text:
#!/usr/bin/env julia

function between(a, b)
    [n for n in a:b]
end

println(between(1, 4)) # [1, 2, 3, 4]
println(between(-2, 2)) # [-2, -1, 0, 1, 2]

max_stars_repo_name: joaopalmeiro/codewars

max_issues_repo_name: joaopalmeiro/codewars

max_forks_repo_name: joaopalmeiro/codewars


text:
### Wrapping elliptic group from ECC
using ECC

struct EllipticGroup{T} <: CyclicGroup
    G::T
    q
    t::Int
end

EllipticGroup{T}(sec,G::EllipticGroup{T}) where T <: S256Point = EllipticGroup(sec2point(sec),G.q,G.t)
EllipticGroup(sec,G::EllipticGroup) = typeof(G)(sec,G)

binary(G::EllipticGroup{T}) where T <: S256Point = point2sec(G.G)
value(G::EllipticGroup{T}) where T <: S256Point = G.G.𝑥.𝑛
security(G::EllipticGroup) = G.t
order(G::EllipticGroup) = G.q

function Scep256k1Group()
    G = ECC.G # G the scep256k1 generator point
    t = 256
    q = ECC.N
    EllipticGroup(G,q,t)
end

^(G::EllipticGroup,n) = EllipticGroup(n*G.G,G.q,G.t)

# When more elliptic groups appears then check for the elements should become important.
import Base.*
*(X::EllipticGroup,Y::EllipticGroup) = EllipticGroup(X.G + Y.G,X.q,X.t)

import Base.==
==(X::EllipticGroup,Y::EllipticGroup) = X.G==Y.G


text:
"""
# module Matrix

# Examples

```jldoctest
julia>
```
"""
module Matrix

export issquare

function issquare(A::AbstractMatrix)::Bool
    m, n = size(A)
    m == n ? true : false
end

end


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

JLLWrappers.@generate_main_file_header("CUTENSOR_CUDA102")
JLLWrappers.@generate_main_file("CUTENSOR_CUDA102", UUID("3383e110-d2c8-5588-847f-51a846eee08b"))
end  # module CUTENSOR_CUDA102_jll


text:
mutable struct Kahan{T <: AbstractFloat}
	sum::T
	c::T
end

Kahan(::Type{T}) where {T <: AbstractFloat} = Kahan{T}(convert(T, 0), convert(T, 0))

#import NumericExtensions.add!
function add!(k::Kahan{T}, x::T) where {T <: AbstractFloat}
	y = x - k.c
	t = k.sum + y
	k.c = (t - k.sum) - y
	k.sum = t
end

import Base.sum
sum(k::Kahan) = k.sum

mutable struct MeanCounter{T <: AbstractFloat}
	n::Int64
	mean::Kahan{T}
end

MeanCounter(::Type{T}) where {T <: AbstractFloat} = MeanCounter{T}(0, Kahan(T))

function add!(m::MeanCounter{T}, x::T) where {T <: AbstractFloat}
	m.n += 1
	add!(m.mean, (x - sum(m.mean)) / m.n)
	return m.mean
end

mean(m::MeanCounter) = sum(m.mean)

export Kahan, MeanCounter, add!, sum, mean


text:
@testset "Read RAWGO" begin
        path = "./test/data/C2FEN02312/scenario_3/"
        raw_path = joinpath(path, "case.raw")
        network = PowerFlowNetwork(path, "RAWGO")
        @test network.baseMVA == 100
        @test size(network.bus, 1) == 2312
        @test nbus(raw_path, "RAWGO") == nbus(network)
        @test nbranch(raw_path, "RAWGO") == nbranch(network)
        @test nbranch(raw_path, "RAWGO", distinct_pair=true) == nbranch(network; distinct_pair=true)
        @test ngen(raw_path, "RAWGO") >= ngen(network)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
myfunc(x::Int, y::Int, z::Int) = x * y * z

function check_getindex(x, ans)
    @test x == ans
    for i in eachindex(x)
        @test (@inferred x[i]) == ans[i]
    end
end

@testset "LazyLayeredArray" begin
    # layer = 1
    data1 = [1,2,3]
    layer1 = LayeredArray(data1)
    check_getindex(@inferred(2 .* layer1)::LazyLayeredArray{1}, 2 * data1)
    check_getindex(@inferred(layer1 ./ 2)::LazyLayeredArray{1}, data1 / 2)
    check_getindex(@inferred(layer1 .+ layer1)::LazyLayeredArray{1}, data1 + data1)
    check_getindex(@inferred(layer1 .- layer1)::LazyLayeredArray{1}, data1 - data1)
    check_getindex(@inferred(layer1 .* layer1)::LazyLayeredArray{1}, data1 .* data1)
    check_getindex(@inferred((layer1 .+ 2*layer1))::LazyLayeredArray{1}, @. data1 + 2*data1)
    check_getindex(@inferred((layer1.*layer1 .+ layer1.*layer1))::LazyLayeredArray{1}, 2 * (data1 .* data1))
    check_getindex(@inferred((layer1.*layer1 .- layer1.*layer1))::LazyLayeredArray{1}, 0 * (data1 .* data1))
    ## with broadcasted object
    check_getindex(@inferred((1 .+ 1) .+ layer1)::LazyLayeredArray{1}, (1 .+ 1) .+ data1)
    check_getindex(@inferred(([1,2,3] .+ 1) .* layer1)::LazyLayeredArray{1}, Ref([1,2,3] .+ 1) .* data1)
    # custom type
    data1′ = [2,4,6]
    layer1′ = MyType(data1′)
    check_getindex(@inferred(2 .* layer1′)::LazyLayeredArray{1}, 2 * data1′)
    check_getindex(@inferred(layer1′ ./ 2)::LazyLayeredArray{1}, data1′ / 2)
    check_getindex(@inferred(layer1′ .+ layer1′)::LazyLayeredArray{1}, data1′ + data1′)
    check_getindex(@inferred(layer1′ .- layer1′)::LazyLayeredArray{1}, data1′ - data1′)
    check_getindex(@inferred(layer1′ .* layer1′')::LazyLayeredArray{1}, data1′ * data1′')
    check_getindex(@inferred((layer1′ .+ 2layer1′))::LazyLayeredArray{1}, @. data1′ + 2*data1′)
    check_getindex(@inferred((layer1′.*layer1′' .+ layer1′.*layer1′'))::LazyLayeredArray{1}, 2 * (data1′ * data1′'))
    check_getindex(@inferred((layer1′.*layer1′' .- layer1′.*layer1′'))::LazyLayeredArray{1}, 0 * (data1′ * data1′'))
    # layer = 1 and 2
    data2 = [4,5,6]
    layer2 = LayeredArray{2}(data2)
    check_getindex(@inferred(layer1 .* layer2)::LazyLayeredArray{2}, copy(Broadcast.broadcasted(*, Ref(data1), data2)))
    check_getindex(@inferred(layer2 .* layer1)::LazyLayeredArray{2}, copy(Broadcast.broadcasted(*, data2, Ref(data1))))
    check_getindex(@inferred(layer1 .+ layer2)::LazyLayeredArray{2}, copy(Broadcast.broadcasted((x,y) -> x .+ y, Ref(data1), data2)))
    check_getindex(@inferred(layer2 .+ layer1)::LazyLayeredArray{2}, copy(Broadcast.broadcasted((x,y) -> x .+ y, data2, Ref(data1))))
    check_getindex(@inferred(layer1 .- layer2)::LazyLayeredArray{2}, copy(Broadcast.broadcasted((x,y) -> x .- y, Ref(data1), data2)))
    check_getindex(@inferred(layer2 .- layer1)::LazyLayeredArray{2}, copy(Broadcast.broadcasted((x,y) -> x .- y, data2, Ref(data1))))
    @test_throws Exception layer1 + layer2
    @test_throws Exception layer2 + layer1
    @test_throws Exception layer1 - layer2
    @test_throws Exception layer2 - layer1
    # layer = 0 and 1
    check_getindex(@inferred(layer1′ .* layer1)::LazyLayeredArray{1}, data1′ .* data1)
    check_getindex(@inferred(layer1 .* layer1′)::LazyLayeredArray{1}, data1 .* data1′)
    check_getindex(@inferred(layer1′ ./ layer1)::LazyLayeredArray{1}, data1′ ./ data1)
    check_getindex(@inferred(layer1 ./ layer1′)::LazyLayeredArray{1}, data1 ./ data1′)
    check_getindex(@inferred(layer1′ .+ layer1)::LazyLayeredArray{1}, data1′ .+ data1)
    check_getindex(@inferred(layer1 .+ layer1′)::LazyLayeredArray{1}, data1 .+ data1′)
    check_getindex(@inferred(layer1′ .- layer1)::LazyLayeredArray{1}, data1′ .- data1)
    check_getindex(@inferred(layer1 .- layer1′)::LazyLayeredArray{1}, data1 .- data1′)
    # layer = 0 and 2
    check_getindex(@inferred(layer1′ .* layer2)::LazyLayeredArray{2}, copy(Broadcast.broadcasted(*, Ref(data1′), data2)))
    check_getindex(@inferred(layer2 .* layer1′)::LazyLayeredArray{2}, copy(Broadcast.broadcasted(*, data2, Ref(data1′))))
    @test_throws Exception layer1′ + layer2
    @test_throws Exception layer2 + layer1′
    @test_throws Exception layer1′ - layer2
    @test_throws Exception layer2 - layer1′
    # custom functions
    check_getindex(@inferred(broadcast(myfunc, layer1, 2, layer1′))::LazyLayeredArray{1}, [myfunc(data1[i], 2, data1′[i]) for i in eachindex(data1, data1′)])
    check_getindex(@inferred(broadcast(myfunc, layer1′, 2, layer1′'))::LazyLayeredArray{1}, [myfunc(data1′[i], 2, data1′[j]) for i in eachindex(data1′), j in eachindex(data1′)])
    check_getindex(@inferred(broadcast(myfunc, layer1, 2, layer2))::LazyLayeredArray{2}, [[myfunc(x1, 2, x2) for x1 in data1] for x2 in data2])
    # other complicated cases
    check_getindex(@inferred(broadcast(sum, layer1 .* layer2))::LazyLayeredArray{2}, [sum([x1 * x2 for x1 in data1]) for x2 in data2])
    check_getindex(@inferred(broadcast(sum, layer1′ .* layer2))::LazyLayeredArray{2}, [sum([x0 * x2 for x0 in data1′]) for x2 in data2])
    # recursive version
    d0 = [layer1′, layer1′, layer1′]
    d1 = [layer1, layer1, layer1]
    x0 = LayeredArray{2}(d0)
    x1 = LayeredArray{2}(d1)
    check_getindex(@inferred(broadcast(sum, @. 2 * x0 * x1))::LazyLayeredArray{2}, [sum([(2 * d0[i][j] * d1[i][j]) for j in 1:length(d0[i])]) for i in 1:length(d0)])
end


text:
function _sum_dimensions(h)
    dim = 0
    for hh in values(h)
        dim += dimension(hh)
    end
    return dim
end


text:
using Markdown
using InteractiveUtils
using Flux: onehot, onecold
# MAIN FUNCTIONS
function get_img_paths(data_dir,extra_format = "")
    """
    Get all images in folder
    """
    img_path = []
    for i in readdir(data_dir)
        if split(i, ".")[end] == "jpg" || split(i, ".")[end] == "png" || split(i, ".")[end] == extra_format
            push!(img_path, joinpath(data_dir, i))
        end
    end
    img_path
end

function load_classes(data_dir)
    """
    - Identify the number of classes in the data. Return an array with all image paths and class labels
    - Does not actually load the images. Just returns paths
    - Assuming heirarchy of classes to be:
        - Main Folder
            - class1
                -img1
                -img2...
            - class2
                -img1
                -img2...
    """
    classes = readdir(main_path)
    @info length(classes), classes
    total_files = sum([length(readdir(main_path*classes[x])) for x in 1:length(classes)] )
    @info "Total images: $total_files"
    labels, paths = Array{String}(undef, total_files), Array{String}(undef, total_files)
    current_index = 1
    for path in classes
        current_path = readdir(main_path*path)
        current_length = length(current_path)

        labels[current_index:current_length+current_index-1] .= path
        # @info paths[current_index:current_length]
        paths[current_index:current_length+current_index-1] = main_path*path.*current_path
        current_index += current_length
    end

    labels, paths
end

function labelFromPattern(paths, labeller)
    """
    - Takes paths generated from load_classes as input
    - Takes a function defining a custom label
    - eg: labeller(x) = split(x, "/")[-1]
    - Returns labels based on function specified
    """
    labeller.(paths)
end

"""
- A function to show a custom labeller function as an example
- This splits the path based on "/" and returns the lower case of the last element
"""
examplelabeller(x) = lowercase(split(x,"/")[end])

"""
- Input path -> Read an image -> Convert to array -> Resize to given dimensions -> Convert to float
"""
loadim(path, dims = (128, 128)) = float.(permutedims(channelview(imresize(load(path), dims...)), (2, 3, 1)))

function load_train_imgs(labels, train_img_paths, dims = (128,128))
    """
    - Load all images
    - Resize them to given size
    - one hot encode labels
    """
    train_imgs = zeros()
    imgs_loaded = 0
    Threads.@spawn for i in train_img_paths
        push!(train_imgs, loadim(i, dims...))
        imgs_loaded += 1
        if imgs_loaded % 1000 == 0
            @info "$imgs_loaded Images have been loaded"
        end
    end

end

# EXAMPLE WORKFLOW

main_path = "/home/subhaditya/Desktop/Datasets/asl-downloaded/"
labels, paths = load_classes(main_path)
@info labelFromPattern(paths, examplelabeller)[1:3]


text:
print("Testing ABCDQR_simula: ")
###################################
nt = 1000
seed = 99
m = ABCDQR_simula(nt,seed)
y = m["y"]

if sum(y.^2) > 1e-6
    println("OK")
end


text:
function run_small_instances()
    mdp_params = Dict(pairs( (n_edges = 3, c_init = 2, demand = Float64[1,1,1], selling_horizon_end = [40,45,50], actions = 15:5:90, objective=:revenue)))
    #mdp_params = Dict(pairs( (n_edges = 2, c_init = 1, demand = Float64[1,1], selling_horizon_end = [20,25], actions = 15:5:30)))
    # mcts_params = Dict(solver= MCTSSolver, n_iterations=1000, depth=30, exploration_constant=40.0, reuse_tree=true)
    mcts_params = Dict(pairs( (solver= DPWSolver, n_iterations=50, depth=30, exploration_constant=40.0, enable_state_pw = true, keep_tree=true, show_progress=false)))
    exp_params = Dict(pairs( (n_runs = 20, vi=true, save=:stats) ))
    params = Dict(:mdp=>mdp_params, :mcts=>mcts_params, :exp=>exp_params)

end


text:
abstract type AbstractPiecewise{T} end

domain(p::AbstractPiecewise) = (first(nodes(p)), last(nodes(p)))
function domain(p::AbstractPiecewise, i)
    1 <= i <= length(p) || throw(BoundsError(p, i))
    n = nodes(p)
    @inbounds begin
        return (n[i], n[i+1])
    end
end

Base.eltype(t::AbstractPiecewise{T}) where {T} = eltype(T)
Base.eltype(::Type{<:AbstractPiecewise{T}}) where {T} = eltype(T)

# Generic implementation: arbitrary list of elements on the intervals between the nodes
struct Piecewise{T,S<:AbstractVector{<:Real}} <: AbstractPiecewise{T}
    nodes::S
    elements::Vector{T}
    function Piecewise(nodes::S, elements::Vector{T}) where {S<:AbstractVector,T}
        @assert length(nodes) == length(elements) + 1
        @assert eltype(nodes) <: Real
        if nodes isa AbstractRange
            @assert step(nodes) > 0
        else
            @assert issorted(nodes)
        end
        return new{T, S}(nodes, elements)
    end
end

nodes(p::Piecewise) = p.nodes
elements(p::Piecewise) = p.elements
nodevalues(p::Piecewise) = Base.Generator(p, nodes(p))

Base.length(p::Piecewise) = length(p.elements)

@inline Base.getindex(p::Piecewise, i) = getindex(p.elements, i)
@inline Base.setindex!(p::Piecewise, v, i) = setindex!(p.elements, v, i)

function (p::Piecewise)(x)
    x₀ = first(nodes(p))
    x₁ = last(nodes(p))
    (x < x₀ || x > x₁) && throw(DomainError())
    x == x₀ && return p[1](x)
    x == x₁ && return p[length(p)](x)
    if nodes isa AbstractRange
        ir = (x-x₀)/step(nodes(p))
        i = floor(Int, ir)
        return p[i](x)
    else
        j = findlast(<=(x), nodes(p))
        @assert j !== nothing && j !== length(nodes(p))
        return p[j](x)
    end
end

Base.copy(p::Piecewise) = Piecewise(copy(nodes(p)), map(copy, elements(p)))

Base.similar(p::Piecewise) = Piecewise(nodes(p), map(similar, elements(p)))

function LinearAlgebra.rmul!(p::Piecewise, α)
    for i = 1:length(p)
        rmul!(p[i], α)
    end
    return p
end

function LinearAlgebra.lmul!(α, p::Piecewise)
    for i = 1:length(p)
        lmul!(α, p[i])
    end
    return p
end

function LinearAlgebra.mul!(pdst::Piecewise, α, psrc::Piecewise)
    @assert nodes(pdst) == nodes(psrc)
    for i = 1:length(psrc)
        mul!(pdst[i], α, psrc[i])
    end
    return pdst
end

function LinearAlgebra.mul!(pdst::Piecewise, psrc::Piecewise, α)
    @assert nodes(pdst) == nodes(psrc)
    for i = 1:length(psrc)
        mul!(pdst[i], psrc[i], α)
    end
    return pdst
end

function LinearAlgebra.axpy!(α, px::Piecewise, py::Piecewise)
    @assert nodes(px) == nodes(py)
    for i = 1:length(px)
        axpy!(α, px[i], py[i])
    end
    return py
end
function LinearAlgebra.axpby!(α, px::Piecewise, β, py::Piecewise)
    @assert nodes(px) == nodes(py)
    for i = 1:length(px)
        axpby!(α, px[i], β, py[i])
    end
    return py
end

Base.zero(p::Piecewise) = Piecewise(nodes(p), zero.(p.elements))
Base.one(p::Piecewise) = Piecewise(nodes(p), one.(p.elements))

function Base.:+(p1::AbstractPiecewise, p2::AbstractPiecewise)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), [p1[i] + p2[i] for i = 1:length(p1)])
end

function Base.:-(p1::AbstractPiecewise, p2::AbstractPiecewise)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), [p1[i] - p2[i] for i = 1:length(p1)])
end

Base.:-(p::AbstractPiecewise) = Piecewise(nodes(p), map(-, elements(p)))

Base.:*(p::AbstractPiecewise, a) = mul!(similar(p), p, a)
Base.:*(a, p::AbstractPiecewise) = mul!(similar(p), a, p)
Base.:/(p::AbstractPiecewise, a) = mul!(similar(p), p, inv(a))
Base.:\(a, p::AbstractPiecewise) = mul!(similar(p), inv(a), p)

for f in (:conj, :adjoint, :transpose, :real, :imag)
    @eval Base.$f(p::AbstractPiecewise) = Piecewise(nodes(p), map($f, elements(p)))
end

function Base.:*(p1::AbstractPiecewise, p2::AbstractPiecewise; kwargs...)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), [*(p1[i], p2[i]; kwargs...) for i = 1:length(p1)])
end

function LinearAlgebra.mul!(p::Piecewise, p1::AbstractPiecewise, p2::AbstractPiecewise,
                                α = true, β = false)
    @assert nodes(p) == nodes(p1) == nodes(p2)
    for i = 1:length(p)
        mul!(p[i], p1[i], p2[i], α, β)
    end
    return p
end

differentiate(p::AbstractPiecewise) = Piecewise(nodes(p), map(differentiate, elements(p)))
function integrate(p::AbstractPiecewise, interval = domain(p))
    @assert interval == domain(p)
    s = integrate(p[1], (nodes(p)[1], nodes(p)[2]))
    for i = 2:length(p)
        s += integrate(p[i], (nodes(p)[i], nodes(p)[i+1]))
    end
    return s
end

function localdot(p1::AbstractPiecewise, p2::AbstractPiecewise)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), [localdot(p1[i], p2[i]) for i = 1:length(p1)])
end

LinearAlgebra.tr(p::AbstractPiecewise) = Piecewise(nodes(p), map(tr, elements(p)))


text:
function Base.show(io::IO, b::T) where {T<:Union{BasicBlock, CompiledBasicBlock}}

    sep() = println(io, "---------------")
    purity_str(instruction) = is_pure(instruction) ? " Pure " : "Impure"

    function print_with_padding(num, max_number)
        lpad("$num", ndigits(max_number))
    end

    println(io, "BasicBlock")
    sep()
    println(io, "")
    println(io, "Instructions")
    sep()
    i_count = length(b.instructions)
    for (i, entry) in enumerate(b.instructions)
        println(io, "$(print_with_padding(i, i_count)) | $(purity_str(entry)) | $(entry.call)")
    end
    sep()
    println(io, "")
    println(io, "Inputs")
    sep()
    for (i, sym) in enumerate(b.inputs)
        println(io, "$i | $sym")
    end
    sep()
    println(io, "")
    println(io, "Variables - Local Mutables")
    sep()
    for v in b.mutable_locals
        println(io, "$v")
    end
    sep()
    println(io, "")
    println(io, "Variables - Local Immutables")
    sep()
    for (sym, value) in b.immutable_locals
        println(io, "$sym | $value")
    end
    sep()
    println(io, "")
    println(io, "Outputs")
    sep()
    for (i, sym) in enumerate(b.outputs)
        println(io, "$i | $sym")
    end
    sep()
end

function Base.show(io::IO, c::Const{T}) where {T}
    print(io, "$(c.val)::$T")
end

Base.show(io::IO, v::Variable) = print(io, "$(v.symbol)")
Base.show(io::IO, ssa::SSAValue) = print(io, "%$(ssa.address)")
Base.show(io::IO, call::CallUnary) = print(io, "call($(call.op), $(call.arg))")
Base.show(io::IO, call::CallBinary) = print(io, "call($(call.op), $(call.arg1), $(call.arg2))")
Base.show(io::IO, call::Call3Args) = print(io, "call($(call.op), $(call.arg1), $(call.arg2), $(call.arg3))")
Base.show(io::IO, call::Call4Args) = print(io, "call($(call.op), $(call.arg1), $(call.arg2), $(call.arg3), $(call.arg4))")
Base.show(io::IO, call::Call5Args) = print(io, "call($(call.op), $(call.arg1), $(call.arg2), $(call.arg3), $(call.arg4), $(call.arg5))")
Base.show(io::IO, call::Call6Args) = print(io, "call($(call.op), $(call.arg1), $(call.arg2), $(call.arg3), $(call.arg4), $(call.arg5), $(call.arg6))")
Base.show(io::IO, op::Op) = print(io, "$(op.op)")
Base.show(io::IO, op::Assignment) = print(io, "$(op.lhs) = $(op.rhs)")


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame);

# Define the Stan language model

stan10_4sl = "
data {
    int N;
    int<lower=0, upper=1> L[N];
    vector[N] P;
    vector[N] C;

    int<lower=1, upper=N> N_chimps;
    int<lower=1, upper=N_chimps> chimp[N];
}
parameters {
    real a_chimp[N_chimps];
    real bp;
    real bpc;
}
model {
    vector[N] p;
    target += normal_lpdf(a_chimp | 0, 10);
    target += normal_lpdf(bp | 0, 10);
    target += normal_lpdf(bpc | 0, 10);
    for (i in 1:N) p[i] = a_chimp[chimp[i]] + (bp + bpc * C[i]) * P[i];
    target += binomial_logit_lpmf(L | 1, p);
}
generated quantities {
    vector[N] log_lik;
    {
        vector[N] p;
        for(n in 1:N) {
          p[n] = a_chimp[chimp[n]] + (bp + bpc * C[n]) * P[n];
          log_lik[n] = binomial_logit_lpmf(L[n] | 1, p[n]);
        }
    }
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m10_4sl = SampleModel("m10.4sl", stan10_4sl);

# Input data for cmdstan

m10_4_data = Dict("N" => size(df, 1), "N_chimps" => length(unique(df[!, :actor])),
"chimp" => df[!, :actor], "L" => df[!, :pulled_left],
"P" => df[!, :prosoc_left], "C" => df[!, :condition]);

# Sample using cmdstan

rc10_4sl = stan_sample(m10_4sl, data=m10_4_data);

# Result rethinking

rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
bp    0.84 0.26  0.43  1.26  2271    1
bpC  -0.13 0.29 -0.59  0.34  2949    1

a[1] -0.74 0.27 -1.16 -0.31  3310    1
a[2] 10.88 5.20  4.57 20.73  1634    1
a[3] -1.05 0.28 -1.52 -0.59  4206    1
a[4] -1.05 0.28 -1.50 -0.60  4133    1
a[5] -0.75 0.27 -1.18 -0.32  4049    1
a[6]  0.22 0.27 -0.22  0.65  3877    1
a[7]  1.81 0.39  1.22  2.48  3807    1
";

# Update sections

if success(rc10_4sl)
    part10_4sl = read_samples(m10_4sl, :particles)
    part10_4sl |> display
end


text:
"""
`prev_quarter(q::TimeType = now())`

Returns Date identifying last day of the previous quarter
"""
function prev_quarter(q::TimeType = now())
    Dates.Date(lastdayofquarter(firstdayofquarter(q)-Dates.Day(1)))
end

"""
`next_quarter(q::TimeType = now())`

Returns Date identifying last day of the next quarter
"""
function next_quarter(q::TimeType = now())
    Dates.Date(lastdayofquarter(lastdayofquarter(q)+Dates.Day(1)))
end

"""
`get_quarter_ends(start_date::Date,end_date::Date)`

Returns an Array of quarter end dates between `start_date` and `end_date`.
"""
function get_quarter_ends(start_date::Dates.Date,end_date::Dates.Date)
    map(lastdayofquarter, collect(start_date:Dates.Month(3):end_date))
end

"""
`quartertodate(string::String)`

Convert `string` in the form "YYqX", "YYYYqX", or "YYYY-qX" to a Date of the end of
the indicated quarter. "X" is in `{1,2,3,4}` and the case of "q" is ignored.
"""
function quartertodate(string::String)
    if occursin(r"^[0-9]{2}[qQ][1-4]$", string)
        year = "20"*string[1:2]
        quarter = string[end]
    elseif occursin(r"^[0-9]{4}[qQ][1-4]$", string)
        year = string[1:4]
        quarter = string[end]
    elseif occursin(r"^[0-9]{4}-[qQ][1-4]$", string)
        year = string[1:4]
        quarter = string[end]
    else
        throw(Meta.ParseError("Invalid format: $string"))
    end

    year = parse(Int, year)
    quarter = parse(Int, quarter)
    month = 3*quarter
    day = 1

    return lastdayofquarter(Date(year, month, day))
end

"""
`datetoquarter(date::Date)`

Convert `string` in the form "YYqX", "YYYYqX", or "YYYY-qX" to a Date of the end of
the indicated quarter. "X" is in `{1,2,3,4}` and the case of "q" is ignored.

Return an integer from the set `{1,2,3,4}`, corresponding to one of the quarters in a year given a Date object.
"""
function datetoquarter(date::Dates.Date)
    month = Dates.month(date)
    if month in 1:3
        return 1
    elseif month in 4:6
        return 2
    elseif month in 7:9
        return 3
    elseif month in 10:12
        return 4
    else
        throw("Must provide a date object with a valid month (in 1:12)")
    end
end

"""
```
function vinttodate(vint)
```

Return the string given by data_vintage(m), which is in the format YYYYMMDD, to a Date object.
"""
function vinttodate(vint::String)
    year  = Meta.parse("20"*vint[1:2])
    month = Meta.parse(vint[3:4])
    day   = Meta.parse(vint[5:6])
    return Date(year, month, day)
end

"""
`subtract_quarters(t1::Date, t0::Date)`

Compute the number of quarters between t1 and t0, including t0 and excluding t1.
"""
function subtract_quarters(t1::Dates.Date, t0::Dates.Date)
    days = t1 - t0
    quarters = round(days.value / 365.25 * 4.0)
    return convert(Int, quarters)
end

"""
`format_dates!(col, df)`

Change column `col` of dates in `df` from String to Date, and map any dates given in the
interior of a quarter to the last day of the quarter.
"""
function format_dates!(col::Symbol, df::DataFrame)
    df[col] = Dates.Date.(df[col])
    map!(lastdayofquarter, df[col], df[col])
end

"""
```
missing2nan(a::Array)
```

Convert all elements of Union{X, Missing.Missing} or Missing.Missing to type Float64.
"""
function missing2nan(a::Array)
    a_new = Array{Float64}(undef, size(a))
    for i in eachindex(a)
        a_new[i] = ismissing(a[i]) ? NaN : a[i]
    end
    return a_new
end

# Temp function to ensure proper transition into 0.7
function nan2missing!(df::DataFrame)
    for col in names(df)
        if col != :date
            df[col] = convert(Vector{Union{Missing, Float64}}, df[col])
            df[col] = replace(x -> isnan(x) ? missing : x, df[col])
        end
    end
end

"""
```
na2nan!(df::DataFrame)
```

Convert all NAs in a DataFrame to NaNs.
"""
function na2nan!(df::DataFrame)
    for col in names(df)
        if typeof(df[col])==Vector{Date}
            nothing
        else
            df[ismissing.(df[col]), col] = NaN
        end
    end
end

"""
```
na2nan!(df::Array)
```

Convert all NAs in an Array to NaNs.
"""
function na2nan!(v::Array)
    for i = 1:length(v)
        v[i] = isnan(v[i]) ?  NaN : v[i]
    end
end

"""
```
missing_cond_vars!(m, df; cond_type = :none)
```

Make conditional period variables not in `cond_semi_names(m)` or
`cond_full_names(m)` missing if necessary.
"""
function missing_cond_vars!(m::AbstractDSGEModel, df::DataFrame; cond_type::Symbol = :none)
    if cond_type in [:semi, :full]
        # Get appropriate
        cond_names = if cond_type == :semi
            cond_semi_names(m)
        elseif cond_type == :full
            cond_full_names(m)
        end

        # Make non-conditional variables missing
        cond_names_missing = setdiff(names(df), [cond_names; :date])
        for var_name in cond_names_missing
            df[var_name] = convert(Vector{Union{Missing, eltype(df[var_name])}}, df[var_name])
            df[df[:date] .>= date_forecast_start(m), var_name] = missing
        end

        # Warn if any conditional variables are missing
        for var in cond_names
            if any(ismissing.(df[df[:date] .>= date_forecast_start(m), var]))
                @warn "Missing some conditional observations for " * string(var)
            end
        end
    end
end

"""
```
get_data_filename(m, cond_type)
```

Returns the data file for `m`, which depends on `data_vintage(m)`, and if
`cond_type in [:semi, :full]`, also on `cond_vintage(m)` and `cond_id(m)`.
"""
function get_data_filename(m::AbstractDSGEModel, cond_type::Symbol)
    filestrings = ["data"]

    # If writing conditional data, append conditional vintage and ID to filename
    if cond_type in [:semi, :full]
        push!(filestrings, "cdid=" * lpad(string(cond_id(m)), 2, string(0)))
        push!(filestrings, "cdvt=" * cond_vintage(m))
    end

    push!(filestrings, "dsid=" * lpad(string(data_id(m)), 2, string(0)))
    push!(filestrings, "vint=" * data_vintage(m))
    filename = join(filestrings, "_")

    return inpath(m, "data", filename * ".csv")
end

"""
```
iterate_quarters(start::Date, quarters::Int)
```

Returns the date corresponding to `start` + `quarters` quarters.

### Inputs
- `start`: starting date
- `quarters`: number of quarters to iterate forward or backward
"""
function iterate_quarters(start::Dates.Date, quarters::Int)

    next = start
    if quarters < 0
        for n = 1:-quarters
            next = Dates.toprev(next) do x
                Dates.lastdayofquarter(x) == x
            end
        end
    elseif quarters > 0
        for n = 1:quarters
            next = Dates.tonext(next) do x
                Dates.lastdayofquarter(x) == x
            end
        end
    end

    next
end

"""
```
reconcile_column_names(a::DataFrame, b::DataFrame)
```

adds columns of missings to a and b so that both have
the same set of column names.
"""
function reconcile_column_names(a::DataFrame, b::DataFrame)
    new_a_cols = setdiff(names(b), names(a))
    new_b_cols = setdiff(names(a), names(b))
    for col in new_a_cols
        a[col] = fill(missing, size(a, 1))
    end
    for col in new_b_cols
        b[col] = fill(missing, size(b, 1))
    end
    return a, b
end


text:
module Client
export query

using JSON, HTTP
using ..Types

function query(client::GraphQLClient, options::QueryOptions = QueryOptions())
  response = HTTP.request("POST", client.url, [], options.query)
  return JSON.parse(String(response.body))
end

end


text:
const testGridLength = 100
const tolerance = 1e-10
testgridArchetype = range(0, stop=1, length=testGridLength)
testGrids = [testgridArchetype, collect(testgridArchetype)]

func1(x,t) = first(promote(zero(x), zero(t)))
func2(x,t) = first(promote(one(x), one(t)))
func3(x,t) = x*t
func4(x,t) = sin(x)-cos(t)
testFuncs = [func1, func2, func3, func4]

@testset "space_trace_*, grid=$testgrid" for testgrid in testGrids
    @testset "Positivity, testFunction $func" for func in testFuncs
        U = [func(x,t) for x in testgrid, t in testgrid]
        @test WN.space_trace_norm(U, testgrid, testGridLength, testGridLength-1) >= 0

        v = [func(x, testgrid[testGridLength-1]) for x in testgrid]
        @test WN.space_trace_dist(U, v, testgrid, testGridLength, testGridLength-1) >= 0
    end
end

@testset "space_dist, grid=$testgrid" for testgrid in testGrids
    @testset "Identity, function=$func" for func in testFuncs
        U = [func(x,t) for x in testgrid, t in testgrid]
        @test isapprox(WN.space_dist(U, U, testgrid, testgrid[2] - testgrid[1]), 0)
    end
end


text:
using Test

using LinearAlgebra
using SymbolicWedderburn
using SymbolicWedderburn.FiniteFields
using LinearAlgebra

using GroupsCore
using PermutationGroups
using Cyclotomics

include("smallgroups.jl")

@testset "Characters" begin
    import SymbolicWedderburn.Characters
    include("gf.jl")
    include("eigenspacedecomposition.jl")
    include("ccmatrix.jl")
    include("dixon.jl")
    include("characters.jl")
end
include("projections.jl")
include("sa_basis.jl")
include("action_permutation.jl")
include("action_linear.jl")
include("action_dihedral.jl")

if VERSION >= v"1.6.0" && !haskey(ENV, "CI")
    @testset "Examples" begin
        using Pkg
        Pkg.activate(joinpath(@__DIR__, "..", "examples"))
        Pkg.instantiate()
        include("../examples/ex_C2_linear.jl")
        include("../examples/ex_S4.jl")
        include("../examples/ex_motzkin.jl")
        include("../examples/ex_robinson_form.jl")
    end
end


text:
function norm_int(cur_func::Function)
  cur_integral = nothing

  try
    cur_integral = quadgk(
      cur_func, integral_zero, integral_one
    )[1]
  catch cur_error
    isa(cur_error, DomainError) || rethrow(cur_error)
    cur_integral = NaN
  end

  cur_integral
end

function norm_int(cur_func::SymEngine.Basic)
  [1, cur_func]
end