text:
"""
Provides a collection of utility functions and types that are used in other submodules.
"""
module Utilities

using Base.Meta, Compat

# Logging output.

const __log__ = Ref(true)
"""
    logging(flag::Bool)

Enable or disable logging output for [`log`](@ref) and [`warn`](@ref).
"""
logging(flag::Bool) = __log__[] = flag

"""
Format and print a message to the user.
"""
log(msg) = __log__[] ? print_with_color(:magenta, STDOUT, "Documenter: ", msg, "\n") : nothing

# Print logging output to the "real" STDOUT.
function log(doc, msg)
    __log__[] && print_with_color(:magenta, STDOUT, "Documenter: ", msg, "\n")
    return nothing
end

debug(msg) = print_with_color(:green, " ?? ", msg, "\n")

"""
    warn(file, msg)
    warn(msg)

Format and print a warning message to the user. Passing a `file` will include the filename
where the warning was raised.
"""
function warn(file, msg)
    if __log__[]
        msg = string(" !! ", msg, " [", file, "]\n")
        print_with_color(:red, STDOUT, msg)
    else
        nothing
    end
end
warn(msg) = __log__[] ? print_with_color(:red, STDOUT, " !! ", msg, "\n") : nothing

function warn(file, msg, err, ex, mod)
    if __log__[]
        warn(file, msg)
        print_with_color(:red, STDOUT, "\nERROR: $err\n\nexpression '$ex' in module '$mod'\n\n")
    else
        nothing
    end
end

function warn(doc, page, msg, err)
    file = page.source
    print_with_color(:red, STDOUT, " !! Warning in $(file):\n\n$(msg)\n\nERROR: $(err)\n\n")
end

# Nullable regex matches.

wrapnothing(T, ::Void) = Nullable{T}()
wrapnothing(T, value)  = Nullable(value)

nullmatch(r::Regex, str::AbstractString) = wrapnothing(RegexMatch, match(r, str))

getmatch(n::Nullable{RegexMatch}, i) = get(n)[i]

# Directory paths.

"""
Returns the current directory.
"""
function currentdir()
    d = Base.source_dir()
    d === nothing ? pwd() : d
end

"""
Returns the path to the Documenter `assets` directory.
"""
assetsdir() = normpath(joinpath(dirname(@__FILE__), "..", "..", "assets"))

cleandir(d::AbstractString) = (isdir(d) && rm(d, recursive = true); mkdir(d))

"""
Find the path of a file relative to the `source` directory. `root` is the path
to the directory containing the file `file`.

It is meant to be used with `walkdir(source)`.
"""
srcpath(source, root, file) = normpath(joinpath(relpath(root, source), file))

# Slugify text.

"""
Slugify a string into a suitable URL.
"""
function slugify(s::AbstractString)
    s = replace(s, r"\s+", "-")
    s = replace(s, r"^\d+", "")
    s = replace(s, r"&", "-and-")
    s = replace(s, r"[^\p{L}\p{P}\d\-]+", "")
    s = strip(replace(s, r"\-\-+", "-"), '-')
end
slugify(object) = string(object) # Non-string slugifying doesn't do anything.

# Parse code blocks.

"""
Returns a vector of parsed expressions and their corresponding raw strings.

The keyword argument `skip = N` drops the leading `N` lines from the input string.
"""
function parseblock(code::AbstractString, doc, page; skip = 0)
    # Drop `skip` leading lines from the code block. Needed for deprecated `{docs}` syntax.
    code = string(code, '\n')
    code = last(split(code, '\n', limit = skip + 1))
    # Check whether we have windows-style line endings.
    local offset = contains(code, "\n\r") ? 2 : 1
    local strlen = length(code)
    local results = []
    local cursor = 1
    while cursor < strlen
        # Check for keywords first since they will throw parse errors if we `parse` them.
        local line = match(r"^(.+)$"m, SubString(code, cursor, strlen)).captures[1]
        local keyword = Symbol(strip(line))
        (ex, ncursor) =
            if haskey(Docs.keywords, keyword)
                (QuoteNode(keyword), cursor + length(line) + offset)
            else
                try
                    parse(code, cursor)
                catch err
                    Utilities.warn(doc, page, "Failed to parse expression.", err)
                    break
                end
            end
        push!(results, (ex, SubString(code, cursor, ncursor - 1)))
        cursor = ncursor
    end
    results
end
isassign(x) = isexpr(x, :(=), 2) && isa(x.args[1], Symbol)

# Checking arguments.

"""
Prints a formatted warning to the user listing unrecognised keyword arguments.
"""
function check_kwargs(kws)
    isempty(kws) && return
    out = IOBuffer()
    println(out, "Unknown keywords:\n")
    for (k, v) in kws
        println(out, "  ", k, " = ", v)
    end
    warn(takebuf_string(out))
end

# Finding submodules.

typealias ModVec Union{Module, Vector{Module}}

"""
Returns the set of submodules of a given root module/s.
"""
function submodules(modules::Vector{Module})
    out = Set{Module}()
    for each in modules
        union!(out, submodules(each))
    end
    out
end
function submodules(root::Module, out = Set([root]))
    for name in names(root, true)
        if isdefined(root, name) && !isdeprecated(root, name)
            object = getfield(root, name)
            if isvalidmodule(root, object)
                push!(out, object)
                submodules(object, out)
            end
        end
    end
    out
end
isvalidmodule(a::Module, b::Module) = a !== b && b !== Main
isvalidmodule(a, b)                 = false

# Compat for `isdeprecated` which does not exist in Julia 0.4.
if isdefined(Base, :isdeprecated)
    isdeprecated(m, s) = Base.isdeprecated(m, s)
else
    isdeprecated(m, s) = ccall(:jl_is_binding_deprecated, Cint, (Any, Any), m, s) != 0
end

## objects
## =======

import Base.Docs: Binding

if VERSION < v"0.5.0-dev"
    @eval function Base.call(::Type{Binding}, m::Module, v::Symbol)
        m = module_name(m) === v ? module_parent(m) : m
        m = Base.binding_module(m, v)
        $(Expr(:new, :Binding, :m, :v))
    end
end

"""
Represents an object stored in the docsystem by its binding and signature.
"""
immutable Object
    binding   :: Binding
    signature :: Type

    function Object(b::Binding, signature::Type)
        m = module_name(b.mod) === b.var ? module_parent(b.mod) : b.mod
        new(Binding(m, b.var), signature)
    end
end

function splitexpr(x::Expr)
    isexpr(x, :macrocall) ? splitexpr(x.args[1]) :
    isexpr(x, :.)         ? (x.args[1], x.args[2]) :
    error("Invalid @var syntax `$x`.")
end
splitexpr(s::Symbol) = :(current_module()), quot(s)
splitexpr(other)     = error("Invalid @var syntax `$other`.")

"""
    object(ex, str)

Returns a expression that, when evaluated, returns an [`Object`](@ref) representing `ex`.
"""
function object(ex::Union{Symbol, Expr}, str::AbstractString)
    binding   = Expr(:call, Binding, splitexpr(Docs.namify(ex))...)
    signature = Base.Docs.signature(ex)
    isexpr(ex, :macrocall, 1) && !endswith(str, "()") && (signature = :(Union{}))
    Expr(:call, Object, binding, signature)
end

function object(qn::QuoteNode, str::AbstractString)
    if haskey(Base.Docs.keywords, qn.value)
        binding = Expr(:call, Binding, Main, qn)
        Expr(:call, Object, binding, Union{})
    else
        error("'$(qn.value)' is not a documented keyword.")
    end
end

function Base.print(io::IO, obj::Object)
    print(io, obj.binding)
    print_signature(io, obj.signature)
end
print_signature(io::IO, signature::Union) = nothing
print_signature(io::IO, signature)        = print(io, '-', signature)

## docs
## ====

"""
    docs(ex, str)

Returns an expression that, when evaluated, returns the docstrings associated with `ex`.
"""
function docs end

# Macro representation changed between 0.4 and 0.5.
if VERSION < v"0.5-"
    function docs(ex::Union{Symbol, Expr}, str::AbstractString)
        :(Base.Docs.@doc $ex)
    end
else
    function docs(ex::Union{Symbol, Expr}, str::AbstractString)
        isexpr(ex, :macrocall, 1) && !endswith(rstrip(str), "()") && (ex = quot(ex))
        :(Base.Docs.@doc $ex)
    end
end
docs(qn::QuoteNode, str::AbstractString) = :(Base.Docs.@doc $(qn.value))

"""
Returns the category name of the provided [`Object`](@ref).
"""
doccat(obj::Object) = startswith(string(obj.binding.var), '@') ?
    "Macro" : doccat(obj.binding, obj.signature)

function doccat(b::Binding, ::Union)
    if b.mod === Main && haskey(Base.Docs.keywords, b.var)
        "Keyword"
    elseif startswith(string(b.var), '@')
        "Macro"
    else
        doccat(getfield(b.mod, b.var))
    end
end

doccat(b::Binding, ::Type)  = "Method"

doccat(::Function) = "Function"
doccat(::DataType) = "Type"
doccat(::Module)   = "Module"
doccat(::Any)      = "Constant"

"""
    filterdocs(doc, modules)

Remove docstrings from the markdown object, `doc`, that are not from one of `modules`.
"""
function filterdocs(doc::Markdown.MD, modules::Set{Module})
    if isempty(modules)
        # When no modules are specified in `makedocs` then don't filter anything.
        Nullable(doc)
    else
        if haskey(doc.meta, :module)
            doc.meta[:module] ∈ modules ? Nullable(doc) : Nullable{Markdown.MD}()
        else
            if haskey(doc.meta, :results)
                out = []
                results = []
                for (each, result) in zip(doc.content, doc.meta[:results])
                    r = filterdocs(each, modules)
                    if !isnull(r)
                        push!(out, get(r))
                        push!(results, result)
                    end
                end
                if isempty(out)
                    Nullable{Markdown.MD}()
                else
                    md = Markdown.MD(out)
                    md.meta[:results] = results
                    Nullable(md)
                end
            else
                out = []
                for each in doc.content
                    r = filterdocs(each, modules)
                    isnull(r) || push!(out, get(r))
                end
                isempty(out) ? Nullable{Markdown.MD}() : Nullable(Markdown.MD(out))
            end
        end
    end
end
# Non-markdown docs won't have a `.meta` field so always just accept those.
filterdocs(other, modules::Set{Module}) = Nullable(other)

"""
Does the given docstring represent actual documentation or a no docs error message?
"""
nodocs(x)      = contains(stringmime("text/plain", x), "No documentation found.")
nodocs(::Void) = false

header_level{N}(::Markdown.Header{N}) = N

# Finding URLs -- based partially on code from the main Julia repo in `base/methodshow.jl`.

function url(repo, file)
    file = abspath(file)
    remote = getremote(dirname(file))
    isempty(repo) && (repo = "https://github.com/$remote/tree/{commit}{path}")
    commit, root = cd(dirname(file)) do
        readchomp(`git rev-parse HEAD`), readchomp(`git rev-parse --show-toplevel`)
    end
    if startswith(file, root)
        _, path = split(file, root; limit = 2)
        repo = replace(repo, "{commit}", commit)
        repo = replace(repo, "{path}", path)
        Nullable{Compat.String}(repo)
    else
        Nullable{Compat.String}()
    end
end

url(remote, repo, doc) = url(remote, repo, doc.data[:module], doc.data[:path], linerange(doc))

# Correct file and line info only available from this version onwards.
if VERSION >= v"0.5.0-dev+3442"
    function url(remote, repo, mod, file, linerange)
        remote = getremote(dirname(file))
        isabspath(file) && isempty(remote) && isempty(repo) && return Nullable{Compat.String}()
        # Format the line range.
        line = format_line(linerange)
        # Macro-generated methods such as those produced by `@deprecate` list their file as
        # `deprecated.jl` since that is where the macro is defined. Use that to help
        # determine the correct URL.
        if inbase(mod) || !isabspath(file)
            base = "https://github.com/JuliaLang/julia/tree"
            dest = "base/$file#$line"
            Nullable{Compat.String}(
                if isempty(Base.GIT_VERSION_INFO.commit)
                    "$base/v$VERSION/$dest"
                else
                    commit = Base.GIT_VERSION_INFO.commit
                    "$base/$commit/$dest"
                end
            )
        else
            commit, root = cd(dirname(file)) do
                readchomp(`git rev-parse HEAD`), readchomp(`git rev-parse --show-toplevel`)
            end
            if startswith(file, root)
                if isempty(repo)
                    repo = "https://github.com/$remote/tree/{commit}{path}#{line}"
                end

                _, path = split(file, root; limit = 2)
                repo = replace(repo, "{commit}", commit)
                repo = replace(repo, "{path}", path)
                repo = replace(repo, "{line}", line)

                Nullable{Compat.String}(repo)
            else
                Nullable{Compat.String}()
            end
        end
    end
