text:
# DXF Parser

export Parser, parse, debugparser, showstate
export DXFParseError, IncompleteDXFInput, UnexpectedDXFInput
export DXFDocument, HeaderVariable, DXFPoint, DXFEntity, DXFBlock
export DXFSection, sectiontype, sections, section
export tracing_parser

abstract type DXFParseError <: Base.Exception end

Base.@kwdef mutable struct Parser
    groups::Vector{DXFGroup}
    # index is the position in groups of the DXFGroup currently being
    # considered:
    index::Int = 1
    # pending is a stack of "opening" groups and reductions:
    pending::Vector = []
    # starts is a stack of indices into pending of pending "start"
    # DXFGroups for which we have not yet encountered the matching
    # "end" group.
    starts::Vector{Int} = Vector{Int}()

    # Debugging flags
    trace_parser_actions = false  # See @tracePA, start, reduce
    showstate_in_parse = false    # Call showstate in each iteration of parse loop
    trace_disambiguation = false  # trace disambiguation of parseraction methods
end

function Parser(groups::Vector{DXFGroup})
    Parser(groups=groups)
end

function showstate(parser::Parser)
    showstate(stdout, parser)
    nothing
end

function showstate(io::IO, parser::Parser)
    level = 1
    println("\n\nParser has $(length(parser.starts)) starts,")
    println("$(length(parser.pending)) pending tokens.")
    println("input index: $(parser.index)/$(length(parser.groups))")
    println("starts: $(parser.starts)")
    for i in 1:length(parser.pending)
        t = parser.pending[i]
        isstart = level <= length(parser.starts) && parser.starts[level] == i
        spaces = repeat("  ", level)
        prefix = if isstart "[" else " " end
        @printf("%8d%s%s %s\n", i, spaces, prefix, summary(t))
        if isstart
            level += 1
        end
    end
    nothing
end

"""
    debugparser(::Parser, rt=false) do ... end

If a DXFParseError occurs during the do block then call `showstate` to
display the `pending` buffer of the parser.  If rt is true then
rethrow the DXFParseError.

All other exceptions are rethrown.
"""
function debugparser(body, parser; io=stderr, rt=false)
    try
        body()
    catch e
        if e isa DXFParseError
            showerror(io, e)
            showstate(io, parser)
            if rt
                rethrow(e)
            end
        else
            rethrow(e)
        end
    end
end

struct UnexpectedDXFInput <: DXFParseError
    parser::Parser
    pending::DXFObject
    current::DXFObject
    message::String
end

function Base.showerror(io::IO, e::UnexpectedDXFInput)
    print(io, "UnexpectedDXFInput, $(e.current): $(e.message)")
end

function Parser(path::String)
    Parser(read_dxf_file(path))
end

function Parser(io::IO)
    parser(read_dxf_file(io))
end

"""
    current(::Parser)
Return the object at the top of the `pending` stack.
"""
function current(parser::Parser)
    last(parser.pending)
end

"""
    pendingindex(::Parser)
Return the index of the topmost start token.
"""
function pendingindex(parser::Parser)
    last(parser.starts)
end

"""
    pending(::Parser)
Return the topmost start token from the `pending` stack.
"""
function pending(parser::Parser)
    parser.pending[pendingindex(parser)]
end

"""
    shift(::Parser)
Shift the next input token from `parser.groups` to `parser.pending`.
"""
function shift(parser::Parser)
    push!(parser.pending, parser.groups[parser.index])
    parser.index += 1
    nothing
end

START_GROUP_TYPES = [
]

"""
    start(::Parser)
Declare that the topmost token of `pending` is a start token
by pushing its index onto the `start` stack.
"""
function start(parser::Parser)
    if parser.trace_parser_actions
        println("start $(last(parser.pending))")
    end
    push!(parser.starts, lastindex(parser.pending))
end

function reduce(parser::Parser, constructor, fromend=0)
    grabbed = grab(parser; fromend)
    reduction = constructor(grabbed)
    insert!(parser.pending,
            lastindex(parser.pending) + 1 - fromend,
            reduction)
    if parser.trace_parser_actions
        println("Reduced $(length(grabbed)) to $reduction")
    end
end

"""
    grab(::Parser)
Remove the elements pf the pending stack from the topmost start token
to the topmost token and return them.

`fromend` is a non-negative number counting from `lastindex`..
"""
function grab(parser::Parser; fromend=0)
    range = pop!(parser.starts):(lastindex(parser.pending) - fromend)
    elts = parser.pending[range]
    deleteat!(parser.pending, range)
    return elts
end

"""
    lookahead(::Parser)
return the next input token or nothing.
A second argument can be passed to specify an offset for additional
lookahead.
"""
function lookahead(parser::Parser, index=0)
    i = parser.index + index
    if i > length(parser.groups)
        nothing
    else
        parser.groups[i]
    end
end

# Supertype for DXFObjects that have a contents shot:
abstract type DXFContentsObject <: DXFObject end

function Base.iterate(o::DXFContentsObject)
    iterate(o.contents)
end

function Base.iterate(o::DXFContentsObject, state)
    iterate(o.contents, state)
end

function Base.summary(io::IO, o::DXFContentsObject)
    print(io, "$(typeof(o)) $(o.contents[1].value)) with $(length(o.contents)) elements")
end

function Base.getproperty(o::DXFContentsObject, p::Symbol)
    if hasfield(typeof(o), p)
        return getfield(o, p)
    end
    if p == :line
        # The line number of a DXFObject is the line number of its
        # first group:
        return o.contents[1].line
    end
    # To get the standard missing field error, since there doesn't
    # seem to be a defined Exception type for this yet:
    getfield(o, p)
end

function Base.propertynames(o::DXFContentsObject, private::Bool)
    [ :line,
      fieldnames(typeof(o, private))... ]
end

"""
DocumentStart is the start token at the beginning of the document.
"""
struct DocumentStart <: DXFObject
end

function Base.getproperty(ds::DocumentStart, p::Symbol)
    if p == :value
        return nothing
    end
    getfield(ds, p)
end

function Base.propertynames(ds::DocumentStart, private::Bool)
    [ :value,
      fieldnames(typeof(ds, private))... ]
end

function Base.summary(io::IO, ds::DocumentStart)
    print(io, "$(typeof(ds))")
end

function parse(groups::Vector{DXFGroup})
    parse(Parser(groups))
end

struct IncompleteDXFInput <: DXFParseError
    parser::Parser
end

function Base.showerror(io::IO, e::IncompleteDXFInput)
    remaining = length(e.parser.groups) - e.parser.index
    print(io, "parse ended prematurely: starts: $(length(e.parser.starts)), remaining input:$remaining")
end

function parse(parser::Parser)
    push!(parser.pending, DocumentStart())
    start(parser)
    while length(parser.starts) > 0 && parser.index <= length(parser.groups)
        shift(parser)
        if parser.showstate_in_parse
            showstate(parser)
        end
        # See if any reductions can be performed
        parseraction(parser)
    end
    if length(parser.starts) > 0 || parser.index <= length(parser.groups)
        remaining = length(parser.groups) - parser.index
        throw(IncompleteDXFInput(parser))
    end
    return parser
end

function tracePA_(frame)
    ip = 3:4   # interesting parameters to show
    psig(t) = fieldtypes(t)[ip]
    @assert frame.func == :parseraction
    println("tracePA: $(frame.func) " *
        "$(psig(frame.linfo.def.sig)) " *
        "$(psig(frame.linfo.specTypes))")
end

macro tracePA(parser_)
    quote
        if ($(esc(parser_))).trace_parser_actions
            tracePA_(first(stacktrace()))
        end
    end
end

"""
Each parseraction is a shift/reduce rule of our DXF parser.
"""
function parseraction end

function parseraction(parser::Parser)
    if length(parser.starts) <= 0
        return
    end
    p = pending(parser)
    c = current(parser)
    if parser.trace_parser_actions
        println("parseraction for $p $c")
    end
    try
        parseraction(parser, p, c)
    catch e
        if !(e isa Base.MethodError)
            rethrow(e)
        end
        if e.f != parseraction
            rethrow(e)
        end
        if length(e.args) == 1
            # Do not handle the one argument call to parseraction
            rethrow(e)
        end
        msg = IOBuffer()
        showerror(msg, e)
        msg = String(take!(msg))
        # I wish MethodError had subtypes to distinguish "no
        # applicable method" from "ambiguous methods".
        if !occursin("is ambiguous", msg)
            rethrow(e)
        end
        best = nothing
        for m in methods(parseraction)
            if length(m.sig.parameters) != 4
                continue
            end
            if !all(isa.(e.args, m.sig.parameters[2:length(m.sig.parameters)]))
                # m is not applicable
                continue
            end
            if best == nothing
                best = m
                continue
            end
            # proper subtype?
            better(spec1, spec2) = (spec1 <: spec2) && !(spec2 <: spec1)
            for pos in ((1, 3, 2).+1)   # m.sig[1] is the type of the function.
                mspec = m.sig.parameters[pos]
                bspec = best.sig.parameters[pos]
                if better(bspec, mspec)
                    break
                end
                if better(mspec, bspec)
                    if parser.trace_disambiguation
                        println("@@ $pos: $mspec of $m preferred to\n      $bspec of $best")
                    end
                    best = m
                    break
                end
            end
        end
        if best == nothing
            rethrow(e)
        end
        invoke(e.f,
               Tuple{best.sig.parameters[2:length(best.sig.parameters)]...},
               e.args...)
    end
end

"""
The default behavior is to shift the new token into pending.
Shifting of input DXFGroups happens in `parser`, whiuch shifts
all new tokens before calling `parseraction`.
"""
function parseraction(parser::Parser, pending::DXFGroup, current::DXFObject)
    @tracePA(parser)
    # Just shift, which happened in the caller.
end

# Document

"""
DocumentStart represents an entire DXF document.
"""
struct DXFDocument <: DXFContentsObject
    contents::Vector
end

function sections(doc::DXFDocument)
    sections = []
    for s in doc
        if !isa(s, DXFSection)
            continue
        end
        n = sectiontype(s)
        if n != nothing
            push!(sections, n)
        end
    end
    sections
end

function section(doc::DXFDocument, sect::String)::Union{Nothing, DXFSection}
    for e in doc.contents
        if e isa DXFSection
            if sectiontype(e) == sect
                return e
            end
        end
    end
    return nothing
end

function parseraction(parser::Parser, pending::DocumentStart, current::EntityType_EOF)
    @tracePA(parser)
    reduce(parser, DXFDocument)
    parseraction(parser)
end

# Sections and Entities

struct DXFSection <: DXFContentsObject
    contents::Vector
end

function Base.summary(io::IO, section::DXFSection)
    println(io, "$(typeof(section)) $(sectiontype(section)) with $(length(section.contents)) elements")
end

function sectiontype(sec::DXFSection)
    if sec.contents[2] isa Name
        sec.contents[2].value
    else
        nothing
    end
end

# A section can only follow the document start token or another sectioon:
function parseraction(parser::Parser, ::DocumentStart, ::DXFSection) @tracePA(parser) end
function parseraction(parser::Parser, ::DXFSection, ::DXFSection) @tracePA(parser) end

function parseraction(parser::Parser, pending::DocumentStart, current::EntityType_SECTION)
    @tracePA(parser)
    start(parser)
end

function parseraction(parser::Parser, pending::DXFGroup, current::EntityType_SECTION)
    @tracePA(parser)
    throw(UnexpectedDXFInput(parser, pending, current,
                             "EntityType is unexpected context"))
end

function parseraction(parser::Parser, pending::EntityType_SECTION, current::EntityType_ENDSEC)
    @tracePA(parser)
    reduce(parser, DXFSection)
    parseraction(parser)
end

# Non-section entities

# It seems that one entity ends where the next one begins.

struct DXFEntity <: DXFContentsObject
    contents::Vector
end

function parseraction(parser::Parser, pending::EntityType_SECTION, current::EntityType)
    @tracePA(parser)
    start(parser)
end

function parseraction(parser::Parser, pending::EntityType, current::EntityType)
    @tracePA(parser)
    # `current` starts a new entity.  First we must finish the previous one:
    reduce(parser, DXFEntity, 1)
    start(parser)
end

function parseraction(parser::Parser, pending::EntityType, current::EntityType_ENDSEC)
    @tracePA(parser)
    # ENDSEC closes the currently open entity and the section that
    # contains it.
    reduce(parser, DXFEntity, 1)
    parseraction(parser)
end

# Header Variables

struct HeaderVariable <: DXFObject
    name::HeaderVariableName
    value::DXFObject

    function HeaderVariable(contents)
        @assert length(contents) == 2 "length $contents ($(length(contents))) != 2"
        new(contents...)
    end
end

function Base.summary(io::IO, v::HeaderVariable)
    print(io, "$(typeof(v)) $(v.name) = $(v.value)")
end

### Maybe make this pending::EntityType_SECTION and test that section is HEADER.
function parseraction(parser::Parser, pending::DXFGroup, current::HeaderVariableName)
    @tracePA(parser)
    start(parser)
end

function parseraction(parser::Parser, pending::HeaderVariableName, current::DXFObject)
    @tracePA(parser)
    reduce(parser, HeaderVariable)
    parseraction(parser)
end

function parseraction(parser::Parser, pending::EntityType_SECTION, current::HeaderVariable)
    @tracePA(parser)
    @assert parser.pending[pendingindex(parser)] == pending
end

# Blocks

struct DXFBlock <: DXFContentsObject
    contents::Vector{DXFObject}
end

function parseraction(parser::Parser, pending::EntityType_BLOCK, current::DXFObject)
    @tracePA(parser)
end

function parseraction(parser::Parser, pending::EntityType_BLOCK, current::EntityType_ENDBLK)
    @tracePA(parser)
    reduce(parser, DXFBlock)
    parseraction(parser)
end

# Points

struct DXFPoint <: DXFObject
    pointX::PointX
    pointY::PointY
    pointZ::Union{Nothing, PointZ}

    DXFPoint(contents) = DXFPoint(contents...)

    function DXFPoint(x::PointX, y::PointY, z::PointZ)
        @assert groupcode(x) + 10 == groupcode(y) "$x $(groupcode(x)), $y $(groupcode(y))"
        @assert groupcode(x) + 20 == groupcode(z) "$x $(groupcode(x)), $z $(groupcode(z))"
        new(x, y, z)
    end

    function DXFPoint(x::PointX, y::PointY)
        @assert groupcode(x) + 10 == groupcode(y) "$x $(groupcode(x)), $y $(groupcode(y))"
        new(x, y, nothing)
    end
end

# We need some simple way to distinguish points of different types.
function groupcode(point::DXFPoint)
    groupcode(point.pointX)
end

function Base.summary(io::IO, p::DXFPoint)
    val(x) = if x isa DXFGroup x.value else nothing end
    print(io, "$(typeof(p)) $(val(p.pointX)), $(val(p.pointY)), $(val(p.pointZ))")
end

function parseraction(parser::Parser, pending::DXFGroup, current::PointX)
    @tracePA(parser)
    start(parser)
    if lookahead(parser) isa PointY
        shift(parser)
    end
    if lookahead(parser) isa PointZ
        shift(parser)
    end
    reduce(parser, DXFPoint)
    parseraction(parser)
end

function parseraction(parser::Parser, pending::PointX, current::DXFGroup)
    @tracePA(parser)
    throw(UnexpectedDXFInput(parser, pending, current,
                             "Expected PointY or PointZ after PointX, got $current"))
end

### For now at least DXFPoint can be contained by a DXFSection
function parseraction(parser::Parser, ::EntityType_SECTION, ::DXFPoint)
    @tracePA(parser)
end

struct ADGroup <: DXFContentsObject
    contents
end

function parseraction(parser::Parser, pending::DXFGroup, current::ADGroupStartEnd)
    @tracePA(parser)
    if length(current.value) > 1 && current.value[1] == '{'
        start(parser)
    else
        #=
        # We see ADGroupStartEnd groups that start with neither '{' nor '}'.
        # Until we understand more, just shift them.
        # We might eventually consaume the next entity to produce
        # an ADGroup with a single element, but until we know what that
        # content element can be we can't tell when to reduce.
        throw(UnexpectedDXFInput(parser, pending, current,
                           "non-opening ADGroupStartEnd in unexpected context"))
        =#
    end
end

function parseraction(parser::Parser, pending::ADGroupStartEnd, current::ADGroupStartEnd)
    @tracePA(parser)
    if length(current.value) == 1 && current.value[1] == '}'
        reduce(parser, ADGroup)
        parseraction(parser)
    else
        throw(UnexpectedDXFInput(parser, pending, current,
                           "non-closing ADGroupStartEnd in unexpected context"))
    end
end

function parseraction(parser::Parser, pending::DXFGroup, current::ADGroup)
    @tracePA(parser)
end


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 QuantumCircuits

include("QCircuits/QCircuits.jl")
include("QML/QML.jl")
include("Execute/Execute.jl")
include("Simulation/Simulation.jl")

using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.Registers
using QuantumCircuits.QCircuits.Circuit

export QuantumCircuit, QuantumDevice, QuantumGate, add!, QCircuit,
       QuantumRegister, ClassicalRegister, tomatrix, decompose, simplify, @gate, @circ

end # module

ext: jl


text:
@testset "$(rpad("Splitting Tableaus",80))" begin

    # instatiate all splitting tableaus
    @test typeof(TableauStrangA())       <: TableauSplitting{Float64}
    @test typeof(TableauStrangB())       <: TableauSplitting{Float64}

    @test typeof(TableauLieA())          <: TableauSplittingNS{Float64}
    @test typeof(TableauLieB())          <: TableauSplittingNS{Float64}
    @test typeof(TableauStrang())        <: TableauSplittingNS{Float64}
    @test typeof(TableauMcLachlan2())    <: TableauSplittingNS{Float64}
    @test typeof(TableauMcLachlan4())    <: TableauSplittingNS{Float64}

    @test typeof(TableauTripleJump())    <: TableauSplittingSS{Float64}
    @test typeof(TableauSuzukiFractal()) <: TableauSplittingSS{Float64}

end


text:
using CSV, CategoricalArrays, Plots; pyplot()

df = CSV.read("../data/companyData.csv")
mktCap = reshape(df.MarketCap, 5, 3)
years  = levels(df.Year)

areaplot(years, mktCap,
	c=[:blue :red :green], labels=["A" "B" "C"],
	xlims=(minimum(years),maximum(years)), ylims=(0,6.5),
	legend=:topleft, xlabel="Years", ylabel="MarketCap")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# ---
# title: 216. Combination Sum III
# id: problem216
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Backtracking
# link: <https://leetcode.com/problems/combination-sum-iii/description/>
# hidden: true
# ---
#
# Find all valid combinations of `k` numbers that sum up to `n` such that the
# following conditions are true:
#
#   * Only numbers `1` through `9` are used.
#   * Each number is used **at most once**.
#
# Return _a list of all possible valid combinations_. The list must not contain
# the same combination twice, and the combinations may be returned in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: k = 3, n = 7
#     Output: [[1,2,4]]
#     Explanation:
#     1 + 2 + 4 = 7
#     There are no other valid combinations.
#
# **Example 2:**
#
#
#
#     Input: k = 3, n = 9
#     Output: [[1,2,6],[1,3,5],[2,3,4]]
#     Explanation:
#     1 + 2 + 6 = 9
#     1 + 3 + 5 = 9
#     2 + 3 + 4 = 9
#     There are no other valid combinations.
#
#
# **Example 3:**
#
#
#
#     Input: k = 4, n = 1
#     Output: []
#     Explanation: There are no valid combinations. [1,2,1] is not valid because 1 is used twice.
#
#
# **Example 4:**
#
#
#
#     Input: k = 3, n = 2
#     Output: []
#     Explanation: There are no valid combinations.
#
#
# **Example 5:**
#
#
#
#     Input: k = 9, n = 45
#     Output: [[1,2,3,4,5,6,7,8,9]]
#     Explanation:
#     1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 = 45
#     There are no other valid combinations.
#
#
#
#
# **Constraints:**
#
#   * `2 <= k <= 9`
#   * `1 <= n <= 60`
#
#
## @lc code=start
using LeetCode

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


text:
"""
    wrap_to_load!(updates, node, nodehash;
                  cachedir=DEFAULT_CACHE_DIR,
                  compression=DEFAULT_COMPRESSION)

Generates a new dispatch node that corresponds to `node::DispatchNode`
and which loads a file from the `cachedir` cache directory whose name and extension
depend on `nodehash` and `compression` and contents are the output of `node`.
The generated node is added to `updates` which maps `node` to the generated node.
"""
function wrap_to_load!(updates::Dict{DispatchNode, DispatchNode},
                       node::DispatchNode,
                       nodehash::String;
                       cachedir::String=DEFAULT_CACHE_DIR,
                       compression::String=DEFAULT_COMPRESSION)
    # Defing loading wrapper function
    """
    Simple load wrapper.
    """
    function loading_wrapper()
        _cachedir = abspath(joinpath(expanduser(cachedir), DEFAULT_HASHCACHE_DIR))
        if compression != "none"
            extension = ".$compression"
            operation = "LOAD-UNCOMPRESS"
        else
            extension = ".bin"
            operation = "LOAD"
        end
        filepath = joinpath(_cachedir, nodehash * extension)
        decompressor = get_compressor(compression, "decompress")
        label = _labelize(node, nodehash)
        @debug "[$nodehash][$(label)] $operation (compression=$compression)"
        if isfile(filepath)
            result = open(decompressor, filepath, "r") do fid
                        deserialize(fid)
                    end
        else
            throw(ErrorException("Cache file $filepath is missing."))
        end
        return result
    end

    # Add wrapped node to updates (no arguments to update :)
    newnode = Op(loading_wrapper)
    newnode.label = _labelize(node, nodehash)
    push!(updates, node => newnode)
    return nothing
end

"""
    wrap_to_store!(graph, node, nodehash;
                   cachedir=DEFAULT_CACHE_DIR,
                   compression=DEFAULT_COMPRESSION,
                   skipcache=false)

Generates a new `Op` node that corresponds to `node::DispatchNode`
and which stores the output of the execution of `node` in a file whose
name and extension depend on `nodehash` and `compression`. The generated
node is added to `updates` which maps `node` to the generated node. The
node output is stored in `cachedir`. The caching is skipped if `skipcache`
is `true`.
"""
function wrap_to_store!(updates::Dict{DispatchNode, DispatchNode},
                        node::DispatchNode,
                        nodehash::String;
                        cachedir::String=DEFAULT_CACHE_DIR,
                        compression::String=DEFAULT_COMPRESSION,
                        skipcache::Bool=false)
    # Defing exec-store wrapper function
    """
    Simple exec-store wrapper.
    """
    function exec_store_wrapper(args...; kwargs...)
        _cachedir = abspath(joinpath(expanduser(cachedir), DEFAULT_HASHCACHE_DIR))
        if compression != "none" && !skipcache
            operation = "EXEC-STORE-COMPRESS"
        elseif compression == "none" && !skipcache
            operation = "EXEC-STORE"
        else
            operation = "EXEC *ONLY*"
        end
        extension = ifelse(compression != "none", ".$compression", ".bin")
        filepath = joinpath(_cachedir, nodehash * extension)
        compressor = get_compressor(compression, "compress")

        # Get calculation result
        result = _run_node(node, args...; kwargs...)

        # Store result
        label = _labelize(node, nodehash)
        @debug "[$nodehash][$(label)] $operation (compression=$compression)"
        if !skipcache
            if !isfile(filepath)
                open(compressor, filepath, "w") do fid
                    serialize(fid, result)
                end
            else
                @debug "`-->[$nodehash][$(node.label)] * SKIPPING $operation"
            end
        end
        return result
    end

    # Update arguments and keyword arguments of node using the updates;
    # the latter should contain at this point only solved nodes (so the
    # dependencies of the current node should be good.
    newnode = Op(exec_store_wrapper)
    newnode.label = _labelize(node, nodehash)
    newnode.args = _arguments(node, updates)
    newnode.kwargs = _kwarguments(node, updates)
    # Add wrapped node to updates
    push!(updates, node=>newnode)
    return nothing
end

# Small wrapper that executes a node
_run_node(node::Op, args...; kwargs...) = node.func(args...; kwargs...)

_run_node(node::IndexNode, args...; kwargs...) = getindex(args..., node.index)

_run_node(node::CollectNode, args...; kwargs...) = vcat(args...)

_run_node(node::DataNode, args...; kwargs...) = identity(args...)

# Small wrapper that gets the label for a wrapped node
_labelize(node::Op, args...) = get_label(node)

_labelize(node::CollectNode, args...) = get_label(node)

_labelize(node::IndexNode, nodehash::String) = "IndexNode_$nodehash"

_labelize(node::DataNode, nodehash::String) = "DataNode_$nodehash"

# Small wrapper that generates new arguments for a wrapped node
_arguments(node::Op, updates) = map(node.args) do arg
    ifelse(arg isa DispatchNode, get(updates, arg, arg), arg)
end

_arguments(node::IndexNode, updates) =
    tuple(get(updates, node.node, node.node))

_arguments(node::CollectNode, updates) =
    Tuple(get(updates, n, n) for n in node.nodes)

_arguments(node::DataNode, updates) =
    tuple(ifelse(node.data isa DispatchNode,
                 get(updates, node.data, node.data),
                 node.data))

# Small wrapper that generates new keyword arguments for a wrapped node
_kwarguments(node::Op, updates) = pairs(
    NamedTuple{(node.kwargs.itr...,)}(
        ((map(node.kwargs.data) do kwarg
                 ifelse(kwarg isa DispatchNode, get(updates, kwarg, kwarg), kwarg)
             end)...,)))

_kwarguments(node::DispatchNode, updates) = pairs(NamedTuple())


