text:
# to get rid of eventually
const Columns = StructVector

# IndexedTable-like API

"""
    colnames(itr)

Returns the names of the "columns" in `itr`.

# Examples:

    colnames(1:3)
    colnames(Columns([1,2,3], [3,4,5]))
    colnames(table([1,2,3], [3,4,5]))
    colnames(Columns(x=[1,2,3], y=[3,4,5]))
    colnames(table([1,2,3], [3,4,5], names=[:x,:y]))
    colnames(ndsparse(Columns(x=[1,2,3]), Columns(y=[3,4,5])))
    colnames(ndsparse(Columns(x=[1,2,3]), [3,4,5]))
    colnames(ndsparse(Columns(x=[1,2,3]), [3,4,5]))
    colnames(ndsparse(Columns([1,2,3], [4,5,6]), Columns(x=[6,7,8])))
    colnames(ndsparse(Columns(x=[1,2,3]), Columns([3,4,5],[6,7,8])))

"""
function colnames end

Base.@pure colnames(t::AbstractVector) = (1,)
columns(v::AbstractVector) = v

Base.@pure colnames(t::Columns) = fieldnames(eltype(t))
Base.@pure colnames(t::Columns{<:Pair}) = colnames(t.first) => colnames(t.second)

"""
    columns(itr, select::Selection = All())

Select one or more columns from an iterable of rows as a tuple of vectors.

`select` specifies which columns to select. Refer to the [`select`](@ref) function for the
available selection options and syntax.

`itr` can be `NDSparse`, `Columns`, `AbstractVector`, or their distributed counterparts.

# Examples

    t = table(1:2, 3:4; names = [:x, :y])

    columns(t)
    columns(t, :x)
    columns(t, (:x,))
    columns(t, (:y, :x => -))
"""
function columns end

columns(c::Columns) = fieldarrays(c)
columns(c::Columns{<:Tuple}) = Tuple(fieldarrays(c))
columns(c::Columns{<:Pair}) = c.first => c.second

"""
    ncols(itr)

Returns the number of columns in `itr`.

# Examples

    ncols([1,2,3]) == 1
    ncols(rows(([1,2,3],[4,5,6]))) == 2
"""
function ncols end
ncols(c::Columns{T, C}) where {T, C} = fieldcount(C)
ncols(c::Columns{<:Pair}) = ncols(c.first) => ncols(c.second)
ncols(c::AbstractArray) = 1

summary(c::Columns{D}) where {D<:Tuple} = "$(length(c))-element Columns{$D}"

_sizehint!(c::Columns, n::Integer) = (foreach(x->_sizehint!(x,n), columns(c)); c)

function _strip_pair(c::Columns{<:Pair})
    f, s = map(columns, fieldarrays(c))
    (f isa AbstractVector) && (f = (f,))
    (s isa AbstractVector) && (s = (s,))
    Columns((f..., s...))
end

# fused indexing operations
# these can be implemented for custom vector types like PooledVector where
# you can get big speedups by doing indexing and an operation in one step.

@inline copyelt!(a, i, j) = (@inbounds a[i] = a[j])
@inline copyelt!(a, i, b, j) = (@inbounds a[i] = b[j])
@inline copyelt!(a::PooledArray, i, j) = (a.refs[i] = a.refs[j])

# row operations

copyrow!(I::Columns, i, src) = foreach(c->copyelt!(c, i, src), columns(I))
copyrow!(I::Columns, i, src::Columns, j) = foreach((c1,c2)->copyelt!(c1, i, c2, j), columns(I), columns(src))
copyrow!(I::AbstractArray, i, src::AbstractArray, j) = (@inbounds I[i] = src[j])
pushrow!(to::Columns, from::Columns, i) = foreach((a,b)->push!(a, b[i]), columns(to), columns(from))
pushrow!(to::AbstractArray, from::AbstractArray, i) = push!(to, from[i])

# test that the row on the right is "as of" the row on the left, i.e.
# all columns are equal except left >= right in last column.
# Could be generalized to some number of trailing columns, but I don't
# know whether that has applications.
@generated function row_asof(c::Columns{D,C}, i, d::Columns{D,C}, j) where {D,C}
    N = length(C.parameters)
    if N == 1
        ex = :(!isless(getfield(fieldarrays(c),1)[i], getfield(fieldarrays(d),1)[j]))
    else
        ex = :(isequal(getfield(fieldarrays(c),1)[i], getfield(fieldarrays(d),1)[j]))
    end
    for n in 2:N
        if N == n
            ex = :(($ex) && !isless(getfield(fieldarrays(c),$n)[i], getfield(fieldarrays(d),$n)[j]))
        else
            ex = :(($ex) && isequal(getfield(fieldarrays(c),$n)[i], getfield(fieldarrays(d),$n)[j]))
        end
    end
    ex
end

# map

"""
    map_rows(f, c...)

Transform collection `c` by applying `f` to each element. For multiple collection arguments, apply `f`
elementwise. Collect output as `Columns` if `f` returns
`Tuples` or `NamedTuples` with constant fields, as `Array` otherwise.

# Examples

    map_rows(i -> (exp = exp(i), log = log(i)), 1:5)
"""
function map_rows(f, iters...)
    collect_columns(f(i...) for i in zip(iters...))
end

# 1-arg case
map_rows(f, iter) = collect_columns(f(i) for i in iter)

## Special selectors to simplify column selector

"""
    All(cols::Union{Symbol, Int}...)

Select the union of the selections in `cols`. If `cols == ()`, select all columns.

# Examples

    t = table([1,1,2,2], [1,2,1,2], [1,2,3,4], [0, 0, 0, 0], names=[:a,:b,:c,:d])
    select(t, All(:a, (:b, :c)))
    select(t, All())
"""
struct All{T}
    cols::T
end

All(args...) = All(args)

"""
    Not(cols::Union{Symbol, Int}...)

Select the complementary of the selection in `cols`. `Not` can accept several arguments,
in which case it returns the complementary of the union of the selections.

# Examples

    t = table([1,1,2,2], [1,2,1,2], [1,2,3,4], names=[:a,:b,:c], pkey = (:a, :b))
    select(t, Not(:a))
    select(t, Not(:a, (:a, :b)))
"""
struct Not{T}
    cols::T
end

Not(args...) = Not(All(args))

"""
    Keys()

Select the primary keys.

# Examples

    t = table([1,1,2,2], [1,2,1,2], [1,2,3,4], names=[:a,:b,:c], pkey = (:a, :b))
    select(t, Keys())
"""
struct Keys; end

"""
    Between(first, last)

Select the columns between `first` and `last`.

# Examples

    t = table([1,1,2,2], [1,2,1,2], 1:4, 'a':'d', names=[:a,:b,:c,:d])
    select(t, Between(:b, :d))
"""
struct Between{T1 <: Union{Int, Symbol}, T2 <: Union{Int, Symbol}}
    first::T1
    last::T2
end

const SpecialSelector = Union{Not, All, Keys, Between, Function, Regex, Type}

hascolumns(t, s) = true
hascolumns(t, s::Symbol) = s in colnames(t)
hascolumns(t, s::Int) = s in 1:length(columns(t))
hascolumns(t, s::Tuple) = all(hascolumns(t, x) for x in s)
hascolumns(t, s::Not) = hascolumns(t, s.cols)
hascolumns(t, s::Between) = hascolumns(t, s.first) && hascolumns(t, s.last)
hascolumns(t, s::All) = all(hascolumns(t, x) for x in s.cols)
hascolumns(t, s::Type) = any(x -> eltype(x) <: s, columns(t))

lowerselection(t, s)                     = s
lowerselection(t, s::Union{Int, Symbol}) = colindex(t, s)
lowerselection(t, s::Tuple)              = map(x -> lowerselection(t, x), s)
lowerselection(t, s::Not)                = excludecols(t, lowerselection(t, s.cols))
lowerselection(t, s::Keys)               = lowerselection(t, IndexedTables.pkeynames(t))
lowerselection(t, s::Between)            = Tuple(colindex(t, s.first):colindex(t, s.last))
lowerselection(t, s::Function)           = colindex(t, Tuple(filter(s, collect(colnames(t)))))
lowerselection(t, s::Regex)              = lowerselection(t, x -> occursin(s, string(x)))
lowerselection(t, s::Type)               = Tuple(findall(x -> eltype(x) <: s, columns(t)))

function lowerselection(t, s::All)
    s.cols == () && return lowerselection(t, valuenames(t))
    ls = (isa(i, Tuple) ? i : (i,) for i in lowerselection(t, s.cols))
    ls |> Iterators.flatten |> union |> Tuple
end

### Iteration API

# For `columns(t, names)` and `rows(t, ...)` to work, `t`
# needs to support `colnames` and `columns(t)`

Base.@pure function colindex(t, col::Tuple)
    fns = colnames(t)
    map(x -> _colindex(fns, x), col)
end

Base.@pure function colindex(t, col)
    _colindex(colnames(t), col)
end

function colindex(t, col::SpecialSelector)
    colindex(t, lowerselection(t, col))
end

function _colindex(fnames::Union{Tuple, AbstractArray}, col, default=nothing)
    if isa(col, Int) && 1 <= col <= length(fnames)
        return col
    elseif isa(col, Symbol)
        idx = something(findfirst(isequal(col), fnames), 0)
        idx > 0 && return idx
    elseif isa(col, Pair{<:Any, <:AbstractArray})
        return 0
    elseif isa(col, Tuple)
        return 0
    elseif isa(col, Pair{Symbol, <:Pair}) # recursive pairs
        return _colindex(fnames, col[2])
    elseif isa(col, Pair{<:Any, <:Any})
        return _colindex(fnames, col[1])
    elseif isa(col, AbstractArray)
        return 0
    end
    default !== nothing ? default : error("column $col not found.")
end

# const ColPicker = Union{Int, Symbol, Pair{Symbol=>Function}, Pair{Symbol=>AbstractVector}, AbstractVector}
column(c, x) = columns(c)[colindex(c, x)]

# optimized method
@inline function column(c::Columns, x::Union{Int, Symbol})
    getfield(fieldarrays(c), x)
end

column(t, a::AbstractArray) = a
column(t, a::Pair{Symbol, <:AbstractArray}) = column(t, a[2])
column(t, a::Pair{Symbol, <:Pair}) = rows(t, a[2]) # renaming a selection
column(t, a::Pair{<:Any, <:Any}) = map(a[2], rows(t, a[1]))
column(t, s::SpecialSelector) = rows(t, lowerselection(t, s))

function columns(c, sel::Union{Tuple, SpecialSelector})
    which = lowerselection(c, sel)
    cnames = colnames(c, which)
    if all(x->isa(x, Symbol), cnames)
        tuplewrap = namedtuple(cnames...)∘tuple
    else
        tuplewrap = tuple
    end
    tuplewrap((rows(c, w) for w in which)...)
end

"""
`columns(itr, which)`

Returns a vector or a tuple of vectors from the iterator.

"""
columns(t, which) = column(t, which)

function colnames(c, cols::Union{Tuple, AbstractArray})
    map(x->colname(c, x), cols)
end

colnames(c, cols::SpecialSelector) = colnames(c, lowerselection(c, cols))

function colname(c, col)
    if isa(col, Union{Int, Symbol})
        col == 0 && return 0
        i = colindex(c, col)
        return colnames(c)[i]
    elseif isa(col, Pair{<:Any, <:Any})
        return col[1]
    elseif isa(col, Tuple)
        #ns = map(x->colname(c, x), col)
        return 0
    elseif isa(col, SpecialSelector)
        return 0
    elseif isa(col, AbstractVector)
        return 0
    end
    error("column named $col not found")
end

"""
    rows(itr, select = All())

Select one or more fields from an iterable of rows as a vector of their values.  Refer to
the [`select`](@ref) function for selection options and syntax.

`itr` can be [`NDSparse`](@ref), `StructArrays.StructVector`, `AbstractVector`, or their distributed counterparts.

# Examples

    t = table([1,2],[3,4], names=[:x,:y])
    rows(t)
    rows(t, :x)
    rows(t, (:x,))
    rows(t, (:y, :x => -))
"""
function rows end