else
    url(remote, repo, mod, file, line) = Nullable{Compat.String}()
end

function getremote(dir::AbstractString)
    remote =
        try
            cd(() -> readchomp(`git config --get remote.origin.url`), dir)
        catch err
            ""
        end
    match  = Utilities.nullmatch(isdefined(Base, :LibGit2) ?
        Base.LibGit2.GITHUB_REGEX : Pkg.Git.GITHUB_REGEX, remote)
    if isnull(match)
        travis = get(ENV, "TRAVIS_REPO_SLUG", "")
        isempty(travis) ? "" : travis
    else
        getmatch(match, 1)
    end
end

function inbase(m::Module)
    if m ≡ Base
        true
    else
        parent = module_parent(m)
        parent ≡ m ? false : inbase(parent)
    end
end

# Find line numbers.
# ------------------

linerange(doc) = linerange(doc.text, doc.data[:linenumber])

function linerange(text, from)
    lines = sum([isodd(n) ? newlines(s) : 0 for (n, s) in enumerate(text)])
    return lines > 0 ? (from:(from + lines + 1)) : (from:from)
end

function format_line(range::Range)
    local top = format_line(first(range))
    return length(range) <= 1 ? top : string(top, '-', format_line(last(range)))
end
format_line(line::Integer) = string('L', line)

newlines(s::AbstractString) = count(c -> c === '\n', s)
newlines(other) = 0

unwrap(f, x::Nullable) = isnull(x) ? nothing : f(get(x))

# Output redirection.
# -------------------

"""
Call a function and capture all `STDOUT` and `STDERR` output.

    withoutput(f) --> (result, success, backtrace, output)

where

  * `result` is the value returned from calling function `f`.
  * `success` signals whether `f` has thrown an error, in which case `result` stores the
    `Exception` that was raised.
  * `backtrace` a `Vector{Ptr{Void}}` produced by `catch_backtrace()` if an error is thrown.
  * `output` is the combined output of `STDOUT` and `STDERR` during execution of `f`.

"""
function withoutput(f)
    # Save the default output streams.
    local stdout = STDOUT
    local stderr = STDERR

    # Redirect both the `STDOUT` and `STDERR` streams to a single `Pipe` object.
    local pipe = Pipe()
    Base.link_pipe(pipe; julia_only_read = true, julia_only_write = true)
    redirect_stdout(pipe.in)
    redirect_stderr(pipe.in)

    # Bytes written to the `pipe` are captured in `output` and converted to a `String`.
    local output = UInt8[]

    # Run the function `f`, capturing all output that it might have generated.
    # Success signals whether the function `f` did or did not throw an exception.
    result, success, backtrace =
        try
            f(), true, Vector{Ptr{Void}}()
        catch err
            err, false, catch_backtrace()
        finally
            # Force at least a single write to `pipe`, otherwise `readavailable` blocks.
            println()
            # Restore the original output streams.
            redirect_stdout(stdout)
            redirect_stderr(stderr)
            # NOTE: `close` must always be called *after* `readavailable`.
            append!(output, readavailable(pipe))
            close(pipe)
        end
    return result, success, backtrace, chomp(Compat.String(output))
end

"""
    issubmodule(sub, mod)

Checks whether `sub` is a submodule of `mod`. A module is also considered to be
its own submodule.

E.g. `A.B.C` is a submodule of `A`, `A.B` and `A.B.C`, but it is not a submodule
of `D`, `A.D` nor `A.B.C.D`.
"""
function issubmodule(sub, mod)
    if (sub === Main) && (mod !== Main)
        return false
    end
    (sub === mod) || issubmodule(module_parent(sub), mod)
end

include("DOM.jl")
include("MDFlatten.jl")

end


text:
using CellSegmentation
using Test

@testset "CellSegmentation.jl" begin
    # Write your tests here.
end

ext: jl


text:
struct Point
    x::Int
    y::Int
    z::Int
end

const empty_points = Vector{Point}()

@enum Color begin
    red
    green
    yellow
    blue
    orange
    phosphorescent_green
    phosphorescent_yellow
    white
    clear_green
    colorless
end


text:
push!(LOAD_PATH,"../src/")
using DickeModel
using Documenter,Bibliography
using Distributed
using Plots
ENV["GKSwstype"]="nul"
addprocs(2,exeflags=`--project=docs/`)
@everywhere push!(LOAD_PATH,"../src/")
@everywhere using DickeModel
#bib = CitationBibliography(Base.Filesystem.joinpath(@__DIR__,"refs.bib"))
include("citations.jl")

bib= CitationBibliography(import_bibtex(Base.Filesystem.joinpath(@__DIR__,"refs.bib")))

#Documenter.Selectors.runner(::Type{Documenter.Expanders.ExampleBlocks}, x, page, doc)= (page.mapping[x] = Documents.RawNode(:hey, x.code))
#Documenter.Selectors.runner(::Type{Documenter.Expanders.SetupBlocks}, x, page, doc)= (page.mapping[x] = Documents.RawNode(:hey, x.code))

makedocs(bib,
         format = Documenter.HTML(
               prettyurls = get(ENV, "CI", nothing) == "true",
               assets = ["assets/customcss.css"]
           ),
         sitename = "DickeModel.jl",
         modules  = [DickeModel.ClassicalDicke, DickeModel.ClassicalSystems,
                    DickeModel.DickeBCE,  DickeModel.UPOs,  DickeModel.ClassicalLMG,  DickeModel.PhaseSpaces,
                    DickeModel.TWA.Weyl,
                    DickeModel.TWA, DickeModel.EnergyShellProjections],
         pages=[
                "DickeModel.jl" => "index.md",
                "Examples" => [
                    "ClassicalDicke" => "ClassicalDickeExamples.md",
                    "DickeBCE (Quantum Dicke)" => "DickeBCEExamples.md",
                    "EnergyShellProjections" => "EnergyShellProjectionsExamples.md",
                    "TWA" => "TWAExamples.md",
                    "UPOs" => "UPOsExamples.md",
                    "ClassicalLMG" => "ClassicalLMGExamples.md",
                ],
                "Documentation" => [
                    "ClassicalDicke" => "ClassicalDicke.md",
                    "DickeBCE" => "DickeBCE.md",
                    "EnergyShellProjections" => "EnergyShellProjections.md",
                    "TWA" => "TWA.md",
                    "UPOs" => "UPOs.md",
                    "ClassicalLMG" => "ClassicalLMG.md",
                    "ClassicalSystems" => "ClassicalSystems.md",
                    "PhaseSpaces" => "PhaseSpaces.md"
                ],
                "References" => "references.md"
               ])

deploydocs(;
    repo="github.com/saulpila/DickeModel.jl",
)


text:
include("laubloomis.jl")

using TaylorModels
using TaylorModels: validated_integ

# Equations of motion
# We write the function such that the operations are either unary or binary:
@taylorize function laubloomis!(t, x, dx)
    dx[1] = 1.4*x[3] - 0.9*x[1]
    dx[2] = 2.5*x[5] - 1.5*x[2]
    dx[3] = 0.6*x[7] - 0.8*(x[2]*x[3])
    dx[4] = 2 - 1.3*(x[3]*x[4])
    dx[5] = 0.7*x[1] - (x[4]*x[5])
    dx[6] = 0.3*x[1] - 3.1*x[6]
    dx[7] = 1.8*x[6] - 1.6*(x[2]*x[7])
    return dx
end

"""
    laubloomis_TMJets(; [t0], [T], [W], [abs_tol], [orderT], [orderQ],
                        [maxsteps], [property])

Build and run the Laub-Loomis model.

### Input

- `t0`       -- (optional, default: `0.0`) initial time
- `T`        -- (optional, default: `20.0`) time horizon
- `W`        -- (optional, default: `0.01`) width of the initial states
- `abs_tol`  -- (optional, default: `1e-10`) absolute tolerance used for time step
- `orderT`   -- (optional, default: `7`) order of the Taylor model in t
- `orderQ`   -- (optional, default: `2`) order of the Taylor model for Jet transport
                variables
- `maxsteps` -- (optional, default: `200`) use this maximum number of steps in
                the validated integration
- `property` -- (optional, default: `(t,x) -> x[4] < 4.5`) safe states property
"""
function laubloomis_TMJets(; t0=0.0, T=20.0, W=0.01, abs_tol=1e-10,
                             orderT=7, orderQ=2, maxsteps=1000,
                             property=(t,x) -> x[4] < 4.5)

    # center of initial conditions
    q0 = [1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45]

    # initial box (around `q0`) of the initial conditions
    δq0 = IntervalBox(-W..W, 7)

    # set variables
    set_variables("δ", numvars=length(q0), order=2orderQ)

    # TODO: wrap as a Reachability algorithm
    tTM, xTM = validated_integ(laubloomis!, q0, δq0, t0, T, orderQ, orderT,
                               abs_tol, maxsteps=maxsteps,
                               check_property=property)
    return tTM, xTM
end

#=
# return `true` if the specification is specification is satisfied and `false`
# otherwise
function check_property(xTM)
    satisfied = true
    for ind in eachindex(xTM[:])
        # check if solution crosses unsafe set (x4 ≥ 4.5)
        if (xTM[ind][4] >= 4.5)
            satisfied = false
            break
        end
    end
    return satisfied
end

using Plots

# plot solution and spec
function plot_solution(xTM)
    p = plot(xTM[:], color=:lightblue, legend=:false)
    plot!(p, [-2.5, 2.5], [2.75, 2.75], color=:red)
    xlims!(p, -2.5,2.5)
    ylims!(p, -3,3)
    return p
end

julia> data = [IntervalBox(IntervalArithmetic.Interval(tTM[i], tTM[i+1]), xTM[i][4]) for i in 1:length(tTM)-1];
julia> plot(data[:], legend=false, xlab="time", ylab="x4")

=#

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
# ---
# title: 685. Redundant Connection II
# id: problem685
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Tree, Depth-first Search, Union Find, Graph
# link: <https://leetcode.com/problems/redundant-connection-ii/description/>
# hidden: true
# ---
#
# In this problem, a rooted tree is a **directed** graph such that, there is
# exactly one node (the root) for which all other nodes are descendants of this
# node, plus every node has exactly one parent, except for the root node which
# has no parents.
#
# The given input is a directed graph that started as a rooted tree with N nodes
# (with distinct values 1, 2, ..., N), with one additional directed edge added.
# The added edge has two different vertices chosen from 1 to N, and was not an
# edge that already existed.
#
# The resulting graph is given as a 2D-array of `edges`. Each element of `edges`
# is a pair `[u, v]` that represents a **directed** edge connecting nodes `u`
# and `v`, where `u` is a parent of child `v`.
#
# Return an edge that can be removed so that the resulting graph is a rooted
# tree of N nodes. If there are multiple answers, return the answer that occurs
# last in the given 2D-array.
#
# **Example 1:**
#
#
#
#     Input: [[1,2], [1,3], [2,3]]
#     Output: [2,3]
#     Explanation: The given directed graph will be like this:
#       1
#      / \
#     v   v
#     2-->3
#
#
# **Example 2:**
#
#
#
#     Input: [[1,2], [2,3], [3,4], [4,1], [1,5]]
#     Output: [4,1]
#     Explanation: The given directed graph will be like this:
#     5 <- 1 -> 2
#          ^    |
#          |    v
#          4 <- 3
#
#
# **Note:**
#
# * The size of the input 2D-array will be between 3 and 1000.
# * Every integer represented in the 2D-array will be between 1 and N, where N is the size of the input array.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
const CLOSING_BRACKET = Dict{Char,Char}(
    '(' => ')',
    '[' => ']',
    '{' => '}',
)
const ISOPENER = Dict(k=>k in keys(CLOSING_BRACKET) for k in [keys(CLOSING_BRACKET)..., values(CLOSING_BRACKET)...])

function matching_brackets(s::String)
    desired_close::Array{Char} = []
    for c in s
        isopener = get(ISOPENER,c,nothing)
        if !isnothing(isopener)
            if isopener
                push!(desired_close, CLOSING_BRACKET[c])
            else
                (isempty(desired_close) || c != pop!(desired_close)) && return false
            end
        end
    end
    isempty(desired_close)
end


text:
"""
    TrelloCred(;api_key, api_token)

A credentials and identity object for Trello.
If the api_key or api_token are not provided, they will be read from
the `TRELLO_API_KEY` and `TRELLO_API_TOKEN` environment variables respectively.

Follow the instructions in the [Trello Documentation](https://developers.trello.com/docs/api-introduction#section--a-name-auth-authentication-and-authorization-a-) to aquire your credentials.
"""
struct TrelloCred
    api_key::String
    api_token::String
end

function TrelloCred(;api_key=ENV["TRELLO_API_KEY"], api_token=ENV["TRELLO_API_TOKEN"])
    return TrelloCred(api_key, api_token)
end

Base.show(io::IO, ::TrelloCred) = println(io, "TrelloCred(<secrets>)")


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module ComplexGates

using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.Gates
using QuantumCircuits.QCircuits.Gates: Parameter
using QuantumCircuits.QCircuits.Registers

import QuantumCircuits.QCircuits.Gates: getqubits, getqubitsids, appendparams!,
       appendRandParams!, setparameters!, getArgs, decompose, bindparameters!

export U4

const U4_params = 4*3 + 3

"Universary two qubits gate."
struct U4 <: QuantumGate
    qubit1::Qubit
    qubit2::Qubit
    params::Vector{Parameter}
end
U4(qubit1::Qubit, qubit2::Qubit) = U4(qubit1, qubit2, [ParameterT(rand() * 2π) for i in 1:U4_params])

Base.:(==)(g1::U4, g2::U4) = g1.qubit1 == g2.qubit1 && g1.qubit2 == g2.qubit2 && g1.params == g2.params
Base.hash(g::U4, h::UInt) = hash((g.qubit1, g.qubit2, g.params), h)

getArgs(g::U4) = (getid(g.qubit1), getid(g.qubit2), params)

"Return the qubits on which operate the gate"
getqubits(gate::U4) = (gate.qubit1, gate.qubit2)

"Return the qubits ids on which operate the gate"
getqubitsids(gate::U4) = (getid(gate.qubit1), getid(gate.qubit2))

Base.show(io::IO, gate::U4) = print(io, "U4($(gate.qubit1), $(gate.qubit2))")

function appendparams!(param, gate::U4)
    for p in gate.params
        @assert !isnothing(p) "Unable to append param without set it first."
        appendparams!(param, p)
    end
end

# function bindparameters!(gate::U4)
#     for p in gate.params
#         @assert p != nothing "Unable to append param without set it first."
#         appendparams!(param, p)
#     end
#     gate.θ = getvalue(gate.θ)
#     gate.ϕ = getvalue(gate.ϕ)
#     gate.λ = getvalue(gate.λ)
# end
#

appendRandParams!(param, gate::U4) = append!(param, rand(U4_params) * 2π)

Base.length(::U4) = U4_params

function setparameters!(gate::U4, params)
    @assert length(params) == U4_params "U4 has $(U4_params) parameters but $(length(param)) was provided."
    for i in 1:U4_params
        gate.params[i] = ParameterT(params[i])
    end
end

function decompose(gate::U4)
    return [
        U3(gate.qubit1, gate.params[1], gate.params[2], gate.params[3]),
        U3(gate.qubit2, gate.params[4], gate.params[5], gate.params[6]),
        CX(gate.qubit2, gate.qubit1),
        Rz(gate.qubit1, gate.params[7]),
        Ry(gate.qubit2, gate.params[8]),
        CX(gate.qubit1, gate.qubit2),
        Ry(gate.qubit2, gate.params[9]),
        CX(gate.qubit2, gate.qubit1),
        U3(gate.qubit1, gate.params[10], gate.params[11], gate.params[12]),
        U3(gate.qubit2, gate.params[13], gate.params[14], gate.params[15])
    ]
end

end  # module ComplexGates


text:
using Documenter
using Erdos

# Copy the list of features from the README to the
# Getting Started page in docs, to avoid duplication.
readme = readlines("README.md")
indexbase = readlines("docs/src/indexbase.md")
idx_features = findfirst(==("## Features"), readme) + 1
open("docs/src/index.md", "w") do f
    for l in indexbase
        println(f, l)
    end
    for l in readme[idx_features:end]
        println(f, l)
    end
end

cp("datasets/README.md", "docs/src/datasets.md", force=true)

makedocs(
    modules     = [Erdos],
    format      = Documenter.HTML(),
    sitename    = "Erdos",
    doctest     = false,
    pages       = Any[
        "Getting Started" => "index.md",
        "Graph Types" => "graph_types.md",
        "Basic Interface" => "core.md",
        "Edge And Vertex Maps" => "maps.md",
        "Network Types" => "network_types.md",
        "Property Interface for Networks" => "network_interface.md",
        "Operators" => "operators.md",
        "Traversals" => "traversals.md",
        "Distance" => "distance.md",
        "Shortest Paths" => "shortestpaths.md",
        "Linear Algebra" => "linalg.md",
        "Reading / Writing Graphs" => "persistence.md",
        "Graph Datasets" => "datasets.md",
        "Deterministic Graphs" => "deterministicgraphs.md",
        "Random Graphs" => "randomgraphs.md",
        "Centrality Measures" => "centrality.md",
        "Community Structures" => "community.md",
        "Flow and Cut" => "flow.md",
        "Matching" => "matching.md",
        "Dismantling" => "dismantling.md",
        "Spanning Trees" => "spanningtrees.md",
    ]
)

deploydocs(
    repo   = "github.com/CarloLucibello/Erdos.jl.git",
    target = "build"
)


text:
abstract type ModelSelectionInput end

struct SimulatedModelSelectionInput{AF<:AbstractFloat, CUD<:ContinuousUnivariateDistribution} <: ModelSelectionInput
	M::Int
	n_particles::Int
	threshold_schedule::AbstractArray{AF,1}
	model_prior::DiscreteUnivariateDistribution
	parameter_priors::AbstractArray{Array{CUD, 1},1}
	distance_simulation_input::AbstractArray{DistanceSimulationInput,1}
	max_iter::Int
end

struct EmulatedModelSelectionInput{AF<:AbstractFloat, CUD<:ContinuousUnivariateDistribution,
	    ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection, ET<:AbstractEmulatorTraining} <: ModelSelectionInput
	M::Int
	n_particles::Int
	threshold_schedule::AbstractArray{AF,1}
	model_prior::DiscreteUnivariateDistribution
	parameter_priors::AbstractArray{Array{CUD,1},1}
	emulator_training_input::AbstractArray{EmulatorTrainingInput{ET},1}
	emulator_retraining::ER
	emulated_particle_selection::EPS
	max_batch_size::Int
	max_iter::Int
end

mutable struct CandidateModelTracker{AF<:AbstractFloat}
	n_accepted::Int
	n_tries::Int
	population::AbstractArray{AF,2}
	distances::AbstractArray{AF,1}
	weight_values::AbstractArray{AF,1}
end

function CandidateModelTracker(n_params::Int)
	return CandidateModelTracker(0, 0, zeros(0, n_params), zeros(0), zeros(0))
end

abstract type ModelSelectionTracker end

mutable struct SimulatedModelSelectionTracker <: ModelSelectionTracker
	M::Int64
	n_particles::Int64
	threshold_schedule::AbstractArray{Float64,1}
	model_prior::DiscreteUnivariateDistribution
	smc_trackers::AbstractArray{SimulatedABCSMCTracker,1}
	max_iter::Int
end

mutable struct EmulatedModelSelectionTracker{CUD<:ContinuousUnivariateDistribution, ET,
        ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection} <: ModelSelectionTracker
	M::Int64
	n_particles::Int64
	threshold_schedule::AbstractArray{Float64,1}
	model_prior::DiscreteUnivariateDistribution
	smc_trackers::AbstractArray{EmulatedABCSMCTracker{CUD, ET, ER, EPS},1}
	max_batch_size::Int64
	max_iter::Int
end

"""
	ModelSelectionOutput

Contains results of a model selection computation, including which models are best supported by the observed data and the parameter poseteriors at each population for each model.

# Fields
- `M::Int64`: The number of models.
- `n_accepted::AbstractArray{AbstractArray{Int64,1},1}`: The number of parameters accepted by each model at each population. `n_accepted[i][j]` contains the number of acceptances for model `j` at population `i`.
- `threshold_schedule::AbstractArray{Float64,1}`: A set of maximum distances from the summarised model output to summarised observed data for a parameter vector to be included in the posterior.
- `smc_outputs::AbstractArray{ABCSMCOutput,1}`: A ['SimulatedABCSMCOutput']@(ref) or ['EmulatedABCSMCOutput']@(ref) for each model. Use to find details of the ABC results at each population.
- `completed_all_populations::Bool`: Indicates whether the algorithm completed all the populations successfully. A successful population is one where at least one model accepts at least one particle.
"""
struct ModelSelectionOutput
	M::Int64
	n_accepted::AbstractArray{AbstractArray{Int64,1},1}
	threshold_schedule::AbstractArray{Float64,1}
	smc_outputs::AbstractArray{ABCSMCOutput,1}
	completed_all_populations::Bool
end


text:
@testset "338.counting-bits.jl" begin
    @test count_bits(40) ==
          [0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3,
           4, 3, 4, 4, 5, 1, 2, 2, 3, 2, 3, 3, 4, 2]
end


text:
"""
$(SIGNATURES)

Extends the tracking function of Tracking.jl to support phased antenna arrays.
"""
function init_tracking(
        ::Val{N},
        system,
        inits,
        sample_freq,
        interm_freq,
        sat_prn;
        pll_bandwidth = 18Hz,
        dll_bandwidth = 1Hz,
        min_integration_time = 0.5ms,
        max_integration_time = 1ms,
        carrier_loop_func = Tracking.init_3rd_order_bilinear_loop_filter,
        code_loop_func = Tracking.init_2nd_order_bilinear_loop_filter,
        cn0_update_time = 20ms
    ) where N
    code_shift = Tracking.CodeShift{3}(system, sample_freq, 0.5) # 3: Early, Prompt, Late
    dopplers = Tracking.Dopplers(inits)
    phases = Tracking.Phases(inits)
    carrier_loop = carrier_loop_func(pll_bandwidth)
    code_loop = code_loop_func(dll_bandwidth)
    correlator_outputs = init_correlator_outputs(Val(N), code_shift)
    data_bits = Tracking.DataBits(system)
    last_valid_correlator_outputs = copy(correlator_outputs)
    cn0_state = Tracking.CN0State(cn0_update_time)
    Tracking.req_signal_and_track(correlator_outputs, last_valid_correlator_outputs, system, sample_freq, interm_freq, inits, dopplers, phases, code_shift, carrier_loop, code_loop, sat_prn, min_integration_time, max_integration_time, 0, data_bits, cn0_state)
end

function init_correlator_outputs(::Val{A}, code_shift::Tracking.CodeShift{N}) where {N,A}
    zeros(SMatrix{N, A, ComplexF64})
end

Base.@propagate_inbounds function dump!(output, signal::Matrix, output_idx, sample, code_carrier)
    for ant_idx = 1:size(output, 1)
        @fastmath output[output_idx,ant_idx] += signal[sample,ant_idx] * code_carrier
    end
end


text:
struct AngularDistance{M,R1,S} <: TaskMap{M,R1,S}
    reference_angle::SVector{1,S} # M = S1
end
AngularDistance{M,R1}(c::SVector{n,S}) where {M,R1,S,n} = AngularDistance{M,R1,S}(c)

struct AngularDistanceT{M,R1,S} <: TaskMapT{M,R1,S}
    base_map::fulltype(AngularDistance{M,R1,S})
end
AngularDistanceT{M,R1}(c::SVector{n,S}) where {M,R1,S,n} =
    AngularDistanceT{M,R1,S}(AngularDistance{M,R1,S}(c))

# S1
domain_coord_rep(::AngularDistance{S1,R1}) = ChartRep()

function task_map_chart_emb(::ChartRep, ::EmbRep, xm, task_map::AngularDistance{S1,R1},
        ::Chart{<:AngleChart,S1})
    Δθ = xm - task_map.reference_angle
    wrapToPi.(Δθ)
end


text:
@testset "Simulator" begin

    T = Float64
    # vis = Visualizer()
    # open(vis)
    # include(joinpath(module_dir(), "src", "dynamics", "quadruped", "visuals.jl"))

    s = get_simulation("quadruped", "flat_2D_lc", "flat")

    ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
        load_type = :split_traj_alt))

    H = 10
    H = ref_traj.H
    h = ref_traj.h
    q0_sim = SVector{s.model.dim.q}(deepcopy(ref_traj.q[1]))
    q1_sim = SVector{s.model.dim.q}(deepcopy(ref_traj.q[2]))

    p = open_loop_policy(ref_traj.u; N_sample = 1)

    sim_b = simulator(s, q0_sim, q1_sim, h, H, p=p, ip_type=:interior_point,
        sim_opts = SimulatorOptions(warmstart=true))

    @time simulate!(sim_b)
    qerr, uerr, γerr, berr = tracking_error(ref_traj, sim_b.traj, 1, idx_shift = 1)
    @test qerr < 0.0555 * 1.5
    @test uerr < 0.0001 * 1.5
    @test γerr < 0.0491 * 1.5
    @test berr < 0.0871 * 1.5
    qerr > 0.0555 * 1.2 && @warn "simulator regression: q"
    uerr > 0.0001 * 1.2 && @warn "simulator regression: u"
    γerr > 0.0491 * 1.2 && @warn "simulator regression: γ"
    berr > 0.0871 * 1.2 && @warn "simulator regression: b"

    # @btime simulate!(deepcopy(sim_b))
    # plot_surface!(vis, s.env, ylims=[-0.5, 0.5])
    # anim = visualize_meshrobot!(vis, s.model, sim_b.traj, α=0.3, sample=1, name=:basic)
    # plot(hcat(sim_b.traj.q...)', legend=false)
end


text:
# Generic data structures to represent meshes read from some input file

# Elements
struct RawElements
    numbers::Vector{Int}
    topology::Matrix{Int}
end
RawElements(;numbers, topology) = RawElements(numbers, topology)

getnumbers(elements::RawElements) = elements.numbers
gettopology(elements::RawElements) = elements.topology

getnumelements(elements::RawElements) = length(getnumbers(elements))

# Nodes
struct RawNodes
    numbers::Vector{Int}
    coordinates::Matrix{Float64}
end
RawNodes(;numbers, coordinates) = RawNodes(numbers, coordinates)

getnumbers(nodes::RawNodes) = nodes.numbers
Ferrite.getcoordinates(nodes::RawNodes) = nodes.coordinates

getnumnodes(nodes::RawNodes) = length(getnumbers(nodes))
getdim(nodes::RawNodes) = size(getcoordinates(nodes),1)
getcoordinate(nodes::RawNodes, number) = getcoordinates(nodes)[:,number]

# Complete mesh
struct RawMesh
    elementsdicts::Dict{String, RawElements}
    nodes::RawNodes
    elementsets::Dict{String, Vector{Int}}
    nodesets::Dict{String, Vector{Int}}
end
RawMesh(;elements,nodes,elementsets,nodesets) = RawMesh(elements,nodes,elementsets,nodesets)

Ferrite.getnodes(mesh::RawMesh) = mesh.nodes
getelementsdicts(mesh::RawMesh) = mesh.elementsdicts
Ferrite.getnodesets(mesh::RawMesh) = mesh.nodesets
getelementsets(mesh::RawMesh) = mesh.elementsets

getdim(mesh::RawMesh) = getdim(getnodes(mesh))
getnumelements(mesh::RawMesh) = sum(getnumelements.(values(getelementsdicts(mesh))))

# Verify that mesh is ok according to what is supported
function checkmesh(mesh::RawMesh)
    # Check that no node numbers are missing
    nodes = getnodes(mesh)
    minimum(getnumbers(nodes)) == 1 || throw(InvalidFileContent("Node numbering must start with 1"))
    maximum(getnumbers(nodes)) == getnumnodes(nodes) || throw(InvalidFileContent("No node numbers may be skipped"))
    # Check that no element numbers are missing
    eldicts = getelementsdicts(mesh)
    minelnum = 2
    maxelnum = 0
    for (_, elems) in eldicts
        minelnum = min(minelnum, minimum(getnumbers(elems)))
        maxelnum = max(maxelnum, maximum(getnumbers(elems)))
    end
    minelnum == 1 || throw(InvalidFileContent("Element numbering must start with 1"))
    maxelnum == getnumelements(mesh) || throw(InvalidFileContent("No element numbers may be skipped"))
end


text:
__precompile__(true)
module JDF

using Blosc: Blosc
using DataFrames
using Missings: Missings
using BufferedStreams
#using RLEVectors
using WeakRefStrings, TimeZones

using StatsBase: rle, inverse_rle, countmap, sample

import Base: show, getindex, setindex!, eltype, names

using Base:size#, @v_str, >=, include, VERSION

using Serialization: serialize, deserialize
using CSV

import DataFrames: nrow, ncol

if VERSION >= v"1.3.0-rc1"
    import Base.Threads: @spawn
else
    macro spawn(_)
        println("JDF: parallel save/load do not work in < Julia 1.3")
    end
end

function __init__()
    if VERSION >= v"1.3.0-rc1"
        Blosc.set_num_threads(Threads.nthreads())
    else
        Blosc.set_num_threads(isdefined(Sys, :CPU_CORES) ? Sys.CPU_CORES : Sys.CPU_THREADS)
    end
end

export savejdf, loadjdf, ssavejdf, sloadjdf
export column_loader, column_loader!
export type_compress!, type_compress
export compress_then_write
export JDFFile, @jdf_str, jdfmetadata, metadata, nrow, ncol, size, names
export IsBitsType, eachcol, some_elm, getindex

"""
    jdf"path/to/JDFfile.jdf"

    JDFFile("path/to/JDFfile.jdf")

Define a JDF file, which you can apply `nrow`, `ncol`, `names` and `size`.

## Example
using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]
nrow(jdf"plsdel.jdf") # 3
ncol(jdf"plsdel.jdf") # 2
size(jdf"plsdel.jdf") # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

# clean up
rm("plsdel.jdf", force = true, recursive = true)
"""
struct JDFFile{T <: AbstractString}
    path::T
end

"""
    jdf"path/to/JDFfile.jdf"

    JDFFile("path/to/JDFfile.jdf")

Define a JDF file, which you can apply `nrow`, `ncol`, `names` and `size`.

## Example
using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]
nrow(jdf"plsdel.jdf") # 3
ncol(jdf"plsdel.jdf") # 2
size(jdf"plsdel.jdf") # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

# clean up
rm("plsdel.jdf", force = true, recursive = true)
"""
macro jdf_str(path)
    return :(JDFFile($path))
end

include("type-writer-loader/Bool.jl")
include("type-writer-loader/Char.jl")
include("type-writer-loader/categorical-arrays.jl")
include("type-writer-loader/Missing.jl")
include("type-writer-loader/Nothing.jl")
include("type-writer-loader/String.jl")
include("type-writer-loader/StringArray.jl")
include("type-writer-loader/Symbol.jl")
include("type-writer-loader/ZonedDateTime.jl")

include("column_loader.jl")
include("compress_then_write.jl")

include("loadjdf.jl")
include("savejdf.jl")
include("type_compress.jl")

include("metadata.jl")
include("eachcol.jl")
include("dataframe-syntax.jl")

# Blosc.set_num_threads(6)

end # module


text:
__precompile__()
module TestUtilities
using OpenEphysLoader
using Base: ImmutableDict

using Test

export filecontext, test_fields

function filecontext(file_reader::Function, args...; kwargs...)
    wrapped_reader(::AbstractString, io::IOStream) = file_reader(io)
    pathiocontext(wrapped_reader, args...; kwargs...)
end

function pathiocontext(
    file_reader::Function,
    file_writer::Function,
    args...;
    kwargs...
)
    path, io = mktemp()
    try
        file_writer(io, args...; kwargs...)
        close(io)
        io = open(path)
        file_reader(path, io)
    finally
        close(io)
        rm(path)
    end
end

@generated function test_fields(a::Any, args...; kwargs...)
    fields = fieldnames(a)
    nfield = length(fields)
    testarr = Vector{Expr}(undef, nfield)
    for (i, field) in enumerate(fields)
        check_ex = :(get(kwarg_dict, $(QuoteNode(field)), true))
        fldtype = fieldtype(a, field)
        moddefined = parentmodule(fldtype) == OpenEphysLoader
        if moddefined
            test_ex = :(test_fields(a.$(field), args[j]...))
        else
            test_ex = :(@test a.$(field) == args[j])
        end
        testarr[i] = quote
            if $(check_ex)
                $(test_ex)
                j += 1
            end
        end
    end
    return quote
        @assert $(nfield) == length(args) + length(kwargs) "Args and kwargs must match number of fields"
        j = 1
        kwarg_dict = Dict(kwargs)
        $(testarr...)
    end
end

end


text:
abstract type AbstractDiagonalOperator{T <: Number} end
abstract type AbstractGeometricOperator{T <: Number} end
Base.length(D::AbstractDiagonalOperator) = length(D.u_cache)
Base.eltype(::AbstractDiagonalOperator{T}) where {T} = T
Base.size(G::AbstractGeometricOperator, inds...) = size(G.u_cache, inds...)

struct DiagonalOperator{T} <: AbstractDiagonalOperator{T}
    ω_vec::Tuple
    u_cache::Vector{T}
end

function DiagonalOperator(funcs...)
    num_type = typeof(funcs[1](0.0))
    DiagonalOperator{num_type}(funcs, zeros(num_type, length(funcs)))
end

DiagonalOperator(funcs::Vector{T}) where {T} = DiagonalOperator(funcs...)

function (D::DiagonalOperator)(t)
    for i in eachindex(D.ω_vec)
        D.u_cache[i] = D.ω_vec[i](t)
    end
    Diagonal(D.u_cache)
end

struct DiagonalFunction{T} <: AbstractDiagonalOperator{T}
    func
    u_cache
end

function DiagonalFunction(func)
    tmp = func(0.0)
    DiagonalFunction{eltype(tmp)}(func, zero(tmp))
end

function (D::DiagonalFunction)(t)
    D.u_cache .= D.func(t)
    Diagonal(D.u_cache)
end

struct GeometricOperator{T} <: AbstractGeometricOperator{T}
    funcs::Tuple
    u_cache::Matrix{T}
end

function GeometricOperator(funcs...)
    dim = (sqrt(1 + 8 * length(funcs)) - 1) / 2
    if !isinteger(dim)
        throw(ArgumentError("Invalid input length."))
    else
        dim = Int(dim)
    end
    num_type = typeof(funcs[1](0.0))
    GeometricOperator{num_type}(funcs, zeros(num_type, dim + 1, dim + 1))
end

GeometricOperator(funcs::Vector{T}) where {T} = GeometricOperator(funcs...)

function (G::GeometricOperator)(t)
    len = size(G.u_cache, 1)
    for j = 1:len
        for i = (j + 1):len
            G.u_cache[i, j] = G.funcs[i - 1 + (j - 1) * len](t)
        end
    end
    Hermitian(G.u_cache, :L)
end

struct ZeroGeometricOperator{T} <: AbstractGeometricOperator{T}
    size::Int
end

(G::ZeroGeometricOperator{T})(s) where {T} = Diagonal(zeros(T, G.size))

"""
$(TYPEDEF)

Defines a time dependent Hamiltonian in adiabatic frame.

# Fields

$(FIELDS)
"""
struct AdiabaticFrameHamiltonian{T} <: AbstractDenseHamiltonian{T}
    "Geometric part"
    geometric::AbstractGeometricOperator
    "Adiabatic part"
    diagonal::AbstractDiagonalOperator
    "Size of the Hamiltonian"
    size::Any
end

function AdiabaticFrameHamiltonian(D::AbstractDiagonalOperator{T}, G::AbstractGeometricOperator) where {T}
    l = length(D)
    AdiabaticFrameHamiltonian{complex(T)}(G, D, (l, l))
end

"""
    function AdiabaticFrameHamiltonian(ωfuns, geofuns)

Constructor of adiabatic frame Hamiltonian. `ωfuns` is a 1-D array of functions which specify the eigen energies (in `GHz`) of the Hamiltonian. `geofuns` is a 1-D array of functions which specifies the geometric phases of the Hamiltonian. `geofuns` can be thought as a flattened lower triangular matrix (without diagonal elements) in column-major order.
"""
function AdiabaticFrameHamiltonian(ωfuns, geofuns)
    if isa(ωfuns, AbstractVector)
        D = DiagonalOperator(ωfuns)
    elseif isa(ωfuns, Function)
        D = DiagonalFunction(ωfuns)
    end
    l = length(D)
    T = eltype(D)

    if isnothing(geofuns) || isempty(geofuns)
        G = ZeroGeometricOperator{T}(l)
    else
        G = GeometricOperator(geofuns)
        if size(G, 1) != l
            error("Diagonal and geometric operators do not match in size.")
        end
    end

    AdiabaticFrameHamiltonian(D, G)
end

function (H::AdiabaticFrameHamiltonian)(tf::Real, s::Real)
    ω = 2π * H.diagonal(s)
    off = H.geometric(s) / tf
    ω + off
end

get_cache(H::AdiabaticFrameHamiltonian{T}) where {T} = zeros(T, size(H))

"""
    function evaluate(H::AdiabaticFrameHamiltonian, s, tf)

Evaluate the adiabatic frame Hamiltonian at (unitless) time `s`, with total annealing time `tf` (in the unit of ``ns``). The final result is given in unit of ``GHz``.
"""
function evaluate(H::AdiabaticFrameHamiltonian, s, tf)
    ω = H.diagonal(s)
    off = H.geometric(s) / tf
    ω + off
end

function (h::AdiabaticFrameHamiltonian)(
    du::Matrix{T},
    u::Matrix{T},
    tf::Real,
    s::Real,
) where {T <: Number}
    ω = h.diagonal(s)
    du .= -2.0im * π * (ω * u - u * ω)
    G = h.geometric(s)
    du .+= -1.0im * (G * u - u * G) / tf
end


text:
type Model
    l1
    l2
    l3
    ls::Vector
end

function Model{T}(wordembeds::Matrix{T})
    l1 = Lookup(wordembeds)
    l2 = GatedLinear(T,300,300)
    l3 = GatedLinear(T,900,300)
    ls = [Linear(T,600,300), Linear(T,300,1)]
    Model(l1, l2, l3, ls)
end

type Model2
    f
    ls::Vector
end

function Model2{T}(wordembeds::Matrix{T})
    x = Var()
    y = Lookup(wordembeds)(x)
    y = window(y, (300,), pads=(100,), strides=(100,))
    y = GatedLinear(T,300,300)(y)
    y = window(y, (900,), pads=(300,), strides=(300,))
    y = GatedLinear(T,900,300)(y)
    y = window(y, (900,), pads=(300,), strides=(300,))
    y = GatedLinear(T,900,300)(y)
    y = window(y, (900,), pads=(300,), strides=(300,))
    y = GatedLinear(T,900,300)(y)
    y = window(y, (900,), pads=(300,), strides=(300,))
    y = GatedLinear(T,900,300)(y)
    y = window(y, (900,), pads=(300,), strides=(300,))
    y = GatedLinear(T,900,300)(y)
    #y = window(y, (900,), pads=(300,), strides=(300,))
    #y = GatedLinear(T,900,300)(y)
    #y = window(y, (900,), pads=(300,), strides=(300,))
    #y = GatedLinear(T,900,300)(y)

    y = pairwise(y, y)
    f = Graph(y, x)
    ls = [Linear(T,600,300), Linear(T,300,1)]
    Model2(f, ls)
end

function (m::Model)(x::Var, y=nothing)
    n = length(x.data)
    x = m.l1(x)
    x = window(x, (300,), pads=(100,), strides=(100,))
    x = m.l2(x)
    for i = 1:length(n)
        y == nothing || (x = dropout(x, 0.5))
        x = window(x, (900,), pads=(300,), strides=(300,))
        x = m.l3(x)
    end
    x = pairwise(x, x)
    n = size(x.data, 2)
    x = reshape(x, size(x.data,1), n*n)
    x = m.ls[1](x)
    x = tanh(x)
    x = m.ls[2](x)
    x = reshape(x, n, n)
    y == nothing ? argmax(x.data,1) : crossentropy(y, x)
end

function (m::Model2)(x::Var, y=nothing)
    x = m.f(x)
    n = size(x.data, 2)
    x = reshape(x, size(x.data,1), n*n)
    x = m.ls[1](x)
    x = tanh(x)
    x = m.ls[2](x)
    x = reshape(x, n, n)
    y == nothing ? argmax(x.data,1) : crossentropy(y, x)
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using Hiccup
using Test

@tags br, link

# hiccup div conflicts with main div, so use this as compromise
ediv = Hiccup.div

@test occursin("class=\"class1 class2\"", sprint(Hiccup.render, Node(:img, "#id.class1.class2", Dict(:src=>"http://www.com"))))

classMatching = ((".section-title", "section-title"),
                 (".test", "test"),
                 (".test1.test2", "test1 test2"),
                 (".-test", "-test"),
                 (".-test1.-test2", "-test1 -test2"),
                 ("#.test", "test"),
                 ("#.test1-hyphen", "test1-hyphen"),
                 ("#.test1-hyphen.test2", "test1-hyphen test2"),
                 ("#id.test1.test2", "test1 test2"),
                 ("#id.test1.-test2", "test1 -test2"),
                 ("#id.test1.-test-2", "test1 -test-2"),
                 ("#id.test-hyphen", "test-hyphen"),
                 ("#id.-test", "-test"),
                 ("#id.-test.test2", "-test test2"),
                 ("#id.-test.-test2", "-test -test2"))
for (in, expected) in classMatching
  @test occursin(expected, sprint(Hiccup.render, Hiccup.div(in, "contents")))
end

# tests for void tags
@test string(br()) == "
"
@test string(img(".image-test", [])) == ""
@test occursin(
  "/>",
  string(link(Dict(:rel => "stylesheet", :href => "test.css"))))
@test_throws ArgumentError img(strong(".test", "test"))

# tests for normal tags
@test string(ediv(ediv(ediv()))) == "

"