text:
# ---
# title: 1221. Split a String in Balanced Strings
# id: problem1221
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String, Greedy
# link: <https://leetcode.com/problems/split-a-string-in-balanced-strings/description/>
# hidden: true
# ---
#
# _Balanced_  strings are those who have equal quantity of 'L' and 'R'
# characters.
#
# Given a balanced string `s` split it in the maximum amount of balanced
# strings.
#
# Return the maximum amount of splitted balanced strings.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "RLRRLLRLRL"
#     Output: 4
#     Explanation: s can be split into "RL", "RRLL", "RL", "RL", each substring contains same number of 'L' and 'R'.
#
#
# **Example 2:**
#
#
#
#     Input: s = "RLLLLRRRLR"
#     Output: 3
#     Explanation: s can be split into "RL", "LLLRRR", "LR", each substring contains same number of 'L' and 'R'.
#
#
# **Example 3:**
#
#
#
#     Input: s = "LLLLRRRR"
#     Output: 1
#     Explanation: s can be split into "LLLLRRRR".
#
#
# **Example 4:**
#
#
#
#     Input: s = "RLRRRLLRLL"
#     Output: 2
#     Explanation: s can be split into "RL", "RRRLLRLL", since each substring contains an equal number of 'L' and 'R'
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 1000`
#   * `s[i] = 'L' or 'R'`
#
#
## @lc code=start
using LeetCode

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


text:
using Distributed

@everywhere function host_id()
    return [gethostname() myid()]
end

vs = Any[]
for worker in workers()
    s = @spawnat worker host_id()
    push!(vs, s)
end

for s in vs
    a = fetch(s)
    println("HOST: $(a[1]),  ID: $(a[2])")
end

max_stars_repo_name: Soyukke/julia-paralell-examples

max_issues_repo_name: Soyukke/julia-paralell-examples

max_forks_repo_name: Soyukke/julia-paralell-examples


text:
#
# Define a global problem and ist constructors
#
# ---

"""
    Problem{T}

Describe the problem that should be optimized by stating all properties, that do not change
during an optimization or that are dependent of a certain solver.

The parameter `T` can be used to distinguish problems with different representations
or implementations.
The default parameter [`AllocatingEvaluation`](@ref), which might be slower but easier to use.
The usually faster parameter value is [`MutatingEvaluation`](@ref)

See [`Options`](@ref) for the changing and solver dependent properties.
"""
abstract type Problem{T} end

"""
    AbstractEvaluationType

An abstract type to specify the kind of evaluation a [`Problem`](@ref) supports.
"""
abstract type AbstractEvaluationType end

"""
    AllocatingEvaluation <: AbstractEvaluationType

A parameter for a [`Problem`](@ref) indicating that the problem uses functions that
allocate memory for their result, i.e. they work out of place.
"""
struct AllocatingEvaluation <: AbstractEvaluationType end

"""
    MutatingEvaluation

A parameter for a [`Problem`](@ref) indicating that the problem uses functions that
do not allocate memory but work on their input, i.e. in place.
"""
struct MutatingEvaluation <: AbstractEvaluationType end

"""
    get_cost(p, x)

evaluate the cost function `F` stored within a [`Problem`](@ref) at the point `x`.
"""
function get_cost(p::Problem, x)
    return p.cost(p.M, x)
end


text:
#import Base.start, Base.next, Base.done
#type Combinations
  #from::Int64
  #choose::Int64
  #total::Int64
  #Combinations(from, choose) = new(from, choose, binomial(from, choose))
#end
#typealias CombinationState Tuple{Array{Int64,1},Int64}
#
#start(x::Combinations) = (collect(1:x.choose),1)
#done(x::Combinations, s::CombinationState) = s[2] >= x.total
#
#inc(c, i) = [c[1:i-1]; (c[i]+1):(c[i]+1+length(c)-i)]
##
#function next_i(x::Combinations, c, i)
  #if c[i] < x.from - (i - x.choose)
    #inc(c, i)
  #else
    #next_i(x, c, i-1)
  #end
#end
#
#function next(x::Combinations, s::CombinationState)
  #ith::Array{Int64,1} = next_i(x, s[1], x.choose)
  #(ith, (ith, (s[2]+1)::Int64))
#end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
# List of available algorithms for AIBECS
algs = [CTKAlg]

@testset "Solvers" begin
    nt = length(T_all)
    n = nt * nb
    x₀ = p₀.xgeo * ones(n)
    @testset "$(string(alg))" for alg in algs
        @test alg <: DiffEqBase.AbstractSteadyStateAlgorithm
        @testset "$(nx₀)x₀, $(np₀)p₀" for nx₀ in 1:2, np₀ in 1:2
            prob = SteadyStateProblem(F, ∇ₓF, nx₀ * x₀, np₀ * p₀)
            s = solve(prob, alg())
            @test s isa SteadyStateSolution
            @test prob isa SteadyStateProblem
            @test norm(s.u) / norm(F(s.u, np₀ * p₀)) > ustrip(upreferred(1e5u"Myr"))
        end
    end
end


text:
using BSON, Dates, Distributions, FileIO, Plots, StaticArrays

using Crux, Flux, POMDPs, POMDPGym, POMDPSimulators, POMDPPolicies

using ParallelTreeSampling

include("utils.jl")

unicodeplots()

# Create a larger GridWorldMDP.
# Use a different set of rewards, especially with a relatively large negative
# reward far out from the center of the grid.
# Note that the default is (4,3)=>-10.0, (4,6)=>-5.0, (9,3)=>10.0, (8,8)=>3.0.
tprob = 0.7
rewards = Dict(GWPos(4,3)=>-10.0, GWPos(4,6)=>-5.0, GWPos(9,3)=>10.0,
               GWPos(8,8)=>3.0, GWPos(16,18)=>-30.0)
# Add small costs for other grids.
Random.seed!(0)
rand_costs = Dict(GWPos(i,j) => rand()
                  for i in 1:20, j in 1:20 if !haskey(rewards, GWPos(i,j)))
mdp = GridWorldMDP(size=(20,20), rewards=rewards, costs=rand_costs, cost_penalty=1.0)

# Learn a policy that solves it.
function train_dqn(mdp::GridWorldMDP)
    N = 5_000_000
    policy = DiscreteNetwork(
        Chain(x -> (x .- 10f0) ./ 10f0, Dense(2, 32, relu), Dense(32, 4)),
        [:up, :down, :left, :right])
    policy = solve(DQN(π=policy, S=state_space(mdp), N=N, ΔN=4,
                       buffer_size=10_000, log=(;period=round(Int, N/10))), mdp)
    atable = Dict(s=>action(policy, [s...]) for s in states(mdp.g))
    BSON.@save "demo/policies/gridworld_large_policy_table.bson" atable
end
# train_dqn(mdp)

atable = BSON.load("demo/policies/gridworld_large_policy_table.bson")[:atable]
function policy_fn(s)
    if 1 <= s[1] <= mdp.g.size[1] && 1 <= s[2] <= mdp.g.size[2]
        return atable[s][1]
    end
    xs = POMDPs.actions(mdp, s)
    return xs[rand(1:end)]
end

px = GenericDiscreteNonParametric(
    [GWPos(0,0),
     GWPos(1,0), GWPos(0,1), GWPos(-1,0), GWPos(0,-1),
     GWPos(2,0), GWPos(0,2), GWPos(-2,0), GWPos(0,-2)],
    [0.600,
     0.075, 0.075, 0.075, 0.075,
     0.025, 0.025, 0.025, 0.025])
function cost_fn(rmdp::RMDP, s, sp)
    amdp = rmdp.amdp
    cost = get(amdp.costs, s, 0)
    reward = POMDPs.reward(amdp, s)
    return reward < 0.0 ? -reward : amdp.cost_penalty * cost
end
rmdp = RMDP(amdp=mdp, π=FunctionPolicy(policy_fn), cost_fn=cost_fn, disturbance_type=:noise)

POMDPs.actions(mdp::RMDP) = px
POMDPs.actionindex(mdp::RMDP, x) = findfirst(px.support .== x)

fixed_s = GWPos(10,10)

N = 500_000
# base_N = 10_000_000
base_N = 1_000_000
c = 0.0
α = 0.001; β = 0.1; γ = 0.3;
vloss = 0.0
is_baseline = false

path = "data"

alpha_list = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5]
date_str = Dates.format(Dates.now(), "mmdd-HHMM")

function baseline(trial=0)
    baseline_out = run_baseline(rmdp, fixed_s, px; N=base_N)
    filename = string("gridworld_large_baseline_$(date_str)_$(base_N)",
                      (trial > 0 ? string("_", trial) : "") ,".jld2")
    save(joinpath(path, filename),
         Dict("risks:" => baseline_out[1], "states:" => baseline_out[2]))

    d = Dict([(alpha, Dict()) for alpha in alpha_list])
    for alpha in alpha_list
        m = eval_metrics(baseline_out[1]; alpha)
        d[alpha][:N] = length(baseline_out[1]); d[alpha][:mean] = round(m.mean, digits=3);
        d[alpha][:var] = round(m.var, digits=3); d[alpha][:cvar] = round(m.cvar, digits=3);
        d[alpha][:worst] = round(m.worst, digits=3);
    end
    return d
end

function mcts(trial=0)
    nominal_distrib_fn = (mdp, s) -> px
    mcts_out, planner = run_mcts(
        rmdp, fixed_s, nominal_distrib_fn; N=N, c=c, vloss=vloss, α=α, β=β, γ=γ)
    mcts_info = mcts_out[4]
    filename = string("gridworld_large_mcts_$(date_str)_$(N)",
                      (trial > 0 ? string("_", trial) : "") ,".jld2")
    save(joinpath(path, filename),
         Dict("risks:" => mcts_out[1], "states:" => mcts_out[2],
              "weights:" => mcts_out[3], "tree:" => mcts_info[:tree]))
    search_t = round(mcts_info[:search_time_s]; digits=3)

    d = Dict([(alpha, Dict()) for alpha in alpha_list])
    for alpha in alpha_list
        m = eval_metrics(mcts_out[1]; weights=exp.(mcts_out[3]), alpha=alpha)
        d[alpha][:N] = length(mcts_out[1]); d[alpha][:mean] = round(m.mean, digits=3);
        d[alpha][:var] = round(m.var, digits=3); d[alpha][:cvar] = round(m.cvar, digits=3);
        d[alpha][:worst] = round(m.worst, digits=3);
    end
    return d, search_t
end

metrics = []
times = []
num_trials = 3
for trial in 1:num_trials
    Random.seed!(trial)
    if is_baseline
        push!(metrics, baseline(trial))
    else
        metric_d, search_t = mcts(trial)
        push!(metrics, metric_d)
        push!(times, search_t)
    end
end

if is_baseline
    println("Baseline metrics")
else
    println("TIS metrics: N=$(N), c=$(c), α=$(α), β=$(β), γ=$(γ)")
end

for (idx, metric) in enumerate(metrics)
    println("Run $(idx)")
    for alpha in alpha_list
        m_d = metric[alpha]
        println(string("  [Alpha=$(alpha)] N: $(m_d[:N]), Mean: $(m_d[:mean]), ",
                       "VaR: $(m_d[:var]), CVaR: $(m_d[:cvar]), Worst: $(m_d[:worst])"))
    end
end
println("Times: $(times)")


text:
mutable struct XORSpace{T<:Integer}
    dimension::Int
    basis::Vector{T}
    function XORSpace{T}(n) where T
        return new(0, zeros(T, n))
    end
end
function appropriate_type(n)
    if n < 8sizeof(Int)
        return Int
    elseif n < 64
        return Int64
    elseif n < 128
        return Int128
    else
        return BigInt
    end
end
function XORSpace(n)
    return XORSpace{appropriate_type(n)}(n)
end
function Base.push!(xs::XORSpace{T}, x::T) where T
    for i in eachindex(xs.basis)
        if !iszero(x & (one(T) << (i - 1)))
            if iszero(xs.basis[i])
                xs.dimension += 1
                xs.basis[i] = x
                break
            else
                x = xor(x, xs.basis[i])
            end
        end
    end
    return xs
end
function orthogonal_complement(xs::XORSpace{T}) where T
    r = Vector{T}(undef, length(xs.basis) - xs.dimension)
    k = 0
    for i in eachindex(xs.basis)
        if iszero(xs.basis[i])
            e = one(T) << (i - 1)
            x = e
            for j in 1:(i - 1)
                if !iszero(xs.basis[j] & e)
                    x |= (one(T) << (j - 1))
                end
            end
            k += 1
            r[k] = x
        end
    end
    @assert k == length(r)
    return r
end
function xor_complement(x, n, ::Type{T} = appropriate_type(n)) where T
    xs = XORSpace{T}(n)
    for xi in x
        push!(xs, xi)
    end
    return orthogonal_complement(xs)
end


text:
mutable struct Gradient{T}
    L :: T
end

mutable struct Fista{T <: Real, AType <: AbstractArray{T}}
    xk :: AType
    xk_tmp :: AType
    tk :: T
    L :: T
    function Fista(A::AbstractArray{T}, L::Real, tk::Real=1) where {T <: Real}
        new{T, typeof(A)}(similar(A), copy(A), T(tk), T(L))
    end
end

mutable struct MFista{T <: Real, AType <: AbstractArray{T}}
    zk :: AType
    xk :: AType
    xk_tmp :: AType
    β :: T
    L :: T
    F :: Function
    F_tmp :: T
    function MFista(A::AbstractArray{T}, L::Real,
                    F::Function, β::Real=1) where {T <: Real}
        new{T, typeof(A)}(similar(A), similar(A), copy(A),
                          T(β), T(L), F, T(Inf))
    end
end

function Base.similar(t::NTuple{N,M where M}) where N
    ([ similar(t[i]) for i in 1:N]...,)
end

function optim_restart(method)
    nothing
end

function optim_restart(method::Fista{T,AType}) where {T <: Real,
                                                AType <: AbstractArray{T}}
    method.tk = 1
    nothing
end

function optim_restart(method::NTuple{N,M where M}) where N
    for i in 1:N
        optim_restart(method[i])
    end
    nothing
end

function gradient_step!(A, ∇A, prox!, method::Nothing)
    A
end

function gradient_step!(A, ∇A::Nothing, prox!, method)
    A
end

function gradient_step!(A::AbstractArray{T},
                        ∇A::AbstractArray{T},
                        prox!::Function,
                        method::Gradient) where {T <: Real}
    A .-= (1/method.L).*∇A
    prox!(A)
end

function gradient_step!(A::AbstractArray{T},
                        ∇A::AbstractArray{T},
                        prox!::Function,
                        method::Fista{T,AType}) where {T <: Real,
                                                       AType <: AbstractArray{T}}
    method.xk .= A .- (1/method.L).*∇A
    prox!(method.xk)
    tk = (1+sqrt(1+4*method.tk^2))/2
    β = (method.tk-1)/tk
    method.tk = tk
    A .= method.xk .+ β.*(method.xk .- method.xk_tmp)
    method.xk_tmp .= method.xk
end

function gradient_step!(A::AbstractArray{T},
                        ∇A::AbstractArray{T},
                        prox!::Function,
                        method::MFista{T,AType}) where {T <: Real,
                                                        AType <: AbstractArray{T}}
    # Different than original Monotone Fista
    method.zk .= A .- (1/method.L).*∇A
    prox!(method.zk)
    method.xk .= method.zk
    Fk = method.F(method.xk)
    if Fk <= method.F_tmp
        method.F_tmp = Fk
        A .= method.zk .+ method.β.*(method.zk .- method.xk_tmp)
        method.xk_tmp .= method.zk
    else
        A .= method.xk_tmp
    end
end

function proximal_step!(A::AbstractArray{T},
                        ∇A::AbstractArray{T},
                        env,
                        uf,
                        compute_fwd!::Function,
                        compute_grad!::Function,
                        method,
                        prox!::Union{Nothing,Function}=nothing) where {T <: Real}
    compute_grad!(env, A, ∇A, uf)
    proxop! = if prox! == nothing; identity else prox! end
    gradient_step!(A, ∇A, proxop!, method)
    uf = compute_fwd!(env, A)
    return uf
end

function proximal_step!(A::NTuple{N,AbstractArray{T}},
                        ∇A::NTuple{N,Union{Nothing,AbstractArray{T}}},
                        env,
                        uf,
                        compute_fwd!::Function,
                        compute_grad!::Function,
                        method::NTuple{N,M where M},
                        prox!::Union{Nothing,
                                     NTuple{N,Function}}=nothing) where {N, T <: Real}
    compute_grad!(env, A, ∇A, uf)
    for i in 1:N
        proxop! = if prox! == nothing; identity else prox![i] end
        gradient_step!(A[i], ∇A[i], proxop!, method[i])
    end
    uf = compute_fwd!(env, A)
    return uf
end

function proximal_iter!(A, ∇A, env,
                        compute_fwd!::Function,
                        compute_err!::Function,
                        compute_grad!::Function,
                        method, n::Int;
                        prox! = nothing, errors=[], print_errors=true,
                        restart=false)
    uf = compute_fwd!(env, A)
    err = compute_err!(env, uf)
    push!(errors, err)
    if print_errors
        println(err)
    end
    for k in 1:n
        err_tmp = err
        uf = proximal_step!(A, ∇A, env, uf,
                            compute_fwd!, compute_grad!, method, prox!)
        err = compute_err!(env, uf)
        push!(errors, err)
        if restart && err > err_tmp
            optim_restart(method)
        end
        if print_errors
            println(err)
        end
    end
    errors
end


text:
function solve_naive(input)
    earliest = earlier = parse(Int64, input[1])
    buses = parse.(Int64, filter(x -> x != "x", split(input[2], ',')))
    idx = findfirst(id -> earlier % id == 0, buses)
    while isnothing(idx)
        earlier += 1
        idx = findfirst(id -> earlier % id == 0, buses)
    end
    buses[idx] * (earlier - earliest)
end

# Taken from https://rosettacode.org/wiki/Chinese_remainder_theorem#Julia
function chineseremainder(n::Array, a::Array)
    Π = prod(n)
    mod(sum(ai * invmod(Π ÷ ni, ni) * (Π ÷ ni) for (ni, ai) in zip(n, a)), Π)
end

function solve_with_flashbacks(input)
    # (Flashbacks OMA) Esto sale con teorema chino del resto...
    buses = map(x -> x == "x" ? 0 : parse(Int64, x), split(input[2], ','))
    n = []
    a = []
    for i in 0:length(buses) - 1
        if buses[i + 1] > 0
            push!(n, -i)
            push!(a, buses[i + 1])
        end
    end
    chineseremainder(a, n)
end

@assert solve_naive(readlines("input/test/test_day13.txt")) == 295 ["test case 1 failed"]
@assert solve_with_flashbacks(readlines("input/test/test_day13.txt")) == 1068781 ["test case 2 failed"]

println("First star: " * string(solve_naive(readlines("input/day13.txt"))))
println("Second star: " * string(solve_with_flashbacks(readlines("input/day13.txt"))))


text:
"""
    expspace(x_min::Float64, x_max::Float64)::Array{Float64}, x_scale::Float64,n_x::Int64)
This generates a vector of length n_x between x_min and x_max,
where the spacing is controlled by x_scale.  When x_scale == 1,
this is equivalent to linspace.  When x_scale > 1, the points
are closer together near x_min.  When 0 < x_scale < 1,
the points are closer together near x_max.
"""
function expspace(x_min::Float64, x_max::Float64, x_scale::Float64,n_x::Int64)::Array{Float64}
    x_vec=zeros(n_x)
    if n_x>1
        step=(x_max-x_min)/(n_x-1)^x_scale
    else
        step=0.0
    end

    for i=1:n_x
        x_vec[i]=x_min+step*(i-1)^x_scale
    end

    return x_vec
end

"""
    tauchen(N::Integer=18, ρ::Real=0.96, σ::Real=0.045, μ::Real=0.0, n_std::Integer=4)
Tauchen's (1996) method for approximating AR(1) process with a finite markov chain
```math
    y_t = \\mu + \\rho y_{t-1} + \\epsilon_t
