text:
abstract type ResamplingStrategy{T} end

function Base.getindex(sampling::ResamplingStrategy, idx::Int)
    if idx > length(sampling)
        throw(ArgumentError("index must satisfy 1 <= idx <= length(sampling)"))
    end

    return (getTrainData(sampling, idx), getTestData(sampling, idx))
end

Base.iterate(sample::ResamplingStrategy, state=1) = state > length(sample) ? nothing : (sample[state], state + 1)

get(sampling::ResamplingStrategy) = [sample for sample in sampling]


text:
using SafeTestsets

@safetestset "Types                                                                           " begin include("types_tests.jl") end
@safetestset "Utils                                                                           " begin include("utils_tests.jl") end
@safetestset "Data Series                                                                     " begin include("dataseries_tests.jl") end
@safetestset "Time Series                                                                     " begin include("timeseries_tests.jl") end
@safetestset "Solution                                                                        " begin include("solution_tests.jl") end

ext: jl


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, 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(...)

A variant of `modify!` which modifies a copy of the passed data set.

See [`modify!`](@ref)
"""
modify(origninal_ds::AbstractDataset, @nospecialize(args...); threads::Bool = true) = modify!(copy(origninal_ds), args...; threads = threads)

modify!(ds::Dataset; threads::Bool = true) = parent(ds)

"""
    modify!(ds::AbstractDataset, args...; [threads = true])

Modify columns of a data set. The `args` arguments must be in the form of `cols => fun => newnames`. The `fun` function will be called on passed `cols`, with the excpetion of two special functions: `byrow` and `splitter`. `fun` assumes a single column as its input, thus passing multiple columns will be broadcasted, i.e. `cols => fun` will be translated to `col1=>fun`, `col2=>fun`, .... When `newname` is not provided `modify!` modifies the passed column.

When a grouped data set is passed to `modify!`, the operation is done on each group of observations.

each `args` can be constructed based on columns in the original data set or the columns which have been created before it.

# Special functions

`byrow` and `splitter` are two special functions which can be passed as `fun`.

`byrow` can accept multiple columns as input and does a given operation on each row of the data set. When a single column is passed to `byrow`, `modify!` modifies the passed column, however, when multiple columns are passed, `byrow` applies the row-wise operation on them and creates a new column.

`splitter` splits a column of tuples to multiple columns. When `splitter` is set as `fun` the `newnames` must be given.

# Using multivariate functions

To pass multiple columns to a `fun` function which operates on multiple inputs, the columns must be passed as tuple of column names, or column indices.

See [`modify`](@ref)

# Examples

```jldoctest
julia> ds = Dataset(x1 = 1:5, x2 = [-2, -1, missing, 1, 2],
                    x3 = [0.0, 0.1, 0.2, missing, 0.4])
5×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Float64?
─────┼──────────────────────────────
   1 │        1        -2       0.0
   2 │        2        -1       0.1
   3 │        3   missing       0.2
   4 │        4         1   missing
   5 │        5         2       0.4

julia> modify!(ds, 2:3 => IMD.sum)
5×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Float64?
─────┼──────────────────────────────
   1 │        1         0       0.7
   2 │        2         0       0.7
   3 │        3         0       0.7
   4 │        4         0       0.7
   5 │        5         0       0.7

julia> modify!(ds, :x1 => x -> x .- mean(x))
5×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Float64?  Int64?    Float64?
─────┼──────────────────────────────
   1 │     -2.0         0       0.7
   2 │     -1.0         0       0.7
   3 │      0.0         0       0.7
   4 │      1.0         0       0.7
   5 │      2.0         0       0.7

julia> body = Dataset(weight = [78.5, 59, 80], height = [160, 171, 183])
3×2 Dataset
 Row │ weight    height
     │ identity  identity
     │ Float64?  Int64?
─────┼────────────────────
   1 │     78.5       160
   2 │     59.0       171
   3 │     80.0       183

julia> modify!(body, :height => byrow(x -> (x/100)^2) => :BMI, [1, 3] => byrow(/) => :BMI)
3×3 Dataset
 Row │ weight    height    BMI
     │ identity  identity  identity
     │ Float64?  Int64?    Float64?
─────┼──────────────────────────────
   1 │     78.5       160   30.6641
   2 │     59.0       171   20.1771
   3 │     80.0       183   23.8884

julia> body = Dataset(weight = [78.5, 59, 80], height = [160, 171, 183])
3×2 Dataset
 Row │ weight    height
     │ identity  identity
     │ Float64?  Int64?
─────┼────────────────────
   1 │     78.5       160
   2 │     59.0       171
   3 │     80.0       183

julia> modify!(body, (:weight, :height)=> cor)
3×3 Dataset
 Row │ weight    height    cor_weight_height
     │ identity  identity  identity
     │ Float64?  Int64?    Float64?
─────┼───────────────────────────────────────
   1 │     78.5       160          0.0890411
   2 │     59.0       171          0.0890411
   3 │     80.0       183          0.0890411
```
"""
function modify!(ds::AbstractDataset, @nospecialize(args...); threads::Bool = true)
    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, threads)
    else
        _modify(ds, normalize_modify_multiple!(idx_cpy, index(ds), args...))
    end
end

# we must take care of all possible types, because, fallback is slow
# _is_scalar(::T, sz) where T <: Number = true
# _is_scalar(::Missing, sz) = true
# _is_scalar(::T, sz) where T <: Tuple = true
# _is_scalar(::TimeType, sz) = true
# _is_scalar(::T, sz) where T <: AbstractString = true
_is_scalar(x, sz) = true
_is_scalar(x::T, sz) where T <: AbstractVector = length(x) != sz

# # TODO can we memorise this and avoid calling it repeatedly in a sesssion
# _is_scalar_barrier(::Val{T}) where T = hasmethod(size, (T,))
#
# function _is_scalar(_res::T, sz) where T
#      resize_col = false
#     if _is_scalar_barrier(Val(T))
#         if size(_res) == () || size(_res,1) != sz
#             resize_col = true
#         end
#     else
#         resize_col = true
#     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!(allocatecol(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!(allocatecol(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!(allocatecol(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 might be a problem in the `byrow` usage, 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)
    @_threadsfor 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)
    @_threadsfor 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, threads)
    if (mssecond isa Base.Callable) && !(mslast isa MultiCol)
        T = _check_the_output_type(ds, msfirst=>mssecond=>mslast)
        _res = allocatecol(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), threads)
        else
            _modify_grouped_fill_one_col!(_res, _columns(ds)[msfirst], mssecond, index(ds).starts, index(ds).ngroups[], nrow(ds), threads)
        end
        ds[!, mslast] = _res
    elseif (mssecond isa Expr)  && mssecond.head == :BYROW
        #TODO we should think about how to pass threads here
        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, threads)
    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, threads)
    end
    return parent(ds)
end


text:
using VegaGraphs
using Test

@testset "VegaGraphs.jl" begin
     x = 2
    y = 2
    @test VegaGraphs.sum_values(x,y) == 4
end

max_stars_repo_name: jpbrasile/VegaGraphs

max_issues_repo_name: jpbrasile/VegaGraphs

max_forks_repo_name: jpbrasile/VegaGraphs


text:
# ------------------------------------------------------------------------------------------
# # Julia 中的线性代数
# > Based on work by Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)
# > with edits from Jane Herriman
#
# ## 目录
#  - [基础线性代数操作](#基础线性代数操作)
#     - [乘法](#乘法)
#     - [转置](#转置)
#     - [转置的乘法](#转置的乘法)
#     - [解线性方程组](#解线性方程组)
#  - [特殊的矩阵结构](#特殊的矩阵结构)
#      - [大规模问题](#大规模问题)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 基础线性代数操作
# 定义一个随机矩阵
# ------------------------------------------------------------------------------------------

A = rand(1:4,3,3)

# ------------------------------------------------------------------------------------------
# 定义一个元素全为 1 的向量
# ------------------------------------------------------------------------------------------

x = fill(1.0, (3))

# ------------------------------------------------------------------------------------------
# 注意 $A$ 的类型为 `Array{Int64,2}`，而 $x$ 的类型为 `Array{Float64,1}`。
#
# Julia 定义 `Array{Type,1}` 的别名为向量 `Vector{Type}`，`Array{Type,2}` 的别名为矩阵 `Matrix{Type}` 。
#
# 许多线性代数的基础操作和其他语言一样
#
# ### 乘法
# ------------------------------------------------------------------------------------------

b = A*x

# ------------------------------------------------------------------------------------------
# ### 转置
# 就像在其他语言中 `A'` 表示对 `A` 进行共轭转置
# ------------------------------------------------------------------------------------------

A'

# ------------------------------------------------------------------------------------------
# 我们还可以通过 `transpose` 函数获得转置矩阵
# ------------------------------------------------------------------------------------------

transpose(A)

# ------------------------------------------------------------------------------------------
# ### 转置的乘法
# Julia 中某些情况下可以省略 `*` 号
# ------------------------------------------------------------------------------------------

A'A

# ------------------------------------------------------------------------------------------
# ### 解线性方程组
# 用方阵 $A$ 表示的线性方程组 $Ax=b$ 可以用左除运算符（函数）`\` 求解
# ------------------------------------------------------------------------------------------

A\b

# ------------------------------------------------------------------------------------------
# ## 特殊的矩阵结构
#
# 矩阵结构在线性代数中非常重要。
# 接触一下大一些的线型系统就可以看到矩阵结构有*多*重要了。
#
# 用线性代数标准包 `LinearAlgebra` 可以获得结构化的矩阵（structured matrices）：
# ------------------------------------------------------------------------------------------

using LinearAlgebra

n = 1000
A = randn(n,n);

# ------------------------------------------------------------------------------------------
# Julia 可以推断特殊矩阵结构，比如判断对称矩阵
# ------------------------------------------------------------------------------------------

Asym = A + A'
issymmetric(Asym)

# ------------------------------------------------------------------------------------------
# 但有时候浮点错误会比较麻烦
# ------------------------------------------------------------------------------------------

Asym_noisy = copy(Asym)
Asym_noisy[1,2] += 5eps()

issymmetric(Asym_noisy)

# ------------------------------------------------------------------------------------------
# 幸运的是我们可以通过如 `Diagonal`，`Triangular`，`Symmetric`，`Hermitian`，`Tridiagonal`
# 和 `SymTridiagonal` 这样的函数来明确地创建矩阵
# ------------------------------------------------------------------------------------------

Asym_explicit = Symmetric(Asym_noisy);

# ------------------------------------------------------------------------------------------
# 我们来看看 Julia 计算 `Asym`，`Asym_noisy` 和 `Asym_explicit` 的特征值各要花多少时间
# ------------------------------------------------------------------------------------------

@time eigvals(Asym);

@time eigvals(Asym_noisy);

@time eigvals(Asym_explicit);

# ------------------------------------------------------------------------------------------
#
# 本例中，使用 `Symmetric()` 处理 `Asym_noisy` 后让计算效率提高了约5倍
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 大规模问题
# 使用 `Tridiagonal` 和 `SymTridiagonal` 类型储存三对角矩阵（tridiagonal matrices）
# 让处理大规模的三对角矩阵问题变为可能。
#
# 以下问题如果使用稠密的 `Matrix` 类型储存，在个人计算机上是无法进行求解的。
# ------------------------------------------------------------------------------------------

n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
module Interpret

using MLStyle

function interpret(expr)
    ctx = Dict{Symbol, Any}()
    interpret!(ctx, expr)
end

function interpret!(ctx , s::Symbol)
    # We'd use `get` for this, but `eval(s)` needs to be lazy
    haskey(ctx, s) && return ctx[s]
    eval(s)
end

interpret!(ctx, ::LineNumberNode) = nothing
interpret!(ctx, n::Number) = n
interpret!(ctx , x) = x

function interpret!(ctx, expr:: Expr)
    # @show ctx, expr
    r(x) = interpret!(ctx,x)

    result = @match expr begin
        Expr(:block, xs...) => for x∈xs r(x) end
        :($f($(args...)))   => r(f)(map(r,args)...)
        :($k = $v)          => (ctx[k] = r(v))
        :(return $x)        => (return x)
        x                   => error("Not yet implemented")
    end

    return result
end

end # module


text:
using Pkg
pkg"activate ."

using Genie, HTTP
using Genie.Router, Genie.Responses

route("/headers") do
  setheaders("X-Foo-Bar" => "Baz")

  "Foo bar baz"
end
route("/headers", method = OPTIONS) do
  setheaders(["X-Foo-Bar" => "Baz", "Access-Control-Allow-Methods" => "GET, POST, OPTIONS"])
  setstatus(200)
end

Genie.AppServer.startup(verbose = true)

response = HTTP.request("GET", "http://localhost:8000/headers") # unhandled, should get default response
@show response

response = HTTP.request("OPTIONS", "http://localhost:8000/headers") # handled
@show response

exit(0)


text:
using FractionalDiffEq, Plots

τ=3.1416; h=0.01; α=0.4; tspan = (0, 70)
x0(t) = [sin(t)*cos(t); sin(t)*cos(t); cos(t)^2-sin(t)^2; cos(t)^2-sin(t)^2]
A=[0 0 1 0; 0 0 0 1; 0 -2 0 0; -2 0 0 0]
B=[0 0 0 0; 0 0 0 0; -2 0 0 0; 0 -2 0 0]
f=[0; 0; 0; 0]

prob = FDDEMatrixProblem(α, τ, A, B, f, x0, tspan)
sol=solve(prob, h, MatrixForm())
plot(sol[:, 1], sol[:, 3])


text:
@testset "Test Constructor for MIMO" begin
    y = rand(2, 100)
    u = rand(3, 100)
    ts = 0.1

    iddata = IdData(y, u, ts)

    @test iddata.t_start == 0
    @test iddata.time_unit == ""
    @test iddata.input_unit == ["", "", ""]
    @test iddata.output_unit == ["", ""]
    @test iddata.input_name == ["u1", "u2", "u3"]
    @test iddata.output_name == ["y1", "y2"]

    @test_throws DimensionMismatch IdData(y[:, 1:end-1], u, ts)
end

@testset "Test Constructor for SISO" begin
    y = rand(100)
    u = rand(100)
    ts = 0.1

    iddata = IdData(y, u, ts)

    @test iddata.t_start == 0
    @test iddata.time_unit == ""
    @test iddata.input_unit == [""]
    @test iddata.output_unit == [""]
    @test iddata.input_name == ["u"]
    @test iddata.output_name == ["y"]

    @test_throws DimensionMismatch IdData(y[1:end-1], u, ts)
end

@testset "Test Constructor for MISO" begin
    y = rand(100)
    u = rand(2, 100)
    ts = 0.1

    iddata = IdData(y, u, ts)

    @test iddata.t_start == 0
    @test iddata.time_unit == ""
    @test iddata.input_unit == ["", ""]
    @test iddata.output_unit == [""]
    @test iddata.input_name == ["u1", "u2"]
    @test iddata.output_name == ["y"]

    @test_throws DimensionMismatch IdData(y[1:end-1], u, ts)
end

@testset "Test Constructor for SIMO" begin
    y = rand(2, 100)
    u = rand(100)
    ts = 0.1

    iddata = IdData(y, u, ts)

    @test iddata.t_start == 0
    @test iddata.time_unit == ""
    @test iddata.input_unit == [""]
    @test iddata.output_unit == ["", ""]
    @test iddata.input_name == ["u"]
    @test iddata.output_name == ["y1", "y2"]

    @test_throws DimensionMismatch IdData(y[:, 1:end-1], u, ts)
end

@testset "Test Get Methods" begin
    y = rand(100)
    u = rand(100)
    ts = 0.1
    t_start = 10

    iddata = IdData(y, u, ts, t_start=t_start)

    @test get_time_vector(iddata) ≈ Vector(10.0:0.1:19.9)
    @test get_output(iddata) ≈ reshape(y, 1, 100)
    @test get_input(iddata) ≈ reshape(u, 1, 100)
end


text:
@testset "Classifier" begin
    # * Setup the problem
    x = [[0.0, 0.0] [1.0, 1.0]]
    y = [-1.0, 1.0]
    svm = LSSVC()

    x_test = [2.0, 2.0]
    x_test = reshape(x_test, 2, :)
    fitted = svmtrain(svm, x, y)
    result = svmpredict(svm, fitted, x_test)
    true_result = [1.0]

    @test all(result .== true_result)

    # ! Test for multiple prediction points
    x_test = [[2.0, 2.0] [3.0, 3.0] [-1.0, -1.0]]
    result = svmpredict(svm, fitted, x_test)
    true_result = [1.0, 1.0, -1.0]

    @test all(result .== true_result)

    # ! Change the hyperparameters
    svm = LSSVC(; γ=5.0, σ=0.5)
    fitted = svmtrain(svm, x, y)
    result = svmpredict(svm, fitted, x_test)
    true_result = [1.0, 1.0, -1.0]

    @test all(result .== true_result)
end

@testset "Regression" begin
    # * Setup the problem
    X_train = [[0.0, 0.0] [2.0, 2.0]]
    y = [0.5, 2.5]
    svr = LSSVR()
    fitted = svmtrain(svr, X_train, y)

    # Create some testing data
    X_test = [1.0, 1.0]
    X_test = reshape(X_test, 2, :)
    y_test = [1.5]
    result = svmpredict(svr, fitted, X_test)
    @test isapprox(result, y_test)
end

function train_and_predict(x, xtest, y, kernel, deg)
    svm = LSSVC(; kernel=kernel, degree=deg)
    fitted = svmtrain(svm, x, y)
    result = svmpredict(svm, fitted, xtest)

    return result
end

@testset "Kernels" begin
    # * Setup the problem
    x = [[0.0, 0.0] [1.0, 1.0]]
    y = [-1.0, 1.0]
    kernels = Dict(:linear => 1, :poly => 2)

    x_test = [2.0, 2.0]
    x_test = reshape(x_test, 2, :)
    true_result = [1.0]

    for (k, v) in kernels
        result = train_and_predict(x, x_test, y, k, v)
        @test all(result .== true_result)
    end
end

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
import Images
using Compat; import Compat.String

# Create a cone in 3d that changes color over time
sz = [201, 301, 31]
center = [(s+1)>>1 for s in sz]  # ceil(Int, sz/2)
C3 = Bool[(i-center[1])^2+(j-center[2])^2 <= k^2 for i = 1:sz[1], j = 1:sz[2], k = sz[3]:-1:1]
cmap1 = round(UInt32, linspace(0,255,60))
cmap = Array(UInt32, length(cmap1))
for i = 1:length(cmap)
    cmap[i] = cmap1[i]<<16 + cmap1[end-i+1]<<8 + cmap1[i]
end
C4 = Array(UInt32, sz..., length(cmap))
for i = 1:length(cmap)
    C4[:,:,:,i] = C3*cmap[i]
end
img = Images.Image(C4, @compat Dict("spatialorder" => ["x", "y", "z"], "timedim" => 4, "colorspace" => "RGB24", "pixelspacing" => [1,1,3]))


text:
module DA_Nlayer_functionTests

using Test
using LightPropagation

# CW comparison of N-layer cylinder to semi-infinite
# note that for large layer thicknesses and radii many more besselroots are needed

### test at multiple SDS

# ρ = 0.0 cm
cylinder_data = Nlayer_cylinder(ρ = 0.0, a = 100.0, l = [10.0, 10.0, 10.0, 10.0], bessels = besselroots[1:10000])
@test fluence_DA_semiinf_CW(0.0, 0.1, 10.0) ≈ fluence_DA_Nlay_cylinder_CW(cylinder_data)

# ρ = 1.0 cm
cylinder_data = Nlayer_cylinder(ρ = 1.0, a = 100.0, l = [10.0, 10.0, 10.0, 10.0], bessels = besselroots[1:10000])
@test fluence_DA_semiinf_CW(1.0, 0.1, 10.0) ≈ fluence_DA_Nlay_cylinder_CW(cylinder_data)

# ρ = 4.0 cm
cylinder_data = Nlayer_cylinder(ρ = 4.0, a = 100.0, l = [10.0, 10.0, 10.0, 10.0], bessels = besselroots[1:10000])
@test fluence_DA_semiinf_CW(4.0, 0.1, 10.0) ≈ fluence_DA_Nlay_cylinder_CW(cylinder_data)

# ρ = 8.0 cm; for long SDS you need more besselroots or decrease the dimensions of cylinder
cylinder_data = Nlayer_cylinder(ρ = 8.0, a = 100.0, l = [10.0, 10.0, 10.0, 10.0], bessels = besselroots[1:10000])
@test isapprox(fluence_DA_semiinf_CW(8.0, 0.1, 10.0), fluence_DA_Nlay_cylinder_CW(cylinder_data), rtol=1e-6)

### test that 2, 3, 4 layers equal with same optical properties
cylinder_4 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0, 10.0, 10.0], μa = [0.1, 0.1, 0.1, 0.1], l = [10.0, 10.0, 10.0, 10.0], n_med = [1.0, 1.0, 1.0, 1.0], bessels = besselroots[1:10000])
cylinder_3 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0, 10.0], μa = [0.1, 0.1, 0.1], l = [10.0, 10.0, 10.0], n_med = [1.0, 1.0, 1.0], bessels = besselroots[1:10000])
cylinder_2 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], bessels = besselroots[1:10000])
@test fluence_DA_Nlay_cylinder_CW(cylinder_4) ≈  fluence_DA_Nlay_cylinder_CW(cylinder_3)  ≈ fluence_DA_Nlay_cylinder_CW(cylinder_2)

### test that arbitrary layers are equal
cylinder_5 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0, 10.0, 10.0, 10.0], μa = [0.1, 0.1, 0.1, 0.1, 0.1], l = [2.0, 2.0, 2.0, 5.0, 5.0], n_med = [1.0, 1.0, 1.0, 1.0, 1.0], bessels = besselroots[1:10000])
cylinder_8 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0], μa = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1], l = [2.0, 2.0, 2.0, 5.0, 5.0, 3.0, 3.0, 3.0], n_med = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], bessels = besselroots[1:10000])

@test fluence_DA_Nlay_cylinder_CW(cylinder_4) ≈  fluence_DA_Nlay_cylinder_CW(cylinder_5) ≈  fluence_DA_Nlay_cylinder_CW(cylinder_8)
### test high scatterring
# SDS = 1.0, μsp = 60.0, once again need more besselroots if you want higher scattering coefficients
cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [60.0, 60.0, 60.0, 60.0], μa = [0.1, 0.1, 0.1, 0.1], l = [2.0, 2.0, 2.0, 2.0], n_med = [1.0, 1.0, 1.0, 1.0], a = 20.0, bessels = besselroots[1:10000])
@test fluence_DA_semiinf_CW(1.0, 0.1, 60.0) ≈ fluence_DA_Nlay_cylinder_CW(cylinder_data)

### test for varying index of refractions
cyl_2 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.5, 1.5], n_ext = 1.4, a = 20.0, bessels = besselroots[1:10000])
cyl_3 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0, 10.0], μa = [0.1, 0.1, 0.1], l = [10.0, 10.0, 10.0], n_med = [1.5, 1.5, 1.5], n_ext = 1.4, a = 20.0, bessels = besselroots[1:10000])
cyl_4 = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0, 10.0, 10.0], μa = [0.1, 0.1, 0.1, 0.1], l = [10.0, 10.0, 10.0, 10.0], n_med = [1.5, 1.5, 1.5, 1.5], n_ext = 1.4, a = 20.0, bessels = besselroots[1:10000])
@test fluence_DA_Nlay_cylinder_CW(cyl_4) ≈ fluence_DA_Nlay_cylinder_CW(cyl_3) ≈ fluence_DA_Nlay_cylinder_CW(cyl_2) ≈ fluence_DA_semiinf_CW(1.0, 0.1, 10.0, n_med = 1.5, n_ext = 1.4)

#### tests for CW fluence in the bottom layer GN
@test fluence_DA_Nlay_cylinder_CW(0.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [1.0, 1.0], 10.0, 2.0, besselroots[1:10000]) ≈ fluence_DA_slab_CW(0.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 2.0, xs = 50)
@test fluence_DA_Nlay_cylinder_CW(1.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [1.0, 1.0], 10.0, 2.0, besselroots[1:10000]) ≈ fluence_DA_slab_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 2.0, xs = 50)
@test fluence_DA_Nlay_cylinder_CW(1.6, [0.21, 0.21], [12.1, 12.1], 1.0, [1.0, 1.0], [1.0, 1.0], 10.0, 2.0, besselroots[1:10000]) ≈ fluence_DA_slab_CW(1.6, 0.21, 12.1; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 2.0, xs = 50)

@test fluence_DA_Nlay_cylinder_CW(0.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [1.5, 1.5], 10.0, 3.0, besselroots[1:10000]) ≈ fluence_DA_slab_CW(0.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 3.0, z = 3.0, xs = 50)

@test fluence_DA_Nlay_cylinder_CW(0.0, [0.1, 0.1, 0.1, 0.1, 0.1, 0.1], [10.0, 10.0, 10.0, 10.0, 10.0, 10.0], 1.0, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [0.5, 0.5, 0.5, 0.5, 0.5, 0.5], 10.0, 3.0, besselroots[1:10000]) ≈ fluence_DA_slab_CW(0.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 3.0, z = 3.0, xs = 50)
@test fluence_DA_Nlay_cylinder_CW(0.0, [0.1, 0.1, 0.1, 0.1, 0.1, 0.1], [10.0, 10.0, 10.0, 10.0, 10.0, 10.0], 1.0, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 10.0, 6.0, besselroots[1:10000]) ≈ fluence_DA_slab_CW(0.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 6.0, z = 6.0, xs = 50)

### test that scalar and vector ρ matches (they call different routines)
data = Nlayer_cylinder(ρ = 1.0:0.1:2.0, bessels = besselroots[1:1000])
@test fluence_DA_Nlay_cylinder_CW(data) ≈ map(ρ -> fluence_DA_Nlay_cylinder_CW(Nlayer_cylinder(ρ = ρ, bessels = besselroots[1:1000])), 1.0:0.1:2.0)

#########
# Time-Domain
##########

## test reflectance in first layer
# test long times at SDS = 0.5 cm
t = range(0.03, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(a = 8.0, ρ = 0.5, bessels = besselroots[1:600])

si = fluence_DA_semiinf_TD(t, 0.5, 0.1, 10.0)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ si
@test a96 ≈ si

# test long times at SDS = 3.0
t = range(0.04, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(a = 8.0, ρ = 3.0, l = [4.0, 4.0, 4.0, 5.0], bessels = besselroots[1:600])

si = fluence_DA_semiinf_TD(t, 3.0, 0.1, 10.0)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ si
@test a96 ≈ si

# test high scattering at SDS = 1.5 cm
t = range(0.04, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(a = 8.0, ρ = 1.5, l = [4.0, 4.0, 4.0, 5.0], μsp = [60.0, 60.0, 60.0, 60.0], bessels = besselroots[1:600])

si = fluence_DA_semiinf_TD(t, 1.5, 0.1, 60.0)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ si
@test a96 ≈ si

## in general the error in the TD is like 1e-14..

## test transmittance in bottom layer
t = range(0.03, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(μa = [0.1, 0.1], μsp = [10.0, 10.0], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 1.0], a = 8.0, ρ = 0.5, z = 2.0, bessels = besselroots[1:600])

slab = fluence_DA_slab_TD(t, 0.5, 0.1, 10.0; s = 2.0, z = 2.0, xs = 50)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ slab
@test a96 ≈ slab

t = range(0.03, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(μa = [0.1, 0.1], μsp = [10.0, 10.0], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.5, 1.5], a = 8.0, ρ = 0.5, z = 3.0, bessels = besselroots[1:600])

slab = fluence_DA_slab_TD(t, 0.5, 0.1, 10.0; s = 3.0, z = 3.0, xs = 50)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ slab
@test a96 ≈ slab

t = range(0.03, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(μa = [0.1, 0.1], μsp = [10.0, 10.0], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.5, 1.5], a = 8.0, ρ = 1.5, z = 3.0, bessels = besselroots[1:600])

slab = fluence_DA_slab_TD(t, 1.5, 0.1, 10.0; s = 3.0, z = 3.0, xs = 50)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ slab
@test a96 ≈ slab

t = range(0.03, 8.0, length = 60)
cylinder_data = Nlayer_cylinder(μa = [0.1, 0.1], μsp = [10.0, 10.0], n_ext = 1.3, n_med = [1.2, 1.2], l = [1.0, 1.0], a = 8.0, ρ = 0.5, z = 2.0, bessels = besselroots[1:600])

slab = fluence_DA_slab_TD(t, 0.5, 0.1, 10.0; s = 2.0, z = 2.0, xs = 50, n_ext = 1.3, n_med = 1.2)

a72 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 72)
a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, N = 96)

@test a72 ≈ slab
@test a96 ≈ slab

### Flux tests

# CW
cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0, bessels = besselroots[1:10000])
@test flux_DA_Nlay_cylinder_CW(cylinder_data) ≈ flux_DA_semiinf_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0)

cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [2.0, 2.0], z = 4.0, n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0, bessels = besselroots[1:10000])
@test flux_DA_Nlay_cylinder_CW(cylinder_data) ≈ flux_DA_slab_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 4.0, z = 4.0, xs = 15)

# TD
t = 1.0
cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0, bessels = besselroots[1:10000])
@test flux_DA_Nlay_cylinder_TD(t, cylinder_data) ≈ flux_DA_semiinf_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0)

cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [2.0, 2.0], z = 4.0, n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0, bessels = besselroots[1:10000])
@test flux_DA_Nlay_cylinder_TD(t, cylinder_data) ≈ flux_DA_slab_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 4.0, z = 4.0, xs = 15)

t = 0.5:0.1:1.5
cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0, bessels = besselroots[1:10000])
@test flux_DA_Nlay_cylinder_TD(t, cylinder_data) ≈ flux_DA_semiinf_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0)

cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [2.0, 2.0], z = 4.0, n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0, bessels = besselroots[1:10000])
@test flux_DA_Nlay_cylinder_TD(t, cylinder_data) ≈ flux_DA_slab_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 4.0, z = 4.0, xs = 15)

end # module

max_stars_repo_name: heltonmc/LightPropagation

max_issues_repo_name: heltonmc/LightPropagation


text:
using Patchwork
using FactCheck

import Patchwork: count, TextNode, Elem

ul(x...) = Elem(:xhtml, :ul, x)
li(x...) = Elem(:xhtml, :li, x)

facts("Testing Nodes") do
    context("testing counts") do
        @fact count(TextNode("a")) --> 0
        @fact count(ul()) --> 0
        @fact count(li("x")) --> 1
        @fact count(ul(li("x"))) --> 2
        @fact count(ul(li("x"), li("x"))) --> 4
    end
end

facts("Testing nested node creation") do
  node = Elem(:div, TextNode("a"))
  @fact count(node) --> 1
end


text:
"""
Supertype for all metadata wrappers.

These allow tracking the contents and origin of metadata.
This can facilitate conversion between metadata types (for saving
a file to a differenet format) or simply saving data back to the
same file type with identical metadata.
"""
abstract type Metadata{K,V} <: AbstractDict{K,V} end

(::Type{T})() where T <: Metadata = T(Dict())

val(metadata::Metadata) = metadata.val

Base.get(metadata::Metadata, args...) = get(val(metadata), args...)
Base.getindex(metadata::Metadata, args...) = getindex(val(metadata), args...)
Base.setindex!(metadata::Metadata, args...) = setindex!(val(metadata), args...)
Base.keys(metadata::Metadata) = keys(val(metadata))
Base.iterate(metadata::Metadata, args...) = iterate(val(metadata), args...)
Base.IteratorSize(::Metadata{M}) where M = Base.IteratorSize(M)
Base.IteratorEltype(::Metadata{M}) where M = Base.IteratorEltype(M)
Base.eltype(::Metadata{M}) where M = eltype(M)
Base.length(metadata::Metadata) = length(val(metadata))

"""
Metadata wrappers to be attached to `Dimension`s.
"""
abstract type DimMetadata{K,V} <: Metadata{K,V} end

"""
Metadata wrappers to be attached to `AbstractGeoArrays`.
"""
abstract type ArrayMetadata{K,V} <: Metadata{K,V} end

"""
Metadata wrappers to be attached to `AbstractGeoStack`.
"""
abstract type StackMetadata{K,V} <: Metadata{K,V} end


text:
if Pkg.installed("ZipFile") == nothing
    Pkg.add("ZipFile")
end

using Mimi
using ZipFile
using Compat

function unzip(inputfilename, outputpath=pwd())
    r = ZipFile.Reader(inputfilename)
    try
        for f in r.files
            outpath = joinpath(outputpath, f.name)
            if isdirpath(outpath)
                mkpath(outpath)
            else
                Base.open(outpath, "w") do io
                    write(io, readstring(f))
                end
            end
        end
    finally
        close(r)
    end
    nothing
end

#list of URLs of branches of packages to test
dependencies = [
"https://github.com/davidanthoff/fund.jl/archive/master.zip"
]

#list of failed tests to build as you go
errors = []
#make a temporary directory to run the tests in
tmp_path = joinpath(dirname(@__FILE__),"tmp_testing/")
mkdir(tmp_path)

#loop through each dependent package
for url in dependencies
  #download the package
  zip_name = chomp(basename(url))
  zip_file_path = joinpath(tmp_path, zip_name)
  download(url, zip_file_path)
  unzip(zip_file_path, tmp_path)
  rm(zip_file_path)
  #find the name of the unzipped package (this only works if the zip archive only has one directory, the package)
  package_name = readdir(tmp_path)[1]

  #test the package
  try
    process = string(tmp_path, package_name, "/test/runtests.jl")
    run(`julia $process`)
  catch e
    append!(errors, [(package_name, e)])
  end
  #delete current package before testing next one
  rm(joinpath(tmp_path, package_name), recursive=true)
end

#remove the temporary directory
rm(tmp_path, recursive=true)

#report the errors that occurred
num_errors = length(errors)
error_message = "Failed tests: $num_errors"

for (package_name, error) in errors
  error_message = string(error_message, "\n", "error in $package_name:", error)
end

if num_errors > 0
  error(error_message)
else
  println("All dependency tests passed.")
end


text:
#==============================================================================#
# TranscribeService.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/transcribe-2017-10-26.normal.json
#==============================================================================#

__precompile__()

module TranscribeService

using AWSCore

"""
    using AWSSDK.TranscribeService.create_vocabulary
    create_vocabulary([::AWSConfig], arguments::Dict)
    create_vocabulary([::AWSConfig]; VocabularyName=, LanguageCode=, Phrases=)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "CreateVocabulary", arguments::Dict)
    transcribe([::AWSConfig], "CreateVocabulary", VocabularyName=, LanguageCode=, Phrases=)

# CreateVocabulary Operation

Creates a new custom vocabulary that you can use to change the way Amazon Transcribe handles transcription of an audio file.

# Arguments

## `VocabularyName = ::String` -- *Required*
The name of the vocabulary. The name must be unique within an AWS account. The name is case-sensitive.

## `LanguageCode = "en-US" or "es-US"` -- *Required*
The language code of the vocabulary entries.

## `Phrases = [::String, ...]` -- *Required*
An array of strings that contains the vocabulary entries.

# Returns

`CreateVocabularyResponse`

# Exceptions

`BadRequestException`, `LimitExceededException`, `InternalFailureException` or `ConflictException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/CreateVocabulary)
"""
@inline create_vocabulary(aws::AWSConfig=default_aws_config(); args...) = create_vocabulary(aws, args)

@inline create_vocabulary(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "CreateVocabulary", args)

@inline create_vocabulary(args) = create_vocabulary(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.delete_vocabulary
    delete_vocabulary([::AWSConfig], arguments::Dict)
    delete_vocabulary([::AWSConfig]; VocabularyName=)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "DeleteVocabulary", arguments::Dict)
    transcribe([::AWSConfig], "DeleteVocabulary", VocabularyName=)

# DeleteVocabulary Operation

Deletes a vocabulary from Amazon Transcribe.

# Arguments

## `VocabularyName = ::String` -- *Required*
The name of the vocabulary to delete.

# Exceptions

`NotFoundException`, `LimitExceededException`, `BadRequestException` or `InternalFailureException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/DeleteVocabulary)
"""
@inline delete_vocabulary(aws::AWSConfig=default_aws_config(); args...) = delete_vocabulary(aws, args)

@inline delete_vocabulary(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "DeleteVocabulary", args)

@inline delete_vocabulary(args) = delete_vocabulary(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.get_transcription_job
    get_transcription_job([::AWSConfig], arguments::Dict)
    get_transcription_job([::AWSConfig]; TranscriptionJobName=)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "GetTranscriptionJob", arguments::Dict)
    transcribe([::AWSConfig], "GetTranscriptionJob", TranscriptionJobName=)

# GetTranscriptionJob Operation

Returns information about a transcription job. To see the status of the job, check the `TranscriptionJobStatus` field. If the status is `COMPLETED`, the job is finished and you can find the results at the location specified in the `TranscriptionFileUri` field.

# Arguments

## `TranscriptionJobName = ::String` -- *Required*
The name of the job.

# Returns

`GetTranscriptionJobResponse`

# Exceptions

`BadRequestException`, `LimitExceededException`, `InternalFailureException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/GetTranscriptionJob)
"""
@inline get_transcription_job(aws::AWSConfig=default_aws_config(); args...) = get_transcription_job(aws, args)

@inline get_transcription_job(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "GetTranscriptionJob", args)

@inline get_transcription_job(args) = get_transcription_job(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.get_vocabulary
    get_vocabulary([::AWSConfig], arguments::Dict)
    get_vocabulary([::AWSConfig]; VocabularyName=)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "GetVocabulary", arguments::Dict)
    transcribe([::AWSConfig], "GetVocabulary", VocabularyName=)

# GetVocabulary Operation

Gets information about a vocabulary.

# Arguments

## `VocabularyName = ::String` -- *Required*
The name of the vocabulary to return information about. The name is case-sensitive.

# Returns

`GetVocabularyResponse`

# Exceptions

`NotFoundException`, `LimitExceededException`, `InternalFailureException` or `BadRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/GetVocabulary)
"""
@inline get_vocabulary(aws::AWSConfig=default_aws_config(); args...) = get_vocabulary(aws, args)

@inline get_vocabulary(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "GetVocabulary", args)

@inline get_vocabulary(args) = get_vocabulary(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.list_transcription_jobs
    list_transcription_jobs([::AWSConfig], arguments::Dict)
    list_transcription_jobs([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "ListTranscriptionJobs", arguments::Dict)
    transcribe([::AWSConfig], "ListTranscriptionJobs", <keyword arguments>)

# ListTranscriptionJobs Operation

Lists transcription jobs with the specified status.

# Arguments

## `Status = "IN_PROGRESS", "FAILED" or "COMPLETED"`
When specified, returns only transcription jobs with the specified status.

## `JobNameContains = ::String`
When specified, the jobs returned in the list are limited to jobs whose name contains the specified string.

## `NextToken = ::String`
If the result of the previous request to `ListTranscriptionJobs` was truncated, include the `NextToken` to fetch the next set of jobs.

## `MaxResults = ::Int`
The maximum number of jobs to return in the response. If there are fewer results in the list, this response contains only the actual results.

# Returns

`ListTranscriptionJobsResponse`

# Exceptions

`BadRequestException`, `LimitExceededException` or `InternalFailureException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/ListTranscriptionJobs)
"""
@inline list_transcription_jobs(aws::AWSConfig=default_aws_config(); args...) = list_transcription_jobs(aws, args)

@inline list_transcription_jobs(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "ListTranscriptionJobs", args)

@inline list_transcription_jobs(args) = list_transcription_jobs(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.list_vocabularies
    list_vocabularies([::AWSConfig], arguments::Dict)
    list_vocabularies([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "ListVocabularies", arguments::Dict)
    transcribe([::AWSConfig], "ListVocabularies", <keyword arguments>)

# ListVocabularies Operation

Returns a list of vocabularies that match the specified criteria. If no criteria are specified, returns the entire list of vocabularies.

# Arguments

## `NextToken = ::String`
If the result of the previous request to `ListVocabularies` was truncated, include the `NextToken` to fetch the next set of jobs.

## `MaxResults = ::Int`
The maximum number of vocabularies to return in the response. If there are fewer results in the list, this response contains only the actual results.

## `StateEquals = "PENDING", "READY" or "FAILED"`
When specified, only returns vocabularies with the `VocabularyState` field equal to the specified state.

## `NameContains = ::String`
When specified, the vocabularies returned in the list are limited to vocabularies whose name contains the specified string. The search is case-insensitive, `ListVocabularies` will return both "vocabularyname" and "VocabularyName" in the response list.

# Returns

`ListVocabulariesResponse`

# Exceptions

`BadRequestException`, `LimitExceededException` or `InternalFailureException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/ListVocabularies)
"""
@inline list_vocabularies(aws::AWSConfig=default_aws_config(); args...) = list_vocabularies(aws, args)

@inline list_vocabularies(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "ListVocabularies", args)

@inline list_vocabularies(args) = list_vocabularies(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.start_transcription_job
    start_transcription_job([::AWSConfig], arguments::Dict)
    start_transcription_job([::AWSConfig]; TranscriptionJobName=, LanguageCode=, MediaFormat=, Media=, <keyword arguments>)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "StartTranscriptionJob", arguments::Dict)
    transcribe([::AWSConfig], "StartTranscriptionJob", TranscriptionJobName=, LanguageCode=, MediaFormat=, Media=, <keyword arguments>)

# StartTranscriptionJob Operation

Starts an asynchronous job to transcribe speech to text.

# Arguments

## `TranscriptionJobName = ::String` -- *Required*
The name of the job. You can't use the strings "." or ".." in the job name. The name must be unique within an AWS account.

## `LanguageCode = "en-US" or "es-US"` -- *Required*
The language code for the language used in the input media file.

## `MediaSampleRateHertz = ::Int`
The sample rate, in Hertz, of the audio track in the input media file.

## `MediaFormat = "mp3", "mp4", "wav" or "flac"` -- *Required*
The format of the input media file.

## `Media = ["MediaFileUri" =>  ::String]` -- *Required*
An object that describes the input media for a transcription job.

## `OutputBucketName = ::String`
The location where the transcription is stored.

If you set the `OutputBucketName`, Amazon Transcribe puts the transcription in the specified S3 bucket. When you call the [GetTranscriptionJob](@ref) operation, the operation returns this location in the `TranscriptFileUri` field. The S3 bucket must have permissions that allow Amazon Transcribe to put files in the bucket. For more information, see [Permissions Required for IAM User Roles](https://docs.aws.amazon.com/transcribe/latest/dg/access-control-managing-permissions.html#auth-role-iam-user).

If you don't set the `OutputBucketName`, Amazon Transcribe generates a pre-signed URL, a shareable URL that provides secure access to your transcription, and returns it in the `TranscriptFileUri` field. Use this URL to download the transcription.

## `Settings = [ ... ]`
A `Settings` object that provides optional settings for a transcription job.
```
 Settings = [
        "VocabularyName" =>  ::String,
        "ShowSpeakerLabels" =>  ::Bool,
        "MaxSpeakerLabels" =>  ::Int,
        "ChannelIdentification" =>  ::Bool
    ]
```

# Returns

`StartTranscriptionJobResponse`

# Exceptions

`BadRequestException`, `LimitExceededException`, `InternalFailureException` or `ConflictException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/StartTranscriptionJob)
"""
@inline start_transcription_job(aws::AWSConfig=default_aws_config(); args...) = start_transcription_job(aws, args)

@inline start_transcription_job(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "StartTranscriptionJob", args)

@inline start_transcription_job(args) = start_transcription_job(default_aws_config(), args)

"""
    using AWSSDK.TranscribeService.update_vocabulary
    update_vocabulary([::AWSConfig], arguments::Dict)
    update_vocabulary([::AWSConfig]; VocabularyName=, LanguageCode=, Phrases=)

    using AWSCore.Services.transcribe
    transcribe([::AWSConfig], "UpdateVocabulary", arguments::Dict)
    transcribe([::AWSConfig], "UpdateVocabulary", VocabularyName=, LanguageCode=, Phrases=)

# UpdateVocabulary Operation

Updates an existing vocabulary with new values. The `UpdateVocabulary` operation overwrites all of the existing information with the values that you provide in the request.

# Arguments

## `VocabularyName = ::String` -- *Required*
The name of the vocabulary to update. The name is case-sensitive.

## `LanguageCode = "en-US" or "es-US"` -- *Required*
The language code of the vocabulary entries.

## `Phrases = [::String, ...]` -- *Required*
An array of strings containing the vocabulary entries.

# Returns

`UpdateVocabularyResponse`

# Exceptions

`BadRequestException`, `LimitExceededException`, `InternalFailureException`, `NotFoundException` or `ConflictException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/transcribe-2017-10-26/UpdateVocabulary)
"""
@inline update_vocabulary(aws::AWSConfig=default_aws_config(); args...) = update_vocabulary(aws, args)

@inline update_vocabulary(aws::AWSConfig, args) = AWSCore.Services.transcribe(aws, "UpdateVocabulary", args)

@inline update_vocabulary(args) = update_vocabulary(default_aws_config(), args)

end # module TranscribeService

#==============================================================================#
# End of file
#==============================================================================#


text:
using Statistics

function format(x)
    if abs(x) < 999
        return @sprintf("%12.7f", x)
    else
        return @sprintf("%12.5e", x)
    end
end

import Base.write
"""

```
write(R::ComplexMixtures.Result, filename::String, solute::Selection, solvent::Selection)
```

Function to write the final results to output files as simple tables
that are human-readable and easy to analyze with other software

If the solute and solvent selections are provides, pass on the atom
names.

"""
write(R::Result, filename::String, solute::Selection, solvent::Selection) =
    write(R, filename, solute_names = solute.names, solvent_names = solvent.names)

"""

```
write(R::ComplexMixtures.Result, filename::String;
      solute_names::Vector{String} = ["nothing"],
      solvent_names::Vector{String} = ["nothing"])
```

Optional passing of atom names.

"""
function write(
    R::Result,
    filename::String;
    solute_names::Vector{String} = ["nothing"],
    solvent_names::Vector{String} = ["nothing"],
)

    # Names of output files containing atomic contibutions
    atom_contrib_solvent =
        FileOperations.remove_extension(filename) *
        "-ATOM_CONTRIB_SOLVENT." *
        FileOperations.file_extension(filename)
    atom_contrib_solute =
        FileOperations.remove_extension(filename) *
        "-ATOM_CONTRIB_SOLUTE." *
        FileOperations.file_extension(filename)

    #
    # GMD computed with minimum distance
    #

    # Conversion factor for volumes (as KB integrals), from A^3 to cm^3/mol
    mole = 6.022140857e23
    convert = mole / 1.e24

    output = open(filename, "w")
    println(output, @sprintf("#"))
    println(output, @sprintf("# Output of ComplexMixtures - MDDF"))
    println(output, @sprintf("# Trajectory files and weights:"))
    for i = 1:length(R.files)
        println(output, "#  $(R.files[i]) - w = $(R.weights[i])")
    end
    println(output, @sprintf("#"))
    println(
        output,
        @sprintf(
            "# Density of solvent in simulation box (sites/A^3): %15.8f",
            R.density.solvent
        )
    )
    println(
        output,
        @sprintf(
            "# Density of solvent in bulk (estimated) (sites/A^3): %15.8f",
            R.density.solvent_bulk
        )
    )
    println(
        output,
        @sprintf(
            "# Molar volume of solvent in simulation (cc/mol): %15.8f",
            convert / R.density.solvent
        )
    )
    println(
        output,
        @sprintf(
            "# Molar volume of solvent in bulk (estimated) (cc/mol): %15.8f",
            convert / R.density.solvent_bulk
        )
    )
    println(output, @sprintf("#"))
    println(output, @sprintf("# Number of atoms solute: %i9", size(R.solute_atom, 2)))
    println(
        output,
        @sprintf("# Number of atoms of the solvent: %i9", size(R.solvent_atom, 2))
    )
    println(output, @sprintf("#"))

    if R.options.usecutoff
        ibulk = setbin(R.options.dbulk, R.options.binstep)
    else
        ibulk = round(Int, R.nbins - 1 / R.options.binstep)
    end
    bulkerror = Statistics.mean(R.mddf[ibulk:R.nbins])
    sdbulkerror = Statistics.std(R.mddf[ibulk:R.nbins])
    println(output, "#")
    println(
        output,
        @sprintf(
            "# Average and standard deviation of bulk-gmd: %12.5f +/- %12.5f",
            bulkerror,
            sdbulkerror
        )
    )
    println(output, "#")

    # Output table

    println(output, "# COLUMNS CORRESPOND TO: ")
    println(output, "#       1  Minimum distance to solute (dmin)")
    println(
        output,
        "#       2  GMD distribution (md count normalized by md count of random-solute distribution)",
    )
    println(
        output,
        "#       3  Kirwood-Buff integral (cc/mol) computed [(1/bulkdensity)*(col(6)-col(7))].",
    )
    println(output, "#       4  Minimum distance site count for each dmin.")
    println(
        output,
        "#       5  Minimum distance site count for each dmin for random solute distribution.",
    )
    println(
        output,
        "#       6  Cumulative number of molecules within dmin in the simulation.",
    )
    println(
        output,
        "#       7  Cumulative number of molecules within dmin for random solute distribution.",
    )
    println(output, "#       8  Volume of the shell of distance dmin and width binstep.")
    println(output, "#")
    println(
        output,
        "#   1-DISTANCE         2-GMD      3-KB INT    4-MD COUNT  5-COUNT RAND      6-SUM MD    7-SUM RAND   8-SHELL VOL",
    )

    for i = 1:R.nbins
        line = "  " * format(R.d[i])                                 #  1-DISTANCE
        line = line * "  " * format(R.mddf[i])                         #  2-GMD
        line = line * "  " * format(R.kb[i])                           #  3-KB INT
        line = line * "  " * format(R.md_count[i])                     #  4-MD COUNT
        line = line * "  " * format(R.md_count_random[i])              #  5-COUNT RAND
        line = line * "  " * format(R.sum_md_count[i])                 #  6-SUM MD
        line = line * "  " * format(R.sum_md_count_random[i])          #  7-SUM RAND
        line = line * "  " * format(R.volume.shell[i])                 #  8-SHELL VOL
        println(output, line)
    end
    close(output)

    # Writting gmd per atom contributions for the solvent

    output = open(atom_contrib_solvent, "w")
    println(output, "#")
    println(output, "# Solvent atomic contributions to total MDDF. ")
    println(output, "#")
    println(output, "# Trajectory files: ", R.files)
    println(output, "#")
    println(output, "# Atoms: ")
    for i = 1:size(R.solvent_atom, 2)
        if solvent_names[1] == "nothing"
            println(output, @sprintf("# %9i", i))
        else
            println(output, @sprintf("# %9i %5s", i, solvent_names[i]))
        end
    end
    println(output, "#")
    string = "#   DISTANCE     GMD TOTAL"
    for i = 1:size(R.solvent_atom, 2)
        string = string * @sprintf("  %12i", i)
    end
    println(output, string)
    for i = 1:R.nbins
        string = format(R.d[i])
        string = string * "  " * format(R.mddf[i])
        for j = 1:size(R.solvent_atom, 2)
            string = string * "  " * format(R.solvent_atom[i, j])
        end
        println(output, string)
    end
    close(output)

    # Writting gmd per atom contributions for the solute

    output = open(atom_contrib_solute, "w")
    println(output, "#")
    println(output, "# Solute atomic contributions to total MDDF. ")
    println(output, "#")
    println(output, "# Trajectory files: ", R.files)
    println(output, "#")
    println(output, "# Atoms:")
    for i = 1:size(R.solute_atom, 2)
        if solute_names[1] == "nothing"
            println(output, @sprintf("# %9i", i))
        else
            println(output, @sprintf("# %9i %5s", i, solute_names[i]))
        end
    end
    println(output, "#")
    string = "#   DISTANCE      GMD TOTAL"
    for i = 1:size(R.solute_atom, 2)
        string = string * @sprintf("  %12i", i)
    end
    println(output, string)
    for i = 1:R.nbins
        string = format(R.d[i]) * "  " * format(R.mddf[i])
        for j = 1:size(R.solute_atom, 2)
            string = string * "  " * format(R.solute_atom[i, j])
        end
        println(output, string)
    end
    close(output)

    # Write final messages with names of output files and their content

    println()
    println(" OUTPUT FILES: ")
    println()
    println(" Wrote solvent atomic GMD contributions to file: ", atom_contrib_solvent)
    println(" Wrote solute atomic GMD contributions to file: ", atom_contrib_solute)
    println()
    println(" Wrote main output file: ", filename)
    println()

    return nothing
end


text:
function debugging end
enable_debug() = @eval debugging() = true
disable_debug() = @eval debugging() = false
# enable_debug()
disable_debug()

macro _assert(args...)
    call = Expr(
        :macrocall,
        getfield(ArgCheck, Symbol("@check")),  # injecting the callable for esc
        __source__,
        args...,
    )
    ex = Expr(:block, __source__, call)
    quote
        if $debugging()
            $ex
        end
        nothing
    end |> esc
end

pause() = ccall(:jl_cpu_pause, Cvoid, ())

function halve(xs::AbstractVector)
    mid = length(xs) ÷ 2
    left = @view xs[firstindex(xs):firstindex(xs)-1+mid]
    right = @view xs[firstindex(xs)+mid:end]
    return (left, right)
end

struct Partitioned{T}
    xs::T
    nchunks::Int
    chunksize::Int
end

Base.length(p::Partitioned) = p.nchunks
Base.eachindex(p::Partitioned) = 1:length(p)

function Base.getindex(p::Partitioned, i::Integer)
    xs = p.xs
    nchunks = p.nchunks
    chunksize = p.chunksize
    nextras = length(xs) - chunksize * nchunks
    nsmalls = nchunks - nextras
    # Produce small chunks first since they are the dependencies of latter chunks:
    if i <= nsmalls
        a = chunksize * (i - 1)
        b = a + chunksize - 1
    else
        a = chunksize * nsmalls + (chunksize + 1) * (i - 1 - nsmalls)
        b = a + (chunksize + 1) - 1
    end
    @view xs[begin+a:begin+b]
end

Base.iterate(p::Partitioned, i = 1) = i <= length(p) ? (p[i], i + 1) : nothing

function partitionto(xs, nchunks::Integer)
    @argcheck length(xs) >= nchunks > 0
    chunksize = length(xs) ÷ nchunks
    return Partitioned(xs, nchunks, chunksize)
end

function consecutive(xs)
    ys = Iterators.accumulate(xs; init = ()) do acc, x
        if acc isa Tuple{Any,Any}
            (acc[2], x)
        else
            (acc..., x)
        end
    end
    return Iterators.filter(x -> x isa Tuple{Any,Any}, ys)
end

donothing() = nothing

waitall(tasks) = waitall(donothing, tasks)
function waitall(onerror, tasks)
    err = sync_end!(onerror, tasks)
    err === nothing || throw(err)
    return
end

function sync_end!(onerror, tasks, err = nothing)
    called = false
    for t in tasks
        try
            wait(t)
        catch e
            called || onerror()
            called = true
            if err === nothing
                err = CompositeException()
            end
            push!(err, e)
        end
    end
    return err
end

function _dac_spawn_foreach(f, xs, onerror)
    if length(xs) == 1
        try
            f(xs[1])
        catch
            onerror()
            rethrow()
        end
    else
        l, r = halve(xs)
        t = Threads.@spawn _dac_spawn_foreach(f, r, onerror)
        try
            _dac_spawn_foreach(f, l, onerror)
        finally
            wait(t)
        end
    end
end

function dac_spawn_foreach(f::F, xs; onerror = donothing) where {F}
    isempty(xs) || _dac_spawn_foreach(f, xs, onerror)
    return
end

const spawn_foreach = dac_spawn_foreach


text:
@cdef :PyList_New PyPtr (Py_ssize_t,)
@cdef :PyList_Append Cint (PyPtr, PyPtr)
@cdef :PyList_AsTuple PyPtr (PyPtr,)

const PyList_Type__ref = Ref(PyNULL)
PyList_Type() = pyglobal(PyList_Type__ref, :PyList_Type)

PyList_Check(o) = Py_TypeCheckFast(o, Py_TPFLAGS_LIST_SUBCLASS)

PyList_CheckExact(o) = Py_TypeCheckExact(o, PyList_Type())

PyList_From(xs::Union{Tuple,AbstractVector}) = PyList_FromIter(xs)

PyList_FromIter(xs) = begin
    r = PyList_New(0)
    isnull(r) && return PyNULL
    try
        for (i, x) in enumerate(xs)
            xo = PyObject_From(x)
            isnull(xo) && (Py_DecRef(r); return PyNULL)
            err = PyList_Append(r, xo)
            Py_DecRef(xo)
            ism1(err) && (Py_DecRef(r); return PyNULL)
        end
        return r
    catch err
        Py_DecRef(r)
        PyErr_SetString(PyExc_Exception(), "Julia error: $err")
        return PyNULL
    end
end


text:
#=
small_scripts:
- Julia version: 1.5.1
- Author: Berk
- Date: 2020-11-04
=#

"""
    data_to_Dict(data::Union{Dict, DataFrame, DataFrameRow}, model::JuMP.Model)

Turns data into a Dict, with JuMP.VariableRef subs.
"""
function data_to_Dict(data::Union{Dict, DataFrame, DataFrameRow}, model::JuMP.Model)
    if data isa Dict
        return Dict(fetch_variable(model, key) => val for (key, val) in data)
    else
        colnames = names(data)
        if size(data, 1) == 1
            return Dict(fetch_variable(model, name) => data[!, name][1] for name in colnames)
        else
            return Dict(fetch_variable(model, name) => data[!, name] for name in colnames)
        end
    end
end

""" Wrapper around Iterators.flatten for variables. """
flat(arr) = Array{Any, 1}(collect(Iterators.flatten(arr)))

""" For substitution into expressions. IMPORTANT. """
function substitute(e::Expr, pair)
    MacroTools.prewalk(e) do s
        s == pair.first && return pair.second
        s
    end
end

""" power function from gams """
power(var, num) = var^num

""" sqr function from gams. """
sqr(var) = var^2

""" Checks outer-boundedness of values of a Dict. """
function check_bounds(bounds::Dict)
    if any(isinf.(Iterators.flatten(values(bounds))))
        throw(OCTHaGOnException("Unbounded variables in model!"))
    else
        return
    end
end

"""
    get_varmap(expr_vars::Array, vars::Array)

Helper function to map vars to flatvars.
Arguments:
    flatvars is a flattened Array{JuMP.VariableRef}
    vars is the unflattened version, usually derived from an Expr.
Returns:
    Dict of ID maps
