text:
##
## Join / merge
##

# Like similar, but returns a nullable array
similar_nullable{T}(dv::AbstractArray{T}, dims::Union{Int, Tuple{Vararg{Int}}}) =
    NullableArray{T}(dims)

similar_nullable{T<:Nullable}(dv::AbstractArray{T}, dims::Union{Int, Tuple{Vararg{Int}}}) =
    NullableArray{eltype(T)}(dims)

similar_nullable{T}(dv::CategoricalArray{T}, dims::Union{Int, Tuple{Vararg{Int}}}) =
    NullableCategoricalArray{T}(dims)

similar_nullable{T}(dv::NullableCategoricalArray{T}, dims::Union{Int, Tuple{Vararg{Int}}}) =
    NullableCategoricalArray{T}(dims)

# helper structure for DataTables joining
immutable DataTableJoiner{DT1<:AbstractDataTable, DT2<:AbstractDataTable}
    dtl::DT1
    dtr::DT2
    dtl_on::DT1
    dtr_on::DT2
    on_cols::Vector{Symbol}

    function DataTableJoiner(dtl::DT1, dtr::DT2, on::Union{Symbol,Vector{Symbol}})
        on_cols = isa(on, Symbol) ? [on] : on
        new(dtl, dtr, dtl[on_cols], dtr[on_cols], on_cols)
    end
end

DataTableJoiner{DT1<:AbstractDataTable, DT2<:AbstractDataTable}(dtl::DT1, dtr::DT2, on::Union{Symbol,Vector{Symbol}}) =
    DataTableJoiner{DT1,DT2}(dtl, dtr, on)

# helper map between the row indices in original and joined table
immutable RowIndexMap
    "row indices in the original table"
    orig::Vector{Int}
    "row indices in the resulting joined table"
    join::Vector{Int}
end

Base.length(x::RowIndexMap) = length(x.orig)

# composes the joined data table using the maps between the left and right
# table rows and the indices of rows in the result
function compose_joined_table(joiner::DataTableJoiner, kind::Symbol,
                              left_ixs::RowIndexMap, leftonly_ixs::RowIndexMap,
                              right_ixs::RowIndexMap, rightonly_ixs::RowIndexMap)
    @assert length(left_ixs) == length(right_ixs)
    # compose left half of the result taking all left columns
    all_orig_left_ixs = vcat(left_ixs.orig, leftonly_ixs.orig)

    ril = length(right_ixs)
    lil = length(left_ixs)
    loil = length(leftonly_ixs)
    roil = length(rightonly_ixs)

    if loil > 0
        # combine the matched (left_ixs.orig) and non-matched (leftonly_ixs.orig) indices of the left table rows
        # preserving the original rows order
        all_orig_left_ixs = similar(left_ixs.orig, lil + loil)
        @inbounds all_orig_left_ixs[left_ixs.join] = left_ixs.orig
        @inbounds all_orig_left_ixs[leftonly_ixs.join] = leftonly_ixs.orig
    else
        # the result contains only the left rows that are matched to right rows (left_ixs)
        all_orig_left_ixs = left_ixs.orig # no need to copy left_ixs.orig as it's not used elsewhere
    end
    # permutation to swap rightonly and leftonly rows
    right_perm = vcat(1:ril, ril+roil+1:ril+roil+loil, ril+1:ril+roil)
    if length(leftonly_ixs) > 0
        # compose right_perm with the permutation that restores left rows order
        right_perm[vcat(right_ixs.join, leftonly_ixs.join)] = right_perm[1:ril+loil]
    end
    all_orig_right_ixs = vcat(right_ixs.orig, rightonly_ixs.orig)

    # compose right half of the result taking all right columns excluding on
    dtr_noon = without(joiner.dtr, joiner.on_cols)

    nrow = length(all_orig_left_ixs) + roil
    @assert nrow == length(all_orig_right_ixs) + loil
    ncleft = ncol(joiner.dtl)
    cols = Vector{Any}(ncleft + ncol(dtr_noon))
    _similar = kind == :inner ? similar : similar_nullable
    for (i, col) in enumerate(columns(joiner.dtl))
        cols[i] = _similar(col, nrow)
        fillcolumn!(cols[i], col, all_orig_left_ixs)
    end
    for (i, col) in enumerate(columns(dtr_noon))
        cols[i+ncleft] = _similar(col, nrow)
        fillcolumn!(cols[i+ncleft], col, all_orig_right_ixs)
        permute!(cols[i+ncleft], right_perm)
    end
    res = DataTable(cols, vcat(names(joiner.dtl), names(dtr_noon)))

    if length(rightonly_ixs.join) > 0
        # some left rows are nulls, so the values of the "on" columns
        # need to be taken from the right
        for (on_col_ix, on_col) in enumerate(joiner.on_cols)
            # fix the result of the rightjoin by taking the nonnull values from the right table
            offset = nrow - length(rightonly_ixs.orig)
            fillcolumn!(res[on_col], joiner.dtr_on[on_col_ix], rightonly_ixs.orig, offset)
        end
    end
    return res
end

function fillcolumn!{T1, T2}(dtcol::AbstractVector{T1}, refcol::AbstractVector{T2},
                             indices::Vector{Int}, offset::Int=0)
    @inbounds for (j, k) in enumerate(indices)
        dtcol[j+offset] = refcol[k]
    end
end

# map the indices of the left and right joined tables
# to the indices of the rows in the resulting table
# if `nothing` is given, the corresponding map is not built
function update_row_maps!(left_table::AbstractDataTable,
                          right_table::AbstractDataTable,
                          right_dict::RowGroupDict,
                          left_ixs::Union{Void, RowIndexMap},
                          leftonly_ixs::Union{Void, RowIndexMap},
                          right_ixs::Union{Void, RowIndexMap},
                          rightonly_mask::Union{Void, Vector{Bool}})
    # helper functions
    @inline update!(ixs::Void, orig_ix::Int, join_ix::Int, count::Int = 1) = nothing
    @inline function update!(ixs::RowIndexMap, orig_ix::Int, join_ix::Int, count::Int = 1)
        n = length(ixs.orig)
        resize!(ixs.orig, n+count)
        ixs.orig[n+1:end] = orig_ix
        append!(ixs.join, join_ix:(join_ix+count-1))
        ixs
    end
    @inline update!(ixs::Void, orig_ixs::AbstractArray, join_ix::Int) = nothing
    @inline function update!(ixs::RowIndexMap, orig_ixs::AbstractArray, join_ix::Int)
        append!(ixs.orig, orig_ixs)
        append!(ixs.join, join_ix:(join_ix+length(orig_ixs)-1))
        ixs
    end
    @inline update!(ixs::Void, orig_ixs::AbstractArray) = nothing
    @inline update!(mask::Vector{Bool}, orig_ixs::AbstractArray) = (mask[orig_ixs] = false)

    # iterate over left rows and compose the left<->right index map
    right_dict_cols = ntuple(i -> right_dict.dt[i], ncol(right_dict.dt))
    left_table_cols = ntuple(i -> left_table[i], ncol(left_table))
    next_join_ix = 1
    for l_ix in 1:nrow(left_table)
        r_ixs = findrows(right_dict, left_table, right_dict_cols, left_table_cols, l_ix)
        if isempty(r_ixs)
            update!(leftonly_ixs, l_ix, next_join_ix)
            next_join_ix += 1
        else
            update!(left_ixs, l_ix, next_join_ix, length(r_ixs))
            update!(right_ixs, r_ixs, next_join_ix)
            update!(rightonly_mask, r_ixs)
            next_join_ix += length(r_ixs)
        end
    end
end

# map the row indices of the left and right joined tables
# to the indices of rows in the resulting table
# returns the 4-tuple of row indices maps for
# - matching left rows
# - non-matching left rows
# - matching right rows
# - non-matching right rows
# if false is provided, the corresponding map is not built and the
# tuple element is empty RowIndexMap
function update_row_maps!(left_table::AbstractDataTable,
                          right_table::AbstractDataTable,
                          right_dict::RowGroupDict,
                          map_left::Bool, map_leftonly::Bool,
                          map_right::Bool, map_rightonly::Bool)
    init_map(dt::AbstractDataTable, init::Bool) = init ?
        RowIndexMap(sizehint!(Vector{Int}(), nrow(dt)),
                    sizehint!(Vector{Int}(), nrow(dt))) : nothing
    to_bimap(x::RowIndexMap) = x
    to_bimap(::Void) = RowIndexMap(Vector{Int}(), Vector{Int}())

    # init maps as requested
    left_ixs = init_map(left_table, map_left)
    leftonly_ixs = init_map(left_table, map_leftonly)
    right_ixs = init_map(right_table, map_right)
    rightonly_mask = map_rightonly ? fill(true, nrow(right_table)) : nothing
    update_row_maps!(left_table, right_table, right_dict, left_ixs, leftonly_ixs, right_ixs, rightonly_mask)
    if map_rightonly
        rightonly_orig_ixs = find(rightonly_mask)
        rightonly_ixs = RowIndexMap(rightonly_orig_ixs,
                                    collect(length(right_ixs.orig) +
                                            (leftonly_ixs === nothing ? 0 : length(leftonly_ixs)) +
                                            (1:length(rightonly_orig_ixs))))
    else
        rightonly_ixs = nothing
    end

    return to_bimap(left_ixs), to_bimap(leftonly_ixs), to_bimap(right_ixs), to_bimap(rightonly_ixs)
end

"""
Join two DataTables

```julia
join(dt1::AbstractDataTable,
     dt2::AbstractDataTable;
     on::Union{Symbol, Vector{Symbol}} = Symbol[],
     kind::Symbol = :inner)
```

### Arguments

* `dt1`, `dt2` : the two AbstractDataTables to be joined

### Keyword Arguments

* `on` : a Symbol or Vector{Symbol}, the column(s) used as keys when
  joining; required argument except for `kind = :cross`

* `kind` : the type of join, options include:

  - `:inner` : only include rows with keys that match in both `dt1`
    and `dt2`, the default
  - `:outer` : include all rows from `dt1` and `dt2`
  - `:left` : include all rows from `dt1`
  - `:right` : include all rows from `dt2`
  - `:semi` : return rows of `dt1` that match with the keys in `dt2`
  - `:anti` : return rows of `dt1` that do not match with the keys in `dt2`
  - `:cross` : a full Cartesian product of the key combinations; every
    row of `dt1` is matched with every row of `dt2`

For the three join operations that may introduce missing values (`:outer`, `:left`,
and `:right`), all columns of the returned data table will be nullable.

### Result

* `::DataTable` : the joined DataTable

### Examples

```julia
name = DataTable(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
job = DataTable(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

join(name, job, on = :ID)
join(name, job, on = :ID, kind = :outer)
join(name, job, on = :ID, kind = :left)
join(name, job, on = :ID, kind = :right)
join(name, job, on = :ID, kind = :semi)
join(name, job, on = :ID, kind = :anti)
join(name, job, kind = :cross)
```

"""
function Base.join(dt1::AbstractDataTable,
                   dt2::AbstractDataTable;
                   on::Union{Symbol, Vector{Symbol}} = Symbol[],
                   kind::Symbol = :inner)
    if kind == :cross
        (on == Symbol[]) || throw(ArgumentError("Cross joins don't use argument 'on'."))
        return crossjoin(dt1, dt2)
    elseif on == Symbol[]
        throw(ArgumentError("Missing join argument 'on'."))
    end

    joiner = DataTableJoiner(dt1, dt2, on)

    if kind == :inner
        compose_joined_table(joiner, kind, update_row_maps!(joiner.dtl_on, joiner.dtr_on,
                                                            group_rows(joiner.dtr_on),
                                                            true, false, true, false)...)
    elseif kind == :left
        compose_joined_table(joiner, kind, update_row_maps!(joiner.dtl_on, joiner.dtr_on,
                                                            group_rows(joiner.dtr_on),
                                                            true, true, true, false)...)
    elseif kind == :right
        compose_joined_table(joiner, kind, update_row_maps!(joiner.dtr_on, joiner.dtl_on,
                                                            group_rows(joiner.dtl_on),
                                                            true, true, true, false)[[3, 4, 1, 2]]...)
    elseif kind == :outer
        compose_joined_table(joiner, kind, update_row_maps!(joiner.dtl_on, joiner.dtr_on,
                                                            group_rows(joiner.dtr_on),
                                                            true, true, true, true)...)
    elseif kind == :semi
        # hash the right rows
        dtr_on_grp = group_rows(joiner.dtr_on)
        # iterate over left rows and leave those found in right
        left_ixs = Vector{Int}()
        sizehint!(left_ixs, nrow(joiner.dtl))
        dtr_on_grp_cols = ntuple(i -> dtr_on_grp.dt[i], ncol(dtr_on_grp.dt))
        dtl_on_cols = ntuple(i -> joiner.dtl_on[i], ncol(joiner.dtl_on))
        @inbounds for l_ix in 1:nrow(joiner.dtl_on)
            if findrow(dtr_on_grp, joiner.dtl_on, dtr_on_grp_cols, dtl_on_cols, l_ix) != 0
                push!(left_ixs, l_ix)
            end
        end
        return joiner.dtl[left_ixs, :]
    elseif kind == :anti
        # hash the right rows
        dtr_on_grp = group_rows(joiner.dtr_on)
        # iterate over left rows and leave those not found in right
        leftonly_ixs = Vector{Int}()
        sizehint!(leftonly_ixs, nrow(joiner.dtl))
        dtr_on_grp_cols = ntuple(i -> dtr_on_grp.dt[i], ncol(dtr_on_grp.dt))
        dtl_on_cols = ntuple(i -> joiner.dtl_on[i], ncol(joiner.dtl_on))
        @inbounds for l_ix in 1:nrow(joiner.dtl_on)
            if findrow(dtr_on_grp, joiner.dtl_on, dtr_on_grp_cols, dtl_on_cols, l_ix) == 0
                push!(leftonly_ixs, l_ix)
            end
        end
        return joiner.dtl[leftonly_ixs, :]
    else
        throw(ArgumentError("Unknown kind of join requested: $kind"))
    end
end

function crossjoin(dt1::AbstractDataTable, dt2::AbstractDataTable)
    r1, r2 = size(dt1, 1), size(dt2, 1)
    cols = Any[[repeat(c, inner=r2) for c in columns(dt1)];
               [repeat(c, outer=r1) for c in columns(dt2)]]
    colindex = merge(index(dt1), index(dt2))
    DataTable(cols, colindex)
end


text:
@testset "Destroy Function Tests" begin

  @test isdefined(Julz, :destroy) == true

  cd("dummy") do

    file_type = "type"

    file_name = "foo"

    src_file = "src/$(EnglishText.pluralize(file_type))/$file_name.jl"
    test_file = "test/$(EnglishText.pluralize(file_type))/$(file_name)_test.jl"

    rm(src_file, force=true)
    rm(test_file, force=true)

    Julz.generate(file_type, file_name)

    @test isfile(src_file)
    @test isfile(test_file)

    Julz.destroy(file_type, file_name)

    @test !isfile(src_file)
    @test !isfile(test_file)

  end

end

ext: jl


text:
using Oceananigans.Architectures
using Oceananigans.Architectures: architecture, arch_array
using Oceananigans.Grids: interior_parent_indices, topology
using Oceananigans.Utils: heuristic_workgroup
using KernelAbstractions: @kernel, @index
using IterativeSolvers, SparseArrays, LinearAlgebra
using CUDA, CUDA.CUSPARSE

mutable struct HeptadiagonalIterativeSolver{G, R, L, D, M, P, PM, PS, I, T, F}
                       grid :: G
               problem_size :: R
        matrix_constructors :: L
                   diagonal :: D
                     matrix :: M
             preconditioner :: P
      preconditioner_method :: PM
    preconditioner_settings :: PS
           iterative_solver :: I
                  tolerance :: T
                previous_Δt :: F
         maximum_iterations :: Int
end

"""
    HeptadiagonalIterativeSolver(coeffs;
                                 grid,
                                 iterative_solver = cg,
                                 maximum_iterations = prod(size(grid)),
                                 tolerance = 1e-13,
                                 reduced_dim = (false, false, false),
                                 placeholder_timestep = -1.0,
                                 preconditioner_method = :Default,
                                 preconditioner_settings = nothing)

`HeptadiagonalIterativeSolver` is a framework to solve the problem `A * x = b`
(provided that `A` is a symmetric matrix).

The solver relies on sparse version of the matrix `A` which are defined by the
field matrix_constructors.

In particular, given coefficients `Ax`, `Ay`, `Az`, `C`, `D`, the solved problem will be

To have the equation solved on Center, Center, Center, the coefficients should be specified
as follows:

- `Ax` -> Face, Center, Center
- `Ay` -> Center, Face, Center
- `Az` -> Center, Center, Face
- `C`  -> Center, Center, Center
- `D`  -> Center, Center, Center

```julia
Axᵢ₊₁ ηᵢ₊₁ + Axᵢ ηᵢ₋₁ + Ayⱼ₊₁ ηⱼ₊₁ + Ayⱼ ηⱼ₋₁ + Azₖ₊₁ ηₖ₊₁ + Azₖ ηⱼ₋₁
- 2 ( Axᵢ₊₁ + Axᵢ + Ayⱼ₊₁ + Ayⱼ + Azₖ₊₁ + Azₖ ) ηᵢⱼₖ
+   ( Cᵢⱼₖ + Dᵢⱼₖ/Δt^2 ) ηᵢⱼₖ = b
```

`solver.matrix` is precomputed with a value of `Δt = -1.0`

The sparse matrix `A` can be constructed with

- `SparseMatrixCSC(constructors...)` for CPU
- `CuSparseMatrixCSC(constructors...)` for GPU

The constructors are calculated based on the pentadiagonal coeffients passed as an input
(`matrix_from_coefficients`).

The diagonal term `- Az / (g * Δt²)` is added later on during the time stepping
to allow for variable time step. It is updated only when the previous time step
is different (`Δt_previous != Δt`).

Preconditioning is done through the incomplete LU factorization.

It works for GPU, but it relies on serial backward and forward substitution which are very
heavy and destroy all the computational advantage, therefore it is switched off until a
parallel backward/forward substitution is implemented. It is also updated based on the
matrix when `Δt != Δt_previous`

The iterative_solver used can is to be chosen from the IterativeSolvers.jl package.
The default solver is a Conjugate Gradient (cg)

```julia
solver = HeptadiagonalIterativeSolver((Ax, Ay, Az, C, D), grid = grid)
```
"""
function HeptadiagonalIterativeSolver(coeffs;
                                      grid,
                                      iterative_solver = cg,
                                      maximum_iterations = prod(size(grid)),
                                      tolerance = 1e-13,
                                      reduced_dim = (false, false, false),
                                      placeholder_timestep = -1.0,
                                      preconditioner_method = :Default,
                                      preconditioner_settings = nothing)

    arch = architecture(grid)

    matrix_constructors, diagonal, problem_size = matrix_from_coefficients(arch, grid, coeffs, reduced_dim)

    # for the moment, placeholder preconditioner and matrix are calculated using a "placeholder" timestep of - 1
    # They will be recalculated before the first time step of the simulation

    placeholder_constructors = deepcopy(matrix_constructors)
    update_diag!(placeholder_constructors, arch, problem_size, diagonal, 1.0)

    placeholder_matrix = arch_sparse_matrix(arch, placeholder_constructors)

    settings       = validate_settings(Val(preconditioner_method), arch, preconditioner_settings)
    reduced_matrix = arch_sparse_matrix(arch, speye(eltype(grid), 2))
    preconditioner = build_preconditioner(Val(preconditioner_method), reduced_matrix, settings)

    return HeptadiagonalIterativeSolver(grid,
                                        problem_size,
                                        matrix_constructors,
                                        diagonal,
                                        placeholder_matrix,
                                        preconditioner,
                                        preconditioner_method,
                                        settings,
                                        iterative_solver,
                                        tolerance,
                                        placeholder_timestep,
                                        maximum_iterations)
end

function matrix_from_coefficients(arch, grid, coeffs, reduced_dim)
    Ax, Ay, Az, C, D = coeffs

    Ax = arch_array(CPU(), Ax)
    Ay = arch_array(CPU(), Ay)
    Az = arch_array(CPU(), Az)
    C  = arch_array(CPU(), C)

    N = size(grid)

    topo = topology(grid)

    dims = validate_laplacian_direction.(N, topo, reduced_dim)
    Nx, Ny, Nz = N = validate_laplacian_size.(N, dims)
    M    = prod(N)
    diag = arch_array(arch, zeros(eltype(grid), M))

    # the following coefficients are the diagonals of the sparse matrix:
    #  - coeff_d is the main diagonal (coefficents of ηᵢⱼₖ)
    #  - coeff_x are the coefficients in the x-direction (coefficents of ηᵢ₋₁ⱼₖ and ηᵢ₊₁ⱼₖ)
    #  - coeff_y are the coefficients in the y-direction (coefficents of ηᵢⱼ₋₁ₖ and ηᵢⱼ₊₁ₖ)
    #  - coeff_z are the coefficients in the z-direction (coefficents of ηᵢⱼₖ₋₁ and ηᵢⱼₖ₊₁)
    #  - periodic boundaries are stored in coeff_bound_

    # position of diagonals for coefficients pos[1] and their boundary pos[2]
    posx = (1, Nx-1)
    posy = (1, Ny-1) .* Nx
    posz = (1, Nz-1) .* Nx .* Ny

    coeff_d       = zeros(eltype(grid), M)
    coeff_x       = zeros(eltype(grid), M - posx[1])
    coeff_y       = zeros(eltype(grid), M - posy[1])
    coeff_z       = zeros(eltype(grid), M - posz[1])
    coeff_bound_x = zeros(eltype(grid), M - posx[2])
    coeff_bound_y = zeros(eltype(grid), M - posy[2])
    coeff_bound_z = zeros(eltype(grid), M - posz[2])

    # initializing elements which vary during the simulation (as a function of Δt)
    loop! = _initialize_variable_diagonal!(Architectures.device(arch), heuristic_workgroup(N...), N)
    event = loop!(diag, D, N; dependencies=Event(Architectures.device(arch)))
    wait(event)

    # filling elements which stay constant in time
    fill_core_matrix!(coeff_d , coeff_x, coeff_y, coeff_z, Ax, Ay, Az, C, N, dims)
    if dims[1]
        fill_boundaries_x!(coeff_d, coeff_bound_x, Ax, N, topo[1])
    end
    if dims[2]
        fill_boundaries_y!(coeff_d, coeff_bound_y, Ay, N, topo[2])
    end
    if dims[3]
        fill_boundaries_z!(coeff_d, coeff_bound_z, Az, N, topo[3])
    end

    sparse_matrix = spdiagm(0=>coeff_d,
                      posx[1]=>coeff_x,      -posx[1]=>coeff_x,
                      posx[2]=>coeff_bound_x,-posx[2]=>coeff_bound_x,
                      posy[1]=>coeff_y,      -posy[1]=>coeff_y,
                      posy[2]=>coeff_bound_y,-posy[2]=>coeff_bound_y,
                      posz[1]=>coeff_z,      -posz[1]=>coeff_z,
                      posz[2]=>coeff_bound_z,-posz[2]=>coeff_bound_z)

    ensure_diagonal_elements_are_present!(sparse_matrix)

    matrix_constructors = constructors(arch, sparse_matrix)

    return matrix_constructors, diag, N
end

@kernel function _initialize_variable_diagonal!(diag, D, N)
    i, j, k = @index(Global, NTuple)
    t  = i + N[1] * (j - 1 + N[2] * (k - 1))
    diag[t] = D[i, j, k]
end

function fill_core_matrix!(coeff_d, coeff_x, coeff_y, coeff_z, Ax, Ay, Az, C, N, dims)
    Nx, Ny, Nz = N
    for k = 1:Nz, j = 1:Ny, i = 1:Nx
        t          = i +  Nx * (j - 1 + Ny * (k - 1))
        coeff_d[t] = C[i, j, k]
    end
    if dims[1]
        for k = 1:Nz, j = 1:Ny, i = 1:Nx-1
            t             = i +  Nx * (j - 1 + Ny * (k - 1))
            coeff_x[t]    = Ax[i+1, j, k]
            coeff_d[t]   -= coeff_x[t]
            coeff_d[t+1] -= coeff_x[t]
        end
    end
    if dims[2]
        for k = 1:Nz, j = 1:Ny-1, i = 1:Nx
            t              = i +  Nx * (j - 1 + Ny * (k - 1))
            coeff_y[t]     = Ay[i, j+1, k]
            coeff_d[t]    -= coeff_y[t]
            coeff_d[t+Nx] -= coeff_y[t]
        end
    end
    if dims[3]
        for k = 1:Nz-1, j = 1:Ny, i = 1:Nx
            t                 = i +  Nx * (j - 1 + Ny * (k - 1))
            coeff_z[t]        = Az[i, j, k+1]
            coeff_d[t]       -= coeff_z[t]
            coeff_d[t+Nx*Ny] -= coeff_z[t]
        end
    end
end

# No-flux boundary conditions are implied in the construction of the matrix, so only
# periodic boundary conditions have to be filled in
#
# In case of a periodic boundary condition we have to modify the diagonal
# as well as the off-diagonals. As an example, for x-periodic boundary conditions
# we have to modify the following diagonal elements
#
# row number (1  + Nx * (j - 1 + Ny * (k - 1))) => corresponding to i = 1 , j = j and k = k
# row number (Nx + Nx * (j - 1 + Ny * (k - 1))) => corresponding to i = Nx, j = j and k = k
#
# Since zero-flux BC were implied, we have to also have to add the coefficients corresponding to i-1 and i+1
# (respectively). Since the off-diagonal elements are symmetric we can fill it in only once

 @inline fill_boundaries_x!(coeff_d, coeff_bound_x, Ax, N, ::Type{Bounded}) = nothing
 @inline fill_boundaries_x!(coeff_d, coeff_bound_x, Ax, N, ::Type{Flat})    = nothing
function fill_boundaries_x!(coeff_d, coeff_bound_x, Ax, N, ::Type{Periodic})
    Nx, Ny, Nz = N
    for k = 1:Nz, j = 1:Ny
        tₘ = 1  + Nx * (j - 1 + Ny * (k - 1))
        tₚ = Nx + Nx * (j - 1 + Ny * (k - 1))
        coeff_bound_x[tₘ] = Ax[1, j, k]
        coeff_d[tₘ]      -= coeff_bound_x[tₘ]
        coeff_d[tₚ]      -= coeff_bound_x[tₘ]
    end
end

 @inline fill_boundaries_y!(coeff_d, coeff_bound_y, Ay, N, ::Type{Bounded}) = nothing
 @inline fill_boundaries_y!(coeff_d, coeff_bound_y, Ay, N, ::Type{Flat})    = nothing
function fill_boundaries_y!(coeff_d, coeff_bound_y, Ay, N, ::Type{Periodic})
    Nx, Ny, Nz = N

    for k = 1:Nz, i = 1:Nx
        tₘ = i + Nx * (1 - 1 + Ny * (k - 1))
        tₚ = i + Nx * (Ny- 1 + Ny * (k - 1))
        coeff_bound_y[tₘ] = Ay[i, 1, k]
        coeff_d[tₘ]      -= coeff_bound_y[tₘ]
        coeff_d[tₚ]      -= coeff_bound_y[tₘ]
    end
end

 @inline fill_boundaries_z!(coeff_d, coeff_bound_z, Az, N, ::Type{Bounded}) = nothing
 @inline fill_boundaries_z!(coeff_d, coeff_bound_z, Az, N, ::Type{Flat})    = nothing
function fill_boundaries_z!(coeff_d, coeff_bound_z, Az, N, ::Type{Periodic})
    Nx, Ny, Nz = N
    for j = 1:Ny, i = 1:Nx
        tₘ = i + Nx * (j - 1 + Ny * (1 - 1))
        tₚ = i + Nx * (j - 1 + Ny * (Nz- 1))
        coeff_bound_z[tₘ] = Az[i, j, 1]
        coeff_d[tₘ]      -= coeff_bound_z[tₘ]
        coeff_d[tₚ]      -= coeff_bound_z[tₘ]
    end
end

function solve!(x, solver::HeptadiagonalIterativeSolver, b, Δt)

    arch = architecture(solver.matrix)

    # update matrix and preconditioner if time step changes
    if Δt != solver.previous_Δt
        constructors = deepcopy(solver.matrix_constructors)
        update_diag!(constructors, arch, solver.problem_size, solver.diagonal, Δt)
        solver.matrix = arch_sparse_matrix(arch, constructors)

        solver.preconditioner = build_preconditioner(Val(solver.preconditioner_method),
                                                     solver.matrix,
                                                     solver.preconditioner_settings)

        solver.previous_Δt = Δt
    end

    #q = solver.iterative_solver(solver.matrix, b, maxiter=solver.maximum_iterations, reltol=solver.tolerance, Pl=solver.preconditioner, verbose=true)
    #q = solver.iterative_solver(solver.matrix, b, maxiter=solver.maximum_iterations, reltol=solver.tolerance, verbose=true)
    q = solver.iterative_solver(solver.matrix, b, maxiter=solver.maximum_iterations, reltol=solver.tolerance, Pl=solver.preconditioner)

    set!(x, reshape(q, solver.problem_size...))
    fill_halo_regions!(x)

    return
end

function Base.show(io::IO, solver::HeptadiagonalIterativeSolver)
    print(io, "Matrix-based iterative solver with: \n")
    print(io, "├── Problem size: "  , solver.problem_size, '\n')
    print(io, "├── Grid: "  , solver.grid, '\n')
    print(io, "├── Solution method: ", solver.iterative_solver, '\n')
    print(io, "└── Preconditioner: ", solver.preconditioner_method)
    return nothing
end


text:
using Pkg, BinDeps, PyCall
using PyCall: python, conda, Conda
using BinDeps: generate_steps, getallproviders, lower, PackageManager

!Sys.islinux() && exit()

if isnothing(Sys.which("sudo")) # in docker
    try run(`apt update`) catch end
    try run(`yum update`) catch end
end

@BinDeps.setup

wget = library_dependency("wget")
git = library_dependency("git")
make = library_dependency("make")
cmake = library_dependency("cmake")
gcc = library_dependency("gcc")
jre = library_dependency("jre")

common = Dict("wget" => wget, "git" => git, "make" => make)
provides(AptGet, Dict(common..., "g++" => gcc, "default-jre" => jre, "cmake" => cmake))
provides(Yum, Dict(common..., "gcc-c++" => gcc,  "java-1.8.0-openjdk-headless" => jre, "cmake3" => cmake))

for dep in bindeps_context.deps
    dp, opts = getallproviders(dep, PackageManager)[1]
    cmd = lower(generate_steps(dep, dp, opts)).steps[1]
    i = findfirst(x -> x == "install", cmd.exec)
    insert!(cmd.exec, i + 1, "-y")
    println(cmd)
    try run(cmd) catch end
end

for pkg in ["scikit-learn", "pandas", "h2o", "lightgbm",  "xgboost", "unidecode", "pyarrow",
            "thundersvm", "tensorflow-ranking", "tensorflow", "treelite", "treelite_runtime"]
    # "catboost==0.16.5"
    try
        run(`$python -m pip install $pkg`)
    catch e
        println(e)
    end
end

buildsh = joinpath(@__DIR__, "build.sh")
ENV["PYTHON"] = python
ENV["JULIA_DEPOT_PATH"] = DEPOT_PATH[1]
run(`bash $buildsh`)


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

JLLWrappers.@generate_main_file_header("EarCut")
JLLWrappers.@generate_main_file("EarCut", UUID("5ae413db-bbd1-5e63-b57d-d24a61df00f5"))
end  # module EarCut_jll


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# only returns real roots
function quadraticroots(a::T, b::T, c::T) where {T<:Real}
    temp = b^2 - 4 * a * c

    if temp < zero(T)
        return nothing # no real roots so no ray cylinder intersection
    end

    radical = sqrt(temp)
    if b >= zero(T)
        radicalterm = -b - radical
        x1 = radicalterm / (2a)
        x2 = (2c) / radicalterm
    else
        radicalterm = -b + radical
        x1 = (2c) / radicalterm
        x2 = radicalterm / (2a)
    end

    return SVector{2,T}(x1, x2)
end

replprint(a) = show(IOContext(stdout), "text/plain", a)

@inline samepoint(pt1, pt2) = isapprox(pt1, pt2, atol = 1e-10)

password() = String(collect(rand(('1':'9'..., 'A':'Z'..., 'a':'z'...)) for i in 1:14))

# 1 if val>= 0 and -1 otherwise. Unlike Math.Sign which returns 0 if val==0

@inline sign(val::T) where {T<:Real} = val >= zero(T) ? 1 : -1

@inline function NaNsafeatan(x::T, y::T)::T where {T<:Real}
    if y == zero(T)
        if x == zero(T)
            return zero(T)
        elseif x > zero(T)
            return T(π / 2)
        else
            return -T(π / 2)
        end
    else
        return atan(x, y)
    end
end

@inline function NaNsafeasin(x::T)::T where {T<:Real}
    if x == zero(T)
        return zero(T)
    elseif x == one(T)
        return T(π / 2)
    elseif x == -one(T)
        return -T(π / 2)
    else
        return asin(x)
    end
end

@inline function NaNsafeacos(x::T)::T where {T<:Real}
    if x == zero(T)
        return T(π / 2)
    elseif x == one(T)
        return zero(T)
    elseif x == -one(T)
        return T(π)
    else
        return acos(x)
    end
end


text:
###################
# address schemas #
###################

abstract type AddressSchema end

struct StaticAddressSchema <: AddressSchema
    keys::Set{Symbol}
end

Base.keys(schema::StaticAddressSchema) = schema.keys

struct VectorAddressSchema <: AddressSchema end
struct SingleDynamicKeyAddressSchema <: AddressSchema end
struct DynamicAddressSchema <: AddressSchema end
struct EmptyAddressSchema <: AddressSchema end
struct AllAddressSchema <: AddressSchema end

export AddressSchema
export StaticAddressSchema # hierarchical
export VectorAddressSchema # hierarchical
export SingleDynamicKeyAddressSchema # hierarchical
export DynamicAddressSchema # hierarchical
export EmptyAddressSchema
export AllAddressSchema

######################
# abstract selection #
######################

"""
    abstract type Selection end

Abstract type for selections of addresses.

All selections implement the following methods:

    Base.in(addr, selection)

Is the address selected?

    Base.getindex(selection, addr)

Get the subselection at the given address.

    Base.isempty(selection)

Is the selection guaranteed to be empty?

    get_address_schema(T)

Return a shallow, compile-time address schema, where `T` is the concrete type of the selection.
"""
abstract type Selection end

Base.in(addr, ::Selection) = false
Base.getindex(::Selection, addr) = EmptySelection()

export Selection

##########################
# hierarchical selection #
##########################

"""
    abstract type HierarchicalSelection <: Selection end

Abstract type for selections that have a notion of sub-selections.

    get_subselections(selection::HierarchicalSelection)

Return an iterator over pairs of addresses and subselections at associated addresses.
"""
abstract type HierarchicalSelection <: Selection end

export HierarchicalSelection
export get_subselections

###################
# empty selection #
###################

"""
    struct EmptySelection <: Selection end

A singleton type for a selection that is always empty.
"""
struct EmptySelection <: Selection end
get_address_schema(::Type{EmptySelection}) = EmptyAddressSchema()
Base.isempty(::EmptySelection) = true

export EmptySelection

#################
# all selection #
#################

"""
    struct AllSelection <: Selection end

A singleton type for a selection that contains all choices at or under an address.
"""
struct AllSelection <: Selection end
get_address_schema(::Type{AllSelection}) = AllAddressSchema()
Base.isempty(::AllSelection) = false # it is not guaranteed to be empty
Base.in(addr, ::AllSelection) = true
Base.getindex(::AllSelection, addr) = AllSelection()

export AllSelection

########################
# complement selection #
########################

struct ComplementSelection <: Selection
    complement::Selection
end
get_address_schema(::Type{ComplementSelection}) = DynamicAddressSchema()
Base.isempty(::ComplementSelection) = false # it is not guaranteed to be empty
Base.in(addr, selection::ComplementSelection) = !(addr in selection.complement)
function Base.getindex(selection::ComplementSelection, addr)
    ComplementSelection(selection.complement[addr])
end

"""
    comp_selection = complement(selection::Selection)

Return a selection that is the complement of the given selection.

An address is in the selection if it is not in the complement selection.
"""
function complement(selection::Selection)
    ComplementSelection(selection)
end

export ComplementSelection, complement

####################
# static selection #
####################

# R is a tuple of symbols..
# T is a tuple of symbols
# U the tuple type of subselections

"""
    struct StaticSelection{T,U} <: HierarchicalSelection .. end

A hierarchical selection whose keys are among its type parameters.
"""
struct StaticSelection{T,U} <: HierarchicalSelection
    subselections::NamedTuple{T,U}
end

function Base.isempty(selection::StaticSelection{T,U}) where {T,U}
    length(R) == 0 && all(isempty(node) for node in selection.subselections)
end

function get_address_schema(::Type{StaticSelection{T,U}}) where {T,U}
    keys = Set{Symbol}()
    for (key, _) in zip(T, U.parameters)
        push!(keys, key)
    end
    StaticAddressSchema(keys)
end

get_subselections(selection::StaticSelection) = pairs(selection.subselections)

function static_getindex(selection::StaticSelection, ::Val{A}) where {A}
    selection.subselections[A]
end

# TODO do we no longer need static_in?

function Base.getindex(selection::StaticSelection, addr::Symbol)
    if haskey(selection.subselections, addr)
        selection.subselections[addr]
    else
        EmptySelection()
    end
end

function Base.getindex(selection::StaticSelection, addr::Pair)
    (first, rest) = addr
    subselection = selection.subselections[first]
    subselection[rest]
end

function Base.in(addr::Symbol, selection::StaticSelection{T,U}) where {T,U}
    addr in T && selection.subselections[addr] == AllSelection()
end

function Base.in(addr::Pair, selection::StaticSelection)
    (first, rest) = addr
    if haskey(selection.subselections, first)
        subselection = selection.subselections[first]
        in(subselection, rest)
    else
        false
    end
end

function StaticSelection(other::HierarchicalSelection)
    keys_and_subselections = collect(get_subselections(other))
    if length(keys_and_subselections) > 0
        (keys, subselections) = collect(zip(keys_and_subselections...))
    else
        (keys, subselections) = ((), ())
    end
    types = map(typeof, subselections)
    StaticSelection{keys,Tuple{types...}}(NamedTuple{keys}(subselections))
end

export StaticSelection

#####################
# dynamic selection #
#####################

"""
    struct DynamicSelection <: HierarchicalSelection .. end

A hierarchical, mutable, selection with arbitrary addresses.

Can be mutated with the following methods:

    Base.push!(selection::DynamicSelection, addr)

Add the address and all of its sub-addresses to the selection.

Example:
```julia
selection = select()
@assert !(:x in selection)
push!(selection, :x)
@assert :x in selection
```

    set_subselection!(selection::DynamicSelection, addr, other::Selection)

Change the selection status of the given address and its sub-addresses that defined by `other`.

Example:
```julia
selection = select(:x)
@assert :x in selection
subselection = select(:y)
set_subselection!(selection, :x, subselection)
@assert (:x => :y) in selection
@assert !(:x in selection)
```

Note that `set_subselection!` does not copy data in `other`, so `other` may be mutated by a later calls to `set_subselection!` for addresses under `addr`.
"""
struct DynamicSelection <: HierarchicalSelection
    # note: only store subselections for which isempty = false
    subselections::Dict{Any,Selection}
end

function Base.isempty(selection::DynamicSelection)
    isempty(selection.subselections)
end

DynamicSelection() = DynamicSelection(Dict{Any,Selection}())

get_address_schema(::Type{DynamicSelection}) = DynamicAddressSchema()

function Base.in(addr, selection::DynamicSelection)
    if haskey(selection.subselections, addr)
        selection.subselections[addr] == AllSelection()
    else
        false
    end
end

function Base.in(addr::Pair, selection::DynamicSelection)
    (first, rest) = addr
    if haskey(selection.subselections, first)
        subselection = selection.subselections[first]
        @assert !isempty(subselection)
        rest in subselection
    else
        false
    end
end

function Base.getindex(selection::DynamicSelection, addr)
    if haskey(selection.subselections, addr)
        selection.subselections[addr]
    else
        EmptySelection()
    end
end

function Base.getindex(selection::DynamicSelection, addr::Pair)
    (first, rest) = addr
    if haskey(selection.subselections, first)
        subselection = selection.subselections[first]
        @assert !isempty(subselection)
        getindex(subselection, rest)
    else
        EmptySelection()
    end
end

function Base.push!(selection::DynamicSelection, addr)
    selection.subselections[addr] = AllSelection()
end

function Base.push!(selection::DynamicSelection, addr::Pair)
    (first, rest) = addr
    if haskey(selection.subselections, first)
        subselection = selection.subselections[first]
    else
        subselection = DynamicSelection()
        selection.subselections[first] = subselection
    end
    push!(subselection, rest)
end

function set_subselection!(selection::DynamicSelection, addr, other::Selection)
    selection.subselections[addr] = other
end

function set_subselection!(selection::DynamicSelection, addr::Pair, other::Selection)
    (first, rest) = addr
    if haskey(selection.subselections, first)
        subselection = selection.subselections[first]
    else
        subselection = DynamicSelection()
        selection.subselections[first] = subselection
    end
    set_subselection!(subselection, rest, other)
end

get_subselections(selection::DynamicSelection) = selection.subselections

"""
    selection = select(addrs...)

Return a selection containing a given set of addresses.

Examples:
```julia
selection = select(:x, "foo", :y => 1 => :z)
selection = select()
selection = select(:x => 1, :x => 2)
```
"""
function select(addrs...)
    selection = DynamicSelection()
    for addr in addrs
        push!(selection, addr)
    end
    selection
end

"""
    selection = selectall()

Construct a selection that includes all random choices.
"""
function selectall()
    AllSelection()
end

export DynamicSelection
export select, selectall, set_subselection!

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
import Base: merge, hcat, vcat, map

###### merge ####################

function _merge_outer(::Type{IndexType}, ta1::TimeArray{T,N,D}, ta2::TimeArray{T,M,D}, padvalue, meta) where {IndexType,T,N,M,D}
    timestamps, new_idx1, new_idx2 = sorted_unique_merge(IndexType, timestamp(ta1), timestamp(ta2))
    vals = fill(convert(T, padvalue), (length(timestamps), length(colnames(ta1)) + length(colnames(ta2))))
    insertbyidx!(vals, values(ta1), new_idx1)
    insertbyidx!(vals, values(ta2), new_idx2, size(values(ta1), 2))
    TimeArray(timestamps, vals, [colnames(ta1); colnames(ta2)], meta; unchecked = true)
end

"""
    merge(ta1::TimeArray{T}, ta2::TimeArray{T}, [tas::TimeArray{T}...];
          method = :inner, colnames = [...], padvalue = NaN)

Merge several `TimeArray`s along with the time index.

## Argument

- `method::Symbol`: `:inner`, `:outer`, `:left` or `:right`.
"""
function merge(ta1::TimeArray{T,N,D}, ta2::TimeArray{T,M,D};
               method::Symbol = :inner, colnames::Vector = Symbol[], meta = nothing,
               padvalue=NaN) where {T,N,M,D}

    if colnames isa Vector{<:AbstractString}
        @warn "`merge(...; colname::Vector{<:AbstractString})` is deprecated, " *
              "use `merge(...; colnames=Symbol.(colnames))` instead."
        colnames = Symbol.(colnames)
    end

    meta = if _meta(ta1) == _meta(ta2) && meta ≡ nothing
        _meta(ta1)
    elseif typeof(_meta(ta1)) <: AbstractString && typeof(_meta(ta2)) <: AbstractString && meta ≡ nothing
        string(_meta(ta1), "_", _meta(ta2))
    else
        meta
    end

    if method == :inner

        idx1, idx2 = overlap(timestamp(ta1), timestamp(ta2))
        vals = [values(ta1[idx1]) values(ta2[idx2])]
        ta = TimeArray(timestamp(ta1[idx1]), vals, [_colnames(ta1); _colnames(ta2)], meta; unchecked = true)

    elseif method == :left

        new_idx2, old_idx2 = overlap(timestamp(ta1), timestamp(ta2))
        right_vals = fill(convert(T, padvalue), (length(ta1), length(_colnames(ta2))))
        insertbyidx!(right_vals, values(ta2), new_idx2, old_idx2)
        ta = TimeArray(timestamp(ta1), [values(ta1) right_vals], [_colnames(ta1); _colnames(ta2)], meta; unchecked = true)

    elseif method == :right

        ta = merge(ta2, ta1, method = :left; padvalue = padvalue)
        ncol2 = length(_colnames(ta2))
        vals = [values(ta)[:, (ncol2+1):end] values(ta)[:, 1:ncol2]]
        ta = TimeArray(timestamp(ta), vals, [_colnames(ta1); _colnames(ta2)], meta; unchecked = true)

    elseif method == :outer

        ta = if (length(timestamp(ta1)) + length(timestamp(ta2))) > typemax(Int32)
            _merge_outer(Int64, ta1, ta2, padvalue, meta)
        else
            _merge_outer(Int32, ta1, ta2, padvalue, meta)
        end

    else
        throw(ArgumentError(
            "merge method must be one of :inner, :left, :right, :outer"))
    end

    return rename!(ta, colnames)

end

merge(x::TimeArray{T}, y::TimeArray{T}, z::TimeArray{T}, a::TimeArray{T}...;
      colnames::Vector = Symbol[], kw...) where {T} =
    rename!(merge(merge(x, y; kw...), z, a...; kw...), colnames)

# hcat ##########################

function hcat(x::TimeArray, y::TimeArray)
    tsx = timestamp(x)
    tsy = timestamp(y)

    if length(tsx) != length(tsx) || tsx != tsy
        throw(DimensionMismatch(
            "timestamps not consistent, please checkout `merge`."))
    end

    m = ifelse(meta(x) == meta(y), meta(x), nothing)

    TimeArray(tsx, [values(x) values(y)], [colnames(x); colnames(y)], m)
end

hcat(x::TimeArray, y::TimeArray, zs::Vararg{TimeArray}) =
    hcat(hcat(x, y), zs...)

# collapse ######################

# accessors functions
# https://github.com/JuliaLang/julia/blob/b617e8d3a77e49cd5625ca663af88e40f7796f15/stdlib/Dates/src/accessors.jl#L50

# `year` doesn't need this wrapper
for (F, T, P) ∈ ((:quarter,     :TimeType,               :Quarter),
                 (:month,       :TimeType,               :Month),
                 (:week,        :TimeType,               :Week),
                 (:day,         :TimeType,               :Day),
                 (:hour,        :(Union{DateTime,Time}), :Hour),
                 (:minute,      :(Union{DateTime,Time}), :Minute),
                 (:second,      :(Union{DateTime,Time}), :Second),
                 (:millisecond, :(Union{DateTime,Time}), :Millisecond),
                 (:microsecond, :Time,                   :Microsecond),
                 (:nanosecond,  :Time,                   :Nanosecond))
    if F === :quarter
        (VERSION < v"1.6") && continue
        F = :(Dates.quarter)
    end
    @eval let
        global collapse
        wrapper(x::$T) = floor(x, $P(1))
        collapse(ta::TimeArray, ::typeof($F), f::Function, g::Function = f; kw...) =
            collapse(ta, wrapper, f, g; kw...)
    end
end

function collapse(ta::TimeArray, period::Function, timestamp::Function,
                  value::Function = timestamp)

    isempty(ta) && return ta

    m, n = length(ta), length(colnames(ta))
    ts   = _timestamp(ta)
    val  = _values(ta)
    idx  = UnitRange{Int}[]
    sizehint!(idx, m)

    t₀ = period(ts[1])
    j  = 1
    for i in 1:m-1
        t₁ = period(ts[i+1])
        t₀ == t₁ && continue
        push!(idx, j:i)
        j = i + 1
        t₀ = t₁
    end
    push!(idx, j:m)

    ts′  = [timestamp(@view(ts[i])) for i ∈ idx]
    val′ = if n == 1
        [value(@view(val[i])) for i ∈ idx]
    else
        [value(@view(val[i, k])) for i ∈ idx, k ∈ 1:n]
    end

    TimeArray(ts′, val′, colnames(ta), meta(ta))
end

collapse(ta::TimeArray, p::Period, t::Function, v::Function = t; kw...) =
    collapse(ta, x -> floor(x, p), t, v; kw...)

# vcat ######################

"""
    $(SIGNATURES)

Concatenate two ``TimeArray`` into single object.

If there are duplicated timestamps, we will keep order as the function input.

```julia-repl
julia> a = TimeArray([Date(2015, 10, 1), Date(2015, 10, 2), Date(2015, 10, 3)], [1, 2, 3]);

julia> b = TimeArray([Date(2015, 10, 2), Date(2015, 10, 3)], [4, 5]);

julia> [a; b]
5×1 TimeArray{Int64,1,Date,Array{Int64,1}} 2015-10-01 to 2015-10-03
│            │ A     │
├────────────┼───────┤
│ 2015-10-01 │ 1     │
│ 2015-10-02 │ 2     │
│ 2015-10-02 │ 4     │
│ 2015-10-03 │ 3     │
│ 2015-10-03 │ 5     │
```
"""
function vcat(tas::TimeArray...)
    # Check all meta fields are identical.
    prev_meta = meta(tas[1])
    for ta in tas
        if meta(ta) != prev_meta
            throw(ArgumentError("metadata doesn't match"))
        end
    end

    # Check column names are identical.
    prev_colnames = colnames(tas[1])
    for ta in tas
        if colnames(ta) != prev_colnames
            throw(ArgumentError("column names don't match"))
        end
    end

    # Concatenate the contents.
    ts  = vcat([timestamp(ta) for ta in tas]...)
    val = vcat([values(ta) for ta in tas]...)

    order = sortperm(ts)
    if ndims(tas[1]) == 1 # Check for 1D to ensure values remains a 1D vector.
        return TimeArray(ts[order], val[order], prev_colnames, prev_meta)
    else
        return TimeArray(ts[order], val[order, :], prev_colnames, prev_meta)
    end
end

# map ######################

@generated function map(f, ta::TimeArray{T,N}) where {T,N}
    input_val  = (N == 1) ? :(values(ta)[i]) : :(vec(values(ta)[i, :]))
    output_val = (N == 1) ? :(vals[i]) : :(vals[i, :])

    output_vals = (N == 1) ? :(vals[order]) : :(vals[order, :])

    quote
        ts   = similar(timestamp(ta))
        vals = similar(values(ta))

        for i in eachindex(ta)
            @inbounds ts[i], $output_val = f(timestamp(ta)[i], $input_val)
        end

        order = sortperm(ts)
        TimeArray(ts[order], $output_vals, colnames(ta), meta(ta))
    end
end


text:
struct BenchmarkResult{TA, TB, TC}
    df::DataFrames.DataFrame
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -20.345339284625858 <= q <= 16.479136897808104)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.24719801624822457 + 0.29438950671798914*1/(1 + exp(-(0.7291821855933653 + 0.701331929085816*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.07780867500186073*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.19203457942980284*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.7065670431246684*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + -0.7293079760660532*1/(1 + exp(-(0.05096633991424904 + 0.31501139336769945*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + -0.3152903616610754*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.8990116916055966*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.8373103942919853*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.8601808654748551*1/(1 + exp(-(0.18255055897839467 + -0.7759686413687503*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.7813463540455929*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + 0.215367127826533*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.9673448256726829*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.9720568561360068*1/(1 + exp(-(0.5199334546165302 + 0.4081309459110076*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.6995024782698374*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.3131098472642604*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.4135863095985335*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))))))))) + 1/(1 + exp(-(-0.49402133773229506 + 0.2465550492028581*1/(1 + exp(-(0.7291821855933653 + 0.701331929085816*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.07780867500186073*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.19203457942980284*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.7065670431246684*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.17626257051164496*1/(1 + exp(-(0.05096633991424904 + 0.31501139336769945*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + -0.3152903616610754*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.8990116916055966*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.8373103942919853*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.8688281780619751*1/(1 + exp(-(0.18255055897839467 + -0.7759686413687503*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.7813463540455929*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + 0.215367127826533*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.9673448256726829*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.3904178182819744*1/(1 + exp(-(0.5199334546165302 + 0.4081309459110076*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.6995024782698374*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.3131098472642604*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.4135863095985335*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))))))))) + 1/(1 + exp(-(-0.6759660684871824 + 0.9311431633625067*1/(1 + exp(-(0.7291821855933653 + 0.701331929085816*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.07780867500186073*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.19203457942980284*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.7065670431246684*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + -0.2015839225926066*1/(1 + exp(-(0.05096633991424904 + 0.31501139336769945*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + -0.3152903616610754*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.8990116916055966*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.8373103942919853*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.01676614484829253*1/(1 + exp(-(0.18255055897839467 + -0.7759686413687503*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.7813463540455929*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + 0.215367127826533*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.9673448256726829*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + -0.07861065087286123*1/(1 + exp(-(0.5199334546165302 + 0.4081309459110076*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.6995024782698374*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.3131098472642604*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.4135863095985335*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))))))))) + 1/(1 + exp(-(-0.08101830092609674 + 0.5708385401734484*1/(1 + exp(-(0.7291821855933653 + 0.701331929085816*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.07780867500186073*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.19203457942980284*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.7065670431246684*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + -0.9629021826429458*1/(1 + exp(-(0.05096633991424904 + 0.31501139336769945*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + -0.3152903616610754*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.8990116916055966*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.8373103942919853*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + 0.4780289946412233*1/(1 + exp(-(0.18255055897839467 + -0.7759686413687503*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.7813463540455929*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + 0.215367127826533*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + 0.9673448256726829*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))))))) + -0.25224290978348307*1/(1 + exp(-(0.5199334546165302 + 0.4081309459110076*1/(1 + exp(-(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])))) + 0.6995024782698374*1/(1 + exp(-(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])))) + -0.3131098472642604*1/(1 + exp(-(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])))) + -0.4135863095985335*1/(1 + exp(-(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using NPZ
using JuMP
using Ipopt
using LinearAlgebra
using WoodburyMatrices
using .Threads
using ProgressMeter
import Base: isless

struct TestResult
    i::Int64
    f::Float64
end
isless(a::TestResult, b::TestResult) = isless(a.f, b.f)

function main()
    T = Float64
    X1 = npzread("Kernels/G1.npy")
    X2 = npzread("Kernels/G2.npy")
    X3 = npzread("Kernels/G4.npy")
    nfactor = npzread("Testing/nfactor.npy")

    ϵ = 1e-2
    λ = 1.0
    lnorm = -log(ϵ)*(2*size(X1,2))

    X = [X3 zeros(T, size(X3));
         zeros(T, size(X3)) X3;
         nfactor*X1 nfactor*X2]
    Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    X ./= Xn

    n = size(X3, 1)
    ndas = size(X1, 1)

    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        l = μ.^2 .*diag(X*(Mc\X'))
        return (M, Mc, l)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, λ, lnorm, n)
        M, Mc, l = Ml(μ)
        l2 = l[1:n] .+ l[(n+1):2*n] # coherence of point sensors only; sum across channels
        return -logdet(M) + λ*maximum(l2) - lnorm
    end

    f(μ) = f(μ, λ, lnorm, n)

    function testf(u, i)
        u_test = copy(u)
        u_test[i] = 1.0
        u_test[i+64] = 1.0
        return TestResult(i, f(u_test))
    end

    u = vcat(zeros(64+64), ones(size(X1,1)))
    to_test = Set(1:64)
    selection_order = Int[]
    f_vals = T[]
    p = Progress(32*65, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

    while length(to_test) > 0
        test_vec = TestResult[]
        for i in collect(to_test)
            push!(test_vec, testf(u, i))
            next!(p)
        end
        best = minimum(test_vec)
        delete!(to_test, best.i)
        push!(selection_order, best.i)
        push!(f_vals, best.f)
        u[best.i] = 1.0
        u[best.i+64] = 1.0
    end

    npzwrite("IncoherenceResults/seq_order.npy", selection_order)
    npzwrite("IncoherenceResults/seq_vals.npy", f_vals)

end

main()
# npzwrite("res_mixed.npy", value.(x))
# ;

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
using Documenter, ProgressLogging

makedocs(;
    modules=[ProgressLogging],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/ProgressLogging.jl/blob/{commit}{path}#L{line}",
    sitename="ProgressLogging.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/tkf/ProgressLogging.jl",
)


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    testField(x,y) = Point(-x, 2y) # x'(t) = -x, y'(t) = 2y
    x =  -2:0.1:4
    y =  -2:0.1:2
    u2(x,y) = -x
    v2(x,y) = 2y
    z = [log10(sqrt(u2(x,y)^2 + v2(x,y)^2)) for x in x, y in y]

    fig = Figure(resolution = (700, 450), fontsize = 18, font = "sans")
    ax = fig[1, 1] = Axis(fig, xlabel = "x", ylabel = "y")
    fs = heatmap!(ax, x, y, z, colormap = Reverse(:plasma))
    streamplot!(ax, testField, x, y, colormap = Reverse(:plasma),
        gridsize= (32,32), arrow_size = 10)
    #limits!(ax, -2, 4, -2, 2)
    cbar = fig[1,2] = Colorbar(fig, fs, label = "log10[(u²+v²)½]",
        labelsize = 14, ticklabelsize=14)
    cbar.width = 20
    colgap!(fig.layout, 5)
    fig
    save(joinpath(@__DIR__, "output", "streamplotAndHeatmap.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# Useful functions for checking proof details of the impossibility proof

"""
Builds the B matrix from the proof
"""
function Bmat(b1,b2)
    k = length(b1)
    @assert(isodd(k))
    r = round(Int64,(k+1)/2)
    b_half = [reverse(b2[r:end]); b1[r:end]]
    Db = diagm(b_half)
    Dk = zeros(k+1,k+1)
    for i = 1:r
        Dk[i,i] = k-i+1
        Dk[k-i+2,k-i+2] = k-i+1
    end
    row = collect(k:-1:0.0)
    r2 = collect(0.0:1:k)
    Top = repeat(row', r)
    Bottom = repeat(r2',r)
    R = [Top; Bottom]
    # @show size(Dk), size(Db), size(R), r, k

    B = Dk - Db*R

    return Dk, Db, R, B

end

"""
Computes primal variables corresponding to a complete hypergraph in the
case of using standard baseline scores for the proof.

These will often be optimal, but not if there is significant imbalance between
the two class sizes
"""
function Xstanard(k,n1,n2)
    n = n1+n2
    Xs = zeros(k+1)
    for i = 1:k+1
        t = i-1
        Xs[i] = binomial(n1,t)*binomial(n2,k-t)/binomial(n,k)
    end

    return Xs
end

"""
Optimal primal solution for but for asymptotic baselines
"""
function Xasym(k,alph)
    Xa = zeros(k+1)
    for i = 1:k+1
        t = i-1
        Xa[i] = alph^t*(1-alph)^(k-t)*binomial(k,t)
    end

    return Xa
end

"""
Use Gurobi to get the solution numerically
"""
function PrimalDualXY(gurobi_env,B,k)

    A = [B -ones(k+1,1)]
    p = [ones(1,k+1) 0]
    m = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(gurobi_env),
        "OutputFlag" => 0, "FeasibilityTol"=>1e-9,"OptimalityTol"=>1e-9))
    @variable(m, x[1:k+2])
    @constraint(m, y, A*x .>= zeros(k+1,1))
    p = [ones(1,k+1) 0]
    @constraint(m, a, p*x .== 1)
    @constraint(m, z, x .>= zeros(k+2,1))

    @objective(m, Max, x[k+2])
    JuMP.optimize!(m)

    X = JuMP.value.(x)
    Y= JuMP.dual.(y)
    gam = X[end]

    return X[1:end-1], Y, gam
end

"""
Get the optimal dual variables for balanced classes.
"""
function BalancedDual(k)
    L = zeros(k+1)
    r = round(Int64,(k+1)/2)
    for i = 1:r
        L[i] = k/(k-i+1)
        L[k-i+2] = k/(k-i+1)
    end
    delta = sum(L)
    Ybal = L/delta
    return Ybal, delta
end

"""
Get optimal dual variables for any class size.
"""
function OptimalDual(b1,b2)
    k = length(b1)
    @assert(isodd(k))
    Ybal, delta = BalancedDual(k)
    r = round(Int64,(k+1)/2)
    s1 = 0
    s2 = 0
    for i = 1:r
        t = k-i+1
        s1 += b2[t]*(k/t - 1)
        s2 += (2-k/t)*b2[t]
    end

    y2k = 2/delta*s1/(1-s2)
    Ydual = zeros(2*r)
    Ydual[1] = y2k
    Ydual[k+1] = 2/delta - y2k

    for i = 2:r
        t = k-i+1
        Ydual[i] = 2/delta*(k/t-1) + (2-k/t)*y2k
        Ydual[k-i+2] = 2*k/(delta*t) - Ydual[i]
    end
    return Ydual
end

"""
For conceptual simplicity, store b and y values in a double-index
fashion
"""
function y_double(b1,b2)
    k = length(b1)
    b = zeros(2,k)
    y = zeros(2,k)
    @assert(isodd(k))

    Y = OptimalDual(b1,b2)
    r = round(Int64,(k+1)/2)

    for i = 1:r
        t = k-i+1
        y[2,t] = Y[i]
        y[1,t] = Y[k-i+2]
    end

    return y
end

# Convert from double indexing to linear
function tt(c,t)
    @assert(c == 1 || c == 2)

    if c == 2
        i = k-t+1
    else
        i = t+1
    end
    return i
end

# Now we have some code for the monotonic LP proof
function monotonic_lp_jump(b1,b2,k)
    """
    Maximize epsilon so that relative affinity scores satisfy
    h(t)/b(t) > epsilon + h(t-1)/b(t-1) for top r = (k+1)/2 scores
    """
    m = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(gurobi_env), "OutputFlag" => 0))

    @variable(m, x[1:k+1])  # x[i+1] = proportion of hyperedges of type i
    @variable(m, y)         # variable we maximize to check beta value
    # @variable(m, d1)        # denominator for class 1 affinities
    # @variable(m, d2)        # denominator for class 2 affinities
    for i = 1:k+1
        @constraint(m, x[i] >= 0)
    end
    @constraint(m, 1 == sum(x[i] for i = 1:k+1))

    # Constraint for whether top J are a beta factor above baseline
    r = round(Int64,(k+1)/2)
    for t = r:k
        @constraint(m, (t)*x[t+1]/b1[t] - (t-1)*x[t]/b1[t-1]>= y)
        @constraint(m, (t)*x[k-t+1]/b2[t] - (t-1)*x[k-t+2]/b2[t-1] >= y)
    end

    @objective(m, Max, y)

    JuMP.optimize!(m)
    Y = JuMP.value.(y)
    X = JuMP.value.(x)

    return X,Y
end

function monotonic_LP_B(b1,b2,k)
    B = zeros(k+1,k+1)
    r = round(Int64,(k+1)/2)
    for i = 1:r
        t = k-i+1
        B[i,i] = t/(b2[t])
        B[i,i+1] = -(t-1)/(b2[t-1])
        B[k-i+2,k-i+2] = t/(b1[t])
        B[k-i+2,k-i+1] = -(t-1)/(b1[t-1])
    end
    return B
end

# Now we have some code for the monotonic LP proof
function monotonic_lp_jump_even_k(b1,b2,k)
    """
    Maximize epsilon so that relative affinity scores satisfy
    h(t)/b(t) > epsilon + h(t-1)/b(t-1) for top r = (k+1)/2 scores
    """
    m = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(gurobi_env), "OutputFlag" => 0))

    @variable(m, x[1:k+1])  # x[i+1] = proportion of hyperedges of type i
    @variable(m, y)         # variable we maximize to check beta value
    for i = 1:k+1
        @constraint(m, x[i] >= 0)
    end
    @constraint(m, 1 == sum(x[i] for i = 1:k+1))

    # Constraint for whether top J are a beta factor above baseline
    half = round(Int64,k/2)
    for t = r:k
        @constraint(m, (t)*x[t+1]/b1[t] - (t-1)*x[t]/b1[t-1]>= y)
        @constraint(m, (t)*x[k-t+1]/b2[t] - (t-1)*x[k-t+2]/b2[t-1] >= y)
    end

    @objective(m, Max, y)

    JuMP.optimize!(m)
    Y = JuMP.value.(y)
    X = JuMP.value.(x)

    return X,Y
end

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
# This file was generated, do not modify it. # hide
test = 1:1000
train = last(test)+1:nrows(Xs);


text:
if basename(pwd()) == "aoc"
    cd("2017/14")
end

include("../10/10.jl")

function tobits(str::AbstractString)
    collect(
        b == '1'
        for c in str
        for b in bitstring(parse(UInt8, c, base = 16))[5:8]
    )
end

@assert tobits("a0c2017") == [
    1, 0, 1, 0,
    0, 0, 0, 0,
    1, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 0,
    0, 0, 0, 1,
    0, 1, 1, 1
]

function loadsquares(input::AbstractString)
    squares = zeros(Bool, 128, 128)
    for row in 0:127
        for (col, b) in enumerate(tobits(densehash(input * "-" * string(row))))
            squares[row+1, col] = b
        end
    end
    squares
end

part1(input::AbstractString) = count(loadsquares(input))

@assert part1("flqrgnkx") == 8108
part1("xlqgujun")

function part2(input::AbstractString)
    squares = loadsquares(input)
    stack = CartesianIndex[]
    total = 0
    index = findfirst(squares)
    while !isnothing(index)
        total += 1
        squares[index] = false
        append!(stack, neighbours(index))
        while !isempty(stack)
            index = popfirst!(stack)
            if squares[index]
                squares[index] = false
                append!(stack, neighbours(index))
            end
        end
        index = findfirst(squares)
    end
    total
end

function neighbours(coord)
    [
        coord + c
        for c in [
            CartesianIndex(0, 1),
            CartesianIndex(0, -1),
            CartesianIndex(1, 0),
            CartesianIndex(-1, 0)
        ]
        if c[1] + coord[1] in 1:128 && c[2] + coord[2] in 1:128
    ]
end

@assert part2("flqrgnkx") == 1242
part2("xlqgujun")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
"""
A `PiecewiseDict` has a dictionary for each piece in a partition. Its representation
is a `BlockVector` containing the expansions of all dictionaries combined.
"""
struct PiecewiseDict{P <: Partition,DICTS,S,T} <: CompositeDict{S,T}
    dicts       ::  DICTS
    offsets     ::  Vector{Int}
    partition   ::  P

    function PiecewiseDict{P,DICTS,S,T}(dicts, partition) where {P <: Partition,DICTS,S,T}
        offsets = compute_offsets(dicts)
        new(dicts, offsets, partition)
    end
end

# Make a PiecewiseDict by scaling one set to each of the elements of the partition
function PiecewiseDict(set::Dictionary1d, partition::Partition, n = ones(Int,length(partition))*length(set))
    dicts = [rescale(resize(set, n[i]), support(partition, i)) for i in 1:length(partition)]
    PiecewiseDict(dicts, partition)
end

# Make a PiecewiseDict from a list of dicts and a given partition
function PiecewiseDict(dicts, partition::Partition)
    S = reduce(promote_type, map(domaintype, dicts))
    T = reduce(promote_type, map(codomaintype, dicts))
    # Make sure that the dicts are an indexable list of Dictionary's
    @assert indexable_list(dicts, Dictionary)
    # TODO: We should check that the supports of the dicts match the partition pieces

    PiecewiseDict{typeof(partition),typeof(dicts),S,T}(dicts, partition)
end

# Construct a piecewise set from a list of dicts in 1d
function PiecewiseDict(dicts::Array{S}) where {S <: Dictionary1d}
    # Is this not overly restrictive? A redefinition of the support of a piecewise set could allow disjointed intervals.
    for i in 1:length(dicts)-1
        @assert supremum(support(dicts[i])) ≈ infimum(support(dicts[i+1]))
    end
    points = map(i->infimum(support(i)), dicts)
    push!(points, supremum(support(dicts[end])))
    partition = PiecewiseInterval(points)
    PiecewiseDict(dicts, partition)
end

function PiecewiseDict(set::Dictionary, partition::Partition)
     dicts = [rescale(set, support(partition, i)) for i in 1:length(partition)]
    PiecewiseDict(dicts, partition)
end

name(dict::PiecewiseDict) = "Piecewise dictionary"

dictionaries(set::PiecewiseDict) = set.dicts

partition(set::PiecewiseDict) = set.partition

function stencilarray(dict::PiecewiseDict)
    A = Any[]
    push!(A, "(")
    i = 1
    for s in dictionaries(dict)
        i != 1 && push!(A, ", ")
        push!(A, s)
        i += 1
    end
    push!(A, ")")
    A
end

similardictionary(set::PiecewiseDict, dicts) = PiecewiseDict(dicts, partition(set))

# The set is orthogonal, biorthogonal, etcetera, if all its subsets are.
for op in (:isbasis, :isframe)
    @eval $op(s::PiecewiseDict) = mapreduce($op, &, components(s))
end

# The set is orthogonal, biorthogonal, etcetera, if all its subsets are.
for op in (:isorthogonal, :isbiorthogonal, :isbasis, :isframe)
    @eval $op(s::PiecewiseDict, m::Measure) =
        (@warn "definition unclear"; mapreduce(x->$op(x, m), &, components(s)))# or take intersection of measure and support of dictpiece
end
isorthonormal(s::PiecewiseDict, m::Measure) = false

for op in (:support,)
    @eval $op(set::PiecewiseDict) = $op(partition(set))
    @eval $op(set::PiecewiseDict, idx::Int) = $op(set, multilinear_index(set, idx))
    @eval $op(set::PiecewiseDict, idx) = $op(component(set, outerindex(idx)), innerindex(idx))
end

# The set has a grid and a transform if all its subsets have it
# Disable for now, until grids can be collected into a MultiGrid or something
#for op in (:hasinterpolationgrid, :hastransform)
#    @eval $op(s::PiecewiseDict) = mapreduce($op, &, components(s))
#end

# We have to override getindex for CompositeDict's, because getindex for a
# CompositeDict calls getindex on one of the underlying sets. However, calling
# a function of the underlying set currently does not return zero for an argument
# that is outside its support. In order to get that behaviour, we need to retain
# the PiecewiseDict.
# Perhaps this should change, and any function should be zero outside its support.
getindex(set::PiecewiseDict, i, j) = sub(set, (i,j))

unsafe_eval_element(set::PiecewiseDict, idx::MultilinearIndex, x) =
    x ∈ set.partition[outerindex(idx)] ? unsafe_eval_element( component(set, outerindex(idx)), innerindex(idx), x) : zero(eltype(x))

function eval_expansion(set::PiecewiseDict, x)
    i = partition_index(set, x)
    eval_expansion(component(set, i), x)
end

# TODO: improve, by subdividing the given grid according to the subregions of the piecewise set
evaluation(::Type{T}, dict::PiecewiseDict, gb::GridBasis, grid::AbstractGrid; options...) where {T} =
    ArrayOperator(evaluation_matrix(T, dict, grid), dict, gb) * LinearizationOperator{T}(dict)

for op in (:differentiation, :antidifferentiation)
    @eval function $op(::Type{T}, src::PiecewiseDict, dest::PiecewiseDict, order; options...) where {T}
        @assert ncomponents(src) == ncomponents(dest)
        # TODO: improve the type of the array elements below
        BlockDiagonalOperator{T}(map( (el_s,el_d) -> $op(T, el_s, el_d, order; options...), components(src), components(dest)), src, dest)
    end
end

"""
Make a piecewise function set from a one-dimensional function set, by inserting
the point x in the support of the original set. This yields a PiecewiseDict on
the partition [left(set), x, right(set)].

The original set is duplicated and rescaled to the two subintervals.
"""
function split_interval(set::Dictionary1d, x)
    @assert infimum(support(set)) < x < supremum(support(set))
    points = [infimum(support(set)), x, supremum(support(set))]
    PiecewiseDict(set, PiecewiseInterval(points))
end

split_interval(set::PiecewiseDict, x) = split_interval(set, partition_index(partition(set), x), x)

function split_interval(set::PiecewiseDict, i::Int, x)
    part = partition(set)
    @assert infimum(support(part, i)) < x < supremum(support(part, i))

    part2 = split_interval(part, i, x)
    two_sets = split_interval(component(set, i), x)
    PiecewiseDict(insert_at(components(set), i, components(two_sets)), part2)
end

# Compute the coefficients in an expansion that results from splitting the given
# expansion at a point x.
function split_interval_expansion(set::Dictionary1d, coefficients, x)
    T = eltype(coefficients)
    pset = split_interval(set, x)
    z = zeros(pset)
    pset1 = component(pset, 1)
    pset2 = component(pset, 2)
    # We will manually evaluate the current function at the approximation grids
    # of the newly created sets with smaller support, and use those function values
    # to reconstruct the original function on each subinterval.
    A1 = approximation(pset1)
    A2 = approximation(pset2)
    E1 = evaluation(T, set, grid(src(A1)))
    E2 = evaluation(T, set, grid(src(A2)))
    z1 = A1*(E1*coefficients)
    z2 = A2*(E2*coefficients)
    coefficients!(z, 1, z1)
    coefficients!(z, 2, z2)
    # We return the new set and its expansion coefficients
    pset, z
end

function split_interval_expansion(set::PiecewiseDict, coefficients::BlockVector, x)
    part = partition(set)
    i = partition_index(part, x)
    set_i = component(set, i)
    coef_i = view(coefficients, Block(i))
    split_set, split_coef = split_interval_expansion(set_i, coef_i, x)
    # We compute the types of the individual sets and their coefficients
    # in a hacky way to help inference further on. TODO: fix, because this
    # violates encapsulation and it assumes homogeneous elements
    S = eltype(set.dicts)
    C = eltype(coefficients)

    # Now we want to replace the i-th set by the two new sets, and same for the coefficients
    # Technicalities arise when i is 1 or i equals the ncomponents of the set
    local dicts, coefs
    old_dicts = components(set)
    old_coef = components(coefficients)
    if i > 1
        if i < ncomponents(set)
            # We retain the old elements before and after the new ones
            dicts = S[old_dicts[1:i-1]..., component(split_set, 1), component(split_set, 2), old_dicts[i+1:end]...]
            coefs = C[old_coef[1:i-1]..., component(split_coef, 1), component(split_coef, 2), old_coef[i+1:end]...]
        else
            # We replace the last element, so no elements come after the new ones
            dicts = S[old_dicts[1:i-1]..., component(split_set, 1), component(split_set, 2)]
            coefs = C[old_coef[1:i-1]..., component(split_coef, 1), component(split_coef, 2)]
        end
    else
        # Here i==1, so there are no elements before the new ones
        dicts = S[component(split_set, 1), component(split_set, 2), old_dicts[i+1:end]...]
        coefs = C[component(split_coef, 1), component(split_coef, 2), old_coef[i+1:end]...]
    end
    PiecewiseDict(dicts), BlockVector(coefs)
end

# TODO: add the measure argument here
gram(::Type{T}, dict::PiecewiseDict; options...) where {T} =
    BlockDiagonalOperator(DictionaryOperator{T}[gram(T, component(dict,i); options...) for i in 1:ncomponents(dict)], dict, dict)


text:
using PyCall

pyQt = pyimport("PyQt5.QtWidgets")
pyQtCore = pyimport("PyQt5.QtCore")

# @pydef  mutable struct msgCpResultShowCls <: pyQtCore.QObject
#       function a( self,checked::Bool)
#             msg = pyQt.QMessageBox(btn)
#             msg.setWindowTitle("提示消息：")
#             msg.setText("Copy Completed")
#             msg.setInformativeText("")
#             msg.setIcon(pyQt.QMessageBox.Information )
#             a1= msg.exec()
#             println( "Ok ",a1 ,b.sender().text())
#       end
# end

function showresult()
        msg = pyQt.QMessageBox(nothing)
        msg.setWindowTitle("提示消息：")
        msg.setText("拷贝完毕")
        msg.setInformativeText(" ")
        msg.setIcon(pyQt.QMessageBox.Information )
        a1= msg.exec()
end


text:
using AbstractTrees

mutable struct BinaryNode{T}
    data::T
    parent::BinaryNode{T}
    left::BinaryNode{T}
    right::BinaryNode{T}

    # Root constructor
    BinaryNode{T}(data) where T = new{T}(data)
    # Child node constructor
    BinaryNode{T}(data, parent::BinaryNode{T}) where T = new{T}(data, parent)
end
BinaryNode(data) = BinaryNode{typeof(data)}(data)

function leftchild(data, parent::BinaryNode)
    !isdefined(parent, :left) || error("left child is already assigned")
    node = typeof(parent)(data, parent)
    parent.left = node
end
function rightchild(data, parent::BinaryNode)
    !isdefined(parent, :right) || error("right child is already assigned")
    node = typeof(parent)(data, parent)
    parent.right = node
end


text:
## Spectral clustering on artists
using DelimitedFiles
using SparseArrays
names = readlines("../2-unit-1-demos/artistsim.names")
data = readdlm("../2-unit-1-demos/artistsim.smat")
xy = readdlm("artistsim.xy")
A = sparse(Int.(data[2:end,1]).+1,
           Int.(data[2:end,2]).+1,
           Int.(data[2:end,3]),
           Int(data[1,1]),Int(data[1,1]))
Asym = max.(A,A') # make the graph undirected
n = size(Asym,1)
#Acc = largest_component(Asym)[1] # fix a small technical issue.
# only these nodes are not connected, so we can just remove them.
# I found this with,
# @show findall(largest_component(Asym)[2] .== false)
notconnected = [ 39896, 39897, 42145, 42146, 44404, 44405]
connected = map(x -> !in(x,Set(notconnected)), 1:n)
Acc = Asym[connected,connected]
n = size(Acc,1)
##
using GraphRecipes
using LinearAlgebra
using Plots
ei, ej = findnz(triu(A,1))
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=2, linecolor=1, linealpha=0.01, linewidth=0.5,
  markeralpha=0.2,
  axis_buffer=0.02, background=nothing)
#savefig("artistsim-plot.png")
##
using MatrixNetworks
pr = pagerank(A,0.85)
##
p=sortperm(pr)
scatter(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0)
## Now show the entire graph
p=sortperm(pr)
ei, ej = findnz(triu(A,1))
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.01, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0)
##
savefig("artistsim-pagerank.png")


text:
module mbDTW
export dtwDist, dtwPath

using mbUtil

function dtwMatrix(A,B,d=(x,y)->norm(x-y))
	m,n = length(A),length(B)
	M = Array{Float64}(m,n)

	M[1,1] = d(A[1],B[1])
	for i in 2:m M[i,1] = M[i-1,1] + d(A[i],B[1]) end
	for j in 2:n M[1,j] = M[1,j-1] + d(A[1],B[j]) end

	for i in 2:m, j in 2:n
		M[i,j] = d(A[i],B[j]) + min(M[i-1, j-1], M[i-1, j], M[i, j-1])
	end

	return M
end

dtwDist = last ∘ dtwMatrix

function dtwPath(A,B)
	M = dtwMatrix(A,B)
	m,n = size(M)

	P = [(m,n)]

	function locmin(M,i,j)
		i≤1 && return (i,j-1)
		j≤1 && return (i-1,j)
		d = [(i-1, j-1), (i, j-1), (i-1, j)]
		rtn = d[indmin(map(i->M[i...], d))]
		return rtn
	end

	while P[1] ≠ (1,1)
		unshift!(P, locmin(M, P[1]...))
	end

	return P
end

end #module


text:
using HDF5
using DelimitedFiles

function list_h5_files(foldername::String; prefix::String="data_")
    path     = abspath(foldername)
    allfiles = readdir(path)

    Ns = length(prefix)

    its_names = Tuple[]
    # append only the files whose names start with the given prefix
    for file in allfiles
        try
            if (file[1:Ns] == prefix && file[end-2:end] == ".h5")
                # extract iteration
                it_str = file[Ns+1:end-3]
                fullname = joinpath(path, file)
                # add to list of tuples with iteration and name
                push!(its_names, (parse(Int64, it_str), fullname))
            end
        catch ex
            if isa(ex, BoundsError)
                # probably triggered by string comparison; do nothing
            else
                throw(ex)
            end
        end
    end

    # sort according to iteration
    sort!(its_names)
    # and extract the list of filenames and iterations
    filenames = [name for (it, name) in its_names]
    its       = [it for (it, name) in its_names]

    (its, filenames)
end

function L2_cmh_t(dir_c, dir_m, dir_h)
    # x grid
    xc = h5read(dir_c * "/x.h5", "x")
    xm = h5read(dir_m * "/x.h5", "x")
    xh = h5read(dir_h * "/x.h5", "x")
    # z grid
    zc = h5read(dir_c * "/z.h5", "z")
    zm = h5read(dir_m * "/z.h5", "z")
    zh = h5read(dir_h * "/z.h5", "z")

    # make sure that we can inject points from the medium and high resolution
    # grids in the coarse grid without interpolation
    @assert xc ≈ xm[1:2:end] ≈ xh[1:4:end]
    @assert zc ≈ zm[1:2:end] ≈ zh[1:4:end]

    # dx of coarse resolution
    dxc = xc[2] - xc[1]
    # dz of coarse resolution
    dzc = zc[2] - zc[1]

    # list all available iterations (and corresponding files)
    (its_c, all_filenames_c) = list_h5_files(dir_c, prefix="data_")
    (its_m, all_filenames_m) = list_h5_files(dir_m, prefix="data_")
    (its_h, all_filenames_h) = list_h5_files(dir_h, prefix="data_")

    # we want to compare the common timesteps. we assume here that there is a
    # factor of 2 between the lowest resolution and the corresponding higher
    # resolution ones
    filenames_c = all_filenames_c[:]
    filenames_m = all_filenames_m[1:2:end]
    filenames_h = all_filenames_h[1:4:end]

    Nf = length(filenames_c)
    @assert length(filenames_m) == length(filenames_h) == Nf

    tt     = zeros(Nf)
    L2_cmt = zeros(Nf)
    L2_mht = zeros(Nf)

    # initiate the grid function for the outgoing norm
    l2_out_cmt_x = zeros(length(xc))
    l2_out_mht_x = zeros(length(xc))

    for it in 1:Nf
        file_c = filenames_c[it]
        file_m = filenames_m[it]
        file_h = filenames_h[it]

        ψc  = h5read(file_c, "ψ")
        ψvc = h5read(file_c, "ψv")
        ϕc  = h5read(file_c, "ϕ")
        tc  = h5readattr(file_c, "./")["time"]

        ψm  = h5read(file_m, "ψ")
        ψvm = h5read(file_m, "ψv")
        ϕm  = h5read(file_m, "ϕ")
        tm  = h5readattr(file_m, "./")["time"]

        ψh  = h5read(file_h, "ψ")
        ψvh = h5read(file_h, "ψv")
        ϕh  = h5read(file_h, "ϕ")
        th  = h5readattr(file_h, "./")["time"]

        # make sure we're comparing the same timestep
        @assert tc ≈ tm ≈ th

        # compute the differences between resolutions (projected onto the coarsest grid)
        ψcm  = ψc  -  ψm[1:2:end, 1:2:end]
        ψvcm = ψvc - ψvm[1:2:end, 1:2:end]
        ϕcm  = ϕc  -  ϕm[1:2:end, 1:2:end]
        ψmh  =  ψm[1:2:end, 1:2:end]  - ψh[1:4:end, 1:4:end]
        ψvmh = ψvm[1:2:end, 1:2:end] - ψvh[1:4:end, 1:4:end]
        ϕmh  =  ϕm[1:2:end, 1:2:end]  - ϕh[1:4:end, 1:4:end]

        tt[it]     = tc

        # define the timestep; needed for the sum in u
        dt0 = 0.25*minimum([dxc, dzc])

        # get the outgoing part of the norm
        l2_out_cmt_x += dt0*dzc*sum( ψvcm.*ψvcm + ϕcm.*ϕcm, dims=2)
        l2_out_mht_x += dt0*dzc*sum( ψvmh.*ψvmh + ϕmh.*ϕmh, dims=2)

        # sum the outgoing and ingoing norms to get the complete one
        L2_cmt[it] = sqrt( sum(dzc*dxc*( ψcm.*ψcm ))) + maximum(sqrt.(l2_out_cmt_x))
        L2_mht[it] = sqrt( sum(dzc*dxc*( ψmh.*ψmh ))) + maximum(sqrt.(l2_out_mht_x))
    end

    tt, L2_cmt, L2_mht
end

# maximum times we double resolution
Nmax = 5

# base lowest resolution
Nx = 17
Nz = 16

root_dir  = "/home/thanasis/repos/BondiToy/examples/run00/"
#"/home/mzilhao/dev/julia/BondiToy/examples/run00/"
toy_model = "SH_smooth_B1/"
    #"WH_smooth_B0/"

# we need 3 different resolutions to build the L2 norm that is used in the self
# convergence ratio
for n in 0:1:Nmax-2
    dir_c = joinpath(root_dir, toy_model, "data_$((Nx-1)*2^n + 1)_$(Nz*2^n)")
    dir_m = joinpath(root_dir, toy_model, "data_$((Nx-1)*2^(n+1) + 1)_$(Nz*2^(n+1))")
    dir_h = joinpath(root_dir, toy_model, "data_$((Nx-1)*2^(n+2) + 1)_$(Nz*2^(n+2))")

    tt, L2_cmt, L2_mht = L2_cmh_t(dir_c, dir_m, dir_h)

    data_dir2 = joinpath(root_dir, toy_model, "norms_self_2")
    mkpath(data_dir2)

    outfile  = joinpath(data_dir2, "L2_$(n)_$(n+1)_c$(n).dat")
    open(outfile, "w") do io
        println(io, "# t      |      L2")
        writedlm(io, [tt L2_cmt])
    end

    outfile  = joinpath(data_dir2, "L2_$(n+1)_$(n+2)_c$(n).dat")
    open(outfile, "w") do io
        println(io, "# t      |      L2")
        writedlm(io, [tt L2_mht])
    end
end

max_stars_repo_name: ThanasisGiannakopoulos/BondiToy

max_issues_repo_name: ThanasisGiannakopoulos/BondiToy

max_forks_repo_name: ThanasisGiannakopoulos/BondiToy


text:
"""
    JDF.metadata(indir)

Load the metadata associated with the JDF in `indir`
"""
jdfmetadata(indir) = begin
    open(joinpath(indir,"metadata.jls")) do io
        deserialize(io)
    end
end

"""
    JDF.metadata(indir)

Load the metadata associated with the JDF in `indir`
"""
metadata(jdf::JDFFile) = jdfmetadata(path(jdf))

"""
    JDF.names(indir)

Load the column names associated with the JDF in `indir`

# Examples
```julia
using JDF, DataFrames
savejdf(DataFrame(a = 1:3, b = 1:3), "plsdel.jdf")
JDF.names("plsdel.jdf")
```
"""
names(jdf::JDFFile) = metadata(jdf).names
# names(indir) = jdfmetadata(indir).names

"""
    JDF.nrow(indir)

Returns number of rows of a JDF
"""
nrow(jdf::JDFFile) = metadata(jdf).rows
# nrow(indir) = metadata(indir).rows

"""
    JDF.ncol(indir)

Returns number of rows of a JDF
"""
ncol(jdf::JDFFile) = metadata(jdf).names |> length
# ncol(indir) = metadata(indir).names |> length

"""
    JDF.size(indir)

Returns the JDF's `size`
"""
# JDF.size(indir) = begin
#     m = metadata(indir)
#     (m.rows, length(m.names))
# end
JDF.size(jdf::JDFFile) = begin
     m = metadata(jdf)
     (m.rows, length(m.names))
end

# JDF.size(path::AbstractString) = JDF.size(JDFFile(path))

# JDF.size(indir, v) = JDF.size(indir, Val(v))
JDF.size(jdf::JDFFile, val) = JDF.size(jdf, Val(val))

"""
    JDF.size(indir, 1)

Returns the number of rows
"""
JDF.size(jdf::JDFFile, ::Val{1}) = nrow(jdf)

"""
    JDF.size(indir, 2)

Returns the number of columns
"""
#JDF.size(indir, ::Val{2}) = ncol(indir)
JDF.size(jdf::JDFFile, ::Val{2}) = ncol(jdf)


text:
function trace_info(t::Trace)
    v = collect(values(t))
    address = [n.address for n in v]
    dist = [n.dist for n in v]
    value = [n.value for n in v]
    logprob = [n.logprob for n in v]
    logprob_sum = [n.logprob_sum for n in v]
    observed = [n.observed for n in v]
    interpretation = [n.interpretation for n in v]
    last_interpretation = [n.last_interpretation for n in v]
    (
        address,
        dist,
        value,
        logprob,
        logprob_sum,
        observed,
        interpretation,
        last_interpretation
    )
end

@doc raw"""
    function to_table(t :: Trace)

Turns a trace into a juliadb table. Does not store parent / child relationships.
"""
function to_table(t :: Trace)
    pk = 1:length(t.trace)
    info = trace_info(t)
    table(
        (
            pk=pk,
            address=info[1],
            dist=info[2],
            value=info[3],
            logprob=info[4],
            logprob_sum=info[5],
            observed=info[6],
            interpretation=info[7],
            last_interpretation=info[8]
        );
        pkey=:pk
    )
end

@doc raw"""
    function save(t :: Trace, f)

Saves a trace to disk at the filepath `f`.
"""
function save(t :: Trace, f)
    if endswith(f, ".jdb")
        trace_table = to_table(t)
        JuliaDB.save(trace_table, f)
    else
        error(
            "File must end with .jdb for saving in JuliaDB format."
        )
    end
    f
end

@doc raw"""
    function save(r :: SamplingResults, f)

Saves a `SamplingResults` to disk in the directory `f`.
"""
function save(r :: SamplingResults, f)
    mkpath(f)
    pk = 1:length(r.traces)
    r_table = table(
        (pk = pk,
        log_weight = r.log_weights,
        return_value = r.return_values,
        trace = r.traces);
        pkey=:pk
    )
    JuliaDB.save(r_table, joinpath(f, "results.jdb"))
    open(joinpath(f, "metadata.txt"), "w") do file
        write(file, "Interpretation: $(string(r.interpretation))")
    end
end

@doc raw"""
    function load_jdb(f) :: Trace

Loads a serialized juliadb table from file `f` and converts it into a trace.
"""
function load_jdb(f) :: Trace
    trace_table = JuliaDB.load(f)
    t = trace()
    for ix in 1:length(trace_table)
        row = trace_table[ix]
        address = row.address
        t[address] = node(
            row.value,
            row.address,
            row.dist,
            row.observed,
            row.interpretation
        )
    end
    logprob!(t)
    t
end

@doc raw"""
    function load_csm(f) :: SamplingResults

Loads a `SamplingResults` from a directory. The directory contains `metadata.txt`,
which currently stores the interpretation of the `SamplingResults` (i.e., what kind of
algorithm generated those results), and a file `results.jdb`, which is a JuliaDB table
of the results.
"""
function load_csm(f) :: NonparametricSamplingResults
    interpretation = open(joinpath(f, "metadata.txt"), "r") do f
        line = readline(f)
        interpretation_string = split(line, ": ")[end]
        eval(Meta.parse(interpretation_string))
    end
    results = NonparametricSamplingResults{typeof(interpretation)}(interpretation, Float64[], [], Trace[])
    results_table = JuliaDB.load(joinpath(f, "results.jdb"))
    for ix in 1:length(results_table)
        row = results_table[ix]
        push!(results.log_weights, row.log_weight)
        push!(results.return_values, row.return_value)
        push!(results.traces, row.trace)
    end
    results
end

@doc raw"""
    function load(f)

Load a `Trace` or `SamplingResults` object from file. The file must be saved in JuliaDB
format with ending `.jdb`, which will be interpreted as a single saved trace, or must be
a directory with ending `.csm`, which will be interpreted as a `SamplingResults` object.
"""
function load(f)
    if endswith(f, ".jdb")
        load_jdb(f)
    elseif endswith(f, ".csm")
        load_csm(f)
    else
        error(
            "Can load only .jdb (JuliaDB) files."
        )
    end
end

export save, load, to_table, to_dataframe


text:
using Pkg
using ABCDGraphGenerator
using Random

@info "Usage: julia abcd_sampler.jl config_filename"
@info "For the syntax of config_filename see example_config.toml file"

filename = ARGS[1]
conf = Pkg.TOML.parsefile(filename)
isempty(conf["seed"]) || Random.seed!(parse(Int, conf["seed"]))

μ = haskey(conf, "mu") ? parse(Float64, conf["mu"]) : nothing
ξ = haskey(conf, "xi") ? parse(Float64, conf["xi"]) : nothing
if !(isnothing(μ) || isnothing(ξ))
    throw(ArgumentError("inconsistent data: only μ or ξ may be provided"))
end

n = parse(Int, conf["n"])

τ₁ = parse(Float64, conf["t1"])
d_min = parse(Int, conf["d_min"])
d_max = parse(Int, conf["d_max"])
d_max_iter = parse(Int, conf["d_max_iter"])
@info "Expected value of degree: $(ABCDGraphGenerator.get_ev(τ₁, d_min, d_max))"
degs = ABCDGraphGenerator.sample_degrees(τ₁, d_min, d_max, n, d_max_iter)
open(io -> foreach(d -> println(io, d), degs), conf["degreefile"], "w")

τ₂ = parse(Float64, conf["t2"])
c_min = parse(Int, conf["c_min"])
c_max = parse(Int, conf["c_max"])
c_max_iter = parse(Int, conf["c_max_iter"])
@info "Expected value of community size: $(ABCDGraphGenerator.get_ev(τ₂, c_min, c_max))"
coms = ABCDGraphGenerator.sample_communities(τ₂, c_min, c_max, n, c_max_iter)
open(io -> foreach(d -> println(io, d), coms), conf["communitysizesfile"], "w")

isCL = parse(Bool, conf["isCL"])
islocal = haskey(conf, "islocal") ? parse(Bool, conf["islocal"]) : false

p = ABCDGraphGenerator.ABCDParams(degs, coms, μ, ξ, isCL, islocal)
edges, clusters = ABCDGraphGenerator.gen_graph(p)
open(conf["networkfile"], "w") do io
    for (a, b) in sort!(collect(edges))
        println(io, a, "\t", b)
    end
end
open(conf["communityfile"], "w") do io
    for (i, c) in enumerate(clusters)
        println(io, i, "\t", c)
    end
end