```
where \$\\epsilon_t \\sim N (0, \\sigma^2)\$
##### Arguments
- `N::Integer`: Number of points in markov process
- `ρ::Real` : Persistence parameter in AR(1) process
- `σ::Real` : Standard deviation of random component of AR(1) process
- `μ::Real(0.0)` : Mean of AR(1) process
- `n_std::Integer(4)` : The number of standard deviations to each side the process
  should span
##### Returns
- `grid, Π, Π_cdf` : The grid, the transition matrix and CDF of transition matrix
##### Notes
- This function has been modified by QuantEcon.jl
"""
function tauchen(N::Integer=9, ρ::Real=0.96, σ::Real=0.045, μ::Real=0.0, n_std::Integer=4)
    if N > 1
        a_bar = n_std * sqrt(σ^2 / (1 - ρ^2))
        y = range(-a_bar,stop= a_bar,length= N)
        d = y[2] - y[1]
    else
        a_bar = 0.0
        y = range(-a_bar,stop= a_bar,length= N)
    end
    # Construct grid

    grid = y .+ μ / (1 - ρ)

    # Get transition probabilities
    if N > 1
        if ρ > 0.0
            Π = zeros(N, N)
            for row = 1:N
    # Do end points first
            Π[row, 1] = cdf(Normal(),(y[1] - ρ*y[row] + d/2) / σ)
            Π[row, N] = 1 - cdf(Normal(),(y[N] - ρ*y[row] - d/2) / σ)

    # fill in the middle columns
                for col = 2:N-1
                Π[row, col] = (cdf(Normal(),(y[col] - ρ*y[row] + d/2) / σ) -
                           cdf(Normal(),(y[col] - ρ*y[row] - d/2) / σ))
                end
            end
            Π_cdf=cumsum(Π,dims=2)
        else
            Π = zeros(N)
            Π[1] = cdf(Normal(),(y[1] + d/2) / σ)
            Π[N] = 1 - cdf(Normal(),(y[N]  - d/2) / σ)
            for col=2:N-1
                Π[col] = (cdf(Normal(),(y[col]  + d/2) / σ) -
                       cdf(Normal(),(y[col] - d/2) / σ))
            end
            Π_cdf=cumsum(Π,dims=1)

        end
    else
        Π=[1.0]
        Π_cdf =[1.0]
    end

    return grid, Π, Π_cdf
end


text:
"""
    find_val_in_df(df::DataFrame,column_of_reference::Symbol,reference::String,value_to_return::Symbol)
Take DataFrame(df) Look in Column (column_of_reference) for the reference value (reference) and return in same row the value in column (value_to_return)
"""
function find_val_in_df(df::DataFrame,
                    column_of_reference::Symbol,
                    reference::String,
                    value_to_return::Symbol
                    )
                    @warn "find_val_in_df deprecated"
    return df[findfirst(df[column_of_reference].==reference),value_to_return]
end

"""
    find_val_in_df(df::DataFrame,column_of_reference::Symbol,reference::String,value_to_return::String)
Take DataFrame(df) Look in Column (column_of_reference) for the reference value (reference) and return corresponding value in column (value_to_return)
"""
function find_val_in_df(df::DataFrame,
                    column_of_reference::Symbol,
                    reference::String,
                    value_to_return::String
                    )
                    @warn "find_val_in_df deprecated"
    return find_val_in_df(df,column_of_reference,reference,Symbol(value_to_return))
end

#Use getindex to return all rows of the DataFrame fulfilling that in column `col` value `val` is found
function Base.getindex(df::DataFrame, col_and_val::Tuple{Symbol,Any})
    return df[findall(df[col_and_val[1]].==col_and_val[2]), :]
end

#Use getindex to return all rows of the DataFrame fulfilling that in column `col` value `val` is found and the column `colon_ind`
function Base.getindex(df::DataFrame, col_and_val::Tuple{Symbol,Any}, colon_ind::Symbol)
    return df[findall(df[col_and_val[1]].==col_and_val[2]), colon_ind]
end

"""
    check_column(df::DataFrame, names_array::Array{Symbol,1})
check if the columns provided in `names_array` exist in the DataFrame `df`
throw an error if they don't
"""
function check_column(df::DataFrame, names_array::Array{Symbol,1})
    for name in names_array
        #Check existance of necessary column
        name in names(df) || throw(@error "No column called `$name` in $(repr(df))")
    end
end

"""
     map_set_in_df(df::DataFrame,column_of_reference::Symbol,reference::String,set_to_return::Symbol)
  Take DataFrame(`df`) Look in Column (`column_of_reference`) for all cases that match the reference value (`reference`) and return the corresponding sets in Column (`set_to_return`)
"""
function map_set_in_df(df::DataFrame,
                    column_of_reference::Symbol,
                    reference::String,
                    set_to_return::Symbol
                    )
                    @warn "mat_set_in_df deprecated"
    return df[df[column_of_reference].==reference,set_to_return]
end

"""
    getindex(variable::OptVariable,index_set::Array)
Get the variable data from the specific Scenario by indicating the `var_name` e.g. "COST" and the `index_set` like `[:;"EUR";"pv"]`
"""
function get_cep_variable_value(variable::OptVariable,
                                index_set::Array
                                )
                                @warn "get_cep_variable_value deprecated"
    index_num=[]
    for i in  1:length(index_set)
        if index_set[i]==Colon()
            push!(index_num,Colon())
        elseif typeof(index_set[i])==Int || typeof(index_set[i])==UnitRange{Int}
            push!(index_num,index_set[i])
        else
            new_index_num=findfirst(variable.axes[i].==index_set[i])
            if new_index_num==[]
                throw(@error("$(index_set[i]) not in indexset #$i of Variable $var_name"))
            else
                push!(index_num,new_index_num)
            end
        end
    end
    return getindex(variable.data,Tuple(index_num)...)
end

"""
    get_cep_variable_value(scenario::Scenario,var_name::String,index_set::Array)
Get the variable data from the specific Scenario by indicating the `var_name` e.g. "COST" and the `index_set` like `[:;"EUR";"pv"]`
"""
function get_cep_variable_value(scenario::Scenario,
                                var_name::String,
                                index_set::Array
                                )
                                @warn "get_cep_variable_value deprecated"
    return get_cep_variable_value(scenario.opt_res.variables[var_name], index_set)
end

"""
    get_cep_variable_set(variable::OptVariable,num_index_set::Int)
Get the variable set from the specific variable and the `num_index_set` like 1
"""
function get_cep_variable_set(variable::OptVariable,
                              num_index_set::Int
                              )
                              @warn "get_cep_variable_set deprecated"
    return variable.axes[num_index_set]
end

"""
    get_cep_variable_set(scenario::Scenario,var_name::String,num_index_set::Int)
Get the variable set from the specific Scenario by indicating the `var_name` e.g. "COST" and the `num_index_set` like 1
"""
function get_cep_variable_set(scenario::Scenario,
                              var_name::String,
                              num_index_set::Int
                              )
                              @warn "get_cep_variable_set deprecated"
    return  get_cep_variable_set(scenario.opt_res.variables[var_name], num_index_set)
end

"""
    get_cep_design_variables(opt_result::OptResult)
Returns all design variables in this opt_result matching the type "dv"
Additionally you can add capacity factors, which scale the design variables by multiplying it with the value in the Dict
"""
function get_cep_design_variables(opt_result::OptResult)
  design_variables=get_cep_variables(opt_result, "dv")
  return design_variables
end

"""
    get_cep_slack_variables(opt_result::OptResult)
Returns all slack variables in this `opt_result` mathing the type "sv"
"""
function get_cep_slack_variables(opt_result::OptResult)
  return get_cep_variables(opt_result, "sv")
end

"""
    get_cep_variables(opt_result::OptResult, variable_type::String)
Returns all variables which types match the String of `variable_type`
"""
function get_cep_variables(opt_result::OptResult, variable_type::String)
  variables=Dict{String,Any}()
  for (key,val) in opt_result.variables
      if val.type==variable_type
          variables[key]=val
      end
  end
  if isempty(variables)
      throw(@error("$variable_type-Variable not provided in $(opt_result.descriptor)"))
  else
      return variables
  end
end

"""
    set_opt_config_cep(opt_data::OptDataCEP; kwargs...)
kwargs can be whatever you need to run the run_opt
it can hold
  -  `fixed_design_variables`: Dictionary{String,Any}
  -  `transmission`: true or false
  -  `generation`: true or false
  -  `storage_p`: true or false
  -  `storage_e`: true or false
  -  `existing_infrastructure`: true or false
  -  `descritor`: a String like "kmeans-10-co2-500" to describe this CEP-Model
  -  `first_stage_vars`: a Dictionary containing the OptVariables from a previous run
The function also checks if the provided data matches your kwargs options (e.g. it let's you know if you asked for transmission, but you have no tech with it in your data)
Returning Dictionary with the variables as entries
"""
function set_opt_config_cep(opt_data::OptDataCEP
                            ;kwargs...)
  # Create new Dictionary and set possible unique categories to false to later check wrong setting
  config=Dict{String,Any}("transmission"=>false, "storage_e"=>false, "storage_p"=>false, "generation"=>false)
  # Check the existence of the categ (like generation or storage - see techs.csv) and write it into Dictionary
  for categ in unique(getfield.(opt_data.techs[:], :categ))
    config[categ]=true
  end
  # Loop through the kwargs and write them into Dictionary
  for kwarg in kwargs
    # Check for false combination
    if String(kwarg[1]) in keys(config)
      if config[String(kwarg[1])]==false && kwarg[2]
        throw(@error("Option "*String(kwarg[1])*" cannot be selected with input data provided for "*opt_data.region))
      end
    end
    config[String(kwarg[1])]=kwarg[2]
  end

  # Return Directory with the information
  return config
end

"""
    set_opt_config_cep!(config::Dict{String,Any}; kwargs...)
add or replace items to an existing config:
- `fixed_design_variables`: `Dict{String,OptVariable}``
- `slack_cost`: Number
"""
function set_opt_config_cep!(config::Dict{String,Any}
                            ;kwargs...)
  # Loop through the kwargs and add them to Dictionary
  for kwarg in kwargs
    config[String(kwarg[1])]=kwarg[2]
  end

  # Return Directory with the information
  return config
end

"""
    check_opt_data_cep(opt_data::OptDataCEP)
