text:
#------------------------------------------------------------------------------

# if the field is finite, enable optimization
function select_divisor(coeffs::Vector{Vector{CoeffFF}}, ch)
    Base.MultiplicativeInverses.UnsignedMultiplicativeInverse(ch)
end

function select_divisor(coeffs::Vector{Vector{CoeffQQ}}, ch)
    ch
end

#------------------------------------------------------------------------------

# Normalize `row` by first coefficient
#
# Finite field magic specialization
function normalize_sparse_row!(row::Vector{CoeffFF}, magic)
    pinv = invmod(row[1], magic.divisor) % magic
    @inbounds for i in 2:length(row)
        # row[i] *= pinv
        row[i] = (row[i] * pinv) % magic
    end
    @inbounds row[1] = one(row[1])
    row
end

# Normalize `row` by first coefficient
#
# Rational field specialization
function normalize_sparse_row!(row::Vector{CoeffQQ}, ch::UInt64)
    pinv = inv(row[1])
    @inbounds for i in 2:length(row)
        row[i] = row[i] * pinv
    end
    row[1] = one(row[1])
    row
end

#------------------------------------------------------------------------------

# reduces row by mul*cfs modulo ch at indices positions
#
# Finite field magic specialization
function reduce_by_pivot!(row::Vector{CoeffFF}, indices::Vector{Int},
                            cfs::Vector{CoeffFF},
                            magic::Base.MultiplicativeInverses.UnsignedMultiplicativeInverse{UInt64})

    # mul = -densecoeffs[i]
    # actually.. not bad!
    @inbounds mul = magic.divisor - row[indices[1]]
    m2  = magic.divisor^2

    # length(row) / length(indices) varies from 10 to 100
    @inbounds for j in 1:length(indices)
        idx = indices[j]
        # x = row[idx] + mul*cfs[j]
        # row[idx] = x - Base.ashr_int(x, 63) & m2
        # row[idx] = x % magic
        # row[idx] = (row[idx] + mul*cfs[j]) % ch
        row[idx] = (row[idx] + mul*cfs[j]) % magic
        # row[idx] = (row[idx] + mul*cfs[j]) & magic.divisor
    end
end

# reduces row by mul*cfs modulo ch at indices positions
#
# Rational field specialization
function reduce_by_pivot!(row::Vector{CoeffQQ}, indices::Vector{Int},
                            cfs::Vector{CoeffQQ}, ch::UInt64)

    # mul = -densecoeffs[i]
    # actually.. not bad!

    mul = -row[indices[1]]

    # length(row) / length(indices) varies from 10 to 100
    @inbounds for j in 1:length(indices)
        idx = indices[j]
        row[idx] = row[idx] + mul*cfs[j]
    end

    nothing
end

#------------------------------------------------------------------------------

# zero entries of densecoeffs and load coefficients cfsref to indices rowexps
#
# Finite field specialization
function load_indexed_coefficients!(densecoeffs::Vector{CoeffFF}, rowexps, cfsref)
    @inbounds for i in 1:length(densecoeffs)
        densecoeffs[i] = CoeffFF(0)
    end
    @inbounds for j in 1:length(rowexps)
        densecoeffs[rowexps[j]] = cfsref[j]
    end
end

# zero entries of densecoeffs and load coefficients cfsref to indices rowexps
#
# Finite field specialization
function load_indexed_coefficients!(densecoeffs::Vector{CoeffQQ}, rowexps, cfsref)
    densecoeffs .= CoeffQQ(0)
    @inbounds for j in 1:length(rowexps)
        densecoeffs[rowexps[j]] = cfsref[j]
    end
end

#------------------------------------------------------------------------------

function reduce_dense_row_by_known_pivots_sparse!(
            densecoeffs::Vector{C}, matrix::MacaulayMatrix{C}, basis::Basis{C},
            pivs::Vector{Vector{Int}}, startcol::Int, tmp_pos::Int, magic;
            exact_colmap::Bool=false) where {C<:Coeff}

    ncols = matrix.ncols
    nleft = matrix.nleft

    # new row nonzero elements count
    k = 0
    uzero = C(0)

    # new pivot index
    np = -1

    for i in startcol:ncols

        # @inbounds if densecoeffs[i] != uzero
        #    densecoeffs[i] = densecoeffs[i] % magic
        # end

        # if row element zero - no reduction
        @inbounds if densecoeffs[i] == uzero
            continue
        end

        # TODO: check this first?
        if !isassigned(pivs, i) || (tmp_pos != -1 && tmp_pos == i)
            if np == -1
                np = i
            end
            k += 1
            continue
        end

        # exponents of reducer row at column i
        reducerexps = pivs[i]

        if exact_colmap # if exact mapping in matrix.low2coef
            @inbounds cfs = matrix.coeffs[tmp_pos] # controversial
        elseif i <= nleft # if reducer is from upper
            @inbounds cfs = basis.coeffs[matrix.up2coef[i]]
        else # of lower part of matrix
            @inbounds cfs = matrix.coeffs[matrix.low2coef[i]]
        end

        # reduce_by_pivot!(densecoeffs, reducerexps, cfs, trick_ch)
        reduce_by_pivot!(densecoeffs, reducerexps, cfs, magic)
    end

    # TODO
    newrow = Vector{Int}(undef, k)
    newcfs = Vector{C}(undef, k)

    # all reduced !
    if k == 0
        return true, newrow, newcfs
    end

    # store new row in sparse format
    # where k - number of structural nonzeros in new reduced row, k > 0
    j = 1
    @inbounds for i in np:ncols # from new pivot
        # densecoeffs[i] = densecoeffs[i] % magic
        @inbounds if densecoeffs[i] != uzero
            newrow[j] = i
            newcfs[j] = densecoeffs[i]
            j += 1
        end
        #=
        if densecoeffs[i] != uzero
            densecoeffs[i] = densecoeffs[i] % magic
        end
        if densecoeffs[i] != uzero
           newrow[j] = i
           newcfs[j] = densecoeffs[i]
           j += 1
       end
       =#
    end

    return false, newrow, newcfs
end

#------------------------------------------------------------------------------

function exact_sparse_rref_par!(matrix::MacaulayMatrix{C}, basis::Basis{C}) where {C<:Coeff}
    ncols  = matrix.ncols
    nlow   = matrix.nlow
    nright = matrix.nright
    nleft  = matrix.nleft

    magic = select_divisor(matrix.coeffs, basis.ch)

    # known pivots
    # no_copy
    pivs = Vector{Vector{Int}}(undef, ncols)
    @inbounds for i in 1:matrix.nup
        pivs[i] = matrix.uprows[i]
    end

    # CHANGED in order to prevent bug
    # when several rows in the matrix are equal
    l2c_tmp = Vector{Int}(undef, max(ncols, matrix.nlow))
    @inbounds for i in 1:nlow
        l2c_tmp[matrix.lowrows[i][1]] = matrix.low2coef[i]
    end
    # CHANGED
    # no_copy
    rowidx2coef = matrix.low2coef
    matrix.low2coef = l2c_tmp

    # unknown pivots
    # (not discovered yet)
    # we will modify them inplace when reducing by pivs
    upivs = matrix.lowrows

    densecoeffs = zeros(C, ncols)

    lowpivs = Vector{Vector{Int64}}(undef, nlow)

    for i in 1:nlow
        # select next row to be reduced
        # npiv ~ exponents
        rowexps = upivs[i]

        # corresponding coefficients from basis
        # (no need to copy here)
        cfsref  = basis.coeffs[rowidx2coef[i]]

        k = 0

        # we load coefficients into dense array
        # into rowexps indices
        # TODO: move this
        # TODO
        load_indexed_coefficients!(densecoeffs, rowexps, cfsref)

        # reduce it with known pivots from matrix.uprows
        # first nonzero in densecoeffs is at startcol position
        startcol = rowexps[1]
        zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1, magic)
        # if fully reduced
        zeroed && continue

        matrix.low2coef[i] = i

        lowpivs[i] = newrow
        # matrix coeffs sparsely stores coefficients of new row
        matrix.coeffs[i] = newcfs

        # normalize if needed
        if !isone( matrix.coeffs[i][1] )
            normalize_sparse_row!(matrix.coeffs[i], magic)
        end
    end

    newpivs = 0

    lowcol2coef = Vector{Int}(undef, ncols)

    for i in 1:nlow
        if !isassigned(lowpivs, i)
            continue
        end

        # select next row to be reduced
        # npiv ~ exponents
        rowexps = lowpivs[i]

        # corresponding coefficients from basis
        # (no need to copy here)
        cfsref  = matrix.coeffs[i]

        # we load coefficients into dense array
        # into rowexps indices
        # TODO: move this
        # TODO
        load_indexed_coefficients!(densecoeffs, rowexps, cfsref)

        startcol = rowexps[1]
        zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, i, magic, exact_colmap=true)

        zeroed && continue

        newpivs += 1

        matrix.lowrows[newpivs] = newrow

        # matrix coeffs sparsely stores coefficients of new row
        matrix.coeffs[i] = newcfs
        # add new pivot at column index newrow[1]
        #  (which is the first nnz column of newrow)
        pivs[newrow[1]]  = newrow
        lowcol2coef[newrow[1]] = i

        # normalize if needed
        if !isone( matrix.coeffs[i][1] )
            normalize_sparse_row!(matrix.coeffs[i], magic)
        end
    end

    matrix.low2coef = lowcol2coef

    # a row to be reduced for each column
    resize!(matrix.lowrows, matrix.nright)

    # shrink matrix
    matrix.npivots = matrix.nrows = matrix.size = newpivs
    resize!(matrix.lowrows, newpivs)
end

#------------------------------------------------------------------------------

function exact_sparse_rref!(matrix::MacaulayMatrix{C}, basis::Basis{C}) where {C<:Coeff}
    ncols  = matrix.ncols
    nlow   = matrix.nlow
    nright = matrix.nright
    nleft  = matrix.nleft

    magic = select_divisor(matrix.coeffs, basis.ch)

    # known pivots
    # no_copy
    pivs = Vector{Vector{Int}}(undef, ncols)
    @inbounds for i in 1:matrix.nup
        pivs[i] = matrix.uprows[i]
    end

    # CHANGED in order to prevent bug
    # when several rows in the matrix are equal
    l2c_tmp = Vector{Int}(undef, max(ncols, matrix.nlow))
    @inbounds for i in 1:nlow
        l2c_tmp[matrix.lowrows[i][1]] = matrix.low2coef[i]
    end
    # CHANGED
    # no_copy
    rowidx2coef = matrix.low2coef
    matrix.low2coef = l2c_tmp

    # unknown pivots
    # (not discovered yet)
    # we will modify them inplace when reducing by pivs
    upivs = matrix.lowrows

    densecoeffs = zeros(C, ncols)

    #=
    @warn "before reducing low"
    dump(matrix, maxdepth=5)
    @warn "lowrow2coef" rowidx2coef
    =#

    for i in 1:nlow
        # select next row to be reduced
        # npiv ~ exponents
        rowexps = upivs[i]

        # corresponding coefficients from basis
        # (no need to copy here)
        cfsref  = basis.coeffs[rowidx2coef[i]]

        # we load coefficients into dense array
        # into rowexps indices
        # TODO: move this
        load_indexed_coefficients!(densecoeffs, rowexps, cfsref)

        # reduce it with known pivots from matrix.uprows
        # first nonzero in densecoeffs is at startcol position
        startcol = rowexps[1]
        # zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1)
        zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1, magic)
        # @warn "reduced " zeroed newrow newcfs
        # if fully reduced
        zeroed && continue

        # matrix coeffs sparsely stores coefficients of new row
        matrix.coeffs[i] = newcfs
        # add new pivot at column index newrow[1]
        #  (which is the first nnz column of newrow)
        pivs[newrow[1]]  = newrow
        # set ref to coefficient to matrix
        # guaranteed to be from lower part
        matrix.low2coef[newrow[1]] = i

        # normalize if needed
        if !isone( matrix.coeffs[i][1] )
            normalize_sparse_row!(matrix.coeffs[i], magic)
        end
    end

    #=
    @warn "after reducing low"
    dump(matrix, maxdepth=5)
    println("PIVS: ", pivs)
    =#

    # number of new pivots
    newpivs = 0

    # a row to be reduced for each column
    resize!(matrix.lowrows, matrix.nright)

    # interreduce new pivots..
    # .. for each right (non-pivotal) column
    for i in 1:nright
        k = ncols - i + 1
        if isassigned(pivs, k)

            if k <= nleft
                cfsref = basis.coeffs[matrix.up2coef[k]]
            else # of lower part of matrix
                cfsref = matrix.coeffs[matrix.low2coef[k]]
            end

            startcol = pivs[k][1]

            @assert length(cfsref) == length(pivs[k])

            load_indexed_coefficients!(densecoeffs, pivs[k], cfsref)

            newpivs += 1

            zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, startcol, magic)

            # update row and coeffs
            matrix.lowrows[newpivs] = newrow
            matrix.coeffs[matrix.low2coef[k]] = newcfs
            matrix.low2coef[k] = matrix.low2coef[k]
            pivs[k] = matrix.lowrows[newpivs]
        end
    end

    # shrink matrix
    matrix.npivots = matrix.nrows = matrix.size = newpivs
    resize!(matrix.lowrows, newpivs)
end

function linear_algebra!(matrix::MacaulayMatrix, basis::Basis, linalg::Val{:exact}, rng)
    resize!(matrix.coeffs, matrix.nlow)
    exact_sparse_rref!(matrix, basis)
end

function linear_algebra!(matrix::MacaulayMatrix, basis::Basis, linalg::Val{:prob}, rng)
    resize!(matrix.coeffs, matrix.nlow)
    randomized_sparse_rref!(matrix, basis, rng)
end

#------------------------------------------------------------------------------

function nblocks_in_randomized(nlow::Int)
    floor(Int, sqrt(nlow / 3)) + 1
end

function randomized_sparse_rref!(matrix::MacaulayMatrix{C}, basis::Basis{C}, rng) where {C<:Coeff}
    ncols  = matrix.ncols
    nlow   = matrix.nlow
    nright = matrix.nright
    nleft  = matrix.nleft

    magic = select_divisor(matrix.coeffs, basis.ch)

    # known pivots
    # no_copy
    pivs = Vector{Vector{Int}}(undef, ncols)
    @inbounds for i in 1:matrix.nup
        pivs[i] = matrix.uprows[i]
    end

    # CHANGED in order to prevent bug
    # when several rows in the matrix are equal
    l2c_tmp = Vector{Int}(undef, max(ncols, matrix.nlow))
    @inbounds for i in 1:nlow
        l2c_tmp[matrix.lowrows[i][1]] = matrix.low2coef[i]
    end
    # CHANGED
    # no_copy
    rowidx2coef = matrix.low2coef
    matrix.low2coef = l2c_tmp

    # unknown pivots
    # (not discovered yet)
    # we will modify them inplace when reducing by pivs
    upivs = matrix.lowrows

    nblocks = nblocks_in_randomized(nlow)
    rem = nlow % nblocks == 0 ? 0 : 1
    rowsperblock = div(nlow, nblocks) + rem

    # @warn "" nlow nblocks rem rowsperblock

    densecoeffs = zeros(UInt64, ncols)
    mulcoeffs = zeros(UInt64, rowsperblock)

    #=
    @warn "before reducing low"
    dump(matrix, maxdepth=5)
    @warn "lowrow2coef" rowidx2coef
    =#

    for i in 1:nblocks
        nrowsupper = nlow > i * rowsperblock ? i * rowsperblock : nlow
        nrowstotal = nrowsupper - (i-1)*rowsperblock

        # @warn "cycle" nrowsupper nrowstotal

        nrowstotal == 0 && continue # end right here

        ctr = 0
        while ctr < nrowstotal

            @inbounds for j in 1:nrowstotal
                # TODO: fixed generator
                mulcoeffs[j] = rand(rng, UInt64) % magic
            end

            densecoeffs .= UInt64(0)
            startcol = ncols

            @inbounds for k in 1:nrowstotal
                rowidx = (i-1)*rowsperblock + k

                rowexps = upivs[rowidx]
                cfsref  = basis.coeffs[rowidx2coef[rowidx]]

                startcol = min(startcol, rowexps[1])

                @inbounds for l in 1:length(rowexps)
                    ridx = rowexps[l]
                    densecoeffs[ridx] = (densecoeffs[ridx] + mulcoeffs[k]*cfsref[l]) % magic
                end
            end

            # reduce it with known pivots from matrix.uprows
            # first nonzero in densecoeffs is at startcol position
            zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1, magic)
            # @warn "reduced " zeroed newrow newcfs
            # if fully reduced

            # this is pure insight

            if zeroed
                ctr = nrowstotal
                break
            end

            absolute_i = (i-1)*rowsperblock + ctr + 1

            # matrix coeffs sparsely stores coefficients of new row
            matrix.coeffs[absolute_i] = newcfs
            # add new pivot at column index newrow[1]
            #  (which is the first nnz column of newrow)
            pivs[newrow[1]]  = newrow
            # set ref to coefficient to matrix
            # guaranteed to be from lower part
            matrix.low2coef[newrow[1]] = absolute_i

            # normalize if needed
            if matrix.coeffs[absolute_i][1] != 1
                normalize_sparse_row!(matrix.coeffs[absolute_i], magic)
            end

            ctr += 1

        end
    end

    #=
    @warn "after reducing low"
    dump(matrix, maxdepth=5)
    println("PIVS: ", pivs)
    =#

    # number of new pivots
    newpivs = 0

    # a row to be reduced for each column
    resize!(matrix.lowrows, matrix.nright)

    # interreduce new pivots..
    # .. for each right (non-pivotal) column
    densecfs = zeros(UInt64, ncols)

    for i in 1:nright
        k = ncols - i + 1
        if isassigned(pivs, k)
            densecfs .= UInt64(0)

            if k <= nleft
                cfsref = basis.coeffs[matrix.up2coef[k]]
            else # of lower part of matrix
                cfsref = matrix.coeffs[matrix.low2coef[k]]
            end

            @inbounds startcol = pivs[k][1]

            @assert length(cfsref) == length(pivs[k])

            @inbounds for j in 1:length(pivs[k])
                densecfs[pivs[k][j]] = cfsref[j]
            end
            newpivs += 1

            zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecfs, matrix, basis, pivs, startcol, startcol, magic)

            # update row and coeffs
            matrix.lowrows[newpivs] = newrow
            matrix.coeffs[matrix.low2coef[k]] = newcfs
            matrix.low2coef[k] = matrix.low2coef[k]
            pivs[k] = matrix.lowrows[newpivs]
        end
    end

    # shrink matrix
    matrix.npivots = matrix.nrows = matrix.size = newpivs
    resize!(matrix.lowrows, newpivs)
end

#------------------------------------------------------------------------------

function interreduce_matrix_rows!(matrix::MacaulayMatrix{C}, basis::Basis{C}) where {C<:Coeff}
    resize!(matrix.lowrows, matrix.ncols)
    resize!(matrix.up2coef, matrix.ncols)
    resize!(matrix.low2coef, matrix.ncols)
    resize!(matrix.coeffs, matrix.ncols)

    magic = select_divisor(matrix.coeffs, basis.ch)

    # same pivs as for rref
    # pivs: column idx --> vector of present columns
    pivs = Vector{Vector{Int}}(undef, matrix.ncols)
    @inbounds for i in 1:matrix.nrows
        pivs[matrix.uprows[i][1]] = matrix.uprows[i]
        matrix.low2coef[matrix.uprows[i][1]] = i
        matrix.coeffs[i] = copy(basis.coeffs[matrix.up2coef[i]])
    end

    densecfs = zeros(C, matrix.ncols)

    k = matrix.nrows

    @inbounds for i in 1:matrix.ncols
        l = matrix.ncols - i + 1
        if isassigned(pivs, l)
            # densecfs .= UInt64(0)
            cfs = matrix.coeffs[matrix.low2coef[l]]
            reducexps = pivs[l]

            startcol = reducexps[1]

            @inbounds for j in 1:length(reducexps)
                densecfs[reducexps[j]] = cfs[j]
            end

            zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecfs, matrix, basis, pivs, l, l, magic)

            # TODO: in `densecfs` zeroe only indices from `newrow`
            matrix.lowrows[k] = newrow
            matrix.coeffs[matrix.low2coef[l]] = newcfs
            pivs[l] = matrix.lowrows[k]
            k -= 1

            for idx in newrow
                densecfs[idx] = zero(densecfs[idx])
            end
        end
    end

    # hmm
    # TODO
    matrix.npivots = matrix.nrows
end

#------------------------------------------------------------------------------

function convert_hashes_to_columns!(
            matrix::MacaulayMatrix, symbol_ht::MonomialHashtable)

    # col2hash = matrix.col2hash
    hdata    = symbol_ht.hashdata
    load     = symbol_ht.load

    # monoms from symbolic table represent one column in the matrix

    col2hash = Vector{Int}(undef, load - 1)
    j = 1
    # number of pivotal cols
    k = 0
    @inbounds for i in symbol_ht.offset:load
        # column to hash index
        col2hash[j] = i
        j += 1

        # meaning the column is pivoted
        if hdata[i].idx == 2
            k += 1
        end
    end

    # sort columns
    # TODO
    sort_columns_by_hash!(col2hash, symbol_ht)

    matrix.nleft = k
    # -1 as long as hashtable load is always 1 more than actual
    matrix.nright = load - matrix.nleft - 1

    # store the other direction of mapping,
    # hash -> column
    for k in 1:length(col2hash)
        hdata[col2hash[k]].idx = k
    end

    nterms = 0
    for k in 1:matrix.nup
        row = matrix.uprows[k]

        @inbounds for j in 1:length(row)
            row[j] = hdata[row[j]].idx
        end

        # TODO: not needed for now
        nterms += length(row)
    end

    for k in 1:matrix.nlow
        row = matrix.lowrows[k]

        @inbounds for j in 1:length(row)
            row[j] = hdata[row[j]].idx
        end

        # TODO
        nterms += length(row)
    end

    matrix.ncols = matrix.nleft + matrix.nright

    @assert matrix.nleft + matrix.nright == symbol_ht.load-1 == matrix.ncols
    @assert matrix.nlow + matrix.nup == matrix.nrows

    matrix.col2hash = col2hash
end

function exact_sparse_linear_algebra_isgroebner!(matrix::MacaulayMatrix{C}, basis::Basis{C}) where {C<:Coeff}
    ncols  = matrix.ncols
    nlow   = matrix.nlow
    nright = matrix.nright
    nleft  = matrix.nleft

    magic = select_divisor(matrix.coeffs, basis.ch)

    # known pivots
    # no_copy
    pivs = Vector{Vector{Int}}(undef, ncols)
    @inbounds for i in 1:matrix.nup
        pivs[i] = matrix.uprows[i]
    end

    # CHANGED in order to prevent bug
    # when several rows in the matrix are equal
    l2c_tmp = Vector{Int}(undef, max(ncols, matrix.nlow))
    @inbounds for i in 1:nlow
        l2c_tmp[matrix.lowrows[i][1]] = matrix.low2coef[i]
    end
    # CHANGED
    # no_copy
    rowidx2coef = matrix.low2coef
    matrix.low2coef = l2c_tmp

    # unknown pivots
    # (not discovered yet)
    # we will modify them inplace when reducing by pivs
    upivs = matrix.lowrows

    densecoeffs = zeros(C, ncols)

    #=
    @warn "before reducing low"
    dump(matrix, maxdepth=5)
    @warn "lowrow2coef" rowidx2coef
    =#

    for i in 1:nlow
        # select next row to be reduced
        # npiv ~ exponents
        rowexps = upivs[i]

        # corresponding coefficients from basis
        # (no need to copy here)
        cfsref  = basis.coeffs[rowidx2coef[i]]

        k = 0

        # we load coefficients into dense array
        # into rowexps indices
        # TODO: move this

        load_indexed_coefficients!(densecoeffs, rowexps, cfsref)

        # reduce it with known pivots from matrix.uprows
        # first nonzero in densecoeffs is at startcol position
        startcol = rowexps[1]
        # zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1)
        zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1, magic)
        # @warn "reduced " zeroed newrow newcfs
        # if fully reduced
        zeroed && continue

        return false
    end
    return true
end

#------------------------------------------------------------------------------

function exact_sparse_rref_nf!(
        matrix::MacaulayMatrix{C},
        tobereduced::Basis{C},
        basis::Basis{C}) where {C<:Coeff}

    ncols  = matrix.ncols
    nlow   = matrix.nlow
    nright = matrix.nright
    nleft  = matrix.nleft

    magic = select_divisor(matrix.coeffs, basis.ch)

    # known pivots
    # no_copy
    pivs = Vector{Vector{Int}}(undef, ncols)
    @inbounds for i in 1:matrix.nup
        pivs[i] = matrix.uprows[i]
    end

    # CHANGED in order to prevent bug
    # when several rows in the matrix are equal
    l2c_tmp = Vector{Int}(undef, max(ncols, matrix.nlow))
    @inbounds for i in 1:nlow
        l2c_tmp[matrix.lowrows[i][1]] = matrix.low2coef[i]
    end
    # CHANGED
    # no_copy
    rowidx2coef = matrix.low2coef
    matrix.low2coef = l2c_tmp

    # unknown pivots
    # (not discovered yet)
    # we will modify them inplace when reducing by pivs
    upivs = matrix.lowrows

    densecoeffs = zeros(C, ncols)

    #=
    @warn "before reducing low"
    dump(matrix, maxdepth=5)
    @warn "lowrow2coef" rowidx2coef
    =#

    # @error "pivs" pivs

    for i in 1:nlow
        # select next row to be reduced
        # npiv ~ exponents
        rowexps = upivs[i]

        # corresponding coefficients from basis
        # (no need to copy here)
        cfsref  = tobereduced.coeffs[rowidx2coef[i]]

        k = 0

        # we load coefficients into dense array
        # into rowexps indices
        # TODO: move this

        load_indexed_coefficients!(densecoeffs, rowexps, cfsref)

        # @error "row " i densecoeffs cfsref

        # reduce it with known pivots from matrix.uprows
        # first nonzero in densecoeffs is at startcol position
        startcol = rowexps[1]
        # zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1)
        zeroed, newrow, newcfs = reduce_dense_row_by_known_pivots_sparse!(densecoeffs, matrix, basis, pivs, startcol, -1, magic)
        # @warn "reduced " zeroed
        # if fully reduced
        zeroed && continue

        # matrix coeffs sparsely stores coefficients of new row
        matrix.coeffs[i] = newcfs
        # add new pivot at column index newrow[1]
        #  (which is the first nnz column of newrow)
        matrix.lowrows[i]  = newrow
        # set ref to coefficient to matrix
        # guaranteed to be from lower part
        matrix.low2coef[i] = i
    end
    matrix.npivots = matrix.nrows = matrix.size = matrix.nlow
end

function exact_sparse_linear_algebra_nf!(
        matrix::MacaulayMatrix,
        tobereduced::Basis,
        basis::Basis)

    resize!(matrix.coeffs, matrix.nlow)
    exact_sparse_rref_nf!(matrix, tobereduced, basis)
end

#------------------------------------------------------------------------------

function convert_matrix_rows_to_basis_elements!(
            matrix::MacaulayMatrix, basis::Basis,
            ht::MonomialHashtable, symbol_ht::MonomialHashtable)

    # we mutate basis array directly by adding new elements

    check_enlarge_basis!(basis, matrix.npivots)
    rows = matrix.lowrows
    crs = basis.ndone

    # @warn "New pivots" matrix.npivots

    @inbounds for i in 1:matrix.npivots
        @inbounds colidx = rows[i][1]

        insert_in_basis_hash_table_pivots(rows[i], ht, symbol_ht, matrix.col2hash)
        # TODO : a constant

        # an interesing way to find coefficients
        basis.coeffs[crs + i] = matrix.coeffs[matrix.low2coef[colidx]]
        basis.gens[crs + i] = matrix.lowrows[i]
    end

    basis.ntotal += matrix.npivots
end

function convert_matrix_rows_to_basis_elements_use_symbol!(
            matrix::MacaulayMatrix, basis::Basis)

    check_enlarge_basis!(basis, matrix.npivots)

    crs = basis.ndone
    rows = matrix.lowrows

    @inbounds for i in 1:matrix.npivots
        row = rows[i]
        colidx = row[1]

        @inbounds for j in 1:length(row)
            row[j] = matrix.col2hash[row[j]]
        end

        basis.coeffs[crs + i] = matrix.coeffs[matrix.low2coef[colidx]]
        basis.gens[crs + i] = row
    end
end

#------------------------------------------------------------------------------

function convert_nf_rows_to_basis_elements!(
        matrix::MacaulayMatrix, basis::Basis,
        ht::MonomialHashtable, symbol_ht::MonomialHashtable)

    check_enlarge_basis!(basis, matrix.npivots)

    @inbounds for i in 1:matrix.npivots
        basis.ndone += 1
        basis.nlead += 1
        basis.nonred[basis.nlead] = basis.ndone
        if isassigned(matrix.coeffs, i)
            row = matrix.lowrows[i]
            insert_in_basis_hash_table_pivots(row, ht, symbol_ht, matrix.col2hash)

            colidx = row[1]
            basis.coeffs[basis.ndone] = matrix.coeffs[i]
            basis.gens[basis.ndone] = row
        else
            # TODO
            empty!(basis.coeffs[basis.ndone])
            empty!(basis.gens[basis.ndone])
        end
    end
end


text:
# push!(LOAD_PATH, "/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/")
include("EbitCloudSpatialExtend.jl")
include("EbitODEMessages.jl")
include("EbitSolver.jl")
include("EbitODEServer.jl")

module EbitServer
using ArgParse
using DifferentialEquations

using Main.EbitSolver
using Main.EbitODEServer
using Main.EbitODEMessages

function parse_commandline(ARGS::Vector{String})
    s = ArgParseSettings(description="Kairos ODE Server", version = "0.0.1")
    @add_arg_table s begin
        "--port", "-p"
            help = "Port used to connect to server."
            arg_type = Int
            default = 2000
    end

    return parse_args(ARGS,s)
end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    arguments = parse_commandline(ARGS)
    Main.EbitODEServer.start_ode_server(arguments["port"])
    return 0
end

end

ext: jl


text:
Random.seed!(123)

n = 2000
h = readdlm("h.txt")
A = Conv(Float64,(n,),h[:])
m = size(A,1)
#A = hcat([[zeros(i);h;zeros(n-1-i)] for i = 0:n-1]...) # Equivalent Full Matrix

x_star = Vector(sprandn(n,0.5))
b = A*x_star+20*randn(m)

f = Translate(SqrNormL2(), -b)
lam = 1e-1*norm(A'*b, Inf)
g = NormL1(lam)

## fast FBS/Adaptive
x0 = zeros(n)
@time it, xFBS, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, fast = true, tol = 1e-8)
@test it < 1000

# ZeroFPR/Adaptive
x0 = zeros(n)
@time it, xZ, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, tol = 1e-8)

@test norm(xFBS-xZ) < 1e-3
@test it < 70

# PANOC/Adaptive
x0 = zeros(n)
@time it, xP, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, tol = 1e-8)

@test norm(xP-xZ) < 1e-3
@test it < 140

#####################
# Complex Variables #
#####################

Random.seed!(123)
m, n = 10,100

A = randn(m,n)+im*randn(m,n)
U,S,V = svd(A)
#create ill conditioned matrix
S[floor(Int,0.5*min(m,n)):end] .= 0.
A = U*diagm(0 => S)*V'

x_star = sprandn(n,0.5).+im.*sprandn(n,0.5)
b = A*x_star+20*(randn(m).+im*randn(m))

f = Translate(SqrNormL2(), -b)
f2 = LeastSquares(A, b)
lam = 1e-2*norm(A'*b, Inf)
g = NormL1(lam)

x0 = zeros(n)+im*zeros(n)
it, x_star, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, verbose = 0, tol = 1e-9)

# fast FBS/Adaptive
x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, fast = true, tol = 1e-6, maxit = 30000)

@test norm(x-x_star) < 1e-4
@test it < 25000

# ZeroFPR/Adaptive
x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, tol = 1e-6)

@test norm(x-x_star) < 1e-4
@test it < 700

# PANOC/Adaptive
x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, tol = 1e-6)

@test norm(x-x_star) < 1e-4
@test it < 1500


text:
using Random, Statistics, Distributions, GLM, ProgressMeter

const SIMULATIONS = 10_000

const N = 10:5:100
const T = 120:120:1200
const K = 1

pvalues = zeros((SIMULATIONS, length(N), length(T))) .* NaN
tstats = zeros((SIMULATIONS, length(N), length(T))) .* NaN

progress = Progress(SIMULATIONS*length(N)*length(T), desc="Running Simulations")
@inbounds for i in 1:SIMULATIONS
    @inbounds for (n,N) in enumerate(N), (t,T) in enumerate(T)
        next!(progress)
        Rt = rand(Normal(0, 1/252), (T,N))
        ft = rand(Normal(), (T,K))

        beta = (Rt'Rt)^-1 * Rt'ft
        mu = ft - Rt*beta

        gamma = zeros(T,2)*NaN
        X = [ones(N) beta]
        for t in 1:T
            gamma[t,:] = (X'X)^-1 * X'Rt[t,:]
        end

        tstat = mean(gamma[:,2]) / (std(gamma[:,2])/sqrt(T))
        pvalue = cdf(TDist(T-1), tstat)
        pvalues[i,n,t] = pvalue
        tstats[i,n,t] = tstat
    end
end

max_stars_repo_name: dysonance/Preferences

max_issues_repo_name: dysonance/Preferences

max_forks_repo_name: dysonance/Preferences


text:
#############################################################################
#############################################################################
#
# This file contains units tests for polynomial factorization
#
#############################################################################
#############################################################################

"""
Test factorization of polynomials.
"""
function mod_factor_test_poly(;N::Int = 5, seed::Int = 0, primes::Vector{Int} = [3,11,17,19,67,101,151])
    Random.seed!(seed)
    for prime in primes
        print("\ndoing prime = $prime \t")
        for _ in 1:N
            print(".")
            p = rand(PolynomialModP,prime)
            factorization = factor(p)
            pr = expand_factorization(factorization)
            @assert p-pr == 0
        end
    end
    println("\nfactor_test_poly - PASSED")
end

mod_factor_test_poly()


text:
function _default_exchange_logic!(iteration_index::Int64, game_world::VLGameWorld,
    order_array::Array{VLOrderModel,1}, current_asset_price_array::Array{Float64,2})::Array{Float64,2}

    try

        # initialize -
        agent_array = game_world.agent_array
        number_of_assets = size(current_asset_price_array, 2)
        λ = game_world.liquidity_parameter_array
        number_of_orders = length(order_array)
        order_book = zeros(Int64, number_of_assets, number_of_orders)

        # build agent table -
        agent_table = Dict{UUID,VLAgentModel}()
        for agent in agent_array
            agent_id = agent.agent_id
            agent_table[agent_id] = agent
        end

        # build order lookup table -
        order_agent_lookup_table = Dict{UUID,VLAgentModel}()
        for order in order_array

            # get info about the order -
            order_id = order.order_id           # get the order id -
            agent_id = order.agent_id           # get the id of the agent who entered this order -
            agent_model = agent_table[agent_id] # get the model for the agent who entered this order -

            # link the order_id and the agent model -
            order_agent_lookup_table[order_id] = agent_model
        end

        # let process the orders - first, build the order book
        for order_index = 1:number_of_orders

            # grab the order -
            order_model = order_array[order_index]
            asset_index = order_model.asset_index   # what asset index is this order?
            order_quantity = order_model.quantity   # what is the order quantity?
            order_action = order_model.action       # what is the order action?

            # add order to the order book for this asset -
            if (order_action == 1)
                order_book[asset_index, order_index] = order_quantity
            elseif (order_action == -1)
                order_book[asset_index, order_index] = -1 * order_quantity
            else
                order_book[asset_index, order_index] = 0
            end
        end

        # generate total order for each asset -
        total_order_quantity = sum(order_book, dims=2)

        # update the price for each asset -
        new_asset_price_array = copy(current_asset_price_array)
        for asset_index = 1:number_of_assets
            current_price = current_asset_price_array[iteration_index - 1, asset_index]
            new_asset_price_array[iteration_index, asset_index] = current_price + (1 / λ[asset_index]) * total_order_quantity[asset_index]
        end

        # default: assume all orders have been filled, update agents
        for order_model in order_array

            order_id = order_model.order_id                         # get the order id -
            agent_model = order_agent_lookup_table[order_id]        # get the agent who made this order -

            # update -
            agent_model.agent_update_logic(agent_model, order_model, new_asset_price_array, iteration_index)
        end

        return new_asset_price_array
    catch error
        throw(error)
    end
end


text:
# using CSV, DataFrames, RCall, StatsBase
#
# csvpath = "C:/Users/dzj/Downloads/Lending Club/csv/"
# dirs = readdir(csvpath)
#
# files = map(dirs) do dir1
#     res = readdir(csvpath * dir1)
#     joinpath(csvpath, dir1, res[1])
# end
#
#
# function fread(path1)
#     R"library(data.table);a = fread($path1 , na.strings = NULL)"
#     @rget a
#     a
# end
#
# data = mapreduce(vcat, files) do file
#     println(file)
#     fread(file)
# end

# CSV.read("C:/Users/dzj/Downloads/Lending Club/csv/LoanStats3b.csv/LoanStats3b.csv"; datarow = 2)
# path1 = "C:/Users/dzj/Downloads/Lending Club/csv/LoanStats3a.csv/LoanStats3a.csv"
# df = fread("C:/Users/dzj/Downloads/Lending Club/csv/LoanStats3a.csv/LoanStats3a.csv")

# function c()
#     local i = 1
#     colwise(df) do col
#         if typeof(col) == DataArrays.DataArray{String, 1}
#             print(i)
#         end
#         i = i + 1
#     end
# end
#
# c()

# function c1()
#     local i = 1
#     map(df[22]) do df1
#         i
#         if typeof(df1) != String
#             print(df1)
#         end
#     end
# end


text:
#=
BFSLevelsMultiSourceTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BFSLevelsMultiSource.jl")
matrix =  GBMatrix([[false, false, false, true, false, false, false] [true, false, false, false, false, false, false] [false, false, false, true, false, true, true] [true, false, false, false, false, false, true] [false, true, false, false, false, false, true] [false, false, true, false, true, false, false] [false, true, false, false, false, false, false]])
S = GBMatrix([[true, false, false] [false, false, false] [false, true, false] [false, false, true] [false, false, false] [false, false, false] [false, false, false]])
@test ms_bfsl(matrix, S, 7, 3) == GBMatrix([[1, 0, 2] [2, 0, 3] [3, 1, 2] [2, 0, 1] [3, 0, 4] [4, 2, 3] [3, 0, 4]])


text:
function fit_mle_all_slices(hmm::PeriodicHMM, observations, n2t::AbstractArray{Int};
    robust = false,
    smooth = true, window = -15:15, kernel = :step,
    history = false,
    kwargs...)

    hmm = copy(hmm)
    N = size(observations, 1)
    @argcheck size(observations, 1) == size(n2t, 1)

    K, T = size(hmm, 1), size(hmm, 3)
    α = hcat([vec(sum(hmm.A[:, :, t], dims = 1) / K) for t = 1:T]...)
    n_in_t = [findall(n2t .== t) for t = 1:T] #
    hist = Vector{HMMBase.EMHistory}(undef, T)
    cycle = CyclicArray(1:T, "1D")
    for t = 1:T
        n_in_t_extanded = sort(vcat([n_in_t[tt] for tt in cycle[[t - 12, t - 6, t, t + 6, t + 12]]]...)) # extend dataset
        hist[t] = fit_mle_B_slice!(@view(α[:, t]), @view(hmm.B[:, t]), observations[n_in_t_extanded, :]; kwargs...)
        # hist[t] = fit_mle_B_slice!(@view(α[:,t]), @view(hmm.B[:,t]), observations[n_in_t[t],:]; kwargs...)
    end

    LL = zeros(N, K)

    if smooth == true
        smooth_B = lagou(hmm.B, dims = 2, window = window, kernel = kernel)
        smooth_α = lagou(α, dims = 2, window = window, kernel = kernel)

        # evaluate likelihood for each type k
        loglikelihoods!(LL, smooth_B, observations, n2t)
        [LL[n, k] += log(smooth_α[k, n2t[n]]) for k = 1:K, n = 1:N]
    else
        # evaluate likelihood for each type k
        loglikelihoods!(LL, hmm.B, observations, n2t)
        [LL[n, k] += log(α[k, n2t[n]]) for k = 1:K, n = 1:N]
    end
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    fit_mle_A_from_slice!(hmm.a, hmm.A, LL, n2t, robust = robust)

    return history ? (hmm, hist) : hmm
end

function fit_mle_all_slices(hmm::HierarchicalPeriodicHMM, observations, n2t::AbstractArray{Int};
    robust = false,
    smooth = true, window = -15:15, kernel = :step,
    history = false,
    kwargs...)

    hmm = copy(hmm)

    N, K, T, size_memory, D = size(observations, 1), size(hmm, 1), size(hmm, 3), size(hmm, 4), size(hmm, 2)
    @argcheck size(observations, 1) == size(n2t, 1)

    # assign category for observation depending in the past observations
    memory = Int(log(size_memory) / log(2))
    lag_cat = conditional_to(observations, memory)

    α = hcat([vec(sum(hmm.A[:, :, t], dims = 1) / K) for t = 1:T]...)
    n_in_t = [findall(n2t .== t) for t = 1:T] #
    hist = Vector{SHHMM.EMHistory}(undef, T)
    cycle = CyclicArray(1:T, "1D")
    for t = 1:T
        n_in_t_extanded = sort(vcat([n_in_t[tt] for tt in cycle[[t - 12, t - 7, t, t + 6, t + 13]]]...)) # extend dataset
        hist[t] = fit_mle_B_slice!(@view(α[:, t]), @view(hmm.B[:, t, :, :]), observations[n_in_t_extanded, :], lag_cat[n_in_t_extanded, :]; kwargs...)
        # hist[t] = fit_mle_B_slice!(α[:,t], @view(hmm.B[:,t,:,:]), observations[n_in_t[t],:], lag_cat[n_in_t[t],:]; kwargs...)
    end

    LL = zeros(N, K)

    if smooth == true
        smooth_B = lagou(hmm.B, dims = 2, window = window, kernel = kernel)
        smooth_α = lagou(α, dims = 2, window = window, kernel = kernel)
        loglikelihoods!(LL, smooth_B, observations, n2t, lag_cat)
        [LL[n, k] += log(smooth_α[k, n2t[n]]) for k = 1:K, n = 1:N]
    else
        loglikelihoods!(LL, hmm.B, observations, n2t, lag_cat)
        [LL[n, k] += log(α[k, n2t[n]]) for k = 1:K, n = 1:N]
    end
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    fit_mle_A_from_slice!(hmm.a, hmm.A, LL, n2t, robust = robust)

    return history ? (hmm, hist) : hmm
end

function loglikelihoods!(LL::AbstractMatrix, B::AbstractArray{F,2} where {F<:MultivariateDistribution}, observations, n2t::AbstractArray{Int})
    N, K = size(observations, 1), size(B, 1)
    @argcheck size(LL) == (N, K)

    for n in OneTo(N)
        t = n2t[n]
        for k in OneTo(K)
            LL[n, k] = logpdf(B[k, t], observations[n, :])
        end
    end
end

function loglikelihoods!(LL::AbstractMatrix, B::AbstractArray{F,2} where {F<:UnivariateDistribution}, observations, n2t::AbstractArray{Int})
    N, K = size(observations, 1), size(B, 1)
    @argcheck size(LL) == (N, K)

    for n in OneTo(N)
        t = n2t[n]
        for k in OneTo(K)
            LL[n, k] = logpdf(B[k, t], observations[n])
        end
    end
end

function loglikelihoods!(LL::AbstractMatrix, B::AbstractArray{F,4} where {F<:UnivariateDistribution}, observations, n2t::AbstractArray{Int}, lag_cat::Matrix{Int})
    N, K, D = size(observations, 1), size(B, 1), size(observations, 2)
    @argcheck size(LL) == (N, K)

    for n in OneTo(N)
        t = n2t[n]
        for k in OneTo(K)
            LL[n, k] = logpdf(product_distribution(B[CartesianIndex.(k, t, 1:D, lag_cat[n, :])]), observations[n, :])
        end
    end
end

function fit_mle_A_from_slice!(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL, n2t; robust = false)

    N, K = size(LL)
    T = size(A, 3)
    c = zeros(N)
    γ = zeros(N, K)

    # get posterior of each category
    c[:] = logsumexp(LL, dims = 2)
    γ[:, :] = exp.(LL .- c)
    a[:] = γ[1, :] # initial date
    max_aposteriori = argmaxrow(γ)
    max_aposteriori_next = lead(max_aposteriori)

    Q = zeros(K, K, T)
    for n = 1:N
        t = n2t[n]
        if max_aposteriori_next[n] === missing
            continue
        else
            Q[max_aposteriori[n], max_aposteriori_next[n], t] += 1
        end
    end
    robust && (Q .+= eps())
    [A[:, :, t] = Q[:, :, t] ./ sum(Q[:, :, t], dims = 2) for t = 1:T]
end

function fit_mle_B_slice!(α::AbstractVector, B::AbstractVector{F} where {F<:Distribution}, observations;
    rand_ini = true,
    n_random_ini = 10, display_random = false,
    Dirichlet_α = 0.8, Dirichlet_categories = 0.85,
    ref_station = 1, kwargs...)
    if rand_ini == true
        α[:], B[:], h = fit_em_multiD_rand(α, B, observations; n_random_ini = n_random_ini, Dirichlet_α = Dirichlet_α, Dirichlet_categories = Dirichlet_categories, display_random = display_random, kwargs...)
    else
        h = fit_em_multiD!(α, B, observations; kwargs...)
    end
    sort_wrt_ref!(α, B, ref_station)
    h
end

function fit_mle_B_slice!(α::AbstractVector, B::AbstractArray{F,3} where {F<:Bernoulli},
    observations::Matrix{Int64}, lag_cat::Matrix{Int};
    rand_ini = true,
    n_random_ini = 10, display_random = false,
    Dirichlet_α = 0.8, Dirichlet_categories = 0.85,
    ref_station = 1, kwargs...)

    K, size_memory = size(B, 1), size(B, 3)
    idx_j = idx_observation_of_past_cat(lag_cat, K, size_memory)

    if rand_ini == true
        α[:], B[:], h = fit_em_multiD_rand(α, B, observations, lag_cat, idx_j; n_random_ini = n_random_ini, Dirichlet_α = Dirichlet_α, display_random = display_random, kwargs...)
    else
        h = fit_em_multiD!(α, B, observations, lag_cat, idx_j; kwargs...)
    end
    sort_wrt_ref!(α, B, ref_station)
    h
end

##

function sort_wrt_ref!(hmm::HierarchicalPeriodicHMM, ref_station)
    K, T = size(hmm.B, 1), size(hmm.B, 2)
    sorting = [[succprob(hmm.B[k, t, ref_station, 1]) for k = 1:K] for t = 1:T] # 1 is by convention the driest category i.e. Y|d....d
    new_order = sortperm.(sorting, rev = true)
    [hmm.B[:, t, :, :] = hmm.B[new_order[t], t, :, :] for t = 1:T]
    [hmm.A[:, :, t] = hmm.A[new_order[t], new_order[t], t] for t = 1:T]
end

function sort_wrt_ref!(hmm::PeriodicHMM, ref_station)
    # TODO
    K = size(α, 1)
    # sorting = [[succprob(hmm.B[k, t, ref_station, 1]) for k in 1:K] for t in 1:T] # 1 is by convention the driest category i.e. Y|d....d
    # [hmm.B[:,:,:,:] = hmm.B[sortperm(sorting[t], rev = true),t,:,:] for t in 1:T]
    # [hmm.A[:,:,t] = hmm.B[sortperm(sorting[t], rev = true), sortperm(sorting[t], rev = true),t] for t in 1:T]
end

function sort_wrt_ref!(α::AbstractVector, B::AbstractArray{F,3}, ref_station) where {F<:Bernoulli}
    K = size(α, 1)
    sorting = [succprob(B[k, ref_station, 1]) for k = 1:K] # 1 is by convention the driest category i.e. Y|d....d
    B[:, :, :] = B[sortperm(sorting, rev = true), :, :]
    α[:] = α[sortperm(sorting, rev = true)]
end

function sort_wrt_ref!(B::AbstractArray{F,3}, ref_station) where {F<:Bernoulli}
    K = size(B, 1)
    sorting = [succprob(B[k, ref_station, 1]) for k = 1:K] # 1 is by convention the driest category i.e. Y|d....d
    B[:, :, :] = B[sortperm(sorting, rev = true), :, :]
end

function sort_wrt_ref!(α::AbstractVector, B::AbstractVector{Product{Discrete,F,Vector{F}}}, ref_station) where {F<:Bernoulli}
    K = size(α, 1)
    sorting = [succprob(B[k].v[ref_station]) for k = 1:K]
    B[:] = B[sortperm(sorting, rev = true)]
    α[:] = α[sortperm(sorting, rev = true)]
end

function sort_wrt_ref!(B::AbstractVector{Product{Discrete,F,Vector{F}}}, ref_station) where {F<:Bernoulli}
    K = size(α, 1)
    sorting = [succprob(B[k].v[ref_station]) for k = 1:K]
    B[:] = B[sortperm(sorting, rev = true)]
end

function sort_wrt_ref!(α::AbstractVector, B::AbstractVector{Product{Discrete,F,Vector{F}}}, ref_station) where {F<:Categorical}
    K = size(α, 1)
    sorting = [probs(B[k].v[ref_station])[1] - 0.01 * probs(B[k].v[ref_station])[4] for k = 1:K] #add the extra term to remove ties

    B[:] = B[sortperm(sorting)]
    α[:] = α[sortperm(sorting)]
end

function sort_wrt_ref!(B::AbstractVector{Product{Discrete,F,Vector{F}}}, ref_station) where {F<:Categorical}
    K = size(α, 1)
    sorting = [probs(B[k].v[ref_station])[1] - 0.01 * probs(B[k].v[ref_station])[4] for k = 1:K] #add the extra term to remove ties
    B[:] = B[sortperm(sorting)]
end

##

function fit_em_multiD_rand(α::AbstractVector, B::AbstractArray{F,3} where {F<:Bernoulli},
    observations::Matrix{Int}, lag_cat::Matrix{Int}, idx_j::Vector{Vector{Vector{Int}}};
    n_random_ini = 10, Dirichlet_α = 0.8, display_random = false, kwargs...)

    D = size(observations, 2)
    K = size(α, 1)
    size_memory = size(B, 3)

    h = fit_em_multiD!(α, B, observations, lag_cat, idx_j; kwargs...)
    log_max = h.logtots[end]
    α_max, B_max = copy(α), copy(B)
    h_max = h
    (display_random == :iter) && println("random IC 1: logtot = $(h.logtots[end])")
    for i = 1:(n_random_ini-1)
        B[:, :, :] = random_product_Bernoulli(D, K, size_memory)
        α[:] = rand(Dirichlet(K, Dirichlet_α))
        h = fit_em_multiD!(α, B, observations, lag_cat, idx_j; kwargs...)
        (display_random == :iter) && println("random IC $(i+1): logtot = $(h.logtots[end])")
        if h.logtots[end] > log_max
            log_max = h.logtots[end]
            h_max = h
            α_max[:], B_max[:] = copy(α), copy(B)
        end
    end
    return α_max, B_max, h_max
end

function fit_em_multiD_rand(α::AbstractVector, B::AbstractVector{F} where {F<:Distribution},
    observations;
    n_random_ini = 10, Dirichlet_α = 0.8, Dirichlet_categories = 0.85, display_random = false, kwargs...)

    D = size(observations, 2)
    K = size(α, 1)

    h = fit_em_multiD!(α, B, observations; kwargs...)
    log_max = h.logtots[end]
    α_max, B_max = copy(α), copy(B)
    h_max = h
    (display_random == :iter) && println("random IC 1: logtot = $(h.logtots[end])")
    for i = 1:(n_random_ini-1)
        random_product!(B, Dirichlet_categories) # Randomize B
        α[:] = rand(Dirichlet(K, Dirichlet_α)) # Randomize α
        h = fit_em_multiD!(α, B, observations; kwargs...)
        (display_random == :iter) && println("random IC $(i+1): logtot = $(h.logtots[end])")
        if h.logtots[end] > log_max
            log_max = h.logtots[end]
            h_max = h
            α_max[:], B_max[:] = copy(α), copy(B)
        end
    end
    return α_max, B_max, h_max
end

function random_product!(B::AbstractVector{F}, Dirichlet_categories::Real) where {F<:Product{Discrete,U,Vector{U}}} where {U<:Categorical}
    K = size(B, 1)
    D = length(B[1].v)
    B[:] = random_product_Categorical(D, K, Dirichlet_categories)
end

function random_product!(B::AbstractVector{F}, Dirichlet_categories::Real) where {F<:Product{Discrete,U,Vector{U}}} where {U<:Bernoulli}
    K = size(B, 1)
    D = length(B[1].v)
    B[:] = random_product_Bernoulli(D, K)
end

function random_product_Categorical(D::Int, K::Int, Dirichlet_categories::Float64)
    [product_distribution([Categorical(rand(Dirichlet(K, Dirichlet_categories), D)[:, j]) for j = 1:D]) for k = 1:K]
end

function random_product_Bernoulli(D::Int, K::Int)
    [Product(Bernoulli.(rand(D))) for k = 1:K]
end

function random_product_Bernoulli(D::Int, K::Int, size_memory::Int)
    [Bernoulli(rand()) for k = 1:K, j = 1:D, m = 1:size_memory]
end

function idx_observation_of_past_cat(lag_cat, K, size_memory)
    # Matrix(T,D) of vector that give the index of data of same past.
    # ie. size_memory = 1 (no memory) -> every data is in category 1
    # ie size_memory = 2 (memory on previous day) -> idx_tj[t,j][1] = vector of index of data where previous day was dry, idx_tj[t,j][2] = index of data where previous day was wet
    D = size(lag_cat, 2)
    idx_j = Vector{Vector{Vector{Int}}}(undef, D)
    for i in OneTo(K)
        for j = 1:D
            idx_j[j] = [findall(lag_cat[:, j] .== m) for m = 1:size_memory]
            ##
        end
    end
    return idx_j
end

function fit_em_multiD!(α::AbstractVector, B::AbstractArray{F,3} where {F<:Bernoulli},
    observations::Matrix{Int}, lag_cat::Matrix{Int}, idx_j::Vector{Vector{Vector{Int}}};
    display = :none, maxiter = 100, tol = 1e-3, robust = false)

    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, D, size_memory = size(observations, 1), size(B, 1), size(B, 2), size(B, 3)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates

    LL = zeros(N, K)
    γ = similar(LL)
    c = zeros(N)

    # Initial parameters already in α, B

    # E-step
    # evaluate likelihood for each type k
    @inbounds for k in OneTo(K), n in OneTo(N)
        LL[n, k] = logpdf(product_distribution(B[CartesianIndex.(k, 1:D, lag_cat[n, :])]), observations[n, :])
    end
    [LL[:, k] .+= log(α[k]) for k = 1:K]
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
    # get posterior of each category
    c[:] = logsumexp(LL, dims = 2)
    γ[:, :] = exp.(LL .- c)

    # Loglikelihood
    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter

        # M-step
        # with γ in hand, maximize (update) the parameters
        α[:] = mean(γ, dims = 1)
        for k in OneTo(K)
            for j = 1:D
                for m = 1:size_memory
                    if sum(γ[idx_j[j][m], k]) > 0
                        B[k, j, m] = fit_mle(Bernoulli, observations[idx_j[j][m], j], γ[idx_j[j][m], k])
                    else
                        B[k, j, m] = Bernoulli(eps())
                    end
                end
            end
        end

        # E-step
        # evaluate likelihood for each type k
        @inbounds for k in OneTo(K), n in OneTo(N)
            LL[n, k] = logpdf(product_distribution(B[CartesianIndex.(k, 1:D, lag_cat[n, :])]), observations[n, :])
        end
        [LL[:, k] .+= log(α[k]) for k = 1:K]
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        # get posterior of each category
        c[:] = logsumexp(LL, dims = 2)
        γ[:, :] = exp.(LL .- c)

        # Loglikelihood
        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history
end

function fit_em_multiD!(α::AbstractVector, B::AbstractVector{F} where {F<:Distribution},
    observations;
    display = :none, maxiter = 200, tol = 1e-3, robust = false)

    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K = size(observations, 1), size(B, 1)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates

    LL = zeros(N, K)
    γ = similar(LL)
    c = zeros(N)

    # Initial parameters already in α, B

    # E-step
    # evaluate likelihood for each type k
    [LL[:, k] = logpdf(B[k], permutedims(observations)) .+ log(α[k]) for k = 1:K]
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
    # get posterior of each category
    c[:] = logsumexp(LL, dims = 2)
    γ[:, :] = exp.(LL .- c)

    # Loglikelihood
    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    # ensure that for categorical data the correct number of categorie is taken. Else nothing
    if typeof(B[1].v[1]) == Categorical{Float64,Vector{Float64}}
        C = (ncategories(B[1].v[1]),)
    else
        C = ()
    end

    for it = 1:maxiter

        # M-step
        # with γ in hand, maximize (update) the parameters
        α[:] = mean(γ, dims = 1)
        for k in OneTo(K)
            if sum(γ[:, k]) > 0
                B[k] = fit_mle(B[k], C..., permutedims(observations), γ[:, k])
            end
        end

        # E-step
        # evaluate likelihood for each type k
        [LL[:, k] = logpdf(B[k], permutedims(observations)) .+ log(α[k]) for k = 1:K]
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        # get posterior of each category
        c[:] = logsumexp(LL, dims = 2)
        γ[:, :] = exp.(LL .- c)

        # Loglikelihood
        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history
end


text:
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../test/useful_routines.jl")

function cell_interface_error(quad, cellmap, facescale, distance_to_interface)
    err = 0.0
    for (idx, (p, w)) in enumerate(quad)
        x = cellmap(p)
        err += (distance_to_interface(x)[1])^2 * facescale[idx] * w
    end
    return err
end

function interface_error(levelset, levelsetcoeffs, cutmesh, numqp, distance_to_interface)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    cellids = findall(CutCell.cell_sign(cutmesh) .== 0)

    err = 0.0
    for cellid in cellids
        quad = interfacequads[+1, cellid]
        if length(quad) > 0
            normals = CutCell.interface_normals(interfacequads, cellid)
            cellmap = CutCell.cell_map(cutmesh, cellid)
            facescale = CutCell.scale_area(cellmap, normals)

            err += cell_interface_error(quad, cellmap, facescale, distance_to_interface)
        end
    end
    return sqrt(err)
end

function shrinking_corner_distance(x,xc,speed,time)
    dx = speed*time
    return corner_distance_function(x,xc-[dx,dx])
end

function final_interface_error(
    levelset,
    mesh,
    speed,
    dt,
    nsteps,
    numqp,
    xc,
)

    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> shrinking_corner_distance(x, xc, speed, 0.0),
        mesh,
    )
    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    levelsetspeed = speed * ones(length(levelsetcoeffs))

    for t = 1:nsteps
        paddedmesh = CutCell.BoundaryPaddedMesh(cutmesh, 1)
        refseedpoints, spatialseedpoints, seedcellids =
            CutCell.seed_zero_levelset(2, levelset, levelsetcoeffs, cutmesh)
        paddedlevelset = CutCell.BoundaryPaddedLevelset(
            paddedmesh,
            refseedpoints,
            spatialseedpoints,
            seedcellids,
            levelset,
            levelsetcoeffs,
            cutmesh,
            1e-10,
        )

        levelsetcoeffs =
            CutCell.step_first_order_levelset(paddedlevelset, levelsetspeed, dt)
        cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    end
    dx = speed*nsteps*dt
    cornerlocation = xc - [dx,dx]
    exact_interface_length = sum(cornerlocation)

    err = interface_error(
        levelset,
        levelsetcoeffs,
        cutmesh,
        numqp,
        x -> shrinking_corner_distance(x,xc,speed,nsteps*dt),
    )/exact_interface_length

    return err
end

function grid_size(mesh)
    w = CutCell.widths(mesh)
    nn = CutCell.nodes_per_mesh_side(mesh)
    return w ./ (nn .- 1)
end

function time_step_size(speed, mesh; CFL = 0.5)
    dx = minimum(grid_size(mesh))
    return CFL * dx / abs(speed)
end

function numsteps(stoptime,dt)
    r = stoptime/dt
    @assert isinteger(r)
    return round(Int,r)
end

function error_for_numelmts(nelmts)
    x0 = [0.0, 0.0]
    L, W = 1.0, 1.0
    numghostlayers = 1
    polyorder = 2
    numqp = 5

    xc = [0.87, 0.87]
    speed = 1.0
    stoptime = 0.2

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)

    mesh = CutCell.Mesh(x0, [L, W], [nelmts, nelmts], basis)
    dt = time_step_size(speed, mesh)
    nsteps = numsteps(stoptime,dt)

    return final_interface_error(
        levelset,
        mesh,
        speed,
        dt,
        nsteps,
        numqp,
        xc,
    )
end

nelmts = [5,10,20,40,80]
# err = error_for_numelmts(80)
dx = 1.0 ./ nelmts
err = [error_for_numelmts(ne) for ne in nelmts]

# using CSV,DataFrames
# df = DataFrame([dx,err],["Element Size", "Error"])
# CSV.write("examples/levelset-propagate/corner-convergence.csv",df)


text:
"""
    function initiate!(pop::simple_sites, c::Int64; overwrite=false)

Iniate a population of random sequences.

# Arguments
- `pop::simple_sites`: population that is to be filled.
- `c::Int64=1`: number of species to create.
"""
function initiate!(pop::simple_sites, c::Int64=1)
    # Create new population
    N_sub = pop.N ÷ c
    rest = pop.N - N_sub * c
    for i in 1:c
        push!(pop.seqs, rand(collect(1:pop.n), pop.l))
        push!(pop.freqs, N_sub)
    end
    pop.freqs[1] += rest
    nothing
end

"""
    function initiate!(pop::simple_sites, c::Int64; overwrite=false)

Iniate a population of random sequences.

# Arguments
- `pop::simple_sites`: population that is to be filled.
- `c::Int64=1`: number of species to create.
"""
function initiate!(pop::simple_sites, emat::Array{Float64, 2})
    if size(emat) != (pop.n, pop.l)
        throw(DimensionMismatch("Energy matrix needs to have dimensions (n, l)"))
    end
    push!(pop.seqs, argmin.(emat[:, i] for i in 1:pop.l))
    push!(pop.freqs, pop.N)
    nothing
end

"""
    function initiate!(
        pop::driver_trailer,
        c::Int64;
        driver::Array{Int64, 1}=Int64[]
        )

Initiates population.

Creates a population of type `driver_trailer`, where a driver sequence can be given.
Given driver sequence can be shorter or longer than initial binding site length. Sequence
will be either filled with random letters or truncated to the right length.

# Arguments
- `pop::driver_trailer`: population that is to be filled.
- `c::Int64=1`: number of species to create.
- `driver::Array{Int64, 1}=Int64[]`: Optional initial driver sequence.
"""
function initiate!(
    pop::driver_trailer,
    c::Int64=1;
    driver::Array{Int64, 1}=Int64[]
    )

    pop.driver[:] = make_driver(driver, pop.m, pop.L)
    # Create new population
    N_sub = pop.N ÷ c
    rest = pop.N - N_sub * c
    for i in 1:c
        push!(pop.seqs, rand(collect(1:pop.n), pop.L))
        push!(pop.freqs, N_sub)
    end
    pop.freqs[1] += rest
    nothing
end

"""
    function initiate!(
        pop::driver_trailer,
        emat::Array{Float64, 2};
        driver::Array{Int64, 1}=Int64[]
        )

Initiates an optimal population.

Creates a population of type `driver_trailer`, where a driver sequence can be given.
Given driver sequence can be shorter or longer than initial binding site length. Optimal sequence
will chosen as the one which minimizes energy.

# Arguments
- `pop::driver_trailer`: population that is to be filled.
- `emat::Array{Float64,2}` : Energy matrix to choose optimal sequence
- `driver::Array{Int64, 1}=Int64[]`: Optional initial driver sequence.
"""
function initiate!(
    pop::driver_trailer,
    emat::Array{Float64, 2};
    driver::Array{Int64, 1}=Int64[]
    )

    if size(emat) != (pop.n, pop.m)
        throw(DimensionMismatch("Energy matrix needs to have dimensions (n, m)"))
    end

    pop.driver[:] = make_driver(driver, pop.m, pop.L)
    push!(pop.freqs, pop.N)

    push!(pop.seqs, argmin.([emat[:, j] for j in pop.driver]))

end

"""
    function initiate!(
        pop::driver_trailer_l,
        c::Int64;
        driver::Array{Int64, 1}=Int64[]
        )

Iniate a population of c species of and a driver sequence, which can be given.

Creates a population of type `driver_trailerPl`, where a driver sequence can be given.
Given driver sequence can be shorter or longer than initial binding site length. Sequence
will be either filled with random letters or truncated to the right length.

# Arguments
- `pop::driver_trailer`: population that is to be filled.
- `c::Int64=1`: number of species to create.
- `driver::Array{Int64, 1}=Int64[]`: Optional initial driver sequence.
"""
function initiate!(
    pop::driver_trailer_l,
    c::Int64=1;
    driver::Array{Int64, 1}=Int64[],
    )

    pop.driver[:] = make_driver(driver, pop.m, pop.L)
    N_sub = pop.N ÷ c
    rest = pop.N - N_sub * c
    for i in 1:c
        push!(pop.seqs, rand(collect(1:pop.n), pop.L))
        push!(pop.freqs, N_sub)
    end
    pop.freqs[1] += rest

    pop.l = ones(Int64, c) * pop.l_0
end

"""
    function initiate!(
        pop::driver_trailer_l,
        emat::Array{Float64, 2};
        driver::Array{Int64, 1}=Int64[]
        )

Iniate a population of c species of and a driver sequence, which can be given.

Creates a population of type `driver_trailerPl`, where a driver sequence can be given.
Given driver sequence can be shorter or longer than initial binding site length. Sequence
will be either filled with random letters or truncated to the right length.

# Arguments
- `pop::driver_trailer`: population that is to be filled.
- `emat::Array{Float64,2}` : Energy matrix to choose optimal sequence
- `driver::Array{Int64, 1}=Int64[]`: Optional initial driver sequence.
"""
function initiate!(
    pop::driver_trailer_l,
    emat::Array{Float64,2};
    driver::Array{Int64, 1}=Int64[],
    )

    pop.driver[:] = make_driver(driver, pop.m, pop.L)
    push!(pop.freqs, pop.N)
    push!(pop.seqs, argmin.([emat[:, j] for j in pop.driver]))

    pop.l = Int64[pop.l_0]
end

"""
    function initiate!(
        pop::mono_pop;
        driver::Array{Int64, 1}=Int64[]
        )

Iniate a monomorphic population and a driver sequence, which can be given.

Creates a population of type `mono_pop`, where a driver sequence can be given.
Given driver sequence can be shorter or longer than initial binding site length. Sequence
will be either filled with random letters or truncated to the right length.

# Arguments
- `pop::driver_trailer`: population that is to be filled.
- `driver::Array{Int64, 1}=Int64[]`: Optional initial driver sequence.
"""
function initiate!(
    pop::mono_pop;
    driver::Array{Int64, 1}=Int64[],
    )
    pop.driver = make_driver(driver, pop.m, pop.l)
    pop.seqs =  rand(collect(1:pop.n), pop.l)
    nothing
end

"""
    function initiate!(
        pop::mono_pop
        emat::Array{Float64, 2};
        driver::Array{Int64, 1}=Int64[]
        )

Iniate a monomorphic population and a driver sequence, which can be given.

Creates a population of type `mono_pop`, where a driver sequence can be given.
Given driver sequence can be shorter or longer than initial binding site length. Sequence
will be either filled with random letters or truncated to the right length.

# Arguments
- `pop::driver_trailer`: population that is to be filled.
- `emat::Array{Float64,2}` : Energy matrix to choose optimal sequence
- `driver::Array{Int64, 1}=Int64[]`: Optional initial driver sequence.
"""
function initiate!(
    pop::mono_pop,
    emat::Array{Float64,2};
    driver::Array{Int64, 1}=Int64[],
    )
    if size(emat) != (pop.n, pop.m)
        throw(DimensionMismatch("Energy matrix needs to have dimensions (n, m)"))
    end
    pop.driver = make_driver(driver, pop.m, pop.l)
    pop.seqs =  argmin.([emat[:, j] for j in pop.driver])
    nothing
end

"""
    make_driver(given_driver, m, L)

Create a driver sequence.
"""
function make_driver(given_driver, m, L)
    driver = zeros(Int64, L)
    if isempty(given_driver)
        driver[:] = rand(collect(1:m), L)
    else
        if any(given_driver .> m)
            throw(ArgumentError("Driver sequence has elements outside of alphabet."))
        elseif any(given_driver .< 0)
            throw(ArgumentError("Don't include negative elements in driver sequence."))
        end
        if length(given_driver) == L
            driver[:] = given_driver
        elseif length(given_driver) < L
            driver[1:length(given_driver)] = given_driver
            driver[length(given_driver) + 1:end] = rand(collect(1:m), L - length(given_driver))
        else
            driver[:] = given_driver[1:L]
        end
    end
    return driver
end


text:
@with_kw mutable struct EpisodicSafetyGym{S, A, O} <: POMDP{S,A,O}
    pomdp::POMDP{S,A,O}
    λ::Float32 = 0f0
    terminate_next::Bool = false
    isfailure = false
end

POMDPs.initialstate(mdp::EpisodicSafetyGym) = initialstate(mdp.pomdp)
POMDPs.initialobs(mdp::EpisodicSafetyGym, s) = initialobs(mdp.pomdp, s)

POMDPs.actions(mdp::EpisodicSafetyGym) = actions(mdp.pomdp)
POMDPs.actionindex(mdp::EpisodicSafetyGym, a) = actionindex(mdp.pomdp)

function POMDPs.isterminal(mdp::EpisodicSafetyGym, s)
    if mdp.terminate_next
        mdp.terminate_next = false
        return true
    else
        isterminal(mdp.pomdp, s)
    end
end
POMDPs.discount(mdp::EpisodicSafetyGym) = discount(mdp.pomdp)

POMDPs.observation(mdp::EpisodicSafetyGym, s) = observation(mdp.pomdp, s)

render(mdp::EpisodicSafetyGym, s, a = nothing; kwargs...) = render(mdp.pomdp, s, a; kwargs...)
render(mdp::EpisodicSafetyGym; kwargs...) = render(mdp.pomdp; kwargs...)

isfailure(mdp::EpisodicSafetyGym, s) = mdp.isfailure

function POMDPs.gen(mdp::EpisodicSafetyGym, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    mdp.isfailure = false
    sp, o, r = gen(mdp.pomdp, s, a, rng, info=info)
    if r > 0.5 # Reached a goal
        mdp.terminate_next = true
    end
    if haskey(info, "cost") && info["cost"] >0
        r -= mdp.λ*info["cost"]
        mdp.terminate_next = true
        mdp.isfailure = true
    end
    return (sp=sp, o=o, r=r)
end

max_forks_repo_name: ancorso/POMDPGym


text:
# Autogenerated wrapper script for muparser_jll for aarch64-apple-darwin
export libmuparser

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("muparser")
JLLWrappers.@declare_library_product(libmuparser, "@rpath/libmuparser.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libmuparser,
        "lib/libmuparser.2.3.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Device query/info/open/close

struct LJDevice
    handle::Cint
end

function lsdev(dt::LJDeviceType=dtANY, ct::LJConnectionType=ctANY; quiet = false)

    found = Ref{Int32}()
    devtypes = Vector{Int32}(undef, LJM.LIST_ALL_SIZE)
    ctypes = Vector{Int32}(undef, LJM.LIST_ALL_SIZE)
    serials = Vector{Int32}(undef, LJM.LIST_ALL_SIZE)
    ips = Vector{Int32}(undef, LJM.LIST_ALL_SIZE)

    LJM.ListAll(dt,ct,found,devtypes,ctypes,serials,ips) |> errorcheck

    found[] == 0 && return nothing
    quiet && return found[]

    for i in 1:found[]
        println("\n[Device #", i, "]")
        println("\tDevice type:     ", LJDeviceType(devtypes[i]))
        println("\tConnection type: ", LJConnectionType(ctypes[i]))
        println("\tSerial number:   ", serials[i])
        println("\tIP Address:      ", ips[i])
    end
end

function Base.open(dt::LJDeviceType, ct::LJConnectionType)
    handle = Ref{Cint}()
    LJM.Open(dt, ct, "ANY", handle) |> errorcheck
    return LJDevice(handle[])
end

function Base.close(dev::LJDevice)
    LJM.Close(dev.handle) |> errorcheck
    dev = nothing
    return
end


text:
include("AbstractTypes/AbstractType.jl")

mutable struct Load <: PassiveAgent
    priceCap::Real # EUR/MWh
    hourlyLoad::Real # load at a specific hour
    representativeLoad::Matrix{Float64} #load factor at each hour in a day

    TimeSeries::Function

    function Load(priceCap::Real,hourlyLoad::Real,representativeLoad::Matrix{Float64})
        self = new()
        self.priceCap = priceCap
        self.hourlyLoad = hourlyLoad
        self.representativeLoad = representativeLoad

        # This function changes load along with hour
        self.TimeSeries = function(day::Int64,hour::Int64)
            self.hourlyLoad = self.representativeLoad[day,hour]
        end
        return self
    end
end


text:
mutable struct MCEuropeanEngine{S <: AbstractBlackScholesProcess, RSG <: AbstractRandomSequenceGenerator} <: MCVanillaEngine{S, RSG}
  process::S
  timeSteps::Int
  timeStepsPerYear::Int
  requiredSamples::Int
  maxSamples::Int
  requiredTolerance::Float64
  brownianBridge::Bool
  seed::Int
  # mcSimulation::MCSimulation{RSG, T}
  antitheticVariate::Bool
  rsg::RSG
end

function MCEuropeanEngine(process::S; timeSteps::Int = -1, timeStepsPerYear::Int = -1, brownianBridge::Bool = false, antitheticVariate::Bool = false,
                          requiredSamples::Int = -1, requiredTolerance::Float64 = -1.0, maxSamples::Int = typemax(Int), seed::Int = 1,
                          rsg::RSG = InverseCumulativeRSG(seed)) where {S <: AbstractBlackScholesProcess, RSG <: AbstractRandomSequenceGenerator}
  # build mc sim
  # mcSim = MCSimulation{RSG, SingleVariate}(antitheticVariate, false, rsg, SingleVariate())

  return MCEuropeanEngine{S, RSG}(process, timeSteps, timeStepsPerYear, requiredSamples, maxSamples, requiredTolerance, brownianBridge, seed, antitheticVariate, rsg)
end

struct EuropeanPathPricer{P <: PlainVanillaPayoff} <: AbstractPathPricer
  payoff::P
  discount::Float64
end

EuropeanPathPricer(optionType::OT, strike::Float64, disc::Float64) where {OT <: OptionType} = EuropeanPathPricer(PlainVanillaPayoff{OT}(optionType, strike), disc)

(pricer::EuropeanPathPricer)(path::Path) = pricer.payoff(path[end]) * pricer.discount

function time_grid(pe::MCVanillaEngine, opt::VanillaOption)
  lastExerciseDate = opt.exercise.dates[end]
  t = get_time(pe.process, lastExerciseDate)
  if pe.timeSteps != -1
    return TimeGrid(t, pe.timeSteps)
  elseif pe.timeStepsPerYear != -1
    steps = round(Int, floor(pe.timeStepsPerYear * t))
    return TimeGrid(t, max(steps, 1))
  else
    error("time steps not specified")
  end
end

function path_generator(pe::MCVanillaEngine, opt::VanillaOption)
  dimensions = get_factors(pe.process)
  grid = time_grid(pe, opt)
  init_sequence_generator!(pe.rsg, dimensions * (length(grid.times) - 1))

  return PathGenerator(pe.process, grid, pe.rsg, pe.brownianBridge)
end

function path_pricer(pe::MCEuropeanEngine, opt::VanillaOption)
  payoff = opt.payoff
  process = pe.process

  return EuropeanPathPricer(payoff.optionType, payoff.strike, discount(process.riskFreeRate, time_grid(pe, opt)[end]))
end

function _calculate!(pe::MCVanillaEngine, opt::VanillaOption)
  mcSim = MCSimulation(pe, false, opt, SingleVariate())
  _calculate!(mcSim, pe.requiredTolerance, pe.requiredSamples, pe.maxSamples)
  opt.results.value = stats_mean(mcSim.mcModel.sampleAccumulator)
  return pe, opt
end


text:
"""
    hfk5z(rh, dh, date1, date2)

Transform a Hipparcos star position into FK5 J2000.0, assuming
zero Hipparcos proper motion.

### Given ###

- `rh`: Hipparcos RA (radians)
- `dh`: Hipparcos Dec (radians)
- `date1`, `date2`: TDB date (Note 1)

### Returned (all FK5, equinox J2000.0, date date1+date2) ###

- `r5`: RA (radians)
- `d5`: Dec (radians)
- `dr5`: FK5 RA proper motion (rad/year, Note 4)
- `dd5`: Dec proper motion (rad/year, Note 4)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

3. The FK5 to Hipparcos transformation is modeled as a pure rotation
   and spin;  zonal errors in the FK5 catalogue are not taken into
   account.

4. It was the intention that Hipparcos should be a close
   approximation to an inertial frame, so that distant objects have
   zero proper motion;  such objects have (in general) non-zero
   proper motion in FK5, and this function returns those fictitious
   proper motions.

5. The position returned by this function is in the FK5 J2000.0
   reference system but at date date1+date2.

6. See also eraFk52h, eraH2fk5, eraFk5zhz.

### Called ###

- `eraS2c`: spherical coordinates to unit vector
- `eraFk5hip`: FK5 to Hipparcos rotation and spin
- `eraRxp`: product of r-matrix and p-vector
- `eraSxp`: multiply p-vector by scalar
- `eraRxr`: product of two r-matrices
- `eraTrxp`: product of transpose of r-matrix and p-vector
- `eraPxp`: vector product of two p-vectors
- `eraPv2s`: pv-vector to spherical
- `eraAnp`: normalize angle into range 0 to 2pi

### Reference ###

- F.Mignard & M.Froeschle, 2000, Astron.Astrophys. 354, 732-739.

"""
function hfk5z(rh, dh, date1, date2)
    r5 = Ref(0.0)
    d5 = Ref(0.0)
    dr5 = Ref(0.0)
    dd5 = Ref(0.0)
    ccall((:eraHfk5z, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
          rh, dh, date1, date2, r5, d5, dr5, dd5)
    r5[], d5[], dr5[], dd5[]
end

"""
    h2fk5(ra, dec, dra, ddec, px, rv)

Transform Hipparcos star data into the FK5 (J2000.0) system.

### Given (all Hipparcos, epoch J2000.0) ###

- `rh`: RA (radians)
- `dh`: Dec (radians)
- `drh`: Proper motion in RA (dRA/dt, rad/Jyear)
- `ddh`: Proper motion in Dec (dDec/dt, rad/Jyear)
- `pxh`: Parallax (arcsec)
- `rvh`: Radial velocity (km/s, positive = receding)

### Returned (all FK5, equinox J2000.0, epoch J2000.0) ###

- `r5`: RA (radians)
- `d5`: Dec (radians)
- `dr5`: Proper motion in RA (dRA/dt, rad/Jyear)
- `dd5`: Proper motion in Dec (dDec/dt, rad/Jyear)
- `px5`: Parallax (arcsec)
- `rv5`: Radial velocity (km/s, positive = receding)

### Notes ###

1. This function transforms Hipparcos star positions and proper
   motions into FK5 J2000.0.

2. The proper motions in RA are dRA/dt rather than
   cos(Dec)*dRA/dt, and are per year rather than per century.

3. The FK5 to Hipparcos transformation is modeled as a pure
   rotation and spin;  zonal errors in the FK5 catalog are not
   taken into account.

4. See also eraFk52h, eraFk5hz, eraHfk5z.

### Called ###

- `eraStarpv`: star catalog data to space motion pv-vector
- `eraFk5hip`: FK5 to Hipparcos rotation and spin
- `eraRv2m`: r-vector to r-matrix
- `eraRxp`: product of r-matrix and p-vector
- `eraTrxp`: product of transpose of r-matrix and p-vector
- `eraPxp`: vector product of two p-vectors
- `eraPmp`: p-vector minus p-vector
- `eraPvstar`: space motion pv-vector to star catalog data

### Reference ###

- F.Mignard & M.Froeschle, Astron. Astrophys. 354, 732-739 (2000).

"""
function h2fk5(ra, dec, dra, ddec, px, rv)
    r = Ref(0.0)
    d = Ref(0.0)
    dr = Ref(0.0)
    dd = Ref(0.0)
    p = Ref(0.0)
    v = Ref(0.0)
    ccall((:eraH2fk5, liberfa), Cvoid,
            (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
            ra, dec, dra, ddec, px, rv, r, d, dr, dd, p, v)
    r[], d[], dr[], dd[], p[], v[]
end


text:
using LinearAlgebra

export OrthoMADS

"""
    OrthoMADS()

Return an empty OrthoMADS object.

OrthoMADS implements the poll set generation as described in Audet and Le Digabel 2015
Section 3.4.
"""
struct OrthoMADS <: AbstractPoll end

"""
    GenerateDirections(p::AbstractProblem, o::OrthoMADS)::Matrix

Generates the poll directions, using the poll set generation
as described in Audet and Le Digabel 2015 Section 3.4.
"""
(GenerateDirections(p::AbstractProblem, o::OrthoMADS)::Matrix) =
    GenerateDirections(p.N, o)
function GenerateDirections(N::Int64, ::OrthoMADS)::Matrix
    dirs_on_unit_sphere = GenerateDirectionsOnUnitSphere(N)
    H = HouseholderTransform(dirs_on_unit_sphere)
	return hcat(H, -H)
end

"""
    GenerateDirectionsOnUnitSphere(N::Int64)

Return an `N` length normalized direction on the unit sphere.
"""
function GenerateDirectionsOnUnitSphere(N::Int64)
    dir = randn(N)
    return dir ./ norm(dir)
end

"""
    HouseholderTransform(q)

Apply the Householder transformation to the vector `q`.
"""
function HouseholderTransform(q)
    nq = norm(q)
    v = q./nq
    return nq^2 .* (I - 2*v*v')
end


text:
module testMoments

using Base.Test
using DataFrames
using Datetime
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running moments tests\n")

#########################
## test portfolio mean ##
#########################

pf = AssetMgmt.Portfolio(ones(4, 8)/8)
mus = DataFrame(rand(1, 8))

## test getPMean for portfolios
AssetMgmt.getPMean(pf, mus)

## test getPMean for investments
vals = rand(8, 4)
valsDf = DataFrame(AssetMgmt.makeWeights(vals))
invs = AssetMgmt.Investments(valsDf, [1:8])

mus = DataFrame(rand(8, 4))

AssetMgmt.getPMean(invs, mus)

#############################
## test portfolio variance ##
#############################

pf = AssetMgmt.Portfolio(ones(1, 4)/4)
covMatr = cov(Timematr(rand(50, 4)))
AssetMgmt.getPVar(pf, covMatr)

valsDf = DataFrame(ones(10, 4)/4)
invs = AssetMgmt.Investments(valsDf, [1:10])
AssetMgmt.getPVar(invs, covMatr)

end


text:
module TestPackage

# Write your package code here.
function testFunction(number)
    return 2*number
end

export testFunction

end


text:
using Turf

fc = GeoJSON.parsefile(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/nearest-point/points.geojson")

function get_point_with_property(features, prop)
    point = nothing
    for (i, feature) in enumerate(features)
        if haskey(feature.properties, prop)
            point = features[i]
            deleteat!(features, i)
        end
    end
    return point
end

destination = get_point_with_property(fc.features, "tag")

nearest = nearestpoint(destination.geometry, [x.geometry for x in fc.features])

index = findfirst(x -> x.geometry.coordinates == nearest.coordinates, fc.features)

fc.features[index].properties = Dict("marker-color" => "#fbb000", "tag" => "nearest-point") # yellow marker
push!(fc.features, destination)

result = geojson(fc)

open(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/nearest-point/points.result.geojson", "w") do file
    write(file, result)
end


text:
using Test
using LinearAlgebra
using DataStructure101
const DS = DataStructure101

@testset "Graph" begin
    @testset "Di Adjacency Matrix" begin
        @testset "Basic operate" begin
            n = 10

            g = DS.Graph(n, :di_matrix)
            row_str = "\t false" ^ n
            data_str = "$row_str\n" ^ n
            @test repr(g) == "Graph(\n$data_str)"

            @test g isa DS.AbstractGraph
            @test g isa DS.DiAdjacencyMatrix
            @test DS.nv(g) == n
            @test DS.ne(g) == 0

            DS.relate!(g, 2, 5)
            @test DS.ne(g) == 1
            @test DS.neighbor(g, 2) == Int[5]

            g = DS.Graph(n, :w_di_matrix, random=false)

            row_str = "\t 0.0" ^ n
            data_str = "$row_str\n" ^ n
            @test repr(g) == "Graph(\n$data_str)"

            @test g isa DS.AbstractGraph
            @test g isa DS.WeightedDiAdjacencyMatrix
            @test DS.nv(g) == n
            @test DS.ne(g) == 0

            DS.relate!(g, 2, 5, 10.5)
            @test DS.ne(g) == 1
            @test DS.neighbor(g, 2) == Int[5]

            w = zeros(Float64, n, n)
            w[2, 5] = 10.5
            @test DS.weight(g) == w
        end

        @testset "Random Walk" begin
            n = 100
            city = DS.Graph(n, :w_di_matrix)
            p = DS.probability(city)

            test_prob = []
            for i=1:n push!(test_prob, sum(p[i, 1:end]) ≈ 1.0) end
            @test all(test_prob)

            # Random Walk Begin
            x = collect(1.0:n)
            x = DS.random_walk(city, x, 10000)
            @test p' * x ≈ x
        end

    end

    @testset "Di AdjacencyList" begin
        g = DS.Graph(10, :di_list)
        @test g isa DS.AbstractGraph
        @test g isa DS.DiAdjacencyList

        @test g[8] == 8
        @test g[11] == -1

        DS.relate!(g, 10, 9)
        DS.relate!(g, 10, 8)

        DS.relate!(g, 9, 7)
        DS.relate!(g, 9, 6)
        DS.relate!(g, 8, 5)
        DS.relate!(g, 8, 4)

        DS.relate!(g, 7, 3)
        DS.relate!(g, 7, 2)
        DS.relate!(g, 6, 1)

        @test DS.nv(g) == 10
        @test DS.ne(g) == 9
        @test DS.degree(g, 8) == 2
        @test DS.degree(g, 6) == 1

        @test DS.neighbor(g, 10) == [9, 8]

        @test repr(g) == "DiAdjacencyList(" *
            "1[], 2[], 3[], 4[], 5[], " *
            "6[1, ], 7[3, 2, ], 8[5, 4, ], 9[7, 6, ], 10[9, 8, ], )"

        DS.unrelate!(g, 6, 1)
        @test repr(g) == "DiAdjacencyList(" *
            "1[], 2[], 3[], 4[], 5[], " *
            "6[], 7[3, 2, ], 8[5, 4, ], 9[7, 6, ], 10[9, 8, ], )"
    end

    @testset "WeightedAdjacencyList" begin
        g = DS.Graph(10, :w_list)
        @test g isa DS.WeightedAdjacencyList
        @test g isa DS.AbstractGraph

        @test g[8] == 8
        @test g[11] == -1

        DS.relate!(g, 10, 9, 9.0)
        DS.relate!(g, 10, 8, 8.0)

        DS.relate!(g, 9, 7, 7.0)
        DS.relate!(g, 9, 6, 6.0)
        DS.relate!(g, 8, 5, 5.0)
        DS.relate!(g, 8, 4, 4.0)

        DS.relate!(g, 7, 3, 3.0)
        DS.relate!(g, 7, 2, 2.0)
        DS.relate!(g, 6, 1, 1.0)

        @test DS.nv(g) == 10
        @test DS.ne(g) == 9
        @test DS.degree(g, 8) == 17
        @test DS.degree(g, 6) == 7

        @test DS.neighbor(g, 10) == [9, 8]

        @test repr(g) == "WeightedAdjacencyList{Float64}(" *
            "1[6(1.0), ], 2[7(2.0), ], 3[7(3.0), ], 4[8(4.0), ], " *
            "5[8(5.0), ], 6[9(6.0), 1(1.0), ], " *
            "7[9(7.0), 3(3.0), 2(2.0), ], 8[10(8.0), 5(5.0), 4(4.0), ], " *
            "9[10(9.0), 7(7.0), 6(6.0), ], 10[9(9.0), 8(8.0), ], )"

        DS.unrelate!(g, 6, 1)
        @test DS.ne(g) == 8
        @test DS.degree(g, 1) == 0
        @test DS.degree(g, 6) == 6
        @test repr(g) == "WeightedAdjacencyList{Float64}(" *
            "1[], 2[7(2.0), ], 3[7(3.0), ], 4[8(4.0), ], " *
            "5[8(5.0), ], 6[9(6.0), ], " *
            "7[9(7.0), 3(3.0), 2(2.0), ], 8[10(8.0), 5(5.0), 4(4.0), ], " *
            "9[10(9.0), 7(7.0), 6(6.0), ], 10[9(9.0), 8(8.0), ], )"
    end
end


text:
macro checked_lib(libname, path)
        (dlopen_e(path) == C_NULL) && error("Unable to load \n\n$libname ($path)\n\nPlease re-run Pkg.build(package), and restart Julia.")
        quote const $(esc(libname)) = $path end
    end
@checked_lib libsundials_nvecserial "/usr/local/lib/libsundials_nvecserial.dylib"
@checked_lib libsundials_cvode "/usr/local/lib/libsundials_cvode.dylib"
@checked_lib libsundials_cvodes "/usr/local/lib/libsundials_cvodes.dylib"
@checked_lib libsundials_kinsol "/usr/local/lib/libsundials_kinsol.dylib"
@checked_lib libsundials_ida "/usr/local/lib/libsundials_ida.dylib"
@checked_lib libsundials_idas "/usr/local/lib/libsundials_idas.dylib"


text:
# pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
# const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
# includet(joinpath(LOOPVECBENCHDIR, "driver.jl"))

pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
include(joinpath(LOOPVECBENCHDIR, "benchmarkflops.jl"))
include(joinpath(LOOPVECBENCHDIR, "plotbenchmarks.jl"))

using Distributed

addprocs(9);

@everywhere begin
    pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
    const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
    include(joinpath(LOOPVECBENCHDIR, "benchmarkflops.jl"))
end

gemm_future = @spawnat 2 benchmark_gemm(2:256);
AtmulB_future = @spawnat 3 benchmark_AtmulB(2:256);
dot_future = @spawnat 4 benchmark_dot(2:256);
selfdot_future = @spawnat 5 benchmark_selfdot(2:256);
gemv_future = @spawnat 6 benchmark_gemv(2:256);
dot3_future = @spawnat 7 benchmark_dot3(2:256);
sse_future = @spawnat 8 benchmark_sse(2:256);
exp_future = @spawnat 9 benchmark_exp(2:256);
aplusBc_future = @spawnat 10 benchmark_aplusBc(2:256);

gemm_bench = fetch(gemm_future)
AtmulB_bench = fetch(AtmulB_future)
dot_bench = fetch(dot_future)
selfdot_bench = fetch(selfdot_future)
gemv_bench = fetch(gemv_future)
dot3_bench = fetch(dot3_future)
sse_bench = fetch(sse_future)
exp_bench = fetch(exp_future)
aplusBc_bench = fetch(aplusBc_future)

plot(gemm_bench)
plot(AtmulB_bench)
plot(dot_bench)
plot(selfdot_bench)
plot(gemv_bench)
plot(dot3_bench)
plot(sse_bench)
plot(exp_bench)
plot(aplusBc_bench)

save(joinpath("~/Pictures", "bench_gemm_v3.png"), plot(gemm_bench));
save(joinpath("~/Pictures", "bench_AtmulB_v3.png"), plot(AtmulB_bench));
save(joinpath("~/Pictures", "bench_dot_v3.png"), plot(dot_bench));
save(joinpath("~/Pictures", "bench_selfdot_v3.png"), plot(selfdot_bench));
save(joinpath("~/Pictures", "bench_gemv_v3.png"), plot(gemv_bench));
save(joinpath("~/Pictures", "bench_dot3_v3.png"), plot(dot3_bench));
save(joinpath("~/Pictures", "bench_sse_v3.png"), plot(sse_bench));
save(joinpath("~/Pictures", "bench_exp_v3.png"), plot(exp_bench));
save(joinpath("~/Pictures", "bench_aplusBc_v3.png"), plot(aplusBc_bench));