rows(x::AbstractVector) = x
rows(cols::Tup) = Columns(cols)

rows(t, which...) = rows(columns(t, which...))

# Replace empty Columns object with one of correct length and eltype
replace_placeholder(t, ::Columns{Tuple{}}) = fill(Tuple(), length(t))
replace_placeholder(t, ::Columns{NamedTuple{(), Tuple{}}}) = fill(NamedTuple(), length(t))
replace_placeholder(t, cols) = cols

_cols_tuple(xs::Columns) = columns(xs)
_cols_tuple(xs::AbstractArray) = (xs,)
concat_cols(xs, ys) = rows(concat_tup(_cols_tuple(xs), _cols_tuple(ys)))

## Mutable Columns Dictionary

mutable struct ColDict{T}
    pkey::Vector{Int}
    src::T
    names::Vector
    columns::Vector
    copy::Union{Nothing, Bool}
end

"""
    d = ColDict(t)

Create a mutable dictionary of columns in `t`.

To get the immutable iterator of the same type as `t`
call `d[]`
"""
function ColDict(t; copy=nothing)
    cnames = colnames(t)
    if cnames isa AbstractArray
        cnames = Base.copy(cnames)
    end
    ColDict(Int[], t, convert(Array{Any}, collect(cnames)), Any[columns(t)...], copy)
end

Base.keys(d::ColDict) = d.names
Base.values(d::ColDict) = d.columns

function Base.getindex(d::ColDict{<:Columns})
    Columns(Tuple(d.columns); names=d.names)
end

Base.getindex(d::ColDict, key) = rows(d[], key)
Base.getindex(d::ColDict, key::AbstractArray) = key

function Base.setindex!(d::ColDict, x, key::Union{Symbol, Integer})
    k = _colindex(d.names, key, 0)
    col = d[x]
    if k == 0
        push!(d.names, key)
        push!(d.columns, col)
    elseif k in d.pkey
        # primary key column has been modified.
        # copy the table as this results in a shuffle
        if d.copy === nothing
            d.copy = true
        end
        d.columns[k] = col
    else
        d.columns[k] = col
    end
end

Base.@deprecate set!(d::ColDict, key, x) setindex!(d, x, key)

transform!(d::ColDict, changes::Pair...) = transform!(d, changes)

function transform!(d::ColDict, changes)
    Base.foreach(changes) do (key, val)::Pair
        d[key] = val
    end
end

function Base.haskey(d::ColDict, key)
    _colindex(d.names, key, 0) != 0
end

function Base.insert!(d::ColDict, index::Integer, (key, col)::Pair)
    if haskey(d, key)
        error("Key $key already exists. Use dict[key] = col instead of inserting.")
    else
        insert!(d.names, index, key)
        insert!(d.columns, index, rows(d.src, col))
        for (i, pk) in enumerate(d.pkey)
            if pk >= index
                d.pkey[i] += 1 # moved right
            end
        end
    end
end

function Base.insert!(d::ColDict, index::Integer, newcols)
    for new::Pair in newcols
        insert!(d, index, new)
        index += 1
    end
end

Base.insert!(d::ColDict, index::Integer, newcols::Pair...) = insert!(d, index, newcols)

function insertafter!(d::ColDict, i, args...)
    k = _colindex(d.names, i, 0)
    if k == 0
        error("$i not found. Cannot insert column after $i")
    end
    insert!(d, k+1, args...)
end

function insertbefore!(d::ColDict, i, args...)
    k = _colindex(d.names, i, 0)
    if k == 0
        error("$i not found. Cannot insert column after $i")
    end
    insert!(d, k, args...)
end

function rename!(d::ColDict, (col, newname)::Pair)
    k = _colindex(d.names, col, 0)
    if k == 0
        error("$col not found. Cannot rename it.")
    end
    d.names[k] = newname
end

@deprecate rename!(t::ColDict, col::Union{Symbol, Integer}, newname) rename!(t, col => newname)

rename!(t::ColDict, changes) = Base.foreach(change::Pair -> rename!(t, change), changes)
rename!(t::ColDict, changes::Pair...) = rename!(t, changes)

function Base.push!(d::ColDict, (key, x)::Pair)
    push!(d.names, key)
    push!(d.columns, rows(d.src, x))
end

function _cols(expr)
    if expr.head == :call
        dict = :(dict = ColDict($(expr.args[2])))
        expr.args[2] = :dict
        quote
            let $dict
                $expr
                dict[]
            end
        end |> esc
    else
        error("This form of @cols is not implemented. Use `@cols f(t,args...)` where `t` is the collection.")
    end
end

macro cols(expr)
    _cols(expr)
end

# Modifying columns

"""
    transform(t::Table, changes::Pair...)

Transform columns of `t`. For each pair `col => value` in `changes` the column `col` is replaced
by the `AbstractVector` `value`. If `col` is not an existing column, a new column is created.

# Examples:

    t = table([1,2], [3,4], names=[:x, :y])

    # change second column to [5,6]
    transform(t, 2 => [5,6])
    transform(t, :y => :y => x -> x + 2)

    # add [5,6] as column :z
    transform(t, :z => 5:6)
    transform(t, :z => :y => x -> x + 2)

    # replacing the primary key results in a re-sorted copy
    t = table([0.01, 0.05], [1,2], [3,4], names=[:t, :x, :y], pkey=:t)
    t2 = transform(t, :t => [0.1,0.05])

    # the column :z is not part of t so a new column is added
    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    pushcol(t, :z => [1//2, 3//4])
"""
transform(t, args...) = @cols transform!(t, args...)

@deprecate setcol(t, args::Pair...) transform(t, args...)
@deprecate setcol(t, key::Union{Int, Symbol}, val) transform(t, key => val)
@deprecate setcol(t, args) transform(t, args)

@deprecate pushcol(t, args::Pair...) transform(t, args...)
@deprecate pushcol(t, key::Union{Int, Symbol}, val) transform(t, key => val)
@deprecate pushcol(t, args) transform(t, args)

@deprecate popcol(t, args...) select(t, Not(args...))
@deprecate popcol(t) select(t, Not(ncols(t)))

"""
    insertcols(t, position::Integer, map::Pair...)

For each pair `name => col` in `map`, insert a column `col` named `name` starting at `position`.
Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    insertcol(t, 2, :w => [0,1])
"""
insertcols(t, i::Integer, args...) = @cols insert!(t, i, args...)

@deprecate insertcol(t, i, name, x) insertcols(t, i, name => x)

"""
    insertcolsafter(t, after, map::Pair...)

For each pair `name => col` in `map`, insert a column `col` named `name` after `after`.
Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    insertcolsafter(t, :t, :w => [0,1])
"""
insertcolsafter(t, after, args...) = @cols insertafter!(t, after, args...)

@deprecate insertcolafter(t, i, name, x) insertcolsafter(t, i, name => x)

"""
insertcolsbefore(t, before, map::Pair...)

For each pair `name => col` in `map`, insert a column `col` named `name` before `before`.
Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    insertcolsbefore(t, :x, :w => [0,1])
"""
insertcolsbefore(t, before, args...) = @cols insertbefore!(t, before, args...)

@deprecate insertcolbefore(t, i, name, x) insertcolsbefore(t, i, name => x)

"""
    rename(t, map::Pair...)

For each pair `col => newname` in `map`, set `newname` as the new name for column `col` in `t`.
Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], names=[:t, :x])
    rename(t, :t => :time)
"""
rename(t, args...) = @cols rename!(t, args...)

@deprecate renamecol(t, args...) rename(t, args...)

## Utilities for mapping and reduction with many functions / OnlineStats

@inline _apply(f::OnlineStat, g, x) = (fit!(g, x); g)
@inline _apply(f::Tup, y::Tup, x::Tup) = _apply(astuple(f), astuple(y), astuple(x))
@inline _apply(f::Tuple, y::Tuple, x::Tuple) = map(_apply, f, y, x)
@inline _apply(f::NamedTuple, y::NamedTuple, x::NamedTuple) = map(_apply, f, y, x)
@inline _apply(f, y, x) = f(y, x)
@inline _apply(f::Tup, x::Tup) = _apply(astuple(f), astuple(x))
@inline _apply(f::NamedTuple, x::NamedTuple) = map(_apply, f, x)
@inline _apply(f::Tuple, x::Tuple) = map(_apply, f, x)
@inline _apply(f, x) = f(x)

@inline init_first(f, x) = x
@inline init_first(f::OnlineStat, x) = (g=copy(f); fit!(g, x); g)
@inline init_first(f::Tup, x::Tup) = map(init_first, f, x)

# Initialize functions to apply and input vectors

function init_inputs(f, x, isvec) # normal functions
    f, x
end

nicename(f::Function) = typeof(f).name.mt.name
nicename(f) = Symbol(last(split(string(f), ".")))
nicename(o::OnlineStat) = Symbol(typeof(o).name.name)

init_funcs(f, isvec) = init_funcs((f,), isvec)

function init_funcs(f::Tup, isvec)
    if isa(f, NamedTuple)
        return init_funcs((map(Pair, fieldnames(typeof(f)), f)...,), isvec)
    end

    funcmap = map(f) do g
        if isa(g, Pair)
            name = g[1]
            if isa(g[2], Pair)
                sel, fn = g[2]
            else
                sel = nothing
                fn = g[2]
            end
            (name, sel, fn)
        else
            (nicename(g), nothing, g)
        end
    end

    ns = map(x->x[1], funcmap)
    ss = map(x->x[2], funcmap)
    fs = map(map(x->x[3], funcmap)) do f
        f
    end

    NamedTuple{(ns...,)}((fs...,)), ss
end

function init_inputs(f::Tup, input, isvec)
    if isa(f, NamedTuple)
        return init_inputs((map(Pair, fieldnames(typeof(f)), f)...,), input, isvec)
    end
    fs, selectors = init_funcs(f, isvec)

    xs = map(s->s === nothing ? input : rows(input, s), selectors)

    ns = fieldnames(typeof(fs))
    NT = namedtuple(ns...)

    # functions and input
    NT((fs...,)), rows(NT((xs...,)))
end

# utils
refs(v::Columns) = Columns(map(refs, fieldarrays(v)))
compact_mem(v::Columns) = replace_storage(compact_mem, v)


text:
import Base: adjoint, getindex
import Plots: plot

abstract type Basis end

abstract type periodic end
abstract type dirichlet end

"""
    BSplineBasis{T}(basis_functions, knot_vector, interval, order, N)

    BSpline basis on `interval` with `N` points and order `k`. Where `T` describes the boundary condition.
"""
struct BSplineBasis{T} <: Basis where {T}
    B #need to define proper type, contains the basis functions
    t::AbstractVector # knotvector
    interval::AbstractInterval # domain
    k::Int # order
    N::Int
end

function Base.show(io::IO, B::BSplineBasis{T}) where {T}
    println(io, T, " BSpline basis of order ", B.k, " on interval ", B.interval, " with ", B.N, " nodes.")
    nothing
end

getindex(B::BSplineBasis, x::Real, j::Integer) = B.B[j](x)
getindex(B::BSplineBasis, x::Real, j::Colon) = [B.B[i](x) for i in 1:length(B.B)]
getindex(B::BSplineBasis, x::AbstractVector, j::Integer) = B.B[j].(x)
getindex(B::BSplineBasis, x::AbstractVector, j::Colon) = [B.B[i](x[j]) for j in 1:length(x), i in 1:length(B.B)]
# Calculate the Spline Basis of order k
# For the Splines, we use CompactBases.jl. For the future, we should write our own libary with a better interface.
"""
    PeriodicBSpline(basis_functions, knot_vector, interval, order, N)

    Construct a `periodic` BSpline basis on `interval` with `N` points and order `k`.
"""
function PeriodicBSpline(interval::AbstractInterval, k::Int, N::Int)
    a = interval.left
    b = interval.right
    Δx = (b - a)/N

    # Create a knotsequence for the cardinal B-Splines
    t_ = range(a-k*Δx, b+k*Δx, step=Δx)
    t = ArbitraryKnotSet(k, t_)
    B = BSpline(t)

    basis_fct = []
    # Make the functions periodic
    for i in k+1:2*k-1
        push!(basis_fct, x -> B[x, i] + B[x, N+i])
    end

    for i in 2*k:N+k
        push!(basis_fct, x -> B[x, i])
    end

    # In order to have the same ordering as in theory
    basis_fct = circshift(basis_fct, -1)

    return BSplineBasis{periodic}(SVector(basis_fct...), t_, interval, k, N)
end

# For the second derivative, the last basis functions has problems with the BC
function derivative(basis::BSplineBasis{periodic})
    k = basis.k-1
    t_d = ArbitraryKnotSet(k, basis.t)
    B_d = BSpline(t_d)

    interval = basis.interval
    a = interval.left
    b = interval.right
    N = basis.N
    Δx = (b - a)/basis.N

    basis_fct = []
    for i in k+2:2*k
        push!(basis_fct, x -> 1/Δx * (B_d[x, i] + B_d[x, N+i]))
    end

    for i in 2*k+1:N+k+1
        push!(basis_fct, x -> 1/Δx * B_d[x, i])
    end

    return BSplineBasis{periodic}(SVector(basis_fct...), basis.t, interval, k, N)
end

function DirichletBSpline(interval, k::Int, N::Int)
    a = interval.left
    b = interval.right
    Δx = (b - a)/N

    t = LinearKnotSet(k, a, b, N)
    B = BSpline(t)

    basis_fct = []

    for i in 1:size(B)[2]
        push!(basis_fct, x -> B[x, i])
    end

    return BSplineBasis{dirichlet}(SVector(basis_fct...), t.t, interval, k, N)
end

function derivative(basis::BSplineBasis{dirichlet})

    interval = basis.interval
    a = interval.left
    b = interval.right
    N = basis.N
    Δx = (b - a)/basis.N

    k = basis.k-1
    t_d = LinearKnotSet(k, a, b, N)
    B_d = BSpline(t_d)

    basis_fct = []

    for i in 1:size(B_d)[2]
        push!(basis_fct, x -> k/(t_d[i+k] - t_d[i]) * B_d[x, i])
    end

    return BSplineBasis{dirichlet}(SVector(basis_fct...), t_d.t, interval, k, N)
end

"""
    adjoint(b::Basis)

    Create a basis that consits of the derivative of basis functions in `b`.
"""
function adjoint(basis::BSplineBasis)
    return derivative(basis)
end

# We can delete this later on
function plot(basis::Basis)
    a = basis.interval.left
    b = basis.interval.right
    xx = range(a, b, length=1000)
    fig = plot()
    for i in 1:length(basis.B)
        plot!(fig, xx, basis[xx, i])
    end

    return fig
end

# You can do this much more elegant, using flattening...
abstract type TensorProductBasis{N} end

# Just collect the 1-D bases as a tuple
"""
    BSplineTensorProductBasis{N,T}(Bases::SArray)

    BSpline basis on a tensor-product domain, consisting of a vector of one-dimensional bases defined on each direction of the tensor-product domain.
    Where `N` describes the dimension and `T` the boundary condition.
"""
struct BSplineTensorProductBasis{N, T} <: TensorProductBasis{N}
    B::SArray{Tuple{N}, BSplineBasis{T}}
end

function Base.show(io::IO, B::BSplineTensorProductBasis{N,T}) where {N,T}
    domain = B[1].interval
    for i in 2:N
        domain = domain × B[i].interval
    end

    println(io, T, " BSpline basis on tensor-product domain ", domain, ", consisting of: ")
    for i in 1:N
        show(io::IO, B.B[i])
    end
    nothing
end

getindex(B::BSplineTensorProductBasis{N, T}, j::Integer) where {N, T} = B.B[j]
getindex(B::BSplineTensorProductBasis{N, T}, j::Colon) where {N, T} = B.B

function PeriodicBSplineTensorProductBasis(domain::ProductDomain, k::NTuple{d,Int}, N::NTuple{d,Int}) where {d}
    basis = []
    for i in 1:d
        push!(basis, PeriodicBSpline(domain.domains[i], k[i], N[i]))
    end
    BSplineTensorProductBasis{d, periodic}(SVector(basis...))
end

function DirichletBSplineTensorProductBasis(domain::ProductDomain, k::NTuple{d,Int}, N::NTuple{d,Int}) where {d}
    basis = []
    for i in 1:d
        push!(basis, DirichletBSpline(domain.domains[i], k[i], N[i]))
    end
    BSplineTensorProductBasis{d, dirichlet}(SVector(basis...))
end

function tensorproduct(a::BSplineBasis{T}, b::BSplineBasis{T}) where T
    Bs = SVector(a, b)
    BSplineTensorProductBasis{2, T}(SVector(a, b))
end

function tensorproduct(a::BSplineTensorProductBasis{d, T}, b::BSplineBasis{T}) where {d, T}
    B = SVector(a.B..., b)
    BSplineTensorProductBasis{d+1, T}(B)
end

function tensorproduct(a::BSplineBasis{T}, b::BSplineTensorProductBasis{d, T}) where {d, T}
    B = SVector(a, b.B...)
    BSplineTensorProductBasis{d+1, T}(B)
end

function tensorproduct(a::BSplineTensorProductBasis{d1, T}, b::BSplineTensorProductBasis{d2, T}) where {d1, d2, T}
    B = SVector(a.B..., b.B...)
    BSplineTensorProductBasis{d1+d2, T}(B)
end

⊗(a::Union{TensorProductBasis, Basis}, b::Union{TensorProductBasis, Basis}) = tensorproduct(a, b)

function derivative(basis::BSplineTensorProductBasis{d, T}) where {d, T}
    b = derivative(basis.B[1])
    for i in 2:d
        bn = derivative(basis.B[i])
        b = b ⊗ bn
    end

    return b
end

function adjoint(basis::TensorProductBasis)
    return derivative(basis)
end

ext: jl


text:
# by Lazaro Alonso
using GLMakie
GLMakie.activate!()
let
    n = 50
    x = y = z = LinRange(-5,5,n)
    vals = zeros(n,n,n)
    for i in 1:n, j in 1:n, k in 1:n
        vals[i,j,k] = 0.5*i^2  + j^2 + 2*k^2
    end
    cmap = :Spectral_11
    fig, _ = contour(x, y, z, vals, colormap = cmap, levels=[100,2500,7000],
         algorithm=:mip, alpha = 0.05, ambient = Vec3f0(0.65, 0.65, 0.65))
    fig
end

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
mutable struct Cell <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    status::Bool
end

"""
``` julia
game_of_life(;
    rules::Tuple = (2, 3, 3, 3),
    dims = (100, 100),
    metric = :chebyshev
)
```
Same as in [Conway's game of life](@ref).
"""
function game_of_life(;
    rules::Tuple = (2, 3, 3, 3),
    dims = (100, 100),
    metric = :chebyshev
)
    space = GridSpace(dims; metric = metric)
    properties = Dict(:rules => rules)
    model = ABM(Cell, space; properties = properties)
    idx = 1
    for x in 1:dims[1]
        for y in 1:dims[2]
            add_agent_pos!(Cell(idx, (x, y), false), model)
            idx += 1
        end
    end
    return model, dummystep, game_of_life_model_step!
end

function game_of_life_model_step!(model)
    new_status = fill(false, nagents(model))
    for agent in allagents(model)
        nlive = nlive_neighbors(agent, model)
        if agent.status == true && (nlive ≤ model.rules[4] && nlive ≥ model.rules[1])
            new_status[agent.id] = true
        elseif agent.status == false && (nlive ≥ model.rules[3] && nlive ≤ model.rules[4])
            new_status[agent.id] = true
        end
    end

    for k in keys(model.agents)
        model.agents[k].status = new_status[k]
    end
end

function nlive_neighbors(agent, model)
    neighbor_positions = nearby_positions(agent, model)
    all_neighbors = Iterators.flatten(ids_in_position(np,model) for np in neighbor_positions)
    sum(model[i].status == true for i in all_neighbors)
end


text:
"MetaFormatter for ConsoleLogger Inspired by (i.e., shamelessly copy-pasted from) PowerModelsDistribution v0.11"
function _pmdse_metafmt(level::Logging.LogLevel, _module, group, id, file, line)
    @nospecialize
    color = Logging.default_logcolor(level)
    prefix = "$(_module) | " * (level == Logging.Warn ? "Warning" : string(level)) * " ] :"
    suffix = ""
    Logging.Info <= level < Logging.Warn && return color, prefix, suffix
    _module !== nothing && (suffix *= "$(_module)")
    if file !== nothing
        _module !== nothing && (suffix *= " ")
        suffix *= Base.contractuser(file)
        if line !== nothing
            suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
        end
    end
    !isempty(suffix) && (suffix = "@ " * suffix)

    return color, prefix, suffix
end

"Sets loglevel for PMDSE to :Warn, silencing Info. PMD starts from :Error"
function silence!()
    set_logging_level!(:Warn)
end

"Resets the log level to Info"
function reset_logging_level!()
    Logging.global_logger(_LOGGER)

    return
end

"Restores the global logger to its default state (before PMDSE was loaded)"
function restore_global_logger!()
    Logging.global_logger(_DEFAULT_LOGGER)

    return
end

"Sets the logging level for PMDSE: :Info, :Warn, :Error"
function set_logging_level!(level::Symbol)
    Logging.global_logger(_make_filtered_logger(getfield(Logging, level)))

    return
end

"Helper function to create the filtered logger for PMDSE"
function _make_filtered_logger(level)
    LoggingExtras.EarlyFilteredLogger(_LOGGER) do log
        if log._module == PowerModelsDistributionStateEstimation && log.level < level
            return false
        else
            return true
        end
    end
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function altaz2hadec(alt::T, az::T, lat::T) where {T<:AbstractFloat}
    # Convert to radians.
    alt_r = deg2rad(alt)
    az_r = deg2rad(az)
    lat_r = deg2rad(lat)
    sin_alt_r, cos_alt_r = sincos(alt_r)
    sin_az_r, cos_az_r = sincos(az_r)
    sin_lat_r, cos_lat_r = sincos(lat_r)
    # Find local hour angle (in degrees, from 0. to 360.).
    ha = rad2deg(atan(-sin_az_r * cos_alt_r,
                      -cos_az_r * sin_lat_r * cos_alt_r +
                      sin_alt_r * cos_lat_r))
    ha = mod(ha, 360)
    # Find declination (positive if north of Celestial Equator, negative if
    # south)
    sindec = sin_lat_r * sin_alt_r + cos_lat_r * cos_alt_r * cos_az_r
    dec = rad2deg(asin(sindec))  # convert dec to degrees
    return ha, dec
end

