text:
# "compatibility matrix"
# M₀[α, β] = 1 if any only if:
#     deg(β ∈ graph) ≥ deg(α ∈ subgraph)
#      and
#     species(α ∈ subgraph) == species(β ∈ graph)
function compatibility_matrix(subgraph::MetaGraph, subgraph_species::Array{Symbol, 1},
                              graph::MetaGraph, graph_species::Array{Symbol, 1},
                              disconnected_component::Bool)::Array{Bool, 2}
    # allocate M. rows correspond to subgraph nodes, columns to graph nodes.
    M₀ = zeros(Bool, nv(subgraph), nv(graph))

    # Get adjacency matrices and 4th/8th degree path matrices
    adjmat_S = adjacency_matrix(subgraph)
    adjmat_G = adjacency_matrix(graph)
    deg_S = adjmat_S^2
    deg_G = adjmat_G^2
    path4_S = deg_S^2
    path4_G = deg_G^2
    path8_S = path4_S^2
    path8_G = path4_G^2

    if !disconnected_component # search for substructures
        @inbounds for α ∈ 1:nv(subgraph) # Loop over rows (subgraph nodes)
            @inbounds for β ∈ 1:nv(graph) # Loop over columns (graph nodes)
                # Record Bool for each (i,j): true if atom species match, graph node degree is sufficient, and 4 and 8 length self-paths are sufficient.
                M₀[α, β] = subgraph_species[α] == graph_species[β] && deg_G[β, β] ≥ deg_S[α, α] && path4_G[β, β] ≥ path4_S[α, α] && path8_G[β, β] ≥ path8_S[α, α]
            end
        end
    else # search only for exact, isolated matches (no substructures)
        @inbounds for α ∈ 1:nv(subgraph) # Loop over rows (subgraph nodes)
            @inbounds for β ∈ 1:nv(graph) # Loop over columns (graph nodes)
                # Record Bool for each (i,j): true if atom species match, and graph node degree matches.
                M₀[α, β] = subgraph_species[α] == graph_species[β] && deg_G[β, β] == deg_S[α, α]
            end
        end
    end
    return M₀
end

# list of nodes β ∈ graph that could possibly correpond with node α ∈ subgraph
function candidate_list(M::Array{Bool, 2}, α::Int)::Array{Int, 1}
    return [β for β ∈ 1:size(M, 2) if M[α, β]]
end

# does node α have possible candidate matches in the graph?
function has_candidates(M::Array{Bool, 2}, α::Int)::Bool
    @inbounds for β ∈ 1:size(M, 2) # loop over graph nodes
        if M[α, β]
            return true
        end
    end
    # if made it this far, no graph node could possibly correspond to α ∈ subgraph
    return false
end

function is_isomorphism(M::Array{Bool, 2})::Bool
    # (1) each row of M, corresponding to a node α ∈ subgraph, contains exactly one 1.
    #     i.e., every subgraph node has exactly one correspondence
    @inbounds for α ∈ 1:size(M, 1)
        if sum(M[α, :]) != 1
            return false
        end
    end
    # (2) no column of M, corresponding to a node β ∈ graph, contains more than one 1.
    #     i.e., a graph node does not correspond to more than 1 subgraph node.
    @inbounds for β ∈ 1:size(M, 2)
        if sum(M[:, β]) > 1
            return false
        end
    end
    # if made it this far, it is indeed an isomorphism.
    return true
end

# idea here:
#   if any subgraph node α has no possible correspondence w/ a node β in the graph, no point in continuing
#   return true iff M has no empty candidate lists for subgraph nodes.
function possibly_contains_isomorphism(M::Array{Bool, 2})::Bool
    @inbounds for α ∈ 1:size(M, 1) # loop over subgraph nodes
        if ! has_candidates(M, α)
            return false # subgraph node α cannot be assigned!
        end
    end
    return true # M may be the intersection of one or more solutions.
end

function prune!(M::Array{Bool, 2}, subgraph::MetaGraph, graph::MetaGraph)
    pruned = true # to enter while loop
    while pruned
        pruned = false
        @inbounds for α ∈ 1:size(M, 1) # loop thru subgraph nodes
            # get neighbors of node α
            neighbors_of_α = neighbors(subgraph, α)
            # loop thru candidate matches β ∈ graph for this subgraph node α
            @inbounds for β ∈ candidate_list(M, α)
                neighbors_of_β = neighbors(graph, β)
                # now, suppose α ∈ subgraph and β ∈ graph correspond...
                @inbounds for x ∈ neighbors_of_α
                    # if there is no neighbor of β that could correspond to x, neighbor of α, then, contradiction.
                    if isdisjoint(candidate_list(M, x), neighbors_of_β)
                        M[α, β] = false
                        pruned = true
                    end
                end
            end
        end
    end
end

function assign_correspondence!(M::Array{Bool, 2}, α::Int, β::Int)
    M[α, :] .= false # zero out row of subgraph node
    M[:, β] .= false # zero out column of graph node
    M[α, β] = true # assign correspondence at intersection
end

# soln:
#    soln[α ∈ subgraph] = β ∈ graph where α corresponds to β
function depth_first_search!(α::Int, subgraph::MetaGraph, graph::MetaGraph,
                            M::Array{Bool, 2}, soln::Array{Int, 1},
                            β_mapped::Array{Bool, 1}, solns::Array{Array{Int, 1}, 1})
    # if reached here from previous solution, exit.
    if α > size(M, 1)
        return nothing
    end

    # loop thru un-assigned graph nodes β that could possibly correspond to subnode α
    @inbounds for β ∈ candidate_list(M, α)
        # if βraph is already mapped, not a viable solution.
        #   (not sure if necessary, i.e. if M already knows this)
        if β_mapped[β]
            continue
        end

        # make a copy so we can restore later.
        M′ = deepcopy(M)

        # explore scenario where α ∈ subgraph corresponds to β ∈ graph
        assign_correspondence!(M′, α, β)
        soln[α] = β
        β_mapped[β] = true

        # prune tree
        prune!(M′, subgraph, graph)

        # if we reached bottom of tree, iso-morphism is found!
        if α == size(M′, 1)
            # do we hv to check if it's a sol'n or is it guarenteed? why prune then?
            if is_isomorphism(M′)
                push!(solns, deepcopy(soln))
            end
            # don't return b/c we need to look at other candidates
        end

        if M′[α, β] && possibly_contains_isomorphism(M′)
            # we've assigned α, go deeper in the depth first search
            depth_first_search!(α + 1, subgraph, graph, M′, soln, β_mapped, solns)
        end

        β_mapped[β] = false
        soln[α] = 0
    end
end

@doc raw"""
returns an array of arrays, each containing one unique subgraph isomorphism
"""
function find_subgraph_isomorphisms(subgraph::MetaGraph, subgraph_species::Array{Symbol, 1},
                                    graph::MetaGraph, graph_species::Array{Symbol, 1}, disconnected_component::Bool=false)
    # store list of solutions here
    solns = Array{Array{Int, 1}, 1}()
    # encodes an isomorhism. maps α ∈ subgraph --> β ∈ graph
    soln = [0 for _ = 1:nv(subgraph)]
    # tell us which β ∈ graph are mapped already.
    #   entry β true iff β mapped
    β_mapped = [false for _ = 1:nv(graph)]
    # initial compatability matrix based on degrees of nodes and species
    M₀ = compatibility_matrix(subgraph, subgraph_species, graph, graph_species, disconnected_component)
    depth_first_search!(1, subgraph, graph, M₀, soln, β_mapped, solns)
    return solns
end


text:
if ~(VERSION < v"1.1-")
    include("interval_special_functions.jl")
end

include("arithmetic.jl")
include("convex_increasing.jl")
include("concave_increasing.jl")
include("mixed_convexity.jl")
include("extrema.jl")
include("division.jl")
include("multiplication.jl")
include("power.jl")
include("other.jl")
include("activation_functions.jl")
include("set_bounds.jl")
include("comparison.jl")
include("trilinear.jl")
include("apriori_mult.jl")

#joinpath(@__DIR__, "no_gradient", "arithmetic.jl")
#joinpath(@__DIR__, "no_gradient", "convex_increasing.jl")
#joinpath(@__DIR__, "no_gradient", "mixed_convexity.jl")

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using .PyPlot
export plotarea, plotarea_dem

"""
    plotit(x, y, dem)
    plotit(x, y, waterflows_output)

Plot DEM, uparea, flow-dir
"""
plotit(x, y, dem) = plotit(x, y, waterflows(dem))
function plotit(x, y, waterflows_output::Tuple)
    area, slen, dir, nout, nin, pits  = waterflows_output

    fig, axs = subplots(3,1)
    sca(axs[1])
    mycontour(x, y, dem)
    colorbar()
    sca(axs[2])
    heatmap(x, y, log10.(area))
    sca(axs[3])
    plotdir(x, y, dir)
end

function plotdir(x, y, dir)
    vecfield = dir2vec.(dir)
    vecfieldx = [v[1] for v in vecfield]
    vecfieldy = [v[2] for v in vecfield]
    quiver(repeat(x,1, length(y)), repeat(y,length(x),1), vecfieldx, vecfieldy)
end

pits2inds(pits) = ([p.I[1] for p in pits],
                   [p.I[2] for p in pits])
pits2vecs(x, y, pits) = (x[[p.I[1] for p in pits if p!=CartesianIndex(-1,-1)]],
                           y[[p.I[2] for p in pits if p!=CartesianIndex(-1,-1)]])

"""
    plotarea(x, y, area, pits; prefn=log10)

Plot uparea, or another variable.  If no pits should be plotted, use `[]`.
"""
function plotarea(x, y, area, pits; prefn=log10, cbar=true)
    px, py = pits2vecs(x, y, pits)
    fig, axs = subplots()
    heatmap(x, y, prefn.(area), cbar=cbar)
    scatter(px, py, 1, "r")
    fig
end

function plotarea_dem(x, y, dem, area, pits; levels=50, threshold=1/100, prefn=log10)
    fig, axs = subplots()
    mycontour(x, y, dem, levels=levels)

    threshold *= prod(size(area))
    area = copy(area)
    area[area.<threshold] .= NaN
    px, py = pits2vecs(x, y, pits)

    heatmap(x, y, prefn.(area))
    scatter(px, py, 1, "r")
end

function plotbnds(x,y,bnds)
    for b in bnds
        px, py = pits2vecs(x, y, b)
        scatter(px, py, 1, "g")
    end
end

function plotlakedepth(x, y, dem)
    area, slen, dir, nout, nin, pits  = waterflows(dem)
    demf = fill_dem(dem, pits, dir)
    heatmap(x, y, demf.-dem)
end

