text:
module CHOLMOD

export                                  # types
 CholmodDense,
 CholmodFactor,
 CholmodSparse,
 CholmodTriplet,

 CholmodSparse!,                        # destructive constructors
 CholmodDense!,

 etree

using Base.LinAlg.UMFPACK               # for decrement, increment, etc.

import Base: (*), convert, copy, ctranspose, eltype, findnz, getindex, hcat,
             isvalid, nfilled, show, size, sort!, transpose, vcat

import ..LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
                 cholfact, cholfact!, copy, det, diag,
                 full, logdet, norm, scale, scale!, solve, sparse

include("cholmod_h.jl")

macro isok(A)
    :($A==CHOLMOD_TRUE || throw(CholmodException()))
end

const chm_ver    = Array(Cint, 3)
if dlsym(dlopen("libcholmod"), :cholmod_version) != C_NULL
    ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), chm_ver)
else
    ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), chm_ver)
end
const chm_com_sz = ccall((:jl_cholmod_common_size,:libsuitesparse_wrapper),Int,())
const chm_com    = fill(0xff, chm_com_sz)
const chm_l_com  = fill(0xff, chm_com_sz)
## chm_com and chm_l_com must be initialized at runtime because they contain pointers
## to functions in libc.so, whose addresses can change
function cmn(::Type{Int32})
    if isnan(reinterpret(Float64,chm_com[1:8])[1])
        @isok ccall((:cholmod_start, :libcholmod), Cint, (Ptr{Uint8},), chm_com)
    end
    chm_com
end
function cmn(::Type{Int64})
    if isnan(reinterpret(Float64,chm_l_com[1:8])[1])
        @isok ccall((:cholmod_l_start, :libcholmod), Cint, (Ptr{Uint8},), chm_l_com)
    end
    chm_l_com
end

typealias CHMITypes Union(Int32,Int64)
typealias CHMVTypes Union(Complex64, Complex128, Float32, Float64)

type CholmodException <: Exception end

## macro to generate the name of the C function according to the integer type
macro chm_nm(nm,typ) string("cholmod_", eval(typ) == Int64 ? "l_" : "", nm) end

### A way of examining some of the fields in chm_com
### Probably better to make this a Dict{ASCIIString,Tuple} and
### save the offsets and the lengths and the types.  Then the names can be checked.
type ChmCommon
    dbound::Float64
    maxrank::Int
    supernodal_switch::Float64
    supernodal::Int32
    final_asis::Int32
    final_super::Int32
    final_ll::Int32
    final_pack::Int32
    final_monotonic::Int32
    final_resymbol::Int32
    prefer_zomplex::Int32               # should always be false
    prefer_upper::Int32
    print::Int32                        # print level. Default: 3
    precise::Int32                      # print 16 digits, otherwise 5
    nmethods::Int32                     # number of ordering methods
    selected::Int32
    postorder::Int32
    itype::Int32
    dtype::Int32
end

### These offsets should be reconfigured to be less error-prone in matches
const chm_com_offsets = Array(Int, length(ChmCommon.types))
ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper),
      Void, (Ptr{Int},), chm_com_offsets)
const chm_final_ll_inds = (1:4) + chm_com_offsets[7]
const chm_prt_inds = (1:4) + chm_com_offsets[13]
const chm_ityp_inds = (1:4) + chm_com_offsets[18]

### there must be an easier way but at least this works.
function ChmCommon(aa::Array{Uint8,1})
    typs = ChmCommon.types
    sz = map(sizeof, typs)
    args = map(i->reinterpret(typs[i], aa[chm_com_offsets[i] + (1:sz[i])])[1], 1:length(sz))
    eval(Expr(:call, unshift!(args, :ChmCommon), Any))
end

function set_chm_prt_lev(cm::Array{Uint8}, lev::Integer) # can probably be removed
    cm[(1:4) + chm_com_offsets[13]] = reinterpret(Uint8, [int32(lev)])
end

## cholmod_dense pointers passed to or returned from C functions are of Julia type
## Ptr{c_CholmodDense}.  The CholmodDense type contains a c_CholmodDense object and other
## fields then ensure the memory pointed to is freed when it should be and not before.
type c_CholmodDense{T<:CHMVTypes}
    m::Int
    n::Int
    nzmax::Int
    lda::Int
    xpt::Ptr{T}
    zpt::Ptr{Void}
    xtype::Cint
    dtype::Cint
end

type CholmodDense{T<:CHMVTypes}
    c::c_CholmodDense
    mat::Matrix{T}
end

if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later
    type c_CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}
        n::Int
        minor::Int
        Perm::Ptr{Ti}
        ColCount::Ptr{Ti}
        IPerm::Ptr{Ti}        # this pointer was added in verison 2.1.0
        nzmax::Int
        p::Ptr{Ti}
        i::Ptr{Ti}
        x::Ptr{Tv}
        z::Ptr{Void}
        nz::Ptr{Ti}
        next::Ptr{Ti}
        prev::Ptr{Ti}
        nsuper::Int
        ssize::Int
        xsize::Int
        maxcsize::Int
        maxesize::Int
        super::Ptr{Ti}
        pi::Ptr{Ti}
        px::Ptr{Ti}
        s::Ptr{Ti}
        ordering::Cint
        is_ll::Cint
        is_super::Cint
        is_monotonic::Cint
        itype::Cint
        xtype::Cint
        dtype::Cint
    end

    type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}
        c::c_CholmodFactor{Tv,Ti}
        Perm::Vector{Ti}
        ColCount::Vector{Ti}
        IPerm::Vector{Ti}
        p::Vector{Ti}
        i::Vector{Ti}
        x::Vector{Tv}
        nz::Vector{Ti}
        next::Vector{Ti}
        prev::Vector{Ti}
        super::Vector{Ti}
        pi::Vector{Ti}
        px::Vector{Ti}
        s::Vector{Ti}
    end

    function CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodFactor{Tv,Ti}})
        cfp = unsafe_load(cp)
        Perm = pointer_to_array(cfp.Perm, (cfp.n,), true)
        ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true)
        IPerm = pointer_to_array(cfp.IPerm, (cfp.IPerm == C_NULL ? 0 : cfp.n + 1,), true)
        p = pointer_to_array(cfp.p, (cfp.p == C_NULL ? 0 : cfp.n + 1,), true)
        i = pointer_to_array(cfp.i, (cfp.i == C_NULL ? 0 : cfp.nzmax,), true)
        x = pointer_to_array(cfp.x, (cfp.x == C_NULL ? 0 : max(cfp.nzmax,cfp.xsize),), true)
        nz = pointer_to_array(cfp.nz, (cfp.nz == C_NULL ? 0 : cfp.n,), true)
        next = pointer_to_array(cfp.next, (cfp.next == C_NULL ? 0 : cfp.n + 2,), true)
        prev = pointer_to_array(cfp.prev, (cfp.prev == C_NULL ? 0 : cfp.n + 2,), true)
        super = pointer_to_array(cfp.super, (cfp.super == C_NULL ? 0 : cfp.nsuper + 1,), true)
        pi = pointer_to_array(cfp.pi, (cfp.pi == C_NULL ? 0 : cfp.nsuper + 1,), true)
        px = pointer_to_array(cfp.px, (cfp.px == C_NULL ? 0 : cfp.nsuper + 1,), true)
        s = pointer_to_array(cfp.s, (cfp.s == C_NULL ? 0 : cfp.ssize + 1,), true)
        cf = CholmodFactor{Tv,Ti}(cfp, Perm, ColCount, IPerm, p, i, x, nz, next, prev,
                                  super, pi, px, s)
        c_free(cp)
        cf
    end
else
    type c_CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}
        n::Int
        minor::Int
        Perm::Ptr{Ti}
        ColCount::Ptr{Ti}
        nzmax::Int
        p::Ptr{Ti}
        i::Ptr{Ti}
        x::Ptr{Tv}
        z::Ptr{Void}
        nz::Ptr{Ti}
        next::Ptr{Ti}
        prev::Ptr{Ti}
        nsuper::Int
        ssize::Int
        xsize::Int
        maxcsize::Int
        maxesize::Int
        super::Ptr{Ti}
        pi::Ptr{Ti}
        px::Ptr{Ti}
        s::Ptr{Ti}
        ordering::Cint
        is_ll::Cint
        is_super::Cint
        is_monotonic::Cint
        itype::Cint
        xtype::Cint
        dtype::Cint
    end

    type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}
        c::c_CholmodFactor{Tv,Ti}
        Perm::Vector{Ti}
        ColCount::Vector{Ti}
        p::Vector{Ti}
        i::Vector{Ti}
        x::Vector{Tv}
        nz::Vector{Ti}
        next::Vector{Ti}
        prev::Vector{Ti}
        super::Vector{Ti}
        pi::Vector{Ti}
        px::Vector{Ti}
        s::Vector{Ti}
    end

    function CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodFactor{Tv,Ti}})
        cfp = unsafe_load(cp)
        Perm = pointer_to_array(cfp.Perm, (cfp.n,), true)
        ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true)
        p = pointer_to_array(cfp.p, (cfp.p == C_NULL ? 0 : cfp.n + 1,), true)
        i = pointer_to_array(cfp.i, (cfp.i == C_NULL ? 0 : cfp.nzmax,), true)
        x = pointer_to_array(cfp.x, (cfp.x == C_NULL ? 0 : max(cfp.nzmax,cfp.xsize),), true)
        nz = pointer_to_array(cfp.nz, (cfp.nz == C_NULL ? 0 : cfp.n,), true)
        next = pointer_to_array(cfp.next, (cfp.next == C_NULL ? 0 : cfp.n + 2,), true)
        prev = pointer_to_array(cfp.prev, (cfp.prev == C_NULL ? 0 : cfp.n + 2,), true)
        super = pointer_to_array(cfp.super, (cfp.super == C_NULL ? 0 : cfp.nsuper + 1,), true)
        pi = pointer_to_array(cfp.pi, (cfp.pi == C_NULL ? 0 : cfp.nsuper + 1,), true)
        px = pointer_to_array(cfp.px, (cfp.px == C_NULL ? 0 : cfp.nsuper + 1,), true)
        s = pointer_to_array(cfp.s, (cfp.s == C_NULL ? 0 : cfp.ssize + 1,), true)
        cf = CholmodFactor{Tv,Ti}(cfp, Perm, ColCount, p, i, x, nz, next, prev,
                                  super, pi, px, s)
        c_free(cp)
        cf
    end
end

type c_CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}
    m::Int
    n::Int
    nzmax::Int
    ppt::Ptr{Ti}
    ipt::Ptr{Ti}
    nzpt::Ptr{Void}
    xpt::Ptr{Tv}
    zpt::Ptr{Void}
    stype::Cint
    itype::Cint
    xtype::Cint
    dtype::Cint
    sorted::Cint
    packed::Cint
end

type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}
    c::c_CholmodSparse{Tv,Ti}
    colptr0::Vector{Ti}
    rowval0::Vector{Ti}
    nzval::Vector{Tv}
end

type c_CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes}
    m::Int
    n::Int
    nzmax::Int
    nnz::Int
    i::Ptr{Ti}
    j::Ptr{Ti}
    x::Ptr{Tv}
    z::Ptr{Void}
    stype:Cint
    itype::Cint
    xtype::Cint
    dtype::Cint
end

type CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes}
    c::c_CholmodTriplet{Tv,Ti}
    i::Vector{Ti}
    j::Vector{Ti}
    x::Vector{Tv}
end

eltype{T<:CHMVTypes}(A::CholmodDense{T}) = T
eltype{T<:CHMVTypes}(A::CholmodFactor{T}) = T
eltype{T<:CHMVTypes}(A::CholmodSparse{T}) = T
eltype{T<:CHMVTypes}(A::CholmodTriplet{T}) = T

## The CholmodDense! constructor does not copy the contents, which is generally what you
## want as most uses of CholmodDense objects are read-only.
function CholmodDense!{T<:CHMVTypes}(aa::VecOrMat{T}) # uses the memory from Julia
    m = size(aa,1); n = size(aa,2)
    CholmodDense(c_CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, aa),
                                   C_NULL, xtyp(T), dtyp(T)),
                 length(size(aa)) == 2 ? aa : reshape(aa, (m,n)))
end

## The CholmodDense constructor copies the contents
function CholmodDense{T<:CHMVTypes}(aa::VecOrMat{T})
    m = size(aa,1); n = size(aa,2)
    acp = length(size(aa)) == 2 ? copy(aa) : reshape(copy(aa), (m,n))
    CholmodDense(c_CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, acp),
                                   C_NULL, xtyp(T), dtyp(T)), acp)
end

function CholmodDense{T<:CHMVTypes}(c::Ptr{c_CholmodDense{T}})
    cp = unsafe_load(c)
    if cp.lda != cp.m || cp.nzmax != cp.m * cp.n
        throw(DimensionMismatch("overallocated cholmod_dense returned object of size $(cp.m) by $(cp.n) with leading dim $(cp.lda) and nzmax $(cp.nzmax)"))
    end
    ## the true in the call to pointer_to_array means Julia will free the memory
    val = CholmodDense(cp, pointer_to_array(cp.xpt, (cp.m,cp.n), true))
    c_free(c)
    val
end
CholmodDense!{T<:CHMVTypes}(c::Ptr{c_CholmodDense{T}}) = CholmodDense(c) # no distinction

function isvalid{T<:CHMVTypes}(cd::CholmodDense{T})
    bool(ccall((:cholmod_check_dense, :libcholmod), Cint,
               (Ptr{c_CholmodDense{T}}, Ptr{Uint8}), &cd.c, cmn(Int32)))
end

function chm_eye{T<:Union(Float64,Complex128)}(m::Integer, n::Integer, t::T)
    CholmodDense(ccall((:cholmod_eye, :libcholmod), Ptr{c_CholmodDense{T}},
                       (Int, Int, Cint, Ptr{Uint8}),
                       m, n,xtyp(T),cmn(Int32)))
end
chm_eye(m::Integer, n::Integer) = chm_eye(m, n, 1.)
chm_eye(n::Integer) = chm_eye(n, n, 1.)

function chm_ones{T<:Union(Float64,Complex128)}(m::Integer, n::Integer, t::T)
    CholmodDense(ccall((:cholmod_ones, :libcholmod), Ptr{c_CholmodDense{T}},
                       (Int, Int, Cint, Ptr{Uint8}),
                       m, n, xtyp(T), cmn(Int32)))
end
chm_ones(m::Integer, n::Integer) = chm_ones(m, n, 1.)

function chm_zeros{T<:Union(Float64,Complex128)}(m::Integer, n::Integer, t::T)
    CholmodDense(ccall((:cholmod_zeros, :libcholmod), Ptr{c_CholmodDense{T}},
                       (Int, Int, Cint, Ptr{Uint8}),
                       m, n, xtyp(T), cmn(Int32)))
end
chm_zeros(m::Integer, n::Integer) = chm_zeros(m, n, 1.)

function chm_print{T<:CHMVTypes}(cd::CholmodDense{T}, lev::Integer, nm::ASCIIString)
    cm = cmn(Int32)
    orig = cm[chm_prt_inds]
    cm[chm_prt_inds] = reinterpret(Uint8, [int32(lev)])
    @isok ccall((:cholmod_print_dense, :libcholmod), Cint,
                   (Ptr{c_CholmodDense{T}}, Ptr{Uint8}, Ptr{Uint8}),
                   &cd.c, nm, cm)
    cm[chm_prt_inds] = orig
end
chm_print(cd::CholmodDense, lev::Integer) = chm_print(cd, lev, "")
chm_print(cd::CholmodDense) = chm_print(cd, int32(4), "")
show(io::IO,cd::CholmodDense) = chm_print(cd, int32(4), "")

function copy{Tv<:CHMVTypes}(B::CholmodDense{Tv})
    CholmodDense(ccall((:cholmod_copy_dense,:libcholmod), Ptr{c_CholmodDense{Tv}},
                       (Ptr{c_CholmodDense{Tv}},Ptr{Uint8}), &B.c, cmn(Int32)))
end

function norm{Tv<:CHMVTypes}(D::CholmodDense{Tv},p::Real=1)
    ccall((:cholmod_norm_dense, :libcholmod), Float64,
          (Ptr{c_CholmodDense{Tv}}, Cint, Ptr{Uint8}),
          &D.c, p == 1 ? 1 :(p == Inf ? 1 : throw(ArgumentError("p must be 1 or Inf"))),cmn(Int32))
end

function CholmodSparse!{Tv<:CHMVTypes,Ti<:CHMITypes}(colpt::Vector{Ti},
                                                     rowval::Vector{Ti},
                                                     nzval::Vector{Tv},
                                                     m::Integer,
                                                     n::Integer,
                                                     stype::Signed)
    bb = colpt[1]
    if bb != 0 && bb != 1 throw(DimensionMismatch("colpt[1] is $bb, must be 0 or 1")) end
    if any(diff(colpt) .< 0) throw(ArgumentError("elements of colpt must be non-decreasing")) end
    if length(colpt) != n + 1 throw(DimensionMismatch("length(colptr) = $(length(colpt)), should be $(n+1)")) end
    if bool(bb)                         # one-based
        decrement!(colpt)
        decrement!(rowval)
    end
    nz = colpt[end]
    if length(rowval) != nz || length(nzval) != nz
        throw(DimensionMismatch("length(rowval) = $(length(rowval)) and length(nzval) = $(length(nzval)) should be $nz"))
    end
    if any(rowval .< 0) || any(rowval .>= m)
        throw(ArgumentError("all elements of rowval must be in the range [0,$(m-1)]"))
    end
    it = ityp(Ti)
    sort!(CholmodSparse(c_CholmodSparse{Tv,Ti}(m,n,int(nz),convert(Ptr{Ti},colpt),
                                               convert(Ptr{Ti},rowval), C_NULL,
                                               convert(Ptr{Tv},nzval), C_NULL,
                                               int32(stype), ityp(Ti),
                                               xtyp(Tv),dtyp(Tv),
                                               CHOLMOD_FALSE,CHOLMOD_TRUE),
                        colpt,rowval,nzval))
end
function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(colpt::Vector{Ti},
                                                    rowval::Vector{Ti},
                                                    nzval::Vector{Tv},
                                                    m::Integer,
                                                    n::Integer,
                                                    stype::Signed)
    CholmodSparse!(copy(colpt),copy(rowval),copy(nzval),m,n,stype)
end
function CholmodSparse!{Tv<:CHMVTypes,Ti<:CHMITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed)
    CholmodSparse!(A.colptr,A.rowval,A.nzval,size(A,1),size(A,2),stype)
end
function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed)
    CholmodSparse!(copy(A.colptr),copy(A.rowval),copy(A.nzval),size(A,1),size(A,2),stype)
end
function CholmodSparse(A::SparseMatrixCSC)
    stype = ishermitian(A) ? 1 : 0
    CholmodSparse(stype > 0 ? triu(A) : A, stype)
end
function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodSparse{Tv,Ti}})
    csp = unsafe_load(cp)
    colptr0 = pointer_to_array(csp.ppt, (csp.n + 1,), true)
    nnz = int(colptr0[end])
    cms = CholmodSparse{Tv,Ti}(csp, colptr0,
                               pointer_to_array(csp.ipt, (nnz,), true),
                               pointer_to_array(csp.xpt, (nnz,), true))
    c_free(cp)
    cms
end
CholmodSparse!{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodSparse{Tv,Ti}}) = CholmodSparse(cp)
CholmodSparse{Tv<:CHMVTypes}(D::CholmodDense{Tv}) = CholmodSparse(D,1) # default Ti is Int

function CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes}(tp::Ptr{c_CholmodTriplet{Tv,Ti}})
    ctp = unsafe_load(tp)
    i = pointer_to_array(ctp.i, (ctp.nnz,), true)
    j = pointer_to_array(ctp.j, (ctp.nnz,), true)
    x = pointer_to_array(ctp.x, (ctp.x == C_NULL ? 0 : ctp.nnz), true)
    ct = CholmodTriplet{Tv,Ti}(ctp, i, j, x)
    c_free(tp)
    ct
end

function chm_rdsp(fnm::String)
    fd = ccall(:fopen, Ptr{Void}, (Ptr{Uint8},Ptr{Uint8}), fnm, "r")
    res = ccall((:cholmod_read_sparse,:libcholmod), Ptr{c_CholmodSparse{Float64,Cint}},
                (Ptr{Void},Ptr{Uint8}),fd,cmn(Cint))
    ccall(:fclose, Cint, (Ptr{Void},), fd) # should do this in try/finally/end
    CholmodSparse(res)
end

for Ti in (:Int32,:Int64)
    @eval begin
        function (*){Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},
                                    B::CholmodSparse{Tv,$Ti})
            CholmodSparse(ccall((@chm_nm "ssmult" $Ti
                                 , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{c_CholmodSparse{Tv,$Ti}},
                                 Cint,Cint,Cint,Ptr{Uint8}), &A.c,&B.c,0,true,true,cmn($Ti)))
        end
        function A_mul_Bc{Tv<:Union(Float32,Float64)}(A::CholmodSparse{Tv,$Ti},
                                                      B::CholmodSparse{Tv,$Ti})
            cm = cmn($Ti)
            aa = Array(Ptr{c_CholmodSparse{Tv,$Ti}}, 2)
            if !is(A,B)
                aa[1] = ccall((@chm_nm "transpose" $Ti
                               ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                              (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                              &B.c,cm)
                aa[2] = ccall((@chm_nm "ssmult" $Ti
                               ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                              (Ptr{c_CholmodSparse{Tv,$Ti}},
                               Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                              &A.c, aa[1], cmn($Ti))
                @isok ccall((@chm_nm "free_sparse" $Ti
                                ,:libcholmod), Cint,
                               (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{Uint8}), aa, cm)
                CholmodSparse(aa[2])
            end
            ## The A*A' case is handled by cholmod_aat. Strangely the matrix returned by
            ## cholmod_aat is not marked as symmetric.  The code following the call to
            ## cholmod_aat is to create the symmetric-storage version of the result then
            ## transpose it to provide sorted columns.  The result is stored in the upper
            ## triangle
            aa[1] = ccall((@chm_nm "aat" $Ti
                           , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                          (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Void}, Int, Cint, Ptr{Uint8}),
                          &A.c, C_NULL, 0, 1, cm)
            ## Create the lower triangle unsorted
            aa[2] = ccall((@chm_nm "copy" $Ti
                           , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                          (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Cint, Ptr{Uint8}),
                          aa[1], -1, 1, cm)
            @isok ccall((@chm_nm "free_sparse" $Ti
                            , :libcholmod), Cint,
                           (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{Uint8}), aa, cm)
            aa[1] = aa[2]
            r = unsafe_load(aa[1])
            ## Now transpose the lower triangle to the upper triangle to do the sorting
            rpt = ccall((@chm_nm "allocate_sparse" $Ti
                         ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}},
                        (Csize_t,Csize_t,Csize_t,Cint,Cint,Cint,Cint,Ptr{Cuchar}),
                        r.m,r.n,r.nzmax,r.sorted,r.packed,-r.stype,r.xtype,cm)
            @isok ccall((@chm_nm "transpose_sym" $Ti
                            ,:libcholmod),Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{$Ti},
                            Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                           aa[1],1,C_NULL,rpt,cm)
            @isok ccall((@chm_nm "free_sparse" $Ti
                            , :libcholmod), Cint,
                           (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{Uint8}), aa, cm)
            CholmodSparse(rpt)
        end
        function Ac_mul_B{Tv<:Union(Float32,Float64)}(A::CholmodSparse{Tv,$Ti},
                                                      B::CholmodSparse{Tv,$Ti})
            cm = cmn($Ti)
            aa = Array(Ptr{c_CholmodSparse{Tv,$Ti}}, 2)
            aa[1] = ccall((@chm_nm "transpose" $Ti
                           ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                          (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                          &B.c,cm)
            if is(A,B)
                Ac = CholmodSparse(aa[1])
                return A_mul_Bc(Ac,Ac)
            end
            aa[2] = ccall((@chm_nm "ssmult" $Ti
                           ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                          (Ptr{c_CholmodSparse{Tv,$Ti}},
                           Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                          aa[1],&B.c,cm)
            @isok ccall((@chm_nm "free_sparse" $Ti
                            , :libcholmod), Cint,
                           (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{Uint8}), aa, cm)
            CholmodSparse(aa[2])
        end
        function CholmodDense{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            CholmodDense(ccall((@chm_nm "sparse_to_dense" $Ti
                                ,:libcholmod), Ptr{c_CholmodDense{Tv,$Ti}},
                               (Ptr{c_CholmodSparse{Tv,Ti}},Ptr{Uint8}),
                               &A.c,cmn($Ti)))
        end
        function CholmodSparse{Tv<:CHMVTypes}(D::CholmodDense{Tv},i::$Ti)
            CholmodSparse(ccall((@chm_nm "dense_to_sparse" $Ti
                                 ,:libcholmod), Ptr{c_CholmodSparse{Tv,Ti}},
                                (Ptr{c_CholmodDense{Tv,Ti}},Ptr{Uint8}),
                                &D.c,cmn($Ti)))
        end
        function CholmodSparse{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti})
            if bool(L.c.is_ll)
                return CholmodSparse(ccall((@chm_nm "factor_to_sparse" $Ti
                                            ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                           (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{Uint8}),
                                           &L.c,cmn($Ti)))
            end
            cm = cmn($Ti)
            Lcll = ccall((@chm_nm "copy_factor" $Ti
                          ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}},
                         (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{Uint8}),
                         &L.c,cm)
            @isok ccall((@chm_nm "change_factor" $Ti
                            ,:libcholmod), Cint,
                           (Cint,Cint,Cint,Cint,Cint,Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{Uint8}),
                           L.c.xtype,true,L.c.is_super,true,true,Lcll,cm)
            val = CholmodSparse(ccall((@chm_nm "factor_to_sparse" $Ti
                                       ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                      (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{Uint8}),
                                      Lcll,cmn($Ti)))
            @isok ccall((@chm_nm "free_factor" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{Ptr{c_CholmodFactor{Tv,$Ti}}},Ptr{Uint8}),
                           [Lcll],cm)
            val
        end
        function CholmodSparse{Tv<:CHMVTypes,Ti<:$Ti}(T::CholmodTriplet{Tv,Ti})
            CholmodSparse(ccall((@chm_nm "triplet_to_sparse" $Ti
                                 ,:libcholmod), Ptr{c_CholmodSparse{Tv,Ti}},
                               (Ptr{c_CholmodTriplet{Tv,Ti}},Ptr{Uint8}),
                               &T.c,cmn($Ti)))
        end
        function CholmodTriplet{Tv<:CHMVTypes,Ti<:$Ti}(A::CholmodSparse{Tv,Ti})
            CholmodTriplet(ccall((@chm_nm "sparse_to_triplet" $Ti
                                 ,:libcholmod), Ptr{c_CholmodTriplet{Tv,Ti}},
                               (Ptr{c_CholmodSparse{Tv,Ti}},Ptr{Uint8}),
                               &A.c,cmn($Ti)))
        end
        function isvalid{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti})
            bool(ccall((@chm_nm "check_factor" $Ti
                        ,:libcholmod), Cint,
                       (Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                       &L.c, cmn($Ti)))
        end
        function isvalid{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            bool(ccall((@chm_nm "check_sparse" $Ti
                        ,:libcholmod), Cint,
                       (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                       &A.c, cmn($Ti)))
        end
        function isvalid{Tv<:CHMVTypes}(T::CholmodTriplet{Tv,$Ti})
            bool(ccall((@chm_nm "check_triplet" $Ti
                        ,:libcholmod), Cint,
                       (Ptr{c_CholmodTriplet{Tv,$Ti}}, Ptr{Uint8}),
                       &T.c, cmn($Ti)))
        end
        function cholfact{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}, ll::Bool)
            cm = cmn($Ti)
            ## may need to change final_asis as well as final_ll
            if ll cm[chm_final_ll_inds] = reinterpret(Uint8, [one(Cint)]) end
            Lpt = ccall((@chm_nm "analyze" $Ti
                         ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}},
                        (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}), &A.c, cm)
            @isok ccall((@chm_nm "factorize" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}},
                            Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                           &A.c, Lpt, cm)
            CholmodFactor(Lpt)
        end
        function cholfact{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},beta::Tv,ll::Bool)
            cm = cmn($Ti)
            ## may need to change final_asis as well as final_ll
            if ll cm[chm_final_ll_inds] = reinterpret(Uint8, [one(Cint)]) end
            Lpt = ccall((@chm_nm "analyze" $Ti
                         ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}},
                        (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}), &A.c, cm)
            @isok ccall((@chm_nm "factorize_p" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Tv}, Ptr{Cint}, Csize_t,
                            Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                           &A.c, &beta, C_NULL, zero(Csize_t), Lpt, cmn($Ti))
            CholmodFactor(Lpt)
        end
        function chm_analyze{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            CholmodFactor(ccall((@chm_nm "analyze" $Ti
                                 ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}},
                                (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}), &A.c, cmn($Ti)))
        end
        function cholfact!{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti},A::CholmodSparse{Tv,$Ti})
            @isok ccall((@chm_nm "factorize" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}},
                            Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                           &A.c, &L.c, cmn($Ti))
            L
        end
        function cholfact!{Tv<:Float64}(L::CholmodFactor{Tv,$Ti},A::CholmodSparse{Tv,$Ti},
                                        beta::Tv)
            @isok ccall((@chm_nm "factorize_p" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Tv}, Ptr{Cint}, Csize_t,
                            Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                           &A.c, &beta, C_NULL, zero(Csize_t), &L.c, cmn($Ti))
            L
        end
        function chm_pack!{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti})
            @isok ccall((@chm_nm "pack_factor" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}),
                           &L.c,cmn($Ti))
            L
        end
        function chm_print{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti},lev,nm)
            cm = cmn($Ti)
            orig = cm[chm_prt_inds]
            cm[chm_prt_inds] = reinterpret(Uint8, [int32(lev)])
            @isok ccall((@chm_nm "print_factor" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{Uint8}, Ptr{Uint8}),
                           &L.c, nm, cm)
            cm[chm_prt_inds] = orig
        end
        function chm_print{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},lev,nm)
            cm = cmn($Ti)
            orig = cm[chm_prt_inds]
            cm[chm_prt_inds] = reinterpret(Uint8, [int32(lev)])
            @isok ccall((@chm_nm "print_sparse" $Ti
                           ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}, Ptr{Uint8}),
                           &A.c, nm, cm)
            cm[chm_prt_inds] = orig
        end
        function chm_scale!{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},
                                           S::CholmodDense{Tv},
                                           typ::Integer)
            @isok ccall((@chm_nm "scale" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodDense{Tv}},Cint,Ptr{c_CholmodSparse{Tv,$Ti}},
                            Ptr{Uint8}), &S.c, typ, &A.c, cmn($Ti))
            A
        end
        function chm_sdmult{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},
                                           trans::Bool,
                                           alpha::Real,
                                           beta::Real,
                                           X::CholmodDense{Tv})
            m,n = size(A)
            nc = trans ? m : n
            nr = trans ? n : m
            if nc != size(X,1)
                throw(DimensionMismatch("Incompatible dimensions, $nc and $(size(X,1)), in chm_sdmult"))
            end
            aa = float64([alpha, 0.])
            bb = float64([beta, 0.])
            Y = CholmodDense(zeros(Tv,nr,size(X,2)))
            @isok ccall((@chm_nm "sdmult" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Ptr{Cdouble},Ptr{Cdouble},
                            Ptr{c_CholmodDense{Tv}}, Ptr{c_CholmodDense{Tv}}, Ptr{Uint8}),
                           &A.c,trans,aa,bb,&X.c,&Y.c,cmn($Ti))
            Y
        end
        function chm_speye{Tv<:CHMVTypes,Ti<:$Ti}(m::Ti, n::Ti, x::Tv)
            CholmodSparse(ccall((@chm_nm "speye" $Ti
                                 , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Int, Int, Cint, Ptr{Uint8}),
                                m, n, xtyp{Tv}, cmn($Ti)))
        end
        function chm_spzeros{Tv<:Union(Float64,Complex128)}(m::$Ti, n::$Ti, nzmax::$Ti, x::Tv)
            CholmodSparse(ccall((@chm_nm "spzeros" $Ti
                         , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Int, Int, Int, Ptr{Uint8}),
                                m, n, nzmax, xtyp{Tv}, cmn($Ti)))
        end
## add chm_xtype and chm_pack
        function copy{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti})
            CholmodFactor(ccall((@chm_nm "copy_factor" $Ti
                                 ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}},
                                (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{Uint8}), &L.c, cmn($Ti)))
        end
        function copy{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            CholmodSparse(ccall((@chm_nm "copy_sparse" $Ti
                                 ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{Uint8}), &A.c, cmn($Ti)))
        end
        function copy{Tv<:CHMVTypes}(T::CholmodTriplet{Tv,$Ti})
            CholmodTriplet(ccall((@chm_nm "copy_triplet" $Ti
                                 ,:libcholmod), Ptr{c_CholmodTriplet{Tv,$Ti}},
                                (Ptr{c_CholmodTriplet{Tv,$Ti}},Ptr{Uint8}), &T.c, cmn($Ti)))
        end
        function ctranspose{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            CholmodSparse(ccall((@chm_nm "transpose" $Ti
                                 ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{Uint8}),
                                &A.c, 2, cmn($Ti)))
        end
        function etree{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            tr = Array($Ti,size(A,2))
            @isok ccall((@chm_nm "etree" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{$Ti},Ptr{Uint8}),
                           &A.c,tr,cmn($Ti))
            tr
        end
        function hcat{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},B::CholmodSparse{Tv,$Ti})
            ccall((@chm_nm "horzcat" $Ti
                   , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                  (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Ptr{Uint8}),
                  &A.c,&B.c,true,cmn($Ti))
        end
        function nfilled{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            ccall((@chm_nm "nnz" $Ti
                   ,:libcholmod), Int, (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{Uint8}),&A.c,cmn($Ti))
        end
        function norm{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},p::Real)
            ccall((@chm_nm "norm_sparse" $Ti
                   , :libcholmod), Float64,
                  (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{Uint8}),
                  &A.c,p == 1 ? 1 :(p == Inf ? 1 : throw(ArgumentError("p must be 1 or Inf"))),cmn($Ti))
        end
        function solve{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti},
                                      B::CholmodDense{Tv}, typ::Integer)
            CholmodDense(ccall((@chm_nm "solve" $Ti
                                ,:libcholmod), Ptr{c_CholmodDense{Tv}},
                               (Cint, Ptr{c_CholmodFactor{Tv,$Ti}},
                                Ptr{c_CholmodDense{Tv}}, Ptr{Uint8}),
                               typ, &L.c, &B.c, cmn($Ti)))
        end
        function solve{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti},
                                      B::CholmodSparse{Tv,$Ti},
                                      typ::Integer)
            CholmodSparse(ccall((@chm_nm "spsolve" $Ti
                                 ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Cint, Ptr{c_CholmodFactor{Tv,$Ti}},
                                 Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                                typ, &L.c, &B.c, cmn($Ti)))
        end
        function sort!{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            @isok ccall((@chm_nm "sort" $Ti
                            ,:libcholmod), Cint,
                           (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Uint8}),
                           &A.c, cmn($Ti))
            A
        end
        function copysym{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            CholmodSparse(ccall((@chm_nm "copy" $Ti
                                 ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Cint,Ptr{Uint8}),
                                &A.c,0,1,cmn($Ti)))
        end
        function transpose{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti})
            CholmodSparse(ccall((@chm_nm "transpose" $Ti
                                 ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}},
                                (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{Uint8}),
                                &A.c, 1, cmn($Ti)))
        end
        function vcat{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},B::CholmodSparse{Tv,$Ti})
            ccall((@chm_nm "vertcat" $Ti
                   , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}},
                  (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Ptr{Uint8}),
                  &A.c,&B.c,true,cmn($Ti))
        end
    end
end
(*){Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::CholmodDense{Tv}) = chm_sdmult(A,false,1.,0.,B)
(*){Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::VecOrMat{Tv}) = chm_sdmult(A,false,1.,0.,CholmodDense(B))

(\){T<:CHMVTypes}(L::CholmodFactor{T},B::CholmodDense{T}) = solve(L,B,CHOLMOD_A)
(\){T<:CHMVTypes}(L::CholmodFactor{T},B::VecOrMat{T}) = solve(L,CholmodDense!(B),CHOLMOD_A).mat
function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti})
    solve(L,B,CHOLMOD_A)
end
function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti})
    sparse!(solve(L,CholmodSparse(B),CHOLMOD_A))
end

function A_mul_Bt{Tv<:Union(Float32,Float64),Ti<:CHMITypes}(A::CholmodSparse{Tv,Ti},
                                                            B::CholmodSparse{Tv,Ti})
    A_mul_Bc(A,B) # in the unlikely event of writing A*B.' instead of A*B'
end

Ac_ldiv_B{T<:CHMVTypes}(L::CholmodFactor{T},B::CholmodDense{T}) = solve(L,B,CHOLMOD_A)
Ac_ldiv_B{T<:CHMVTypes}(L::CholmodFactor{T},B::VecOrMat{T}) = solve(L,CholmodDense!(B),CHOLMOD_A).mat
function Ac_ldiv_B{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti})
    solve(L,B,CHOLMOD_A)
end
function Ac_ldiv_B{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti})
    sparse!(solve(L,CholmodSparse(B),CHOLMOD_A))
end

function Ac_mul_B{Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::CholmodDense{Tv})
    chm_sdmult(A,true,1.,0.,B)
end
function Ac_mul_B{Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::VecOrMat{Tv})
    chm_sdmult(A,true,1.,0.,CholmodDense(B))
end

function At_mul_B{Tv<:Union(Float32,Float64),Ti<:CHMITypes}(A::CholmodSparse{Tv,Ti},
                                                            B::CholmodSparse{Tv,Ti})
    Ac_mul_B(A,B) # in the unlikely event of writing A.'*B instead of A'*B
end