# test escapes
@test string(Node(:pre, "
fish & chips
")) ==
  "
<p>fish &amp; chips</p>
"

@test string(Node(:a, "link", href="http://example.com/test?a&b")) ==
  "link"


text:
defaulttol(a) = eps(real(float(one(eltype(a)))))^(1/2)

"""
 Take the L2 Tikhonov regularised inverse of a matrix `m`.

 The regularisation parameter is the larger of `delta` (the optional argument that defaults
 to zero) and square root of machine epsilon. The inverse is done using an SVD.
 """
 function reginv(a::AbstractMatrix, δ = defaulttol(a))
     U, S, V = svd(a)
     Sinv = inv.(hypot.(S, δ))
     return  V * Diagonal(Sinv) * U'
 end

 function posreginv(a::AbstractMatrix, δ = defaulttol(a))
     # assumes a is positive definite (or close to)
     D, V = eigen(Hermitian(a))
     Dinv = inv.(hypot.(D, δ))
     return  V * Diagonal(Dinv) * V'
 end

 # computes solution C of A * C + C * A = B for A hermitian
  function symm_sylvester_reg(A::AbstractMatrix, B::AbstractMatrix, δ = zero(eltype(A)))
      n = LinearAlgebra.checksquare(A)
      Λ, U = eigen(Hermitian(A))
      UdCU = U' * B * U
      for j = 1:n
          for i = 1:n
              UdCU[i,j] /= sqrt((Λ[i] + Λ[j])^2 + δ^2)
          end
      end
      C = U*UdCU*U'
      return C
  end


text:
module MPSKit
    using TensorKit,KrylovKit,Parameters, OptimKit, FastClosures
    using Base.Threads,FLoops,Transducers
    using Base.Iterators
    using RecipesBase

    using LinearAlgebra:diag,Diagonal;
    import LinearAlgebra

    #bells and whistles for mpses
    export InfiniteMPS,FiniteMPS,MPSComoving,PeriodicArray,MPSMultiline
    export transfer_left,transfer_right
    export leftorth,rightorth,leftorth!,rightorth!,poison!,uniform_leftorth,uniform_rightorth
    export r_LL,l_LL,r_RR,l_RR,r_RL,r_LR,l_RL,l_LR #should be properties

    #useful utility functions?
    export add_util_leg,max_Ds,virtualspace, recalculate!

    #hamiltonian things
    export Hamiltonian,Operator,Cache
    export MPOHamiltonian,InfiniteMPO,MPOMultiline
    export ac_prime,c_prime,environments,ac2_prime,expectation_value,effective_excitation_hamiltonian
    export leftenv,rightenv

    #algos
    export find_groundstate!, find_groundstate, Vumps, Dmrg, Dmrg2, GradDesc, Idmrg1, Idmrg2, GradientGrassmann
    export leading_boundary
    export excitations,FiniteExcited,QuasiparticleAnsatz
    export marek_gap, correlation_length
    export timestep!,timestep,Tdvp,Tdvp2,make_time_mpo,WI,WII
    export splitham,infinite_temperature, entanglement_spectrum, transfer_spectrum, variance
    export changebonds!,changebonds,VumpsSvdCut,OptimalExpand,SvdCut,UnionTrunc,RandExpand
    export entropy
    export dynamicaldmrg
    export fidelity_susceptibility
    export approximate!,approximate
    export periodic_boundary_conditions
    export exact_diagonalization

    @deprecate params(args...) environments(args...)

    #default settings
    module Defaults
        const eltype = ComplexF64
        const maxiter = 100
        const tolgauge = 1e-14
        const tol = 1e-12
        const verbose = true
        _finalize(iter,state,opp,envs) = (state,envs);
        import KrylovKit: GMRES
        const solver = GMRES(tol=1e-12, maxiter=100)
    end

    include("utility/periodicarray.jl")
    include("utility/utility.jl") #random utility functions
    export entanglementplot, transferplot
    include("utility/plotting.jl")

    #maybe we should introduce an abstract state type
    include("states/abstractmps.jl")
    include("states/transfer.jl") # mps transfer matrices
    include("states/infinitemps.jl")
    include("states/multiline.jl")
    include("states/finitemps.jl")
    include("states/comoving.jl")
    include("states/orthoview.jl")
    include("states/quasiparticle_state.jl")
    include("states/ortho.jl")

    abstract type Operator end
    abstract type Hamiltonian <: Operator end

    include("operators/mpohamiltonian/mpohamiltonian.jl") #the mpohamiltonian objects
    include("operators/umpo.jl")

    abstract type Cache end #cache "manages" environments

    include("environments/FinEnv.jl")
    include("environments/abstractinfenv.jl")
    include("environments/permpoinfenv.jl")
    include("environments/mpohaminfenv.jl")
    include("environments/overlapenv.jl")
    include("environments/qpenv.jl")
    include("environments/idmrgenv.jl")

    abstract type Algorithm end

    include("algorithms/derivatives.jl")
    include("algorithms/expval.jl")
    include("algorithms/toolbox.jl")
    include("algorithms/grassmann.jl")

    include("algorithms/changebonds/optimalexpand.jl")
    include("algorithms/changebonds/vumpssvd.jl")
    include("algorithms/changebonds/svdcut.jl")
    include("algorithms/changebonds/randexpand.jl")

    include("algorithms/timestep/tdvp.jl")
    include("algorithms/timestep/timeevmpo.jl")

    include("algorithms/groundstate/vumps.jl")
    include("algorithms/groundstate/idmrg.jl")
    include("algorithms/groundstate/dmrg.jl")
    include("algorithms/groundstate/gradient_grassmann.jl")

    include("algorithms/propagator/corvector.jl")

    include("algorithms/excitation/quasiparticleexcitation.jl")
    include("algorithms/excitation/dmrgexcitation.jl")

    include("algorithms/statmech/vumps.jl")
    include("algorithms/statmech/gradient_grassmann.jl")

    include("algorithms/fidelity_susceptibility.jl")

    include("algorithms/approximate/vomps.jl")
    include("algorithms/approximate/fvomps.jl")
    include("algorithms/approximate/idmrg.jl")

    include("algorithms/ED.jl")

    include("algorithms/unionalg.jl")

end


text:
function assemble_robj(data, program, bb, primitive, pre_fun, post_fun)
    transp = get(data, :transparency, Node(false))
    overdraw = get(data, :overdraw, Node(false))
    pre = if pre_fun != nothing
        _pre_fun = GLAbstraction.StandardPrerender(transp, overdraw)
        function ()
            _pre_fun()
            pre_fun()
        end
    else
        GLAbstraction.StandardPrerender(transp, overdraw)
    end
    robj = RenderObject(data, program, pre, nothing, bb, nothing)
    post = if haskey(data, :instances)
        GLAbstraction.StandardPostrenderInstanced(data[:instances], robj.vertexarray, primitive)
    else
        GLAbstraction.StandardPostrender(robj.vertexarray, primitive)
    end
    robj.postrenderfunction = if post_fun != nothing
        () -> begin
            post()
            post_fun()
        end
    else
        post
    end
    robj
end

function assemble_shader(data)
    shader = data[:shader]
    delete!(data, :shader)
    glp = get(data, :gl_primitive, GL_TRIANGLES)
    return assemble_robj(
        data, shader, FRect3D(), glp,
        get(data, :prerender, nothing),
        get(data, :postrender, nothing)
    )
end

points2f0(positions::Vector{T}, range::AbstractRange) where {T} = Point2f0[Point2f0(range[i], positions[i]) for i=1:length(range)]

"""
Converts index arrays to the OpenGL equivalent.
"""
to_index_buffer(x::GLBuffer) = x
to_index_buffer(x::TOrSignal{Int}) = x
to_index_buffer(x::VecOrSignal{UnitRange{Int}}) = x
to_index_buffer(x::TOrSignal{UnitRange{Int}}) = x
"""
For integers, we transform it to 0 based indices
"""
to_index_buffer(x::AbstractVector{I}) where {I <: Integer} = indexbuffer(Cuint.(x .- 1))
function to_index_buffer(x::Node{<: AbstractVector{I}}) where I <: Integer
    indexbuffer(lift(x-> Cuint.(x .- 1), x))
end

"""
If already GLuint, we assume its 0 based (bad heuristic, should better be solved with some Index type)
"""
function to_index_buffer(x::VectorTypes{I}) where I <: Union{GLuint, LineFace{GLIndex}}
    indexbuffer(x)
end

to_index_buffer(x) = error(
    "Not a valid index type: $(typeof(x)).
    Please choose from Int, Vector{UnitRange{Int}}, Vector{Int} or a signal of either of them"
)


text:
# ## Non-zero inflated, continuous. e.g. LogNormal
# function expert_ll_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:NonZIContinuousExpert}

#     expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, yl) - logcdf.(d, yu))
#     expert_ll = (tu == 0.) ? -Inf : expert_ll

#     return expert_ll
# end

# ## Zero inflated, continuous. e.g. ZILogNormal
# function expert_ll_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:ZIContinuousExpert}

#     expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, yl) - logcdf.(d, yu))
#     expert_ll = (tu == 0.) ? -Inf : expert_ll

#     return expert_ll
# end

# ## Non-zero inflated, discrete. e.g. Poisson
# function expert_ll_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:NonZIDiscreteExpert}

#     expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, ceil.(yl) .- 1) - logcdf.(d, yu))

#     return expert_ll
# end

# ## Zero inflated, discrete. e.g. ZIPoisson
# function expert_ll_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:ZIDiscreteExpert}

#     expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, ceil.(yl) .- 1) - logcdf.(d, yu))

#     return expert_ll
# end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# Duck typing, v1 and v2 have to implement addition/substraction and scalar multiplication
function midpoint(v1, v2, dom::Domain, tol)
    # There has to be a midpoint
    @assert in(v2,dom) != in(v1,dom)
    if in(v2,dom)
        min=v1
        max=v2
    else
        min=v2
        max=v1
    end
    mid = NaN
    while sum(abs.(max-min)) > tol
        step = (max-min)/2
        mid = min+step
        in(mid,dom) ? max=mid : min=mid
    end
    mid
end

## Avoid ambiguity (because everything >=2D is tensor but 1D is not)
function boundary(g::ProductGrid{TG,T,N},dom::Domain{<:Number}) where {TG,T,N}
    println("This method being called means there is a 1D ProductGrid.")
end

function boundary(g::MaskedGrid{G,M},dom::Domain{<:Number}) where {G,M}
  # TODO merge supergrid?
    boundary(grid(g),dom)
end

"""
    boundary(g::AbstractGrid{TG,T},dom::Domain{N},tol=1e-12)

Create a grid on the boundary of the domain.
The grid determines the resolution of the boundary.
"""
function boundary(g::ProductGrid{TG,T},dom::EuclideanDomain{N},tol=1e-12) where {TG,N,T}
    # Initialize neighbours
    CartesianNeighbours = CartesianIndices(ntuple(k->-1:1,Val(N)))
    periodic = map(isperiodic, elements(g))
    midpoints = eltype(g)[]
    # for each element
    for i in eachindex(g)
        # for all neighbours
        for neighbourindex in CartesianNeighbours
            neighbour = ModCartesianIndicesBase.add_offset_mod((neighbourindex+i).I, ntuple(k->1, Val(N)), size(g), periodic)
            # check if any are on the other side of the boundary
            try
                if in(g[i],dom) != in(g[neighbour...],dom)
                    # add the midpoint to the grid
                    push!(midpoints, midpoint(g[i],g[neighbour...],dom,tol))
                end
            catch y
                isa(y,BoundsError) || rethrow(y)
            end
        end
    end
    ScatteredGrid(midpoints)
end

function boundary(g::AbstractGrid{T},dom::Domain{T},tol=1e-12) where {T <: Number}
    midpoints = T[]
    # for each element
    for i in eachindex(g)
        # check if any are on the other side of the boundary
        try
            if in(g[i],dom) != in(g[i+1],dom)
                # add the midpoint to the grid
                push!(midpoints, midpoint(g[i],g[i+1],dom,tol))
            end
        catch y
            isa(y,BoundsError) || rethrow(y)
        end
    end
    ScatteredGrid(midpoints)
end

function boundary(g::MaskedGrid{G,M},dom::EuclideanDomain{N}) where {G,M,N}
    boundary(supergrid(g),dom)
end

"""
A Masked grid that contains the elements of grid that are on the boundary of the domain
"""
function boundary_grid(grid::AbstractGrid, domain::Domain)
    mask = boundary_mask(grid, domain);
    MaskedGrid(grid,mask,domain);
end

isperiodics(g::ProductGrid) = map(isperiodic,elements(g))
isperiodics(g::AbstractGrid) = ntuple(k->false, Val(dimension(g)))

function boundary_mask(grid::AbstractGrid, domain::Domain, periodic=isperiodics(grid))
    S = size(grid)
    m = BitArray(undef,S)
    m[:] .= 0
    for i in eachindex(grid)
        if grid[i]∈domain
            t = true
            for bi in ModCartesianIndicesBase.nbindexlist(i, S, periodic)
                if !(grid[bi]∈domain)
                    t = false
                    break
                end
            end
            m[i] = !t
        end
    end
    m
end


text:
using Molly

using Base.Threads
using Statistics
using Test

@testset "Lennard-Jones energy conservation" begin
    temp = 1.0u"K"
    n_steps = 10_000
    box_size = SVector(50.0, 50.0, 50.0)u"nm"
    n_atoms = 2_000
    atom_mass = 40.0u"u"
    simulator = VelocityVerlet(dt=0.005u"ps")

    parallel_list = nthreads() > 1 ? (false, true) : (false,)
    lj_potentials = (
        LennardJones(cutoff=DistanceCutoff(        3.0u"nm"), nl_only=false),
        LennardJones(cutoff=ShiftedPotentialCutoff(3.0u"nm"), nl_only=false),
        LennardJones(cutoff=ShiftedForceCutoff(    3.0u"nm"), nl_only=false),
    )

    for parallel in parallel_list
        @testset "$lj_potential" for lj_potential in lj_potentials
            s = System(
                atoms=[Atom(charge=0.0, mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms],
                general_inters=(lj_potential,),
                coords=place_atoms(n_atoms, box_size, 0.6u"nm"),
                velocities=[velocity(atom_mass, temp) for i in 1:n_atoms],
                box_size=box_size,
                loggers=Dict("coords" => CoordinateLogger(100),
                                "energy" => TotalEnergyLogger(100)),
            )

            E0 = total_energy(s)
            @time simulate!(s, simulator, n_steps; parallel=parallel)

            ΔE = total_energy(s) - E0
            @test abs(ΔE) < 2e-2u"kJ * mol^-1"

            Es = s.loggers["energy"].energies
            maxΔE = maximum(abs.(Es .- E0))
            @test maxΔE < 2e-2u"kJ * mol^-1"

            @test abs(Es[end] - Es[1]) < 2e-2u"kJ * mol^-1"

            final_coords = last(s.loggers["coords"].coords)
            @test all(all(c .> 0.0u"nm") for c in final_coords)
            @test all(all(c .< box_size) for c in final_coords)
        end
    end
end


text:
# kuramoto_7-27-19.jl

# Starting off the day strong with some code from my other project,
# grok-the-dot.

function f(x)
  return 2x+1
end

function g(x)
  return x^2
end

function h(x)
  return 4x^2 + 4x + 11   # (2x + 1)^2 = 4x^2 + 4x + 11
end

G = collect(1:1_000_000)  # You can use underscores to
                          # separate out numbers for easier
                          # reading. It's lowkey one of my
                          # favorite features.

@time h.(G)
@time g.(f.(G))

@time map(λ -> λ^2, f.(G))
@time map(λ -> λ^2, map(λ -> 2λ+1, G))

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


text:
using PowerDynamics
using OrdinaryDiffEq

##
dir =  @__DIR__
include("$dir/real_comps.jl")

## parameters from Raphael

omega_0 = 0.
VBase = 320e3
WBase = 1e9
ratio = 10
kappa = atan(ratio)
lengthLine = 100.
res = 0.03 * lengthLine * WBase / VBase^2 # VBase^2/WBase
adm = 1 / (res + im * ratio * res)

## other components

L1 = LineReal(; from = 1, to = 2, Y = adm)
L2 = LineReal(; from = 2, to = 3, Y = adm)
S = SlackReal(; U = complex(1.)) # 1pu
T = ThirdOrderReal(H = 5.0, P = 1.0, D = 0.1, Ω = 50, E_f = 1., X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)
PQ = PQReal(; P = -0.5, Q=-0.1)

##

pg = PowerGrid([S, T, PQ], [L1, L2])

op = find_operationpoint(pg, sol_method=:rootfind)

##

rpg = rhs(pg)

ode = ODEProblem(rpg, op.vec, (0., 10.))

##

using ModelingToolkit

# https://mtk.sciml.ai/stable/tutorials/modelingtoolkitize/

mtk_sys = modelingtoolkitize(ode)

jac = generate_jacobian(mtk_sys)[1]

J = eval(jac)

J(op.vec, nothing, 0.) |> eigvals

##
using ForwardDiff: jacobian

constraint_idx = findall( 0 .== [rpg.mass_matrix[i, i] for i in 1:systemsize(pg)])

f!(dx, x) = rpg(dx, x, nothing, 0.)
f(x) = (dx = similar(x); rpg(dx, x, nothing, 0.); dx)

J = jacobian(f, op.vec)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
module Help

export help, apropos, @help

MODULE_DICT   = nothing
FUNCTION_DICT = nothing

function clear_cache()
    global MODULE_DICT   = nothing
    global FUNCTION_DICT = nothing
end

function decor_help_desc(func::String, mfunc::String, desc::String)
    sd = convert(Array{ByteString,1}, split(desc, '\n'))
    for i = 1:length(sd)
        if beginswith(sd[i], func)
            sd[i] = mfunc * sd[i][length(func)+1:end]
        else
            break
        end
    end
    return join(sd, '\n')
end

function helpdb_filename()
    root = "$JULIA_HOME/../share/julia"
    file = "helpdb.jl"
    for loc in [Base.locale()]
        fn = joinpath(root, loc, file)
        if isfile(fn)
            return fn
        end
    end
    joinpath(root, file)
end

function init_help()
    global MODULE_DICT, FUNCTION_DICT
    if FUNCTION_DICT == nothing
        info("Loading help data...")
        helpdb = evalfile(helpdb_filename())
        MODULE_DICT = Dict()
        FUNCTION_DICT = Dict()
        for (mod,func,desc) in helpdb
            if !isempty(mod)
                mfunc = mod * "." * func
                desc = decor_help_desc(func, mfunc, desc)
            else
                mfunc = func
            end
            if !haskey(FUNCTION_DICT, mfunc)
                FUNCTION_DICT[mfunc] = {}
            end
            push!(FUNCTION_DICT[mfunc], desc)
            if !haskey(MODULE_DICT, func)
                MODULE_DICT[func] = {}
            end
            if !in(mod, MODULE_DICT[func])
                push!(MODULE_DICT[func], mod)
            end
        end
    end
end

function help(io::IO)
    init_help()
    print(io, """

     Welcome to Julia. The full manual is available at

        http://docs.julialang.org

     To get help, try help(function), help("@macro"), or help("variable").
     To search all help text, try apropos("string").
    """)
end

function print_help_entries(io::IO, entries)
    first = true
    for desc in entries
        if !first
            println(io)
        end
        println(io, strip(desc))
        first = false
    end
end

func_expr_from_symbols(s::Vector{Symbol}) = length(s) == 1 ? s[1] : Expr(:., func_expr_from_symbols(s[1:end-1]), Expr(:quote, s[end]))

function help(io::IO, fname::String, obj=0)
    init_help()
    found = false
    if haskey(FUNCTION_DICT, fname)
        print_help_entries(FUNCTION_DICT[fname])
        found = true
    elseif haskey(MODULE_DICT, fname)
        allmods = MODULE_DICT[fname]
        alldesc = {}
        for mod in allmods
            mfname = isempty(mod) ? fname : mod * "." * fname
            if isgeneric(obj)
                mf = eval(func_expr_from_symbols(map(symbol, split(mfname, "."))))
                if mf === obj
                    append!(alldesc, FUNCTION_DICT[mfname])
                    found = true
                end
            else
                append!(alldesc, FUNCTION_DICT[mfname])
                found = true
            end
        end
        found && print_help_entries(io, alldesc)
    elseif haskey(FUNCTION_DICT, "Base." * fname)
        print_help_entries(io, FUNCTION_DICT["Base." * fname])
        found = true
    end
    if !found
        if isa(obj, DataType)
            print(io, "DataType   : ")
            writemime(io, "text/plain", obj)
            println(io)
            println(io, "  supertype: ", super(obj))
            if obj.abstract
                st = subtypes(obj)
                if length(st) > 0
                    print(io, "  subtypes : ")
                    showcompact(io, st)
                    println(io)
                end
            end
            if length(obj.names) > 0
                println(io, "  fields   : ", obj.names)
            end
        elseif isgeneric(obj)
            writemime(io, "text/plain", obj); println()
        else
            println(io, "No help information found.")
        end
    end
end

apropos(s::String) = apropos(STDOUT, s)
function apropos(io::IO, txt::String)
    init_help()
    n = 0
    r = Regex("\\Q$txt", Base.PCRE.CASELESS)
    for (func, entries) in FUNCTION_DICT
        if ismatch(r, func) || any(e->ismatch(r,e), entries)
            for desc in entries
                nl = search(desc,'\n')
                if nl != 0
                    println(io, desc[1:(nl-1)])
                else
                    println(io, desc)
                end
            end
            n+=1
        end
    end
    if n == 0
        println(io, "No help information found.")
    end
end

function help(io::IO, f::Function)
    if is(f,help)
        return help(io)
    end
    help(io, string(f), f)
end

help(io::IO, t::DataType) = help(io, string(t.name), t)
help(io::IO, t::Module) = help(io, string(t))

function help(io::IO, x)
    show(io, x)
    t = typeof(x)
    if isa(t,DataType)
        println(io, " is of type")
        help(io, t)
    else
        println(io, " is of type $t")
    end
end

help(args...) = help(STDOUT, args...)

# check whether an expression is a qualified name, e.g. Base.FFTW.FORWARD
isname(n::Symbol) = true
isname(ex::Expr) = ((ex.head == :. && isname(ex.args[1]) && isname(ex.args[2]))
                    || (ex.head == :quote && isname(ex.args[1])))

macro help(ex)
    if !isa(ex, Expr) || isname(ex)
        return Expr(:call, :help, esc(ex))
    elseif ex.head == :macrocall && length(ex.args) == 1
        # e.g., "julia> @help @printf"
        return Expr(:call, :help, string(ex.args[1]))
    else
        return Expr(:macrocall, symbol("@which"), esc(ex))
    end
end

end # module

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
module BayesianSurvivalAnalysis

using Distributions: Weibull, Normal, MvNormal, pdf, cdf, logpdf, Truncated, Exponential
using LinearAlgebra: I, dot
using StatsBase: sample, mean
using SpecialFunctions: gamma

# Load files
#include("./src/utils.jl")
include("./src/simulate.jl")
include("./src/estimate.jl")

# Export
#export gibbs,
#metropolis,
#simulate,
#probability,
#latent,
#effective_sample_size

end


text:
# Conversion functions for date and time objects in the Python datetime
# module and the Julia Dates module.
import Compat.Dates

# Unfortunately, the Python C API (in Python/Include/datetime.h) is somewhat
# painful to call from Julia because it consists mainly of macros that
# translate to lookups in datastructures that are loaded at runtime by
# a PyDateTime_IMPORT macro.  So, we have to mirror those declarations here.
# Fortunately, they don't seem to have changed much since Python 2.7, with
# the biggest difference being the use of a 64-bit hash type.

struct PyDateTime_CAPI
    # type objects:
    DateType::PyPtr
    DateTimeType::PyPtr
    TimeType::PyPtr
    DeltaType::PyPtr
    TZInfoType::PyPtr

    # singletons:
    @static if pyversion >= v"3.7"
        TimeZone_UTC::PyPtr
    end

    # function pointers:
    Date_FromDate::Ptr{Cvoid}
    DateTime_FromDateAndTime::Ptr{Cvoid}
    Time_FromTime::Ptr{Cvoid}
    Delta_FromDelta::Ptr{Cvoid}
    @static if pyversion >= v"3.7"
        TimeZone_FromTimeZone::Ptr{Cvoid}
    end
    DateTime_FromTimestamp::Ptr{Cvoid}
    Date_FromTimestamp::Ptr{Cvoid}
end

struct PyDateTime_Delta{H} # H = Clong in Python 2, Py_hash_t in Python 3
    # PyObject_HEAD (for non-Py_TRACE_REFS build):
    ob_refcnt::Int
    ob_type::PyPtr
    hashcode::H
    days::Cint
    seconds::Cint
    microseconds::Cint
end

# the DateTime, Date, and Time objects are a struct with fields:
#     ob_refcnt::Int
#     ob_type::PyPtr
#     hashcode::Py_hash_t
#     hastzinfo::UInt8
#     unsigned char data[LEN]
#     tzinfo::PyPtr
# where LEN = 4 for Date, 6 for Time, and 10 for DateTime.  We
# will access this via raw Ptr{UInt8} loads, with the following offset
# for data:
const PyDate_HEAD = sizeof(Int)+sizeof(PyPtr)+sizeof(Py_hash_t)+1

# called from __init__
const DateType = Ref{PyPtr}(0)
const DateTimeType = Ref{PyPtr}(0)
const DeltaType = Ref{PyPtr}(0)
const Date_FromDate = Ref{Ptr{Cvoid}}(0)
const DateTime_FromDateAndTime = Ref{Ptr{Cvoid}}(0)
const Delta_FromDelta = Ref{Ptr{Cvoid}}(0)
function init_datetime()
    # emulate PyDateTime_IMPORT:
    PyDateTimeAPI = unsafe_load(@pycheckn ccall((@pysym :PyCapsule_Import),
                                     Ptr{PyDateTime_CAPI}, (Ptr{UInt8}, Cint),
                                     "datetime.datetime_CAPI", 0))
    DateType[] = PyDateTimeAPI.DateType
    DateTimeType[] = PyDateTimeAPI.DateTimeType
    DeltaType[] = PyDateTimeAPI.DeltaType
    Date_FromDate[] = PyDateTimeAPI.Date_FromDate
    DateTime_FromDateAndTime[] = PyDateTimeAPI.DateTime_FromDateAndTime
    Delta_FromDelta[] = PyDateTimeAPI.Delta_FromDelta
end

PyObject(d::Dates.Date) =
    PyObject(@pycheckn ccall(Date_FromDate[], PyPtr,
                             (Cint, Cint, Cint, PyPtr),
                             Dates.year(d), Dates.month(d), Dates.day(d),
                             DateType[]))

PyObject(d::Dates.DateTime) =
    PyObject(@pycheckn ccall(DateTime_FromDateAndTime[], PyPtr,
                             (Cint, Cint, Cint, Cint, Cint, Cint, Cint,
                              PyPtr, PyPtr),
                             Dates.year(d), Dates.month(d), Dates.day(d),
                             Dates.hour(d), Dates.minute(d), Dates.second(d),
                             Dates.millisecond(d) * 1000,
                             pynothing[], DateTimeType[]))

PyDelta_FromDSU(days, seconds, useconds) =
    PyObject(@pycheckn ccall(Delta_FromDelta[], PyPtr,
                             (Cint, Cint, Cint, Cint, PyPtr),
                             days, seconds, useconds,
                             1, DeltaType[]))

PyObject(p::Dates.Day) = PyDelta_FromDSU(Dates.value(p), 0, 0)

function PyObject(p::Dates.Second)
    # normalize to make Cint overflow less likely
    s = Dates.value(p)
    d = div(s, 86400)
    s -= d * 86400
    PyDelta_FromDSU(d, s, 0)
end

function PyObject(p::Dates.Millisecond)
    # normalize to make Cint overflow less likely
    ms = Dates.value(p)
    s = div(ms, 1000)
    ms -= s * 1000
    d = div(s, 86400)
    s -= d * 86400
    PyDelta_FromDSU(d, s, ms * 1000)
end

PyDate_Check(o::PyObject) = pyisinstance(o, DateType[])
PyDateTime_Check(o::PyObject) = pyisinstance(o, DateTimeType[])
PyDelta_Check(o::PyObject) = pyisinstance(o, DeltaType[])

function pydate_query(o::PyObject)
    if PyDate_Check(o)
        return PyDateTime_Check(o) ? Dates.DateTime : Dates.Date
    elseif PyDelta_Check(o)
        return Dates.Millisecond
    else
        return Union{}
    end
end

function convert(::Type{Dates.DateTime}, o::PyObject)
    if PyDate_Check(o)
        dt = convert(Ptr{UInt8}, o.o) + PyDate_HEAD
        if PyDateTime_Check(o)
            Dates.DateTime((UInt(unsafe_load(dt,1))<<8)|unsafe_load(dt,2), # Y
                           unsafe_load(dt,3), unsafe_load(dt,4), # month, day
                           unsafe_load(dt,5), unsafe_load(dt,6), # hour, minute
                           unsafe_load(dt,7), # second
                           div((UInt(unsafe_load(dt,8)) << 16) |
                               (UInt(unsafe_load(dt,9)) << 8) |
                               unsafe_load(dt,10), 1000)) # μs ÷ 1000
        else
            Dates.DateTime((UInt(unsafe_load(dt,1))<<8)|unsafe_load(dt,2), # Y
                           unsafe_load(dt,3), unsafe_load(dt,4)) # month, day
        end
    else
        throw(ArgumentError("unknown DateTime type $o"))
    end
end

function convert(::Type{Dates.Date}, o::PyObject)
    if PyDate_Check(o)
        dt = convert(Ptr{UInt8}, o.o) + PyDate_HEAD
        Dates.Date((UInt(unsafe_load(dt,1)) << 8) | unsafe_load(dt,2), # Y
                   unsafe_load(dt,3), unsafe_load(dt,4)) # month, day
    else
        throw(ArgumentError("unknown Date type $o"))
    end
end

function delta_dsμ(o::PyObject)
    PyDelta_Check(o) || throw(ArgumentError("$o is not a timedelta instance"))
    p = unsafe_load(convert(Ptr{PyDateTime_Delta{Py_hash_t}}, o.o))
    return (p.days, p.seconds, p.microseconds)
end

# Should we throw an InexactError when converting a period
# that is not an exact multiple of the resulting unit?   For now,
# follow the lead of Dates and truncate; see Julia issue #9169.

function convert(::Type{Dates.Millisecond}, o::PyObject)
    (d,s,μs) = delta_dsμ(o)
    return Dates.Millisecond((86400d + s) * 1000 + div(μs, 1000))
end

function convert(::Type{Dates.Second}, o::PyObject)
    (d,s,μs) = delta_dsμ(o)
    return Dates.Second(86400d + s + div(μs, 1000000))
end

function convert(::Type{Dates.Day}, o::PyObject)
    (d,s,μs) = delta_dsμ(o)
    return Dates.Day(d + div(s + div(μs, 1000000), 86400))
end


text:
using Statistics

@testsuite "statistics" AT->begin
    @testset "std" begin
        @test compare(std, AT, rand(10))
        @test compare(std, AT, rand(10,1,2))
        @test compare(std, AT, rand(10,1,2); corrected=true)
        @test compare(std, AT, rand(10,1,2); dims=1)
    end

    @testset "var" begin
        @test compare(var, AT, rand(10))
        @test compare(var, AT, rand(10,1,2))
        @test compare(var, AT, rand(10,1,2); corrected=true)
        @test compare(var, AT, rand(10,1,2); dims=1)
        @test compare(var, AT, rand(10,1,2); dims=[1])
        @test compare(var, AT, rand(10,1,2); dims=(1,))
        @test compare(var, AT, rand(10,1,2); dims=[2,3])
        @test compare(var, AT, rand(10,1,2); dims=(2,3))
    end

    @testset "mean" begin
        @test compare(mean, AT, rand(2,2))
        @test compare(mean, AT, rand(2,2); dims=2)
        @test compare(mean, AT, rand(2,2,2); dims=[1,3])
        @test compare(x->mean(sin, x), AT, rand(2,2))
        @test compare(x->mean(sin, x; dims=2), AT, rand(2,2))
        @test compare(x->mean(sin, x; dims=[1,3]), AT, rand(2,2,2))
    end

    @testset "cov" begin
        s = 100
        @test compare(cov, AT, rand(s))
        @test compare(cov, AT, rand(Complex{Float64}, s))
        @test compare(cov, AT, rand(s, 2))
        @test compare(cov, AT, rand(Complex{Float64}, s, 2))
        @test compare(cov, AT, rand(s, 2); dims=2)
        @test compare(cov, AT, rand(Complex{Float64}, s, 2); dims=2)
        @test compare(cov, AT, rand(1:100, s))
    end

    @testset "cor" begin
        s = 100
        @test compare(cor, AT, rand(s))
        @test compare(cor, AT, rand(Complex{Float64}, s))
        @test compare(cor, AT, rand(s, 2))
        @test compare(cor, AT, rand(Complex{Float64}, s, 2))
        @test compare(cor, AT, rand(s, 2); dims=2)
        @test compare(cor, AT, rand(Complex{Float64}, s, 2); dims=2)
        @test compare(cor, AT, rand(1:100, s))
    end
end


text:
fp = dirname(@__FILE__)

@testset "Impulse responses to structural shocks identified by a DSGE-VECM" begin
    matdata = load(joinpath(fp, "../../../reference/dsgevecm_lambda_irfs.jld2"))
    ŷ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                               vec(matdata["DD"]), matdata["MM"],
                               matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                               Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                               matdata["cointvec"], Int(matdata["qahead"]),
                               vec(matdata["XXpred"]);
                               test_shocks = matdata["Shocks"])

    @test ŷ ≈ matdata["yypred"]

    Random.seed!(1793) # drawing shocks, so this shouldn't be the same
    ŷ1 = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                vec(matdata["DD"]), matdata["MM"],
                                matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                                matdata["cointvec"],
                                Int(matdata["qahead"]), vec(matdata["XXpred"]); draw_shocks = true)
    Random.seed!(1793) # test re-seeding works
    ŷ2 = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                vec(matdata["DD"]), matdata["MM"],
                                matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                                matdata["cointvec"],
                                Int(matdata["qahead"]), vec(matdata["XXpred"]); draw_shocks = true)
    Random.seed!(1793) # flipping shocks should yield the same thing
    @info "The following warning is expected."
    ŷ3 = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                vec(matdata["DD"]), matdata["MM"],
                                matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                                matdata["cointvec"],
                                Int(matdata["qahead"]), vec(matdata["XXpred"]); draw_shocks = true,
                                flip_shocks = true)

    @test !(ŷ1 ≈ matdata["yypred"])
    @test ŷ1 ≈ ŷ2
    @test ŷ1 ≈ ŷ3

    # Test drawing shocks
    Random.seed!(1793)
    ŷ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                               vec(matdata["DD"]), matdata["MM"],
                               matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                               Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                               matdata["cointvec"], 5, vec(matdata["XXpred"]);
                               draw_shocks = true)
    @test @test_matrix_approx_eq ŷ matdata["ypred_drawshocks"]

    # Now try passing in X̂ = zeros(...)
    test_shocks = matdata["Shocks"]
    test_shocks[:, 2:end] .= 0.
    ŷ1_fcast = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                      vec(matdata["DD"]), matdata["MM"],
                                      matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                      Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                                      matdata["cointvec"], 5;
                                      test_shocks = test_shocks)

    # Compute forecast with no shock
    n_obs = Int(matdata["nvar"])
    p = 4
    baseline_ŷ = zeros(n_obs, 5) # baseline ŷ
    baseline_X̂ = zeros(length(matdata["XXpred"]))
    for t = 1:5
        out = vec(baseline_X̂' * matdata["cct_sim"])
        baseline_ŷ[:, t] = out
        addcoint = baseline_X̂[1:3] + matdata["cointvec"] * out
        baseline_X̂ .= vcat(addcoint, 1., out, baseline_X̂[3 + 2:Int(matdata["k"]) - n_obs])
    end
    ŷ1 = ŷ1_fcast  - baseline_ŷ

    # Now use deviations keyword
    ŷ2 = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                               vec(matdata["DD"]), matdata["MM"],
                               matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                               Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                               matdata["cointvec"], 5;
                               deviations = true,
                               test_shocks = test_shocks)

    # Replicate using VECM IRFs
    ZZ = matdata["ZZ"][1:n_obs, :]
    DD = matdata["DD"][1:n_obs]
    MM = matdata["MM"][1:n_obs, :]
    a0_m = convert(Matrix{Float64},
                   dropdims(impulse_responses(matdata["TTT"], matdata["RRR"], ZZ, DD, MM,
                                              sqrt.(matdata["QQ"]), 1,
                                              accumulate = false); dims = 2)')
    rotation, r_a = qr(a0_m)
    rotation = sign.(diag(r_a)) .* rotation'
    Σ_chol = cholesky(matdata["sig_sim"]).L * rotation

    expY = zeros(n_obs, p + 5)
    expY[:, p + 1] = Σ_chol * test_shocks[:, 1]
    addcoint = zeros(3)
    for t = 2:5
        addcoint += matdata["cointvec"] * expY[:, p + t - 1]
        xT = vcat(addcoint, vec(expY[:, p + t - 1:-1:t])) # stacks addcoint w/ lag yₜ₋₁, ..., yₜ₋ₚ
        expY[:, p + t] = vec(xT' * matdata["cct_sim"][vcat(1:3, 5:32), :])
    end
    expY = expY[:, p + 1:end]

    @test @test_matrix_approx_eq expY ŷ1
    @test @test_matrix_approx_eq expY ŷ2

    # Check IRFs to individual exogenous structural shocks
    ŷ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                               vec(matdata["DD"]), matdata["MM"],
                               matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                               Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                               matdata["cointvec"], 5, vec(matdata["XXpred"]))
    @test @test_matrix_approx_eq ŷ matdata["ypred_indivshocks"]

    # Now use deviations keyword
    ŷ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                               vec(matdata["DD"]), matdata["MM"],
                               matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                               Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                               matdata["cointvec"], 5, deviations = true)
    neg_ŷ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                   vec(matdata["DD"]), matdata["MM"],
                                   matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                   Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                                   matdata["cointvec"], 5, deviations = true,
                                   flip_shocks = true)
    @test @test_matrix_approx_eq neg_ŷ -ŷ

    # Replicate using VECM irfs
        ZZ = matdata["ZZ"][1:n_obs, :]
        DD = matdata["DD"][1:n_obs]
        MM = matdata["MM"][1:n_obs, :]
        a0_m = convert(Matrix{Float64},
                       dropdims(impulse_responses(matdata["TTT"], matdata["RRR"], ZZ, DD, MM,
                                                  sqrt.(matdata["QQ"]), 1,
                                                  accumulate = false); dims = 2)')
        rotation, r_a = qr(a0_m)
        rotation = sign.(diag(r_a)) .* rotation'
        Σ_chol = cholesky(matdata["sig_sim"]).L * rotation

    for i in 1:size(matdata["RRR"], 2)
        shocks = zeros(n_obs)
        shocks[i] = -sqrt(matdata["QQ"][i, i])
        expY = zeros(n_obs, p + 5)
        expY[:, p + 1] = Σ_chol * shocks
        addcoint = zeros(3)
        for t = 2:5
            addcoint += matdata["cointvec"] * expY[:, p + t - 1]
            xT = vcat(addcoint, vec(expY[:, p + t - 1:-1:t])) # stacks addcoint w/ lag yₜ₋₁, ..., yₜ₋ₚ
            expY[:, p + t] = vec(xT' * matdata["cct_sim"][vcat(1:3, 5:32), :])
        end
        expY = expY[:, p + 1:end]

        @test @test_matrix_approx_eq expY ŷ[:, :, i]
    end
end

nothing


text:
module TensorBoardLogger

using ProtoBuf: readproto, writeproto, ProtoType
using CRC32c

using ImageCore: colorview, channelview
using ColorTypes: Colorant, Gray, GrayA, RGB, RGBA
using FileIO: FileIO, @format_str, Stream, save, load

# hasproperty is not defined before 1.2. This is Compat.hasproperty
if VERSION < v"1.2.0-DEV.272"
    using ProtoBuf: hasproperty
end

#TODO: Is there a more lightweight package for compmuting an histogram?
using StatsBase: Histogram, fit

using  Base.CoreLogging: CoreLogging, AbstractLogger, LogLevel, Info,
    handle_message, shouldlog, min_enabled_level, catch_exceptions

export TBLogger, reset!, set_step!, increment_step!, set_step_increment!
export log_histogram, log_value, log_vector, log_text, log_image, log_images,
       log_audio, log_audios, log_graph, log_embeddings, log_custom_scalar
export map_summaries

export ImageFormat, L, CL, LC, LN, NL, NCL, NLC, CLN, LCN, HW, WH, HWC, WHC,
       CHW, CWH,HWN, WHN, NHW, NWH, HWCN, WHCN, CHWN, CWHN, NHWC, NWHC, NCHW, NCWH

# Custom Scalar Plugin
export tb_multiline, tb_margin

# Wrapper types
export TBText, TBVector, TBHistogram, TBImage, TBImages, TBAudio, TBAudios

# Protobuffer definitions for tensorboard
include("protojl/tensorboard/tensorboard.jl")
using .tensorboard: Summary_Value, GraphDef, Summary, Event, SessionLog_SessionStatus, SessionLog
using .tensorboard: TensorShapeProto_Dim, TensorShapeProto, TextPluginData
using .tensorboard: TensorProto, SummaryMetadata, SummaryMetadata_PluginData, _DataType

include("PNG.jl")
using .PNGImage
include("ImageFormat.jl")

# Plugin name used to store julia-specific metadata.
const TB_PLUGIN_JLARRAY_NAME = "_jl_tbl_array_sz"

include("TBLogger.jl")
include("utils.jl")  # CRC Utils
include("event.jl")
include("Loggers/base.jl")
include("Loggers/LogValue.jl")
include("Loggers/LogText.jl")
include("Loggers/LogHistograms.jl")
include("Loggers/LogAudio.jl")
include("Loggers/LogEmbeddings.jl")

# Custom Scalar Plugin
include("Loggers/LogCustomScalar.jl")

include("logger_dispatch.jl")
include("logger_dispatch_overrides.jl")

# TODO move it: It depends on TBImage for the explicit interface...
# maybe we could split the explicit and the @log interfaces?
include("Loggers/LogImage.jl")

include("Deserialization/deserialization.jl")
include("Deserialization/lookahead.jl")
include("Deserialization/helpers.jl")
include("Deserialization/audio.jl")
include("Deserialization/histograms.jl")
include("Deserialization/images.jl")
include("Deserialization/tensor.jl")

using Requires
function __init__()
    @require LightGraphs="093fc24a-ae57-5d10-9952-331d41423f4d" begin
        using .LightGraphs
        include("Optional/LogGraph.jl")
    end
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include("Optional/Plots.jl")
    end
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
        include("Optional/PyPlot.jl")
    end
    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
        include("Optional/Tracker.jl")
    end
    @require ValueHistories="98cad3c8-aec3-5f06-8e41-884608649ab7" begin
        include("Optional/ValueHistories.jl")
    end
end

end # module