Check the consistency of the data
"""
function check_opt_data_cep(opt_data::OptDataCEP)
  # Check lines
  # Only when Data provided
  if !isempty(opt_data.lines)
    # Check existence of start and end node
    for tech in axes(opt_data.lines,"tech")
      for node in getfield.(opt_data.lines[tech,:],:node_end)
        if !(node in axes(opt_data.nodes, "node"))
          throw(@error("Node "*node*" set as ending node, but not included in nodes-Data"))
        end
      end
    end
  end
end

"""
    get_total_demand(cep::OptModelCEP, ts_data::ClustData)
Return the total demand by multiplying demand with deltas and weights for the OptModel CEP
"""
function get_total_demand(cep::OptModelCEP,
                          ts_data::ClustData)
  ## DATA ##
  set=cep.set
  #ts          Dict( tech-node ): t x k
  ts=ts_data.data
  #ts_weights: k - weight of each period:
  ts_weights=ts_data.weights
  #ts_deltas:  t x k - Δt of each segment x period
  ts_deltas=ts_data.delta_t
  total_demand=0
  for node in set["nodes"]
    for t in set["time_T"]
      for k in set["time_K"]
        total_demand+=ts["el_demand-"*node][t,k]*ts_deltas[t,k]*ts_weights[k]
      end
    end
  end

  return total_demand
end

"""
    get_cost_series(cep_data::OptDataCEP,clust_res::ClustResultBest, opt_res::OptResult)
Return an array for the time series of costs in all the impact dimensions and the set of impacts
"""
function get_cost_series(nodes::DataFrame,
                        var_costs::DataFrame,
                       clust_res::ClustResultBest,
                       set::Dict{String,Array},
                       variables::Dict{String,OptVariable})
  ## DATA ##
  # ts_ids:   n_clustered periods
  ts_ids=clust_res.best_ids
  #ts_weights: k - weight of each period:
  ts_weights=clust_res.best_results.weights
  #ts_deltas:  t x k - Δt of each segment x period
  ts_deltas=clust_res.best_results.delta_t

  #emision at each period-step
  cost_ts=zeros(length(ts_ids)+1,length(set["impact"]))
  #At the beginning yearly fixed and cap costs
  cost_ts[1,:]=sum(get_cep_variable_value(variables["COST"],["cap_fix",:,:]),dims=2)

  for n in 1:length(ts_ids)
        var_cost=zeros(length(set["impact"]))
        i=1
        for impact in set["impact"]
          for tech in set["tech"]
            for node in set["nodes"]
              var_cost[i] += find_cost_in_df(var_costs,nodes,tech,node,impact)*  sum(get_cep_variable_value(variables["GEN"],["el",tech,:,ts_ids[n],node])' * ts_deltas[:,ts_ids[n]])
            end
          end
          i+=1
        end
        cost_ts[n+1,:]=cost_ts[n,:]+var_cost
  end
  return cost_ts, set["impact"]
end

"""
    get_cost_series(cep_data::OptDataCEP,scenario::Scenario)
Return an array for the time series of costs in all the impact dimensions and the set of impacts
"""
function get_cost_series(cep_data::OptDataCEP,scenario::Scenario)
  return get_cost_series(cep_data.nodes,cep_data.var_costs,scenario.clust_res, scenario.opt_res.model_set,scenario.opt_res.variables)
end

"""
    get_met_cap_limit(cep::OptModelCEP, opt_data::OptDataCEP, variables::Dict{String,OptVariable})
Return the technologies that meet the capacity limit
"""
function get_met_cap_limit(cep::OptModelCEP, opt_data::OptDataCEP, variables::Dict{String,Any})
  ## DATA ##
  # Set
  set=cep.set
  # nodes with limits
  nodes=opt_data.nodes

  met_cap_limit=Array{String,1}()
  for tech in set["tech_cap"]
    for node in set["nodes"]
      #Check if the limit is reached in any capacity at any node
      if sum(variables["CAP"][tech,:,node]) == nodes[tech,node].power_lim
        #Add this technology and node to the met_cap_limit Array
        push!(met_cap_limit,tech*"-"*node)
      end
    end
  end
  # If the array isn't empty throw an error (as limits are only for numerical speedup)
  if !isempty(met_cap_limit)
    #TODO change to warning
    throw( @error "Limit is reached for techs $met_cap_limit")
  end
  return met_cap_limit
end


text:
using Lucas, Base.Test

# strcut

@test strcut("ello", 1) == "…"
@test strcut("ello", 2) == "e…"
@test strcut("ello", 3) == "el…"
@test strcut("ello", 4) == "ello"

@test strcut("elo", 1, dots=false) == "e"
@test strcut("elo", 2, dots=false) == "el"
@test strcut("elo", 3, dots=false) == "elo"

@test strcut("ที่รัก", 1) == "…"
@test strcut("ที่รัก", 2) == "ที่…"
@test strcut("ที่รัก", 3) == "ที่รัก"
@test strcut("ที่รักa",3) == "ที่รั…"

@test strcut("ที่รัก", 1, dots=false) == "ที่"
@test strcut("ที่รัก", 2, dots=false) == "ที่รั"
@test strcut("ที่รัก", 3, dots=false) == "ที่รัก"
@test strcut("ที่รักa",3, dots=false) == "ที่รัก"


text:
@testset "Quaternion Midpoint" begin
	# Test quaternion sqrt
	q0 = rand(4)
	q0 ./= norm(q0)
	sq0 = sqrt_quat(q0)
	Q0 = Quaternion(q0...)
	sQ0 = sqrt(Q0)
	sq0_ = [sQ0.s, sQ0.v1, sQ0.v2, sQ0.v3]
	@test norm(sq0 - sq0_) < 1e-10

	q0 = [1,0,0,0]
	q0 ./= norm(q0)
	sq0 = sqrt_quat(q0)
	Q0 = Quaternion(q0...)
	sQ0 = sqrt(Q0)
	sq0_ = [sQ0.s, sQ0.v1, sQ0.v2, sQ0.v3]
	@test norm(sq0 - sq0_) < 1e-10

	# Test midpoint
	q0 = rand(4)
	q0 ./= norm(q0)
	ψ = 0.05 * rand(3)
	qψ = cayley_map(ψ)
	sqψ = sqrt_quat(qψ)
	q1 = L_multiply(qψ) * q0
	q2 = L_multiply(sqψ) * L_multiply(sqψ) * q0
	@test norm(q1 - q2) < 1e-10
	@test norm(R_multiply(q0)' * q1 - qψ) < 1e-10

	qmid = L_multiply(sqψ) * q0
	qmid_ = L_multiply(sqrt_quat(R_multiply(q0)' * q1)) * q0
	@test norm(qmid - qmid_) < 1e-10
end


text:
using SwitchOnSafety
using Documenter, Literate

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

const EXAMPLES = [
    "AJPR14e54.jl",
]

for example in EXAMPLES
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR)
    Literate.notebook(example_filepath, OUTPUT_DIR)
    Literate.script(example_filepath, OUTPUT_DIR)
end

makedocs(
    sitename = "SwitchOnSafety",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    # See https://github.com/jump-dev/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Index" => "index.md",
        "Invariant Sets" => "invariant.md",
        "Joint Spectral Radius" => "jsr.md",
        "Examples" => Any[
            "AJPR14e54" => "generated/AJPR14e54.md",
        ]
    ]
)

deploydocs(
    repo   = "github.com/blegat/SwitchOnSafety.jl.git",
    push_preview = true,
)


text:
@defstruct ChannelPoolingLayer Layer (
  name :: AbstractString = "channel-pooling",
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) > 0),
  (tops :: Vector{Symbol} = Symbol[], length(tops) == length(bottoms)),
  (kernel :: Int = 1, kernel > 0),
  (stride :: Int = 1, stride > 0),
  (pad :: NTuple{2, Int} = (0,0), all([pad...] .>= 0)),
  (channel_dim :: Int = -2, channel_dim != 0),
  pooling :: PoolingFunction = Pooling.Max(),
)
@characterize_layer(ChannelPoolingLayer,
  can_do_bp => true,
)

type ChannelPoolingLayerState <: LayerState
  layer      :: ChannelPoolingLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}

  op_dims    :: Vector{Int}
  etc        :: Any
end

function setup_etc(backend::CPUBackend, layer::ChannelPoolingLayer, inputs, blobs)
  if isa(layer.pooling, Pooling.Max)
    masks = Array(Array, length(inputs))
    for i = 1:length(inputs)
      masks[i] = Array(Csize_t, size(blobs[i]))
    end
    etc = masks
  elseif isa(layer.pooling, Pooling.Mean)
    integrals = Array(Array, length(inputs))
    for i = 1:length(inputs)
      integrals[i] = Array(eltype(inputs[i]), size(inputs[i])[1:end-1])
    end
    etc = integrals
  else
    etc = nothing
  end
  return etc
end

function setup(backend::Backend, layer::ChannelPoolingLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  pooled_chann_all = Array(Int, length(inputs))
  blobs = Array(Blob, length(inputs))
  blobs_diff = Array(Blob, length(inputs))
  op_dims = Array(Int, length(inputs))

  for i = 1:length(inputs)
    dim_total = ndims(inputs[i])
    op_dim = layer.channel_dim < 0 ? layer.channel_dim + dim_total+1 : layer.channel_dim
    @assert 1 <= op_dim <= dim_total
    @assert op_dim != dim_total

    op_dims[i] = op_dim

    dims = [size(inputs[i])...]
    pool_dim = dims[op_dim]
    pooled_dim = round(Int, ceil(float(pool_dim + layer.pad[1]+layer.pad[2] - layer.kernel) / layer.stride)) + 1

    # make sure the last pooling is not purely pooling padded area
    if ((pooled_dim-1)*layer.stride >= pool_dim + layer.pad[1])
      pooled_dim -= 1
    end
    pooled_chann_all[i] = pooled_dim

    output_dims = copy(dims)
    output_dims[op_dim] = pooled_dim
    output_dims = tuple(output_dims...)

    data_type = eltype(inputs[i])
    blobs[i] = make_blob(backend, data_type, output_dims)
    if isa(diffs[i], NullBlob)
      blobs_diff[i] = NullBlob()
    else
      blobs_diff[i] = make_blob(backend, data_type, output_dims)
    end
  end

  etc = setup_etc(backend, layer, inputs, blobs)
  state = ChannelPoolingLayerState(layer, blobs, blobs_diff, op_dims, etc)
end
function shutdown_etc(backend::CPUBackend, state::ChannelPoolingLayerState)
end
function shutdown(backend::Backend, state::ChannelPoolingLayerState)
  map(destroy, state.blobs)
  map(destroy, state.blobs_diff)
  shutdown_etc(backend, state)
end

function forward(backend::CPUBackend, state::ChannelPoolingLayerState, inputs::Vector{Blob})
  forward(backend, state.layer.pooling, state, inputs)
end
function forward(backend::CPUBackend, pool::StdPoolingFunction,
    state::ChannelPoolingLayerState, inputs::Vector{Blob})

  for i = 1:length(inputs)
    input = inputs[i].data
    output = state.blobs[i].data

    dims_in = split_dims(input, state.op_dims[i])
    dims_out = split_dims(output, state.op_dims[i])

    if isa(pool, Pooling.Max)
      max_channel_pooling_forward(reshape(input,dims_in), reshape(output,dims_out), reshape(state.etc[i],dims_out), state.layer)
    elseif isa(pool, Pooling.Mean)
      mean_channel_pooling_forward(reshape(input,dims_in), reshape(output,dims_out), state.etc[i], state.layer)
    end
  end
end

function backward(backend::CPUBackend, state::ChannelPoolingLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  backward(backend, state.layer.pooling, state, inputs, diffs)
end

function backward(backend::CPUBackend, pool::StdPoolingFunction, state::ChannelPoolingLayerState,
    inputs::Vector{Blob}, diffs::Vector{Blob})

  for i = 1:length(inputs)
    diff = diffs[i]

    if !isa(diff, NullBlob)
      dims_in = split_dims(inputs[i], state.op_dims[i])
      dims_out = split_dims(state.blobs[i], state.op_dims[i])

      if isa(pool, Pooling.Max)
        max_channel_pooling_backward(reshape(diff.data,dims_in), reshape(state.blobs_diff[i].data,dims_out),
            reshape(state.etc[i],dims_out), state.layer)
      elseif isa(pool, Pooling.Mean)
        mean_channel_pooling_backward(reshape(diff.data,dims_in), reshape(state.blobs_diff[i].data,dims_out), state.layer)
      end
    end
  end
end


text:
export Sites

"""
    Sites{L <: SiteLabel, T, N} <: AbstractSites{L, T, N}

Lattice Sites are Julia Arrays with certain Label, it is able to use array
interface.

TODO: type promote to normal arrays when mixed with normal arrays.
"""
mutable struct Sites{L <: SiteLabel, T, N} <: AbstractSites{L, T, N}
    data::Array{T, N}
end

# Enable type inference
Sites(::Type{L}, data::Array{T, N}) where {L, T, N} = Sites{L, T, N}(data)

Sites(::Type{L}, shape::Int...) where L <: SiteLabel = Sites(L, shape)
Sites(::Type{L}, shape::Tuple) where L <: SiteLabel =
    Sites{L, eltype(L), length(shape)}(fill(down(L), shape))

data(s::Sites) = s.data

# use array interface
import Base: eltype, length, ndims, size, eachindex,
    getindex, setindex!, stride, strides, copy
import Compat: axes

eltype(x::Sites{L, T}) where {L, T} = T
length(x::Sites) = length(x.data)
ndims(x::Sites) = ndims(x.data)
size(x::Sites) = size(x.data)
size(x::Sites, n::Integer) = size(x.data, n)
axes(x::Sites) = axes(x.data)
axes(x::Sites, d::Integer) = axes(x.data, d)
eachindex(x::Sites) = eachindex(x.data)
stride(x::Sites, k::Integer) = stride(x.data, k)
strides(x::Sites) = strides(x.data)
getindex(x::Sites, index::Integer...) = getindex(x.data, index...)
getindex(x::Sites, index::NTuple{N, T}) where {N, T <: Integer} = getindex(x.data, index...)
setindex!(x::Sites, val, index::Integer...) = setindex!(x.data, val, index...)
setindex!(x::Sites, val, index::NTuple{N, T}) where {N, T <: Integer} = setindex!(x.data, val, index...)

@inline function copy(b::Sites{L, T, N}) where {L, T, N}
    return Sites{L, T, N}(copy(b.data))
end

sitetype(::Type{Sites{L}}) where L = L
sitetype(x::Sites{L}) where L = L

export reset!, set!

@inline function reset!(b::Sites{L, T, N}) where {L, T, N}
    fill!(b.data, down(L))
    return b
end

@inline function set!(s::Sites{Bit, T, N}, n::Int) where {T, N}
    for i in eachindex(s)
        @inbounds s[i] = ((n>>(i-1)) & 0x1)
    end
    return s
end

@inline function set!(s::Sites{Spin, T, N}, n::Int) where {T, N}
    for i in eachindex(s)
        @inbounds s[i] = 2 * ((n>>(i-1)) & 0x1) - 1
    end
    return s
end

@inline function set!(s::Sites{Half, T, N}, n::Int) where {T, N}
    for i in eachindex(s)
        @inbounds s[i] = ((n>>(i-1)) & 0x1) - 0.5
    end
    return s
end

#######################
# randomize interface
#######################

function rand!(rng::AbstractRNG, b::Sites{L, T, N}) where {L, T, N}
    rand!(rng, b.data, [up(L), down(L)])
    return b
end

rand!(b::Sites{L, T, N}) where {L <: SiteLabel, T, N} = rand!(GLOBAL_RNG, b)
rand(rng::AbstractRNG, ::Type{L}, shape::Dims) where {L <: SiteLabel} =
    rand!(rng, Sites(L, shape))

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

@inline function flip!(b::Sites{L, T, N}, index::Integer...) where {L, T, N}
    if b[index...] == up(L)
        b[index...] = down(L)
    else
        b[index...] = up(L)
    end
    return b
end

@inline function randflip!(b::Sites{L, T, N}) where {L, T, N}
    offset = rand(1:length(b))
    return flip!(b, offset)
end

# carry bit
function carrybit!(a::Sites{L}) where L
    for i in eachindex(a)
        if a[i] == up(L)
            a[i] = down(L)
        else
            a[i] = up(L)
            break
        end
    end
end

# TODO: whether this should have side effect?
# will side effect affects performance (extra instance)?
import Base: <<
function <<(a::Sites, b::Int)
    for i = 1:b
        carrybit!(a)
    end
    return a
end

##############
# Conversions
##############

import Base: convert

function convert(::Type{Integer}, x::Sites{L, T, N}) where {L, T, N} end

for IntType in (:Int8, Int16, Int32, Int64, Int128, BigInt)
@eval begin
        function convert(::Type{$IntType}, x::Sites{Bit, T, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(Compat.InexactError(:convert, $IntType, x))
            end

            sum($IntType(each) << (i-1) for (i, each) in enumerate(x))
        end

        function convert(::Type{$IntType}, x::Sites{Spin, T, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(Compat.InexactError(:convert, $IntType, x))
            end

            sum($IntType(div(each+1, 2)) << (i-1) for (i, each) in enumerate(x))
        end

        function convert(::Type{$IntType}, x::Sites{Half, T, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(Compat.InexactError(:convert, $IntType, x))
            end

            sum($IntType(each+0.5) << (i-1) for (i, each) in enumerate(x))
        end
    end
end


text:
"""
setupitpab(window::T, N_itp_samples::Int, domain_percentage::T;
   N_fit_positions::Int = 15,
   p0 = [0.5; 0.0],
   p_lb = [0.1; -5.0],
   p_ub = [0.6; 5.0],
   max_iters = 5000,
   xtol_rel = 1e-5,
   ftol_rel = 1e-5,
   maxtime = Inf,
   optim_algorithm = :LN_BOBYQA) where T <: Real

window ∈ (0,1)
optim_algorithm can be :GN_ESCH, :GN_ISRES, :LN_BOBYQA, :GN_DIRECT_L
"""
function setupitpab(window::T, N_itp_samples::Int, domain_percentage::T;
    N_fit_positions::Int = 15,
    p0 = [0.5; 0.0],
    p_lb = [0.1; -5.0],
    p_ub = [0.6; 5.0],
    max_iters = 5000,
    xtol_rel = 1e-5,
    ftol_rel = 1e-5,
    maxtime = Inf,
    optim_algorithm = :LN_BOBYQA) where T <: Real

    # get piece-wise linear monotone maps.
    infos, zs, p_range = MonotoneMaps.getendomorphismpiecewiselinear(zero(T), one(T),
    window;
    N_itp_samples = N_itp_samples,
    domain_percentage = domain_percentage)

    # get compact sigmoid parameters fitted to each of the piece-wise linear maps.
    # costfuncs, minxs, rets = MonotoneMaps.getcompactsigmoidparameters(infos;
    # p0 = p0, p_lb = p_lb, p_ub = p_ub, optim_algorithm = optim_algorithm,
    # max_iters = max_iters, N_fit_positions = N_fit_positions, xtol_rel = xtol_rel,
    # ftol_rel = ftol_rel, maxtime = maxtime)
    costfuncs, minxs, rets = MonotoneMaps.getcompactsigmoidparameters(infos;
    N_fit_positions = N_fit_positions, max_iters = max_iters,
    xtol_rel = xtol_rel, ftol_rel = ftol_rel, maxtime = maxtime,
    p0 = p0, p_lb = p_lb, p_ub = p_ub, optim_algorithm = optim_algorithm)
    #qs = collect( tt->MonotoneMaps.evalcompositelogisticprobit(tt, minxs[i][1], minxs[i][2]) for i = 1:length(minxs) )

    Δp = p_range[2]-p_range[1]
    itp_range = p_range[1]:Δp:p_range[end]

    # itp.
    a_samples = collect( minxs[i][1] for i = 1:length(minxs) )
    b_samples = collect( minxs[i][2] for i = 1:length(minxs) )

    a_itp = Interpolations.interpolate(a_samples, Interpolations.BSpline(Interpolations.Cubic(Interpolations.Line(Interpolations.OnGrid()))))
    a_sitp = Interpolations.scale(a_itp, itp_range)
    a_setp = Interpolations.extrapolate(a_sitp, Interpolations.Flat()) # zero outside interp range.

    b_itp = Interpolations.interpolate(b_samples, Interpolations.BSpline(Interpolations.Cubic(Interpolations.Line(Interpolations.OnGrid()))))
    b_sitp = Interpolations.scale(b_itp, itp_range)
    b_setp = Interpolations.extrapolate(b_sitp, Interpolations.Flat()) # zero outside interp range.

    return a_setp, b_setp,
        minxs, rets # debug
end


text:
###
  Name     : c8_25_cluster.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 3/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using Clustering
srand(12345)
d = 10
n = 500
x = rand(d, n)
S = -pairwise(Euclidean(), x, x)
# set diagonal value to median value
S = S - diagm(diag(S)) + median(S)*eye(size(S,1))
R = affinityprop(S)


text:
include("trackingcell.jl");
include("lineage.jl");

data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF";
ret_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200316_result";
img = load(File(format"TIFF","$data_dir/20200316/HE7-11-1-80uw-PWM_1_s3.ome.tiff"));
img_edge, mask_edge, mask_clear = split_cell_LoG(img, size(img)[3]÷20);
save("$ret_dir/img_edge.tiff", img_edge);

nucleus_pos = [ component_centroids(label_components(mask_clear[:,:,i])) for i in 1:size(mask_clear)[3] ];

paths = generate_path(nucleus_pos);

for i in 1:length(paths)
	save("$ret_dir/s3_c$i.tiff", N0f16.(export_cell(paths[i], img, mask_edge)));
end


text:
abstract AbstractIndicatorMatrix

###Factory function. Currently implemented is :stationary
indicator{T<:AbstractFloat}(s::Symbol,nproposals::Int,nsamples::Int,::Type{T},args...) = _indicator(Val{s},nproposals,nsamples,T,args...)

indicatorsamples(indicator::AbstractIndicatorMatrix) = indicator.samples

function accepted(indicator::AbstractIndicatorMatrix)
    acc = 0
    @simd for i=2:length(indicator.samples)
        @inbounds indicator.samples[i]!=indicator.samples[i-1]?acc+=1:nothing
    end
    acc
end

numsamples(indicator::AbstractIndicatorMatrix) = length(indicator.samples) - 1


text:
module VisualTests

using Gadfly, Interpolations
const nx = 10
xg = 1:nx
xf = -1.2nx:.1:3.3nx

f1(x) = sin(2pi/nx * (x-1))

y1 = f1.(xg)

p = plot()

if true

Btypes = (Periodic, Flat, Line, Free, Reflect)
Gtypes = (OnCell, OnGrid)
degrees = (
    Constant(), Linear(),
    [Quadratic(T(G())) for T in Btypes, G in Gtypes]...,
    [Cubic(T(G())) for T in Btypes[1:end-1], G in Gtypes]...,  # no Reflect for Cubic
)
Etypes = (Flat, Line, Reflect, Periodic)

for deg in degrees, ET in Etypes
    itp = extrapolate(interpolate(y1, BSpline(deg)), ET())

    stuff = Any[]

    push!(stuff, layer(x=xg,y=y1,Geom.point,Theme(default_color=colorant"green")))
    push!(stuff, layer(x=xf,y=[itp[x] for x in xf],Geom.path,Theme(point_size=2px)))
    title = "$deg, $(nameof(ET))"
    push!(stuff, Guide.title(title))
    display(plot(stuff...))
end

nx = 10
xg = 1:nx
ny = 6
yg = 1:ny

f2(x,y) = sin(2pi/nx * (x-1)) * cos(2pi/ny * (y-1))
zg = f2.(xg, yg')

xf = -1:.1:nx+1
yf = -1:.1:ny+1

itp2 = extrapolate(interpolate(zg, BSpline(Quadratic(Flat(OnCell())))), Line())

display(plot(
    layer(x=xf,y=yf,z=[itp2[x,y] for x in xf, y in yf], Geom.contour),
    Guide.title("Quadratic(Flat(OnCell)), Line")
))

end

end


text:
export
    Quat,
    get_axis,
    get_rotation_angle,

    RPY

immutable Quat
    x::Float64
    y::Float64
    z::Float64
    w::Float64
end

Base.length(::Quat) = 4
Base.copy(a::Quat) = Quat(a.x, a.y, a.z, a.w)
Base.convert(::Type{Vector{Float64}}, a::Quat) = [a.x, a.y, a.z, a.w]
function Base.convert{R<:Real}(::Type{Quat}, a::AbstractArray{R})
    @assert(length(a) == 4)
    Quat(a[1], a[2], a[3], a[4])
end
Base.show(io::IO, q::Quat) = @printf(io, "QUAT({%6.3f, %6.3f, %6.3f}, %6.3f)", q.x, q.y, q.z, q.w)

Base.abs(q::Quat) = sqrt(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w)
function Base.norm(q::Quat)
    n = abs(q)
    Quat(q.x/n, q.y/n, q.z/n, q.w/n)
end

function get_axis(q::Quat)

    θ = 2*atan2(sqrt(q.x*q.x + q.y*q.y + q.z*q.z), q.w)
    h = sin(θ/2)

    x = q.x / h
    y = q.y / h
    z = q.z / h

    VecE3(x, y, z)
end
get_rotation_angle(q::Quat) = 2.0*acos(q.w)

"""
Roll Pitch Yaw
"""
immutable RPY
    r::Float64
    p::Float64
    y::Float64
end
function Base.convert(::Type{RPY}, q::Quat)

    q2 = norm(q)
    x = q.x
    y = q.y
    z = q.z
    w = q.w

    roll  = atan2(y*z+w*x, 0.5-x^2-y^2)
    pitch = asin(-2*(x*z + w*y))
    yaw   = atan2(x*y+w*z, 0.5-y^2-z^2)

    RPY(roll, pitch, yaw)
end
function Base.convert(::Type{Matrix{Float64}}, quat::Quat)

    # convert it to a rotation matrix

    w = quat.w
    x = quat.x
    y = quat.y
    z = quat.z

    n = w * w + x * x + y * y + z * z
    s = (n == 0.0) ? 0.0 : 2.0 / n
    wx = s * w * x; wy = s * w * y; wz = s * w * z
    xx = s * x * x; xy = s * x * y; xz = s * x * z
    yy = s * y * y; yz = s * y * z; zz = s * z * z

    [ 1 - (yy + zz)         xy - wz          xz + wy;
           xy + wz     1 - (xx + zz)         yz - wx;
           xz - wy          yz + wx     1 - (xx + yy) ]
end


text:
using Test
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
using LinearAlgebraicRepresentation.Arrangement
Lar = LinearAlgebraicRepresentation
Lara = LinearAlgebraicRepresentation.Arrangement

@testset "Pipeline Arrangement 1" begin
	V = [0.0 1.0; 1.0 1.0; 3.0 1.0; 4.0 1.0; 2.0 2.0; 2.0 0.0];
	EV = SparseArrays.sparse(Array{Int8, 2}([
		[1 0 0 1 0 0] #1->1,4
		[0 1 1 0 0 0] #2->2,3
		[0 0 0 0 1 1] #3->5,6
		[0 1 0 0 1 0] #4->2,5
		[0 0 1 0 0 1] #5->3,6
	]))
	bigPI = Lar.spaceindex((convert(Lar.Points, V'), Lar.cop2lar(EV)));

	@testset "intersect_edges" begin
		@test Lara.intersect_edges(V, EV[1, :], EV[3, :]) == [([2.0 1.0], 0.5)];
		@test Lara.intersect_edges(V, EV[1, :], EV[3, :]) == Lara.intersect_edges(V, EV[2, :], EV[3, :])
		@test Lara.intersect_edges(V, EV[1, :], EV[2, :]) == [([1.0 1.0], 0.25); ([3.0 1.0], 0.75)];
		@test Lara.intersect_edges(V, EV[2, :], EV[1, :]) == Lara.intersect_edges(V, EV[4, :], EV[5, :]);
		@test Lara.intersect_edges(V, EV[1, :], EV[4, :]) == [([1.0 1.0], 0.25)];
		@test Lara.intersect_edges(V, EV[2, :], EV[4, :]) == [([1.0 1.0], 0.0)];

	end

	@testset "merge_vertices" begin
		V1, EV1 = Lara.frag_edge(V, EV, 1, bigPI);
		ans1 = SparseArrays.sparse(Array{Int8, 2}([
				[1 0 1 0 0] #1->1,3
				[0 0 1 0 1] #2->3,5
				[0 0 0 1 1] #3->4,5
				[0 1 0 1 0] #4->2,4
			]));
		ans2 = SparseArrays.sparse(Array{Int8, 2}([
				[1 0 1] #1->1,3
				[0 1 1] #2->2,3
			]));
		@test Lara.merge_vertices!(V1, EV1)[1] == [0.0 1.0; 4.0 1.0; 1.0 1.0; 3.0 1.0; 2.0 1.0];
		@test Lara.merge_vertices!(V1, EV1)[2] == ans1;

		@test Lara.merge_vertices!(V1, EV1, [[-1]], 1)[1] == [0.0 1.0; 4.0 1.0; 2.0 1.0];
		@test Lara.merge_vertices!(V1, EV1, [[-1]], 1)[2] == ans2;
	end

	@testset "Complete Pipeline 1" begin
		Vans = [0.0 1.0; 4.0 1.0; 1.0 1.0; 3.0 1.0; 2.0 1.0; 2.0 2.0; 2.0 0.0];
		EVans = SparseArrays.sparse(Array{Int8, 2}([
				[1 0 1 0 0 0 0] # 1->3
				[0 0 1 0 1 0 0] # 3->4
				[0 0 0 1 1 0 0] # 4->5
				[0 1 0 1 0 0 0] # 2->4
				[0 0 0 0 1 1 0] # 5->6
				[0 0 0 0 1 0 1] # 5->7
				[0 0 1 0 0 1 0] # 3->6
				[0 0 0 1 0 0 1] # 4->7
			]));
		@test Lara.planar_arrangement_1(V, EV)[1] == Vans
		@test Lara.planar_arrangement_1(V, EV)[2] == EVans
	end
end

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

@testset "Biconnected Components" begin
	@testset "Normal Execution" begin
		EV = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 0 0 0 0 0 0 0 0] #1 -> 1,2		|
			[0 1 1 0 0 0 0 0 0 0] #2 -> 2,3		|
			[1 0 1 0 0 0 0 0 0 0] #3 -> 3,1		|
			[0 0 1 1 0 0 0 0 0 0] #4 -> 3,4		 |
			[0 0 1 0 0 1 0 0 0 0] #5 -> 3,6		 |
			[0 0 1 0 0 0 1 0 0 0] #6 -> 3,7		-
			[0 0 0 1 1 0 0 0 0 0] #7 -> 4,5		 |
			[0 0 0 0 1 1 0 0 0 0] #8 -> 5,6		 |
			[0 0 0 0 0 0 0 1 1 0] #9 -> 8,9		  |
			[0 0 0 0 0 0 0 0 1 1] #10-> 9,10	  |
			[0 0 0 0 0 0 0 1 0 1] #11-> 8,10	  |
		]));
		bicomp = sort(Lara.biconnected_components(EV));

		@test length(Lara.biconnected_components(EV)) == 3;
		@test sort(bicomp[1]) == [1, 2, 3];
		@test sort(bicomp[2]) == [4, 5, 7, 8];
		@test sort(bicomp[3]) == [9, 10, 11];
	end

	@testset "Void Components" begin
		EV = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 0 0] #1 -> 1,2
			[1 0 1 0] #2 -> 1,3
			[1 0 0 1] #3 -> 1,4
		]));

		@test length(Lara.biconnected_components(EV)) == 0;
	end

	@testset "Couple of Components" begin
		EV = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 0] #1 -> 1,2  |
			[1 1 0] #2 -> 1,2  |
			[1 0 1] #3 -> 1,2
		]));

		@test length(Lara.biconnected_components(EV)) == 1;
		@test sort(Lara.biconnected_components(EV)[1]) == [1, 2];
	end
end

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

@testset "Pipeline Arrangement 2" begin

	@testset "Containment Boxes on Nested Triangular" begin
		W = [
			0.0 0.0; 8.0 0.0; 4.0 6.0;
			2.0 2.0; 6.0 2.0; 4.0 5.0;
			8.0 1.0; 5.0 6.0; 8.0 6.0;
			3.0 3.0; 5.0 3.0; 4.0 4.0];
		copEV = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 0 0 0 0 0 0 0 0 0 0] # 1 -> 1,2
			[0 1 1 0 0 0 0 0 0 0 0 0] # 2 -> 2,3
			[1 0 1 0 0 0 0 0 0 0 0 0] # 3 -> 3,1
			[0 0 0 1 1 0 0 0 0 0 0 0] # 4 -> 4,5
			[0 0 0 0 1 1 0 0 0 0 0 0] # 5 -> 5,6
			[0 0 0 1 0 1 0 0 0 0 0 0] # 6 -> 4,6
			[0 0 0 0 0 0 1 1 0 0 0 0] # 7 -> 7,8
			[0 0 0 0 0 0 0 1 1 0 0 0] # 8 -> 8,9
			[0 0 0 0 0 0 1 0 1 0 0 0] # 9 -> 7,9
			[0 0 0 0 0 0 0 0 0 1 1 0] #10 -> 10,11
			[0 0 0 0 0 0 0 0 0 0 1 1] #11 -> 11,12
			[0 0 0 0 0 0 0 0 0 1 0 1] #12 -> 10,12
		]));

		EVs = [
			SparseArrays.sparse(Array{Int8, 2}([
				[-1  1 0 0 0 0 0 0 0 0 0 0] #1 -> 1,2
				[ 0 -1 1 0 0 0 0 0 0 0 0 0] #2 -> 2,3
				[-1  0 1 0 0 0 0 0 0 0 0 0] #3 -> 3,1
			])),
			SparseArrays.sparse(Array{Int8, 2}([
				[0 0 0 -1  1 0 0 0 0 0 0 0] #4 -> 4,5
				[0 0 0  0 -1 1 0 0 0 0 0 0] #5 -> 5,6
				[0 0 0 -1  0 1 0 0 0 0 0 0] #6 -> 4,6
			])),
			SparseArrays.sparse(Array{Int8, 2}([
				[0 0 0 0 0 0 -1  1 0 0 0 0] #7 -> 7,8
				[0 0 0 0 0 0  0 -1 1 0 0 0] #8 -> 8,9
				[0 0 0 0 0 0 -1  0 1 0 0 0] #9 -> 7,9
			])),
			SparseArrays.sparse(Array{Int8, 2}([
				[0 0 0 0 0 0 0 0 0 -1  1 0] #10 -> 10,11
				[0 0 0 0 0 0 0 0 0  0 -1 1] #11 -> 11,12
				[0 0 0 0 0 0 0 0 0 -1  0 1] #12 -> 10,12
			]))
		];

		shells = [
			SparseArrays.sparse(Array{Int8}([-1; -1;  1])),
			SparseArrays.sparse(Array{Int8}([-1; -1;  1])),
			SparseArrays.sparse(Array{Int8}([-1; -1;  1])),
			SparseArrays.sparse(Array{Int8}([-1; -1;  1])),
		];

		bicon_comps = [[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]];
		shell_bboxes = [
			([0.0 0.0], [8.0 6.0]),
			([2.0 2.0], [6.0 5.0]),
			([5.0 1.0], [8.0 6.0]),
			([3.0 3.0], [5.0 4.0])
		];
		precontainment = sparse(Array{Int8,2}([
			[0 0 0 0]
			[1 0 0 0]
			[1 0 0 0]
			[1 1 0 0]
		]));
		containment_graph = sparse(Array{Int8,2}([
			[0 0 0 0]
			[1 0 0 0]
			[0 0 0 0]
			[1 1 0 0]
		]));
		transitive = sparse(Array{Int8,2}([
			[0 0 0 0]
			[1 0 0 0]
			[0 0 0 0]
			[0 1 0 0]
		]));

		@test Lara.pre_containment_test(shell_bboxes) == precontainment;
		@test Lara.prune_containment_graph(3, W, EVs, shells, containment_graph) == containment_graph;

		Lar.Arrangement.transitive_reduction!(containment_graph)

		@test containment_graph == transitive;
	end

	@testset "Pipeline Arrangement 2 over K^4" begin

		# K^4 + edge
		W = [0.0 0.0; 4.0 0.0; 2.0 3.0; 2.0 1.5];
		copEV = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 0 0] #1 -> 1,2
			[0 1 1 0] #2 -> 2,3
			[1 0 1 0] #3 -> 3,1
			[1 0 0 1] #4 -> 1,2
			[0 1 0 1] #5 -> 2,3
			[0 0 1 1] #6 -> 3,1
		]));

		copFE = SparseArrays.sparse(Array{Int8, 2}([
			[1 0 0 1 1 0] #1 -> 1,4,5
			[0 1 0 0 1 1] #2 -> 2,5,6
			[0 0 1 1 0 1] #3 -> 3,4,6
			[1 1 1 0 0 0] #4 -> 1,2,3 External
		]));

		bicon_comps = [[1, 2, 3, 4, 5, 6]];

		@testset "External Cycle" begin
			@test Lara.get_external_cycle(W, copEV, copFE) == 4;
			@test typeof(Lara.get_external_cycle(W, copEV, copFE[1:3,:])) == Nothing
		end

		@testset "Biconnected Components" begin
			@test length(Lara.biconnected_components(copEV)) == 1;
			@test sort(Lara.biconnected_components(copEV)[1]) == bicon_comps[1];
		end

		@testset "Complete Component Graph" begin
			EVs = [SparseArrays.sparse(Array{Int8, 2}([
				[-1  1  0  0] #1 -> 1,2
				[ 0 -1  1  0] #2 -> 2,3
				[-1  0  1  0] #3 -> 3,1
				[-1  0  0  1] #4 -> 1,2
				[ 0 -1  0  1] #5 -> 2,3
				[ 0  0 -1  1] #6 -> 3,1
			]))];

			boundaries = [SparseArrays.sparse(Array{Int8, 2}([
				[ 1  0  0 -1  1  0] #1 -> 1,4,5
				[ 0  1  0  0 -1  1] #2 -> 2,5,6
				[ 0  0 -1  1  0 -1] #3 -> 3,4,6
			]))];

			@test Lara.componentgraph(W, copEV, bicon_comps)[1] == 1;
			@test size(Lara.componentgraph(W, copEV, bicon_comps)[2]) == (1, 1);
			@test Lara.componentgraph(W, copEV, bicon_comps)[3] == W;
			@test Lara.componentgraph(W, copEV, bicon_comps)[4] == EVs;
			@test Lara.componentgraph(W, copEV, bicon_comps)[5] == boundaries;
			@test Lara.componentgraph(W, copEV, bicon_comps)[6] == [sparse(Array{Int8}([-1; -1; 1; 0; 0; 0]))];
			@test Lara.componentgraph(W, copEV, bicon_comps)[7] == [([0.0 0.0], [4.0 3.0])];
		end
	end

	@testset "Pipeline Arrangement 2 over planar_tile" begin
		W = [
			0.0 0.0; 0.4 0.0; 0.5 0.0; 0.6 0.0; 1.0 0.0;
			0.0 0.4; 1.0 0.4;
			0.0 0.5; 1.0 0.5;
			0.0 0.6; 1.0 0.6;
			0.0 1.0; 0.4 1.0; 0.5 1.0; 0.6 1.0; 1.0 1.0;
			0.5 0.5
		];
		copEV = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 0 0 0  0 0  0 0  0 0  0 0 0 0 0  0]	#1  ->  1  2
			[1 0 0 0 0  1 0  0 0  0 0  0 0 0 0 0  0]	#2  ->  1  6
			[0 1 0 0 0  1 0  0 0  0 0  0 0 0 0 0  0]	#3  ->  2  6

			[0 0 0 1 1  0 0  0 0  0 0  0 0 0 0 0  0]	#4  ->  4  5
			[0 0 0 1 0  0 1  0 0  0 0  0 0 0 0 0  0]	#5  ->  4  7
			[0 0 0 0 1  0 1  0 0  0 0  0 0 0 0 0  0]	#6  ->  5  7

			[0 0 0 0 0  0 0  0 0  1 0  1 0 0 0 0  0]	#7  -> 10 12
			[0 0 0 0 0  0 0  0 0  1 0  0 1 0 0 0  0]	#8  -> 10 13
			[0 0 0 0 0  0 0  0 0  0 0  1 1 0 0 0  0]	#9  -> 12 13

			[0 0 0 0 0  0 0  0 0  0 1  0 0 0 1 0  0]	#10 -> 11 15
			[0 0 0 0 0  0 0  0 0  0 1  0 0 0 0 1  0]	#11 -> 11 16
			[0 0 0 0 0  0 0  0 0  0 0  0 0 0 1 1  0]	#12 -> 15 16

			[0 0 1 0 0  0 0  1 0  0 0  0 0 0 0 0  0]	#13 ->  3  8
			[0 0 1 0 0  0 0  0 1  0 0  0 0 0 0 0  0]	#14 ->  3  9
			[0 0 1 0 0  0 0  0 0  0 0  0 0 0 0 0  1]	#15 ->  3 17
			[0 0 0 0 0  0 0  1 0  0 0  0 0 0 0 0  1]	#16 ->  8 17
			[0 0 0 0 0  0 0  1 0  0 0  0 0 1 0 0  0]	#17 ->  8 14
			[0 0 0 0 0  0 0  0 1  0 0  0 0 0 0 0  1]	#18 ->  9 17
			[0 0 0 0 0  0 0  0 1  0 0  0 0 1 0 0  0]	#19 ->  9 14
			[0 0 0 0 0  0 0  0 0  0 0  0 0 1 0 0  1]	#20 -> 14 17
		]));

		copFE = SparseArrays.sparse(Array{Int8, 2}([
			[1 1 1  0 0 0  0 0 0  0 0 0  0 0 0 0 0 0 0 0]	#1 -> 1,2,3
			[0 0 0  1 1 1  0 0 0  0 0 0  0 0 0 0 0 0 0 0]	#2 -> 4,5,6
			[0 0 0  0 0 0  1 1 1  0 0 0  0 0 0 0 0 0 0 0]	#3 -> 7,8,9
			[0 0 0  0 0 0  0 0 0  1 1 1  0 0 0 0 0 0 0 0]	#4 -> 10,11,12
			[0 0 0  0 0 0  0 0 0  0 0 0  1 0 1 1 0 0 0 0]	#5 -> 13,15,16
			[0 0 0  0 0 0  0 0 0  0 0 0  0 1 1 0 0 1 0 0]	#6 -> 14,15,18
			[0 0 0  0 0 0  0 0 0  0 0 0  0 0 0 1 1 0 0 1]	#7 -> 16,17,20
			[0 0 0  0 0 0  0 0 0  0 0 0  0 0 0 0 0 1 1 1]	#8 -> 18,19,20
			[1 1 1  1 1 1  1 1 1  1 1 1  1 1 0 0 1 0 1 0]	#9 -> external
		]));

		bicon_comps = [
			[ 1,  2,  3],
			[ 4,  5,  6],
			[ 7,  8,  9],
			[10, 11, 12],
			[13, 14, 15, 16, 17, 18, 19, 20]
		];

		bboxes = [
			([0.0 0.0], [0.4 0.4])
			([0.6 0.0], [1.0 0.4])
			([0.0 0.6], [0.4 1.0])
			([0.6 0.6], [1.0 1.0])
			([0.0 0.0], [1.0 1.0])
		];

		@testset "External Cycle" begin
			@test Lara.get_external_cycle(W, copEV, copFE) == 9;
			@test typeof(Lara.get_external_cycle(W, copEV, copFE[1:8,:])) == Nothing
		end

		@testset "Biconnected Components" begin
			@test length(Lara.biconnected_components(copEV)) == 5;
			@test sort(sort(Lara.biconnected_components(copEV))[1]) == bicon_comps[1];
			@test sort(sort(Lara.biconnected_components(copEV))[2]) == bicon_comps[2];
			@test sort(sort(Lara.biconnected_components(copEV))[3]) == bicon_comps[3];
			@test sort(sort(Lara.biconnected_components(copEV))[4]) == bicon_comps[4];
			@test sort(sort(Lara.biconnected_components(copEV))[5]) == bicon_comps[5];
		end

		@testset "Containment Boxes" begin
			pre_containment = SparseArrays.sparse(Array{Int8, 2}([
				[0 0 0 0 1]
				[0 0 0 0 1]
				[0 0 0 0 1]
				[0 0 0 0 1]
				[0 0 0 0 0]
			]));
			@test Lara.pre_containment_test(bboxes) == pre_containment;
		end

		@testset "Complete Component Graph" begin
			shells = [
				sparse(Array{Int8}([-1;  1; -1])),
				sparse(Array{Int8}([-1;  1; -1])),
				sparse(Array{Int8}([ 1; -1;  1])),
				sparse(Array{Int8}([ 1; -1;  1])),
				sparse(Array{Int8}([ 1; -1;  0;  0;  1;  0; -1; 0]))
			];

			@test Lara.componentgraph(W, copEV, bicon_comps)[1] == 5;
			@test Lara.componentgraph(W, copEV, bicon_comps)[2] == SparseArrays.spzeros(Int8, 5, 5);
			@test Lara.componentgraph(W, copEV, bicon_comps)[3] == W;
			# @test abs(Lara.componentgraph(W, copEV, bicon_comps)[4]) == EV;
			# @test abs(Lara.componentgraph(W, copEV, bicon_comps)[5]) == boundaries;
			@test Lara.componentgraph(W, copEV, bicon_comps)[6] == shells;
			@test Lara.componentgraph(W, copEV, bicon_comps)[7] == bboxes;
		end

	end
end


text:
mutable struct Cascade{T<:Sampling}
    start::T
    stop::T
    steps::Vector{T}
    permutation::Vector
    correlation::AbstractMatrix
    ispermuted::Bool
    iscorrelated::Bool
end

function Cascade( ; start, stop, steps, correlation, permutation, ispermuted, iscorrelated)
    return Cascade(start, stop, steps, permutation, correlation, ispermuted, iscorrelated)
end

"""
    collect_data(x::Cascade)