"""
function get_varmap(expr_vars::Array, vars::Array)
    length(flat(expr_vars)) >= length(vars) || throw(OCTHaGOnException(string("Insufficiently many input
        variables declared in ", vars, ".")))
    unique(vars) == vars || throw(OCTHaGOnException(string("Nonunique variables among ", vars, ".")))
    if expr_vars == vars
        return collect(1:length(vars))
    end
    varmap = Tuple[(0,0) for i=1:length(vars)]
    for i = 1:length(expr_vars)
        if expr_vars[i] isa JuMP.VariableRef
            try
                varmap[findall(x -> x == expr_vars[i], vars)[1]] = (i,0)
            catch
                throw(OCTHaGOnException(string("Scalar variable ", expr_vars[i], " was not properly declared in vars.")))
            end
        else
            for j=1:length(expr_vars[i])
                try
                    varmap[findall(x -> x == expr_vars[i][j], vars)[1]] = (i,j)
                catch
                    continue
                end
            end
        end
    end
    length(varmap) == length(vars) || throw(OCTHaGOnException(string("Could not properly map
                                            expr_vars: ", expr_vars,
                                            "to vars: ", vars, ".")))
    return varmap
end

get_varmap(expr_vars::Nothing, vars::Array) = nothing

""" Returns the mapping from flattened expr\\_vars to vars. """
function get_datamap(expr_vars::Array, vars::Array)
    datamap = []
    flatvars = flat(expr_vars)
    for var in vars
        idx = findall(x -> x == var, flatvars)
        if length(idx) == 1
            push!(datamap, idx[1])
        else
            throw(OCTHaGOnException("There was an issue with getting data mapping."))
        end
    end
    return datamap
end

"""Returns the relevant index ranges for variables in expr\\_vars."""
function get_var_ranges(expr_vars::Array)
    var_ranges = []
    count = 0
    for varlist in expr_vars
        if varlist isa VariableRef
            count += 1
            push!(var_ranges, count)
        else
            push!(var_ranges, (count + 1 : count + length(varlist)))
            count += length(varlist)
        end
    end
    return var_ranges
end

"""
    zeroarray(varmap::Array)

Creates a template array for deconstruct function.
"""
function zeroarray(var_ranges::Array)
    arr = []
    for i in var_ranges
        if i isa UnitRange || i isa Array
            push!(arr, zeros(length(i)))
        elseif i isa Int64
            push!(arr, 0)
        end
    end
    return arr
end

"""
    deconstruct(data::DataFrame, vars::Array, varmap::Array)

Takes in data for input into a Function, and rips it apart into appropriate arrays.
"""
function deconstruct(data::DataFrame, vars::Array, expr_vars::Array, varmap::Array)
    n_samples, n_vars = size(data)
    zeroarr = zeroarray(get_var_ranges(expr_vars))
    arrs = [];
    stringvars = string.(vars)
    for i = 1:n_samples
        narr = deepcopy(zeroarr)
        for j = 1:length(varmap)
            if varmap[j] isa Tuple && varmap[j][2] != 0
                narr[varmap[j][1]][varmap[j][2]] = data[i, stringvars[j]]
            else
                narr[varmap[j][1]] = data[i, stringvars[j]]
            end
        end
        push!(arrs, narr)
    end
    return arrs
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CLIME_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CLIME")
JLLWrappers.@generate_main_file("CLIME", UUID("3c6ae550-c37b-5556-a07e-d40b4910cf1c"))
end  # module CLIME_jll


text:
# funexpr.jl - extract arguments and (sanitized) expression of a function body.
#
# Here sanitization means removing things like LineNumberNode and replacing
# hard-to-consume nodes (e.g. `GloabalRef` and `return`) with their simple
# counterparts (e.g. expression with `.` and variable node).

sanitize(x) = x
sanitize(ex::Expr) = sanitize(ExH(ex))
sanitize(ex::LineNumberNode) = nothing
sanitize(m::Module) = Symbol(string(m))
sanitize(ex::ExH{:line}) = nothing
sanitize(ex::ExH{:return}) = sanitize(ex.args[1])
sanitize(x::Core.SSAValue) = Symbol("SSAValue_$(x.id)")
# sanitize(ex::ExH{:macrocall}) = nothing  # note: ignoring macros, experimental

function sanitize(ex::ExH{:block})
    sanitized_args = [sanitize(arg) for arg in ex.args]
    new_args = filter(arg -> arg != nothing, sanitized_args)
    return length(new_args) == 1 ? new_args[1] : Expr(ex.head, new_args...)
end

function sanitize(ex::ExH{:quote})
    return length(ex.args) == 1 ? QuoteNode(ex.args[1]) : ex
end

function sanitize(ex::ExH{H}) where H
    sanitized_args = [sanitize(arg) for arg in ex.args]
    new_args = filter(arg -> arg != nothing, sanitized_args)
    return Expr(H, new_args...)
end

function sanitize(ref::GlobalRef)
    return Expr(:., Symbol(string(ref.mod)), QuoteNode(ref.name))
end

## recover lowered

const RECOVER_LOWERED_RULES = [
    :(Base.broadcast(_f, _xs...)) => :(_f.(_xs...)),
    :(Base.broadcast(_m._f, _xs...)) => :(_m._f.(_xs...)),
    :(A_mul_B(_x, _y)) => :(_x * _y),
    :(A_mul_Bt(_x, _y)) => :(_x * _y'),
    :(At_mul_B(_x, _y)) => :(_x' * _y),
    :(A_mul_Bc(_x, _y)) => :(_x * _y'),
    :(Ac_mul_B(_x, _y)) => :(_x' * _y),
    :(Base.literal_pow(Main.:^, _x, _v)) => :(_x ^ _v),
    :(Core.apply_type(Base.Val, _x)) => :_x,
    # :(Core.getfield(_x, _y)) => :(_x._y),
    :(Core.getfield(_x, $(QuoteNode(:_y)))) => :(_x._y),
]

"""
Try to recover an expression from a lowered form. Example:

    ex = (Main.sum)((Base.literal_pow)(Main.^, (Base.broadcast)(Main.-, (Main.predict)(W, b, x), y), (Core.apply_type)(Base.Val, 2)))
"""
recover_lowered(x) = x
recover_lowered(ex::Expr) = recover_lowered(ExH(ex))

function recover_lowered(ex::ExH{:block})
    recovered_args = [recover_lowered(arg) for arg in ex.args]
    new_args = filter(arg -> arg != nothing, recovered_args)
    return length(new_args) == 1 ? new_args[1] : Expr(ex.head, new_args...)
end

function recover_lowered(ex::ExH{:call})
    # recover arguments
    recovered_args = [recover_lowered(arg) for arg in ex.args[2:end]]
    new_args = filter(arg -> arg != nothing, recovered_args)
    ex = Expr(:call, ex.args[1], new_args...)
    # check patterns
    for (pat, rpat) in RECOVER_LOWERED_RULES
        rex = tryrewrite(ex, pat, rpat)
        if rex != nothing
            ex = get(rex)
            break
        end
    end
    ex = canonical_calls(@__MODULE__, ex) |> subs_bcast_with_dot
    return ex
end

function recover_lowered(ex::ExH{H}) where H
    recovered_args = [recover_lowered(arg) for arg in ex.args]
    new_args = filter(arg -> arg != nothing, recovered_args)
    return Expr(H, new_args...)
end

function recover_lowered_rec(ex)
    ex_old = ex
    ex = recover_lowered(ex)
    # TODO: create a function for this
    while string(ex) != string(ex_old)  # no more changes
        ex_old = ex
        ex = recover_lowered(ex)
    end
    return ex
end

## funexpr

function replace_slots(ex::Expr, slotnames::Vector)
    new_args = Array{Any}(undef, length(ex.args))
    for (i, arg) in enumerate(ex.args)
        if isa(arg, Slot)
            new_args[i] = slotnames[arg.id]
        elseif isa(arg, Expr)
            new_args[i] = replace_slots(arg, slotnames)
        else
            new_args[i] = arg
        end
    end
    new_ex = Expr(ex.head, new_args...)
    return new_ex
end

function arg_names(sig::Expr)
    # note: ignoring type parameters, may not always be the right thing
    while sig.head == :where
        sig = sig.args[1]
    end
    return [isa(arg, Symbol) ? arg : arg.args[1] for arg in sig.args[2:end]]
end

function arg_types(sig::Expr)
    # note: ignoring type parameters, may not always be the right thing
    while sig.head == :where
        sig = sig.args[1]
    end
    return [isa(arg, Symbol) ? Any : arg.args[2] for arg in sig.args[2:end]]
end

function to_expr(src::CodeInfo)
    if isa(src.code, Array{Any,1})
        slotnames = [Symbol(name) for name in Base.sourceinfo_slotnames(src)]
        body = Expr(:body, src.code...)
        result = replace_slots(body, slotnames)
        return result
    else
        error("Can't convert CodeInfo to expression: CodeInfo is compressed: $src")
    end
end

function concretise_types(code::Expr, types::NTuple{N, DataType}) where N
    sig_types = arg_types(code.args[1])
    st = Dict(zip(sig_types, types))
    return subs(code, st)
end

"""
Replace all calls to an inner constructor with the corresponding outer constructor
"""
function replace_inner_constr(f, ex::Expr)
    f_name = Meta.parse(string(f))
    constr = :($f_name(_xs...))
    ex = rewrite_all(ex, :(new{_T1, _T2, _T3}(_xs...)), constr)
    ex = rewrite_all(ex, :(new{_T1, _T2}(_xs...)), constr)
    ex = rewrite_all(ex, :(new{_T}(_xs...)), constr)
    ex = rewrite_all(ex, :(new(_xs...)), constr)
    return ex
end

function funexpr(f::Union{Function, DataType, UnionAll}, types::NTuple{N,DataType}) where N
    method = get_method(f, types)
    file = string(method.file)
    linestart = method.line
    try
        ex, _ = get_source_at(file, linestart)
        ex.head == :toplevel && throw(LoadError(file, Int(linestart), "Bad code found"))
        ex = concretise_types(ex, types)
        ex = replace_inner_constr(f, ex)
        return arg_names(ex.args[1]), sanitize(ex.args[2])
    catch err
        if isa(err, LoadError)
            code = code_lowered(f, types)[1]
            args = convert(Vector{Symbol}, code.slotnames[2:end])
            ex = to_expr(code) |> sanitize |> recover_lowered_rec
            # ex = subs(ex, Dict(:new => parse(string(f))))
            return args, ex
        else
            rethrow(err)
        end
    end
end

func_expr = funexpr

function get_or_generate_argnames(f, types)
    try
        args, _ = funexpr(f, types)
        return args
    catch
        return [Symbol("arg$i") for i=1:length(types)]
    end
end


text:
"""
    align(x::Union{Int,Ptr}, [n])

Return aligned memory address with minimum increment. `align` assumes `n` is a
power of 2.
"""
function align end
@inline align(x::Integer) = (x + REGISTER_SIZE-1) & -REGISTER_SIZE
@inline align(x::Ptr{T}, args...) where {T} = reinterpret(Ptr{T}, align(reinterpret(UInt, x), args...))
@inline align(x::Integer, n) = (nm1 = n - 1; (x + nm1) & -n)
@inline align(x::Integer, ::Type{T}) where {T} = align(x, REGISTER_SIZE ÷ sizeof(T))
@inline align(x::Integer, ::Type{T}) where {T<:Union{Float64,Int64,UInt64}} = align(x, REGISTER_SIZE >>> 3)
@inline align(x::Integer, ::Type{T}) where {T<:Union{Float32,Int32,UInt32}} = align(x, REGISTER_SIZE >>> 2)
@inline align(x::Integer, ::Type{T}) where {T<:Union{Float16,Int16,UInt16}} = align(x, REGISTER_SIZE >>> 1)
@inline align(x::Integer, ::Type{T}) where {T<:Union{Int8,UInt8}} = align(x, REGISTER_SIZE)

# @generated align(::Val{L}, ::Type{T}) where {L,T} = align(L, T)
aligntrunc(x::Integer, n) = x & -n
aligntrunc(x::Integer) = aligntrunc(x, REGISTER_SIZE)
aligntrunc(x::Integer, ::Type{T}) where {T} = aligntrunc(x, REGISTER_SIZE ÷ sizeof(T))
alignment(x::Integer, N = 64) = reinterpret(Int, x) % N

function valloc(N::Int, ::Type{T} = Float64) where {T}
    # We want alignment to both vector and cacheline-sized boundaries
    a = max(REGISTER_SIZE, L₁CACHE.linesize)
    reinterpret(Ptr{T}, align(reinterpret(UInt,Libc.malloc(sizeof(T)*N + a - 1)), a))
end


text:
using SumOfSquares
using MosekTools
using DynamicPolynomials

function sos_lower_bound(factory, sparsity::Sparsity)
    @polyvar x y z
    K = @set 1-x^2-y^2>=0 && 1-y^2-z^2>=0
    p = 1 - x*y - y*z
    model = Model(factory)
    @variable(model, t)
    @objective(model, Max, t)
    @constraint(model, p - t in SOSCone(), domain = K, sparsity = sparsity)
    optimize!(model)
    return objective_value(model)
end

const factory = Mosek.Optimizer

# use chordal sparsity
sparse_value = sos_lower_bound(factory, VariableSparsity())

# no chordal sparsity
dense_value = sos_lower_bound(factory, NoSparsity())

println("Sparse and dense equivalent?")
println("Absolute different of objective value: ", abs(dense_value - sparse_value))


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Cbc
using CSV
using DataFrames

@testset "BenchmarkRunner" begin
    @info "Building training data..."
    instances = [build_knapsack_file_instance(), build_knapsack_file_instance()]
    stats = parallel_solve!(LearningSolver(Cbc.Optimizer), instances)
    @test length(stats) == 2
    @test stats[1] !== nothing
    @test stats[2] !== nothing

    benchmark = BenchmarkRunner(
        solvers = Dict(
            "baseline" => LearningSolver(Cbc.Optimizer, components = []),
            "ml-exact" => LearningSolver(Cbc.Optimizer),
            "ml-heur" => LearningSolver(Cbc.Optimizer, mode = "heuristic"),
        ),
    )
    @info "Fitting..."
    fit!(benchmark, instances)

    @info "Benchmarking..."
    parallel_solve!(benchmark, instances, n_trials = 2)

    csv_filename = tempname()
    write_csv!(benchmark, csv_filename)
    @test isfile(csv_filename)
    csv = DataFrame(CSV.File(csv_filename))
    @test size(csv)[1] == 12
end


text:
using QCMaterial
using Documenter

makedocs(;
    modules=[QCMaterial],
    authors="Rihito Sakurai <sakurairihito@gmail.com> and contributors",
    repo="https://github.com/sakurairihito/QCMaterialNew/blob/{commit}{path}#L{line}",
    sitename="QCMaterialNew",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/sakurairihito/QCMaterialNew",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Chapter1"=>[
            "Introduction" => "chapter1/intro.md",
            "About Goma-chan" => "chapter1/goma.md",
        ],
        "Chapter2"=>"chapter2/azarashi.md",
    ],
)

deploydocs(;
    devbranch="vqs",
    target="build",
    repo="github.com/sakurairihito/QCMaterialNew",
    versions=nothing
)

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
# This file was generated, do not modify it. # hide
scitype(X) <: input_scitype(booster)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@info "Caesar.jl is loading tools using ImageDraw.jl."

using .ImageDraw
using .Colors

export makeImage!

function makeImage!(pc::Caesar._PCL.PointCloud,
                    x_domain::Tuple{<:Real,<:Real}=(-1000,1000),
                    y_domain::Tuple{<:Real,<:Real}=x_domain;
                    pose=nothing,
                    ppose=nothing,
                    rows::Integer=1000,
                    cols::Integer=rows,
                    color::C=Gray(0.1),
                    trajCol=Gray(1.0),
                    img::AbstractMatrix{<:Colorant} = Gray.(zeros(rows,cols)),
                    circle_size::Real=1,
                    drawkws... ) where {C <: Colorant}
  #

  x_range = (x_domain[2]-x_domain[1])
  y_range = (y_domain[2]-y_domain[1])
  gridsize_x = x_range/rows
  gridsize_y = y_range/cols

  for pt in pc.points
    _x = trunc(Int, ( pt.x-x_domain[1])/gridsize_x ) #  + rows/2
    _y = trunc(Int, (-pt.y-y_domain[1])/gridsize_y ) #  + cols/2
    draw!( img, Ellipse(CirclePointRadius(_x, _y, circle_size)), color; drawkws... )  #; thickness = 1, fill = true)) )
  end

  # draw pose location
  if !(pose isa Nothing)
    _x = trunc(Int, ( pose[1]-x_domain[1])/gridsize_x ) #  + rows/2
    _y = trunc(Int, (-pose[2]-y_domain[1])/gridsize_y ) #  + cols/2

    draw!( img, Ellipse(CirclePointRadius(_x, _y, 4)), trajCol )

    if !(ppose isa Nothing)
      __x = trunc(Int, ( pose[1]-x_domain[1])/gridsize_x ) #  + rows/2
      __y = trunc(Int, (-pose[2]-y_domain[1])/gridsize_y ) #  + cols/2

      draw!( img, LineSegment(Point(_x, _y), Point(__x, __y)), trajCol )
    end
  end

  return img
end


text:
using WAV
using SpeechAugment: addEcho

T60 = 0.2
# read a wav file
wav,fs = wavread("xxx.wav")
# add reverberation effect on audio
echowav = addEcho(wav, fs, T60,
                  (3, 4, 2.5),      # room size
                  (1, 2, 1.2),      # sound src
                  (1.6, 2.6, 1.0)); # micphone

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
module LinSpaces2

export linspace2

import Base:
    getindex, unsafe_getindex, length,
    isempty, step,
    first, last, start, done, next

immutable LinSpace2{T<:AbstractFloat} <: Range{T}
    _start::T
    _stop::T
    len::T
    mult::T
end

function linspace2{T<:AbstractFloat}(start::T, stop::T, len::T)
    # Setup.  TODO: Need to copy some additional sanity checks from the real
    # LinSpace here.
    N1 = len-1
    mult = 1/N1
    # Compute values at ends according to formula used in getindex().
    # Unfortunately it's possible to have N1*mult == prevfloat(one(T))
    # in addition to the expected N1*mult == 1.  In IEEE round-to-nearest mode,
    # we can prove that these are the only possibilites and that optionally
    # modifying the stored end points by 1ulp is enough to bring the *computed*
    # end points back to the desired values.
    _start = N1*mult*start == start ? start : nextfloat(start)
    _stop  = N1*mult*stop  == stop  ? stop  : nextfloat(stop)

    LinSpace2(_start, _stop, len, mult)
end

function linspace2{T<:AbstractFloat}(start::T, stop::T, len::Real)
    T_len = convert(T, len)
    T_len == len || throw(InexactError())
    linspace2(start, stop, T_len)
end
linspace2(start::Real, stop::Real, len::Real=50) =
    linspace2(promote(AbstractFloat(start), AbstractFloat(stop))..., len)

@inline getindex{T}(r::LinSpace2{T}, i::Integer) = (checkbounds(r, i); unsafe_getindex(r, i))
@inline unsafe_getindex{T}(r::LinSpace2{T}, i::Integer) =
    (r.len-i)*r.mult*r._start + (i-1)*r.mult*r._stop

isempty{T}(r::LinSpace2{T}) = length(r) == 0

step{T}(r::LinSpace2{T}) = (r._stop - r._start)*r.mult

length{T}(r::LinSpace2{T}) = Integer(r.len + signbit(r.len - 1))

first{T}(r::LinSpace2{T}) = unsafe_getindex(r,1)
# FIXME: last() is wrong because it passes a T rather than an Integer, but for
# some mysterious reason this breakage makes the @simd benchmark much faster.
last{T}(r::LinSpace2{T}) = unsafe_getindex(r,r.len)

start{T}(r::LinSpace2{T}) = 1
done{T}(r::LinSpace2{T}, i::Int) = length(r) < i
@inline next{T}(r::LinSpace2{T}, i::Int) = (unsafe_getindex(r,i), i+1)

end # module