text:
module HGLET

include("utils.jl")

using ..GraphSignal, ..GraphPartition, ..BasisSpecification, ..GHWT, SparseArrays

include("common.jl")

export HGLET_Synthesis, HGLET_jkl, HGLET_Analysis, HGLET_Analysis_All,
HGLET_BestBasis, HGLET_GHWT_BestBasis, HGLET_GHWT_Synthesis, HGLET_GHWT_BestBasis_minrelerror

"""
function HGLET_Synthesis(dvec::Vector{Float64}, GP::GraphPart, BS::BasisSpec,
                         G::GraphSig; gltype::Symbol = :L)

### Input Arguments
 * `dvec`: the expansion coefficients corresponding to the chosen basis
 * `GP`: a GraphPart object
 * `BS`: a BasisSpec object
 * `G`: a GraphSig object
 * `gltype`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used

### Output Argument
 * `f`: the reconstructed signal
 * `GS`: the reconstructed GraphSig object
"""
function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
                         G::GraphSig; gltype::Symbol = :L)
    # Preliminaries
    W = G.W
    jmax = size(GP.rs,2)

    # Fill in the appropriate entries of dmatrix
    dmatrix = dvec2dmatrix(dvec,GP,BS)
    f = dmatrix[:,jmax,:]

    # Perform the signal synthesis from the given coefficients
    for j = jmax:-1:1
        regioncount = count(!iszero, GP.rs[:,j]) - 1
        for r = 1:regioncount
            # the index that marks the start of the region
            rs1 = GP.rs[r,j]

            # the index that is one after the end of the region
            rs3 = GP.rs[r+1,j]

            # the number of points in the current region
            n = rs3 - rs1

            # only proceed forward if coefficients do not exist
            if (j == jmax || count(!iszero, dmatrix[rs1:rs3-1,j+1,:]) == 0) && count(!iszero, dmatrix[rs1:rs3-1,j,:]) > 0

                if n == 1
                    f[rs1,:] = dmatrix[rs1,j,:]
                elseif n > 1
                    indrs = GP.ind[rs1:rs3-1]
                    W_temp = W[indrs,indrs]
                    D = Diagonal(vec(sum(W_temp, dims = 1)))
                    D2 =Diagonal(vec(sum(W_temp, dims = 1)).^(-0.5))
                    normalizep = false # a flag for normalizing L for Lsym/Lrw
                    if minimum(diag(D)) > 10^3*eps() # avoiding 1/small entries
                        normalizep = true
                    end
                    # compute the GL eigenvectors via svd.
                    if ( gltype == :Lrw || gltype == :Lsym ) && normalizep
                        v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
                    else
                        v,_,_ = svd(D - Matrix(W_temp))
                        if gltype != :L
                            @warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(gltype))
                        end
                    end
                    v = v[:,end:-1:1] # reorder the ev's in the decreasing ew's

                    # standardize the eigenvector signs
                    for col = 1:n
                        row = 1
                        standardized = false
                        while !standardized
                            if v[row,col] > 10^3*eps()
                                standardized = true
                            elseif v[row,col] < -10^3*eps()
                                v[:,col] = - v[:,col]
                                standardized = true
                            else
                                row += 1
                            end
                        end
                    end

                    # reconstruct the signal
                    if gltype == :Lrw && normalizep
                        f[rs1:rs3-1,:] = D2*v*dmatrix[rs1:rs3-1,j,:]
                    else
                        f[rs1:rs3-1,:] = v*dmatrix[rs1:rs3-1,j,:]
                    end
                end
            end
        end
    end

    # put the reconstructed values in the correct order
    f[GP.ind,:] = f

    # creat a GraphSig object with the reconstructed data
    GS = deepcopy(G)
    replace_data!(GS,f)

    return f, GS
end

"""
function HGLET_jkl(GP::GraphPart, drow::Int, dcol::Int)

    Generate the (j,k,l) indices for the HGLET basis vector corresponding to
        the coefficient dmatrix(drow,dcol)

### Input Arguments
 * `GP`:    a GraphPart object
 * `drow`:  the row of the expansion coefficient
 * `dcol`:  the column of the expansion coefficient

### Output Argument
 * `j`: the level index of the expansion coefficient
 * `k`: the subregion index of the expansion coefficient
 * `l`: the eigenvector index of the expansion coefficient
"""
function HGLET_jkl(GP::GraphPart, drow::Int, dcol::Int)
    j = dcol - 1

    k = findfirst(GP.rs[:,dcol] .> drow)
    k = k-2

    l = drow - GP.rs[k+1,dcol]

    return j,k,l
end

"""
function HGLET_Analysis(G::GraphSig, GP::GraphPart, gltype::Symbol = :L)

    For a GraphSig object 'G', generate the matrix of HGLET expansion
    coefficients corresponding to the eigenvectors of specified version of
    the graph Laplacian matrix, i.e., either L, Lrw, or Lsym

### Input Arguments
 * `G`:  a GraphSig object
 * `GP`: a GraphPart object
 * `gltype`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used

### Output Argument
 * `dmatrix`: the matrix of expansion coefficients using the specific GL matrix
"""
function HGLET_Analysis(G::GraphSig, GP::GraphPart, gltype::Symbol = :L)
    # Preliminaries
    W = G.W
    ind = GP.ind
    rs = GP.rs
    N = size(G.W,1)
    jmax = size(rs,2)
    fcols = size(G.f,2)
    dmatrix = zeros(N,jmax,fcols)
    dmatrix[:,jmax,:] = G.f[ind,:]

    # Perform the HGLET analysis, i.e., generating the HGLET coefficients
    for j = jmax-1:-1:1
        regioncount = count(!iszero, rs[:,j]) - 1
        for r = 1:regioncount
            # the index that marks the start of the region
            rs1 = rs[r,j]

            # the index that is one after the end of the region
            rs3 = rs[r+1,j]

            # the number of points in the current region
            n = rs3 - rs1

            if n == 1
                dmatrix[rs1,j,:] = G.f[ind[rs1],:]
            elseif n > 1
                indrs = ind[rs1:rs3-1]
                W_temp = W[indrs,indrs]
                D = Diagonal(vec(sum(W_temp, dims = 1)))
                D2 = Diagonal(vec(sum(W_temp, dims = 1)).^(-0.5))
                normalizep = false # a flag for normalizing L for Lsym and Lrw
                if minimum(diag(D)) > 10^3*eps() # avoiding 1/small entries
                    normalizep = true
                end

                # compute the GL eigenvectors via svd
                if ( gltype == :Lrw || gltype == :Lsym ) && normalizep
                    v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
                else
                    v,_,_ = svd(D - Matrix(W_temp))
                    if gltype != :L
                        @warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(gltype))
                    end
                end
                v = v[:,end:-1:1] # reorder the ev's in the decreasing ew's

                # standardize the eigenvector signs
                for col = 1:n
                    row = 1
                    standardized = false
                    while !standardized
                        if v[row,col] > 10^3*eps()
                            standardized = true
                        elseif v[row,col] < -10^3*eps()
                            v[:,col] = - v[:,col]
                            standardized = true
                        else
                            row += 1
                        end
                    end
                end

                # obtain the expansion coeffcients
                if gltype == :Lrw && normalizep
                    dmatrix[rs1:rs3-1,j,:] = v'*(D.^0.5)*G.f[indrs,:]
                else
                    dmatrix[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
                end
            end
        end
    end
    return dmatrix
end

"""
function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)

    For a GraphSig object 'G', generate the 3 matrices of HGLET expansion
    coefficients corresponding to the eigenvectors of L, Lrw and Lsym

### Input Arguments
 * `G`:  a GraphSig object
 * `GP`: a GraphPart object

### Output Argument
 * `dmatrixH`:        the matrix of expansion coefficients for L
 * `dmatrixHrw`:      the matrix of expansion coefficients for Lrw
 * `dmatrixHsym`:     the matrix of expansion coefficients for Lsym
"""
function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
# Preliminaries
    W = G.W
    ind = GP.ind
    rs = GP.rs
    # This is the type of graph Laplacian matrix used for graph partitioning;
    # Not the eigenvector dictionary computed for the HGLET.
    method = GP.method
    N = size(G.W,1)
    jmax = size(rs,2)
    fcols = size(G.f,2)
    dmatrixH = zeros(N,jmax,fcols)
    dmatrixH[:,jmax,:] = G.f[ind,:]
    dmatrixHrw = deepcopy(dmatrixH)
    dmatrixHsym = deepcopy(dmatrixH)

# Perform the transform ==> eigenvectors of L
    for j = jmax-1:-1:1
        regioncount = count(!iszero, rs[:,j]) - 1
        for r = 1:regioncount
            # the index that marks the start of the region
            rs1 = rs[r,j]

            # the index that is one after the end of the region
            rs3 = rs[r+1,j]

            # the number of points in the current region
            n = rs3-rs1

            if n == 1
                dmatrixH[rs1,j,:] = G.f[ind[rs1],:]

            elseif n > 1
                indrs = ind[rs1:rs3-1]

                # compute the eigenvectors of L ==> svd(L)
                v,_,_ = svd(Diagonal(vec(sum(W[indrs,indrs],dims = 1)))
                        - Matrix(W[indrs,indrs]))
                v = v[:,end:-1:1]

                # standardize the eigenvector signs
                for col = 1:n
                    row = 1
                    standardized = false
                    while !standardized
                        if v[row,col] > 10^3*eps()
                            standardized = true
                        elseif v[row,col] < -10^3*eps()
                            v[:,col] = - v[:,col]
                            standardized = true
                        else
                            row += 1
                        end
                    end
                end

                # obtain the expansion coefficients
                dmatrixH[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
            end
        end
    end

    for j = jmax-1:-1:1
        regioncount = count(!iszero, rs[:,j])-1
        for r = 1:regioncount
            # the index that marks the start of the region
            rs1 = rs[r,j]

            # the index that is one after the end of the region

            rs3 = rs[r+1,j]

            # the number of points in the current region
            n = rs3 - rs1

            if n == 1
                dmatrixHrw[rs1,j,:] = G.f[ind[rs1],:]
                dmatrixHsym[rs1,j,:] = G.f[ind[rs1],:]

            elseif n > 1
                indrs = ind[rs1:rs3-1]

                # compute the eigenvectors
                if minimum(sum(W[indrs,indrs],dims = 1)) > 10^3*eps()
                    useLrw = true

                    ### eigenvectors of L_rw ==> svd(L_sym)
                    W_temp = W[indrs,indrs]
                    D = Diagonal(vec(sum(W_temp,dims = 1)))
                    D2 =Diagonal(vec(sum(W_temp,dims = 1)).^(-0.5))
                    v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
                    v = v[:,end:-1:1]
                else
                    useLrw = false

                    ### eigenvectors of L ==> svd(L)
                    W_temp = W[indrs,indrs]
                    D = Diagonal(vec(sum(W_temp,dims = 1)))
                    v,_,_ = svd(D - Matrix(W_temp))
                    v = v[:,end:-1:1]
                end

                #standardized the eigenvector signs
                for col = 1:n
                    row = 1
                    standardized = false
                    while !standardized
                        if v[row,col] > 10^3*eps()
                            standardized = true
                        elseif v[row,col] < -10^3*eps()
                            v[:,col] = - v[:,col]
                            standardized = true
                        else
                            row += 1
                        end
                    end
                end

                # obtain the expansion coefficients for L_sym
                dmatrixHsym[rs1:rs3-1,j,:] = v'*G.f[indrs,:]

                # obtain the expansion coeffcients for L_rw
                if useLrw
                    dmatrixHrw[rs1:rs3-1,j,:] = v'*(Diagonal(vec(sum(W_temp,dims = 1)).^0.5))*G.f[indrs,:]
                else
                    dmatrixHrw[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
                end

            end
        end
    end

    return dmatrixH, dmatrixHrw, dmatrixHsym

end

"""
function HGLET_BestBasis(GP::GraphPart;
    dmatrix::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    gltype::Symbol = :L, cfspec::Any = 0.1, flatten::Any = 1)

    Select the HGLET best basis from the input matrix of expansion coefficients

### Input Arguments
* `GP`:      a GraphPart object
* `dmatrix`: the matrix of HGLET expansion coefficients
* `gltype`:  the type of graph Laplacian matrix used for HGLET dictionary
    (default = :L)
* `cfspec`: the cost functional specification to be used: see utils.jl
    for the detail. If it's a number, say, p, then the l^p norm is used.
    If it's a function, then that function is used. Default is 0.1, i.e., l^0.1
* `flatten`: the method for flattening vector-valued data to scalar-valued data;
    If it's a number, say, p, then the sum of the l^p norm is computed.
    There are many other options, such as :ash, :entropy, etc. See utils.jl
    for the detail. Default is 1, i.e, the sum of the l^1 norm of the coefs.

### Output Argument
* `dvec`: the vector of expansion coefficients corresponding to the bestbasis
* `BS`:   a BasisSpec object which specifies the best basis
"""
function HGLET_BestBasis(GP::GraphPart;
    dmatrix::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    gltype::Symbol = :L, cfspec::Any = 0.1, flatten::Any = 1)

    # the cost functional to be used
    costfun = cost_functional(cfspec)

    # constants and dmatrix cleanup
    N, jmax, fcols = size(dmatrix)
    dmatrix[abs.(dmatrix) .< 10^2*eps()] .= 0

    # flatten dmatrix for more than one input signals
    if fcols > 1
        dmatrix0 = deepcopy(dmatrix)
        dmatrix = dropdims(dmatrix_flatten(dmatrix, flatten), dims = 3)
    end

    #
    # Find the HGLET best basis now!
    #

    # allocate/initialize ==> order matters here
    dvec = dmatrix[:,jmax]
    levlist = jmax*ones(Int,N)

    # set the tolerance
    tol = 10^4*eps()

    # perform the basis search
    for j = jmax:-1:1
        regioncount = count(!iszero, GP.rs[:,j]) - 1
        for r = 1:regioncount
            indr = GP.rs[r,j]:(GP.rs[r+1,j]-1)
            ### compute the cost of the current best basis
            costBB = costfun(dvec[indr])

            ### compute the cost of the HGLET coefficients
            costNew = costfun(dmatrix[indr,j])
            # change the best basis if the new cost is less expensive
            if costBB >= costNew - tol
                costBB, dvec[indr], levlist[indr] = BBchange(costNew, dmatrix[indr,j], j)
            end
        end
    end

    levlist = collect(enumerate(levlist))

    BS = BasisSpec(levlist, c2f = true, description = "HGLET Best Basis")
    #levlist2levlengths!(GP, BS)

    # if we flattened dmatrix, then "unflatten" the expansion coefficients
    if fcols > 1
        dvec = dmatrix2dvec(dmatrix0, GP, BS)
    end

    return dvec, BS
end

"""
function HGLET_GHWT_BestBasis(GP::GraphPart;
    dmatrixH::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    cfspec::Any = 0.1,flatten::Any = 1)

    Select the best basis from several matrices of expansion coefficients

### Input Arguments
* `dmatrixH`:    the matrix of HGLET expansion coefficients ==> eigenvectors of L
* `dmatrixHrw`:  the matrix of HGLET expansion coefficients ==> eigenvectors of Lrw
* `dmatrixHsym`: the matrix of HGLET expansion coefficients ==> eigenvectors of Lsym
* `dmatrixG`:    the matrix of GHWT expansion coefficients
* `GP`:          a GraphPart object
* `cfspec`: the cost functional specification to be used: see utils.jl
    for the detail. If it's a number, say, p, then the l^p norm is used.
    If it's a function, then that function is used. Default is 0.1, i.e., l^0.1
* `flatten`: the method for flattening vector-valued data to scalar-valued data;
    If it's a number, say, p, then the sum of the l^p norm is computed.
    There are many other options, such as :ash, :entropy, etc. See utils.jl
    for the detail. Default is 1, i.e, the sum of the l^1 norm of the coefs.

### Output Argument
* `dvec`:     the vector of expansion coefficients corresponding to the bestbasis
* `BS`:       a BasisSpec object which specifies the best basis
* `trans`:    specifies which transform was used for that portion of the signal:
    00 = HGLET with L
    01 = HGLET with Lrw
    10 = HGLET with Lsym
    11 = GHWT
"""
function HGLET_GHWT_BestBasis(GP::GraphPart;
    dmatrixH::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    cfspec::Any = 0.1,flatten::Any = 1)
    # specify transform codes
    transHsym = [true false]
    transG = [true true]
    transHrw = [false true]
    transH = [false false]

    # the cost functional to be used
    costfun = cost_functional(cfspec)

    # constants and dmatrix cleanup
    if !isempty(dmatrixHsym)
        N, jmax, fcols = size(dmatrixHsym)
        dmatrixHsym[abs.(dmatrixHsym).<10^2*eps()] .= 0
    end
    if !isempty(dmatrixG)
        N, jmax, fcols = size(dmatrixG)
        dmatrixG[abs.(dmatrixG).<10^2*eps()] .= 0
    end
    if !isempty(dmatrixHrw)
        N, jmax, fcols = size(dmatrixHrw)
        dmatrixHrw[abs.(dmatrixHrw).<10^2*eps()] .= 0
    end
    if !isempty(dmatrixH)
        N, jmax, fcols = size(dmatrixH)
        dmatrixH[abs.(dmatrixH).<10^2*eps()] .= 0
    end

    # flatten dmatrix
    if fcols > 1
        if !isempty(dmatrixHsym)
            dmatrix0Hsym = deepcopy(dmatrixHsym)
            dmatrixHsym = dropdims(dmatrix_flatten(dmatrixHsym,flatten),dims = 3)
        end
        if !isempty(dmatrixG)
            dmatrix0G = deepcopy(dmatrixG)
            dmatrixG = dropdims(dmatrix_flatten(dmatrixG,flatten),dims = 3)
        end
        if !isempty(dmatrixHrw)
            dmatrix0Hrw = deepcopy(dmatrixHrw)
            dmatrixHrw = dropdims(dmatrix_flatten(dmatrixHrw,flatten),dims = 3)
        end
        if !isempty(dmatrixH)
            dmatrix0H = deepcopy(dmatrixH)
            dmatrixH = dropdims(dmatrix_flatten(dmatrixH,flatten),dims = 3)
        end
    end

    # Find the HGLET/GHWT best basis

    # allocate/initialize ==> order matters here
    if !isempty(dmatrixHsym)
        dvec = dmatrixHsym[:,jmax]
        trans = repeat(transHsym,N,1)
    end
    if !isempty(dmatrixG)
        dvec = dmatrixG[:,jmax]
        trans = repeat(transG,N,1)
    end
    if !isempty(dmatrixHrw)
        dvec = dmatrixHrw[:,jmax]
        trans = repeat(transHrw,N,1)
    end
    if !isempty(dmatrixH)
        dvec = dmatrixH[:,jmax]
        trans = repeat(transH,N,1)
    end
    levlist = jmax*ones(Int,N)

    # set the tolerance
    tol = 10^4*eps()

    # perform the basis search
    for j = jmax:-1:1
        regioncount = count(!iszero, GP.rs[:,j]) - 1
        for r = 1:regioncount
            indr = GP.rs[r,j]:(GP.rs[r+1,j]-1)
            ### compute the cost of the current best basis
            costBB = costfun(dvec[indr])

            ### compute the cost of the HGLET-Lsym coefficients
            if !isempty(dmatrixHsym)
                costNew = costfun(dmatrixHsym[indr,j])
                # change the best basis if the new cost is less expensive
                if costBB >= costNew - tol
                    costBB, dvec[indr], levlist[indr], trans[indr,:] = BBchange(costNew, dmatrixHsym[indr,j],j,transHsym)
                end
            end

            ### compute the cost of the GHWT coefficients
            if !isempty(dmatrixG)
                costNew = costfun(dmatrixG[indr,j])
                # change the best basis if the new cost is less expensive
                if costBB >= costNew - tol
                    costBB, dvec[indr], levlist[indr], trans[indr,:] = BBchange(costNew, dmatrixG[indr,j],j,transG)
                end
            end

            ### compute the cost of the HGLET-Lrw coefficients
            if !isempty(dmatrixHrw)
                costNew = costfun(dmatrixHrw[indr,j])
                # change the best basis if the new cost is less expensive
                if costBB >= costNew - tol
                    costBB, dvec[indr], levlist[indr], trans[indr,:] = BBchange(costNew, dmatrixHrw[indr,j],j,transHrw)
                end
            end

            ### compute the cost of the HGLET-L coefficients
            if !isempty(dmatrixH)
                costNew = costfun(dmatrixH[indr,j])
                # change the best basis if the new cost is less expensive
                if costBB >= costNew - tol
                    _, dvec[indr],levlist[indr],trans[indr,:] = BBchange(costNew, dmatrixH[indr,j],j,transH)
                end
            end
        end
    end

    transfull = deepcopy(trans)
    #trans = trans[levlist.!=0,:]
    levlist = collect(enumerate(levlist))

    BS = BasisSpec(levlist,c2f = true, description = "HGLET-GHWT Best Basis")
    #levlist2levlengths!(GP,BS)

    # if we flattened dmatrix, then "unflatten" the expansion coefficients
    if fcols > 1
        # create vectors of coefficients (which are zero if the transform's coefficients were not included as function inputs)
        if !isempty(dmatrixH)
            dvecH = dmatrix2dvec(dmatrix0H,GP,BS)
        else
            dvecH = zeros(N,fcols)
        end
        if !isempty(dmatrixHrw)
            dvecHrw = dmatrix2dvec(dmatrix0Hrw,GP,BS)
        else
            dvecHrw = zeros(N,fcols)
        end
        if !isempty(dmatrixHsym)
            dvecHsym = dmatrix2dvec(dmatrix0Hsym,GP,BS)
        else
            dvecHsym = zeros(N,fcols)
        end
        if !isempty(dmatrixG)
            dvecG = dmatrix2dvec(dmatrix0G,GP,BS)
        else
            dvecG = zeros(N,fcols)
        end

        dvec = dvecHsym.*(transfull[:,1].*(.~transfull[:,2]))
        + dvecG.*(transfull[:,1].*transfull[:,2])
        + dvecHrw.*((.~transfull[:,1]).*transfull[:,2])
        + dvecH.*((.~transfull[:,1]).*(.~transfull[:,2]))
    end

return dvec, BS, transfull

end

function BBchange(costNew::Float64, dvec::Vector{Float64}, j::Int)
    #change to the new best basis
    costBB = costNew

    n = length(dvec)

    levlist = fill(j, n)

    return costBB, dvec, levlist
end

function BBchange(costNew::Float64, dvec::Vector{Float64}, j::Int, trans::Array{Bool,2})
    #change to the new best basis
    costBB = costNew

    n = length(dvec)

    levlist = fill(j, n)

    trans = repeat(trans,n,1)

    return costBB, dvec, levlist, trans
end

"""
function HGLET_GHWT_Synthesis(dvec::Matrix{Float64},GP::GraphPart,BS::BasisSpec,trans::Array{Bool,2},G::GraphSig)

    Given a vector of HGLET & GHWT expansion coefficients, info about the
    graph partitioning, and the choice of basis and corresponding transforms,
    reconstruct the signal.

### Input Arguments
* `dvec`:   the expansion coefficients corresponding to the chosen basis
* `GP`:     a GraphPart object
* `BS`:     a BasisSpec object
* `trans`:  a specification of the transforms used for the HGLET-GHWT hybrid transform
    00 = HGLET with L
    01 = HGLET with Lrw
    10 = HGLET with Lsym
    11 = GHWT
* `G`:      a GraphSig object

### Output Argument
* `f`:      the reconstructed signal
* `GS`:     the reconstructed GraphSig object
"""
function HGLET_GHWT_Synthesis(dvec::Matrix{Float64},GP::GraphPart,BS::BasisSpec,trans::Array{Bool,2},G::GraphSig)
    # fill out trans
    #_,_,transfull = BSfull(GP,BS,trans)
    transfull = trans
    # decompose dvec into GHWT and HGLET components
    dvecHsym= dvec.*(transfull[:,1].*(.~transfull[:,2]))
    dvecG   = dvec.*(transfull[:,1].*transfull[:,2])
    dvecHrw = dvec.*((.~transfull[:,1]).*transfull[:,2])
    dvecH   = dvec.*((.~transfull[:,1]).*(.~transfull[:,2]))

    # Synthesize using the transforms separately

    fH,_ = HGLET_Synthesis(dvecH, GP, BS, G)
    fHrw,_ = HGLET_Synthesis(dvecHrw, GP, BS, G, gltype = :Lrw)
    fHsym,_ = HGLET_Synthesis(dvecHsym, GP, BS, G, gltype = :Lsym)
    fG = ghwt_synthesis(dvecG, GP, BS)

    f = fH + fHrw + fHsym + fG

    GS = deepcopy(G)
    replace_data!(GS,f)

    return f, GS

end

"""
function HGLET_GHWT_BestBasis_minrelerror(GP::GraphPart, G::GraphSig;
    dmatrixH::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    compare::Bool = true)

    Find the best basis for approximating the signal 'G' by performing the
        best basis search with a range of tau-measures as cost functionals
        (tau = 0.1,0.2,...,1.9) and minimizing the relative error.

### Input argument:
* `dmatrixH`: the matrix of HGLET expansion coefficients ==> eigenvectors of L
* `dmatrixHrw`: the matrix of HGLET expansion coefficients ==> eigenvectors of Lrw
* `dmatrixHsym`: the matrix of HGLET expansion coefficients ==> eigenvectors of Lsym
* `dmatrixG`: the matrix of GHWT expansion coefficients
* `GP`: a GraphPart object
* `G`: the GraphSig object
* `compare`: if it is false, don't compare the hybrid best basis to the GHWT fine-to-coarse best basis

### Output argument:
* `dvec`: the vector of expansion coefficients corresponding to the bestbasis
* `BS`: a BasisSpec object which specifies the best basis
* `trans`: specifies which transform was used for that portion of the signal
    00 = HGLET with L
    01 = HGLET with Lrw
    10 = HGLET with Lsym
    11 = GHWT
* `tau`: the tau that yields the smallest relative error

"""
function HGLET_GHWT_BestBasis_minrelerror(GP::GraphPart, G::GraphSig;
    dmatrixH::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
    compare::Bool = true)

    dvec, BS, trans, tau = 0, 0, 0, 0 # predefine
    sumrelerror = Inf
    for tau_temp = 0.1:0.1:1.9
        # we are only considering the GHWT
        if isempty(dmatrixH) && isempty(dmatrixHrw) && isempty(dmatrixHsym) && !isempty(dmatrixG)
            dvec_temp, BS_temp = ghwt_bestbasis(dmatrixG, GP, cfspec = tau_temp)
            trans_temp = trues(length(BS_temp.levlist),2)
            orthbasis = true
        else
            dvec_temp, BS_temp, trans_temp = HGLET_GHWT_BestBasis(GP,
                dmatrixH = dmatrixH, dmatrixG = dmatrixG,
                dmatrixHrw = dmatrixHrw, dmatrixHsym = dmatrixHsym,
                cfspec = tau_temp)

            # check whether any HGLET Lrw basis vectors are in the best basis
            orthbasis = true
            rows = size(trans_temp,1)
            for row = 1:rows
                if !trans_temp[row,1] && trans_temp[row,2]
                    orthbasis = false
                    break
                end
            end
        end

        # compute the relative errors
        if orthbasis
            relerror_temp = orth2relerror(dvec_temp[:])
        else
            B,_ = HGLET_GHWT_Synthesis(Matrix{Float64}(I,length(dvec_temp),length(dvec_temp)),
                GP, BS_temp, trans_temp, G)
            relerror_temp = nonorth2relerror(dvec_temp[:],B)
        end
        sumrelerror_temp = sum(relerror_temp)

        # consider the GHWT fine-to-coarse best basis
        if compare == true && !isempty(dmatrixG)
            dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrixG, GP, cfspec = tau_temp)
            sumrelerror_f2c = sum(orth2relerror(dvec_f2c[:]))
            if sumrelerror_f2c < sumrelerror_temp
                sumrelerror_temp = sumrelerror_f2c
                dvec_temp = copy(dvec_f2c)
                BS_temp = deepcopy(BS_f2c)
                trans_temp = repeat([true true], length(BS_f2c.levlist),1)
            end
        end

        # compare to the current lowest sum of relative errors
        if  sumrelerror_temp < sumrelerror
            dvec = copy(dvec_temp)
            BS = deepcopy(BS_temp)
            trans = copy(trans_temp)
            sumrelerror = sumrelerror_temp
            tau = tau_temp
        end
    end
    return dvec, BS, trans, tau
end

end # end of module HGLET


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames
using Random

# 画面レイアウト
win = Window("Bingo Machine")
vbox = Box(:v)
push!(win, vbox)

## メニューバー
menu_bar = MenuBar()
file = MenuItem("File")
filemenu = Menu(file)
menu_history = MenuItem("History")
menu_quit = MenuItem("Quit")
push!(vbox, menu_bar)
push!(menu_bar, file)
push!(filemenu, menu_history)
push!(filemenu, menu_quit)

## 数字表示部分
number_display = Label("")
set_gtk_property!(number_display, :use_markup, true)
hbox = Box(:h)
push!(vbox, number_display)
push!(vbox, hbox)

## 数字を更新、リセットするためのボタン
next_button = Button("Next")
reset_button = Button("Reset")
push!(hbox, next_button)
push!(hbox, reset_button)

set_gtk_property!(vbox, :expand, number_display, true)
set_gtk_property!(hbox, :expand, next_button, true)
set_gtk_property!(hbox, :expand, reset_button, true)

maxnum = 75
number_list = shuffle(1:maxnum)
number_history = []

function next_number()
    if length(number_list) != 0
        n = pop!(number_list)
        push!(number_history, n)
        GAccessor.markup(number_display, "$n")
    else
        GAccessor.markup(number_display, "Finish!")
    end

    return nothing
end

function reset_number()
    global number_list = shuffle(1:maxnum)
    global number_history = []
    set_gtk_property!(number_display, :label, "")

    return nothing
end

function show_history()
    win = Window("Number history")
    history = Label("" * join(number_history, " ") * "")
    set_gtk_property!(history, :use_markup, true)
    push!(win, history)
    GAccessor.line_wrap(history, true)
    showall(win)

    return nothing
end

showall(win)

signal_connect(x -> next_number(), next_button, "clicked")
signal_connect(x -> reset_number(), reset_button, "clicked")
signal_connect((x,y) -> show_history(), menu_history, :activate, Nothing, (), false)

signal_connect(menu_quit, :activate) do w
  if !isinteractive()
    Gtk.gtk_quit()
  else
    exit()
  end
end

if !isinteractive()
    c = Condition()
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    wait(c)
end

ext: jl

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
# Tests for CartesianGrid
using DiffEqJump, LightGraphs
using Test, Random

io = IOBuffer()
rng = MersenneTwister()
dims = (4,3,2)
sites = rand(1:prod(dims), 10)
num_samples = 10^5
rel_tol = 0.01
grids = [DiffEqJump.CartesianGridRej(dims), DiffEqJump.CartesianGridIter(dims), LightGraphs.grid(dims)]
for grid in grids
    show(io, "text/plain", grid)
    @test String(take!(io)) !== nothing
    @test DiffEqJump.num_sites(grid) == prod(dims)
    @test DiffEqJump.outdegree(grid, 1) == 3
    @test DiffEqJump.outdegree(grid, 4) == 3
    @test DiffEqJump.outdegree(grid, 17) == 4
    @test DiffEqJump.outdegree(grid, 21) == 3
    @test DiffEqJump.outdegree(grid, 6) == 5
    for site in sites
        d = Dict{Int,Int}()
        for i in 1:num_samples
            nb = DiffEqJump.rand_nbr(rng, grid, site)
            nb in keys(d) ? d[nb] += 1 : d[nb] = 1
        end
        for val in values(d)
            @test abs(val/num_samples - 1/DiffEqJump.outdegree(grid,site)) < rel_tol
        end
    end
end

# setup
rel_tol = 0.05
num_samples = 10^4
dims = (3,3,3)
g = CartesianGridRej(dims)
num_nodes = DiffEqJump.num_sites(g)
num_species = 3
reactstoch = [[1 => 1, 2 => 1],[3 => 1]]
netstoch = [[1 => -1, 2 => -1, 3 => 1],[1 => 1, 2 => 1, 3 => -1]]
rates = [0.1, 1.]
num_rxs = length(rates)
ma_jumps = MassActionJump(rates, reactstoch, netstoch)
u = ones(Int, num_species, num_nodes)
rng = MersenneTwister()

# Tests for RxRates
rx_rates = DiffEqJump.RxRates(num_nodes, ma_jumps)
show(io, "text/plain", rx_rates)
for site in 1:num_nodes
    DiffEqJump.update_rx_rates!(rx_rates, 1:num_rxs, u, site)
    rx_props = [DiffEqJump.evalrxrate(u[:, site], rx, ma_jumps) for rx in 1:num_rxs]
    rx_probs = rx_props/sum(rx_props)
    d = Dict{Int,Int}()
    for i in 1:num_samples
        rx = DiffEqJump.sample_rx_at_site(rx_rates, site, rng)
        rx in keys(d) ? d[rx] += 1 : d[rx] = 1
    end
    for (k,v) in d
        @test abs(v/num_samples - rx_probs[k]) < rel_tol
    end
end

# Tests for HopRatesGraphDsi
hopping_constants = ones(num_species, num_nodes)
hop_rates = DiffEqJump.HopRatesGraphDsi(hopping_constants)
show(io, "text/plain", hop_rates)
spec_probs = ones(num_species)/num_species

for site in 1:num_nodes
    DiffEqJump.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
    num_nbs = DiffEqJump.outdegree(g, site)
    target_probs = ones(num_nbs)/num_nbs
    d1 = Dict{Int,Int}()
    d2 = Dict{Int,Int}()
    for i in 1:num_samples
        spec, target = DiffEqJump.sample_hop_at_site(hop_rates, site, rng, g)
        d1[spec] = get(d1, spec, 0) + 1
        d2[target] = get(d2, target, 0) + 1
    end
    @test maximum(abs.(collect(values(d1))/num_samples - spec_probs)) < rel_tol
    @test maximum(abs.(collect(values(d2))/num_samples - target_probs)) < rel_tol
end

# Tests for HopRatesGraphDsij
hop_constants = Matrix{Vector{Float64}}(undef, num_species, num_nodes)
for ci in CartesianIndices(hop_constants)
    (species, site) = Tuple(ci)
    hop_constants[ci] = repeat([1.0], DiffEqJump.outdegree(g, site))
end
spec_probs = ones(num_species)/num_species
hop_rates_structs = [DiffEqJump.HopRatesGraphDsij(hop_constants), DiffEqJump.HopRatesGridDsij(hop_constants, g)]

for hop_rates in hop_rates_structs
    show(io, "text/plain", hop_rates)
    for site in 1:num_nodes
        DiffEqJump.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
        num_nbs = DiffEqJump.outdegree(g, site)
        target_probs = ones(num_nbs)/num_nbs
        d1 = Dict{Int,Int}()
        d2 = Dict{Int,Int}()
        for i in 1:num_samples
            spec, target = DiffEqJump.sample_hop_at_site(hop_rates, site, rng, g)
            d1[spec] = get(d1, spec, 0) + 1
            d2[target] = get(d2, target, 0) + 1
        end
        @test maximum(abs.(collect(values(d1))/num_samples - spec_probs)) < rel_tol
        @test maximum(abs.(collect(values(d2))/num_samples - target_probs)) < rel_tol
    end
end

# Tests for HopRatesGraphDsLij
species_hop_constants = ones(num_species)
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
for site in 1:num_nodes
    site_hop_constants[site] = repeat([1.0], DiffEqJump.outdegree(g, site))
end
spec_probs = ones(num_species)/num_species
hop_rates_structs = [DiffEqJump.HopRatesGraphDsLij(species_hop_constants, site_hop_constants), DiffEqJump.HopRatesGridDsLij(species_hop_constants, site_hop_constants, g)]

for hop_rates in hop_rates_structs
    show(io, "text/plain", hop_rates)
    for site in 1:num_nodes
        DiffEqJump.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
        num_nbs = DiffEqJump.outdegree(g, site)
        target_probs = ones(num_nbs)/num_nbs
        d1 = Dict{Int,Int}()
        d2 = Dict{Int,Int}()
        for i in 1:num_samples
            spec, target = DiffEqJump.sample_hop_at_site(hop_rates, site, rng, g)
            d1[spec] = get(d1, spec, 0) + 1
            d2[target] = get(d2, target, 0) + 1
        end
        @test maximum(abs.(collect(values(d1))/num_samples - spec_probs)) < rel_tol
        @test maximum(abs.(collect(values(d2))/num_samples - target_probs)) < rel_tol
    end
end

# Tests for HopRatesGraphDsiLij
species_hop_constants = ones(num_species, num_nodes)
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
for site in 1:num_nodes
    site_hop_constants[site] = repeat([1.0], DiffEqJump.outdegree(g, site))
end
spec_probs = ones(num_species)/num_species
hop_rates_structs = [DiffEqJump.HopRatesGraphDsiLij(species_hop_constants, site_hop_constants), DiffEqJump.HopRatesGridDsiLij(species_hop_constants, site_hop_constants, g)]

for hop_rates in hop_rates_structs
    show(io, "text/plain", hop_rates)
    for site in 1:num_nodes
        DiffEqJump.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
        num_nbs = DiffEqJump.outdegree(g, site)
        target_probs = ones(num_nbs)/num_nbs
        d1 = Dict{Int,Int}()
        d2 = Dict{Int,Int}()
        for i in 1:num_samples
            spec, target = DiffEqJump.sample_hop_at_site(hop_rates, site, rng, g)
            d1[spec] = get(d1, spec, 0) + 1
            d2[target] = get(d2, target, 0) + 1
        end
        @test maximum(abs.(collect(values(d1))/num_samples - spec_probs)) < rel_tol
        @test maximum(abs.(collect(values(d2))/num_samples - target_probs)) < rel_tol
    end
end

@test String(take!(io)) !== nothing


text:
println("start run")

using JuMP, EAGO

m = Model(EAGO.Optimizer)
set_optimizer_attribute(m, "mul_relax_style", 1)
set_optimizer_attribute(m, "verbosity", 1)
set_optimizer_attribute(m, "output_iterations", 1000)
set_optimizer_attribute(m, "iteration_limit", 10000000)
set_optimizer_attribute(m, "cut_max_iterations", 2)

# OBBT depth 0 -> 20... increases number of iterations...
set_optimizer_attribute(m, "obbt_depth", 8)
set_optimizer_attribute(m, "obbt_repetitions", 2)

# ----- Variables ----- #
x_Idx = Any[2, 3, 4]
@variable(m, x[x_Idx])
set_lower_bound(x[2], 1.0e-6)
set_upper_bound(x[2], 1.0)
set_lower_bound(x[3], 1.0e-6)
set_upper_bound(x[3], 1.0)
set_lower_bound(x[4], 1.0e-6)
set_upper_bound(x[4], 1.0)

# ----- Constraints ----- #
s = time()
#=
@NLobjective(m, Min, (15.3261663216011*x[2]+23.2043471859416*x[3]+6.69678129464404*x[4])*log(2.1055*x[2]+3.1878*x[3] + 0.92*x[4]) +
                      1.04055250396734*x[2]-2.24199441248417*x[3]+3.1618173099828*x[4]+6.4661663216011*x[2]*log(x[2]/(2.1055*x[2]+3.1878*x[3]+0.92*x[4])) +
                      12.2043471859416*x[3]*log(x[3]/(2.1055*x[2]+3.1878*x[3]+0.92*x[4]))+0.696781294644034*x[4]*log(x[4]/(2.1055*x[2]+3.1878*x[3]+0.92*x[4])) +
                      9.86*x[2]*log(x[2]/(1.972*x[2]+2.4*x[3]+1.4*x[4]))+12*x[3]*log(x[3]/(1.972*x[2]+2.4*x[3]+1.4*x[4]))+7*x[4]*log(x[4]/(1.972*x[2]+2.4*x[3] +
                      1.4*x[4]))+(1.972*x[2]+2.4*x[3]+1.4*x[4])*log(1.972*x[2]+2.4*x[3]+1.4*x[4])+1.972*x[2]*log(x[2]/(1.972*x[2]+0.283910843616504*x[3] +
                      3.02002220174195*x[4]))+2.4*x[3]*log(x[3]/(1.45991339466884*x[2]+2.4*x[3]+0.415073537580851*x[4]))+1.4*x[4]*log(x[4]/(0.602183324335333*x[2] +
                      0.115623371371275*x[3]+1.4*x[4]))-17.2981663216011*x[2]*log(x[2])-25.6043471859416*x[3]*log(x[3])-8.09678129464404*x[4]*log(x[4]))
=#
#@NLobjective(m, Min, (x[3]^3 - x[3]^3)*(x[3]^3 - x[3]^3))
@NLobjective(m, Min, x[2]*x[3]*x[4] + x[2]*x[4])
@constraint(m, e2, x[2]+x[3]+x[4] == 1.0)

m = m 		 # model get returned when including this script.

optimize!(m)
@show termination_status(m)
sout = s - time()
@show solve_time(m)
@show sout

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
# loaders.jl
# utilities for loading in data / images

using FileIO
using Colors

function ir_load_brainweb_t1_256()
	file = dirname(pathof(MIRT)) * "/../data/mri/brainweb_t1_256.gif"
	return Float32.(Gray.(load(file)))' # transpose!
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module NetworkSecurityGroupClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_NetworkSecurityGroupsApi.jl")
include("api_SecurityRulesApi.jl")

# export models
export convert, ApplicationGatewayBackendAddress
export convert, ApplicationGatewayBackendAddressPool
export convert, ApplicationGatewayBackendAddressPoolPropertiesFormat
export convert, ApplicationSecurityGroup
export convert, ApplicationSecurityGroupPropertiesFormat
export convert, BackendAddressPool
export convert, BackendAddressPoolPropertiesFormat
export convert, CloudError
export convert, CloudErrorBody
export convert, CustomDnsConfigPropertiesFormat
export convert, DdosSettings
export convert, Delegation
export convert, FlowLog
export convert, FlowLogFormatParameters
export convert, FlowLogPropertiesFormat
export convert, FrontendIPConfiguration
export convert, FrontendIPConfigurationPropertiesFormat
export convert, IPAllocationMethod
export convert, IPConfiguration
export convert, IPConfigurationProfile
export convert, IPConfigurationProfilePropertiesFormat
export convert, IPConfigurationPropertiesFormat
export convert, IPVersion
export convert, InboundNatRule
export convert, InboundNatRulePropertiesFormat
export convert, IpTag
export convert, LoadBalancerBackendAddress
export convert, LoadBalancerBackendAddressPropertiesFormat
export convert, NetworkInterface
export convert, NetworkInterfaceDnsSettings
export convert, NetworkInterfaceIPConfiguration
export convert, NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties
export convert, NetworkInterfaceIPConfigurationPropertiesFormat
export convert, NetworkInterfacePropertiesFormat
export convert, NetworkInterfaceTapConfiguration
export convert, NetworkInterfaceTapConfigurationPropertiesFormat
export convert, NetworkSecurityGroup
export convert, NetworkSecurityGroup2
export convert, NetworkSecurityGroupListResult
export convert, NetworkSecurityGroupPropertiesFormat
export convert, NetworkSecurityGroupPropertiesFormat2
export convert, PrivateEndpoint
export convert, PrivateEndpointProperties
export convert, PrivateLinkServiceConnection
export convert, PrivateLinkServiceConnectionProperties
export convert, PrivateLinkServiceConnectionState
export convert, ProvisioningState
export convert, PublicIPAddress
export convert, PublicIPAddressDnsSettings
export convert, PublicIPAddressPropertiesFormat
export convert, PublicIPAddressSku
export convert, Resource
export convert, ResourceNavigationLink
export convert, ResourceNavigationLinkFormat
export convert, RetentionPolicyParameters
export convert, Route
export convert, RouteNextHopType
export convert, RoutePropertiesFormat
export convert, RouteTable
export convert, RouteTablePropertiesFormat
export convert, SecurityRule
export convert, SecurityRule2
export convert, SecurityRuleAccess
export convert, SecurityRuleAccess2
export convert, SecurityRuleDirection
export convert, SecurityRuleDirection2
export convert, SecurityRuleListResult
export convert, SecurityRulePropertiesFormat
export convert, SecurityRulePropertiesFormat2
export convert, ServiceAssociationLink
export convert, ServiceAssociationLinkPropertiesFormat
export convert, ServiceDelegationPropertiesFormat
export convert, ServiceEndpointPolicy
export convert, ServiceEndpointPolicyDefinition
export convert, ServiceEndpointPolicyDefinitionPropertiesFormat
export convert, ServiceEndpointPolicyPropertiesFormat
export convert, ServiceEndpointPropertiesFormat
export convert, SubResource
export convert, Subnet
export convert, SubnetPropertiesFormat
export convert, TagsObject
export convert, TrafficAnalyticsConfigurationProperties
export convert, TrafficAnalyticsProperties
export convert, TransportProtocol
export convert, VirtualNetworkTap
export convert, VirtualNetworkTapPropertiesFormat

# export operations
export convert, NetworkSecurityGroupsApi, SecurityRulesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
@time using Plots
@time using LsqFit
@time using Statistics
# %%
@inline function getindexbc(grid, i)
    if checkbounds(Bool, grid, i)
        @inbounds return grid[i]
    else
        return 0
    end
end
@inline iscritical(i::Int, hc::Int) = i >= hc
getindexbc(grid) = i -> getindexbc(grid)
function evolveat(
        grid::Array{Int,2},
        hood::Array{CartesianIndex{2},1},
        index::CartesianIndex,
        hc::Int,
        topple::Array{Int},
        step::Int
    )
    res = grid[index]
    if res >= hc
        res -= 4
        @inbounds topple[step] += 1
    end
    c = count(h -> iscritical(getindexbc(grid, h + index), hc), hood)
    res = res + c
    return res
end
function evolve!(
        grid::Array{Int,2},
        simgrid::Array{Int,2},
        hood::Array{CartesianIndex{2},1},
        hc::Int,
        topple::Array{Int,1},
        step::Int
    )
    for index in CartesianIndices(grid)
        @inbounds simgrid[index] = evolveat(grid, hood, index, hc, topple, step)
    end
end
function ready_next!(grid::Array{Int,2}, hc::Int, topple::Array{Int,1}, step::Int, max_steps::Int)
    if count(iscritical.(grid, hc)) == 0
        if step == max_steps
            return false, step
        else
            xrand, yrand = rand.(Base.OneTo.(size(grid)))
            grid[xrand, yrand] += 1
            step += 1
        end
    end
    return true, step
end

function simulate(grid::Array{Int}, hc::Int, max_steps::Int)
    hood = CartesianIndex.([(-1, 0), (1, 0), (0, -1), (0, 1)])
    topple = fill(0, max_steps)
    step = 1
    st = true
    simgrid = similar(grid)
    results = Array{Int}[]
    while st
        evolve!(grid, simgrid, hood, hc, topple, step)
        grid, simgrid = simgrid, grid
        st, step = ready_next!(grid, hc, topple, step, max_steps)
    end
    return results, topple
end
# %%

hc = 4
results, topple = simulate(rand(1:hc - 1, 20, 20), 4, 10000);

# %%

function freqtable(k::Array{T}) where T
    d = Dict{T,Int}()
    for i ∈ k
        if i != 0
            d[i] = get(d, i, 0) + 1
        end
    end
    return d
end
function getlogdata(topple_freq)
    x = log.(keys(topple_freq))
    perm = sortperm(x)
    y = log.(values(topple_freq))[perm]
    x = x[perm]
    return x, y
end
function bin(xs, ys, length)
    binx = range(x[1], stop=x[end], length=length)
    _step = step(binx)
    biny = [
        mean(ys[findall(x -> val <= x < (val + _step), xs)])
        for val in binx
    ]
    binx = (binx .+ _step)[1:end - 1]
    return binx, biny
end

x, y = getlogdata(freqtable(topple))
binx, biny = bin(x, y, 100)
plot(binx, biny, st=:scatter)
# %%
@. model(x, p) = p[1] * x + p[2]
fit = curve_fit(model, x[1:45], y[1:45], [0.,0.])
println(fit.param)
xp = 0:0.01:maximum(x)
yp = model(xp, fit.param)
plot!(xp, yp)


text:
export Pose2Point2, PackedPose2Point2

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

"""
    $TYPEDEF

Bearing and Range constraint from a Pose2 to Point2 variable.
"""
struct Pose2Point2{T <: IIF.SamplableBelief} <: IIF.AbstractRelativeMinimize
    Zij::T
    partial::Tuple{Int,Int}
    # empty constructor
    # Pose2Point2{T}() where {T} = new{T}()
    # # regular constructor
    # Pose2Point2{T}(x1::T) where {T <: IIF.SamplableBelief} = new{T}(x1)
end
# convenience and default constructor
Pose2Point2(x1::T=MvNormal(zeros(2),LinearAlgebra.diagm([0.01;0.01]))) where {T <: IIF.SamplableBelief} = Pose2Point2{T}(x1, (1,2))

# prescribed sampling function
function getSample(cfo::CalcFactor{<:Pose2Point2}, N::Int=1)
  smpls = zeros(2, N)
  smpls[1:2,:] .= rand(cfo.factor.Zij, N)

  return (smpls,)
end

function IIF.getParametricMeasurement(s::Pose2Point2{<:MvNormal})

  return IIF.getParametricMeasurement(s.Zij)

end

# define the conditional probability constraint
function (cfo::CalcFactor{<:Pose2Point2})(meas,
                                          wXi,
                                          wLj )
  #

  wLj_pred = SE2(wXi)*SE2([meas;0.0])
  return  wLj .- se2vee(wLj_pred)[1:2]
end

## Serialization support

mutable struct PackedPose2Point2 <: IncrementalInference.PackedInferenceType
    Zij::String
    PackedPose2Point2() = new()
    PackedPose2Point2(s1::AS) where {AS <: AbstractString} = new(string(s1))
end

function convert(::Type{PackedPose2Point2}, obj::Pose2Point2{T}) where {T <: IIF.SamplableBelief}
  return PackedPose2Point2(convert(PackedSamplableBelief, obj.Zij))
end

# TODO -- should not be resorting to string, consider specialized code for parametric distribution types and KDEs
function convert(::Type{Pose2Point2}, packed::PackedPose2Point2)
  Pose2Point2(convert(SamplableBelief, packed.Zij))
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the Lattice-Boltzmann equations for the D2Q9 scheme

equations = LatticeBoltzmannEquations2D(Ma=0.05, Re=2000)

initial_condition = initial_condition_couette_unsteady
boundary_conditions = (
                       x_neg=boundary_condition_periodic,
                       x_pos=boundary_condition_periodic,
                       y_neg=boundary_condition_wall_noslip,
                       y_pos=boundary_condition_couette,
                      )

surface_flux = flux_lax_friedrichs
solver = DGSEM(3, surface_flux)

coordinates_min = (0, 0)
coordinates_max = (1, 1)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                periodicity=(true, false),
                n_cells_max=10_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 40.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=1000,
                                   save_final_restart=true)

# Custom solution variables: normalize velocities by reference speed `u0`
@inline function macroscopic_normalized(u, equations::LatticeBoltzmannEquations2D)
  macroscopic = cons2macroscopic(u, equations)
  rho, v1, v2, p = macroscopic

  # Use `typeof(macroscopic)` to avoid having to explicitly add `using StaticArrays`
  convert(typeof(macroscopic), (rho, v1/equations.u0, v2/equations.u0, p))
end
Trixi.varnames(::typeof(macroscopic_normalized), equations::LatticeBoltzmannEquations2D) = ("rho", "v1_normalized", "v2_normalized", "p")

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=macroscopic_normalized)

stepsize_callback = StepsizeCallback(cfl=1.0)

collision_callback = LBMCollisionCallback()

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback,
                        collision_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# Direct product of different sectors
#==============================================================================#
const SectorTuple = Tuple{Vararg{Sector}}

struct ProductSector{T<:SectorTuple} <: Sector
    sectors::T
end
_sectors(::Type{Tuple{}}) = ()
_sectors(::Type{T}) where {T<:SectorTuple} =
    (Base.tuple_type_head(T),  _sectors(Base.tuple_type_tail(T))...)

Base.IteratorSize(::Type{SectorValues{ProductSector{T}}}) where {T<:SectorTuple} =
    Base.IteratorSize(Base.Iterators.product(map(values, _sectors(T))...))
Base.@pure Base.size(::SectorValues{ProductSector{T}}) where {T<:SectorTuple} =
    map(s->length(values(s)), _sectors(T))
Base.@pure Base.length(P::SectorValues{<:ProductSector}) = *(size(P)...)

function Base.iterate(::SectorValues{ProductSector{T}}, args...) where {T<:SectorTuple}
    next = iterate(product(values.(_sectors(T))...), args...)
    next === nothing && return nothing
    val, state = next
    return ProductSector{T}(val), state
end
function Base.getindex(P::SectorValues{ProductSector{T}}, i::Int) where {T<:SectorTuple}
    Base.IteratorSize(P) isa IsInfinite &&
        throw(ArgumentError("cannot index into infinite product sector"))
    ProductSector{T}(getindex.(values.(_sectors(T)), Tuple(CartesianIndices(size(P))[i])))
end
function findindex(P::SectorValues{ProductSector{T}}, c::ProductSector{T}) where
                                                                    {T<:SectorTuple}
    Base.IteratorSize(P) isa IsInfinite &&
        throw(ArgumentError("cannot index into infinite product sector"))
    LinearIndices(size(P))[CartesianIndex(findindex.(values.(_sectors(T)), c.sectors))]
end

ProductSector{T}(args...) where {T<:SectorTuple} = ProductSector{T}(args)
Base.convert(::Type{ProductSector{T}}, t::Tuple) where {T<:SectorTuple} =
    ProductSector{T}(convert(T, t))

Base.one(::Type{ProductSector{T}}) where {G<:Sector, T<:Tuple{G}} = ProductSector((one(G),))
Base.one(::Type{ProductSector{T}}) where {G<:Sector, T<:Tuple{G,Vararg{Sector}}} =
    one(G) × one(ProductSector{Base.tuple_type_tail(T)})

Base.conj(p::ProductSector) = ProductSector(map(conj, p.sectors))
function ⊗(p1::P, p2::P) where {P<:ProductSector}
    if FusionStyle(P) isa Abelian
        (P(first(product(map(⊗, p1.sectors, p2.sectors)...))),)
    else
        return SectorSet{P}(product(map(⊗,p1.sectors,p2.sectors)...))
    end
end

Nsymbol(a::P, b::P, c::P) where {P<:ProductSector} =
    prod(map(Nsymbol, a.sectors, b.sectors, c.sectors))
function Fsymbol(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:ProductSector}
    if FusionStyle(P) isa Abelian || FusionStyle(P) isa SimpleNonAbelian
        return prod(map(Fsymbol, a.sectors, b.sectors, c.sectors,
                                    d.sectors, e.sectors, f.sectors))
    else
        # TODO: DegenerateNonAbelian case, use kron ?
        throw(MethodError(Fsymbol,(a,b,c,d,e,f)))
    end
end
function Rsymbol(a::P, b::P, c::P) where {P<:ProductSector}
    if FusionStyle(P) isa Abelian || FusionStyle(P) isa SimpleNonAbelian
        return prod(map(Rsymbol, a.sectors, b.sectors, c.sectors))
    else
        # TODO: DegenerateNonAbelian case, use kron ?
        throw(MethodError(Rsymbol,(a,b,c)))
    end
end
function Asymbol(a::P, b::P, c::P) where {P<:ProductSector}
    if FusionStyle(P) isa Abelian || FusionStyle(P) isa SimpleNonAbelian
        return prod(map(Asymbol, a.sectors, b.sectors, c.sectors))
    else
        # TODO: DegenerateNonAbelian case, use kron ?
        throw(MethodError(Asymbol,(a,b,c)))
    end
end
function Bsymbol(a::P, b::P, c::P) where {P<:ProductSector}
    if FusionStyle(P) isa Abelian || FusionStyle(P) isa SimpleNonAbelian
        return prod(map(Bsymbol, a.sectors, b.sectors, c.sectors))
    else
        # TODO: DegenerateNonAbelian case, use kron ?
        throw(MethodError(Bsymbol,(a,b,c)))
    end
end
frobeniusschur(p::ProductSector) = prod(map(frobeniusschur, p.sectors))

fusiontensor(a::ProductSector{T}, b::ProductSector{T}, c::ProductSector{T},
                v::Nothing = nothing) where {T<:SectorTuple} =
    _kron(fusiontensor(a.sectors[1], b.sectors[1], c.sectors[1]),
            fusiontensor(ProductSector(tail(a.sectors)), ProductSector(tail(b.sectors)),
                            ProductSector(tail(c.sectors))))

fusiontensor(a::ProductSector{T}, b::ProductSector{T}, c::ProductSector{T},
                v::Nothing = nothing) where {T<:Tuple{<:Sector}} =
    fusiontensor(a.sectors[1], b.sectors[1], c.sectors[1])

FusionStyle(::Type{<:ProductSector{T}}) where {T<:SectorTuple} =
    Base.:&(map(FusionStyle, _sectors(T))...)
BraidingStyle(::Type{<:ProductSector{T}}) where {T<:SectorTuple} =
    Base.:&(map(BraidingStyle, _sectors(T))...)
Base.isreal(::Type{<:ProductSector{T}}) where {T<:SectorTuple} = all(isreal, _sectors(T))

fermionparity(P::ProductSector) = mapreduce(fermionparity, xor, P.sectors)

dim(p::ProductSector) = *(dim.(p.sectors)...)

Base.isequal(p1::ProductSector, p2::ProductSector) = isequal(p1.sectors, p2.sectors)
Base.hash(p::ProductSector, h::UInt) = hash(p.sectors, h)
Base.isless(p1::ProductSector{T}, p2::ProductSector{T}) where {T} =
    isless(reverse(p1.sectors), reverse(p2.sectors))

# Default construction from tensor product of sectors
#-----------------------------------------------------
×(S1, S2, S3, S4...) = ×(×(S1, S2), S3, S4...)

×(S1::Sector, S2::Sector) = ProductSector((S1, S2))
×(P1::ProductSector, S2::Sector) = ProductSector(tuple(P1.sectors..., S2))
×(S1::Sector, P2::ProductSector) = ProductSector(tuple(S1, P2.sectors...))
×(P1::ProductSector, P2::ProductSector) =
    ProductSector(tuple(P1.sectors..., P2.sectors...))

×(G1::Type{ProductSector{Tuple{}}},
                    G2::Type{ProductSector{T}}) where {T<:SectorTuple} = G2
×(G1::Type{ProductSector{T1}},
                    G2::Type{ProductSector{T2}}) where {T1<:SectorTuple,T2<:SectorTuple} =
    tuple_type_head(T1) × (ProductSector{tuple_type_tail(T1)} × G2)
×(G1::Type{ProductSector{Tuple{}}}, G2::Type{<:Sector}) =
    ProductSector{Tuple{G2}}
×(G1::Type{ProductSector{T}}, G2::Type{<:Sector}) where {T<:SectorTuple} =
    Base.tuple_type_head(T) × (ProductSector{Base.tuple_type_tail(T)} × G2)
×(G1::Type{<:Sector}, G2::Type{ProductSector{T}}) where {T<:SectorTuple} =
    ProductSector{Base.tuple_type_cons(G1,T)}
×(G1::Type{<:Sector}, G2::Type{<:Sector}) = ProductSector{Tuple{G1,G2}}

function Base.show(io::IO, P::ProductSector)
    sectors = P.sectors
    compact = get(io, :typeinfo, nothing) === typeof(P)
    sep = compact ? ", " : " × "
    print(io,"(")
    for i = 1:length(sectors)
        i == 1 || print(io, sep)
        io2 = compact ? IOContext(io, :typeinfo => typeof(sectors[i])) : io
        print(io2, sectors[i])
    end
    print(io,")")
end

function Base.show(io::IO, ::Type{ProductSector{T}}) where {T<:SectorTuple}
    sectors = T.parameters
    print(io,"(")
    for i = 1:length(sectors)
        i == 1 || print(io, " × ")
        print(io, sectors[i])
    end
    print(io,")")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceIPConfigurationPropertiesFormat <: SwaggerModel
    virtualNetworkTaps::Any # spec type: Union{ Nothing, Vector{VirtualNetworkTap} } # spec name: virtualNetworkTaps
    applicationGatewayBackendAddressPools::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayBackendAddressPool} } # spec name: applicationGatewayBackendAddressPools
    loadBalancerBackendAddressPools::Any # spec type: Union{ Nothing, Vector{BackendAddressPool2} } # spec name: loadBalancerBackendAddressPools
    loadBalancerInboundNatRules::Any # spec type: Union{ Nothing, Vector{InboundNatRule2} } # spec name: loadBalancerInboundNatRules
    privateIPAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddress
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, IPAllocationMethod } # spec name: privateIPAllocationMethod
    privateIPAddressVersion::Any # spec type: Union{ Nothing, IPVersion } # spec name: privateIPAddressVersion
    subnet::Any # spec type: Union{ Nothing, Subnet } # spec name: subnet
    primary::Any # spec type: Union{ Nothing, Bool } # spec name: primary
    publicIPAddress::Any # spec type: Union{ Nothing, PublicIPAddress } # spec name: publicIPAddress
    applicationSecurityGroups::Any # spec type: Union{ Nothing, Vector{ApplicationSecurityGroup} } # spec name: applicationSecurityGroups
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    privateLinkConnectionProperties::Any # spec type: Union{ Nothing, NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties } # spec name: privateLinkConnectionProperties

    function NetworkInterfaceIPConfigurationPropertiesFormat(;virtualNetworkTaps=nothing, applicationGatewayBackendAddressPools=nothing, loadBalancerBackendAddressPools=nothing, loadBalancerInboundNatRules=nothing, privateIPAddress=nothing, privateIPAllocationMethod=nothing, privateIPAddressVersion=nothing, subnet=nothing, primary=nothing, publicIPAddress=nothing, applicationSecurityGroups=nothing, provisioningState=nothing, privateLinkConnectionProperties=nothing)
        o = new()
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("virtualNetworkTaps"), virtualNetworkTaps)
        setfield!(o, Symbol("virtualNetworkTaps"), virtualNetworkTaps)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("applicationGatewayBackendAddressPools"), applicationGatewayBackendAddressPools)
        setfield!(o, Symbol("applicationGatewayBackendAddressPools"), applicationGatewayBackendAddressPools)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("loadBalancerBackendAddressPools"), loadBalancerBackendAddressPools)
        setfield!(o, Symbol("loadBalancerBackendAddressPools"), loadBalancerBackendAddressPools)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("loadBalancerInboundNatRules"), loadBalancerInboundNatRules)
        setfield!(o, Symbol("loadBalancerInboundNatRules"), loadBalancerInboundNatRules)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("privateIPAddress"), privateIPAddress)
        setfield!(o, Symbol("privateIPAddress"), privateIPAddress)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("privateIPAddressVersion"), privateIPAddressVersion)
        setfield!(o, Symbol("privateIPAddressVersion"), privateIPAddressVersion)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("primary"), primary)
        setfield!(o, Symbol("primary"), primary)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("publicIPAddress"), publicIPAddress)
        setfield!(o, Symbol("publicIPAddress"), publicIPAddress)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("applicationSecurityGroups"), applicationSecurityGroups)
        setfield!(o, Symbol("applicationSecurityGroups"), applicationSecurityGroups)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(NetworkInterfaceIPConfigurationPropertiesFormat, Symbol("privateLinkConnectionProperties"), privateLinkConnectionProperties)
        setfield!(o, Symbol("privateLinkConnectionProperties"), privateLinkConnectionProperties)
        o
    end
end # type NetworkInterfaceIPConfigurationPropertiesFormat

const _property_map_NetworkInterfaceIPConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("virtualNetworkTaps")=>Symbol("virtualNetworkTaps"), Symbol("applicationGatewayBackendAddressPools")=>Symbol("applicationGatewayBackendAddressPools"), Symbol("loadBalancerBackendAddressPools")=>Symbol("loadBalancerBackendAddressPools"), Symbol("loadBalancerInboundNatRules")=>Symbol("loadBalancerInboundNatRules"), Symbol("privateIPAddress")=>Symbol("privateIPAddress"), Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("privateIPAddressVersion")=>Symbol("privateIPAddressVersion"), Symbol("subnet")=>Symbol("subnet"), Symbol("primary")=>Symbol("primary"), Symbol("publicIPAddress")=>Symbol("publicIPAddress"), Symbol("applicationSecurityGroups")=>Symbol("applicationSecurityGroups"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("privateLinkConnectionProperties")=>Symbol("privateLinkConnectionProperties"))
const _property_types_NetworkInterfaceIPConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("virtualNetworkTaps")=>"Vector{VirtualNetworkTap}", Symbol("applicationGatewayBackendAddressPools")=>"Vector{ApplicationGatewayBackendAddressPool}", Symbol("loadBalancerBackendAddressPools")=>"Vector{BackendAddressPool2}", Symbol("loadBalancerInboundNatRules")=>"Vector{InboundNatRule2}", Symbol("privateIPAddress")=>"String", Symbol("privateIPAllocationMethod")=>"IPAllocationMethod", Symbol("privateIPAddressVersion")=>"IPVersion", Symbol("subnet")=>"Subnet", Symbol("primary")=>"Bool", Symbol("publicIPAddress")=>"PublicIPAddress", Symbol("applicationSecurityGroups")=>"Vector{ApplicationSecurityGroup}", Symbol("provisioningState")=>"ProvisioningState", Symbol("privateLinkConnectionProperties")=>"NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties")
Base.propertynames(::Type{ NetworkInterfaceIPConfigurationPropertiesFormat }) = collect(keys(_property_map_NetworkInterfaceIPConfigurationPropertiesFormat))
Swagger.property_type(::Type{ NetworkInterfaceIPConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfaceIPConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ NetworkInterfaceIPConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_NetworkInterfaceIPConfigurationPropertiesFormat[property_name]

function check_required(o::NetworkInterfaceIPConfigurationPropertiesFormat)
    true
end

function validate_property(::Type{ NetworkInterfaceIPConfigurationPropertiesFormat }, name::Symbol, val)
end


text:
module Encoref
using DataFrames, JSON
using StringDistances, Statistics, DataStructures, Random, Pipe
import ParquetFiles
import Base: isless

include("utils.jl")
include("matchtypes.jl")
include("constructors.jl")
include("steps/init.jl")
include("steps/extend.jl")
include("steps/match.jl")
include("steps/integrate.jl")
include("steps/common.jl")

export CorefSystem, coref

end


text:
import PredictMD
import Test

if PredictMD.is_travis_ci()
    Test.@testset "Testing print_list_of_package_imports()" begin
        PredictMD.print_list_of_package_imports()
    end
end


text:
"""

`MayOptimize.CodeChoice` is the abstract parent type for the different types
implemented by the package `MayOptimize`.  Derived types are:

- `MayOptimize.Standard` is a singleton type intended to specify that a
  standard Julia method must be used.  This type is purposely not supported by
  `MayOptimize` macros as they should not be used with it.  This type is not
  exported by default.

- [`OptimLevel`](@ref) is an abstract type whose sub-types are used to compile
  Julia methods with specific optimizations thanks to the macros of the
  `MayOptimize` package.

"""
abstract type CodeChoice end

"""

`OptimLevel` is the abstract parent type for the different optimization levels
implemented by the macros of the package `MayOptimize`.  Derived types are:

- `Debug` for debugging or reference code that performs bound checking and no
  vectorization.

- `InBounds` for code that assumes valid indices and thus avoids
  bound checking.

- `Vectorize` for code that assumes valid indices and requires vectorization.

See macros [`@maybe_inbounds`] and [`@maybe_vectorized`] for examples.

"""
abstract type OptimLevel <: CodeChoice  end # reference code for debugging
abstract type Debug      <: OptimLevel  end # reference code for debugging
abstract type InBounds   <: OptimLevel  end # assume in-bounds
abstract type Vectorize  <: InBounds    end # to vectorize (also assume in-bounds)
@doc @doc(OptimLevel) Debug
@doc @doc(OptimLevel) InBounds
@doc @doc(OptimLevel) Vectorize

"""

`OptimLevel.AVX` is an abstract optimization level.  It is intended to use the
`@avx` macro from `LoopVectorization` if available, or `@simd` otherwise (this
is the reason to inherit from [`Vectorize`](@ref)).

For now, this type is not exported by default and using it with the
`MayOptimize` macros behaves as if [`Vectorize`](@ref) optimization level has
been selected.  Specialized methods have to be written to exploit the power of
`@avx` which can deal with several loop levels (unlike `@simd` which is
restricted to the innermost loop level).

"""
abstract type AVX <: Vectorize end

"""

`MayOptimize.Standard` is a singleton type derived from
`MayOptimize.CodeChoice` and intended to specify that standard Julia methods
should be used.  This type is purposely not supported by `MayOptimize` macros
as they should not be used with it.

This type is not exported by default.

"""
struct Standard <: CodeChoice end

max_stars_repo_name: emmt/ConditionallyOptimize


text:
struct AddVariableSpec
    variable_name::Symbol
    binary::Bool
    expression_name::Union{Nothing, Symbol}
    sign::Float64
    devices_filter_func::Union{Nothing, Function}
    initial_value_func::Union{Nothing, Function}
    lb_value_func::Union{Nothing, Function}
    ub_value_func::Union{Nothing, Function}
end

"""
Construct AddVariableSpec.

Accepts a single variable_name or a vector variable_names. One must be passed but not both.
"""
function AddVariableSpec(;
    variable_name,
    binary,
    expression_name = nothing,
    sign = 1.0,
    devices_filter_func = nothing,
    initial_value_func = nothing,
    lb_value_func = nothing,
    ub_value_func = nothing,
)
    return AddVariableSpec(
        variable_name,
        binary,
        expression_name,
        sign,
        devices_filter_func,
        initial_value_func,
        lb_value_func,
        ub_value_func,
    )
end

function AddVariableSpec(
    ::Type{<:T},
    ::Type{<:U},
    ::PSIContainer,
) where {T <: VariableType, U <: PSY.Component}
    error("AddVariableSpec is not implemented for $T / $U")
end

"""
Add variables to the PSIContainer for any component.
"""
function add_variables!(
    psi_container::PSIContainer,
    ::Type{T},
    devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}},
) where {T <: VariableType, U <: PSY.Component}
    _add_variables!(psi_container, devices, AddVariableSpec(T, U, psi_container))
end

"""
Add variables to the PSIContainer for a service.
"""
function add_variables!(
    psi_container::PSIContainer,
    ::Type{T},
    service::U,
    contributing_devices::Vector{V},
) where {T <: VariableType, U <: PSY.Reserve, V <: PSY.Device}
    spec = AddVariableSpec(T, psi_container, service)
    _add_variables!(psi_container, contributing_devices, spec)
end

function _add_variables!(
    psi_container::PSIContainer,
    devices::Union{Vector{T}, IS.FlattenIteratorWrapper{T}},
    spec::AddVariableSpec,
) where {T <: PSY.Component}
    variable_name = spec.variable_name
    binary = spec.binary
    expression_name = spec.expression_name
    sign = spec.sign
    initial_value_func = spec.initial_value_func
    lb_value_func = spec.lb_value_func
    ub_value_func = spec.ub_value_func

    filter_func = nothing
    if isnothing(spec.devices_filter_func)
        if T <: PSY.Device
            filter_func = x -> PSY.get_available(x)
        end
    else
        filter_func = spec.devices_filter_func
    end
    if !isnothing(filter_func)
        devices = filter!(filter_func, collect(devices))
    end

    add_variable!(
        psi_container,
        devices,
        variable_name,
        binary,
        expression_name,
        sign;
        initial_value = initial_value_func,
        lb_value = lb_value_func,
        ub_value = ub_value_func,
    )
end

@doc raw"""
Adds a variable to the optimization model and to the affine expressions contained
in the psi_container model according to the specified sign. Based on the inputs, the variable can
be specified as binary.

# Bounds

``` lb_value_function <= varstart[name, t] <= ub_value_function ```

If binary = true:

``` varstart[name, t] in {0,1} ```

# LaTeX

``  lb \ge x^{device}_t \le ub \forall t ``

``  x^{device}_t \in {0,1} \forall t iff \text{binary = true}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* devices : Vector or Iterator with the devices
* var_name::Symbol : Base Name for the variable
* binary::Bool : Select if the variable is binary
* expression_name::Symbol : Expression_name name stored in psi_container.expressions to add the variable
* sign::Float64 : sign of the addition of the variable to the expression_name. Default Value is 1.0

# Accepted Keyword Arguments
* ub_value : Provides the function over device to obtain the value for a upper_bound
* lb_value : Provides the function over device to obtain the value for a lower_bound. If the variable is meant to be positive define lb = x -> 0.0
* initial_value : Provides the function over device to obtain the warm start value

"""
function add_variable!(
    psi_container::PSIContainer,
    devices::D,
    var_name::Symbol,
    binary::Bool,
    expression_name::Union{Nothing, Symbol} = nothing,
    sign::Float64 = 1.0;
    kwargs...,
) where {D <: Union{Vector{<:PSY.Component}, IS.FlattenIteratorWrapper{<:PSY.Component}}}
    @assert !isempty(devices)
    time_steps = model_time_steps(psi_container)
    variable = add_var_container!(
        psi_container,
        var_name,
        [PSY.get_name(d) for d in devices],
        time_steps,
    )

    lb_f = get(kwargs, :lb_value, nothing)
    init_f = get(kwargs, :initial_value, nothing)
    ub_f = get(kwargs, :ub_value, nothing)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        variable[name, t] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "$(var_name)_{$(name), $(t)}",
            binary = binary
        )

        !isnothing(ub_f) && JuMP.set_upper_bound(variable[name, t], ub_f(d))
        !isnothing(lb_f) && !binary && JuMP.set_lower_bound(variable[name, t], lb_f(d))
        !isnothing(init_f) && JuMP.set_start_value(variable[name, t], init_f(d))

        if !(isnothing(expression_name))
            bus_number = PSY.get_number(PSY.get_bus(d))
            add_to_expression!(
                get_expression(psi_container, expression_name),
                bus_number,
                t,
                variable[name, t],
                sign,
            )
        end
    end

    return
end

@doc raw"""
Adds a bounds to a variable in the optimization model.

# Bounds

``` bounds.min <= varstart[name, t] <= bounds.max  ```

# LaTeX

``  x^{device}_t >= bound^{min;} \forall t ``

``  x^{device}_t <= bound^{max} \forall t ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* bounds::DeviceRangeConstraintInfo : contains names and vector of min / max
* var_type::AbstractString : type of the variable
* T: type of the device

"""
function set_variable_bounds!(
    psi_container::PSIContainer,
    bounds::Vector{DeviceRangeConstraintInfo},
    var_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Component}
    var = get_variable(psi_container, var_type, T)
    for t in model_time_steps(psi_container), bound in bounds
        _var = var[get_component_name(bound), t]
        JuMP.set_upper_bound(_var, bound.limits.max)
        JuMP.set_lower_bound(_var, bound.limits.min)
    end
end


text:
"""
Provides some useful functions for dealing with tree models and tree approximation using nested distance
Needs more ideas!

"""

function summary(tree0::Tree)
    println()
    println("Tree with $(maximum(stage(tree0))+1) stages and branching structure $(tree0.name[6:end]) .")
    println("The tree has $(length(tree0.state)) nodes")
    println("If applicable, the tree has a nested distance to the fitted paths of $(round(tree0.dist,digits=5)).")
end

function splitStructure(tree0::Tree)
    return(tree0.name[5:end])
end

################################
# Provide Julia wrapper functions to stochastic mortality models (hrStoMoMore)
# Provide plots of sample paths?
# Provide counter of number of nodes
# Provide estimated run time of tree_approximation! ?

# Provide shading of likely tree scenarios
# Provide scaled plots to eliminated trend? --> nicer plot, better shading available

# make use of the idea in bushiness to find better tree structure ???
# helper function to save and read trees for documentation
# create size and length methods for tree structures


text:
immutable StreamInfoBlock
    minblocksize::Int32
    maxblocksize::Int32
    minframesize::Int32
    maxframesize::Int32
    samplerate::Int32
    channels::Int32
    bitspersample::Int32
    totalsamples::Int64
end

function Base.show(io::IO, ib::StreamInfoBlock)
    println("StreamInfo metadata block")
    println(io," minblocksize: ", ib.minblocksize)
    println(io," maxblocksize: ", ib.maxblocksize)
    println(io," minframesize: ", ib.minframesize)
    println(io," maxframesize: ", ib.maxframesize)
    println(io," samplerate: ", ib.samplerate)
    println(io," channels: ", ib.channels)
    println(io," bitspersample: ", ib.bitspersample)
    println(io," totalsamples: ", ib.totalsamples)
    println(io)
end

immutable VorbisCommentEntry
    len::UInt32
    entry::Ptr{UInt8}
end

Base.show(io::IO,e::VorbisCommentEntry) =
    println(io, bytestring(e.entry,e.len))

immutable FLACVorbisComment
    vendor::VorbisCommentEntry
    n::UInt32
    comments::Ptr{VorbisCommentEntry}
end

function Base.show(io::IO,c::FLACVorbisComment)
    show(io,c.vendor)
    for cc in pointer_to_array(c.comments,c.n)
        show(io,cc)
    end
    println(io)
end

@doc "standard metadata callback function" ->
function mcallback(d::Ptr{Void}, mp::Ptr{Int32}, client::Ptr{Void})
    gv = pointer_to_array(mp,3)
    show(gv)
    println()
    println(" type: ", gv[1])
    println(" islast: ", Bool(gv[2]))
    println(" length: ", gv[3])
    if gv[1] == 0
        show(unsafe_load(reinterpret(Ptr{StreamInfoBlock},mp + 4*sizeof(Int32))))
    end
    if gv[1] == 4
        show(unsafe_load(reinterpret(Ptr{FLACVorbisComment},mp+4*sizeof(Int32))))
    end
    nothing
end

const mcallback_c = cfunction(mcallback, Void, (Ptr{Void}, Ptr{Int32}, Ptr{Void}))

@doc "error callback"->
function ecallback(d::Ptr{Void}, status::Int32, client::Ptr{Void})
    error("Got error callback with status = $status")
end

const ecallback_c = cfunction(ecallback, Void,
                              (Ptr{Void}, Int32, Ptr{Void}))

@doc "write callback" ->
function wcallback(d::Ptr{Void}, frame::Ptr{Int32},
                   buffer::Ptr{Ptr{Int32}}, client::Ptr{Void})
    fr = pointer_to_array(frame,7)
    println("Frame")
    println(" blocksize: ", fr[1])
    println(" samplerate: ", fr[2])
    println(" channels: ", fr[3])
    println(" channel assignment: ", fr[4])
    println(" bits per sample: ", fr[5])
    println(" number type: ", fr[6])
    println(" frame or sample number: ", fr[7])
    zero(Int32)
end

const wcallback_c = cfunction(wcallback, Int32,
                              (Ptr{Void}, Ptr{Int32},
                               Ptr{Ptr{Int32}}, Ptr{Void}))

function flacopen(s::AbstractString)
    dd = ccall((:FLAC__stream_decoder_new,flac),Ptr{Void},())
    ccall((:FLAC__stream_decoder_set_md5_checking,flac),Void,(Ptr{Void},Bool),dd,true)
    ccall((:FLAC__stream_decoder_set_metadata_respond_all,flac),Bool,(Ptr{Void},),dd)
    err = ccall((:FLAC__stream_decoder_init_file,flac),Int32,
                (Ptr{Void}, Ptr{UInt8}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
                dd, s, wcallback_c, mcallback_c, ecallback_c)
    dd
end

setmd5(dd::Ptr{Void}) =
    ccall((:FLAC__stream_decoder_set_md5_checking,flac),Void,(Ptr{Void},Bool),dd,true)

setrespondall(dd::Ptr{Void}) = ccall((:FLAC__stream_decoder_set_metadata_respond_all,flac),Bool,(Ptr{Void},),dd)

process_single(dd::Ptr{Void}) = ccall((:FLAC__stream_decoder_process_single,flac),Bool,(Ptr{Void},),dd)

process_metadata(dd::Ptr{Void}) = ccall((:FLAC__stream_decoder_process_until_end_of_metadata,flac),Bool,(Ptr{Void},),dd)

process_stream(dd::Ptr{Void}) = ccall((:FLAC__stream_decoder_process_until_end_of_stream,flac),Bool,(Ptr{Void},),dd)

type StreamEncoder  # type not immutable so that finalizer can be applied
    v::Ptr{Void}
end

Base.unsafe_convert(::Type{Ptr{Void}},en::StreamEncoder) = en.v

function StreamEncoder()
    en = StreamEncoder(ccall((:FLAC__stream_encoder_new,flac),Ptr{Void},()))
    finalizer(en,x->ccall((:FLAC__stream_encoder_delete,flac),Void,(Ptr{Void},),x.v))
    en
end

for (nm,typ) in (("verify",:Bool),
                 ("compression_level",:Cuint), # a value between 0 (fastest, least compression) to 8 (most compression) - default 5
                 ("channels",:Cuint),          # default 2
                 ("bits_per_sample",:Cuint),   # default 16
                 ("sample_rate",:Cuint),       # default 44100
                 ("total_samples_estimate",:Cuint), # default 0
                 ("blocksize",:Cuint))         # default 0 - encoder estimates
    @eval begin
        function $(symbol(string("set_",nm)))(en::StreamEncoder,val)
            ccall(($(string("FLAC__stream_encoder_set_",nm)),flac),Bool,(Ptr{Void},$typ),en,val)
        end
        function $(symbol(string("get_",nm)))(en::StreamEncoder)
            ccall(($(string("FLAC__stream_encoder_get_",nm)),flac),$typ,(Ptr{Void},),en)
        end
    end
end
get_state(en::StreamEncoder) = ccall((:FLAC__stream_encoder_get_state,flac),Int32,(Ptr{Void},),en)

function init_file(en::StreamEncoder,fnm::ByteString)
    ec = ccall((:FLAC__stream_encoder_init_file,flac),UInt32,
               (Ptr{Void},Ptr{Uint8},Ptr{Void},Ptr{Void}),
               en,fnm,C_NULL,C_NULL)
    ec == 0 || error("Error code $ec from stream_encoder_init_file")
    en
end

function process_interleaved(en::StreamEncoder,buf::Vector{Int32})
    nsamp = div(length(buf),get_channels(en))
    ccall((:FLAC__stream_encoder_process_interleaved,flac),Bool,
          (Ptr{Void},Ptr{Int32},Uint32),en,buf,nsamp) ||
         error("process_interleaved failed: encoder_state is $(get_state(en))")
    nothing
end

function flacwrite(fnm::AbstractString,hdr::WAVHeader,mm::Vector{Int16})
    en = StreamEncoder()
    set_compression_level(en,8)
    set_total_samples_estimate(en, hdr.dsiz >> 2)
    set_channels(en,hdr.channels)
    set_sample_rate(en,hdr.rate)
    set_bits_per_sample(en,hdr.samplebits)
    get_state(en) == 1 || error("encoder state before init_file is not 1")
    init_file(en,fnm)
    get_state(en) == 0 || error("encoder state after init_file is $(get_state(en))")
    ints_per_block = 2048
    buf = Array(Int32,ints_per_block)
    indsm = 1:ints_per_block
    for blk in 1:div(length(mm),ints_per_block)
        process_interleaved(en,copy!(buf,sub(mm,indsm)))
        indsm += ints_per_block
        rem(blk,1000) == 0 && println(lpad(string(blk),8))
    end
    indsm = (indsm.start):length(mm)
    process_interleaved(en,convert(Vector{Int32},sub(mm,indsm)))
    ccall((:FLAC__stream_encoder_finish,flac),Bool,(Ptr{Void},),en)
    show(en)
end

function Base.show(io::IO,en::StreamEncoder)
    println("FLAC Stream Encoder")
    println(" channels: ", Int(get_channels(en)))
    println(" bits per sample: ", Int(get_bits_per_sample(en)))
    println(" sample rate: ", Int(get_sample_rate(en)))
    println(" blocksize: ", Int(get_blocksize(en)))
    println(" encoder state: ", get_state(en))
    println()
end

function readinf(nm::ByteString)
    res = Dict{ASCIIString,ByteString}()
    if isfile(nm)
        open(nm,"r") do io
            for ln in eachline(io)
                vals = map(strip,split(split(chomp(ln),"#")[1],"="))
                if length(vals) == 2 && length(vals[2]) > 0
                    res[vals[1]] = strip(vals[2],'\'')
                end
            end
        end
    end
    res
end

function process_dir(dnm::ByteString)
    cd(dnm) do
        count = 0
        for nm in readdir()
            bb,ee = splitext(nm)
            if ee == ".wav"
                meta = readinf(string(bb,".inf"))
                println(meta)
                hdr,mm = WAVopen(nm)
                flacwrite(string(bb,".flac"),hdr,mm)
                count += 1
            end
        end
    end
    count
end


text:
function get_hp(obj::ParameterConfiguration)
    return jcall(obj, "getHP", jdouble, ())
end

function get_parameter_name(obj::ParameterConfiguration)
    return jcall(obj, "getParameterName", JString, ())
end


text:
# Run some non-test portion of standard script here.  See also:
# https://github.com/travis-ci/travis-build/blob/master/lib/travis/build/script/julia.rb

using Pkg

# Manually installing other CI dependencies
packages = []
if get(ENV, "TRAVIS", "") == "true"
    append!(packages, ["Coverage"])
end
specs = PackageSpec.(packages)
@info string("Installing:\n", join(specs, "\n"))
Pkg.add(specs)


text:
abstract type IndexEncoding end

Base.eltype(encoding::IndexEncoding) = eltype(encoding.indices)
Base.collect(encoding::IndexEncoding) = collect(encoding.indices)

abstract type TopologyClass end

struct Line <: TopologyClass end
struct Triangle <: TopologyClass end

parametric_dimension(::Type{Line}) = 2
parametric_dimension(::Type{Triangle}) = 3

struct Strip{C<:TopologyClass,I} <: IndexEncoding
    indices::I
end

topology_class(::Type{<:Strip{C}}) where {C} = C

Strip{C}(indices) where {C} = Strip{C,typeof(indices)}(indices)

const LineStrip = Strip{Line}
const TriangleStrip = Strip{Triangle}

struct Fan{C<:TopologyClass,I} <: IndexEncoding
    indices::I
end

topology_class(::Type{<:Fan{C}}) where {C} = C

Fan{C}(indices) where {C} = Fan{C,typeof(indices)}(indices)

const TriangleFan = Fan{Triangle}

@auto_hash_equals struct IndexList{C<:TopologyClass,Dim,T} <: IndexEncoding
    indices::Vector{SVector{Dim,T}}
end

topology_class(::Type{<:IndexList{C}}) where {C} = C

const LineList = IndexList{Line}
const TriangleList = IndexList{Triangle}

function IndexList{C}(indices) where {C<:TopologyClass}
    IndexList{C,parametric_dimension(C),eltype(eltype(indices))}(indices)
end

function TriangleList(encoding::TriangleFan, T=Int)
    indices = encoding.indices
    origin = first(indices)
    list = map(2:lastindex(indices) - 1) do i
        SVector{3,T}(origin, indices[i], indices[i + 1])
    end
    TriangleList(list)
end

function TriangleList(encoding::TriangleStrip, T=Int)
    indices = encoding.indices
    list = map(1:lastindex(indices) - 2) do i
        SVector{3,T}(indices[i], indices[i + 1 + (i - 1) % 2], indices[i + 2 - (i - 1) % 2])
    end
    TriangleList(list)
end

function LineList(encoding::LineStrip, T=Int)
    indices = encoding.indices
    LineList([SVector{2,T}(indices[i], indices[i+1]) for i in 1:length(indices) - 1])
end

abstract type MeshEncoding end

"""
Mesh encoded by vertex data and vertex indices.
"""
struct MeshVertexEncoding{I<:IndexEncoding,T}
    encoding::I
    vertex_data::Vector{T}
end

MeshVertexEncoding(vertex_data, ::Type{C}) where {C<:TopologyClass} = MeshVertexEncoding(Strip{C}(1:length(vertex_data)), vertex_data)

function MeshVertexEncoding(encoding::IndexEncoding, vertex_data)
    MeshVertexEncoding{typeof(encoding),eltype(vertex_data)}(encoding, vertex_data)
end


text:
#Function for calulating black body decay rate for 87Rb Rydberg states
#23/07/2017

function Rb87blackbody(neff,temp,ll,jj)
    #Numbers from Beterov et al. (2009)
    if ll==0
        AA = 0.134
        BB = 0.251
        CC = 2.567
        DD = 4.426
    elseif ll==1&&jj==0.5
        AA = 0.053
        BB = 0.128
        CC = 2.183
        DD = 3.989
    elseif ll==1&&jj==1.5
        AA = 0.046
        BB = 0.109
        CC = 2.085
        DD = 3.901
    elseif ll==2&&jj==1.5
        AA = 0.033
        BB = 0.084
        CC = 1.912
        DD = 3.716
    elseif ll==2&&jj==2.5
        AA = 0.032
        BB = 0.082
        CC = 1.898
        DD = 3.703
    end

    term3 = AA/(neff^DD)
    term4 = 2.14e10/(exp(315780*BB/((neff^CC)*temp))-1)
    BBdecayrate = term3*term4
    return BBdecayrate
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
abstract type AbstractStridedBitPointer{N,C,B,R,X,O} <: AbstractStridedPointer{Bool,N,C,B,R,X,O} end
struct StridedBitPointer{N,C,B,R,X,O} <: AbstractStridedBitPointer{N,C,B,R,X,O}
    p::Ptr{UInt8}
    strd::X
    offsets::O
end

# @inline function ptr_index(ptr::AbstractStridedBitPointer, i::MM{1})
#     Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >>> 3
# end
# @inline function ptr_index(ptr::AbstractStridedBitPointer, i::MM{2})
#     Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >>> 3
# end
# @inline function bitload(ptr::AbstractStridedBitPointer, i::MM{4})
#     Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >>> 3
# end
@inline function ptr_index(ptr::AbstractStridedBitPointer, i::MM{8})
    Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >>> 3
end
@inline function ptr_index(ptr::AbstractStridedBitPointer, i::MM{16})
    Base.unsafe_convert(Ptr{UInt16}, ptr.ptr), i.i >>> 3
end
@inline function ptr_index(ptr::AbstractStridedBitPointer, i::MM{32})
    Base.unsafe_convert(Ptr{UInt32}, ptr.ptr), i.i >>> 3
end
@inline function ptr_index(ptr::AbstractStridedBitPointer, i::MM{64})
    Base.unsafe_convert(Ptr{UInt64}, ptr.ptr), i.i >>> 3
end
@inline function bitload(ptr::AbstractStridedBitPointer, i::MM{W}) where {W}
    ptr, ind = ptr_index(ptr, i)
    Mask{W}(vload(ptr, ind))
end
@inline bitload(ptr::AbstractStridedBitPointer, i, ::Union{Unsigned,Mask}) = bitload(ptr, i)
@inline bitload(ptr::AbstractStridedBitPointer, i::Integer) = getindexzerobased(bitload(ptr, MM{8}(i)), i & 7)

# @inline function vstore!(ptr::AbstractStridedBitPointer, m::Mask{8}, i::Integer)
    # vstore!(Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), (m.u % Bool), i)
# end
@inline function bitstore!(ptr::AbstractStridedBitPointer, m::Mask{W}, i::MM{W}) where {W}
    ptr, ind = ptr_index(ptr, i)
    vstore!(ptr, m.u, ind)
end
@inline function bitstore!(ptr::AbstractStridedBitPointer, m::Mask{W}, i::MM{W}, mask::Mask{W}) where {W}
    ptr, ind = ptr_index(ptr, i)
    vstore!(ptr, m.u, ind)
end

@generated function bitstore!(
    ptr::Ptr{T}, v::Mask{W,U}, mask::Mask{W,U}
) where {W,T,U<:Unsigned}
    @assert isa(Aligned, Bool)
    ptyp = JuliaPointerType
    mtyp_input = LLVM_TYPES[U]
    mtyp_trunc = "i$W"
    decls = String[]
    instrs = String[]
    align = sizeof(U)
    push!(instrs, "%ptr = inttoptr $ptyp %0 to <$W x i1>*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_input %2 to <$W x i1>")
    else
        push!(instrs, "%vtrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%masktrunc = trunc $mtyp_input %2 to $mtyp_trunc")
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$W x i1>")
    end
    push!(decls,
        "declare void @llvm.masked.store.v$(W)i1(<$W x i1>, <$W x i1>*, i32, <$W x i1>)"
    )
    push!(instrs,
        "call void @llvm.masked.store.v$(W)i1(<$W x i1> %v, <$W x i1>* %ptr, i32 $align, <$W x i1> %mask)"
    )
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, $U, $U},
            ptr, v.u, mask.u)
    end
end
@generated function bitstore!(
    ptr::Ptr{T}, v::Mask{W,U}, ind::I, mask::Mask{W,U}
) where {W,T,I<:Integer,U<:Unsigned}
    @assert isa(Aligned, Bool)
    ptyp = JuliaPointerType
    mtyp_input = LLVM_TYPES[U]
    mtyp_trunc = "i$W"
    decls = String[]
    instrs = String[]
    align = sizeof(U)
    push!(instrs, "%ptr = inttoptr $ptyp %0 to i8*")
    push!(instrs, "%offsetptri8 = getelementptr inbounds i8, i8* %ptr, i$(8sizeof(I)) %2")
    push!(instrs, "offsetptr = bitcast i* %offsetptri8 to <$W x i1>*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_input %3 to <$W x i1>")
    else
        push!(instrs, "%vtrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%masktrunc = trunc $mtyp_input %3 to $mtyp_trunc")
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$W x i1>")
    end
    push!(decls,
        "declare void @llvm.masked.store.v$(W)i1(<$W x i1>, <$W x i1>*, i32, <$W x i1>)"
    )
    push!(instrs,
        "call void @llvm.masked.store.v$(W)i1(<$W x i1> %v, <$W x i1>* %offsetptr, i32 $align, <$W x i1> %mask)"
    )
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, $U, $I, $U},
            ptr, v.u, ind, mask.u)
    end
end

@inline vload(ptr::AbstractStridedBitPointer, i::Tuple) = bitload(ptr, offset(ptr, vadd_fast(i, ptr.offsets)))
@inline vload(ptr::AbstractStridedBitPointer, i::Tuple, ::Mask) = vload(ptr, i)
# @inline function vload(bptr::PackedStridedBitPointer{1}, (i,j)::Tuple{MM{W},<:Any}) where {W}
#     j = vadd_fast(j, bptr.offsets[2])
#     s = bptr.strides[1]
#     # shift = vmul_fast(s, j) & (W - 1)
#     U = mask_type(Val{W}())
#     UW = widen(U)
#     indbits = vadd_fast(vadd_fast(i.i, bptr.offsets[1]), vmul_fast(j,s))
#     ptr, ind = ptr_index(bptr, MM{W}(indbits))
#     u = vload(Base.unsafe_convert(Ptr{UW}, gepbyte(ptr, ind)))
#     shift = indbits & 7
#     # @show ind, shift, u
#     Mask{W}((u >>> shift) % U)
# end

# @inline getind(a::PackedStridedBitPointer{0}) = a.offsets[1]
# @inline getind(a::PackedStridedBitPointer{1}) = vadd_fast(a.offsets[1], vmul_fast(a.offsets[2],a.strides[1]))
# @inline Base.:(≥)(a::PackedStridedBitPointer, b::PackedStridedBitPointer) = getind(a) ≥ getind(b)
# @inline Base.:(≤)(a::PackedStridedBitPointer, b::PackedStridedBitPointer) = getind(a) ≤ getind(b)
# @inline Base.:(>)(a::PackedStridedBitPointer, b::PackedStridedBitPointer) = getind(a) > getind(b)
# @inline Base.:(<)(a::PackedStridedBitPointer, b::PackedStridedBitPointer) = getind(a) < getind(b)
# @inline Base.:(==)(a::PackedStridedBitPointer, b::PackedStridedBitPointer) = getind(a) == getind(b)

# @inline function vstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, (i,j)::Tuple{MM{W},<:Integer}) where {W}
#     j -= 1
#     s = bptr.strides[1]
#     shift = (s * j) & (W - 1)
#     U = mask_type(Val{W}())
#     UW = widen(U)
#     ptr, ind = ptr_index(bptr, MM{W}(i.i - 1 + j*s))
#     um = ((vload(Base.unsafe_convert(Ptr{UW}, gep(ptr, ind))) )
#     u = (v.u % UW) << shift
#     # @show ind, shift, u
#     vstore!(Base.unsafe_convert(Ptr{UW}, gep(ptr, ind)), u | um)

# end
# @inline vstore!(bptr::PackedStridedBitPointer{1}, v::Vec{W,Bool}, i::Tuple{MM{W},<:Integer}) where {W} = vstore!(bptr, tomask(v), i)
# @inline vstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, i::Tuple{MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, v, i)
# @inline vstore!(bptr::PackedStridedBitPointer{1}, v::Vec{W,Bool}, i::Tuple{MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, tomask(v), i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, i::Tuple{MM{W},<:Integer}) where {W} = vstore!(bptr, v, i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, i::Tuple{MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, v, i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::Vec{W,Bool}, i::Tuple{MM{W},<:Integer}) where {W} = vstore!(bptr, tomask(v), i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::Vec{W,Bool}, i::Tuple{MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, tomask(v), i)

@inline vstore!(ptr::AbstractStridedBitPointer, v::Mask, i::Tuple) = bitstore!(ptr, v, offset(ptr, vadd_fast(i, ptr.offsets)))
@inline vstore!(ptr::AbstractStridedBitPointer, v::Mask, i::Tuple, u::AbstractMask) = bitstore!(ptr, v, offset(ptr, vadd_fast(i, ptr.offsets)), tomask(u))
@inline vnoaliasstore!(ptr::AbstractStridedBitPointer, v::Mask, i::Tuple) = bitstore!(ptr, v, offset(ptr, vadd_fast(i, ptr.offsets)))
@inline vnoaliasstore!(ptr::AbstractStridedBitPointer, v::Mask, i::Tuple, u::AbstractMask) = bitstore!(ptr, v, offset(ptr, vadd_fast(i, ptr.offsets)), tomask(u))
@inline vstore!(ptr::AbstractStridedBitPointer, v::Vec{<:Any,Bool}, i::Tuple) = bitstore!(ptr, tomask(v), offset(ptr, vadd_fast(i, ptr.offsets)))
@inline vstore!(ptr::AbstractStridedBitPointer, v::Vec{<:Any,Bool}, i::Tuple, u::AbstractMask) = bitstore!(ptr, tomask(v), offset(ptr, vadd_fast(i, ptr.offsets)), tomask(u))
@inline vnoaliasstore!(ptr::AbstractStridedBitPointer, v::Vec{<:Any,Bool}, i::Tuple) = bitstore!(ptr, tomask(v), offset(ptr, vadd_fast(i, ptr.offsets)))
@inline vnoaliasstore!(ptr::AbstractStridedBitPointer, v::Vec{<:Any,Bool}, i::Tuple, u::AbstractMask) = bitstore!(ptr, tomask(v), offset(ptr, vadd_fast(i, ptr.offsets)), tomask(u))