cholfact{T<:CHMVTypes}(A::CholmodSparse{T},beta::T) = cholfact(A,beta,false)
cholfact(A::CholmodSparse) = cholfact(A,false)
cholfact(A::SparseMatrixCSC,ll::Bool) = cholfact(CholmodSparse(A),ll)
cholfact(A::SparseMatrixCSC) = cholfact(CholmodSparse(A),false)
function cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::CholmodSparse{T},beta::Number)
    cholfact!(L,A,convert(T,beta))
end
function cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::SparseMatrixCSC{T},beta::Number)
    cholfact!(L,CholmodSparse(A),convert(T,beta))
end
cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::SparseMatrixCSC{T}) = cholfact!(L,CholmodSparse(A))

chm_analyze(A::SparseMatrixCSC) = chm_analyze(CholmodSparse(A))

chm_print(A::CholmodSparse, lev::Integer) = chm_print(A, lev, "")
chm_print(A::CholmodFactor, lev::Integer) = chm_print(L, lev, "")

function chm_scale!{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T},typ::Integer)
    chm_scale!(A,CholmodDense(b),typ)
    A
end
chm_scale{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T},typ::Integer) = chm_scale!(copy(A),b,typ)

chm_speye(m::Integer, n::Integer) = chm_speye(m, n, 1., 1) # default element type is Float32
chm_speye(n::Integer) = chm_speye(n, n, 1.)             # default shape is square

chm_spzeros(m::Integer,n::Integer,nzmax::Integer) = chm_spzeros(m,n,nzmax,1.)

function scale!{T<:CHMVTypes}(b::Vector{T}, A::CholmodSparse{T})
    chm_scale!(A,CholmodDense(b),CHOLMOD_ROW)
    A
end
scale{T<:CHMVTypes}(b::Vector{T}, A::CholmodSparse{T}) = scale!(b,copy(A))
function scale!{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T})
    chm_scale!(A,CholmodDense(b),CHOLMOD_COL)
    A
end
scale{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T}) = scale!(copy(A), b)

norm(A::CholmodSparse) = norm(A,1)

show(io::IO,L::CholmodFactor) = chm_print(L,int32(4),"")
show(io::IO,A::CholmodSparse) = chm_print(A,int32(4),"")

size(B::CholmodDense) = size(B.mat)
size(B::CholmodDense,d) = size(B.mat,d)
size(A::CholmodSparse) = (int(A.c.m), int(A.c.n))
function size(A::CholmodSparse, d::Integer)
    d == 1 ? A.c.m : (d == 2 ? A.c.n : 1)
end
size(L::CholmodFactor) = (n = int(L.c.n); (n,n))
size(L::CholmodFactor,d::Integer) = d < 1 ? throw(BoundsError("Dimension $d out of range")) : (d <= 2 ? int(L.c.n) : 1)

function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},
                                            B::SparseMatrixCSC{Tv,Ti},typ::Integer)
    sparse!(solve(L,CholmodSparse(B),typ))
end
function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti})
    sparse!(solve(L,CholmodSparse(B),CHOLMOD_A))
end
function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti})
    solve(L,B,CHOLMOD_A)
end
function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti})
    solve(L,B,CHOLMOD_A)
end
function solve{T<:CHMVTypes}(L::CholmodFactor{T},B::Vector{T},typ::Integer)
    solve(L,CholmodDense!(B),typ).mat[:]
end
function solve{T<:CHMVTypes}(L::CholmodFactor{T},B::Matrix{T},typ::Integer)
    solve(L,CholmodDense!(B),typ).mat
end
solve{T<:CHMVTypes}(L::CholmodFactor{T},B::CholmodDense{T}) = solve(L,B,CHOLMOD_A)

function findnz{Tv,Ti}(A::CholmodSparse{Tv,Ti})
    jj = similar(A.rowval0)             # expand A.colptr0 to a vector of indices
    for j in 1:A.c.n, k in (A.colptr0[j]+1):A.colptr0[j+1]
        jj[k] = j
    end

    ind = similar(A.rowval0)
    ipos = 1
    count = 0
    for k in 1:length(A.nzval)
        if A.nzval[k] != 0
            ind[ipos] = k
            ipos += 1
            count += 1
        else
            println("warning: sparse matrix contains explicitly stored zeros")
        end
    end
    ind = ind[1:count]                  # ind is the indices of nonzeros in A.nzval
    (increment!(A.rowval0[ind]), jj[ind], A.nzval[ind])
end

findnz(L::CholmodFactor) = findnz(CholmodSparse(L))

function diag{Tv}(A::CholmodSparse{Tv})
    minmn = minimum(size(A))
    res = zeros(Tv,minmn)
    cp0 = A.colptr0
    rv0 = A.rowval0
    anz = A.nzval
    for j in 1:minmn, k in (cp0[j]+1):cp0[j+1]
        if rv0[k] == j-1
            res[j] += anz[k]
        end
    end
    res
end

function diag{Tv}(L::CholmodFactor{Tv})
    res = zeros(Tv,L.c.n)
    xv  = L.x
    if bool(L.c.is_super)
        nec = decrement!(diff(L.super))  # number of excess columns per supernode
        dstride = increment!(diff(L.pi)) # stride of diagonal elements (nrow + 1)
        px = L.px
        pos = 1
        for i in 1:length(nec)
            base = px[i] + 1
            res[pos] = xv[base]
            pos += 1
            for j in 1:nec[i]
                res[pos] = xv[base + j*dstride[i]]
                pos += 1
            end
        end
    else
        c0 = L.p
        r0 = L.i
        xv = L.x
        for j in 1:L.c.n
            jj = c0[j]+1
            assert(r0[jj] == j-1)
            res[j] = xv[jj]
        end
    end
    res
end

function logdet{Tv,Ti}(L::CholmodFactor{Tv,Ti},sub)
    res = zero(Tv)
    for d in diag(L)[sub] res += log(abs(d)) end
    bool(L.c.is_ll) ? 2res : res
end
logdet(L::CholmodFactor) = logdet(L, 1:L.c.n)
det(L::CholmodFactor) = exp(logdet(L))

full(A::CholmodSparse) = CholmodDense(A).mat
full(A::CholmodDense) = A.mat
function sparse(A::CholmodSparse)
    if bool(A.c.stype) return sparse!(copysym(A)) end
    SparseMatrixCSC(A.c.m, A.c.n, increment(A.colptr0), increment(A.rowval0), copy(A.nzval))
end
function sparse!(A::CholmodSparse)
    SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval)
end
sparse(L::CholmodFactor) = sparse!(CholmodSparse(L))
sparse(D::CholmodDense) = sparse!(CholmodSparse(D))
sparse(T::CholmodTriplet) = sparse!(CholmodSparse(T))
end #module

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "MHD Multicomponent" begin

  @trixi_testset "elixir_mhdmulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [4.30029920e-02, 4.29875057e-02, 2.57471806e-02, 1.62185617e-01, 1.74536934e-02, 1.74545523e-02, 2.68731904e-02, 1.36464770e-15, 1.21243408e-02, 2.42486817e-02],
      linf = [3.13715220e-01, 3.03783978e-01, 2.15002288e-01, 9.04249573e-01, 9.39809810e-02, 9.47028202e-02, 1.52772540e-01, 8.24570183e-15, 7.87460541e-02, 1.57492108e-01])
  end

  @trixi_testset "elixir_mhdmulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
      l2   = [4.25140100e-02, 4.25084765e-02, 2.38571816e-02, 1.15557310e-01, 1.63671989e-02, 1.63673846e-02, 2.58192976e-02, 2.34524850e-04, 1.08353280e-02, 2.16706560e-02],
      linf = [2.34580322e-01, 2.34685766e-01, 1.19133412e-01, 5.33372522e-01, 6.17688441e-02, 6.17909503e-02, 9.60217532e-02, 4.44580923e-03, 6.18887680e-02, 1.23777536e-01])
  end

  @trixi_testset "elixir_mhdmulti_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_convergence.jl"),
      l2   = [0.0003827494352191891, 0.0003827494352191414, 0.0005161253331853331, 0.0005704446900280551, 0.000443264790332803, 0.0004432647903328268, 0.0005073530103410933, 0.00038526540117132834, 7.438157880192512e-5, 0.00014876315760385024, 0.0002975263152077005],
      linf = [0.0013927512912668982, 0.0013927512912661499, 0.002597102723833799, 0.0016567175206610996, 0.0020103639634265758, 0.002010363963426798, 0.0025762126248859716, 0.0011800383681934948, 0.00028863972195533605, 0.0005772794439106721, 0.0011545588878213442])
  end

  @trixi_testset "elixir_mhdmulti_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_rotor.jl"),
      l2   = [7.17665548e-01, 7.22954510e-01, 0.00000000e+00, 8.34512452e-01, 4.74961443e-02, 7.91093545e-02, 0.00000000e+00, 1.98437701e-03, 1.91605624e-01, 9.58028120e-02],
      linf = [1.11390901e+01, 1.10181050e+01, 0.00000000e+00, 1.24464832e+01, 6.35439011e-01, 1.06796978e+00, 0.00000000e+00, 6.41596561e-02, 3.49017111e+00, 1.74508556e+00],
      tspan = (0.0, 0.01))
end

end

end # module

ext: jl


text:
using HackerNewsAnalysis

fetch_top_stories()
#=
julia> fetch_top_stories()
484-element JSON3.Array{Int64,Base.CodeUnits{UInt8,JSON3.VectorString{Array{UInt8,1}}},Array{UInt64,1}}:
 21676252
 21676933
 21677389
 21674752
 21675456
 21676027
        ⋮
=#

fetch_top_stories() |> first
#=
julia> fetch_top_stories() |> first
21676252
=#

first(fetch_top_stories())
#=
julia> first(fetch_top_stories())
21676252
=#

fetch_top_stories() |> first |> fetch_story
#=
Story("kkm", 384, 538, 1575218518, 21676252, "The world needs more search engines", [21677682, 21676507, 21676549, 21676457, 21676352, 21677053, 21676790, 21676616, 21676315, 21676385  …  21676442, 21676505, 21676858, 21676598, 21676736, 21677029, 21676919, 21676853, 21676397, 21676703], "https://www.0x65.dev/blog/2019-12-01/the-world-needs-cliqz-the-world-needs-more-search-engines.html")
=#


text:
__precompile__()

module Structs

export IntP, FloatP
export DmdStruct

const IntP = Int64
const FloatP = Float64

mutable struct DmdStruct
    facecon::Array{IntP,3}; # faces-to-DG nodes connectivities (used to obtain DG nodes on faces)
    eblks::Array{IntP,2};  # blocks of elements for parallel computation
    fblks::Array{IntP,2};  # blocks of faces for parallel computation
    nbsd::Array{IntP,2};   # neighboring subdomains (neighbors)
    elemsend::Array{IntP,2};  # elements to be sent to neighbors
    elemrecv::Array{IntP,2};  # elements to be received from neighbors
    elemsendpts::Array{IntP,2}; # markers for elements to be sent to neighbors
    elemrecvpts::Array{IntP,2}; # markers for elements to be received from neighbors
    elempart::Array{IntP,2};    # element partitions
    elempartpts::Array{IntP,2}; # markers for element partitions
    facepart::Array{IntP,2};    # element partitions
    facepartpts::Array{IntP,2}; # markers for element partitions
    facepartbnd::Array{IntP,2};
    t2f::Array{IntP,2};
    f::Array{IntP,2};
    #elem2cpu::Array{IntP,2};    # element partitions
    DmdStruct() = new();
end

# mutable struct AppStruct
#     codename::String; # Exasim
#     version::String;  # Set version
#     appname::String;  # application name
#     platform::String; # CPU or GPU
#     cpucompiler::String; # Path to CPU compiler
#     mpicompiler::String; # Path to MPI compiler
#     gpucompiler::String; # Path to GPU compiler
#     mpirun::String;      # Path to MPI run command and MPI run options
#     cpuflags::String;    # options for CPU compiler
#     gpuflags::String;    # options for GGU compiler
#     # Flux::String;
#     # Source::String;
#     # Mass::String;
#     # Avfield::String;
#     # Output::String;
#     # Ubou::String;
#     # Fbou::String;
#     # Uhat::String;
#     # Fhat::String;
#     # Stab::String;
#     pdemodel::String;# used to indicate PDE model
#     pdemodelfile::String;# PDE model file name
#
#     preprocessmode::IntP; # preprocessing mode
#     mpiprocs::IntP; # number of MPI ranks
#     nd::IntP; # physical dimension
#     nc::IntP; # number of compoments of (u, q, p)
#     ncu::IntP;# number of compoments of (u)
#     ncq::IntP;# number of compoments of (q)
#     ncp::IntP;# number of compoments of (p)
#     nco::IntP;# number of compoments of (odg)
#     nch::IntP;# number of compoments of (uhat)
#     ncx::IntP;# number of compoments of (xdg)
#     ncw::IntP;# number of compoments of (wdg)
#     nce::IntP;# number of compoments of (output field)
#     neb::IntP;# number of element blocks for parallel computation
#     nfb::IntP;# number of face blocks for parallel computation
#     elemtype::IntP; # type of elements
#     nodetype::IntP; # type of nodes
#     hybrid::IntP; # discretization method
#     tdep::IntP; # flag for steady-state or time-dependent problem
#     wave::IntP; # flag for wave problem
#     linearproblem::IntP; # flag for linear problem
#     debugmode::IntP; # flag for debug mode
#     stgNmode::IntP; # number of synthetic turbulence generation modes
#
#     porder::IntP; # polymnomial degree
#     pgauss::IntP; # Gauss quadrature polynomial degree
#     temporalscheme::IntP; # temporal scheme
#     torder::IntP; # temporal order of accuracy
#     nstage::IntP; # number of RK stages
#     convStabMethod::IntP; # flag for convection stabilization method
#     diffStabMethod::IntP; # flag for diffusion stabilization method
#     rotatingFrame::IntP;  # flag for rotating frame
#     viscosityModel::IntP; # flag for viscosity model
#     SGSmodel::IntP; # flag for SGS model
#     ALE::IntP; # flag for ALE
#     AV::IntP; # flag for artificial viscosity
#     nonlinearsolver::IntP; # flag for nonlinear solver (Newton default)
#     linearsolver::IntP; # flag for linear solver (GMRES default)
#     NLiter::IntP; # maximum number of nonlinear iterations
#     linearsolveriter::IntP; # maximum number of linear iterations
#     GMRESrestart::IntP; # number of GMRES restarts
#     GMRESortho::IntP; # GMRES orthogonalization method
#     preconditioner::IntP; # flag for preconditioner method
#     precMatrixType::IntP; # flag for type of preconditioner matrix
#     NLMatrixType::IntP;
#     runmode::IntP; # flag for run mode
#     tdfunc::IntP; # flag for time-dependent function associated with time-derivative
#     source::IntP; # flag for source-term function
#
#     matvecorder::IntP; # flag for order of accuracy for matrix-vector multiplication
#     RBdim::IntP; # reduced-basis dimension
#     saveSolFreq::IntP; # flag for how often the solution saved in binary files
#     saveSolOpt::IntP; # option for how the solution be saved: 0 -> u only, 1 -> u and q
#     timestepOffset::IntP; # for restarting the simulation from the saved solution
#     saveSolBouFreq::IntP; # how often the solution be saved on a particular boundary
#     ibs::IntP; # the boundary on which the solution be saved
#     compudgavg::IntP; # flag if time-average solution is computed
#     extFhat::IntP;
#     extUhat::IntP;
#
#     time::FloatP; # starting time (usually 0, however >0 if restarting from the saved solution)
#     NLparam::FloatP;
#     NLtol::FloatP; # nonlinear solver tolerance
#     linearsolvertol::FloatP; # linear solver tolerance
#     matvectol::FloatP; # matrix-vector product tolerance
#
#     flag::Array{IntP,2};   # flag parameters
#     problem::Array{IntP,2};# problem parameters
#     boundaryconditions::Array{IntP,2};# a list of boundary condition numbers
#     stgib::Array{IntP,2};  # synthetic turbulence
#
#     uinf::Array{FloatP,2};    # boundary data
#     dt::Array{FloatP,2};      # time steps
#     factor::Array{FloatP,2};  # factors
#     physicsparam::Array{FloatP,2}; # physical parameters
#     solversparam::Array{FloatP,2}; # solvers parameters
#     tau::Array{FloatP,2}; # stabilization parameters
#     stgdata::Array{FloatP,2}; # synthetic turbulence
#     stgparam::Array{FloatP,2}; # synthetic turbulence
#
#     gmsh::String;
#     metis::String;
#     paraview::String;
#     visfilename::String;
#     visscalars;
#     visvectors;
#     viselem;
#     AppStruct() = new();
# end
#
# function InitializeAppStruct(app,version)
#     app.codename = "Exasim";
#     app.version = version;
#     app.appname = "app";
#     app.platform = "cpu";
#     app.cpucompiler = "g++";
#     app.mpicompiler = "mpicxx";
#     app.gpucompiler = "nvcc";
#     app.mpirun = "mpirun";
#     #app.cpuflags = "-O2 -Wall -ldl -lm -lblas -llapack";
#     app.cpuflags = "-O2 -ldl -lm -lblas -llapack";
#     app.gpuflags = "-lcudart -lcublas";
#     app.pdemodelfile = "";
#     # app.Source = "";
#     # app.Mass = "";
#     # app.Avfield = "";
#     # app.Output = "";
#     # app.Ubou = "";
#     # app.Fbou = "";
#     # app.Uhat = "";
#     # app.Fhat = "";
#     # app.Stab = "";
#     app.preprocessmode = 1;
#     app.mpiprocs = 1;
#     app.nd = 1;
#     app.nc = 1;
#     app.ncu = 1;
#     app.ncq = 0;
#     app.ncp = 0;
#     app.nco = 0;
#     app.nch = 1;
#     app.ncx = 1;
#     app.ncw = 0;
#     app.nce = 0;
#     app.neb = 512*8;
#     app.nfb = 512*32;
#     app.elemtype = 1;
#     app.nodetype = 1;
#     app.pdemodel="ModelD";
#     app.hybrid = 0;
#     app.tdep = 0;
#     app.wave = 0;
#     app.linearproblem = 0;
#     app.debugmode = 0;
#     app.stgNmode = 0;
#     app.porder = 1;
#     app.pgauss = 2;
#     app.temporalscheme = 0;
#     app.torder = 1;
#     app.nstage = 1;
#     app.convStabMethod = 0;
#     app.diffStabMethod = 0;
#     app.rotatingFrame = 0;
#     app.viscosityModel = 0;
#     app.SGSmodel = 0;
#     app.ALE = 0;
#     app.AV = 0;
#     app.nonlinearsolver = 0;
#     app.linearsolver = 0;
#     app.NLiter = 20;
#     app.linearsolveriter = 200;
#     app.GMRESrestart = 25;
#     app.GMRESortho = 0;
#     app.preconditioner = 0;
#     app.precMatrixType = 0;
#     app.NLMatrixType = 0;
#     app.runmode = 0;
#     app.tdfunc = 1;
#     app.source = 1;
#     app.matvecorder = 1;
#     app.RBdim = 5;
#     app.saveSolFreq = 1;
#     app.saveSolOpt = 1;
#     app.timestepOffset = 0;
#     app.saveSolBouFreq = 0;
#     app.ibs = 0;
#     app.compudgavg = 0;
#     app.extFhat = 0;
#     app.extUhat = 0;
#
#     app.time = 0.0;
#     app.NLparam = 0.0;
#     app.NLtol = 1e-6;
#     app.linearsolvertol = 1e-3;
#     app.matvectol = 1e-6;
#
#     app.flag = [0 0];
#     app.problem = [0 0];
#     app.boundaryconditions = [0 0];
#     app.stgib = [0 0];
#
#     app.tau = reshape([1.0],1,1); # stabilization parameters
#     app.uinf = [0.0 0.0]; # freestream values
#     app.dt = [0.0 0.0];  # time steps
#     app.factor = [0.0 0.0];  # factors
#     app.physicsparam = [0.0 0.0]; # physical parameters
#     app.solversparam = [0.0 0.0]; # solvers parameters
#     app.stgdata = [0.0 0.0]; # synthetic turbulence
#     app.stgparam = [0.0 0.0]; # synthetic turbulence
#
#     app.gmsh = "";
#     app.metis = "";
#     app.paraview = "";
#     app.visfilename = "";
#     app.visscalars = [];
#     app.visvectors = [];
#     app.viselem = [];
#     return app;
# end
#
# struct MasterStruct
#     nd::IntP; # physical dimension
#     porder::IntP; # polynomial degree
#     pgauss::IntP; # order of Gauss quadrature
#     elemtype::IntP; # type of elements
#     nodetype::IntP; # type of nodes
#     npe::IntP; # number of polynomials per element
#     npf::IntP; # number of polynomials per face
#     nge::IntP; # number of Gauss points per element
#     ngf::IntP; # number of Gauss points per face
#     np1d::IntP; # number of polynomials per 1D element
#     ng1d::IntP; # number of Gauss points per 1D element
#     perm::Array{IntP,2}; # to get face nodes from element nodes
#     shapegt::Array{FloatP,3}; # element shape functions at Gauss points (transpose)
#     shapegw::Array{FloatP,3}; # element shape functions at Gauss points multiplied by Gauss weights
#     shapfgt::Array{FloatP,3}; # face shape functions at Gauss points (transpose)
#     shapfgw::Array{FloatP,3}; # face shape functions at Gauss points multiplied by Gauss weights
#     shapent::Array{FloatP,3}; # element shape functions at nodes (transpose)
#     shapen::Array{FloatP,3};  # element shape functions at nodes
#     shapfnt::Array{FloatP,3}; # face shape functions at nodes (transpose)
#     shapfn::Array{FloatP,3};  # face shape functions at nodes
#     xpe::Array{FloatP,2}; # nodal points on master element
#     gpe::Array{FloatP,2}; # gauss points on master element
#     gwe::Array{FloatP,2}; # gauss weighs on master element
#     xpf::Array{FloatP,2}; # nodal points on master face
#     gpf::Array{FloatP,2}; # gauss points on master face
#     gwf::Array{FloatP,2}; # gauss weighs on master face
#     shap1dgt::Array{FloatP,3}; # 1D shape functions at Gauss points (transpose)
#     shap1dgw::Array{FloatP,3}; # 1D shape functions at Gauss points multiplied by Gauss weights
#     shap1dnt::Array{FloatP,3}; # 1D shape functions at nodes (transpose)
#     shap1dn::Array{FloatP,3}; # 1D shape functions at nodes
#     xp1d::Array{FloatP,2}; # node points on 1D element
#     gp1d::Array{FloatP,2}; # gauss points on 1D element
#     gw1d::Array{FloatP,2}; # gauss weights on 1D element
#     telem::Array{IntP,2}; # element
#     tface::Array{IntP,2}; # face
# end

# mutable struct MeshStruct
#     p::Array{FloatP,2};       # points of a linear mesh
#     t::Array{IntP,2};         # elements of a linear mesh
#     f::Array{IntP,2};         # faces of a linear mesh
#     dgnodes::Array{FloatP,3}; # spatial nodes of a high-order mesh
#     tprd::Array{IntP,2};      # elements for periodic conditions
#     boundaryexpr;             # expressions to determine boundaries
#     periodicexpr;             # expressions to determine periodic boundaries
#     curvedboundary::Array{IntP,2}; # flags for curved boundaries
#     curvedboundaryexpr;       # expressions to determine curved boundaries
#     MeshStruct() = new();
# end
#
# mutable struct SolStruct
#     UDG::Array{FloatP,3}; # intial solution (u, q) at DG nodes
#     ODG::Array{FloatP,3}; # auxilary field at DG nodes
#     WDG::Array{FloatP,3}; # dw/dt = u (wave problem) at DG nodes
#     Uout::Array{FloatP,3}; # numerical solution (u, q, w)
#     SolStruct() = new();
# end

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
a = parse(Float64, readline(stdin))
b = parse(Float64, readline(stdin))
println(stdout, "Got values: $a, $b")
if b > a
    println(stderr, "Wrong values: ", b, ">", a)
    exit(1)
end
println(stdout, "log(", a, "-", b, ")=", log(a-b))
exit(0)

max_stars_repo_path: Chapter03/01. Managing streams, and reading and writing files/iotest.jl


text:
function bs(val::Float64)
	if isfinite(val)
		out = @sprintf("% .3e", val)
	elseif isnan(val)
		out = "       NaN"
	else
		out = "       Inf"
	end
	out
end

function createBLOracles(trf,trl,numTP, tef, tel, L2reg::Bool; outputLevel::Int64=0)

	outputLevel > 0  && println("Fraction of ones in training set: $(trl.numPlus/length(trl))")
	outputLevel > 0  && println("Fraction of ones in testing  set: $(tel.numPlus/length(tel))")

	trft=trf'
	gradientOracle(W,index) = BL_get_f_g(trft, trl,W,index)
	gradientOracle(W) = BL_get_f(trf, trl,W)
	testFunction(W) = BL_for_output(tef, tel,W)

	function outputsFunction(W)
		(f,pcc, mcc, tp, tn, fp, fn) = testFunction(W)
		yo = gradientOracle(W)
		ResultFromOO(bs(yo)*" "*bs(f)*" "*bs(pcc)*" "*bs(mcc),[f, pcc, mcc, tp, tn, fp, fn, yo])
	end

	restoreGradient(cs,indices) = BL_restore_gradient(trft, trl,cs,indices)

	if L2reg
		mygradientOracle(a) = L2regGradient(gradientOracle, 1/numTP, a)
		mygradientOracle(a, b) = L2regGradient(gradientOracle, 1/numTP, a, b)
		csDataType = Matrix{Float64}
	else
		mygradientOracle=gradientOracle
		csDataType = Matrix{Float64}
	end

	(mygradientOracle, size(trf)[2], numTP, csDataType, "BL", Outputter(outputsFunction,  "     f-train       f         pcc        mcc",8), L2reg)
end

function createMLOracles(trf,trl,numTP,numClasses, tef, tel, L2reg::Bool; outputLevel::Int64=0)
	trft=trf'
	gradientOracle(W,index) = ML_get_f_g(trft, trl,W,index)
	gradientOracle(W) = ML_get_f(trf, trl,W)
	testFunction(W) = ML_for_output(tef, tel,W)

	function outputsFunction(W)
		(f,pcc) = testFunction(W)
		yo = gradientOracle(W)
		ResultFromOO(bs(yo)*" "*bs(f) *" "*bs(pcc), [f, pcc, yo])
	end

	restoreGradient(cs,index) = ML_restore_gradient(trft, trl,cs,index)
	if L2reg
		mygradientOracle(a) = L2regGradient(gradientOracle, 1/numTP, a)
		mygradientOracle(a, b) = L2regGradient(gradientOracle, 1/numTP, a, b)
		csDataType=Matrix{Float64}
	else
		mygradientOracle=gradientOracle
		myrestoreGradient=restoreGradient
		csDataType=Matrix{Float64}
	end

	(mygradientOracle, size(trf)[2]*numClasses,  numTP, csDataType, "ML", Outputter(outputsFunction,  "     f-train       f         pcc    ",3), L2reg,restoreGradient)
end

function createSBLOracles(trf,trl,numTP, tef, tel, L2reg::Bool; outputLevel::Int64=0)

	outputLevel > 0  && println("Fraction of ones in training set: $(trl.numPlus/length(trl))")
	outputLevel > 0  && println("Fraction of ones in testing  set: $(tel.numPlus/length(tel))")

	trft=trf'
	gradientOracle(W,index) = SBL_get_f_g(trft, trl,W,index)
	gradientOracle(W) = SBL_get_f(trf, trl,W)
	testFunction(W) = SBL_for_output(tef, tel,W)

	function outputsFunction(W)
		(f,pcc, mcc, tp, tn, fp, fn) = testFunction(W)
		yo = gradientOracle(W)
		ResultFromOO(bs(yo)*" "*bs(f)*" "*bs(pcc)*" "*bs(mcc),[f, pcc, mcc, tp, tn, fp, fn, yo])
	end

	restoreGradient(cs,indices) = SBL_restore_gradient(trft, trl,cs,indices)

	if L2reg
		mygradientOracle(a) = L2regGradient(gradientOracle, 1/numTP, a)
		mygradientOracle(a, b) = L2regGradient(gradientOracle, 1/numTP, a, b)
		csDataType = Matrix{Float64}
	else
		mygradientOracle=gradientOracle
		csDataType = SparseMatrixCSC{Float64,Int64}
	end

	(mygradientOracle, size(trf)[2], numTP, csDataType, "SBL", Outputter(outputsFunction,  "    f-train       f         pcc        mcc",8), L2reg)
end


text:
struct FileLocInfo
    filepath::String
    line::Int
    # 0 if unknown
    column::Int
    # The line at which the current context starts, 0 if unknown
    defline::Int
    # typemax(int) if unknown
    endline::Int
end

struct BufferLocInfo
    data::String
    line::Int
    # 0 if unknown
    column::Int
    defline::Int
    # typemax(int) if unknown
    endline::Int
end

function loc_for_fname(file::String, line::Integer, defline::Integer, endline::Integer)
    if startswith(file, "REPL[")
        hist_idx = parse(Int,string(file)[6:end-1])
        isdefined(Base, :active_repl) || return nothing, ""
        hp = Base.active_repl.interface.modes[1].hist
        return BufferLocInfo(hp.history[hp.start_idx+hist_idx], line, 0, defline, endline)
    else
        for path in SEARCH_PATH
            fullpath = joinpath(path,string(file))
            if isfile(fullpath)
                return FileLocInfo(fullpath, line, 0, defline, endline)
            end
        end
    end
    return nothing
end

function locinfo(frame::Frame)
    if frame.framecode.scope isa Method
        meth = frame.framecode.scope
        def_file, def_line = JuliaInterpreter.whereis(meth)
        current_file, current_line = JuliaInterpreter.whereis(frame)
        n_stmts = length(frame.framedata.ssavalues)
        total_lines = JuliaInterpreter.linenumber(frame, 1n_stmts) - def_line + 1
        # We currently cannot see a difference between f(x) = x and
        # function f(x)
        #     x
        # end
        # If we could, we would do the += 1 below only in the second case
        # This means that we now miss printing the end for cases like the second
        # (one line bodies)
        total_lines == 1 || (total_lines += 1)
        end_line = def_line + total_lines - 1
        if def_file != current_file || def_line > current_line
            def_line = 0 # We are not sure where the context start in cases like these, could be improved?
            end_line = typemax(Int)
        end
        return loc_for_fname(current_file, current_line, def_line, end_line)
    else
        println("not yet implemented")
    end
end

# Used for the tests
const _print_full_path = Ref(true)

function locdesc(frame::Frame)
    sprint() do io
        if frame.framecode.scope isa Method
            locdesc(io, frame.framecode)
        else
            println(io, "not yet implemented")
        end
    end
end

function locdesc(io, framecode::FrameCode)
    meth = framecode.scope
    argnames = framecode.src.slotnames[2:meth.nargs]
    spectypes = Any[Any for i=1:length(argnames)]
    print(io, meth.name,'(')
    first = true
    for (argname, argT) in zip(argnames, spectypes)
        first || print(io, ", ")
        first = false
        print(io, argname)
        !(argT === Any) && print(io, "::", argT)
    end
    path = _print_full_path[] ? meth.file : string(basename(String(meth.file)), ":", meth.line)
    path = CodeTracking.replace_buildbot_stdlibpath(String(path))
    print(io, ") at ", path)
end


text:
using JuLIP, StaticArrays, LinearAlgebra
using ACE: State, filter
using JuLIP.Potentials: neigsz
# using ACE: BondEnvelope, filter, State, CylindricalBondEnvelope

# TODO: make this type-stable
struct BondsIterator
   at
   nlist_bond
   nlist_env
   filter
end

"""
* rcutbond: include all bonds (i,j) such that rij <= rcutbond
* `rcutenv`: include all bond environment atoms k such that `|rk - mid| <= rcutenv`
* `filter` : `filter(X) == true` if particle `X` is to be included; `false` if to be discarded from the environment
"""
function bonds(at::Atoms, rcutbond, rcutenv, filter)
   nlist_bond = neighbourlist(at, rcutbond; recompute=true, storelist=false)
   nlist_env = neighbourlist(at, rcutenv; recompute=true, storelist=false)
   return BondsIterator(at, nlist_bond, nlist_env, filter)
end

function Base.iterate(iter::BondsIterator, state=(1,0))
   i, q = state
   # store temporary arrays for those...
   Js, Rs = neigs(iter.nlist_bond, i)

   # nothing left to do
   if i >= length(iter.at) && q >= length(Js)
      return nothing
   end

   # increment:
   #   case 1: we haven't yet exhausted the current neighbours.
   #           just increment the q index pointing into Js
   if q < length(Js)
      q += 1
      # here we could build in a rule to skip any pair for which we don't
      # want to do the computation.

    #  case 2: if i < length(at) but q >= length(Js) then we need to
    #          increment the i index and get the new neighbours
   elseif i < length(iter.at)
      i += 1
      Js, Rs = neigs(iter.nlist_bond, i)
      q = 1
   else
      return nothing
   end

   j = Js[q]   # index of neighbour (in central cell)
   rr0 = rrij = Rs[q]  # position of neighbour (in shifted cell) relative to i
   # ssj = Rs[q] - iter.at.X[j]   # shift of atom j into shifted cell

   # now we construct the environment
   Js_e, Rs_e, Zs_e = _get_bond_env(iter, i, j, rrij)

   return (i, j, rrij, Js_e, Rs_e, Zs_e), (i, q)
end

function _get_bond_env(iter::BondsIterator, i, j, rrij)
   # TODO: store temporary arrays
   Js_i, Rs_i, Zs_i = neigsz(iter.nlist_env, iter.at, i)

   rri = iter.at.X[i]
   rrmid = rri + 0.5 * rrij
   Js = Int[]; sizehint!(Js,  length(Js_i) ÷ 4)
   Rs = typeof(rrij)[]; sizehint!(Rs,  length(Js_i) ÷ 4)
   Zs = AtomicNumber[]; sizehint!(Zs,  length(Js_i) ÷ 4)

   ŝ = rrij/norm(rrij)

   # find the bond and remember it;
   # TODO: this could now be integrated into the second loop
   q_bond = 0
   for (q, rrq) in enumerate(Rs_i)
      # rr = rrq + rri - rrmid
      if rrq ≈ rrij   # TODO: replace this with checking for j and shift!
         @assert Js_i[q] == j
         q_bond = q
         break
      end
   end
   if q_bond == 0
      error("the central bond neigbour atom j was not found")
   end

   # now add the environment
   for (q, rrq) in enumerate(Rs_i)
      # skip the central bond
      if q == q_bond; continue; end
      # add the rest provided they fall within the provided filter
      rr = rrq + rri - rrmid
      z = dot(rr, ŝ)
      r = norm(rr - z * ŝ)
      if iter.filter(r, z)
         push!(Js, Js_i[q])
         push!(Rs, rr)
         push!(Zs, Zs_i[q])
      end
   end

   return Js, Rs, Zs
end


text:
using OpenBCI

boardIP = "192.168.1.2"  // change for your board
myIP = "192.168.1.1"     // change for your host
idfile = "patientdata.json"

OpenBCI.makeganglionbdfplus("testfile.bdf", boardIP, myIP, 30, idfile=idfile)

true


text:
# General idea here: An ILPOptimizer is a lightweight type used to create the a full model.
# An `ILPHolder` is the struct that holds most of the metadata about ILP construction.
# An `ILPOptimizer` is used to construct an `ILPHolder`.

# ILP Optimizer
abstract type ILPType end

struct Static <: ILPType end
struct Synchronous <: ILPType end
struct Asynchronous <: ILPType end

struct ILPOptimizer{T, U <: ILPType} <: AbstractOptimizer{T}
    # PMM to DRAM ratio - depending on `T`
    x::T
end

# Hijack `_getratio` to forward to the correct field.
_getratio(I::ILPOptimizer) = I.x

function (O::ILPOptimizer{Rational{Int64}})(data, backend::nGraph.Backend)
    bounds = Profiler.allocation_bounds(data)

    top = max(
        bounds.upper_bound / (getratio(O) + 1) / 1E6,
    )

    x = fill(round(Int, top), size(nodes(data)))
    println("Trying to use $(maximum(x)) MB of memory")
    return ILPHolder(O, x, backend; defrag = !iszero(_numerator(O)))
end

function (O::ILPOptimizer{Int64})(data, backend::nGraph.Backend)
    x = fill(round(Int, getlimit(O) / 1E6), size(nodes(data)))
    println("Trying to use $(maximum(x)) MB of memory")
    return ILPHolder(O, x, backend)
end

Static(x::T) where {T <: OPTIM_TYPES} = ILPOptimizer{T, Static}(x)
Synchronous(x::T) where {T <: OPTIM_TYPES} = ILPOptimizer{T, Synchronous}(x)
Asynchronous(x::T) where {T <: OPTIM_TYPES} = ILPOptimizer{T, Asynchronous}(x)

#####
##### TensorMeta
#####

struct TensorMeta
    graph::MetaGraph

    # Nodes using this tensor
    users::Vector{XNode}

    # Look-up a node wrapper, get the node that serves as a reference for this
    reference_map::Dict{XNode, XNode}

    # Flag to indicate that this is a fixed tensor - can decrease variable generation.
    isfixed::Bool
end

#####
##### ILPHolder
#####

mutable struct ILPHolder{T <: ILPType}
    dram_limits::Vector{Int}

    descriptors::Dict{XTensor{XNode}, TensorMeta}
    async_move_vars::Dict{XNode, Vector{JuMP.VariableRef}}

    # Need to key this with string names instead of XNode because of recompiling the ngraph
    # function and reporifling.
    #
    # New XNodes get created that don't hash the same as the old ones.
    node_to_limit_index::Dict{String, Int}

    # Bandwidths
    read_bandwidth::Int64
    write_bandwidth::Int64
    read_bandwidth_async::Int64
    write_bandwidth_async::Int64

    # Flag to determine if we need to defrag
    defrag::Bool
end

# Build an ILP Holder that will constain most of the metadata for IL: Construction.
function ILPHolder(O::ILPOptimizer{<:Any, T}, limits, backend; defrag = true) where {T}
    return ILPHolder{T}(
        limits,
        Dict{XTensor{XNode}, TensorMeta}(),
        Dict{XNode, Vector{JuMP.VariableRef}}(),
        Dict{String, Int}(),
        _bw_remote_local_sync(backend),
        _bw_local_remote_sync(backend),
        _bw_remote_local_async(backend),
        _bw_local_remote_async(backend),
        defrag,
    )
end

name(::Type{Static}) = "static"
name(::Type{Synchronous}) = "synchronous"
name(::Type{Asynchronous}) = "asynchronous"
name(::ILPOptimizer{<:Any, T}) where {T} = name(T)

# Struct to be passed around since all these items are generally used together anyways.
mutable struct Frame{T <: ILPHolder}
    modeltype::T
    model::JuMP.Model
    profile_data::FunctionData

    # Function parameters that are in the "local" memory and thus should be counted for the
    # total memory consumption.
    local_args::Vector{XTensor{XNode}}
end

Frame(modeltype, model::JuMP.Model, profile_data::FunctionData) =
    Frame(modeltype, model, profile_data, XTensor{XNode}[])

limit(F::Frame) = limit(F.modeltype)
JuMP.optimize!(F::Frame) = optimize!(F.model)

#####
##### Model Types
#####

# Implementations
include("tensor_graphs.jl")
include("formulation.jl")
include("defrag.jl")
include("configure.jl")
include("inspect.jl")
include("compare.jl")

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
function simulated_annealing(fcn::Function,
                             x0::Array,
                             args...;
                             neighbor!::Function   = Optim.default_neighbor!,
                             iterations::Int       = 1000,
                             store_trace::Bool     = false,
                             show_trace::Bool      = false,
                             extended_trace::Bool  = false,
                             temperature::Function = Optim.log_temperature,
                             kwargs...)

    Optim.optimize(fcn, x0,
                   method = SimulatedAnnealing(neighbor = neighbor!,temperature = temperature),
                   iterations = iterations, store_trace = store_trace, show_trace = show_trace,
                   extended_trace = extended_trace)

end

function log_temperature(t::Real; initial_temperature::Real = 1.0)
    return initial_temperature/log(1+t)
end

function exponential_temperature(t::Real; initial_temperature::Real = 1.0, α::Real = .99)
    if !( 0.0 < α < 1.0)
        error("α must be in (0,1)")
    end
    return initial_temperature*α^t
end

function linear_temperature(t::Real;  initial_temperature::Real = 1, β = .005)
    return max(initial_temperature - t*β, 0.0)
end


text:
@testset "Test NLP Evals" begin

    x0 = zeros(2)

    nlp = ADNLPModel(x -> norm(x.^2), x0)
    max_nlp = init_max_counters(obj = 2)
    nlp_stop_evals = NLPStopping(nlp, max_cntrs = max_nlp)

    slack = SlackModel(nlp)
    max_slack = init_max_counters(obj = 2)
    slack_stop_evals = NLPStopping(nlp, max_cntrs = max_slack)

    nls = ADNLSModel(x -> x.^2, x0, 2)
    max_nls = init_max_counters_NLS(obj = 2, residual = 1)
    nls_stop_evals = NLPStopping(nls, max_cntrs = max_nls)

    @test typeof(nlp_stop_evals.pb.counters) == Counters
    @test nlp_stop_evals.meta.max_cntrs[:neval_obj] == 2
    @test typeof(nls_stop_evals.pb.counters) == NLSCounters
    @test nls_stop_evals.meta.max_cntrs[:neval_obj] == 2
    @test nls_stop_evals.meta.max_cntrs[:neval_residual] == 1

    Stopping._resources_check!(nlp_stop_evals, x0)
    @test nlp_stop_evals.meta.resources == false
    Stopping._resources_check!(slack_stop_evals, x0)
    @test slack_stop_evals.meta.resources == false
    Stopping._resources_check!(nls_stop_evals, x0)
    @test nls_stop_evals.meta.resources == false

    nlp_stop_evals.meta.max_cntrs[:neval_obj] = 0
    Stopping._resources_check!(nlp_stop_evals, x0)
    @test nlp_stop_evals.meta.resources == false
    nls_stop_evals.meta.max_cntrs[:neval_residual] = -1
    Stopping._resources_check!(nls_stop_evals, x0)
    @test nls_stop_evals.meta.resources == true

    max_evals!(nlp_stop_evals, 10)
    @test nlp_stop_evals.meta.max_cntrs[:neval_obj] == 10
    @test nlp_stop_evals.meta.max_cntrs[:neval_sum] == typemax(Int)

    max_evals!(nlp_stop_evals, allevals = 10, obj = 2)
    @test nlp_stop_evals.meta.max_cntrs[:neval_cons] == 10
    @test nlp_stop_evals.meta.max_cntrs[:neval_sum]  == typemax(Int)
    @test nlp_stop_evals.meta.max_cntrs[:neval_obj] == 2

end


text:
abstract type AbstractMPS end

"""
    length(::MPS/MPO)

The number of sites of an MPS/MPO.
"""
Base.length(m::AbstractMPS) = m.length

"""
    ITensors.data(::MPS/MPO)

The Vector storage of an MPS/MPO.

This is mostly for internal usage, please let us
know if there is functionality not available for
MPS/MPO you would like.
"""
data(m::AbstractMPS) = m.data

leftlim(m::AbstractMPS) = m.llim

rightlim(m::AbstractMPS) = m.rlim

function setleftlim!(m::AbstractMPS, new_ll::Int)
  m.llim = new_ll
end

function setrightlim!(m::AbstractMPS, new_rl::Int)
  m.rlim = new_rl
end

isortho(m::AbstractMPS) = leftlim(m)+1 == rightlim(m)-1

function orthocenter(m::T) where {T<:AbstractMPS}
  !isortho(m) && error("$T has no well-defined orthogonality center")
  return leftlim(m)+1
end

Base.getindex(M::AbstractMPS,
              n::Integer) = getindex(data(M),n)

function Base.setindex!(M::AbstractMPS,
                        T::ITensor,
                        n::Integer)
  (n <= leftlim(M)) && setleftlim!(M,n-1)
  (n >= rightlim(M)) && setrightlim!(M,n+1)
  setindex!(data(M),T,n)
end

Base.copy(m::AbstractMPS) = typeof(m)(length(m),
                                      copy(data(m)),
                                      leftlim(m),
                                      rightlim(m))

Base.similar(m::AbstractMPS) = typeof(m)(length(m),
                                         similar(data(m)),
                                         0,
                                         length(m))

Base.deepcopy(m::AbstractMPS) = typeof(m)(length(m),
                                          deepcopy(data(m)),
                                          leftlim(m),
                                          rightlim(m))

Base.eachindex(m::AbstractMPS) = 1:length(m)

Base.iterate(M::AbstractMPS) = iterate(data(M))

Base.iterate(M::AbstractMPS, state) = iterate(data(M), state)

"""
    dag(m::MPS)

    dag(m::MPO)

Hermitian conjugation of a matrix product state or operator `m`.
"""
function dag(m::AbstractMPS)
  N = length(m)
  mdag = typeof(m)(N)
  for i in eachindex(m)
    mdag[i] = dag(m[i])
  end
  return mdag
end

"""
    prime!(m::MPS, args...; kwargs...)

    prime!(m::MPO, args...; kwargs...)

Prime all ITensors of an MPS/MPO in-place.
"""
function prime!(M::AbstractMPS, args...; kwargs...)
  for i in eachindex(M)
    M[i] = prime(M[i], args...; kwargs...)
  end
  return M
end

"""
    prime(m::MPS, args...; kwargs...)

    prime(m::MPO, args...; kwargs...)

Prime all ITensors of an MPS/MPO.
"""
prime(M::AbstractMPS, vargs...) = prime!(copy(M), vargs...)

"""
    noprime!(m::MPS, args...; kwargs...)

    noprime!(m::MPO, args...; kwargs...)

Set the prime level to zero for all ITensors of an MPS/MPO,
in-place.
"""
function noprime!(M::AbstractMPS, args...; kwargs...)
  for i in eachindex(M)
    M[i] = noprime(M[i], args...; kwargs...)
  end
  return M
end

"""
    noprime(m::MPS, args...; kwargs...)

    noprime(m::MPO, args...; kwargs...)

Set the prime level to zero for all ITensors of an MPS/MPO.
"""
noprime(M::AbstractMPS, vargs...) = noprime!(copy(M), vargs...)

function primelinkinds!(M::AbstractMPS, plinc::Integer = 1)
  for i ∈ eachindex(M)[1:end-1]
    l = linkind(M,i)
    prime!(M[i], plinc, l)
    prime!(M[i+1], plinc, l)
  end
end

function simlinkinds!(M::AbstractMPS)
  for i ∈ eachindex(M)[1:end-1]
    isnothing(commonind(M[i],M[i+1])) && continue
    l = linkind(M,i)
    l̃ = sim(l)
    replaceind!(M[i],l,l̃)
    replaceind!(M[i+1],l,dag(l̃))
  end
end

"""
    maxlinkdim(M::MPS)

    maxlinkdim(M::MPO)

Get the maximum link dimension of the MPS or MPO.
"""
function maxlinkdim(M::AbstractMPS)
  md = 0
  for b ∈ eachindex(M)[1:end-1]
    md = max(md,dim(linkind(M,b)))
  end
  md
end

function Base.show(io::IO, M::AbstractMPS)
  print(io,"$(typeof(M))")
  (length(M) > 0) && print(io,"\n")
  for (i, A) ∈ enumerate(data(M))
    if order(A) != 0
      println(io,"[$i] $(inds(A))")
    else
      println(io,"[$i] ITensor()")
    end
  end
end

function linkind(M::AbstractMPS, j::Int)
  N = length(M)
  j ≥ length(M) && error("No link index to the right of site $j (length of MPS is $N)")
  li = commonind(M[j],M[j+1])
  if isnothing(li)
    error("linkind: no MPS link index at link $j")
  end
  return li
end

function plussers(::Type{T},
                  left_ind::Index,
                  right_ind::Index,
                  sum_ind::Index) where {T<:Array}
  total_dim    = dim(left_ind) + dim(right_ind)
  total_dim    = max(total_dim, 1)
  # TODO: I am not sure if we should be using delta
  # tensors for this purpose? I think we should consider
  # not allowing them to be made with different index sizes
  #left_tensor  = δ(left_ind, sum_ind)
  left_tensor  = diagITensor(1.0,left_ind, sum_ind)
  right_tensor = ITensor(right_ind, sum_ind)
  for i in 1:dim(right_ind)
    right_tensor[right_ind(i), sum_ind(dim(left_ind) + i)] = 1
  end
  return left_tensor, right_tensor
end

"""
    add(A::MPS, B::MPS; kwargs...)
    +(A::MPS, B::MPS; kwargs...)

    add(A::MPO, B::MPO; kwargs...)
    +(A::MPO, B::MPO; kwargs...)

Add two MPS/MPO with each other, with some optional
truncation.
"""
function Base.:+(A::T, B::T; kwargs...) where {T <: AbstractMPS}
    N = length(A)
    length(B) != N && throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
    orthogonalize!(A, 1; kwargs...)
    orthogonalize!(B, 1; kwargs...)
    C = similar(A)
    rand_plev = 13124
    lAs = [linkind(A, i) for i in 1:N-1]
    prime!(A, rand_plev, "Link")

    first  = Vector{ITensor{2}}(undef,N-1)
    second = Vector{ITensor{2}}(undef,N-1)
    for i in 1:N-1
        lA = linkind(A, i)
        lB = linkind(B, i)
        r  = Index(dim(lA) + dim(lB), tags(lA))
        f, s = plussers(typeof(data(A[1])), lA, lB, r)
        first[i]  = f
        second[i] = s
    end
    C[1] = A[1] * first[1] + B[1] * second[1]
    for i in 2:N-1
        C[i] = dag(first[i-1]) * A[i] * first[i] + dag(second[i-1]) * B[i] * second[i]
    end
    C[N] = dag(first[N-1]) * A[N] + dag(second[N-1]) * B[N]
    prime!(C, -rand_plev, "Link")
    truncate!(C; kwargs...)
    return C
end

add(A::T, B::T;
    kwargs...) where {T <: AbstractMPS} = +(A, B; kwargs...)

"""
    sum(A::Vector{MPS}; kwargs...)

    sum(A::Vector{MPO}; kwargs...)

Add multiple MPS/MPO with each other, with some optional
truncation.
"""
function Base.sum(A::Vector{T};
                  kwargs...) where {T <: AbstractMPS}
  length(A) == 0 && return T()
  length(A) == 1 && return A[1]
  length(A) == 2 && return +(A[1], A[2]; kwargs...)
  nsize = isodd(length(A)) ? (div(length(A) - 1, 2) + 1) : div(length(A), 2)
  newterms = Vector{T}(undef, nsize)
  np = 1
  for n in 1:2:length(A) - 1
    newterms[np] = +(A[n], A[n+1]; kwargs...)
    np += 1
  end
  if isodd(length(A))
    newterms[nsize] = A[end]
  end
  return sum(newterms; kwargs...)
end

function orthogonalize!(M::AbstractMPS,
                        j::Int;
                        kwargs...)
  while leftlim(M) < (j-1)
    (leftlim(M) < 0) && setleftlim!(M, 0)
    b = leftlim(M)+1
    linds = uniqueinds(M[b],M[b+1])
    L,R = factorize(M[b], linds)
    M[b] = L
    M[b+1] *= R

    setleftlim!(M,b)
    if rightlim(M) < leftlim(M)+2
      setrightlim!(M, leftlim(M)+2)
    end
  end

  N = length(M)

  while rightlim(M) > (j+1)
    (rightlim(M) > (N+1)) && setrightlim!(M,N+1)
    b = rightlim(M)-2
    rinds = uniqueinds(M[b+1],M[b])
    L,R = factorize(M[b+1], rinds)
    M[b+1] = L
    M[b] *= R

    setrightlim!(M,b+1)
    if leftlim(M) > rightlim(M)-2
      setleftlim!(M, rightlim(M)-2)
    end
  end
end

function NDTensors.truncate!(M::AbstractMPS;
                             kwargs...)
  N = length(M)

  # Left-orthogonalize all tensors to make
  # truncations controlled
  orthogonalize!(M,N)

  # Perform truncations in a right-to-left sweep
  for j in reverse(2:N)
    rinds = uniqueinds(M[j],M[j-1])
    U,S,V = svd(M[j],rinds;kwargs...)
    M[j] = U
    M[j-1] *= (S*V)
    setrightlim!(M,j)
  end

end

NDTensors.contract(A::AbstractMPS,
                   B::AbstractMPS;
                   kwargs...) = *(A, B; kwargs...)

"""
    *(x::Number, M::MPS)

    *(x::Number, M::MPO)

Scale the MPS or MPO by the provided number.

Note: right now this just naively scales the
middle tensor. In the future, the plan would be
to scale the tensors between the left limit
and right limit by `x^(1/N)` where `N` is the distance
between the left limit and right limit.
"""
function Base.:*(x::Number, M::AbstractMPS)
  N = deepcopy(M)
  c = div(length(N), 2)
  N[c] .*= x
  return N
end

Base.:-(M::AbstractMPS) = Base.:*(-1,M)

hasqns(M::AbstractMPS) = hasqns(M[1])

function flux(M::AbstractMPS)::QN
  hasqns(M) || error("MPS or MPO does not conserve QNs")
  q = QN()
  for j=M.llim+1:M.rlim-1
    q += flux(M[j])
  end
  return q
end

totalqn(M::AbstractMPS) = flux(M)

@doc """
    orthogonalize!(M::MPS, j::Int; kwargs...)

    orthogonalize!(M::MPO, j::Int; kwargs...)

Move the orthogonality center of the MPS
to site `j`. No observable property of the
MPS will be changed, and no truncation of the
bond indices is performed. Afterward, tensors
`1,2,...,j-1` will be left-orthogonal and tensors
`j+1,j+2,...,N` will be right-orthogonal.
""" orthogonalize!

@doc """
    truncate!(M::MPS; kwargs...)

    truncate!(M::MPO; kwargs...)

Perform a truncation of all bonds of an MPS/MPO,
using the truncation parameters (cutoff,maxdim, etc.)
provided as keyword arguments.
""" truncate!

@deprecate orthoCenter(args...;
                       kwargs...) orthocenter(args...; kwargs...)

import .NDTensors.store

@deprecate store(m::AbstractMPS) data(m)

@deprecate replacesites!(args...;
                         kwargs...) ITensors.replacesiteinds!(args...; kwargs...)

@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...)

@deprecate applympo(args...; kwargs...) contract(args...; kwargs...)

@deprecate errorMPOprod(args...;
                        kwargs...) error_contract(args...;
                                                  kwargs...)

@deprecate error_mpoprod(args...;
                         kwargs...) error_contract(args...;
                                                   kwargs...)

@deprecate error_mul(args...;
                     kwargs...) error_contract(args...;
                                               kwargs...)

@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...)

import Base.sum

@deprecate sum(A::AbstractMPS,
               B::AbstractMPS; kwargs...) add(A, B; kwargs...)

@deprecate multmpo(args...;
                   kwargs...) contract(args...; kwargs...)

@deprecate set_leftlim!(args...;
                        kwargs...) ITensors.setleftlim!(args...;
                                                        kwargs...)

@deprecate set_rightlim!(args...;
                         kwargs...) ITensors.setrightlim!(args...;
                                                          kwargs...)

@deprecate tensors(args...;
                   kwargs...) ITensors.data(args...; kwargs...)

@deprecate primelinks!(args...;
                       kwargs...) ITensors.primelinkinds!(args...;
                                                          kwargs...)

@deprecate simlinks!(args...;
                     kwargs...) ITensors.simlinkinds!(args...;
                                                      kwargs...)

@deprecate mul(A::AbstractMPS,
               B::AbstractMPS;
               kwargs...) contract(A, B; kwargs...)


text:
module CalculatorController
  using JSON
  using TypedPolynomials
  using StaticArrays
  using LinearAlgebra
  using Optim
  using ReverseDiff
  using DiffResults
  using BenchmarkTools

  #if !(@isdefined coeffs)
      # Configuration constants
      const PATH = "../../electrode/clients/lib/efield/comsol/data/fit_coeffs.json"
      const KEY_ORDER = ["LP", "UP", "LW", "LE", "UE", "UW"]
      const PARAM_ORDER = ["bias", "dipole", "angle", "dEdx", "dEdy", "nux", "nuy", "d2Edx2", "d2Edy2"]

      # Optimization constants
      const weights = SA_F64[1, 0, 1, 1, 1, 0, 0, 1, 1]
      const Vmax = SA_F64[5000, 5000, 5000, 5000, 5000, 5000]
      const Vmin = SA_F64[-5000, -5000, -5000, -5000, -5000, -5000]

      # Physical constants
      const amu = 1.66054e-27 # kg
      const DVcmToJ = 3.33564e-28 # J

      # Physical parameters for KRb
      const D_POLY_FIT = SA_F64[0, 0.0561, -0.00337, 1.03812e-4, -1.42283e-6, 5.34778e-9]
      const m = 127*amu

      # Polynomial orders
      const xorder = 7
      const yorder = 7

      # Location to evaluate polynomials
      const xval = 0.0
      const yval = 0.0

      # Load JSON with polynomial coefficients
      const coeffs = JSON.parsefile(PATH)

      # Construct potential as polynomial
      @polyvar x y
      const monomials = [10^(i + j) * x^i * y^j for i in 0:xorder for j in 0:yorder]
      const Vpoly = [dot(convert(Vector{Float64}, coeffs[k]), monomials) for k in KEY_ORDER]

      # Compute polynomials corresponding to electric field and its derivatives
      const Epoly = SMatrix{6,2}([differentiate(-v, var) for v in Vpoly, var in [x, y]])
      const dEdxpoly = SMatrix{6,2}([differentiate.(e, x) for e in Epoly])
      const dEdypoly = SMatrix{6,2}([differentiate.(e, y) for e in Epoly])
      const d2Edx2poly = SMatrix{6,2}([differentiate.(de, x) for de in dEdxpoly])
      const d2Edy2poly = SMatrix{6,2}([differentiate.(de, y) for de in dEdypoly])

      # Evaluate polynomials corresponding to electric field and its derivatives
      const ExV = subs.(Epoly, (x, y) => (xval, yval))
      const dEdxV = subs.(dEdxpoly, (x, y) => (xval, yval))
      const dEdyV = subs.(dEdypoly, (x, y) => (xval, yval))
      const d2Edx2V = subs.(d2Edx2poly, (x, y) => (xval, yval))
      const d2Edy2V = subs.(d2Edy2poly, (x, y) => (xval, yval))

      # Construct dipole moment and its derivatives as polynomials
      @polyvar E
      const Dpoly = dot(D_POLY_FIT, [(E/1000)^i for i in 0:length(D_POLY_FIT)-1])
      const dDdEpoly = differentiate(Dpoly, E)
      const d2DdE2poly = differentiate(dDdEpoly, E)
 # end

  # Compute the parameters of the field given electrode potentials
  function get_params(V)
      # Evaluate the electric field and its derivatives given electrode potentials
      Ev = sum(ExV .* V, dims=1)
      dEdxv = sum(dEdxV .* V, dims=1)
      dEdyv = sum(dEdyV .* V, dims=1)
      d2Edx2v = sum(d2Edx2V .* V, dims=1)
      d2Edy2v = sum(d2Edy2V .* V, dims=1)

      # Compute the electric field magnitude, direction and gradient
      bias = norm(Ev)
      angle = atan(Ev[1], Ev[2]) * 180/π
      dEdx = dot(Ev, dEdxv)/bias
      dEdy = dot(Ev, dEdyv)/bias
      d2Edx2 = (dot(dEdxv, dEdxv) + dot(d2Edx2v, Ev) - dEdx^2)/bias
      d2Edy2 = (dot(dEdyv, dEdyv) + dot(d2Edy2v, Ev) - dEdy^2)/bias

      # Evaluate the dipole moment at the bias field and its derivatives
      dipole = subs(Dpoly, E => bias)
      dDdEv = subs(dDdEpoly, E => bias)
      d2DdE2v = subs(d2DdE2poly, E => bias)

      # Compute the trap frequency
      biaspoly = (-dipole/bias*0 + 2*dDdEv + bias*d2DdE2v)
      kx = (dipole/bias + dDdEv) * dot(Ev, d2Edx2v) + dEdx^2 * biaspoly
      ky = (dipole/bias + dDdEv) * dot(Ev, d2Edy2v) + dEdy^2 * biaspoly
      # TODO: Correct these formulae after testing program
      νx = -sign(kx)*sqrt(abs(100*DVcmToJ*kx) / (2 * π * m))
      νy = -sign(ky)*sqrt(abs(100*DVcmToJ*ky) / (2 * π * m))
      out = SVector(bias, dipole, angle, dEdx, dEdy, νx, νy, d2Edx2, d2Edy2)

      return out
  end

  function barrier(x, xmin, xmax)
      return 1E0 .* (-log.(max.(x .- xmin, nextfloat(0.0))) -log.(max.(xmax .- x, nextfloat(0.0))))
  end
  # Compute the least squares error
  function err(p, V, w)
      return dot(w, (get_params(V) .- p).^2) + sum(barrier(V, Vmin, Vmax))
  end

  # Optimize electrode potentials for parameters p
  function opt(p; V0=nothing, w=weights)
      if V0 === nothing
          V0 = randn(6)
      end
      f(V) = err(p, V, w)
      tape = ReverseDiff.GradientTape(f, V0)
      compiled_tape = ReverseDiff.compile(tape)
      result = DiffResults.GradientResult(V0)
      function fg!(F, G, V)
          ReverseDiff.gradient!(result, compiled_tape, V)
          G .= DiffResults.gradient(result)
          return DiffResults.value(result)
      end
      opts = Optim.Options(x_tol=1E-8, f_tol=1E-8, g_tol=1E-8)
      optimize(Optim.only_fg!(fg!), V0, BFGS(), opts)
      # TODO: Why doesn't this work?
      #optimize(Optim.only_fg!(fg!), Vmin, Vmax, V0, Fminbox(BFGS()), opts)
  end

  function opt_json(p; V0=nothing, w=weights)
    optres = opt(p; V0=V0, w=w)
    iter = optres.iterations
    bestV = max.(min.(optres.minimizer, Vmax), Vmin)
    p2 = get_params(bestV)
    Vdict = Dict(k => bestV[i] for (i, k) in enumerate(KEY_ORDER))
    pdict = Dict(k => p2[i] for (i, k) in enumerate(PARAM_ORDER))
    edict = Dict("err" => err(p, bestV, w), "iter" => iter, "barrier" => sum(barrier(bestV, Vmin, Vmax)))
    return merge(Vdict, pdict, edict)
  end

  function params_json(V)
    p = get_params(V)
    return Dict(k => p[i] for (i, k) in enumerate(PARAM_ORDER))
  end

  # Benchmark the optimization
  # goal = [12446.91, 0, 45, 250, 63.83, 0, 0]
  # println("\nBenchmark:")
  # @benchmark opt(x) setup=(x = randn(length(goal)))
end