text:
# a helper function that checks if there is enough memory for the output data frame
#  If type is not Number, probably something is wrong about setting the variables and it is better to be conservative. here 10^7 threshhold is arbitarary
_check_allocation_limit(T, rows, cols) = T !== Union{} && T <: Number ? sizeof(T)*rows*cols / Base.Sys.total_memory() : rows*cols/10^7

_default_renamecolid_function_withoutid(x, y) = "_c" * string(x)
_default_renamecolid_function_withid(x, y) = identity(string(values(x)))
_default_renamerowid_function(x) = identity(x)
# handling simplest case
function _simple_ds_transpose!(outx, inx, i)
    @views copy!(outx[i,:], inx)
end

function _generate_col_row_names(renamecolid, renamerowid, ids, dsnames; max_length = 0)
    local new_col_names
    try
        new_col_names = map(x -> renamecolid(x, dsnames), ids)
    catch e
        if (e isa MethodError)
            new_col_names = map(renamecolid, ids)
        else
            rethrow(e)
        end
    end
    row_names = map(renamerowid, dsnames)
    row_names = allowmissing(row_names)
    r_n_l = length(row_names)
    if r_n_l < max_length
        resize!(row_names, max_length)
        for i in r_n_l+1:max_length
            row_names[i] = missing
        end
    end
    (new_col_names, row_names)
end

function _simple_transpose_ds_generate(T, in_cols, row_names_length, new_col_names, variable_name, threads)
    outputmat = Matrix{T}(undef, row_names_length, length(new_col_names))
    if threads
        Threads.@threads for i in 1:length(in_cols)
            _simple_ds_transpose!(outputmat, in_cols[i], i)
        end
    else
        for i in 1:length(in_cols)
            _simple_ds_transpose!(outputmat, in_cols[i], i)
        end
    end
    outputmat

end

function _find_id_unique_values(ds, ididx::MultiColumnIndex, perms; mapformats = true)
    groups, gslots, ngroups = _gather_groups(ds[perms, ididx], :, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64), mapformats = mapformats)
    res = falses(nrow(ds))
    seen_groups = falses(ngroups)
    @inbounds for i in 1:length(res)
        !seen_groups[groups[i]] ? (seen_groups[groups[i]] = true; res[i] = true) : nothing
    end
    return groups, res
end
_find_id_unique_values(ds, ididx::ColumnIndex, perms; mapformats = true) = _find_id_unique_values(ds, [ididx], perms; mapformats = mapformats)

"""
    transpose(ds::AbstractDataset, cols;
        id = nothing,
        renamecolid = (x -> "_c" * string(x)),
        renamerowid = identity,
        variable_name = "_variables_",
        default = missing,
        mapformats = true,
        threads = true)

transpose `ds[!, cols]`. When `id` is set, the values of `ds[!, id]` will be used to label the columns in the new data set. The function uses the `renamecolid` function to generate the new columns labels. The `renamerowid` function is applied to stringified names of `ds[!, cols]` and these are attached to the output as a new column with the label `variable_name`. When a grouped dataset (e.g. by using `groupby!(ds, gcols)`, ,`groupby(ds, gcols)`, or `gatherby(ds, gcols)`) passed as the first argument the transposing is done within each group constructed by grouping columns. If the number of rows in a group is smaller than other groups, the extra columns for that group in the output data frame is filled with `missing` by default, however, the default value can be changed by passing `default = ` argument.

When `cols` is a Tuple of column indices, the transposing is done for each set of indices and at the end all transposed columns are horizontally concatenated. In this case, by default, the `variable_name` keyword argument is set to `nothing`.

* `renamecolid`: When `id` is not set, the argument to `renamecolid` must be an `Int`. And when `id` is set, the `renamecolid` will be applied to each row of `ds[!, id]` as Tuple.
* When `id` is set, `renamecolid` is defined as `x -> identity(string(values(x)))`
* By default, `transpose` uses the formatted value for the id variables, to change this the `mapformats = false` can be used.
* When `threads = true`, `transpose` uses all available cores to `Julia` to do the computations.

# Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,4], x2 = [1,4,9,16])
4×2 Dataset
 Row │ x1        x2
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         1
   2 │        2         4
   3 │        3         9
   4 │        4        16

julia> transpose(ds, [:x1,:x2])
2×5 Dataset
 Row │ _variables_  _c1       _c2       _c3       _c4
     │ identity     identity  identity  identity  identity
     │ String?      Int64?    Int64?    Int64?    Int64?
─────┼─────────────────────────────────────────────────────
   1 │ x1                  1         2         3         4
   2 │ x2                  1         4         9        16

julia> pop = Dataset(country = ["c1","c1","c2","c2","c3","c3"],
                             sex = repeat(["male", "female"],3),
                             pop_2000 = [100, 120, 150, 155, 170, 190],
                             pop_2010 = [110, 120, 155, 160, 178, 200],
                             pop_2020 = [115, 130, 161, 165, 180, 203])
6×5 Dataset
 Row │ country   sex       pop_2000  pop_2010  pop_2020
     │ identity  identity  identity  identity  identity
     │ String?   String?   Int64?    Int64?    Int64?
─────┼──────────────────────────────────────────────────
   1 │ c1        male           100       110       115
   2 │ c1        female         120       120       130
   3 │ c2        male           150       155       161
   4 │ c2        female         155       160       165
   5 │ c3        male           170       178       180
   6 │ c3        female         190       200       203

julia> groupby!(pop, :country);

julia> transpose(pop, r"pop_",
                       id = :sex, variable_name = "year",
                       renamerowid = x -> replace(x, "pop_" => ""),
                       renamecolid = x -> x * "_pop")
9×4 Dataset
 Row │ country   year      male_pop  female_pop
     │ identity  identity  identity  identity
     │ String?   String?   Int64?    Int64?
─────┼──────────────────────────────────────────
   1 │ c1        2000           100         120
   2 │ c1        2010           110         120
   3 │ c1        2020           115         130
   4 │ c2        2000           150         155
   5 │ c2        2010           155         160
   6 │ c2        2020           161         165
   7 │ c3        2000           170         190
   8 │ c3        2010           178         200
   9 │ c3        2020           180         203
```
"""
Base.transpose(::Dataset, cols; [id , renamecolid , renamerowid , variable_name, default, threads, mapformats])

function ds_transpose(ds, cols::Union{Tuple, MultiColumnIndex}; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = "_variables_", threads = true, mapformats = true)
    if cols isa Tuple
        tcols = [cols[j] isa ColumnIndex ? index(ds)[[cols[j]]] : multiple_getindex(index(ds), cols[j]) for j in 1:length(cols)]
    else
        tcols = [multiple_getindex(index(ds), cols)]
    end
    max_num_col = maximum(length, tcols)
    if variable_name isa AbstractString || variable_name isa Symbol || variable_name === nothing
        var_name = repeat([variable_name], length(tcols))
    elseif variable_name isa AbstractVector
        var_name = variable_name
        @assert length(var_name) == length(tcols)
    else
        throw(ArgumentError("`variable_name` must be a string, symbol, nothing, or a vector of them"))
    end
    if id !== nothing
        ididx = multiple_getindex(index(ds), id)

        if renamecolid === nothing
            renamecolid = _default_renamecolid_function_withid
        end
        ids_refs, unique_loc  = _find_id_unique_values(ds, ididx, _get_perms(ds); mapformats = mapformats)

        if length(ididx) == 1
            unique_ids = getindex(ds, view(_get_perms(ds), unique_loc), ididx[1]; mapformats = mapformats)
        else
            #TODO not very good way to do this
            unique_ids = Tables.rowtable(Dataset([getindex(ds, view(_get_perms(ds), unique_loc), ididx[k], mapformats = mapformats) for k in 1:length(ididx)], :auto, copycols = false))
        end
        @assert (size(unique_ids,1)) == nrow(ds) "Duplicate ids are not allowed."
    end

    local newds
    for j in 1:length(tcols)
        sel_cols = tcols[j]
        ECol = view(_columns(ds), sel_cols)
        T = mapreduce(eltype, promote_type, ECol)
        # make sure the new columns have missing support by default
        T = Union{Missing, T}

        if id === nothing
            if renamecolid === nothing
                renamecolid = _default_renamecolid_function_withoutid
            end
            new_col_names, row_names = _generate_col_row_names(renamecolid, renamerowid, 1:nrow(ds), names(ds)[sel_cols], max_length = max_num_col)
        else

            new_col_names, row_names = _generate_col_row_names(renamecolid, renamerowid, unique_ids, names(ds)[sel_cols], max_length = max_num_col)
        end

        outputmat = _simple_transpose_ds_generate(T, ECol, max_num_col, new_col_names, variable_name, threads)
        if j == 1
            newds = Dataset(outputmat, new_col_names)
            if var_name[j] !== nothing
                new_var_label = Symbol(var_name[j])
                newds = insertcols!(newds, 1,  new_var_label => row_names, unsupported_copy_cols = false)
            end
        else
            ds2 = Dataset(outputmat, new_col_names)
            u = add_names(index(newds), index(ds2), makeunique=true)
            for i in 1:length(u)
                newds[!, u[i]] = ds2[!, i].val
            end
            if var_name[j] !== nothing
                new_var_label = Symbol(var_name[j])
                insertcols!(newds, j, new_var_label => row_names, unsupported_copy_cols = false, makeunique = true)
            end
        end
    end
    newds
end

# groupby case
function _fill_onecol_for_tr!(y, x, ntimes, perms)
    for i in 1:length(perms)
        fill!(view(y, (i-1)*ntimes+1:(i*ntimes)), x[perms[i]])
    end
end
function _fill_onecol_for_tr_threaded!(y, x, ntimes, perms)
    Threads.@threads for i in 1:length(perms)
        fill!(view(y, (i-1)*ntimes+1:(i*ntimes)), x[perms[i]])
    end
end

function _fill_row_names!(res, row_names, ntimes)
    n = length(row_names)
    for i in 1:ntimes
        @views copy!(res[(i-1)*n+1:i*n], row_names)
    end
    res
end

function _fill_gcol!(res, ds, gcolindex, colsidx_length, perms, nrows, threads)
    ntimes = colsidx_length
    totalrow = nrows * ntimes
    for i in 1:length(gcolindex)
        _tmp = allocatecol(ds[!,gcolindex[i]].val, totalrow)
        push!(res, _tmp)
        if DataAPI.refpool(res[i]) !== nothing
            if threads
                _fill_onecol_for_tr_threaded!(res[i].refs, DataAPI.refarray(_columns(ds)[gcolindex[i]]), ntimes, perms)
            else
                _fill_onecol_for_tr!(res[i].refs, DataAPI.refarray(_columns(ds)[gcolindex[i]]), ntimes, perms)
            end
        else
            if threads
                _fill_onecol_for_tr_threaded!(res[i], _columns(ds)[gcolindex[i]], ntimes, perms)
            else
                _fill_onecol_for_tr!(res[i], _columns(ds)[gcolindex[i]], ntimes, perms)
            end
        end
    end
    res
end

function _fill_col_val_f_barrier_threaded!(res, xvals, perms, ntimes, max_num_col, ds_n_row, j)
    Threads.@threads for i in 1:ds_n_row
        res[(i-1)*max_num_col+j] = xvals[perms[i]]
    end
end
function _fill_col_val_f_barrier!(res, xvals, perms, ntimes, max_num_col, ds_n_row, j)
    for i in 1:ds_n_row
        res[(i-1)*max_num_col+j] = xvals[perms[i]]
    end
end

function _fill_col_val!(res, in_cols, ntimes, max_num_col, ds_n_row, perms, threads)
    for j in 1:ntimes
        if threads
            _fill_col_val_f_barrier_threaded!(res, in_cols[j], perms, ntimes, max_num_col, ds_n_row, j)
        else
            _fill_col_val_f_barrier!(res, in_cols[j], perms, ntimes, max_num_col, ds_n_row, j)
        end
    end
end

function fast_stack_gcols(T, ds, in_cols, colsidx_length, gcolsidx, threads)
    # construct group columns
    g_array = AbstractArray[]
    _fill_gcol!(g_array, parent(ds), gcolsidx, colsidx_length, _get_perms(ds), nrow(ds), threads)
    ds1 = Dataset(g_array, _names(ds)[gcolsidx], copycols = false)
    ds1
end

function _extend_repeat_row_names!(row_names, max_num_col)
    r_n_l = length(row_names)
    if r_n_l < max_num_col
        resize!(row_names, max_num_col)
        for i in r_n_l+1:max_num_col
            row_names[i] = missing
        end
    end
    row_names
end

function _obtain_maximum_groups_size(starts, nrows)
    maxvalue = nrows - starts[end] + 1
    for i in (length(starts) - 1):-1:1
        diffvalue = starts[i + 1] - starts[i]
        if diffvalue > maxvalue
            maxvalue = diffvalue
        end
    end
    return maxvalue
end

function _fill_one_col_transpose!(outputmat, xval, starts, perms, n_row_names, row, ngrps)
    for g in 1:ngrps
        lo = starts[g]
        g == ngrps ? hi = length(xval) : hi = starts[g+1] - 1
        cnt = 1
        for i in lo:hi
            outputmat[cnt][(g - 1) * n_row_names + row] = xval[perms[lo + cnt - 1]]
            cnt += 1
        end
    end
end
function _fill_one_col_transpose_threaded!(outputmat, xval, starts, perms, n_row_names, row, ngrps)
    Threads.@threads for g in 1:ngrps
        lo = starts[g]
        g == ngrps ? hi = length(xval) : hi = starts[g+1] - 1
        cnt = 1
        for i in lo:hi
            outputmat[cnt][(g - 1) * n_row_names + row] = xval[perms[lo + cnt - 1]]
            cnt += 1
        end
    end
end

function _fill_one_col_transpose_id!(outputmat, xval, starts, perms, n_row_names, _is_cell_filled, ids, row, ngrps)
    for g in 1:ngrps
        lo = starts[g]
        g == ngrps ? hi = length(xval) : hi = starts[g+1] - 1
        counter = 1
        for i in lo:hi
            cnt = ids[i]
            _row_ = (g - 1) * n_row_names + row
            if _is_cell_filled[_row_, cnt]
                throw(AssertionError("Duplicate id within a group is not allowed"))
            else
                outputmat[cnt][_row_] = xval[perms[lo + counter - 1]]
                _is_cell_filled[_row_, cnt] = true
                counter += 1
            end
        end
    end
end
function _fill_one_col_transpose_id_threaded!(outputmat, xval, starts, perms, n_row_names, _is_cell_filled, ids, row, ngrps)
    Threads.@threads for g in 1:ngrps
        lo = starts[g]
        g == ngrps ? hi = length(xval) : hi = starts[g+1] - 1
        counter = 1
        for i in lo:hi
            cnt = ids[i]
            _row_ = (g - 1) * n_row_names + row
            if _is_cell_filled[_row_, cnt]
                throw(AssertionError("Duplicate id within a group is not allowed"))
            else
                outputmat[cnt][_row_] = xval[perms[lo + counter - 1]]
                _is_cell_filled[_row_, cnt] = true
                counter += 1
            end
        end
    end
end

function update_outputmat!(outputmat, x, starts, perms, n_row_names, threads)
    if threads
        for j in 1:length(x)
            _fill_one_col_transpose_threaded!(outputmat, x[j], starts, perms, n_row_names, j, length(starts))
        end
    else
        for j in 1:length(x)
            _fill_one_col_transpose!(outputmat, x[j], starts, perms, n_row_names, j, length(starts))
        end
    end
end

function update_outputmat!(outputmat, x, starts, perms, ids, n_row_names, _is_cell_filled, threads)
    if threads
        for j in 1:length(x)
            _fill_one_col_transpose_id_threaded!(outputmat, x[j], starts, perms, n_row_names, _is_cell_filled, ids, j, length(starts))
        end
    else
        for j in 1:length(x)
            _fill_one_col_transpose_id!(outputmat, x[j], starts, perms, n_row_names, _is_cell_filled, ids, j, length(starts))
        end
    end
end

function _preallocate_outputmat!(res, n1, n2, fillval, threads, ::Val{CT}) where CT
    @_threadsfor threads for j in 1:n2
        res[j] = fill!(_our_vect_alloc(CT, n1), fillval)
    end
    res
end

function _fill_outputmat_withoutid(T, in_cols, ds, starts, perms, new_col_names, row_names_length, threads; default_fill = missing)

    @assert _check_allocation_limit(nonmissingtype(T), row_names_length*_ngroups(ds), length(new_col_names)) < 1.0 "The output data frame is huge and there is not enough resource to allocate it."
    CT = promote_type(T, typeof(default_fill))
    # outputmat = [__fill!(_our_vect_alloc(CT, row_names_length*_ngroups(ds)), default_fill) for _ in 1:length(new_col_names)]
    outputmat = Vector{typeof(_our_vect_alloc(CT, 0))}(undef, length(new_col_names))
    _preallocate_outputmat!(outputmat, row_names_length*_ngroups(ds), length(new_col_names), default_fill, threads, Val(CT) )
    update_outputmat!(outputmat, in_cols, starts, perms, row_names_length, threads)

    outputmat
end

function _fill_outputmat_withid(T, in_cols, ds, starts, perms, ids, new_col_names, row_names_length, threads; default_fill = missing)

    @assert _check_allocation_limit(nonmissingtype(T), row_names_length*_ngroups(ds), length(new_col_names)) < 1.0 "The output data frame is huge and there is not enough resource to allocate it."
    CT = promote_type(T, typeof(default_fill))
    # outputmat = [fill!(_our_vect_alloc(CT, row_names_length*_ngroups(ds)), default_fill) for _ in 1:length(new_col_names)]
    outputmat = Vector{typeof(_our_vect_alloc(CT, 0))}(undef, length(new_col_names))
    _preallocate_outputmat!(outputmat, row_names_length*_ngroups(ds), length(new_col_names), default_fill, threads, Val(CT) )

    _is_cell_filled = zeros(Bool, row_names_length*_ngroups(ds), length(new_col_names))

    update_outputmat!(outputmat, in_cols, starts, perms, ids, row_names_length, _is_cell_filled, threads)

    outputmat
end

function ds_transpose(ds::Union{Dataset, GroupBy, GatherBy}, cols::Union{Tuple, MultiColumnIndex}, gcols::MultiColumnIndex; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = "_variables_", default_fill = missing, threads = true, mapformats = true)
    if cols isa Tuple
        tcols = [cols[j] isa ColumnIndex ? index(ds)[[cols[j]]] : multiple_getindex(index(ds), cols[j]) for j in 1:length(cols)]
    else
        tcols = [multiple_getindex(index(ds), cols)]
    end
    max_num_col = maximum(length, tcols)
    gcolsidx = gcols
    if variable_name isa AbstractString || variable_name isa Symbol || variable_name === nothing
        var_name = repeat([variable_name], length(tcols))
    elseif variable_name isa AbstractVector
        var_name = variable_name
        @assert length(var_name) == length(tcols)
    else
        throw(ArgumentError("`variable_name` must be a string, symbol, nothing, or a vector of them"))
    end
    local outds

    need_fast_stack = false
    if _ngroups(ds) == nrow(ds)
        need_fast_stack = true
    end
    if need_fast_stack && id === nothing
        if renamecolid === nothing
            renamecolid = _default_renamecolid_function_withoutid
        end
        # fast_stack path, while keeping the row order consistent
        for j in 1:length(tcols)
            sel_cols = tcols[j]
            ECol = view(_columns(ds), sel_cols)

            T = mapreduce(eltype, promote_type, ECol)
            T = Union{T, Missing}

            if j == 1
                outds = fast_stack_gcols(T, ds, ECol, max_num_col, gcolsidx, threads)
                for j in 1:length(gcolsidx)
                    setformat!(outds, j => getformat(parent(ds), gcolsidx[j]))
                end
            end
            if var_name[j] !== nothing
                _repeat_row_names = allowmissing(PooledArray(renamerowid.(names(ds)[sel_cols])))
                _extend_repeat_row_names!(_repeat_row_names, max_num_col)
                _repeat_row_names.refs = repeat(_repeat_row_names.refs, nrow(ds))
                new_var_label = Symbol(var_name[j])
                insertcols!(outds, length(gcolsidx)+j, new_var_label => _repeat_row_names, unsupported_copy_cols = false, makeunique = true)
            end
            res = _our_vect_alloc(Union{Missing, T}, nrow(ds) * max_num_col)
            _fill_col_val!(res, ECol, length(sel_cols), max_num_col, nrow(ds), _get_perms(ds), threads)
            local new_col_id
            try
                new_col_id = Symbol(renamecolid(1, names(ds)[sel_cols]))
            catch e
                if (e isa MethodError)
                    new_col_id = Symbol(renamecolid(1))
                else
                    rethrow(e)
                end
            end
            insertcols!(outds, ncol(outds)+1, new_col_id => res, unsupported_copy_cols = false, makeunique = true)
        end
        return outds

    end
    if id !== nothing
        ididx = multiple_getindex(index(ds), id)
        if renamecolid === nothing
            renamecolid = _default_renamecolid_function_withid
        end
        ids_refs, unique_loc  = _find_id_unique_values(parent(ds), ididx, _get_perms(ds); mapformats = mapformats)

        # we assume the unique function keep the same order as original data, which is the case sofar
        if length(ididx) == 1
            unique_ids = getindex(parent(ds), view(_get_perms(ds), unique_loc), ididx[1], mapformats = mapformats)
        else
            #TODO not very good way to do this
            unique_ids = Tables.rowtable(Dataset([getindex(parent(ds), view(_get_perms(ds), unique_loc), ididx[k], mapformats = mapformats) for k in 1:length(ididx)], :auto, copycols = false))
        end

        out_ncol = length(unique_ids)
    end

    for j in 1:length(tcols)
        sel_cols = tcols[j]
        ECol = view(_columns(ds), sel_cols)

        T = mapreduce(eltype, promote_type, ECol)
        T = Union{T, Missing}

        if id === nothing
            if renamecolid === nothing
                renamecolid = _default_renamecolid_function_withoutid
            end

            out_ncol = _obtain_maximum_groups_size(view(_group_starts(ds), 1:_ngroups(ds)), nrow(ds))
            new_col_names, row_names = _generate_col_row_names(renamecolid, renamerowid, 1:out_ncol, names(ds)[sel_cols], max_length = max_num_col)
            outputmat = _fill_outputmat_withoutid(T, ECol, ds, view(_group_starts(ds), 1:_ngroups(ds)), _get_perms(ds), new_col_names, max_num_col, threads; default_fill = default_fill)
        else

            new_col_names, row_names = _generate_col_row_names(renamecolid, renamerowid, unique_ids, names(ds)[sel_cols], max_length = max_num_col)
            outputmat = _fill_outputmat_withid(T, ECol, ds, view(_group_starts(ds), 1:_ngroups(ds)),  _get_perms(ds), ids_refs, new_col_names, max_num_col, threads; default_fill = default_fill)
        end
        # rows_with_group_info = _find_group_row(gds)
        new_var_label = Symbol(var_name[j])
        if j == 1
            g_array = AbstractArray[]
            _fill_gcol!(g_array, parent(ds), gcolsidx, max_num_col, view(_get_perms(ds), view(_group_starts(ds), 1:_ngroups(ds))), _ngroups(ds), threads)
            outds = Dataset(g_array, _names(ds)[gcolsidx], copycols = false)
            for j in 1:length(gcolsidx)
                setformat!(outds, j => getformat(parent(ds), gcolsidx[j]))
            end
        end

        # _repeat_row_names = Vector{eltype(row_names)}(undef, _ngroups(ds)*length(colsidx))
        # _fill_row_names!(_repeat_row_names, row_names, _ngroups(ds))
        if var_name[j] !== nothing
            _repeat_row_names = allowmissing(PooledArray(row_names))
            _repeat_row_names.refs = repeat(_repeat_row_names.refs, _ngroups(ds))
            insertcols!(outds, length(gcolsidx)+j, new_var_label => _repeat_row_names, unsupported_copy_cols = false, makeunique = true)
        end
        outds2 = Dataset(outputmat, new_col_names, copycols = false)

         for j in 1:ncol(outds2)
            push!(_columns(outds), _columns(outds2)[j])
        end
        merge!(index(outds), index(outds2), makeunique = true)
    end
    outds
end

function Base.transpose(ds::AbstractDataset, cols::MultiColumnIndex; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = "_variables_", default = missing, threads = true, mapformats = true)
    if !isgrouped(ds)
        ds_transpose(ds, cols; id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, threads = threads, mapformats = mapformats)
    else
        ds_transpose(ds, cols, _groupcols(ds); id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, threads = threads, default_fill = default, mapformats = mapformats)
    end
end

Base.transpose(ds::Union{GroupBy, GatherBy}, cols::MultiColumnIndex; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = "_variables_", default = missing, threads = true, mapformats = true) =
    ds_transpose(ds, cols, _groupcols(ds); id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, threads = threads, default_fill = default, mapformats = mapformats)

Base.transpose(ds::Union{GatherBy, GroupBy, AbstractDataset}, col::ColumnIndex; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = "_variables_", default = missing, threads = true, mapformats = true) =
    transpose(ds, [col]; id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, default = default, threads = threads, mapformats = mapformats)