"""
    altaz2hadec(alt, az, lat) -> ha, dec

### Purpose ###

Convert Horizon (Alt-Az) coordinates to Hour Angle and Declination.

### Explanation ###

Can deal with the NCP singularity.  Intended mainly to be used by program
`hor2eq`.

### Arguments ###

Input coordinates may be either a scalar or an array, of the same dimension.

* `alt`: local apparent altitude, in degrees, scalar or array.
* `az`: the local apparent azimuth, in degrees, scalar or vector, measured
  *east* of *north*!!!  If you have measured azimuth west-of-south (like the book
  Meeus does), convert it to east of north via: `az = (az + 180) % 360`.
* `lat`: the local geodetic latitude, in degrees, scalar or array.

`alt` and `az` can be given as a 2-tuple `(alt, az)`.

### Output ###

2-tuple `(ha, dec)`

* `ha`: the local apparent hour angle, in degrees.  The hour angle is the time
  that right ascension of 0 hours crosses the local meridian.  It is
  unambiguously defined.
* `dec`: the local apparent declination, in degrees.

The output coordinates are always floating points and have the same type (scalar
or array) as the input coordinates.

### Example ###

Arcturus is observed at an apparent altitude of 59d,05m,10s and an azimuth
(measured east of north) of 133d,18m,29s while at the latitude of +43.07833
degrees.  What are the local hour angle and declination of this object?

```jldoctest
julia> using AstroLib

julia> ha, dec = altaz2hadec(ten(59,05,10), ten(133,18,29), 43.07833)
(336.6828582472844, 19.182450965120402)
```
The widely available XEPHEM code gets:

    Hour Angle = 336.683
    Declination = 19.1824

### Notes ###

`hadec2altaz` converts Hour Angle and Declination to Horizon (Alt-Az)
coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
altaz2hadec(alt::Real, az::Real, lat::Real) =
    altaz2hadec(promote(float(alt), float(az), float(lat))...)

altaz2hadec(altaz::Tuple{Real, Real}, lat::Real) = altaz2hadec(altaz..., lat)

function altaz2hadec(alt::AbstractArray{R}, az::AbstractArray{<:Real},
                     lat::AbstractArray{<:Real}) where {R<:Real}
    @assert length(alt) == length(az) == length(lat)
    typealt = float(R)
    ha  = similar(alt, typealt)
    dec = similar(alt, typealt)
    for i in eachindex(alt)
        ha[i], dec[i] = altaz2hadec(alt[i], az[i], lat[i])
    end
    return ha, dec
end


text:
abstract type StencilOrder end
struct StencilOrder4 <: StencilOrder end
struct StencilOrder8 <: StencilOrder end
show(io::IO, s::StencilOrder4) = write(io, "O4")
show(io::IO, s::StencilOrder8) = write(io, "O8")

# order 4 deriv :  a1*(f3-f2) + a2*(f4-f1)
@inline dfdr4(a1,a2,f1,f2,f3,f4) = a1*(f3 - f2) + a2*(f4 - f1)

# order 8 deriv : a1*(f5-f4) + a2*(f6-f3) + a3*(f7-f2) + a4*(f8-f1)
@inline dfdr8(a1,a2,a3,a4,f1,f2,f3,f4,f5,f6,f7,f8) = a1*(f5-f4) + a2*(f6-f3) + a3*(f7-f2) + a4*(f8-f1)

function stencilcoeffs(stenciltype::Symbol, Q::Int, D::Int, ndim::Int, dtmod::Real, velmax::Real, mindel::Real, T::Type)
    if stenciltype == :Fornberg
        return stencilcoeffs_fornberg(Q,D,T)
    end
    if stenciltype == :Nihei
        return stencilcoeffs_nihei(Q,ndim,velmax,mindel,dtmod,T)
    end
    throw(ArgumentError("Unknonw stenciltype: $(stenciltype)"))
end

courantnumber(a::Array{T,1}, ndim) where {T<:Real} = T(1.0/(sqrt(ndim)*sum(abs.(a))))

function stencilcoeffs_fornberg(Q,D,T)
    a = D == 1 ? Array{T}(undef, div(Q,2)) : Array{T}(undef, div(Q,2)+1)
    if Q == 4 && D == 1 # assumes staggered
        a[1] = 9.0/8.0
        a[2] = -1.0/24.0
    elseif Q == 8 && D == 1 # assumes staggered
        a[1] = +1225.0/1024.0
        a[2] = -245.0/3072.0
        a[3] = +49.0/5120.0
        a[4] = -5.0/7168.0
    elseif Q == 8 && D == 2 # assumes centered
        a[1] = -205.0/72.0 # centered coeff
        a[2] = 8.0/5.0
        a[3] = -1.0/5.0
        a[4] = 8.0/315.0
        a[5] = -1.0/560.0
    else
        throw(ArgumentError("Unknown combo of accruacy order=$(Q) and derivative order=$(D)"))
    end
    a
end

function stencilcoeffs_nihei(Q,ndim,velmax,mindel,dtmod,T;nkdx=500,pts_per_wavelength=2.8,rtol=1e-9)
    a = stencilcoeffs_fornberg(Q, T)
    csafe = 0.95*courantnumber(a, ndim) # Approximate and safe courant number of O(8)
    c = velmax*dtmod/mindel/.99
    @assert dtmod > 0 && c < csafe
    kdx_max = 2pi/pts_per_wavelength
    dkdx = kdx_max/nkdx

    # Newton's method
    itmax=10
    aold = copy(a)
    for it = 1:itmax
        g = stencilcoeffs_nihei_gradient(a,Q,c,nkdx,dkdx)
        H = stencilcoeffs_nihei_hessian(a,Q,c,nkdx,dkdx)
        a -= H\g
        phi = stencilcoeffs_nihei_cost(a,Q,c,nkdx,dkdx)
        if norm(a-aold) < rtol
            break
        end
        copyto!(aold,a)
    end
    a
end

function stencilcoeffs_nihei_gradient(a, Q, c, nkdx, dkdx)
    Q2 = round(Int,Q/2)
    grad = zeros(Q2)
    for ikdx = 1:nkdx
        kdx = ikdx*dkdx

        A = 2.0/(kdx*c)
        E = 0.0
        for r = 1:Q2
            E += c*(a[r]*sin(kdx*(r - 0.5)))
        end
        vnorm = A*asin(E)  # normalized velocity
        diff_vnorm =  vnorm - 1.0  # normalized velocity misfit
        for p = 1:Q2
            B = c*sin(kdx*(p - 0.5))
            grad[p] += diff_vnorm*A*B/sqrt(1.0 - E^2)  # Gradient for a[p]
        end
    end
    grad
end

function stencilcoeffs_nihei_hessian(a, Q, c, nkdx, dkdx)
    Q2 = round(Int,Q/2)
    hess = zeros(Q2, Q2)
    for ikdx = 1:nkdx
        kdx = convert(Float64, ikdx)*dkdx

            A = 2.0/(kdx*c)
        E = 0.0
        for r = 1:Q2
            E += c*(a[r]*sin(kdx*(r - 0.5)))
        end
        vnorm = A*asin(E)  # normalized velocity
        diff_vnorm =  vnorm - 1.0  # normalized velocity misfit

        for p = 1:Q2
            B = c*sin(kdx*(p - 0.5))
            for q = 1:Q2
                D = c*sin(kdx*(q - 0.5))
                hess[p,q] += diff_vnorm*A*B*D*E/(1.0 - E^2)^1.5 + (A^2)*B*D/(1.0 - E^2)  # Hpq
            end
        end
    end
    hess
end

function stencilcoeffs_nihei_cost(a, Q, c, nkdx, dkdx)
    Q2 = round(Int,Q/2)
    obj = 0.0
    for ikdx = 1:nkdx
        kdx = convert(Float64, ikdx)*dkdx

        A = 2.0/(kdx*c)
        E = 0.0
        for r = 1:Q2
            E += c*(a[r]*sin(kdx*(r - 0.5)))
        end
        vnorm = A*asin(E)  # normalized velocity
        obj += 0.5*(vnorm - 1.0)^2  # L2 norm objective function of normalized velocity misfit
    end
    return obj
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
"""
Alternate positional forms

"""
module Alternate

# For reference only in deference to the form below with keyword step
@static if false
"""
    range( length, start = 1, stop = length; step = nothing )

Alternate one, two, and three argument range argument `range`. Compact.
"""
Base.range(length, start = 1, stop = nothing) =
    start === 1 && stop === nothing ?
        Base.OneTo(length) :
        Base.range(start, stop; length)
"""
```julia
julia> Base.range(5)
Base.OneTo(5)

julia> Base.range(5, 2)
2:6

julia> Base.range(5, 2, 11)
2.0:2.25:11.0
```
"""

"""
    range( length, start = 1, stop = length; step = nothing )

Alternate one, two, and three argument range argument `range`
"""
function Base.range(length, start = 1, stop = nothing; step = nothing)
    if step === nothing
        if start === 1 && stop === length
            Base.OneTo(length)
        else
            Base._range(start, step, stop, length)
        end
    elseif stop === nothing
        Base._range(start, step, stop, length)
    else
        throw(ArgumentError("Cannot specify all of length = $length, start = $start, stop = $stop, and step = $step"))
    end
end
"""
```julia
julia> Base.range(5)
Base.OneTo(5)

julia> Base.range(5, 2)
2:6

julia> Base.range(5, 2, 11)
2.0:2.25:11.0
```
"""
end

end


text:
module SAT
using Random
using ExtractMacro
using Printf

export solve, energy,
        CNF, randomcnf, readcnf, writecnf

include("cnf.jl")
include("belief_propagation.jl")
include("gradient_descent.jl")

end # module


text:
#=
in this example, we just put a batch of the data into the
network at once over and over again to optimizes the
parameters of network.
=#

using Gadfly
using Random

# first class data and label
input1 = 3*randn(2,2000);
label1 = zeros(2,2000);
label1[1,:] .= 1.0;

# second class data and label
t = 0:pi/1000:2*pi;
input2 = zeros(2,length(t));
label2 = zeros(2,length(t));
label2[2,:] .= 1.0;
for i = 1:length(t)
    r = randn() + 16.0;
    input2[1,i] = r * cos(t[i]);
    input2[2,i] = r * sin(t[i]);
end

input = hcat(input1,input2)
label = hcat(label1,label2)

batch = 2048
epoch = 1000
lrate = 1e-4
topology = [2, 32,32,32, 2]
operator = [relu, relu, relu, relu]
indexs   = collect(1:batch)

mlpmodel = MLP(topology, operator)
paramter = paramsof(mlpmodel)
lossval  = zeros(epoch,1)
tic = time()
for i = 1:epoch
    ids  = shuffle!(indexs)
    outs = forward(mlpmodel, Variable(input[:,ids[1:batch]]))
    outs = softmax(outs; dims=1)
    COST = crossEntropyCost(outs, Variable(label[:,ids[1:batch]]))
    backward()
    update!(paramter,lrate)
    zerograds!(paramter)
    lossval[i] = COST
end
toc = time()
println("\n------- sgd optimizer loss: ", lossval[end])
l1 = layer(y=log.(lossval), Geom.line,Theme(default_color=colorant"cyan"))

mlpmodel  = MLP(topology, operator)
paramter  = paramsof(mlpmodel)
optimizer = Momentum(paramter,lr=lrate,inertia=0.9)
lossval   = zeros(epoch,1)
tic = time()
for i = 1:epoch
    ids  = shuffle!(indexs)
    outs = forward(mlpmodel, Variable(input[:,ids[1:batch]]))
    outs = softmax(outs; dims=1)
    COST = crossEntropyCost(outs, Variable(label[:,ids[1:batch]]))
    backward()
    update!(optimizer,paramter)
    zerograds!(paramter)
    lossval[i] = COST
end
toc = time()
println("-- Momentum optimizer loss: ", lossval[end])
l2 = layer(y=log.(lossval), Geom.line,Theme(default_color=colorant"yellow"))

mlpmodel  = MLP(topology, operator)
paramter  = paramsof(mlpmodel)
optimizer = Adam(paramter,lr=1e-3)
lossval   = zeros(epoch,1)
tic = time()
for i = 1:epoch
    ids  = shuffle!(indexs)
    outs = forward(mlpmodel, Variable(input[:,ids[1:batch]]))
    outs = softmax(outs; dims=1)
    COST = crossEntropyCost(outs, Variable(label[:,ids[1:batch]]))
    backward()
    update!(optimizer,paramter)
    zerograds!(paramter)
    lossval[i] = COST
end
toc = time()
println("------ Adma optimizer loss: ", lossval[end])
l3 = layer(y=log.(lossval), Geom.line,Theme(default_color=colorant"red"))

plot(l1,l2,l3)

mlpmodel  = Chain(
    Maxout(2,32,k=3),
    Residual(Dense(32,32),Dense(32,32),Linear(32,32)),
    Residual(Dense(32,32),Dense(32,32),Linear(32,32)),
    Residual(Dense(32,32),Dense(32,32),Linear(32,32)),
    Dense(32,2)
)

epoch = 70
parameter = paramsof(mlpmodel)
lossval  = zeros(epoch,1)
for i = 1:epoch
    ids  = shuffle!(indexs)
    outs = forward(mlpmodel, Variable(input[:,ids[1:batch]]))
    outs = softmax(outs; dims=1)
    COST = crossEntropyCost(outs, Variable(label[:,ids[1:batch]]))
    backward()
    update!(parameter,lrate)
    zerograds!(parameter)
    lossval[i] = COST
end
println(" sgd optimizer loss: ", lossval[end])

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
"""
    growGenerator(gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

Creates a CGP tree using the \"grow\" method.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG)`: random number generator for random number consistency
    along an experiment.

# User Arguments
None

# Returns
An instance of type `CGPGenotype` created by the \"grow\" method.
"""
function growGenerator(gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)
    functionSet = gpExperimentInfo._functionSet
    terminalSet = deepcopy(gpExperimentInfo._terminalSet)
    maxDepth = gpExperimentInfo._maxTreeDepth

    function growRecursiveMethod(maxDepthAux)
        expr = Array{Node}(undef, 0)
        prob = length(terminalSet)/(length(terminalSet) + length(functionSet))

        if maxDepthAux == 0 || rand(rng) < prob
            push!(expr, rand(rng, terminalSet))
        else
            push!(expr, rand(rng, functionSet))

            arity = getFunctionArity(expr[end])

            for i = 1:arity
                append!(expr, growRecursiveMethod(maxDepthAux-1))
            end
        end

        return expr
    end # function

    return CGPGenotype(growRecursiveMethod(maxDepth))
end # function

"""
    fullGenerator(gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

Creates a CGP tree using the \"full\" method.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG)`: random number generator for random number consistency
    along an experiment.

# User Arguments
None

# Returns
An instance of type `CGPGenotype` created by the \"full\" method.
"""
function fullGenerator(gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)
    functionSet = gpExperimentInfo._functionSet
    terminalSet = deepcopy(gpExperimentInfo._terminalSet)
    maxDepth = gpExperimentInfo._maxTreeDepth

    function fullRecursiveMethod(maxDepthAux)

        expr = Array{Node}(undef, 0)

        if maxDepthAux == 0
            push!(expr, rand(rng, terminalSet))
        else

            push!(expr, rand(rng, functionSet))

            arity = getFunctionArity(expr[end])

            for i = 1:arity
                append!(expr, fullRecursiveMethod(maxDepthAux-1))
            end
        end

        return expr
    end # function

    return CGPGenotype(fullRecursiveMethod(maxDepth))
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
function maketable1(;trialrunonly::Bool=trialrunonly)::Nothing
  #load the data we will work with
  panel::DataFrame = constructportfolios(refreshportfolios=false, trialrunonly=trialrunonly)
  table1panelD(panel)

  return nothing
end

#NOTE: We require
function table1panelD(panel::DataFrame)
  local cqports::CQuartileConstructions = leverageconstructions()
  Fgrps = Fgroups(cqports)

  fss::Vector{FactorSpecification} = Vector{FactorSpecification}()
  sizehint!(fss, length(T1_PANEL_D_ROWNAME) * length(Fgrps))

  #make the regression specifications
  @time for rname ∈ T1_PANEL_D_ROWNAME
    for grp ∈ Fgrps
      push!(fss, FactorSpecification(panel, T1_PANEL_D_RET[rname], :date, :permno,
        Ffactors=nothing,
        Fcharacteristics = [grp],
        Fcontrols = T1_PANEL_D_CONTROLS[rname],
        name=Symbol("Row$(rname)_Col$(grp)"),
        skiptimeseriesdfcheck=true)) #only doing crossectional regressions
    end
  end

  #WARNING: When BLAS paralllelization issues are solved replace with parallel=PARALLEL[]
  frs::Vector{FactorRegression} = (fs->characteristicFM(fs, parallel=false)).(fss)

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Mimi
using OptiMimi
using Test

using ForwardDiff

include("testlp.jl")

# Differentiation-free optimization
# Create quadratic component
@defcomp quad1 begin
    regions = Index()

    # The x-value of the maximum of the quadratic
    maximum = Parameter(index=[regions])

    # The x-value to evaluate the quadratic
    input = Parameter(index=[regions])

    # The y-value of the quadratic at the x-value
    value = Variable(index=[regions])

    function run_timestep(p, v, d, t)
        v.value[:] = -(p.input - p.maximum).^2
    end
end

# Prepare model
model1 = Model()
set_dimension!(model1, :time, 1)
set_dimension!(model1, :regions, 2)
add_comp!(model1, quad1)
set_param!(model1, :quad1, :maximum, [2., 10.])
set_param!(model1, :quad1, :input, [0., 0.])

objective1(model::Model) = sum(model[:quad1, :value])

optprob = problem(model1, [:quad1], [:input], [0.], [100.0], objective1)
(maxf, maxx) = solution(optprob, () -> [0., 0.])

@test maxf ≈ 0 atol=0.01
@test maxx[1] ≈ 2 atol=0.01
@test maxx[2] ≈ 10 atol=0.01

# Automatic differentiation

# Create quadratic component
@defcomp quad2 begin
    regions = Index()

    # The x-value of the maximum of the quadratic
    maximum = Parameter(index=[regions])

    # The x-value to evaluate the quadratic
    input = Parameter(index=[regions])

    # The y-value of the quadratic at the x-value
    value = Variable(index=[regions])

    function run_timestep(p, v, d, t)
        v.value[:] = -(p.input - p.maximum).^2
    end
end

# Prepare model
model2 = Model(Number)
set_dimension!(model2, :time, 1)
set_dimension!(model2, :regions, 2)
add_comp!(model2, quad2)
set_param!(model2, :quad2, :maximum, convert(Array{Number,1}, [2., 10.]))
set_param!(model2, :quad2, :input, convert(Array{Number,1}, [0., 0.]))

objective2(model::Model) = sum(model[:quad2, :value])

# Test the translation to a simple objective function
uo = unaryobjective(model2, [:quad2], [:input], objective2)
guo(xx) = ForwardDiff.gradient(uo, xx)
guos = guo([0., 0.])
@test guos[1] ≈ 4 atol=0.01
@test guos[2] ≈ 20 atol=0.01

# Optimize
optprob = problem(model2, [:quad2], [:input], [0.], [100.0], objective2)
(maxf, maxx) = solution(optprob, () -> [0., 0.])

@test maxf ≈ 0 atol=0.01
@test maxx[1] ≈ 2 atol=0.01
@test maxx[2] ≈ 10 atol=0.01

include("test_dupover.jl")


text:
using DataFrames

df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"])
n, m = size(df)
println(n,m)
#println(df)
#println(df[1,:])

for i =1:m
    println(df[1,:][i][1])
end
#df2 = DataFrame([1,"L"])
push!(df, [1,"F"])
println(df)

a = DataFrame()
for i=1:5
    a[i] = []
end
push!(a, [1, "F", 2, 3, 5])
print(a)

max_stars_repo_name: giladroyz/MWEM-project

max_issues_repo_name: giladroyz/MWEM-project

max_forks_repo_name: giladroyz/MWEM-project


text:
@testset "collect_pair" begin

# test of collect on pair rdds

sc = SparkContext(master="local")

nums1 = parallelize(sc, 1:2)
nums2 = parallelize(sc, 11:12)
pairs = cartesian(nums1, nums2)
values = collect(pairs)

@test values == [(1,11), (1,12), (2, 11), (2, 12)]

close(sc)

end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

# Reordering analysis: Feasibility, and where and what to reorder.

# ISSUE: We assume that there are NO aliases between any two arrays. We should
# use both static alias analysis and dynamic alias check to get rid of the
# assumption.

@doc """
A map from a permutation vector to its inverse.
"""
const inverse_color_map = Dict{Int, Int}(
    ROW_PERM     => ROW_INV_PERM,
    ROW_INV_PERM => ROW_PERM,
    COL_PERM     => COL_INV_PERM,
    COL_INV_PERM => COL_PERM
)

@doc """
A vertex in the inter-dependence graph. Each vertex represents an array's
row or column permutation vector. Color indicates the vector.
"""
type InterDependenceGraphVertex
    array      :: Symexpr
    row_perm   :: Bool
    color      :: Int
    neighbours :: Set{Tuple{InterDependenceGraphVertex, Bool}}

    InterDependenceGraphVertex(_array, _row_perm) = new(_array, _row_perm, NOT_PERM_YET,
                               Set{Tuple{InterDependenceGraphVertex, Bool}}())
end

@doc """
Inter-dependence graph. It separately represents rows and columns of arrays. It
will be colored starting from the seed.
"""
type InterDependenceGraph
    rows             :: Dict{Symexpr, InterDependenceGraphVertex}
    columns          :: Dict{Symexpr, InterDependenceGraphVertex}
    seed             :: Sym
    perm_restriction :: Int

    InterDependenceGraph(_seed) = new(
        Dict{Symexpr, InterDependenceGraphVertex}(),
        Dict{Symexpr, InterDependenceGraphVertex}(),
        _seed, PERM_VECTORS_ARE_INDEPENDENT)
end

@doc """
The CallSites' extra field for reordering.
"""
type ReorderingExtra
    seed                     :: Sym
    decider_ast              :: Expr
    matrix2mknob             :: Dict{Symexpr, Symbol}
    decider_bb               :: Any # BasicBlock
    decider_stmt_idx         :: Int
    bb                       :: Any # BasicBlock
    stmt_idx                 :: StatementIndex
    prev_stmt_idx            :: StatementIndex
    inter_dependence_graph   :: InterDependenceGraph

    # Scratch variables.
    live_in_before_prev_expr :: Set{Sym}
    live_in_before_expr      :: Set{Sym}

    ReorderingExtra(_seed, _decider_ast, _matrix2mknob) = new(_seed, _decider_ast,
             _matrix2mknob, nothing, 0, nothing, 0, 0, InterDependenceGraph(_seed),
             Set{Sym}(), Set{Sym}())
end

@doc """
Argument types of setindex!.

A matrix has the following access patterns: [:, j], [i, :], [:, :], and [i, j]
A vector has the following access patterns: [i], and [:].
"""
const setindex!_arg_types1 = (AbstractMatrix, Any, GlobalRef(Main, :(:)), Number)
const setindex!_arg_types2 = (AbstractMatrix, Any, Number, GlobalRef(Main, :(:)))
const setindex!_arg_types3 = (AbstractMatrix, Any, GlobalRef(Main, :(:)), GlobalRef(Main, :(:)))
const setindex!_arg_types4 = (AbstractMatrix, Any, Number, Number)
const setindex!_arg_types5 = (AbstractVector, Any, Number)
const setindex!_arg_types6 = (AbstractVector, Any, GlobalRef(Main, :(:)))

@doc """
For setindex!(A, x, I), map from the setindex! skeleton to the relation
between A and  x. If a dimension of A should not be permuted, then map to
NEVER_PERM.
"""
const setindex!2relations = Dict(
    setindex!_arg_types1 => (ROW_ROW,    NEVER_PERM   ),
    setindex!_arg_types2 => (NEVER_PERM, COLUMN_ROW   ),
    setindex!_arg_types3 => (ROW_ROW,    COLUMN_COLUMN),
    setindex!_arg_types4 => (NEVER_PERM, NEVER_PERM   ),
    setindex!_arg_types5 => (NEVER_PERM,              ),
    setindex!_arg_types6 => (ROW_ROW,                 ),
)

@doc """
Argument types of getindex.

A matrix has the following access patterns: [:, j], [i, :], [:, :], and [i, j]
A vector has the following access patterns: [i], and [:].
"""
const getindex_arg_types1  = (AbstractMatrix, GlobalRef(Main, :(:)), Number)
const getindex_arg_types2  = (AbstractMatrix, Number, GlobalRef(Main, :(:)))
const getindex_arg_types3  = (AbstractMatrix, GlobalRef(Main, :(:)), GlobalRef(Main, :(:)))
const getindex_arg_types4  = (AbstractMatrix, Number, Number)
const getindex_arg_types5  = (AbstractVector, Number)
const getindex_arg_types6  = (AbstractVector, GlobalRef(Main, :(:)))

@doc """
For getindex(A, I), map from the getindex skeleton to the relation
between A and getindex(A, I). If a dimension of A should not be permuted,
then map to NEVER_PERM.
"""
const getindex2relations = Dict(
    getindex_arg_types1  => (ROW_ROW,    NEVER_PERM   ),
    getindex_arg_types2  => (NEVER_PERM, COLUMN_ROW   ),
    getindex_arg_types3  => (ROW_ROW,    COLUMN_COLUMN),
    getindex_arg_types4  => (NEVER_PERM, NEVER_PERM   ),
    getindex_arg_types5  => (NEVER_PERM,              ),
    getindex_arg_types6  => (ROW_ROW,                 )
)

@doc """"
Build an edge between the two vertices for the inter-dependence graph.
Inverse indicates if their permutation vectors are inverse to each other.
"""
function build_edge(
    vertex1 :: InterDependenceGraphVertex,
    vertex2 :: InterDependenceGraphVertex,
    inverse :: Bool
)
    push!(vertex1.neighbours, (vertex2, inverse))
    push!(vertex2.neighbours, (vertex1, inverse))
end

@doc """
Build a vertex for the array's row or column permutation vector in the
inter-dependence graph.
"""
function build_vertex(
    array    :: Symexpr,
    row_perm :: Bool,
    graph    :: InterDependenceGraph
)
    if row_perm
        if !haskey(graph.rows, array)
            vertex = InterDependenceGraphVertex(array, row_perm)
            graph.rows[array] = vertex
            return vertex
        end
        return graph.rows[array]
    else
        if !haskey(graph.columns, array)
            vertex = InterDependenceGraphVertex(array, row_perm)
            graph.columns[array] = vertex
            return vertex
        end
        return graph.columns[array]
    end
end

@doc """
Build two vertices and and an edge between them in the inter-dependence graph.
"""
function build_vertices_and_edge(
    array1   :: Any,
    array2   :: Any,
    relation :: Int,
    graph    :: InterDependenceGraph
)
    if relation == ROW_ROW
        vertex1 = build_vertex(array1, true, graph)
        vertex2 = build_vertex(array2, true, graph)
        build_edge(vertex1, vertex2, false)
    elseif relation == COLUMN_COLUMN
        vertex1 = build_vertex(array1, false, graph)
        vertex2 = build_vertex(array2, false, graph)
        build_edge(vertex1, vertex2, false)
    elseif relation == COLUMN_ROW_INVERSE
        vertex1 = build_vertex(array1, false, graph)
        vertex2 = build_vertex(array2, true, graph)
        build_edge(vertex1, vertex2, true)
    else
        assert(relation == COLUMN_ROW)
        vertex1 = build_vertex(array1, false, graph)
        vertex2 = build_vertex(array2, true, graph)
        build_edge(vertex1, vertex2, false)
    end
end

@doc """
Build inter-dependence graph with the inter-dependence information drawn from
the AST of a function call to setindex!.
"""
function build_inter_dependence_graph_for_setindex!(
    ast       :: Any,
    arg_types :: Tuple,
    args      :: Vector,
    relations :: Tuple,
    graph     :: InterDependenceGraph,
)
    assert(arg_types[1] <: AbstractMatrix || arg_types[1] <: AbstractVector)

    A = get_symexpr(args[1])

    # Build dependence between A and the AST, which represents the result of
    # setindex!(). Since setindex!(A, x, I) returns the updated A, the relations
    # between them are simply ROW_ROW and COLUMN_COLUMN.
    build_vertices_and_edge(A, ast, ROW_ROW, graph)
    if arg_types[1] <: AbstractMatrix
        build_vertices_and_edge(A, ast, COLUMN_COLUMN, graph)
    end

    # Now build dependence between A and x.
    x = get_symexpr(args[2])
    for i in 1 : length(relations)
        relation = relations[i]
        if relation == NEVER_PERM
            # This is not actually a relation between A and x. Instead, it says
            # the first/second dimension of A (A's rows/columns) cannot be permuted.
            if i == 1
                graph.rows[A].color    = NEVER_PERM
            else
                graph.columns[A].color = NEVER_PERM
            end
        else
            assert(relation == ROW_ROW || relation == COLUMN_ROW || relation == COLUMN_COLUMN)
            build_vertices_and_edge(A, x, relation, graph)
        end
    end
end

@doc """
Build inter-dependence graph with the inter-dependence information drawn from
the AST of a function call to getindex.
"""
function build_inter_dependence_graph_for_getindex(
    ast       :: Any,
    arg_types :: Tuple,
    args      :: Vector,
    relations :: Tuple,
    graph     :: InterDependenceGraph,
)
    assert(arg_types[1] <: AbstractMatrix || arg_types[1] <: AbstractVector)

    A = get_symexpr(args[1])

    # Build dependence between A and the AST, which represents the result of
    # getindex().
    for i in 1 : length(relations)
        relation = relations[i]
        if relation == NEVER_PERM
            # This is not actually a relation between A and x. Instead, it says
            # the first/second dimension of A (A's rows/columns) cannot be permuted.
            if i == 1
                graph.rows[A].color    = NEVER_PERM
            else
                graph.columns[A].color = NEVER_PERM
            end
        else
            assert(relation == ROW_ROW || relation == COLUMN_ROW || relation == COLUMN_COLUMN)
            build_vertices_and_edge(A, ast, relation, graph)
        end
    end
end

@doc """
Build inter-dependence graph with the inter-dependence information drawn from
the AST of a function call, if the call is setindex! or getindex.
"""
function build_inter_dependence_graph_for_set_get_index(
    ast             :: Any,
    module_name     :: AbstractString,
    function_name   :: AbstractString,
    arg_types       :: Tuple,
    args            :: Vector,
    graph           :: InterDependenceGraph,
    handle_getindex :: Bool
)
    if module_name != "Main" || (function_name != "setindex!" && function_name != "getindex")
        return false
    end

    index2relations = handle_getindex ? getindex2relations : setindex!2relations
    for (skeleton, relations) in index2relations
        # We cannot use match_skeletons(arg_types, skeleton), since arg_types
        # contains a GlobalRef type, not a GlobalRef(Main, :(:)) (which is not a
        # type, but an instantiation of the type).
        matched = true
        if length(args) == length(skeleton)
            for i in 1 : length(skeleton)
                if skeleton[i] == GlobalRef(Main, :(:))
                    if args[i] != GlobalRef(Main, :(:))
                        matched = false
                        break
                    end
                else
                    if !(arg_types[i] <: skeleton[i])
                        matched = false
                        break
                    end
                end
            end
        else
            matched = false
        end

        if matched
            if handle_getindex
                build_inter_dependence_graph_for_getindex(ast, arg_types, args, relations, graph)
            else
                build_inter_dependence_graph_for_setindex!(ast, arg_types, args, relations, graph)
            end
            return true
        end
    end
    return false
end

@doc """
Build inter-dependence graph with the inter-dependence information drawn from
the AST of a function call, if the call is setindex! or getindex.
"""
function build_inter_dependence_graph_for_set_get_index(
    ast             :: Any,
    module_name     :: AbstractString,
    function_name   :: AbstractString,
    arg_types       :: Tuple,
    args            :: Vector,
    graph           :: InterDependenceGraph,
    handle_getindex :: Bool
)
    if module_name != "Main" || (function_name != "setindex!" && function_name != "getindex")
        return false
    end

    index2relations = handle_getindex ? getindex2relations : setindex!2relations
    for (skeleton, relations) in index2relations
        # We cannot use match_skeletons(arg_types, skeleton), since arg_types
        # contains a GlobalRef type, not a GlobalRef(Main, :(:)) (which is not a
        # type, but an instantiation of the type).
        matched = true
        if length(args) == length(skeleton)
            for i in 1 : length(skeleton)
                if skeleton[i] == GlobalRef(Main, :(:))
                    if args[i] != GlobalRef(Main, :(:))
                        matched = false
                        break
                    end
                else
                    if !(arg_types[i] <: skeleton[i])
                        matched = false
                        break
                    end
                end
            end
        else
            matched = false
        end

        if matched
            if handle_getindex
                build_inter_dependence_graph_for_getindex(ast, arg_types, args, relations, graph)
            else
                build_inter_dependence_graph_for_setindex!(ast, arg_types, args, relations, graph)
            end
            return true
        end
    end
    return false
end

@doc """
Get the symbol or expr representing the array. If array_index is something like
"0.1", it means the return result of the AST is a tuple, and element 1 of the tuple
is needed; in this case, we return a manually made expression
Expr(:TupleElement, ast, 1).
"""
function get_array(
    array_index :: Union{Int, AbstractString},
    ast         :: Expr,
    args        :: Vector,
)
    if typeof(array_index) <: Int
        a = (array_index == 0) ? ast :
             (typeof(args[array_index]) == SymbolNode ?
              args[array_index].name : args[array_index])
        return a
    else
        assert(array_index[1] == '0' && array_index[2] == '.')
        index = parse(Int, array_index[3:end])
        return Expr(:TupleElement, ast, index)
    end
end

@doc """
Build inter-dependence graph with the inter-dependence information drawn from
the AST of a function call.
"""
function build_inter_dependence_graph_for_call(
    ast           :: Any,
    module_name   :: AbstractString,
    function_name :: AbstractString,
    arg_types     :: Tuple,
    args          :: Vector,
    call_sites    :: CallSites
)
    # Recursively handle each argument
    for arg in args
        build_inter_dependence_graph(arg, call_sites)
    end

    if function_name == "tuple" || function_name == "apply_type"
        # TODO: for each element of a tuple, build an interdependence between
        # it and an Expr(TupleElement, result ast, index of the element)
        return nothing
    end

    # Now handle the call itself
    all_numbers, some_arrays = numbers_or_arrays(ast.typ, arg_types)
    if all_numbers || !some_arrays
        # The function call's result and arguments are all numbers, or
        # some are non-numbers (like Range{UInt64}) but not regular arrays.
        # No arrays to care about.
        return nothing
    end

    if function_name == ""
        throw(UnresolvedFunction(ast.head, args[1]))
    end

    graph = call_sites.extra.inter_dependence_graph

    # Special handling to setindex! and getindex: so far, we do not generate code
    # for individual array element access. Thus we need to indicate that the
    # index of the array element is not subject to reordering (Set the rows
    # or columns as NEVER_PERMUTE).
    # Try setindex!
    if build_inter_dependence_graph_for_set_get_index(ast, module_name, function_name, arg_types, args, graph, false)
        return nothing
    end

    # Try getindex
    if build_inter_dependence_graph_for_set_get_index(ast, module_name, function_name, arg_types, args, graph, true)
        return nothing
    end

    fd = look_for_function_description(module_name, function_name, arg_types)
    if fd == nothing
        throw(UndescribedFunction(module_name, function_name, arg_types, ast))
    end
    if !fd.distributive
        throw(NonDistributiveFunction(module_name, function_name, arg_types))
    end

    for (array_index1, array_index2, relation) in fd.IA
        array1 = get_array(array_index1, ast, args)
        array2 = get_array(array_index2, ast, args)
        build_vertices_and_edge(array1, array2, relation, graph)
    end

    return nothing
end

@doc """
Build inter-dependence graph with the inter-dependence information drawn from
the AST (of a statement or part of a statement).
"""
function build_inter_dependence_graph(
    ast        :: Any,
    call_sites :: CallSites
)
    if ast == call_sites.extra.decider_ast
        assert(call_sites.extra.decider_bb == nothing)
        assert(call_sites.extra.decider_stmt_idx == 0)
        call_sites.extra.decider_bb       = call_sites.extra.bb
        call_sites.extra.decider_stmt_idx = call_sites.extra.stmt_idx
    end

    local asttyp = typeof(ast)
    if asttyp <: Tuple || asttyp <: Array
        for i = 1:length(ast)
            build_inter_dependence_graph(ast[i], call_sites)
        end
        return nothing
    elseif asttyp == Expr
        symbol_info = call_sites.symbol_info
        head        = ast.head
        args        = ast.args
        if head == :return
            return build_inter_dependence_graph(args, call_sites)
        elseif head == :gotoifnot
            if_clause = args[1]
            return build_inter_dependence_graph(if_clause, call_sites)
        elseif head == :line
            return nothing
        elseif head == :call || head == :call1
            if type_of_ast_node(args[end], symbol_info) == FunctionKnob
                # Ignore the fknob added to the call.
                num_args = length(args) - 2
            else
                num_args = length(args) - 1
            end
            module_name, function_name = resolve_call_names(args)
            arg_types                  = ntuple(i-> type_of_ast_node(args[i+1],
                                                symbol_info), num_args)
            arguments                  = args[2 : num_args + 1]
            return build_inter_dependence_graph_for_call(ast, module_name,
                                function_name, arg_types, arguments, call_sites)
        elseif head == :(=)
            prev_stmt_idx = call_sites.extra.prev_stmt_idx
            if prev_stmt_idx > 0
                # Test if the previous and the current expr compose of a pattern
                # like this:
                #   GenSym(1) = (top(indexed_next))(GenSym(0),x,#s73)
                #   L = (top(getfield))(GenSym(1),1)
                # where GenSym(0) is a tuple of SparseMatrixCSC.
                prev_expr                = call_sites.extra.bb.statements[prev_stmt_idx].expr
                cur_expr                 = ast
                live_in_before_prev_expr = call_sites.extra.live_in_before_prev_expr
                live_in_before_expr      = call_sites.extra.live_in_before_expr
                pattern_skeleton1        = (:(=), Tuple{SparseMatrixCSC,Int}, Expr(:call, TopNode(:indexed_next), Tuple{SparseMatrixCSC, SparseMatrixCSC}, Int, Int))
                pattern_skeleton2        = (:(=), SparseMatrixCSC, Expr(:call, TopNode(:getfield), :f1, 1))
                if match_skeletons(prev_expr, pattern_skeleton1,  call_sites, live_in_before_prev_expr, prev_expr, cur_expr, false) &&
                   match_skeletons(cur_expr, pattern_skeleton2,  call_sites, live_in_before_expr, prev_expr, cur_expr, false)
                    graph       = call_sites.extra.inter_dependence_graph
                    arr         = cur_expr.args[1]          # L
                    tup         = prev_expr.args[2].args[2] # GenSym(0)
                    index       = prev_expr.args[2].args[3] # x
                    tup_element = Expr(:TupleElement, tup, index)
                    build_vertices_and_edge(arr, tup_element, ROW_ROW, graph)
                    build_vertices_and_edge(arr, tup_element, COLUMN_COLUMN, graph)
                    return nothing
                end
            end

            rhs_type = type_of_ast_node(args[2], symbol_info)
            if !(rhs_type <: Vector) && !(rhs_type <: AbstractMatrix)
                # Look into rhs to see if its sub-expressions have anything to
                # do with inter-dependence. For example, if rhs is dot(x, y),
                # although rhs type is not an array, x and y has ROW_ROW
                # inter-dependence
                return build_inter_dependence_graph(args[2], call_sites)
            end
            lhs_type = type_of_ast_node(args[1], symbol_info)
            module_name   = ""
            function_name = ":="
            arg_types     = (lhs_type, rhs_type)
            arguments     = args
            return build_inter_dependence_graph_for_call(ast, module_name,
                                function_name, arg_types, arguments, call_sites)
        else
            throw(UnhandledExpr(head, args))
        end
    elseif asttyp == SymbolNode  || asttyp == Symbol   || asttyp == GenSym ||
           asttyp == LabelNode   || asttyp == GotoNode || asttyp == LineNumberNode ||
           asttyp == ASCIIString || asttyp == LambdaStaticData ||
           asttyp <: Number      || asttyp == NewvarNode || asttyp == GlobalRef ||
           asttyp == QuoteNode
        return nothing
    else
        throw(UnknownASTDistributivity(ast, typeof(ast)))
    end
    return nothing
end

@doc """ A handler to visit_expressions(). """
build_inter_dependence_graph(ast, call_sites :: CallSites, top_level_number, is_top_level, read) =
    build_inter_dependence_graph(ast, call_sites)

@doc """
Color the inter-dependence graph, starting with given vertex, which has been
colored itself.
"""
function color_inter_dependence_graph(
    liveness :: Liveness,
    graph    :: InterDependenceGraph,
    from     :: InterDependenceGraphVertex,
    visited  :: Set{InterDependenceGraphVertex}
)
    if in(from, visited)
        return
    end
    push!(visited, from)
    for (vertex, inverse) in from.neighbours
        color = inverse ? inverse_color_map[from.color] : from.color
        if vertex.color != NOT_PERM_YET
            if vertex.color != color
                # Already colored, but in a different color. A conflict exists
                # But this can indicate that some colors are equal. For example,
                # from p = A * p, we know p.rowPermutation = A.rowPermutation=
                # A.colInversePermutation. Thus p's rowPermutation can have two
                # results (A.rowPermutation and A.colInversePermutation) but that
                # just means that A's rows and columns must be permuted in an
                # inverse relationship: P1 * A * P2, where P1 must equal P2'.
                # This is a nice constraint that compiler should tell libraries
                # so that rows and columns of A won't be reordered independently.
                if vertex.color == NEVER_PERM
                    vertex2index = sort_inter_dependence_graph_vertices(graph)
                    dprintln(1, 0, "\nReordering cannot proceed: vertex ",
                             vertex2index[vertex], " in the inter-dependence ",
                             "graph cannot be reordered, as it might have array ",
                             "element access that the current implementation ",
                             "does not support. However, its neighbour, vertex ",
                             vertex2index[from], " requires its color to be ",
                             permutation_color_to_str(color), ".")
                    dprintln(1, 0, "\nInter-dependence graph:")
                    dprintln(1, 1, liveness, graph)
                    throw(ConflictPermutation(from, vertex, color))
                else
                    dprintln(1, 0, "\nPermutation constraints found: ",
                             permutation_color_to_str(vertex.color),
                             " must equal ",
                             permutation_color_to_str(color))
                    if (vertex.color == ROW_PERM && color == COL_INV_PERM) ||
                       (vertex.color == COL_INV_PERM && color == ROW_PERM) ||
                       (vertex.color == ROW_INV_PERM && color == COL_PERM) ||
                       (vertex.color == COL_PERM && color == ROW_INV_PERM)
                            graph.perm_restriction |= PERM_VECTORS_ARE_INVERSE
                    elseif (vertex.color == ROW_PERM && color == COL_PERM) ||
                           (vertex.color == COL_PERM && color == ROW_PERM) ||
                           (vertex.color == ROW_INV_PERM && color == COL_INV_PERM) ||
                           (vertex.color == COL_INV_PERM && color == ROW_INV_PERM)
                            graph.perm_restriction |= PERM_VECTORS_ARE_SAME
                    else
                        # There is no other known case so far
                        assert(false)
                    end
                end
            end
        else
            vertex.color = color
        end
        color_inter_dependence_graph(liveness, graph, vertex, visited)
    end
end

@doc """
Color the inter-dependence graph, starting with the seed in it.
"""
function color_inter_dependence_graph(
    liveness :: Liveness,
    graph    :: InterDependenceGraph
)
    seed                      = graph.seed
    seed_rows_vertex          = graph.rows[seed]
    seed_columns_vertex       = graph.columns[seed]
    seed_rows_vertex.color    = ROW_PERM
    seed_columns_vertex.color = COL_PERM
    visited                   = Set{InterDependenceGraphVertex}()
    color_inter_dependence_graph(liveness, graph, seed_rows_vertex, visited)
    color_inter_dependence_graph(liveness, graph, seed_columns_vertex, visited)
end

@doc """
Add to the new expression a constant that indicates which permutation
vector to use for reordering
"""
function add_permutation_vector(
    new_expr :: Expr,
    perm     :: Int
)
    assert(perm == NOT_PERM_YET || perm == ROW_PERM || perm == ROW_INV_PERM ||
           perm == COL_PERM || perm == COL_INV_PERM)
    push!(new_expr.args,
          perm == NOT_PERM_YET ? GlobalRef(Sparso, :NOT_PERM_YET) :
          perm == ROW_PERM ? GlobalRef(Sparso, :ROW_PERM) :
          perm == ROW_INV_PERM ? GlobalRef(Sparso, :ROW_INV_PERM) :
          perm == COL_PERM ? GlobalRef(Sparso, :COL_PERM) :
          GlobalRef(Sparso, :COL_INV_PERM))
end

@doc """
Add to the new expression a matrix knob
"""
function add_matrix_knob(
    new_expr :: Expr,
    mknob    :: Symbol
)
    push!(new_expr.args, mknob)
end

@doc """
Add to the new expression an array for reordering or reverse reordering.
"""
function add_array(
    new_expr      :: Expr,
    A             :: Sym,
    symbol_info   :: Sym2TypeMap,
    graph         :: InterDependenceGraph,
    matrices_done :: Bool,
    matrix2mknob  :: Dict{Symexpr, Symbol}
)
    if haskey(graph.rows, A) # If no key, the array does not appear in the loop at all
        if !matrices_done && type_of_ast_node(A, symbol_info) <: AbstractSparseMatrix
            Pr = graph.rows[A].color
            Pc = graph.columns[A].color
            if Pr != NOT_PERM_YET || Pc != NOT_PERM_YET
                push!(new_expr.args, A)
                add_permutation_vector(new_expr, Pr)
                add_permutation_vector(new_expr, Pc)
                add_matrix_knob(new_expr, matrix2mknob[A])
            end
        elseif  matrices_done && type_of_ast_node(A, symbol_info) <: Vector
            Pr = graph.rows[A].color
            if Pr != NOT_PERM_YET
                push!(new_expr.args, A)
                add_permutation_vector(new_expr, Pr)
            end
        end
    end
end

@doc """
Add to the new expression for reordering all the arrays live out of the
statement that contains the reordering decider.
This function should be called twice. The first time matrices_done = false, and
the second time matrices_done = true. They will add matrices and vectors,
respectively.
"""
function add_arrays_to_reorder(
    new_expr      :: Expr,
    decider_stmt  :: Statement,
    symbol_info   :: Sym2TypeMap,
    liveness      :: Liveness,
    graph         :: InterDependenceGraph,
    FAR           :: Vector, # Vector{Symbol},
    matrices_done :: Bool,
    matrix2mknob  :: Dict{Symexpr, Symbol}
)
    live_out = LivenessAnalysis.live_out(decider_stmt, liveness)
    def      = LivenessAnalysis.def(decider_stmt, liveness)
#    use      = LivenessAnalysis.use(decider_stmt, liveness)

    # FAR (the symbols defined and used in the decider function call) have already
    # been reordered during the execution of that call. We expect a decide call
    # f happens in this way:
    #   f(...)  or  x = f(...)
    # In the secon case, the defininition of x has also been reordered.
    # Reorder all the live-out symbols except them.
    # ISSUE: how to make sure f() is the only call happens in the statement it
    # is called? What if for a statement like  x= f(...) * g(...)?
    for A in setdiff(live_out, union(FAR, def))
        add_array(new_expr, A, symbol_info, graph, matrices_done, matrix2mknob)
    end
end

@doc """
Add to the new expression for reverse reodering all the arrays live out of
from_bb and live into to_bb.
This function should be called twice. The first time matrices_done = false, and
the second time matrices_done = true. They will add matrices and vectors,
respectively.
"""
function add_arrays_to_reversely_reorder(
    new_expr      :: Expr,
    from_bb       :: BasicBlock,
    to_bb         :: BasicBlock,
    symbol_info   :: Sym2TypeMap,
    liveness      :: Liveness,
    graph         :: InterDependenceGraph,
    matrices_done :: Bool,
    matrix2mknob  :: Dict{Symexpr, Symbol}
)
    live_out = LivenessAnalysis.live_out(from_bb, liveness)
    live_in  = LivenessAnalysis.live_in(to_bb, liveness)
    for A in intersect(live_out, live_in)
        add_array(new_expr, A, symbol_info, graph, matrices_done, matrix2mknob)
    end
end

@doc """
Add to the new expression for reverse reordering all the arrays live into stmt.
This function should be called twice. The first time matrices_done = false, and
the second time matrices_done = true. They will add matrices and vectors,
respectively.
"""
function add_arrays_to_reversely_reorder(
    new_expr      :: Expr,
    stmt          :: Statement,
    symbol_info   :: Sym2TypeMap,
    liveness      :: Liveness,
    graph         :: InterDependenceGraph,
    matrices_done :: Bool,
    matrix2mknob  :: Dict{Symexpr, Symbol}
)
    live_in  = LivenessAnalysis.live_in(stmt, liveness)
    for A in live_in
        add_array(new_expr, A, symbol_info, graph, matrices_done, matrix2mknob)
    end
end

@doc """
Create reorder actions for the loop region.
"""
function create_reorder_actions(
    actions     :: Vector{Action},
    region      :: Region,
    symbol_info :: Sym2TypeMap,
    liveness    :: Liveness,
    graph       :: InterDependenceGraph,
    cfg         :: CFG,
    call_sites  :: CallSites,
    FAR         :: Vector, # Vector{Symbol},
    fknob       :: Symbol,
)
    decider_bb       = call_sites.extra.decider_bb
    decider_stmt_idx = call_sites.extra.decider_stmt_idx
    matrix2mknob     = call_sites.extra.matrix2mknob

    # Create an action that would insert new statements before the region
    # loop's head block (If it is a loop region), or at the beginning (if it is
    # a function region). There are two new statements: one is to set
    # the reordering decision maker, with the restrictioins of permutations,
    # and the other is to initialize reordering_status.
    region               = call_sites.region
    action_before_region = (typeof(region) == FunctionRegion) ?
        InsertBeforeOrAfterStatement(Vector{Statement}(), true, region.entry, 1) :
        InsertBeforeLoopHead(Vector{Statement}(), region.loop, true)
    push!(actions, action_before_region)

    restriction = graph.perm_restriction
    stmt        = Expr(:call, GlobalRef(Sparso, :set_reordering_decision_maker), fknob, restriction)
    push!(action_before_region.new_stmts,  Statement(0, stmt))

    reordering_status = new_symbol("reordering_status")
    stmt              = Expr(:(=), reordering_status,
                              Expr(:call, TopNode(:vect), false,
                                    GlobalRef(Main, :C_NULL),
                                    GlobalRef(Main, :C_NULL),
                                    GlobalRef(Main, :C_NULL),
                                    GlobalRef(Main, :C_NULL),
                                    0.0))
    push!(action_before_region.new_stmts,  Statement(0, stmt))

    # In the region, after the statement that contains the reordering decision maker,
    # insert a statement to reorder other arrays.
    inside_region_action = InsertBeforeOrAfterStatement(Vector{Statement}(),
                            false, decider_bb, decider_stmt_idx)
    push!(actions, inside_region_action)

    stmt = Expr(:call, GlobalRef(Sparso, :reordering),
                 fknob, reordering_status)
    push!(inside_region_action.new_stmts,  Statement(0, stmt))

    decider_stmt = decider_bb.statements[decider_stmt_idx]
    add_arrays_to_reorder(stmt, decider_stmt, symbol_info, liveness, graph, FAR, false, matrix2mknob)
    push!(stmt.args, QuoteNode(:__delimitor__))
    add_arrays_to_reorder(stmt, decider_stmt, symbol_info, liveness, graph, FAR, true, matrix2mknob)

    # At each region exit, insert reverse reordering of array.
    # Create statements that will delete the function knob at region exits
    if (typeof(region) == FunctionRegion)
        for pred_bb in region.exit.preds
            len = length(pred_bb.statements)
            assert(len > 0)
            last_stmt = pred_bb.statements[len]
            last_expr = last_stmt.expr
            assert(typeof(last_expr) <: Expr && last_expr.head == :return)

            action  = InsertBeforeOrAfterStatement(Vector{Statement}(), true, pred_bb, len)
            push!(actions, action)

            stmt = Expr(:call, GlobalRef(Sparso, :reverse_reordering),
                         reordering_status)
            push!(action.new_stmts,  Statement(0, stmt))

            add_arrays_to_reversely_reorder(stmt, last_stmt, symbol_info, liveness, graph, false, matrix2mknob)
            push!(stmt.args, QuoteNode(:__delimitor__))
            add_arrays_to_reversely_reorder(stmt, last_stmt, symbol_info, liveness, graph, true, matrix2mknob)
        end
    else
        for exit in region.exits
            action = InsertOnEdge(Vector{Statement}(), exit.from_bb, exit.to_bb)
            push!(actions, action)

            stmt = Expr(:call, GlobalRef(Sparso, :reverse_reordering),
                         reordering_status)
            push!(action.new_stmts,  Statement(0, stmt))

            add_arrays_to_reversely_reorder(stmt, exit.from_bb, exit.to_bb, symbol_info, liveness, graph, false, matrix2mknob)
            push!(stmt.args, QuoteNode(:__delimitor__))
            add_arrays_to_reversely_reorder(stmt, exit.from_bb, exit.to_bb, symbol_info, liveness, graph, true, matrix2mknob)
        end
    end
end

@doc """
Perform analyses for reordering. Write the intended transformation into actions.
"""
function reordering(
    actions     :: Vector{Action},
    region      :: Region,
    symbol_info :: Sym2TypeMap,
    liveness    :: Liveness,
    cfg         :: CFG,
    call_sites  :: CallSites
)
    old_extra   = call_sites.extra
    try
        fknob = call_sites.extra.reordering_decider_fknob
        if !in(fknob, call_sites.extra.function_knobs)
            # The AST that can decides reordering was probably removed due to,
            # for example, parent tree being changed.
            return actions
        end
        decider          = call_sites.extra.fknob2expr[fknob]
        FAR              = call_sites.extra.reordering_FAR
        seed             = FAR[1]
        matrix2mknob     = call_sites.extra.matrix2mknob
        call_sites.extra = ReorderingExtra(seed, decider, matrix2mknob)

        # An empty inter-dependence graph has been built. Build a row and a
        # column vertex for the seed.
        graph = call_sites.extra.inter_dependence_graph
        build_vertex(seed, true, graph)
        build_vertex(seed, false, graph)

        recursive  = false
        visit_expressions(region, cfg, call_sites, recursive, build_inter_dependence_graph)

        color_inter_dependence_graph(liveness, graph)

        dprintln(1, 0, "\nColored inter-dependence graph:")
        dprintln(1, 1, liveness, graph)

        new_actions = Vector{Action}()
        create_reorder_actions(new_actions, region, symbol_info, liveness, graph,
                               cfg, call_sites, FAR, fknob)

        dprintln(1, 0, "\nReordering actions to take:\n", new_actions)

        for action in new_actions
            push!(actions, action)
        end
    catch ex
        # In case any exception happen in the printing, try
        try
            Libc.flush_cstdio()
            flush(STDOUT)
            dprintln(1, 0, "Exception! Sparse Accelerator skips reordering the loop.")
            dprintln(1, 1, ex)
            dprintln(1, 0, "********************************************************************************")
            Libc.flush_cstdio()
            flush(STDOUT)
        catch
            # Do nothing
        end
    finally
        call_sites.extra = old_extra
        return actions
    end
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
using POMDPModelChecking
using POMDPModels
using POMDPs
# discrete MDP 1
ns = 4
na = 6
R = zeros(ns, na)
T = zeros(ns, na, ns) # sp, a, s

# from state 1
T[2,1, 1] = 1.0  # if go, move to 2
T[1,2:na,1] .= 1.0 # else stay

# from state 2
T[1,2,2] = 0.7 # if (s2, safe) end up in 1
T[3,2,2] = 0.3 # if (s2, safe) end up in 3

T[3,3,2] = 0.5
T[4,3,2] = 0.5
T[2,1,2] = 1.
T[2,4:end,2] .= 1.

# from state 3
T[3,:,3] .= 1

# from state 4
T[4,:,4] .= 1.
T[:,4,4] .= 0.
T[1,4,4] = 1.

mdp = TabularMDP(T, R, 0.95)

# product mdp

property = "F G succ & (G!fail)"
ltl2tgba(property, "discrete.hoa") # should translate to DRA
automata = hoa2rabin("discrete.hoa")

pmdp = ProductMDP(mdp, automata)

b0 = initialstate_distribution(pmdp)

# state_space = states(pmdp)
# action_space = actions(pmdp)
# γ = discount(pmdp)
# n_states(pmdp) == length(state_space)
# n_actions(pmdp) == length(action_space)
# statetype(pmdp) == ProductState{Int64, Int64}
# actiontype(pmdp) == Int64

# test_stateindexing(pmdp)
# test_transition(pmdp)

for s in states(pmdp)
    for a in actions(pmdp, s)
        d = transition(pmdp, s, a)
        tsum = 0.
        for (sp, p) in weighted_iterator(d)
            tsum += p
        end
        if !(tsum ≈ 1.0)
            println("error for s ", s, " a ", a, " sums to ", tsum, " next states ", d.vals)
        end
    end
end

s = ProductState(2, 1)
a = 3
d = transition(pmdp, s, a)
for (sp, p) in weighted_iterator(d)
    println((sp.s, sp.q), " ", p)
end
l = labels(mdp, s.s)

d = transition(pmdp.mdp, s.s, a)
for (sp, p) in weighted_iterator(d)
    println(sp, " ", p)
end

MEC = maximal_end_components(pmdp)
state_space = states(pmdp);
MEC = Set(state_space[i] for i in MEC)
us = accepting_states(pmdp)


text:
# ---
# title: 155. Min Stack
# id: problem155
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: Stack, Design
# link: <https://leetcode.com/problems/min-stack/description/>
# hidden: true
# ---
#
# Design a stack that supports push, pop, top, and retrieving the minimum
# element in constant time.
#
#   * push(x) -- Push element x onto stack.
#   * pop() -- Removes the element on top of the stack.
#   * top() -- Get the top element.
#   * getMin() -- Retrieve the minimum element in the stack.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["MinStack","push","push","push","getMin","pop","top","getMin"]
#     [[],[-2],[0],[-3],[],[],[],[]]
#
#     **Output**
#     [null,null,null,null,-3,null,0,-2]
#
#     **Explanation**
#     MinStack minStack = new MinStack();
#     minStack.push(-2);
#     minStack.push(0);
#     minStack.push(-3);
#     minStack.getMin(); // return -3
#     minStack.pop();
#     minStack.top();    // return 0
#     minStack.getMin(); // return -2
#
#
#
#
# **Constraints:**
#
#   * Methods `pop`, `top` and `getMin` operations will always be called on **non-empty** stacks.
#
#
## @lc code=start
using LeetCode
Base.@kwdef struct MinStack{V <: Number}
    nums::Vector{V} = Int[]
    min_nums::Vector{V} = Int[]
end

function push!(stack::MinStack, val::Int)::Nothing
    Base.push!(stack.nums, val)
    if isempty(stack.min_nums) || val <= stack.min_nums[end]
        Base.push!(stack.min_nums, val)
    end
    return
end

function pop!(stack::MinStack)::Nothing
    if Base.pop!(stack.nums) == stack.min_nums[end]
        Base.pop!(stack.min_nums)
    end
    return
end

top(stack::MinStack)::Int = stack.nums[end]

get_min(stack::MinStack)::Int = stack.min_nums[end]
## @lc code=end
## @lc test=start
@testset "155.min-stack.jl" begin
    @testset "155-test case 1" begin
        stack = MinStack()
        @test push!(stack, -2) === nothing
        @test push!(stack, 0) === nothing
        @test push!(stack, -3) === nothing
        @test get_min(stack) == -3
        @test pop!(stack) === nothing
        @test top(stack) == 0
        @test get_min(stack) == -2
    end
    @testset "155-test case 2" begin
        stack = MinStack()
        push!(stack, 2) === nothing
        @test push!(stack, 0) === nothing
        @test push!(stack, 3) === nothing
        @test push!(stack, 0) === nothing
        @test get_min(stack) == 0
        @test pop!(stack) === nothing
        @test get_min(stack) == 0
        @test pop!(stack) === nothing
        @test get_min(stack) == 0
        @test pop!(stack) === nothing
        @test get_min(stack) == 2
    end
end
## @lc test=end

## @lc info=start
# link: [solution to 155](https://leetcode-cn.com/problems/min-stack/solution/155-zui-xiao-zhan-ji-lu-zui-xiao-zhi-jul-nlr6/)
## @lc info=end