function heatmap(x, y, mat; cbar=true)
    dx2 = (x[2]-x[1])/2
    imshow(Array(mat'), origin="lower", extent=(x[1]-dx2,x[end]+dx2,y[1]-dx2,y[end]+dx2))
    cbar && colorbar()
end

function mycontour(x, y, mat; levels=50, cbar=true)
    contour(x, y, Array(mat'), levels, linewidths=1)
    cbar && colorbar()
end

function plot_catchments(x, y, c, lim=10)
    c = copy(c)
    tokeep = []
    for cc =1:maximum(c)
        inds = c.==cc
        if sum(inds)<lim
            c[inds] .= 0
        else
            push!(tokeep, inds)
        end
    end
    for (i,inds) in enumerate(tokeep)
        c[inds] .= i
    end
    heatmap(x,y,c)
end


text:
struct PipelineInputAssemblyStateCreateInfo
    mHandleRef::Ref{vk.VkPipelineInputAssemblyStateCreateInfo}

    function PipelineInputAssemblyStateCreateInfo(;
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineInputAssemblyStateCreateFlags = vk.VkFlags(0),
        topology::vk.VkPrimitiveTopology = vk.VK_PRIMITIVE_TOPOLOGY_TRIANGLE_LIST,
        primitiveRestartEnable::vk.VkBool32 = VkExt.VK_FALSE
    )

        this = new(vk.VkPipelineInputAssemblyStateCreateInfo(
            vk.VK_STRUCTURE_TYPE_PIPELINE_INPUT_ASSEMBLY_STATE_CREATE_INFO, #sType::VkStructureType
            pNext, #::Ptr{Cvoid}
            flags, #::VkPipelineInputAssemblyStateCreateFlags
            topology, #::VkPrimitiveTopology
            primitiveRestartEnable #::VkBool32
        ))
    end
end

function handleRef(this::PipelineInputAssemblyStateCreateInfo)::Ref{vk.VkPipelineInputAssemblyStateCreateInfo}
    this.mHandleRef
end


text:
# verify that a LinearMap makes it through the chain

using Pseudospectra, LinearMaps, Test

@testset "Linear map" begin

    A = LinearMap(Pseudospectra.convdiff_fd(10))

    # d,nconv,niter,nmult,resid = Pseudospectra.xeigs(A,I; nev=20, ncv=60,
    #                                                 which=:SR )

    opts = Dict{Symbol,Any}(:ARPACK_plots => false,:npts=>20)
    opts[:arpack_opts] = ArpackOptions{eltype(A)}(which=:SR,nev=6,ncv=60)

    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
end


text:
function CallTest()
    UnitRange{Int64}(2,2)
end

step_through(enter_call_expr(:($(CallTest)())))

# Properly handle :meta annotations
function MetaTest()
    @Base._pure_meta
    0
end

step_through(enter_call_expr(:($(MetaTest)())))

# Test Vararg handling
function VATest(x...)
    x
end
callVA() = VATest()

step_through(enter_call_expr(:($(VATest)())))

# Test Loops
function LoopTest()
    x = Int[]
    for i = 1:2
        push!(x, i)
    end
    x
end

step_through(enter_call_expr(:($(LoopTest)())))

# Test continue
function ContinueTest()
    x = Int[]
    for i = 1:3
        if true
            push!(x, i)
            continue
        end
        error("Fell through")
    end
    x
end

step_through(enter_call_expr(:($(ContinueTest)())))

#foo() = 1+1
function foo(n)
    x = n+1
    ((BigInt[1 1; 1 0])^x)[2,1]
end

step_through(enter_call_expr(:($(foo)(20))))

# Make sure that Symbols don't get an extra QuoteNode
function foo_sym()
    x = :ok
    typeof(x)
end

@test step_through(enter_call_expr(:($(foo_sym)()))) == Symbol

# Make sure evalling "new" works with symbols

function new_sym()
  a = :a
  () -> a
end

step_through(enter_call_expr(:($new_sym())))


text:
# abstract type AbstractPlanarPoint{T} <: StaticArrays.FieldVector{2, T} end
# abstract type AbstractPlanarVector{T} <: StaticArrays.FieldVector{2, T} end
#
# struct PlanarPoint{T} <: AbstractPlanarPoint{T}
#     u::T
#     v::T
# end

"""
    struct CartesianPoint{T} <: AbstractCoordinatePoint{T, Cartesian}

Describes a three-dimensional point in Cartesian coordinates.

## Fields
* `x`: x-coordinate (in m).
* `y`: y-coordinate (in m).
* `z`: z-coordinate (in m).

See also [`CylindricalPoint`](@ref).
"""
struct CartesianPoint{T} <: AbstractCoordinatePoint{T, Cartesian}
    x::T
    y::T
    z::T
end

zero(PT::Type{<:AbstractCoordinatePoint{T}}) where {T} = PT(zero(T),zero(T),zero(T))

"""
    struct CylindricalPoint{T} <: AbstractCoordinatePoint{T, Cylindrical}

Describes a three-dimensional point in cylindrical coordinates.

## Fields
* `r`: Radius (in m).
* `φ`: Polar angle (in rad).
* `z`: `z`-coordinate (in m).

!!! note
    `φ == 0` corresponds to the `x`-axis in the Cartesian coordinate system.

See also [`CartesianPoint`](@ref).
"""
struct CylindricalPoint{T} <: AbstractCoordinatePoint{T, Cylindrical}
    r::T
    φ::T
    z::T
    CylindricalPoint{T}(r::T, φ::T, z::T) where {T} = new(r, mod(φ,T(2π)), z)
    CylindricalPoint{T}(r::Real, φ::Real, z::Real) where {T} = new(T(r), mod(T(φ),T(2π)), T(z))
end

function CylindricalPoint(pt::CartesianPoint{T})::CylindricalPoint{T} where {T}
    return CylindricalPoint{T}(hypot(pt.x, pt.y), atan(pt.y, pt.x), pt.z)
end

function CartesianPoint(pt::CylindricalPoint{T})::CartesianPoint{T} where {T}
    sφ::T, cφ::T = sincos(pt.φ)
    return CartesianPoint{T}(pt.r * cφ, pt.r * sφ, pt.z)
end

@inline CylindricalPoint(pt::CylindricalPoint) = pt
@inline CartesianPoint(pt::CartesianPoint) = pt

@inline _convert_point(pt::AbstractCoordinatePoint, ::Type{Cylindrical}) = CylindricalPoint(pt)
@inline _convert_point(pt::AbstractCoordinatePoint, ::Type{Cartesian}) = CartesianPoint(pt)

# function _Δφ(φ1::T, φ2::T)::T where {T}
#     δφ = mod(φ2 - φ1, T(2π))
#     min(δφ, T(2π) - δφ)
# end
#
# _φNear(φ::Real, φMin::T, φMax::T) where {T} = _Δφ(T(φ),φMin) ≤ _Δφ(T(φ),φMax) ? φMin : φMax


text:
#=
interp1.jl
emulate matlab's interp1()
=#

export interp1

using Plots: scatter, plot!, gui
using Interpolations
using Test: @test, @test_throws, @inferred

"""
    yi = interp1(x, y, xi ; how=Gridded(Linear()), extrap=0)

1D interpolation of `y = f(x)` at points `xi`

In:
- `x::AbstractVector{<:Real}`
- `y::AbstractVector{<:Number}`

Option:
- `how::Interpolations.InterpolationType` default `Gridded(Linear())`
- `extrap::Any` how to extrapolate, e.g., `Flat()`; default `0`
other options from Interpolations.jl are `Line()` `Periodic()` `Reflect()` `Throw()`

Output is same size as input `xi`
"""
function interp1(x::AbstractVector{<:Real}, y::AbstractVector{<:Number}, xi ;
		how::Interpolations.InterpolationType = Gridded(Linear()),
		extrap::Any = 0)
#		extrap::Union{<:Interpolations.BoundaryCondition, <:Number} = 0 # fails?

	fun = interpolate((x,), y, how)
	if !isnothing(extrap)
		fun = extrapolate(fun, extrap)
	end
	fun.(xi)
end

"""
    interp1(:test)
self test
"""
function interp1(test::Symbol)
	@test (@inferred interp1(1:4, 2:5, 1.5)) == 2.5
	@test (@inferred interp1(1:2, [3im,4im], 1.5)) == 3.5im
	x = LinRange(-1, 1, 21)
	y = sin.(x*π)
	xi = LinRange(-1, 1, 101) * 2
	@test_throws BoundsError interp1(x, y, xi, extrap=nothing)
	@test_throws BoundsError interp1(x, y, xi, extrap=Throw())
	scatter(x, y, label="")
	y0 = @inferred interp1(x, y, xi) # 0 extrapolation
	plot!(xi, y0, label="0", ylim=[-2,2])
	for ex in (Flat(), Line(), Periodic(), Reflect())
		yy = @inferred interp1(x, y, xi, extrap=ex)
		plot!(xi, yy, label="$ex")
	end
	gui()
	true
end


text:
using BAT
using StatsBase
using Distributions
using Test

@testset "BATHist from StatsBase.Hist" begin
    statsbase_hist_1d = fit(Histogram, randn(10_000), closed=:left)
    statsbase_hist_2d = fit(Histogram, (randn(10_000), randn(10_000)), closed=:left)

    bathist_1d = BATHistogram(statsbase_hist_1d)
    bathist_2d = BATHistogram(statsbase_hist_2d)

    @test isa(bathist_1d, BATHistogram)
    @test isa(bathist_2d, BATHistogram)
    @test isa(bathist_2d.h, StatsBase.Histogram)
end

likelihood = let D = 5
    params -> begin
        r = logpdf(MvNormal(zeros(D), ones(D)), params.θ)
        r2 = logpdf(Normal(5, 1.5), params.ϕ)
        LogDVal(r+r2)
    end
end;

prior = BAT.NamedTupleDist(
    θ = MvNormal(zeros(5), ones(5)),
    ϕ = Normal(1, 0.2)
)

posterior = PosteriorDensity(likelihood, prior);
algorithm = MCMCSampling(mcalg = MetropolisHastings(), nchains = 4, nsteps = 10^4)
shaped_samples = bat_sample(posterior, algorithm).result
unshaped_samples = BAT.unshaped.(shaped_samples)

#TODO: what about :θ ? it refers to a multidimensional parameter
@testset "asindex" begin
    @test BAT.asindex(prior, :θ) == 1
    @test BAT.asindex(prior, 1) == 1

    @test BAT.asindex(shaped_samples, :θ) == 1
    @test BAT.asindex(shaped_samples, 1) == 1

    @test_throws ArgumentError BAT.asindex(unshaped_samples, :θ) == 1
    @test BAT.asindex(unshaped_samples, 1) == 1
end

@testset "BATHist from samples" begin
    @test isa(BATHistogram(shaped_samples, :ϕ), BATHistogram)
    @test isa(BATHistogram(shaped_samples, 2), BATHistogram)

    @test_throws ArgumentError BATHistogram(unshaped_samples, :ϕ), BATHistogram
    @test isa(BATHistogram(unshaped_samples, 2), BATHistogram)

    @test isa(BATHistogram(shaped_samples, (:θ, :ϕ)), BATHistogram)
    @test isa(BATHistogram(shaped_samples, (1, 2)), BATHistogram)
end

@testset "BATHist from prior" begin
    @test isa(BATHistogram(prior, :ϕ), BATHistogram)
    @test isa(BATHistogram(prior, 2), BATHistogram)

    @test isa(BATHistogram(prior, (:θ, :ϕ)), BATHistogram)
    @test isa(BATHistogram(prior, (1, 2)), BATHistogram)
end


text:
# Julia version 1.0

function partitions(cards, subtotal)

    m = 0
    # Hit
    for i = 1:10
        if (cards[i]>0)
            total = subtotal+i
            if (total < 21)
                # Stand
                m += 1
                # Hit again
                cards[i] -= 1
                m += partitions(cards, total)
                cards[i] += 1
            elseif (total==21)
                # Stand; hit again is an automatic bust
                m += 1
                break
            end
        end
    end
    return m

end

#deck = Array{Int,10}(4,4,4,4,4,4,4,4,4,16)
deck = Int8[4,4,4,4,4,4,4,4,4,16]
d = 0

for i = 1:10
    # Dealer showing
    deck[i] -= 1
    p = 0
    for j = 1:10
        deck[j] -= 1
        p += partitions(deck, j)
        deck[j] += 1
    end
    println("Dealer showing ", i-1," partitions = ",p)
    global d += p
    deck[i] += 1
end

println("Total partitions = ",d)

max_stars_repo_name: octonion/recursion

max_issues_repo_name: octonion/recursion

max_forks_repo_name: octonion/recursion


text:
import Base: +, -

const addition_types = Dict((:guess, :guess) => :guess,
                            (:exact, :exact) => :exact,
                            (:exact, :lower) => :lower,
                            (:exact, :upper) => :upper,
                            (:lower, :exact) => :lower,
                            (:upper, :exact) => :upper,
                            (:lower, :lower) => :lower,
                            (:upper, :upper) => :upper)

@generated function +(lhs::Sigmoid{N, ES, lhs_mode}, rhs::Sigmoid{N, ES, rhs_mode}) where {N, ES, lhs_mode, rhs_mode}
  if (ES == 0)
    breakdown = :(@breakdown lhs arithmetic; @breakdown rhs arithmetic)
    sub_algorithm = arithmetic_sub
    add_algorithm = arithmetic_add
    build_algorithm = build_arithmetic
  else
    breakdown = :(@breakdown lhs; @breakdown rhs)
    sub_algorithm = numeric_sub
    add_algorithm = numeric_add
    build_algorithm = build_numeric
  end

  #dealing with modes for addition.
  haskey(addition_types, (lhs_mode, rhs_mode)) || throw(ArgumentError("incompatible types passed to addition function!"))
  mode = addition_types[(lhs_mode, rhs_mode)]

  S = Sigmoid{N, ES, mode}
  quotemode = QuoteNode(mode)

  quote
    #correctly set the mode.
    mode = $quotemode

    #adding infinities is infinite.
    if !isfinite(lhs)
        if !isfinite(rhs)
            if (mode == :guess)
                throw(NaNError(+, [$S(Inf), $S(Inf)]))
            elseif (lhs_mode == :exact) || (rhs_mode == :exact)
                return Sigmoid{N,ES,:exact}(Inf)
            else
                return Sigmoid{N,ES,lhs_mode}(Inf)
            end
        end
        return Sigmoid{N,ES,lhs_mode}(Inf)
    end

    isfinite(rhs) || return Sigmoid{N,ES,rhs_mode}(Inf)

    #adding zeros is identity.
    (reinterpret((@UInt), lhs) == zero(@UInt)) && return reinterpret($S, rhs)
    (reinterpret((@UInt), rhs) == zero(@UInt)) && return reinterpret($S, lhs)

    (reinterpret((@UInt), lhs) == -reinterpret((@UInt), rhs)) && return zero($S)

    $breakdown

    if (lhs_sgn != rhs_sgn)

      (dif_sgn, dif_exp, dif_frc) = $sub_algorithm(lhs_sgn, lhs_exp, lhs_frc, rhs_sgn, rhs_exp, rhs_frc)

      __round($build_algorithm($S, dif_sgn, dif_exp, dif_frc))
    else
      (sum_exp, sum_frc) = $add_algorithm(lhs_exp, lhs_frc, rhs_exp, rhs_frc)

      __round($build_algorithm($S, lhs_sgn, sum_exp, sum_frc))
    end
  end
end

-(operand::Sigmoid{N, ES, mode}) where {N, ES, mode} = reinterpret(Sigmoid{N, ES, mode}, -@s(operand))
-(lhs::Sigmoid{N, ES, mode}, rhs::Sigmoid{N, ES, mode}) where {N, ES, mode} = lhs + (-rhs)

function arithmetic_sub(lhs_sgn, lhs_exp, lhs_frc, rhs_sgn, rhs_exp, rhs_frc)
  ## assign top and bottom fractions, ascertain the final sign, and set the
  ## exponential realm.

  if (lhs_exp > rhs_exp)
    top_frc = lhs_frc
    bot_frc = rhs_frc >> (lhs_exp - rhs_exp)
    diff_sgn = lhs_sgn
    diff_exp = lhs_exp
  elseif (rhs_exp > lhs_exp)
    top_frc = rhs_frc
    bot_frc = lhs_frc >> (rhs_exp - lhs_exp)
    diff_sgn = rhs_sgn
    diff_exp = rhs_exp
  elseif (lhs_frc > rhs_frc)
    top_frc = lhs_frc
    bot_frc = rhs_frc
    diff_sgn = lhs_sgn
    diff_exp = lhs_exp
  else
    top_frc = rhs_frc
    bot_frc = lhs_frc
    diff_sgn = rhs_sgn
    diff_exp = rhs_exp
  end

  #subtract the two from each other, with carry check.
  diff_frc = top_frc - bot_frc

  #underflow check.
  if (diff_frc > top_frc)
    #decrease the exponent - note that this can't be less than zero.
    diff_exp -= 1
    diff_frc = (diff_frc << 1)
  else
    #do shifting
    shiftcount = min(leading_zeros(diff_frc), diff_exp)
    diff_exp -= shiftcount
    diff_frc <<= shiftcount
  end

  #return the base sign, the adjusted exponent, and the result fraction.
  (diff_sgn, diff_exp, diff_frc)
end

function numeric_sub(lhs_sgn, lhs_exp, lhs_frc, rhs_sgn, rhs_exp, rhs_frc)
  ## assign top and bottom fractions, ascertain the final sign, and set the
  ## exponential realm.

  if (lhs_exp > rhs_exp)
    top_frc = lhs_frc
    sub_exp_shift = lhs_exp - rhs_exp
    bot_frc = (rhs_frc >> sub_exp_shift) | ((@signbit) >> (sub_exp_shift - 1))
    #set the carry and the guard
    carry_lost = false
    #infer sign and exponent from the top.
    diff_sgn = lhs_sgn
    diff_exp = lhs_exp
  elseif (rhs_exp > lhs_exp)
    top_frc = rhs_frc
    sub_exp_shift = rhs_exp - lhs_exp
    bot_frc = (lhs_frc >> sub_exp_shift) | ((@signbit) >> (sub_exp_shift - 1))
    #set the carry and the guard
    carry_lost = false
    #infer sign and exponent from the top.
    diff_sgn = rhs_sgn
    diff_exp = rhs_exp
  elseif (lhs_frc > rhs_frc)
    top_frc = lhs_frc
    bot_frc = rhs_frc
    #set the carry and the guard
    carry_lost = true
    diff_sgn = lhs_sgn
    diff_exp = lhs_exp
  else
    top_frc = rhs_frc
    bot_frc = lhs_frc
    #set the carry and the guard
    carry_lost = true
    #infer sign and exponent from the top.
    diff_sgn = rhs_sgn
    diff_exp = rhs_exp
  end

  #subtract the two from each other, with carry check.
  diff_frc = top_frc - bot_frc

  #underflow check.
  carry_lost |= (diff_frc > top_frc)

  #reassign the exponent if we've lost the carry bit
  if (carry_lost)
    sub_shift = leading_zeros(diff_frc) + 1
    diff_exp -= sub_shift
    diff_frc <<= sub_shift
  end

  #return the base sign, the adjusted exponent, and the result fraction.
  (diff_sgn, diff_exp, diff_frc)
end

function numeric_add(lhs_exp, lhs_frc, rhs_exp, rhs_frc)

  carry = 1

  if (lhs_exp > rhs_exp)
    add_exp_shift = lhs_exp - rhs_exp
    rhs_frc >>= add_exp_shift
    rhs_frc |= (@signbit) >> (add_exp_shift - 1)
    sum_exp = lhs_exp
    carry = 1
  elseif (lhs_exp < rhs_exp)
    add_exp_shift = rhs_exp - lhs_exp
    lhs_frc >>= add_exp_shift
    lhs_frc |= (@signbit) >> (add_exp_shift - 1)
    sum_exp = rhs_exp
    carry = 1
  else
    sum_exp = lhs_exp
    carry = 2
  end

  #add the two together, with carry check.
  sum_frc = rhs_frc + lhs_frc

  #carry check.
  carry += (sum_frc < rhs_frc)

  if carry > 1
    sum_exp += 1
    #move the fraction over one, and append the top bit.
    sum_frc >>= 1
    (carry == 3) && (sum_frc |= (@signbit))
  end

  (sum_exp, sum_frc)
end

function arithmetic_add(lhs_exp, lhs_frc, rhs_exp, rhs_frc)
  if (lhs_exp > rhs_exp)
    rhs_frc >>= (lhs_exp - rhs_exp)
    sum_exp = lhs_exp
  else
    lhs_frc >>= (rhs_exp - lhs_exp)
    sum_exp = rhs_exp
  end
  #add the two together, with carry check.
  sum_frc = rhs_frc + lhs_frc

  #carry check.
  if (sum_frc < rhs_frc)
    #increase the exponent.
    sum_exp += 1
    #move the fraction over one, and append the top bit.
    sum_frc = (sum_frc >> 1) | @signbit
  end

  (sum_exp, sum_frc)
end

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
@doc raw"""
```
t = find_threshold(Balanced(), histogram, edges)
```

Computes the threshold value using Yen's maximum corrilation criterion for
bileval thresholding.

# Output

Returns a real number `t` that specifies the threshold.

# Details

The algorithm computes the corilation between the two sides of the threshold and
attempts to maximize the corilation. The corrilation is defined as ``C(X)=-\ln{
\sum\limits_{i \ge 0} p_i^2}``.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia

using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img, 256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Yen(), counts[1:end], edges)
```

# Reference

1. Yen JC, Chang FJ, Chang S (1995), "A New Criterion for Automatic Multilevel Thresholding", IEEE Trans. on Image Processing 4 (3): 370-378, [doi:10.1109/83.366472](https://doi.org/10.1109/83.366472)
"""
function find_threshold(algorithm::Yen, histogram::AbstractArray, edges::AbstractRange)
    total = sum(histogram)
    m = length(histogram)
    p = zeros(m)

    # calulate probability
    for i in eachindex(p)
        p[i] = histogram[i]/total
    end

    # setup sums
    Pₛ = 0
    Gₛ = 0
    G′ₛ = sum(p.^2)
    max_TC = typemin(Float64)
    t = 0
    for s in eachindex(p[1:m-1])
        # update sums
        Pₛ += p[s]
        Gₛ += p[s]^2
        G′ₛ -= p[s]^2

        # calculate TC
        TC = -log(Gₛ * G′ₛ) + 2 * log(Pₛ * (1 - Pₛ))
        if TC > max_TC
            max_TC = TC
            t = s
        end
    end

    return edges[t]
end


text:
module Advent

export
    # input
    read_charmatrix,

    # site
    save_input

include("parse.jl")
include("site.jl")

end


text:
module PedCar

#=
Computationally efficient implementation of an autonomous driving scenario modeled as an MDP
This problem involves three agents:
    - an ego vehicle
    - an other car
    - a pedestrian
=#

using StaticArrays
using GridInterpolations
using POMDPs
using POMDPModelTools
using ParticleFilters
using AutomotiveDrivingModels
using AutoViz
using AutomotivePOMDPs
#using AutomotiveSensors
using Parameters
using Random
using LinearAlgebra

export
    PedCarMDP,
    PedCarMDPState,
    PedCarMDPAction,
    init_car_grid,
    init_ped_grid,
    init_l_grid,
    init_v_grid,
    car_transition,
    car_interpolation,
    car_reset,
    ped_transition,
    ped_interpolation,
    ped_reset,
    ego_transition,
    crash,
    ind2state,
    init_transition!,
    init_collision_checker!,
    animate_states,
    state2scene,
    get_mdp_state

include("mdp_type.jl")
include("spaces.jl")
include("interpolation.jl")
include("driver_models_helpers.jl")
include("transition.jl")
include("helpers.jl")
include("reward.jl")
include("render_helpers.jl")

#export
#    init_car_belief,
#    init_ped_belief,
#    ParticleOverlay,
#    CarState,
#    CarStateEstimator,
#    PedStateEstimator

#include("belief_updater.jl")

end # module


text:
function CircularOrbit(arg0::Orbit)
    return CircularOrbit((Orbit,), arg0)
end

function CircularOrbit(arg0::PVCoordinates, arg1::Frame, arg2::AbsoluteDate, arg3::jdouble)
    return CircularOrbit((PVCoordinates, Frame, AbsoluteDate, jdouble), arg0, arg1, arg2, arg3)
end

function CircularOrbit(arg0::TimeStampedPVCoordinates, arg1::Frame, arg2::jdouble)
    return CircularOrbit((TimeStampedPVCoordinates, Frame, jdouble), arg0, arg1, arg2)
end

function CircularOrbit(arg0::jdouble, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::PositionAngle, arg13::Frame, arg14::AbsoluteDate, arg15::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return CircularOrbit((jdouble, jdouble, jdouble, jdouble, PositionAngle, Frame, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function CircularOrbit(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::PositionAngle, arg7::Frame, arg8::AbsoluteDate, arg9::jdouble)
    return CircularOrbit((jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, PositionAngle, Frame, AbsoluteDate, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function add_kepler_contribution(obj::CircularOrbit, arg0::PositionAngle, arg1::jdouble, arg2::Vector{jdouble})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, jdouble, Vector{jdouble}), arg0, arg1, arg2)
end

function eccentric_to_mean(::Type{CircularOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(CircularOrbit, "eccentricToMean", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function eccentric_to_true(::Type{CircularOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(CircularOrbit, "eccentricToTrue", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_a(obj::CircularOrbit)
    return jcall(obj, "getA", jdouble, ())
end

function get_a_dot(obj::CircularOrbit)
    return jcall(obj, "getADot", jdouble, ())
end

function get_alpha(obj::CircularOrbit, arg0::PositionAngle)
    return jcall(obj, "getAlpha", jdouble, (PositionAngle,), arg0)
end

function get_alpha_dot(obj::CircularOrbit, arg0::PositionAngle)
    return jcall(obj, "getAlphaDot", jdouble, (PositionAngle,), arg0)
end

function get_alpha_e(obj::CircularOrbit)
    return jcall(obj, "getAlphaE", jdouble, ())
end

function get_alpha_e_dot(obj::CircularOrbit)
    return jcall(obj, "getAlphaEDot", jdouble, ())
end

function get_alpha_m(obj::CircularOrbit)
    return jcall(obj, "getAlphaM", jdouble, ())
end

function get_alpha_m_dot(obj::CircularOrbit)
    return jcall(obj, "getAlphaMDot", jdouble, ())
end

function get_alpha_v(obj::CircularOrbit)
    return jcall(obj, "getAlphaV", jdouble, ())
end

function get_alpha_v_dot(obj::CircularOrbit)
    return jcall(obj, "getAlphaVDot", jdouble, ())
end

function get_circular_ex(obj::CircularOrbit)
    return jcall(obj, "getCircularEx", jdouble, ())
end

function get_circular_ex_dot(obj::CircularOrbit)
    return jcall(obj, "getCircularExDot", jdouble, ())
end

function get_circular_ey(obj::CircularOrbit)
    return jcall(obj, "getCircularEy", jdouble, ())
end

function get_circular_ey_dot(obj::CircularOrbit)
    return jcall(obj, "getCircularEyDot", jdouble, ())
end

function get_e(obj::CircularOrbit)
    return jcall(obj, "getE", jdouble, ())
end

function get_e_dot(obj::CircularOrbit)
    return jcall(obj, "getEDot", jdouble, ())
end

function get_equinoctial_ex(obj::CircularOrbit)
    return jcall(obj, "getEquinoctialEx", jdouble, ())
end

function get_equinoctial_ex_dot(obj::CircularOrbit)
    return jcall(obj, "getEquinoctialExDot", jdouble, ())
end

function get_equinoctial_ey(obj::CircularOrbit)
    return jcall(obj, "getEquinoctialEy", jdouble, ())
end

function get_equinoctial_ey_dot(obj::CircularOrbit)
    return jcall(obj, "getEquinoctialEyDot", jdouble, ())
end

function get_hx(obj::CircularOrbit)
    return jcall(obj, "getHx", jdouble, ())
end

function get_hx_dot(obj::CircularOrbit)
    return jcall(obj, "getHxDot", jdouble, ())
end

function get_hy(obj::CircularOrbit)
    return jcall(obj, "getHy", jdouble, ())
end

function get_hy_dot(obj::CircularOrbit)
    return jcall(obj, "getHyDot", jdouble, ())
end

function get_i(obj::CircularOrbit)
    return jcall(obj, "getI", jdouble, ())
end

function get_i_dot(obj::CircularOrbit)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_le(obj::CircularOrbit)
    return jcall(obj, "getLE", jdouble, ())
end

function get_le_dot(obj::CircularOrbit)
    return jcall(obj, "getLEDot", jdouble, ())
end

function get_lm(obj::CircularOrbit)
    return jcall(obj, "getLM", jdouble, ())
end

function get_lm_dot(obj::CircularOrbit)
    return jcall(obj, "getLMDot", jdouble, ())
end

function get_lv(obj::CircularOrbit)
    return jcall(obj, "getLv", jdouble, ())
end

function get_lv_dot(obj::CircularOrbit)
    return jcall(obj, "getLvDot", jdouble, ())
end

function get_right_ascension_of_ascending_node(obj::CircularOrbit)
    return jcall(obj, "getRightAscensionOfAscendingNode", jdouble, ())
end

function get_right_ascension_of_ascending_node_dot(obj::CircularOrbit)
    return jcall(obj, "getRightAscensionOfAscendingNodeDot", jdouble, ())
end

function get_type(obj::CircularOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function has_derivatives(obj::Orbit)
    return jcall(obj, "hasDerivatives", jboolean, ())
end

function interpolate(obj::CircularOrbit, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", CircularOrbit, (AbsoluteDate, Stream), arg0, arg1)
end

function mean_to_eccentric(::Type{CircularOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(CircularOrbit, "meanToEccentric", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function shifted_by(obj::CircularOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", CircularOrbit, (jdouble,), arg0)
end

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

function true_to_eccentric(::Type{CircularOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(CircularOrbit, "trueToEccentric", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end


text:
"""
Returns a matrix of marginal damages per one ton of additional emissions of the specified gas in the specified year.
"""
function getmarginaldamages(; emissionyear=2010, parameters = nothing, yearstoaggregate = 1000, gas = :C)

    # Calculate number of years to run the models
    yearstorun = min(1050, getindexfromyear(emissionyear) + yearstoaggregate)

    # Get default FUND model
    m1 = MimiFUND.get_model(nsteps = yearstorun, params = parameters)

    # Get model to add marginal emissions to
    m2 = MimiFUND.get_model(nsteps = yearstorun, params = parameters)
    add_comp!(m2, Mimi.adder, :marginalemission, before = :climateco2cycle, first = 1951)
    addem = zeros(yearstorun)
    addem[getindexfromyear(emissionyear)-1:getindexfromyear(emissionyear) + 8] .= 1.0
    set_param!(m2, :marginalemission, :add, addem)

    # Reconnect the appropriate emissions in the marginal model
    if gas == :C
        connect_param!(m2, :marginalemission, :input, :emissions, :mco2)
        connect_param!(m2, :climateco2cycle, :mco2, :marginalemission, :output, repeat([missing], yearstorun + 1))
    elseif gas == :CH4
        connect_param!(m2, :marginalemission, :input, :emissions, :globch4)
        connect_param!(m2, :climatech4cycle, :globch4, :marginalemission, :output, repeat([missing], yearstorun + 1))
    elseif gas == :N2O
        connect_param!(m2, :marginalemission, :input, :emissions, :globn2o)
        connect_param!(m2, :climaten2ocycle, :globn2o, :marginalemission, :output, repeat([missing], yearstorun + 1))
    elseif gas == :SF6
        connect_param!(m2, :marginalemission, :input, :emissions,:globsf6)
        connect_param!(m2, :climatesf6cycle, :globsf6, :marginalemission, :output, repeat([missing], yearstorun + 1))
    else
        error("Unknown gas.")
    end

    # Run both models
    run(m1; ntimesteps = yearstorun)
    run(m2; ntimesteps = yearstorun)

    # Get damages
    damage1 = m1[:impactaggregation, :loss]
    damage2 = m2[:impactaggregation, :loss]

    # Calculate the marginal damage between run 1 and 2 for each year/region
    marginaldamages = (damage2 .- damage1) / 10000000.0

    return marginaldamages
end


text:
__precompile__()

module Ev3dev

import Base: read, parse, write, call

immutable Brick
    root_path::AbstractString

    Brick(root_path::AbstractString="/") = new(root_path)
end

include("attributes.jl")
include("classes.jl")
include("devices.jl")
include("higher_level_commands.jl")

export Brick,
       LargeMotor,
       MediumMotor,
       UltrasoundSensor,
       ColorSensor,
       TouchSensor,
       GyroSensor,
       find_devices,
       find_device_at_address,
       scaled_values,
       command,
       run_continuous,
       servo_absolute,
       servo_relative,
       stop,
       Degrees,
       Radians

end


text:
function estimation_step(obj::KalmanEstimator, arg0::ObservedMeasurement)
    return jcall(obj, "estimationStep", Vector{AbstractIntegratedPropagator}, (ObservedMeasurement,), arg0)
end

function get_current_date(obj::KalmanEstimator)
    return jcall(obj, "getCurrentDate", AbsoluteDate, ())
end

function get_current_measurement_number(obj::KalmanEstimator)
    return jcall(obj, "getCurrentMeasurementNumber", jint, ())
end

function get_estimated_measurements_parameters(obj::KalmanEstimator)
    return jcall(obj, "getEstimatedMeasurementsParameters", ParameterDriversList, ())
end

function get_orbital_parameters_drivers(obj::KalmanEstimator, arg0::jboolean)
    return jcall(obj, "getOrbitalParametersDrivers", ParameterDriversList, (jboolean,), arg0)
end

function get_physical_estimated_covariance_matrix(obj::KalmanEstimator)
    return jcall(obj, "getPhysicalEstimatedCovarianceMatrix", RealMatrix, ())
end

function get_physical_estimated_state(obj::KalmanEstimator)
    return jcall(obj, "getPhysicalEstimatedState", RealVector, ())
end

function get_propagation_parameters_drivers(obj::KalmanEstimator, arg0::jboolean)
    return jcall(obj, "getPropagationParametersDrivers", ParameterDriversList, (jboolean,), arg0)
end

function process_measurements(obj::KalmanEstimator, arg0::Iterable)
    return jcall(obj, "processMeasurements", Vector{AbstractIntegratedPropagator}, (Iterable,), arg0)
end

function set_observer(obj::KalmanEstimator, arg0::KalmanObserver)
    return jcall(obj, "setObserver", void, (KalmanObserver,), arg0)
end


text:
"""
	isrooted(tree)

Return `true` if the `tree` is rooted.
"""
isrooted(tree::AbstractTree) = tree.rooted

"""
	root!(tree)

Set a `tree` as rooted. The origin node of the tree is set as the root.
"""
function root!(tree::AbstractTree)
	tree.rooted = true

	return nothing
end

"""
	root!(tree, node)

Set `node` as the root of a `tree`.
"""
function root!(tree::AbstractTree, p::AbstractNode)
	tree.origin = p
	root!(tree)

	return nothing
end

"""
	unroot!(tree)

Set a `tree` as unrooted.
"""
function unroot!(tree::AbstractTree)
	tree.rooted = false

	return nothing
end

"""
	origin(tree)

Return the "origin" of a `tree`: The node that serves as the default point of entry to the tree.
"""
origin(tree::AbstractTree) = tree.origin

"""
	origin!(tree, node)

Set `node` as the "origin" of `tree`. The origin of the tree serves as a "pseudoroot": the default point of entry to the tree.
"""
function origin!(tree::AbstractTree, p::AbstractNode)
	tree.origin = p
end

"""
	label(tree)

Get the label of a `tree`.
"""
label(tree::AbstractTree) = tree.label

"""
	label(tree, v)

Set the string `v` as the label of a `tree`.
"""
function label!(tree::AbstractTree, v::String)
	tree.label = v

	return nothing
end

"""
	annotation(tree, k)

Get the value of the annotation with key `k` in `tree`.
"""
annotation(tree::AbstractTree, k::AbstractString) = tree.annotations[k]

"""
	annotate!(tree, key, value)

Set the `value` for an annotation `key` in `tree`.

Annotation keys can only be `AbstractString`s.
"""
function annotate!(tree::AbstractTree, k::AbstractString, v::Any)
	tree.annotations[k] = v

	return nothing
end

"""
	deannotate!(tree, key)

Remove an annotation `key` from a `tree`.
"""
function deannotate!(tree::AbstractTree, k::AbstractString)
	delete!(tree.annotations, k)

	return nothing
end

function _setspecies!(p::AbstractNode, q::AbstractNode, dir::SpeciesDirectory)
	if label(q) ∈ dir
		q.species = dir[label(q)]
	end
	for link in q.links
		link.to == p && continue
		_setspecies!(q, link.to, dir)
	end
end

"""
	setspecies!(tree [, species_directory])

Set the species in `tree` by matching node labels to species names in the species directory of the tree.
"""
setspecies!(tree::AbstractTree) =
	_setspecies!(tree.origin, tree.origin, tree.dir)


text:
using IntervalMatrices, Test, LinearAlgebra, SparseArrays

using IntervalMatrices: _truncated_exponential_series,
                        horner,
                        scale_and_square,
                        correction_hull,
                        input_correction
include("models.jl")

include("constructors.jl")
include("arithmetic.jl")
include("setops.jl")
include("exponential.jl")
include("affine.jl")


text:
module BenchWords
using BenchmarkTools
include("../../examples/words.jl")

suite = BenchmarkGroup()

function makewords(n)
    genpos() = 2 ^ rand(0:5) + rand(0:10)
    io = IOBuffer()
    for _ in 1:n
        for _ in 1:genpos()
            write(io, ' ')
        end
        write(io, String(rand('a':'z', genpos())))
    end
    return String(take!(io))
end

str = makewords(2^14)

for n in [1, 2, 4]
    b = length(str) ÷ n
    suite["nthreads=$n"] = @benchmarkable countwords($str; basesize=$b)
end

end  # module
BenchWords.suite


text:
# Automatically generated test file
function super_exponential1()

g = [
1.000000000000000,
1.064494458917860,
1.284025416687741,
1.755054656960298,
2.718281828459045,
4.770733181967603,
9.487735836358526,
21.380942759123343,
54.598150033144236,
157.984985495187459,
518.012824668341977,
1924.651131676947216,
8103.083927575384223]

ζ_vec = [2.09751e6, 2.19867e5, 23047.2, 2415.87, 253.239, 26.5452, 2.78255, 0.291675, 0.0305742, 0.00320488]

I_sols = [
[1, 13],
[1, 12, 13],
[1, 11, 12, 13],
[1, 10, 11, 12, 13],
[1, 8, 10, 11, 12, 13],
[1, 7, 9, 10, 11, 12, 13],
[1, 6, 8, 9, 10, 11, 12, 13],
[1, 5, 7, 8, 9, 10, 11, 12, 13],
[1, 4, 6, 7, 8, 9, 10, 11, 12, 13],
[1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13],
[1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]]

f_sols = [
40316842.674393124878407,
2097512.634645521640778,
136141.474100708961487,
10875.019243597984314,
837.684379041194916,
79.105205863714218,
9.202944725751877,
1.286040574312210,
0.211227089166641,
0.041290149092674,
0.004006102681160,
0.000000044703484]

return g, ζ_vec, I_sols, f_sols

end


text:
module GrammaticalEvolutions

using ExprRules
using StatsBase

using ExprOptimization: ExprOptAlgorithm, ExprOptResult, BoundedPriorityQueue, enqueue!
import ExprOptimization: optimize

export GrammaticalEvolution

const OPERATORS = [:reproduction, :crossover, :mutation]

abstract type SelectionMethod end
abstract type MutationMethod end
abstract type TrackingMethod end

"""
    GrammaticalEvolution

Grammatical Evolution.
# Arguments
- `grammar::Grammar`: grammar
- `typ::Symbol`: start symbol
- `pop_size::Int`: population size
- `iterations::Int`: number of iterations
- `init_gene_length::Int`: initial length of genotype integer array
- `max_gene_length::Int`: maximum length of genotype integer array
- `max_depth::Int`: maximum depth of derivation tree
- `p_reproduction::Float64`: probability of reproduction operator
- `p_crossover::Float64`: probability of crossover operator
- `p_mutation::Float64`: probability of mutation operator
- `select_method::SelectionMethod`: selection method (default: tournament selection)
- `mutate_method::InitializationMethod`: mutation method (default: multi-mutate)
- `track_method::TrackingMethod`: additional tracking, e.g., track top k exprs (default: no additional tracking)
"""
struct GrammaticalEvolution <: ExprOptAlgorithm
    pop_size::Int
    iterations::Int
    init_gene_length::Int
    max_gene_length::Int
    max_depth::Int
    p_operators::Weights
    select_method::SelectionMethod
    mutate_method::MutationMethod
    track_method::TrackingMethod

    function GrammaticalEvolution(
        grammar::Grammar,
        typ::Symbol,
        pop_size::Int,                          #population size
        iterations::Int,                        #number of generations
        init_gene_length::Int,                  #initial length of genotype Int vector
        max_gene_length::Int,                   #maximum length of genotype Int vector
        max_depth::Int,                         #maximum depth of derivation tree
        p_reproduction::Float64,                #probability of reproduction operator
        p_crossover::Float64,                   #probability of crossover operator
        p_mutation::Float64;                    #probability of mutation operator
        select_method::SelectionMethod=TournamentSelection(),   #selection method
        mutate_method::MutationMethod=MultiMutate(grammar, typ),
        track_method::TrackingMethod=NoTracking())   #tracking method

        p_operators = Weights([p_reproduction, p_crossover, p_mutation])
        new(pop_size, iterations, init_gene_length, max_gene_length, max_depth, p_operators,
            select_method, mutate_method, track_method)
    end
end

"""
    TournamentSelection

Tournament selection method with tournament size k.
"""
struct TournamentSelection <: SelectionMethod
    k::Int
end
TournamentSelection() = TournamentSelection(2)

"""
    TruncationSelection

Truncation selection method keeping the top k individuals
"""
struct TruncationSelection <: SelectionMethod
    k::Int
end
TruncationSelection() = TruncationSelection(100)

"""
    RandomMutation

Randomly change each entry in integer array with probability p_mutate
"""
struct RandomMutation <: MutationMethod
    p_mutate::Float64
end
RandomMutation() = RandomMutation(0.1)

"""
    GeneDuplication

Pick a random segment of the gene and duplicate it at the end.
"""
struct GeneDuplication <: MutationMethod end

"""
    GenePruning

With probability p_prune, decode the gene using grammar and start symbol typ and discard the genes not used.
"""
struct GenePruning <: MutationMethod
    p_prune::Float64
    grammar::Grammar
    typ::Symbol
end
GenePruning(grammar::Grammar, typ::Symbol) = GenePruning(0.1, grammar, typ)

"""
    MultiMutate

Apply a sequence of mutation operators.
"""
struct MultiMutate <: MutationMethod
    Ms::Vector{MutationMethod}
end
MultiMutate(grammar::Grammar, typ::Symbol) = MultiMutate([
    RandomMutation(),
    GeneDuplication(),
    GenePruning(grammar, typ)])

"""
    NoTracking

No additional tracking of expressions.
"""
struct NoTracking <: TrackingMethod end

"""
    TopKTracking

Track the top k expressions.
"""
struct TopKTracking <: TrackingMethod
    k::Int
    q::BoundedPriorityQueue{RuleNode,Float64}

    function TopKTracking(k::Int)
        q = BoundedPriorityQueue{RuleNode,Float64}(k,Base.Order.Reverse) #lower is better
        obj = new(k, q)
        obj
    end
end

"""
    optimize(p::GrammaticalEvolution, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)

Grammatical Evolution algorithm with parameters p, grammar 'grammar', start symbol typ, and loss function 'loss'.  Loss function has the form: los::Float64=loss(node::RuleNode, grammar::Grammar).

See: Ryan, Collins, O'Neil, "Grammatical Evolution: Evolving Programs for an Arbitrary Language",
    in European Conference on Genetic Programming, Spring, 1998, pp. 83-96.
"""
function optimize(p::GrammaticalEvolution, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)
    grammatical_evolution(p, grammar, typ, loss; kwargs...)
end

"""
    grammatical_evolution(p::GrammaticalEvolution, grammar::Grammar, typ::Symbol, loss::Function)

Grammatical Evolution algorithm with parameters p, grammar 'grammar', start symbol typ, and loss function 'loss'.  Loss function has the form los::Float64=loss(node::RuleNode, grammar::Grammar).

See: Ryan, Collins, O'Neil, "Grammatical Evolution: Evolving Programs for an Arbitrary Language",
    in European Conference on Genetic Programming, Spring, 1998, pp. 83-96.
"""
function grammatical_evolution(p::GrammaticalEvolution, grammar::Grammar, typ::Symbol,
    loss::Function)
    iseval(grammar) && error("Grammatical Evolution does not support _() functions in the grammar")

    pop0 = initialize(p.pop_size, p.init_gene_length)
    pop1 = [Int[] for i=1:p.pop_size]
    losses0 = Vector{Union{Float64,Missing}}(missing,p.pop_size)
    losses1 = Vector{Union{Float64,Missing}}(missing,p.pop_size)

    best_tree, best_loss = evaluate!(p, grammar, typ, loss, pop0, losses0, RuleNode(0), Inf)
    for iter = 1:p.iterations
        fill!(losses1, missing)
        i = 0
        while i < p.pop_size
            op = sample(OPERATORS, p.p_operators)
            if op == :reproduction
                ind1,j = select(p.select_method, pop0, losses0)
                pop1[i+=1] = ind1
                losses1[i] = losses0[j]
            elseif op == :crossover
                ind1,_ = select(p.select_method, pop0, losses0)
                ind2,_ = select(p.select_method, pop0, losses0)
                child = crossover(ind1, ind2)
                limit_length!(child, p.max_gene_length)
                pop1[i+=1] = child
            elseif op == :mutation
                ind1,_ = select(p.select_method, pop0, losses0)
                child1 = mutation(p.mutate_method, ind1)
                limit_length!(child1, p.max_gene_length)
                pop1[i+=1] = child1
            end
        end
        pop0, pop1 = pop1, pop0
        losses0, losses1 = losses1, losses0
        best_tree, best_loss = evaluate!(p, grammar, typ, loss, pop0, losses0, best_tree, best_loss)
    end
    alg_result = Dict{Symbol,Any}()
    _add_result!(alg_result, p.track_method)
    ExprOptResult(best_tree, best_loss, get_executable(best_tree, grammar), alg_result)
end

"""
    _add_result!(d::Dict{Symbol,Any}, t::NoTracking)

Add tracking results to alg_result.  No op for NoTracking.
"""
_add_result!(d::Dict{Symbol,Any}, t::NoTracking) = nothing
"""
    _add_result!(d::Dict{Symbol,Any}, t::TopKTracking)

Add tracking results to alg_result.
"""
function _add_result!(d::Dict{Symbol,Any}, t::TopKTracking)
    d[:top_k] = collect(t.q)
    d
end

"""
    initialize(pop_size::Int, len::Int)

Randomly initialize population.
"""
initialize(pop_size::Int, len::Int) = [rand(Int, len) for i = 1:pop_size]

"""
    limit_length!(ind::Vector{Int}, max_length::Int)

Limit the length of the gene ind to max_length.
"""
function limit_length!(ind::Vector{Int}, max_length::Int)
    resize!(ind, min(length(ind), max_length))
    ind
end

"""
    select(p::TournamentSelection, pop::Vector{RuleNode}, losses::Vector{Union{Float64,Missing}})

Tournament selection.
"""
function select(p::TournamentSelection, pop::Vector{Vector{Int}}, losses::Vector{Union{Float64,Missing}})
    ids = StatsBase.sample(1:length(pop), p.k; replace=false, ordered=true)
    i = ids[1] #assumes pop is sorted
    pop[i], i
end

"""
    select(p::TruncationSelection, pop::Vector{RuleNode}, losses::Vector{Union{Float64,Missing}})

Truncation selection.
"""
function select(p::TruncationSelection, pop::Vector{Vector{Int}}, losses::Vector{Union{Float64,Missing}})
    i = rand(1:p.k)  #assumes pop is sorted
    pop[i], i
end

"""
    crossover(a::Vector{Int}, b::Vector{Int})

Crossover genetic operator.  Pick a random crossover point from 'a', then pick a crossover point from 'b', splice first part of 'a' with second part of 'b'
"""
function crossover(a::Vector{Int}, b::Vector{Int})
    i = rand(1:length(a))
    return vcat(a[1:i], b[i+1:end])
end

"""
    mutation(::RandomMutation, a::Vector{Int})

Mutation the gene using RandomMutation.  Randomly change each entry in integer array with probability p_mutate
"""
function mutation(p::RandomMutation, a::Vector{Int})
    child = deepcopy(a)
    for i in eachindex(child)
        if rand() < p.p_mutate
            child[i] = rand(Int)
        end
    end
    return child
end

"""
    mutation(p::GeneDuplication, a::Vector{Int})

Mutate the gene using GeneDuplication.  Pick a random segment of the gene and duplicate it at the end.
"""
function mutation(p::GeneDuplication, a::Vector{Int})
    child = a
    n = length(a)
    i, j = rand(1:n), rand(1:n)
    interval = min(i,j) : max(i,j)
    return vcat(child, child[interval])
end

"""
    mutation(p::GenePruning, a::Vector{Int})

Mutate the gene using GenePruning.  With probability p_prune, decodes the gene and discards the genes not used.
"""
function mutation(p::GenePruning, a::Vector{Int})
    child = deepcopy(a)
    if rand() < p.p_prune
        c = decode(child, p.grammar, p.typ).n_rules_applied
        if c < length(a)
            child = child[1:c]
        end
    end
    return child
end

"""
    mutation(p::MultiMutate, a::Vector{Int})

Mutate the gene using MultiMutate.  Apply a sequence of mutation operators.
"""
function mutation(p::MultiMutate, child::Vector{Int})
    for m in p.Ms
        child = mutation(m, child)
    end
    return child
end

"""
    evaluate!(p::GrammaticalEvolution, grammar::Grammar, typ::Symbol, loss::Function, pop::Vector{RuleNode}, losses::Vector{Union{Float64}}, best_tree::RuleNode, best_loss::Float64)

Evaluate the loss function for population and sort.  Update the globally best tree, if needed.
"""
function evaluate!(p::GrammaticalEvolution, grammar::Grammar, typ::Symbol, loss::Function,
                   pop::Vector{Vector{Int}}, losses::Vector{Union{Float64,Missing}},
                   best_tree::RuleNode, best_loss::Float64)

    # Pre-compute indics that are missing to make
    # the loop multi-threading friendly.
    idcs_missing = filter(i -> ismissing(losses[i]), eachindex(pop))
    Threads.@threads for i in idcs_missing
        decoded = decode(pop[i], grammar, typ)
        losses[i] = depth(decoded.node) > p.max_depth ?  Inf : loss(decoded.node, grammar)
    end

    perm = sortperm(losses)
    pop[:], losses[:] = pop[perm], losses[perm]
    if losses[1] < best_loss
        best_loss = losses[1]
        best_tree = decode(pop[1], grammar, typ).node
    end
    _update_tracker!(p.track_method, pop, losses, grammar, typ)

    (best_tree, best_loss)
end

"""
    _update_tracker!(t::NoTracking, pop::Vector{Vector{Int}},
        losses::Vector{Union{Float64,Missing}})

Update the tracker.  No op for NoTracking.
"""
function _update_tracker!(t::NoTracking, pop::Vector{Vector{Int}},
    losses::Vector{Union{Float64,Missing}}, grammar::Grammar, typ::Symbol)
    nothing
end
"""
    _update_tracker!(t::TopKTracking, pop::Vector{Vector{Int}},
        losses::Vector{Union{Float64,Missing}})

Update the tracker.  Track top k expressions.
"""
function _update_tracker!(t::TopKTracking, pop::Vector{Vector{Int}},
    losses::Vector{Union{Float64,Missing}}, grammar::Grammar, typ::Symbol)
    n = 0
    for i = 1:length(pop)
        tree = decode(pop[i], grammar, typ).node
        r = enqueue!(t.q, tree, losses[i])
        r >= 0 && (n += 1) #no clash
        n >= t.k && break
    end
end

"""
    DecodedExpression

Results of a decode operation.  The decoded expression tree is stored in node and the number of rule is in n_rule_applied.
"""
struct DecodedExpression
    node::RuleNode
    n_rules_applied::Int
end

"""
    decode(x::Vector{Int}, grammar::Grammar, typ::Symbol, c_max=1000, c=0)

Decode an integer array (genotype) to a derivation tree (phenotype) using given grammar and start symbol typ.  Unlimited wraps.
"""
function decode(x::Vector{Int}, grammar::Grammar, typ::Symbol, c_max=1000, c=0)
    node, c = _decode(x, grammar, typ, c_max, c)
    DecodedExpression(node, c)
end
function _decode(x::Vector{Int}, grammar::Grammar, typ::Symbol, c_max=1000, c=0)
    types = grammar[typ]
    if length(types) > 1
        g = x[mod1(c+=1, length(x))]
        rule = types[mod1(g, length(types))]
    else
        rule = types[1]
    end
    node = RuleNode(rule)
    childtypes = child_types(grammar, node)
    if !isempty(childtypes) && c < c_max
        for ctyp in childtypes
            cnode, c = _decode(x, grammar, ctyp, c_max, c)
            push!(node.children, cnode)
        end
    end
    return (node, c)
end

end #module


text:
NORM_TYPE = Inf

function lb_s_thres(iter::Class_iterate,dir::Class_point,pars::Class_parameters)
  ex = pars.ls.fraction_to_boundary_predict_exp
  #alpha_P_tilde = simple_max_step(iter.point.s,dir.s,zeros(length(dir.s)))
  #x_thres = max.(2.0 * LinearAlgebra.norm(dir.x,NORM_TYPE)^(2 + ex), LinearAlgebra.norm(dir.x,NORM_TYPE)^(1.0 + ex) + LinearAlgebra.norm(dir.y,NORM_TYPE)^0.5 * LinearAlgebra.norm(dir.x,NORM_TYPE))
  x_thres = LinearAlgebra.norm(dir.x,NORM_TYPE)^ex #+ LinearAlgebra.norm(dir.y,NORM_TYPE) + get_delta(iter)

  return min.(iter.point.s, LinearAlgebra.norm(dir.x,NORM_TYPE) * x_thres)
end

function lb_s_predict(iter::Class_iterate,dir::Class_point,pars::Class_parameters)
  #@show iter.frac_bd_predict
  return iter.frac_bd_predict .* lb_s_thres(iter,dir,pars)
end

function lb_y(iter::Class_iterate,dir::Class_point,pars::Class_parameters)
  #alpha_P_tilde = simple_max_step(iter.point.s,dir.s,zeros(length(dir.s)))
  return iter.frac_bd .*  iter.point.y * min(1.0, LinearAlgebra.norm(dir.x,NORM_TYPE) )
end

function lb_s(iter::Class_iterate,dir::Class_point,pars::Class_parameters)
    #alpha_P_tilde = simple_max_step(iter.point.s,dir.s,zeros(length(dir.s)))
    lb_s = iter.frac_bd .* lb_s_thres(iter,dir,pars)
    #min.(iter.point.s, LinearAlgebra.norm(dir.x,NORM_TYPE)^2)
    #@show lb_s
    return lb_s
end

function simple_max_step(val::Array{Float64,1}, dir::Array{Float64,1}, lb::Array{Float64,1})
    ratio = maximum( [1.0; -dir ./ (val - lb) ] ) #; -q * dir.y ./ iter.point.y] );
    max_step = 1.0 / ratio
    #@show max_step
    return max_step
end


text:
module DVIDArrays

include("type.jl")
include("backends.jl")

end # end of DVIDArrays module


text:
module DatasetModule

import ..ProgramConstantsModule: BATCH_DIM, FEATURE_DIM

struct Dataset{T<:Real}
    X::AbstractMatrix{T}
    y::AbstractVector{T}
    n::Int
    nfeatures::Int
    weighted::Bool
    weights::Union{AbstractVector{T},Nothing}
    varMap::Array{String,1}
end

"""
    Dataset(X::AbstractMatrix{T}, y::AbstractVector{T};
            weights::Union{AbstractVector{T}, Nothing}=nothing,
            varMap::Union{Array{String, 1}, Nothing}=nothing)

Construct a dataset to pass between internal functions.
"""
function Dataset(
    X::AbstractMatrix{T},
    y::AbstractVector{T};
    weights::Union{AbstractVector{T},Nothing}=nothing,
    varMap::Union{Array{String,1},Nothing}=nothing,
) where {T<:Real}
    Base.require_one_based_indexing(X, y)
    n = size(X, BATCH_DIM)
    nfeatures = size(X, FEATURE_DIM)
    weighted = weights !== nothing
    if varMap === nothing
        varMap = ["x$(i)" for i in 1:nfeatures]
    end

    return Dataset{T}(X, y, n, nfeatures, weighted, weights, varMap)
end

end


text:
function convert2weightedGraph(
    network::AbstractGraph,
	agents::AbstractArray
    )

	src_network = deepcopy(network)

	weighted_net = SimpleWeightedGraph(nv(src_network))

	for v in vertices(src_network)
		for neighbor in neighbors(src_network, v)
			weighted_distance = sum(agents[v].cultureVector .- agents[neighbor].cultureVector)
			add_edge!(weighted_net, v, neighbor, weighted_distance)
			rem_edge!(src_network, v, neighbor)
		end
	end

	return weighted_net
end

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
# Module that wraps running an experiment into a nice package
module CautiousEngine

using Logging
using LoggingExtras
using BSON
using Dates
using Distributed
using Base.Threads
using DataFrames
using LinearAlgebra
using Base.Threads
using GaussianProcesses
using Random
using Distributions
using StatsBase
using Serialization
using MAT
using Printf
using IntervalSets
using SparseArrays
using Plots
using Spot

# Lazy sets, the hot new way to work with polytopes
using LazySets

include("metadata.jl")
include("parameters.jl")
include("validation.jl")
include("e2e.jl")

# My own packages
include("GPBounding/GPBounding.jl")
using .GPBounding

include("abstraction/discretization.jl")
include("abstraction/images.jl")
include("abstraction/intersection.jl")
include("abstraction/regions.jl")
include("abstraction/refinement.jl")

include("imdp/dfa.jl")
include("imdp/imdp.jl")
include("imdp/multimodal.jl")
include("imdp/pimdp.jl")
include("imdp/transitions.jl")

include("online/online_control.jl")
include("online/run_simulation.jl")

include("regression/data.jl")
include("regression/gps.jl")
include("regression/error_bounds.jl")

include("synthesis/verification.jl")
include("synthesis/synthesis.jl")

include("visualization/visualize.jl")

export SystemParameters,
       DataParameters,
       ExperimentParameters

export general_label_fcn

export perform_synthesis_from_result_dirs,
       generate_transition_bounds,
       generate_linear_truth,
       verify_experiment,
       generate_result_dir_name

## TODO: Where tf to put this
"""
Wrapper function for unpacking the experiment parameters.
"""
function generate_training_data(params::ExperimentParameters)

    x_train, y_train = generate_training_data(params.system_params.unknown_dynamics_fcn, params.domain, params.data_params.data_num,
                                              known_fn=params.system_params.known_dynamics_fcn, random_seed=params.random_seed,
                                              n_dims_out=params.system_params.n_dims_out, process_dist=params.system_params.process_noise_dist,
                                              measurement_dist=params.system_params.measurement_noise_dist)

    return x_train, y_train
end

end


text:
"""
    function datatype(str::String)
This function evaluates an input string as a DataType if it is defined.
Otherwise, it will return false.
See: (thread on discourse.julialang.org)[https://discourse.julialang.org/t/parse-string-to-datatype/7118/9]
"""
function datatype(str::String)
    type = :($(Symbol(str)))
    return isdefined(SLiDE, type) ? eval(type) : nothing
end

"""
    Base.broadcastable(x::InvertedIndex{T}) where {T<:Any}
(!!!!) add docs here.
"""
Base.broadcastable(x::InvertedIndex{T}) where {T <: Any} = [x];

"""
    Base.split(x::Missing)
    Base.split(x::Number)
Extends `split` to ignore missing fields.
"""
Base.split(str::Missing) = str
Base.split(str::Missing, splitter::Any) = str
Base.split(x::Symbol) = Symbol.(split(string(x)))

"""
    Base.strip(x::Missing)
    Base.strip(x::Number)
Extends `strip` to ignore missing fields and numbers.
"""
Base.strip(x::String) = replace(x, r"^\s*\"*|\"*\s*$" => "")
Base.strip(x::Missing) = x
Base.strip(x::Number) = x

"""
    Base.lowercase(x::Symbol)
Extends `lowercase` to handle other data types.
"""
Base.lowercase(x::Int) = x
Base.lowercase(x::Symbol) = Symbol(lowercase(string(x)))
Base.lowercase(x::Missing) = missing

"""
    Base.titlecase(x::Symbol)
Extends `titlecase` to handle other data types.
"""
Base.titlecase(x::Int) = x
Base.titlecase(x::Symbol) = Symbol(titlecase(string(x)))
Base.titlecase(x::Missing) = missing

"""
    Base.uppercase(x::Symbol)
Extends `uppercase` to handle other data types.
"""
Base.uppercase(x::Int) = x
Base.uppercase(x::Symbol) = Symbol(uppercase(string(x)))
Base.uppercase(x::Missing) = missing

"""
    Base.uppercasefirst(x::Symbol)
Extends `uppercasefirst` to handle other data types.
"""
Base.uppercasefirst(x::Int) = x
Base.uppercasefirst(x::Symbol) = Symbol(uppercasefirst(string(x)))
Base.uppercasefirst(x::Missing) = missing

"""
    Base.occursin(x::Symbol, y::Symbol)
    Base.occursin(x::String, y::Symbol)
Extends `occursin` to work for symbols. Potentially helpful for DataFrame columns.
"""
Base.occursin(x::Symbol, y::Symbol) = occursin(string(x), y)
Base.occursin(x::String, y::Symbol) = occursin(x, string(y))

"""
    convert_type(::Type{T}, x::Any)
    convert_type(::Dict{Any,Any}, df::DataFrame, value_col::Symbol; kwargs...)
Converts `x` into the specified `Type{xT}`.

Consider extending [convert](https://docs.julialang.org/en/v1/base/base/#Base.convert)
function (!!!!)

# Arguments
- `::Type{T}`: target DataType.
- `x<:Any`: value to convert.

# Keywords
Options available when converting a DataFrame into a dictionary of keys pointing to a value:
- `drop_cols = []`: Columns not to include in
    the keys. By default, no columns are dropped.
- `value_col::Symbol = :end`: If converting

# Returns
Data in specified type
"""
convert_type(::Type{T}, x::Any) where T <: AbstractString = string(x)
convert_type(::Type{T}, x::Dates.Date) where T <: Integer = Dates.year(x)

convert_type(::Type{Map}, x::Group) = Map(x.file, [x.from], x.to, [x.input], x.output, :inner)
convert_type(::Type{CSVInput}, x::DataInput) = CSVInput(x.name, x.descriptor)
convert_type(::Type{CSVInput}, x::SetInput) = CSVInput(x.name, "set")

convert_type(::Type{T}, x::AbstractString) where T <: AbstractString = string(strip(x))

function convert_type(::Type{T}, x::AbstractString) where T <: Integer
    return convert_type(T, convert_type(Float64, x))
end

function convert_type(::Type{T}, x::AbstractString) where T <: Real
    return parse(T, reduce(replace, ["," => "", "\"" => ""], init=x))
end

convert_type(::Type{T}, x::Symbol) where T <: Real = convert_type(T, convert_type(String, x))

function convert_type(::Type{DataFrame}, x::Array{Any,2})
    col = Symbol.(replace(x[1,:], "" => missing))
    data = x[2:end,:]
    return DataFrame(data, col, makeunique = true)
end

convert_type(::Type{DataFrame}, lst::Array{Dict{Any,Any},1}) = [DataFrame.(lst)...;]

function convert_type(::Type{DataFrame}, arr::JuMP.Containers.DenseAxisArray; cols=[])
    cols = ensurearray(cols)

    val = JuMP.value.(arr.data)
    ind = permute(arr.axes);
    val = collect(Iterators.flatten(val));

    df = hcat(DataFrame(ensuretuple.(ind)), DataFrame([val], [:value]))
    return edit_with(df, Rename.(propertynames(df)[1:length(cols)], cols))
end

function convert_type(::Type{Dict}, df::DataFrame; drop_cols=[], value_col::Symbol=:Float)
    # Find and save the column containing values and that/those containing keys.
    # If no value column indicator is specified, find the first DataFrame column of floats.
    value_col == :Float && (value_col = find_oftype(df, AbstractFloat)[1])
    key_cols = setdiff(propertynames(df), convert_type.(Symbol, ensurearray(drop_cols)), [value_col])
    ONEKEY = length(key_cols) == 1

    d = Dict((ONEKEY ? row[key_cols[1]] : (row[key_cols]...,)) => row[value_col]
        for row in eachrow(df))
    return d
end

convert_type(::Type{DataType}, x::AbstractString) = datatype(x)

convert_type(::Type{Array{T}}, x::Any) where T <: Any = convert_type.(T, x)
convert_type(::Type{Array{T,1}}, x::Any) where T <: Any = convert_type.(T, x)

# !!!! WARNING: DEPRECIATED.
convert_type(::Type{Array}, d::Dict) = [collect(values(d))...;]

convert_type(::Type{Array{Tuple}}, df::DataFrame) = ensuretuple.(eachrow(df))
convert_type(::Type{T}, x::Missing) where T <: Real = x;
convert_type(::Type{T}, x::Missing) where T <: AbstractString = x
convert_type(::Type{Any}, x::AbstractString) = "missing" == lowercase(x) ? missing : x
convert_type(::Type{Any}, x::Any) = x

convert_type(::Type{T}, x::Any) where T = T(x)

convert_type(::Type{Bool}, x::AbstractString) = lowercase(x) == "true" ? true : false

convert_type(::Type{UnitRange}, x::AbstractArray) = x[1]:x[end]
convert_type(::Type{UnitRange}, x::Any) = convert_type(UnitRange, ensurearray(x))

# [@printf("%-8s %s\n", T, fieldpropertynames(T)[T.types .== Any]) for T in subtypes(Edit) if Any in T.types]

"""
    isarray(x::Any)
Returns true/false if the the DataType or object is an array.
"""
isarray(::Type{Array{T,1}}) where T <: Any = true
isarray(x::Array{T,1}) where T <: Any = true
isarray(::Any) = false

"""
    ensurearray(x::Any)
# Returns
- `x`::AbstractArray
"""
ensurearray(x::Array{T,1}) where T <: Any = x
ensurearray(x::Tuple{Vararg{Any}}) = collect(x)
ensurearray(x::UnitRange) = collect(x)
ensurearray(x::Base.ValueIterator) = [collect(x)...;]
ensurearray(x::DataFrameRow) = ensurearray(values(x))
ensurearray(x::Any) = [x]

"""
    ensuretuple(x::Any)
# Returns
- `x`::Tuple
"""
ensuretuple(x::Tuple{Vararg{Any}}) = x
ensuretuple(x::DataFrameRow) = Tuple(x)
ensuretuple(x::Any) = tuple(x)
# !!!! add examples showing dimensionality to the documentation.

"""
    ensurefinite(x::Float64)
# Returns
- `x`, replacing `NaN` or `Inf` values with `0.0`.
"""
ensurefinite(x::Float64) = (isnan(x) || x==Inf) ? 0.0 : x

"""
"""
istype(df::DataFrame, T::DataType) = broadcast(<:, findtype(df), T)

findtype(df::DataFrame) = eltype.(eachcol(dropmissing(df)))
findtype(x::AbstractArray) = eltype(skipmissing(x))

"""
    find_oftype(df::DataFrame, T::DataType)
    find_oftype(df::Dict, T::DataType)
Returns DataFrame column names of the specified type.
"""
find_oftype(df::DataFrame, T::DataType) = propertynames(df)[istype(df, T)]
find_oftype(df::DataFrame, T::InvertedIndex{DataType}) = propertynames(df)[.!istype(df, T.skip)]

function find_oftype(d::Dict, T::DataType)
    return Dict(k => v for (k, v) in d if any(broadcast(<:, typeof.(ensurearray(v)), T)))
end

"""
    add_permutation!(set, x)
This function adds a permutation of existing set keys to the input dictionary.
If the dictionary does not contiain all of the sets specified in `x`,
the function will produce an error.

# Arguments
- `set::Dict` dictionary to update with permutations
- `x::Tuple{Symbol,1}`: set keys to permute
"""
function add_permutation!(set::Dict, x::Tuple)
    if !(x in keys(set))
        missing_keys = setdiff(ensurearray(x), keys(set))
        if !isempty(missing_keys)
            @error("Cannot create a composite $x. Key(s) $missing_keys missing from set.")
        end
        set[x] = sort(permute([[set[k] for k in x]...,]))
        end
    return set[x]
end

"""
# Returns
- `val::Array{Symbol,1}` of input DataFrame propertynames indicating values, which are
    defined as columns that DO contain `AbstractFloat` or `Bool` DataTypes.
"""
findvalue(df::DataFrame) = [find_oftype(df, AbstractFloat); find_oftype(df, Bool)]

"""
# Returns
- `idx::Array{Symbol,1}` of input DataFrame propertynames indicating indices, which are
    defined as columns that do NOT contain `AbstractFloat` or `Bool` DataTypes.
"""
findindex(df::DataFrame) = setdiff(propertynames(df), findvalue(df))

"""
# Returns
- `utx::Array{Symbol,1}` of input DataFrame propertynames with "units" in the name.
"""
findunits(df::DataFrame) = propertynames_with(df, :units)
# This function is a bit niche, but will be used heavily in the EEM.
# Should maybe keep it internal.

"""
"""
function propertynames_with(df::DataFrame, id::Symbol)
    col = propertynames(df)
    return col[occursin.(id,col)]
end

"""
Appends inputs, maintaining type of first input (currently only works for symbols)
"""
append(x::Array{Symbol,1}) = Symbol(Symbol.(x[1:end-1], :_)..., x[end])
append(x1::Symbol, x2::Any) = Symbol(x1,:_,x2)
append(x1::Symbol, x2::Vararg{Any,N}) where N = append(Symbol.([x1; ensurearray(x2)]))

append(x::Tuple) = append(ensurearray(x))

append(x::Array{String,1}) = string(string.(x[1:end-1], :_)..., x[end])
append(x1::String, x2::Any) = string(x1,:_,x2)
append(x1::String, x2::Vararg{Any,N}) where N = append(string.([x1; ensurearray(x2)]))

append(x1::Nothing, x2::Vararg{Any,N}) where N = append(x2)
append(x::Any) = x
# !!!! clean up this function.

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
abstract type GTR <: NASM end

_πA(mod::GTR) = mod.πA
_πC(mod::GTR) = mod.πC
_πG(mod::GTR) = mod.πG
_πT(mod::GTR) = mod.πT

_α(mod::GTR) = mod.α
_β(mod::GTR) = mod.β
_γ(mod::GTR) = mod.γ
_δ(mod::GTR) = mod.δ
_ϵ(mod::GTR) = mod.ϵ


text:
using HDF5, Logging, BenchmarkTools, LinearAlgebra, Statistics
import TrajectoryOptimization: get_time

IMAGE_DIR = joinpath(TrajectoryOptimization.root_dir(),"development","paper","images")

function run_step_size_comparison(model, obj, U0, group::String, Ns;
        integrations::Vector{Symbol}=[:rk3,:snopt,:ipopt],
        dt_truth=1e-3,opts=opts,infeasible=false,X0=Matrix{Float64}(undef,0,1),benchmark=false)
    solver = Solver(model, obj, integration=:rk3, N=size(U0,2))
    if isempty(X0)
        if infeasible
                X0 = line_trajectory(solver)
        else
            X0 = rollout(solver,U0)
        end
    end

    solver_truth = Solver(model, obj, dt=dt_truth)
    X_truth, U_truth = get_dircol_truth(solver_truth,X0,U0,group)[2:3]
    interp(t) = TrajectoryOptimization.interpolate_trajectory(solver_truth, X_truth, U_truth, t)

    h5open("data.h5","cw") do file
        N_group =  group * "/N_plots"
        if exists(file, N_group)
            g_parent = file[N_group]
        else
            g_parent = g_create(file, N_group)
        end

        if has(g_parent, "Ns")
            o_delete(g_parent, "Ns")
        end
        g_parent["Ns"] = Ns
    end

    # run_Ns(model, obj, X0, U0, Ns, :hermite_simpson, interp, group, opts=opts, infeasible=infeasible, benchmark=benchmark)
    for integration in integrations
        println("Starting Integration :$integration")
        run_Ns(model, obj, X0, U0, Ns, integration, interp, group, opts=opts, infeasible=infeasible, benchmark=benchmark)
    end

end

function get_dircol_truth(solver::Solver,X0,U0,group)#::Tuple{Solver, Matrix{Float64}, Matrix{Float64}}
    read_dircol_truth = false
    if isfile("data.h5")
        file = h5open("data.h5","r")
        if exists(file,group)
            if exists(file,group * "/dircol_truth")
                if exists(file,group * "/dircol_truth/dt")
                    if read(file,group * "/dircol_truth/dt")[1] == solver.dt
                        read_dircol_truth = true
                    end
                end
            end
        end
    else
        file = h5open("data.h5","cw")
    end
    if read_dircol_truth
        @info "Reading Dircol Truth from file"
        dt = Float64(read(file,group * "/dircol_truth/dt")[1])
        solver_truth = Solver(solver)
        X_truth = read(file,group * "/dircol_truth/X")
        U_truth = read(file,group * "/dircol_truth/U")
        close(file)
    else
        @info "Running Dircol Truth..."
        close(file)
        solver_truth, res_truth = run_dircol_truth(solver, X0, U0, group)
        X_truth = Array(res_truth.X)
        U_truth = Array(res_truth.U)
    end

    return solver_truth, X_truth, U_truth
end

function run_dircol_truth(solver::Solver, X0, U0, group::String)
    println("Solving DIRCOL \"truth\"")
    solver_truth = Solver(solver)
    res_truth, stat_truth = solve_dircol(solver_truth, X0, U0, method=:hermite_simpson)

    println("Writing results to file")
    h5open("data.h5","cw") do file
        group *= "/dircol_truth"
        if exists(file, group)
            o_delete(file, group)
        end
        g_truth = g_create(file, group)
        g_truth["X"] = Array(res_truth.X)
        g_truth["U"] = Array(res_truth.U)
        g_truth["dt"] = solver.dt
    end

    return solver_truth, res_truth, stat_truth
end

function run_Ns(model, obj, X0, U0, Ns, integration, interp::Function, group::String; infeasible=false, opts=SolverOptions(), benchmark=false)
    err, err_final, stats = run_Ns(model, obj, X0, U0, Ns, interp, integration=integration, infeasible=infeasible, opts=opts, benchmark=benchmark)
    save_Ns(err, err_final, stats, group, String(integration))

    return err, err_final, stats
end

function run_Ns(model, obj, X0, U0, Ns, interp; integration=:rk3, infeasible=false, opts=SolverOptions(),benchmark=false)
    num_N = length(Ns)

    err = zeros(num_N)
    err_final = zeros(num_N)
    stats = Array{Dict,1}(undef,num_N)
    disable_logging(Logging.Info)
    for (i,N) in enumerate(Ns)
        println("Solving with $N knot points")

        if integration in [:snopt, :ipopt]
            solver = Solver(model, obj, N=N, opts=opts)
            dircol_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)
            res,stat = solve_dircol(solver, X0, U0, method=:hermite_simpson, options=dircol_options, nlp=integration)
            if benchmark
                b = @benchmark solve_dircol($solver, $X0, $U0, method=:hermite_simpson, options=$dircol_options)
                times = filter(isfinite,b.times)
                stat["runtime"] = time(median(b))/1e9
                stat["std"] = std(times)/1e9
                stat["samples"] = length(times)
            else
                stat["std"] = 0
                stat["samples"] = 1
            end
            t = get_time(solver)
            Xi,Ui = interp(t)
            err[i] = norm(Xi-res.X)/N
            err_final[i] = norm(res.X[:,N] - obj.xf)
            X = Array(res.X)
        else
            solver = Solver(model,obj,N=N,opts=opts,integration=integration)
            if infeasible
                res,stat = solve(solver,X0,U0)
                if benchmark
                    b = @benchmark solve($solver,$X0,$U0)
                    times = filter(isfinite,b.times)
                    stat["runtime"] = time(median(b))/1e9
                    stat["std"] = std(times)/1e9
                    stat["samples"] = length(times)
                else
                    stat["std"] = 0
                    stat["samples"] = 1
                end
            else
                res,stat = solve(solver,U0)
                if benchmark
                    b = @benchmark solve($solver,$U0)
                    times = filter(isfinite,b.times)
                    stat["runtime"] = time(median(b))/1e9
                    stat["std"] = std(times)/1e9
                    stat["samples"] = length(times)
                else
                    stat["std"] = 0
                    stat["samples"] = 1
                end
            end
            X = to_array(res.X)
            t = get_time(solver)
            Xi,Ui = interp(t)
            err[i] = norm(Xi-to_array(res.X))/N
            err_final[i] = norm(res.X[N] - obj.xf)
        end

        stats[i] = stat
    end
    disable_logging(Logging.Debug)

    return err, err_final, stats
end

function save_Ns(err, err_final, stats, group::String, name::String)
    # Save to h5 file
    h5open("data.h5","cw") do file
        group *= "/N_plots/"
        if exists(file, group)
            g_parent = file[group]
        else
            g_parent = g_create(file, group)
        end

        if has(g_parent, name)
            o_delete(g_parent, name)
        end
        g_create(g_parent, name)
        g = g_parent[name]

        g["runtime"] = [stat["runtime"] for stat in stats]
        g["error"] = err
        g["error_final"] = err_final
        g["iterations"] = [stat["iterations"] for stat in stats]
        g["std"] = [stat["std"] for stat in stats]
        g["samples"] = [stat["samples"] for stat in stats]
        if ~isempty(stats[1]["c_max"])
            g["c_max"] = [stat["c_max"][end] for stat in stats]
        end
    end
end

function plot_stat(stat::String, group, plot_names::Vector{Symbol}; kwargs...)
    plot_names = [string(name) for name in plot_names]
    plot_stat(stat, group, plot_names; kwargs...)
end

function plot_stat(stat::String, group, plot_names::Vector{String}=["midpoint", "rk3", "rk3_foh", "rk4", "ipopt", "snopt"]; kwargs...)
    fid = h5open("data.h5","r")
    file_names = names(fid[group * "/N_plots"])
    use_names = intersect(file_names, plot_names)
    close(fid)

    if :label in keys(kwargs) && length(kwargs[:label]) == length(use_names)
        label = kwargs[:label]
    else
        label = use_names
    end

    Ns, data = load_data(stat, use_names, group)
    plot_vals(Ns, data, label, stat; kwargs...)
end

function plot_vals(Ns,vals,labels,name::String; kwargs...)
    if :color in keys(kwargs) && length(kwargs[:color]) == length(vals)
        colors = kwargs[:color]
    else
        colors = 1:length(vals)
    end
    p = plot(Ns,vals[1], label=labels[1], marker=:circle, ylabel=name, xlabel="Number of Knot Points",
        color=colors[1],markerstrokecolor=colors[1]; kwargs...)
    for (val,label,color) in zip(vals[2:end],labels[2:end],colors[2:end],labels[2:end])
        if ~isempty(val)
            plot!(Ns,val,label=label,marker=:circle,color=color,markerstrokecolor=color)
        end
    end
    p
end

function save_data(group)
    all_err = [err_mid, err_rk3, err_foh, err_rk4]
    all_eterm = [eterm_mid, eterm_rk3, eterm_foh, eterm_rk4]
    all_stats = [stats_mid, stats_rk3, stats_foh, stats_rk4]
    all_names = ["midpoint", "rk3", "rk3_foh", "rk4"]
    h5open("data.h5","cw") do file

        # Create "N_plots" group
        group *= "/N_plots"
        if exists(file, group)
            g_parent = file[group]
        else
            g_parent = g_create(file, group)
        end

        # Store Ns
        if has(g_parent, "Ns")
            o_delete(g_parent, "Ns")
        end
        g_parent["Ns"] = Ns

        for name in all_names
            if has(g_parent,name)
                o_delete(g_parent, name)
            end
            g_create(g_parent, name)
        end
        gs = [g_parent[name] for name in all_names]

        for i = 1:4
            g = gs[i]
            g["runtime"] = [stat["runtime"] for stat in all_stats[i]]
            g["error"] = all_err[i]
            g["error_final"] = all_eterm[i]
            g["iterations"] = [stat["iterations"] for stat in all_stats[i]]
            if ~isempty(all_stats[i][1]["c_max"])
                g["c_max"] = [stat["c_max"][end] for stat in all_stats[i]]
            end
        end
    end
end

function load_data(stat::String, names::Vector{String}, group)
    data = [load_data(stat, name, group)[2] for name in names]
    Ns = load_data(stat,names[1], group)[1]
    return Ns,data
end

function load_data(stat::String, name::String, group)
    h5open("data.h5","r") do file
        g_parent = file[group * "/N_plots"]
        Ns = read(g_parent, "Ns")
        data = read(g_parent[name], stat)
        return vec(Ns), data
    end
end

function constraint_plot(solver,U0::Trajectory, nesterov::Bool=true; kwargs...)
    solver.opts.verbose = false
    solver.opts.square_root = false
    solver.opts.use_nesterov = false
    res0, stats0 = solve(solver,U0)

    solver.opts.square_root = true
    res1, stats1 = solve(solver,U0)

    solver.opts.use_nesterov = true
    res2, stats2 = solve(solver,U0)

    plot(stats0["c_max"],yscale=:log10,label="normal",ylabel="max constraint violation",xlabel="iteration",linewidth=2)
    plot!(stats1["c_max"],label="square root",linewidth=2)
    plot!(stats2["c_max"],label="nesterov",linewidth=2; kwargs...)
end

function constraint_plot(solver,X0::Matrix,U0::Matrix, nesterov::Bool=true; kwargs...)
    solver.opts.verbose = false
    solver.opts.square_root = false
    solver.opts.use_nesterov = false
    res0, stats0 = solve(solver,X0,U0)

    solver.opts.square_root = true
    res1, stats1 = solve(solver,X0,U0)

    if nesterov
        solver.opts.use_nesterov = true
        res2, stats2 = solve(solver,X0,U0)
    end

    p = plot(stats0["c_max"],yscale=:log10,label="normal",ylabel="max constraint violation",xlabel="iteration",linewidth=2)
    plot!(stats1["c_max"],label="square root",linewidth=2)
    if nesterov
        plot!(stats2["c_max"],label="nesterov",linewidth=2; kwargs...)
    end
    p
end


text:
abstract type AbstractStorageForm end

abstract type BookKeepingForm <: AbstractStorageForm end

include("storage/storage_variables.jl")
include("storage/book_keeping_constraints.jl")


text:
# ---
# title: 590. N-ary Tree Postorder Traversal
# id: problem590
# author: zhwang
# date: 2022-03-12
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/n-ary-tree-postorder-traversal/description/>
# hidden: true
# ---
#
# Given an n-ary tree, return the _postorder_ traversal of its nodes' values.
#
# _Nary-Tree input serialization  is represented in their level order traversal,
# each group of children is separated by the null value (See examples)._
#
#
#
# **Follow up:**
#
# Recursive solution is trivial, could you do it iteratively?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/narytreeexample.png)
#
#
#
#     Input: root = [1,null,3,2,4,null,5,6]
#     Output: [5,6,3,2,4,1]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/08/sample_4_964.png)
#
#
#
#     Input: root = [1,null,2,3,4,5,null,null,6,7,null,8,null,9,10,null,null,11,null,12,null,13,null,null,14]
#     Output: [2,6,14,11,7,3,12,8,4,13,9,10,5,1]
#
#
#
#
# **Constraints:**
#
#   * The height of the n-ary tree is less than or equal to `1000`
#   * The total number of nodes is between `[0, 10^4]`
#
#
## @lc code=start
using LeetCode

## @lc code=end

## @lc test=start
@testset "590.n-ary-tree-postorder-traversal.jl" begin
    @test
end
## @lc test=end

## @lc info=start
# link: [solution to 590]()
## @lc info=end


text:
abstract type HeAttributes end

"""
He is a mutable composite type that is describe the attributes of the data contained
in the Helium format file..

"""
mutable struct HeInfo  <: HeAttributes
    fileName::String
    matType::DataType
    numCols::Int64
    numRows::Int64
    hasColNames::Bool
    hasRowNames::Bool
    strMiss::String
    sep::String
    skipCol::Int64
end


text:
using Gadfly, Compat

set_default_plot_size(6inch, 3inch)

plot((x,y) -> x*exp(-(x-(round(Int, x)))^2-y^2), -8., 8, -2., 2)