text:
using BinDeps
using Compat

@BinDeps.setup

ecos = library_dependency("ecos", aliases=["libecos"])

if is_apple()
    using Homebrew
    provides( Homebrew.HB, "ecos", ecos, os = :Darwin )
end

version = "2.0.5"
win_version = "2.0.2"
provides(Sources, URI("https://github.com/ifa-ethz/ecos/archive/v$version.tar.gz"),
    [ecos], os = :Unix, unpacked_dir="ecos-$version")

prefix = joinpath(BinDeps.depsdir(ecos),"usr")
srcdir = joinpath(BinDeps.depsdir(ecos),"src","ecos-$version")

provides(Binaries, URI("https://cache.julialang.org/https://bintray.com/artifact/download/tkelman/generic/ecos-$win_version.7z"),
    [ecos], unpacked_dir="usr/bin$(Sys.WORD_SIZE)", os = :Windows,
    SHA="b90254220a9a63cba08700f3664519d360f45d363454e5c107e6f30e144a60a1")

# We'll keep this around for emergencies, but OSX users should be able to use Homebrew
provides(SimpleBuild,
    (@build_steps begin
        GetSources(ecos)
        CreateDirectory(joinpath(prefix,"lib"))
        FileRule(joinpath(prefix,"lib","libecos.dylib"),@build_steps begin
            ChangeDirectory(srcdir)
            `make shared`
            `mv libecos.dylib $prefix/lib`
        end)
    end),[ecos], os = :Darwin)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(ecos)
        CreateDirectory(joinpath(prefix,"lib"))
        FileRule(joinpath(prefix,"lib","libecos.so"),@build_steps begin
            ChangeDirectory(srcdir)
            `make shared`
            `mv libecos.so $prefix/lib`
        end)
    end),[ecos], os = :Unix)

@BinDeps.install Dict(:ecos => :ecos)


text:
mutable struct MatA
    size::Int
    values::Matrix{TensorDec}
end
copy(A::MatA) = MatA(A.size,reshape([copy(A.values[j,i]) for i in 1:3 for j in 1:3],(3,3)))

function to_dense(A::MatA)::Matrix{Float64}
    d = Matrix{Float64}(undef,3*A.size^2,3*A.size^2)
    for i in 1:3, j in 1:3
        d[(j-1)*A.size^2+1:j*A.size^2,(i-1)*A.size^2+1:i*A.size^2] = to_dense(A.values[j,i])
    end
    return d
end

mutable struct MatB
    size::Int
    values::Vector{TensorDec}
end
copy(B::MatB) = MatB(B.size,[copy(B.values[i]) for i in 1:3])

function to_dense(B::MatB)::Matrix{Float64}
    d = Matrix{Float64}(undef,B.size^2,3*B.size^2)
    for i in 1:3
        d[1:B.size^2,(i-1)*B.size^2+1:i*B.size^2] = to_dense(B.values[i])
    end
    return d
end

mutable struct MatC
    size::Int
    values::Vector{TensorDec}
end
copy(C::MatC) = MatC(C.size,[copy(C.values[i]) for i in 1:3])

function to_dense(C::MatC)::Matrix{Float64}
    d = Matrix{Float64}(undef,3*C.size^2,C.size^2)
    for i in 1:3
        d[(i-1)*C.size^2+1:i*C.size^2,1:C.size^2] = to_dense(C.values[i])
    end
    return d
end

function mul(A::MatA,B::MatA,atol::Float64)::MatA
    if A.size != B.size
        throw("sizes of A and B not compatible")
    end

    C = MatA(A.size,Matrix{TensorDec}(undef,3,3))

    # Distribute error threshold
    atol2 = atol/9.0
    # atol2 = atol

    for i in 1:3, j in 1:3
        C.values[i,j] = create_zeros(A.size)
        for k in 1:3
            axpy(1.,mul(A.values[i,k],B.values[k,j],atol2),C.values[i,j],atol2)
            # reduce(C.values[i,j],atol2)
            # reduce_abs(C.values[i,j],atol2)
        end
    end
    return C
end

function mul(A::MatA,B::MatC,atol::Float64)::MatC
    if A.size != B.size
        throw("sizes of A and B not compatible")
    end

    C = MatC(A.size,Vector{TensorDec}(undef,3))

    # Distribute error threshold
    atol2 = atol/3.0
    # atol2 = atol

    for i in 1:3
        C.values[i] = create_zeros(A.size)
    end

    for i in 1:3
        for j in 1:3
            axpy(1.,mul(A.values[i,j],B.values[j],atol2),C.values[i],atol2)
            # reduce(C.values[i],atol2)
            # reduce_abs(C.values[i],atol2)
        end
    end
    return C
end

function mul(A::MatB,B::MatA,atol::Float64)::MatB
    if A.size != B.size
        throw("sizes of A and B not compatible")
    end

    C = MatB(A.size,Vector{TensorDec}(undef,3))

    # Distribute error threshold
    atol2 = atol/3.0
    # atol2 = atol

    for i in 1:3
        C.values[i] = create_zeros(A.size)
    end

    for i in 1:3
        for j in 1:3
            axpy(1.,mul(A.values[j],B.values[j,i],atol2),C.values[i],atol2)
            # reduce(C.values[i],atol2)
            # reduce_abs(C.values[i],atol2)
        end
    end

    return C
end

function mul(A::MatB,B::MatC,atol::Float64)::TensorDec
    if A.size != B.size
        throw("sizes of A and B not compatible")
    end

    C = create_zeros(A.size)
    for i in 1:3
        axpy(1.,mul(A.values[i],B.values[i],atol),C,atol)
        # reduce(C,atol)
        # reduce_abs(C,atol)
    end

    return C
end

mutable struct NSFMatCoefs
    size::Int
    A::MatA
    B::MatB
    C::MatC
    T::TensorDec
    NSFMatCoefs(size) = new(size)
    NSFMatCoefs(size,A,B,C,T) = new(size,A,B,C,T)
end

function to_dense(X::NSFMatCoefs)::Matrix{Float64}
    d = Matrix{Float64}(undef,4*X.size^2,4*X.size^2)
    d[1:3*X.size^2,1:3*X.size^2] = to_dense(X.A)
    d[3*X.size^2+1:4*X.size^2,1:3*X.size^2] = to_dense(X.B)
    d[1:3*X.size^2,3*X.size^2+1:4*X.size^2] = to_dense(X.C)
    d[3*X.size^2+1:4*X.size^2,3*X.size^2+1:4*X.size^2] = to_dense(X.T)
    return d
end

function get_rank(X::NSFMatCoefs)::Matrix{Int}
    ranks = Matrix{Int}(undef,4,4)
    for i in 1:3, j in 1:3
        ranks[j,i] = X.A.values[j,i].rank
    end
    for i in 1:3
        ranks[4,i] = X.B.values[i].rank
        ranks[i,4] = X.C.values[i].rank
    end
    ranks[4,4] = X.T.rank
    return ranks
end

function get_NSFMatCoefs(X::TensorDec,m::Int,H::Vector{Float64},G::Vector{Float64},atol::Float64)::NSFMatCoefs

    # Distribute error threshold
    # atol2 = atol/16.0
    atol2 = atol/1e1
    println("atol2 = $(atol2)")

    ut = [dwtmat(X.size,X.u[i],m,H,G) for i in 1:X.rank]
    vt = [dwtmat(X.size,X.v[i],m,H,G) for i in 1:X.rank]

    index = Dict((1,1)=>1,(2,1)=>2,(1,2)=>3,(2,2)=>4)

    # Create 4-by-4 block matrix
    mat = Matrix{TensorDec}(undef,4,4)
    for i in 1:2, j in 1:2, k in 1:2, l in 1:2
        mat[2*(k-1)+l,2*(i-1)+j] = TensorDec(X.size÷2,X.rank,X.rank,copy(X.sig),
            [ut[p][index[(k,i)]] for p in 1:X.rank],[vt[p][index[(l,j)]] for p in 1:X.rank])
        normalize(mat[2*(k-1)+l,2*(i-1)+j],atol2)
    end

    # Reduce rank of mat[1,1], then reduce ranks in other blocks wrt the largest s-value of mat[1,1]
    reduce(mat[1,1],atol2)
    maxsval = get_max_sval(mat[1,1])
    for i in 1:4, j in 1:4
        if i != 1 || j != 1
            reduce(mat[i,j],atol2,maxsig=maxsval)
        end
    end

    A = MatA(X.size÷2,Matrix{TensorDec}(undef,3,3))
    for i in 1:3, j in 1:3
        A.values[j,i] = mat[j,i]
    end
    B = MatB(X.size÷2,[mat[4,i] for i in 1:3])
    C = MatC(X.size÷2,[mat[i,4] for i in 1:3])
    T = mat[4,4]

    Y = NSFMatCoefs(X.size÷2,A,B,C,T)
end

mutable struct NSFMatrix
    size::Int
    moments::Int
    scales::Int
    coefs::Vector{NSFMatCoefs}
end

# """
# Compute inverse of MatA
# """
# function get_inv_MatA(A::MatA,atol::Float64;max_iter=10)::MatA

#     # Adense = to_dense(A)
#     # println("is A is_diagally_dominant? ",is_diagally_dominant(Adense))

#     atol2 = atol/9
#     # atol2 = atol
#     # DA =  [get_inv_diag(A.values[i,i],atol2) for i in 1:3]
#     DA =  [get_inv_diag(A.values[i,i],1e-5) for i in 1:3]

#     atol2 = atol/9
#     # atol2 = atol
#     T = MatA(A.size,Matrix{TensorDec}(undef,3,3))
#     for i in 1:3, j in 1:3
#         if i == j
#             T.values[j,i] = mul(DA[j],get_offdiag(A.values[j,i]),atol2)
#             scale(T.values[j,i],-1.)
#         else
#             T.values[j,i] = mul(DA[j],A.values[j,i],atol2)
#             scale(T.values[j,i],-1.)
#         end
#     end

#     B = copy(T)
#     for i in 1:3
#         aipy(1.,B.values[i,i])
#     end

#     T = mul(T,T,atol)
#     est = sum(sum(T.values[j,i].sig[1:T.values[j,i].rank]) for i in 1:3, j in 1:3)
#     println("est = ", est)
#     for _ in 1:max_iter
#         if est < atol
#             @goto converged
#         else
#             S = copy(T)
#             for i in 1:3
#                 aipy(1.0,S.values[i,i])
#             end
#             B = mul(B,S,atol2)
#             T = mul(T,T,atol2)
#             est = sum(sum(T.values[j,i].sig[1:T.values[j,i].rank]) for i in 1:3, j in 1:3)
#             println("est = ",est)
#         end
#     end
#     @warn "maximum iterations exceeded!"
#     @label converged

#     for i in 1:3, j in 1:3
#         B.values[j,i] = mul(B.values[j,i],DA[i],atol2)
#         reduce(B.values[j,i],atol2)
#         # fastchol(B.values[j,i],atol=atol2)
#     end

#     return B
# end

"""
Compute inverse of MatA via Schulz iterations
"""
function invMatA(A::MatA,atol::Float64;max_iter=50)::Union{Nothing,MatA}
    # Select α
    # atol2 = min(1e-3,atol*1e2)
    atol2 = 1e-4
    α = 1/snorm(A,atol2)

    # Initialize
    G0 = copy(A)
    for i in 1:3, j in 1:3
        scale(G0.values[i,j],α)
    end

    atol2 = 1e-1
    err = 1.0

    # Main iteration
    for _ in 1:max_iter

        println("atol2 = $(atol2)")

        # Check for termination
        # err_ = maximum(norm(diagvec(GA.values[j,j]).-1.0,Inf) for j in 1:3)
        err_ = max(maximum([norm(sum(diagmul(G0.values[i,j],A.values[j,i]) for j in 1:3).-1.0,Inf) for i in 1:3]),
            maximum([norm(sum(diagmul(A.values[i,j],G0.values[j,i]) for j in 1:3).-1.0,Inf) for i in 1:3]))
        # err_ = opnorm(to_dense(GA)-Matrix(1.0I,3*A.size^2,3*A.size^2),2)

        println("invMatA err is $(err_)")

        if err_ < atol
            println("reach asked error threshold $(atol), return")
            return G0
        end

        println("err reduction is $(abs(err_-err)/err)")
        atol2 = min(atol2,err_/10)
        if err_ >= err || abs(err_-err) <= 0.5*err
            atol2 /= 2
        end
        if err_ < 1e-3
            atol2 = max(err_*err_,atol/1e2)
        end

        err = err_

        # # Compute G0*A*G0 directly
        # B = mul(mul(G0,A,atol2),G0,atol2)
        # for i in 1:3, j in 1:3
        #     scale(G0.values[i,j],2.0)
        # end

        # for i in 1:3, j in 1:3
        #     axpy(-1.0,B.values[i,j],G0.values[i,j],atol2)
        #     reduce(G0.values[i,j],atol2)
        # end

        ############################################################################################################
        # Compute diagonal blocks of 2*G0-G0*A*G0
        B = MatA(A.size,Matrix{TensorDec}(undef,3,3))
        for i in 1:3
            B.values[i,i] = copy(G0.values[i,i])
            scale(B.values[i,i],2.0)
            for k in 1:3
                for l in 1:3
                    axpy(-1.0,mul(mul(G0.values[i,k],A.values[k,l],atol2),G0.values[l,i],atol2),B.values[i,i],atol2)
                end
            end
            reduce(B.values[i,i],atol2)
        end

        # Compute off-diagonal blocks of 2G0-G0*A*G0
        # Use max svalue of the diagonal blocks for truncation
        maxsval = max(get_max_sval(B.values[1,1]),
            get_max_sval(B.values[2,2]),
            get_max_sval(B.values[3,3]))
        println("invMatA: maximim svalue is $(maxsval)")

        for i in 1:3
            for j in 1:3
                if i != j
                    B.values[i,j] = copy(G0.values[i,j])
                    scale(B.values[i,j],2.0)
                    for k in 1:3
                        for l in 1:3
                            axpy(-1.0,mul(mul(G0.values[i,k],A.values[k,l],atol2),G0.values[l,j],atol2),B.values[i,j],atol2)
                        end
                    end
                    reduce(B.values[i,j],atol2,maxsig=maxsval)
                end
            end
        end

        for i in 1:3
            for j in 1:3
                G0.values[i,j] = copy(B.values[i,j])
            end
        end
        ##############################################################################################################

    end

    # No convergence
    @warn("invMatA maximum number of iteration reached")
end

function invMatAr(A::MatA,atol::Float64,rank::Matrix{Int64};max_iter=50)::Union{Nothing,MatA}
    if size(rank) != (3,3)
        throw("expected a rank matrix of size (3,3), got $(size(rank))")
    end

    # Select α
    # atol2 = min(1e-3,atol*1e2)
    atol2 = 1e-4
    α = 1/snorm(A,atol2)

    # Initialize
    G0 = copy(A)
    for i in 1:3, j in 1:3
        scale(G0.values[i,j],α)
    end

    atol2 = 1e-1
    err = 1.0

    # Main iteration
    for _ in 1:max_iter

        println("atol2 = $(atol2)")

        # Check for termination
        # Compute maximum diagonal of I-G0*A and check for termination
        err_ = max(maximum([norm(sum(diagmul(G0.values[i,j],A.values[j,i]) for j in 1:3).-1.0,Inf) for i in 1:3]),
            maximum([norm(sum(diagmul(A.values[i,j],G0.values[j,i]) for j in 1:3).-1.0,Inf) for i in 1:3]))
        # err_ = maximum(norm(diagvec(GA.values[j,j]).-1.0,Inf) for j in 1:3)
        # err_ = opnorm(to_dense(GA)-Matrix(1.0I,3*A.size^2,3*A.size^2),2)

        println("invMatA err is $(err_)")

        if err_ < atol
            println("reach asked error threshold $(atol), return")
            return G0
        end

        println("err reduction is $(abs(err_-err)/err)")
        atol2 = min(atol2,err_/2)
        if err_ >= err || abs(err_-err) <= 0.5*err
            atol2 /= 2
        end
        err = err_

        #
        B = mul(mul(G0,A,atol2),G0,atol2)

        for i in 1:3, j in 1:3
            scale(G0.values[i,j],2.0)
        end

        for i in 1:3, j in 1:3
            axpy(-1.0,B.values[i,j],G0.values[i,j],atol2)
            reduce(G0.values[i,j],atol2)
            droprank(G0.values[i,j],rank[i,j])
        end
    end

    # No convergence
    @warn("invMatAr maximum number of iteration reached")
end

"""
Second version if invMatAr
"""
function invMatAr2(A::MatA,atol::Float64;max_iter=50)::Union{Nothing,MatA}
    # Select α
    # atol2 = min(1e-3,atol*1e2)
    atol2 = 1e-4
    α = 1/snorm(A,atol2)

    # Initialize
    G0 = copy(A)
    for i in 1:3, j in 1:3
        scale(G0.values[i,j],α)
    end

    atol2 = 1e-1
    err = 1.0

    # Main iteration
    for _ in 1:max_iter

        println("atol2 = $(atol2)")

        # Compute maximum diagonal of I-G0*A and check for termination
        @time err_ = max(maximum([norm(sum(diagmul(G0.values[i,j],A.values[j,i]) for j in 1:3).-1.0,Inf) for i in 1:3]),
            maximum([norm(sum(diagmul(A.values[i,j],G0.values[j,i]) for j in 1:3).-1.0,Inf) for i in 1:3]))
        # err_ = opnorm(to_dense(GA)-Matrix(1.0I,3*A.size^2,3*A.size^2),2)

        println("invMatA err is $(err_)")

        if err_ < atol
            println("reach asked error threshold $(atol), return")
            return G0
        end

        println("err reduction is $(abs(err_-err)/err)")
        atol2 = min(atol2,err_/2)
        if err_ >= err || abs(err_-err) <= 0.5*err
            atol2 /= 2
        end
        err = err_

        # Compute diagonal blocks of 2*G0-G0*A*G0
        B = MatA(A.size,Matrix{TensorDec}(undef,3,3))
        for i in 1:3
            B.values[i,i] = copy(G0.values[i,i])
            scale(B.values[i,i],2.0)
            for k in 1:3
                for l in 1:3
                    axpy(-1.0,mul(mul(G0.values[i,k],A.values[k,l],atol2),G0.values[l,i],atol2),B.values[i,i],atol2)
                end
            end
            reduce(B.values[i,i],atol2)
        end

        # Compute off-diagonal blocks of G0*A*G0
        # Use maximum diagonal rank for truncation
        rank = maximum([B.values[i,i].rank for i in 1:3])
        println("invMatA: maximim rank is $(rank)")
        for i in 1:3
            for j in 1:3
                B.values[i,j] = copy(G0.values[i,j])
                scale(B.values[i,j],2.0)
                for k in 1:3
                    for l in 1:3
                        axpy(-1.0,mul(mul(G0.values[i,k],A.values[k,l],atol2),G0.values[l,j],atol2),B.values[i,j],atol2)
                    end
                end
                reduce(B.values[i,j],atol2)
                droprank(B.values[i,j],rank)
            end
        end

        for i in 1:3
            for j in 1:3
                G0.values[i,j] = copy(B.values[i,j])
            end
        end
    end

    # No convergence
    @warn("invMatAr2 maximum number of iteration reached")
end

# """
# Estimate the spectral norm (square root of the largest eigenvalue of A^H*A) of MatA using
# random power iteration
# """
function snorm(A::MatA,atol::Float64;max_iter::Int=500)::Union{Float64,Nothing}
    x = [rand(A.size,A.size),rand(A.size,A.size),rand(A.size,A.size)]
    xnorm = norm(x)
    s = xnorm

    for _ in 1:max_iter
        x /= xnorm
        s_ = s

        x = [sum(apply(A.values[i,j],x[j]) for j in 1:3) for i in 1:3]
        x = [sum(apply(A.values[i,j],x[j],TRANS="T") for j in 1:3) for i in 1:3]

        xnorm = norm(x)
        s = xnorm

        if xnorm <= atol || abs(s-s_) <= s_*atol
            return s
        end

    end

    # No convergence
    @warn("snorm maximum number of iteration reached")
end

function factorize(X::NSFMatCoefs,atol::Float64)::NSFMatCoefs
    Y = NSFMatCoefs(X.size)

    # atol2 = 9/16*atol
    atol2 = atol
    # Y.A = get_inv_MatA(X.A,atol2)
    Y.A = invMatA(X.A,atol2)
    maxsval = max(get_max_sval(Y.A.values[1,1]),
        get_max_sval(Y.A.values[2,2]),
        get_max_sval(Y.A.values[2,2]))
    println("factorize: max svalue is $(maxsval)")

    Y.B = copy(X.B)
    for j in 1:3
        reduce(Y.B.values[j],atol2,maxsig=maxsval)
    end

    # atol2 = 3/16*atol
    atol2 = atol
    Y.C = mul(Y.A,X.C,atol2)
    for j in 1:3
        reduce(Y.C.values[j],atol2,maxsig=maxsval)
    end

    Y.T = copy(X.T)
    # atol2 = atol/16
    atol2 = atol
    axpy(-1.,mul(X.B,Y.C,atol2),Y.T,atol2/1e5)
    reduce(Y.T,atol2,maxsig=maxsval)
    # reduce_abs(Y.T,atol2)

    return Y
end

function factorize2(X::NSFMatCoefs,atol::Float64,rank::Matrix{Int64})::NSFMatCoefs
    Y = NSFMatCoefs(X.size)

    # atol2 = 9/16*atol
    atol2 = atol
    # Y.A = get_inv_MatA(X.A,atol2)
    # Y.A = invMatA(X.A,atol2)
    Y.A = invMatAr(X.A,atol2,3*rank[1:3,1:3])
    # Y.A = invMatAr2(X.A,atol2)
    # for i in 1:3
    #     for j in 1:3
    #         droprank(Y.A.values[i,j],rank[i,j])
    #     end
    # end

    Y.B = copy(X.B)
    # for i in 1:3
    #     droprank(Y.B.values[i],rank[4,i])
    # end

    # atol2 = 3/16*atol
    atol2 = atol
    Y.C = mul(Y.A,X.C,atol2)
    # for i in 1:3
    #     droprank(Y.C.values[i],rank[i,4])
    # end

    Y.T = copy(X.T)
    # atol2 = atol/16
    atol2 = atol
    axpy(-1.,mul(X.B,Y.C,atol2),Y.T,atol2)
    reduce(Y.T,atol2)
    # droprank(Y.T,rank[4,4])

    return Y
end

function get_NSFMatrix(X::TensorDec,
        moments::Int,
        scales::Int,
        atol::Float64,
        )::NSFMatrix

    atol2 = atol

    H,G = get_qmf(moments)
    Y = NSFMatrix(X.size,moments,scales,Vector{NSFMatCoefs}(undef,scales))
    T = X
    for j in 1:scales
        Y.coefs[j] = get_NSFMatCoefs(T,moments,H,G,atol2)
        T = Y.coefs[j].T
        atol2 /= 16.
    end
    return Y
end

# function get_NSFMatrix_decomp(X::TensorDec,
#         moments::Int,
#         scales::Int,
#         atol::Float64,
#         )::NSFMatrix
#     atol2 = atol
#     H,G = get_qmf(moments)
#     Y = NSFMatrix(X.size,moments,scales,Vector{NSFMatCoefs}(undef,scales))
#     T = X

#     # At first scale compute decomposition with truncation of s-values
#     println("construct decomposition of scale 1")
#     println("rank of T is $(T.rank)")
#     println("getting NSFMatCoefs ")
#     @time coefs = get_NSFMatCoefs(T,moments,H,G,atol2)
#     println("done")
#     println("compute factorization")
#     @time Y.coefs[1] = factorize(coefs,atol2)
#     println("done")
#     T = Y.coefs[1].T

#     rank = get_rank(Y.coefs[1])
#     atol2 *= 2

#     # Start from second scale, compute decomposition with truncation of ranks
#     for j in 2:scales
#         println("construct decomposition of scale $(j)")
#         println("rank of T is $(T.rank)")
#         println("getting NSFMatCoefs ")
#         @time coefs = get_NSFMatCoefs(T,moments,H,G,atol2)
#         println("done")

#         println("compute factorization")
#         @time Y.coefs[j] = factorize2(coefs,atol2,rank)
#         println("done")
#         T = Y.coefs[j].T
#         # atol2 *= 2
#         # atol2 /= 16.
#         # atol2 /= 4

#     end
#     return Y
# end

function get_NSFMatrix_decomp(X::TensorDec,
        moments::Int,
        scales::Int,
        atol::Float64,
        )::NSFMatrix
    atol2 = atol
    H,G = get_qmf(moments)
    Y = NSFMatrix(X.size,moments,scales,Vector{NSFMatCoefs}(undef,scales))
    T = X

    # Start from second scale, compute decomposition with truncation of ranks
    for j in 1:scales

        # if j==scales
        #     atol2 /= 100.
        # end

        println("construct decomposition of scale $(j)")
        println("rank of T is $(T.rank)")
        println("getting NSFMatCoefs ")
        @time coefs = get_NSFMatCoefs(T,moments,H,G,atol2)
        println("done")

        println("compute factorization")
        @time Y.coefs[j] = factorize(coefs,atol2)
        println("done")
        T = Y.coefs[j].T
        # atol2 *= 2
        # atol2 /= 2.
        atol2 /= 4.

    end
    return Y
end

"""
Apply
"""
function apply(A::NSFMatrix,x::NSFVector)::NSFVector
    throw("unimplemented")
end

"""
Forward-backward solve M * x = b
"""
function fbsub(M::NSFMatrix,
    b::NSFVector,
    atol::Float64,
    )::NSFVector

    println("fbsub using atol = $(atol)")

    if M.size != b.size || M.moments != b.moments || M.scales != b.scales
        throw(DimensionMismatch("size of M not compatible with size of b"))
    end

    n = M.size
    moments = M.moments
    scales = M.scales

    H,G = get_qmf(moments)

    println("forward sub....")

    # Forward substitution L * y = b
    y = NSFVector(n,moments,scales,Vector{NSFVecCoefs}(undef,scales))
    n ÷= 2
    proj = NSFVecCoefs(zeros(Float64,n,n),
        zeros(Float64,n,n),
        zeros(Float64,n,n),
        zeros(Float64,n,n))

    for j in 1:scales

        y.coefs[j] = NSFVecCoefs()

        rhs = [b.coefs[j].A-proj.A,b.coefs[j].B-proj.B,b.coefs[j].C-proj.C]

        y.coefs[j].A = apply(M.coefs[j].A.values[1,1],rhs[1])+
            apply(M.coefs[j].A.values[1,2],rhs[2])+
            apply(M.coefs[j].A.values[1,3],rhs[3])

        y.coefs[j].B = apply(M.coefs[j].A.values[2,1],rhs[1])+
            apply(M.coefs[j].A.values[2,2],rhs[2])+
            apply(M.coefs[j].A.values[2,3],rhs[3])

        y.coefs[j].C = apply(M.coefs[j].A.values[3,1],rhs[1])+
            apply(M.coefs[j].A.values[3,2],rhs[2])+
            apply(M.coefs[j].A.values[3,3],rhs[3])

        if j < scales
            coefs = apply(M.coefs[j].B.values[1],y.coefs[j].A)+
                    apply(M.coefs[j].B.values[2],y.coefs[j].B)+
                    apply(M.coefs[j].B.values[3],y.coefs[j].C)+
                    proj.T
            A,B,C,T = dwtmat(n,coefs,moments,H,G)
            proj = NSFVecCoefs(A,B,C,T)
            n ÷= 2
        end
    end

    # Solve the last scale
    rhs = b.coefs[scales].T-proj.T-
        apply(M.coefs[scales].B.values[1],y.coefs[scales].A)-
        apply(M.coefs[scales].B.values[2],y.coefs[scales].B)-
        apply(M.coefs[scales].B.values[3],y.coefs[scales].C)
    y.coefs[scales].T = solve(M.coefs[scales].T,rhs,atol)

    println("backward sub....")

    # Backward substitution U * x = y
    x = NSFVector(M.size,moments,scales,Vector{NSFVecCoefs}(undef,scales))
    x.coefs[scales] = NSFVecCoefs()
    x.coefs[scales].T = y.coefs[scales].T
    for j in scales:-1:1
        x.coefs[j].A = y.coefs[j].A-apply(M.coefs[j].C.values[1],x.coefs[j].T)
        x.coefs[j].B = y.coefs[j].B-apply(M.coefs[j].C.values[2],x.coefs[j].T)
        x.coefs[j].C = y.coefs[j].C-apply(M.coefs[j].C.values[3],x.coefs[j].T)

        if j > 1
            x.coefs[j-1] = NSFVecCoefs()
            x.coefs[j-1].T = idwtmat(n,moments,H,G,x.coefs[j].A,x.coefs[j].B,x.coefs[j].C,x.coefs[j].T)
            n *= 2
        end
    end

    return x
end

ext: jl

max_stars_repo_name: xinshuoyang/FactNSF

max_issues_repo_name: xinshuoyang/FactNSF

max_forks_repo_name: xinshuoyang/FactNSF


text:
@testset "Boruvka" begin

g4 = complete_graph(4)

distmx = [
    0  1  5  6
    1  0  4  10
    5  4  0  3
    6  10  3  0
]

vec_mst = Vector{Edge}([Edge(1, 2), Edge(3, 4), Edge(2, 3)])
cost_mst = sum(distmx[src(e),dst(e)] for e in vec_mst)

max_vec_mst = Vector{Edge}([Edge(1, 4), Edge(2, 4), Edge(1, 3)])
cost_max_vec_mst = sum(distmx[src(e),dst(e)] for e in max_vec_mst)

for g in testgraphs(g4)
# Testing Boruvka's algorithm
    res1 = boruvka_mst(g,distmx)
    g1t = SimpleGraph(res1.mst)
    @test res1.weight == cost_mst
    # acyclic graphs have n - c edges
    @test nv(g1t) - length(connected_components(g1t)) == ne(g1t)
    @test nv(g1t) == nv(g)

    res2 = boruvka_mst(g,distmx,minimize=false)
    g2t = SimpleGraph(res2.mst)
    @test res2.weight == cost_max_vec_mst
    @test nv(g2t) - length(connected_components(g2t)) == ne(g2t)
    @test nv(g2t) == nv(g)

end
#second test
distmx_sec = [
    0       0       0.26    0       0.38    0       0.58    0.16
    0       0       0.36    0.29    0       0.32    0       0.19
    0.26    0.36    0       0.17    0       0       0.4     0.34
    0       0.29    0.17    0       0       0       0.52    0
    0.38    0       0       0       0       0.35    0.93    0.37
    0       0.32    0       0       0.35    0       0       0.28
    0.58    0       0.4     0.52    0.93    0       0       0
    0.16    0.19    0.34    0       0.37    0.28    0       0
]

gx = SimpleGraph(distmx_sec)

vec2 = Vector{Edge}([Edge(1, 8), Edge(2, 8), Edge(3, 4), Edge(5, 6), Edge(6, 8), Edge(3, 7), Edge(1, 3)])
weight_vec2 = sum(distmx_sec[src(e),dst(e)] for e in vec2)

max_vec2 = Vector{Edge}([Edge(1, 7), Edge(2, 3), Edge(3, 7), Edge(4, 7), Edge(5, 7), Edge(5, 6), Edge(5, 8)])
weight_max_vec2 = sum(distmx_sec[src(e),dst(e)] for e in max_vec2)

for g in testgraphs(gx)

    res3 = boruvka_mst(g, distmx_sec)
    g3t = SimpleGraph(res3.mst)
    @test res3.weight == weight_vec2
    @test nv(g3t) - length(connected_components(g3t)) == ne(g3t)
    @test nv(g3t) == nv(gx)

    res4 = boruvka_mst(g, distmx_sec,minimize = false)
    g4t = SimpleGraph(res4.mst)
    @test res4.weight == weight_max_vec2
    @test nv(g4t) - length(connected_components(g4t)) == ne(g4t)
    @test nv(g4t) == nv(gx)

end

#third test with two components

distmx_third = [
    0       0       0.26    0       0.38    0       0.58    0.16        0       0       0       0
    0       0       0.36    0.29    0       0.32    0       0.19        0       0       0       0
    0.26    0.36    0       0.17    0       0       0.4     0.34        0       0       0       0
    0       0.29    0.17    0       0       0       0.52    0           0       0       0       0
    0.38    0       0       0       0       0.35    0.93    0.37        0       0       0       0
    0       0.32    0       0       0.35    0       0       0.28        0       0       0       0
    0.58    0       0.4     0.52    0.93    0       0       0           0       0       0       0
    0.16    0.19    0.34    0       0.37    0.28    0       0           0       0       0       0
    0       0       0       0       0       0       0       0           0       1       5       6
    0       0       0       0       0       0       0       0           1       0       4       10
    0       0       0       0       0       0       0       0           5       4       0       3
    0       0       0       0       0       0       0       0           6       10      3       0
]

gd = SimpleGraph(distmx_third)

vec3 = Vector{Edge}([Edge(1, 8), Edge(2, 8), Edge(3, 4),
                    Edge(5, 6), Edge(6, 8), Edge(3, 7),
                    Edge(9, 10), Edge(11,12),
                    Edge(1, 3), Edge(10,11)])
weight_vec3 = sum(distmx_third[src(e),dst(e)] for e in vec3)

max_vec3 = Vector{Edge}([Edge(1, 7), Edge(2, 3), Edge(3, 7),
                        Edge(4, 7), Edge(5, 7), Edge(5, 6),
                        Edge(5, 8), Edge(9,12), Edge(10,12), Edge(9,11)])
weight_max_vec3 = sum(distmx_third[src(e),dst(e)] for e in max_vec3)

for g in testgraphs(gd)

    res5 = boruvka_mst(g, distmx_third)
    g5t = SimpleGraph(res5.mst)
    @test res5.weight == weight_vec3
    @test nv(g5t) - length(connected_components(g5t)) == ne(g5t)
    @test nv(g5t) == nv(gd)

    res6 = boruvka_mst(g, distmx_third,minimize = false)
    g6t = SimpleGraph(res6.mst)
    @test res6.weight == weight_max_vec3
    @test nv(g6t) - length(connected_components(g6t)) == ne(g6t)
    @test nv(g6t) == nv(gd)

end

end


text:
@with_kw struct lower
  a = nothing
  b = nothing
  c = nothing
  d = nothing
end

@with_kw struct upper
  a = nothing
  b = nothing
  c = nothing
  d = nothing
end

export lower, upper


text:
function truncate_batch(limits::Matrix{Float64}, X::AbstractMatrix)
    X1 = deepcopy(X)
    d, n = size(X1)
    for i in 1:n, j in 1:d
        X1[j, i] = max(limits[j, 1], X1[j, i])
        X1[j, i] = min(limits[j, 2], X1[j, i])
    end
    return X1
end

"""
    PolynomialThompsonSampling(d::Int64, num_arms::Int64, initial_batches::Int64,
                               retrain_freq::Vector{Int64}; <keyword arguments>)

Construct a Thompson sampling policy that uses a [`PartitionedBayesPM`](@ref) to
model the expected rewards.

# Arguments

- `d::Int64`: The number of features.
- `num_arms::Int64`: The number of available actions.
- `inital_batches::Int64`: The number of batches to sample before training the polnomial
    models.
- `retrain::Vector{Int64}`: The frequency (in terms of batches) at which the partition/basis
    selection is retrained from scratch.
- `α::Float64=1.0`: Thompson sampling inflation. `α > 1` and increasing alpha increases the
    amount of exploration.
- `Jmax::Int64=3`: The maximum degree of any polynomial region.
- `Pmax::Int64=100`: The maximum number of features in any polynomial region.
- `Kmax::Int64=500`: The maximum number of regions in the partition.
- `λ::Float64=1.0`: Prior scaling.
- `shape0::Float64=1e-3`: Inverse-gamma prior shape hyperparameter.
- `scale0::Float64=1e-3`: Inverse-gamma prior scale hyperparameter.
- `ratio::Float64=1.0`: Polynomial degrees are reduced until `size(X, 2) < ratio * length(tpbasis(d, J))`.
- `tol::Float64=1e-4`: The required increase in the model evidence to accept a split.
- `verbose_retrain::Bool=true`: Print details of the partition search.
"""
mutable struct PolynomialThompsonSampling <: AbstractPolicy
    t::Int64
    batches::Int64
    data::BanditDataset
    arms::Vector{PartitionedBayesPM}
    initial_batches::Int64
    α::Float64
    retrain::Vector{Int64}
    limits::Matrix{Float64}
    limits_cache::Matrix{Float64}

    Jmax::Int64
    Pmax::Int64
    Kmax::Int64
    λ::Float64
    shape0::Float64
    scale0::Float64
    ratio::Float64
    tol::Float64
    verbose_retrain::Bool

    function PolynomialThompsonSampling(
        d::Int64,
        num_arms::Int64,
        initial_batches::Int64,
        retrain::Vector{Int64};
        α::Float64=1.0,
        Jmax::Int64=3,
        Pmax::Int64=100,
        Kmax::Int64=500,
        λ::Float64=1.0,
        shape0::Float64=1e-3,
        scale0::Float64=1e-3,
        ratio::Float64=1.0,
        tol::Float64=1e-4,
        verbose_retrain::Bool=false,
    )
        limits = repeat([0.0 0.0], d, 1)
        limits_cache = deepcopy(limits)
        data = BanditDataset(d)
        arms = Vector{PartitionedBayesPM}(undef, num_arms)
        return new(
            0,
            0,
            data,
            arms,
            initial_batches,
            α,
            retrain,
            limits,
            limits_cache,
            Jmax,
            Pmax,
            Kmax,
            λ,
            shape0,
            scale0,
            ratio,
            tol,
            verbose_retrain,
        )
    end
end

function (pol::PolynomialThompsonSampling)(X::AbstractMatrix)
    n = size(X, 2)
    num_arms = length(pol.arms)
    actions = zeros(Int64, n)

    # Check if inital batches have been completed
    if pol.batches <= pol.initial_batches
        for i in 1:n
            actions[i] = (pol.t + i) % num_arms + 1
        end
        return actions
    end

    X1 = truncate_batch(pol.limits, X)
    thompson_samples = zeros(num_arms)
    for i in 1:n
        for a in 1:num_arms
            shape, scale = shape_scale(pol.arms[a])
            x = X1[:, i:i]
            k = locate(pol.arms[a].P, x)[1]
            varsim = rand(InverseGamma(shape, scale))
            β = pol.arms[a].models[k].lm.β
            Σ = pol.α * varsim * pol.arms[a].models[k].lm.Σ
            βsim = rand(MvNormal(β[:, 1], Σ))
            z = expand(
                x,
                pol.arms[a].models[k].basis,
                pol.arms[a].P.regions[k];
                J=pol.arms[a].models[k].J,
            )
            thompson_samples[a] = βsim' * z[:, 1]
        end
        actions[i] = argmax(thompson_samples)
    end
    return actions
end

function update!(
    pol::PolynomialThompsonSampling,
    X::AbstractMatrix,
    a::AbstractVector{<:Int},
    r::AbstractMatrix,
)
    pol.t += size(X, 2)
    pol.batches += 1
    append_data!(pol.data, X, a, r)

    for d in 1:size(X, 1)
        lower = minimum(X[d, :])
        upper = maximum(X[d, :])
        pol.limits_cache[d, 1] = min(lower, pol.limits_cache[d, 1])
        pol.limits_cache[d, 2] = max(upper, pol.limits_cache[d, 2])
    end

    if pol.batches < pol.initial_batches
        return nothing
    end
    if pol.batches in pol.retrain || pol.batches == pol.initial_batches
        for a in 1:length(pol.arms)
            pol.limits = deepcopy(pol.limits_cache)
            Xa, ra = arm_data(pol.data, a)
            pol.arms[a] = PartitionedBayesPM(
                Xa,
                ra,
                pol.limits;
                Jmax=pol.Jmax,
                Pmax=pol.Pmax,
                Kmax=pol.Kmax,
                λ=pol.λ,
                shape0=pol.shape0,
                scale0=pol.scale0,
                ratio=pol.ratio,
                tol=pol.tol,
                verbose=pol.verbose_retrain,
            )
        end
    else
        X1 = truncate_batch(pol.limits, X)
        for i in unique(a)
            fit!(pol.arms[i], X1[:, a .== i], r[:, a .== i])
        end
    end
end

function predict(policy::PolynomialThompsonSampling, X::AbstractMatrix, a::Int64)
    X1 = truncate_batch(policy.limits, X)
    return policy.arms[a](X1)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: ec2_instance_connect
using AWS.Compat
using AWS.UUIDs

"""
    send_sshpublic_key(availability_zone, instance_id, instance_osuser, sshpublic_key)
    send_sshpublic_key(availability_zone, instance_id, instance_osuser, sshpublic_key, params::Dict{String,<:Any})

Pushes an SSH public key to a particular OS user on a given EC2 instance for 60 seconds.

# Arguments
- `availability_zone`: The availability zone the EC2 instance was launched in.
- `instance_id`: The EC2 instance you wish to publish the SSH key to.
- `instance_osuser`: The OS user on the EC2 instance whom the key may be used to
  authenticate as.
- `sshpublic_key`: The public key to be published to the instance. To use it after
  publication you must have the matching private key.

"""
send_sshpublic_key(AvailabilityZone, InstanceId, InstanceOSUser, SSHPublicKey; aws_config::AbstractAWSConfig=global_aws_config()) = ec2_instance_connect("SendSSHPublicKey", Dict{String, Any}("AvailabilityZone"=>AvailabilityZone, "InstanceId"=>InstanceId, "InstanceOSUser"=>InstanceOSUser, "SSHPublicKey"=>SSHPublicKey); aws_config=aws_config)
send_sshpublic_key(AvailabilityZone, InstanceId, InstanceOSUser, SSHPublicKey, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = ec2_instance_connect("SendSSHPublicKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AvailabilityZone"=>AvailabilityZone, "InstanceId"=>InstanceId, "InstanceOSUser"=>InstanceOSUser, "SSHPublicKey"=>SSHPublicKey), params)); aws_config=aws_config)


text:
using HPAT

@acc hpat function q26(category, item_count, num_centroids, file_name)
    store_sales = DataSource(DataTable{:ss_item_sk=Int64,:ss_customer_sk=Int64}, HDF5, file_name)
    item = DataSource(DataTable{:i_item_sk=Int64,:i_category=ASCIIString,:i_class_id=Int64}, HDF5, file_name)

    sale_items = join(store_sales, item, :ss_item_sk==:i_item_sk, :ss_item_sk)

    sale_items = sale_items[:i_category==category]

    customer_i_class = aggregate(sale_items, :ss_customer_sk, :ss_item_count = length(:ss_item_sk),
                                                              :id1 = sum(:i_class_id==1),
                                                              :id2 = sum(:i_class_id==2),
                                                              :id3 = sum(:i_class_id==3),
                                                              :id4 = sum(:i_class_id==4),
                                                              :id5 = sum(:i_class_id==5),
                                                              :id6 = sum(:i_class_id==6),
                                                              :id7 = sum(:i_class_id==7),
                                                              :id8 = sum(:i_class_id==8),
                                                              :id9 = sum(:i_class_id==9),
                                                              :id10 = sum(:i_class_id==10),
                                                              :id11 = sum(:i_class_id==11),
                                                              :id12 = sum(:i_class_id==12),
                                                              :id13 = sum(:i_class_id==13),
                                                              :id14 = sum(:i_class_id==14),
                                                              :id15 = sum(:i_class_id==15))

    customer_i_class = customer_i_class[:ss_item_count>item_count]

    model = Kmeans(customer_i_class, num_centroids)
    return model
end

println(q26("food", 20, 10, "data.hdf5"))


text:
using FFMPEG

const av_load_path = joinpath(dirname(@__FILE__), "ffmpeg")

const avcodec_dir    = joinpath(av_load_path, "AVCodecs",   "v$(FFMPEG._avcodec_version().major)")
const avformat_dir   = joinpath(av_load_path, "AVFormat",   "v$(FFMPEG._avformat_version().major)")
const avutil_dir     = joinpath(av_load_path, "AVUtil",     "v$(FFMPEG._avutil_version().major)")
const swscale_dir    = joinpath(av_load_path, "SWScale",    "v$(FFMPEG._swscale_version().major)")
const avdevice_dir   = joinpath(av_load_path, "AVDevice",   "v$(FFMPEG._avdevice_version().major)")
const avfilter_dir   = joinpath(av_load_path, "AVFilters",  "v$(FFMPEG._avfilter_version().major)")
#const avresample_dir = joinpath(av_load_path, "AVResample", "v$(_avresample_version().major)")
#const swresample_dir = joinpath(av_load_path, "SWResample", "v$(_swresample_version().major)")


text:
using NullableArrays
using DataArrays

srand(1)
A = rand(5_000_000)
B = rand(Bool, 5_000_000)
mu_A = mean(A)
X = NullableArray(A)
Y = NullableArray(A, B)
D = DataArray(A)
E = DataArray(A, B)

f(x) = 5 * x
f{T<:Number}(x::Nullable{T}) =
    ifelse(isnull(x), Nullable{typeof(5 * x.value)}(), Nullable(5 * x.value))

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

function profile_reduce_methods()
    A = rand(5_000_000)
    B = rand(Bool, 5_000_000)
    X = NullableArray(A)
    Y = NullableArray(A, B)
    D = DataArray(A)
    E = DataArray(A, B)

    profile_mapreduce(A, X, Y, D, E)
    profile_reduce(A, X, Y, D, E)

    for method in (
        sum,
        prod,
        minimum,
        maximum,
    )
        (method)(A)
        (method)(X)
        (method)(D)
        println("Method: $method(A) (0 missing entries)")
        print("  for Array{Float64}:          ")
        @time((method)(A))
        print("  for NullableArray{Float64}:  ")
        @time((method)(X))
        print("  for DataArray{Float64}:      ")
        @time((method)(D))
        println()

        (method)(f, A)
        (method)(f, X)
        (method)(f, D)
        println("Method: $method(f, A) (0 missing entries)")
        print("  for Array{Float64}:          ")
        @time((method)(f, A))
        print("  for NullableArray{Float64}:  ")
        @time((method)(f, X))
        print("  for DataArray{Float64}:      ")
        @time((method)(f, D))
        println()
    end

    for method in (
        sum,
        prod,
        minimum,
        maximum,
    )
        (method)(Y)
        println("Method: $method(A) (~half missing entries, skip=false)")
        print("  for NullableArray{Float64}:  ")
        @time((method)(Y))
        (method)(E)
        print("  for DataArray{Float64}:      ")
        @time((method)(E))
        println()

        (method)(f, Y)
        println("Method: $method(f, A) (~half missing entries, skip=false)")
        print("  for NullableArray{Float64}:  ")
        @time((method)(f, Y))
        if in(method, (sum, prod))
            (method)(f, E)
            print("  for DataArray{Float64}:      ")
            @time((method)(f, E))
        else
            println("  $method(f, A::DataArray) currently incurs error")
        end
        println()
    end

    for method in (
        sum,
        prod,
        minimum,
        maximum,
    )
        (method)(Y, skipnull=true)
        println("Method: $method(A) (~half missing entries, skip=true)")
        print("  for NullableArray{Float64}:  ")
        @time((method)(Y, skipnull=true))
        (method)(E, skipna=true)
        print("  for DataArray{Float64}:      ")
        @time((method)(E, skipna=true))
        println()

        (method)(f, Y, skipnull=true)
        println("Method: $method(f, A) (~half missing entries, skip=true)")
        print("  for NullableArray{Float64}:  ")
        @time((method)(f, Y, skipnull=true))
        (method)(f, E, skipna=true)
        print("  for DataArray{Float64}:      ")
        @time((method)(f, E, skipna=true))
        println()
    end

    for method in (
        sumabs,
        sumabs2
    )
        (method)(A)
        (method)(X)
        (method)(D)
        println("Method: $method(A) (0 missing entries)")
        print("  for Array{Float64}:          ")
        @time((method)(A))
        print("  for NullableArray{Float64}:  ")
        @time((method)(X))
        print("  for DataArray{Float64}:      ")
        @time((method)(D))
        println()
    end

    for method in (
        sumabs,
        sumabs2
    )
        (method)(Y)
        println("Method: $method(A) (~half missing entries, skip=false)")
        print("  for NullableArray{Float64}:  ")
        @time((method)(Y))
        (method)(E)
        print("  for DataArray{Float64}:      ")
        @time((method)(E))
        println()
    end

    for method in (
        sumabs,
        sumabs2
    )
        (method)(Y, skipnull=true)
        println("Method: $method(A) (~half missing entries, skip=true)")
        print("  for NullableArray{Float64}:  ")
        @time((method)(Y, skipnull=true))
        (method)(E, skipna=true)
        print("  for DataArray{Float64}:      ")
        @time((method)(E, skipna=true))
        println()
    end
end

function profile_mapreduce(A, X, Y, D, E)
    println("Method: mapreduce(f, op, A) (0 missing entries)")
    mapreduce(f, Base.(:+), A)
    print("  for Array{Float64}:          ")
    @time(mapreduce(f, Base.(:+), A))
    mapreduce(f, Base.(:+), X)
    print("  for NullableArray{Float64}:  ")
    @time(mapreduce(f, Base.(:+), X))
    mapreduce(f, Base.(:+), D)
    print("  for DataArray{Float64}:      ")
    @time(mapreduce(f, Base.(:+), D))
    println()

    println("Method: mapreduce(f, op, A) (~half missing entries, skip=false)")
    mapreduce(f, Base.(:+), Y)
    print("  for NullableArray{Float64}:  ")
    @time(mapreduce(f, Base.(:+), Y))
    mapreduce(f, Base.(:+), E)
    print("  for DataArray{Float64}:      ")
    @time(mapreduce(f, Base.(:+), E))
    println()

    println("Method: mapreduce(f, op, A) (~half missing entries, skip=true)")
    mapreduce(f, Base.(:+), Y, skipnull=true)
    print("  for NullableArray{Float64}:  ")
    @time(mapreduce(f, Base.(:+), Y, skipnull=true))
    mapreduce(f, Base.(:+), E, skipna=true)
    print("  for DataArray{Float64}:      ")
    @time(mapreduce(f, Base.(:+), E, skipna=true))
    println()
end

function profile_reduce(A, X, Y, D, E)
    println("Method: reduce(f, op, A) (0 missing entries)")
    reduce(Base.(:+), A)
    print("  for Array{Float64}:          ")
    @time(reduce(Base.(:+), A))
    reduce(Base.(:+), X)
    print("  for NullableArray{Float64}:  ")
    @time(reduce(Base.(:+), X))
    reduce(Base.(:+), D)
    print("  for DataArray{Float64}:      ")
    @time(reduce(Base.(:+), D))
    println()

    println("Method: reduce(f, op, A) (~half missing entries, skip=false)")
    reduce(Base.(:+), Y)
    print("  for NullableArray{Float64}:  ")
    @time(reduce(Base.(:+), Y))
    reduce(Base.(:+), E)
    print("  for DataArray{Float64}:      ")
    @time(reduce(Base.(:+), E))
    println()

    println("Method: reduce(f, op, A) (~half missing entries, skip=true)")
    reduce(Base.(:+), Y, skipnull=true)
    print("  for NullableArray{Float64}:  ")
    @time(reduce(Base.(:+), Y, skipnull=true))
    reduce(Base.(:+), E, skipna=true)
    print("  for DataArray{Float64}:      ")
    @time(reduce(Base.(:+), E, skipna=true))
    println()
end

# # NullableArray vs. DataArray comparison
function profile_skip(skip::Bool)
    println("Comparison of skipnull/skipNA methods")
    println()
    println("f := identity, op := +")
    println("mapreduce(f, op, X; skipnull/skipNA=$skip) (0 missing entries)")

    mapreduce(identity, +, X, skipnull=skip)
    print("  for NullableArray{Float64}:  ")
    @time(mapreduce(identity, +, X, skipnull=skip))

    mapreduce(identity, +, D, skipna=skip)
    print("  for DataArray{Float64}:      ")
    @time(mapreduce(identity, +, D, skipna=skip))

    println()
    println("reduce(op, X; skipnull/skipNA=$skip) (0 missing entries)")
    reduce(+, X, skipnull=skip)
    print("  for NullableArray{Float64}:  ")
    @time(reduce(+, X, skipnull=skip))

    reduce(+, D, skipna=skip)
    print("  for DataArray{Float64}:      ")
    @time(reduce(+, D, skipna=skip))

    println()
    println("mapreduce(f, op, X; skipnull/skipNA=$skip) (~half missing entries)")
    mapreduce(identity, +, Y, skipnull=skip)
    print("  for NullableArray{Float64}:  ")
    @time(mapreduce(identity, +, Y, skipnull=skip))

    mapreduce(identity, +, E, skipna=skip)
    print("  for DataArray{Float64}:      ")
    @time(mapreduce(identity, +, E, skipna=skip))

    println()
    println("reduce(op, X; skipnull/skipNA=$skip) (~half missing entries)")
    reduce(+, Y, skipnull=skip)
    print("  for NullableArray{Float64}:  ")
    @time(reduce(+, Y, skipnull=skip))

    reduce(+, E, skipna=true)
    print("  for DataArray{Float64}:      ")
    @time(reduce(+, E, skipna=true))
    nothing
end

function profile_skip_impl()
    println("Comparison of internal skip methods:")
    println("mapreduce_impl_skipnull(f, op, X) (0 missing entries)")
    NullableArrays.mapreduce_impl_skipnull(identity, +, X)
    print("  for NullableArray{Float64}:  ")
    @time(NullableArrays.mapreduce_impl_skipnull(identity, +, X))

    DataArrays.mapreduce_impl_skipna(identity, +, D)
    print("  for DataArray{Float64}:      ")
    @time(DataArrays.mapreduce_impl_skipna(identity, +, D))

    println()
    println("mapreduce_impl_skipnull(f, op, X) (~half missing entries)")
    NullableArrays.mapreduce_impl_skipnull(identity, +, Y)
    print("  for NullableArray{Float64}:  ")
    @time(NullableArrays.mapreduce_impl_skipnull(identity, +, Y))

    DataArrays.mapreduce_impl_skipna(identity, +, E)
    print("  for DataArray{Float64}:      ")
    @time(DataArrays.mapreduce_impl_skipna(identity, +, E))
    nothing
end


text:
module WriteGmsh

include("properties.jl")
include("writeGmshv2.jl")

#from properties
export PhysicalNames, push!, createElementDict, getGmshElementType, Element
#from writeGmshv2
export writeGmshv2
end # module


text:
@doc doc"""
    schreier_sims!(sc::StabilizerChain)
Complete the `sc`, i.e. run the _full, deterministic_ Schreier-Sims algorithm on `sc`.
"""
function schreier_sims!(sc::StabilizerChain)
    depth = length(sc)
    while depth ≥ 1
        @label outer_loop

        _, S, Δ = sc[depth]
        @debug "going through orbit $(collect(Δ))"
        for β in Δ
            uᵦ = Δ[β]
            @debug "Next point in orbit" β uᵦ
            for g in S
                @debug "Scheier generator for" g
                schreier_generator = uᵦ*g*getinv(Δ, β^g)
                isone(schreier_generator) && continue
                found_new_generator = false
                @debug "a new non-trivial schreier generator" schreier_generator
                h, new_depth = sift(schreier_generator, sc)
                @debug "sifts to $h at depth $new_depth"
                if new_depth ≤ length(sc)
                    @debug "its a new strong generator at depth $new_depth !"
                    found_new_generator = true
                elseif !isone(h)
                    @debug "$new_depth > $(length(sc)): extending the chain..."
                    found_new_generator = true
                    push!(sc, firstmoved(h))
                    @debug "New base" sc.base
                end

                if found_new_generator
                    for l in (depth+1):new_depth # h fixes all points <= depth
                        @debug "recomputing Schreier trees at level $l ∈ $((depth+1):new_depth)"
                        push!(sc, h, l, recompute=true) # recompute Schreier trees
                        @debug "...done! new transversal is" sc.transversals[l]
                    end
                    @debug "restarting the procedure at depth" new_depth
                    depth = new_depth
                    @goto outer_loop
                end
            end
        end
        depth -= 1
        @debug "Checked all gens, decreasing depth to" depth
    end
    return sc
end

@doc doc"""
    schreier_sims(gens::Vector{perm}, B::Vector{<:Integer}) → StabilizerChain
Complete `gens` to strong generator set by including `B` as (partial) base.
"""
function schreier_sims(gens::AbstractVector{Perm{I}}, B::AbstractVector{I}=I[]) where I<:Integer
    sc = StabilizerChain(gens, B)
    schreier_sims!(sc)
    return sc.base, sc.sgs, sc.transversals
end


text:
#=
Run the method on all 1dsfm datasets! Save a .json file for viz, and also save the
eigenvalues, eigenvectors, and timing data for future analysis.
=#

include("../set_loadpath.jl")
using LinAlgUtils
using JacobianUtils

# JLD requires that packages whose types are being serialized be included, even
# if they aren't explicitly used otherwise in the file.
using SparseArrays
using LinearAlgebra
using Rotations

import Glob
using JSON
using JLD

########
# MAIN #
########
num_vecs = 20

dataset_root = "dataset/1dsfm"

# output location
outdir_root = "output/1dsfm"

# bin file location
convert_to_bal = normpath(joinpath(@__DIR__, "bundle2bal.jl"))
compute_jacobian = normpath(joinpath(@__DIR__, "..", "build", "bin", "compute_jacobian"))

problems = Glob.glob("$dataset_root/*")
for problem_file in problems

    dataset = splitext(basename(problem_file))[1]
    # find all of the individual problems in this dataset
    # only get post-bundle adjustment problems.
    outdir = "$outdir_root/$dataset"

    # check that the outdir exists. If not, make it.
    if !isdir(outdir)
        mkpath(outdir)
    end
    results_file = "$outdir/results.jld"
    if ~isfile(results_file)
        save(results_file, "all_runs", [])
    end

    println("Beginning computation on dataset: $dataset")

    problem_shortname = dataset
    bal_file = "$outdir/$problem_shortname.txt"
    json_file = "$outdir/$problem_shortname.json"
    jacobian_file = "$outdir/$problem_shortname.jac.h5"

    # check: does the bal file already exist? If so, skip this problem.
    if isfile(bal_file)
        println("Saved results found for $dataset/$problem_shortname. Skipping.")
        continue
    end

    # convert the bundle problem to a bal problem.
    println("Converting .out to BAL for problem $dataset/$problem_shortname")
    t_conv = @elapsed run(`julia $convert_to_bal "$problem_file" "$bal_file"`)
    println("  ---> that took $t_conv seconds.")

    # extract the jacobian from a problem instance
    println("Computing jacobian for problem $dataset/$problem_shortname")
    t_jac = @elapsed run(`$compute_jacobian "$bal_file" "$jacobian_file"`)
    println("  ---> that took $t_jac seconds.")

    # get the eigenvectors from this jacobian
    println("Computing eigenvectors for problem $dataset/$problem_shortname")
    problem = readproblem(jacobian_file)
    t_eig = @elapsed λ, X, residuals = compute_covariance_evecs(problem, num_vecs, method=:dropping)
    rm(jacobian_file)
    println("  ---> that took $t_eig seconds.")

    println("Saving output for problem $dataset/$problem_shortname")

    # TEMP: camera center calculations so we don't have to figure out how to do
    # this in Javascript
    camera_centers = zeros(Float64, 3*problem.num_cameras)
    for k = 0:problem.num_cameras-1
        rot = problem.camera_parameters[9*k+1:9*k+3]
        θ = norm(rot)
        rot_aa = AngleAxis(θ, rot[1], rot[2], rot[3])
        t = problem.camera_parameters[9*k+4:9*k+6]
        camera_centers[3*k+1 : 3*k+3] = -inv(rot_aa) * t
    end

    # save the vis information to a json file
    # stuff this into a json file
    json_data = Dict()
    json_data["num_cameras"] = problem.num_cameras
    json_data["num_points"] = problem.num_points
    json_data["camera_parameters"] = problem.camera_parameters
    json_data["point_parameters"] = problem.point_parameters
    json_data["camera_centers"] = camera_centers
    evecs = []
    for i=1:size(X,2)
        push!(evecs, X[:,i])
    end
    json_data["eigenvectors"] = evecs
    json_data["eigenvalues"] = λ
    open(json_file, "w") do f
        write(f, "sample_scene = ")
        JSON.print(f, json_data)
    end

    # prepare a structure with all of the results data
    data = Dict(
        "dataset_name"         => dataset,
        "problem_name"         => problem_shortname,
        "num_cameras"          => problem.num_cameras,
        "num_points"           => problem.num_points,
        "num_observations"     => problem.num_observations,
        "eigenvectors"         => evecs,
        "eigenvalues"          => λ,
        "computation_time"     => t_eig,
        "residuals"            => residuals,
        "source_problem_file"  => problem_file
    )
    # read the records file, append the new data, and push back to file
    all_records = load(results_file, "all_runs")
    rm(results_file)
    push!(all_records, data)
    save(results_file, "all_runs", all_records)

    println("")

end

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
init(𝒜::ASB07_decomposed, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒜, 𝑆, copy(𝑂))

function init!(𝒜::ASB07_decomposed, 𝑆::AbstractSystem, 𝑂::Options)
    𝑂copy = copy(𝑂)
    𝑂copy[:n] = statedim(𝑆)

    𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy)

    𝑂validated[:partition] = 𝒜.options[:partition]

    # compute blocks
    𝑂validated[:blocks] =
        compute_blocks(𝒜.options[:vars], 𝑂validated[:partition])

    return 𝑂validated
end


text:
using SparseArrays # blockdiag() function
using LinearAlgebra # Identity Matrix
function sdp_ymat(lines, Ybus)

nbus = max(maximum([lines[l].fbus for l in lineset]), maximum([lines[l].tbus for l in lineset]))
e = (k) -> sparse(Matrix(I, nbus, nbus)[:,k])
Yk_small = (k) -> e(k)*e(k)'*Ybus

Yk = (k) -> (1/2)*
    [real(Yk_small(k) + conj.(Yk_small(k)')) imag(conj.(Yk_small(k)') - Yk_small(k));
    imag(Yk_small(k) - conj.(Yk_small(k)')) real(Yk_small(k) + conj.(Yk_small(k)'))]

Yk_ = (k) -> -(1/2)*
    [imag(Yk_small(k) + conj.(Yk_small(k)')) real(Yk_small(k) - conj.(Yk_small(k)'));
    real(conj.(Yk_small(k)') - Yk_small(k)) imag(Yk_small(k) + conj.(Yk_small(k)'))]

Mk = (k) -> blockdiag( e(k)*e(k)', e(k)*e(k)' )

gl = (l) -> real( 1 / (lines[l].r+im*lines[l].x)); # Real part of line admittance
bl = (l) -> imag( 1 / (lines[l].r+im*lines[l].x)); # Imaginary part of line admittance

tau = (l) -> ifelse(lines[l].tap == 0, 1, lines[l].tap);
theta = (l) -> lines[l].shft
gbcosft = (l) -> gl(l)*cos(theta(l)) + bl(l)*cos(theta(l)+pi/2)
gbsinft = (l) -> gl(l)*sin(theta(l)) + bl(l)*sin(theta(l)+pi/2)
gbcostf = (l) -> gl(l)*cos(-theta(l)) + bl(l)*cos(-theta(l)+pi/2)
gbsintf = (l) -> gl(l)*sin(-theta(l)) + bl(l)*sin(-theta(l)+pi/2)

Ylineft = (l) -> 0.5*(
    sparse(
        [lines[l].fbus, lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].fbus+nbus],
        [lines[l].fbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [gl(l)/(tau(l)^2), -gbcosft(l)/tau(l), gbsinft(l)/tau(l), gl(l)/(tau(l)^2), -gbsinft(l)/tau(l), -gbcosft(l)/tau(l)],
                2*nbus, 2*nbus)
    +
    sparse(
        [lines[l].fbus, lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].fbus+nbus],
        [lines[l].fbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [gl(l)/(tau(l)^2), -gbcosft(l)/tau(l), gbsinft(l)/tau(l), gl(l)/(tau(l)^2), -gbsinft(l)/tau(l), -gbcosft(l)/tau(l)],
            2*nbus,2*nbus)')

Y_lineft = (l) -> 0.5*(
    sparse(
        [lines[l].fbus, lines[l].fbus,  lines[l].fbus,  lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].fbus+nbus ],
        [lines[l].fbus, lines[l].tbus,  lines[l].tbus+nbus, lines[l].fbus+nbus, lines[l].tbus,  lines[l].tbus+nbus ],
        [-(bl(l)+lines[l].b/2)/(tau(l)^2), gbsinft(l)/tau(l), gbcosft(l)/tau(l), -(bl(l)+lines[l].b/2)/(tau(l)^2), -gbcosft(l)/tau(l), gbsinft(l)/tau(l)],
            2*nbus,2*nbus)
    +
    sparse(
        [lines[l].fbus, lines[l].fbus, lines[l].fbus,  lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].fbus+nbus ],
        [lines[l].fbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [-(bl(l)+lines[l].b/2)/(tau(l)^2), gbsinft(l)/tau(l), gbcosft(l)/tau(l), -(bl(l)+lines[l].b/2)/(tau(l)^2), -gbcosft(l)/tau(l), gbsinft(l)/tau(l)],
            2*nbus,2*nbus)')

Ylinetf = (l) -> 0.5*(
    sparse(
        [lines[l].fbus, lines[l].fbus,  lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus ],
        [lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus,  lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus ],
        [-gbcostf(l)/tau(l),   -gbsintf(l)/tau(l),    gbsintf(l)/tau(l),  -gbcostf(l)/tau(l), gl(l), gl(l)],
            2*nbus,2*nbus)
    +
    sparse(
        [lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [-gbcostf(l)/tau(l), -gbsintf(l)/tau(l), gbsintf(l)/tau(l), -gbcostf(l)/tau(l), gl(l), gl(l)],
            2*nbus,2*nbus)')

Y_linetf = (l) -> 0.5*(
    sparse(
        [lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [gbsintf(l)/tau(l),   -gbcostf(l)/tau(l),    gbcostf(l)/tau(l),     gbsintf(l)/tau(l),     -(bl(l)+lines[l].b/2), -(bl(l)+lines[l].b/2)],
            2*nbus,2*nbus)
    +
    sparse(
        [lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
        [gbsintf(l)/tau(l),   -gbcostf(l)/tau(l),    gbcostf(l)/tau(l),     gbsintf(l)/tau(l),     -(bl(l)+lines[l].b/2),  -(bl(l)+lines[l].b/2)],
                2*nbus,2*nbus)')

YL = (l) -> sparse(
    [lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus+nbus],
    [lines[l].fbus, lines[l].tbus, lines[l].fbus+nbus, lines[l].tbus+nbus, lines[l].fbus, lines[l].tbus, lines[l].fbus+nbus, lines[l].tbus+nbus],
    [1, -1, 1, -1, -1, 1, -1, 1],
        2*nbus,2*nbus) * lines[l].r * ( gl(l)^2 + bl(l)^2 )

YL_ = (l) -> (sparse(
    [lines[l].fbus, lines[l].fbus, lines[l].fbus+nbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus, lines[l].tbus+nbus, lines[l].tbus+nbus],
    [lines[l].fbus, lines[l].tbus, lines[l].fbus+nbus, lines[l].tbus+nbus, lines[l].fbus, lines[l].tbus, lines[l].fbus+nbus, lines[l].tbus+nbus],
    [1, -1, 1, -1, -1, 1, -1, 1],
        2*nbus,2*nbus) * lines[l].x * ( gl(l)^2 + bl(l)^2 )
    -
    sparse(
    [lines[l].fbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
    [lines[l].fbus, lines[l].fbus+nbus, lines[l].tbus, lines[l].tbus+nbus],
    [1, 1, 1, 1],
        2*nbus, 2*nbus) * lines[l].b / 2)

    return Yk, Yk_, Mk, Ylineft, Ylinetf, Y_lineft, Y_linetf, YL, YL_
end

max_stars_repo_name: JipKim/ACOPF

max_issues_repo_name: JipKim/ACOPF

max_forks_repo_name: JipKim/ACOPF


text:
@testset "MultiobjCoverProblem" begin
    @testset "empty" begin
        sm = SetMosaic(Set{Symbol}[])
        msm = mask(sm, [Set{Symbol}()])
        problem = MultiobjCoverProblem(msm)
        @test nvars(problem) == 0
        @test_skip nmasks(problem) == 1
        @test_throws DimensionMismatch score([1.0], problem)
        @test OESC.miscover_score(Float64[], problem) == (0.0, 0.0)
        @test score(Float64[], problem) == (0.0, 0.0, 0.0, 0.0)
        @test aggscore(Float64[], problem) == 0.0

        res = optimize(problem)
        @test isempty(res)
        @test nvars(res) == 0
        @test nsolutions(res) == 0
        @test best_index(res) == 0
    end

    @testset "[a]" begin # FIXME take element detection probability into account
        # empty mask problem
        empty_prob = MultiobjCoverProblem(mask(SetMosaic([Set([:a])]), [Set{Symbol}()]))
        @test nvars(empty_prob) == 0
        @test OESC.miscover_score(Float64[], empty_prob) == (0.0, 0.0)
        @test score(Float64[], empty_prob) == (0.0, 0.0, 0.0, 0.0)
        @test aggscore(Float64[], empty_prob) == 0.0
        empty_res = optimize(empty_prob)
        @test isempty(empty_res)
        @test nvars(empty_res) == 0
        @test nsolutions(empty_res) == 0

        # enabled because sel prob is high, and although :a is all elements it needs to be covered
        prob1_mosaic = mask(SetMosaic([Set([:a])], Set([:a])), [Set([:a])], min_nmasked=1)
        en1_prob = MultiobjCoverProblem(prob1_mosaic, CoverParams(sel_tax=-log(0.5), uncovered_factor=1.0))
        @test nvars(en1_prob) == 1
        @test_throws DimensionMismatch score(Float64[], en1_prob)
        @test OESC.miscover_score([0.0], en1_prob) == (1.0, 0.0)
        @test OESC.miscover_score([1.0], en1_prob) == (0.0, 0.0)
        en1_score = score([1.0], en1_prob)
        @test en1_score == (en1_prob.var_scores[1], 0.0, 0.0, 0.0)
        @test score([0.0], en1_prob) == (0.0, 0.0, 1.0, 0.0)
        @test en1_score[1] < 1.0 # score[3] at [0.0], required for enabling
        @test score([0.5], en1_prob) == (0.5*en1_prob.var_scores[1], 0.0, 0.5, 0.0)
        @test aggscore([1.0], en1_prob) < aggscore([0.9], en1_prob) <
              aggscore([0.5], en1_prob) < aggscore([0.0], en1_prob)
        en1_res = optimize(en1_prob, MultiobjOptimizerParams(ϵ=0.001))
        @test nsolutions(en1_res) > 0
        @test nvars(en1_res) == 1
        @test best_varweights(en1_res) == ones(Float64, 1)

        # disabled because :a is all elements and selection probability low
        dis1_prob = MultiobjCoverProblem(prob1_mosaic, CoverParams(sel_tax=-log(0.01), uncovered_factor=1.0))
        @test nvars(dis1_prob) == 1
        @test dis1_prob.var_scores[1] > en1_prob.var_scores[1]
        dis1_score = score([1.0], dis1_prob)
        @test dis1_score == (dis1_prob.var_scores[1], 0.0, 0.0, 0.0)
        @test score([0.0], dis1_prob) == (0.0, 0.0, 1.0, 0.0)
        @test dis1_score[1] > 1.0 # score[3] at [0.0], required for disabling
        @test aggscore([1.0], dis1_prob) > aggscore([0.5], dis1_prob) > aggscore([0.0], dis1_prob)
        dis1_res = optimize(dis1_prob)
        @test nsolutions(dis1_res) > 0
        @test nvars(dis1_res) == 1
        @test best_varweights(dis1_res) == zeros(Float64, 1)

        # enabled because there is :a and :b
        en2_prob = MultiobjCoverProblem(mask(SetMosaic([Set([:a])], Set([:a, :b])),
                                             [Set([:a])], min_nmasked=1),
                                        CoverParams(sel_tax=-log(0.9)))
        @test nvars(en2_prob) == 1
        en2_res = optimize(en2_prob)
        @test nsolutions(en2_res) > 0
        @test nvars(en2_res) == 1
        @test best_varweights(en2_res) == ones(Float64, 1)
    end

    @testset "[a b] [c d] [a b c d], mask=[a b]" begin # FIXME take weights into account
        sm = SetMosaic([Set([:a, :b]), Set([:c, :d]), Set([:a, :b, :c])])
        sm_ab = mask(sm, [Set([:a, :b])], min_nmasked=1)

        problem_def = MultiobjCoverProblem(sm_ab) # cd set is dropped

        @test OESC.var_scores([0.0, 0.0], problem_def) == (0.0, 0.0)
        @test OESC.var_scores([1.0, 0.0], problem_def)[2] == 0.0
        @test OESC.var_scores([0.0, 1.0], problem_def)[2] == 0.0
        @test OESC.var_scores([1.0, 1.0], problem_def)[2] > 0

        @test OESC.miscover_score([0.0, 0.0], problem_def) == (2.0, 0.0)
        @test OESC.miscover_score([0.25, 0.0], problem_def) == (1.5, 0.0)
        @test OESC.miscover_score([0.75, 0.0], problem_def) == (0.5, 0.0)
        @test OESC.miscover_score([1.0, 0.0], problem_def) == (0.0, 0.0)
        @test OESC.miscover_score([0.0, 1.0], problem_def) == (0.0, 1.0)
        @test OESC.miscover_score([1.0, 1.0], problem_def) == (0.0, 1.0)

        res_def = optimize(problem_def)
        @test nsolutions(res_def) > 0
        @test best_varweights(res_def) ≈ [1.0, 0.0] atol=0.01

        problem_no_penalty = MultiobjCoverProblem(sm_ab, CoverParams(setXset_factor=0.0, covered_factor=0.0, sel_tax=0.0))
        @test nvars(problem_no_penalty) == 2
        res_no_penalty = optimize(problem_no_penalty, MultiobjOptimizerParams(ϵ=[0.001, 0.001], WeightDigits=nothing))
        @test nvars(res_no_penalty) == 2
        @test best_varweights(res_no_penalty) ≈ [1.0, 1.0] atol=0.01
    end

    @testset "[a b] [b c] [a b c], mask=[b]" begin # FIXME take weights into account
        sm = SetMosaic([Set([:a, :b]), Set([:b, :c]), Set([:a, :b, :c])])
        sm_b = mask(sm, [Set([:b])], min_nmasked=1)

        prob_hi_sXs = MultiobjCoverProblem(sm_b, CoverParams(setXset_factor=10.0, sel_tax=-log(1E-25)))
        @test nvars(prob_hi_sXs) == 3
        res_ignore_overlap = optimize(prob_hi_sXs, MultiobjOptimizerParams(ϵ=0.001))
        @test nvars(res_ignore_overlap) == 3
        @test best_varweights(res_ignore_overlap) ≈ [0.0, 0.0, 0.0] atol=0.01

        problem_b = MultiobjCoverProblem(sm_b, CoverParams(setXset_factor=1.0, uncovered_factor=1.0, sel_tax=-log(0.5)))
        res_b = optimize(problem_b)
        @test nsolutions(res_b) > 0
        @test nvars(res_b) == 3
        # [:a :b] is as good as [:b :c]
        @test aggscore([0.0, 1.0, 0.0], problem_b) < aggscore([0.0, 0.0, 1.0], problem_b)
        @test aggscore([0.0, 1.0, 0.0], problem_b) < aggscore([0.0, 0.0, 0.0], problem_b)
        @test aggscore([0.0, 1.0, 0.0], problem_b) == aggscore([1.0, 0.0, 0.0], problem_b)
        @test max(best_varweights(res_b)[1], best_varweights(res_b)[2]) ≈ 1.0 atol=0.01
        @test best_varweights(res_b)[3] ≈ 0.0 atol = 0.01

        @test OESC.miscover_score([1.0, 0.0, 0.0], problem_b) == (0.0, 1.0)
        @test OESC.miscover_score([0.0, 1.0, 0.0], problem_b) == (0.0, 1.0)
        @test OESC.miscover_score([0.0, 0.5, 0.0], problem_b) == (0.5, 0.5)
        @test OESC.miscover_score([0.0, 0.0, 1.0], problem_b) == (0.0, 2.0)
        @test OESC.miscover_score([1.0, 0.0, 1.0], problem_b) == (0.0, 2.0)
        @test OESC.miscover_score([1.0, 1.0, 1.0], problem_b) == (0.0, 2.0)
        @test OESC.miscover_score([1.0, 1.0, 0.5], problem_b) == (0.0, 2.0)
    end

    @testset "[a b d] [b c d] [c] [d] [a b c d e] [c d e], mask=[a b c]" begin # FIXME take weights into account
        sm = SetMosaic([Set([:a, :b, :d]), Set([:b, :c, :d]), Set([:c]), Set([:d]),
                        Set([:a, :b, :c, :d, :e]), Set([:c, :d, :e, :f])],
                        Set([:a, :b, :c, :d, :e, :f]))
        sm_abc = mask(sm, [Set([:a, :b, :c])], min_nmasked=1)

        # low prior probability to select sets, high probability to miss active element, so select abc
        prob_hi_sXs = MultiobjCoverProblem(sm_abc, CoverParams(setXset_factor=10.0, uncovered_factor=1.0, sel_tax=-log(0.1)))
        @test nvars(prob_hi_sXs) == 5 # d is out

        @test aggscore([0.0, 0.0, 0.0, 1.0, 0.0], prob_hi_sXs) < aggscore([0.0, 0.0, 0.0, 0.0, 0.0], prob_hi_sXs)
        @test aggscore([0.0, 0.0, 0.0, 1.0, 0.0], prob_hi_sXs) < aggscore([0.0, 1.0, 0.0, 0.0, 0.0], prob_hi_sXs)
        res_ignore_overlap = optimize(prob_hi_sXs, MultiobjOptimizerParams(ϵ=[0.001, 0.001]))
        @test aggscore([0.0, 0.0, 0.0, 1.0, 0.0], prob_hi_sXs) < aggscore([1.0, 0.0, 0.0, 0.0, 0.0], prob_hi_sXs)
        @test best_varweights(res_ignore_overlap) ≈ [0.0, 0.0, 0.0, 1.0, 0.0]

        # higher prior probability to select sets, lower probability to miss active element, so select a and b
        problem_low_penalty = MultiobjCoverProblem(sm_abc, CoverParams(setXset_factor=1.0, uncovered_factor=1.0, sel_tax=-log(0.9)))
        @test aggscore([1.0, 0.0, 1.0, 0.0, 0.0], problem_low_penalty) < aggscore([1.0, 1.0, 0.0, 0.0, 0.0], problem_low_penalty)
        @test aggscore([1.0, 0.0, 1.0, 0.0, 0.0], problem_low_penalty) < aggscore([0.0, 0.0, 0.0, 1.0, 0.0], problem_low_penalty)
        res_low_penalty = optimize(problem_low_penalty, MultiobjOptimizerParams(ϵ=[0.001, 0.001]))
        @test best_varweights(res_low_penalty) ≈ [1.0, 0.0, 1.0, 0.0, 0.0]
    end

    @testset "multimask: [a b d] [b c d] [c] [d] [a b c d e] [c d e], mask=[[a b c] [b e]]" begin # FIXME take weights into account
        sm = SetMosaic([Set([:a, :b, :d]), Set([:b, :c, :d]), Set([:c]), Set([:d]),
                        Set([:a, :b, :c, :d, :e]), Set([:c, :d, :e, :f])],
                        Set([:a, :b, :c, :d, :e, :f]))
        sm_abc = mask(sm, [Set([:a, :b, :c]), Set([:b, :e])], min_nmasked=1)

        # lower prior probability to select sets, high overlap penalty, so select abd and c FIXME update desc
        prob_hi_sXs = MultiobjCoverProblem(sm_abc, CoverParams(setXset_factor=1.0, uncovered_factor=1.0, covered_factor=0.5, sel_tax=-log(0.6)))
        @test_skip nmasks(prob_hi_sXs) == 2
        @test nvars(prob_hi_sXs) == 5 # d is out: abd bcd c abcde cdef

        @test OESC.miscover_score([0.0, 0.0, 0.0, 0.0, 0.0], prob_hi_sXs) == (5.0, 0.0)
        @test OESC.miscover_score([1.0, 0.0, 0.0, 0.0, 0.0], prob_hi_sXs) == (2.0, 1.0 + 0.9 * 2.0)
        @test OESC.miscover_score([0.0, 1.0, 0.0, 0.0, 0.0], prob_hi_sXs) == (2.0, 1.0 + 0.9 * 2.0)
        @test OESC.miscover_score([0.0, 0.0, 1.0, 0.0, 0.0], prob_hi_sXs) == (4.0, 0.0 + 0.9 * 1.0)

        @test aggscore([1.0, 0.0, 1.0, 0.0, 0.0], prob_hi_sXs) < aggscore([0.0, 0.0, 0.0, 0.0, 0.0], prob_hi_sXs)
        @test aggscore([0.0, 0.0, 0.0, 1.0, 0.0], prob_hi_sXs) < aggscore([0.0, 0.0, 0.0, 0.0, 0.0], prob_hi_sXs)
        @test aggscore([1.0, 0.0, 1.0, 0.0, 0.0], prob_hi_sXs) > aggscore([0.0, 0.0, 0.0, 1.0, 0.0], prob_hi_sXs)
        res_hi_sXs = optimize(prob_hi_sXs, MultiobjOptimizerParams(ϵ=[0.001, 0.001]))
        @test best_varweights(res_hi_sXs) ≈ [0.0, 0.0, 0.0, 1.0, 0.0] atol=1E-2

        # higher prior probability to select sets, no overlap penalty, so select abd, bcd, c and abcde + cdef
        prob_low_sXs = MultiobjCoverProblem(sm_abc, CoverParams(setXset_factor=0.01, covered_factor=1E-5, uncovered_factor=1.0, sel_tax=-log(0.9)))
        @test_skip nmasks(prob_low_sXs) == 2
        @test aggscore([1.0, 0.0, 1.0, 0.0, 0.0], prob_low_sXs) < aggscore([0.0, 0.0, 0.0, 0.0, 0.0], prob_low_sXs)
        @test aggscore([1.0, 0.0, 1.0, 1.0, 0.0], prob_low_sXs) < aggscore([1.0, 0.0, 1.0, 0.0, 0.0], prob_low_sXs)
        @test aggscore([1.0, 0.0, 1.0, 1.0, 0.0], prob_low_sXs) < aggscore([0.0, 0.0, 1.0, 1.0, 0.0], prob_low_sXs)
        @test aggscore([1.0, 1.0, 1.0, 1.0, 0.0], prob_low_sXs) < aggscore([0.0, 0.0, 1.0, 1.0, 0.0], prob_low_sXs)
        @test aggscore([1.0, 1.0, 1.0, 1.0, 0.0], prob_low_sXs) < aggscore([1.0, 0.0, 1.0, 1.0, 0.0], prob_low_sXs)
        res_low_sXs = optimize(prob_low_sXs, MultiobjOptimizerParams(ϵ=[0.001, 0.001]))
        @test best_varweights(res_low_sXs) ≈ [1.0, 1.0, 1.0, 1.0, 0.0] atol=1E-2# [1.0, 0.0, 1.0, 1.0, 0.0] atol=1E-2
    end

    @testset "WeightedSetMosaic-based problem" begin
        @testset "empty" begin
            sm = SetMosaic(Set{Symbol}[])
            wsm = assignweights(sm, [Dict{Int, Float64}()])
            problem = MultiobjCoverProblem(wsm)
            @test nvars(problem) == 0
            @test_skip nexperiments(problem) == 1
            @test_throws DimensionMismatch score([1.0], problem)
            @test score(Float64[], problem) == (0.0, 0.0, 0.0, 0.0)
            @test aggscore(Float64[], problem) == 0.0

            res = optimize(problem)
            @test isempty(res)
            @test nvars(res) == 0
            @test nsolutions(res) == 0
            @test best_index(res) == 0
        end

        @testset "[a]" begin # FIXME take element detection probability into account
            # empty mask problem
            empty_prob = MultiobjCoverProblem(assignweights(SetMosaic([Set([:a])]), [Dict{Int, Float64}()]))
            @test nvars(empty_prob) == 0
            @test score(Float64[], empty_prob) == (0.0, 0.0, 0.0, 0.0)
            @test aggscore(Float64[], empty_prob) == 0.0
            empty_res = optimize(empty_prob)
            @test isempty(empty_res)
            @test nvars(empty_res) == 0
            @test nsolutions(empty_res) == 0

            # enabled because sel prob is high, and although :a is all elements it needs to be covered
            prob1_mosaic = assignweights(SetMosaic([Set([:a])], Set([:a])), [Dict(1 => -1.0)])
            en1_prob = MultiobjCoverProblem(prob1_mosaic, CoverParams(sel_tax=-log(0.5), uncovered_factor=1.0))
            @test nvars(en1_prob) == 1
            @test_throws DimensionMismatch score(Float64[], en1_prob)
            en1_score = score([1.0], en1_prob)
            @test en1_score == (en1_prob.var_scores[1], 0.0, 0.0, 0.0)
            @test score([0.0], en1_prob) == (0.0, 0.0, 1.0, 0.0)
            @test en1_score[1] < 1.0 # score[3] at [0.0], required for enabling
            @test score([0.5], en1_prob) == (0.5*en1_prob.var_scores[1], 0.0, 0.5, 0.0)
            @test aggscore([1.0], en1_prob) < aggscore([0.9], en1_prob) <
                aggscore([0.5], en1_prob) < aggscore([0.0], en1_prob)
            en1_res = optimize(en1_prob, MultiobjOptimizerParams(ϵ=0.001))
            @test nsolutions(en1_res) > 0
            @test nvars(en1_res) == 1
            @test best_varweights(en1_res) == ones(Float64, 1)

            # disabled because the weight is zero
            dis1_prob = MultiobjCoverProblem(prob1_mosaic, CoverParams(sel_tax=-log(0.01), uncovered_factor=1.0))
            @test nvars(dis1_prob) == 1
            @test dis1_prob.var_scores[1] > en1_prob.var_scores[1]
            dis1_score = score([1.0], dis1_prob)
            @test dis1_score == (dis1_prob.var_scores[1], 0.0, 0.0, 0.0)
            @test score([0.0], dis1_prob) == (0.0, 0.0, 1.0, 0.0)
            @test dis1_score[1] > 1.0 # score[3] at [0.0], required for disabling
            @test aggscore([1.0], dis1_prob) > aggscore([0.5], dis1_prob) > aggscore([0.0], dis1_prob)
            dis1_res = optimize(dis1_prob)
            @test nsolutions(dis1_res) > 0
            @test nvars(dis1_res) == 1
            @test best_varweights(dis1_res) == zeros(Float64, 1)

            # enabled because the weight is negative
            prob2_mosaic = assignweights(SetMosaic([Set([:a])], Set([:a])), [Dict(1 =>  -1)])
            en2_prob = MultiobjCoverProblem(prob2_mosaic, CoverParams(sel_tax=-log(0.9)))
            @test nvars(en2_prob) == 1
            en2_res = optimize(en2_prob)
            @test nsolutions(en2_res) > 0
            @test nvars(en2_res) == 1
            @test best_varweights(en2_res) == ones(Float64, 1)
        end

        @testset "[a b] [c d] [a b c d], mask=[a b]" begin # FIXME take weights into account
            sm = SetMosaic([Set([:a, :b]), Set([:c, :d]), Set([:a, :b, :c])])
            sm_ab = assignweights(sm, [Dict(1 => -5, 3 => -3)])

            problem_def = MultiobjCoverProblem(sm_ab) # cd set is dropped

            @test OESC.var_scores([0.0, 0.0], problem_def) == (0.0, 0.0)
            @test OESC.var_scores([1.0, 0.0], problem_def)[2] == 0.0
            @test OESC.var_scores([0.0, 1.0], problem_def)[2] == 0.0
            @test OESC.var_scores([1.0, 1.0], problem_def)[2] > 0

            res_def = optimize(problem_def)
            @test nsolutions(res_def) > 0
            @test best_varweights(res_def) ≈ [1.0, 0.0] atol=0.01

            problem_no_penalty = MultiobjCoverProblem(sm_ab, CoverParams(setXset_factor=0.0, covered_factor=0.0, sel_tax=0.0))
            @test nvars(problem_no_penalty) == 2
            res_no_penalty = optimize(problem_no_penalty, MultiobjOptimizerParams(ϵ=[0.001, 0.001], WeightDigits=nothing))
            @test nvars(res_no_penalty) == 2
            @test best_varweights(res_no_penalty) ≈ [1.0, 1.0] atol=0.01
        end
    end
end


text:
using DataFrames
using Rmath
using Distributions
using GLM
using StatsBase
using Winston
using DataArrays

module RobustStats
using DataFrames
using Rmath
using Winston
using Distributions
using StatsBase
using GLM
using DataArrays
import Base.show
import DataFrames.complete_cases
import DataFrames.deleterows!

export
    outOutput,
    idealfOutput,
    testOutput,
    tmean,
    winval,
    winmean,
    winvar,
    trimse,
    trimci,
    stein1,
    stein2,
    idealf,
    pbvar,
    bivar,
    tauloc,
    tauvar,
    outbox,
    akerd,
    akerd_C,
    rdplot,
    sint,
    binomci,
    acbinomci,
    near,
    msmedse,
    stein1_tr,
    stein2_tr,
    sintv2,
    seq,
    cnorm,
    hpsi,
    onestep,
    onesampb,
    mom,
    momci,
    trimpb,
    trimcibt,
    indt,
    pcorb,
    bootse,
    indirectTest,
    yuend,
    t1way

include("types.jl")
include("utilis.jl")
include("show.jl")
include("functions.jl")
include("data.jl")
end


text:
using SafeTestsets

@safetestset "Analytical 2x2 example test" begin include("twobytwo_tests.jl") end
@safetestset "QL/CE2 and NL/GQL(M)/GCE2(M) tests" begin include("qlnl_tests.jl") end
@safetestset "GQL(Λ)/GCE2(Λ) tests" begin include("gql_tests.jl") end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

"""
TODO
"""
module App

## INCLUDES
include("libs.jl")
include("shader.jl")
include("compileAndLink.jl")

## COMPILE C File
const compileAndLink = isdefined(:USE_COMPILE)
if compileAndLink compileWithGCC() end

WIREFRAME = false
TEXTUREMODE = true
LIGHTMODE = true
FRUSTUM_CULLING = true
HIDE_UNSEEN_CUBES = true
RENDER_METHOD = 1
mychunk = nothing
program_chunks = 0
SCENE = 2
fstm = nothing
planeData = nothing
boxData = nothing
CHUNK_SIZE = 64
current_program = 0
load_once = true

"""
TODO
"""
function use_program(program, f)
  if program != current_program
    glUseProgram(program)
    glCheckError("glUseProgram")
  end
  f()
  if program != current_program
    glUseProgram(current_program)
    glCheckError("glUseProgram")
  end
end

"""
sets a mode in a shader.
"""
function setMode(program, name, value, mode="")
  use_program(program, () -> begin
    l = glGetUniformLocation(program, name)
    if l>-1
      if isa(value, Integer) glUniform1i(l, value)
      elseif isa(value, AbstractFloat) glUniform1f(l, value)
      elseif isa(value, AbstractArray) glUniform1fv(l, 1, value)
      else warn("MODE("*displayInYellow(mode)*") for "*typeof(value)*" is not implemented yet.")
      end
    end
  if mode != "" println("MODE(",displayInYellow(mode),"): ",displayInYellow(value)) end
  end)
end

"""
TODO
"""
function setFrustumCulling(load=true)
  global load_once
  if load || load_once
    SetCamera(fstm, Vec3f(CAMERA.position), Vec3f(CAMERA.position+forward(CAMERA)), Vec3f(0,1,0))
    load_once = false
  end

  if FRUSTUM_CULLING checkInFrustum(mychunk, fstm)
  else showAll(mychunk)
  end

  if HIDE_UNSEEN_CUBES hideUnseen(mychunk) end

  if load
    update(mychunk)
    global BLOCK_COUNT = mychunk.fileredCount
    upload(chunkData, :instances, getData(mychunk))
    if FRUSTUM_CULLING
      upload(planeData, :vertices, getVertices(fstm))
      upload(boxData, :vertices, getBox(fstm))
    end
  end
end

"""
TODO
"""
function chooseRenderMethod(method=RENDER_METHOD)

  if method == 1 name="INSTANCES of POINTS + GEOMETRY SHADER"
  elseif method == 2  name="INSTANCES of TRIANGLES"
  elseif method == 3 name="INSTANCES of TRIANGLES + INDICIES"
  elseif method == 4 name="TRIANGLES + INDICIES"
  elseif method == 5 name="TRIANGLES"
  elseif method == 6 name="TRIANGLES"
  else name="NOT DEFINED"
  end

  global RENDER_METHOD = method

  println("METHOD(",method,"): ",displayInYellow(name))

  if name == "NOT DEFINED" return end

  info("Create chunk & create blocks...")

  ## BLOCKS
  global mychunk
  clean(mychunk) # free memory

  mychunk = Chunk(CHUNK_SIZE)

  if SCENE == 0 createSingle(mychunk)
  elseif SCENE == 1 createExample(mychunk)
  else createLandscape(mychunk)
  end

  setFrustumCulling(false)

  update(mychunk)

  global BLOCK_COUNT = mychunk.fileredCount
  #global VERT_COUNT = mychunk.verticesCount

  info("Upload data...")
  glCheckError("Upload data")

  if method == 1 linkData(chunkData, :instances=>getData(mychunk))
  elseif method == 2 linkData(chunkData, :vertices=>(DATA_CUBE,3), :instances=>getData(mychunk))
  elseif method == 3 linkData(chunkData, :vertices=>(DATA_CUBE_VERTEX,3), :indicies=>(DATA_CUBE_INDEX,1,GL_ELEMENT_ARRAY_BUFFER, true), :instances=>getData(mychunk))
  elseif method == 4 linkData(chunkData, :vertices=>(DATA_CUBE_VERTEX,3), :indicies=>(DATA_CUBE_INDEX,1,GL_ELEMENT_ARRAY_BUFFER, true))
  elseif method > 4 linkData(chunkData, :vertices=>(DATA_CUBE,3))
  end

  global program_chunks

  if program_chunks != 0
    println("Unbind & Delete previous program")
    glUseProgram(0)
    glDeleteProgram(program_chunks)
    glCheckError("glDeleteProgram")
  end

  (INST_VSH, INST_VSH_GSH, INST_GSH, INST_FSH, VSH_TEXTURE, VSH, FSH, GSH) = loadShaders()

  if method == 1 program_chunks = createShaderProgram(INST_VSH_GSH, INST_FSH, INST_GSH)
  elseif method > 1 && method <= 3 program_chunks = createShaderProgram(INST_VSH, INST_FSH)
  elseif method == 4 program_chunks = createShaderProgram(VSH_TEXTURE, INST_FSH)
  elseif method == 5 program_chunks = createShaderProgram(VSH_TEXTURE, INST_FSH)
  elseif method > 5 program_chunks = createShaderProgram(VSH, FSH)
  end

  info("Compile Shader Programs...")

  if program_chunks <= 0 error("No Program") end
  glUseProgram(program_chunks)
  glCheckError("glUseProgram")

  setAttributes(chunkData, program_chunks)
  setMatrix(program_chunks, "iMVP", CAMERA.MVP)
  glCheckError("setMatrix")

  global location_position = glGetUniformLocation(program_chunks, "iPosition")
  global location_texindex = glGetUniformLocation(program_chunks, "iTexIndex")

  setMode(program_chunks, "iUseLight", LIGHTMODE)
  setMode(program_chunks, "iUseTexture", TEXTUREMODE)
end

"""
TODO
"""
function checkForUpdate()
  global keyPressed, keyValue, cam_updated, FRUSTUM_KEY, ALL_KEY, CAMERA, mychunk, chunkData, fstm, planeData

  if keyPressed

    if keyValue == 80 #p
      setPosition(CAMERA,[0f0,0,0])

    elseif keyValue == 81 #q
      global WIREFRAME=!WIREFRAME

    elseif keyValue == 84 #t
      global TEXTUREMODE=!TEXTUREMODE
      setMode(program_chunks, "iUseTexture", TEXTUREMODE, "TEXTURE")

    elseif keyValue == 76 #l
      global LIGHTMODE=!LIGHTMODE
      setMode(program_chunks, "iUseLight", LIGHTMODE, "LIGHT")

    elseif keyValue == 82 #r
      chooseRenderMethod()

    elseif (keyValue >= 290 && keyValue <= 301) # F1 - F12
      chooseRenderMethod(keyValue - 289)

    elseif keyValue == 49 #1
      global CHUNK_SIZE = 1
      chooseRenderMethod()

    elseif keyValue == 50 #2
      global CHUNK_SIZE = 4
      chooseRenderMethod()

    elseif keyValue == 51 #3
      global CHUNK_SIZE = 8
      chooseRenderMethod()

    elseif keyValue == 52 #4
      global CHUNK_SIZE = 16
      chooseRenderMethod()

    elseif keyValue == 53 #5
      global CHUNK_SIZE = 24
      chooseRenderMethod()

    elseif keyValue == 54 #6
      global CHUNK_SIZE = 32
      chooseRenderMethod()

    elseif keyValue == 55 #7
      global CHUNK_SIZE = 40
      chooseRenderMethod()

    elseif keyValue == 56 #8
      global CHUNK_SIZE = 48
      chooseRenderMethod()

    elseif keyValue == 57 #9
      global CHUNK_SIZE = 56
      chooseRenderMethod()

    elseif keyValue == 48 #0
      global CHUNK_SIZE = 64
      chooseRenderMethod()

    elseif keyValue == 45 #ß
      global CHUNK_SIZE = 72
      chooseRenderMethod()

    elseif keyValue == 61 #´
      global CHUNK_SIZE = 96
      chooseRenderMethod()

    elseif keyValue == 96 #^
      global CHUNK_SIZE = 128
      chooseRenderMethod()

    elseif keyValue == 66 #b
      global SCENE = 0
      chooseRenderMethod()
    elseif keyValue == 78 #n
      global SCENE = 1
      chooseRenderMethod()
    elseif keyValue == 77 #m
      global SCENE = 2
      chooseRenderMethod()

    elseif keyValue == 70 #f
      global FRUSTUM_CULLING = !FRUSTUM_CULLING
      setFrustumCulling()
    elseif keyValue == 86 #v
      setFrustumCulling()
    elseif keyValue == 79 #o
      global HIDE_UNSEEN_CUBES = !HIDE_UNSEEN_CUBES
      setFrustumCulling()
    end
  end

  if keyPressed keyPressed=false end
end

"""
TODO
"""
useProgram(program) = begin global current_program = program; glUseProgram(program) end

"""
TODO
"""
setMatrix(program, name, m) = begin const cm = SMatrix{4,4,Float32}(m); glUniformMatrix4fv(glGetUniformLocation(program, name), 1, false, cm) end

"""
TODO
"""
function setMVP(program, mvp)
  use_program(program, () -> begin
    setMatrix(program, "iMVP", mvp)
    glCheckError("setMatrix")
  end)
end

#setMVP(mvp) = setMatrix(program, "iMVP", mvp)

## PROGRAM

"""
TODO
"""
function run()

println("---------------------------------------------------------------------")
println("Start Program @ ", Dates.time())
versioninfo()

# OS X-specific GLFW hints to initialize the correct version of OpenGL
GLFW.Init()

# Create a windowed mode window and its OpenGL context
global window = GLFW.CreateWindow(WIDTH, HEIGHT, "OpenGL Example")

# Make the window's context current
GLFW.MakeContextCurrent(window)

GLFW.SetWindowSize(window, WIDTH, HEIGHT) # Seems to be necessary to guarantee that window > 0
rezizeWindow(WIDTH,HEIGHT)

# Window settings
GLFW.SwapInterval(0) # intervall between canvas images (min. 2 images)

# Graphcis Settings
GLFW.WindowHint(GLFW.OPENGL_DEBUG_CONTEXT,true) #show debug
#GLFW.WindowHint(GLFW.SAMPLES,4)

# OpenGL Version
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR,4)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR,6)

GLFW.SetCursorPosCallback(window, OnCursorPos)
GLFW.SetKeyCallback(window, OnKey)
GLFW.SetMouseButtonCallback(window, OnMouseKey)

#setEventCallbacks(OnCursorPos,OnKey,OnMouseKey)

GLFW.ShowWindow(window)

glinfo = createcontextinfo()

println("OpenGL $(displayInRed(glinfo[:gl_version]))")
println("GLSL $(glinfo[:glsl_version])")
println("Vendor $(displayInRed(glinfo[:gl_vendor]))")
println("Renderer $(displayInRed(glinfo[:gl_renderer]))")
println("---------------------------------------------------------------------")

## CAMERA

setPosition(CAMERA,[0f0,0,0])
setProjection(CAMERA, projection_perspective(FOV, RATIO, CLIP_NEAR, CLIP_FAR))

global fstm = Frustum()
SetFrustum(fstm, FOV, RATIO, CLIP_NEAR, 1000f0)
SetCamera(fstm, Vec3f(CAMERA.position), Vec3f(CAMERA.position+forward(CAMERA)), Vec3f(0,1,0))

Update(CAMERA)

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

program = 0

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

global chunkData = MeshData()
global planeData = MeshData()
global boxData = MeshData()

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

## TEXTURES

uploadTexture("blocks.png")

## LOAD DEFAULT
chooseRenderMethod()

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

#linkData(chunkData, :vertices=>(DATA_PLANE_VERTEX,3), :indicies=>(DATA_PLANE_INDEX,1,GL_ELEMENT_ARRAY_BUFFER, true), :instances=>getData(mychunk))
linkData(planeData,  :vertices=>getVertices(fstm))
linkData(boxData,  :vertices=>getBox(fstm))

#chunkData.arrays[:vertices].count
#n = length(cubeVertices_small) / 3

#function compileShaderPrograms()
#  global program_chunks, program_normal

global program_normal = createShaderProgram(VSH, FSH) #, createShader(GSH, GL_GEOMETRY_SHADER)

setAttributes(planeData, program_normal)
setAttributes(boxData, program_normal)

setMVP(program_normal, CAMERA.MVP)

#end

#compileShaderPrograms()

global location_position = -1
global location_texindex = -1

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

function updateBlocks()
  #const mmvp = SMatrix{4,4,Float32}(CAMERA.MVP)
  #setMVP(CAMERA.MVP)
  #glUniform3fv(location_shift, 1, CAMERA.position)
  #glUniform3fv(location_shift, 1, shiftposition)
  #for b in blocks; b.mvp=mmvp*b.model; end
end

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

glEnable(GL_DEPTH_TEST)
glEnable(GL_BLEND)
glEnable(GL_CULL_FACE)
#glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
#glBlendEquation(GL_FUNC_ADD)
#glFrontFace(GL_CCW)
glCullFace(GL_BACK)
#glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)#GL_FILL,GL_LINE
glClearColor(0.0, 0.0, 0.0, 1.0)

# Loop until the user closes the window
render = function(x)
  #mvp = mmvp*MMatrix{4,4,Float32}(unsafe_wrap(Array, mvp, (4,4)))
  #setMVP(CAMERA.MVP*translation([c.x,c.y,c.z]))

  #glUniformMatrix4fv(location_mvp, 1, false, x.mvp)
  #glDrawElements(GL_TRIANGLES, chunkData.draw.count, GL_UNSIGNED_INT, C_NULL )
  #glDrawArrays(GL_TRIANGLES, 0, chunkData.draw.count)
  nothing
end

#=
if use_geometry_shader
  const loopBlocks() = render(mychunk.childs[1])
else
  if compileAndLink
    objptr = createLoop(1,refblocks,render) #compileAndLink
    const loopBlocks() = loopByObject(objptr) #compileAndLink
  else
    const loopBlocks() = for b in mychunk.childs; render(b); end
  end
end
=#

cam_updated=false

const SLEEP=0 #1f0/200

i=0

"""
TODO
"""
while !GLFW.WindowShouldClose(window)
  showFrames()
  UpdateCounters()

  if OnUpdate(CAMERA)
    setMVP(program_chunks, CAMERA.MVP)
    setMVP(program_normal, CAMERA.MVP)
    cam_updated=true
  end

  checkForUpdate()
  if cam_updated cam_updated=false end

  # Pulse the background
  #c=0.5 * (1 + sin(i * 0.01)); i+=1
  #glClearColor(c, c, c, 1.0)
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

  #print("loopBlocks "); @time
  #loopBlocks()

  glPolygonMode(GL_FRONT_AND_BACK, WIREFRAME ? GL_LINE : GL_FILL)

  """
  TODO
  """
  if isValid(mychunk)
    useProgram(program_chunks)
    #glCheckError("useProgram")
    #glPolygonMode(GL_FRONT_AND_BACK, WIREFRAME ? GL_LINE : GL_FILL)
    glBindVertexArray(chunkData.vao)
    #glCheckError("glBindVertexArray bind")

    if RENDER_METHOD == 1 glDrawArraysInstanced(GL_POINTS, 0, 1, mychunk.fileredCount) #GL_TRIANGLE_STRIP
    elseif RENDER_METHOD == 2 glDrawArraysInstanced(GL_TRIANGLES, 0, chunkData.draw.count, mychunk.fileredCount)
    elseif RENDER_METHOD == 3 glDrawElementsInstanced(GL_TRIANGLES, chunkData.draw.count, GL_UNSIGNED_INT, C_NULL, mychunk.fileredCount)
    #glDrawElementsInstancedBaseVertex(GL_TRIANGLES, chunkData.draw.count / 6, GL_UNSIGNED_INT, C_NULL, mychunk.count, 0)
    elseif RENDER_METHOD == 4
      #* thats slow! (glDrawElements ~60 fps, glDrawElementsInstanced ~ 200 fps !!!)
      for b in getFilteredChilds(mychunk)
        if location_texindex > -1 glUniform1f(location_texindex, b.typ) end
        if location_position > -1 glUniform3fv(location_position, 1, b.pos) end
        glDrawElements(GL_TRIANGLES, chunkData.draw.count, GL_UNSIGNED_INT, C_NULL )
        #glCheckError("glDrawElements")
      end
    elseif RENDER_METHOD == 5
      glDrawArrays(GL_TRIANGLES, 0, chunkData.draw.count)
    end
    glBindVertexArray(0)
    #glCheckError("glBindVertexArray unbind")
  end

  #################################

  useProgram(program_normal)
  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)

  #setMode(program_normal, "color", Vec3f(1,0,0))

  glBindVertexArray(boxData.vao)
  glDrawArrays(GL_TRIANGLES, 0, boxData.draw.count)
  glBindVertexArray(0)

  #setMode(program_normal, "color", Vec3f(1,1,0))

  glBindVertexArray(planeData.vao)
  glDrawArrays(GL_TRIANGLES, 0, planeData.draw.count)
  glBindVertexArray(0)

  #################################

  # Swap front and back buffers
  GLFW.SwapBuffers(window)
  # Poll for and process events
  GLFW.PollEvents()

  if SLEEP>0 Libc.systemsleep(SLEEP) end
end

GLFW.DestroyWindow(window)
GLFW.Terminate()

end

end

"""
TODO
"""
function main()
  App.run()
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using Printf
@show VERSION

function calc_julia()
    total = 0

    for i in 1:10^5
        for j in 1:10^4
            tmp = if j % 4 == 0
                     i + j
                  elseif j % 4 == 1
                     i - j
                  elseif j % 4 == 2
                     i * j
                  else
                    i ÷ j
                  end
            total += tmp
        end
    end

    return total
end

println()
@time total = calc_julia()
@time total = calc_julia()
@time total = calc_julia()

@printf("total = %d\n", total)


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

using Actors

import EasyPkg
@EasyPkg.using_BaseTest

@testset "Basic Actor operations" begin
    echo_actor = @actor begin
        while true
            const sender, msg = receive()
            tell(sender, msg)
        end
    end

    try
        @testset "tell() and receive()" begin
            tell(echo_actor, "foo")
            tell(echo_actor, "bar")

            @test receive() == (echo_actor => "foo")
            @test receive() == (echo_actor => "bar")
        end

        @testset "ask()" begin
            @test ask(echo_actor, "baz") == "baz"
        end

    finally
        if VERSION.minor <= 4
            # Currently doesn't work right for Julia >= v0.5:
            @testset "kill()" begin
                try
                    kill(echo_actor)
                catch exc
                    @test exc == InterruptException()
                end
            end
        end
    end
end


text:
using HybridSystems
include(joinpath(dirname(dirname(pathof(HybridSystems))), "examples", "cruise_control.jl"));

if true
    D = 1.0
    U = 1.0
    v_shift = 2.0
    vmin = -1.0
    vmax = 2.0
    v = (1.0,)
    m = 1.0
    m0 = 1.0
    h = 0.5
    kd = 1/2
    ks = 1/2
    Δv = 5.0
else
    function constant(scaling)
        va = 15.6 * scaling
        vb = 24.5 * scaling
        vc = 29.5 * scaling
        v = (va, vb, vc)
        U = 4.0 * scaling
        D = 0.5 * scaling
        vmin = 5.0 * scaling
        return D, U, vmin, v
    end
    D, U, vmin, v = constant(1.0)
    vmax = 35.0
    m = 1000
    m0 = 500
    h = 0.8
    kd = 4600
    ks = 4500
    Δv = (v[1] - vmin) / 2
end

T = 2
N = 1 + (T+1) * length(v)
function system(M, T; v=v, N = 1 + (T+1) * length(v), sym=false, vmax = vmax)
    H = h * T
    #return cruise_control_example(N, M, vmin = vmin, vmax=vmax, v=v, U=U, H=H, D=D, T=T, sym=sym, m0 = m0)
    return cruise_control_example(N, M, vmin = vmin, vmax=vmax, v=v,
        U=U, H=h*T, D=D, T=T, sym=sym, m0 = m0, ks=ks, m=m, kd=kd)
end
_vec(M, d, v, u) = [repeat([d, v], M); v; u]
# |string_elongation| < D
# |acceleration| < U
sym_rect(M) = _vec(M, D, 1.0, U) # _vec(M, D, Δv, U)

using SetProg
using SwitchOnSafety
const SOS = SwitchOnSafety;

using MathOptInterface
const MOI = MathOptInterface

function symsolve(M, set_variable; volume_heuristic=set -> L1_heuristic(set, sym_rect(M)))
    hs = system(M, T, N=1, sym=true)
    # Shift interval 5 <= v <= 15.6 -> -5.3 <= v <= 5.3
    sets = invariant_sets(hs, sdp_factory, [set_variable], volume_heuristic=volume_heuristic)
    [SOS.SetProg.Sets.Translation(set, _vec(M, 0, 2.0, 0)) for set in sets]
end

function fullsolve(T, M, setvar::Function;
                   volume_heuristic=set -> L1_heuristic(set, sym_rect(M)),
                   oneshot = false,
                   onlyone = false)
    hs = system(M, T)
    function hv(v)
        h = zeros(statedim(hs, 1))
        for i in 1:M
            h[2i] = (vmin .+ v) / 2
        end
        h[2M+1] = (vmin .+ v) / 2
        h
    end
    habc = SOS.SetProg.InteriorPoint.(hv.(v))
    ha = habc[1]
    @show nstates(hs)
    hi = fill(ha, nstates(hs))
    set_variables = map(setvar, hi)
    if oneshot
        λ = Dict(t => 1.0 for t in transitions(hs))
    else
        λ = Dict(t => 1.0 for t in transitions(hs) if source(hs, t) == target(hs, t))
    end
    if oneshot
        return invariant_sets(hs, sdp_factory, set_variables, λ = λ, volume_heuristic=volume_heuristic);
    else
        sets = Vector{SwitchOnSafety.SetProg.Sets.AbstractSet{Float64}}(undef, nstates(hs))
        for i in (onlyone ? (1:1) : (1:nstates(hs)))
            println("Computing set $i")
            invariant_sets!(sets, i:i, hs, sdp_factory, set_variables, λ = λ, enabled=1:i,
                            volume_heuristic=volume_heuristic)
        end
        return sets
    end
end


text:
@testset "py2jl mod gen" begin
    py2jl"""
def sumBy(f, seq):
    s = 0
    for each in seq:
        s = s + f(each)
    return s

result = sumBy(lambda x: 1 + x, [100, 200])
    """
    @testset "transpiled codes run in the pure Julia runtime" begin
        @test result == 302
        @test sumBy(identity, [1, 2, 3]) === 6
    end
end


text:
using Femshop2
using Test

@testset "Femshop2.jl" begin
    # Write your tests here.
end