function Base.transpose(ds::AbstractDataset, cols::Tuple; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = nothing, default = missing, threads = true, mapformats = true)
    if !isgrouped(ds)
        ds_transpose(ds, cols; id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, threads = threads, mapformats = mapformats)
    else
        ds_transpose(ds, cols, _groupcols(ds); id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, threads = threads, default_fill = default, mapformats = mapformats)
    end
end

Base.transpose(ds::Union{GroupBy, GatherBy}, cols::Tuple; id = nothing, renamecolid = nothing, renamerowid = _default_renamerowid_function, variable_name = nothing, default = missing, threads = true, mapformats = true) =
    ds_transpose(ds, cols, _groupcols(ds); id = id, renamecolid = renamecolid, renamerowid = renamerowid, variable_name = variable_name, threads = threads, default_fill = default, mapformats = mapformats)

#### flatten

"""
    flatten(ds::AbstractDataset, cols; mapformats = false, threads = true)

When columns `cols` of data set `ds` have iterable elements that define
`length` (for example a `Vector` of `Vector`s), return a `Dataset` where each
element of each `col` in `cols` is flattened, meaning the column corresponding
to `col` becomes a longer vector where the original entries are concatenated.
Elements of row `i` of `ds` in columns other than `cols` will be repeated
according to the length of `ds[i, col]`. These lengths must therefore be the
same for each `col` in `cols`, or else an error is raised. Note that these
elements are not copied, and thus if they are mutable changing them in the
returned `Dataset` will affect `ds`.

When `mapformats = true`, the function uses the formatted values of `cols`.

`cols` can be any column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).

To turn off multithreaded computations pass `threads = false`.

See [`flatten!`](@ref)

# Examples

```jldoctest
julia> ds1 = Dataset(a = [1, 2], b = [[1, 2], [3, 4]], c = [[5, 6], [7, 8]])
2×3 Dataset
 Row │ a         b         c
     │ identity  identity  identity
     │ Int64?    Array…?   Array…?
─────┼──────────────────────────────
   1 │        1  [1, 2]    [5, 6]
   2 │        2  [3, 4]    [7, 8]

julia> flatten(ds1, :b)
4×3 Dataset
 Row │ a         b         c
     │ identity  identity  identity
     │ Int64?    Int64?    Array…?
─────┼──────────────────────────────
   1 │        1         1  [5, 6]
   2 │        1         2  [5, 6]
   3 │        2         3  [7, 8]
   4 │        2         4  [7, 8]

julia> flatten(ds1, [:b, :c])
4×3 Dataset
 Row │ a         b         c
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         1         5
   2 │        1         2         6
   3 │        2         3         7
   4 │        2         4         8

julia> ds2 = Dataset(a = [1, 2], b = [("p", "q"), ("r", "s")])
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Tuple…?
─────┼──────────────────────
   1 │        1  ("p", "q")
   2 │        2  ("r", "s")

julia> flatten(ds2, :b)
4×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    String?
─────┼────────────────────
   1 │        1  p
   2 │        1  q
   3 │        2  r
   4 │        2  s

julia> ds3 = Dataset(a = [1, 2], b = [[1, 2], [3, 4]], c = [[5, 6], [7]])
2×3 Dataset
 Row │ a         b         c
     │ identity  identity  identity
     │ Int64?    Array…?   Array…?
─────┼──────────────────────────────
   1 │        1  [1, 2]    [5, 6]
   2 │        2  [3, 4]    [7]

julia> flatten(ds3, [:b, :c])
ERROR: ArgumentError: Lengths of iterables stored in columns :b and :c are not the same in row 2

julia> ds = Dataset(x=1:3, y=["ab,cd", "e", missing], z=[[1,2], 2, 3])
3×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    String?   Any
─────┼──────────────────────────────
   1 │        1  ab,cd     [1, 2]
   2 │        2  e         2
   3 │        3  missing   3

julia> fmt(x) = split(x, ",")
fmt (generic function with 2 methods)

julia> fmt(::Missing) = missing
fmt (generic function with 2 methods)

julia> setformat!(ds, :y => fmt)
3×3 Dataset
 Row │ x         y                              z
     │ identity  fmt                            identity
     │ Int64?    String?                        Any
─────┼───────────────────────────────────────────────────
   1 │        1  SubString{String}["ab", "cd"]  [1, 2]
   2 │        2  SubString{String}["e"]         2
   3 │        3  missing                        3

julia> flatten(ds, :y)
7×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    Char?     Any
─────┼──────────────────────────────
   1 │        1  a         [1, 2]
   2 │        1  b         [1, 2]
   3 │        1  ,         [1, 2]
   4 │        1  c         [1, 2]
   5 │        1  d         [1, 2]
   6 │        2  e         2
   7 │        3  missing   3

julia> flatten(ds, :y, mapformats = true)
4×3 Dataset
 Row │ x         y           z
     │ identity  identity    identity
     │ Int64?    SubStrin…?  Any
─────┼────────────────────────────────
   1 │        1  ab          [1, 2]
   2 │        1  cd          [1, 2]
   3 │        2  e           2
   4 │        3  missing     3

julia> flatten(ds, 2:3, mapformats = true)
4×3 Dataset
 Row │ x         y           z
     │ identity  identity    identity
     │ Int64?    SubStrin…?  Int64?
─────┼────────────────────────────────
   1 │        1  ab                 1
   2 │        1  cd                 2
   3 │        2  e                  2
   4 │        3  missing            3
```
"""
flatten(ds, cols)

"""
    flatten!(ds, cols; mapformats = false, threads = true)

Variant of `flatten` that does flatten `ds` in-place.
"""
flatten!

function _ELTYPE(x)
    eltype(x)
end
function _ELTYPE(x::Missing)
    Missing
end

function _LENGTH(x)::Int
    length(x)
end

function _LENGTH(x::Missing)::Int
    1
end

function flatten!(ds::Dataset,
                 cols::Union{ColumnIndex, MultiColumnIndex}; mapformats = false, threads = true)
     _check_consistency(ds)

     idxcols = index(ds)[cols]
     isempty(idxcols) && return ds
     col1 = first(idxcols)
     all_idxcols = Any[]
     if mapformats
         f_fmt = getformat(ds, col1)
         push!(all_idxcols, byrow(ds, f_fmt, col1, threads = threads))
     else
        push!(all_idxcols, _columns(ds)[col1])
    end
    lengths = byrow(Dataset(all_idxcols, [:x], copycols = false), _LENGTH, 1, threads = threads, allowmissing = false)
    if length(idxcols) > 1
        for col in 2:length(idxcols)
             if mapformats
                 f_fmt = getformat(ds, idxcols[col])
                 push!(all_idxcols, byrow(ds, f_fmt, idxcols[col]), threads = threads)
             else
                 push!(all_idxcols, _columns(ds)[idxcols[col]])
             end
             v = all_idxcols[col]
             if any(x -> _LENGTH(x[1]) != x[2], zip(v, lengths))
                 r = findfirst(x -> x != 0, _LENGTH.(v) .- lengths)
                 colnames = _names(ds)
                 throw(ArgumentError("Lengths of iterables stored in columns :$(colnames[col1]) " *
                                     "and :$(colnames[idxcols[col]]) are not the same in row $r"))
             end
         end
     end
     r_index = _create_index_for_repeat(lengths, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64))
     _permute_ds_after_sort!(ds, r_index, check = false, cols = Not(cols), threads = threads)
     if threads
         new_total = hp_sum(lengths)
     else
         new_total = sum(lengths)
     end
     if length(idxcols) > 1
         sort_permute_idxcols = sortperm(idxcols)
         idxcols_sorted = idxcols[sort_permute_idxcols]
     else
         sort_permute_idxcols = [1]
         idxcols_sorted = idxcols
     end
     cumsum!(lengths, lengths)
     for col in 1:length(idxcols_sorted)
         col_to_flatten = all_idxcols[sort_permute_idxcols[col]]

         T = mapreduce(_ELTYPE, promote_type, col_to_flatten)
         _res = allocatecol(T, new_total)
         _fill_flatten!(_res, col_to_flatten, lengths, threads = threads)
         if length(idxcols) == ncol(ds)
             _columns(ds)[idxcols_sorted[col]] = _res
         else
             deleteat!(_columns(ds), idxcols_sorted[col])
             insert!(_columns(ds), idxcols_sorted[col], _res)
         end
         removeformat!(ds, idxcols_sorted[col])
     end
     _reset_grouping_info!(ds)
     _modified(_attributes(ds))
     ds
end

function flatten(ds::AbstractDataset,
                 cols::Union{ColumnIndex, MultiColumnIndex}; mapformats = false, threads = true)
     _check_consistency(ds)

     idxcols = index(ds)[cols]
     isempty(idxcols) && return copy(ds)
     col1 = first(idxcols)
     all_idxcols = Any[]
     if mapformats
         f_fmt = getformat(ds, col1)
         push!(all_idxcols, byrow(ds, f_fmt, col1, threads = threads))
     else
        push!(all_idxcols, _columns(ds)[col1])
    end
    lengths = byrow(Dataset(all_idxcols, [:x], copycols = false), _LENGTH, 1, threads = threads, allowmissing = false)
    if length(idxcols) > 1
        for col in 2:length(idxcols)
             if mapformats
                 f_fmt = getformat(ds, idxcols[col])
                 push!(all_idxcols, byrow(ds, f_fmt, idxcols[col]), threads = threads)
             else
                 push!(all_idxcols, _columns(ds)[idxcols[col]])
             end
             v = all_idxcols[col]
             if any(x -> _LENGTH(x[1]) != x[2], zip(v, lengths))
                 r = findfirst(x -> x != 0, _LENGTH.(v) .- lengths)
                 colnames = _names(ds)
                 throw(ArgumentError("Lengths of iterables stored in columns :$(colnames[col1]) " *
                                     "and :$(colnames[idxcols[col]]) are not the same in row $r"))
             end
         end
     end
     if threads
         new_total = hp_sum(lengths)
     else
         new_total = sum(lengths)
     end
     new_ds = similar(ds[!, Not(cols)], new_total)
     for name in _names(new_ds)
         col_name = index(ds)[name]
        repeat_lengths_v2!(new_ds[!, name].val, _columns(ds)[col_name], lengths)
     end
     if length(idxcols) > 1
         sort_permute_idxcols = sortperm(idxcols)
         idxcols_sorted = idxcols[sort_permute_idxcols]
     else
         sort_permute_idxcols = [1]
         idxcols_sorted = idxcols
     end
     cumsum!(lengths, lengths)
     for col in 1:length(idxcols_sorted)
         col_to_flatten = all_idxcols[sort_permute_idxcols[col]]

         T = mapreduce(_ELTYPE, promote_type, col_to_flatten)
         _res = allocatecol(T, new_total)
         _fill_flatten!(_res, col_to_flatten, lengths, threads = threads)
         insertcols!(new_ds, idxcols_sorted[col], _names(ds)[idxcols_sorted[col]] => _res, unsupported_copy_cols = false)
     end
     for j in setdiff(1:ncol(ds), idxcols)
         setformat!(new_ds, j=>getformat(ds, j))
     end
     setinfo!(new_ds, _attributes(ds).meta.info[])
     _reset_grouping_info!(new_ds)
     new_ds
end

function _fill_flatten!_barrier(_res, val, lo)
    if ismissing(val)
        _res[lo] = missing
    else

        cnt = 0
        for j in val
            _res[lo+cnt] = j
            cnt += 1
        end
    end
end

function _fill_flatten!(_res, col_to_flatten, lengths; threads = false)
    @_threadsfor threads for i in 1:length(col_to_flatten)
        i == 1 ? lo = 1 : lo = lengths[i-1]+1
        _fill_flatten!_barrier(_res, col_to_flatten[i], lo)
    end
end


text:
module STL

struct Vertex
    x::Float32
    y::Float32
    z::Float32
end

const Normal = Vertex

struct Triangle
    normal::Normal
    v1::Vertex
    v2::Vertex
    v3::Vertex
end

function parse(path::AbstractString)
    open(path) do stl
        skip(stl, 80)  # skip header

        trianglecount = read(stl, UInt32)
        triangles = Array{Triangle, 1}(undef, trianglecount)

        for i in 1:trianglecount
            normal = Normal(read(stl, Float32), read(stl, Float32), read(stl, Float32))
            v1 = Vertex(read(stl, Float32), read(stl, Float32), read(stl, Float32))
            v2 = Vertex(read(stl, Float32), read(stl, Float32), read(stl, Float32))
            v3 = Vertex(read(stl, Float32), read(stl, Float32), read(stl, Float32))
            triangles[i] = Triangle(normal, v1, v2, v3)

            skip(stl, 2)  # skip attribute byte count
        end

        @assert eof(stl)
        return triangles
    end
end

export parse

end # module

ext: jl

max_stars_repo_name: SimonDanisch/stl-benchmark

max_issues_repo_name: SimonDanisch/stl-benchmark

max_forks_repo_name: SimonDanisch/stl-benchmark


text:
# Extract data from a 2D/3D DG solution and prepare it for visualization as a heatmap/contour plot.
#
# Returns a tuple with
# - x coordinates
# - y coordinates
# - nodal 2D data
# - x vertices for mesh visualization
# - y vertices for mesh visualization
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_data_2d(center_level_0, length_level_0, leaf_cells, coordinates, levels, ndims,
                     unstructured_data, n_nodes,
                     grid_lines=false, max_supported_level=11, nvisnodes=nothing,
                     slice_axis=:z, slice_axis_intercept=0)
  # Determine resolution for data interpolation
  max_level = maximum(levels)
  if max_level > max_supported_level
    error("Maximum refinement level $max_level is higher than " *
          "maximum supported level $max_supported_level")
  end
  max_available_nodes_per_finest_element = 2^(max_supported_level - max_level)
  if nvisnodes === nothing
    max_nvisnodes = 2 * n_nodes
  elseif nvisnodes == 0
    max_nvisnodes = n_nodes
  else
    max_nvisnodes = nvisnodes
  end
  nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes)
  resolution = nvisnodes_at_max_level * 2^max_level
  nvisnodes_per_level = [2^(max_level - level)*nvisnodes_at_max_level for level in 0:max_level]
  # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate
  # level-0-cell) that contains the number of visualization nodes for any
  # refinement level to visualize on an equidistant grid

  if ndims == 3
    (unstructured_data, coordinates, levels,
        center_level_0) = unstructured_2d_to_3d(unstructured_data,
        coordinates, levels, length_level_0, center_level_0, slice_axis,
        slice_axis_intercept)
  end

  # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]²
  n_elements = length(levels)
  normalized_coordinates = similar(coordinates)
  for element_id in 1:n_elements
    @views normalized_coordinates[:, element_id] .= (
          (coordinates[:, element_id] .- center_level_0) ./ (length_level_0 / 2 ))
  end

  # Interpolate unstructured DG data to structured data
  (structured_data =
      unstructured2structured(unstructured_data, normalized_coordinates,
                              levels, resolution, nvisnodes_per_level))

  # Interpolate cell-centered values to node-centered values
  node_centered_data = cell2node(structured_data)

  # Determine axis coordinates for contour plot
  xs = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[1]
  ys = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[2]

  # Determine element vertices to plot grid lines
  if grid_lines
    mesh_vertices_x, mesh_vertices_y = calc_vertices(coordinates, levels, length_level_0)
  else
    mesh_vertices_x = Matrix{Float64}(undef, 0, 0)
    mesh_vertices_y = Matrix{Float64}(undef, 0, 0)
  end

  return xs, ys, node_centered_data, mesh_vertices_x, mesh_vertices_y
end

# Change order of dimensions (variables are now last) and convert data to `solution_variables`
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_unstructured_data(u, semi, solution_variables)
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)

  if solution_variables === cons2cons
    raw_data = u
    n_vars = size(raw_data, 1)
  else
    # Reinterpret the solution array as an array of conservative variables,
    # compute the solution variables via broadcasting, and reinterpret the
    # result as a plain array of floating point numbers
    raw_data = Array(reinterpret(eltype(u),
           solution_variables.(reinterpret(SVector{nvariables(equations),eltype(u)}, u),
                      Ref(equations))))
    n_vars = size(raw_data, 1)
  end

  unstructured_data = Array{Float64}(undef,
                                     ntuple((d) -> nnodes(solver), ndims(equations))...,
                                     nelements(solver, cache), n_vars)
  for variable in 1:n_vars
    @views unstructured_data[.., :, variable] .= raw_data[variable, .., :]
  end

  return unstructured_data
end


text:
# Override mul! method.
#

"""
    mul!(C, A, B, α, β)
BLIS-based matrix multiply.
Generic strides are directly supported.
"""
mul!(C::StridedMatrix{TC},
     A::StridedMatrix{TA},
     B::StridedMatrix{TB},
     α::Tα,
     β::Tβ) where {TA<:BliCompatibleType,
                   TB<:BliCompatibleType,
                   TC<:BliCompatibleType,
                   Tα<:BliCompatibleType,
                   Tβ<:BliCompatibleType} = begin

    oA = BliObj(A)
    oB = BliObj(B)
    oC = BliObj(C)
    oα = BliObj(α)
    oβ = BliObj(β)

    ObjectBackend.bli_gemm!(oα, oA, oB, oβ, oC)
    C
end

"""
    mul!(C, A, B)
BLIS-based matrix multiply.
Generic strides are directly supported.
"""
mul!(C::StridedMatrix{TC},
     A::StridedMatrix{TA},
     B::StridedMatrix{TB}
    ) where {TA<:BliCompatibleType,
             TB<:BliCompatibleType,
             TC<:BliCompatibleType} = mul!(C, A, B, TA(1.0), TC(0.0))

#--------------------------------------------------
# Instantiate mul! separately to avoid conflicting
#  with Julia's intrinsic.

mul!(C::StridedMatrix{TC},
     A::StridedMatrix{TA},
     B::Adjoint{TB, <:StridedMatrix{TB}},
     α::Tα,
     β::Tβ) where {TA<:BliCompatibleType,
                   TB<:BliCompatibleType,
                   TC<:BliCompatibleType,
                   Tα<:BliCompatibleType,
                   Tβ<:BliCompatibleType} = begin

    oA = BliObj(A)
    oB = BliObj(B')
    oC = BliObj(C)
    oα = BliObj(α)
    oβ = BliObj(β)

    ObjectBackend.bli_obj_set_conjtrans!(BLIS_CONJ_TRANSPOSE, oB.obj)
    ObjectBackend.bli_gemm!(oα, oA, oB, oβ, oC)
    C
end

# mul!(C::StridedMatrix{TC},
#      A::Adjoint{TA, <:StridedMatrix{TA}},
#      B::StridedMatrix{TB},
#      α::Tα,
#      β::Tβ) where {TA<:BliCompatibleType,
#                    TB<:BliCompatibleType,
#                    TC<:BliCompatibleType,
#                    Tα<:BliCompatibleType,
#                    Tβ<:BliCompatibleType} = begin
#
#     oA = BliObj(A')
#     oB = BliObj(B)
#     oC = BliObj(C)
#     oα = BliObj(α)
#     oβ = BliObj(β)
#     @show oα
#
#     ObjectBackend.bli_obj_set_conjtrans!(BLIS_CONJ_TRANSPOSE, oA.obj)
#     ObjectBackend.bli_gemm!(oα, oA, oB, oβ, oC)
#     C
# end

mul!(C::StridedMatrix{TC},
     A::Adjoint{TA, <:StridedMatrix{TA}},
     B::Adjoint{TB, <:StridedMatrix{TB}},
     α::Tα,
     β::Tβ) where {TA<:BliCompatibleType,
                   TB<:BliCompatibleType,
                   TC<:BliCompatibleType,
                   Tα<:BliCompatibleType,
                   Tβ<:BliCompatibleType} = begin

    oA = BliObj(A')
    oB = BliObj(B')
    oC = BliObj(C)
    oα = BliObj(α)
    oβ = BliObj(β)

    ObjectBackend.bli_obj_set_conjtrans!(BLIS_CONJ_TRANSPOSE, oA.obj)
    ObjectBackend.bli_obj_set_conjtrans!(BLIS_CONJ_TRANSPOSE, oB.obj)
    ObjectBackend.bli_gemm!(oα, oA, oB, oβ, oC)
    C
end

mul!(C::StridedMatrix{TC},
     A::Adjoint{TA, <:StridedMatrix{TA}},
     B::Adjoint{TB, <:StridedMatrix{TB}}
    ) where {TA<:BliCompatibleType,
             TB<:BliCompatibleType,
             TC<:BliCompatibleType} = mul!(C, A, B, TA(1.0), TC(0.0))

mul!(C::StridedMatrix{TC},
     A::StridedMatrix{TA},
     B::Adjoint{TB, <:StridedMatrix{TB}}
    ) where {TA<:BliCompatibleType,
             TB<:BliCompatibleType,
             TC<:BliCompatibleType} = mul!(C, A, B, TA(1.0), TC(0.0))

# A'B instantiated here as mul!(C, A', B, α, β)
#  has additional ambiguity.
mul!(C::StridedMatrix{TC},
     A::Adjoint{TA, <:StridedMatrix{TA}},
     B::StridedMatrix{TB}
    ) where {TA<:BliCompatibleType,
             TB<:BliCompatibleType,
             TC<:BliCompatibleType} = begin
    oA = BliObj(A')
    oB = BliObj(B)
    oC = BliObj(C)
    oα = BliObj(1.0)
    oβ = BliObj(0.0)

    ObjectBackend.bli_obj_set_conjtrans!(BLIS_CONJ_TRANSPOSE, oA.obj)
    ObjectBackend.bli_gemm!(oα, oA, oB, oβ, oC)
    C
end


text:
"""
    Ground(;name)

Zero magnetic potential.
"""
function Ground(;name)
    @named port = PositiveMagneticPort()
    eqs = [port.V_m ~ 0]
    ODESystem(eqs, t, [], [], systems=[port], name=name)
end

"""
    Idle(;name)

Idle running branch.
"""
function Idle(;name)
    @named two_port = TwoPort()
    @unpack Phi = two_port
    eqs = [
        Phi ~ 0,
    ]
    extend(ODESystem(eqs, t, [], [], systems=[], name=name), two_port)
end

"""
    Short(;name)

Short cut branch.
"""
function Short(;name)
    @named two_port = TwoPort()
    @unpack V_m = two_port
    eqs = [
        V_m ~ 0,
    ]
    extend(ODESystem(eqs, t, [], [], systems=[], name=name), two_port)
end

"""
    Crossing(;name)

Crossing of two branches.

This is a simple crossing of two branches. The ports port_p1 and port_p2 are connected, as well as port_n1 and port_n2.
"""
function Crossing(;name)
    @named port_p1 = PositiveMagneticPort()
    @named port_p2 = PositiveMagneticPort()
    @named port_n1 = NegativeMagneticPort()
    @named port_n2 = NegativeMagneticPort()
    eqs = [
        connect(port_p1, port_p2),
        connect(port_n1, port_n2),
    ]
    ODESystem(eqs, t, [], [], systems=[port_p1, port_p2, port_n1, port_n2], name=name)
end

"""
    ConstantPermeance(;name, G_m=1.0)

Constant permeance.

# Parameters:
- `G_m`: [H] Magnetic permeance
"""
function ConstantPermeance(;name, G_m=1.0)
    @named two_port = TwoPort()
    @unpack V_m, Phi = two_port
    @parameters G_m=G_m
    eqs = [
        Phi ~ G_m * V_m,
    ]
    extend(ODESystem(eqs, t, [], [G_m], name=name), two_port)
end

"""
    ConstantReluctance(;name, R_m=1.0)

Constant reluctance.

# Parameters:
- `R_m`: [H^-1] Magnetic reluctance
"""
function ConstantReluctance(;name, R_m=1.0)
    @named two_port = TwoPort()
    @unpack V_m, Phi = two_port
    @parameters R_m=R_m
    eqs = [
        V_m ~ Phi * R_m,
    ]
    extend(ODESystem(eqs, t, [], [R_m], name=name), two_port)
end

"""
    ElectroMagneticConverter(;name, N, Phi_start=0.0)

Ideal electromagnetic energy conversion.

The electromagnetic energy conversion is given by Ampere's law and Faraday's law respectively
V_m = N * i
N * dΦ/dt = -v

# Parameters:
- `N`: Number of turns
- `Phi_start`: [Wb] Initial magnetic flux flowing into the port_p
"""
function ElectroMagneticConverter(;name, N, Phi_start=0.0)
    @named port_p = PositiveMagneticPort()
    @named port_n = NegativeMagneticPort()
    @named p = Pin()
    @named n = Pin()

    sts = @variables v(t) i(t) V_m(t) Phi(t)=Phi_start
    pars = @parameters N=N
    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
        V_m ~ port_p.V_m - port_n.V_m
        0 ~ port_p.Phi + port_n.Phi
        Phi ~ port_p.Phi
        #converter equations:
        V_m ~ i * N # Ampere's law
        D(Phi) ~ -v / N # Faraday's law
    ]
    ODESystem(eqs, t, sts, pars, systems=[port_p, port_n, p, n], name=name)
end

"""
    EddyCurrent(;name, rho=0.098e-6, l=1, A=1, Phi_start=0.0)

For modelling of eddy current in a conductive magnetic flux tube.

# Parameters:
- `rho`: [ohm * m] Resistivity of flux tube material (default: Iron at 20degC)
- `l`: [m] Average length of eddy current path
- `A`: [m^2] Cross sectional area of eddy current path
- `Phi_start`: [Wb] Initial magnetic flux flowing into the port_p
"""
function EddyCurrent(;name, rho=0.098e-6, l=1, A=1, Phi_start=0.0)
    @named two_port = TwoPort(Phi_start=Phi_start)
    @unpack V_m, Phi = two_port
    @parameters R = rho * l / A # Electrical resistance of eddy current path
    eqs = [
        D(Phi) ~ V_m * R,
    ]
    extend(ODESystem(eqs, t, [], [R], name=name), two_port)
end


text:
using StateSpace
using Distributions
using Test
using Statistics
using LinearAlgebra
using Random

function random_cov(n :: Int)
    A = rand(n, n)
    A = A + A'
    A + n * Matrix(I,n,n)
end

include("test_dispatch.jl")
include("test_KF.jl")
include("test_EKF.jl")
include("test_UKF.jl")
include("test_EnKF.jl")
include("test_particle_filter.jl")


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ b58fd1ff-e776-4721-8ea4-dd0d52f7de2e
using Pkg

# ╔═╡ 36d2a4e8-f81f-11eb-3c33-bb29b29df2fb
Pkg.activate("../")# The FTensorKit.jl enviroment

# ╔═╡ ebea242e-b9c0-49b7-a3cb-452b8d515831
using FTensorKit,Viznet,Compose,LinearAlgebra,OMEinsum,LightGraphs,DataStructures

# ╔═╡ 5f9fd5d2-a0f9-4b68-9412-73a0d1e548d1
# Project preliminary
# see: https://arxiv.org/abs/2006.08289

# ╔═╡ 08bcba33-afed-43b8-a5c6-eeeb0fd0cfd5
lattice = Cubic{OBC}(2,3,4); # create a lattice (contains graph and coordinates)

# ╔═╡ 288502fa-a0b3-432a-94e3-ffa35845490a
lattice2D = PerspectiveMap3D(lattice,0.17,0.19); # map this lattice to 2D

# ╔═╡ febaf067-af8c-4628-bf31-8df9fb227b2c
begin
	# set brush
	linebrush = bondstyle(:default, stroke("grey"))
	nodebrush = nodestyle(:default, stroke("lightblue"),fill("lightblue"),r=0.011)
	gatesbrush = nodestyle(:square, stroke("lightgreen"),fill("lightgreen");r=0.005) # gate brush
	textbrush = textstyle(:default)
end;

# ╔═╡ b696ca61-56cc-4bfd-9672-849a60b091a9
canvas() do
	# plot 2D version for the original lattice
	l2D = lattice2D.locations
	g = lattice2D.graph
    for i =1:size(l2D)[1]
        nodebrush >> Tuple(l2D[i,:])
    end
    for edge in edges(g)
        linebrush >> (Tuple(l2D[src(edge),:]),Tuple(l2D[dst(edge),:]))
    end
end |> x-> compose(context(0.0, 0.0, 1.0, 1.0), x) |> SVG("origin.svg")

# ╔═╡ 4588198b-afdd-47a3-9027-ff83f0a898fe
# Contraction EinCode from this graph
# TODO: As there is no information about input tensor, this eincode does not process legs order in each tensor.
eincode(lattice2D.graph)

# ╔═╡ e1e37180-d93c-4fa8-9998-71749d52d52f
# Now we are going to insert Swap Gates for Fermions into this Network

# ╔═╡ 0517aed7-3de5-4199-99a6-2a3febc49094
gates_locations = intersects2D(lattice2D);
# contains gates locations for plot
# A good way to check if the gates are generated correctly.

# ╔═╡ df3c1ea2-a013-4cd8-9736-15fb29c051e4
canvas() do
	# plot 2D version for the original lattice
	l2D = lattice2D.locations
	g = lattice2D.graph
    for i =1:size(l2D)[1]
        nodebrush >> Tuple(l2D[i,:])
    end
    for edge in edges(g)
        linebrush >> (Tuple(l2D[src(edge),:]),Tuple(l2D[dst(edge),:]))
    end
	for i in gates_locations
		gatesbrush >> Tuple(i)
	end
end |> x-> compose(context(0.0, 0.0, 1.0, 1.0), x) |> SVG("gates.svg")

# ╔═╡ a453055a-26a3-4a4a-827b-0426387ae7ec
# Transform a normal 2D tensor network into Fermionic Form by adding Fermionic SwapGates.
ftn = FTN(lattice2D);

# ╔═╡ 29edc4e4-bcdf-4bb0-9cb2-207b4d234b18
canvas() do
	# plot 2D version for the Fermionic Tensor Network
	l2D = ftn.locations
	g = ftn.graph
    for i =iter_origin(ftn)
        nodebrush >> Tuple(l2D[i,:])
		textbrush >> (Tuple(l2D[i,:]),"v$i")
    end

	edges = [edge for edge in edges(g)]
    for i = 1:ne(g)
		edge = edges[i]
        linebrush >> (Tuple(l2D[src(edge),:]),Tuple(l2D[dst(edge),:]))
		textbrush >> (Tuple((l2D[src(edge),:]+l2D[dst(edge),:])./2),"b$i")
    end
	for i in iter_gates(ftn)
		gatesbrush >> Tuple(l2D[i,:])
		textbrush >> (Tuple(l2D[i,:]),"$i")
	end
end  |> x-> compose(context(0.0, 0.0, 1.0, 1.0), x) |> SVG("labels.svg")

# ╔═╡ d077e71f-e121-45b4-aa12-d0db624816a3
# Contraction EinCode from this new graph
# TODO: As there is no information about input tensor, this eincode does not process legs order in each tensor.
typeof(eincode(ftn.graph))

# ╔═╡ f5a9c980-d67a-44ab-ba26-928e581550db
exit()

# ╔═╡ Cell order:
# ╠═5f9fd5d2-a0f9-4b68-9412-73a0d1e548d1
# ╠═b58fd1ff-e776-4721-8ea4-dd0d52f7de2e
# ╠═36d2a4e8-f81f-11eb-3c33-bb29b29df2fb
# ╠═ebea242e-b9c0-49b7-a3cb-452b8d515831
# ╠═08bcba33-afed-43b8-a5c6-eeeb0fd0cfd5
# ╠═288502fa-a0b3-432a-94e3-ffa35845490a
# ╠═febaf067-af8c-4628-bf31-8df9fb227b2c
# ╠═b696ca61-56cc-4bfd-9672-849a60b091a9
# ╠═4588198b-afdd-47a3-9027-ff83f0a898fe
# ╠═e1e37180-d93c-4fa8-9998-71749d52d52f
# ╠═0517aed7-3de5-4199-99a6-2a3febc49094
# ╠═df3c1ea2-a013-4cd8-9736-15fb29c051e4
# ╠═a453055a-26a3-4a4a-827b-0426387ae7ec
# ╠═29edc4e4-bcdf-4bb0-9cb2-207b4d234b18
# ╠═d077e71f-e121-45b4-aa12-d0db624816a3
# ╠═f5a9c980-d67a-44ab-ba26-928e581550db

max_stars_repo_name: qiyang-ustc/FTensorKit

max_issues_repo_name: qiyang-ustc/FTensorKit

max_forks_repo_name: qiyang-ustc/FTensorKit


text:
@testset "Basic reading" begin
    recordings_directory = joinpath(@__DIR__, "recordings")
    mkpath(recordings_directory)
    BrokenRecord.configure!(; path = recordings_directory)

    anonymous_auth = FHIRClientXML.AnonymousAuth()
    oauth2_auth = FHIRClientXML.OAuth2()
    FHIRClientXML.set_token!(oauth2_auth, "helloworld")
    jwt_auth = FHIRClientXML.JWTAuth()
    FHIRClientXML.set_token!(jwt_auth, "helloworld")
    username_password_auth_1 = FHIRClientXML.UsernamePassAuth("helloworld")
    FHIRClientXML.set_password!(username_password_auth_1, "helloworld")
    username_password_auth_2 = FHIRClientXML.UsernamePassAuth(; username = "helloworld")
    FHIRClientXML.set_password!(username_password_auth_2, "helloworld")
    auths = [
        anonymous_auth,
        oauth2_auth,
        jwt_auth,
        username_password_auth_1,
        username_password_auth_2,
    ]

    for (i, auth) in enumerate(auths)
        base_url = FHIRClientXML.BaseURL("https://hapi.fhir.org/baseR4")
        client = FHIRClientXML.Client(base_url, auth)
        @test FHIRClientXML.get_base_url(client) == base_url
        search_request_path = "/Patient?given=Jason&family=Argonaut"
        response_search_results_bundle = BrokenRecord.playback(() -> FHIRClientXML.request_xml(client, "GET", search_request_path), "recording-$(i)-1.bson")
        patient_id = response_search_results_bundle["Bundle"]["entry"][1]["resource"]["Patient"]["id"][:value]
        patient_request = "/Patient/$(patient_id)"
        patients = [
            BrokenRecord.playback(() -> FHIRClientXML.request_xml(client, "GET", patient_request), "recording-$(i)-2.bson")
            BrokenRecord.playback(() -> FHIRClientXML.request_xml(client, "GET", patient_request; body = Dict()), "recording-$(i)-3.bson")
            BrokenRecord.playback(() -> FHIRClientXML.request_xml(client, "GET", patient_request; query = Dict{String, String}()), "recording-$(i)-4.bson")
            BrokenRecord.playback(() -> FHIRClientXML.request_xml(client, "GET", patient_request; body = Dict(), query = Dict{String, String}()), "recording-$(i)-5.bson")
        ]
        for patient in patients
            patient_resource = patient["Patient"]
            @test patient isa AbstractDict
            @test patient_resource isa AbstractDict
            @test patient_resource["name"]["use"][:value] == "usual"
            @test patient_resource["name"]["text"][:value] == "Jason Argonaut"
            @test patient_resource["name"]["family"][:value] == "Argonaut"
            @test patient_resource["name"]["given"][:value] == "Jason"
            @test patient_resource["birthDate"][:value] == "1985-08-01"
        end
        Base.shred!(auth)
        Base.shred!(client)
    end

    for i in 1:length(auths)
        Base.shred!(auths[i])
    end
end


text:
#EasyPlot Tools to assist with themes
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
#Different from the one in base.jl (used only in WfrmAttributes)
const ColorRef_WFA = Union{Symbol, Int, Colorant}

#Low-level waveform attributes (compared to LineAttributes & GlyphAttributes)
#NOTE: line & glyph values are "untangled" here.
mutable struct WfrmAttributes
	linestyle::Symbol
	linewidth::Float64 #[0, 10]
	linecolor::ColorRef_WFA

	glyphshape::Symbol
	glyphsize::Float64 #[0, 10]
	glyphlinecolor::ColorRef_WFA
	#glyphlinewidth::Float64 #Use linewidth
	glyphfillcolor::ColorRef_WFA
end

#==Selecting colors
===============================================================================#
getcolor(::Nothing, v) = getcolor(COLORSCHEME_DEFAULT, v) #If no theme selected
#getcolor(t::Theme, v) = getcolor(t.colorscheme, v)

#==Computing final waveform attributes
===============================================================================#
#=TODO
Probably more efficient to run through this algorithm *before* calling
EasyPlot.addwfrm.  So the _addwfrm algorithms should use WfrmAttributes

Problem?:
Cannot mirror glyph.color from line.color when glyph.color==nothing
because we assign line.color AFTER we split up a mulit-dim waveform into multiple
waveforms.
=#

#Compute final attributes, given the theme:
#Assertion: no field in result will equal "nothing".
#TODO: Obtain other defaults from Theme
function WfrmAttributes(t::Theme, la::LineAttributes, ga::GlyphAttributes;
	resolvecolors::Bool = true)
	local linestyle, glyphlinecolor
	haslinestyle = (la.style != nothing)
	hasglyph = ga.shape != :none #Alt w/nothing: (ga.shape != nothing)
	hasglyphcolor = (ga.color != nothing)

	linewidth = (nothing == la.width) ? (1) : la.width
	linecolor = (nothing == la.color) ? (:default) : la.color

	glyphshape = hasglyph ? ga.shape : (:none)
	glyphsize = (nothing == ga.size) ? (1) : ga.size

	glyphfillcolor = hasglyphcolor ? ga.color : (:transparent)

	if hasglyph && !haslinestyle
		linestyle = :none
	elseif !haslinestyle
		linestyle = :solid
	else
		linestyle = la.style
	end

	if hasglyph && hasglyphcolor
		glyphlinecolor = glyphfillcolor
	else
		glyphlinecolor = linecolor #Needs something
	end

	if resolvecolors
		linecolor = getcolor(t.colorscheme, linecolor)
		glyphlinecolor = getcolor(t.colorscheme, glyphlinecolor)
		glyphfillcolor = getcolor(t.colorscheme, glyphfillcolor)
	end
	return WfrmAttributes(linestyle, linewidth, linecolor,
		glyphshape, glyphsize, glyphlinecolor, glyphfillcolor
	)
end

#Compute final, concrete values for marker lines:
function resolve_markerline(t::Theme, la::LineAttributes)
	linestyle = (nothing == la.style) ? (:solid) : la.style
	linewidth = (nothing == la.width) ? (1) : la.width
	linecolor = (nothing == la.color) ? (:default) : la.color
	linecolor = getcolor(t.colorscheme, linecolor)
	return LineAttributes(style=linestyle, width=linewidth, color=linecolor)
end

#Last line


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/25/2014

using SISLES
using ArgParse
using HDF5, JLD
using PyPlot
#using RiskRatio
using PyCall

# Examples
#
# to get help
# $ julia simulate.jl -h
# $ julia simulate.jl run -h
# $ julia simulate.jl compare -h
# $ julia simulate.jl validate -h
# $ julia simulate.jl generate -h
# $ julia simulate.jl plot -h
#
# run/plot simulation. sample is generated randomly.
# $ julia simulate.jl run
# $ julia simulate.jl plot
#
# read 11th sample from sample files. run/plot simulation with TCAS
# $ julia simulate.jl run -t -r -n 11
# $ julia simulate.jl plot
#
# read 12th sample from sample files. run/plot simulation with/out TCAS
# $ julia simulate.jl compare -n 12
# $ julia simulate.jl plot -c
#
# read 23th sample from Lincoln Lab sample files. run/plot simulation
# $ julia simulate.jl validate -n 23
# $ julia simulate.jl plot -d
#
# generate 100 samples with 50 seconds transitions for each sample to file. "initial.txt" and "transition.txt" files are generated under current directory.
# $ julia simulate.jl generate -n 100
#
# generate 10 samples with 20 seconds transitions for each sample to file. "init.txt" and "tran.txt" files are generated under current directory.
# $ julia simulate.jl generate --init_file "init.txt" --ninit 10 --tran_file "tran.txt" --ntran 20 --param_file "../Encounter/CorrAEMImpl/params/cor.txt"

function initialize_simulation(; bValidate = false, bReadSampleFromFile = false, initial_sample_filename = "initial.txt", transition_sample_filename = "transition.txt")

    SISLES_PATH = Pkg.dir("SISLES", "src")
    number_of_aircraft = 2

    if bValidate
        aem = LLAEM(number_of_aircraft, "$SISLES_PATH/Encounter/LLAEMImpl/data/cor_enc.txt", 50, ["$SISLES_PATH/Encounter/LLAEMImpl/data/cor_ac1.txt", "$SISLES_PATH/Encounter/LLAEMImpl/data/cor_ac2.txt"], [50, 50])
    else
        if bReadSampleFromFile
            aem = CorrAEM("$SISLES_PATH/Encounter/CorrAEMImpl/params/cor.txt", initial_sample_filename, transition_sample_filename)
            #validate(aem)
        else
            aem = CorrAEM("$SISLES_PATH/Encounter/CorrAEMImpl/params/cor.txt")
        end
    end

    pr_1 = SimplePilotResponse()
    pr_2 = SimplePilotResponse()

    #adm_1 = SimpleADM()
    #adm_2 = SimpleADM()
    adm_1 = SimpleADM(number_of_substeps = 1)
    adm_2 = SimpleADM(number_of_substeps = 1)

    as = AirSpace(number_of_aircraft)

    sr_1 = SimpleTCASSensor(1)
    sr_2 = SimpleTCASSensor(2)

    cas_1 = SimpleTCAS()
    cas_2 = SimpleTCAS()

    sim = TCASSimulator()

    sim.parameters.em = aem
    sim.parameters.pr = [pr_1, pr_2]
    sim.parameters.dm = [adm_1, adm_2]
    sim.parameters.wm = as
    sim.parameters.sr = [sr_1, sr_2]
    sim.parameters.cas = [cas_1, cas_2]
    sim.parameters.number_of_aircraft = number_of_aircraft

    return sim
end

function run_simulation(; bTCAS = false, bReadSampleFromFile = false, initial_sample_filename = "initial.txt", transition_sample_filename = "transition.txt", sample_number = 1)

    sim = initialize_simulation(bReadSampleFromFile = bReadSampleFromFile, initial_sample_filename = initial_sample_filename, transition_sample_filename = transition_sample_filename)

    aem = sim.parameters.em

    AC1_trajectory_ = Vector{Float64}[]
    AC2_trajectory_ = Vector{Float64}[]

    adm_1, adm_2 = sim.parameters.dm

    addObserver(adm_1, "adm_simulate", x -> push!(AC1_trajectory_, x))
    addObserver(adm_2, "adm_simulate", x -> push!(AC2_trajectory_, x))

    if bReadSampleFromFile
        simulate(sim, bTCAS = bTCAS, sample_number = sample_number)
    else
        simulate(sim, bTCAS = bTCAS)
    end

    labels = ["A, L, chi(1: front, 2: back), beta(deg), C1, C2, hmd(ft), vmd(ft)", "time(sec), x_1(ft), y_1(ft), h_1(ft), x_2(ft), y_2(ft), h_2(ft)"]

    initial = [aem.A, aem.L, aem.geometry_at_TCA[1], aem.geometry_at_TCA[2], aem.C[1], aem.C[2], aem.geometry_at_TCA[3], aem.geometry_at_TCA[4]]

    AC1_trajectory = zeros(length(AC1_trajectory_), 4)
    AC2_trajectory = zeros(length(AC2_trajectory_), 4)

    for i = 1:length(AC1_trajectory_)
        AC1_trajectory[i, :] = AC1_trajectory_[i]
        AC2_trajectory[i, :] = AC2_trajectory_[i]
    end

    SimulationResult = Any[labels, initial, AC1_trajectory, AC2_trajectory]

    save("result.jld", "data", SimulationResult)
end

# Multiple comparisons

function run_multiple_tcas_sims(; initial_sample_filename = "initial.txt", transition_sample_filename = "transition.txt", sample_number = 10)

    sim = initialize_simulation(bReadSampleFromFile = true, initial_sample_filename = initial_sample_filename, transition_sample_filename = transition_sample_filename)

    aem = sim.parameters.em

    AC1_trajectory_ML = []
    AC2_trajectory_ML = []

    AC1_trajectory_NC = []
    AC2_trajectory_NC = []

    AC1_trajectory_ = Vector{Float64}[]
    AC2_trajectory_ = Vector{Float64}[]

    adm_1, adm_2 = sim.parameters.dm

    addObserver(adm_1, "adm_simulate", x -> push!(AC1_trajectory_, x))
    addObserver(adm_2, "adm_simulate", x -> push!(AC2_trajectory_, x))

    for i = [1:sample_number]
      simulate(sim, bTCAS = false, sample_number = sample_number)
      AC1_trajectory = zeros(length(AC1_trajectory_), 4)
      AC2_trajectory = zeros(length(AC2_trajectory_), 4)

      for i = 1:length(AC1_trajectory_)
        AC1_trajectory[i, :] = AC1_trajectory_[i]
        AC2_trajectory[i, :] = AC2_trajectory_[i]
      end

      push!(AC1_trajectory_NC, AC1_trajectory)
      push!(AC2_trajectory_NC, AC2_trajectory)

      AC1_trajectory_ = Vector{Float64}[]
      AC2_trajectory_ = Vector{Float64}[]

    simulate(sim, bTCAS = true, sample_number = sample_number)

    AC1_trajectory_tcas = zeros(length(AC1_trajectory_), 4)
    AC2_trajectory_tcas = zeros(length(AC2_trajectory_), 4)

    for i = 1:length(AC1_trajectory_)
        AC1_trajectory_tcas[i, :] = AC1_trajectory_[i]
        AC2_trajectory_tcas[i, :] = AC2_trajectory_[i]
    end

    push!(AC1_trajectory_ML, AC1_trajectory_tcas)
    push!(AC2_trajectory_ML, AC2_trajectory_tcas)

    AC1_trajectory_tcas = zeros(length(AC1_trajectory_), 4)
    AC2_trajectory_tcas = zeros(length(AC2_trajectory_), 4)

    labels = ["A, L, chi(1: front, 2: back), beta(deg), C1, C2, hmd(ft), vmd(ft)", "time(sec), x_1(ft), y_1(ft), h_1(ft), x_2(ft), y_2(ft), h_2(ft)"]

    initial = [aem.A, aem.L, aem.geometry_at_TCA[1], aem.geometry_at_TCA[2], aem.C[1], aem.C[2], aem.geometry_at_TCA[3], aem.geometry_at_TCA[4]]

    SimulationResult = Any[labels, initial, AC1_trajectory_NC, AC2_trajectory_NC, AC1_trajectory_ML, AC2_trajectory_ML]

    save("mresult.jld", "data", SimulationResult)

    end
end

# Multiple comparisons

function run_simulation_for_comparison(; initial_sample_filename = "initial.txt", transition_sample_filename = "transition.txt", sample_number = 1)

    sim = initialize_simulation(bReadSampleFromFile = true, initial_sample_filename = initial_sample_filename, transition_sample_filename = transition_sample_filename)

    aem = sim.parameters.em

    AC1_trajectory_ = Vector{Float64}[]
    AC2_trajectory_ = Vector{Float64}[]

    adm_1, adm_2 = sim.parameters.dm

    addObserver(adm_1, "adm_simulate", x -> push!(AC1_trajectory_, x))
    addObserver(adm_2, "adm_simulate", x -> push!(AC2_trajectory_, x))

    simulate(sim, bTCAS = false, sample_number = sample_number)

    AC1_trajectory = zeros(length(AC1_trajectory_), 4)
    AC2_trajectory = zeros(length(AC2_trajectory_), 4)

    for i = 1:length(AC1_trajectory_)
        AC1_trajectory[i, :] = AC1_trajectory_[i]
        AC2_trajectory[i, :] = AC2_trajectory_[i]
    end

    AC1_trajectory_ = Vector{Float64}[]
    AC2_trajectory_ = Vector{Float64}[]

    simulate(sim, bTCAS = true, sample_number = sample_number)

    AC1_trajectory_tcas = zeros(length(AC1_trajectory_), 4)
    AC2_trajectory_tcas = zeros(length(AC2_trajectory_), 4)

    for i = 1:length(AC1_trajectory_)
        AC1_trajectory_tcas[i, :] = AC1_trajectory_[i]
        AC2_trajectory_tcas[i, :] = AC2_trajectory_[i]
    end

    labels = ["A, L, chi(1: front, 2: back), beta(deg), C1, C2, hmd(ft), vmd(ft)", "time(sec), x_1(ft), y_1(ft), h_1(ft), x_2(ft), y_2(ft), h_2(ft)"]

    initial = [aem.A, aem.L, aem.geometry_at_TCA[1], aem.geometry_at_TCA[2], aem.C[1], aem.C[2], aem.geometry_at_TCA[3], aem.geometry_at_TCA[4]]

    SimulationResult = Any[labels, initial, AC1_trajectory, AC2_trajectory, AC1_trajectory_tcas, AC2_trajectory_tcas]

    save("result.jld", "data", SimulationResult)
end

function run_simulation_for_validation(; sample_number = 1)

    sim = initialize_simulation(bValidate = true)

    aem = sim.parameters.em

    AC1_trajectory_ = Vector{Float64}[]
    AC2_trajectory_ = Vector{Float64}[]

    adm_1, adm_2 = sim.parameters.dm

    addObserver(adm_1, "adm_simulate", x -> push!(AC1_trajectory_, x))
    addObserver(adm_2, "adm_simulate", x -> push!(AC2_trajectory_, x))

    simulate(sim, sample_number = sample_number)

    AC1_trajectory = zeros(length(AC1_trajectory_), 4)
    AC2_trajectory = zeros(length(AC2_trajectory_), 4)

    for i = 1:length(AC1_trajectory_)
        AC1_trajectory[i, :] = AC1_trajectory_[i]
        AC2_trajectory[i, :] = AC2_trajectory_[i]
    end

    AC1_trajectory_ll = getTrajectory(aem, 1)
    AC2_trajectory_ll = getTrajectory(aem, 2)

    labels = ["A, L, chi(1: front, 2: back), beta(deg), C1, C2, hmd(ft), vmd(ft)", "time(sec), x_1(ft), y_1(ft), h_1(ft), x_2(ft), y_2(ft), h_2(ft)"]

    initial = [aem.A, aem.L, aem.geometry_at_TCA[1], aem.geometry_at_TCA[2], aem.C[1], aem.C[2], aem.geometry_at_TCA[3], aem.geometry_at_TCA[4]]

    SimulationResult = Any[labels, initial, AC1_trajectory, AC2_trajectory, AC1_trajectory_ll, AC2_trajectory_ll]

    save("result.jld", "data", SimulationResult)
end

function generate_samples_to_file(initial_sample_filename = "initial.txt",
    number_of_initial_samples = 1, transition_sample_filename = "transition.txt",
    number_of_transition_samples = 50,
    parameter_file = Pkg.dir("SISLES/src/Encounter/CorrAEMImpl/params/cor.txt"))

    aem = CorrAEM(parameter_file, initial_sample_filename, number_of_initial_samples,
        transition_sample_filename, number_of_transition_samples)

    Encounter.generateEncountersToFile(aem)
end

function plot_result()

    SimulationResult = load("result.jld", "data")

    labels, initial, AC1_trajectory, AC2_trajectory = SimulationResult

    println(labels[1])
    println(initial')

    index_tca = find(x -> abs(x - 40) < 0.001, AC1_trajectory[:, 1])

    figure(1)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    axis("equal")
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2"], loc = 0)

    figure(2)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    axis("equal")
    l = initial[7] * 1.2
    axis([-l, l, -l, l])
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2"], loc = 0)

    figure(3)
    plot(AC1_trajectory[:, 1] , AC1_trajectory[:, 4], AC2_trajectory[:, 1], AC2_trajectory[:, 4])
    plot(AC1_trajectory[index_tca, 1], AC1_trajectory[index_tca, 4], "*", AC2_trajectory[index_tca, 1], AC2_trajectory[index_tca, 4], "*")
    grid(true)
    xlabel("time(sec)")
    ylabel("h(ft)")
    legend(["AC1", "AC2"], loc = 0)
end

function plot_result_for_comparison()

    SimulationResult = load("result.jld", "data")

    labels, initial, AC1_trajectory, AC2_trajectory, AC1_trajectory_tcas, AC2_trajectory_tcas = SimulationResult

    println(labels[1])
    println(initial')

    index_tca = find(x -> abs(x - 40) < 0.001, AC1_trajectory[:, 1])
    index_tca_tcas = find(x -> abs(x - 40) < 0.001, AC1_trajectory_tcas[:, 1])

    figure(1)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory_tcas[:, 3], AC1_trajectory_tcas[:, 2], "--", AC2_trajectory_tcas[:, 3], AC2_trajectory_tcas[:, 2], "--")
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    plot(AC1_trajectory_tcas[index_tca_tcas, 3], AC1_trajectory_tcas[index_tca_tcas, 2], "x", AC2_trajectory_tcas[index_tca_tcas, 3], AC2_trajectory_tcas[index_tca_tcas, 2], "x")
    axis("equal")
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2", "AC1_TCAS", "AC2_TCAS"], loc = 0)

    figure(2)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory_tcas[:, 3], AC1_trajectory_tcas[:, 2], "--", AC2_trajectory_tcas[:, 3], AC2_trajectory_tcas[:, 2], "--")
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    plot(AC1_trajectory_tcas[index_tca_tcas, 3], AC1_trajectory_tcas[index_tca_tcas, 2], "x", AC2_trajectory_tcas[index_tca_tcas, 3], AC2_trajectory_tcas[index_tca_tcas, 2], "x")
    axis("equal")
    l = initial[7] * 1.2
    axis([-l, l, -l, l])
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2", "AC1_TCAS", "AC2_TCAS"], loc = 0)

    figure(3)
    h1 = AC1_trajectory_tcas[1, 4] - AC1_trajectory[1, 4]
    h2 = AC2_trajectory_tcas[1, 4] - AC2_trajectory[1, 4]
    plot(AC1_trajectory[:, 1] , AC1_trajectory[:, 4] + h1, AC2_trajectory[:, 1], AC2_trajectory[:, 4] + h2)
    plot(AC1_trajectory_tcas[:, 1], AC1_trajectory_tcas[:, 4], "--", AC2_trajectory_tcas[:, 1], AC2_trajectory_tcas[:, 4], "--")
    plot(AC1_trajectory[index_tca, 1], AC1_trajectory[index_tca, 4] + h1, "*", AC2_trajectory[index_tca, 1], AC2_trajectory[index_tca, 4] + h2, "*")
    plot(AC1_trajectory_tcas[index_tca_tcas, 1], AC1_trajectory_tcas[index_tca_tcas, 4], "x", AC2_trajectory_tcas[index_tca_tcas, 1], AC2_trajectory_tcas[index_tca_tcas, 4], "x")
    grid(true)
    xlabel("time(sec)")
    ylabel("h(ft)")
    legend(["AC1", "AC2", "AC1_TCAS", "AC2_TCAS"], loc = 0)
end

function plot_risk_ratios()

    # Generate risk ratios

    figure(1)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory_ll[:, 3], AC1_trajectory_ll[:, 2], "--", AC2_trajectory_ll[:, 3], AC2_trajectory_ll[:, 2], "--")
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    plot(AC1_trajectory_ll[index_tca_ll, 3], AC1_trajectory_ll[index_tca_ll, 2], "x", AC2_trajectory_ll[index_tca_ll, 3], AC2_trajectory_ll[index_tca_ll, 2], "x")
    axis("equal")
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2", "AC1_LL", "AC2_LL"], loc = 0)

    figure(2)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory_ll[:, 3], AC1_trajectory_ll[:, 2], "--", AC2_trajectory_ll[:, 3], AC2_trajectory_ll[:, 2], "--")
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    plot(AC1_trajectory_ll[index_tca_ll, 3], AC1_trajectory_ll[index_tca_ll, 2], "x", AC2_trajectory_ll[index_tca_ll, 3], AC2_trajectory_ll[index_tca_ll, 2], "x")
    axis("equal")
    l = initial[7] * 1.2
    axis([-l, l, -l, l])
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2", "AC1_LL", "AC2_LL"], loc = 0)

    figure(3)
    plot(AC1_trajectory[:, 1] , AC1_trajectory[:, 4], AC2_trajectory[:, 1], AC2_trajectory[:, 4])
    plot(AC1_trajectory_ll[:, 1], AC1_trajectory_ll[:, 4], "--", AC2_trajectory_ll[:, 1], AC2_trajectory_ll[:, 4], "--")
    plot(AC1_trajectory[index_tca, 1], AC1_trajectory[index_tca, 4], "*", AC2_trajectory[index_tca, 1], AC2_trajectory[index_tca, 4], "*")
    plot(AC1_trajectory_ll[index_tca_ll, 1], AC1_trajectory_ll[index_tca_ll, 4], "x", AC2_trajectory_ll[index_tca_ll, 1], AC2_trajectory_ll[index_tca_ll, 4], "x")
    grid(true)
    xlabel("time(sec)")
    ylabel("h(ft)")
    legend(["AC1", "AC2", "AC1_LL", "AC2_LL"], loc = 0)
end

function plot_result_for_validation()

    SimulationResult = load("result.jld", "data")

    labels, initial, AC1_trajectory, AC2_trajectory, AC1_trajectory_ll, AC2_trajectory_ll = SimulationResult

    println(labels[1])
    println(initial')

    index_tca = find(x -> abs(x - 40) < 0.001, AC1_trajectory[:, 1])
    index_tca_ll = find(x -> abs(x - 40) < 0.001, AC1_trajectory_ll[:, 1])

    figure(1)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory_ll[:, 3], AC1_trajectory_ll[:, 2], "--", AC2_trajectory_ll[:, 3], AC2_trajectory_ll[:, 2], "--")
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    plot(AC1_trajectory_ll[index_tca_ll, 3], AC1_trajectory_ll[index_tca_ll, 2], "x", AC2_trajectory_ll[index_tca_ll, 3], AC2_trajectory_ll[index_tca_ll, 2], "x")
    axis("equal")
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2", "AC1_LL", "AC2_LL"], loc = 0)

    figure(2)
    plot(AC1_trajectory[:, 3] , AC1_trajectory[:, 2], AC2_trajectory[:, 3], AC2_trajectory[:, 2])
    plot(AC1_trajectory_ll[:, 3], AC1_trajectory_ll[:, 2], "--", AC2_trajectory_ll[:, 3], AC2_trajectory_ll[:, 2], "--")
    plot(AC1_trajectory[index_tca, 3], AC1_trajectory[index_tca, 2], "*", AC2_trajectory[index_tca, 3], AC2_trajectory[index_tca, 2], "*")
    plot(AC1_trajectory_ll[index_tca_ll, 3], AC1_trajectory_ll[index_tca_ll, 2], "x", AC2_trajectory_ll[index_tca_ll, 3], AC2_trajectory_ll[index_tca_ll, 2], "x")
    axis("equal")
    l = initial[7] * 1.2
    axis([-l, l, -l, l])
    grid(true)
    xlabel("y(ft)")
    ylabel("x(ft)")
    legend(["AC1", "AC2", "AC1_LL", "AC2_LL"], loc = 0)

    figure(3)
    plot(AC1_trajectory[:, 1] , AC1_trajectory[:, 4], AC2_trajectory[:, 1], AC2_trajectory[:, 4])
    plot(AC1_trajectory_ll[:, 1], AC1_trajectory_ll[:, 4], "--", AC2_trajectory_ll[:, 1], AC2_trajectory_ll[:, 4], "--")
    plot(AC1_trajectory[index_tca, 1], AC1_trajectory[index_tca, 4], "*", AC2_trajectory[index_tca, 1], AC2_trajectory[index_tca, 4], "*")
    plot(AC1_trajectory_ll[index_tca_ll, 1], AC1_trajectory_ll[index_tca_ll, 4], "x", AC2_trajectory_ll[index_tca_ll, 1], AC2_trajectory_ll[index_tca_ll, 4], "x")
    grid(true)
    xlabel("time(sec)")
    ylabel("h(ft)")
    legend(["AC1", "AC2", "AC1_LL", "AC2_LL"], loc = 0)
end

function parse_commandline()

    println("Collision Avoidance System Safety Assessment Tool -- CASSATT \r\n")

    settings = ArgParseSettings()

    settings.description = "Airspace Encounter Simulation"

    @add_arg_table settings begin
        "run"
            help = "simulate trajectories"
            action = :command
        "multiple"
            help = "run multiple simulations"
            action = :command
        "compare"
            help = "compare simulated trajectories with/out TCAS"
            action = :command
        "validate"
            help = "compare simulated trajectories to those from Lincoln Lab"
            action = :command
        "generate"
            help = "generate samples to file"
            action = :command
        "plot"
            help = "draw trajectories"
            action = :command
    end

    @add_arg_table settings ["multiple"] begin
        "--init_file"
          help = "initial sample file"
          arg_type = AbstractString
          default = "initial.txt"
      "--tran_file"
          help = "transition sample file"
          arg_type = AbstractString
          default = "transition.txt"
               "--sample", "-n"
            help = "sample number"
            arg_type = Int
            default = 1
    end
    @add_arg_table settings["run"] begin
        "--tcas", "-t"
            help = "simulate with TCAS"
            action = :store_true
        "--read_from_file", "-r"
            help = "read samples from file"
            action = :store_true
        "--init_file"
            help = "initial sample file"
            arg_type = AbstractString
            default = "initial.txt"
        "--tran_file"
            help = "transition sample file"
            arg_type = AbstractString
            default = "transition.txt"
        "--sample", "-n"
            help = "sample number"
            arg_type = Int
            default = 1
    end

    @add_arg_table settings["compare"] begin
        "--init_file"
            help = "initial sample file"
            arg_type = AbstractString
            default = "initial.txt"
        "--tran_file"
            help = "transition sample file"
            arg_type = AbstractString
            default = "transition.txt"
        "--sample", "-n"
            help = "sample number"
            arg_type = Int
            default = 1
    end

    @add_arg_table settings["validate"] begin
        "--sample", "-n"
            help = "sample number"
            arg_type = Int
            default = 1
    end

    @add_arg_table settings["generate"] begin
        "--init_file"
            help = "initial sample file"
            arg_type = AbstractString
            default = "initial.txt"
        "--ninit", "-n"
            help = "number of initial samples"
            arg_type = Int
            default = 1
        "--tran_file"
            help = "transition sample file"
            arg_type = AbstractString
            default = "transition.txt"
        "--ntran"
            help = "number of transition samples"
            arg_type = Int
            default = 50
        "--param_file"
            help = "encounter parameter file"
            arg_type = AbstractString
            default = Pkg.dir("SISLES/src/Encounter/CorrAEMImpl/params/cor.txt")
    end

    @add_arg_table settings["plot"] begin
        "--compare", "-c"
            help = "compare simulated trajectories with/out TCAS"
            action = :store_true
        "--validate", "-d"
            help = "compare simulated trajectories to those from Lincoln Lab"
            action = :store_true
    end

    return parse_args(settings)
end

function main()

    parsed_args = parse_commandline()

    if parsed_args["%COMMAND%"] == "run"
        run_simulation(bTCAS = parsed_args["run"]["tcas"], bReadSampleFromFile = parsed_args["run"]["read_from_file"], initial_sample_filename = parsed_args["run"]["init_file"], transition_sample_filename = parsed_args["run"]["tran_file"], sample_number = parsed_args["run"]["sample"])

    elseif parsed_args["%COMMAND%"] == "compare"
        run_simulation_for_comparison(initial_sample_filename = parsed_args["compare"]["init_file"], transition_sample_filename = parsed_args["compare"]["tran_file"], sample_number = parsed_args["compare"]["sample"])

    elseif parsed_args["%COMMAND%"] == "validate"
        run_simulation_for_validation(sample_number = parsed_args["validate"]["sample"])

    elseif parsed_args["%COMMAND%"] == "generate"
        generate_samples_to_file(parsed_args["generate"]["init_file"], parsed_args["generate"]["ninit"], parsed_args["generate"]["tran_file"], parsed_args["generate"]["ntran"], parsed_args["generate"]["param_file"])
    elseif parsed_args["%COMMAND%"] == "multiple"
        run_multiple_tcas_sims()
    elseif parsed_args["%COMMAND%"] == "plot"
        if parsed_args["plot"]["compare"]
            plot_result_for_comparison()
        elseif parsed_args["plot"]["validate"]
            plot_result_for_validation()
        else
            plot_result()
        end

        readline()
    end
end

main()


text:
function decomposeX_balance_allfactors(A,B,k,c1,c2,c3)
  nA = size(A,1)
  nB = size(B,1)

  u = ones(nA,1)./nA
  v = ones(nB,1)./nB
  du = zeros(k)
  dv = zeros(k)
  du[end] = 1
  dv[end] = 1

  eA = ones(nA,1)
  eB = ones(nB,1)

  U = zeros(nA,k)
  U[:,end] = eA
  for i = k-1:-1:2
    U[:,i] = A*U[:,i+1]
    mui = maximum(U[:,i])
    U[:,i] = U[:,i]./mui
    du[i] = du[i+1]*mui
  end
  U[:,1] = u
  rksums = zeros(k-1)
  du[1] = 1
  for i = 1:k-1
    rksums[i] = du[1]*sum(U[:,1])
    U[:,1] = A*U[:,1]
    mui = maximum(U[:,1])
    U[:,1] = U[:,1]./mui
    du[1] = du[1]*mui
  end
  usums = vec(sum(U,dims=1))
  usums = usums.*du
  # println("done with U")

  V = zeros(nB,k)
  V[:,end] = eB
  for i = k-1:-1:2
    V[:,i] = B*V[:,i+1]
    mvi = maximum(V[:,i])
    V[:,i] = V[:,i]./mvi
    dv[i] = dv[i+1]*mvi
  end
  V[:,1] = v
  hksums = zeros(k-1)
  dv[1] = 1
  for i = 1:k-1
    hksums[i] = dv[1]*sum(V[:,1])
    V[:,1] = B*V[:,1]
    mvi = maximum(V[:,1])
    V[:,1] = V[:,1]./mvi
    dv[1] = dv[1]*mvi
  end
  vsums = vec(sum(V,dims=1))
  vsums = vsums.*dv

  rhos = copy(du)
  gams = copy(dv)

    # create W
    W_prev = 1*rhos[1]*gams[1]
    W1 = 1
    W2 = 1

    scalerow = 1
    scalecol = 1

    for i = 2:k
      localk = i-1
      rk = zeros(localk)
      rk[1] = rksums[localk]
      lt = length(2:localk)
      rk[2:localk] = usums[end-lt+1:end]

      hk = zeros(localk)
      hk[1] = hksums[localk]
      hk[2:localk] = vsums[end-lt+1:end]

      Du = diagm(rhos[i]./rhos[1:i-1])
      Dv = diagm(gams[i]./gams[1:i-1])

      v1 = Dv*hk
      v2 = Du*rk

      W_curr = vcat(hcat(c1*W_prev,scalecol*c2*W_prev*v1),
               hcat(scalerow*c2*v2'*W_prev,scalerow*scalecol*c3*v2'*W_prev*v1))

#      W_curr = vcat(hcat(c1*W_prev,scalecol*c2*W_prev*Dv*hk),
#                    hcat(scalerow*c2*rk'*Du*W_prev,scalerow*scalecol*c3*rk'*Du*W_prev*Dv*hk))

      W1 = vcat(hcat(W_prev,zeros(size(W_prev,1),localk)),
                  hcat(zeros(1,size(W_prev,2)),rk'))
      W2 = vcat(hcat(c1*I(size(W_prev,1)),c2*hk),
                  hcat(c2*W_prev,c3*W_prev*hk))

      W_curr = W_curr/W_curr[end,end]

      W_prev = W_curr
    end
  return U,V,W_prev,W1,W2
end

# old version of decomposeX
function decomposeX(A,B,k,c1,c2,c3)
  println("decomposeX")
  nA = size(A,1)
  nB = size(B,1)

  u = ones(nA,1)./nA
  v = ones(nB,1)./nB

  eA = ones(nA,1)
  eB = ones(nB,1)

  U = zeros(nA,k)
  U[:,end] = eA
  for i = k-1:-1:2
    U[:,i] = A*U[:,i+1]
  end
  U[:,1] = u
  rksums = zeros(k-1)
  for i = 1:k-1
    rksums[i] = sum(U[:,1])
    U[:,1] = A*U[:,1]
  end
  usums = sum(U,1)
# println("done with U")
  V = zeros(nB,k)
  V[:,end] = eB
  for i = k-1:-1:2
    V[:,i] = B*V[:,i+1]
  end
  V[:,1] = v
  hksums = zeros(k)
  for i = 1:k-1
    hksums[i] = sum(V[:,1])
    V[:,1] = B*V[:,1]
  end
  vsums = sum(V,1)
# println("done with U")
  # create W
  W_prev = 1
  W1 = 1
  W2 = 1
  @show usums
  @show vsums

  for i = 2:k
# println("i is $i out of $k")
    localk = i-1
    rk = zeros(localk)
    rk[1] = rksums[localk]
    lt = length(2:localk)
    rk[2:localk] = usums[end-lt+1:end]
    @show rk

    # rk[1] = sum(A^(localk-1)*u)
    # for j = 2:localk
    #   rk[j] = sum(A^(localk-j)*eA)
    # end

    hk = zeros(localk)
    hk[1] = hksums[localk]
    hk[2:localk] = vsums[end-lt+1:end]
    @show hk
    # hk[1] = sum(B^(localk-1)*v)
    # for j = 2:localk
    #   hk[j] = sum(B^(localk-j)*eB)
    # end

    # println("Wcalc start")
    W_curr = vcat(hcat(c1*W_prev,c2*W_prev*hk),
                  hcat(c2*rk'*W_prev,c3*rk'*W_prev*hk))

    W1 = vcat(hcat(W_prev,zeros(size(W_prev,1),localk)),
                hcat(zeros(1,size(W_prev,2)),rk'))
    W2 = vcat(hcat(c1*I(size(W_prev,1)),c2*hk),
                hcat(c2*W_prev,c3*W_prev*hk))

    W_prev = W_curr
    # println("Wcalc end")
  end
  return U,V,W_prev,W1,W2
end

function decomposeX_normalized(A,B,k,c1,c2,c3)
  nA = size(A,1)
  nB = size(B,1)

  u = ones(nA,1)./nA
  v = ones(nB,1)./nB

  eA = ones(nA,1)
  eB = ones(nB,1)

  U = zeros(nA,k)
  U[:,end] = eA
  U[:,end] = U[:,end]./sum(U[:,end])
  for i = k-1:-1:2
    U[:,i] = A*U[:,i+1]
    U[:,i] = U[:,i]./sum(U[:,i])
  end
  U[:,1] = u
  rksums = zeros(k-1)
  for i = 1:k-1
    rksums[i] = sum(U[:,1])
    U[:,1] = A*U[:,1]
  end
  U[:,1] = U[:,1]./sum(U[:,1])
  usums = sum(U,1)
# println("done with U")
  V = zeros(nB,k)
  V[:,end] = eB
  V[:,end] = V[:,end]./sum(V[:,end])
  for i = k-1:-1:2
    V[:,i] = B*V[:,i+1]
    V[:,i] = V[:,i]./sum(V[:,i])
  end
  V[:,1] = v
  hksums = zeros(k)
  for i = 1:k-1
    hksums[i] = sum(V[:,1])
    V[:,1] = B*V[:,1]
  end
  V[:,1] = V[:,1]./sum(V[:,1])
  vsums = sum(V,1)
# println("done with U")
  # create W
  W_prev = 1
  W1 = 1
  W2 = 1
  @show usums
  @show vsums
  @show hksums
  @show rksums
  for i = 2:k
# println("i is $i out of $k")
    localk = i-1
    rk = zeros(localk)
    rk[1] = rksums[localk]
    lt = length(2:localk)
    rk[2:localk] = usums[end-lt+1:end]
    # rk[1] = sum(A^(localk-1)*u)
    # for j = 2:localk
    #   rk[j] = sum(A^(localk-j)*eA)
    # end

    hk = zeros(localk)
    hk[1] = hksums[localk]
    hk[2:localk] = vsums[end-lt+1:end]
    # hk[1] = sum(B^(localk-1)*v)
    # for j = 2:localk
    #   hk[j] = sum(B^(localk-j)*eB)
    # end

    # println("Wcalc start")
    W_curr = vcat(hcat(c1*W_prev,c2*W_prev*hk),
                  hcat(c2*rk'*W_prev,c3*rk'*W_prev*hk))

    W1 = vcat(hcat(W_prev,zeros(size(W_prev,1),localk)),
                hcat(zeros(1,size(W_prev,2)),rk'))
    W2 = vcat(hcat(c1*I(size(W_prev,1)),c2*hk),
                hcat(c2*W_prev,c3*W_prev*hk))

    W_prev = W_curr
    # println("Wcalc end")
  end
  return U,V,W_prev,W1,W2
end

function decomposeX_with_uv(A,B,k,c1,c2,c3,u,v)
  nA = size(A,1)
  nB = size(B,1)

  # u = ones(nA,1)./nA
  # v = ones(nB,1)./nB

  eA = ones(nA,1)
  eB = ones(nB,1)

  U = zeros(nA,k)
  U[:,end] = eA
  for i = k-1:-1:2
    U[:,i] = A*U[:,i+1]
  end
  U[:,1] = u
  rksums = zeros(k-1)
  for i = 1:k-1
    rksums[i] = sum(U[:,1])
    U[:,1] = A*U[:,1]
  end
  usums = sum(U,1)
# println("done with U")
  V = zeros(nB,k)
  V[:,end] = eB
  for i = k-1:-1:2
    V[:,i] = B*V[:,i+1]
  end
  V[:,1] = v
  hksums = zeros(k)
  for i = 1:k-1
    hksums[i] = sum(V[:,1])
    V[:,1] = B*V[:,1]
  end
  vsums = sum(V,1)
# println("done with U")
  # create W
  W_prev = 1
  W1 = 1
  W2 = 1
  for i = 2:k
# println("i is $i out of $k")
    localk = i-1
    rk = zeros(localk)
    rk[1] = rksums[localk]
    lt = length(2:localk)
    rk[2:localk] = usums[end-lt+1:end]
    # rk[1] = sum(A^(localk-1)*u)
    # for j = 2:localk
    #   rk[j] = sum(A^(localk-j)*eA)
    # end

    hk = zeros(localk)
    hk[1] = hksums[localk]
    hk[2:localk] = vsums[end-lt+1:end]
    # hk[1] = sum(B^(localk-1)*v)
    # for j = 2:localk
    #   hk[j] = sum(B^(localk-j)*eB)
    # end

    # println("Wcalc start")
    W_curr = vcat(hcat(c1*W_prev,c2*W_prev*hk),
                  hcat(c2*rk'*W_prev,c3*rk'*W_prev*hk))

    W1 = vcat(hcat(W_prev,zeros(size(W_prev,1),localk)),
                hcat(zeros(1,size(W_prev,2)),rk'))
    W2 = vcat(hcat(c1*I(size(W_prev,1)),c2*hk),
                hcat(c2*W_prev,c3*W_prev*hk))

    W_prev = W_curr
    # println("Wcalc end")
  end
  return U,V,W_prev,W1,W2
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME