text:
#cd("D:\\OneDrive\\Current\\meanVar\\Code\\hox_krox")
cd("/home/crackauc/.julia/RAMeanVar/Code/hox_krox")
include("bp_sharpening_experiment_function.jl")
cd("data")
cd("noise_noise_plot")
cd("cyp_knocks")
############################# Solo Run

# 100.0 is sharp, 1.0 is unsharp
Vbp = 100.0#1.0
kdeg = 3.3e4
ϵₕₖ=0.30
sol = boundary_sharpening(Vbp=Vbp,ϵₕₖ=ϵₕₖ)
using JLD
#save("wildtype.jld", "u", sol.u, "t", sol.t)
save("bp_$(Int(Vbp))_epshk_$(ϵₕₖ).jld", "u", sol.u, "t", sol.t)

############################### Multi Run

Vbp = 5.0
ϵₕₖ = 0.3
using JLD
for Vbp in Vbps
  for i in 1:10
    sol = boundary_sharpening(Vbp=Vbp,ϵₕₖ=ϵₕₖ)
    save("bp_$(Vbp)_$i.jld", "u", sol.u, "t", sol.t)
  end
end

kdeg=3.3e4
using JLD
for i in 1:10
  sol = boundary_sharpening(kdeg=kdeg)
  save("cyp_$(kdeg)_$i.jld", "u", sol.u, "t", sol.t)
end

Vbps = 5.0:10.0:105.0
ϵₕₖ = 0.3

for Vbp in Vbps
    sol = boundary_sharpening(Vbp=Vbp,ϵₕₖ=ϵₕₖ)
    save("bp_$(Int(Vbp))_epshk_$(ϵₕₖ).jld", "u", sol.u, "t", sol.t)
end

Vbps = 5.0:10.0:105.0
ϵₕₖs = 0.1:0.025:0.3
for Vbp in Vbps, ϵₕₖ in ϵₕₖs
  sol = boundary_sharpening(Vbp=Vbp,ϵₕₖ=ϵₕₖ)
  save("bp_$(Int(Vbp))_epshk_$(ϵₕₖ).jld", "u", sol.u, "t", sol.t)
end

max_stars_repo_name: ChrisRackauckas/RAMeanVar

max_issues_repo_name: ChrisRackauckas/RAMeanVar

max_forks_repo_name: ChrisRackauckas/RAMeanVar


text:
function cse(q::Query)
    gr = query2ev(q)
    q′ = ev2query(gr)
    q′
end

immutable EvNode
    q::Query
    hc::UInt
    self::Int
    pred::Int
    args::Vector{Int}
end

const NO_EV_ARGS = Int[]

immutable EvGraph
    q::Query
    nodes::Vector{EvNode}
    hc2idxs::Dict{UInt,Vector{Int}}
end

function query2ev(q::Query)
    gr = EvGraph(q, EvNode[], Dict{UInt64,Vector{Int}}())
    query2ev!(gr, q, 0)
    gr
end

function query2ev!(gr::EvGraph, q::Query, pred::Int)
    if isa(q.sig, ItSig) && domain(input(q)) == domain(output(q))
        return pred
    end
    if isa(q.sig, ComposeSig)
        self = pred
        for arg in q.args
            self = query2ev!(gr, arg, self)
        end
        return self
    end
    args =
        if !isempty(q.args)
            Int[query2ev!(gr, arg, pred) for arg in q.args]
        else
            NO_EV_ARGS
        end
    hc = hash((q.sig, pred, args))
    if hc in keys(gr.hc2idxs)
        for idx in gr.hc2idxs[hc]
            n = gr.nodes[idx]
            if n.q.sig == q.sig && n.pred == pred
                return idx
            end
        end
    else
        gr.hc2idxs[hc] = Int[]
    end
    self = length(gr.nodes) + 1
    n = EvNode(q, hc, self, pred, args)
    push!(gr.nodes, n)
    push!(gr.hc2idxs[hc], self)
    self
end

ev2query(gr::EvGraph) =
    ev2query(gr, Dict{Tuple{Int,Int},Query}(), length(gr.nodes), 0)

function ev2query(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, self::Int, tail::Int)
    if self == 0
        @assert tail == 0
        dom = domain(input(gr.q))
        return ItQuery(dom)
    end
    n = gr.nodes[self]
    if self == tail
        dom = domain(output(n.q))
        return ItQuery(dom)
    end
    if (self, tail) in keys(refs)
        return refs[(self, tail)]
    end
    q = ev2query(gr, refs, n, n.q.sig)
    if n.pred != tail
        prq = ev2query(gr, refs, n.pred, tail)
        q = prq >> q
    end
    q
end

ev2query(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, n::EvNode, ::AbstractPrimitive) =
    n.q

function ev2query(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, n::EvNode, sig::ConnectSig)
    args = Query[ev2query(gr, refs, arg, n.pred) for arg in n.args]
    return ConnectQuery(sig.refl, args...)
end

function ev2query(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, n::EvNode, ::FrameSig)
    args = Query[ev2query(gr, refs, arg, n.pred) for arg in n.args]
    return FrameQuery(args...)
end

function ev2query(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, n::EvNode, ::GivenSig)
    args = Query[ev2query(gr, refs, arg, n.pred) for arg in n.args]
    return GivenQuery(args...)
end

function ev2query(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, n::EvNode, ::RecordSig)
    dups = finddups(gr, refs, n.args, n.pred)
    if isempty(dups)
        args = Query[ev2query(gr, refs, arg, n.pred) for arg in n.args]
        return RecordQuery(args)
    end
    dupqs = Vector{Query}(length(dups))
    for i = endof(dups):-1:1
        tag = gensym()
        dup = dups[i]
        dupq = ev2query(gr, refs, dup, n.pred) |> decorate(:tag => tag)
        dupqs[i] = dupq
        ref = SlotQuery(tag, output(dupq))
        refs[(dup, n.pred)] = ref
    end
    args = Query[ev2query(gr, refs, arg, n.pred) for arg in n.args]
    q = RecordQuery(args)
    for i = 1:endof(dups)
        dup = dups[i]
        dupq = dupqs[i]
        delete!(refs, (dup, n.pred))
        q = GivenQuery(q, dupq)
    end
    q
end

function finddups(gr::EvGraph, refs::Dict{Tuple{Int,Int},Query}, fs::Vector{Int}, pred::Int)
    if length(fs) <= 1
        return NO_EV_ARGS
    end
    f2k = Dict{Int,Int}()
    stack = Int[]
    for (k, f) in enumerate(fs)
        push!(stack, f)
        while !isempty(stack)
            self = pop!(stack)
            if self == pred || ((self, pred) in keys(refs))
                continue
            end
            if !(self in keys(f2k))
                n = gr.nodes[self]
                if !cheap(n.q.sig)
                    f2k[self] = k
                end
                push!(stack, n.pred)
                if isa(n.q.sig, RecordSig) && n.pred == pred
                    for arg in n.args
                        push!(stack, arg)
                    end
                end
            elseif f2k[self] != k
                f2k[self] = 0
            end
        end
    end
    dups = Int[]
    for (self, k) in f2k
        if k == 0
            push!(dups, self)
        end
    end
    if length(dups) > 1
        dups1 = dups
        dups2 = finddups(gr, refs, dups, pred)
        dups2set = Set(dups2)
        dups = Int[]
        for dup in dups1
            if !(dup in dups2set)
                push!(dups, dup)
            end
        end
        for dup in dups2
            push!(dups, dup)
        end
    end
    dups
end

cheap(sig::AbstractSignature) =
    isa(sig, ConstSig) ||
    isa(sig, NullSig) ||
    isa(sig, ItSig)

ext: jl

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# abstract ScaLapackLite type
abstract type AbstractScaLapackLite end

# ScaLapackLite parameters
mutable struct ScaLapackLiteParams <: AbstractScaLapackLite
    mblocks::ScaInt; nblocks::ScaInt;
    mprocs::ScaInt; nprocs::ScaInt;
    rootproc::ScaInt
    # Constructor
    ScaLapackLiteParams() = new(1, 1, 1, 1, 0)
    ScaLapackLiteParams(mblocks, nblocks, mprocs, nprocs, rootproc) = new(mblocks, nblocks, mprocs, nprocs, rootproc)
end
function Base.:(==)(A::ScaLapackLiteParams, B::ScaLapackLiteParams)
    is_equal_blocks = (A.mblocks == B.mblocks) && (A.nblocks == B.nblocks)
    is_equal_procs = (A.mprocs == B.mprocs) && (A.nprocs == B.nprocs)
    is_equal_rootproc = (A.rootproc == B.rootproc)
    return is_equal_blocks && is_equal_procs && is_equal_rootproc
end

# ScalapackLiteVector type
mutable struct ScaLapackLiteVector <: AbstractScaLapackLite
    params::ScaLapackLiteParams
    x::Vector{T} where {T}
    ScaLapackLiteVector() = new(ScaLapackLiteVector())
    ScaLapackLiteVector(params) = new(params)
    v(params, x) = new(params, x)
end

# ScaLapackLiteMatrix type
# plannings:
#   - add member "localpart" that distributed matrix in each process
mutable struct ScaLapackLiteMatrix <: AbstractScaLapackLite
    params::ScaLapackLiteParams
    X::Matrix{T} where {T}
    ScaLapackLiteMatrix() = new(ScaLapackLiteParams())
    ScaLapackLiteMatrix(params) = new(params)
    ScaLapackLiteMatrix(params, X) = new(params, X)
end

#--- operator overloadings ---#

# perform A'
function Base.:(adjoint)(A::ScaLapackLiteMatrix)
    if length(A.X) > 0
        elty = typeof(A.X[1,1])
        if elty <: Complex
            return ScaLapackLiteMatrix(A.params, conj(transpose(A.X)))
        else
            return ScaLapackLiteMatrix(A.params, transpose(A.X))
        end
    else
        return ScaLapackLiteMatrix(A.params, transpose(A.X))
    end
end

# perform A*B
function Base.:(*)(A::ScaLapackLiteMatrix, B::ScaLapackLiteMatrix)
    # check A and B have same params
    if A.params == B.params
        C = multiple(A.X, B.X, A.params)
        return ScaLapackLiteMatrix(A.params, C)
    else
        throw(DimensionMismatch("ScaLapackLiteMatrix A and B have different ScaLapaclLiteParams"))
    end
end

#--- generate zero and identity matrix ---#
for (elty) in (:Float32, :Float64, :ComplexF32, :ComplexF64)
    @eval begin
        function zero!(A::Matrix{$elty}, m::ScaInt, n::ScaInt, desc::Vector{ScaInt})
            uplo = 'A'
            α = convert($elty, 0); β = convert($elty, 1);
            ia = one; ja = one;
            ScaLapack.pXlaset!(uplo, m, n,
                               α, β,
                               A, ia, ja, desc)
        end
    end
end
for (elty) in (:Float32, :Float64, :ComplexF32, :ComplexF64)
    @eval begin
        function identity!(I::Matrix{$elty}, m::ScaInt, n::ScaInt, desc::Vector{ScaInt})
            uplo = 'A'
            α = convert($elty, 0); β = convert($elty, 1);
            ia = one; ja = one;
            ScaLapack.pXlaset!(uplo, m, n,
                               α, β,
                               I, ia, ja, desc)
        end
    end
end

#--- multiple ---#
for (elty) in (:Float32, :Float64, :ComplexF32, :ComplexF64)
    @eval begin
        function multiple(A::Matrix{$elty}, B::Matrix{$elty}, params::ScaLapackLiteParams)

            # MPI parameters
            mblocks = params.mblocks
            nblocks = params.nblocks
            mprocs = params.mprocs
            nprocs = params.nprocs
            rootproc = params.rootproc
            comm = MPI.COMM_WORLD
            rank = MPI.Comm_rank(comm)

            # copy matrix from master to slave
            A = MPI.bcast(A, rootproc, comm)
            B = MPI.bcast(B, rootproc, comm)
            MPI.Barrier(comm)

            # check both dimensions; cols(A) = rows(B)
            # (mA, nA), (mB, nB) = (mA, nA), (nA, nB)
            if rank == rootproc
                check_multiple(A, B)
            end

            # type of the matrix A
            elty = typeof(A[1,1])

            # create context
            ctxt = ScaLapack.sl_init(mprocs, nprocs)
            ctxt0 = ScaLapack.sl_init(mprocs, nprocs)

            # matrix parameters
            mA, nA = size(A)
            mB, nB = size(B)
            mproc, nproc, myrow, mycol = BLACS.gridinfo(ctxt)
            mxlocr_A = ScaLapack.numroc(mA, mblocks, mproc, zero, mprocs)
            mxlocr_B = ScaLapack.numroc(mB, mblocks, mproc, zero, mprocs)
            mxlocc_A = ScaLapack.numroc(nA, nblocks, nproc, zero, nprocs)
            mxlocc_B = ScaLapack.numroc(nB, nblocks, nproc, zero, nprocs)
            mxllda_A = max(one, mxlocr_A)
            mxllda_B = max(one, mxlocr_B)

            if mproc >= zero && nproc >= zero

                # get array descriptor
                # (mA, nA) x (mB, nB) = (mA, nA) x (nA, nB) = (mA, nB)
                # (mA, nB) --> (mxlocr_A, mxlocc_B) --> mxllda_A
                desc_A = ScaLapack.descinit(mA, nA, mA, nA, zero, zero, ctxt0, mA)
                desc_B = ScaLapack.descinit(mB, nB, mB, nB, zero, zero, ctxt0, mB)
                desc_C = ScaLapack.descinit(mA, nB, mA, nB, zero, zero, ctxt0, mA)
                desc_myA = ScaLapack.descinit(mA, nA, mblocks, nblocks, zero, zero, ctxt, mxllda_A)
                desc_myB = ScaLapack.descinit(mB, nB, mblocks, nblocks, zero, zero, ctxt, mxllda_B)
                desc_myC = ScaLapack.descinit(mA, nB, mblocks, nblocks, zero, zero, ctxt, mxllda_A)

                # generate process matrix
                myA = zeros($elty, mxlocr_A, mxlocc_A)
                myB = zeros($elty, mxlocr_B, mxlocc_B)
                ScaLapack.pXgemr2d!(mA, nA,
                                    A, one, one, desc_A,
                                    myA, one, one, desc_myA,
                                    ctxt)
                ScaLapack.pXgemr2d!(mB, nB,
                                    B, one, one, desc_B,
                                    myB, one, one, desc_myB,
                                    ctxt)

                # perform A*B
                α = convert(elty, one); β = convert(elty, zero);
                myC = zeros($elty, mxlocr_A, mxlocc_B)
                transa = 'n'; transb = 'n';
                ScaLapack.pXgemm!(transa, transb,
                                  mA, nB, nA,
                                  α,
                                  myA, one, one, desc_myA,
                                  myB, one, one, desc_myB,
                                  β,
                                  myC, one, one, desc_myC)

                # merge local matrix to global
                C = Matrix{$elty}(undef, mA, nB)
                ScaLapack.pXgemr2d!(mA, nB,
                                    myC, one, one, desc_myC,
                                    C, one, one, desc_C,
                                    ctxt0)

            end

            BLACS.barrier(ctxt, 'A')
            BLACS.barrier(ctxt0, 'A')
            BLACS.gridexit(ctxt)
            BLACS.gridexit(ctxt0)

            # free
            if rank != rootproc
                A = Matrix{$elty}(undef, zero, zero)
                B = Matrix{$elty}(undef, zero, zero)
                C = Matrix{$elty}(undef, zero, zero)
            end

            return C

        end     # function

    end         # eval
end             # for

#--- find hessenberg matrix ---#
function hessenberg(sllm_A::ScaLapackLiteMatrix, reduced::Bool = true, get_Q::Bool = true)

    # assume: ScaLapackLiteParams must be correctly given for all processes
    params = sllm_A.params

    # MPI parameters
    mblocks = params.mblocks
    nblocks = params.nblocks
    mprocs = params.mprocs
    nprocs = params.nprocs
    rootproc = params.rootproc
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    # copy matrix from master to slave
    sllm_A = MPI.bcast(sllm_A, rootproc, comm)

    # decompose
    A = sllm_A.X

    # type of the matrix A
    elty = typeof(A[1,1])

    # check dimensions
    if rank == rootproc
        check_hessenberg(A)
        if mblocks != nblocks
            throw(DimensionMismatch("process grid must be contained from squared blocks"))
        end
    end

    # create context
    ctxt = ScaLapack.sl_init(mprocs, nprocs)
    ctxt0 = ScaLapack.sl_init(mprocs, nprocs)

    # matrix parameters
    m, n = convert.(ScaInt, size(A))
    mproc, nproc, myrow, mycol = BLACS.gridinfo(ctxt)
    mxlocr = ScaLapack.numroc(m, mblocks, mproc, zero, mprocs)
    mxlocc = ScaLapack.numroc(n, nblocks, nproc, zero, nprocs)
    mxllda = max(one, mxlocr)

    if mproc >= zero && nproc >= zero

        two = convert(ScaInt, 2)
        three = convert(ScaInt, 3)

        # Q matrix is enabled if get_Q = true
        Q = Matrix{elty}(undef, zero, zero)

        # get array descriptor
        desc_A = ScaLapack.descinit(m, n, m, n, zero, zero, ctxt0, m)
        desc_myA = ScaLapack.descinit(m, n, mblocks, nblocks, zero, zero, ctxt, mxllda)

        # prepare local matrix
        myA = zeros(elty, mxllda, mxlocc)
        ScaLapack.pXgemr2d!(m, n,
                            A, one, one, desc_A,
                            myA, one, one, desc_myA,
                            ctxt)

        # MPI params
        numblocks = mblocks
        rsrc_a = zero
        csrc_a = zero

        # find hessenberg matrix
        ilo = one; ihi = m;
        ia = one; ja = one;
        τ = zeros(elty, ScaLapack.numroc(ja+n-two, numblocks, mycol, csrc_a, nprocs))
        ScaLapack.pXgehrd!(m, ilo, ihi,
                           myA, ia, ja, desc_myA,
                           τ)

        if get_Q

            # temporary variables
            desc_Q = ScaLapack.descinit(m, n, m, n, zero, zero, ctxt0, m)
            desc_myQ = ScaLapack.descinit(m, n, mblocks, nblocks, zero, zero, ctxt, mxllda)
            Q = zeros(elty, m, n)
            myQ = zeros(elty, mxllda, mxlocc)
            identity!(myQ, m, n, desc_myQ)

            # generate orthogonal Q matrix; Q = Q * I
            # the orthogonal matrix Q is generated from the non-reduced Hessenberg matrix myA
            # Q = ( I-tau*v(1)*v(1)' )*( I-tau*v(2)*v(2)' )...
            # the v(j) is contained in the lower off-diag excepted the sub-siag of myA
            side = 'L'; trans = 'N';
            ScaLapack.pXYYmhr!(side, trans,
                               m, n, ilo, ihi,
                               myA, ia, ja, desc_myA,
                               τ,
                               myQ, one, one, desc_myQ)
            # merge local matrix to global
            ScaLapack.pXgemr2d!(m, n,
                                myQ, one, one, desc_myQ,
                                Q, one, one, desc_Q,
                                ctxt0)
        end

        if reduced
            # remove non-reduced part of the upper Hessenberg matrix
            uplo = 'L'
            α = convert(elty, zero); β = convert(elty, zero);
            ia = three; ja = one;
            ScaLapack.pXlaset!(uplo, m-two, n-two,
                               α, β,
                               myA, ia, ja, desc_myA)
        end

        # merge local matrix to global
        ScaLapack.pXgemr2d!(m, n,
                            myA, one, one, desc_myA,
                            A, one, one, desc_A,
                            ctxt0)

    end

    BLACS.barrier(ctxt, 'A')
    BLACS.barrier(ctxt0, 'A')
    BLACS.gridexit(ctxt)
    BLACS.gridexit(ctxt0)

    # free
    if rank != rootproc
        A = Matrix{elty}(undef, zero, zero)
        if get_Q; Q = Matrix{elty}(undef, zero, zero); end;
    end

    if get_Q; return (ScaLapackLiteMatrix(params, A), ScaLapackLiteMatrix(params, Q));
    else; return ScaLapackLiteMatrix(params, A); end;

end     # function

#--- Schur decomposition ---#
function schur(sllm_A::ScaLapackLiteMatrix, force_complex::Bool = false, get_T::Bool = true)

    # assume: ScaLapackLiteParams must be correctly given for all processes
    params = sllm_A.params

    # MPI parameters
    mblocks = params.mblocks
    nblocks = params.nblocks
    mprocs = params.mprocs
    nprocs = params.nprocs
    rootproc = params.rootproc
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    # perform find reduced upper hessenberg matrix
    # the arguments: general matrix A, reducing flag, return Q matrix
    (sllm_H, sllm_Q) = ScaLapackLite.hessenberg(sllm_A, true, true)

    # copy matrix from master to slave
    sllm_H = MPI.bcast(sllm_H, rootproc, comm)
    sllm_Q = MPI.bcast(sllm_Q, rootproc, comm)

    # type of the matrix H
    elty = typeof(sllm_H.X[1,1])
    elty_w = elty_e = elty
    if elty == Float32; elty_w = ComplexF32;
    elseif elty == Float64; elty_w = ComplexF64; end
    if force_complex; elty_e = elty_w; end

    # check dimensions
    if rank == rootproc
        check_schur(sllm_H.X)
        if mblocks != nblocks
            throw(DimensionMismatch("process grid must be contained from squared blocks"))
        end
    end

    # convert Float to ComplexF
    if force_complex
        sllm_H.X = convert(Matrix{elty_e}, sllm_H.X)
        sllm_Q.X = convert(Matrix{elty_e}, sllm_Q.X)
    end

    # decompose
    H = sllm_H.X; Q = sllm_Q.X;

    # create context
    ctxt = ScaLapack.sl_init(mprocs, nprocs)
    ctxt0 = ScaLapack.sl_init(mprocs, nprocs)

    # matrix parameters
    m, n = convert.(ScaInt, size(H))
    mproc, nproc, myrow, mycol = BLACS.gridinfo(ctxt)
    mxlocr = ScaLapack.numroc(m, mblocks, mproc, zero, mprocs)
    mxlocc = ScaLapack.numroc(n, nblocks, nproc, zero, nprocs)
    mxllda = max(one, mxlocr)

    if mproc >= zero && nproc >= zero

        two = convert(ScaInt, 2)
        three = convert(ScaInt, 3)

        # get array descriptor
        desc_H = ScaLapack.descinit(m, n, m, n, zero, zero, ctxt0, m)
        desc_Q = ScaLapack.descinit(m, n, m, n, zero, zero, ctxt0, m)
        desc_myH = ScaLapack.descinit(m, n, mblocks, nblocks, zero, zero, ctxt, mxllda)
        desc_myQ = ScaLapack.descinit(m, n, mblocks, nblocks, zero, zero, ctxt, mxllda)

        # generate process matrix
        myH = zeros(elty_e, mxllda, mxlocc)
        ScaLapack.pXgemr2d!(m, n,
                            H, one, one, desc_H,
                            myH, one, one, desc_myH,
                            ctxt)
        myQ = zeros(elty_e, mxllda, mxlocc)
        ScaLapack.pXgemr2d!(m, n,
                            Q, one, one, desc_Q,
                            myQ, one, one, desc_myQ,
                            ctxt)

        # perform Schur decomposition: H = S * T * S^H
        # its eigenvalues ( = diag(T) ) are stored into the w
        # myQ will be updated to the local part of Q * S
        #       because of the: H = Q^H * A * Q = S * T * S^H
        #                         <=> A = (Q*S) * T (Q*S)^H
        # if wantt = true: myA will be overwitten to T
        wantt = get_T; wantz = true;
        ilo = one; ihi = m;
        w = zeros(elty_w, m)
        # for real input
        if elty_e == Float32 || elty_e == Float64
            wr = zeros(elty, m); wi = zeros(elty, m)
            ScaLapack.pXlahqr!(wantt, wantz, m,
                               ilo, ihi, myH, desc_myH,
                               wr, wi,
                               ilo, ihi, myQ, desc_myQ)
            w = wr + im*wi
        # for complex input
        elseif elty_e == ComplexF32 || elty_e == ComplexF64
            ScaLapack.pXlahqr!(wantt, wantz, m,
                               ilo, ihi, myH, desc_myH,
                               w,
                               ilo, ihi, myQ, desc_myQ)
        end

        # merge local Schur vectors matrix to global
        ScaLapack.pXgemr2d!(m, n,
                            myQ, one, one, desc_myQ,
                            Q, one, one, desc_Q,
                            ctxt0)
        # merge local upper triangular matrix T
        if get_T
            ScaLapack.pXgemr2d!(m, n,
                                myH, one, one, desc_myH,
                                H, one, one, desc_H,
                                ctxt0)
        end

    end

    BLACS.barrier(ctxt, 'A')
    BLACS.barrier(ctxt0, 'A')
    BLACS.gridexit(ctxt)
    BLACS.gridexit(ctxt0)

    # free
    if rank != rootproc
        w = Vector{elty_w}(undef, zero)
        Q = Matrix{elty_e}(undef, zero, zero)
        H = Matrix{elty_e}(undef, zero, zero)
    end

    if get_T; return (w, ScaLapackLiteMatrix(params, Q), ScaLapackLiteMatrix(params, H));
    else; return (w, ScaLapackLiteMatrix(params, Q)); end;

end

#--- find eigenvalues ---#
function eigen(sllm_A::ScaLapackLiteMatrix, eidx::Vector{ScaInt} = [], is_leftside::Bool = false)

    # assume: ScaLapackLiteParams must be correctly given for all processes
    params = sllm_A.params

    # MPI parameters
    mblocks = params.mblocks
    nblocks = params.nblocks
    mprocs = params.mprocs
    nprocs = params.nprocs
    rootproc = params.rootproc
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    # perform Schur decomposition
    # the arguments: general matrix A, the flag to force convert to complex, return T matrix
    (eigvals, sllm_Z, sllm_T) = schur(sllm_A, true, true)
    sllm_Z = MPI.bcast(sllm_Z, rootproc, comm)
    sllm_T = MPI.bcast(sllm_T, rootproc, comm)

    # type of the matrix T
    elty = typeof(sllm_T.X[1,1])
    elty_ev = elty;
    if elty == Float32; elty_ev = ComplexF32;
    elseif elty == Float64; elty_ev = ComplexF64; end

    # check dimensions
    if rank == rootproc
        check_eigen(sllm_T.X)
        if mblocks != nblocks
            throw(DimensionMismatch("process grid must be contained from squared blocks"))
        end
    end

    # convert Float to ComplexF
    if elty != elty_ev
        sllm_T.X = convert(Matrix{elty_ev}, sllm_T.X)
        sllm_Z.X = convert(Matrix{elty_ev}, sllm_Z.X)
    end

    # decompose and copy Z to initial input V for pXtrevc
    T = sllm_T.X; Z = sllm_Z.X;

    # create context
    ctxt = ScaLapack.sl_init(mprocs, nprocs)
    ctxt0 = ScaLapack.sl_init(mprocs, nprocs)

    # matrix parameters
    m, n = convert.(ScaInt, size(T))
    mproc, nproc, myrow, mycol = BLACS.gridinfo(ctxt)
    mxlocr = ScaLapack.numroc(m, mblocks, mproc, zero, mprocs)
    mxlocc = ScaLapack.numroc(n, nblocks, nproc, zero, nprocs)
    mxllda = max(one, mxlocr)

    # eigenvector parameters
    mm = m; neigs = convert(ScaInt, length(eidx));
    if neigs > zero
        # specified eigen vectors
        select = Vector{Bool}([ false for j = one:m ])
        select[eidx] .= true
    else
        if neigs < zero
            throw(DimensionMismatch("illegal value; num of specified eigen index = $neigs"))
        else
            # all eigen vectors ( neigs = 0 )
            eidx = Vector{ScaInt}([ j for j = one:m ])
            select = Vector{Bool}([ true for j = one:m ])
        end
    end

    if mproc >= zero && nproc >= zero

        # get array descriptor
        desc_T = ScaLapack.descinit(m, n, m, n, zero, zero, ctxt0, m)
        desc_Z = ScaLapack.descinit(m, n, m, n, zero, zero, ctxt0, m)
        desc_myT = ScaLapack.descinit(m, n, mblocks, nblocks, zero, zero, ctxt, mxllda)
        desc_myZ = ScaLapack.descinit(m, n, mblocks, nblocks, zero, zero, ctxt, mxllda)

        # generate process matrix of T and Z
        myT = zeros(elty_ev, mxllda, mxlocc)
        myZ = zeros(elty_ev, mxllda, mxlocc)
        ScaLapack.pXgemr2d!(m, n,
                            T, one, one, desc_T,
                            myT, one, one, desc_myT,
                            ctxt)
        ScaLapack.pXgemr2d!(m, n,
                            Z, one, one, desc_Z,
                            myZ, one, one, desc_myZ,
                            ctxt)

        # find eigenvectors
        o = zeros(ScaInt, zero)
        O = zeros(elty_ev, zero, zero)
        if is_leftside; side = 'L'; else; side = 'R'; end;
        if neigs == zero; howmny = 'B'; else; howmny = 'S'; end
        # for each side
        nV = zero
        if is_leftside
            nV = ScaLapack.pXtrevc!(side, howmny, select,
                                    m, myT, desc_myT,
                                    myZ, desc_myZ, O, o)
        else
            nV = ScaLapack.pXtrevc!(side, howmny, select,
                                    m, myT, desc_myT,
                                    O, o, myZ, desc_myZ)
        end

        # merge local eigen vectors matrix to global
        ScaLapack.pXgemr2d!(m, n,
                            myZ, one, one, desc_myZ,
                            Z, one, one, desc_Z,
                            ctxt0)

    end

    BLACS.barrier(ctxt, 'A')
    BLACS.barrier(ctxt0, 'A')
    BLACS.gridexit(ctxt)
    BLACS.gridexit(ctxt0)

    # free
    if rank != rootproc
        eigvals = Vector{elty_ev}(undef, zero)
        T = Matrix{elty_ev}(undef, zero, zero)
        Z = Matrix{elty_ev}(undef, zero, zero)
    else
        # store actual size of result
        eigvals = eigvals[eidx]
        Z = Z[:, eidx]
    end

    return (eigvals, Z)

end

#--- ERROR dumpings ---#
function check_multiple(A::Matrix{<:Number}, B::Matrix{<:Number})
    mA, nA = size(A)
    mB, nB = size(B)
    if nA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA, $nA), matrix B has dimensions ($mB, $nB)"))
    end
end
function check_hessenberg(A::Matrix{<:Number})
    mA, nA = size(A)
    if mA != nA
        throw(DimensionMismatch("matrix A has no-squared dimensions ($mA, $nA)"))
    end
end
function check_schur(A::Matrix{<:Number})
    mA, nA = size(A)
    if mA != nA
        throw(DimensionMismatch("matrix A has no-squared dimensions ($mA, $nA)"))
    end
end
function check_eigen(A::Matrix{<:Number})
    mA, nA = size(A)
    if mA != nA
        throw(DimensionMismatch("matrix A has no-squared dimensions ($mA, $nA)"))
    end
end


text:
#= src/pqr.jl

References:

  P. Businger, G.H. Golub. Linear least squares solutions by Householder
    transformations. Numer. Math. 7: 269-276, 1965.

  M. Gu, S.C. Eisenstat. Efficient algorithms for computing a strong
    rank-revealing QR factorization. SIAM J. Sci. Comput. 17 (4): 848-869, 1996.
=#

# PartialQRFactors

mutable struct PartialQRFactors
  Q::Nullable{Matrix}
  R::Nullable{Matrix}
  p::Vector{Int}
  k::Int
  T::Nullable{Matrix}
end
const PQRFactors = PartialQRFactors

function getindex(A::PQRFactors, d::Symbol)
  if     d == :P  return ColumnPermutation(A.p)
  elseif d == :Q  return get(A.Q)
  elseif d == :R  return UpperTrapezoidal(get(A.R))
  elseif d == :T  return get(A.T)
  elseif d == :k  return A.k
  elseif d == :p  return A.p
  else            throw(KeyError(d))
  end
end

# PartialQR

mutable struct PartialQR{T} <: Factorization{T}
  Q::Matrix{T}
  R::Matrix{T}
  p::Vector{Int}
end

conj!(A::PartialQR) = PartialQR(conj!(A.Q), conj!(A.R), A.p)
conj(A::PartialQR) = PartialQR(conj(A.Q), conj(A.R), A.p)

convert(::Type{PartialQR{T}}, A::PartialQR) where {T} =
  PartialQR(convert(Array{T}, A.Q), convert(Array{T}, A.R), A.p)
convert(::Factorization{T}, A::PartialQR) where {T} = convert(PartialQR{T}, A)
convert(::Type{Array}, A::PartialQR) = Matrix(A)
convert(::Type{Array{T}}, A::PartialQR) where {T} = convert(Array{T}, Matrix(A))
convert(::Type{Matrix}, A::PartialQR) = Matrix(A)
convert(::Type{Matrix{T}}, A::PartialQR) where {T} = convert(Array{T}, Matrix(A))
Array(A::PartialQR)  = convert(Array, A)

copy(A::PartialQR) = PartialQR(copy(A.Q), copy(A.R), copy(A.p))

  Matrix(A::PartialQR) = rmul!(A[:Q]*A[:R], A[:P]')
  adjoint(A::PartialQR) = Adjoint(A)
  transpose(A::PartialQR) = Transpose(A)

function getindex(A::PartialQR, d::Symbol)
  if     d == :P  return ColumnPermutation(A.p)
  elseif d == :Q  return A.Q
  elseif d == :R  return UpperTrapezoidal(A.R)
  elseif d == :k  return size(A.Q, 2)
  elseif d == :p  return A.p
  else            throw(KeyError(d))
  end
end

ishermitian(::PartialQR) = false
issym(::PartialQR) = false

isreal(::PartialQR{T}) where {T} = T <: Real

ndims(::PartialQR) = 2

size(A::PartialQR) = (size(A.Q,1), size(A.R,2))
size(A::PartialQR, dim::Integer) =
  dim == 1 ? size(A.Q,1) : (dim == 2 ? size(A.R,2) : 1)

# BLAS/LAPACK multiplication/division routines

  ## left-multiplication

function mul!!(C::AbstractVector{T}, A::PartialQR{T}, B::AbstractVector{T}) where T
    lmul!(A[:P]', B)
    mul!(C, A[:Q], A[:R]*B)
end  # overwrites B
function mul!!(C::AbstractMatrix{T}, A::PartialQR{T}, B::AbstractMatrix{T}) where T
    lmul!(A[:P]', B)
    mul!(C, A[:Q], A[:R]*B)
end  # overwrites B
mul!(C::AbstractVector{T}, A::PartialQR{T}, B::AbstractVector{T}) where {T} =
    mul!!(C, A, copy(B))

mul!(C::AbstractMatrix{T}, A::PartialQR{T}, B::AbstractMatrix{T}) where {T} =
    mul!!(C, A, copy(B))

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!(C::AbstractMatrix{T}, A::PartialQR{T}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        tmp = $Adj(A[:P]) * Bc
        mul!(C, A[:Q], A[:R]*tmp)
      end
      function mul!(C::AbstractVector{T}, Ac::$Adj{T,PartialQR{T}}, B::AbstractVector{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:Q]) * B
        mul!(C, $Adj(A[:R]), tmp)
        lmul!(A[:P], C)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,PartialQR{T}}, B::AbstractMatrix{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:Q]) * B
        mul!(C, $Adj(A[:R]), tmp)
        lmul!(A[:P], C)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,PartialQR{T}}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        A = parent(Ac)
        tmp = $Adj(A[:Q]) * Bc
        mul!(C, $Adj(A[:R]), tmp)
        lmul!(A[:P], C)
      end
    end
  end

  ## right-multiplication

  function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::PartialQR{T}) where T
    mul!(C, A*B[:Q], B[:R])
    rmul!(C, B[:P]')
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,PartialQR{T}}) where T
        B = parent(Bc)
        rmul!(A, B[:P])
        tmp = A * $Adj(B[:R])
        mul!(C, tmp, $Adj(B[:Q]))
      end  # overwrites A
      mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,PartialQR{T}}) where {T} =
        mul!!(C, copy(A), Bc)
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::PartialQR{T}) where T
        tmp = Ac * B[:Q]
        mul!(C, tmp, B[:R])
        rmul!(C, $Adj(B[:P]))
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, Bc::$Adj{T,PartialQR{T}}) where T
        B = parent(Bc)
        tmp = Ac * B[:P]
        tmp = tmp * $Adj(B[:R])
        mul!(C, tmp, $Adj(B[:Q]))
      end
    end
  end

  ## left-division (pseudoinverse left-multiplication)
  function ldiv!(C::AbstractVector{T}, A::PartialQR{T}, B::AbstractVector{T}) where T
    tmp = (A[:R]*A.R')\(A[:Q]'*B)
    mul!(C, A[:R]', tmp)
    lmul!(A[:P], C)
  end

function ldiv!(C::AbstractMatrix{T}, A::PartialQR{T}, B::AbstractMatrix{T}) where T
    tmp = (A[:R]*A.R')\(A[:Q]'*B)
    mul!(C, A[:R]', tmp)
    lmul!(A[:P], C)
end

  # standard operations

  ## left-multiplication

  function *(A::PartialQR{TA}, B::AbstractVector{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(PartialQR{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1))
    mul!(CT, AT, BT)
  end
  function *(A::PartialQR{TA}, B::AbstractMatrix{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(PartialQR{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end
  function *(A::AbstractMatrix{TA}, B::PartialQR{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = (T == TA ? A : convert(Array{T}, A))
    BT = convert(PartialQR{T}, B)
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function *(Ac::$Adj{TA,PartialQR{TA}}, B::AbstractVector{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(PartialQR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1))
        mul!(CT, $Adj(AT), BT)
      end
      function *(A::PartialQR{TA}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(PartialQR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,PartialQR{TA}}, B::AbstractMatrix{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(PartialQR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,PartialQR{TA}}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(PartialQR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
      function *(A::AbstractMatrix{TA}, Bc::$Adj{TB,PartialQR{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialQR{T}, B)
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, B::PartialQR{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialQR{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, Bc::$Adj{TB,PartialQR{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialQR{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
    end
  end

## left-division
function \(A::PartialQR{TA}, B::AbstractVector{TB}) where {TA,TB}
  T = promote_type(TA, TB)
  AT = convert(PartialQR{T}, A)
  BT = (T == TB ? B : convert(Array{T}, B))
  CT = Array{T}(undef, size(A,2))
  ldiv!(CT, AT, BT)
end
function \(A::PartialQR{TA}, B::AbstractMatrix{TB}) where {TA,TB}
  T = promote_type(TA, TB)
  AT = convert(PartialQR{T}, A)
  BT = (T == TB ? B : convert(Array{T}, B))
  CT = Array{T}(undef, size(A,2), size(B,2))
  ldiv!(CT, AT, BT)
end

# factorization routines

  for sfx in ("", "!")
    f = Symbol("pqrfact", sfx)
    g = Symbol("pqrfact_none", sfx)
    h = Symbol("pqr", sfx)
    @eval begin
      function $f(
          trans::Symbol, A::AbstractMatOrLinOp{S}, opts::LRAOptions=LRAOptions(S);
          args...) where S
        chktrans(trans)
        opts = copy(opts; args...)
        chkopts!(opts, A)
        opts.sketch == :none && return $g(trans, A, opts)
        V = idfact(trans, A, opts)
        F = qr!(convert(Matrix, getcols(trans, A, V[:sk])))
        retq = occursin("q", opts.pqrfact_retval)
        retr = occursin("r", opts.pqrfact_retval)
        rett = occursin("t", opts.pqrfact_retval)
        Q = retq ? Matrix(F.Q) : nothing
        R = retr ? pqrr(F.R, V[:T]) : nothing
        T = rett ? V[:T] : nothing
        retq && retr && !rett && return PartialQR(Q, R, V[:p])
        PQRFactors(Q, R, V[:p], V[:k], T)
      end
      $f(trans::Symbol, A, args...; kwargs...) =
        $f(trans, LinOp(A), args...; kwargs...)
      $f(A, args...; kwargs...) = $f(:n, A, args...; kwargs...)

      function $h(trans::Symbol, A, args...; kwargs...)
        defaults = (;pqrfact_retval="qr")
        kwargs = merge(defaults,kwargs)
        F = $f(trans, A, args...; kwargs...)
        F.Q, F.R, F.p
      end
      $h(A, args...; kwargs...) = $h(:n, A, args...; kwargs...)
    end
  end

pqrfact_none!(trans::Symbol, A::AbstractMatrix, opts::LRAOptions) =
  pqrfact_backend!(convert(Matrix, trans == :n ? A : A'), opts)

  pqrfact_none(trans::Symbol, A::AbstractMatrix, opts::LRAOptions) =
    pqrfact_backend!(Matrix(trans == :n ? A : A'), opts)

function pqrr(R::Matrix{S}, T::Matrix{S}) where S
  k, n = size(T)
  n += k
  R_ = Array{S}(undef, k, n)
  R1 = view(R_, :,   1:k)
  R2 = view(R_, :, k+1:n)
  copyto!(R1, R)
  copyto!(R2, T)
  lmul!(UpperTriangular(R1), R2)
  R_
end

## core backend routine: rank-adaptive GEQP3 with determinant maximization
function pqrfact_backend!(A::AbstractMatrix, opts::LRAOptions)
  p, tau, k = geqp3_adap!(A, opts)
  pqrback_postproc(A, p, tau, k, opts)
end

function geqp3_adap!(A::AbstractMatrix{T}, opts::LRAOptions) where T
  m, n = size(A)
  jpvt = collect(BlasInt, 1:n)
  l    = min(m, n)
  k    = (opts.rank < 0 || opts.rank > l) ? l : opts.rank
  tau  = Array{T}(undef, k)
  if k > 0
    k = geqp3_adap_main!(A, jpvt, tau, opts)
  end
  jpvt = convert(Array{Int}, jpvt)
  jpvt, tau, k
end

function geqp3_adap_main!(
    A::AbstractMatrix{T}, jpvt::Vector{BlasInt}, tau::Vector{T},
    opts::LRAOptions) where T<:BlasFloat
  chkstride1(A)
  lda = stride(A, 2)
  n   = length(jpvt)
  k   = length(tau)

  # set block size and allocate work array
  nb      = min(opts.nb, k)
  is_real = T <: Real
  lwork   = 2*n*is_real + (n + 1)*nb
  work    = Array{T}(undef, lwork)

  # initialize column norms
  if is_real
    @inbounds for j = 1:n
      work[j] = work[n+j] = norm(view(A,:,j))
    end
  else
    rwork = Array{eltype(real(zero(T)))}(undef, 2*n)
    @inbounds for j = 1:n
      rwork[j] = rwork[n+j] = norm(view(A,:,j))
    end
  end
  maxnrm = maximum(view(is_real ? work : rwork, 1:n))

  # set pivot threshold
  ptol = max(opts.atol, opts.rtol*maxnrm)

  # block factorization
  j = 1
  fjb = Ref{BlasInt}()
  while j <= k
    jb = BlasInt(min(nb, k-j+1))
    if is_real
      _LAPACK.laqps!(
        BlasInt(j-1), jb, fjb, view(A,:,j:n), view(jpvt,j:n), view(tau,j:k),
        view(work,j:n), view(work,n+j:2*n),
        view(work,2*n+1:2*n+nb), view(work,2*n+jb+1:lwork))
    else
      _LAPACK.laqps!(
        BlasInt(j-1), jb, fjb, view(A,:,j:n), view(jpvt,j:n), view(tau,j:k),
        view(rwork,j:n), view(rwork,n+j:2*n),
        view(work,1:nb), view(work,jb+1:lwork))
    end
    jn = j + fjb[]

    # check for rank termination
    if abs(A[jn-1,jn-1]) <= ptol
      @inbounds for i = j:jn-1
        abs(A[i,i]) <= ptol && return i - 1
      end
    end
    j = jn
  end
  k
end

function pqrback_postproc(
    A::AbstractMatrix{S}, p::Vector{Int}, tau::Vector{S}, k::Integer,
    opts::LRAOptions) where S
  retq = occursin("q", opts.pqrfact_retval)
  retr = occursin("r", opts.pqrfact_retval)
  rett = occursin("t", opts.pqrfact_retval)
  maxdet = 0 < k < size(A,2) && opts.maxdet_tol >= 0
  Q = retq ? LAPACK.orgqr!(A[:,1:k], tau, k) : nothing
  R = retr || rett || maxdet ? triu!(A[1:k,:]) : nothing
  T = rett || maxdet ? maxdet_t(R) : nothing
  if maxdet
    maxdet_swapcols!(Q, R, p, T, opts)
    R = retr ? R : nothing
  end
  retq && retr && !rett && return PartialQR(Q, R, p)
  PQRFactors(Q, R, p, k, T)
end

function maxdet_t(R::AbstractMatrix{S}) where S
  k, n = size(R)
  T = R[:,k+1:n]
  ldiv!(UpperTriangular(view(R,1:k,1:k)), T)
end

## rank-revealing QR determinant maximization
function maxdet_swapcols!(
    Q::Union{Matrix{S}, Nothing}, R::Matrix{S}, p::Vector{Int}, T::Matrix{S},
    opts::LRAOptions) where S
  k, n  = size(R)
  R1    = view(R, :, 1:k)
  work  = Array{S}(undef, max(n, 2*k))
  retq  = occursin("q", opts.pqrfact_retval)
  retr  = occursin("r", opts.pqrfact_retval)
  niter = 0
  while true
    Tmax, idx = findmaxabs(T)
    Tmax <= 1 + opts.maxdet_tol && break
    if niter == opts.maxdet_niter
      if opts.verb
        @warn "iteration limit ($niter) reached in determinant maximization"
      end
      break
    end
    niter += 1
    i, j = Tuple(ind2sub((k, n-k), idx))
    maxdet_update!(R1, p, T, work, i, j, retr)
  end
  niter == 0 && return
  F = qrfact!(R1)
  retq && LAPACK.gemqrt!('R', 'N', F.factors, F.T, Q)
  if retr
    triu!(R1)
    R2 = view(R, :, k+1:n)
    copyto!(R2, T)
    lmul!(UpperTriangular(R1), R2)
  end
end

## column swap update based on Sherman-Morrison
function maxdet_update!(
    R1::AbstractMatrix{S}, p::Vector{Int}, T::Matrix{S}, work::Vector{S},
    i::Integer, j::Integer, retr::Bool) where S
  k, n = size(T)
  n += k
  p[i], p[k+j] = p[k+j], p[i]
  @inbounds @simd for l = 1:k
    work[l] = T[l,j]
    T[l,j]  = 0
  end
  work[i] -= 1
  T[i,j]   = 1
  @inbounds for l = 1:n-k
    work[k+l] = conj(T[i,l])
  end
  BLAS.ger!(-1/(1 + work[i]), view(work,1:k), view(work,k+1:n), T)
  if retr
    mul!(view(work,k+1:2*k), R1, view(work,1:k))
    @inbounds @simd for l = 1:k
      R1[l,i] += work[k+l]
    end
  end
end


text:
import java.io.BufferedReader;
import java.io.ByteArrayInputStream;
import java.io.ByteArrayOutputStream;
import java.io.InputStreamReader;
import java.io.PrintStream;
import java.nio.file.Files;
import java.nio.file.OpenOption;
import java.nio.file.Path;
import java.nio.file.Paths;
import java.nio.file.StandardOpenOption;
import java.security.KeyPair;
import java.security.KeyPairGenerator;
import java.security.PrivateKey;
import java.security.PublicKey;
import java.util.ArrayList;

import accrue.cryptoerase.runtime.CryptoLibrary;

public class Log {
    private final SimpleIRCClient parent;
    private static final OpenOption[] defaultOpenOptions = { StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.CREATE };
    // ANNOTATION: key declaration
    private PrivateKey PRIVKEY(L /parent.clearHistory T){L} privkey = null;
    // ANNOTATION: key declaration
    private PublicKey PUBKEY(L /parent.clearHistory T){L} pubkey = null;

    // ANNOTATION: field erasure type
    private ArrayList{L /parent.clearHistory T} buffer;
    private int historyLength;

    public Log(SimpleIRCClient parent) {
    	this.parent = parent;
    	KeyPair kp = null;
    	try {
    		KeyPairGenerator kpg = KeyPairGenerator.getInstance("RSA");
    		kpg.initialize(2048); // 2048 is the keysize
    		kp = kpg.generateKeyPair();
    	} catch (Exception e) {}

	privkey = kp.getPrivate();
	pubkey = kp.getPublic();
	//    	privkey = (PRIVKEY(L /parent.clearHistory T){L}) kp.getPrivate();
	//    	pubkey = (PUBKEY(L /parent.clearHistory T){L}) kp.getPublic();

    	buffer = new ArrayList();
    	historyLength = 0;
    }

    public int historyLength() {
    	return (this.historyLength < 1) ? 0 : this.historyLength - 1;
    }

    public ArrayList getScreen(int num) {
	try {
	    if (num < 0 || num >= historyLength - 1) {
    		return null;
	    }

	    Path file = Paths.get("irclog." + num, new String[0]);
	    byte[] bytes = Files.readAllBytes(file);
	    ArrayList compressed = CryptoLibrary.decryptStrings(privkey, bytes);
	    return CryptoLibrary.decompressStrings(compressed);
	} catch (Throwable e) {}
	return null;
    }

    public void writeLine(String line) {
	try {
	    int bufferSize = 0;
	    try {
		buffer.add(line);
		bufferSize = buffer.size();
	    } catch (Throwable e) {}
	    // Declassify whether we've received more than 20 messages since the last write
	    // ANNOTATION: declassification
	    boolean writeOut = ({L}) bufferSize >= 20;
	    if (writeOut) {
		ArrayList compressed = null;
		try { compressed = CryptoLibrary.compressStrings(buffer); } catch (Throwable e) {}
		byte[] bytes = CryptoLibrary.encryptStrings(pubkey, compressed);
		/*
		   // INJECTABLE FAULT: attempt to write high data to disk
		   byte[] bytes = new byte[100];
		   bytes[0] = ({H}) ((byte)42);
		*/
		Path file = Paths.get("irclog." + historyLength, new String[0]);
		Files.write(file, bytes, defaultOpenOptions);
		buffer.clear();
		historyLength++;
	    }
	} catch (Throwable e) {}
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(AcuteML.render2file)),NamedTuple{(:id, :age, :field, :GPA, :courses),Tuple{Int64,Int64,String,Float64,Array{String,1}}},typeof(render2file),String,Bool})
    Base.precompile(Tuple{typeof(AcuteML.aml_create),Expr,Array{Union{Expr, Symbol},1},Array{Any,1},Array{Union{Expr, Symbol, Type},1},Array{Union{Expr, Symbol},1},Array{Union{Missing, String},1},Array{Union{Missing, Function, Symbol},1},Array{Union{Missing, Type},1},String,Type,Array{Union{Missing, Function, Symbol},0},Bool,Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Expr})
    Base.precompile(Tuple{typeof(AcuteML.aml_create),Expr,Array{Union{Expr, Symbol},1},Array{Any,1},Array{Union{Expr, Symbol, Type},1},Array{Union{Expr, Symbol},1},Array{Union{Missing, String},1},Array{Union{Missing, Function, Symbol},1},Array{Union{Missing, Type},1},String,Type,Array{Union{Missing, Function, Symbol},0},Bool,Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Symbol})
    Base.precompile(Tuple{typeof(AcuteML.aml_parse),Expr})
    Base.precompile(Tuple{typeof(AcuteML.findfirsttext),Int64,Node})
    Base.precompile(Tuple{typeof(AcuteML.get_arg_xmlcreator),Expr,Bool,Expr,Expr,String,Type,Expr,QuoteNode,Expr})
    Base.precompile(Tuple{typeof(AcuteML.get_arg_xmlcreator),Nothing,Bool,Expr,Expr,String,Type,Expr,QuoteNode,Expr})
    Base.precompile(Tuple{typeof(AcuteML.get_arg_xmlextractor),Nothing,Bool,Expr,Expr,String,Type,Expr,QuoteNode,Expr})
    Base.precompile(Tuple{typeof(AcuteML.get_arg_xmludpater),Nothing,Bool,Expr,Expr,String,Type,Expr,QuoteNode,Expr})
    Base.precompile(Tuple{typeof(AcuteML.get_struct_xmlcreator),Expr,Array{Any,1},Array{Union{Expr, Symbol},1},Nothing,Expr,Array{Expr,1},Array{Expr,1},Nothing})
    Base.precompile(Tuple{typeof(AcuteML.get_struct_xmlextractor),Expr,Array{Any,1},Array{Expr,1},Nothing,Array{Expr,1},Nothing})
    Base.precompile(Tuple{typeof(AcuteML.multiString),Array{String,1}})
    Base.precompile(Tuple{typeof(AcuteML.multiString),Float64})
    Base.precompile(Tuple{typeof(AcuteML.multiString),Int64})
    Base.precompile(Tuple{typeof(AcuteML.multiString),String})
    Base.precompile(Tuple{typeof(AcuteML.nodeparse),Type,Node})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Array{Any,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Array{Float64,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Array{Int64,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Array{String,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Nothing,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,String,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Array{Any,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Array{Float64,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Array{Int64,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Int64,Type{AbsText}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Nothing,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(convert),Type{Union{Missing, Type}},Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(convert),Type{Union{Missing, Type}},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(convert),Type{Union{Missing, Type}},Type{AbsText}})
    Base.precompile(Tuple{typeof(createnode),Type{AbsEmpty},String})
    Base.precompile(Tuple{typeof(createnode),Type{AbsHtml},String})
    Base.precompile(Tuple{typeof(createnode),Type{AbsHtml}})
    Base.precompile(Tuple{typeof(createnode),Type{AbsXml}})
    Base.precompile(Tuple{typeof(findcontent),String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Any,1}},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Any,1}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Float64,1}},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Float64,1}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Int64,1}},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Int64,1}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{String,1}},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Int64},String,Node,Type{AbsText}})
    Base.precompile(Tuple{typeof(findcontent),Type{Nothing},String,Document,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(findcontent),Type{Nothing},String,Node,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(findcontent),Type{String},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Union{Nothing, String}},String,Document,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(findcontent),Type{Union{Nothing, String}},String,Node,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(newTemplate),String,Symbol})
    Base.precompile(Tuple{typeof(newTemplate),String})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Any,1},Array{Node,1},Node,Type})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Any,1},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Any,1},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Float64,1},Array{Node,1},Node,Type})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Float64,1},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Float64,1},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Int64,1},Array{Node,1},Node,Type})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Int64,1},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Int64,1},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{String,1},String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{String,1},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Int64,String,Node,Type{AbsText}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Document,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Document,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Node,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Node,Type{AbsText}})
    isdefined(AcuteML, Symbol("##5#6")) && Base.precompile(Tuple{getfield(AcuteML, Symbol("##5#6")),Symbol})
end


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

__precompile__(true)

module Devices

# package code goes here

end # module


text:
s = """
r = 5:1
"""
msgs = lintstr(s)
@test msgs[1].code == :E433
@test contains(msgs[1].message, "for a decreasing range, use a negative step e.g. 10:-1:1")

s = """
x = [1,2,7,8]
y = [3,4,5,6]
splice!(x, 3:2, y)
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(r::UnitRange)
    return r == 0:-1
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(r::UnitRange)
    a = r[2]
    b = r[3:4]
    @lintpragma("Info type r")
    @lintpragma("Info type a")
    @lintpragma("Info type b")
    (a,b)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(r) == UnitRange")
@test msgs[2].code == :I271
@test contains(msgs[2].message, "typeof(a) == Any")
@test msgs[3].code == :I271
@test contains(msgs[3].message, "typeof(b) == UnitRange")


text:
using SHIPs, LinearAlgebra

trans = PolyTransform(2, 1.0)
fcut = PolyCutoff1s(2, 0.1, 2.0)
basis = SHIPBasis(SparseSHIP(5, 8, 1.0), trans, fcut)

length(basis)

A = SHIPs.Rotations.CoeffArray(20)

for bo = 2:5
   ctr = 0
   println();
   @show bo
   for νz in basis.NuZ[bo]
      ν = νz.ν
      kk = getfield.(basis.KL[1][ν], :k)
      if norm(kk) != 0
         continue
      end
      ll = getfield.(basis.KL[1][ν], :l)
      U = SHIPs.Rotations.basis(A, ll)
      c = [ SHIPs._Bcoeff(ll, mm, A.cg) for mm in SHIPs._mrange(ll) ]
      def = norm(U * (U' * c) - c) > 1e-12
      if def > 1e-12
         @show ll, def
         ctr += 1
         if ctr > 10
            break
         end
      else
         print("✓")
      end
   end
end

using StaticArrays

for bo = 2:5
   println()
   @show bo
   ctr = 0
   cartrg = CartesianIndices( ntuple(_->0:3, bo) )
   for ill in cartrg
      ll = SVector(ill.I...)
      if issorted(ll) && sum(ll) <= 20
         U =  try
            SHIPs.Rotations.basis(A, ll)
         catch
            println()
            @show ll
            continue
         end
         c = [ SHIPs._Bcoeff(ll, mm, A.cg) for mm in SHIPs._mrange(ll) ]
         def = norm(U * (U' * c) - c)
         if def > 1e-12
            @show ll, def
            ctr += 1
            if ctr > 10
               break
            end
         else
            print("✓")
         end
      end
   end
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using LinearAlgebra, IntervalSets, UnPack
import ClimaCore: Domains, Topologies, Meshes, Spaces, Geometry

using Test

@testset "Surface area of a sphere" begin
    FT = Float64
    radius = FT(3)
    ne = 4
    Nq = 4
    domain = Domains.SphereDomain(radius)
    mesh = Meshes.EquiangularCubedSphere(domain, ne)
    grid_topology = Topologies.Topology2D(mesh)
    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    @test sum(ones(space)) ≈ 4pi * radius^2 rtol = 1e-3
end

@testset "Volume of a spherical shell" begin
    FT = Float64
    radius = FT(128)
    zlim = (0, 1)
    helem = 4
    zelem = 10
    Nq = 4

    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = zelem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.SphereDomain(radius)
    horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
    horztopology = Topologies.Topology2D(horzmesh)
    quad = Spaces.Quadratures.GLL{Nq}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)

    # "shallow atmosphere" spherical shell: volume = surface area * height
    @test sum(ones(hv_center_space)) ≈ 4pi * radius^2 * (zlim[2] - zlim[1]) rtol =
        1e-3
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

include("adstring.jl")
export adstring

include("airtovac.jl")
export airtovac

include("aitoff.jl")
export aitoff

include("altaz2hadec.jl")
export altaz2hadec

include("bprecess.jl")
export bprecess

include("calz_unred.jl")
export calz_unred

include("ct2lst.jl")
export ct2lst

include("daycnv.jl")
export daycnv

include("deredd.jl")
export deredd

include("eci2geo.jl")
export eci2geo

include("eqpole.jl")
export eqpole

include("flux2mag.jl")
export flux2mag

include("gal_uvw.jl")
export gal_uvw

include("geo2eci.jl")
export geo2eci

include("geo2geodetic.jl")
export geo2geodetic

include("geo2mag.jl")
export geo2mag

include("geodetic2geo.jl")
export geodetic2geo

include("get_date.jl")
export get_date

include("get_juldate.jl")
export get_juldate

include("gcirc.jl")
export gcirc

include("hadec2altaz.jl")
export hadec2altaz

include("helio_jd.jl")
export helio_jd

include("helio_rv.jl")
export helio_rv

include("jdcnv.jl")
export jdcnv

include("jprecess.jl")
export jprecess

include("juldate.jl")
export juldate

include("kepler_solver.jl")
export kepler_solver

include("lsf_rotate.jl")
export lsf_rotate

include("mag2flux.jl")
export mag2flux

include("mag2geo.jl")
export mag2geo

include("month_cnv.jl")
export month_cnv

include("moonpos.jl")
export moonpos

include("mphase.jl")
export mphase

include("nutate.jl")
export nutate

include("paczynski.jl")
export paczynski

include("planck_freq.jl")
export planck_freq

include("planck_wave.jl")
export planck_wave

include("polrec.jl")
export polrec

include("posang.jl")
export posang

include("precess.jl")
export precess

include("precess_xyz.jl")
export precess_xyz

include("premat.jl")
export premat

include("radec.jl")
export radec

include("recpol.jl")
export recpol

include("rhotheta.jl")
export rhotheta

include("sixty.jl")
export sixty

include("sphdist.jl")
export sphdist

include("sunpos.jl")
export sunpos

include("ten.jl")
export ten, tenv

include("trueanom.jl")
export trueanom

include("vactoair.jl")
export vactoair

include("xyz.jl")
export xyz

include("ydn2md.jl")
export ydn2md

include("ymd2dn.jl")
export ymd2dn


text:
struct Record_1
    a
    b
end

@as_record Record_1

@test @match Record_1(1, 2) begin
    Record_1(1, 2) => true
    _ => false
end

@test @match Record_1(1, 2) begin
    Record_1(a=1) => true
    _ => false
end

@test @match Record_1(1, 2) begin
    Record_1(b = 2) => true
    _ => false
end

@test @match Record_1(1, 2) begin
    Record_1(;a=1) => true
    _ => false
end

@test @match Record_1(1, 2) begin
    Record_1(;b = 2) => true
    _ => false
end

@test @match Record_1(1, 2) begin
    Record_1(;b) => b == 2
    _ => false
end

struct Record_2{A}
    x :: A
    y :: Int
end

@as_record visible in (@__MODULE__) Record_2

@test @match Record_2(1, 2) begin
    Record_2{A}(_) where A => A == typeof(1)
    _ => false
end


text:
using RoadRunnerJulia
#println("in the callAPI file: ", rrlib)

RoadRunnerJulia.setConfigInt("LOADSBMLOPTIONS_CONSERVED_MOIETIES", 1)
#rrlib = Libdl.dlopen("C:/vs_rebuild/install/roadrunner/bin/roadrunner_c_api.dll")
ant_str = """
    const Xo, X1
     Xo -> S1; k1*Xo - k2*S1
     S1 -> S2; k3*S2
     S2 -> X1; k4*S2

     Xo = 1;   X1 = 0
     S1 = 0;   S2 = 0
     k1 = 0.1; k2 = 0.56
     k3 = 1.2; k4 = 0.9
"""
rr = loada(ant_str)
println(RoadRunnerJulia.getConfigInt("LOADSBMLOPTIONS_CONSERVED_MOIETIES"))
RoadRunnerJulia.setConfigInt("sfds", 1)
println(RoadRunnerJulia.getConfigInt("sfds"))
data = simulate(rr, 0, 50, 51)
println(data)
ss = steadyState(rr)
println("this is the steady state value: ", ss)
ssValues = RoadRunnerJulia.computeSteadyStateValues(rr)
# str = RoadRunnerJulia.getListOfConfigKeys()
ids = getFloatingSpeciesIds(rr)
#data = RoadRunnerJulia.simulate(rr)
#println(data)


text:
# Specialized printing/showing methods
#

function Base.print(io::IO, seq::LongSequence{A}; width::Integer = 0) where {A<:Alphabet}
    return _print(io, seq, width, codetype(A()))
end

# Dispatch to generic method in biosequences/printing.jl
function _print(io::IO, seq::LongSequence{<:Alphabet}, width::Integer, ::AlphabetCode)
    return _print(SimpleBuffer(io), seq, width)
end

# Specialized method for ASCII alphabet
function _print(io::IO, seq::LongSequence{<:Alphabet}, width::Integer, ::AsciiAlphabet)
    # If seq is large, always buffer for memory efficiency
    if length(seq) ≥ 4096
        return _print(SimpleBuffer(io), seq, width, AsciiAlphabet())
    end
    if (width < 1) | (length(seq) ≤ width)
        # Fastest option
        return print(io, String(seq))
    else
        # Second fastest option
        buffer = SimpleBuffer(io, padded_length(length(seq), width))
        return _print(buffer, seq, width, AsciiAlphabet())
    end
end

function _print(buffer::SimpleBuffer, seq::LongSequence{<:Alphabet}, width::Integer, ::AsciiAlphabet)
    col = 0
    @inbounds for i in eachindex(seq)
        col += 1
        if (width > 0) & (col > width)
            write(buffer, UInt8('\n'))
            col = 1
        end
        write(buffer, stringbyte(seq[i]))
    end
    close(buffer)
    return nothing
end


text:
#Defines the Navigation system type, including its discrete time behaviour
#include("DiscreteTimeSystem.jl")
#= NavigationFunction should output the "measured state" of the vehicle, based on the
ODE state (simulated measured state) through some filter.
=#

#Define ControlSystem as a subtype of the DiscreteTimeSystem
abstract type AbstractNavigationSystem <: DiscreteTimeSystem

########### Template function ###########################
mutable struct NavigationSystem{T<:Number} <: AbstractNavigationSystem
  Active::Bool
  Ts::T #sample time
  NavigationRelation::Function
  Tprev::T #storage for previous time it was updated, hence must be a mutable type
end
NavigationSystem() = NavigationSystem(false, 0.0 , x-> nothing, 0.0)
NavigationSystem(Ts::Number) = NavigationSystem(false, Float64(Ts), x -> (), 0.0)
NavigationSystem(Ts::Number, csfunc::Function) = NavigationSystem(true, Float64(Ts), csfunc, 0.0)
NavigationSystem(Active::Bool, Ts::Number, csfunc::Function) = NavigationSystem(Active, Float64(Ts), csfunc, 0.0)

#Set the callback affect condition of the ControlSystem
#placeholder!
function affect!(int,dts::AbstractNavigationSystem)
    for cu in user_cache(int)
        cu.u = [dts.NavigationRelation(int.u)...]
    end
  dts.Tprev = int.t
end


text:
using LightGraphs: topological_sort_by_dfs, outneighbors
using MetaGraphs: get_prop

"convert literal+/- to probability value 0/1"
@inline lit2value(l::Lit)::Int = (l > 0 ? 1 : 0)

"""
Learning from data a circuit with several structure learning algorithms
"""
function learn_credal_circuit(data::Union{XData, WXData},s_idm::Float64;
        pseudocount = 1.0, algo = "chow-liu", algo_kwargs=(α=1.0, clt_root="graph_center"))::CredalΔ
    if algo == "chow-liu"
        clt = learn_chow_liu_tree(data; algo_kwargs...)
        pc = compile_cred_circuit_from_clt(clt)
        estimate_credal_parameters(pc, convert(XBatches,data),s_idm; pseudocount = pseudocount)
        pc
    else
        error("Cannot learn a probabilistic circuit with algorithm $algo")
    end
end

"""
Learning from data a circuit with several structure learning algorithms
"""
function learn_probabilistic_circuit(data::Union{XData, WXData};
        pseudocount = 1.0, algo = "chow-liu", algo_kwargs=(α=1.0, clt_root="graph_center"))::ProbΔ
    if algo == "chow-liu"
        clt = learn_chow_liu_tree(data; algo_kwargs...)
        pc = compile_prob_circuit_from_clt(clt)
        estimate_parameters(pc, convert(XBatches,data); pseudocount = pseudocount)
        pc
    else
        error("Cannot learn a probabilistic circuit with algorithm $algo")
    end
end

"Build decomposable probability circuits from Chow-Liu tree"
function compile_prob_circuit_from_clt(clt::CLT)::ProbΔ
    topo_order = Var.(reverse(topological_sort_by_dfs(clt::CLT))) #order to parse the node
    lin = Vector{ProbΔNode}()
    node_cache = Dict{Lit, LogicalΔNode}()
    prob_cache = ProbCache()
    parent = parent_vector(clt)

    prob_children(n)::Vector{<:ProbΔNode{<:node_type(n)}} =
        copy_with_eltype(map(c -> prob_cache[c], n.children), ProbΔNode{<:node_type(n)})

    "default order of circuit node, from left to right: +/1 -/0"

    "compile leaf node into circuits"
    function compile_leaf(ln::Var)
        pos = LiteralNode( var2lit(ln))
        neg = LiteralNode(-var2lit(ln))
        node_cache[var2lit(ln)] = pos
        node_cache[-var2lit(ln)] = neg
        pos2 = ProbLiteral(pos)
        neg2 = ProbLiteral(neg)
        push!(lin, pos2)
        push!(lin, neg2)
        prob_cache[pos] = pos2
        prob_cache[neg] = neg2
    end

    "compile inner disjunction node"
    function compile_⋁inner(ln::Lit, children::Vector{Var})::Vector{⋁Node}
        logical_nodes = Vector{⋁Node}()
        v = lit2value(ln)

        for c in children
            #build logical ciruits
            temp = ⋁Node([node_cache[lit] for lit in [var2lit(c), - var2lit(c)]])
            push!(logical_nodes, temp)
            n = Prob⋁(temp, prob_children(temp))
            prob_cache[temp] = n
            n.log_thetas = zeros(Float64, 2)
            cpt = get_prop(clt, c, :cpt)
            weights = [cpt[(1, v)], cpt[(0, v)]]
            n.log_thetas = log.(weights)
            push!(lin, n)
        end

        return logical_nodes
    end

    "compile inner conjunction node into circuits, left node is indicator, rest nodes are disjunction children nodes"
    function compile_⋀inner(indicator::Lit, children::Vector{⋁Node})
        leaf = node_cache[indicator]
        temp = ⋀Node(vcat([leaf], children))
        node_cache[indicator] = temp
        n = Prob⋀(temp, prob_children(temp))
        prob_cache[temp] = n
        push!(lin, n)
    end

    "compile inner node, 1 inner variable to 2 leaf nodes, 2 * num_children disjunction nodes and 2 conjunction nodes"
    function compile_inner(ln::Var, children::Vector{Var})
        compile_leaf(ln)
        pos⋁ = compile_⋁inner(var2lit(ln), children)
        neg⋁ = compile_⋁inner(-var2lit(ln), children)
        compile_⋀inner(var2lit(ln), pos⋁)
        compile_⋀inner(-var2lit(ln), neg⋁)
    end

    "compile root, add another disjunction node"
    function compile_root(root::Var)
        temp = ⋁Node([node_cache[s] for s in [var2lit(root), -var2lit(root)]])
        n = Prob⋁(temp, prob_children(temp))
        prob_cache[temp] = n
        n.log_thetas = zeros(Float64, 2)
        cpt = get_prop(clt, root, :cpt)
        weights = [cpt[1], cpt[0]]
        n.log_thetas = log.(weights)
        push!(lin, n)
        return n
    end

    function compile_independent_roots(roots::Vector{ProbΔNode})
        temp = ⋀Node([c.origin for c in roots])
        n = Prob⋀(temp, prob_children(temp))
        prob_cache[temp] = n
        push!(lin, n)
        temp = ⋁Node([temp])
        n = Prob⋁{LogicalΔNode}(temp, prob_children(temp))
        prob_cache[temp] = n
        n.log_thetas = [0.0]
        push!(lin, n)
    end

    roots = Vector{ProbΔNode}()
    for id in topo_order
        children = Var.(outneighbors(clt, id))
        if isequal(children, [])
            compile_leaf(id)
        else
            compile_inner(id, children)
        end
        if 0 == parent[id]
            push!(roots, compile_root(id))
        end
    end

    if length(roots) > 1
        compile_independent_roots(roots)
    end

    return lin
end

"Build decomposable credal circuit from Chow-Liu tree"
function compile_cred_circuit_from_clt(clt::CLT)::CredalΔ
    topo_order = Var.(reverse(topological_sort_by_dfs(clt::CLT))) #order to parse the node
    lin = Vector{CredalΔNode}()
    node_cache = Dict{Lit, LogicalΔNode}()
    cred_cache = CredalCache()
    parent = parent_vector(clt)

    cred_children(n)::Vector{<:CredalΔNode{<:node_type(n)}} =
        copy_with_eltype(map(c -> cred_cache[c], n.children), CredalΔNode{<:node_type(n)})

    "default order of circuit node, from left to right: +/1 -/0"

    "compile leaf node into circuits"
    function compile_leaf(ln::Var)
        pos = LiteralNode( var2lit(ln))
        neg = LiteralNode(-var2lit(ln))
        node_cache[var2lit(ln)] = pos
        node_cache[-var2lit(ln)] = neg
        pos2 = CredalLiteral(pos)
        neg2 = CredalLiteral(neg)
        push!(lin, pos2)
        push!(lin, neg2)
        cred_cache[pos] = pos2
        cred_cache[neg] = neg2
    end

    "compile inner disjunction node"
    function compile_⋁inner(ln::Lit, children::Vector{Var})::Vector{⋁Node}
        logical_nodes = Vector{⋁Node}()
        v = lit2value(ln)

        for c in children
            #build logical ciruits
            temp = ⋁Node([node_cache[lit] for lit in [var2lit(c), - var2lit(c)]])
            push!(logical_nodes, temp)
            n = Credal⋁(temp, cred_children(temp))
            cred_cache[temp] = n
            n.log_thetas = zeros(Float64, 2)
            cpt = get_prop(clt, c, :cpt)
            weights = [cpt[(1, v)], cpt[(0, v)]]
            n.log_thetas = log.(weights)
            push!(lin, n)
        end

        return logical_nodes
    end

    "compile inner conjunction node into circuits, left node is indicator, rest nodes are disjunction children nodes"
    function compile_⋀inner(indicator::Lit, children::Vector{⋁Node})
        leaf = node_cache[indicator]
        temp = ⋀Node(vcat([leaf], children))
        node_cache[indicator] = temp
        n = Credal⋀(temp, cred_children(temp))
        cred_cache[temp] = n
        push!(lin, n)
    end

    "compile inner node, 1 inner variable to 2 leaf nodes, 2 * num_children disjunction nodes and 2 conjunction nodes"
    function compile_inner(ln::Var, children::Vector{Var})
        compile_leaf(ln)
        pos⋁ = compile_⋁inner(var2lit(ln), children)
        neg⋁ = compile_⋁inner(-var2lit(ln), children)
        compile_⋀inner(var2lit(ln), pos⋁)
        compile_⋀inner(-var2lit(ln), neg⋁)
    end

    "compile root, add another disjunction node"
    function compile_root(root::Var)
        temp = ⋁Node([node_cache[s] for s in [var2lit(root), -var2lit(root)]])
        n = Credal⋁(temp, cred_children(temp))
        cred_cache[temp] = n
        n.log_thetas = zeros(Float64, 2)
        cpt = get_prop(clt, root, :cpt)
        weights = [cpt[1], cpt[0]]
        n.log_thetas = log.(weights)
        push!(lin, n)
        return n
    end

    function compile_independent_roots(roots::Vector{CredalΔNode})
        temp = ⋀Node([c.origin for c in roots])
        n = Credal⋀(temp, cred_children(temp))
        cred_cache[temp] = n
        push!(lin, n)
        temp = ⋁Node([temp])
        n = Credal⋁{LogicalΔNode}(temp, cred_children(temp))
        cred_cache[temp] = n
        n.log_thetas = [0.0]
        push!(lin, n)
    end

    roots = Vector{CredalΔNode}()
    for id in topo_order
        children = Var.(outneighbors(clt, id))
        if isequal(children, [])
            compile_leaf(id)
        else
            compile_inner(id, children)
        end
        if 0 == parent[id]
            push!(roots, compile_root(id))
        end
    end

    if length(roots) > 1
        compile_independent_roots(roots)
    end

    return lin
end


text:
B_data = [ -2.5   0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0          -2.0         1.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
1.0  -3.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0          -2.0         1.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.5           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           2.0         0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.5          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         2.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.5           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           2.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.5           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   2.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.5          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   2.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.5           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   2.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.5    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   2.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0    0.5   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   2.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.5    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.5           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           2.0         0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.5          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         2.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.5           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0           0.0         0.0           2.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.5           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   2.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.5          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   2.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.5           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0           0.0         0.0           0.0   0.0   0.0   2.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.5    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   2.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0    0.5   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   2.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.5    0.0   0.0
0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -5.0           5.0        -2.0           0.0   0.0   0.0   0.0   0.0   0.0   1.0           0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5          -2.0         0.5           0.0   0.0   0.0   0.0   0.0   0.0   0.0           1.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5          -2.0        -3.5           4.5  -3.0   1.5   0.0   0.0   0.0   0.0           0.0          1.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0          -5.0   5.0  -2.0   0.0   0.0   0.0   0.0           0.0          0.0           1.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0          -0.5  -2.0   0.5   0.0   0.0   0.0   0.0           0.0          0.0           0.0           1.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.5  -2.0  -3.5   4.5  -3.0   1.5   0.0           0.0          0.0           0.0           0.0          1.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0  -5.0   5.0  -2.0   0.0           0.0          0.0           0.0           0.0          0.0           1.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0  -0.5  -2.0   0.5   0.0           0.0          0.0           0.0           0.0          0.0           0.0    1.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.5  -2.0  -3.5   0.0           0.0          0.0           0.0           0.0          0.0           1.5   -3.0   5.5    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   8.75         -4.5         2.25          0.0   0.0   0.0   0.0   0.0   0.0  -5.25         -3.0          0.75          0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.81025e-16   2.0         9.36751e-17   0.0   0.0   0.0   0.0   0.0   0.0   0.75         -3.0         -0.75          0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
-0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0  -1.8735e-15    1.249e-15   2.0           0.0   0.0   0.0   0.0   0.0   0.0  -3.0           7.5         -7.5           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           2.0   0.0   0.0   0.0   0.0   0.0   2.25         -4.5          6.75         -5.25         -3.0          0.75          0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   2.0   0.0   0.0   0.0   0.0   9.36751e-17  -1.8735e-16   2.81025e-16   0.75         -3.0         -0.75          0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   2.0   0.0   0.0   0.0  -6.245e-16     1.249e-15   -1.8735e-15   -3.0           7.5         -7.5           0.0    0.0   0.0    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   2.0   0.0   0.0   0.0           0.0          0.0           2.25         -4.5          6.75         -5.25  -3.0   0.75   0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   2.0   0.0   0.0           0.0          0.0           9.36751e-17  -1.8735e-16   2.81025e-16   0.75  -3.0  -0.75   0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   2.0   0.0           0.0          0.0          -6.245e-16     1.249e-15   -1.8735e-15   -3.0    7.5  -7.5    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.25  -0.5   0.75  -2.5   2.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0         0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.5   -1.0   1.5    1.0  -2.0]


text:
##############################################################################
##
## LSMR GPU
##
##############################################################################
using CuArrays.CUDAnative
import CuArrays: allowscalar
allowscalar(false)

##############################################################################
##
## Conversion FixedEffect between CPU and GPU
##
##############################################################################

# https://github.com/JuliaGPU/CuArrays.jl/issues/306
cuzeros(T::Type, n::Integer) = fill!(CuVector{T}(undef, n), zero(T))
function CuArrays.cu(T::Type, fe::FixedEffect)
	refs = CuArray(fe.refs)
	interaction = cu(T, fe.interaction)
	FixedEffect{typeof(refs), typeof(interaction)}(refs, interaction, fe.n)
end
CuArrays.cu(T::Type, w::Union{Fill, Ones, Zeros}) = fill!(CuVector{T}(undef, length(w)), w[1])
CuArrays.cu(T::Type, w::AbstractVector) = CuVector{T}(w)

##############################################################################
##
## FixedEffectLinearMap on the GPU (code by Paul Schrimpf)
##
##############################################################################

const N_THREADS = 256

function FixedEffectLinearMap{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Type{Val{:lsmr_gpu}}) where {T}
	fes = [cu(T, fe) for fe in fes]
	sqrtw = cu(T, sqrt.(values(weights)))
	colnorm = [colnorm!(cuzeros(T, fe.n), fe.refs, fe.interaction, sqrtw) for fe in fes]
	caches = [cache!(cuzeros(T, length(sqrtw)), fe.interaction, sqrtw, scale, fe.refs) for (fe, scale) in zip(fes, colnorm)]
	return FixedEffectLinearMap{T}(fes, sqrtw, colnorm, caches)
end

function colnorm!(fecoef::CuVector, refs::CuVector, y::CuVector, sqrtw::CuVector)
	nblocks = cld(length(refs), N_THREADS)
	@cuda threads=N_THREADS blocks=nblocks colnorm!_kernel!(fecoef, refs, y, sqrtw)
	fecoef .= sqrt.(fecoef)
end

function colnorm!_kernel!(fecoef, refs, y, sqrtw)
	index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
	stride = blockDim().x * gridDim().x
	@inbounds for i = index:stride:length(y)
		CuArrays.CUDAnative.atomic_add!(pointer(fecoef, refs[i]), abs2(y[i] * sqrtw[i]))
	end
end

function cache!(y::CuVector, interaction::CuVector , sqrtw::CuVector, fecoef::CuVector, refs::CuVector)
	nblocks = cld(length(y), N_THREADS)
	@cuda threads=N_THREADS blocks=nblocks cache_kernel!(y, interaction, sqrtw, fecoef, refs)
	return y
end

function cache_kernel!(y, interaction, sqrtw, fecoef, refs)
	index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
	stride = blockDim().x * gridDim().x
	@inbounds for i = index:stride:length(y)
		y[i] += interaction[i] * sqrtw[i] / fecoef[refs[i]]
	end
end

function mean!(fecoef::CuVector, refs::CuVector, α::Number, y::CuVector, cache::CuVector)
	nblocks = cld(length(y), N_THREADS)
	@cuda threads=N_THREADS blocks=nblocks mean_kernel!(fecoef, refs, α, y, cache)
end

function mean_kernel!(fecoef, refs, α, y, cache)
	index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
	stride = blockDim().x * gridDim().x
	@inbounds for i = index:stride:length(y)
		CuArrays.CUDAnative.atomic_add!(pointer(fecoef, refs[i]), α * y[i] * cache[i])
	end
end

function demean!(y::CuVector, α::Number, fecoef::CuVector, refs::CuVector, cache::CuVector)
	nblocks = cld(length(y), N_THREADS)
	@cuda threads=N_THREADS blocks=nblocks demean_kernel!(y, α, fecoef, refs, cache)
end

function demean_kernel!(y, α, fecoef, refs, cache)
	index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
	stride = blockDim().x * gridDim().x
	@inbounds for i = index:stride:length(y)
		y[i] += α * fecoef[refs[i]] * cache[i]
	end
end

##############################################################################
##
## Implement AbstractFixedEffectSolver interface
##
##############################################################################
struct FixedEffectSolverLSMRGPU{T} <: AbstractFixedEffectSolver{T}
	m::FixedEffectLinearMap{T}
	b::CuVector{T}
	r::CuVector{T}
	x::FixedEffectCoefficients{T}
	v::FixedEffectCoefficients{T}
	h::FixedEffectCoefficients{T}
	hbar::FixedEffectCoefficients{T}
	tmp::Vector{T} # used to convert AbstractVector to Vector{T}
	fes::Vector{<:FixedEffect}
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Type{Val{:lsmr_gpu}}) where {T}
	m = FixedEffectLinearMap{T}(fes, weights, Val{:lsmr_gpu})
	b = cuzeros(T, length(weights))
	r = cuzeros(T, length(weights))
	x = FixedEffectCoefficients([cuzeros(T, fe.n) for fe in fes])
	v = FixedEffectCoefficients([cuzeros(T, fe.n) for fe in fes])
	h = FixedEffectCoefficients([cuzeros(T, fe.n) for fe in fes])
	hbar = FixedEffectCoefficients([cuzeros(T, fe.n) for fe in fes])
	tmp = zeros(T, length(weights))
	FixedEffectSolverLSMRGPU(m, b, r, x, v, h, hbar, tmp, fes)
end

function solve_residuals!(r::AbstractVector, feM::FixedEffectSolverLSMRGPU{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
	copyto!(feM.tmp, r)
	copyto!(feM.r, feM.tmp)
	feM.r .*=  feM.m.sqrtw
	copyto!(feM.b, feM.r)
	fill!(feM.x, 0.0)
	x, ch = lsmr!(feM.x, feM.m, feM.b, feM.v, feM.h, feM.hbar; atol = tol, btol = tol, maxiter = maxiter)
	mul!(feM.r, feM.m, feM.x, -1.0, 1.0)
	feM.r ./=  feM.m.sqrtw
	copyto!(feM.tmp, feM.r)
	copyto!(r, feM.tmp)
	return r, div(ch.mvps, 2), ch.isconverged
end

function solve_residuals!(X::AbstractMatrix, feM::FixedEffectSolverLSMRGPU; kwargs...)
	iterations = Int[]
	convergeds = Bool[]
	for j in 1:size(X, 2)
		_, iteration, converged = solve_residuals!(view(X, :, j), feM; kwargs...)
		push!(iterations, iteration)
		push!(convergeds, converged)
	end
	return X, iterations, convergeds
end

function solve_coefficients!(r::AbstractVector, feM::FixedEffectSolverLSMRGPU{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
	copyto!(feM.tmp, r)
	copyto!(feM.b, feM.tmp)
	feM.b .*= feM.m.sqrtw
	fill!(feM.x, 0.0)
	x, ch = lsmr!(feM.x, feM.m, feM.b, feM.v, feM.h, feM.hbar; atol = tol, btol = tol, maxiter = maxiter)
	for (x, scale) in zip(feM.x.x, feM.m.colnorm)
		x ./=  scale
	end
	x = Vector{eltype(r)}[collect(x) for x in feM.x.x]
	full(normalize!(x, feM.fes; tol = tol, maxiter = maxiter), feM.fes), div(ch.mvps, 2), ch.isconverged
end


text:
mutable struct GravTreeGather{N,T}
	acc::SVector{N,T}
	pot::T
    mass::T
    #nodecenters::Vector{Vector{T}} #for debugging/visualization
    #nodelengths::Vector{Vector{T}} #for debugging/visualization
end

GravTreeGather{N,T}() where {N,T} = GravTreeGather{N,T}(zero(SVector{N,T}),0,0)

#when changing the argument of treewalk, remember to also change the recursive treewalk call inside the function body!
function gravity_treewalk!(ga::GravTreeGather{N,T}, p::SVector{N,T}, node::Node{N,T,D},
	openingangle::T, softening::T, boxsizes::SVector{N,T}) where {N,T,D}
    if isLeaf(node)
        #println("in a leaf node")
        if node.p != nothing
            #println("nonempty leaf")
            ga.mass += node.p.mass
            dx = nearest.(node.p.pos - p, boxsizes)
			r2 = sum(dx.^2)
			#@show r2, node.p.pos, p
			if r2 > 0.0 #exclude self contribution (r2=0)
				r_inv = sqrt(1.0 / (r2 + softening^2))
				r3_inv = r_inv^3
				ga.acc += node.p.mass * r3_inv * dx
				ga.pot += node.p.mass * r_inv
	            #push!(ga.nodecenters, node.p.pos)
	            #push!(ga.nodelengths, zeros(N))
			end
            #@show "it's a particle..." node.center, node.length
            #r2 = get_distance2(node.part, p, boxsizes, true)
        end
    else
        #println("This is a node... check the opening angle")
		#r2 = get_distance2(pos_c, p, boxsizes, true)
		#don't use get_distance2 as we need both r2 and dx (for vec2pix)
		dx = nearest.(node.n.pos - p, boxsizes)
		r2 = sum(dx.^2)
		if r2 > (node.length[1] / openingangle)^2
			#println("skip node ", i)
			ga.mass += node.n.mass
			r_inv = sqrt(1.0 / (r2 + softening^2))
			r3_inv = r_inv^3
			ga.acc += node.n.mass * r3_inv * dx
			ga.pot += node.n.mass * r_inv
			#push!(ga.nodecenters, node.center)
			#push!(ga.nodelengths, node.length)
			#@show "use this node!" node.center, node.length
		else
	        @inbounds for i in 1:2^N
	            #println("open this node")
	            gravity_treewalk!(ga, p, node.child[i], openingangle, softening, boxsizes)
	        end
		end
    end
end

function gravity_bruteforce!(ga::GravTreeGather{N,T}, p::SVector{N,T}, X::Vector{SVector{N,T}},
							 mass::Vector{T}, boxsizes::SVector{N,T}) where {N,T}
	for i in eachindex(X)
		ga.mass += mass[i]
		dx = nearest.(X[i] - p, boxsizes)
		r2 = sum(dx.^2)
		if r2 > 0.0 #exclude self contribution (r2=0)
			r_inv = sqrt(1.0 / r2)
			r3_inv = r_inv^3
			ga.acc += mass[i] * r3_inv * dx
			ga.pot += mass[i] * r_inv
		end
	end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -1.263105417837706 <= q <= 2.4766737340949954)

                     add_NL_constraint(m, :(softplus(0.7635571507341847 + -0.35250618093605723*softplus(0.5051094873509578 + -0.17934163301535788*softplus(-0.331542222824047 + 0.7450661875335634*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9610191438924418*softplus(0.6781688636512797 + -0.22580923517577522*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))) + 0.1709251594821839*softplus(0.058205314222445015 + 0.7954955842799771*softplus(-0.331542222824047 + 0.7450661875335634*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9762293687340509*softplus(0.6781688636512797 + -0.22580923517577522*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))) + softplus(-0.29663194049169395 + 0.013502946151282558*softplus(0.5051094873509578 + -0.17934163301535788*softplus(-0.331542222824047 + 0.7450661875335634*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9610191438924418*softplus(0.6781688636512797 + -0.22580923517577522*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))) + -0.2459643040595232*softplus(0.058205314222445015 + 0.7954955842799771*softplus(-0.331542222824047 + 0.7450661875335634*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9762293687340509*softplus(0.6781688636512797 + -0.22580923517577522*softplus(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*softplus(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions