text:
struct MultiCol
    x
end

struct splitter end

# should we also define byrow as a structure?
function byrow(@nospecialize(f); @nospecialize(args...))
    br = [:($f, $args)]
    br[1].head = :BYROW
    br
end

function _check_ind_and_add!(outidx::Index, val)
    if !haskey(outidx, val)
        push!(outidx, val)
    end
end

# splitting a column to multiple columns
function normalize_modify!(outidx::Index, idx, @nospecialize(sel::Pair{<:ColumnIndex,
                                                                <:Pair{<:typeof(splitter),
                                                                    <:Vector{<:Union{Symbol, AbstractString}}}}))
    src, (fun, dst) = sel
    for i in 1:length(dst)
        _check_ind_and_add!(outidx, Symbol(dst[i]))
    end
    return outidx[src] => fun => MultiCol(Symbol.(dst))

end
function normalize_modify!(outidx::Index, idx, @nospecialize(sel::Pair{<:ColumnIndex,
                                                                <:splitter}
                                                                    ))
    throw(ArgumentError("for `splitter` the destinations must be specified"))

end

# col => fun => dst, the job is to create col => fun => :dst
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                                    <:Pair{<:Union{Base.Callable},
                                                        <:Union{Symbol, AbstractString}}})
                                                        )
    src, (fun, dst) = sel
    _check_ind_and_add!(outidx, Symbol(dst))
    return outidx[src] => fun => Symbol(dst)
end

# (col1, col2) => fun => dst, the job is to create (col1, col2) => fun => :dst
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:NTuple{N, ColumnIndex},
                                                    <:Pair{<:Union{Base.Callable},
                                                        <:Union{Symbol, AbstractString}}})
                                                        ) where N
    src, (fun, dst) = sel
    N < 2 && throw(ArgumentError("For multivariate functions (Tuple of column names), the number of input columns must be greater than 1"))
    _check_ind_and_add!(outidx, Symbol(dst))
    return ntuple(i->outidx[src[i]], N) => fun => Symbol(dst)
end

# col => fun, the job is to create col => fun => :colname
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                                    <:Union{Base.Callable}}))

    src, fun = sel
    return outidx[src] => fun => _names(outidx)[outidx[src]]
end

# (col1, col2) => fun, the job is to create (col1, col2) => fun => :colname
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:NTuple{N, ColumnIndex},
                                                    <:Union{Base.Callable}})) where N

    src, fun = sel
    N < 2 && throw(ArgumentError("For multivariate functions (Tuple of column names), the number of input columns must be greater than 1"))
    col1, col2 = outidx[src[1]], outidx[src[2]]
    var1, var2 = _names(outidx)[col1], _names(outidx)[col2]
    if N > 2
        nname = Symbol(funname(sel.second), "_", var1, "_", var2, "_etc")
    else
        nname = Symbol(funname(sel.second), "_", var1, "_", var2)
    end

    _check_ind_and_add!(outidx, nname)
    return ntuple(i->outidx[src[i]], length(src)) => fun => nname
end

# col => byrow
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                                    <:Vector{Expr}}))
    colsidx = outidx[sel.first]
    if sel.second[1].head == :BYROW
        # TODO needs a better name for destination
        # _check_ind_and_add!(outidx, Symbol("row_", funname(sel.second.args[1])))
        return outidx[colsidx] => sel.second[1] => _names(outidx)[colsidx]
    end
    throw(ArgumentError("only byrow is accepted when using expressions"))
end
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                                    <:Expr}))
    colsidx = outidx[sel.first]
    if sel.second.head == :BYROW
        # TODO needs a better name for destination
        # _check_ind_and_add!(outidx, Symbol("row_", funname(sel.second.args[1])))
        return outidx[colsidx] => sel.second => _names(outidx)[colsidx]
    end
    throw(ArgumentError("only byrow is accepted when using expressions"))
end
# col => byrow => dst
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                        <:Pair{<:Vector{Expr},
                                            <:Union{Symbol, AbstractString}}}))
    colsidx = outidx[sel.first]
    if sel.second.first[1].head == :BYROW
        # TODO needs a better name for destination
        _check_ind_and_add!(outidx, Symbol(sel.second.second))
        return outidx[colsidx] => sel.second.first[1] => Symbol(sel.second.second)
    end
    throw(ArgumentError("only byrow is accepted when using expressions"))
end
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                        <:Pair{<:Expr,
                                            <:Union{Symbol, AbstractString}}}))
    colsidx = outidx[sel.first]
    if sel.second.first.head == :BYROW
        # TODO needs a better name for destination
        _check_ind_and_add!(outidx, Symbol(sel.second.second))
        return outidx[colsidx] => sel.second.first => Symbol(sel.second.second)
    end
    throw(ArgumentError("only byrow is accepted when using expressions"))
end

# cols => fun, the job is to create [col1 => fun => :col1name, col2 => fun => :col2name ...]
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:MultiColumnIndex,
                                                    <:Vector{Expr}}))
    colsidx = outidx[sel.first]
    if sel.second isa AbstractVector && sel.second[1] isa Expr
        if sel.second[1].head == :BYROW
            # TODO needs a better name for destination
            _check_ind_and_add!(outidx, Symbol("row_", funname(sel.second[1].args[1])))
            return outidx[colsidx] => sel.second[1] => Symbol("row_", funname(sel.second[1].args[1]))
        end
    end
    # res = Any[normalize_modify!(outidx, idx, colsidx[1] => sel.second)]
    # for i in 2:length(colsidx)
    #     push!(res, normalize_modify!(outidx, idx, colsidx[i] => sel.second))
    # end
    # return res
end
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:MultiColumnIndex,
                                                    <:Union{Base.Callable, Expr}}))
    colsidx = outidx[sel.first]
    if sel.second isa Expr
        if sel.second.head == :BYROW
            # TODO needs a better name for destination
            _check_ind_and_add!(outidx, Symbol("row_", funname(sel.second.args[1])))
            return outidx[colsidx] => sel.second => Symbol("row_", funname(sel.second.args[1]))
        end
    end
    res = Any[normalize_modify!(outidx, idx, colsidx[1] => sel.second)]
    for i in 2:length(colsidx)
        push!(res, normalize_modify!(outidx, idx, colsidx[i] => sel.second))
    end
    return res
end
# cols => funs which will be normalize as col1=>fun1, col2=>fun2, ...
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:MultiColumnIndex,
                                                    <:Vector{<:Base.Callable}}))
    colsidx = outidx[sel.first]
    if !(length(colsidx) == length(sel.second))
        throw(ArgumentError("The input number of columns and the length of the number of functions should match"))
    end
    res = Any[normalize_modify!(outidx, idx, colsidx[1] => sel.second[1])]
    for i in 2:length(colsidx)
        push!(res, normalize_modify!(outidx, idx, colsidx[i] => sel.second[i]))
    end
    return res
end

function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:ColumnIndex,
                                                    <:Vector{<:Base.Callable}}))
    colsidx = outidx[sel.first]
    normalize_modify!(outidx, idx, colsidx .=> sel.second[i])
    return res
end

# special case cols => byrow(...) => :name
function normalize_modify!(outidx::Index, idx,
    @nospecialize(sel::Pair{<:MultiColumnIndex,
                            <:Pair{<:Vector{Expr},
                                <:Union{Symbol, AbstractString}}}))
    colsidx = outidx[sel.first]
    if sel.second.first[1].head == :BYROW
        _check_ind_and_add!(outidx, Symbol(sel.second.second))
        return outidx[colsidx] => sel.second.first[1] => Symbol(sel.second.second)
    else
        throw(ArgumentError("only byrow operation is supported for cols => fun => :name"))
    end
end
function normalize_modify!(outidx::Index, idx,
    @nospecialize(sel::Pair{<:MultiColumnIndex,
                            <:Pair{<:Expr,
                                <:Union{Symbol, AbstractString}}}))
    colsidx = outidx[sel.first]
    if sel.second.first.head == :BYROW
        _check_ind_and_add!(outidx, Symbol(sel.second.second))
        return outidx[colsidx] => sel.second.first => Symbol(sel.second.second)
    else
        throw(ArgumentError("only byrow operation is supported for cols => fun => :name"))
    end
end

# cols .=> fun .=> dsts, the job is to create col1 => fun => :dst1, col2 => fun => :dst2, ...
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:MultiColumnIndex,
                                                    <:Pair{<:Union{Base.Callable,Vector{Expr}},
                                                        <:AbstractVector{<:Union{Symbol, AbstractString}}}}))
    colsidx = outidx[sel.first]
    if !(length(colsidx) == length(sel.second.second))
        throw(ArgumentError("The input number of columns and the length of the output names should match"))
    end
    res = Any[normalize_modify!(outidx, idx, colsidx[1] => sel.second.first => sel.second.second[1])]
    for i in 2:length(colsidx)
        push!(res, normalize_modify!(outidx, idx, colsidx[i] => sel.second.first => sel.second.second[i]))
    end
    return res
end
# cols .=> fun .=> dsts, the job is to create col1 => fun => :dst1, col2 => fun => :dst2, ...
function normalize_modify!(outidx::Index, idx,
                            @nospecialize(sel::Pair{<:MultiColumnIndex,
                                                    <:Pair{<:Expr,
                                                        <:AbstractVector{<:Union{Symbol, AbstractString}}}}))
    colsidx = outidx[sel.first]
    if !(length(colsidx) == length(sel.second.second))
        throw(ArgumentError("The input number of columns and the length of the output names should match"))
    end
    res = Any[normalize_modify!(outidx, idx, colsidx[1] => sel.second.first => sel.second.second)]
    for i in 2:length(colsidx)
        push!(res, normalize_modify!(outidx, idx, colsidx[i] => sel.second.first => sel.second.second[i]))
    end
    return res
end

function normalize_modify!(outidx::Index, idx, arg::AbstractVector)
    res = Any[]
    for i in 1:length(arg)
        _res = normalize_modify!(outidx::Index, idx::Index, arg[i])
        if _res isa AbstractVector
            for j in 1:length(_res)
                push!(res, _res[j])
            end
        else
            push!(res, _res)
        end
    end
    return res
end

function normalize_modify_multiple!(outidx::Index, idx, @nospecialize(args...))
    res = Any[]
    for i in 1:length(args)
        _res = normalize_modify!(outidx, idx, args[i])
        if typeof(_res) <: Pair
            push!(res, _res)
        else
            for j in 1:length(_res)
                push!(res, _res[j])
            end
        end
    end
    res
end

modify(origninal_ds::AbstractDataset, @nospecialize(args...)) = modify!(copy(origninal_ds), args...)

modify!(ds::Dataset) = parent(ds)
function modify!(ds::AbstractDataset, @nospecialize(args...))
    if ds isa SubDataset
		idx_cpy = copy(index(parent(ds)))
	else
        idx_cpy = Index(copy(index(ds).lookup), copy(index(ds).names), copy(index(ds).format))
    end
    if isgrouped(ds)
        norm_var = normalize_modify_multiple!(idx_cpy, index(ds), args...)
        allnewvars = map(x -> x.second.second, norm_var)
        all_new_var = Symbol[]
        for i in 1:length(allnewvars)
            if typeof(allnewvars[i]) <: MultiCol
                for j in 1:length(allnewvars[i].x)
                    push!(all_new_var, allnewvars[i].x[j])
                end
            else
                push!(all_new_var, allnewvars[i])
            end
        end
        var_index = idx_cpy[unique(all_new_var)]
        any(index(ds).sortedcols .∈ Ref(var_index)) && throw(ArgumentError("the grouping variables cannot be modified, first use `ungroup!(ds)` to ungroup the data set"))
        _modify_grouped(ds, norm_var)
    else
        _modify(ds, normalize_modify_multiple!(idx_cpy, index(ds), args...))
    end
end

# size() is better option for checking if the result is scalar,
# it works for numbers and it won't work for strings and symbols
_is_scalar(::Missing, sz) = true
_is_scalar(::Tuple, sz) = true
function _is_scalar(_res, sz)
     resize_col = false
    try
        size(_res)
        if size(_res) == () || size(_res,1) != sz
            # fill!(Tables.allocatecolumn(typeof(_res), nrow(ds)),
            #                           _res)
            # _res = repeat([_res], nrow(ds))
            resize_col = true
        end
    catch e
        if (e isa MethodError)
            # fill!(Tables.allocatecolumn(typeof(_res), nrow(ds)),
                                      # _res)
            # _res = repeat([_res], nrow(ds))
            resize_col = true
       end

    end
    return resize_col
end

function _resize_result!(ds, _res, newcol)
    resize_col = _is_scalar(_res, nrow(ds))
    if resize_col
        if ds isa SubDataset
            if haskey(index(ds), newcol)
                ds[:, newcol] = fill!(Tables.allocatecolumn(typeof(_res), nrow(ds)), _res)
            elseif !haskey(index(parent(ds)), newcol)
                parent(ds)[!, newcol] = missings(typeof(_res), nrow(parent(ds)))
                _update_subindex!(index(ds), index(parent(ds)), newcol)
                ds[:, newcol] = fill!(Tables.allocatecolumn(typeof(_res), nrow(ds)), _res)
            else
                throw(ArgumentError("modifing a parent's column which doesn't appear in SubDataset is not allowed"))
            end
        else
            ds[!, newcol] = fill!(Tables.allocatecolumn(typeof(_res), nrow(ds)), _res)
        end
    else
        if ds isa SubDataset
            if haskey(index(ds), newcol)
                ds[:, newcol] = _res
            elseif !haskey(index(parent(ds)), newcol)
                parent(ds)[!, newcol] = missings(eltype(_res), nrow(parent(ds)))
                _update_subindex!(index(ds), index(parent(ds)), newcol)
                ds[:, newcol] = _res
            else
                throw(ArgumentError("modifing a parent's column which doesn't appear in SubDataset is not allowed"))
            end
        else
            ds[!, newcol] = _res
        end
    end
end

function _modify_single_var!(ds, _f, x, dst)
    _res = _f(x)
    _resize_result!(ds, _res, dst)
end
function _modify_single_tuple_var!(ds, _f, x, dst)
    _res = _f(x...)
    _resize_result!(ds, _res, dst)
end

# the number of destination can be smaller or greater than the number of elements of Tuple,
function _modify_multiple_out!(ds, x, dst)
    !(nonmissingtype(eltype(x)) <: Tuple) && throw(ArgumentError("to use `splitter`, the source column must be a vector of Tuple"))
    tb = Tables.columntable(x)
    for j in 1:length(dst)
        try
            _resize_result!(ds, Tables.getcolumn(tb, j), dst[j])
        catch
            _resize_result!(ds, missings(nrow(ds)), dst[j])
        end
    end
end

function _modify_f_barrier(ds, msfirst, mssecond, mslast)
    if (mssecond isa Base.Callable) && !(mslast isa MultiCol)
        if msfirst isa NTuple
            _modify_single_tuple_var!(ds, mssecond, ntuple(i -> _columns(ds)[msfirst[i]], length(msfirst)), mslast)
        else
            _modify_single_var!(ds, mssecond, _columns(ds)[msfirst], mslast)
        end
    elseif (mssecond isa Expr) && mssecond.head == :BYROW
        try
            if ds isa SubDataset
                _res = byrow(ds, mssecond.args[1], msfirst; mssecond.args[2]...)
                if haskey(index(ds), mslast)
                    ds[:, mslast] = _res
                elseif !haskey(index(parent(ds)), mslast)
                    parent(ds)[!, mslast] = missings(eltype(_res), nrow(parent(ds)))
                    _update_subindex!(index(ds), index(parent(ds)), mslast)
                    ds[:, mslast] = _res
                else
                    throw(ArgumentError("modifing a parent's column which doesn't appear in SubDataset is not allowed"))
                end
            else
                ds[!, mslast] = byrow(ds, mssecond.args[1], msfirst; mssecond.args[2]...)
            end
        catch e
            if e isa MethodError
                throw(ArgumentError("There is problem in your `byrow`, make sure that the output of `byrow` is a vector"))
            end
            rethrow(e)
        end
    elseif  (mssecond isa Base.Callable) && (mslast isa MultiCol) && (mssecond isa typeof(splitter))
        _modify_multiple_out!(ds, _columns(ds)[msfirst], mslast.x)
    else
        @error "not yet know how to handle this situation $(msfirst => mssecond => mslast)"
    end
end

function _modify(ds, ms)
    needs_reset_grouping = false
    for i in 1:length(ms)
        _modify_f_barrier(ds, ms[i].first, ms[i].second.first, ms[i].second.second)
    end
    return ds
end

function _check_the_output_type(ds::Dataset, ms)
    if ms.first isa Tuple
        CT = return_type(ms.second.first, ntuple(i -> _columns(ds)[ms.first[i]], length(ms.first)))
    else
        CT = return_type(ms.second.first, _columns(ds)[ms.first])
    end
    # TODO check other possibilities:
    # the result can be
    # * AbstractVector{T} where T
    # * Vector{T}
    # * not a Vector
    CT == Union{} && throw(ArgumentError("compiler cannot assess the return type of calling `$(ms.second.first)` on `:$(_names(ds)[ms.first])`, you may want to try using `byrow`"))
    if CT <: AbstractVector
        if hasproperty(CT, :var)
            T = Union{Missing, CT.var.ub}
        else
            T = Union{Missing, eltype(CT)}
        end
    else
        T = Union{Missing, CT}
    end
    T
end

# FIXME notyet complete
# fill _res for grouped data: col => f => :newcol
function _modify_grouped_fill_one_col!(_res, x, _f, starts, ngroups, nrows)
    Threads.@threads for g in 1:ngroups
        lo = starts[g]
        g == ngroups ? hi = nrows : hi = starts[g + 1] - 1
        _tmp_res = _f(view(x, lo:hi))
        resize_col = _is_scalar(_tmp_res, length(lo:hi))
        if resize_col
            fill!(view(_res, lo:hi), _tmp_res)
        else
            copy!(view(_res, lo:hi), _tmp_res)
        end
    end
    _res
end

function _modify_grouped_fill_one_col_tuple!(_res, x, _f, starts, ngroups, nrows)
    Threads.@threads for g in 1:ngroups
        lo = starts[g]
        g == ngroups ? hi = nrows : hi = starts[g + 1] - 1
        _tmp_res = do_call(_f, x, lo:hi)
        resize_col = _is_scalar(_tmp_res, length(lo:hi))
        if resize_col
            fill!(view(_res, lo:hi), _tmp_res)
        else
            copy!(view(_res, lo:hi), _tmp_res)
        end
    end
    _res
end

function _modify_grouped_f_barrier(ds, msfirst, mssecond, mslast)
    if (mssecond isa Base.Callable) && !(mslast isa MultiCol)
        T = _check_the_output_type(ds, msfirst=>mssecond=>mslast)
        _res = Tables.allocatecolumn(T, nrow(ds))
        if msfirst isa Tuple
            _modify_grouped_fill_one_col_tuple!(_res, ntuple(i->_columns(ds)[msfirst[i]], length(msfirst)), mssecond, index(ds).starts, index(ds).ngroups[], nrow(ds))
        else
            _modify_grouped_fill_one_col!(_res, _columns(ds)[msfirst], mssecond, index(ds).starts, index(ds).ngroups[], nrow(ds))
        end
        ds[!, mslast] = _res
    elseif (mssecond isa Expr)  && mssecond.head == :BYROW
        ds[!, mslast] = byrow(ds, mssecond.args[1], msfirst; mssecond.args[2]...)
    elseif (mssecond isa Base.Callable) && (mslast isa MultiCol) && (mssecond isa typeof(splitter))
        _modify_multiple_out!(ds, _columns(ds)[msfirst], mslast.x)
    else
                # if something ends here, we should implement new functionality for it
        @error "not yet know how to handle the situation $(msfirst => mssecond => mslast)"
    end
end

function _modify_grouped(ds, ms)
    needs_reset_grouping = false
    for i in 1:length(ms)
        _modify_grouped_f_barrier(ds, ms[i].first, ms[i].second.first, ms[i].second.second)
    end
    return parent(ds)
end


text:
using PSplines
using Documenter

DocMeta.setdocmeta!(PSplines, :DocTestSetup, :(using PSplines); recursive=true)

makedocs(;
    modules=[PSplines],
    authors="Víctor Amores",
    repo="https://github.com/victorjesusamoresmedianero/PSplines.jl/blob/{commit}{path}#{line}",
    sitename="PSplines.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://victorjesusamoresmedianero.github.io/PSplines.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/victorjesusamoresmedianero/PSplines.jl",
    devbranch = "main",
)

ext: jl


text:
module Calendar

import ICU
using Compat

export CalendarTime,
       format,
       now,
       today,
       parse_date,
       ymd,
       ymd_hms,
       yDhms,

       # tests
       isAM,
       isPM,
       isleapyear,

       # fields
       year,
       month,
       week,
       dayofyear,
       dayofweek,
       day,
       hour,
       hour12,
       minute,
       second,
       tz,
       timezone,

       # deprecated
       am, pm,
       isleap,

       # durations
       CalendarDuration,
       FixedCalendarDuration,
       years,
       months,
       weeks,
       days,
       hours,
       minutes,
       seconds,
       hms,
       hm,

       # ranges
       CalendarTimeRange,

       #Constants
       Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday,
       January, February, March, April, May, June, July, August, September, October, November, December

import Base.show, Base.(+), Base.(-), Base.(<), Base.(==), Base.length,
       Base.colon, Base.getindex, Base.start, Base.next, Base.done, Base.(*), Base.(.*),
       Base.size, Base.step, Base.vcat, Base.isless, Base.hash, Base.isequal

immutable CalendarTime
    millis::Float64
    cal::ICU.ICUCalendar
end

# needed by DataFrames
length(::CalendarTime) = 1

# support hashing
hash(ct::CalendarTime) = bitmix(hash(ct.millis),hash(ct.cal))
isequal(a::CalendarTime, b::CalendarTime) = a.millis == b.millis && a.cal == b.cal

# default timezone
_tz = ICU.getDefaultTimeZone()

const _cal_cache = Dict()
const _tz_cache = Dict()
function _get_cal(tz)
    if !haskey(_cal_cache, tz)
        cal = ICU.ICUCalendar(tz)
        _cal_cache[tz] = cal
        _tz_cache[cal] = tz
    end
    _cal_cache[tz]
end

const _format_cache = Dict()
function _get_format(cal)
    tz = _tz_cache[cal]
    if !haskey(_format_cache, tz)
        _format_cache[tz] = ICU.ICUDateFormat(ICU.UDAT_LONG, ICU.UDAT_MEDIUM, tz)
    end
    _format_cache[tz]
end
function _get_format(pattern, cal)
    tz = _tz_cache[cal]
    k = pattern,tz
    if !haskey(_format_cache, k)
        _format_cache[k] = ICU.ICUDateFormat(pattern, tz)
    end
    _format_cache[k]
end

now(tz=_tz) = CalendarTime(ICU.getNow(), _get_cal(tz))

function today()
    n = now()
    ymd(year(n), month(n), day(n))
end

function ymd_hms(y::Integer, mo::Integer, d::Integer, h::Integer, mi::Integer, s::Real, tz=_tz)
    cal = _get_cal(tz)
    is = trunc(Integer,s)
    ms = rem(s,1)*1e3
    ICU.clear(cal)
    ICU.setDateTime(cal, y, mo, d, h, mi, is)
    CalendarTime(ICU.getMillis(cal) + ms, cal)
end

function ymd_hms(s::AbstractString, tz=_tz)
    re = r"(\d{4})[^\d]+(\d{1,2})[^\d]+(\d{1,2})[^\d]+(\d{1,2})[^\d]+(\d{2})[^\d]+(\d{2})"
    m = match(re, s)
    if m === nothing; return nothing end
    y  = int(m.captures[1])
    mo = int(m.captures[2])
    d  = int(m.captures[3])
    h  = int(m.captures[4])
    mi = int(m.captures[5])
    s  = int(m.captures[6])
    ymd_hms(y, mo, d, h, mi, s, tz)
end

function yDhms(y::Integer, jd::Integer, h::Integer, mi::Integer, s::Real, tz=_tz)
    cal = _get_cal(tz)
    is = trunc(Integer,s)
    ms = rem(s,1)*1e3
    ICU.clear(cal)
    ICU.set(cal, ICU.UCAL_YEAR, y)
    ICU.set(cal, ICU.UCAL_DAY_OF_YEAR, jd)
    ICU.set(cal, ICU.UCAL_HOUR, h)
    ICU.set(cal, ICU.UCAL_MINUTE, mi)
    ICU.set(cal, ICU.UCAL_SECOND, is)
    CalendarTime(ICU.getMillis(cal) + ms, cal)
end

ymd(y::Integer, m::Integer, d::Integer, tz=_tz) = ymd_hms(y, m, d, 0, 0, 0, tz)

tz(t::CalendarTime) = _tz_cache[t.cal]
tz(t::CalendarTime, tz) = CalendarTime(t.millis, _get_cal(tz))
const timezone = tz

for (f,k,o) in [(:year,ICU.UCAL_YEAR,0),
                (:month,ICU.UCAL_MONTH,1),
                (:week,ICU.UCAL_WEEK_OF_YEAR,0),
                (:dayofyear,ICU.UCAL_DAY_OF_YEAR,0),
                (:dayofweek,ICU.UCAL_DAY_OF_WEEK,0),
                (:day,ICU.UCAL_DATE,0),
                (:hour12,ICU.UCAL_HOUR,0),
                (:hour,ICU.UCAL_HOUR_OF_DAY,0),
                (:minute,ICU.UCAL_MINUTE,0),
                (:second,ICU.UCAL_SECOND,0)]

    @eval begin
        function ($f)(t::CalendarTime)
            ICU.setMillis(t.cal, t.millis)
            ICU.get(t.cal, $k) + $o
        end

        function ($f)(t::CalendarTime, val::Integer)
            ICU.setMillis(t.cal, t.millis)
            ICU.set(t.cal, $k, val - $o)
            CalendarTime(ICU.getMillis(t.cal), t.cal)
        end

        @vectorize_1arg CalendarTime $f
    end
end

isleapyear(t::CalendarTime) = isleapyear(year(t))
isleapyear(y::Integer) = (((y % 4 == 0) && (y % 100 != 0)) || (y % 400 == 0))
@vectorize_1arg CalendarTime isleapyear
@deprecate isleap isleapyear

function isPM(t::CalendarTime)
    ICU.setMillis(t.cal, t.millis)
    ICU.get(t.cal, ICU.UCAL_AM_PM) == 1
end
isAM(t::CalendarTime) = !isPM(t)
@vectorize_1arg CalendarTime isAM
@vectorize_1arg CalendarTime isPM

@deprecate am isAM
@deprecate pm isPM

for op in [:<, :(==), :isless]
    @eval ($op)(t1::CalendarTime, t2::CalendarTime) = ($op)(t1.millis, t2.millis)
end

function format(pattern::AbstractString, t::CalendarTime)
    utf8(ICU.format(_get_format(pattern,t.cal), t.millis))
end
format(t::CalendarTime, pattern::AbstractString) = format(pattern, t)

function parse_date(pattern::AbstractString, s::AbstractString, tz::AbstractString)
    try
        cal = _get_cal(tz)
        millis = ICU.parse(_get_format(pattern,cal), s)
        return CalendarTime(millis, cal)
    catch
        error("failed to parse '", s, "' with '", pattern, "'")
    end
end
parse_date(pattern, s) = parse_date(pattern, s, _tz)
parse_date{S<:AbstractString}(pattern::AbstractString, s::AbstractArray{S}, tz::AbstractString) = map(x -> parse_date(pattern, x, tz), s)
parse_date{S<:AbstractString}(pattern::AbstractString, s::AbstractArray{S}) = map(x -> parse_date(pattern, x, _tz), s)
const parse = parse_date

function show(io::IO, t::CalendarTime)
    s = ICU.format(_get_format(t.cal), t.millis)
    print(io, s)
end

abstract AbstractCalendarDuration

immutable CalendarDuration <: AbstractCalendarDuration
    years::Int
    months::Int
    weeks::Int
    millis::Float64
end

immutable FixedCalendarDuration <: AbstractCalendarDuration
    millis::Float64
end

CalendarDuration() = CalendarDuration(0,0,0,0.)

years(x::Integer)  = CalendarDuration(x, 0, 0, 0.)
months(x::Integer) = CalendarDuration(0, x, 0, 0.)
weeks(x::Integer)  = CalendarDuration(0, 0, x, 0.)
for f in [:years,:months,:weeks]
    @eval @vectorize_1arg Integer $f
end
for (f,a) in [(:days,86400e3),(:hours,3600e3),(:minutes,60e3),(:seconds,1e3)]
    @eval ($f)(x::Real) = FixedCalendarDuration(x*($a))
    @eval @vectorize_1arg Real $f
end
hms(h::Real, m::Real, s::Real) = FixedCalendarDuration(1e3(60(60h + m) + s))
hm(h::Real, m::Real) = hms(h, m, 0)

function print_millis(io::IO, millis::Float64, first::Bool)
    negative = millis < 0
    secs = abs(millis*1e-3)
    for (x,tag) in [(86400," day"),(3600," hour"), (60," minute")]
        if secs >= x
            n = floor(Integer,secs/x)
            secs -= n*x

            if first
                if negative write(io,'-') end
                first = false
            else
                write(io, negative ? " - " : " + ")
            end
            print(io, n, tag)
            if n > 1 write(io,'s') end
        end
    end
    if secs > 0
        if first
            if negative write(io,'-') end
            first = false
        else
            write(io, negative ? " - " : " + ")
        end
        print(io, isinteger(secs) ? round(Int,secs) : secs, " second")
        if secs != 1. write(io,'s') end
    end
end

function show(io::IO, d::CalendarDuration)
    first = true
    for f in [:years,:months,:weeks]
        n = d.(f)
        if n != 0
            negative = n < 0
            if first
                if negative write(io,'-') end
                first = false
            else
                write(io, negative ? " - " : " + ")
            end
            print(io, abs(n))
            write(io, ' ')
            sf = string(f)
            write(io, abs(n) > 1 ? sf : sf[1:end-1])
        end
    end
    if d.millis != 0.
        print_millis(io, d.millis, first)
    end
end

function show(io::IO, d::FixedCalendarDuration)
    print_millis(io, d.millis, true)
end

for op in [:+, :-]
    @eval begin
        ($op)(d1::CalendarDuration, d2::CalendarDuration) =
            CalendarDuration($op(d1.years, d2.years),
                             $op(d1.months, d2.months),
                             $op(d1.weeks, d2.weeks),
                             $op(d1.millis, d2.millis))

        ($op)(d::CalendarDuration, f::FixedCalendarDuration) =
            CalendarDuration(d.years,
                             d.months,
                             d.weeks,
                             $op(d.millis, f.millis))

        ($op)(f::FixedCalendarDuration, d::CalendarDuration) =
            CalendarDuration($op(d.years),
                             $op(d.months),
                             $op(d.weeks),
                             $op(f.millis, d.millis))

        ($op)(d1::FixedCalendarDuration, d2::FixedCalendarDuration) =
            FixedCalendarDuration($op(d1.millis, d2.millis))

        function ($op)(t::CalendarTime, d::CalendarDuration)
            ICU.setMillis(t.cal, $op(t.millis, d.millis))
            for (f,v) in [(:years,ICU.UCAL_YEAR),
                          (:months,ICU.UCAL_MONTH),
                          (:weeks,ICU.UCAL_WEEK_OF_YEAR)]
                n = d.(f)
                if n != 0
                    ICU.add(t.cal, v, $op(n))
                end
            end
            CalendarTime(ICU.getMillis(t.cal), t.cal)
        end

        ($op)(t::CalendarTime, d::FixedCalendarDuration) =
            CalendarTime(($op)(t.millis, d.millis), t.cal)

        ($op)(d::AbstractCalendarDuration, t::CalendarTime) = ($op)(t, d)

        @vectorize_2arg @compat(Union{CalendarTime, AbstractCalendarDuration}) $op
        # Fix up cases that @vectorize doesn't cover (without this, the result is of type Array{Any}):
        # I can't figure out how to consolidate this.
        ($op)(x::Array{CalendarTime}, y::Array{CalendarDuration}     ) = reshape(CalendarTime[ ($op)(x[i], y[i]) for i=1:length(x) ], promote_shape(size(x),size(y)))
        ($op)(x::Array{CalendarTime}, y::Array{FixedCalendarDuration}) = reshape(CalendarTime[ ($op)(x[i], y[i]) for i=1:length(x) ], promote_shape(size(x),size(y)))
        ($op)(y::Array{CalendarDuration}     , x::Array{CalendarTime}) = reshape(CalendarTime[ ($op)(x[i], y[i]) for i=1:length(x) ], promote_shape(size(x),size(y)))
        ($op)(y::Array{FixedCalendarDuration}, x::Array{CalendarTime}) = reshape(CalendarTime[ ($op)(x[i], y[i]) for i=1:length(x) ], promote_shape(size(x),size(y)))
    end
end

for op in [:*, :.*]
    @eval begin
        ($op)(d::CalendarDuration, i::Integer) = CalendarDuration(d.years*i, d.months*i, d.weeks*i, d.millis*i)
        ($op)(d::FixedCalendarDuration, x::Real) = FixedCalendarDuration(d.millis*x)
        ($op){T<:AbstractCalendarDuration}(d::Array{T}, i::Integer) =
            reshape([ d[j] * i for j=1:length(d) ], size(d))
        ($op){I<:Integer}(d::AbstractCalendarDuration, i::Array{I}) =
            reshape([ d * i[j] for j=1:length(i) ], size(i))
        ($op)(x, d::AbstractCalendarDuration) = ($op)(d, x)
    end
end

(-)(t1::CalendarTime, t2::CalendarTime) = FixedCalendarDuration(t1.millis - t2.millis)
@vectorize_2arg CalendarTime (-)

(-)(d::CalendarDuration) = CalendarDuration(-d.years, -d.months, -d.weeks, -d.millis)
(-)(d::FixedCalendarDuration) = FixedCalendarDuration(-d.millis)
@vectorize_1arg AbstractCalendarDuration (-)

for op in [:<, :(==)]
    @eval begin
        ($op)(d1::FixedCalendarDuration, d2::FixedCalendarDuration) =
            ($op)(d1.millis, d2.millis)
    end
end

immutable CalendarTimeRange{T<:AbstractCalendarDuration} <: Range{CalendarTime}
    start::CalendarTime
    step::T
    len::Int
end

function show{T}(io::IO, r::CalendarTimeRange{T})
    print(io, r.start, ':', r.step, ':', r.start + (r.len-1)*r.step)
end

function colon(t1::CalendarTime, d::FixedCalendarDuration, t2::CalendarTime)
    n = floor(Integer, (t2.millis - t1.millis)/d.millis) + 1
    CalendarTimeRange(t1, d, n)
end

function colon(t1::CalendarTime, d::CalendarDuration, t2::CalendarTime)
    approx_d = (365.2425*d.years + 30.436875*d.months + 7.*d.weeks)*86400e3 + d.millis
    n = floor(Integer, (t2.millis - t1.millis)/approx_d) + 1
    CalendarTimeRange(t1, d, n)
end

colon(t1::CalendarTime, t2::CalendarTime) = colon(t1, seconds(1), t2)

function getindex(r::CalendarTimeRange, i::Integer)
    if !(1 <= i <= r.len); error(BoundsError); end
    r.start + (i-1)*r.step
end

size(r::CalendarTimeRange) = (r.len,)
length(r::CalendarTimeRange) = r.len
step(r::CalendarTimeRange) = r.step
start(r::CalendarTimeRange) = 0
next{T}(r::CalendarTimeRange{T}, i) = (r.start + i*r.step, i+1)
done(r::CalendarTimeRange, i) = length(r) <= i

function vcat(r::CalendarTimeRange)
    n = length(r)
    a = Array(CalendarTime,n)
    i = 1
    for x in r
        a[i] = x
        i += 1
    end
    return a
end

convert(::Type{Array{CalendarTime,1}}, r::CalendarTimeRange) = vcat(r)

const Sunday = 1
const Monday = 2
const Tuesday = 3
const Wednesday = 4
const Thursday = 5
const Friday = 6
const Saturday = 7

const January = 1
const February = 2
const March = 3
const April = 4
const May = 5
const June = 6
const July = 7
const August = 8
const September = 9
const October = 10
const November = 11
const December = 12

end # module


text:
"""

gradcheck(f, w, x...; gcheck=10, o...)

Numerically checks the gradient of `f(w,x...;o...)` with respect to
its first argument `w`, which could be a Number, Array, Tuple or Dict.
`gcheck` largest entries from each numeric array in gradient
`dw=(grad(f))(w,x...;o...)` is compared to its numerical estimate and
a Warning is issued for comparisons that fail.

"""
function gradcheck(f, w, x...; gcheck=10, o...)
    g = grad(f)
    d = g(w, x...; o...)
    if isa(w, Number)
        gc_number(d, f, w, x...; o...)
    elseif isa(w, KnetArray) || (isa(w, Array) && isbits(eltype(w)))
        gc_array(w, d, f, w, x...; gcheck=gcheck, o...)
    else
        k = gc_indices(w)
        for i in k
            gc_index(w, d, i, f, w, x...; gcheck=gcheck, o...)
        end
    end
end

function gc_index(w, d, i, f, w0, x...; gcheck=10, o...)
    if isa(w[i], Number)
        gc_array(w, d, f, w0, x...; gcheck=1, icheck=i, o...)
    elseif isa(w[i],KnetArray) || (isa(w[i], Array) && isbits(eltype(w[i])))
        gc_array(w[i], d[i], f, w0, x...; gcheck=gcheck, o...)
    else
        k = gc_indices(w[i])
        for j in k
            gc_index(w[i], d[i], j, f, w0, x...; gcheck=gcheck, o...)
        end
    end
end

function gc_array(w, d, f, worig, x...; gcheck=10, icheck=0, o...)
    irange = (icheck > 0 ? (icheck:icheck) :
              length(w) <= gcheck ? (1:length(w)) :
              d == nothing ? rand(1:length(w), gcheck) :
              sortperm(abs(vec(Array(d))),rev=true)[1:gcheck])
    (delta, atol, rtol) = gc_params(typeof(w[first(irange)]))
    for i in irange
        w0 = w[i]
        (w1, w2) = gc_interval(w0, delta)
        w[i] = w1
        f1 = f(worig, x...; o...)
        w[i] = w2
        f2 = f(worig, x...; o...)
        w[i] = w0
        nd = (f2-f1) / (w2-w1)
        di = (d==nothing ? 0 : d[i])
        if !isapprox(di, nd; rtol=rtol, atol=atol)
            warn("d=$di nd=$nd")
        else
            println("gcheck: d=$di nd=$nd")
        end
    end
end

function gc_number(d, f, w, x...; o...)
    (delta, atol, rtol) = gc_params(typeof(w))
    (w1, w2) = gc_interval(w, delta)
    (f1, f2) = (f(w1,x...;o...), f(w2,x...;o...))
    nd = (f2-f1) / (w2-w1)
    if !isapprox(d, nd; rtol=rtol, atol=atol)
        warn("d=$d nd=$nd")
    else
        println("gcheck: d=$d nd=$nd")
    end
end

gc_params(t)=(a=cbrt(eps(t)); (a,a,a))
gc_indices(w)=eachindex(w)
gc_indices(w::Tuple)=(1:length(w))

function gc_interval(w,d)
    w1=w-d/2
    w2=w+d/2
    (w1 < 0 < w) && (w1=zero(w))
    (w2 > 0 > w) && (w2=zero(w))
    return (w1,w2)
end


text:
module CategoricalTimeSeries

include("IntegerIB.jl")
include("SerialDependence.jl")
include("SpectralEnvelope.jl")
include("MotifRecognition.jl")

export IB, search_optima!, brute_optimize!, IB_optimize!, calc_metrics, get_IB_curve, get_y, print_results
export cramer_coefficient, cohen_coefficient, conditional_entropy, LaggedBivariateProbability, H, theils_u, rate_evolution, bootstrap_CI
export spectral_envelope, get_mappings, detrend, smooth, power_spectrum, apply_mappings
export detect_motifs, mapping, find_motifs, plot_motif, expected_matches, least_occurence_threshold

end


text:
"""
ilqr.jl
"""

struct iLQRSolver{IR,Tm,To,Tix,Tiu,Txx,Tuu,Tuud,Tux,Tuxd,Txu,Tx,Tu,TG,T
                  } <: UnconstrainedSolver{T}
    # problem info
    n::Int
    m::Int
    N::Int
    model::Tm
    obj::To
    ix::Tix
    iu::Tiu
    # TODO: mem bottleneck
    # state at each time step: n X N
    X::Vector{Tx}
    X_tmp::Vector{Tx}
    U::Vector{Tu}
    U_tmp::Vector{Tu}
    ts::Vector{T}
    # gains
    # TODO: mem bottleneck
    # state feedback gains (m, n) x N
    K::Vector{Tux}
    K_dense::Tuxd
    # feedforward gains (m) x N
    d::Vector{Tu}
    # discrete dynamics jacobians
    # w.r.t. state (n, n)
    A::Txx
    # w.r.t. control (n, m)
    B::Txu
    # state difference jacobian (n̄, n)
    G::TG
    # quadratic expansion of obj
    E::QuadraticCost{Txx,Tuu,Tux,Tx,Tu,T}
    Qxx::Txx
    Qxx_tmp::Txx
    Quu::Tuu
    Quu_dense::Tuud
    Quu_reg::Tuu
    Qux::Tux
    Qux_tmp::Tux
    Qux_reg::Tux
    Qx::Tx
    Qu::Tu
    # cost-to-go
    P::Txx
    P_tmp::Tux
    p::Tx
    p_tmp::Tu
    ΔV::Vector{T}
    # regularization
    ρ::Vector{T}
    # regularization rate of change
    dρ::Vector{T}
    # gradient
    grad::Vector{T}
    # misc
    opts::SolverOptions{T}
    stats::SolverStats{T}
    logger::SolverLogger
    solver_name::Symbol
end

function iLQRSolver(prob::Problem{IR,T}, opts::SolverOptions, stats::SolverStats) where {IR,T}
    n, m, N = size(prob)
    n̄ = state_diff_size(prob.model)
    M = prob.M
    Md = prob.Md
    V = prob.V
    K = [M(zeros(T, m, n̄)) for k = 1:N-1]
    K_dense = Md(zeros(T, m, n̄))
    d = [V(zeros(T, m)) for k = 1:N-1]
    A = M(zeros(T, n, n))
    B = M(zeros(T, n, m))
    G = M(zeros(T, n̄, m))
    Qxx = M(zeros(T, n, n))
    Qxx_tmp = copy(Qxx)
    Quu = M(zeros(T, m, m))
    Quu_dense = Md(zeros(T, m, m))
    Quu_reg = copy(Quu)
    Qux = M(zeros(T, m, n))
    Qux_tmp = copy(Qux)
    Qux_reg = copy(Qux)
    Qx = V(zeros(T, n))
    Qu = V(zeros(T, m))
    P = M(zeros(T, n, n))
    P_tmp = M(zeros(T, m, n))
    p = V(zeros(T, n))
    p_tmp = V(zeros(T, m))
    ΔV = zeros(T, 2)
    ρ = zeros(T, 1)
    dρ = zeros(T, 1)
    grad = zeros(T, N-1)
    logger = SolverLogging.default_logger(opts.verbose >= 2)
    Tm = typeof(prob.model)
    To = typeof(prob.obj)
    Tix = typeof(prob.ix)
    Tiu = typeof(prob.iu)
    Txx = typeof(Qxx)
    Tuu = typeof(Quu)
    Tuud = typeof(Quu_dense)
    Tux = typeof(Qux)
    Tuxd = typeof(K_dense)
    Txu = typeof(B)
    Tx = typeof(Qx)
    Tu = typeof(Qu)
    TG = typeof(G)
    solver = iLQRSolver{IR,Tm,To,Tix,Tiu,Txx,Tuu,Tuud,Tux,Tuxd,Txu,Tx,Tu,TG,T}(
        prob.n, prob.m, prob.N, prob.model, prob.obj, prob.ix, prob.iu, prob.X,
        prob.X_tmp, prob.U, prob.U_tmp, prob.ts,
        K, K_dense, d, A, B, G, prob.E, Qxx, Qxx_tmp, Quu, Quu_dense,
        Quu_reg, Qux, Qux_tmp, Qux_reg, Qx, Qu, P, P_tmp,
        p, p_tmp, ΔV, ρ, dρ, grad, opts, stats, logger, :iLQR)
    reset!(solver)
    return solver
end

# methods
Base.size(solver::iLQRSolver) = solver.n, solver.m, solver.N

log_level(::iLQRSolver) = InnerLoop

function reset!(solver::iLQRSolver{T}) where T
    # reset regularization
    solver.ρ[1] = 0.0
    solver.dρ[1] = 0.0
    return nothing
end


text:
using Fire
using ReadCoverage

"""
    relcov(path_bam::String, path_bed12::String, out::String; max_depth::Int=0)

Calculates relative gene body coverage from a BAM file.
(Almost same as genebody_coverage.py in RSeQC.)
"""
@main function relcov(path_bam::String, path_bed12::String, out::String; max_depth::Int=0)
    relative_genebodycoverage(path_bam, path_bed12, output_prefix=out, max_depth=max_depth);
end

"""
    abcov(path_bam::String, path_bed12::String, out::String; bin_size::Int=100)

Cacluates absolute gene body coverage from a BAM file
"""
@main function abcov(path_bam::String, path_bed12::String, out::String; bin_size::Int=100)
    absolute_genebodycoverage(path_bam, path_bed12, output_prefix=out, bin_size=bin_size);
end

"""
    coverage(path_bam::String, chrom::String, leftpos::Int64, rightpos::Int64, out::String)

Cacluates absolute gene body coverage from a BAM file
"""
@main function coverage(path_bam::String, chrom::String, leftpos::Int64, rightpos::Int64, out::String)
    readcoverage_bam(path_bam, chrom, leftpos, rightpos, output_prefix=out);
end


text:
function b(t,n,k)
    out = binomial(n,k).*(t.^k).*((1 .- t).^(n-k))
    return out
end


text:
@test_throws ErrorException deploydocs(
    repo = "github.com/JuliaDocs/Documenter.jl.git",
    julia = 0.5,        # must be a string
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using OptiMo, Bazinga, ScSTO
using Printf
using PyPlot, PyCall
using DelimitedFiles

foldername = "/home/albertodm/Documents/ScSTO.jl/demo/"

################################################################################
################################################################################
# unconstrained STO
################################################################################
################################################################################
# time interval
t0 = 0.0
tf = 12.0

# control sequence [nu x N]
uvec = Matrix([repeat([0.0; 1.0], 5, 1); 0.0]')
nu = size(uvec, 1)
N = size(uvec, 2)

# cost matrix
C = Matrix([1.0 0.0 -1.0 0.0; 0.0 1.0 0.0 -1.0])
Q = Matrix(C' * C)

# initial state
x0 = [0.5; 0.7; 1; 1] # nx x 1
nx = length(x0)

# system dynamics
# nx, nu  ->  nx
# nx, nu  -> nx x nx
function dynam(x::Vector{Float64}, u::Vector{Float64})
    n = length(x)
    f = zeros(n)
    f[1] = x[1] - x[1] * x[2] - 0.4 * x[1] * u[1]
    f[2] = -x[2] + x[1] * x[2] - 0.2 * x[2] * u[1]
    return f
end
function d_dynam(x::Vector{Float64}, u::Vector{Float64})
    n = length(x)
    df = zeros(n, n)
    df[1, 1] = 1.0 - x[2] - 0.4 * u[1]
    df[1, 2] = -x[1]
    df[2, 1] = x[2]
    df[2, 2] = -1.0 + x[1] - 0.2 * u[1]
    return df
end

# number of points on the fixed grid
ngrid = 100

# switching cost
swc_grid = [0.0; 0.1; 1.0]
col_grid = [:steelblue, :orangered, :green]
ng = length(swc_grid)

# solver
solver = Bazinga.ZEROFPR( max_iter = 50, verbose = true )

# allocation
nt = 1000
objective = Array{Float64}(undef, ng, 2)
cputime = Array{Float64}(undef, ng, 2)
swdelta = Array{Float64}(undef, N, ng, 2)
swtau = Array{Float64}(undef, N - 1, ng, 2)
swtf = Array{Float64}(undef, ng, 2)
xsim = Array{Float64}(undef, nx, nt, ng, 2)
usim = Array{Float64}(undef, nu, nt, ng, 2)
repo_objf = Array{Any}(undef, ng, 2)

for k = 1:ng
    swc = swc_grid[k]

    # STO problem with switching costs
    prob = scstoproblem(
        x0,
        dynam,
        d_dynam,
        uvec,
        ngrid = ngrid,
        t0 = t0,
        tf = tf,
        Q = Q,
        swc = swc,
    )

    # solve!
    out = solver(prob)
    print(out)

    cputime[k, 1] = out.time
    swdelta[:, k, 1] .= out.x
    swtau[:, k, 1], swtf[k, 1] = gettau(prob, swdelta[:, k, 1])

    xsim[:, :, k, 1], _, objective[k, 1], t = simulate(prob, swtau[:, k, 1], length = nt)
    usim[:, :, k, 1], _ = simulateinput(prob, swtau[:, k, 1], t)

    repo_objf[k, 1] = prob.repo.objf

end

t = collect(range(t0, stop = tf, length = nt))

tmpmat = Array{Float64}(undef, nt, 1+3*ng)
tmpmat[:,1] = t
for k = 1:ng
    tmpmat[:,3*k-1] = xsim[1, :, k, 1]
    tmpmat[:,3*k]   = xsim[2, :, k, 1]
    tmpmat[:,3*k+1] = usim[1, :, k, 1]
end
writedlm( foldername * "fishing_unc_traj.csv",  tmpmat, ',')

#=
figure()
subplot(3, 1, 1)
for k = 1:ng
    plot(t, xsim[1, :, k, 1], c = col_grid[k])
end
ylim(0, 1.75)
xlim(0, 12)
yticks([0; 1])
ylabel(L"x_1")
subplot(3, 1, 2)
for k = 1:ng
    plot(t, xsim[2, :, k, 1], c = col_grid[k])
end
ylabel(L"x_2")
ylim(0, 1.75)
xlim(0, 12)
yticks([0; 1])
subplot(3, 1, 3)
for k = 1:ng
    plot(t, usim[1, :, k, 1], label = "σ = $(swc_grid[k])", c = col_grid[k])
end
legend()
ylim(-0.2, 1.2)
yticks([0; 1])
xlim(0, 12)
ylabel(L"u")
xlabel(L"$\mathrm{Time}\; [s]$")
gcf()
savefig(foldername * "ScSTO.jl/demo/data/fishing_unc_traj.pdf")

objfmin = minimum(repo_objf[1, 1])
for k = 2:ng
    global objfmin = min(objfmin, minimum(repo_objf[k, 1]))
end

figure()
for k = 1:ng
    semilogy(repo_objf[k, 1] .- objfmin)
end
gcf()
savefig(foldername * "ScSTO.jl/demo/data/fishing_unc_objf.pdf")

figure()
for k in 1:ng
    semilogy(repo_objf[k,1] .- minimum(repo_objf[k,1]))
end
gcf()

swnnz = similar(swc_grid)
for k in 1:ng
    swnnz[k] = sum(swdelta[:,k,1] .> 0)
end
swobjective = swc_grid .* swnnz
swobjective .+= objective
=#

################################################################################
################################################################################
# constrained STO
################################################################################
################################################################################
# constraints
# N-1  ->  ncon
# N-1, ncon  ->  N-1
# ncon  ->  ncon
ncon = 3
function constr(tau::Vector{Float64}, c::Vector{Float64})
    c[1] = tau[1]
    c[2] = tau[3]
    c[3] = tau[5]
    return nothing
end
function d_constr(tau::Vector{Float64}, v::Vector{Float64}, jtv::Vector{Float64})
    jtv .= 0.0
    jtv[1] = v[1]
    jtv[3] = v[2]
    jtv[5] = v[3]
    return nothing
end
function p_constr(c::Vector{Float64}, p::Vector{Float64})
    p[1] = max(1.0, min(c[1], 2.0))
    p[2] = max(4.0, min(c[2], 5.0))
    p[3] = max(7.0, min(c[3], 8.0))
    return nothing
end

solver2 = Bazinga.ALPX(
    max_iter = 10,
    max_sub_iter = 10,
    subsolver = :zerofpr,
    verbose = true,
)

for k = 1:ng
    swc = swc_grid[k]

    prob2 = scstoproblem(
        x0,
        dynam,
        d_dynam,
        uvec,
        ngrid = ngrid,
        t0 = t0,
        tf = tf,
        Q = Q,
        swc = swc,
        ncon = ncon,
        constr = constr,
        dconstr = d_constr,
        pconstr = p_constr,
    )

    warmstart!(prob2, swtau[:, k, 1])

    out2 = solver2(prob2)
    print(out2)

    swdelta[:, k, 2] = out2.x
    swtau[:, k, 2], _ = gettau(prob2, swdelta[:, k, 2])
    xsim[:, :, k, 2], _, _, _ = simulate(prob2, swtau[:, k, 2], t)
    usim[:, :, k, 2], _ = simulateinput(prob2, swtau[:, k, 2], t)
end

tmpmat = Array{Float64}(undef, nt, 1+2*3)
k = 2
tmpmat[:,1] = t
tmpmat[:,2] = xsim[1, :, k, 1]
tmpmat[:,3] = xsim[2, :, k, 1]
tmpmat[:,4] = usim[1, :, k, 1]
tmpmat[:,5] = xsim[1, :, k, 2]
tmpmat[:,6] = xsim[2, :, k, 2]
tmpmat[:,7] = usim[1, :, k, 2]
writedlm( foldername * "fishing_con_traj.csv",  tmpmat, ',')

#=
figure()
subplot(3, 1, 1)
k = 2
plot(t, xsim[1, :, k, 1], ls = :dashed)
plot(t, xsim[1, :, k, 2])
ylim(0, 1.75)
xlim(0, 12)
yticks([0; 1])
ylabel(L"x_1")
subplot(3, 1, 2)
plot(t, xsim[2, :, k, 1], ls = :dashed)
plot(t, xsim[2, :, k, 2])
ylabel(L"x_2")
ylim(0, 1.75)
xlim(0, 12)
yticks([0; 1])
subplot(3, 1, 3)
plot(t, usim[1, :, k, 1], ls = :dashed, label = "σ = $(swc_grid[k]) : unc")
plot(t, usim[1, :, k, 2], label = "σ = $(swc_grid[k]) : con")
legend() # loc = "upper right", ncol = 2)
ylim(-0.2, 1.2)
xlim(0, 12)
yticks([0; 1])
ylabel(L"u")
xlabel(L"$\mathrm{Time}\; [s]$")
gcf()
savefig(foldername * "fishing_con_traj.pdf")
=#


text:
# WIP
import Base: ==, hash

"""
struct Scope

The `Scope` tracks a typevar and its 'rules' or 'constraints'.  For example,
it represents `:(T where T <: AbstractFloat)`.
"""
struct Scope
    var::Symbol
    rule::Union{Nothing, Expr}
end
Scope(var::Symbol) = Scope(var, nothing)

as_scope_expr(s::Scope) = s.rule === nothing ? s.var : s.rule

"""
function rule_symbols(s::Scope)

From a [`Scope`](@ref) extract all typevar symbols
"""
function rule_symbols(s::Scope)
    extract_symbols(e::Symbol) = [e]
    extract_symbols(e::Expr) = extract_symbols(Val(e.head), e)
    extract_symbols(::Val{:(.)}, e) = extract_symbols(strip_module(e))
    extract_symbols(::Val{:curly}, e) = mapreduce(extract_symbols,append!, e.args[2:end];init=Symbol[])
    extract_symbols(::Val{:where}, e) = extract_symbols(e.args[1])
    function extract_symbols(::Union{Val{:(>:)},Val{:(<:)}}, e)
        if length(e.args) == 1
            extract_symbols(e.args[1])
        else
            if s.var === e.args[1]
                extract_symbols(e.args[2])
            else
                [extract_symbols(e.args[1]); extract_symbols(e.args[2])...]
            end
        end
    end

    s.rule === nothing && return Symbol[]
    extract_symbols(s.rule)
end

"""
function replace_rule_symbol!(expr::Expr, s, t)

In the given expression, replace all references to the symbol `s` with the
symbol `t`.
"""
function replace_rule_symbol!(expr::Expr, s, t)
    rrs(e::Symbol) = e === s ? t : e
    rrs(e::Expr) = (rrs(Val(e.head), e); e)
    rrs(::Val{:where}, e) = (rrs(e.args[1]); e)
    function rrs(::Val{:curly}, e)
        for (i,a) in enumerate(e.args[2:end])
            e.args[i+1] = rrs(a)
        end
        e
    end
    function rrs(::Union{Val{:(>:)},Val{:(<:)}}, e)
        e.args[1] = rrs(e.args[1])
        if length(e.args) > 1
            e.args[2] = rrs(e.args[2])
        end
        e
    end

    rrs(expr)
end

Base.hash(obj::Scope, h::UInt) = hash(obj.var)
==(o1::Scope, o2::Scope) = o1.var === o2.var
==(o1::Scope, o2::Symbol) = o1.var === o2
==(o1::Symbol, o2::Scope) = o2.var === o1

"""
struct Dependencies

`Dependencies` tracks the set of implementation's expressions that are bound to a
given interface's [`Scope`](@ref), and the dependencies on those bindings, so
every `Scope` has a related `Dependencies` struct (perhaps empty).

The `bindings` field contains the set of implementation's expressions dependent on
a `Scope`, while the `dependencies` field maps the specific interface's scope to the
implementation's signature's expression.
"""
struct Dependencies
    bindings::Vector
    dependencies::Vector{<:Pair}

    Dependencies() = new(Any[], Pair[])
end

"""
    Signature

A structure that holds information about the signature of a method; it includes the
expression for every argument and the return value, along with the evaluated types
for each, and includes the type parameters following the method (the `where` clause)
"""
struct Signature
    mod::Module
    args::Vector{<:SE}
    args_evald::Vector
    retval::SE
    retval_evald
    wparams

    function Signature(mod::Module, args::Vector{<:SE}, retval::SE, wparams)
        newargs,newretval = add_wparams_to_args(wparams, args, retval)

        # @debug "Signature" args=args retval=retval wparams=wparams

        # here we evaluate the args, and then recreate the Expressions
        # for each arg; the idea is to normalize type expressions
        args_evald = [mod.eval(a) for a in newargs]
        newargs_from_types = map(zip(args_evald, args)) do (arg,oldarg)
            e = strip_module(Meta.parse(string(arg)))
            replace_vector_with_array(e)
            e
        end

        newarg_exprs = map(zip(args_evald, args, newargs, newargs_from_types)) do (t,old,a,b)
            e = decide_between_exprs(t,a,b,wparams)
            e = strip_added_where(old, e, wparams)
        end

        # @debug "Signature" newargs=newargs newargs_from_types=newargs_from_types evald=args_evald newarg_exprs=newarg_exprs

        retval_evald = mod.eval(newretval)
        newretval_from_type = strip_module(Meta.parse(string(retval_evald)))
        replace_vector_with_array(newretval_from_type)

        newretval_expr = decide_between_exprs(retval_evald, newretval, newretval_from_type, wparams)

        newretval_expr = strip_module(newretval_expr)
        newretval_expr = strip_added_where(retval, newretval_expr, wparams)

        new(mod,
            SE[na for na in newarg_exprs],
            args_evald,
            newretval_expr,
            retval_evald,
            wparams)
    end
end

function add_wparams_to_args(wparams, args, retval)
    if isempty(wparams)
        newargs = args
        newretval = retval
    else
        newargs = map(args) do arg
            # wparams will include all typevars, but it works
            Expr(:where, [arg; wparams...]...)
        end
        newretval = Expr(:where, [retval; wparams...]...)
    end
    (newargs, newretval)
end

function strip_added_where(old, new::Symbol, wparams)
    # @info "strip_added_where(Any,Symbol)" old=old new=new
    new
end

function strip_added_where(old::Symbol, new::Expr, wparams)
    if new.head === :where
        try_replace_symbol!(new, wparams)
    end
    new
end

#=
    old = :Score
    new = :(Score{I} where I)
    transform to :(Score{I} where I)

    old = :(Score{<:O})
    new = :(Score{var"#s24"} where {O, var"#s24" <: O})
    transform to :(Score{var"#s24"} where {var"#s24" <: O})
=#
function strip_added_where(old::Expr, new::Expr, wparams)
    inwhere(s::Symbol) = Any[s]
    inwhere(e::Expr) = inwhere(Val(e.head), e)
    inwhere(::Union{Val{:(<:)}, Val{:(>:)}}, e) = mapreduce(arg->inwhere(arg), append!, e.args; init=Any[])
    inwhere(::Union{Val{:where},Val{:curly}}, e) = mapreduce(arg->inwhere(arg), append!, e.args[2:end]; init=Any[])

    incurly(x) = Any[]
    incurly(s::Symbol) = Any[s]
    incurly(e::Expr) = incurly(Val(e.head), e)
    incurly(e::QuoteNode) = incurly(e.value)
    incurly(::Val{:(.)}, e) = incurly(e.args[2])
    incurly(::Val{:curly}, e) = mapreduce(arg->incurly(arg), append!, e.args[2:end]; init=Any[])
    incurly(::Union{Val{:(<:)}, Val{:(>:)}}, e) = incurly(e.args[1])

    isrelated(::Val, e, tokeep) = false
    function isrelated(::Union{Val{:(<:)}, Val{:(>:)}}, e, tokeep)
        isrelated(e.args[1], tokeep) || (length(e.args) > 1 && isrelated(e.args[2], tokeep))
    end
    isrelated(::Val{:curly}, e, tokeep) = any(arg->isrelated(arg,tokeep), e.args[2:end])
    isrelated(s::Symbol, tokeep) = s in tokeep
    isrelated(e, tokeep) = isrelated(Val(e.head), e, tokeep)

    extract_related(s::Symbol) = s
    extract_related(e::Expr) = extract_related(Val(e.head), e)
    extract_related(::Union{Val{:(<:)}, Val{:(>:)}}, e::Expr) = e.args[1]

    old == new && return new

    inwhere = old.head === :where ? inwhere(Val(old.head), old) : Any[]
    incurly = old.head === :curly ? incurly(Val(old.head), old) : Any[]

    if new.head === :where
        toremove = Int64[]
        for (i,arg) in enumerate(new.args[2:end])
            if !(isrelated(arg, inwhere)) && (arg in incurly)
                push!(toremove, i+1)
            end
        end
        deleteat!(new.args, toremove)
        if length(new.args) == 1
            new = new.args[1]
        else
            try_replace_symbol!(new, wparams)
        end
    end

    new
end

function try_replace_symbol!(new, wparams) # new is a :where Expr
    # We skip replacing :(T where T) constructs; we
    # only handle :(X{I,O} where O) constructs
    new.args[1] isa Expr || return
    new.args[1].head === :curly || return

    shouldreplace(arg, p::Expr) = shouldreplace(arg, p.args[1])
    shouldreplace(arg, p::Symbol) = arg === p

    function replace_symbol(::Val, e::Expr, arg, sym) # :curly, :(<:), :(>:)
        for i in 1:length(e.args)
            e.args[i] = replace_symbol(e.args[i], arg, sym)
        end
        e
    end
    function replace_symbol(::Val{:curly}, e::Expr, arg, sym)
        for i in 2:length(e.args)
            e.args[i] = replace_symbol(e.args[i], arg, sym)
        end
        e
    end
    replace_symbol(e, arg, sym) = e
    replace_symbol(e::Symbol, arg, sym) = e === arg ? sym : e
    replace_symbol(e::Expr, arg, sym) = replace_symbol(Val(e.head), e, arg, sym)

    for p in wparams
        for i in 2:length(new.args)
            arg = new.args[i]
            if shouldreplace(arg, p)
                replace_symbol(new, arg, gensym())
            end
        end
    end
end

#= We prefer to use the parse(string(eval(expr))) version (the b parameter),
   since it tends to be the most complete.  In certain cases,
   though, (e.g. (T where T)) it removes our typevars and
   convert it to `Any`.  In the first case, see below where the `C` is
   not specified.

   ```
   struct Foo{A,B,C} end

   foo(f::Foo{A,B}) where {A,B} = (A,B)
   ```

   Once we eval and reparse the type, we get the `C` back.

   We have to handle all special cases, such as Vectors, which
   are aliases, and 'decompile' differently depending on whether
   there is a typevar:  e.g. Vector{T} where T => Array{T,1} where T
   while Vector{Int64} => Vector{Int64}.  We convert all Vector
   references to Array references.
=#
function decide_between_exprs(t, a, b, wparams)
    # this handles the case where we have Q where Q <: Tuple, so b == :Tuple but
    # we want to keep Q where Q <: Tuple
    inwparams = in(a isa Expr ? a.args[1] : a, map(x->x isa Expr ? x.args[1] : x, wparams))

    # @debug "mapping newargs" a=a b=b t=t inwparams=inwparams
    if t == Vector
        b
    elseif t <: Vector
        x = copy(b)
        replace_vector_with_array(x)
        x
    elseif t == Any || inwparams
        a
    else
        b
    end

end

replace_vector_with_array(curr) = nothing

function replace_vector_with_array(expr::Expr)
    for (i,arg) in enumerate(expr.args)
        if arg === :Vector
            expr.args[i] = :Array
            if i == 1
                push!(expr.args, 1)
            end
        else
            replace_vector_with_array(arg)
        end
    end
end

body(t::UnionAll) = body(t.body)
body(t::Type) = t

extract_vars(t::Type, a::Vector) = nothing
function extract_vars(t::UnionAll, a::Vector)
    push!(a, t.var.name)
    extract_vars(t.body, a)
end

function extract_vars(::Val{:where}, t, a)
    append!(a, t.args[2:end])
    extract_vars(t.args[1], a)
end

function extract_vars(::Union{Val{:(<:)}, Val{:(>:)}}, t, a)
    for arg in t.args[1:end]
        extract_vars(arg, a)
    end
end

function extract_vars(::Val{:curly}, t, a)
    for arg in t.args[2:end]
        extract_vars(arg, a)
    end
end

extract_vars(t::Expr, a::Vector) = extract_vars(Val(t.head), t, a)
extract_vars(t::Symbol, a::Vector) = begin
    if startswith(string(t), "#")
        push!(a, t)
    end
end
extract_vars(::Number, a::Vector) = nothing
extract_vars(t) = (a = []; extract_vars(t, a); a)

function is_number(s)
    rx = r"^\d+(?:\.\d*)?$"
    match(rx, s) !== nothing
end

function wrap_in_where(expr, wparams)
    isempty(wparams) ? expr : Expr(:where, [expr;unique(wparams)...]...)
end

wheres(scope, sig) = unique(mapreduce(d->[as_scope_expr(x) for x in keys(d)], append!, scope; init=Any[sig.wparams...]))

function eval_expr(mod, expr, params)
    is_number(string(expr)) && return mod.eval(expr)
    if expr isa Expr && length(expr.args) == 1 && (expr.head == :(<:) || expr.head == :(>:))
        s = gensym()
        pushfirst!(expr.args, s)
        whered = wrap_in_where(expr, [s])
    end
    whered = wrap_in_where(expr, params)

    @debug "eval_expr" whered=whered

    evald = try
        mod.eval(whered)
    catch e
        if e isa UndefVarError
            @debug "Handling eval $whered; returning Any"
            Any
        else
            rethrow()
        end
    end

    return evald
end

function eval_in_union(mod, i_expr, s_expr, i_params, s_params, errors)
    i_eval = eval_expr(mod, i_expr, i_params)
    s_eval = eval_expr(mod, s_expr, s_params)

    @debug "eval_in_union" i_eval=i_eval s_eval=s_eval
    if !(s_eval isa Core.TypeofVararg || i_eval isa Core.TypeofVararg) && !(s_eval <: i_eval)
        push!(errors, """Invalid Union:
            \tdefined :  $(s_eval)
            \texpected:  $(i_eval)
            \tdefined :  $(s_expr)
            \texpected:  $(i_expr)""")
        return (false, Nothing, i_expr, s_expr)
    end

    # new_iexpr = Meta.parse(string(i_eval))
    # replace_vector_with_array(new_iexpr)
    # new_sexpr = Meta.parse(string(s_eval))
    # replace_vector_with_array(new_sexpr)

    return (true, s_eval, i_expr, s_expr)
end

ntharg(sig::Signature, n) = n > length(sig.args_evald) ? sig.retval_evald : sig.args_evald[n]

# we added the first where clause to every argument
stripwheres(s::Symbol, wheres) = s
function stripwheres(e::Expr, wheres)
    gf(a) = a
    gf(a::Expr) = a.args[1]

    todel = Int64[]
    e.head === :where || return e
    c = copy(e)
    for (i,arg) in enumerate(c.args[2:end])
        for w in wheres
            if gf(arg) == gf(w)
                push!(todel, i+1)
            end
        end
    end

    deleteat!(c.args, todel)
    length(c.args) == 1 ? c.args[1] : c
end

find_scope(::Val, expr, v) = (expr, v)

function find_scope(::Val{:where}, expr, v)
    push!(v, expr.args[2:end]...)
    find_scope(expr.args[1], v)
end

find_scope(expr) = (expr, [])
find_scope(expr::Expr) = find_scope(Val(expr.head), expr, [])

find_scope(expr, v) = (expr, v)
find_scope(expr::Expr, v) = find_scope(Val(expr.head), expr, v)

to_scope(::Val{:(>:)}, expr) = Scope(expr.args[1], expr)
to_scope(::Val{:(<:)}, expr) = Scope(expr.args[1], expr)
to_scope(expr::Symbol) = Scope(expr)
to_scope(expr::Expr) = to_scope(Val(expr.head), expr)

strip_module(expr::Expr, args1::Expr) = begin
    expr.args[1] = strip_module(Val(args1.head), args1)
    expr
end
strip_module(expr::Expr, args1::Symbol) = expr

strip_module(::Union{Val{:curly}, Val{:where}}, expr) = strip_module(expr, expr.args[1])
strip_module(::Val{:(.)}, expr) = strip_module(expr.args[2])

strip_module(expr::QuoteNode) = expr.value
strip_module(expr::Symbol) = expr
strip_module(expr::Expr) = strip_module(Val(expr.head), expr)

in_curly(::Val, expr) = (expr, [])
in_curly(::Val{:curly}, expr) = (expr.args[1], collect(expr.args[2:end]))
in_curly(expr::Symbol) = (expr, [])
in_curly(expr::Expr) = in_curly(Val(expr.head), expr)

create_scope_dicts(list) = Dict(map(x->to_scope(x)=>Dependencies(), list))

function check_supertype_sig(ifc_expr::Expr, ifc_scope, sig_scope, ifc, sig, sigarg, errors, source)
        # take a look at the supertype
        sigarg === Any && return

        @debug "check_supertype_sig" ifc_expr=ifc_expr sigarg=sigarg

        super_sig = supertype(sigarg)
        super_expr = strip_module(Meta.parse(string(super_sig)))

        @debug "check_supertype_sig" super_expr=super_expr
        # update the wparam context for the signature
        d = last(sig_scope)
        sig_expr, sig_sub_scope = find_scope(super_expr)
        replace_vector_with_array(sig_expr)
        merge!(d, create_scope_dicts(sig_sub_scope))

        bind_sig_symbols(ifc_expr, sig_expr, ifc_scope, sig_scope, ifc, sig, super_sig, errors, source)
end

function collect_unions(expr)
    base, typevars = in_curly(expr)
    base === :Union && return mapreduce(collect_unions, append!, typevars; init=[])
    [expr]
end

function collect_orset_to_check(ifc_expr, sig_expr, errors)
    ifc_unions = collect_unions(ifc_expr)
    sig_unions = collect_unions(sig_expr)

    # handle the `Union{}` case
    if xor(isempty(ifc_unions), isempty(sig_unions))
        return [(ifc_expr, sig_expr)]
    end

    collect(Iterators.product(ifc_unions, sig_unions))
end

function collect_orset_to_check(ifc_expr::Number, sig_expr, errors)
    if ifc_expr != sig_expr
        push!(errors, """Invalid numeric type paramter $sig_expr:
            \tdefined :  $(sig_expr)
            \texpected:  $(ifc_expr)""")
    end
    return []
end

function bind_sig_symbols(ifc_expr::Number, sig_expr, ifc_scope, sig_scope, ifc, sig, sigarg, errors, source)
    if ifc_expr != sig_expr
        push!(errors, """Invalid numeric type paramter $sig_expr:
            \tdefined :  $(sig_expr)
            \texpected:  $(ifc_expr)""")
    end
end

function bind_sig_symbols(ifc_expr::Expr, sig_expr, ifc_scope, sig_scope, ifc, sig, sigarg, errors, source)
    base_ifc_expr, ifc_typevars = in_curly(ifc_expr)
    base_sig_expr, sig_typevars = in_curly(sig_expr)

    @debug "bind_sig_symbols(Expr, Expr)" ifc=ifc_expr sig=sig_expr base_ifc=base_ifc_expr base_sig=base_sig_expr iscope=ifc_scope sscope=sig_scope

    if base_ifc_expr === base_sig_expr
        if length(ifc_typevars) != length(sig_typevars)
            sigparms = [sig.wparams...;[as_scope_expr(x) for x in keys(last(sig_scope))]...]
            ifcparms = [ifc.wparams...;[as_scope_expr(x) for x in keys(last(ifc_scope))]...]
            if !(eval_expr(sig.mod, sig_expr, sigparms) <: eval_expr(ifc.mod, ifc_expr, ifcparms))
                push!(errors, """Invalid typevar length for paramter $sig_expr:
                    \tdefined :  $(sig_expr)
                    \texpected:  $(ifc_expr)""")
                return
            end
        end
        for (ifc_typevar, sig_typevar) in zip(ifc_typevars, sig_typevars)
            parms = [sig.wparams...;[as_scope_expr(x) for x in keys(last(sig_scope))]...]
            evald_sigtypevar = eval_expr(sig.mod, sig_typevar, parms)
            walk_arg_depth(ifc_typevar, sig_typevar, ifc_scope, sig_scope, ifc, sig, evald_sigtypevar, source)
        end
    else
        check_supertype_sig(ifc_expr, ifc_scope, sig_scope, ifc, sig, sigarg, errors, source)
    end
end

function bind_sig_symbols(ifc_symbol::Symbol, sig_expr, ifc_scope, sig_scope, ifc, sig, sigarg, errors, source)
    @debug "bind_sig_symbols(Symbol)" ifc=ifc_symbol sig=sig_expr ifc_scope=ifc_scope
    # this handles the case where the interface is `:(I where I)`
    haskey(ifc_scope[end], ifc_symbol) && return

    for (i, iscope) in reverse(collect(enumerate(ifc_scope)))
        for (k,v) in iscope
            if k.var === ifc_symbol
                # push! binds symbols (the var itself)
                push!(v.bindings, sig_expr)
                # push_dependencies_into_scope! binds the variables in the dependencies (sub/supertyping)
                push_dependencies_into_scope!(k, sig_expr, ifc_scope[i:-1:1])
                @goto jump_out_of_loop;
            end
        end
    end
    @label jump_out_of_loop

    @debug "bind_sig_symbols(Symbol)" iscope=ifc_scope sscope=sig_scope
end

function push_dependencies_into_scope!(scope::Scope, binding, ifc_scope)
    symbols = rule_symbols(scope)
    @debug "push_dependencies_into_scope!" scope=scope symbols=symbols binding=binding ifc_scope=ifc_scope
    !isempty(symbols) || return

    for scp in ifc_scope
        for (k,v) in scp
            for s in symbols
                if s === k.var
                    push!(v.dependencies, scope=>binding)
                end
            end
        end
    end

end

function check_for_mismatched_scoping(mod, ifc, sig, ifc_scope, sig_scope, errors)
    scope_wheres = wheres(ifc_scope, sig)

    scope = pop!(ifc_scope)
    pop!(sig_scope)

    @debug "check_for_mismatched_scoping" scope=scope

    for (k,v) in scope
        types = unique(v.bindings)
        if length(types) > 1
            push!(errors, """Invalid typing--interface typevar bound to multiple types
                \ttypevar:  $(k.var)
                \ttypes  :  $types""")
        end
    end

    # TODO must also check for invalid subtyping based on scoping
    # e.g. var"#s10" <: T => Float64 and T => AbstractFloat
    vars = Symbol[k.var for k in keys(scope)]
    for (_,v) in scope
        if isempty(v.dependencies)
            continue
        end

        for dep in v.dependencies
            scp::Scope = first(dep)
            rule_symbs = rule_symbols(scp)
            rule_binding = last(dep)

            if !all(in(vars), rule_symbs)
                @warn """Cannot currently check depencency $dep in current scope;
                        only $vars are available.  In the future, we must collect
                        bindings and push them into the parent scope."""
                continue
            end

            d = Dict()
            for v in rule_symbs
                try
                    vals = get(scope, v, nothing).bindings
                    if !isempty(vals)
                        push!(d, v=>vals[1])
                    else
                        # TODO push error and continue
                    end
                catch e
                    @debug "$e" scope=scope v=v d=d
                    rethrow()
                end
            end

            new_rule = copy(scp.rule)
            replace_rule_symbol!(new_rule, scp.var, rule_binding)
            for (k, rep) in d
                replace_rule_symbol!(new_rule, k, rep)
            end

            error_rule = copy(new_rule)
            if new_rule.head === :(<:) || new_rule.head === :(>:)
                new_rule.args[1] = wrap_in_where(new_rule.args[1], scope_wheres)
                new_rule.args[2] = wrap_in_where(new_rule.args[2], scope_wheres)
            end

            @debug "check_for_mismatched_scoping" new_rule=new_rule
            # TODO try/catch around this; cannot judge
            try
                if !mod.eval(new_rule)
                    push!(errors, """Invalid typing--interface typevars misbound
                    \trule        :  $(scp.rule)
                    \tbound rule  :  $(error_rule)""")
                end
            catch e
                @debug """Cannot currently check unrealized dependency scopes
                \t$(new_rule)
                """
                # occurs when we are checking against type parameters
                # when we do not have a where for them;
                # FIXME add where clause to dependencies
            end
        end
    end
end

function walk_arg_depth(ifc_expr, sig_expr, ifc_scope, sig_scope, ifc, sig, sig_evald_arg, source)
    bare_ifc_expr, sub_ifc_scope = find_scope(ifc_expr)
    bare_sig_expr, sub_sig_scope = find_scope(sig_expr)

    last_ifc_scope = [as_scope_expr(x) for x in keys(last(ifc_scope))]
    last_sig_scope = [as_scope_expr(x) for x in keys(last(sig_scope))]
    @debug "walk_arg_depth" ifc=ifc_expr sig=sig_expr iscope=ifc_scope sscope=sig_scope

    push!(ifc_scope, create_scope_dicts(sub_ifc_scope))
    push!(sig_scope, create_scope_dicts(sub_sig_scope))

    allerrors = String[]
    orset = collect_orset_to_check(bare_ifc_expr, bare_sig_expr, allerrors)
    @debug "walk_arg_depth" orset=orset sis=sub_ifc_scope sss=sub_sig_scope
    for (i_expr, s_expr) in orset
        errors = String[]
        i_scope = copy(ifc_scope)
        s_scope = copy(sig_scope)

        evald, s_eval, i_expr, s_expr = eval_in_union(sig.mod, i_expr, s_expr,
                            [ifc.wparams...;last_ifc_scope...;sub_ifc_scope...],
                            [sig.wparams...;last_sig_scope...;sub_sig_scope...],
                            allerrors)
        evald || continue

        bind_sig_symbols(i_expr, s_expr, i_scope, s_scope, ifc, sig, s_eval, errors, source)
        check_for_mismatched_scoping(sig.mod, ifc, sig, i_scope, s_scope, errors)
        isempty(errors) || push!(allerrors, join(unique(errors), "\n"))
    end

    pop!(ifc_scope)
    pop!(sig_scope)

    @debug "walk_arg_depth" lenerrors=length(allerrors)
    (length(allerrors) < length(orset)) || isempty(allerrors) || error(string(source, "\n", join(allerrors, "\n")))
end

function check_signature_type_params_against_interface(ifc, sig, source)
    # if the interface has no type parameters, then the previous subtyping checks cover everything
    isempty(ifc.wparams) && return nothing
    ifc_scope = Any[create_scope_dicts(ifc.wparams)]
    sig_scope = Any[create_scope_dicts(sig.wparams)]

    @debug "check_signature_type_params_against_interface" allsig=sig.args allifc=ifc.args sigrv=sig.retval ifcrv=ifc.retval

    for (n, (sig_expr, ifc_expr)) in enumerate(zip([sig.args;sig.retval], [ifc.args;ifc.retval]))
        # check both expressions against each other, and strip only when sig_params does not match ifcparams?
        ifc_expr = stripwheres(ifc_expr, ifc.wparams)
        sig_expr = stripwheres(sig_expr, sig.wparams)

        @debug "check_signature_type_params_against_interface" ifcexpr=ifc_expr sigexpr=sig_expr

        walk_arg_depth(ifc_expr, sig_expr, ifc_scope, sig_scope, ifc, sig, ntharg(sig, n), source)
    end

    @debug "check_signature_type_params_against_interface" ifc_scope=ifc_scope
    errors = []
    check_for_mismatched_scoping(sig.mod, ifc, sig, ifc_scope, sig_scope, errors)
    isempty(errors) || error(string(source, "\n", join(errors, "\n")))
end

function __validate_signature(fkey, fname, sig::Signature, source)
    ifc = get_interface(fkey)
    ifc !== nothing || error("Interface not defined for $fname (impl at $source)")

    length(sig.args_evald) == length(ifc.args_evald) || error("""Invalid argument length for $fname at $source:
    \tdefined :  $(sig.args_evald)
    \texpected:  $(ifc.args_evald)""")

    subtype = map(<:, sig.args_evald, ifc.args_evald)
    all(subtype) || error("""Invalid definition of $fname at $source:
    \tdefined :  $(sig.args_evald)
    \texpected:  $(ifc.args_evald)
    \tdefined :  $(sig.args)
    \texpected:  $(ifc.args)
    \tat      :  $subtype""")

    (sig.retval_evald <: ifc.retval_evald) || error("""Invalid return type for $fname at $source:
    \tdefined :  $(sig.retval_evald)
    \texpected:  $(ifc.retval_evald)
    \tdefined :  $(sig.retval)
    \texpected:  $(ifc.retval)""")

    @debug "__validate_signature" f=fname ifc=ifc sig=sig

    check_signature_type_params_against_interface(ifc, sig, source)
end

function __validate_signature(mod, source, fkey::Symbol, f, args::Vector{<:Vector}, wparams)
    body(t::UnionAll) = body(t.body)
    body(t::Type) = t

    for arglist in args
        @debug "__validate_signature" f=string(f) args=arglist wparams=wparams
        sig = Signature(mod, arglist[2:end], arglist[1], wparams)
        try
            __validate_signature(fkey, Symbol(f), sig, source)
        catch e
            if !(e isa ErrorException)
                @error("Error in function $f at $source")
            end
            rethrow()
        end
    end
end


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

export SimData

"""
    SimData

Data structure containing the time series data for a simulation.

It is a collection of [`TSeries`](@ref) of the same frequency and containing
data for the same range. When used for simulation, the range must include the
initial conditions, the simulation range and the final conditions, although it
could extend beyond that. It must contain time series for all variables and
shocks in the model, although it might contain other time series.

"""
struct SimData{F <: Frequency,C <: AbstractMatrix{Float64}} <: AbstractMatrix{Float64}
    firstdate::MIT{F}
    columns::NamedTuple
    values::C

    # inner constructor enforces constraints.
    function SimData(fd::MIT, names::NTuple{N,Symbol}, values::AbstractMatrix) where {N}
        if N != size(values, 2)
            throw(ArgumentError("Number of names and columns don't match: $N ≠ $(size(values, 2))."))
        end
        columns = NamedTuple{names}([TSeries(fd, view(values, :, i)) for i in 1:N ])
        new{frequencyof(fd),typeof(values)}(fd, columns, values)
    end
end

# External constructors
"""
    SimData(fd, vars, data)

Construct an instance of SimData with the given variable names and data. The
first date is provided in the `fd::MIT` argument, which also contains
information about the frequency. The number of names must match the number of
columns in `data`. The names must be String or Symbol or anything that converts
to Symbol.

"""
SimData(fd, vars=(), data=zeros(0,0)) = SimData(fd, tuple(Symbol.(vars)...), data)
SimData(rng::AbstractUnitRange{<:Integer}) = SimData(first(rng), (), zeros(length(rng),0))

TimeSeriesEcon.firstdate(sd::SimData) = getfield(sd, :firstdate)
TimeSeriesEcon.lastdate(sd::SimData) = (firstdate(sd) - 1) + size(_values(sd), 1)
TimeSeriesEcon.frequencyof(::SimData{F}) where F <: Frequency = F
TimeSeriesEcon.mitrange(sd::SimData) = (firstdate(sd) - 1) .+ Base.axes1(sd)

@inline _columns(sd::SimData) = getfield(sd, :columns)
@inline _names(sd::SimData) = keys(getfield(sd, :columns))
@inline _values(sd::SimData) = getfield(sd, :values)
@inline _values_view(sd::SimData, i1, i2) = view(getfield(sd, :values), i1, i2)
@inline _values_slice(sd::SimData, i1, i2) = getindex(getfield(sd, :values), i1, i2)

Base.size(sd::SimData) = size(_values(sd))
Base.IndexStyle(sd::SimData) = IndexStyle(_values(sd))

# Indexing with integers falls back to AbstractArray
const FallbackType = Union{Integer,Colon,AbstractUnitRange{<:Integer},AbstractVector{<:Integer},CartesianIndex}
Base.getindex(sd::SimData, i1::FallbackType...) = getindex(_values(sd), i1...)
Base.setindex!(sd::SimData, val, i1::FallbackType...) = setindex!(_values(sd), val, i1...)

Base.similar(sd::SimData) = SimData(firstdate(sd), _names(sd), similar(_values(sd)))

Base.:(==)(a::SimData, b::SimData) = frequencyof(a) == frequencyof(b) && firstdate(a) == firstdate(b) && _values(a) == _values(b)

Base.dataids(sd::SimData) = Base.dataids(_values(sd))

export rawdata, colnames
rawdata(sd::SimData) = _values(sd)
colnames(sd::SimData) = keys(_columns(sd))
Base.pairs(sd::SimData) = pairs(_columns(sd))

const ColumnTypes = Union{Symbol,ModelVariable,AbstractString}
# Define dot access to columns
Base.propertynames(sd::SimData) = _names(sd)
Base.getproperty(sd::SimData, col::ColumnTypes) = getproperty(sd, Symbol(col))
function Base.getproperty(sd::SimData, col::Symbol)
    if col in _names(sd)
        return getfield(_columns(sd), col)
    else
        throw(BoundsError(sd, [col,]))
    end
end

# Access to columns by [:xyz] notation
Base.getindex(sd::SimData, col::ColumnTypes) = getproperty(sd, col)
Base.setindex!(sd::SimData, val, col::ColumnTypes) = setproperty!(sd, col, val)

Base.setproperty!(sd::SimData, col::ColumnTypes, val) = setproperty!(sd, Symbol(col), val)
function Base.setproperty!(sd::SimData, col::Symbol, val)
    try
        col = getproperty(sd, col)
    catch BoundsError
        error("Cannot assign new data column this way. Use hcat(sd, column=value, ...) instead")
    end
    if Base.mightalias(col, val)
        val = copy(val)
    end
    if val isa TSeries && frequencyof(val) == frequencyof(sd)
        rng = intersect(mitrange(sd), mitrange(val))
        col[rng] = val[rng].values
    else
        return setindex!(col.values, val, :)
    end
end

@inline Base.getindex(sd::SimData, cols::AbstractVector{<:Number}) = getindex(sd, (_names(sd)[cols])...)
# @inline Base.getindex(sd::SimData, cols::AbstractVector) = getindex(sd, tuple(map(Symbol, cols)...))
@inline Base.getindex(sd::SimData, cols::AbstractVector{<:ModelVariable}) = getindex(sd, tuple(map(Symbol, cols)...))
@inline Base.getindex(sd::SimData, cols::Tuple) = getindex(sd, tuple(map(Symbol, cols)...))

function Base.getindex(sd::SimData, cols::NTuple{N,Symbol}) where N
    SimData(firstdate(sd), tuple(cols...), hcat((getproperty(sd, c) for c in cols)...))
end

function Base.hcat(sd::SimData; KW...)
    l1 = size(sd, 1)
    as_vect(v::Number) = fill(Float64(v), l1)
    as_vect(v) = v
    names = (_names(sd)..., keys(KW)...)
    vals = hcat(_values(sd), (as_vect(v) for v in values(KW))...)
    return SimData(firstdate(sd), names, vals)
end

# Indexing access with MIT

function check_frequency(a, b)
    Fa = frequencyof(a)
    Fb = frequencyof(b)
    Fa == Fb || throw(ArgumentError("wrong frequency: expected $Fa got $Fb."))
end

macro if_same_frequency(a, b, expr)
    :( frequencyof($a) == frequencyof($b) ? $(expr) : throw(ArgumentError("Wrong frequency")) ) |> esc
end

# A single MIT materializes as a NamedTuple (row of the matrix with column names attached to the values)
Base.getindex(sd::SimData, i1::MIT, ::Colon) = sd[i1 - firstdate(sd) + 1, :]
function Base.getindex(sd::SimData, i1::MIT)
    check_frequency(sd, i1)
    if firstdate(sd) <= i1 <= lastdate(sd)
        return @inbounds NamedTuple{_names(sd)}(_values_view(sd, i1 - firstdate(sd) + 1, :))
    else
        throw(BoundsError(sd, i1))
    end
end

# Modifying a row in a table -> one must pass in a vector
Base.setindex!(sd::SimData, val, i1::MIT, ::Colon) = begin sd[i1 - firstdate(sd) + 1, :] = val end
function Base.setindex!(sd::SimData, val::AbstractVector{<:Real}, i1::MIT)
    check_frequency(sd, i1)
    if firstdate(sd) <= i1 <= lastdate(sd)
        row = i1 - firstdate(sd) + 1
        setindex!(_values_view(sd, row, :), val, :)
        return _values(sd)[row,:]
    else
        throw(BoundsError(sd, i1))
    end
end

function Base.setindex!(sd::SimData, val::NamedTuple, i1::MIT)
    check_frequency(sd, i1)
    if firstdate(sd) <= i1 <= lastdate(sd)
        for (n, v) in pairs(val)
            setindex!(getproperty(sd, n), v, i1)
        end
        return sd[i1]
    else
        throw(BoundsError(sd, i1))
    end
    return sd[i1]
end

# A selection of several rows returns a slice from the original SimData
function Base.getindex(sd::SimData, i1::AbstractUnitRange{<:MIT})
    check_frequency(sd, i1)
    if firstdate(sd) <= minimum(i1) <= maximum(i1) <= lastdate(sd)
        return SimData(first(i1), _names(sd), _values_slice(sd, i1 .- firstdate(sd) .+ 1, :))
    else
        throw(BoundsError(sd, i1))
    end
end

function Base.setindex!(sd::SimData, val, i1::AbstractUnitRange{<:MIT})
    check_frequency(sd, i1)
    if firstdate(sd) <= minimum(i1) <= maximum(i1) <= lastdate(sd)
        rows = i1 .- firstdate(sd) .+ 1
    else
        throw(BoundsError(sd, i1))
    end
    if val isa Number
        fill!(_values_view(sd, rows, :), val)
    else
        setindex!(_values_view(sd, rows, :), val, :, :)
    end
    return sd[rows,:]
end

# Base.getindex(sd::SimData, ::Colon, c) = getindex(sd, mitrange(sd), c)

Base.getindex(sd::SimData, r, c) = getproperty(sd, Symbol(c))[r]
Base.getindex(sd::SimData, r, c::AbstractVector) = getindex(sd, r, tuple(map(Symbol, c)...))
# Base.getindex(sd::SimData, r, c::Tuple) = getindex(sd, r, tuple(map(Symbol, c)...))

Base.getindex(sd::SimData, r::MIT, c::NTuple{N,Symbol}) where N = (a = sd[r]; NamedTuple{c}([a[cc] for cc in c]))
function Base.getindex(sd::SimData, r::AbstractUnitRange{<:MIT}, c::NTuple{N,Symbol}) where N
    check_frequency(sd, r)
    if firstdate(sd) <= first(r) <= last(r) <= lastdate(sd)
        col_inds = indexin(c, [_names(sd)...])
        for (cc, cn) in zip(c, col_inds)
            if cn === nothing
                throw(BoundsError(sd, [cc]))
            end
        end
        return SimData(first(r), c, _values_slice(sd, r .- firstdate(sd) .+ 1, col_inds))
    else
        throw(BoundsError(sd, r))
    end
end

Base.setindex!(sd::SimData, v, r, c::ColumnTypes) = setindex!(getproperty(sd, c), v, r)
Base.setindex!(sd::SimData, v, r, c::AbstractVector) = setindex!(sd, v, r, tuple(map(Symbol, c)...))
Base.setindex!(sd::SimData, v, r, c::Tuple) = setindex!(sd, v, r, tuple(map(Symbol, c)...))

function Base.setindex!(sd::SimData, v, r::MIT, c::NTuple{N,Symbol}) where N
    for (cc, vv) in zip(c, v)
        setindex!(sd, vv, r, cc)
    end
    return sd[r, c]
end

function Base.setindex!(sd::SimData, v, r::AbstractUnitRange{<:MIT}, c::NTuple{N,Symbol}) where N
    check_frequency(sd, r)
    if firstdate(sd) <= first(r) <= last(r) <= lastdate(sd)
        cols = indexin(c, [_names(sd)...])
        for (cc, cn) in zip(c, cols)
            if cn === nothing
                throw(BoundsError(sd, [cc]))
            end
        end
        rows = r .- firstdate(sd) .+ 1
        if v isa Number
            v = fill(Float64(v), length(r) * length(c))
        end
        setindex!(_values_view(sd, rows, cols), v, :, :)
    else
        throw(BoundsError(sd, r))
    end
end

#### Pretty printing

sprint_names(names) = length(names) > 10 ? "$(length(names)) variables" : "variables (" * join(names, ",") * ")"
Base.summary(io::IO, sd::SimData) = isempty(sd) ?
        print(IOContext(io, :limit => true), "Empty SimData with ", sprint_names(_names(sd)), " in ", mitrange(sd)) :
        print(IOContext(io, :limit => true), size(sd, 1), '×', size(sd, 2), " SimData with ", sprint_names(_names(sd)), " in ", mitrange(sd))

Base.show(io::IO, ::MIME"text/plain", sd::SimData) = show(io, sd)
function Base.show(io::IO, sd::SimData)
    summary(io, sd)
    isempty(sd) && return
    print(io, ":")
    limit = get(io, :limit, true)
    nval, nsym = size(sd)

    from = firstdate(sd)
    dheight, dwidth = displaysize(io)
    if get(io, :compact, nothing) === nothing
        io = IOContext(io, :compact => true)
    end
    dwidth -= 11

    names_str = let
        names = map(_names(sd)) do n
            sn = string(n)
            return sn[1:min(10, length(sn))]
        end
        reshape([names...], 1, :)
    end
    sd_with_names = [names_str; _values(sd)]
    A = Base.alignment(io, sd_with_names, axes(sd_with_names, 1), 1:nsym, dwidth, dwidth, 2)

    all_cols = true
    if length(A) ≠ nsym
        dwidth = div(dwidth - 1, 2)
        AL = Base.alignment(io, sd_with_names, axes(sd_with_names, 1), 1:nsym, dwidth, dwidth, 2)
        AR = reverse(Base.alignment(io, sd_with_names, axes(sd_with_names, 1), reverse(1:nsym), dwidth, dwidth, 2))
        Linds = [1:length(AL)...]
        Rinds = [nsym - length(AR) + 1:nsym...]
        all_cols = false
    end

    local vdots = "\u22ee"
    local hdots = " \u2026 "
    local ddots = " \u22f1 "

    print_aligned_val(io, v, (al, ar), showsep=true; sep=showsep ? "  " : "") = begin
        sv = sprint(print, v, context=io, sizehint=0)
        if v isa Number
            vl, vr = Base.alignment(io, v)
        else
            if length(sv) > al + ar
                sv = sv[1:al + ar - 1] * '…'
            end
            vl, vr = al, length(sv) - al
        end
        print(io, repeat(" ", al - vl), sv, repeat(" ", ar - vr), sep)
    end

    print_colnames(io, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
        local nLcols = length(Lcols)
        local nRcols = length(Rcols)
        for (i, (col, align)) in enumerate(zip(Lcols, LAligns))
            print_aligned_val(io, names[col], align, i < nLcols)
        end
        nRcols == 0 && return
        print(io, hdots)
        for (i, (col, align)) in enumerate(zip(Rcols, RAligns))
            print_aligned_val(io, names[col], align, i < nRcols)
        end
    end

    print_rows(io, rows, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
        local nLcols = length(Lcols)
        local nRcols = length(Rcols)
        for row in rows
            mit = from + (row - 1)
            print(io, '\n', lpad(mit, 8), " : ")
            for (i, (val, align)) in enumerate(zip(vals[row, Lcols], LAligns))
                print_aligned_val(io, val, align, i < nLcols)
            end
            nRcols == 0 && continue
            print(io, hdots)
            for (i, (val, align)) in enumerate(zip(vals[row, Rcols], RAligns))
                print_aligned_val(io, val, align, i < nRcols)
            end
        end
    end

    print_vdots(io, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
        print(io, '\n', repeat(" ", 11))
        local nLcols = length(Lcols)
        local nRcols = length(Rcols)
        for (i, (col, align)) in enumerate(zip(Lcols, LAligns))
            print_aligned_val(io, vdots, align, i < nLcols)
        end
        nRcols == 0 && return
        print(io, ddots)
        for (i, (col, align)) in enumerate(zip(Rcols, RAligns))
            print_aligned_val(io, vdots, align, i < nRcols)
        end
    end

    names = _names(sd)
    vals = _values(sd)
    if !limit
        print(io, "\n", repeat(" ", 11))
        print_colnames(io, 1:nsym, A)
        print_rows(io, 1:nval, 1:nsym, A)
    elseif nval > dheight - 6 # all rows don't fit
                # unable to show all rows
        if all_cols
            print(io, "\n", repeat(" ", 11))
            print_colnames(io, 1:nsym, A)
            top = div(dheight - 6, 2)
            print_rows(io, 1:top, 1:nsym, A)
            print_vdots(io, 1:nsym, A)
            bot = nval - dheight + 7 + top
            print_rows(io, bot:nval, 1:nsym, A)
        else # not all_cols
            print(io, "\n", repeat(" ", 11))
            print_colnames(io, Linds, AL, Rinds, AR)
            top = div(dheight - 6, 2)
            print_rows(io, 1:top, Linds, AL, Rinds, AR)
            print_vdots(io, Linds, AL, Rinds, AR)
            bot = nval - dheight + 7 + top
            print_rows(io, bot:nval, Linds, AL, Rinds, AR)
        end # all_cols
    else # all rows fit
        if all_cols
            print(io, '\n', repeat(" ", 11))
            print_colnames(io, 1:nsym, A)
            print_rows(io, 1:nval, 1:nsym, A)
        else
            print(io, '\n', repeat(" ", 11))
            print_colnames(io, Linds, AL, Rinds, AR)
            print_rows(io, 1:nval, Linds, AL, Rinds, AR)
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractDomain

An abstract type to aid with custom spatial domain types. Users
can subtype their domains from this type, and implement the methods
in `geotraits/domain.jl`.
"""
abstract type AbstractDomain end

"""
    view(domain, inds)

Return a view of `domain` at given indices `inds`.
"""
Base.view(domain::AbstractDomain, inds::AbstractVector{Int}) =
  DomainView(domain, inds)

#------------------
# IMPLEMENTATIONS
#------------------
"""
    SpatialDomain{T,N}

Spatial domain in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type SpatialDomain{T,N} <: AbstractDomain end

ncoords(::SpatialDomain{T,N}) where {T,N} = N
coordtype(::SpatialDomain{T,N}) where {T,N} = T

include("domains/point_set.jl")
include("domains/regular_grid.jl")
include("domains/structured_grid.jl")


text:
using WebPlayer, FileIO, Images, Blink
using Base.Test, Colors
img = Float64.(Gray.(restrict(restrict(load(homedir() * "/Desktop/pliss.jpg")))))
video = Array{Float64, 3}(size(img)..., 100)
for f = 1:100
    r = rand() * 0.1
    for i = 1:size(img, 1), j = 1:size(img, 2)
        video[i, j, f] = clamp(img[i, j] + r, 0, 1)
    end
end

videos = [video, video];

x = WebPlayer.playvideo(videos, ["test", "test"], width = 200)
#
# w = Window()
# body!(w, x)
# tools(w)


text:
module JSONSchema

using JSON
import HTTP
import URIs

export Schema, validate

include("schema.jl")
include("validation.jl")

export diagnose
function diagnose(x, schema)
    Base.depwarn(
        "`diagnose(x, schema)` is deprecated. Use `validate(schema, x)` instead.",
        :diagnose,
    )
    ret = validate(schema, x)
    return ret === nothing ? nothing : sprint(show, ret)
end

end

max_stars_repo_name: fredo-dedup/JSONSchema

max_issues_repo_name: fredo-dedup/JSONSchema

max_forks_repo_name: fredo-dedup/JSONSchema


text:
# m5.4.AMs.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

# ### snippet 5.1

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

df = DataFrame(
    :A => df[:, :MedianAgeMarriage],
    :M => df[:, :Marriage],
    :D => df[:, :Divorce]
 )

scale!(df, [:M, :A, :D])

# Define the Stan language model

stan5_4_AM = "
data {
    int N;
    vector[N] A;
    vector[N] M;
}
parameters {
    real a;
    real bAM;
    real<lower=0> sigma;
}
model {
    vector[N] mu = a + bAM * M;
    a ~ normal( 0 , 0.2 );
    bAM ~ normal( 0 , 0.5 );
    sigma ~ exponential( 1 );
    A ~ normal( mu , sigma );
}
";

# Define the SampleModel
#tmpdir=ProjDir*"/tmp"
m5_4_AMs = SampleModel("m5.4.AM", stan5_4_AM);

# Input data

m5_4_data = Dict(
    "N" => size(df, 1),
    "M" => df[:, :M_s],
    "A" => df[:, :A_s]
);

# Sample using cmdstan

rc5_4_AMs = stan_sample(m5_4_AMs, data=m5_4_data);

if success(rc5_4_AMs)

    part5_4_AMs = read_samples(m5_4_AMs, :particles)
    part5_4_AMs |> display

end

# End of m5.3.AMs.jl


text:
using Plots

# define the Lorenz attractor
mutable struct Lorenz
    dt; σ; ρ; β; x; y; z
end

function step!(l::Lorenz)
    dx = l.σ*(l.y - l.x)       ; l.x += l.dt * dx
    dy = l.x*(l.ρ - l.z) - l.y ; l.y += l.dt * dy
    dz = l.x*l.y - l.β*l.z     ; l.z += l.dt * dz
end

attractor = Lorenz((dt = 0.02, σ = 10., ρ = 28., β = 8//3, x = 1., y = 1., z = 1.)...)

# initialize a 3D plot with 1 empty series
plt = plot3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
                title = "Lorenz Attractor", marker = 2)

# build an animated gif by pushing new points to the plot, saving every 10th frame
@gif for i=1:1500
    step!(attractor)
    push!(plt, attractor.x, attractor.y, attractor.z)
end every 10

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function test_random_mask_layer(backend, T, eps)
  println("-- Testing RandomMask on $(typeof(backend)){$T}")

  n_inputs = 3
  tensor_dims = [abs(rand(Int)) % 6 + 1 for i = 1:n_inputs]
  dims_all = [tuple((abs(rand(Int, tensor_dims[i])) % 8 + 1) ...) for i = 1:n_inputs]
  println("    > $n_inputs input blobs with tensor dims $tensor_dims")

  ratio = min(abs(rand()) + 0.1, 0.9)
  inputs = Array[rand(T, dims_all[i]) for i = 1:n_inputs]
  input_blobs = Blob[make_blob(backend, inputs[i]) for i = 1:n_inputs]
  diff_blobs = Blob[make_blob(backend, inputs[i]) for i = 1:n_inputs]

  println("    > Setup")
  layer = RandomMaskLayer(bottoms=[symbol("inputs-$i") for i = 1:n_inputs], ratio=ratio)
  state = setup(backend, layer, input_blobs, diff_blobs)

  println("    > Forward")
  forward(backend, state, input_blobs)
  for i = 1:n_inputs
    rand_vals = to_array(state.dropouts[i].rand_vals)
    expected_output = inputs[i] .* (rand_vals .> ratio)
    got_output = to_array(input_blobs[i])
    @test all(abs(got_output - expected_output) .< eps)
  end

  println("    > Backward")
  top_diffs = [rand(T, dims_all[i]) for i = 1:n_inputs]
  for i = 1:n_inputs
    copy!(diff_blobs[i], top_diffs[i])
  end
  backward(backend, state, input_blobs, diff_blobs)
  for i = 1:n_inputs
    rand_vals = to_array(state.dropouts[i].rand_vals)
    got_grad = to_array(diff_blobs[i])
    expected_grad = top_diffs[i] .* (rand_vals .> ratio)
    @test all(abs(got_grad - expected_grad) .< eps)
  end

  shutdown(backend, state)
end
function test_random_mask_layer(backend::Backend)
  test_random_mask_layer(backend, Float64, 1e-10)
  test_random_mask_layer(backend, Float32, 1e-4)
end

if test_cpu
  test_random_mask_layer(backend_cpu)
end
if test_gpu
  test_random_mask_layer(backend_gpu)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "nv_codec_headers"
version = v"10.0.26"

# Collection of sources required to build this package
sources = [
    ArchiveSource("https://github.com/FFmpeg/nv-codec-headers/releases/download/n10.0.26.0/nv-codec-headers-10.0.26.0.tar.gz",
                  "2552723fef0adbed1c327e9cca167226b35d8a5d17fb86c753b71d70b26f141c"),
]

# Bash recipe for building across all platforms
script = raw"""
cd nv-codec-headers-*
make install PREFIX=${prefix}
install_license /usr/share/licenses/MIT
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("include/ffnvcodec/dynlink_cuda.h", :dynlink_cuda_h),
    FileProduct("include/ffnvcodec/dynlink_cuviddec.h", :dynlink_cuviddec_h),
    FileProduct("include/ffnvcodec/dynlink_loader.h", :dynlink_loader_h),
    FileProduct("include/ffnvcodec/dynlink_nvcuvid.h", :dynlink_nvcuvid_h),
    FileProduct("include/ffnvcodec/nvEncodeAPI.h", :nvEncodeAPI_h),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil