text:
module CSparse

    using Metis

    export testmat

    include("utilities.jl")
    include("testmat.jl")
    include("simplicialchol.jl")
    include("trisolvers.jl")

## Replace calls to these simple functions by negation, abs, etc.
JS_FLIP(i::Integer) = -i    # CSparse must flip about -1 because 0 is a valid index
JS_UNFLIP(i::Integer) = (i <= 0) ? -i : i        # or abs(i)
JS_MARKED(w::Vector{Integer},j::Integer) = w[j] <= 0
JS_MARK(w::Vector{Integer},j::Integer) = w[j] = -w[j]

function js_chol{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, S::CholSymb{Ti})
    n = size(A,2)
    cp = S.cp; pinv = S.pinv; parent = S.parent
## maybe change the natural case (i.e. no permutation) to pinv = Array(Ti, 0)
    C  = (pinv == 1:n) ? A : js_symperm(A, pinv)
    Cp = C.colptr; Ci = C.rowval; Cx = C.nzval
    Lp = copy(cp); Li = Array(Ti, cp[n+1]); Lx = Array(Tv, cp[n+1])
    c  = copy(cp[1:n])
    s  = Array(Ti, n)                   # Ti workspace
    x  = Array(Tv, n)                   # Tv workspace
    for k in 1:n                        # compute L[k,:] for L*L' = C
        ## Nonzero pattern of L[k,:]
        top = cs_ereach(C, k, parent, s, c) # find pattern of L(k,:)
        x[k] = 0                            # x[1:k] is now zero
        for p in Cp[k]:(Cp[k+1]-1)          # x = full(triu(C(:,k)))
            if (Ci[p] <= k) x[Ci[p]] = Cx[p] end
        end
        d = x[k]                        # d = C[k,k]
        x[k] = zero(Tv)                 # clear x for k+1st iteration
        ## Triangular solve
        while top <= n             # solve L[0:k-1,0:k-1] * x = C[:,k]
            i = s[top]             # s[top..n-1] is pattern of L[k,:]
            lki = x[i]/Lx[Lp[i]]   # L[k,i] = x[i]/L[i,i]
            x[i] = zero(Tv)        # clear x for k+1st iteration
            for p in (Lp[i]+1):(c[i]-1)
                x[Li[p]] -= Lx[p] * lki
            end
            d -= lki * lki ;            # d = d - L[k,i]*L[k,i]
            p = c[i]
            c[i] += 1
            Li[p] = k                   # store L[k,i] in column i
            Lx[p] = lki
        end
        ## Compute L[k,k]
        if (d <= 0) error("Matrix is not positive definite, detected at row $k") end
        p = c[k]
        c[k] += 1
        Li[p] = k                 # store L[k,k] = sqrt(d) in column k
        Lx[p] = sqrt(d)
    end
    SparseMatrixCSC(n, n, Lp, Li, Lx)
end

## depth-first-search of the graph of a matrix, starting at node j
function js_dfs{T}(j::Integer, G::SparseMatrixCSC{T}, top::Integer,
                   xi::Vector{Integer}, pstack::Vector{Integer},
                   pinv::Vector{Integer})
    head = 0; Gp = G.colptr; Gi = G.rowval
    xi[0] = j                        # initialize the recursion stack
    while (head >= 0)
        j = xi[head]       # get j from the top of the recursion stack
        jnew = pinv[j]
        if !JS_MARKED(Gp, j)
            JS_MARK(Gp, j)              # mark node j as visited
            pstack[head] = (jnew < 0) ? 0 : JS_UNFLIP(Gp[jnew])
        end
        done = 1            # node j done if no unvisited neighbors
        p2 = (jnew < 0) ? 0 : JS_UNFLIP(Gp[jnew+1])
        for p in pstack[head]:(p2-1)    # examine all neighbors of j
            i = Gi[p]                   # consider neighbor node i
            if (JS_MARKED(Gp, i)) continue end # skip visited node i
            pstack[head] = p      # pause depth-first search of node j
            xi[++head] = i        # start dfs at node i
            done = 0              # node j is not done
            break                 # break, to start dfs[i]
        end
        if done                 # depth-first search at node j is done
            head -= 1           # remove j from the recursion stack
            top -= 1
            xi[top] = j                # and place in the output stack
        end
    end
    top
end

# compute the etree of A (using triu(A), or A'A without forming A'A
# A root of the etree (which may be a forest) is indicated by 0
function js_etree{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, ata::Bool)
    m,n = size(A); Ap = A.colptr; Ai = A.rowval
    parent = zeros(Ti, n); w = zeros(Ti, n + (ata ? m : 0))
    ancestor = 0; prev = n              # offsets into w
    if (ata) w[prev + (1:m)] = 0 end
    for k in 1:n
        parent[k] = 0                   # node k has no parent yet
        w[ancestor + k] = 0             # nor does k have an ancestor
        for p in Ap[k]:(Ap[k+1] - 1)
            i = ata ? w[Ai[p] + prev] : Ai[p]
            while i != 0 && i < k
                inext = w[ancestor + i] # inext = ancestor of i
                w[ancestor + i] = k     # path compression
                if (inext == 0) parent[i] = k end # no anc., parent is k
                i = inext
            end
            if (ata) w[Ai[p] + prev] = k end
        end
    end
    parent
end

js_etree(A::SparseMatrixCSC) = js_etree(A, false)

## post order a forest
function js_post{T<:Union(Int32,Int64)}(parent::Vector{T})
    n = length(parent)
    head = zeros(T,n)                   # empty linked lists
    next = zeros(T,n)
    for j in n:-1:1                  # traverse nodes in reverse order
        if (parent[j] == 0) continue end # j is a root
        next[j] = head[parent[j]]      # add j to list of its parent
        head[parent[j]] = j
    end
    return head, next
    stack = zeros(T,n)
    post = zeros(T,n)
    k = 1
    for j in 1:n
        if (parent[j] != 0) continue end # skip j if it is not a root
        k = js_tdfs(j, k, head, next, post, stack) ;
    end
    post
end

## xi[top:n] = nodes reachable from graph of G*P' via nodes in B[:,k]
## xi[n:2n] used as workspace
function js_reach{Tv,Ti}(G::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti},
                         k::Integer, xi::Vector{Ti}, pinv::Vector{Ti})
    n = size(G,2); Bp = B.colptr; Bi = B.rowval; Gp = G.colptr
    top = n
    for p in Bp[k]:(Bp[k+1]-1)
        if !JS_MARKED(Gp, Bi[p])      # start a dfs at unmarked node i
            top = js_dfs(Bi[p], G, top, xi, xi+n, pinv)
        end
    end
    for p in top:n JS_MARK(Gp, xi[p]) end  #restore G
    top
end

# depth-first search and postorder of a tree rooted at node j
function js_tdfs{T<:Union(Int32,Int64)}(j::Integer, k::Integer, head::Vector{T},
                 next::Vector{T}, post::Vector{T}, stack::Vector{T})
    top = 1
    stack[1] = j                        # place j on the stack
    while (top > 0)                     # while (stack is not empty)
        p = stack[top]                  # p = top of stack
        i = head[p]                     # i = youngest child of p
        if (i == 0)
            top -= 1              # p has no unordered children left
            post[k] = p           # node p is the kth postordered node
            p += 1
        else
            head[p] = next[i]           # remove i from children of p
            top += 1
            stack[top] = i              # start dfs on child node i
        end
    end
    k
end

## This is essentially a copy of the transpose function in sparse.jl
function js_transpose{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti})
    m,n = size(A)
    Ap = A.colptr; Ai = A.rowval; Ci = similar(Ai); Ax = A.nzval; Cx = similar(Ax)
    w  = zeros(Ti, m + 1)
    w[1] = one(Ti)
    for i in Ai w[i + 1] += 1 end         # row counts
    Cp = cumsum(w)
    w[:] = Cp[:]
    for j in 1:n, p in Ap[j]:(Ap[j+1] - 1)
        i = Ai[p]; q = w[i]; w[i] += 1
        Ci[q] = j
        Cx[q] = Ax[p]
    end
    SparseMatrixCSC(n, m, Cp, Ci, Cx)
end

## Copied from extras/suitesparse.jl
function _jl_convert_to_0_based_indexing!(S::SparseMatrixCSC)
    for i in 1:length(S.rowval) S.rowval[i] -= 1 end
    for p in 1:length(S.colptr) S.colptr[p] -= 1 end
    S
end

function _jl_convert_to_1_based_indexing!(S::SparseMatrixCSC)
    for i in 1:length(S.rowval) S.rowval[i] += 1 end
    for p in 1:length(S.colptr) S.colptr[p] += 1 end
    S
end

_jl_convert_to_0_based_indexing(S) = _jl_convert_to_0_based_indexing!(copy(S))
_jl_convert_to_1_based_indexing(S) = _jl_convert_to_1_based_indexing!(copy(S))

type cs{Tv<:Union(Float64,Complex128),Ti<:Union(Int32,Int64)} # the CXSparse cs struct
    nzmax::Ti
    m::Ti
    n::Ti
    p::Ptr{Ti}
    i::Ptr{Ti}
    x::Ptr{Tv}
    nz::Ti
end

function cs{Tv<:Union(Float64,Complex128),Ti<:Union(Int32,Int64)}(A::SparseMatrixCSC{Tv,Ti})
    if A.colptr[1] != 0 error("Sparse matrix must be in 0-based indexing") end
    cs{Tv,Ti}(convert(Ti,A.colptr[end]), convert(Ti,A.m), convert(Ti,A.n),
              pointer(A.colptr), pointer(A.rowval), pointer(A.nzval), convert(Ti, -1))
end

## Primary CXSparse-based functions, cs_cholsol, cs_lusol, cs_print, cs_qrsol
for (cholsol, lusol, prt, qrsol, vtyp, ityp) in
    ((:cs_ci_cholsol, :cs_ci_lusol, :cs_ci_print, :cs_ci_qrsol, :Complex128, :Int32),
     (:cs_di_cholsol, :cs_di_lusol, :cs_di_print, :cs_di_qrsol, :Float64, :Int32),
     (:cs_cl_cholsol, :cs_cl_lusol, :cs_cl_print, :cs_cl_qrsol, :Complex128, :Int64),
     (:cs_dl_cholsol, :cs_dl_lusol, :cs_dl_print, :cs_dl_qrsol, :Float64, :Int64))
    @eval begin
        ## A is symmetric but only triu(A) is passed
        function cs_cholsol!(A::SparseMatrixCSC{$vtyp,$ityp}, b::Vector{$vtyp}, order::Integer)
            if !(0 <= order <= 3) error("order = $order is not in the range 0 to 3") end
            m,n = size(A)
            if m != n || n != length(b) error("Dimension mismatch") end
            st = ccall(($(string(cholsol)),"libcxsparse"), $ityp, ($ityp, Ptr{Void}, Ptr{$vtyp}),
                       order, pack(cs(_jl_convert_to_0_based_indexing!(A))).data, b)
            _jl_convert_to_1_based_indexing!(A)
            if st == 0 error("Failure in cholsol") end
            b
        end

        ## A must be square
        function cs_lusol!(A::SparseMatrixCSC{$vtyp,$ityp}, b::Vector{$vtyp}, order::Integer, tol::Real)
            if !(0 <= order <= 3) error("order = $order is not in the range 0 to 3") end
            m,n = size(A)
            if m != n || n != length(b) error("Dimension mismatch") end
            st = ccall(($(string(lusol)),"libcxsparse"), $ityp, ($ityp, Ptr{Void}, Ptr{$vtyp}, Float64),
                       order, pack(cs(_jl_convert_to_0_based_indexing!(A))).data, b, tol)
            _jl_convert_to_1_based_indexing!(A)
            if st == 0 error("Failure in lusol") end
            b
        end

        function cs_print(A::SparseMatrixCSC{$vtyp,$ityp}, brief::Bool)
            ccall(($(string(prt)),"libcxsparse"), Void, (Ptr{Void}, $ityp),
                  pack(cs(_jl_convert_to_0_based_indexing!(A))).data, brief)
            _jl_convert_to_1_based_indexing!(A)
            None
        end

        ## Overdetermined least squares system when m >= n
        function cs_qrsol!(A::SparseMatrixCSC{$vtyp,$ityp}, b::Vector{$vtyp}, order::Integer)
            if (order != 0 && order != 3) error("order = $order but must be 0 or 3") end
            if size(A, 1) != length(b) error("Dimension mismatch") end
            st = ccall(($(string(qrsol)),"libcxsparse"), $ityp, ($ityp, Ptr{Void}, Ptr{$vtyp}),
                       order, pack(cs(_jl_convert_to_0_based_indexing!(A))).data, b)
            _jl_convert_to_1_based_indexing!(A)
            if st == 0 error("Failure in qrsol") end
            b
        end
    end
end

## default ordering for cs_cholsol is 1 (i.e. amd(A))
cs_cholsol!{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T}) = cs_cholsol!(A, b, 1)
function cs_cholsol{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T}, ord::Integer)
    cs_cholsol!(A, copy(b), ord)
end
function cs_cholsol{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T})
    cs_cholsol(A, copy(b), 1)
end

## default ordering for cs_lusol is 2.
## There are other cases with order=1 but I haven't yet deciphered them
cs_lusol!{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T}) = cs_lusol!(A, b, 2, 1.)
function cs_lusol{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T}, ord::Integer)
    cs_lusol!(A, copy(b), ord, ord == 1 ? 0.001 : 1.)
end
function cs_lusol{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T})
    cs_lusol!(A, copy(b), 2, 1.)
end

cs_qrsol!{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T}) = cs_lusol!(A, b, 3)
function cs_qrsol{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T}, ord::Integer)
    cs_qrsol!(A, copy(b), ord)
end
function cs_qrsol{T<:Union(Float64,Complex128)}(A::SparseMatrixCSC{T}, b::Vector{T})
    cs_qrsol(A, copy(b), 3)
end

for (amd, etree, post, counts, norm, vtyp, ityp) in
    ((:cs_ci_amd, :cs_ci_etree, :cs_ci_post, :cs_ci_counts, :cs_ci_norm, :Complex128, :Int32),
     (:cs_di_amd, :cs_di_etree, :cs_di_post, :cs_di_counts, :cs_di_norm, :Float64, :Int32),
     (:cs_cl_amd, :cs_cl_etree, :cs_cl_post, :cs_cl_counts, :cs_cl_norm, :Complex128, :Int64),
     (:cs_dl_amd, :cs_dl_etree, :cs_dl_post, :cs_dl_counts, :cs_dl_norm, :Float64, :Int64))
    @eval begin
        ## Approximate minimal degree ordering
        function cs_amd(A::SparseMatrixCSC{$vtyp,$ityp}, order::Integer)
            if !(0 < order < 4) error("Valid values of order are 1:Chol, 2:LU, 3:QR") end
            ppt   = ccall(($(string(amd)),"libcxsparse"), Ptr{$ityp}, ($ityp, Ptr{Void}),
                          order, pack(cs(_jl_convert_to_0_based_indexing!(A))).data)
            _jl_convert_to_1_based_indexing!(A)
            pointer_to_array(ppt, (size(A,2),)) + 1
        end

        ## cs_counts returns all the information from cs_etree plus the column counts
        function cs_counts(A::SparseMatrixCSC{$vtyp,$ityp}, col::Bool)
            n = size(A, 2)
            cspk = pack(cs(_jl_convert_to_0_based_indexing!(A)))
            etrpt = ccall(($(string(etree)),"libcxsparse"), Ptr{$ityp},
                          (Ptr{Void}, $ityp), cspk.data, col)
            pospt = ccall(($(string(post)),"libcxsparse"), Ptr{$ityp},
                          (Ptr{$ityp}, $ityp), etrpt, n)
            coupt = ccall(($(string(counts)),"libcxsparse"), Ptr{$ityp},
                          (Ptr{Void}, Ptr{$ityp}, Ptr{$ityp}, $ityp),
                          cspk.data, etrpt, pospt, col)
            _jl_convert_to_1_based_indexing!(A)
            (pointer_to_array(etrpt, (n,)) + 1, pointer_to_array(pospt, (n,)) + 1,
             pointer_to_array(coupt, (n,)))
        end

        ## returns the elimination tree and the post-ordering permutation
        function cs_etree(A::SparseMatrixCSC{$vtyp,$ityp}, col::Bool)
            ept = ccall(($(string(etree)),"libcxsparse"), Ptr{$ityp}, (Ptr{Void}, $ityp),
                        pack(cs(_jl_convert_to_0_based_indexing!(A))).data, col)
            _jl_convert_to_1_based_indexing!(A)
            n = size(A, 2)
            popt  = ccall(($(string(post)),"libcxsparse"), Ptr{$ityp},
                          (Ptr{$ityp}, $ityp), ept, n)
            pointer_to_array(ept, (n,)) + 1, pointer_to_array(popt, (n,)) + 1
        end

        ## 1-norm of a sparse matrix (better to use js_norm, this is just for illustration)
        function cs_norm(A::SparseMatrixCSC{$vtyp,$ityp})
            res = ccall(($(string(norm)),"libcxsparse"), Float64, (Ptr{Void},),
                        pack(cs(_jl_convert_to_0_based_indexing!(A))).data)
            _jl_convert_to_1_based_indexing!(A)
            res
        end
    end
end

cs_amd(A::SparseMatrixCSC) = cs_amd(A, 1)
cs_counts(A::SparseMatrixCSC) = cs_counts(A, false)
cs_etree(A::SparseMatrixCSC) = cs_etree(A, false)
cs_print(A::SparseMatrixCSC) = cs_print(A, true)

type cs_symb{Tv<:Union(Float64,Complex128),Ti<:Union(Int32,Int64)} # the CXSparse cs_symbolic struct
    pinv::Ptr{Ti}
    q::Ptr{Ti}
    parent::Ptr{Ti}
    cp::Ptr{Ti}
    leftmost::Ptr{Ti}
    m2::Ti
    lnz::Float64
    unz::Float64
end

type cs_num{Tv<:Union(Float64,Complex128),Ti<:Union(Int32,Int64)} # the CXSparse cs_numeric struct
    L::Ptr{cs{Tv,Ti}}
    U::Ptr{cs{Tv,Ti}}
    pinv::Ptr{Ti}
    B::Ptr{Float64}
end

for (chol, qr, schol, sqr, symperm, vtyp, ityp) in
    ((:cs_ci_chol, :cs_ci_qr, :cs_ci_schol, :cs_ci_sqr, :cs_ci_symperm, :Complex128, :Int32),
     (:cs_di_chol, :cs_di_qr, :cs_di_schol, :cs_di_sqr, :cs_di_symperm, :Float64, :Int32),
     (:cs_cl_chol, :cs_cl_qr, :cs_cl_schol, :cs_cl_sqr, :cs_cl_symperm, :Complex128, :Int64),
     (:cs_dl_chol, :cs_dl_qr, :cs_dl_schol, :cs_dl_sqr, :cs_dl_symperm, :Float64, :Int64))
    @eval begin
        ## Numeric Cholesky factorization
        function cs_chol(A::SparseMatrixCSC{$vtyp,$ityp}, S::cs_symb{$vtyp,$ityp})
            pt  = ccall(($(string(chol)),"libcxsparse"), Ptr{Uint8}, (Ptr{Void}, Ptr{Void}),
                        pack(cs(_jl_convert_to_0_based_indexing!(A))).data, pack(S).data)
            _jl_convert_to_1_based_indexing!(A)
            outtyp = cs_num{$vtyp,$ityp}
            unpack(IOString(pointer_to_array(pt, (sum(map(sizeof, outtyp.types)),))), outtyp)
        end
        ## Symbolic Cholesky factorization
        function cs_schol(A::SparseMatrixCSC{$vtyp,$ityp}, order::Integer)
            if !(0 <= order <= 3) error("order = $order is not in the range 0 to 3") end
            pt  = ccall(($(string(schol)),"libcxsparse"), Ptr{Uint8}, ($ityp, Ptr{Void}),
                        order, pack(cs(_jl_convert_to_0_based_indexing!(A))).data)
            _jl_convert_to_1_based_indexing!(A)
            outtyp = cs_symb{$vtyp,$ityp}
            unpack(IOString(pointer_to_array(pt, (sum(map(sizeof, outtyp.types)),))), outtyp)
        end
        ## Symbolic QR decomposition
        function cs_sqr(A::SparseMatrixCSC{$vtyp,$ityp}, order::Integer)
            if !(0 <= order <= 3) error("order = $order is not in the range 0 to 3") end
            pt  = ccall(($(string(sqr)),"libcxsparse"), Ptr{Uint8}, ($ityp, Ptr{Void}, $ityp),
                        order, pack(cs(_jl_convert_to_0_based_indexing!(A))).data, 1)
            _jl_convert_to_1_based_indexing!(A)
            outtyp = cs_symb{$vtyp,$ityp}
            unpack(IOString(pointer_to_array(pt, (sum(map(sizeof, outtyp.types)),))), outtyp)
        end
    end
end

cs_schol(A::SparseMatrixCSC) = cs_schol(A, 1)
cs_sqr(A::SparseMatrixCSC) = cs_sqr(A, 3)

end


text:
module Interfaces

export ⊕, ⊗, ⋅
export add!, sub!, mul!, div!
export id, rank, dimension, kind, sequence, dims, index, inds
export update!, reset!, register!, prepare!, run!, expand, expand!, decompose, decompose!, regularize, regularize!
export permute, vector, matrix

"Direct sum."
function ⊕ end

"Direct product."
function ⊗ end

"Dot product."
function ⋅ end

"Inplace addition."
function add! end

"Inplace subtraction."
function sub! end

"Inplace multiplication."
function mul! end

"Inplace division."
function div! end

"Inplace Update."
function update! end

"Inplace reset."
function reset! end

"Inplace register."
function register! end

"Inplace prepare."
function prepare! end

"Inplace run."
function run! end

"Id."
function id end

"Rank."
function rank end

"Dimension."
function dimension end

"Kind."
function kind end

"Sequence."
function sequence end

"Dimensions."
function dims end

"Index."
function index end

"Indices."
function inds end

"Get the expansion."
function expand end

"In place expansion."
function expand! end

"Decompose."
function decompose end

"In place decomposition."
function decompose! end

"Regularize."
function regularize end

"In place regularization."
function regularize! end

"Get the permutation."
function permute end

"Vector representation."
function vector end

"Matrix representation."
function matrix end

end # module

ext: jl


text:
#!/usr/bin/env julia

using ArgParse

include("../src/findsl.jl")

function parse_commandline()
     s = ArgParseSettings("This program finds Synthetic Lethal pairs that knock out one or more genes from the essential gene database",
                         version = "0.0.1",
                         add_version = true)

    @add_arg_table s begin
        "--downregulated-cutoff"
            help = "This determines at which level the node is determined to be down regulated."
            arg_type = Float64
            default = 0.99
        "--upregulated-cutoff"
            help = "This determines at which level the node is determined to be upregulated."
            arg_type = Float64
            default = 1.02
        "--upperbound", "-u"
            arg_type = Int
            default = 10
        "--lowerbound", "-l"
            arg_type = Float64
            default = 0.001
        "--tissue-type"
            help = "This will correspond to a tissue type in the file data/expression-data.tsv"
            arg_type = AbstractString
            default = "esophagus"
        "--verbose", "-v"
            help = "This will cause output to be printed to standard out."
            action = :store_true
        "pairwise-interaction-file"
            help = "This is the full path to the pairwise interaction file."
            required = true
        "sl-results-file"
            help = "This file will be generated by this script."
            required = true
        "db-id-file"
            help = "This is the full path to the db_id_to_name_mapping file."
            required = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    if parsed_args["verbose"]
        println("Parsed args:")
        for (arg,val) in parsed_args
            println("  $arg  =>  $val")
        end
    end

    FindSL.run(parsed_args["pairwise-interaction-file"],
               parsed_args["sl-results-file"],
               parsed_args["lowerbound"],
               parsed_args["upperbound"],
               parsed_args["downregulated-cutoff"],
               parsed_args["upregulated-cutoff"],
               parsed_args["pairwise-interaction-file"],
               parsed_args["db-id-file"],
               parsed_args["tissue-type"],
               parsed_args["verbose"])
end

main()

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
function _second_order_sensitivities(loss,prob,alg,sensealg::ForwardDiffOverAdjoint,
                                     args...;kwargs...)
   ForwardDiff.jacobian(prob.p) do p
     x = Zygote.gradient(p) do _p
       loss(concrete_solve(prob,alg,prob.u0,_p,args...;sensealg=sensealg.adjalg,kwargs...))
     end
     first(x)
   end
end

function _second_order_sensitivity_product(loss,v,prob,alg,sensealg::ForwardDiffOverAdjoint,
                                           args...;kwargs...)

   θ = ForwardDiff.Dual.(prob.p,v)
   _loss = p -> loss(concrete_solve(prob,alg,prob.u0,p,args...;sensealg=sensealg.adjalg,kwargs...))
   getindex.(ForwardDiff.partials.(Zygote.gradient(_loss,θ)[1]),1)
end


text:
using MyPkgDemo
using Test

@testset "test for cat_str.jl" begin
 # Write your tests here.
   @test cat_str("a", "b") == "ab"
   @test cat_str("c", "s") == "cs"
   @test cat_str("w", "wsxdc") == "wwsxdc"
   @test cat_str("shfjg", "sdhdh", "ddd") == "shfjgsdhdhddd"
   @test cat_str("d", "w", "x") == "dwx"
end

max_stars_repo_name: junedwx/MyPkgDemo

max_issues_repo_name: junedwx/MyPkgDemo

max_forks_repo_name: junedwx/MyPkgDemo


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

using Test
using Mortar2D: calculate_normals

elements = Dict(
    1 => [1, 2],
    2 => [2, 3])
element_types = Dict(
    1 => :Seg2,
    2 => :Seg2)
X = Dict(
    1 => [7.0, 7.0],
    2 => [4.0, 3.0],
    3 => [0.0, 0.0])
n = calculate_normals(elements, element_types, X)
@test isapprox(n[1], 1/5*[4, -3])
@test isapprox(n[2], sqrt(2)/2*[1, -1])
@test isapprox(n[3], 1/5*[3, -4])


text:
# Configuration for running simple neutral copy error simulation
#=
Recommended command line to run:
> julia -p 4 -L NeutralEvolution.jl nrun.jl examples/example1
=#
export simtype
@everywhere simtype = 3
@everywhere const N = 1        # population size
const mutstddev = 0.04
const ngens = 1000
#const initial_value = exp(5.0)
const initial_value = 1.0
const num_trials=100000
const record_interval = 50
const use_population = true
const log_error = true
const neutral = true


text:
#production and acquisition

function PACost!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    key = (options.postprocess ? "postprocessing" : "evaluating")
    println(key*" windTurbine cost")
    windTurbine!(options,sizes,data,output)
    println(key*" supportStructures cost")
    supportStructures!(options,sizes,data,output)
    println(key*" powerTransmissionSystem cost")
    powerTransmissionSystem!(options,sizes,data,output)
    println(key*" monitoringSystem cost")
    monitoringSystem!(options,sizes,data,output)

    nothing
end

function windTurbine!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.projectManagement == 0
        if options.postprocess
            nothing
        else
            Cwtmat = 3e6*log.(data.turbinePR).-662400
            Cwttrans = data.vesselCost.*data.vesselTimeWT
            for i in 1:sizes.winds,j in 1:sizes.turbines
                output.CAPEX[i,j] += (Cwttrans[i]+Cwtmat[j])*data.turbineQuantity[j]
                output.windTurbine[i,j] += (Cwttrans[i]+Cwtmat[j])*data.turbineQuantity[j]
            end
        end
    else
        error("windTurbine mode not found.")
    end

    nothing

end

function supportStructures!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.supportStructures == 0
        if options.postprocess
            nothing
        else
            Cssmat = zeros(sizes.winds,sizes.turbines)
            Csstrans = data.vesselCost.*data.vesselTimeSS

            for i in 1:sizes.winds,j in 1:sizes.turbines
                Cssmat[i,j] = 339200*data.turbinePR[j]*(1+0.02*(data.WD[i]-8))*(1+8e-7*(data.height[j]*data.diameter[j]*data.diameter[j]/2-1e5))
                output.CAPEX[i,j] += (Csstrans[i]+Cssmat[i,j])*data.turbineQuantity[j]
                output.supportStructures[i,j] += (Csstrans[i]+Cssmat[i,j])*data.turbineQuantity[j]
            end
        end
    else
        error("supportStructures mode not found.")
    end

    nothing

end

function powerTransmissionSystem!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.powerTransmissionSystem == 0
        if options.postprocess
            nothing
        else
            Cofsubs = 583300 .+ 107900*data.IC
            Consubs = Cofsubs/2
            Ccable = data.Ntln.*data.Ctln.*data.onshoreDistance + data.Ntlf.*data.Ctlf.*data.offshoreDistance + data.Cprot

            output.CAPEX .+= transpose(Cofsubs + Consubs)
            output.CAPEX .+= Ccable
            output.powerTransmissionSystem .+= transpose(Cofsubs + Consubs)
            output.powerTransmissionSystem .+= Ccable
        end
    else
        error("powerTransmissionSystem mode not found.")
    end
    nothing

end

function monitoringSystem!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.monitoringSystem == 0
        if options.postprocess
            nothing
        else
            output.CAPEX .+= transpose((data.cmsCost+data.scadaCost).*data.turbineQuantity)
            output.monitoringSystem .+= transpose((data.cmsCost+data.scadaCost).*data.turbineQuantity)
        end
    else
        error("monitoringSystem mode not found.")
    end
    nothing

end

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
using Pkg, DrWatson
#@quickactivate "SR2StanPluto"

function tf1(f::Number; l=12, kwargs...)
    kwargs
 end

tf1(1) |> display
tf1(3) |> display
res = tf1(5; xlims=(6, 9), ylims=(0, 20))
res |> display


text:
# The increments to the `__counter` variable is for testing purposes

function _check_transa(t::Char)
    if !(t in ('C', 'N', 'T'))
        error("transa: is '$t', must be 'N', 'T', or 'C'")
    end
end

mkl_size(t::Char, M::AbstractVecOrMat) = t == 'N' ? size(M) : reverse(size(M))

# Checks sizes for the multiplication C <- A * B
function _check_mat_mult_matvec(C, A, B, tA)
    _size(v::AbstractMatrix) = size(v)
    _size(v::AbstractVector) = (size(v,1), 1)
    _str(v::AbstractMatrix) = string("[", size(v,1), ", ", size(v,2), "]")
    _str(v::AbstractVector) = string("[", size(v,1), "]")
    mA, nA = mkl_size(tA, A)
    mB, nB = _size(B)
    mC, nC = _size(C)
    if nA != mB || mC != mA || nC != nB
        t = ""
        if tA == 'T'; t = ".\'"; end
        if tA == 'C'; t = "\'"; end
        str = string("arrays had inconsistent dimensions for C <- A", t, " * B: ", _str(C), " <- ", _str(A), t, " * ", _str(B))
        throw(DimensionMismatch(str))
    end
end

for (mv, sv, mm, sm, T) in ((:mkl_scscmv, :mkl_scscsv, :mkl_scscmm, :mkl_scscsm, :Float32),
                            (:mkl_dcscmv, :mkl_dcscsv, :mkl_dcscmm, :mkl_dcscsm, :Float64),
                            (:mkl_ccscmv, :mkl_ccscsv, :mkl_ccscmm, :mkl_ccscsm, :ComplexF32),
                            (:mkl_zcscmv, :mkl_zcscsv, :mkl_zcscmm, :mkl_zcscsm, :ComplexF64))
@eval begin
function cscmv!(transa::Char, α::$T, matdescra::String,
                A::SparseMatrixCSC{$T, BlasInt}, x::StridedVector{$T},
                β::$T, y::StridedVector{$T})
    _check_transa(transa)
    _check_mat_mult_matvec(y, A, x, transa)
    __counter[] += 1
    ccall(($(string(mv)), :libmkl_rt), Cvoid,
        (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$T},
         Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
         Ptr{BlasInt}, Ptr{$T}, Ref{$T}, Ptr{$T}),
        transa, A.m, A.n, α,
        matdescra, A.nzval, A.rowval, A.colptr,
        pointer(A.colptr, 2), x, β, y)
    return y
end

function cscmm!(transa::Char, α::$T, matdescra::String,
                A::SparseMatrixCSC{$T, BlasInt}, B::StridedMatrix{$T},
                β::$T, C::StridedMatrix{$T})
    _check_transa(transa)
    _check_mat_mult_matvec(C, A, B, transa)
    mB, nB = size(B)
    mC, nC = size(C)
    __counter[] += 1
    ccall(($(string(mm)), :libmkl_rt), Cvoid,
        (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
         Ref{$T}, Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt},
         Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ref{BlasInt},
         Ref{$T}, Ptr{$T}, Ref{BlasInt}),
        transa, A.m, nC, A.n,
        α, matdescra, A.nzval, A.rowval,
        A.colptr, pointer(A.colptr, 2), B, mB,
        β, C, mC)
    return C
end

function cscsv!(transa::Char, α::$T, matdescra::String,
                A::SparseMatrixCSC{$T, BlasInt}, x::StridedVector{$T},
                y::StridedVector{$T})
    n = checksquare(A)
    _check_transa(transa)
    _check_mat_mult_matvec(y, A, x, transa)
    __counter[] += 1
    ccall(($(string(sv)), :libmkl_rt), Cvoid,
        (Ref{UInt8}, Ref{BlasInt}, Ref{$T}, Ptr{UInt8},
         Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
         Ptr{$T}, Ptr{$T}),
        transa, A.m, α, matdescra,
        A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2),
        x, y)
    return y
end

function cscsm!(transa::Char, α::$T, matdescra::String,
                A::SparseMatrixCSC{$T, BlasInt}, B::StridedMatrix{$T},
                C::StridedMatrix{$T})
    mB, nB = size(B)
    mC, nC = size(C)
    n = checksquare(A)
    _check_transa(transa)
    _check_mat_mult_matvec(C, A, B, transa)
    __counter[] += 1
    ccall(($(string(sm)), :libmkl_rt), Cvoid,
        (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$T},
         Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
         Ptr{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{$T},
         Ref{BlasInt}),
        transa, A.n, nC, α,
        matdescra, A.nzval, A.rowval, A.colptr,
        pointer(A.colptr, 2), B, mB, C,
        mC)
    return C
end

end # @eval
end


text:
using CompScienceMeshes
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)
pVI = point(10,11,12)
pVII = point(10,11,13)
pVIII = point(11,11,12)

Sourcechart = simplex(pI,pII,pIII)
Testchart = simplex(pVIII,pVII,pVI)

Accuracy = 12

function integrand(x,y)
			return(((x-pI)'*(y-pVII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

result = sauterschwabintegral(Sourcechart, Testchart, integrand, Accuracy, Accuracy)
println(result)

#=For those who want to test the sauterschwab_nonparameterized() function,
may uncomment the following three lines=#

#pd = PositiveDistance(Accuracy)
#result2 = sauterschwab_nonparameterized(Sourcechart, Testchart, integrand, pd)
#println(result2)

#=In this case the two charts from above can be used and the order of the points
within the simplex() fucntions can be changed arbitrarily. The user may also
compare the two results and see that both are equal.=#


text:
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

const fontface = cv_fontface("DejaVu Sans")

function setup_star_arc_painters(setup::CV_SceneSetupChain, cut_test=nothing)

    layout = setup.layout
    multiply_pos = CV_MultiplyByFactor(ComplexF64)

    trafo = cv_get_trafo(layout)
    trafo_domain = multiply_pos
    trafo_codomain = w -> trafo(multiply_pos(w))

    common_style = cv_op_over → cv_linewidth(3) → cv_antialias_best
    star_style = common_style → cv_color(1,1,1,0.8)
    arc_style = common_style → cv_color(0,0,0,0.8)
    star_lines, arc_lines = cv_star_arc_lines(
        tuple(LinRange(0.2, 1.5, 6)...), tuple(LinRange(-π/4, π/4, 6)...))

    CLP = CV_LinePainter
    star_painter_domain   = star_style ↦ CLP(trafo_domain,   star_lines)
    arc_painter_domain    = arc_style  ↦ CLP(trafo_domain,   arc_lines )
    star_painter_codomain = star_style ↦ CLP(trafo_codomain, star_lines, false,
        trafo_domain, cut_test)
    arc_painter_codomain  = arc_style  ↦ CLP(trafo_codomain, arc_lines, false,
        trafo_domain, cut_test)

    state_counter = cv_get_state_counter(layout)
    cc_can_domain = cv_get_cc_can_domain(layout)
    cc_can_codomain = cv_get_cc_can_codomain(layout)

    update_painter_func = z -> begin
        multiply_pos.factor = z
        state = state_counter.value

        if state==3 || state==5
            cv_paint(cc_can_domain, arc_painter_domain)
            cv_paint(cc_can_codomain, arc_painter_codomain)
        end

        if state==4 || state==5
            cv_paint(cc_can_domain, star_painter_domain)
            cv_paint(cc_can_codomain, star_painter_codomain)
        end
    end
    return cv_combine(setup; update_painter_func)
end

function setup_axis(setup::CV_SceneSetupChain)
    label_style = cv_black → fontface → cv_fontsize(20)
    ticks1 = "%.0f" ⇒  -2:2
    ticks1h =    "" ⇒  -2.5:2.5
    ticks2 = "%.0f" ⇒  -3:3
    ticks3 = ("-π" ⇒ -pi, "-π/2" ⇒ -pi/2, "0" ⇒ 0.0, "π/2" ⇒ pi/2, "π" ⇒ pi,)

    app_l = CV_TickLabelAppearance(; label_style, tick_length=10)
    app_s = CV_TickLabelAppearance(; label_style, tick_length=6)
    return cv_setup_lr_axis(setup,
        (app_l ↦ ticks1, app_s ↦ ticks1h,), (app_l ↦ ticks1, app_s ↦ ticks1h,),
        (app_l ↦ ticks2,), (app_l ↦ ticks3,))
end

function setup_axis_grid(setup::CV_SceneSetupChain)

    layout = setup.layout
    axis_grid_style = cv_op_over → cv_linewidth(1) → cv_color(0.7, 0.7, 0.7, 0.8)
    ag_domain_painter = CV_GridPainter(
        range(-2.0, stop=2.0, step=0.5), range(-2.0, stop=2.0, step=0.5))
    axis_grid_domain = axis_grid_style ↦ ag_domain_painter
    axis_grid_codomain = axis_grid_style ↦ CV_GridPainter(
        range(-3, stop=3, step=1), range(-π, stop=π, step=π/2))

    cc_can_domain = cv_get_cc_can_domain(layout)
    cc_can_codomain = cv_get_cc_can_codomain(layout)

    update_painter_func = z -> begin
        cv_paint(cc_can_domain, axis_grid_domain)
        cv_paint(cc_can_codomain, axis_grid_codomain)
        return nothing
    end
    return cv_combine(setup; update_painter_func)
end

function get_layout()
    trafo = log

    domain   = CV_Math2DCanvas(-2.5 + 2.5im, 2.5 - 2.5im, 100)
    codomain = CV_Math2DCanvas(-3.5 + (pi+0.1)im, 3.5 - (pi+0.1)im, 100)

    layout1 = CV_StateLayout(CV_2DLayout(), CV_CyclicValue(5))
    layout2 = cv_add(layout1, trafo, domain, codomain)
    return cv_do_lr_layout(layout2, 80)
end

function do_setup(layout, cut_test)
    setup = cv_setup_lr_border(cv_setup_cycle_state(CV_LRSetupChain(layout)))
    layout = setup.layout
    trafo = cv_get_trafo(layout)
    cc_can_domain = cv_get_cc_can_domain(layout)
    cc_can_codomain = cv_get_cc_can_codomain(layout)
    portrait_painter_domain = CV_PortraitPainter(trafo)
    portrait_painter_codomain = CV_PortraitPainter()
    update_painter_func = z -> begin
        cv_paint(cc_can_domain, portrait_painter_domain)
        cv_paint(cc_can_codomain, portrait_painter_codomain)
        return nothing
    end
    setup = cv_combine(setup; update_painter_func)
    setup = setup_axis_grid(setup)
    setup = cv_setup_lr_painters(setup; cut_test,
        portrait_painter_domain=nothing,
        portrait_painter_codomain=nothing)
    setup = setup_star_arc_painters(setup, cut_test)
    setup = setup_axis(setup)
    cv_add_padding!(setup.layout, 30)
    scene = cv_setup_domain_codomain_scene(setup)
    return scene
end

function main()
    cut_test = cv_create_angle_cross_test(pi, 0, Inf; δ=1e-2)

    setup = do_setup(get_layout(), cut_test)

    cv_scene_lr_start(setup; z_start=1.0+0.0im, state_start=3)

    scene = setup.layout

    cv_save_image(cv_get_can_layout(scene), "./Log02.png")

    if CVG !== nothing
        handler = CVG.cvg_visualize(scene)
        CVG.cvg_wait_for_destroy(handler.window)

        CVG.cvg_close(handler); CVG.cv_destroy(scene.layout)
    end
    return nothing
end

main()

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
"""
A single point or an array of points on the cells of a Triangulation
CellField objects can be evaluated efficiently at CellPoint instances.
"""
struct CellPoint{DS} <: CellDatum
  cell_ref_point::AbstractArray{<:Union{Point,AbstractArray{<:Point}}}
  cell_phys_point::AbstractArray{<:Union{Point,AbstractArray{<:Point}}}
  trian::Triangulation
  domain_style::DS
end

function CellPoint(
  cell_ref_point::AbstractArray{<:Union{Point,AbstractArray{<:Point}}},
  trian::Triangulation,
  domain_style::ReferenceDomain)

  cell_map = get_cell_map(trian)
  cell_phys_point = lazy_map(evaluate,cell_map,cell_ref_point)
  CellPoint(cell_ref_point,cell_phys_point,trian,domain_style)
end

function CellPoint(
  cell_phys_point::AbstractArray{<:Union{Point,AbstractArray{<:Point}}},
  trian::Triangulation,
  domain_style::PhysicalDomain)
  cell_map = get_cell_map(trian)
  cell_invmap = lazy_map(inverse_map,cell_map)
  cell_ref_point = lazy_map(evaluate,cell_invmap,cell_phys_point)
  CellPoint(cell_ref_point,cell_phys_point,trian,domain_style)
end

function get_data(f::CellPoint)
  if DomainStyle(f) == ReferenceDomain()
    f.cell_ref_point
  else
    f.cell_phys_point
  end
end

get_triangulation(f::CellPoint) = f.trian
DomainStyle(::Type{CellPoint{DS}}) where DS = DS()

function change_domain(a::CellPoint,::ReferenceDomain,::PhysicalDomain)
  CellPoint(a.cell_ref_point,a.cell_phys_point,a.trian,PhysicalDomain())
end

function change_domain(a::CellPoint,::PhysicalDomain,::ReferenceDomain)
  CellPoint(a.cell_ref_point,a.cell_phys_point,a.trian,ReferenceDomain())
end

# Possibly with a different name
"""
"""
function get_cell_points(trian::Triangulation)
  cell_ref_coords = get_cell_ref_coordinates(trian)
  cell_phys_coords = get_cell_coordinates(trian)
  CellPoint(cell_ref_coords,cell_phys_coords,trian,ReferenceDomain())
end

"""
"""
abstract type CellField <: CellDatum end

function Base.show(io::IO,::MIME"text/plain",f::CellField)
  show(io,f)
  print(io,":")
  print(io,"\n num_cells: $(num_cells(f))")
  print(io,"\n DomainStyle: $(DomainStyle(f))")
  print(io,"\n Triangulation: $(get_triangulation(f))")
  print(io,"\n Triangulation id: $(objectid(get_triangulation(f)))")
end

function CellField(f::Function,trian::Triangulation,domain_style::DomainStyle)
  s = size(get_cell_map(trian))
  cell_field = Fill(GenericField(f),s)
  GenericCellField(cell_field,trian,PhysicalDomain())
end

function CellField(f::Number,trian::Triangulation,domain_style::DomainStyle)
  s = size(get_cell_map(trian))
  cell_field = Fill(ConstantField(f),s)
  GenericCellField(cell_field,trian,domain_style)
end

function CellField(f::AbstractArray{<:Number},trian::Triangulation,domain_style::DomainStyle)
  @check length(f)==num_cells(trian)  """\n
  You are trying to build a CellField from an array of length $(length(f))
  on a Triangulation with $(num_cells(trian)) cells. The length of the given array
  and the number of cells should match.
  """
  cell_field = lazy_map(ConstantField,f)
  GenericCellField(cell_field,trian,domain_style)
end

function CellField(f::CellField,trian::Triangulation,domain_style::DomainStyle)
  change_domain(f,trian,domain_style)
end

function CellField(f,trian::Triangulation)
  CellField(f,trian,ReferenceDomain())
end

function get_normal_vector(trian::Triangulation)
  cell_normal = get_facet_normal(trian)
  @assert ! isa(cell_normal,SkeletonPair)
  GenericCellField(cell_normal,trian,ReferenceDomain())
end

evaluate!(cache,f::Function,x::CellPoint) = CellField(f,get_triangulation(x))(x)

function change_domain(a::CellField,::ReferenceDomain,::PhysicalDomain)
  trian = get_triangulation(a)
  cell_map = get_cell_map(trian)
  cell_invmap = lazy_map(inverse_map,cell_map)
  cell_field_ref = get_data(cell_field)
  cell_field_phys = lazy_map(Broadcasting(∘),cell_field_ref,cell_invmap)
  GenericCellField(cell_field_phys,trian,PhysicalDomain())
end

function change_domain(a::CellField,::PhysicalDomain,::ReferenceDomain)
  trian = get_triangulation(a)
  cell_map = get_cell_map(trian)
  cell_field_phys = get_data(a)
  cell_field_ref = lazy_map(Broadcasting(∘),cell_field_phys,cell_map)
  GenericCellField(cell_field_ref,trian,ReferenceDomain())
end

"""
"""
function change_domain(a::CellField,target_trian::Triangulation,target_domain::DomainStyle)
  change_domain(a,DomainStyle(a),target_trian,target_domain)
end

function change_domain(a::CellField,::ReferenceDomain,trian::Triangulation,::ReferenceDomain)
  trian_a = get_triangulation(a)
  if have_compatible_domains(trian_a,trian)
    return a
  elseif have_compatible_domains(trian_a,get_background_triangulation(trian))
    cell_id = get_cell_to_bgcell(trian)
    @assert ! isa(cell_id,SkeletonPair)
    cell_a_q = lazy_map(Reindex(get_data(a)),cell_id)
    cell_s2q = get_cell_ref_map(trian)
    cell_field = lazy_map(Broadcasting(∘),cell_a_q,cell_s2q)
    GenericCellField(cell_field,trian,ReferenceDomain())
  else
    @unreachable """\n
    We cannot move the given CellField to the reference domain of the requested triangulation.
    Make sure that the given triangulation is either the same as the triangulation on which the
    CellField is defined, or that the latter triangulation is the background of the former.
    """
  end
end

function change_domain(a::CellField,::PhysicalDomain,trian::Triangulation,::PhysicalDomain)
  trian_a = get_triangulation(a)
  if have_compatible_domains(trian_a,trian)
    return a
  elseif have_compatible_domains(trian_a,get_background_triangulation(trian))
    cell_id = get_cell_to_bgcell(trian)
    @assert ! isa(cell_id,SkeletonPair)
    cell_field = lazy_map(Reindex(get_data(a)),cell_id)
    GenericCellField(cell_field,trian,PhysicalDomain())
  else
    @unreachable """\n
    We cannot move the given CellField to the physical domain of the requested triangulation.
    Make sure that the given triangulation is either the same as the triangulation on which the
    CellField is defined, or that the latter triangulation is the background of the former.
    """
  end
end

function change_domain(a::CellField,::PhysicalDomain,trian::Triangulation,::ReferenceDomain)
  a_trian = change_domain(a,trian,PhysicalDomain())
  change_domain(a_trian,ReferenceDomain())
end

function change_domain(a::CellField,::ReferenceDomain,trian::Triangulation,::PhysicalDomain)
  a_phys = change_domain(a,PhysicalDomain())
  change_domain(a_phys,trian,PhysicalDomain())
end

"""
"""
struct GenericCellField{DS} <: CellField
  cell_field::AbstractArray{<:Union{Field,AbstractArray{<:Field}}}
  trian::Triangulation
  domain_style::DS
  function GenericCellField(
    cell_field::AbstractArray{<:Union{Field,AbstractArray{<:Field}}},
    trian::Triangulation,
    domain_style::DomainStyle)

    DS = typeof(domain_style)
    new{DS}(Fields.MemoArray(cell_field),trian,domain_style)
  end
end

get_data(f::GenericCellField) = f.cell_field
get_triangulation(f::GenericCellField) = f.trian
DomainStyle(::Type{GenericCellField{DS}}) where DS = DS()

# Evaluation of CellFields

(a::CellField)(x) = evaluate(a,x)

function evaluate!(cache,f::CellField,x::Point)
  @notimplemented """\n
  Evaluation of a CellField at a given Point is not implemented yet.

  This is a feature that we want to have at some point in Gridap.
  If you are ready to help with this implementation, please contact the
  Gridap administrators.
  """
end

function evaluate!(cache,f::CellField,x::CellPoint)
  _f, _x = _to_common_domain(f,x)
  cell_field = get_data(_f)
  cell_point = get_data(_x)
  lazy_map(evaluate,cell_field,cell_point)
end

function _to_common_domain(f::CellField,x::CellPoint)

  trian_f = get_triangulation(f)
  trian_x = get_triangulation(x)

  if have_compatible_domains(trian_f,trian_x)
    nothing
  elseif have_compatible_domains(trian_f,get_background_triangulation(trian_x))
    nothing
  elseif have_compatible_domains(trian_x,get_background_triangulation(trian_f))
    @unreachable """\n
    CellField objects defined on a sub-triangulation cannot be evaluated
    on the underlying background mesh.

    This happens e.g. when trying to evaluate a CellField defined on a Neumann boundary
    at a CellPoint defined on the underlying background mesh.
    """
  else
    @unreachable """\n
    Your are trying to evaluate a CellField on a CellPoint object defined on incompatible
    triangulations. Verify that either the two objects are defined in the same triangulation
    or that the triangulaiton of the CellField is the background triangulation of the CellPoint.
    """
  end

  f_on_trian_x = change_domain(f,trian_x,DomainStyle(x))
  f_on_trian_x, x
end

# Gradient

function gradient(a::CellField)
  cell_∇a = lazy_map(Broadcasting(∇),get_data(a))
  if DomainStyle(a) == PhysicalDomain()
    g = cell_∇a
  else
    cell_map = get_cell_map(get_triangulation(a))
    g = lazy_map(Broadcasting(push_∇),cell_∇a,cell_map)
  end
  GenericCellField(g,get_triangulation(a),DomainStyle(a))
end

function ∇∇(a::CellField)
  cell_∇∇a = lazy_map(Broadcasting(∇∇),get_data(a))
  if DomainStyle(a) == PhysicalDomain()
    h = cell_∇∇a
  else
    cell_map = get_cell_map(get_triangulation(a))
    h = lazy_map(Broadcasting(push_∇∇),cell_∇∇a,cell_map)
  end
  GenericCellField(h,get_triangulation(a),DomainStyle(a))
end

# This function has to be removed when ∇⋅∇(a) is implemented
laplacian(a::CellField) = tr(∇∇(a))

# Operations between CellField

function evaluate!(cache,k::Operation,a::CellField...)
  _operate_cellfields(k,a...)
end

function evaluate!(cache,k::Operation,a::Union{Function,CellField}...)
  _operate_cellfields(k,_convert_to_cellfields(a...)...)
end

function evaluate!(cache,k::Operation,a::Union{Number,CellField}...)
  _operate_cellfields(k,_convert_to_cellfields(a...)...)
end

function evaluate!(cache,k::Operation,a::Union{AbstractArray{<:Number},CellField}...)
  _operate_cellfields(k,_convert_to_cellfields(a...)...)
end

# Why julia hangs with this method????
#
#function evaluate!(cache,k::Operation,a::Union{Function,Number,AbstractArray,CellField}...)
#  b = _convert_to_cellfields(a...)
#  _operate_cellfields(k,b...)
#end

struct OperationCellField{DS} <: CellField
  op::Operation
  args::Tuple
  trian::Triangulation
  domain_style::DS
  memo::Dict{Any,Any}
  function OperationCellField(op::Operation,args::CellField...)

    @assert length(args) > 0
    trian = get_triangulation(first(args))
    domain_style = DomainStyle(first(args))
    @check all( map(i->DomainStyle(i)==domain_style,args) )
    @check all( map(i->have_compatible_domains(get_triangulation(i),trian),args) )

    if num_cells(trian)>0
      x = _get_cell_points(args...)
      try
         ax = map(i->i(x),args)
         axi = map(first,ax)
         r = Fields.BroadcastingFieldOpMap(op.op)(axi...)
      catch
        @unreachable """\n
        It is not possible to perform operation $(op.op) on the given cell fields.

        See the caught error for more information. (If you are using Visual
          Studio Code REPL you might not see the caught error, please use the
          command-line REPL).
        """
      end
    end

    new{typeof(domain_style)}(op,args,trian,domain_style,Dict())
  end
end

function _get_cell_points(args::CellField...)
  k = findfirst(i->isa(i,CellState),args)
  if k === nothing
    j = findall(i->isa(i,OperationCellField),args)
    if length(j) == 0
      _get_cell_points(first(args))
    else
      _get_cell_points(args[j]...)
    end
  else
    args[k].points
  end
end

function _get_cell_points(a::CellField)
  trian = get_triangulation(a)
  get_cell_points(trian)
end

function _get_cell_points(a::OperationCellField...)
  b = []
  for ai in a
    for i in ai.args
      push!(b,i)
    end
  end
  _get_cell_points(b...)
end

function _get_cell_points(a::OperationCellField)
  _get_cell_points(a.args...)
end

function get_data(f::OperationCellField)
  a = map(get_data,f.args)
  lazy_map(Broadcasting(f.op),a...)
end
get_triangulation(f::OperationCellField) = f.trian
DomainStyle(::Type{OperationCellField{DS}}) where DS = DS()

function evaluate!(cache,f::OperationCellField,x::CellPoint)
  #key = (:evaluate,objectid(x))
  #if ! haskey(f.memo,key)
  #  ax = map(i->i(x),f.args)
  #  f.memo[key] = lazy_map(Fields.BroadcastingFieldOpMap(f.op.op),ax...)
  #end
  #f.memo[key]
  ax = map(i->i(x),f.args)
  lazy_map(Fields.BroadcastingFieldOpMap(f.op.op),ax...)
end

function change_domain(f::OperationCellField,target_trian::Triangulation,target_domain::DomainStyle)
  args = map(i->change_domain(i,target_trian,target_domain),f.args)
  OperationCellField(f.op,args...)
end

function _operate_cellfields(k::Operation,a...)
  b = _to_common_domain(a...)
  OperationCellField(k,b...)
end

function _convert_to_cellfields(a...)
  a1 = filter(i->isa(i,CellField),a)
  a2 = _to_common_domain(a1...)
  target_domain = DomainStyle(first(a2))
  target_trian = get_triangulation(first(a2))
  map(i->CellField(i,target_trian,target_domain),a)
end

function _to_common_domain(a::CellField...)

  # Find a suitable domain style
  if any( map(i->DomainStyle(i)==ReferenceDomain(),a) )
    target_domain = ReferenceDomain()
  else
    target_domain = PhysicalDomain()
  end

  # Find a suitable triangulation
  msg = """\n
  You are trying to operate CellField objects defined on incompatible triangulations.

  Make sure that all CellField objects are defined on the background triangulation
  or that the number of different sub-triangulations is equal to one.

  For instace:

  - 3 cell fields 2, two them on the same Neumann boundary and the other on the background mesh is OK.

  - 2 cell fields defined on 2 different Neumann boundaries is NOT OK.
  """
  trian_candidates = unique(objectid,map(get_triangulation,a))
  if length(trian_candidates) == 1
    target_trian = first(trian_candidates)
  elseif length(trian_candidates) == 2
    trian_a, trian_b = trian_candidates
    if have_compatible_domains(trian_a,trian_b)
      target_trian = trian_a
    elseif have_compatible_domains(trian_a,get_background_triangulation(trian_b))
      target_trian = trian_b
    elseif have_compatible_domains(trian_b,get_background_triangulation(trian_a))
      target_trian = trian_a
    else
      @unreachable msg
    end
  else
    @unreachable msg
  end
  map(i->change_domain(i,target_trian,target_domain),a)
end

# Composition (this replaces the @law macro)

Base.:(∘)(f::Function,g::CellField) = Operation(f)(g)
Base.:(∘)(f::Function,g::Tuple{Vararg{CellField}}) = Operation(f)(g...)
Base.:(∘)(f::Function,g::Tuple{Vararg{Union{AbstractArray{<:Number},CellField}}}) = Operation(f)(g...)
Base.:(∘)(f::Function,g::Tuple{Vararg{Union{Function,CellField}}}) = Operation(f)(g...)

# Define some of the well known arithmetic ops

# Unary ops

for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj)
  @eval begin
    ($op)(a::CellField) = Operation($op)(a)
  end
end

# Binary ops

for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/)
  @eval begin
    ($op)(a::CellField,b::CellField) = Operation($op)(a,b)
    ($op)(a::CellField,b::Number) = Operation($op)(a,b)
    ($op)(a::Number,b::CellField) = Operation($op)(a,b)
    ($op)(a::CellField,b::Function) = Operation($op)(a,b)
    ($op)(a::Function,b::CellField) = Operation($op)(a,b)
    ($op)(a::CellField,b::AbstractArray{<:Number}) = Operation($op)(a,b)
    ($op)(a::AbstractArray{<:Number},b::CellField) = Operation($op)(a,b)
  end
end

dot(::typeof(∇),f::CellField) = divergence(f)
function (*)(::typeof(∇),f::CellField)
  msg = "Syntax ∇*f has been removed, use ∇⋅f (\\nabla \\cdot f) instead"
  error(msg)
end
outer(::typeof(∇),f::CellField) = gradient(f)
outer(f::CellField,::typeof(∇)) = transpose(gradient(f))
cross(::typeof(∇),f::CellField) = curl(f)

"""
    get_physical_coordinate(trian::Triangulation)

In contrast to get_cell_map, the returned object:
- is a [`CellField`](@ref)
- its gradient is the identity tensor
"""
function get_physical_coordinate(trian::Triangulation)
  CellField(_phys_coord,trian)
end

_phys_coord(x) = x

_phys_coord_grad(x) = one(typeof(outer(x,x)))

gradient(::typeof(_phys_coord)) = _phys_coord_grad

# Skeleton related Operations

function Base.getproperty(x::CellField, sym::Symbol)
  if sym in (:⁺,:plus)
    CellFieldAt{:plus}(x)
  elseif sym in (:⁻, :minus)
    CellFieldAt{:minus}(x)
  else
    getfield(x, sym)
  end
end

function Base.propertynames(x::CellField, private=false)
  (fieldnames(typeof(x))...,:⁺,:plus,:⁻,:minus)
end

struct CellFieldAt{T,F} <: CellField
  parent::F
  CellFieldAt{T}(parent::CellField) where T = new{T,typeof(parent)}(parent)
end

get_data(f::CellFieldAt) = get_data(f.parent)
get_triangulation(f::CellFieldAt) = get_triangulation(f.parent)
DomainStyle(::Type{CellFieldAt{T,F}}) where {T,F} = DomainStyle(F)
gradient(a::CellFieldAt{P}) where P = CellFieldAt{P}(gradient(a.parent))
∇∇(a::CellFieldAt{P}) where P = CellFieldAt{P}(∇∇(a.parent))

function CellFieldAt{T}(parent::OperationCellField) where T
  args = map(i->CellFieldAt{T}(i),parent.args)
  OperationCellField(parent.op,args...)
end

function get_normal_vector(trian::SkeletonTriangulation)
  cell_normal_plus = get_facet_normal(trian.plus)
  cell_normal_minus = get_facet_normal(trian.minus)
  plus = GenericCellField(cell_normal_plus,trian,ReferenceDomain())
  minus = GenericCellField(cell_normal_minus,trian,ReferenceDomain())
  SkeletonPair(plus,minus)
end

for op in (:outer,:*,:dot)
  @eval begin
    ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
    ($op)(a::SkeletonPair{<:CellField},b::CellField) = Operation($op)(a,b)
  end
end

function evaluate!(cache,k::Operation,a::CellField,b::SkeletonPair{<:CellField})
  plus = k(a.plus,b.plus)
  minus = k(a.minus,b.minus)
  SkeletonPair(plus,minus)
end

function evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField},b::CellField)
  plus = k(a.plus,b.plus)
  minus = k(a.minus,b.minus)
  SkeletonPair(plus,minus)
end

jump(a::CellField) = a.⁺ - a.⁻
jump(a::SkeletonPair{<:CellField}) = a.⁺ + a.⁻ # a.⁻ results from multiplying by n.⁻. Thus we need to sum.

mean(a::CellField) = Operation(_mean)(a.⁺,a.⁻)
_mean(x,y) = 0.5*x + 0.5*y

# This is the fundamental part to make operations on the skeleton work.

function change_domain(a::CellField,target_trian::SkeletonTriangulation,target_domain::DomainStyle)
  trian_a = get_triangulation(a)
  if have_compatible_domains(trian_a,target_trian)
    return change_domain(a,target_domain)
  elseif have_compatible_domains(trian_a,get_background_triangulation(target_trian))
    # In this case, we can safely take either plus or minus arbitrarily.
    if isa(a,GenericCellField) && isa(get_array(a.cell_field),Fill{<:ConstantField})
      a_on_target_trian = change_domain(a,target_trian.plus,target_domain)
      return GenericCellField(get_data(a_on_target_trian),target_trian,target_domain)
    elseif isa(a,GenericCellField) && isa(get_array(a.cell_field),Fill{<:GenericField{<:Function}})
      a_on_target_trian = change_domain(a,target_trian.plus,target_domain)
      return GenericCellField(get_data(a_on_target_trian),target_trian,target_domain)
    else
      @unreachable """\n
      It is not possible to use the given CellField on a SkeletonTriangulation.
      Make sure that you are specifying which of the two possible traces,
      either plus (aka ⁺) or minus (aka ⁻) you want to use.
      """
    end
  else
    @unreachable """\n
    We cannot move the given CellField to the requested triangulation.
    Make sure that the given CellField is defined on the triangulation you want to work with.
    """
  end
end

function change_domain(a::CellFieldAt,trian::SkeletonTriangulation,target_domain::DomainStyle)
  trian_a = get_triangulation(a)
  if have_compatible_domains(trian_a,get_background_triangulation(trian))
    plus, minus = change_domain_skeleton(a.parent,trian,target_domain)
    if isa(a,CellFieldAt{:plus})
      return plus
    elseif isa(a,CellFieldAt{:minus})
      return minus
    else
      @unreachable
    end
  else
    @unreachable """\n
    It is not allowd to writte `u.⁺` of `u.⁻` for the given CellField.
    Make sure that the CellField `u` is either defined on the background mesh
    or it is a normal vector extracted from a SkeletonTriangulation.
    """
  end
end

function change_domain_skeleton(a::CellField,trian::SkeletonTriangulation,target_domain::DomainStyle)
  a_on_plus_trian = change_domain(a,trian.plus,target_domain)
  a_on_minus_trian = change_domain(a,trian.minus,target_domain)
  plus = GenericCellField(get_data(a_on_plus_trian),trian,target_domain)
  minus = GenericCellField(get_data(a_on_minus_trian),trian,target_domain)
  plus, minus
end

function change_domain(f::OperationCellField,target_trian::SkeletonTriangulation,target_domain::DomainStyle)
  args = map(i->change_domain(i,target_trian,target_domain),f.args)
  OperationCellField(f.op,args...)
end

# Just to provide more meaningful error messages
function (a::SkeletonPair{<:CellField})(x)
  @unreachable """\n
  You are trying to evaluate a CellField on a mesh skeleton but you have not specified which of the
  two sides i.e. plus (aka ⁺) or minus (aka ⁻) you want to select.

  For instance, if you have extracted the normal vector and the cell points from a SkeletonTriangulation

      x = get_cell_points(strian)
      n = get_normal_vector(strian)

  Evaluating `n(x)` is not allowed. You need to call either `n.⁺(x)` or `n.⁻(x)`.
  """
end


text:
module StopWatchesBase

export @stopwatch

using UUIDs: UUID

include("recording.jl")
include("tables.jl")
include("noop.jl")

function __init__()
    init_records()
end

end


text:
# This is a small scale experiment as an alternative to Sage's (and Nemo's)
# parent model to allow for easier low-level operations.
#
# For example, a monomial is represented as a vector (or rather a static vector)
# but to manipulate it, you probably want to know more: what is the term order?
# what variables are module position? This is store in a "parent" variable.
#
# So typically, a monomial will be represented by a structure with two fields:
# the actual exponents, and the parent. But say that you have an array with many
# monomials, all sharing the same parent. Julia will probably box them (because
# of the pointer to the parent) and the information of the parent is very
# repetitive.
#
# So instead, we define an abstract type Context{T}, which is meant to
# "interpret" the elements of type T. And when, we want to store a element
# together with its parent, we enclose is in a Γpair. For low-level code, we
# just carry around the context and never use Γpairs.

abstract type Context{T} end
Base.eltype(::Context{T}) where T = T

struct Γpair0{T, Γ<:Context{T}}
    ctx::Γ
    dat::T
end
const Γpair{Γ, T} = Γpair0{T, Γ}

gpair(ctx::Γ, dat) where {T, Γ<:Context{T}} = Γpair{Γ, T}(ctx, dat)

# Default conversion, but some conversion may require more information about the context.
(ctx::Context{T})(x...) where {T} = eltype(ctx)(x...)

#.. Sorting interface

# To make it possible to sort arrays of type T, define lt(::Context{T}, ::T, ::T)
# Of course, you can also define directly lt(::T, ::T) but it is not always possible.

struct Ordering{Γ<:Context}<:Base.Order.Ordering
    ctx::Γ
end
order(ctx::Γ) where Γ = Ordering{Γ}(ctx)

function Base.Order.lt(order::Ordering{Γ}, a, b) where {Γ<:Context}
    lt(order.ctx, a, b)
end


text:
module Util

"""
create annotations for any struct
"""
macro namedargs(StructDec)
  sstruct = string(StructDec)
  typename = sstruct[findfirst(r"(?<=struct).*(?=\n)", sstruct)]
  typename = strip(typename)
  parametric = findfirst(r"{.*}", typename)
  if parametric == nothing
    parametric = ""
  else
    parametric = "where $(typename[parametric])"
  end
  typename = typename[findfirst(r"\w+", typename)]

  function fieldNameAndType(expr)
    return split(string(expr), "::")
  end

  args = StructDec.args[3].args[2:2:end]
  args = filter(a -> !occursin(r"new\(+", string(a)), args)
  args = map(a -> fieldNameAndType(a), args)
  fnames = [a[1] for a in args]
  ftypes = [a[2] for a in args]

  callArgs = join(["$(a[1])" for a in args], ", ")
  oldConstructorCall = "$typename($callArgs)"

  fn_args = join(["$(a[1])::$(a[2])" for a in args], ", ")
  namedArgsConstructor = "$typename(; $fn_args) $parametric"

  return quote
    $(esc(StructDec))
    $(esc(Meta.parse("$namedArgsConstructor = $oldConstructorCall")))
  end
end

end # module Util

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
# Julia Interface to DSP
# Original interface by: Kibaek Kim; Argonne National Laboratory, 2016
module DspCInterface

#import ..Dsp
import Compat: String, unsafe_wrap
using Pkg
#import JuMP
using MPI
#Pkg.installed("MPI") == nothing || using MPI
#"MPI" in keys(Pkg.installed())  || using MPI
#export DspModel, @dsp_ccall,freeDSP,

###############################################################################
# Help functions
###############################################################################
macro dsp_ccall(func, args...)
    @static if Sys.isunix()
        return esc(quote
            ccall(($func, "libDsp"), $(args...))
        end)
    end
    @static if Sys.iswindows()
        return esc(quote
            ccall(($func, "libDsp"), stdcall, $(args...))
        end)
    end
end

mutable struct DspModel
    p::Ptr{Nothing}

    # Number of blocks
    nblocks::Int

    # solve_type should be one of these:
    # :Dual
    # :Benders
    # :Extensive
    solve_type

    numRows::Int
    numCols::Int
    primVal
    dualVal
    colVal::Vector{Float64}
    rowVal::Vector{Float64}

    # MPI settings
    comm
    comm_size::Int
    comm_rank::Int

    # Array of block ids:
    # The size of array is not necessarily same as nblocks,
    # as block ids may be distributed to multiple processors.
    block_ids::Vector{Integer}

    function DspModel()
        # assign Dsp pointer
        p = @dsp_ccall("createEnv", Ptr{Nothing}, ())
        # initialize variables
        nblocks = 0
        solve_type = :Dual
        numRows = 0
        numCols = 0
        primVal = NaN
        dualVal = NaN
        colVal = Vector{Float64}()
        rowVal = Vector{Float64}()
        comm = nothing
        comm_size = 1
        comm_rank = 0
        block_ids = Vector{Integer}()
        # create DspModel
        dsp = new(p, nblocks, solve_type, numRows, numCols, primVal, dualVal, colVal, rowVal, comm, comm_size, comm_rank, block_ids)
        # with finalizer
        finalizer(freeDSP,dsp)
        # return DspModel
        return dsp
    end
end

function freeDSP(dsp::DspModel)
    if dsp.p == C_NULL
        return
    else
        @dsp_ccall("freeEnv", Nothing, (Ptr{Nothing},), dsp.p)
        dsp.p = C_NULL
    end
    dsp.nblocks = 0
    dsp.solve_type = nothing
    dsp.numRows = 0
    dsp.numCols = 0
    dsp.primVal = NaN
    dsp.dualVal = NaN
    dsp.colVal = Vector{Float64}()
    dsp.rowVal = Vector{Float64}()
    return
end

function freeModel(dsp::DspModel)
    check_problem(dsp)
    @dsp_ccall("freeModel", Nothing, (Ptr{Nothing},), dsp.p)
    dsp.nblocks = 0
    dsp.numRows = 0
    dsp.numCols = 0
    dsp.primVal = NaN
    dsp.dualVal = NaN
    dsp.colVal = Vector{Float64}()
    dsp.rowVal = Vector{Float64}()
end

function check_problem(dsp::DspModel)
    if dsp.p == C_NULL
        error("Invalid DspModel")
    end
end

function readParamFile(dsp::DspModel, param_file::AbstractString)
    check_problem(dsp)
    @dsp_ccall("readParamFile", Nothing, (Ptr{Nothing}, Ptr{UInt8}), dsp.p, param_file);
end

###############################################################################
# Block IDs
###############################################################################
function setBlockIds(dsp::DspModel, nblocks::Integer)
    check_problem(dsp)
    # set number of blocks
    dsp.nblocks = nblocks
    # set MPI settings
    if isdefined(Main,:MPI) && MPI.Initialized()
        dsp.comm = MPI.COMM_WORLD
        dsp.comm_size = MPI.Comm_size(dsp.comm)
        dsp.comm_rank = MPI.Comm_rank(dsp.comm)
    end
    #@show dsp.nblocks
    #@show dsp.comm
    #@show dsp.comm_size
    #@show dsp.comm_rank
    # get block ids with MPI settings
    dsp.block_ids = getBlockIds(dsp)
    #@show dsp.block_ids
    #@show dsp.block_ids
    # send the block ids to Dsp
    @dsp_ccall("setIntPtrParam", Nothing, (Ptr{Nothing}, Ptr{UInt8}, Cint, Ptr{Cint}),
        dsp.p, "ARR_PROC_IDX", convert(Cint, length(dsp.block_ids)), convert(Vector{Cint}, dsp.block_ids .- 1))
end

function getBlockIds(dsp::DspModel)
    check_problem(dsp)
    # processor info
    mysize = dsp.comm_size
    myrank = dsp.comm_rank
    # empty block ids
    proc_idx_set = Int[]
    # DSP is further parallelized with mysize > dsp.nblocks.
    modrank = myrank % dsp.nblocks
    # If we have more than one processor,
    # do not assign a sub-block to the master.
    if mysize > 1
        if myrank == 0
            return proc_idx_set
        end
        # exclude master
        mysize -= 1;
        modrank = (myrank-1) % dsp.nblocks
    end
    # assign sub-blocks in round-robin fashion
    for s = modrank:mysize:(dsp.nblocks-1)
        push!(proc_idx_set, s+1)
    end
    # return assigned block ids
    return proc_idx_set
end

###############################################################################
# Load problems
###############################################################################

function readSmps(dsp::DspModel, filename::AbstractString)
    # Check pointer to TssModel
    check_problem(dsp)
    # read smps files
    @dsp_ccall("readSmps", Nothing, (Ptr{Nothing}, Ptr{UInt8}), dsp.p, convert(Vector{UInt8}, filename))
    # set block Ids
    setBlockIds(dsp, getNumScenarios(dsp))
end

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

for func in [:freeSolver,
             :solveDe,
             :solveBd,
             :solveDd]
    strfunc = string(func)
    @eval begin
        function $func(dsp::DspModel)
            return @dsp_ccall($strfunc, Nothing, (Ptr{Nothing},), dsp.p)
        end
    end
end

for func in [:solveBdMpi, :solveDdMpi]
    strfunc = string(func)
    @eval begin
        function $func(dsp::DspModel, comm)
            return @dsp_ccall($strfunc, Nothing, (Ptr{Nothing}, MPI.CComm), dsp.p, convert(MPI.CComm, comm))
        end
    end
end

#Solve a DSP model
function solve(dsp::DspModel)
    check_problem(dsp)
    if dsp.comm_size == 1
        if dsp.solve_type == :Dual
            solveDd(dsp);
        elseif dsp.solve_type == :Benders
            solveBd(dsp);
        elseif dsp.solve_type == :Extensive
            solveDe(dsp);
        end
    elseif dsp.comm_size > 1
        if dsp.solve_type == :Dual
            solveDdMpi(dsp, dsp.comm);
        elseif dsp.solve_type == :Benders
            solveBdMpi(dsp, dsp.comm);
        elseif dsp.solve_type == :Extensive
            solveDe(dsp);
        end
    end
end

###############################################################################
# Get functions
###############################################################################
for (func,rtn) in [(:getNumScenarios, Cint),
                   (:getTotalNumCols, Cint),
                   (:getTotalNumRows, Cint),
                   (:getStatus, Cint),
                   (:getNumIterations, Cint),
                   (:getNumNodes, Cint),
                   (:getWallTime, Cdouble),
                   (:getPrimalBound, Cdouble),
                   (:getDualBound, Cdouble),
                   (:getNumCouplingRows, Cint)]
    strfunc = string(func)
    @eval begin
        function $func(dsp::DspModel)
            check_problem(dsp)
            return @dsp_ccall($strfunc, $rtn, (Ptr{Nothing},), dsp.p)
        end
    end
end
getSolutionStatus(dsp::DspModel) = getStatus(dsp)
function getNumRows(dsp::DspModel, num::Integer)
    @dsp_ccall("getNumRows", Cint, (Ptr{Nothing}, Cint), dsp.p, num)
end
function getNumCols(dsp::DspModel, num::Integer)
    @dsp_ccall("getNumCols", Cint, (Ptr{Nothing}, Cint), dsp.p, num)
end

function getObjCoef(dsp::DspModel)
    check_problem(dsp)
    num = getTotalNumCols()
    obj = Array{Cdouble}(undef,num)
    @dsp_ccall("getObjCoef", Nothing, (Ptr{Nothing}, Ptr{Cdouble}), dsp.p, obj)
    return obj
end

function getSolution(dsp::DspModel, num::Integer)
    sol = Array{Cdouble}(undef,num)
    @dsp_ccall("getPrimalSolution", Nothing, (Ptr{Nothing}, Cint, Ptr{Cdouble}), dsp.p, num, sol)
    return sol
end
getSolution(dsp::DspModel) = getSolution(dsp, getTotalNumCols(dsp))

function getDualSolution(dsp::DspModel, num::Integer)
    sol = Array{Cdouble}(undef,num)
    @dsp_ccall("getDualSolution", Nothing, (Ptr{Nothing}, Cint, Ptr{Cdouble}), dsp.p, num, sol)
    return sol
end
getDualSolution(dsp::DspModel) = getDualSolution(dsp, getNumCouplingRows(dsp))

###############################################################################
# Set functions
###############################################################################
function setSolverType(dsp::DspModel, solver)
    check_problem(dsp)
    solver_types = [:DualDecomp, :Benders, :ExtensiveForm]
    if solver in solver_types
        dsp.solver = solver
    else
        @warn("Solver type $solver is invalid.")
    end
end

end # end of module