text:
function get_date_range(df::DataFrame, date_col_name::Symbol = :date)
    minimum(df[!, date_col_name]), maximum(df[!, date_col_name])
end

function filter_raw_data(df::DataFrame,
    test_fdate::D, test_tdate::D) where D<:Union{DateTime, ZonedDateTime}

    # filter dataframe by total date range and station Code
    df = @from i in df begin
        @where test_fdate <= i.date <= test_tdate
        @orderby i.date
        @select i
        @collect DataFrame
    end

    df
end

function filter_raw_data(df::DataFrame,
    train_fdate::D, train_tdate::D,
    test_fdate::D, test_tdate::D) where D<:Union{DateTime, ZonedDateTime}

    # filter dataframe by total date range and station Code
    df = @from i in df begin
        @where train_fdate <= i.date <= train_tdate ||
            test_fdate <= i.date <= test_tdate
        @orderby i.date
        @select i
        @collect DataFrame
    end

    df
end

"""
    filter_raw_data(df, stn_cdoe, fdate, tdate)

Used in analysis
"""

function filter_raw_data(df::DataFrame,
    stn_code::Integer,
    fdate::D, tdate::D) where D<:Union{DateTime, ZonedDateTime}

    # filter dataframe by total date range and station Code
    # @where fdate <= i.date <= tdate && i.stationCode == stn_code
    df = @from i in df begin
        @where i.stationCode == stn_code && fdate <= i.date <= tdate
        @orderby i.date
        @select i
        @collect DataFrame
    end

    df
end

function filter_station(df, stn_code)
    stn_df = @from i in df begin
        @where i.stationCode == stn_code
        @select i
        @collect DataFrame
    end

    stn_df
end

"""
    validate_dates(_from_date, _to_date, window_size, df)

Validate dates range
"""
function validate_dates(_from_date::D, _to_date::D,
    window_size::Integer, df::DataFrame) where {I<:Integer, D<:Union{DateTime, ZonedDateTime}}
    # Round date and simple date range validation
    from_date = ceil(max(_from_date, df.date[1]), Dates.Hour(1))
    to_date = floor(min(_to_date, df.date[end]), Dates.Hour(1))
    if from_date > to_date
        throw(ArgumentError("invalid date range: $from_date ~ $to_date"))
    end

    # sample_size + output_size should be avalable
    # .value : to get `Dates` value
    if Dates.Hour(to_date - from_date).value < window_size
        throw(BoundsError(
            "window size($window_size) is smaller than date range: $from_date ~ $to_date"))
    end

    from_date, to_date
end

"""
    window_df(df, sample_size, output_size, offset = 0)

Create list of dataframe with sample_size + output_size window with offset
"""
function window_df(df::DataFrame, sample_size::I, output_size::I,
    offset::I = zero(I)) where I <: Integer

    # start index for window
    # sample_size + hours (hours for Y , < sample_size) should be avalable
    start_idxs = collect(1:(size(df, 1) - sample_size - output_size + 1))
    final_idxs = start_idxs .+ (sample_size + output_size - 1)

    # create Dataframe array
    dfs = DataFrame[]

    # filter dataframe by total date range and station Code
    for (si, fi) in zip(start_idxs, final_idxs)
        push!(dfs, df[(si+offset):(fi+offset), :])
    end

    dfs
end