This function returns a list of all [`Data`](@ref) Types stored in `x`,
ordered [`x.start; x.steps, x.stop]`.
"""
collect_data(x::Cascade) = Vector{}([x.start; x.steps; x.stop])

"""
    collect_correlation(x::Cascade)
This function extends the correlation matrix ``A \\in \\mathbb{R}^{N,N}`` stored in
`x.correlation` by defining ``A' \\in \\mathbb{R}^{N+2,N+2}`` to account for
`x.start` and `x.stop` for consistent dimensionality with the vector returned by
[`collect_data`](@ref).

```math
A' =
\\begin{pmatrix}
1 &   &     \\
    & A &   \\
    &   & 1
\\end{pmatrix}
```

# Returns
- `L::Matrix`:
"""
function collect_correlation(x::Cascade; kwargs...)
    mat = x.correlation
    N = size(mat,1)

    L = 1.0*Matrix(I(N+2))
    L[2:N+1, 2:N+1] .= mat
    return L
end

"""
    collect_permutation(lst)
    collect_permutation(cascade)
"""
collect_permutation(lst::AbstractArray) = [1; get_permutation(lst).+1; length(lst)+2]

function collect_permutation(x::T; kwargs...) where T<:Cascade
    return collect_permutation(get_permutation(x))
end

"""
    collect_ticks(x::Cascade{Data}; kwargs...)
This function returns a list of major y-axis ticks
"""
function collect_ticks(cascade; kwargs...)
    vmin, vmax, vscale = vlim(cascade; kwargs...)
    return collect(vmin:minimum(get_order.([vmin,vmax])):vmax)
end

"""
"""
collect_value(cascade::Cascade) = get_value(collect_data(cascade))

"""
    get_permutation(idx)
    get_permutation(cascade)
"""
function get_permutation(x::Cascade)
    return x.ispermuted ? get_permutation(x.permutation) : collect(1:length(x.steps))
end

function get_permutation(idx::AbstractVector{Int})
    ordered = DataFrames.DataFrame(old=idx,)
    collapsed = DataFrames.DataFrame(old=sort(idx), new=1:length(idx))
    df = DataFrames.leftjoin(collapsed, ordered, on=:old)
    return df[:,:new]
end

"""
    get_correlation(cascade)
"""
function get_correlation(x::Cascade; kwargs...)
    return x.correlation[get_permutation(x; kwargs...), :]
end

"""
    update_stop!(x::)
This function updates `cascade.stop` or the last row of a matrix to equal the negative
of the cumulative sum of all previous values.
"""
function update_stop!(x::Cascade)
    val = collect_value(x)
    set_value!(x.stop, update_stop!(val)[end,:])
    return x
end

function update_stop!(x::Matrix)
    x[end,:] .= -Statistics.cumsum(x; dims=1)[end-1,:]
    return x
end

function update_stop!(df::DataFrames.DataFrame; value=VALUE_COL)
    df[:,value] .= -update_stop!(df[:,value])
    return df
end


text:
function toplanetext(data::String)
   return data
end

formats = Dict{String,Function}(
  "*/*"              =>  toplanetext
  , "text/plain"       =>  toplanetext
  , ""       =>  toplanetext
)


text:
struct MacroContext
    source::LineNumberNode
    module_::Module
end

mapcat(f, xs) = collect(MapCat(f), xs)
firstsomething(f, xs) = foldl(right, xs |> Map(f) |> ReduceIf(!isnothing); init = nothing)
ifnothing(f) = x -> x === nothing ? f() : x
donothing(_...) = nothing

@inline _any() = false
@inline _any(x) = x
@inline _any(x, xs...) = x || _any(xs...)
@inline _any(a, b) = a || b
@inline _any(a, b, c) = a || b || c
@inline _any(a, b, c, d) = a || b || c || d
@inline _any(a, b, c, d, e) = a || b || c || d || e
@inline _any(a, b, c, d, e, f) = a || b || c || d || e || f
@inline _any(a, b, c, d, e, f, g) = a || b || c || d || e || f || g
@inline _any(a, b, c, d, e, f, g, h) = a || b || c || d || e || f || g || h
@inline _any(a, b, c, d, e, f, g, h, i) = a || b || c || d || e || f || g || h || i
@inline _any(a, b, c, d, e, f, g, h, i, j) = a || b || c || d || e || f || g || h || i || j

function flattenblockargs(ex)
    @match ex begin
        Expr(:block, args...) => mapcat(flattenblockargs, args)
        _ => [ex]
    end
end

function is_function(ex)
    @match ex begin
        Expr(:function, _...) => true
        Expr(:(=), Expr(:call, _...), _) => true
        Expr(:(=), Expr(:where, _...), _) => true
        Expr(:do, _...) => true
        Expr(:->, _...) => true
        _ => false
    end
end

is_rebinding_update(ex) =
    ex isa Expr && length(ex.args) == 2 && is_rebinding_update_op(ex.head)

is_dot_update(ex) = ex isa Expr && length(ex.args) == 2 && is_dot_update_op(ex.head)

function is_rebinding_update_op(sym::Symbol)
    s = String(sym)
    endswith(s, "=") || return false
    op = Symbol(s[1:end-1])
    Base.isbinaryoperator(op) || return false
    return !isdefined(Base, sym)  # handle `<=` etc.
end

function is_dot_update_op(sym::Symbol)
    s = String(sym)
    startswith(s, ".") && endswith(s, "=") || return false
    op = Symbol(s[2:end-1])
    return Base.isbinaryoperator(op)
end

is_dot_op(_) = false
function is_dot_op(sym::Symbol)
    s = String(sym)
    startswith(s, ".") || return false
    op = Symbol(s[2:end])
    return Base.isoperator(op)
end

is_dotcall(ex) =
    if isexpr(ex, :call)
        length(ex.args) > 1 && is_dot_op(ex.args[1])
    else
        isexpr(ex, :., 2) && isexpr(ex.args[2], :tuple)
    end

is_dotcall(ex, nargs) =
    if isexpr(ex, :call)
        length(ex.args) == nargs + 1 && is_dot_op(ex.args[1])
    else
        isexpr(ex, :., 2) && isexpr(ex.args[2], :tuple, nargs)
    end


text:
using IsPurelyImmutable
using Test

@testset "is_purely_immutable literals" begin
    @test is_purely_immutable(1)
    @test is_purely_immutable(true)
    @test is_purely_immutable(nothing)

    @test is_purely_immutable((2,3))

    @test is_purely_immutable(1:2:10)

    # Test custom extension for strings & symbols
    @test is_purely_immutable("hi")
    @test is_purely_immutable(:hey)

    @test is_purely_immutable((:hi,"hey"))
    @test is_purely_immutable((a=:hi, b="hey"))

    # Things that aren't pure (because they're mutable)
    @test !is_purely_immutable([])
    @test !is_purely_immutable(Int[1,2,3])
    @test !is_purely_immutable(Dict(1=>2))
    @test !is_purely_immutable(:(2+3))
end

struct EmptyImmutable end
mutable struct EmptyMutable end
struct ImmutableStructMutableFields
    x::Vector
end
mutable struct MutableStructMutableFields
    x::Vector
end
struct RecursivelyImmutable
    x::Union{RecursivelyImmutable,Nothing}
    RecursivelyImmutable(x = nothing) = new(x)
end
struct ImmutableStructAbstractField
    x  # Might be mutable, might be immutable
end
@testset "is_purely_immutable() with custom types" begin
    @test is_purely_immutable(EmptyImmutable())
    @test is_purely_immutable(EmptyMutable())
    @test !is_purely_immutable(ImmutableStructMutableFields([]))
    @test !is_purely_immutable(MutableStructMutableFields([]))
    @test is_purely_immutable(RecursivelyImmutable())
    @test is_purely_immutable(RecursivelyImmutable(RecursivelyImmutable()))

    # Purity of immutable types with Abstract Fields deepends on their value.
    @test is_purely_immutable(ImmutableStructAbstractField(1))
    @test is_purely_immutable(ImmutableStructAbstractField("hi"))
    @test !is_purely_immutable(ImmutableStructAbstractField([]))
end

# Overriding is_purely_immutable of custom type: Mutable Until Shared example
mutable struct MutableStructUntilShared
    v::Int
    mutable::Bool
    MutableStructUntilShared(x, mut=true) = new(x, mut)
end

function Base.setproperty!(m::MutableStructUntilShared, f::Symbol, v)
    if !m.mutable
        error("setproperty! Cannot modify MutableStructUntilShard `m` once it's been" *
              " marked immutable")
    else
        setfield!(m,f,v)
    end
end

IsPurelyImmutable.is_purely_immutable(x::MutableStructUntilShared) = !x.mutable

@testset "is_purely_immutable() for Mutable Until Shared struct" begin
    @test !is_purely_immutable(MutableStructUntilShared(1))
    @test is_purely_immutable(MutableStructUntilShared(1, false))

    x = MutableStructUntilShared(0)
    @test !is_purely_immutable(x)
    x.v = 1

    x.mutable = false  # Mark immutable
    @test_throws ErrorException x.v = 2
    @test is_purely_immutable(x)
end