"""
    window_df(df, sample_size, output_size, _from_date, _to_date)

Create list of dataframe with sample_size + output_size window with date range
"""
function window_df(df::DataFrame, sample_size::I, output_size::I,
    _from_date::D, _to_date::D) where {I<:Integer, D<:Union{DateTime, ZonedDateTime}}

    window_size = sample_size + output_size
    from_date, to_date = validate_dates(_from_date, _to_date, window_size, df)

    # create Dataframe array
    dfs = DataFrame[]

    total_it = Dates.Hour((to_date - Dates.Hour(window_size - 1)) - from_date).value
    p = Progress(total_it, dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    @info "    Construct windows for current station..."

    for _wd_from_date = from_date:Dates.Hour(1):(to_date - Dates.Hour(window_size - 1))
        ProgressMeter.next!(p);
        _wd_to_date = _wd_from_date + Dates.Hour(window_size)
        _df = @from i in df begin
            @where _wd_from_date <= i.date < _wd_to_date
            @orderby i.date
            @select i
            @collect DataFrame
        end

        push!(dfs, _df)
    end

    dfs
end

"""
    window_df(df, sample_size, output_size, _from_date, _to_date, station_code)

Create list of dataframe with sample_size + output_size window with date range and station code

TODO: Too slow!
"""
function window_df(df::DataFrame, sample_size::I, output_size::I,
    _from_date::D, _to_date::D,
    station_code::I) where {I<:Integer, D<:Union{DateTime, ZonedDateTime}}

    window_size = sample_size + output_size
    from_date, to_date = validate_dates(_from_date, _to_date, window_size, df)

    # create Dataframe array
    dfs = DataFrame[]

    total_it = Dates.Hour((to_date - Dates.Hour(window_size - 1)) - from_date).value
    p = Progress(total_it, dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    @info "    Construct Training Set window for station $(station_code)..."

    for _wd_from_date = from_date:Dates.Hour(1):(to_date - Dates.Hour(window_size - 1))
        ProgressMeter.next!(p);
        _wd_to_date = _wd_from_date + Dates.Hour(window_size)
        _df = @from i in df begin
            @where _wd_from_date <= i.date < _wd_to_date && i.stationCode == station_code
            @orderby i.date
            @select i
            @collect DataFrame
        end

        push!(dfs, _df)
    end

    dfs
end

"""
    split_sizes3(total_size, batch_size)

Split `total_size`for train/valid/test set with size checking compared to `batch_size`

# Examples
```julia-repl
julia> split_sizes3(100, 10)
(64, 16, 20)

julia> split_sizes3(10, 20)
ERROR: AssertionError: train_size >= batch_size
Stacktrace:
 [1] split_sizes3(::Int64, ::Int64) at ./REPL[85]:5
 [2] top-level scope at none:0
````
"""
function split_sizes3(total_size::I, batch_size::I) where I <: Integer
    train_size, valid_size, test_size = split_sizes3(total_size)

    # at least contains single batch
    @assert train_size >= batch_size

    train_size, valid_size, test_size
end

"""
    split_sizes3(total_size)

Split `total_size` for train/valid/test set (0.64:0.16:0.20)

# Default ratios
(train + valid) : test = 0.8 : 0.2
train : valid = 0.8 : 0.2
train : valid : test = 0.64 : 0.16 : 0.20

# Examples
```julia-repl
julia> split_sizes3(100)
(64, 16, 20)
```
"""
function split_sizes3(total_size::I,
    train_ratio::AbstractFloat = 0.64,
    valid_ratio::AbstractFloat = 0.16) where I <: Integer

    train_size = round(total_size * train_ratio)
    valid_size = round(total_size * valid_ratio)
    test_size = total_size - (train_size + valid_size)

    @assert train_size + valid_size + test_size == total_size

    I(train_size), I(valid_size), I(test_size)
end

"""
    split_sizes2(total_size, batch_size)

Split `total_size` for train/valid set with size checking compared to `batch_size`

# Default ratios
train : valid = 0.8 : 0.2

# Examples
```julia-repl
julia> split_sizes2(100)
(80, 20)

julia> split_sizes2(10, 20)
ERROR: AssertionError: train_size >= batch_size
Stacktrace:
 [1] split_sizes2(::Int64, ::Int64) at ./REPL[85]:5
 [2] top-level scope at none:0
```
"""
function split_sizes2(total_size::I, batch_size::I) where I <: Integer

    train_size, valid_size = split_sizes2(total_size)

    # at least contains single batch
    @assert train_size >= batch_size

    train_size, valid_size
end

"""
    split_sizes2(total_size)

Split `total_size` for train/valid set (0.8:0.2)

# Default ratios
train : valid = 0.8 : 0.2

# Examples
```julia-repl
julia> split_sizes2(100)
(80, 20)
```
"""
function split_sizes2(total_size::I,
    train_ratio::AbstractFloat = 0.8) where I <: Integer

    train_size = round(total_size * train_ratio)
    valid_size = total_size - train_size

    @assert train_size + valid_size == total_size

    I(train_size), I(valid_size)
end

"""
    remove_sparse_input!(X, Y, ratio = 0.5)

Remove this pair if sparse ratio is greather than given `ratio`.
"""
function remove_sparse_input!(X, Y, missing_ratio=0.5)
    @assert 0.0 <= missing_ratio <= 1.0

    invalid_idxs = UnitRange[]
    size_y = length(Y)
    size_m = length(findall(_y -> ismissing(_y) || _y == zero(_y), Y))

    if (size_m / size_y) >= missing_ratio
        # set X and Y to zeros
        X = fill!(X, zero(X[1, 1]))
        Y = fill!(Y, zero(Y[1, 1]))
    end

    nothing
end

"""
    getX(df::DataFrame, features, sample_size)

Get 2D input X from DataFrame.

return size: (sample_size, length(features))
"""
function getX(df::DataFrame, features::Array{Symbol, 1}, sample_size::I) where I <: Integer
    convert(Matrix, df[1:sample_size, features])
end

"""
    getY(df, ycol, sample_size, output_size)

Get 1D output Y from DataFrame.

return size: (output_size,)
"""
function getY(df::DataFrame, ycol::Symbol, sample_size::I) where I <: Integer
    Array(df[(sample_size+1):end, ycol])
end

"""
    make_pair_DNN(df, ycol, features, sample_size, output_size)

Split single window to serialized X and Y
"""
function make_pair_DNN(df::DataFrame,
    ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, output_size::I, _eltype::DataType=Float32) where I <: Integer

    # get X (2D)
    _X = _eltype.(getX(df, features, sample_size))
    # get Y (1D)
    _Y = _eltype.(getY(df, ycol, sample_size))

    # for validation, sparse check not applied
    # GPU
    _X = vec(_X)
    _Y = _Y
    @assert length(_X) == sample_size * length(features)
    @assert length(_Y) == output_size
    @assert ndims(_X) == 1
    @assert ndims(_Y) == 1

    _X, _Y
end

"""
    make_pair_DNN(df, ycol, features, sample_size, output_size)

Split single window to serialized X and Y
"""
function make_pair_date_DNN(df::DataFrame,
    ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, output_size::I, _eltype::DataType=Float32) where I <: Integer

    # get X (2D)
    _X = _eltype.(getX(df, features, sample_size))
    # get Y (1D)
    _Y = _eltype.(getY(df, ycol, sample_size))
    _date = DateTime.(df[(sample_size+1):end, :date], Local)

    # for validation, sparse check not applied
    # GPU
    _X = vec(_X)
    _Y = _Y
    @assert length(_X) == sample_size * length(features)
    @assert length(_Y) == output_size
    @assert ndims(_X) == 1
    @assert ndims(_Y) == 1

    _X, _Y, _date
end

"""
    make_minibatch_DNN(df, ycol, chnks, features, sample_size, output_size, batch_size, Float64)

Create batch consists of pairs in `df` along with `idx` (row) and `features` (columns)
The batch is higher dimensional input consists of multiple pairs

single pairs to multiple
pairs = [(X_1, Y_1), (X_2, Y_2), ...]

[(input_single, output_single)...] =>
[
    ((input_single, input_single, ..., input_single,),  (output_single, output_single, ...output_single,)), // minibatch
    ((input_single, input_single, ..., input_single,),  (output_single, output_single, ...output_single,)), // minibatch
    ((input_single, input_single, ..., input_single,),  (output_single, output_single, ...output_single,)), // minibatch
]

and make pairs to be column stacked

# Example

batch_size = 5
pairs = [
            ([ 1, 2, 3, 4], [ 5, 6]),
            ([ 7, 8, 9,10], [11,12]),
            ([13,14,15,16], [17,18]),
            ([19,20,21,22], [23,24]),
            ([25,26,27,28], [29,30])]

become

([1,  7, 13, 19, 25;
  2,  8, 14, 20, 26;
  3,  9, 15, 21, 27;
  4, 10, 16, 22, 28],
 [5, 11, 17, 23, 29;
  6, 12, 18, 24, 30])

https://discourse.julialang.org/t/flux-support-for-mini-batches/14718/3
"""
function make_batch_DNN(dfs::Array{DataFrame, 1}, ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, output_size::I, batch_size::I,
    missing_ratio::AbstractFloat=0.5, _eltype::DataType=Float32) where I <: Integer

    # check Dataframe array size be `batch_size`
    @assert length(dfs) <= batch_size

    X = zeros(_eltype, sample_size * length(features), batch_size)
    Y = zeros(_eltype, output_size, batch_size)

    # input_pairs Array{Tuple{Array{Int64,1},Array{Int64,1}},1}
    for (i, df) in enumerate(dfs)
        # get X (2D)
        _X = _eltype.(vec(getX(df, features, sample_size)))
        # get Y (1D)
        _Y = _eltype.(getY(df, ycol, sample_size))

        # fill zeros if data is too sparse (too many misisngs)
        # only in training
        remove_sparse_input!(_X, _Y)

        # make X & Y as column stacked batch
        X[:, i] = _X
        Y[:, i] = _Y
    end

    X, Y
end

function make_batch_date_DNN(dfs::Array{DataFrame, 1}, ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, output_size::I, batch_size::I,
    missing_ratio::AbstractFloat=0.5, _eltype::DataType=Float32) where I <: Integer

    # check Dataframe array size be `batch_size`
    @assert length(dfs) <= batch_size

    X = zeros(_eltype, sample_size * length(features), batch_size)
    Y = zeros(_eltype, output_size, batch_size)
    _dates = Array{DateTime}(undef, output_size, batch_size)

    # input_pairs Array{Tuple{Array{Int64,1},Array{Int64,1}},1}
    for (i, df) in enumerate(dfs)
        # get X (2D)
        _X = _eltype.(vec(getX(df, features, sample_size)))
        # get Y (1D)
        _Y = _eltype.(getY(df, ycol, sample_size))
        _date = DateTime.(df[(sample_size+1):end, :date], Local)

        # fill zeros if data is too sparse (too many misisngs)
        # only in training
        remove_sparse_input!(_X, _Y)

        # make X & Y as column stacked batch
        X[:, i] = _X
        Y[:, i] = _Y
        _dates[:, i] = _date
    end

    X, Y, _dates
end

"""
    getX(df::DataFrame, idxs, features, sample_size)

Get 2D input X from DataFrame.

return size: (sample_size, length(features))
"""
function getX(df::DataFrame, idxs::UnitRange{I}, features::Array{Symbol, 1}) where I <: Integer
    convert(Matrix, df[idxs, features])
end

"""
    getY(df, ycol, sample_size, output_size)

Get 1D output Y from DataFrame.

return size: (output_size,)
"""
function getY(df::DataFrame, idxs::UnitRange{I}, ycol::Symbol) where I <: Integer
    Array(df[idxs, ycol])
end

"""
    make_pair_RNN(df, ycol, features,
    sample_size, kernel_length, output_size,
    _eltype=Float32)

Create batch CNN Input for LSTNet with batch_size == 1
Used in valid/test set
"""
function make_pair_RNN(df::DataFrame,
    ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, kernel_length::I, output_size::I,
    _eltype::DataType=Float32) where I <: Integer

    # O = (W - K + 2P) / S + 1
    # W = Input size, O = Output size
    # K = Kernel size, P = Padding size, S = Stride size
    pad_sample_size = kernel_length - 1

    # W : row in DataFrame
    # H : col in DataFrame
    #X_enc = zeros(_eltype, length(features), pad_sample_size + sample_size, 1, 1)
    X_enc = zeros(_eltype, pad_sample_size + sample_size, length(features), 1, 1)
    # decode array (Array of Array, batch sequences)
    X_dec = similar([zeros(_eltype, output_size, 1)], output_size)
    Y = zeros(_eltype, output_size, 1)

    # get X (2D)
    _X = _eltype.(getX(df, 1:sample_size, features))
    # get Y (1D)
    _Y = _eltype.(getY(df, (sample_size + 1):(sample_size + output_size), ycol))

    # WHCN order, Channel is 1 because this is not an image
    # left zero padding
    X_enc[(pad_sample_size + 1):end, :, 1, 1] = _X
    Y[:, 1] = _Y

    # feed decoder output
    # Ref : Yujin Tang, et. al, Sequence-to-Sequence Model with Attention for Time Series Classification
    for i in 1:output_size
        # first row is zero, but from second row, it feeds output from previous time step
        # this makes to assert sample_size >= output_size
        #_Y = _eltype.(getY(df, (i + 1):(i + output_size - 1), ycol))
        #_Y = view(df, (i + 1):(i + output_size - 1), ycol)
        _X_dec = zeros(_eltype, output_size, 1)
        _X_dec[2:output_size] = view(df, (i + 1):(i + output_size - 1), ycol)
        X_dec[i] = _X_dec
    end

    X_enc, X_dec, Y
end

function make_pair_date_RNN(df::DataFrame,
    ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, kernel_length::I, output_size::I,
    _eltype::DataType=Float32) where I <: Integer

    # O = (W - K + 2P) / S + 1
    # W = Input size, O = Output size
    # K = Kernel size, P = Padding size, S = Stride size
    pad_sample_size = kernel_length - 1
    #X_enc = zeros(_eltype, length(features), pad_sample_size + sample_size, 1, 1)
    X_enc = zeros(_eltype, pad_sample_size + sample_size, length(features), 1, 1)
    # decode array (Array of Array, batch sequences)
    _x = zeros(_eltype, output_size, 1)
    X_dec = similar([zeros(_eltype, output_size, 1)], output_size)
    Y = zeros(_eltype, output_size, 1)

    # get X (2D)
     _X = _eltype.(getX(df, 1:sample_size, features))
    # get Y (1D)
     _Y = _eltype.(getY(df, (sample_size + 1):(sample_size + output_size), ycol))
    _date = DateTime.(df[(sample_size+1):end, :date], Local)

    # WHCN order, Channel is 1 because this is not an image
    # left zero padding
    X_enc[(pad_sample_size + 1):end, :, 1, 1] = _X
    Y[:, 1] = _Y

    # feed decoder output
    # Ref : Yujin Tang, et. al, Sequence-to-Sequence Model with Attention for Time Series Classification
    for i in 1:output_size
        # first row is zero, but from second row, it feeds output from previous time step
        # this makes to assert sample_size >= output_size
        X_dec[i] = vcat(zeros(_eltype, 1, 1), _eltype.(getY(df, (i + 1):(i + output_size - 1), ycol)))
    end

    X_enc, X_dec, Y, _date
end

"""
    make_batch_RNN(dfs, ycol, features,
    sample_size, kernel_length, output_size, batch_size,
    _eltype=Float32)

Create batch CNN Input for LSTNet
Used in train set
"""
function make_batch_RNN(dfs::Array{DataFrame, 1},
    ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, kernel_length::I, output_size::I, batch_size::I,
    _eltype::DataType=Float32) where I <: Integer

    # O = (W - K + 2P) / S + 1
    # W = Input size, O = Output size
    # K = Kernel size, P = Padding size, S = Stride size
    pad_sample_size = kernel_length - 1
    #X_enc = zeros(_eltype, length(features), pad_sample_size + sample_size, 1, batch_size)
    X_enc = zeros(_eltype, pad_sample_size + sample_size, length(features), 1, batch_size)
    # decode array (Array of Array, batch sequences)
    X_dec = similar([zeros(_eltype, output_size, 1)], output_size)
    Y = zeros(_eltype, output_size, batch_size)

    # feed decoder output
    # zero padding on input matrix
    for (i, df) in enumerate(dfs)
        # get X (2D)
        _X = _eltype.(getX(df, 1:sample_size, features))
        #_X = view(df, 1:sample_size, features)
        # get Y (1D)
        _Y = _eltype.(getY(df, (sample_size + 1):(sample_size + output_size), ycol))
        #_Y = view(df, (sample_size + 1):(sample_size + output_size), ycol)

        # WHCN order, Channel is 1 because this is not an image
        # left zero padding
        #X_enc[:, (pad_sample_size + 1):end, 1, i] = transpose(_X)
        #X_enc[:, (pad_sample_size + 1):end, 1, i] = transpose(collect(_X))
        X_enc[(pad_sample_size + 1):end, :, 1, i] = _X
        Y[:, i] = _Y
    end

    # Ref : Yujin Tang, et. al, Sequence-to-Sequence Model with Attention for Time Series Classification
    for i in 1:output_size
        _X_dec = zeros(_eltype, output_size, batch_size)
        for (j, df) in enumerate(dfs)
            # first row is zero, but from second row, it feeds output from previous time step
            # this makes to assert sample_size >= output_size
            #_Y = _eltype.(getY(df, (i + 1):(i + output_size - 1), ycol))

            _X_dec[:, j] = vcat(zeros(_eltype, 1, 1), _eltype.(getY(df, (i + 1):(i + output_size - 1), ycol)))
        end
        X_dec[i] = _X_dec
    end

    X_enc, X_dec, Y
end

function make_batch_date_RNN(dfs::Array{DataFrame, 1},
    ycol::Symbol, features::Array{Symbol, 1},
    sample_size::I, kernel_length::I, output_size::I, batch_size::I,
    _eltype::DataType=Float32) where I <: Integer

    # O = (W - K + 2P) / S + 1
    # W = Input size, O = Output size
    # K = Kernel size, P = Padding size, S = Stride size
    pad_sample_size = kernel_length - 1
    #X_enc = zeros(_eltype, length(features), pad_sample_size + sample_size, 1, batch_size)
    X_enc = zeros(_eltype, pad_sample_size + sample_size, length(features), 1, batch_size)
    # decode array (Array of Array, batch sequences)
    X_dec = similar([zeros(_eltype, output_size, 1)], output_size)
    Y = zeros(_eltype, output_size, batch_size)
    _dates = Array{DateTime}(undef, output_size, batch_size)

    # feed decoder output
    # zero padding on input matrix
    for (i, df) in enumerate(dfs)
        # get X (2D)
        _X = _eltype.(getX(df, 1:sample_size, features))
        #_X = view(df, 1:sample_size, features)
        # get Y (1D)
        _Y = _eltype.(getY(df, (sample_size + 1):(sample_size + output_size), ycol))
        #_Y = view(df, (sample_size + 1):(sample_size + output_size), ycol)
        #_date = DateTime.(df[(sample_size+1):end, :date], Local)

        # WHCN order, Channel is 1 because this is not an image
        # left zero padding
        #X_enc[:, (pad_sample_size + 1):end, 1, i] = transpose(_X)
        # X_enc[:, (pad_sample_size + 1):end, 1, i] = transpose(collect(_X))
        X_enc[(pad_sample_size + 1):end, :, 1, i] = _X
        Y[:, i] = _Y
        _dates[:, i] = DateTime.(df[(sample_size+1):end, :date], Local)
    end

    # Ref : Yujin Tang, et. al, Sequence-to-Sequence Model with Attention for Time Series Classification
    for i in 1:output_size
        _X_dec = zeros(_eltype, output_size, batch_size)
        for (j, df) in enumerate(dfs)
            # first row is zero, but from second row, it feeds output from previous time step
            # this makes to assert sample_size >= output_size
            #_Y = _eltype.(getY(df, (i + 1):(i + output_size - 1), ycol))

            #_X_dec[2:output_size, j] = _Y
            _X_dec[:, j] = vcat(zeros(_eltype, 1, 1), _eltype.(getY(df, (i + 1):(i + output_size - 1), ycol)))
        end
        X_dec[i] = _X_dec
    end

    X_enc, X_dec, Y, _dates
end

"""
    serializeBatch(X, Y)

Reconstruct batch for RNN by serialize batch array.
Output will be arryas of 1D
"""
function serializeBatch(X, Y; dims=2)
    # after CNN, X and Y should have same size
    @assert size(X, 2) == size(Y, 2)

    collect(zip(eachslice(X, dims=dims), eachslice(Y, dims=dims)))
end

"""
    serializeBatch(X)

Convert 2D batch array to array of 1D array
(m x n) array -> n-of m x 1 arrays

i.e.
X = (m x n) -> [(m x 1), ...]
"""
serializeBatch(A; dims=2) = collect(eachslice(A, dims=dims))

"""
    whcn2cnh(X)

Convert 4D WHCN array to CNH

(sample_size, 1, hidCNN, batch_size)
-> (hidCNN, batch_size, sample_size)
"""
function whcn2cnh(x::AbstractArray{N, 4}) where {N<:Number}
    @assert size(x, 2) == 1

    x = reshape(x, size(x, 1), size(x, 3), size(x, 4))

    permutedims(x, [2, 3, 1])
end

"""
    matrix2arrays(X)

2D matrix to array of arrays
"""
function matrix2arrays(x::AbstractArray{N, 2}) where {N<:Number}
    [x[i, :] |> gpu for i in axes(x, 1)]
end

"""
    is_sparse_Y(Y, ratio = 0.5)

Check Y is sparse, 0s are more than `raito`
"""
function is_sparse_Y(Y, missing_ratio=0.5)

    @assert 0.0 <= missing_ratio <= 1.0

    size_y = length(Y)
    size_m = length(findall(_y -> ismissing(_y) || _y == 0.0, Y))

    if (size_m / size_y) >= missing_ratio
        # if use deleteat here, it manipulates pairs while loop and return wrong results
        return true
    end

    false
end

zero2Missing!(df::DataFrame, ycol::Symbol) = replace!(df[!, ycol], 0 => missing, 0.0 => missing, NaN => missing)

"""
    impute!(df, ycol, method; total_mean)

# References
## BCPA
    * Shigeyuki Oba, et al. A BPCA Based Missing Value Imputing Method for Traffic Flow Volume Data
    * Li Qu, et. al. A Bayesian missing value estimation method for gene expression profile data
"""
#=
function impute!(df::DataFrame, ycol::Symbol, method::Symbol; total_mean::AbstractFloat = 0.0, ycols::Array{Symbol, 1} = [ycol], leafsize::Integer = 10)
    # filter missing value
    nomissings = filter(x -> x !== missing, df[!, ycol])

    if method == :mean
        replace!(df[!, ycol], missing => Int(round(total_mean)))
    elseif method == :sample
        _weeks = 53
        arr = collect(skipmissing(float.(df[!, ycol])))
        dfs = []

        for _week in 1:_weeks
            push!(dfs,
                DataFrames.filter(row -> week(row[:date]) == _week && row[ycol] !== missing && row[ycol] !== NaN, df))
        end

        weeks = week.(df[!, :date])
        for (i, _val) in enumerate(df[!, ycol])
            if _val === missing
                _df = dfs[weeks[i]]
                df[i, ycol] = sample(_df[!, ycol])
            end
        end
    #elseif method == :knn
    elseif method == :bpca
        # TODO: variational Bayes (VB) algorithm
        # bpca_init!(df, ycol, total_mean)
    end
end
=#

"""
    padded_push!(df, array, name, train_fdate, train_tdate)

append arrays to dataframe padded by train data

This use simple loop, so it may be slow
"""
padded_push!(df::DataFrame, array::AbstractArray, name::Symbol, train_fdate::ZonedDateTime, train_tdate::ZonedDateTime) =
    padded_push!(df, array, name, DateTime(train_fdate, Local), DateTime(train_tdate, Local))

function padded_push!(df::DataFrame, array::AbstractArray, name::Symbol, train_fdate::DateTime, train_tdate::DateTime)
    # get array filled with zero instead of missing
    # If I fill with missing, zscore fails due to msising value
    # for test_date, I don't need that.
    padded = zeros(size(df, 1))

    dates = train_fdate:Dates.Hour(1):train_tdate

    # should match
    @assert size(array) == size(dates)

    # now use dates as index by key-value structure
    train_values = Dict(dates .=> array)

    # fill value from date
    for (i, row) in enumerate(eachrow(df))
        date = DateTime(row[:date], Local)

        if haskey(train_values, date)
            padded[i] = train_values[date]
        end
    end

    DataFrames.insertcols!(df, 1, name => padded)
end

"""
    construct_annual_table(sea_year, sea_day, train_fdate, train_tdate)

Construct table of seasonality data within train date range
table accepts month, day, hour as keys.
"""
construct_annual_table(sea_year::AbstractVector, sea_day::AbstractVector, train_fdate::ZonedDateTime, train_tdate::ZonedDateTime) =
    construct_annual_table(sea_year, sea_day, DateTime(train_fdate, Local), DateTime(train_tdate, Local))

function construct_annual_table(df::DataFrame, target::Symbol, train_fdate::DateTime, train_tdate::DateTime)
    months = []
    days = []
    hours = []
    sea_years = []
    sea_days = []

    # ndsparse is immutable object, so make it dataframe first then group it.
    # 1. filter dataframe by date
    _df = DataFrames.filter(row -> train_fdate < DateTime(row[:date], Local) < train_tdate, df)

    # 2. Then make month, day, hour columns
    DataFrames.insertcols!(_df, 1, :month => Dates.month.(_df[!, :date]))
    DataFrames.insertcols!(_df, 1, :day => Dates.day.(_df[!, :date]))
    DataFrames.insertcols!(_df, 1, :hour => Dates.hour.(_df[!, :date]))

    # 2. group df using keys (month, day, hour)
    gd = DataFrames.groupby(_df, [:month, :day, :hour])

    # unique
    u_months = []
    u_days = []
    u_hours = []
    u_sea_years_fit = []
    u_sea_years_org = []
    u_sea_days = []

    for gdf in gd
        ugdf = unique(gdf)

        # length of ugdf is same as # of stations
        push!(u_months, ugdf[!, :month][1])
        push!(u_days, ugdf[!, :day][1])
        push!(u_hours, ugdf[!, :hour][1])
        push!(u_sea_years_fit, ugdf[!, Symbol(target, "_year_fit")][1])
        push!(u_sea_years_org, ugdf[!, Symbol(target, "_year_org")][1])
        push!(u_sea_days, ugdf[!, Symbol(target, "_day")][1])
    end

    ndsparse((month = u_months, day = u_days, hour = u_hours),
    (season_year_fit = u_sea_years_fit, season_year_org = u_sea_years_org,
    season_day = u_sea_days,))
end

"""
    padded_push!(df, array, name, train_fdate, train_tdate)

Append arrays to dataframe padded by test data
Assume that train data have been already appended,
so just replace missing to real value by given seaonality table

This use simple loop, so it may be slow
"""
padded_push!(df::DataFrame, target::Symbol,
    year_sea_col::Symbol, day_sea_col::Symbol, day_res_col::Symbol,
    season_table::AbstractNDSparse,
    test_fdate::ZonedDateTime, test_tdate::ZonedDateTime) =
    padded_push!(df, target, year_sea_col, day_sea_col, day_res_col, season_table,
    DateTime(test_fdate, Local), DateTime(test_tdate, Local))

function padded_push!(df::DataFrame, target::Symbol,
    year_sea_col::Symbol, day_sea_col::Symbol, day_res_col::Symbol,
    season_table::AbstractNDSparse,
    test_fdate::DateTime, test_tdate::DateTime)

    # fill value from table and compute residual
    for (i, row) in enumerate(eachrow(df))
        # check date in range
        if test_fdate <= DateTime(row.date, Local) <= test_tdate
            df[i, year_sea_col] = season_table[Dates.month(row.date), Dates.day(row.date), Dates.hour(row.date)].season_year_fit
            df[i, day_sea_col] = season_table[Dates.month(row.date), Dates.day(row.date), Dates.hour(row.date)].season_day
            df[i, day_res_col] = df[i, target] - df[i, year_sea_col] - df[i, day_sea_col]
        end
    end
end

padded_push!(df::DataFrame, target::Symbol,
    year_sea_fit_col::Symbol, year_sea_org_col::Symbol, day_sea_col::Symbol, day_res_col::Symbol,
    season_table::AbstractNDSparse,
    test_fdate::ZonedDateTime, test_tdate::ZonedDateTime) =
    padded_push!(df, target, year_sea_fit_col, year_sea_org_col, day_sea_col, day_res_col, season_table,
    DateTime(test_fdate, Local), DateTime(test_tdate, Local))

function padded_push!(df::DataFrame, target::Symbol,
    year_sea_fit_col::Symbol, year_sea_org_col::Symbol, day_sea_col::Symbol, day_res_col::Symbol,
    season_table::AbstractNDSparse,
    test_fdate::DateTime, test_tdate::DateTime)

    # fill value from table and compute residual
    for (i, row) in enumerate(eachrow(df))
        # check date in range
        if test_fdate <= DateTime(row.date, Local) <= test_tdate
            df[i, year_sea_fit_col] = season_table[Dates.month(row.date), Dates.day(row.date), Dates.hour(row.date)].season_year_fit
            df[i, year_sea_org_col] = season_table[Dates.month(row.date), Dates.day(row.date), Dates.hour(row.date)].season_year_org
            df[i, day_sea_col] = season_table[Dates.month(row.date), Dates.day(row.date), Dates.hour(row.date)].season_day
            df[i, day_res_col] = df[i, target] - df[i, year_sea_fit_col] - df[i, day_sea_col]
        end
    end
end

function compose_seasonality(dates::Array{DateTime, 1}, result, season_table::AbstractNDSparse, fitted=true)
    org_values = Array{Float64}(undef, size(dates, 1))

    if fitted
        for (i, date) in enumerate(dates)
            org_values[i] = result[i] +
                season_table[Dates.month(date), Dates.day(date), Dates.hour(date)].season_year_fit +
                season_table[Dates.month(date), Dates.day(date), Dates.hour(date)].season_day
        end
    else
        for (i, date) in enumerate(dates)
            org_values[i] = result[i] +
                season_table[Dates.month(date), Dates.day(date), Dates.hour(date)].season_year_org +
                season_table[Dates.month(date), Dates.day(date), Dates.hour(date)].season_day
        end
    end

    org_values
end
#=
function bpca_fill(df::DataFrame, ycol_t::Symbol, ycols::Array{Symbol, 1})
    # init values
    # size of given array (N x D)
    N = size(df[!, ycol], 2)
    # Univarate
    D = 1
    # default number of PCA axis
    Q = D - 1 > 0 ? D - 1 : D

    y_est_mean = df[!, ycol_t]
    y_est_zero = df[!, ycol_t]

    # for SVD,
    replace!(y_est_zero, missing => 0)
    # fill by mean values
    replace!(y_est_mean, missing => Int(round(total_mean)))

    mean_arr = map(ycols) do ycol
        mean(skipmissing(float.(df[!, ycol])))
    end

    # covariance
    covy = Statistics.cov(y_est_zero)
    # SVD
    F = LinearAlgebra.cov(cov)

    # PCA
    W = F.U * sqrt.(S)

    # variance
    tau = 1.0 ./( sum(Diagonal(covy)) - sum(Diagonal(F.S)) )

    taumax = 1e10;
    taumin = 1e-10;
    tau = max( min( tau, taumax), taumin );

    galpha0 = 1e-10;
    balpha0 = 1;
    alpha = ( 2 * galpha0 + d) ./ (tau * Diagonal(transpose(W)*W) + 2 * galpha0 / balpha0);

    gmu0  = 0.001;

    btau0 = 1;
    gtau0 = 1e-10;
    SigW = Matrix{Float64}(I, Q, Q);

    epochs = 200

    for i = 1:epochs

        trS = _bpca_no_miss()
    end
end

function _bpca_no_miss(Q::I), tau, W) where I<:Integer
    Iq = Matrix{Float64}(I, Q, Q)
    # E-step (Expectation)
    Rx = Iq .+ tau * transpose(W) * W + SigW;
    Rxinv = inv( Rx );
    dy = df[!, ycol] .- size(df[!, ycol], 2);
    x = tau * Rxinv * transpose(W) * transpose(dy);

    T = transpose(dy) * transpose(x);

    # trS
    sum( sum( dy .* dy ));

    # M-step (Maximization)
end
=#


text:
using StochasticPrograms
using Random

##

default_es_pars = Dict((
    :c_i => .3,
    :c_o => .05,
    :c_sto_op => 0.0001,
    :asset_lifetime => 10.,
    :c_pv => 700.,
    :c_wind => 2000.,
    :c_storage => 600.,
    :c_heat_storage => 30.,
    :c_heatpump => 20.,
    :inv_budget => 500000000.,
    :recovery_time => 72,
    :COP => 0.8,
    :heat_losses => 0.01
))

function simple_flex_sampler(n, F_max, t_max)
    [@scenario t_xi = rand(1:t_max) s_xi = rand([-1, 1]) F_xi = rand() * F_max probability = 1/n
        for i in 1:n]
end

function no_flex_pseudo_sampler()
    [@scenario t_xi = 1 s_xi = 1 F_xi = 0. probability = 1.
    ,]
end

function define_energy_system(pv, wind, demand, heatdemand; p = default_es_pars, strict_flex=false)
    number_of_hours = minimum([length(pv), length(demand), length(wind)])
    c_sto_op = p[:c_sto_op]
    c_i = p[:c_i]
    c_o = p[:c_o]
    recovery_time = p[:recovery_time]
    COP = p[:COP]
    heat_losses = p[:heat_losses]
    energy_system = @stochastic_model begin
        @stage 1 begin
            @parameters begin
                # Expected lifetime of components, years
                asset_lifetime = p[:asset_lifetime]
                # Costs in Euro/1kWp
                c_pv = p[:c_pv]
                c_wind = p[:c_wind]
                c_storage = p[:c_storage]
                c_sto_op = c_sto_op
                c_i = c_i
                c_o = c_o
                c_heat_storage = p[:c_heat_storage]
                c_heatpump = p[:c_heatpump]
                COP = COP
                heat_losses = heat_losses
                # Euro
                inv_budget = p[:inv_budget] # Make the problem bounded
            end
            lifetime_factor = asset_lifetime * 365 * 24 / number_of_hours
            # Component units to be invested in, kWp
            @decision(model, u_pv >= 0)
            @decision(model, u_wind >= 0)
            @decision(model, u_storage >= 0)
            @constraint(model, c_pv * u_pv + c_wind * u_wind + c_storage * u_storage <= inv_budget)
            # Grid connection
            @decision(model, gci[t in 1:number_of_hours] >= 0)
            @decision(model, gco[t in 1:number_of_hours] >= 0)
            # Storage model
            @decision(model, sto_in[t in 1:number_of_hours] >= 0) # into the bus from storage
            @decision(model, sto_out[t in 1:number_of_hours] >= 0)
            @decision(model, sto_soc[t in 1:number_of_hours] >= 0)
            @constraint(model, [t in 1:number_of_hours-1], sto_soc[t+1] == sto_soc[t] - sto_out[t] + sto_in[t])
            @constraint(model, [t in 1:number_of_hours], sto_soc[t] <= u_storage)
            @constraint(model, sto_soc[1] == u_storage / 2)
            @constraint(model, sto_soc[number_of_hours] - sto_out[number_of_hours] + sto_in[number_of_hours] == sto_soc[1])
            # Heat model
            @decision(model, u_heatpump >= 0)
            @decision(model, u_heat_storage >= 0)
            @decision(model, heat_sto_in[t in 1:number_of_hours] >= 0) # from the heat storage
            @decision(model, heat_sto_out[t in 1:number_of_hours] >= 0)
            @decision(model, heat_sto_soc[t in 1:number_of_hours] >= 0)
            @decision(model, heatpumpflow[t in 1:number_of_hours] >= 0)

            @constraint(model, [t in 1:number_of_hours-1], heat_sto_soc[t+1] == heat_sto_soc[t] - heat_sto_out[t] + heat_sto_in[t])
            @constraint(model, [t in 1:number_of_hours], heat_sto_soc[t] <= u_heat_storage)
            @constraint(model, heat_sto_soc[1] == u_heat_storage / 2)
            @constraint(model, heat_sto_soc[number_of_hours] - heat_sto_out[number_of_hours] + heat_sto_in[number_of_hours] == heat_sto_soc[1])
            # Energy balance
            @constraint(model, [t in 1:number_of_hours],
            gci[t] - gco[t] + u_pv * pv[t] + u_wind * wind[t] - demand[t] + sto_in[t] - sto_out[t] - heatpumpflow[t] == 0)
            # Heat balance
            @constraint(model, [t in 1:number_of_hours], -heatdemand[t] - heat_sto_in[t] + heat_sto_out[t] + COP*heatpumpflow[t] - heat_losses*heat_sto_soc[t] == 0)
            # Investment costs
            @objective(model, Min, (u_pv * c_pv + u_wind * c_wind + u_storage * c_storage
            + u_heat_storage * c_heat_storage + u_heatpump * c_heatpump) / lifetime_factor
            + c_i * sum(gci) - c_o * sum(gco) +
            c_sto_op * sum(sto_in) + c_sto_op * sum(sto_out))
        end
        @stage 2 begin
            @parameters begin
                recovery_time = recovery_time
                c_sto_op = c_sto_op
                c_i = c_i
                c_o = c_o
                penalty = 10000
                heat_losses = heat_losses
                COP = COP
            end
            @uncertain t_xi s_xi F_xi # t_xi the time of flexibility demand, s_xi - sign (±1 or 0)
            t_xi_final = t_xi + recovery_time - 1
            @known(model, u_pv)
            @known(model, u_wind)
            @known(model, u_storage)
            @known(model, gci)
            @known(model, gco)
            @known(model, sto_in)
            @known(model, sto_out)
            # Post event components
            # Grid connection
            @recourse(model, gci2[t in 1:recovery_time] >= 0)
            @recourse(model, gco2[t in 1:recovery_time] >= 0)
            # Storage model
            @recourse(model, sto_in2[t in 1:recovery_time] >= 0) # into the bus from storage
            @recourse(model, sto_out2[t in 1:recovery_time] >= 0)
            @recourse(model, sto_soc2[t in 1:recovery_time] >= 0)
            @constraint(model, [t in 1:recovery_time-1], sto_soc2[t+1] == sto_soc2[t] - sto_out2[t] + sto_in2[t])
            @constraint(model, [t in 1:recovery_time], sto_soc2[t] <= u_storage)
            @constraint(model, sto_soc2[1] == sto_soc[t_xi])
            @constraint(model, sto_soc2[recovery_time] - sto_out2[recovery_time] + sto_in2[recovery_time] == sto_soc[t_xi+recovery_time])
            # Heat model
            @recourse(model, heat_sto_in2[t in 1:recovery_time] >= 0) # from the heat storage
            @recourse(model, heat_sto_out2[t in 1:recovery_time] >= 0)
            @recourse(model, heat_sto_soc2[t in 1:recovery_time] >= 0)
            @recourse(model, heatpumpflow2[t in 1:recovery_time] >= 0)
            @constraint(model, [t in 1:recovery_time-1], heat_sto_soc2[t+1] == heat_sto_soc2[t] - heat_sto_out2[t] + heat_sto_in2[t])
            @constraint(model, [t in 1:recovery_time], heat_sto_soc2[t] <= u_heat_storage)
            @constraint(model, heat_sto_soc2[1] == heat_sto_soc[t_xi])
            @constraint(model, heat_sto_soc2[recovery_time] - heat_sto_out2[recovery_time] + heat_sto_in2[recovery_time] == heat_sto_soc[t_xi+recovery_time])

            # Event energy balance
            # The storage and other fast acting components use the recourse variables here.
            # They provide the balance. Grid connection is not allowed, as we are suporting the grid here.
            @constraint(model, gci2[1] - gco2[1] + u_pv * pv[t_xi] + u_wind * wind[t_xi]
             - demand[t_xi] + sto_in2[1] - sto_out2[1]
             - heatpumpflow2[1] + F_xi * s_xi == 0)
            if strict_flex
                @constraint(model, gci2[1] == gci[t_xi])
                @constraint(model, gco2[1] == gco[t_xi])
            end
            ## Utility variables to linearize min|gci[t_xi]-gci2[1]|
            @recourse(model, gi1>=0)
            @recourse(model, gi2>=0)
            @constraint(model, gci2[1]-gci[t_xi] == gi1-gi2)
            @recourse(model, go1>=0)
            @recourse(model, go2>=0)
            @constraint(model, gco2[1]-gco[t_xi] == go1-go2)
            # Post event energy balance
            @constraint(model, [t in 2:recovery_time],
            gci2[t] - gco2[t]
            + u_pv * pv[t + t_xi - 1] + u_wind * wind[t + t_xi - 1]
            - demand[t + t_xi - 1] + sto_in2[t] - sto_out2[t] - heatpumpflow2[t] == 0)
            # Heat balance
            @constraint(model, [t in 1:recovery_time], -heatdemand[t] - heat_sto_in2[t] + heat_sto_out2[t] + COP*heatpumpflow2[t] - heat_losses*heat_sto_soc2[t] == 0)

            @objective(model, Min,
            + c_i * (sum(gci2) - sum(gci[t_xi:t_xi_final]))
            - c_o * (sum(gco2) - sum(gco[t_xi:t_xi_final]))
            + penalty * (gi1 + gi2) + penalty * (go1 + go2)
            + c_sto_op * (sum(sto_in2) + sum(sto_out2)
            - sum(sto_in[t_xi:t_xi_final]) - sum(sto_out[t_xi:t_xi_final])))
        end
    end
    energy_system
end

function fix_investment!(sp, vars, values)
    for (var, val) in zip(vars, values)
        fix(decision_by_name(sp, 1,var), val)
    end
end

function unfix_investment!(sp, vars, )
    for var in vars
        unfix(decision_by_name(sp, 1,var))
    end
end

ext: jl


text:
struct Sample{T, I<:Integer}
  size::Size{I}
  patches::Vector{Patch{T,I}}
  background::T
  set::Set{T}
end

Sample(size::Size{I}, background::T, set::Set{T}) where {I<:Integer, T} =
  Sample(size, Vector{Patch{T,I}}(), background, set)

Sample(size::Size{I}, background::T, set::Vector{T}) where {I<:Integer, T} =
  Sample(size, Vector{Patch{T,I}}(), background, Set(set))

patches(s::Sample) = s.patches
Base.size(s::Sample) = (height(s.size), width(s.size))
Base.size(s::Sample, i::Integer) = i == 1 ? height(s.size) : width(s.size)

"Fill a sample with random patches"
function Base.rand(s::Sample{T,I}; patchsize::Size = Size(8), nbpatches::Integer = 5,
  nbsame::Integer = 3, diff::Integer = 1) where {T, I<:Integer}
  res = deepcopy(s)
  # Gather non-overlapping positions for patches
  patchpos = nonoverlappingrandrect(width(res.size), height(res.size),
                width(patchsize), height(patchsize), nbpatches)
  # Generate actual patches from the positions
  pos = pop!(patchpos)
  patch = rand(Patch, patchsize, s.set, top(pos), left(pos))
  push!(patches(res), patch)
  # Generate same patches
  for i in 1:nbsame-1
    pos = pop!(patchpos)
    newpatch = Patch(copy(content(patch)), set(patch), top(pos), left(pos))
    push!(patches(res), newpatch)
  end
  # Generate altered additional patches
  while(!isempty(patchpos))
    pos = pop!(patchpos)
    newpatch = Patch(copy(content(patch)), set(patch), top(pos), left(pos))
    flip!(newpatch, diff)
    push!(patches(res), newpatch)
  end
  return res
end

"Return an Array representing the sample"
function render(s::Sample)
  res = fill(s.background, (s.size.height, s.size.width))
  for p in s.patches
    paste!(res, p.content, p.top, p.left)
  end
  return res
end


text:
module BenchMisc
using BenchmarkTools
using Statistics
using Gadfly
using MLStyle
using DataFrames
import Match
import Rematch
using ..ArbitrarySampler
using ..Utils

abstract type AbstractUserType end
struct UserTy1{A} <: AbstractUserType
    field1 :: Int
    field2 :: String
    field3 :: A
end
@as_record UserTy1
struct UserTy2 <: AbstractUserType end
@as_record UserTy2

mod′(n) = x -> mod(x, n)

any_user1 = @spec UserTy1(::Int, ::String, _)
any_user  = @spec UserTy2() || $any_user1
specs = [
    :s1 => @spec((1 || 2, _, ::String || ::Symbol)),
    :s2 => @spec([2 || 0, ::Integer, ::Integer, 3 || 0]),
    :s3 => @spec(UserTy1(::Int isa mod′(4), ::String, (7||8 , 8||9))),
    :s4 => @spec("yes" || "no" || "not sure" || _),
    :s5 => @spec(UserTy1(::Int isa mod′(5), "no" || "yes", $any_user)),
    :s6 => @spec([7 || _, (1 || _, $any_user), (_{2})..., ::Int || _]),
    :s7 => @spec((_::String, ("1" || "2" || "5", UserTy2(), ), ::Int isa (x -> 9 + mod(x, 5)))),
    :s8 => @spec((_, 10, _, 20, _)),
    :_  => @spec(_)
]

implementations = [
    :MLStyle => (@λ begin
        (1, _, ::String)              -> 1
        ([2, _, _, 3] && [_, a, b, _] && if a > b end) -> 2
        UserTy1(2, ::String, (7, 9))  -> 3
        "yes"                         -> 4
        UserTy1(2, "no", ::AbstractUserType)      -> 5
        (2, 3, ::Real, 5)                         -> 6
        [7, (1, ::AbstractUserType), _..., ::Int] -> 7
        (::String, ("5", ::UserTy2, ), 10) -> 8
        _                                  -> 9
    end),
    :Rematch => function (x)
        Rematch.@match x begin
            (1, _, _::String)              => 1
            ([2, a, b, 3] where a > b)  => 2
            UserTy1(2, _::String, (7, 9))  => 3
            "yes"                         => 4
            UserTy1(2, "no", _::AbstractUserType)      => 5
            (2, 3, _::Real, 5)                         => 6
            [7, (1, _::AbstractUserType), _..., _::Int] => 7
            (_::String, ("5", _::UserTy2, ), 10) => 8
            _                                  => 9
        end
    end,
    Symbol("Match.jl") => function (x)
        Match.@match x begin
            (1, _, _::String)              => 1
            ([2, a, b, 3], if a > b end)  => 2
            UserTy1(2, _::String, (7, 9))  => 3
            "yes"                         => 4
            UserTy1(2, "no", _::AbstractUserType)      => 5
            (2, 3, _::Real, 5)                         => 6
            [7, (1, _::AbstractUserType), _..., _::Int] => 7
            (_::String, ("5", _::UserTy2, ), 10) => 8
            _                                  => 9
        end
    end
]

records = NamedTuple{(:time_mean, :implementation, :case)}[]
for (spec_id, spec) in specs
    # group_key = string(spec_id)
    # suite[group_key] = BenchmarkGroup()
    for (impl_id, impl_fn) in implementations
        bench′ =
            @benchmark $impl_fn(sample) setup = (sample = $generate($spec)) samples = 2000
        time′ = mean(bench′.times)
        @info :bench (spec_id, impl_id, time′)
        push!(records, (time_mean = time′, implementation = impl_id, case = spec_id))
    end
end

df = DataFrame(records)
@info df

theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size = 5px,
)

report_meantime, df_time = report(
    df, Guide.title("Misc");
    benchfield = :time_mean
)
open("stats/bench-misc.txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/bench-misc.svg", 14inch, 6inch), report_meantime)

end


text:
module Pmodules

export @Parent, @P

""" Find the first undefined identifier in the lineage of another.

The given value should be an iterable of symbols, the "full name" of
an identifier. This function will find and return the first undefined
identifier in the path of the full name, or Nothing if the given identifier
is defined.
"""
function first_undefined(fullid, base_module::Module)
    search_mod = base_module

    # Make sure that the first id in fullid is the base
    # module
    if fullid[1] != nameof(base_module)
        error("Incorrect base module given to undefined search.")
    end

    for (i, modsym) = enumerate(fullid[2:end])
        if isdefined(search_mod, modsym)
            search_mod = Base.getproperty(search_mod, modsym)
        else
            return fullid[1:i + 1]
        end
    end
    return nothing
end

""" Get the object that corresponds to the given identifier.

Simply start in main and search for the given id (which is
probably a module). Kind of makes an assumption that the top of
the path has been included/imported. Generally this should only
be used to get an identifier within the current top-level package.
"""
function get_id(fullid, base_module::Module)
    search_mod = base_module

    # Make sure that the first id in fullid is the base
    # module
    if fullid[1] != nameof(base_module)
        error("Incorrect base module given to undefined search.")
    end

    for sym = fullid[2:end]
        # Might throw if this function was called inappropriately
        search_mod = getproperty(search_mod, sym)
    end
    return search_mod
end

""" Resolve relative identifier from the given base."""
function resolve_relative(maybe_relative, relative_base)
    # If we are absolute, just return the ident
    if maybe_relative[1] != :(.)
        return maybe_relative
    end
    # Otherwise, the number of dots minus 1 is the
    # number of symbols taken off the end of relative base
    count_dots = 1
    for sym = maybe_relative[2:end]
        if sym == :(.)
            count_dots += 1
        end
    end

    return (relative_base[1:end-count_dots+1]..., maybe_relative[count_dots+1:end]...)
end

""" Ensure that the given identifier name will be available for import."""
function ensure_ident(ident_name, calling_module)
    mroot = Base.moduleroot(calling_module)
    # Make sure we are dealing with an internal identifier
    if ident_name[1] ∉ (nameof(mroot), :(.))
        error("Ensuring identifier of a non-internal indentifier." *
              " Ensuring $(ident_name) in $(fullname(calling_module)).")
    end

    ident_name = resolve_relative(ident_name, fullname(calling_module))

    first_undef = first_undefined(ident_name, mroot)

    if first_undef != nothing
        # We need to find a file to include
        root_dir = dirname(pathof(mroot))

        include_pen_dir = joinpath(root_dir, map(String, ident_name[2:end - 1])...)

        # Always search for bath in case of weird edge cases so we can notify the users.

        proposed_leaf = joinpath(include_pen_dir, String(ident_name[end]) * ".jl")
        proposed_parent = joinpath(include_pen_dir, String(ident_name[end]), String(ident_name[end]) * ".jl")

        is_proposed_leaf = isfile(proposed_leaf)
        is_proposed_parent = isfile(proposed_parent)

        if is_proposed_leaf && (!is_proposed_parent)
            include_dir = proposed_leaf
        elseif is_proposed_parent && (!is_proposed_leaf)
            include_dir = proposed_parent
        elseif is_proposed_leaf && is_proposed_parent
            error("Two files exist that would result in the same logical name: $(proposed_leaf) and $(proposed_parent).")
        else
            error("Could not find file that would ensure identifier $(ident_name)")
        end

        include_mod = get_id(first_undef[1:end - 1], mroot)

        Base.include(include_mod, include_dir)
    end
    # Otherwise the identifier is already defined and thus ensured.
end

""" Gives the ultimate name in the path, without the extension (if there is one)."""
finalname(path) = splitext(basename(path))[1]

""" Determine whether the given module is a parent.

Filepath should be an absolute path.
"""
function is_parent(filepath::String, mod_fullname)::Bool
    # TODO: fix this shit
    # Perform validity check
    isabspath(filepath) || error("Given path must be an absolute path.")
    # Get the relevant parts of the filepath
    directory = basename(dirname(filepath))
    filebase = splitext(basename(filepath))[1]

    # This only will occur on a top-level declaration, so the length must be 3
    # (Base, __toplevel__, module name)
    if length(mod_fullname) == 3 && mod_fullname[1:2] == fullname(Base.__toplevel__)
        # Top level Pmodule so it must be a parent
        return true
    end

    if directory == "src"
        # If the above doesn't apply and the directory is src, then the full name with
        # one entry is the parent, and the rest are leaves.
        return length(mod_fullname) == 1
    end

    # Otherwise we can use normal logic
    return directory == filebase && filebase == String(last(mod_fullname))
end

""" Define a parent module in the Pmodule system."""
macro Parent()
    # First check to see that it satisfies the parent module requirement.
    fullmod = fullname(__module__)
    filepath = String(__source__.file)

    # Make sure we are being called from a file
    if ! isfile(filepath)
        error("Pmodule called in non-file context.")
    end

    filename = finalname(filepath)
    directory = finalname(dirname(filepath))

    # Make sure we are called in a parent module.
    # TODO: More sophisticated parent detection?
    if ! is_parent(filepath, fullmod)
        error("Pmodule called in non-parent module.")
    end

    # Find all other julia files in this directory
    fulldir = dirname(filepath)
    # Use only .jl files that are not the current one
    modfiles = [finalname(fp) for fp in readdir(fulldir)
                if (endswith(fp, ".jl") && fp != basename(filepath))]

    # Ensure all identifiers before running the import statement
    return esc(Expr(:block,
        :(import Pmodules),
        [:(Pmodules.ensure_ident($(tuple(fullmod..., Symbol(mname))),
            $(__module__)))
        for mname in modfiles]...))
end

""" Support for ensuring on internal imports.
"""
macro P(ex)
    if isa(ex, Expr)
        if ex.head == :import || ex.head == :using
            # Import expression
            return esc(process_pimport(ex, __source__, __module__))
        else
            error("Unsupported expression in Pmodule system.")
        end
    else
        error("Invalid input to Pmodule macro.")
    end
end

""" Process a Pmodule import.

First determine if the desired module/modules are defined. If they are, simply import
them. If they are not, then determine the closest common ancestor, and include the
direct child of that common ancestor. If there is no common ancestor, assume that the
module is available via the usual import mechanism and import directly.

Ultimately, the import expression is used as is. It is only parsed to perform automatic
inclusion of needed files.
"""
function process_pimport(import_ex::Expr, source::LineNumberNode, module_::Module)
    # First see if the import expression is a list of imports or an import with
    # an identification specifier. If there is more than 1, we don't allow a colon part
    allow_idents = length(import_ex.args) > 1

    # Get all uniquely needed modules
    module_specs = unique([mspec for mexpr = import_ex.args for mspec = process_mexpr(mexpr, allow_idents)])

    # For each module spec, simply remove all of those that do not share a module root with the module root
    # of this module
    rootmod = Base.nameof(Base.moduleroot(module_))
    internal_mspecs = filter(t->t[1] ∈ (rootmod, :(.)), module_specs)
    # Now ensure all identifers import the import expression
    # Assume Pmodules has been imported
    return Expr(:block, :(import Pmodules), [:(Pmodules.ensure_ident($(mname), $(module_))) for mname in internal_mspecs]..., import_ex)
end

""" Process a single part of a module import expr.

The input should one of the arguments to a valid import statement.

The return is a list of module specs that outline all of the required modules
to satisfy this import expression. The module spec is a tuple of symbols that
indicate module-submodule relationship (symbols to the right are submodules
of symbols to the left).
"""
function process_mexpr(ex::Expr, allow_idents::Bool)
    if ex.head == :(.)
        # If the head is :(.), then there will only be a single tuple return
        return [tuple(ex.args...)]
    elseif ex.head == :(:)
        # Check for colon in multiple-argument import
        if allow_idents
            error("Import statement colon identifier specifier in import with multiple arguments.")
        end

        # First we get the base identifier
        base_id = tuple(ex.args[1].args...)

        # For each id on the right side of the colon, we join it with the base module
        # Note that even if the terminal identifier of a right-side part is not a module, that is
        # okay because include logic will make it so the only conflict that will occur is when
        # a name is imported from the same module that the import resides in, which is a corner
        # case that I am okay with.
        return [tuple(base_id..., rhs_id.args...) for rhs_id in ex.args[2:end]]
    else
        error("Invalid import statement.")
    end
end

end # module


text:
# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
function ITensors.OpSum(
  ::Model"fqhe_2b_pot", n1, n2; Ly::Float64, Vs::Array{Float64,1}, prec::Float64
)
  rough_N = round(Int64, 2 * Ly)
  coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs)
  opt = optimize_coefficients(coeff; prec=prec)
  opt = filter_optimized_Hamiltonian_by_first_site(opt)
  #sorted_opt = sort_by_configuration(opt);
  #println(opt)
  return generate_Hamiltonian(opt)
end

function ITensors.MPO(::Model"fqhe_2b_pot", s::CelledVector, n::Int64; Ly, Vs, prec)
  rough_N = round(Int64, 2 * Ly)
  coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs)
  opt = optimize_coefficients(coeff; prec=prec)
  opt = filter_optimized_Hamiltonian_by_first_site(opt)
  range_model = check_max_range_optimized_Hamiltonian(opt)
  while range_model >= rough_N
    rough_N = 2 * (rough_N + 1)
    coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs)
    opt = optimize_coefficients(coeff; prec=prec)
    opt = filter_optimized_Hamiltonian_by_first_site(opt)
    range_model = check_max_range_optimized_Hamiltonian(opt)
  end
  opsum = generate_Hamiltonian(opt)

  return splitblocks(linkinds, MPO(opsum, [s[x] for x in n:(n + range_model)]))
end

#Please contact Loic Herviou before using this part of the code for production
# loic.herviou@epfl.ch
###############################
###Pseudpotentials
###############################
"""
    two_body_factorized_pseudomomentum_auxiliary(x, Ly, m)
Auxiliary function that computes the value of the factorized form of the pseudocoefficient
Inputs: x momentum `2*pi*n/Ly`, Ly the length in the periodic direction and m the moment of the pseudo coefficient
Outputs: a floating value
"""
function two_body_factorized_pseudomomentum_auxiliary(x, Ly, m::Int64)
  return 2^(0.75) * sqrt(2 * pi / Ly) * exp(-x^2) * Hermite_polynomial(2 * x; n=m) /
         sqrtfactorial(m) / pi^0.25
end

"""
    build_two_body_coefficient_pseudopotential_factorized_cylinder(Lx, Ly, N_phi, maximalMoment::Int64; prec=1e-12, fermions = true)
Build the dictionary of the factorized form of the two body pseudopotentials on the cylinder
Inputs: Lx and Ly are the dimensions on the torus (Ly the orbital basis), N_phi the number of orbitals, and m the maximal moment (starting at 0)
Outputs: a dictionary of coefficients with keys the shift of the center of mass, values are either Float64
"""
function build_two_body_coefficient_pseudopotential_factorized_cylinder(
  Lx, Ly, N_phi::Int64, maximalMoment::Int64; prec=1e-12, fermions=true
)
  if fermions
    sg_ph = -1
    admissible_m = 1:2:maximalMoment
  else
    sg_ph = 1
    admissible_m = 0:2:maximalMoment
  end
  coefficients = Dict{Float64,Array{Float64,1}}()
  for k in 0.5:0.5:(N_phi / 2 - 0.5)
    coefficients[k] = zeros(Float64, length(admissible_m))
    for (idx_m, m) in enumerate(admissible_m)
      coefficients[k][idx_m] +=
        (
          two_body_factorized_pseudomomentum_auxiliary(2 * pi * k / Ly, Ly, m) +
          sg_ph * two_body_factorized_pseudomomentum_auxiliary(-2 * pi * k / Ly, Ly, m)
        ) / 2
    end
  end
  return coefficients
end

function build_two_body_coefficient_pseudopotential_cylinder(coeff, Vs, N_phi)
  full_coeff = Dict()
  for j in 0.5:0.5:(N_phi - 1.5)
    for l in mod1(j, 1):1:min(j, N_phi - 1 - j)
      sg_1 = 1
      n_1 = mod(round(Int64, j + l), N_phi)
      n_2 = mod(round(Int64, j - l), N_phi)
      if n_1 > n_2
        sg_1 = -1
        temp = n_2
        n_2 = n_1
        n_1 = temp
      elseif n_1 == n_2
        continue
      end
      for k in mod1(j, 1):1:min(j, N_phi - 1 - j)
        sg_2 = 1
        m_1 = mod(round(Int64, j + k), N_phi)
        m_2 = mod(round(Int64, j - k), N_phi)
        if m_1 > m_2
          sg_2 = -1
          temp = m_2
          m_2 = m_1
          m_1 = temp
        elseif m_1 == m_2
          continue
        end
        if !haskey(full_coeff, (m_1, m_2, n_2, n_1))
          full_coeff[[m_1, m_2, n_2, n_1]] =
            sg_1 * sg_2 * sum([coeff[l][x] * Vs[x] * coeff[k][x] for x in 1:length(Vs)])
        else
          full_coeff[[m_1, m_2, n_2, n_1]] +=
            sg_1 * sg_2 * sum([coeff[l][x] * Vs[x] * coeff[k][x] for x in 1:length(Vs)])
        end
      end
    end
  end
  return full_coeff
end

function build_two_body_coefficient_pseudopotential(;
  r::Float64=1.0,
  Lx::Float64=-1.0,
  Ly::Float64=-1.0,
  N_phi::Int64=10,
  Vs::Array{Float64,1}=[0],
  prec=1e-12,
)
  if Lx != -1
    println("Generating pseudopotential coefficients from Lx")
    Ly = 2 * pi * N_phi / Lx
    r = Lx / Ly
  elseif Ly != -1
    println("Generating pseudopotential coefficients from Ly")
    Lx = 2 * pi * N_phi / Ly
    r = Lx / Ly
  else
    println("Generating pseudopotential coefficients from r")
    Lx = sqrt(2 * pi * N_phi * r)
    Ly = sqrt(2 * pi * N_phi / r)
  end
  println(
    string(
      "Parameters are N_phi=",
      N_phi,
      ", r=",
      round(r; digits=3),
      ", Lx =",
      round(Lx; digits=3),
      " and Ly =",
      round(Ly; digits=3),
    ),
  )
  maximalMoment = 2 * length(Vs) - 1
  temp_coeff = build_two_body_coefficient_pseudopotential_factorized_cylinder(
    Lx, Ly, N_phi, maximalMoment; prec=prec, fermions=true
  )
  return build_two_body_coefficient_pseudopotential_cylinder(temp_coeff, Vs, N_phi)
end

######################
## Generic
######################
function filter_optimized_Hamiltonian_by_first_site(coeff::Dict; n=1)
  res = Dict()
  for (k, v) in coeff
    if k[2] == n
      res[k] = v
    end
  end
  return res
end

function check_max_range_optimized_Hamiltonian(coeff::Dict)
  temp = 0
  for k in keys(coeff)
    temp = max(temp, k[end] - k[2])
  end
  return temp
end

function sort_by_configuration(coeff::Dict)
  res = Dict()
  for (k, v) in coeff
    if !haskey(res, k[1:2:end])
      res[k[1:2:end]] = Dict()
    end
    res[k[1:2:end]][k] = v
  end
  return res
end
"""
    sqrtfactorial(m)
Return `sqrt(m!)`
Inputs: `m` a number
"""
function sqrtfactorial(m)
  res = 1
  for j in 1:m
    res *= sqrt(j)
  end
  return res
end

"""
    Laguerre_polynomial(x; n=0)
Return the Laguerre polynomial of arbitrary degre
Inputs: `x` the input and `n` the degree of the Laguerre polynomial
Output: the value of the Laguerre polynomial
"""
function Laguerre_polynomial(x; n=0)
  if n == 0
    return 1
  elseif n == 1
    return 1 - x
  elseif n == 2
    return x^2 / 2 - 2x + 1
  end
  Lm2 = 1
  Lm1 = 1 - x
  Lm = x^2 / 2 - 2x + 1
  for j in 3:n
    Lm2 = Lm1
    Lm1 = Lm
    Lm = ((2j - 1 - x) * Lm1 - (j - 1) * Lm2) / j
  end
  return Lm
end

"""
    Hermite_polynomial(x; n=0)
Return the Hermite polynomial of arbitrary degre
Inputs: `x` the input and `n` the degree of the Hermite polynomial (probabilist notation)
Output: the value of the Laguerre polynomial
"""
function Hermite_polynomial(x; n=0)
  if n == 0
    return 1
  elseif n == 1
    return x
  elseif n == 2
    return x^2 - 1
  end
  Lm2 = 1
  Lm1 = x
  Lm = x^2 - 1
  for j in 3:n
    Lm2 = Lm1
    Lm1 = Lm
    Lm = x * Lm1 - (j - 1) * Lm2
  end
  return Lm
end

#########################
#Building the Hamiltonian
#########################
function get_perm!(lis, name)
  for j in 1:(length(lis) - 1)
    if lis[j] > lis[j + 1]
      c = lis[j]
      lis[j] = lis[j + 1]
      lis[j + 1] = c
      c = name[j]
      name[j] = name[j + 1]
      name[j + 1] = c
      sg = get_perm!(lis, name)
      return -sg
    end
  end
  return 1
end

function filter_op!(lis, name)
  x = 1
  while x <= length(lis) - 1
    if lis[x] == lis[x + 1]
      if name[x] == "Cdag" && name[x + 1] == "C"
        popat!(lis, x + 1)
        popat!(name, x + 1)
        name[x] = "N"
      else
        print("Wrong order in filter_op")
      end
    end
    x += 1
  end
end

function optimize_coefficients(coeff::Dict; prec=1e-12)
  optimized_dic = Dict()
  for (ke, v) in coeff
    if abs(v) < prec
      continue
    end
    if length(ke) == 4
      name = ["Cdag", "Cdag", "C", "C"]
    elseif length(ke) == 6
      name = ["Cdag", "Cdag", "Cdag", "C", "C", "C"]
    elseif length(ke) == 8
      name = ["Cdag", "Cdag", "Cdag", "Cdag", "C", "C", "C", "C"]
    else
      println("Optimization not implemented")
      continue
    end
    k = Base.copy(ke)
    sg = get_perm!(k, name)
    filter_op!(k, name)
    new_k = [isodd(n) ? name[n ÷ 2 + 1] : k[n ÷ 2] + 1 for n in 1:(2 * length(name))]
    optimized_dic[new_k] = sg * v
  end
  return optimized_dic
end

function generate_Hamiltonian(mpo::OpSum, coeff::Dict; global_factor=1, prec=1e-12)
  for (k, v) in coeff
    if abs(v) > prec
      add!(mpo, global_factor * v, k...)
    end
  end
  return mpo
end

function generate_Hamiltonian(coeff::Dict; global_factor=1, prec=1e-12)
  mpo = OpSum()
  return generate_Hamiltonian(mpo, coeff; global_factor=global_factor, prec=prec)
end


text:
using Pkg
Pkg.activate(".")
using Logging
using DaemonMode

logger = SimpleLogger(stdout, Logging.Error)

with_logger(logger) do
    DaemonMode.serve()
end


text:
module FourRoomsContUtil

using ..Resampling
import ..Resampling.FourRoomsParams
import ..Resampling: AbstractPolicy
using ..Reproduce
import StatsBase
import StatsBase: Weights
using Random
import ..JuliaRL
import ..build_algorithm_dict

include("FourRoomsCont/policies.jl")

export
    TCFourRoomsContAgent, NNFourRoomsContAgent, IntFourRoomsContAgent,
    predict, predict!

include("FourRoomsCont/agent.jl")
include("FourRoomsCont/int_agent.jl")
include("FourRoomsCont/nn_agent.jl")

include("FourRoomsCont/dmu.jl")

export
    env_settings!,
    get_experience,
    policy_settings,
    max_is_dict,
    get_policy

const constructed_policy_settings = Dict(
    "uniform"=>UniformRandomPolicy(4),
    "corners"=>FourCornersVariant(),
    "favor_down"=>FavoredRandomPolicy(4, FourRoomsParams.DOWN, 0.1),
    "random_state_variant"=>RandomStateVariant(),
    "random_state_weight_variant"=>RandomStateWeightVariant()
)

const GVFS = Dict(
    "collide_down"=>GVF(FunctionalCumulant((state_t, action_t, state_tp1, action_tp1, preds_tp1)-> (state_tp1[2] ? 1.0 : 0.0)),
                StateTerminationDiscount(0.9, (state_t, action_t, state_tp1)->(state_tp1[2])),
                FavoredRandomPolicy(4, 3, 1.0)),
    "favored_down"=>GVF(FunctionalCumulant((state_t, action_t, state_tp1, action_tp1, preds_tp1)-> (state_tp1[2] ? 1.0 : 0.0)),
                        StateTerminationDiscount(0.9, (state_t, action_t, state_tp1)->(state_tp1[2])),
                        FavoredRandomPolicy(4, 3, 0.9))
)

get_max_is_ratio(μ::RandomPolicy, π::RandomPolicy) =
    maximum(π.probabilities./μ.probabilities)
get_max_is_ratio(μ::FRStateDependentPolicy, π::RandomPolicy) =
    maximum([maximum(π.probabilities./μ.weights[state]) for state in μ.states])

function env_settings!(s::ArgParseSettings)
    @add_arg_table s begin
        "--policy"
        help = "The policy settings: $(keys(constructed_policy_settings))"
        range_tester=(p)->(p ∈ keys(constructed_policy_settings))
        required=true
        "--gvf"
        help = "The gvf to learn: $(keys(GVFS))"
        range_tester=(gvf)->(gvf ∈ keys(GVFS))
        required=true
    end
end

function get_policy(parsed::Dict)
    return constructed_policy_settings[parsed["policy"]]
end

import ProgressMeter

function get_dmu_mu_dist(env::FourRoomsCont, policy; seed=0, num_runs=1e5, num_steps=1000)
    rng = Random.MersenneTwister(seed)

    heat_map = zeros(size(env)...)

    ProgressMeter.@showprogress 0.1 "Run:" for i in 1:num_runs
        _, state = start!(env; rng=rng)
        action = StatsBase.sample(rng, policy, state[1])
        for i in 1:num_steps
            _, state, _, _ = step!(env, action; rng=rng)
            action = StatsBase.sample(rng, policy, state[1])
        end
        heat_map[Resampling.project(env, state[1])...] += 1
    end
    return heat_map./sum(heat_map)
end

function sample_according_to_dmu(env::FourRoomsCont, policy_str::AbstractString, num_states; rng=Random.GLOBAL_RNG)
    dmu = mu_dmu_dict[policy_str]
    indicies = reshape(CartesianIndices(dmu), :)
    w = StatsBase.Weights(dmu[indicies])

    prj_states = StatsBase.sample(rng, indicies, w, num_states)

    [[rand(rng) + (s[1]-1), rand(rng) + (s[2]-1)] for s in prj_states]
end

end


text:
using MomentTensors, Test, StaticArrays

@testset "Construction" begin
    @testset "Direct" begin
        @test MT(1, 2, 3, 4, 5, 6) isa MT{Float64}
        @test MT{Float32}(1, 2, 3, 4, 5, 6) isa MT{Float32}
        @testset "$T" for T in (Float32, Float64)
            @test eltype(MT{T}(1, 2, 3, 4, 5, 6)) == T
            @test eltype(MT{Complex{T}}(1, 2, 3, 4, 5, 6)) == Complex{T}
        end
        @test MT(1, 2, 3, 4, 5, 6) == MT(@SVector[1, 2, 3, 4, 5, 6])
        @test MT(Int32[1, 2, 3, 4, 5, 6]) ==
            MT{float(Int32)}(@SVector[1.f0, 2.f0, 3.f0, 4.f0, 5.f0, 6.f0])
        @test MT(Float32[1, 2, 3, 4, 5, 6]) == MT{Float32}(1, 2, 3, 4, 5, 6)
        @test MT([11 12 13; 12 22 23; 13 23 33]) == MT(11, 22, 33, 12, 13, 23)
        @test MT(11, 22, 33, 12, 13, 23) == MT([11, 22, 33, 12, 13, 23])
        @test_throws ArgumentError MT([1 2 3; 4 5 6])
        # Warn for asymmetric matrices
        @test_logs (:warn,
            "Supplied tensor is not symmetric; using upper elements") MT([1 2 3
                                                                          4 5 6
                                                                          7 8 9])
        # Do not warn with `warn=false`
        @test_logs MT([1 2 3; 4 5 6; 7 8 9], false)
    end

    @testset "Strike-dip-rake" begin
        @test MT(30, 60, -45, 7e15) ≈
            MT(1e16.*(-0.429, -0.264, 0.693, -0.338, -0.091, -0.029)...) atol=0.01e16
    end
end


text:
# SAM Header
# ==========

struct Header
    metainfo::Vector{MetaInfo}
end

"""
    SAM.Header()

Create an empty header.
"""
function Header()
    return Header(MetaInfo[])
end

function Base.:(==)(a::Header, b::Header)
    return a.metainfo == b.metainfo
end

function Base.copy(header::Header)
    return Header(header.metainfo)
end

function Base.eltype(::Type{Header})
    return MetaInfo
end

function Base.length(header::Header)
    return length(header.metainfo)
end

function Base.iterate(header::Header, i=1)
    if i > length(header.metainfo)
        return nothing
    end
    return header.metainfo[i], i + 1
end

"""
    find(header::Header, key::AbstractString)::Vector{MetaInfo}

Find metainfo objects satisfying `SAM.tag(metainfo) == key`.
"""
function Base.findall(header::Header, key::AbstractString)::Vector{MetaInfo}
    return filter(m -> isequalkey(m, key), header.metainfo)
end

function Base.pushfirst!(header::Header, metainfo::MetaInfo)
    pushfirst!(header.metainfo, metainfo)
    return header
end

function Base.push!(header::Header, metainfo::MetaInfo)
    push!(header.metainfo, metainfo)
    return header
end


text:
# Raw memory management

export Mem, attribute, attribute!, memory_type, is_managed

module Mem

using ..CUDA
using ..CUDA: @enum_without_prefix, CUstream, CUdevice, CuDim3, CUarray, CUarray_format

using Base: @deprecate_binding

# TODO: needs another redesign
#
# - Buffers should be typed, since ArrayBuffers are tied to a format.
#   untyped buffers can use T=UInt8 or T=Nothing and reinterpret the output pointer
#   (this is impossible with textures, though)
# - copyto! methods should take a Buffer so that we can populate the memcpy structs directly
# - allocate buffers on construction, instead of using an alloc method? free on finalizer?
# - have CuArray contain a buffer, use that for dispatch.
# - trait to determine which pointers a buffer can yield?

#
# untyped buffers
#

abstract type Buffer end

# expected interface:
# - similar()
# - ptr, bytesize and ctx fields
# - convert() to Ptr and CuPtr

Base.pointer(buf::Buffer) = buf.ptr

Base.sizeof(buf::Buffer) = buf.bytesize

CUDA.device(buf::Buffer) = device(buf.ctx)

# ccall integration
#
# taking the pointer of a buffer means returning the underlying pointer,
# and not the pointer of the buffer object itself.
Base.unsafe_convert(T::Type{<:Union{Ptr,CuPtr,CuArrayPtr}}, buf::Buffer) = convert(T, buf)

## device buffer

"""
    Mem.DeviceBuffer
    Mem.Device

A buffer of device memory residing on the GPU.
"""
struct DeviceBuffer <: Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int
    ctx::CuContext
end

Base.similar(buf::DeviceBuffer, ptr::CuPtr{Cvoid}=pointer(buf),
             bytesize::Int=sizeof(buf), ctx::CuContext=buf.ctx) =
    DeviceBuffer(ptr, bytesize, ctx)

Base.convert(::Type{<:Ptr}, buf::DeviceBuffer) =
    throw(ArgumentError("cannot take the CPU address of a GPU buffer"))

Base.convert(::Type{CuPtr{T}}, buf::DeviceBuffer) where {T} =
    convert(CuPtr{T}, pointer(buf))

"""
    Mem.alloc(DeviceBuffer, bytesize::Integer)

Allocate `bytesize` bytes of memory on the device. This memory is only accessible on the
GPU, and requires explicit calls to `unsafe_copyto!`, which wraps `cuMemcpy`,
for access on the CPU.
"""
function alloc(::Type{DeviceBuffer}, bytesize::Integer)
    bytesize == 0 && return DeviceBuffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CUDA.CUdeviceptr}()
    CUDA.cuMemAlloc(ptr_ref, bytesize)

    return DeviceBuffer(reinterpret(CuPtr{Cvoid}, ptr_ref[]), bytesize, CuCurrentContext())
end

function free(buf::DeviceBuffer)
    if pointer(buf) != CU_NULL
        CUDA.cuMemFree(buf)
    end
end

## host buffer

"""
    Mem.HostBuffer
    Mem.Host

A buffer of pinned memory on the CPU, possible accessible on the GPU.
"""
struct HostBuffer <: Buffer
    ptr::Ptr{Cvoid}
    bytesize::Int
    ctx::CuContext

    mapped::Bool
end

Base.similar(buf::HostBuffer, ptr::Ptr{Cvoid}=pointer(buf), bytesize::Int=sizeof(buf),
             ctx::CuContext=buf.ctx, mapped::Bool=buf.mapped) =
    HostBuffer(ptr, bytesize, ctx, mapped)

Base.convert(::Type{Ptr{T}}, buf::HostBuffer) where {T} =
    convert(Ptr{T}, pointer(buf))

function Base.convert(::Type{CuPtr{T}}, buf::HostBuffer) where {T}
    if buf.mapped
        pointer(buf) == C_NULL && return convert(CuPtr{T}, CU_NULL)
        ptr_ref = Ref{CuPtr{Cvoid}}()
        CUDA.cuMemHostGetDevicePointer(ptr_ref, pointer(buf), #=flags=# 0)
        convert(CuPtr{T}, ptr_ref[])
    else
        throw(ArgumentError("cannot take the GPU address of a pinned but not mapped CPU buffer"))
    end
end

@deprecate_binding HOSTALLOC_DEFAULT 0 false
const HOSTALLOC_PORTABLE = CUDA.CU_MEMHOSTALLOC_PORTABLE
const HOSTALLOC_DEVICEMAP = CUDA.CU_MEMHOSTALLOC_DEVICEMAP
const HOSTALLOC_WRITECOMBINED = CUDA.CU_MEMHOSTALLOC_WRITECOMBINED

"""
    Mem.alloc(HostBuffer, bytesize::Integer, [flags])

Allocate `bytesize` bytes of page-locked memory on the host. This memory is accessible from
the CPU, and makes it possible to perform faster memory copies to the GPU. Furthermore, if
`flags` is set to `HOSTALLOC_DEVICEMAP` the memory is also accessible from the GPU.
These accesses are direct, and go through the PCI bus.
If `flags` is set to `HOSTALLOC_PORTABLE`, the memory is considered mapped by all CUDA contexts,
not just the one that created the memory, which is useful if the memory needs to be accessed from
multiple devices. Multiple `flags` can be set at one time using a bytewise `OR`:

    flags = HOSTALLOC_PORTABLE | HOSTALLOC_DEVICEMAP

"""
function alloc(::Type{HostBuffer}, bytesize::Integer, flags=0)
    bytesize == 0 && return HostBuffer(C_NULL, 0, CuContext(C_NULL), false)

    ptr_ref = Ref{Ptr{Cvoid}}()
    CUDA.cuMemHostAlloc(ptr_ref, bytesize, flags)

    mapped = (flags & HOSTALLOC_DEVICEMAP) != 0
    return HostBuffer(ptr_ref[], bytesize, CuCurrentContext(), mapped)
end

const HOSTREGISTER_PORTABLE = CUDA.CU_MEMHOSTREGISTER_PORTABLE
const HOSTREGISTER_DEVICEMAP = CUDA.CU_MEMHOSTREGISTER_DEVICEMAP
const HOSTREGISTER_IOMEMORY = CUDA.CU_MEMHOSTREGISTER_IOMEMORY

"""
    Mem.register(HostBuffer, ptr::Ptr, bytesize::Integer, [flags])

Page-lock the host memory pointed to by `ptr`. Subsequent transfers to and from devices will
be faster, and can be executed asynchronously. If the `HOSTREGISTER_DEVICEMAP` flag is
specified, the buffer will also be accessible directly from the GPU.
These accesses are direct, and go through the PCI bus.
If the `HOSTREGISTER_PORTABLE` flag is specified, any CUDA context can access the memory.
"""
function register(::Type{HostBuffer}, ptr::Ptr, bytesize::Integer, flags=0)
    bytesize == 0 && throw(ArgumentError("Cannot register an empty range of memory."))

    CUDA.cuMemHostRegister(ptr, bytesize, flags)

    mapped = (flags & HOSTREGISTER_DEVICEMAP) != 0
    return HostBuffer(ptr, bytesize, CuCurrentContext(), mapped)
end

"""
    Mem.unregister(HostBuffer)

Unregisters a memory range that was registered with [`Mem.register`](@ref).
"""
function unregister(buf::HostBuffer)
    CUDA.cuMemHostUnregister(buf)
end

function free(buf::HostBuffer)
    if pointer(buf) != CU_NULL
        CUDA.cuMemFreeHost(buf)
    end
end

## unified buffer

"""
    Mem.UnifiedBuffer
    Mem.Unified

A managed buffer that is accessible on both the CPU and GPU.
"""
struct UnifiedBuffer <: Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int
    ctx::CuContext
end

Base.similar(buf::UnifiedBuffer, ptr::CuPtr{Cvoid}=pointer(buf),
             bytesize::Int=sizeof(buf), ctx::CuContext=buf.ctx) =
    UnifiedBuffer(ptr, bytesize, ctx)

Base.convert(::Type{Ptr{T}}, buf::UnifiedBuffer) where {T} =
    convert(Ptr{T}, reinterpret(Ptr{Cvoid}, pointer(buf)))

Base.convert(::Type{CuPtr{T}}, buf::UnifiedBuffer) where {T} =
    convert(CuPtr{T}, pointer(buf))

@enum_without_prefix CUDA.CUmemAttach_flags CU_MEM_

"""
    Mem.alloc(UnifiedBuffer, bytesize::Integer, [flags::CUmemAttach_flags])

Allocate `bytesize` bytes of unified memory. This memory is accessible from both the CPU and
GPU, with the CUDA driver automatically copying upon first access.
"""
function alloc(::Type{UnifiedBuffer}, bytesize::Integer,
              flags::CUDA.CUmemAttach_flags=ATTACH_GLOBAL)
    bytesize == 0 && return UnifiedBuffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CuPtr{Cvoid}}()
    CUDA.cuMemAllocManaged(ptr_ref, bytesize, flags)

    return UnifiedBuffer(ptr_ref[], bytesize, CuCurrentContext())
end

function free(buf::UnifiedBuffer)
    if pointer(buf) != CU_NULL
        CUDA.cuMemFree(buf)
    end
end

"""
    prefetch(::UnifiedBuffer, [bytes::Integer]; [device::CuDevice], [stream::CuStream])

Prefetches memory to the specified destination device.
"""
function prefetch(buf::UnifiedBuffer, bytes::Integer=sizeof(buf);
                  device::CuDevice=device(buf), stream::CuStream=CuDefaultStream())
    bytes > sizeof(buf) && throw(BoundsError(buf, bytes))
    CUDA.cuMemPrefetchAsync(buf, bytes, device, stream)
end

@enum_without_prefix CUDA.CUmem_advise CU_MEM_

"""
    advise(::UnifiedBuffer, advice::CUDA.CUmem_advise, [bytes::Integer]; [device::CuDevice])

Advise about the usage of a given memory range.
"""
function advise(buf::UnifiedBuffer, advice::CUDA.CUmem_advise, bytes::Integer=sizeof(buf);
                device::CuDevice=device(buf))
    bytes > sizeof(buf) && throw(BoundsError(buf, bytes))
    CUDA.cuMemAdvise(buf, bytes, advice, device)
end

## array buffer

mutable struct ArrayBuffer{T,N} <: Buffer
    ptr::CuArrayPtr{T}
    dims::Dims{N}
    ctx::CuContext
end

Base.convert(::Type{CuArrayPtr{T}}, buf::ArrayBuffer) where {T} =
    convert(CuArrayPtr{T}, pointer(buf))

function alloc(::Type{<:ArrayBuffer{T}}, dims::Dims{N}) where {T,N}
    format = convert(CUarray_format, eltype(T))

    if N == 2
        width, height = dims
        depth = 0
        @assert 1 <= width "CUDA 2D array (texture) width must be >= 1"
        # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_WIDTH
        @assert 1 <= height "CUDA 2D array (texture) height must be >= 1"
        # @assert height <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_HEIGHT
    elseif N == 3
        width, height, depth = dims
        @assert 1 <= width "CUDA 3D array (texture) width must be >= 1"
        # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_WIDTH
        @assert 1 <= height "CUDA 3D array (texture) height must be >= 1"
        # @assert height <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_HEIGHT
        @assert 1 <= depth "CUDA 3D array (texture) depth must be >= 1"
        # @assert depth <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_DEPTH
    elseif N == 1
        width = dims[1]
        height = depth = 0
        @assert 1 <= width "CUDA 1D array (texture) width must be >= 1"
        # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_WIDTH
    else
        "CUDA arrays (texture memory) can only have 1, 2 or 3 dimensions"
    end

    allocateArray_ref = Ref(CUDA.CUDA_ARRAY3D_DESCRIPTOR(
        width, # Width::Csize_t
        height, # Height::Csize_t
        depth, # Depth::Csize_t
        format, # Format::CUarray_format
        UInt32(CUDA.nchans(T)), # NumChannels::UInt32
        0))

    handle_ref = Ref{CUarray}()
    CUDA.cuArray3DCreate(handle_ref, allocateArray_ref)
    ptr = reinterpret(CuArrayPtr{T}, handle_ref[])

    return ArrayBuffer{T,N}(ptr, dims, CuCurrentContext())
end

function free(buf::ArrayBuffer)
    CUDA.cuArrayDestroy(buf)
end

## convenience aliases

const Device  = DeviceBuffer
const Host    = HostBuffer
const Unified = UnifiedBuffer
const Array   = ArrayBuffer

#
# typed pointers
#

## initialization

"""
    Mem.set!(buf::CuPtr, value::Union{UInt8,UInt16,UInt32}, len::Integer;
             async::Bool=false, stream::CuStream)

Initialize device memory by copying `val` for `len` times. Executed asynchronously if
`async` is true, in which case a valid `stream` is required.
"""
set!

for T in [UInt8, UInt16, UInt32]
    bits = 8*sizeof(T)
    fn_sync = Symbol("cuMemsetD$(bits)")
    fn_async = Symbol("cuMemsetD$(bits)Async")
    @eval function set!(ptr::CuPtr{$T}, value::$T, len::Integer;
                        async::Bool=false, stream::Union{Nothing,CuStream}=nothing)
        if async
          stream===nothing &&
              throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(CUDA, fn_async))(ptr, value, len, stream)
        else
          stream===nothing ||
              throw(ArgumentError("Synchronous memory operations cannot be issues on a stream."))
            $(getproperty(CUDA, fn_sync))(ptr, value, len)
        end
    end
end

## copy operations

for (f, srcPtrTy, dstPtrTy) in (("cuMemcpyDtoH", CuPtr,      Ptr),
                                ("cuMemcpyHtoD", Ptr,        CuPtr),
                                ("cuMemcpyDtoD", CuPtr,      CuPtr),
                               )
    @eval function Base.unsafe_copyto!(dst::$dstPtrTy{T}, src::$srcPtrTy{T}, N::Integer;
                                       stream::Union{Nothing,CuStream}=nothing,
                                       async::Bool=false) where T
        if async
            stream===nothing &&
                throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(CUDA, Symbol(f * "Async")))(dst, src, N*sizeof(T), stream)
        else
            stream===nothing ||
                throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
            $(getproperty(CUDA, Symbol(f)))(dst, src, N*sizeof(T))
        end
        return dst
    end
end

function Base.unsafe_copyto!(dst::CuArrayPtr{T}, doffs::Integer, src::Ptr{T}, N::Integer;
                             stream::Union{Nothing,CuStream}=nothing,
                             async::Bool=false) where T
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        CUDA.cuMemcpyHtoAAsync(dst, doffs, src, N*sizeof(T), stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        CUDA.cuMemcpyHtoA(dst, doffs, src, N*sizeof(T))
    end
end

function Base.unsafe_copyto!(dst::Ptr{T}, src::CuArrayPtr{T}, soffs::Integer, N::Integer;
                             stream::Union{Nothing,CuStream}=nothing,
                             async::Bool=false) where T
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        CUDA.cuMemcpyAtoHAsync(dst, src, soffs, N*sizeof(T), stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        CUDA.cuMemcpyAtoH(dst, src, soffs, N*sizeof(T))
    end
end

Base.unsafe_copyto!(dst::CuArrayPtr{T}, doffs::Integer, src::CuPtr{T}, N::Integer) where {T} =
    CUDA.cuMemcpyDtoA(dst, doffs, src, N*sizeof(T))

Base.unsafe_copyto!(dst::CuPtr{T}, src::CuArrayPtr{T}, soffs::Integer, N::Integer) where {T} =
    CUDA.cuMemcpyAtoD(dst, src, soffs, N*sizeof(T))

Base.unsafe_copyto!(dst::CuArrayPtr, src, N::Integer; kwargs...) =
    Base.unsafe_copyto!(dst, 0, src, N; kwargs...)

Base.unsafe_copyto!(dst, src::CuArrayPtr, N::Integer; kwargs...) =
    Base.unsafe_copyto!(dst, src, 0, N; kwargs...)

function unsafe_copy2d!(dst::Union{Ptr{T},CuPtr{T},CuArrayPtr{T}}, dstTyp::Type{<:Buffer},
                        src::Union{Ptr{T},CuPtr{T},CuArrayPtr{T}}, srcTyp::Type{<:Buffer},
                        width::Integer, height::Integer=1;
                        dstPos::CuDim=(1,1), srcPos::CuDim=(1,1),
                        dstPitch::Integer=0, srcPitch::Integer=0,
                        async::Bool=false, stream::Union{Nothing,CuStream}=nothing) where T
    srcPos = CUDA.CuDim3(srcPos)
    @assert srcPos.z == 1
    dstPos = CUDA.CuDim3(dstPos)
    @assert dstPos.z == 1

    srcMemoryType, srcHost, srcDevice, srcArray = if srcTyp == Host
        CUDA.CU_MEMORYTYPE_HOST,
        src::Ptr,
        0,
        0
    elseif srcTyp == Mem.Device
        CUDA.CU_MEMORYTYPE_DEVICE,
        0,
        src::CuPtr,
        0
    elseif srcTyp == Mem.Unified
        CUDA.CU_MEMORYTYPE_UNIFIED,
        0,
        reinterpret(CuPtr{Cvoid}, src),
        0
    elseif srcTyp == Mem.Array
        CUDA.CU_MEMORYTYPE_ARRAY,
        0,
        0,
        src::CuArrayPtr
    end

    dstMemoryType, dstHost, dstDevice, dstArray = if dstTyp == Host
        CUDA.CU_MEMORYTYPE_HOST,
        dst::Ptr,
        0,
        0
    elseif dstTyp == Mem.Device
        CUDA.CU_MEMORYTYPE_DEVICE,
        0,
        dst::CuPtr,
        0
    elseif dstTyp == Mem.Unified
        CUDA.CU_MEMORYTYPE_UNIFIED,
        0,
        reinterpret(CuPtr{Cvoid}, dst),
        0
    elseif dstTyp == Mem.Array
        CUDA.CU_MEMORYTYPE_ARRAY,
        0,
        0,
        dst::CuArrayPtr
    end

    params_ref = Ref(CUDA.CUDA_MEMCPY2D(
        # source
        (srcPos.x-1)*sizeof(T), srcPos.y-1,
        srcMemoryType, srcHost, srcDevice, srcArray,
        srcPitch,
        # destination
        (dstPos.x-1)*sizeof(T), dstPos.y-1,
        dstMemoryType, dstHost, dstDevice, dstArray,
        dstPitch,
        # extent
        width*sizeof(T), height
    ))
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        CUDA.cuMemcpy2DAsync(params_ref, stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        CUDA.cuMemcpy2D(params_ref)
    end
end

"""
    unsafe_copy3d!(dst, dstTyp, src, srcTyp, width, height=1, depth=1;
                   dstPos=(1,1,1), dstPitch=0, dstHeight=0,
                   srcPos=(1,1,1), srcPitch=0, srcHeight=0,
                   async=false, stream=nothing)

Perform a 3D memory copy between pointers `src` and `dst`, at respectively position `srcPos`
and `dstPos` (1-indexed). Both pitch and destination can be specified for both the source
and destination; consult the CUDA documentation for more details. This call is executed
asynchronously if `async` is set, in which case `stream` needs to be a valid CuStream.
"""
function unsafe_copy3d!(dst::Union{Ptr{T},CuPtr{T},CuArrayPtr{T}}, dstTyp::Type{<:Buffer},
                        src::Union{Ptr{T},CuPtr{T},CuArrayPtr{T}}, srcTyp::Type{<:Buffer},
                        width::Integer, height::Integer=1, depth::Integer=1;
                        dstPos::CuDim=(1,1,1), srcPos::CuDim=(1,1,1),
                        dstPitch::Integer=0, dstHeight::Integer=0,
                        srcPitch::Integer=0, srcHeight::Integer=0,
                        async::Bool=false, stream::Union{Nothing,CuStream}=nothing) where T
    srcPos = CUDA.CuDim3(srcPos)
    dstPos = CUDA.CuDim3(dstPos)

    srcMemoryType, srcHost, srcDevice, srcArray = if srcTyp == Host
        CUDA.CU_MEMORYTYPE_HOST,
        src::Ptr,
        0,
        0
    elseif srcTyp == Mem.Device
        CUDA.CU_MEMORYTYPE_DEVICE,
        0,
        src::CuPtr,
        0
    elseif srcTyp == Mem.Unified
        CUDA.CU_MEMORYTYPE_UNIFIED,
        0,
        reinterpret(CuPtr{Cvoid}, src),
        0
    elseif srcTyp == Mem.Array
        CUDA.CU_MEMORYTYPE_ARRAY,
        0,
        0,
        src::CuArrayPtr
    end

    dstMemoryType, dstHost, dstDevice, dstArray = if dstTyp == Host
        CUDA.CU_MEMORYTYPE_HOST,
        dst::Ptr,
        0,
        0
    elseif dstTyp == Mem.Device
        CUDA.CU_MEMORYTYPE_DEVICE,
        0,
        dst::CuPtr,
        0
    elseif dstTyp == Mem.Unified
        CUDA.CU_MEMORYTYPE_UNIFIED,
        0,
        reinterpret(CuPtr{Cvoid}, dst),
        0
    elseif dstTyp == Mem.Array
        CUDA.CU_MEMORYTYPE_ARRAY,
        0,
        0,
        dst::CuArrayPtr
    end

    params_ref = Ref(CUDA.CUDA_MEMCPY3D(
        # source
        (srcPos.x-1)*sizeof(T), srcPos.y-1, srcPos.z-1,
        0, # LOD
        srcMemoryType, srcHost, srcDevice, srcArray,
        C_NULL, # reserved
        srcPitch, srcHeight,
        # destination
        (dstPos.x-1)*sizeof(T), dstPos.y-1, dstPos.z-1,
        0, # LOD
        dstMemoryType, dstHost, dstDevice, dstArray,
        C_NULL, # reserved
        dstPitch, dstHeight,
        # extent
        width*sizeof(T), height, depth
    ))
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        CUDA.cuMemcpy3DAsync(params_ref, stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        CUDA.cuMemcpy3D(params_ref)
    end
end

#
# auxiliary functionality
#

## memory info

function info()
    free_ref = Ref{Csize_t}()
    total_ref = Ref{Csize_t}()
    CUDA.cuMemGetInfo(free_ref, total_ref)
    return convert(Int, free_ref[]), convert(Int, total_ref[])
end

end # module Mem

"""
    available_memory()

Returns the available_memory amount of memory (in bytes), available for allocation by the CUDA context.
"""
available_memory() = Mem.info()[1]

"""
    total_memory()

Returns the total amount of memory (in bytes), available for allocation by the CUDA context.
"""
total_memory() = Mem.info()[2]

## pointer attributes

"""
    attribute(X, ptr::Union{Ptr,CuPtr}, attr)

Returns attribute `attr` about pointer `ptr`. The type of the returned value depends on the
attribute, and as such must be passed as the `X` parameter.
"""
function attribute(X::Type, ptr::Union{Ptr{T},CuPtr{T}}, attr::CUpointer_attribute) where {T}
    ptr = reinterpret(CuPtr{T}, ptr)
    data_ref = Ref{X}()
    cuPointerGetAttribute(data_ref, attr, ptr)
    return data_ref[]
end

"""
    attribute!(ptr::Union{Ptr,CuPtr}, attr, val)

Sets attribute` attr` on a pointer `ptr` to `val`.
"""
function attribute!(ptr::Union{Ptr{T},CuPtr{T}}, attr::CUpointer_attribute, val) where {T}
    ptr = reinterpret(CuPtr{T}, ptr)
    cuPointerSetAttribute(Ref(val), attr, ptr)
    return
end

@enum_without_prefix CUpointer_attribute CU_

# some common attributes

@enum_without_prefix CUmemorytype CU_
memory_type(x) = CUmemorytype(attribute(Cuint, x, POINTER_ATTRIBUTE_MEMORY_TYPE))

is_managed(x) = convert(Bool, attribute(Cuint, x, POINTER_ATTRIBUTE_IS_MANAGED))

## shared texture/array stuff

function Base.convert(::Type{CUarray_format}, T::Type)
    if T === UInt8
        return CU_AD_FORMAT_UNSIGNED_INT8
    elseif T === UInt16
        return CU_AD_FORMAT_UNSIGNED_INT16
    elseif T === UInt32
        return CU_AD_FORMAT_UNSIGNED_INT32
    elseif T === Int8
        return CU_AD_FORMAT_SIGNED_INT8
    elseif T === Int16
        return CU_AD_FORMAT_SIGNED_INT16
    elseif T === Int32
        return CU_AD_FORMAT_SIGNED_INT32
    elseif T === Float16
        return CU_AD_FORMAT_HALF
    elseif T === Float32
        return CU_AD_FORMAT_FLOAT
    else
        throw(ArgumentError("CUDA does not support texture arrays for element type $T."))
    end
end

nchans(::Type{<:NTuple{C}}) where {C} = C
nchans(::Type) = 1


text:
using HorizonSideRobots

function krest(r::Robot)
    for side in [Nord,West,Sud,Ost]
        movekrest(r,side)
        reversekrest(r,side)
    end
    putmarker!(r)
end

function movekrest(r,side)
    while !isborder(r,side)
        move!(r,side)
        putmarker!(r)
    end
end

function reversekrest(r,side)
    while ismarker(r)==1
        move!(r,HorizonSide((Int(side)+2)%4))
    end
end

krest(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "utils" begin
    include("test_error_log.jl")
    include("test_array_utils.jl")
    include("test_coord_utils.jl")
end