text:
# ---------------------------------------------------------------------------- #
#
#   added_functionality.jl
#
#   Additional functions that make using egads from julia more straightforward.
#
#   egads.jl
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    EG_open( )

Opens an `egads` context.
"""
function EG_open( )
    ptr_obj = Ref{ego}()
    status = ccall((:EG_open, libegads_include), Cint, (Ptr{ego},), ptr_obj)
    return (ptr_obj[], status)
end

"""
    EG_getContext(object::ego)

Returns the `egads` context.
"""

function EG_getContext(object::ego)
    context_ptr = Ref{ego}()
    status = ccall((:EG_getContext, libegads_include), Cint, (ego, Ptr{ego}), object, context_ptr)
    return (context_ptr[], status)
end

"""
    EG_solidBoolean(src::ego, tool::ego, oper::Cint )

Performs a Boolean operation on the `ego` objects `src` and `tool`. `oper`
specifies the type of Boolean operation, such as `FUSION`, `INTERSECTION`, etc.
"""
function EG_solidBoolean(src::ego, tool::ego, oper::Cint )
    model_ptr = Ref{ego}()
    status = ccall((:EG_solidBoolean, libegads_include), Cint, (ego, ego,
                Cint, Ptr{ego}), src, tool, oper, model_ptr)
    return (model_ptr[], status)
end

"""
    EG_makeSolidBody(context::ego, stype::Cint, rvec)

Makes a solid object. `context` is the current egads context, `stype` specifies
the type of object, e.g. CYLINDER, SPHERE, etc. `rvec` is an array which holds
the data defining the object.
"""
function EG_makeSolidBody(context::ego, stype::Cint, rvec)
    body_ptr = Ref{ego}()
    status = ccall((:EG_makeSolidBody, libegads_include), Cint, (ego, Cint,
                Ptr{Cdouble}, Ptr{ego}), context, stype, rvec, body_ptr)
    return (body_ptr[], status)
end

"""
    EG_filletBody(src::ego, nedge::Cint, edges, radius::Cdouble, facemap)

Takes in a body `src` and puts a fillet of radius `radius` around the edges
defined in `edges`.
"""
function EG_filletBody(src::ego, nedge::Cint, edges, radius::Cdouble, facemap)
    result_ptr = Ref{ego}()
    status = ccall((:EG_filletBody, libegads_include), Cint, (ego, Cint, Ptr{ego}, Cdouble, Ptr{ego}, Ptr{Ptr{Cint}}), src, nedge, edges, radius, result_ptr, facemap)
    return (result_ptr[],status)
end

"""
    EG_makeTopology(context::ego, geom::ego, oclass::Cint, mtype::Cint, limits, nChildren::Cint, children, senses)

Makes a topology.
"""
function EG_makeTopology(context::ego, geom::ego, oclass::Cint, mtype::Cint, limits, nChildren::Cint, children, senses)
    topo_ptr = Ref{ego}()
    status = ccall((:EG_makeTopology, libegads_include), Cint, (ego, ego, Cint, Cint, Ptr{Cdouble}, Cint, Ptr{ego}, Ptr{Cint},
                Ptr{ego}), context, geom, oclass, mtype, limits, nChildren, children, senses, topo_ptr)
    return (topo_ptr[], status)
end

"""
    EG_makeGeometry(context::ego, oclass::Cint, mtype::Cint, refGeom::ego, ivec, rvec)

Returns a geometry.
"""
function EG_makeGeometry(context::ego, oclass::Cint, mtype::Cint, refGeom::ego, ivec, rvec)
    geom_ptr = Ref{ego}()
    status = ccall((:EG_makeGeometry, libegads_include), Cint, (ego, Cint, Cint, ego, Ptr{Cint}, Ptr{Cdouble}, Ptr{ego}), context, oclass, mtype, refGeom, ivec, rvec, geom_ptr)
    return (geom_ptr[],status)
end

"""
    getBodyFromModel( topo::ego )

Returns a body (`ego` object) from a model (also an `ego` object). This is
necessary to perform Boolean operations.
"""
function getBodyFromModel( topo::ego )
    geom     = Ref{ego}()
    oclass   = Ref{Int32}(0)
    _type    = Ref{Int32}(0)
    limits   = fill(0.0, 4)
    nchild   = Ref{Int32}(0)
    sense    = Ref{Ptr{Int32}}(0)
    children = Ref{Ptr{ego}}()

    status  = ccall((:EG_getTopology, libegads_include), Cint, (ego, Ptr{ego},
                     Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint},
                     Ptr{Ptr{ego}}, Ptr{Ptr{Cint}}), topo, geom, oclass, _type,
                     limits, nchild, children, sense)

    carray = children[]
    child0 = unsafe_load(carray,1)

    if nchild[] != 1
        println("nchild ", nchild[])
        warn("You are trying to get one body from a multi-body model.")
    end

    return (child0, status)
end

"""
    EG_invEvaluate(geom::ego, xyz, results)

Returns parameters from inverse evaluation.
"""
function EG_invEvaluate(geom::ego, xyz, results)
    param_ptr = Ref{Cdouble}(0.0)
    status = ccall((:EG_invEvaluate, libegads_include), Cint, (ego, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), geom, xyz, param_ptr, results)
    return (param_ptr[], status)
end

"""
    EG_copyObject(object::ego, oform)

Copies the `ego` object `object`.
"""
function EG_copyObject(object::ego, oform)
    copy_ptr = Ref{ego}()
    status = ccall((:EG_copyObject, libegads_include), Cint, (ego, Ptr{Void}, Ptr{ego}), object, oform, copy_ptr)
    return (copy_ptr[],status)
end

"""
    EG_makeFace(object::ego, mtype::Cint, limits)

Makes a face and returns it.
"""
function EG_makeFace(object::ego, mtype::Cint, limits)
    face_ptr = Ref{ego}()
    status = ccall((:EG_makeFace, libegads_include), Cint, (ego, Cint, Ptr{Cdouble}, Ptr{ego}), object, mtype, limits, face_ptr)
    return (face_ptr[],status)
end

"""
    EG_otherCurve(surface::ego, curve::ego, tol::Cdouble)

"""
function EG_otherCurve(surface::ego, curve::ego, tol::Cdouble)
    newcurve_ptr = Ref{ego}()
    status = ccall((:EG_otherCurve, libegads_include), Cint, (ego, ego, Cdouble, Ptr{ego}), surface, curve, tol, newcurve_ptr)
    return (newcurve_ptr[],status)
end

"""
    EG_sewFaces(nobj::Cint, objs, toler::Cdouble, flag::Cint)

Sews the faces in `objs` together and returns a new model.
"""
function EG_sewFaces(nobj::Cint, objs, toler::Cdouble, flag::Cint)
    result_ptr = Ref{ego}()
    status = ccall((:EG_sewFaces, libegads_include), Cint, (Cint, Ptr{ego}, Cdouble, Cint, Ptr{ego}), nobj, objs, toler, flag, result_ptr)
    return (result_ptr[], status)
end

"""
    EG_join( ebody1::ego, ebody2::ego, toler::Cdouble)

Joins two bodies with common faces.
"""
function EG_join( ebody1::ego, ebody2::ego, toler::Cdouble)
    # find faces - 1
    nface1_ptr = Ref{Cint}(0)
    faces1_ptr = Ref{Ptr{ego}}()
    egads.EG_getBodyTopos(ebody1, NULL_E, FACE, nface1_ptr, faces1_ptr)
    nface1 = nface1_ptr[]

    # find faces - 2
    nface2_ptr = Ref{Cint}(0)
    faces2_ptr = Ref{Ptr{ego}}()
    egads.EG_getBodyTopos(ebody2, NULL_E, FACE, nface2_ptr, faces2_ptr)
    nface2 = nface2_ptr[]

    # match faces
    nmatch_ptr = Ref{Cint}(0)
    match_ptr  = Ref{Ptr{Cint}}()
    EG_matchBodyFaces(ebody1, ebody2, toler, nmatch_ptr, match_ptr)
    nmatch = nmatch_ptr[]
    matches = fill( 0, 2*nmatch )
    for ii in 1:nmatch
        matches[2*ii-1] = unsafe_load(match_ptr[],2*ii-1)
        matches[2*ii]   = unsafe_load(match_ptr[],2*ii)
    end

    # delete all matching faces
    efaces = Vector{ ego }( nface1 + nface2 )
    jj = 1
    for ii in 1:nface1
        efaces[jj] = unsafe_load(faces1_ptr[],ii); jj += 1
        for kk in 1:nmatch
            if (ii == matches[2*kk-1])
                jj -= 1
            end
        end
    end
    for ii in 1:nface2
        efaces[jj] = unsafe_load(faces2_ptr[],ii); jj += 1
        for kk in 1:nmatch
            if (ii == matches[2*kk])
                jj -= 1
            end
        end
    end
    jj -= 1

    # sew faces
    (emodel, status) = EG_sewFaces( Cint(jj), efaces, toler, Cint(0) )

    return (emodel, status)

end

"""
    EG_getBodyFaces( ebody::ego, toler::Cdouble)

Get faces from body.
"""
function EG_getBodyFaces( ebody::ego, toler::Cdouble)

    nface_ptr = Ref{Cint}(0)
    faces_ptr = Ref{Ptr{ego}}()
    status = egads.EG_getBodyTopos(ebody, NULL_E, FACE, nface_ptr, faces_ptr)
    nface  = nface_ptr[]

    return (nface, faces_ptr, status)

end

"""
    EG_matchBodyFaces( ebody1::ego, ebody2::ego, toler::Cdouble)

Matches the faces on both bodies.
"""
function EG_matchBodyFaces( ebody1::ego, ebody2::ego, toler::Cdouble)

    # match faces
    nmatch_ptr = Ref{Cint}(0)
    match_ptr  = Ref{Ptr{Cint}}()
    status = EG_matchBodyFaces(ebody1, ebody2, toler, nmatch_ptr, match_ptr)
    nmatch = nmatch_ptr[]
    matches = fill( 0, 2*nmatch )
    for ii in 1:nmatch
        matches[2*ii-1] = unsafe_load(match_ptr[],2*ii-1)
        matches[2*ii]   = unsafe_load(match_ptr[],2*ii)
    end

    return (nmatch, matches, status)

end

"""
    cleanup( status::Cint, context::ego )

Cleans up the context if something went wrong.
"""
function cleanup( status::Cint, context::ego )
    warn("No success, status: ", status )
    status = EG_close(context)
    @printf("EG_close -> status=%d\n", status);
end

"""
    cleanup( context::ego )

Cleans up the context at end of a file.
"""
function cleanup( context::ego )
    status = EG_close(context)
    @printf("EG_close -> status=%d\n", status);
end


text:
const MEMORY_ORDERING_EXPLANATION = """
specified as a symbol (e.g., `:sequentially_consistent`) or a `Val` of a symbol (e.g.,
`Val(:sequentially_consistent)`)
"""

"""
    LLVM.Interop.atomic_pointerref(pointer::LLVMPtr{T}, ordering) -> value::T

Load a `value` from `pointer` with the given memory `ordering` atomically.

`ordering` is a Julia atomic ordering $MEMORY_ORDERING_EXPLANATION.

See also: `getproperty`, `getfield`
"""
atomic_pointerref

"""
    LLVM.Interop.atomic_pointerset(pointer::LLVMPtr{T}, x::T, ordering) -> pointer

Store a value `x` in `pointer` with the given memory `ordering` atomically.

`ordering` is a Julia atomic ordering $MEMORY_ORDERING_EXPLANATION.

See also: `setproperty!`, `setfield!`
"""
atomic_pointerset

"""
    LLVM.Interop.atomic_pointermodify(
        pointer::LLVMPtr{T},
        op,
        x::T,
        ordering,
    ) -> (old => new)::Pair{T,T}

Replace an `old` value stored at `pointer` with a `new` value comped as `new = op(old, x)`
with the given memory `ordering` atomically.  Return a pair `old => new`.

`ordering` is a Julia atomic ordering $MEMORY_ORDERING_EXPLANATION.

See also: `modifyproperty!`, `modifyfield!`
"""
atomic_pointermodify

"""
    LLVM.Interop.atomic_pointerswap(pointer::LLVMPtr{T}, op, new::T, ordering) -> old::T

Replace an `old` value stored at `pointer` with a `new` value with the given memory
`ordering` atomically.  Return the `old` value.

`ordering` is a Julia atomic ordering $MEMORY_ORDERING_EXPLANATION.

See also: `modifyproperty!`, `modifyfield!`
"""
atomic_pointerswap

"""
    LLVM.Interop.atomic_pointerreplace(
        pointer::LLVMPtr{T},
        expected::T,
        desired::T,
        success_ordering,
        fail_ordering,
    ) -> (; old::T, success::Bool)

Try to replace the value of `pointer` from an `expected` value to a `desired` value
atomically with the ordering `success_ordering`.  The property `old` of the returned value
is the value stored in the `pointer`.  The property `success` of the returned value
indicates if the replacement was successful.  The ordering `fail_ordering` specifies the
ordering used for loading the `old` value.

`success_ordering` and `fail_ordering` are Julia atomic orderings
$MEMORY_ORDERING_EXPLANATION.

See also: `replaceproperty!`, `replacefield!`
"""
atomic_pointerreplace

const _llvm_from_julia_ordering = (
    not_atomic = LLVM.API.LLVMAtomicOrderingNotAtomic,
    unordered = LLVM.API.LLVMAtomicOrderingUnordered,
    monotonic = LLVM.API.LLVMAtomicOrderingMonotonic,
    acquire = LLVM.API.LLVMAtomicOrderingAcquire,
    release = LLVM.API.LLVMAtomicOrderingRelease,
    acquire_release = LLVM.API.LLVMAtomicOrderingAcquireRelease,
    sequentially_consistent = LLVM.API.LLVMAtomicOrderingSequentiallyConsistent,
)

_julia_ordering(p) =
    Union{map(x -> p(x) ? Val{x} : Union{}, keys(_llvm_from_julia_ordering))...}

const AllOrdering = _julia_ordering(_ -> true)
const AtomicOrdering = _julia_ordering(!=(:not_atomic))

const LLVMOrderingVal = Union{map(x -> Val{x}, values(_llvm_from_julia_ordering))...}

is_stronger_than_monotonic(order::Symbol) =
    !(order === :monotonic || order === :unordered || order === :not_atomic)

for (julia, llvm) in pairs(_llvm_from_julia_ordering)
    @eval llvm_from_julia_ordering(::Val{$(QuoteNode(julia))}) = Val{$llvm}()
end

"""
    @dynamic_order(order) do order
        ... use order ...
    end

It is expanded to an expression similar to:

    if order === :not_atomic
        let order = Val(:not_atomic)
            ... use order ...
        end
    elseif order === :unordered
        let order = Val(:unordered)
            ... use order ...
        end
    elseif ...
        ...
    else
        throw(ConcurrencyViolationError(...))
    end

This is used for helping the compiler to optimize expressions such as
`atomic_pointerref(ptr, :monotonic)` and also to avoid abstract run-time dispatch.
"""
macro dynamic_order(thunk, order)
    @assert Meta.isexpr(thunk, :->, 2) && Meta.isexpr(thunk.args[1], :tuple, 1)
    ordervar = esc(thunk.args[1].args[1])
    body = esc(thunk.args[2])
    expr = foldr(
        keys(_llvm_from_julia_ordering),
        init = :(throw(ConcurrencyViolationError("invalid atomic ordering: ", order))),
    ) do key, r
        quote
            if order === $(QuoteNode(key))
                let $ordervar = Val{$(QuoteNode(key))}()
                    $body
                end
            else
                $r
            end
        end
    end
    quote
        order = $(esc(order))
        $expr
    end
end

_valueof(::Val{x}) where {x} = x

@inline function atomic_pointerref(pointer, order::Symbol)
    @dynamic_order(order) do order
        atomic_pointerref(pointer, order)
    end
end

@inline function atomic_pointerset(pointer, x, order::Symbol)
    @dynamic_order(order) do order
        atomic_pointerset(pointer, x, order)
    end
end

@generated function atomic_pointerref(ptr::LLVMPtr{T,A}, order::AllOrdering) where {T,A}
    sizeof(T) == 0 && return T.instance
    llvm_order = _valueof(llvm_from_julia_ordering(order()))
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)

        T_ptr = convert(LLVMType, ptr; ctx)

        T_typed_ptr = LLVM.PointerType(eltyp, A)

        # create a function
        param_types = [T_ptr]
        llvm_f, _ = create_function(eltyp, param_types)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)
            ld = load!(builder, typed_ptr)
            ordering!(ld, llvm_order)

            if A != 0
                metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A; ctx)
            end
            alignment!(ld, sizeof(T))

            ret!(builder, ld)
        end

        call_function(llvm_f, T, Tuple{LLVMPtr{T,A}}, :ptr)
    end
end

@generated function atomic_pointerset(
    ptr::LLVMPtr{T,A},
    x::T,
    order::AllOrdering,
) where {T,A}
    if sizeof(T) == 0
        # Mimicking what `Core.Intrinsics.atomic_pointerset` generates.
        # See: https://github.com/JuliaLang/julia/blob/v1.7.2/src/cgutils.cpp#L1570-L1572
        is_stronger_than_monotonic(order) || return :ptr
        return quote
            Core.Intrinsics.fence($(QuoteNode(order)))
            ptr
        end
    end
    llvm_order = _valueof(llvm_from_julia_ordering(order()))
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)
        T_ptr = convert(LLVMType, ptr; ctx)
        T_typed_ptr = LLVM.PointerType(eltyp, A)

        # create a function
        param_types = [T_ptr, eltyp]
        llvm_f, _ = create_function(LLVM.VoidType(ctx), param_types)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)
            val = parameters(llvm_f)[2]
            st = store!(builder, val, typed_ptr)
            ordering!(st, llvm_order)

            if A != 0
                metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A; ctx)
            end
            alignment!(st, sizeof(T))

            ret!(builder)
        end

        call = call_function(llvm_f, Cvoid, Tuple{LLVMPtr{T,A},T}, :ptr, :x)
        quote
            $call
            ptr
        end
    end
end

right(_, r) = r

const binoptable = [
    (:xchg, right, LLVM.API.LLVMAtomicRMWBinOpXchg),
    (:add, +, LLVM.API.LLVMAtomicRMWBinOpAdd),
    (:sub, -, LLVM.API.LLVMAtomicRMWBinOpSub),
    (:and, &, LLVM.API.LLVMAtomicRMWBinOpAnd),
    (:or, |, LLVM.API.LLVMAtomicRMWBinOpOr),
    (:xor, xor, LLVM.API.LLVMAtomicRMWBinOpXor),
    (:max, max, LLVM.API.LLVMAtomicRMWBinOpMax),
    (:min, min, LLVM.API.LLVMAtomicRMWBinOpMin),
    (:umax, max, LLVM.API.LLVMAtomicRMWBinOpUMax),
    (:umin, min, LLVM.API.LLVMAtomicRMWBinOpUMin),
    (:fadd, +, LLVM.API.LLVMAtomicRMWBinOpFAdd),
    (:fsub, -, LLVM.API.LLVMAtomicRMWBinOpFSub),
]

const AtomicRMWBinOpVal = Union{(Val{binop} for (_, _, binop) in binoptable)...}

@generated function llvm_atomic_op(
    binop::AtomicRMWBinOpVal,
    ptr::LLVMPtr{T,A},
    val::T,
    order::LLVMOrderingVal,
) where {T,A}
    Context() do ctx
        T_val = convert(LLVMType, T; ctx)
        T_ptr = convert(LLVMType, ptr; ctx)

        T_typed_ptr = LLVM.PointerType(T_val, A)

        llvm_f, _ = create_function(T_val, [T_ptr, T_val])

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)

            single_threaded = false
            rv = atomic_rmw!(
                builder,
                _valueof(binop()),
                typed_ptr,
                parameters(llvm_f)[2],
                _valueof(order()),
                single_threaded,
            )

            ret!(builder, rv)
        end

        call_function(llvm_f, T, Tuple{LLVMPtr{T,A},T}, :ptr, :val)
    end
end

@inline function atomic_pointermodify(pointer, op::OP, x, order::Symbol) where {OP}
    @dynamic_order(order) do order
        atomic_pointermodify(pointer, op, x, order)
    end
end

@inline function atomic_pointermodify(
    ptr::LLVMPtr{T},
    op,
    x::T,
    ::Val{:not_atomic},
) where {T}
    old = atomic_pointerref(ptr, Val(:not_atomic))
    new = op(old, x)
    atomic_pointerset(ptr, new, Val(:not_atomic))
    return old => new
end

@inline function atomic_pointermodify(
    ptr::LLVMPtr{T},
    ::typeof(right),
    x::T,
    order::AtomicOrdering,
) where {T}
    old = llvm_atomic_op(
        Val(LLVM.API.LLVMAtomicRMWBinOpXchg),
        ptr,
        x,
        llvm_from_julia_ordering(order),
    )
    return old => x
end

const atomictypes = Any[
    Int8,
    Int16,
    Int32,
    Int64,
    Int128,
    UInt8,
    UInt16,
    UInt32,
    UInt64,
    UInt128,
    Float16,
    Float32,
    Float64,
]

for (opname, op, llvmop) in binoptable
    opname === :xchg && continue
    types = if opname in (:min, :max)
        filter(t -> t <: Signed, atomictypes)
    elseif opname in (:umin, :umax)
        filter(t -> t <: Unsigned, atomictypes)
    elseif opname in (:fadd, :fsub)
        filter(t -> t <: AbstractFloat, atomictypes)
    else
        filter(t -> t <: Integer, atomictypes)
    end
    for T in types
        @eval @inline function atomic_pointermodify(
            ptr::LLVMPtr{$T},
            ::$(typeof(op)),
            x::$T,
            order::AtomicOrdering,
        )
            old = llvm_atomic_op($(Val(llvmop)), ptr, x, llvm_from_julia_ordering(order))
            return old => $op(old, x)
        end
    end
end

@inline atomic_pointerswap(pointer, new) = first(atomic_pointermodify(pointer, right, new))
@inline atomic_pointerswap(pointer, new, order) =
    first(atomic_pointermodify(pointer, right, new, order))

@inline function atomic_pointermodify(
    ptr::LLVMPtr{T},
    op,
    x::T,
    order::AllOrdering,
) where {T}
    # Should `fail_order` be stronger?  Ref: https://github.com/JuliaLang/julia/issues/45256
    fail_order = Val(:monotonic)
    old = atomic_pointerref(ptr, fail_order)
    while true
        new = op(old, x)
        (old, success) = atomic_pointerreplace(ptr, old, new, order, fail_order)
        success && return old => new
    end
end

@generated function llvm_atomic_cas(
    ptr::LLVMPtr{T,A},
    cmp::T,
    val::T,
    success_order::LLVMOrderingVal,
    fail_order::LLVMOrderingVal,
) where {T,A}
    llvm_success = _valueof(success_order())
    llvm_fail = _valueof(fail_order())
    Context() do ctx
        T_val = convert(LLVMType, T; ctx)
        T_pointee = T_val
        if T_val isa LLVM.FloatingPointType
            T_pointee = LLVM.IntType(sizeof(T) * 8; ctx)
        end
        T_ptr = convert(LLVMType, ptr; ctx)
        T_success = convert(LLVMType, Ptr{Int8}; ctx)

        T_typed_ptr = LLVM.PointerType(T_pointee, A)
        T_ok_ptr = LLVM.PointerType(convert(LLVMType, Int8; ctx))

        llvm_f, _ = create_function(T_val, [T_ptr, T_val, T_val, T_success])

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)
            ok_ptr = inttoptr!(builder, parameters(llvm_f)[4], T_ok_ptr)

            cmp_int = parameters(llvm_f)[2]
            if T_val isa LLVM.FloatingPointType
                cmp_int = bitcast!(builder, cmp_int, T_pointee)
            end

            val_int = parameters(llvm_f)[3]
            if T_val isa LLVM.FloatingPointType
                val_int = bitcast!(builder, val_int, T_pointee)
            end

            single_threaded = false
            res = atomic_cmpxchg!(
                builder,
                typed_ptr,
                cmp_int,
                val_int,
                llvm_success,
                llvm_fail,
                single_threaded,
            )

            rv = extract_value!(builder, res, 0)
            ok = extract_value!(builder, res, 1)
            ok = zext!(builder, ok, LLVM.Int8Type(ctx))
            store!(builder, ok, ok_ptr)

            if T_val isa LLVM.FloatingPointType
                rv = bitcast!(builder, rv, T_val)
            end

            ret!(builder, rv)
        end

        expr = call_function(
            llvm_f,
            T,
            Tuple{LLVMPtr{T,A},T,T,Ptr{Int8}},
            :ptr,
            :cmp,
            :val,
            :success_ptr,
        )
        quote
            success = Ref{Int8}()
            old = GC.@preserve success begin
                success_ptr = Ptr{Int8}(pointer_from_objref(success))
                $expr
            end
            (; old, success = success[] != zero(Int8))
        end
    end
end

@inline function atomic_pointerreplace(
    pointer,
    expected,
    desired,
    success_order::Symbol,
    fail_order::Symbol,
)
    # This avoids abstract dispatch at run-time but probably too much codegen?
    #=
    @dynamic_order(success_order) do success_order
        @dynamic_order(fail_order) do fail_order
            atomic_pointerreplace(pointer, expected, desired, success_order, fail_order)
        end
    end
    =#

    # This avoids excessive codegen while hopefully imposes no cost when const-prop works:
    so = @dynamic_order(success_order) do success_order
        success_order
    end
    fo = @dynamic_order(fail_order) do fail_order
        fail_order
    end
    return atomic_pointerreplace(pointer, expected, desired, so, fo)
end

@inline function atomic_pointerreplace(
    ptr::LLVMPtr{T},
    expected::T,
    desired::T,
    ::Val{:not_atomic},
    ::Val{:not_atomic},
) where {T}
    old = atomic_pointerref(ptr, Val(:not_atomic))
    if old === expected
        atomic_pointerset(ptr, desired, Val(:not_atomic))
        success = true
    else
        success = false
    end
    return (; old, success)
end

@inline atomic_pointerreplace(
    ptr::LLVMPtr{T},
    expected::T,
    desired::T,
    success_order::_julia_ordering(∉((:not_atomic, :unordered))),
    fail_order::_julia_ordering(∉((:not_atomic, :unordered, :release, :acquire_release))),
) where {T} = llvm_atomic_cas(
    ptr,
    expected,
    desired,
    llvm_from_julia_ordering(success_order),
    llvm_from_julia_ordering(fail_order),
)

ext: jl


text:
using SolitonBEC
using Test

names = [
    "./test_evolution.jl",
]

detected_tests = filter(
    name->startswith(name, "test_") && endswith(name, ".jl"),
    readdir("."))

unused_tests = setdiff(detected_tests, names)
if length(unused_tests) != 0
    error("The following tests are not used:\n", join(unused_tests, "\n"))
end

unavailable_tests = setdiff(names, detected_tests)
if length(unavailable_tests) != 0
    error("The following tests could not be found:\n", join(unavailable_tests, "\n"))
end

for name=names
    if startswith(name, "test_") && endswith(name, ".jl")
        include(name)
    end
end

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC


text:
# Contiguous on-device arrays

export
    CuDeviceArray, CuDeviceVector, CuDeviceMatrix, CuBoundsError, ldg

## construction

"""
    CuDeviceArray(dims, ptr)
    CuDeviceArray{T}(dims, ptr)
    CuDeviceArray{T,A}(dims, ptr)
    CuDeviceArray{T,A,N}(dims, ptr)

Construct an `N`-dimensional dense CUDA device array with element type `T` wrapping a
pointer, where `N` is determined from the length of `dims` and `T` is determined from the
type of `ptr`. `dims` may be a single scalar, or a tuple of integers corresponding to the
lengths in each dimension). If the rank `N` is supplied explicitly as in `Array{T,N}(dims)`,
then it must match the length of `dims`. The same applies to the element type `T`, which
should match the type of the pointer `ptr`.
"""
CuDeviceArray

# NOTE: we can't support the typical `tuple or series of integer` style construction,
#       because we're currently requiring a trailing pointer argument.

struct CuDeviceArray{T,N,A} <: AbstractArray{T,N}
    shape::Dims{N}
    ptr::DevicePtr{T,A}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    CuDeviceArray{T,N,A}(shape::Dims{N}, ptr::DevicePtr{T,A}) where {T,A,N} = new(shape,ptr)
end

const CuDeviceVector = CuDeviceArray{T,1,A} where {T,A}
const CuDeviceMatrix = CuDeviceArray{T,2,A} where {T,A}

# outer constructors, non-parameterized
CuDeviceArray(dims::NTuple{N,<:Integer}, p::DevicePtr{T,A})                where {T,A,N} = CuDeviceArray{T,N,A}(dims, p)
CuDeviceArray(len::Integer,              p::DevicePtr{T,A})                where {T,A}   = CuDeviceVector{T,A}((len,), p)

# outer constructors, partially parameterized
CuDeviceArray{T}(dims::NTuple{N,<:Integer},   p::DevicePtr{T,A}) where {T,A,N} = CuDeviceArray{T,N,A}(dims, p)
CuDeviceArray{T}(len::Integer,                p::DevicePtr{T,A}) where {T,A}   = CuDeviceVector{T,A}((len,), p)
CuDeviceArray{T,N}(dims::NTuple{N,<:Integer}, p::DevicePtr{T,A}) where {T,A,N} = CuDeviceArray{T,N,A}(dims, p)
CuDeviceVector{T}(len::Integer,               p::DevicePtr{T,A}) where {T,A}   = CuDeviceVector{T,A}((len,), p)

# outer constructors, fully parameterized
CuDeviceArray{T,N,A}(dims::NTuple{N,<:Integer}, p::DevicePtr{T,A}) where {T,A,N} = CuDeviceArray{T,N,A}(Int.(dims), p)
CuDeviceVector{T,A}(len::Integer,               p::DevicePtr{T,A}) where {T,A}   = CuDeviceVector{T,A}((Int(len),), p)

## getters

Base.pointer(a::CuDeviceArray) = a.ptr

Base.size(g::CuDeviceArray) = g.shape
Base.length(g::CuDeviceArray) = prod(g.shape)

## conversions

Base.unsafe_convert(::Type{DevicePtr{T,A}}, a::CuDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

# from CuArray
function Base.convert(::Type{CuDeviceArray{T,N,AS.Global}}, a::CuArray{T,N}) where {T,N}
    ptr = Base.unsafe_convert(Ptr{T}, Base.cconvert(Ptr{T}, a))
    CuDeviceArray{T,N,AS.Global}(a.shape, DevicePtr{T,AS.Global}(ptr))
end
cudaconvert(a::CuArray{T,N}) where {T,N} = convert(CuDeviceArray{T,N,AS.Global}, a)

## indexing

# TODO: arrays as allocated by the CUDA APIs are 256-byte aligned. we should keep track of
#       this information, because it enables optimizations like Load Store Vectorization
#       (cfr. shared memory and its wider-than-datatype alignment)

@inline function Base.getindex(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = datatype_align(T)
    Base.unsafe_load(pointer(A), index, Val(align))::T
end

@inline function Base.setindex!(A::CuDeviceArray{T}, x, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = datatype_align(T)
    Base.unsafe_store!(pointer(A), x, index, Val(align))
end

"""
    ldg(A, i)

Index the array `A` with the linear index `i`, but loads the value through the read-only
texture cache for improved cache behavior. You should make sure the array `A`, or any
aliased instance, is not written to for the duration of the current kernel.

This function can only be used on devices with compute capability 3.5 or higher.

See also: [`Base.getindex`](@ref)
"""
@inline function ldg(A::CuDeviceArray{T}, index::Integer) where {T}
    # FIXME: this only works on sm_35+, but we can't verify that for now
    @boundscheck checkbounds(A, index)
    align = datatype_align(T)
    unsafe_cached_load(pointer(A), index, Val(align))::T
end

Base.IndexStyle(::Type{<:CuDeviceArray}) = Base.IndexLinear()

## other

Base.show(io::IO, a::CuDeviceVector) =
    print(io, "$(length(a))-element device array at $(pointer(a))")
Base.show(io::IO, a::CuDeviceArray) =
    print(io, "$(join(a.shape, '×')) device array at $(pointer(a))")

Base.show(io::IO, mime::MIME"text/plain", a::CuDeviceArray) = show(io, a)

@inline function Base.unsafe_view(A::CuDeviceVector{T}, I::Vararg{Base.ViewIndex,1}) where {T}
    ptr = pointer(A) + (I[1].start-1)*sizeof(T)
    len = I[1].stop - I[1].start + 1
    return CuDeviceArray(len, ptr)
end

@inline function Base.iterate(A::CuDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end


text:
module FastLapackInterface

import Base.strides

using LinearAlgebra
using LinearAlgebra: BlasInt, BlasFloat, checksquare, chkstride1
using LinearAlgebra.BLAS: @blasfunc
using LinearAlgebra.LAPACK: chklapackerror

@static if VERSION < v"1.7"
    using LinearAlgebra.LAPACK: liblapack
else
    const liblapack = "libblastrampoline"
end

include("LinSolveAlgo.jl")
export LinSolveWs, linsolve_core!, linsolve_core_no_lu!, lu!
include("QrAlgo.jl")
export QrWs, QrpWs, geqrf_core!, geqp3!, ormqr_core!
include("SchurAlgo.jl")
export DgeesWs, dgees!, DggesWs, dgges!

end #module
#import LinearAlgebra: USE_BLAS64, LAPACKException


text:
using HierarchicalMatrices
using Test, LinearAlgebra, Random, InteractiveUtils

for r in map(BLOCKRANK, subtypes(AbstractFloat))
    @test iseven(r)
end

Random.seed!(0)

@hierarchical UpperTriangularHierarchicalMatrix Matrix (UpperTriangular, LowRankMatrix) Diagonal

for T in (Float32, Float64)
    A = rand(T, 10, 5)
    x = rand(T, 40)

    y = zero(x)
    HierarchicalMatrices.mul!(y, A, x, 1, 1)
    @test norm(y[1:10] - A*x[1:5]) ≤ eps(T)*norm(A*x[1:5])

    fill!(y, 0.0)
    HierarchicalMatrices.mul!(y, A, x, 5, 5, 2, 2)
    @test norm(y[5:2:23] - A*x[5:2:14]) ≤ eps(T)*norm(A*x[5:2:14])

    fill!(y, 0.0)
    HierarchicalMatrices.mul!(y, transpose(A), x, 1, 5, 2, 1)
    @test norm(y[1:5] - A'x[5:2:23]) ≤ eps(T)*norm(A'x[5:2:23])

    fill!(y, 0.0)
    HierarchicalMatrices.mul!(y, transpose(A), x, 6, 3, 1, 3)
    @test norm(y[6:3:18] - A'*x[3:12]) ≤ eps(T)*norm(A'*x[3:12])

    A = map(big, A)
    x = map(big, x)

    y = zero(x)
    HierarchicalMatrices.mul!(y, A, x, 1, 1)
    @test y[1:10] == A*x[1:5]

    fill!(y, 0.0)
    HierarchicalMatrices.mul!(y, A, x, 5, 5, 2, 2)
    @test y[5:2:23] == A*x[5:2:14]

    fill!(y, 0.0)
    HierarchicalMatrices.mul!(y, transpose(A), x, 1, 5, 2, 1)
    @test y[1:5] == A'x[5:2:23]

    fill!(y, 0.0)
    HierarchicalMatrices.mul!(y, transpose(A), x, 6, 3, 1, 3)
    @test y[6:3:18] == A'*x[3:12]

    UTHM = UpperTriangularHierarchicalMatrix(T, 2, 2)

    @test typeof(UTHM.UpperTriangularblocks) == Matrix{UpperTriangular{T,LowRankMatrix{T}}}
end

include("../examples/Kernel.jl")


text:
using Wynn
using SymPy
using Test

@syms x
s = [1, x, x^2/2, x^3/6, x^4/24]

@testset "Exponential EpsilonTable" begin
    @test EpsilonTable(s).etable[0,-1] == 0
    @test EpsilonTable(s).etable[1,-1] == 0
    @test EpsilonTable(s).etable[2,-1] == 0
    @test EpsilonTable(s).etable[3,-1] == 0
    @test EpsilonTable(s).etable[4,-1] == 0
    @test EpsilonTable(s).etable[-1,2] == -1/(x-1)
    @test EpsilonTable(s).etable[0,0]  == 1
    @test EpsilonTable(s).etable[1,0]  == 1+x
    @test EpsilonTable(s).etable[4,0]  == 1 + x + x^2/2 + x^3/6 + x^4/24
    @test EpsilonTable(s).etable[-2,4] == 2 / (x^2 - 2*x + 2)
    @test EpsilonTable(s).etable[0,4]  == (x^2 + 6*x + 12) / (x^2 -6*x + 12)
    @test EpsilonTable(s).etable[1,2]  == simplify(-(x^2 + 4*x + 6)/(2*x - 6))
    @test EpsilonTable(s).etable[-3,4] == 0
    @test EpsilonTable(s).etable[-1,0] == 0
end


text:
module ParserCombinator
  export Success,
         Failure,
         Parser,
         @rule,
         @define_rule,
         or,
         cat,
         repeat,
         chainl,
         map,
         p_string,
         regex

  type Success
    value
    next
    success
    Success(value, next) = new(value, next, true)
  end

  type Failure
    next
    success
    Failure(next) = new(next, false)
  end

  type Parser
    parse
    Parser(parse) = new(parse)
  end

  macro rule(body, name)
    esc(quote
      function $(name)()
        Parser((input) -> ($body)().parse(input))
      end
    end)
  end

  function or(lhs :: Parser, rhs :: Parser)
    Parser() do input
      r = lhs.parse(input)
      if !r.success
        rhs.parse(input)
      else
        r
      end
    end
  end

  function cat(lhs :: Parser, rhs :: Parser)
    Parser() do input
      r1 = lhs.parse(input)
      if r1.success
        r2 = rhs.parse(r1.next)
        if r2.success
          Success((r1.value, r2.value), r2.next)
        else
          Failure(r2.next)
        end
      else
        Failure(r1.next)
      end
    end
  end

  function map(parser :: Parser, f :: Function)
    Parser() do input
      r = parser.parse(input)
      if r.success
        Success(f(r.value), r.next)
      else
        r
      end
    end
  end

  function repeat(parser :: Parser)
    Parser() do input
      rest = input
      values = []
      while true
        r = parser.parse(rest)
        if !r.success
          return Success(values, rest)
        end
        values = vcat(values, r.value)
        rest = r.next
      end
    end
  end

  function chainl(p :: Parser, q :: Parser)
    value = map(cat(p, repeat(cat(q, p))),
                  function(value)
                    x = value[1]
                    xs = value[2]
                    a = x
                    while length(xs) > 0
                      f = xs[1][1]
                      b = xs[1][2]
                      a = f(a, b)
                      xs = xs[2:end]
                    end
                    return a
                  end
                 )
    return value
  end

  function p_string(literal)
    Parser() do input
      range = search(input, literal)
      if range.start != 1
        Failure(input)
      else
        Success(literal, input[length(literal)+1:end])
      end
    end
  end

  function regex(regex :: Regex)
    Parser() do input
      r = match(regex, input)
      if r != nothing
        Success(r.match, input[length(r.match)+1:end])
      else
        Failure(input)
      end
    end
  end

  Base.:+(p1::Parser, p2::Parser) = cat(p1, p2)
  Base.:/(p1::Parser, p2::Parser) = or(p1, p2)
  Base.:>(p::Parser, f::Function) = map(p, f)
end

max_stars_repo_name: kmizu/mlcombinators


text:
include("default_config.jl")
include("custom_salt.jl")
include("min_length.jl")
include("custom_alphabet.jl")
include("default_config_hex.jl")
include("custom_params_hex.jl")


text:
const EMPTY_DATEFORMAT = Dates.DateFormat("")

# Constructors
# independent constructor
function Source(fullpath::Union{AbstractString,IO};

              delim=COMMA,
              quotechar=QUOTE,
              escapechar=ESCAPE,
              null::AbstractString=String(""),

              header::Union{Integer,UnitRange{Int},Vector}=1, # header can be a row number, range of rows, or actual string vector
              datarow::Int=-1, # by default, data starts immediately after header or start of file
              types::Union{Dict{Int,DataType},Dict{String,DataType},Vector{DataType}}=DataType[],
              nullable::Bool=true,
              dateformat::Union{AbstractString,Dates.DateFormat}=EMPTY_DATEFORMAT,

              footerskip::Int=0,
              rows_for_type_detect::Int=100,
              rows::Int=0,
              use_mmap::Bool=true)
    # make sure character args are UInt8
    isascii(delim) || throw(ArgumentError("non-ASCII characters not supported for delim argument: $delim"))
    isascii(quotechar) || throw(ArgumentError("non-ASCII characters not supported for quotechar argument: $quotechar"))
    isascii(escapechar) || throw(ArgumentError("non-ASCII characters not supported for escapechar argument: $escapechar"))
    # isascii(decimal) || throw(ArgumentError("non-ASCII characters not supported for decimal argument: $decimal"))
    # isascii(separator) || throw(ArgumentError("non-ASCII characters not supported for separator argument: $separator"))
    delim = delim % UInt8
    quotechar = quotechar % UInt8
    escapechar = escapechar % UInt8
    decimal = CSV.PERIOD; decimal % UInt8
    separator = CSV.COMMA; separator % UInt8
    dateformat = isa(dateformat,Dates.DateFormat) ? dateformat : Dates.DateFormat(dateformat)
    return CSV.Source(fullpath=fullpath,
                        options=CSV.Options(delim,quotechar,escapechar,separator,
                                    decimal,ascii(null),null!="",dateformat,dateformat==Dates.ISODateFormat,-1,0,1,DataType[]),
                        header=header, datarow=datarow, types=types, nullable=nullable, footerskip=footerskip,
                        rows_for_type_detect=rows_for_type_detect, rows=rows, use_mmap=use_mmap)
end

function Source(;fullpath::Union{AbstractString,IO}="",
                options::CSV.Options=CSV.Options(),

                header::Union{Integer,UnitRange{Int},Vector}=1, # header can be a row number, range of rows, or actual string vector
                datarow::Int=-1, # by default, data starts immediately after header or start of file
                types::Union{Dict{Int,DataType},Dict{String,DataType},Vector{DataType}}=DataType[],
                nullable::Bool=true,

                footerskip::Int=0,
                rows_for_type_detect::Int=100,
                rows::Int=0,
                use_mmap::Bool=true)
    # delim=CSV.COMMA;quotechar=CSV.QUOTE;escapechar=CSV.ESCAPE;separator=CSV.COMMA;decimal=CSV.PERIOD;null="";header=1;datarow=2;types=DataType[];formats=UTF8String[];skipblankrows=true;footerskip=0;rows_for_type_detect=250;countrows=true
    # argument checks
    isa(fullpath,AbstractString) && (isfile(fullpath) || throw(ArgumentError("\"$fullpath\" is not a valid file")))
    isa(header,Integer) && datarow != -1 && (datarow > header || throw(ArgumentError("data row ($datarow) must come after header row ($header)")))

    isa(fullpath,IOStream) && (fullpath = chop(replace(fullpath.name,"<file ","")))

    # open the file for property detection
    if isa(fullpath,IOBuffer)
        source = fullpath
        fullpath = "<IOBuffer>"
    elseif isdefined(fullpath, :gz_file) # hack to detect GZip.Stream
        source = fullpath
        fullpath = fullpath.name
    elseif isa(fullpath,IO)
        source = IOBuffer(readbytes(fullpath))
        fullpath = fullpath.name
    else
        source = IOBuffer(use_mmap ? Mmap.mmap(fullpath) : open(readbytes,fullpath))
    end
    options.datarow != -1 && (datarow = options.datarow)
    options.rows != 0 && (rows = options.rows)
    options.header != 1 && (header = options.header)
    options.types != DataType[] && (types = options.types)
    rows = rows == 0 ? CSV.countlines(source,options.quotechar,options.escapechar) : rows
    # rows == 0 && throw(ArgumentError("No rows of data detected in $fullpath"))
    seekstart(source)

    datarow = datarow == -1 ? (isa(header,Vector) ? 0 : last(header)) + 1 : datarow # by default, data starts on line after header
    rows = rows - datarow + 1 - footerskip # rows now equals the actual number of rows in the dataset

    # figure out # of columns and header, either an Integer, Range, or Vector{String}
    # also ensure that `f` is positioned at the start of data
    cols = 0
    if isa(header,Integer)
        # default header = 1
        if header <= 0
            CSV.skipto!(source,1,datarow,options.quotechar,options.escapechar)
            datapos = position(source)
            cols = length(CSV.readsplitline(source,options.delim,options.quotechar,options.escapechar))
            seek(source, datapos)
            columnnames = ["Column$i" for i = 1:cols]
        else
            CSV.skipto!(source,1,header,options.quotechar,options.escapechar)
            columnnames = [x for x in CSV.readsplitline(source,options.delim,options.quotechar,options.escapechar)]
            cols = length(columnnames)
            datarow != header+1 && CSV.skipto!(source,header+1,datarow,options.quotechar,options.escapechar)
            datapos = position(source)
        end
    elseif isa(header,Range)
        CSV.skipto!(source,1,first(header),options.quotechar,options.escapechar)
        columnnames = [x for x in readsplitline(source,options.delim,options.quotechar,options.escapechar)]
        cols = length(columnnames)
        for row = first(header):(last(header)-1)
            for (i,c) in enumerate([x for x in readsplitline(source,options.delim,options.quotechar,options.escapechar)])
                columnnames[i] *= "_" * c
            end
        end
        datarow != last(header)+1 && CSV.skipto!(source,last(header)+1,datarow,options.quotechar,options.escapechar)
        datapos = position(source)
    else
        CSV.skipto!(source,1,datarow,options.quotechar,options.escapechar)
        datapos = position(source)
        cols = length(readsplitline(source,options.delim,options.quotechar,options.escapechar))
        seek(source,datapos)
        if isempty(header)
            columnnames = ["Column$i" for i = 1:cols]
        else
            length(header) == cols || throw(ArgumentError("length of provided header doesn't match number of columns of data at row $datarow"))
            columnnames = header
        end
    end

    # Detect column types
    columntypes = Array(DataType,cols)
    if isa(types,Vector) && length(types) == cols
        columntypes = types
    elseif isa(types,Dict) || isempty(types)
        poss_types = Array(DataType,min(rows < 0 ? rows_for_type_detect : rows, rows_for_type_detect),cols)
        fill!(poss_types, NullField)
        lineschecked = 0
        while !eof(source) && lineschecked < min(rows < 0 ? rows_for_type_detect : rows, rows_for_type_detect)
            vals = CSV.readsplitline(source,options.delim,options.quotechar,options.escapechar)
            lineschecked += 1
            for i = 1:cols
                i > length(vals) && continue
                poss_types[lineschecked,i] = CSV.detecttype(vals[i],options.dateformat,options.null)
            end
        end
        # detect most common/general type of each column of types
        d = Set{DataType}()
        for i = 1:cols
            for n = 1:lineschecked
                t = poss_types[n,i]
                t == NullField && continue
                push!(d,t)
            end
            columntypes[i] = (isempty(d) || WeakRefString{UInt8} in d ) ? WeakRefString{UInt8} :
                                (Date     in d) ? Date :
                                (DateTime in d) ? DateTime :
                                (Float64  in d) ? Float64 :
                                (Int      in d) ? Int : WeakRefString{UInt8}
            empty!(d)
        end
    else
        throw(ArgumentError("$cols number of columns detected; `types` argument has $(length(types)) entries"))
    end
    if isa(types,Dict{Int,DataType})
        for (col,typ) in types
            columntypes[col] = typ <: AbstractString ? WeakRefString{UInt8} : typ
        end
    elseif isa(types,Dict{String,DataType})
        for (col,typ) in types
            c = findfirst(columnnames, col)
            columntypes[c] = typ <: AbstractString ? WeakRefString{UInt8} : typ
        end
    end
    (any(columntypes .== DateTime) && options.dateformat == EMPTY_DATEFORMAT) && (options.dateformat = Dates.ISODateTimeFormat)
    (any(columntypes .== Date) && options.dateformat == EMPTY_DATEFORMAT) && (options.dateformat = Dates.ISODateFormat)
    if nullable
        columntypes = DataType[T <: Nullable ? T : Nullable{T} for T in columntypes]
    else
        # WeakRefStrings won't be safe if they don't end up in a NullableArray
        columntypes = DataType[T <: WeakRefString ? String : T for T in columntypes]
    end
    seek(source,datapos)
    return Source(Data.Schema(columnnames,columntypes,rows),
                  options,source,datapos,String(fullpath))
end

# construct a new Source from a Sink that has been streamed to (i.e. DONE)
function Source(s::CSV.Sink)
    io = s.data
    if isa(io,IOStream)
        nm = String(chop(replace(io.name,"<file ","")))
        data = IOBuffer(Mmap.mmap(nm))
    else
        seek(io, s.datapos)
        data = IOBuffer(readbytes(io))
        nm = String("")
    end
    seek(data,s.datapos)
    for (i, T) in enumerate(s.schema.types)
        if T <: AbstractString
            s.schema.types[i] = WeakRefString{UInt8}
        end
    end
    return Source(s.schema,s.options,data,s.datapos,nm)
end

# Data.Source interface
Data.reset!(io::CSV.Source) = seek(io.data,io.datapos)
Data.isdone(io::CSV.Source, row, col) = eof(io.data) || (row > io.schema.rows && io.schema.rows > -1)
Data.streamtype{T<:CSV.Source}(::Type{T}, ::Type{Data.Field}) = true
Data.getfield{T}(source::CSV.Source, ::Type{T}, row, col) = get(CSV.parsefield(source.data, T, source.options, row, col))
Data.getfield{T}(source::CSV.Source, ::Type{Nullable{T}}, row, col) = CSV.parsefield(source.data, T, source.options, row, col)
Data.reference(source::CSV.Source{Base.AbstractIOBuffer{Array{UInt8,1}}}) = source.data.data

"""
`CSV.read(fullpath::Union{AbstractString,IO}, sink::Type{T}=DataFrame, args...; kwargs...)` => `typeof(sink)`
`CSV.read(fullpath::Union{AbstractString,IO}, sink::Data.Sink; kwargs...)` => `Data.Sink`

parses a delimited file into a Julia structure (a DataFrame by default, but any `Data.Sink` may be given).

Positional arguments:

* `fullpath`; can be a file name (string) or other `IO` instance
* `sink::Type{T}`; `DataFrame` by default, but may also be other `Data.Sink` types that support streaming via `Data.Field` interface; note that the method argument can be the *type* of `Data.Sink`, plus any required arguments the sink may need (`args...`).
                    or an already constructed `sink` may be passed (2nd method above)

Keyword Arguments:

* `delim::Union{Char,UInt8}`; a single character or ascii-compatible byte that indicates how fields in the file are delimited; default is `UInt8(',')`
* `quotechar::Union{Char,UInt8}`; the character that indicates a quoted field that may contain the `delim` or newlines; default is `UInt8('"')`
* `escapechar::Union{Char,UInt8}`; the character that escapes a `quotechar` in a quoted field; default is `UInt8('\\')`
* `null::String`; an ascii string that indicates how NULL values are represented in the dataset; default is the empty string, `""`
* `header`; column names can be provided manually as a complete Vector{String}, or as an Int/Range which indicates the row/rows that contain the column names
* `datarow::Int`; specifies the row on which the actual data starts in the file; by default, the data is expected on the next row after the header row(s); for a file without column names (header), specify `datarow=1`
* `types`; column types can be provided manually as a complete Vector{DataType}, or in a Dict to reference individual columns by name or number
* `nullable::Bool`; indicates whether values can be nullable or not; `true` by default. If set to `false` and missing values are encountered, a `NullException` will be thrown
* `dateformat::Union{AbstractString,Dates.DateFormat}`; how all dates/datetimes in the dataset are formatted
* `footerskip::Int`; indicates the number of rows to skip at the end of the file
* `rows_for_type_detect::Int=100`; indicates how many rows should be read to infer the types of columns
* `rows::Int`; indicates the total number of rows to read from the file; by default the file is pre-parsed to count the # of rows; `-1` can be passed to skip a full-file scan, but the `Data.Sink` must be setup account for a potentially unknown # of rows
* `use_mmap::Bool=true`; whether the underlying file will be mmapped or not while parsing
* `append::Bool=false`; if the `sink` argument provided is an existing table, `append=true` will append the source's data to the existing data instead of doing a full replace

Note by default, "string" or text columns will be parsed as the [`WeakRefString`](https://github.com/quinnj/WeakRefStrings.jl) type. This is a custom type that only stores a pointer to the actual byte data + the number of bytes.
To convert a `String` to a standard Julia string type, just call `string(::WeakRefString)`, this also works on an entire column.
Oftentimes, however, it can be convenient to work with `WeakRefStrings` depending on the ultimate use, such as transfering the data directly to another system and avoiding all the intermediate copying.

Example usage:
```
julia> dt = CSV.read("bids.csv")
7656334×9 DataFrames.DataFrame
│ Row     │ bid_id  │ bidder_id                               │ auction │ merchandise      │ device      │
├─────────┼─────────┼─────────────────────────────────────────┼─────────┼──────────────────┼─────────────┤
│ 1       │ 0       │ "8dac2b259fd1c6d1120e519fb1ac14fbqvax8" │ "ewmzr" │ "jewelry"        │ "phone0"    │
│ 2       │ 1       │ "668d393e858e8126275433046bbd35c6tywop" │ "aeqok" │ "furniture"      │ "phone1"    │
│ 3       │ 2       │ "aa5f360084278b35d746fa6af3a7a1a5ra3xe" │ "wa00e" │ "home goods"     │ "phone2"    │
...
```
"""
function read(fullpath::Union{AbstractString,IO}, sink=DataFrame, args...; append::Bool=false,
              delim=COMMA,
              quotechar=QUOTE,
              escapechar=ESCAPE,
              null::AbstractString=String(""),
              header::Union{Integer,UnitRange{Int},Vector}=1, # header can be a row number, range of rows, or actual string vector
              datarow::Int=-1, # by default, data starts immediately after header or start of file
              types::Union{Dict{Int,DataType},Dict{String,DataType},Vector{DataType}}=DataType[],
              nullable::Bool=true,
              dateformat::Union{AbstractString,Dates.DateFormat}=EMPTY_DATEFORMAT,

              footerskip::Int=0,
              rows_for_type_detect::Int=100,
              rows::Int=0,
              use_mmap::Bool=true)
    source = Source(fullpath;
                  delim=delim, quotechar=quotechar, escapechar=escapechar, null=null,
                  header=header, datarow=datarow, types=types, nullable=nullable, dateformat=dateformat,
                  footerskip=footerskip, rows_for_type_detect=rows_for_type_detect, rows=rows, use_mmap=use_mmap)
    sink = Data.stream!(source, sink, append, args...)
    Data.close!(sink)
    return sink
end

function read{T}(fullpath::Union{AbstractString,IO}, sink::T; append::Bool=false,
              delim=COMMA,
              quotechar=QUOTE,
              escapechar=ESCAPE,
              null::AbstractString=String(""),
              header::Union{Integer,UnitRange{Int},Vector}=1, # header can be a row number, range of rows, or actual string vector
              datarow::Int=-1, # by default, data starts immediately after header or start of file
              types::Union{Dict{Int,DataType},Dict{String,DataType},Vector{DataType}}=DataType[],
              nullable::Bool=true,
              dateformat::Union{AbstractString,Dates.DateFormat}=EMPTY_DATEFORMAT,

              footerskip::Int=0,
              rows_for_type_detect::Int=100,
              rows::Int=0,
              use_mmap::Bool=true)
    source = Source(fullpath;
                  delim=delim, quotechar=quotechar, escapechar=escapechar, null=null,
                  header=header, datarow=datarow, types=types, nullable=nullable, dateformat=dateformat,
                  footerskip=footerskip, rows_for_type_detect=rows_for_type_detect, rows=rows, use_mmap=use_mmap)
    sink = Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end

read(source::CSV.Source, sink=DataFrame, args...; append::Bool=false) = (sink = Data.stream!(source, sink, append, args...); Data.close!(sink); return sink)
read{T}(source::CSV.Source, sink::T; append::Bool=false) = (sink = Data.stream!(source, sink, append); Data.close!(sink); return sink)


text:
# N = 13195 # 29
 N = 600851475143 # 6857
ubound = Int(ceil(sqrt(N)))
if ubound%2 == 0
	ubound -= 1
end
lpf = 1

function isPrime(n::Int)
	for i in 2:ceil(sqrt(n))
		if n%i==0
			return false
		end
	end
	return true
end

while ubound > 1
	if N%ubound==0 && isPrime(ubound)
		lpf = ubound
		break
	end
	ubound -= 2
end
println(N, '\t', lpf)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
#
# Computes the endpoints for fundamental solution integral by stepping
# along the contour until it finds a point whose value is almost smaller
# than machine epsilon
#

function find_endpoints(a::Float64, b::Float64, x::Float64, y::Float64, derivs::Bool)
    ϵ = 1e-14
    temp1 = sqrt(b^2-a+ZIM)
    tm = sqrt(2.0(b-temp1))
    tp = sqrt(2.0(b+temp1))
    st1 = log(tp)
    st3 = log(tm)
    if imag(st3) ≥ M_PI_2 st3 -= im*M_PI end

    if b ≤ sqrt(a)
        c1 = real(st1)
        stdist = 0.5abs(c1)
        if stdist == 0 stdist = 1.0 end

        # We start with finding lm1
        t = c1 - stdist
        gammaforbidden!(fe_gam, st3, t)
        dlm1 = 1.0
        integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
        if abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
            while abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
                dlm1 *= jump_ratio
                t = c1 - dlm1*stdist
                gammaforbidden!(fe_gam, st3, t)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
        else
            while abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ && dlm1 > ϵ
                dlm1 /= jump_ratio
                t = c1 - dlm1*stdist
                gammaforbidden!(fe_gam, st3, t)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
            if dlm1 < ϵ dlm1 = 0.0 end
            dlm1 *= jump_ratio
        end
        lm1 = c1 - dlm1*stdist

        # lp1 comes next
        t = c1 + stdist
        gammaforbidden!(fe_gam, st3, t)
        dlm1 = 1.0
        integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
        if abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ
            while abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ && dlm1 > ϵ
                dlm1 /= jump_ratio
                t = c1 + dlm1*stdist
                gammaforbidden!(fe_gam, st3, t)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
            if dlm1 < ϵ dlm1 = 0.0 end
            dlm1 *= jump_ratio
        else
            while abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
                dlm1 *= jump_ratio
                t = c1 + dlm1*stdist
                gammaforbidden!(fe_gam, st3, t)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
        end
        lp1 = c1 + dlm1*stdist

        c2 = c1
        lp2 = 0.0
        lm2 = 0.0
    else
        c1,c2 = real(st1),real(st3)
        if c1 ≥ c2
            c1,c2 = c2,c1
        end
        tmin,fmin = locate_minimum(st1, st3, a, b, x, y, derivs)

        f2 = -real(st3) + W2
        g2 = -real(st1) - V4

        stdist = abs(c1-tmin)
        if stdist == 0 stdist = 1.0 end

        # We start with finding lm1
        t = c1 - stdist

        f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
        g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
        fe_gam[1] = complex(t,f*g)

        dlm1 = 1.0
        integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
        if abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
            while abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
                dlm1 *= jump_ratio
                t = c1 - dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
        else
            while abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ && dlm1 > ϵ
                dlm1 /= jump_ratio
                t = c1 - dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
            if dlm1 < ϵ dlm1 = 0.0 end
            dlm1 *= jump_ratio
        end
        lm1 = c1 - dlm1*stdist

        # lp1 comes next
        t = c1 + stdist

        f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
        g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
        fe_gam[1] = complex(t,f*g)

        dlm1 = 1.0
        integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
        if abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ
            while abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ && dlm1 > ϵ
                dlm1 /= jump_ratio
                t = c1 + dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
            if dlm1 < ϵ dlm1 = 0.0 end
            dlm1 *= jump_ratio
        else
            while abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
                dlm1 *= jump_ratio
                t = c1 + dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
        end
        lp1 = c1 + dlm1*stdist

        #
        # Next up: lm2
        #
        stdist = abs(c2-tmin)
        if stdist == 0 stdist = 1.0 end

        t = c2 - stdist

        f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
        g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
        fe_gam[1] = complex(t,f*g)

        dlm1 = 1.0
        integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
        if abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
            while abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
                dlm1 *= jump_ratio
                t = c2 - dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
        else
            while abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ && dlm1 > ϵ
                dlm1 /= jump_ratio
                t = c2 - dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
            if dlm1 < ϵ dlm1 = 0.0 end
            dlm1 *= jump_ratio
        end
        lm2 = c2 - dlm1*stdist

        # Last but not least: lp2
        t = c2 + stdist

        f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
        g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
        fe_gam[1] = complex(t,f*g)

        dlm1 = 1.0
        integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
        if abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ
            while abs(fe_u[1]) < ϵ && abs(fe_ux[1]) < ϵ && abs(fe_uy[1]) < ϵ && dlm1 > ϵ
                dlm1 /= jump_ratio
                t = c2 + dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
            if dlm1 < ϵ dlm1 = 0.0 end
            dlm1 *= jump_ratio
        else
            while abs(fe_u[1]) > ϵ || abs(fe_ux[1]) > ϵ || abs(fe_uy[1]) > ϵ
                dlm1 *= jump_ratio
                t = c2 + dlm1*stdist
                f = M_F_1 * atan( 2 * ( t + f2)) - M_F_2
                g = M_G_1 * atan( -4 * (t + g2)) - M_G_2
                fe_gam[1] = complex(t,f*g)
                integrand!(fe_u, fe_ux, fe_uy, fe_gam, a, b, x, y, derivs)
            end
        end
        lp2 = c2 + dlm1*stdist

        # final cleanup
        if lm2 < c1
            if lm1 < lm2
                lm1 = lm2
            end
            lm2 = c2
        end
        if lp1 > c2
            if lp1 < lp2
                lp2 = lp1
            end
            lp1 = c1
        end
    end
    return lm1,lp1,lm2,lp2,c1,c2,tm,tp
end


text:
export go_board_graph, go_game_graph

# Can't use import or using here without getting a warning about
# missing dependencies. Instead we have to qualify everything used
# from LightGraphs. Let's make that more convenient with an alias.
const LG = LightGraphs

"""
    go_board_graph(n::Integer)

Construct an undirected graph according to the numbering scheme of the
GoGameGraphs package. The graph is returned as a `LightGraphs.SimpleGraph`.
"""
function go_board_graph(n::Integer)
    adjacency_list = Board(n).edges
    return LG.SimpleGraph(sum(length.(adjacency_list)) ÷ 2, adjacency_list)
end

"""
    go_game_graph(board_graph::LightGraphs.SimpleGraph)

Return the directed graph corresponding to playing the game of go on
the undirected graph `board_graph`. Although go is traditionally
played on a square grid, the rules can be applied to any undirected
graph. All the possible board positions act as vertices in the game
graph and the moves as edges. The empty board corresponds to vertex 1.
The result is a sparse directed and strongly connected graph with lots
of cycles. It is returned as a `LightGraphs.SimpleDiGraph`.

The size of the game graph grows exponentially with the number of
vertices in the corresponding board graph. This function supports
board graphs of up to 12 vertices. The size limit can be overridden
but you have been warned.

    go_game_graph(board_graph::LightGraphs.SimpleGraph; allow_self_capture = false)

The keyword argument `allow_self_capture` can be set to `false` to
disallow moves that causes the player to lose his own stones. This
results in a game graph with fewer edges than when self captures are
allowed.

(Side note: all variations of go rules agree on the mechanics of
playing moves but can have different opinions about which moves are
valid. Whether to allow multi-stone self captures is one such
difference. However, all go rules agree that single-stone self
captures are invalid since that would not change the position. In
graph terms that means no self-edges, regardless of the setting of
`allow_self_capture`.)
"""
function go_game_graph(board_graph::LG.SimpleGraph; allow_self_capture = true,
                       override_size_limit = false,
                       boards = nothing)
    if LG.nv(board_graph) > 12 && !override_size_limit
        error("Board has more than 12 vertices, which will produce a very large game graph. If you are confident that your computer can handle that, add the keyword argument `override_size_limit=true`.")
    end
    board_edges = [LG.outneighbors(board_graph, v) for v in LG.vertices(board_graph)]
    game_graph = GameGraph(Board(board_edges),
                           allow_suicide = allow_self_capture)
    forward_adjacency_list = game_graph.edges
    nv = length(forward_adjacency_list)
    ne = sum(length.(forward_adjacency_list))
    backward_adjacency_list = [Int[] for v in 1:nv]
    for v = 1:nv
        for w in forward_adjacency_list[v]
            push!(backward_adjacency_list[w], v)
        end
    end

    if boards isa Vector{Vector{Int}}
        append!(boards, game_graph.board)
    end

    return LG.SimpleDiGraph(ne, forward_adjacency_list, backward_adjacency_list)
end

max_stars_repo_name: GunnarFarneback/GoGameGraphs

max_issues_repo_name: GunnarFarneback/GoGameGraphs

max_forks_repo_name: GunnarFarneback/GoGameGraphs


text:
##################
# All 1D RefCube #
##################
function create_face_quad_rule(quad_rule::QuadratureRule{0,shape,T}, ::Interpolation{1,shape}) where {T,shape<:RefCube}
    w = getweights(quad_rule)
    face_quad_rule = QuadratureRule{1,shape,T}[]

    # Face 1
    new_points = [Vec{1,T}((-one(T),))] # ξ = -1
    push!(face_quad_rule, QuadratureRule{1,shape,T}(w, new_points))
    # Face 2
    new_points = [Vec{1,T}((one(T),))] # ξ = 1
    push!(face_quad_rule, QuadratureRule{1,shape,T}(w, new_points))

    return face_quad_rule
end

function weighted_normal(::Tensor{2,1,T}, ::FaceValues{1,1,T,RefCube}, face::Int) where {T}
    face == 1 && return Vec{1,T}((-one(T),))
    face == 2 && return Vec{1,T}(( one(T),))
    throw(ArgumentError("unknown face number: $face"))
end

##################
# All 2D RefCube #
##################
function create_face_quad_rule(quad_rule::QuadratureRule{1,shape,T}, ::Interpolation{2,shape}) where {T,shape<:RefCube}
    w = getweights(quad_rule)
    p = getpoints(quad_rule)
    n_points = length(w)
    face_quad_rule = QuadratureRule{2,shape,T}[]

    # Face 1
    new_points = [Vec{2,T}((p[i][1], -one(T))) for i in 1:n_points] # ξ = t, η = -1
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))
    # Face 2
    new_points = [Vec{2,T}((one(T), p[i][1])) for i in 1:n_points] # ξ = 1, η = t
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))
    # Face 3
    new_points = [Vec{2,T}((p[i][1], one(T))) for i in 1:n_points] # ξ = t, η = 1
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))
    # Face 4
    new_points = [Vec{2,T}((-one(T), p[i][1])) for i in 1:n_points] # ξ = -1, η = t
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))

    return face_quad_rule
end

function weighted_normal(J::Tensor{2,2}, ::FaceValues{2,2,T,RefCube}, face::Int) where {T}
    @inbounds begin
        face == 1 && return Vec{2}(( J[2,1], -J[1,1]))
        face == 2 && return Vec{2}(( J[2,2], -J[1,2]))
        face == 3 && return Vec{2}((-J[2,1],  J[1,1]))
        face == 4 && return Vec{2}((-J[2,2],  J[1,2]))
    end
    throw(ArgumentError("unknown face number: $face"))
end

#########################
# All RefTetrahedron 2D #
#########################
function create_face_quad_rule(quad_rule::QuadratureRule{1,shape,T}, ::Interpolation{2,shape}) where {T,shape<:RefTetrahedron}
    w = getweights(quad_rule)
    p = getpoints(quad_rule)
    n_points = length(w)
    face_quad_rule = QuadratureRule{2,shape,T}[]

    # Face 1
    new_points = [Vec{2,T}((p[i][1], one(T)-p[i][1])) for i in 1:n_points] # ξ = t, η = 1-t
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))
    # Face 2
    new_points = [Vec{2,T}((zero(T), p[i][1])) for i in 1:n_points] # ξ = 0, η = t
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))
    # Face 3
    new_points = [Vec{2,T}((p[i][1], zero(T))) for i in 1:n_points] # ξ = t, η = 0
    push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))

    return face_quad_rule
end

function weighted_normal(J::Tensor{2,2}, ::FaceValues{2,2,T,RefTetrahedron}, face::Int) where {T}
    @inbounds begin
        face == 1 && return Vec{2}((-(J[2,1] - J[2,2]), J[1,1] - J[1,2]))
        face == 2 && return Vec{2}((-J[2,2], J[1,2]))
        face == 3 && return Vec{2}((J[2,1], -J[1,1]))
    end
    throw(ArgumentError("unknown face number: $face"))
end

##################
# All RefCube 3D #
##################
function create_face_quad_rule(quad_rule::QuadratureRule{2,shape,T}, ::Interpolation{3,shape}) where {T,shape<:RefCube}
    w = getweights(quad_rule)
    p = getpoints(quad_rule)
    n_points = length(w)
    face_quad_rule = QuadratureRule{3,shape,T}[]

    # Face 1
    new_points = [Vec{3,T}((p[i][1], p[i][2], -one(T))) for i in 1:n_points] # ξ = t, η = s, ζ = -1
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 2
    new_points = [Vec{3,T}((p[i][1], -one(T), p[i][2])) for i in 1:n_points] # ξ = t, η = -1, ζ = s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 3
    new_points = [Vec{3,T}((one(T), p[i][1], p[i][2])) for i in 1:n_points] # ξ = 1, η = t, ζ = s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 4
    new_points = [Vec{3,T}((p[i][1], one(T), p[i][2])) for i in 1:n_points] # ξ = t, η = 1, ζ = s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 5
    new_points = [Vec{3,T}((-one(T), p[i][1], p[i][2])) for i in 1:n_points] # ξ = -1, η = t, ζ = s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 6
    new_points = [Vec{3,T}((p[i][1], p[i][2], one(T))) for i in 1:n_points] # ξ = t, η = s, ζ = 1
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))

    return face_quad_rule
end

function weighted_normal(J::Tensor{2,3}, ::FaceValues{3,3,T,RefCube}, face::Int) where {T}
    @inbounds begin
        face == 1 && return J[:,2] × J[:,1]
        face == 2 && return J[:,1] × J[:,3]
        face == 3 && return J[:,2] × J[:,3]
        face == 4 && return J[:,3] × J[:,1]
        face == 5 && return J[:,3] × J[:,2]
        face == 6 && return J[:,1] × J[:,2]
    end
    throw(ArgumentError("unknown face number: $face"))
end

#########################
# All RefTetrahedron 3D #
#########################
function create_face_quad_rule(quad_rule::QuadratureRule{2,shape,T}, ::Interpolation{3,shape}) where {T,shape<:RefTetrahedron}
    w = getweights(quad_rule)
    p = getpoints(quad_rule)
    n_points = length(w)
    face_quad_rule = QuadratureRule{3,shape,T}[]

    # Face 1
    new_points = [Vec{3,T}((p[i][1], p[i][2], zero(T))) for i in 1:n_points] # ξ = t, η = s, ζ = 0
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 2
    new_points = [Vec{3,T}((p[i][1], zero(T), p[i][2])) for i in 1:n_points] # ξ = t, η = 0, ζ = s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 3
    new_points = [Vec{3,T}((p[i][1], p[i][2], one(T)-p[i][1]-p[i][2])) for i in 1:n_points] # ξ = t, η = s, ζ = 1-t-s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))
    # Face 4
    new_points = [Vec{3,T}((zero(T), p[i][1], p[i][2])) for i in 1:n_points] # ξ = 0, η = t, ζ = s
    push!(face_quad_rule, QuadratureRule{3,shape,T}(w, new_points))

    return face_quad_rule
end

function weighted_normal(J::Tensor{2,3}, ::FaceValues{3,3,T,RefTetrahedron}, face::Int) where {T}
    @inbounds begin
        face == 1 && return J[:,2] × J[:,1]
        face == 2 && return J[:,1] × J[:,3]
        face == 3 && return (J[:,1]-J[:,3]) × (J[:,2]-J[:,3])
        face == 4 && return J[:,3] × J[:,2]
    end
    throw(ArgumentError("unknown face number: $face"))
end


text:
# Entropy
# =======

function StatsBase.entropy(table::Probabilities{T,N,A}) where {T,N,A}
    H = zero(T)
    p = gettablearray(table)
    @inbounds for pᵢ in p
        if pᵢ > zero(T)
            H += pᵢ * log(pᵢ)
        end
    end
    -H # Default base: e
end

function StatsBase.entropy(table::Counts{T,N,A}) where {T,N,A}
    H = zero(T)
    total = gettotal(table)
    n = gettablearray(table)
    @inbounds for nᵢ in n
        if nᵢ > zero(T)
            H += nᵢ * log(nᵢ/total)
        end
    end
    -H/total # Default base: e
end

"""
It calculates the Shannon entropy (H) from a table of `Counts` or `Probabilities`.
Use last and optional positional argument to change the base of the log. The default base
is e, so the result is in nats. You can use 2.0 as base to get the result in bits.
"""
function StatsBase.entropy(table::Union{Counts{T,N,A},Probabilities{T,N,A}}, base::Real) where {T,N,A}
    entropy(table) / log(base)
end

# Marginal Entropy
# ----------------

function marginal_entropy(table::Probabilities{T,N,A}, margin::Int) where {T,N,A}
    H = zero(T)
    marginals = getmarginalsarray(table)
    @inbounds for pi in view(marginals, :, margin)
        if pi > zero(T)
            H += pi * log(pi)
        end
    end
    -H # Default base: e
end

function marginal_entropy(table::Counts{T,N,A}, margin::Int) where {T,N,A}
    H = zero(T)
    total = gettotal(table)
    marginals = getmarginalsarray(table)
    @inbounds for ni in view(marginals, :, margin)
        if ni > zero(T)
            H += ni * log(ni/total)
        end
    end
    -H/total # Default base: e
end

"""
It calculates marginal entropy (H) from a table of `Counts` or `Probabilities`. The second
positional argument is used to indicate the magin used to calculate the entropy, e.g. it
estimates the entropy H(X) if marginal is 1, H(Y) for 2, etc.
Use last and optional positional argument to change the base of the log. The default base
is e, so the result is in nats. You can use 2.0 as base to get the result in bits.
"""
function marginal_entropy(table::Union{Counts{T,N,A},Probabilities{T,N,A}},
                          margin::Int, base::Real) where {T,N,A}
    marginal_entropy(table, margin) / log(base)
end

# Kullback-Leibler
# ================

function kullback_leibler(probabilities::Probabilities{T,N,A},
                          background::Array{T,N}) where {T,N,A}
    p = getcontingencytable(probabilities)
    @assert size(background)==size(p) "probabilities and background must have the same size."
    KL = zero(T)
    @inbounds for i in 1:length(p)
        pi = p[i]
        if pi > zero(T)
            KL += pi * log(pi/background[i])
        end
    end
    KL # Default base: e
end

function kullback_leibler(probabilities::Probabilities{T,N,A},
         background::Union{Probabilities{T,N,A},ContingencyTable{T,N,A}}=BLOSUM62_Pi) where {T,N,A}
    kullback_leibler(probabilities, gettablearray(background))
end

"""
It calculates the Kullback-Leibler (KL) divergence from a table of `Probabilities`. The
second positional argument is a `Probabilities` or `ContingencyTable` with the background
distribution. It's optional, the default is the `BLOSUM62_Pi` table.
Use last and optional positional argument to change the base of the log. The default base
is e, so the result is in nats. You can use 2.0 as base to get the result in bits.
"""
kullback_leibler(p::Probabilities{T,1,A}, q, base::Real) where {T,A} = kullback_leibler(p, q)/log(base)
kullback_leibler(p::Probabilities{T,1,A}, base::Real) where {T,A} = kullback_leibler(p)/log(base)

# Mutual Information
# ==================

# It avoids ifelse() because log is expensive (https://github.com/JuliaLang/julia/issues/8869)
@inline function _mi(::Type{T}, pij, pi, pj) where T
    (pij > zero(T)) && (pi > zero(T)) ? T(pij * log(pij/(pi*pj))) : zero(T)
end

function mutual_information(table::Probabilities{T,2,A}) where {T,A}
    MI = zero(T)
    marginals = getmarginalsarray(table)
    p = gettablearray(table)
    N = size(marginals,1)
    @inbounds for j in 1:N
        pj = marginals[j,2]
        if pj > zero(T)
            @inbounds @simd for i in 1:N
                MI += _mi(T, p[i,j], marginals[i,1], pj)
            end
        end
    end
    MI # Default base: e
end

# It avoids ifelse() because log is expensive (https://github.com/JuliaLang/julia/issues/8869)
@inline function _mi(total::T, nij, ni, nj) where T
    (nij > zero(T)) && (ni > zero(T)) ? T(nij * log((total * nij)/(ni * nj))) : zero(T)
end

function mutual_information(table::Counts{T,2,A}) where {T,A}
    MI = zero(T)
    marginals = getmarginalsarray(table)
    n = gettablearray(table)
    total = gettotal(table)
    N = size(marginals,1)
    @inbounds for j in 1:N
        nj = marginals[j,2]
        if nj > zero(T)
            @inbounds @simd for i in 1:N
                MI += _mi(total, n[i,j], marginals[i,1], nj)
            end
        end
    end
    MI/total # Default base: e
end

"""
It calculates Mutual Information (MI) from a table of `Counts` or `Probabilities`.
Use last and optional positional argument to change the base of the log. The default base
is e, so the result is in nats. You can use 2.0 as base to get the result in bits.
Calculation of MI from `Counts` is faster than from `Probabilities`.
"""
function mutual_information(table::Union{Counts{T,N,A}, Probabilities{T,N,A}},
                            base::Real) where {T,N,A}
    mutual_information(table) / log(base)
end

function mutual_information(pxyz::Union{Counts{T,3,A}, Probabilities{T,3,A}}) where {T,A}
    pxy = delete_dimensions(pxyz, 3)
    return(
        marginal_entropy(pxyz, 1) +                 # H(X) +
        marginal_entropy(pxyz, 2) +                 # H(Y) +
        marginal_entropy(pxyz, 3) -                 # H(Z) -
        entropy(pxy) -                              # H(X,Y) -
        entropy(delete_dimensions!(pxy, pxyz, 2)) - # H(X,Z) -
        entropy(delete_dimensions!(pxy, pxyz, 2)) + # H(Y,Z) +
        entropy(pxyz)                               # H(X,Y,Z)
    )
end

# Normalized Mutual Information by Entropy
# ----------------------------------------

"""
It calculates a Normalized Mutual Information (nMI) by Entropy from a table of `Counts` or
`Probabilities`.

`nMI(X, Y) = MI(X, Y) / H(X, Y)`
"""
function normalized_mutual_information(table::Union{Counts{T,N,A},Probabilities{T,N,A}}) where {T,N,A}
    H = entropy(table)
    if H != zero(T)
        MI = mutual_information(table)
        return(T(MI/H))
    else
        return(zero(T))
    end
end


text:
using LinearAlgebra

@doc """
  `block_qx(W::ITensor,ul::reg_form)::Tuple{ITensor,ITensor,Index}`

Perform a block respecting QX decomposition of the operator valued matrix `W`.
The appropriate decomposition, QR, RQ, QL, LQ is selected based on the `reg_form` `ul`
and the `orth` keyword argument.
The new internal `Index` between Q and R/L is modified so that the tags are "Link,qx" instead
"Link,qr" etc. returned by the qr/rq/ql/lq routines.  Q and R/L are also gauge fixed so that
the corner element of R/L is 1.0 and Q⁺Q=d𝕀 where d is the dimensionality of the local
Hilbert space.

# Arguments
- `W` Operator valued matrix for decomposition.
- `ul` upper/lower regular form of `W`. We can auto detect here, but is more efficient if this is done by the higher level calling routines.

# Keywords
- `orth::orth_type = right` : choose `left` or `right` orthogonal (canonical) form
- `epsrr::Float64 = 1e-14` : cutoff for rank revealing QX which removes zero pivot rows and columns.
   All rows with max(abs(R[r,:]))<epsrr are considered zero and removed.  epsrr==0.0 indicates no rank reduction.

# Returns a Tuple containing
- `Q` with orthonormal columns or rows depending on orth=left/right, dimensions: (χ+1)x(χq+1)
- `R` or `L` depending on `ul` with dimensions: (χq+2)x(χ\'+2)
- `iq` the new internal link index between `Q` and `R`/`L` with dimensions χq+2 and tags="Link,qx"

# Example
```julia
julia>using ITensors
julia>using ITensorMPOCompression
julia>N=5; #5 sites
julia>NNN=2; #Include 2nd nearest neighbour interactions
julia>sites = siteinds("S=1/2",N);
#
#  Make a Hamiltonian directly, i.e. not using autoMPO
#
julia>H=make_transIsing_MPO(sites,NNN);
#
#  Use pprint to see the structure for site #2. I = unit operator and S = any
#  other operator
#
julia>pprint(H[2]) #H[1] is a row vector, so let's see what H[2] looks like
I 0 0 0 0
S 0 0 0 0
S 0 0 0 0
0 0 I 0 0
0 S 0 S I
#
#  Now do a block respecting QX decomposition. QL decomposition is chosen because
#  H[2] is in lower regular form and the default ortho direction is left.
#
julia>Q,L,iq=block_qx(H[2]); #Block respecting QL
#
#  The first column of Q is unchanged because it is outside the V-block.
#  Also one column was removed because rank revealing QX is the default algorithm.
#
julia>pprint(Q)
I 0 0 0
S 0 0 0
S 0 0 0
0 I 0 0
0 0 S I
#
#  Similarly L is missing one row due to the rank revealing QX algorithm
#
julia>pprint(L,iq) #we need to tell pprint which index is the row index.
I 0 0 0 0
0 0 I 0 0
0 S 0 S 0
0 0 0 0 I
```
"""
function block_qx(W::ITensor,ul::reg_form=lower;kwargs...)::Tuple{ITensor,ITensor,Index}
  d,n,r,c=parse_links(W)
  return block_qx(W,n,r,c,ul;kwargs...)
end

function block_qx(W_::ITensor,n::Int64,r::Index,c::Index,ul::reg_form=lower;kwargs...)::Tuple{ITensor,ITensor,Index}
  #
  # Copy so that we don't mess up the original MPO
  #
  W=copy(W_)
  #
  # settle the left/right && upper/lower question
  #
  lr::orth_type=get(kwargs, :orth, left)
  ms=matrix_state(ul,lr)
  #
  #  decide some strings and variables based on lr.
  #
  (tln,cr)= lr==left ? ("l=$n",c) : ("l=$(n-1)",r)

  ilw=filterinds(inds(W),tags=tln)[1] #get the link to the next site
  offset=V_offsets(ms)
  V=getV(W,offset) #extract the V block
  il=filterinds(inds(V),tags=tln)[1] #link to next site
  iothers=noncommoninds(inds(V),il) #group all other indices for QX factorization

  if ul==lower
    if lr==left
      Q,RL=ql(V,iothers;positive=true,kwargs...) #block respecting QL decomposition
    else #right
      RL,Q=lq(V,iothers;positive=true,kwargs...) #block respecting LQ decomposition
    end
  else #upper
    if lr==left
      Q,RL=ITensorMPOCompression.qr(V,iothers;positive=true,kwargs...) #block respecting QR decomposition
    else #right
      RL,Q=rq(V,iothers;positive=true,kwargs...) #block respecting RQ decomposition
    end
  end
  set_scale!(RL,Q,offset) #rescale so the L(n,n)==1.0
  ITensors.@debug_check begin
    @assert norm(V-RL*Q)<1e-12 #make sure decomp worked
  end
  qx=String(tags(commonind(RL,Q))[2]) #should be "ql","lq","qr" os "rq"
  replacetags!(Q ,qx,"qx") #releive client code from the burden dealing ql,lq,qr,rq tags
  replacetags!(RL,qx,"qx")
  W=setV(W,Q,ms) #Q is the new V, stuff Q into W. THis can resize W
  RLplus,iqx=growRL(RL,ilw,offset) #Now make a full size version of RL
  ilw=filterinds(W,tags=tln)[1]
  replaceind!(W,ilw,iqx)
  @assert hastags(W,"qx")
  return W,RLplus,iqx
end

function ql!(A::StridedMatrix{<:LAPACK.BlasFloat}, ::NoPivot; blocksize=36)
  tau=similar(A, Base.min(size(A)...))
  x=LAPACK.geqlf!(A, tau)
  #save L from the lower portion of A, before orgql! mangles it!
  nr,nc=size(A)
  mn=Base.min(nr,nc)
  L=similar(A,(mn,nc))
  for r in 1:mn
    for c in 1:r+nc-mn
      L[r,c]=A[r+nr-mn,c]
    end
    for c in r+1+nc-mn:nc
      L[r,c]=0.0
    end
  end
  # Now we need shift the orth vectors from the right side of Q over the left side, before
  if (mn<nc)
      for r in 1:nr
        for c in 1:mn
            A[r,c]=A[r,c+nc-mn];
        end
      end
  end
# A may now have extra columns from mn+1:nc, but they get chopped of in the outer functions
  LAPACK.orgql!(A,tau)
  return  A,L
end

function lq!(A::StridedMatrix{<:LAPACK.BlasFloat}, ::NoPivot; blocksize=36)
  tau=similar(A, Base.min(size(A)...))
  x=LAPACK.gelqf!(A, tau)
  #save L from the lower portion of A, before orgql! mangles it!
  nr,nc=size(A)
  mn=Base.min(nr,nc)
  L=similar(A,(nr,mn))
  for c in 1:mn
    for r in 1:c-1
      L[r,c]=0.0
    end
    for r in c:nr
      L[r,c]=A[r,c]
    end
  end

  LAPACK.orglq!(A,tau)
  # now A is Q but it may nave extra rows that need to be chopped off
  Q=A[1:mn,:]
  return  L,Q
end

function rq!(A::StridedMatrix{<:LAPACK.BlasFloat}, ::NoPivot; blocksize=36)
  tau=similar(A, Base.min(size(A)...))
  x=LAPACK.gerqf!(A, tau)
  #save R from the lower portion of A, before orgql! mangles it!
  nr,nc=size(A)
  mn=Base.min(nr,nc)
  R=similar(A,(nr,mn))
  for c in 1:mn
    for r in 1:c+nr-mn
      R[r,c]=A[r,c+nc-mn]
    end
    for r in c+1+nr-mn:nr
      R[r,c]=0.0
    end
  end
  #
  # If nr>nc we need shift the orth vectors from the bottom of Q up to top before
  # unpacking the reflectors.
  #
  if mn<nr
    for c in 1:nc
      for r in 1:mn
        A[r,c]=A[r+nr-mn,c]
      end
    end
    A=A[1:mn,:] #whack the extra rows in A or orgrq! will complain
  end
  LAPACK.orgrq!(A,tau)
  return  R,A
end

ql!(A::AbstractMatrix) = ql!(A, NoPivot())
lq!(A::AbstractMatrix) = lq!(A, NoPivot())
rq!(A::AbstractMatrix) = rq!(A, NoPivot())

function ql(A::AbstractMatrix{T}, arg...; kwargs...) where T
    Base.require_one_based_indexing(A)
    AA = similar(A, LinearAlgebra._qreltype(T), size(A))
    copyto!(AA, A)
    return ql!(AA, arg...; kwargs...)
end

function lq(A::AbstractMatrix{T}, arg...; kwargs...) where T
  Base.require_one_based_indexing(A)
  AA = similar(A, LinearAlgebra._qreltype(T), size(A))
  copyto!(AA, A)
  return lq!(AA, arg...; kwargs...)
end

function rq(A::AbstractMatrix{T}, arg...; kwargs...) where T
  Base.require_one_based_indexing(A)
  AA = similar(A, LinearAlgebra._qreltype(T), size(A))
  copyto!(AA, A)
  return rq!(AA, arg...; kwargs...)
end

function ql_positive(M::AbstractMatrix)
  sparseQ, L = ql(M)
  Q = convert(Matrix, sparseQ)
  nr,nc = size(Q)
  dc=nc>nr ? nc-nr : 0 #diag is shifted over by dc if nc>nr
  for c in 1:nc
    if c<=nr && real(L[c, c+dc]) < 0.0
      L[c, 1:c+dc] *= -1
      Q[:,c] *= -1
    end
  end
  return (Q, L)
end

function lq_positive(M::AbstractMatrix)
  L, sparseQ = lq(M)
  Q = convert(Matrix, sparseQ)
  nr,nc = size(Q)
  for r in 1:nr
    if r<=nc && real(L[r, r]) < 0.0
      L[r:end, r] *= -1
      Q[r,:] *= -1
    end
  end
  return (L, Q)
end

function rq_positive(M::AbstractMatrix)
  R, sparseQ = rq(M)
  Q = convert(Matrix, sparseQ)
  nr,nc = size(R)
  dr=nr>nc ? nr-nc : 0 #diag is shifted down by dr if nr>nc
  for r in 1:nr
    if r<=nc && real(R[r+dr, r]) < 0.0
      R[1:r+dr, r] *= -1
      Q[r,:] *= -1
    end
  end
  return (R, Q)
end

function ql(T::NDTensors.DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
    positive = get(kwargs, :positive, false)
    # TODO: just call qr on T directly (make sure
    # that is fast)
    if positive
      QM, LM = ql_positive(matrix(T))
    else
      QM, LM = ql(matrix(T))
    end
    # Make the new indices to go onto Q and L
    q, l = inds(T)
    q = dim(q) < dim(l) ? sim(q) : sim(l)
    Qinds = IndsT((ind(T, 1), q))
    Linds = IndsT((q, ind(T, 2)))
    Q = NDTensors.tensor(NDTensors.Dense(vec(Matrix(QM))), Qinds)
    L = NDTensors.tensor(NDTensors.Dense(vec(LM)), Linds)
    return Q, L
  end

  function lq(T::NDTensors.DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
    positive = get(kwargs, :positive, false)
    # TODO: just call qr on T directly (make sure
    # that is fast)
    if positive
      LM, QM = lq_positive(matrix(T))
    else
      LM, QM = lq(matrix(T))
    end

    # Make the new indices to go onto Q and L
    l, q = inds(T)
    #@show inds(T)
    q = dim(q) < dim(l) ? sim(q) : sim(l)
    Qinds = IndsT((q,ind(T, 2)))
    Linds = IndsT((ind(T, 1),q))
    #@show Linds Qinds
    Q = NDTensors.tensor(NDTensors.Dense(vec(Matrix(QM))), Qinds)
    L = NDTensors.tensor(NDTensors.Dense(vec(LM)), Linds)
    return L, Q
  end

  function rq(T::NDTensors.DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
    positive = get(kwargs, :positive, false)
    # TODO: just call qr on T directly (make sure
    # that is fast)
    if positive
      RM, QM = rq_positive(matrix(T))
    else
      RM, QM = rq(matrix(T))
    end

    # Make the new indices to go onto Q and L
    #@show inds(T)
    l, q = inds(T)
    q = dim(q) < dim(l) ? sim(q) : sim(l)
    Qinds = IndsT((q,ind(T, 2)))
    Linds = IndsT((ind(T, 1),q))
    #@show Linds Qinds
    Q = NDTensors.tensor(NDTensors.Dense(vec(Matrix(QM))), Qinds)
    R = NDTensors.tensor(NDTensors.Dense(vec(RM)), Linds)
    return R, Q
  end

# ql decomposition of an order-n dense tensor according to
# positions Lpos and Rpos
function ql(
  T::NDTensors.DenseTensor{<:Number,N,IndsT},
  Lpos::NTuple{NL,Int},
  Rpos::NTuple{NR,Int};
  kwargs...
) where {N,IndsT,NL,NR}
  M = NDTensors.permute_reshape(T, Lpos, Rpos)
  QM, LM = ql(M; kwargs...)
  l = ind(LM, 1)
  # TODO: simplify this by permuting inds(T) by (Lpos,Rpos)
  # then grab Linds,Rinds
  Linds = NDTensors.similartype(IndsT, Val{NL})(ntuple(i -> inds(T)[Lpos[i]], Val(NL)))
  Qinds = NDTensors.push(Linds, l)
  Q = reshape(QM, Qinds)
  Rinds = NDTensors.similartype(IndsT, Val{NR})(ntuple(i -> inds(T)[Rpos[i]], Val(NR)))
  Rinds = NDTensors.pushfirst(Rinds, l)
  L = NDTensors.reshape(LM, Rinds)
  return Q, L
end
# ql decomposition of an order-n block sparse tensor according to
# positions Lpos and Rpos
function ql(
  T::NDTensors.BlockSparseTensor{Elt,N,IndsT},
  Lpos::NTuple{NL,Int},
  Rpos::NTuple{NR,Int};
  kwargs...
) where {Elt,N,IndsT,NL,NR}

  nb=nnzblocks(T)
  Qs = Vector{DenseTensor{Elt,NL+1,NTuple{NL+1, Int64},<:NDTensors.Dense}}(undef, nb)
  Ls = Vector{DenseTensor{Elt,NR+1,NTuple{NR+1, Int64},<:NDTensors.Dense}}(undef, nb)

  #@show Lpos Rpos inds(T) dim(T,1) dim(T,2)
  #
  @show "in ql block sparse" Lpos Rpos
  for (n, b) in enumerate(eachnzblock(T))
    blockT = blockview(T, b) #this is a DenseTensor so we call the dense ql for each block.
    Qs[n],Ls[n] = ql(blockT,Lpos,Rpos;kwargs...) #All reshaping is taken care of in this call.
    @show blockT Qs[n] Ls[n]
  end

  #
  #  We need to figure out what is the index between Q and L this code is stolen from
  #  ITensors/NDTensors/src/blocksparse/linearalgebra.jl line ~107 in
  #  function LinearAlgebra.svd(T::BlockSparseMatrix{ElT}; kwargs...)
  #  I don;t know what it is supposed to be doing
  #
  nb1_lt_nb2 = (
    nblocks(T)[1] < nblocks(T)[2] ||
    (nblocks(T)[1] == nblocks(T)[2] && dim(T, 1) < dim(T, 2))
  )

  if nb1_lt_nb2
    qind = sim(ind(T, 1))
  else
    qind = sim(ind(T, 2))
  end
   # qind may have too many blocks ... really why?
  if nblocks(qind) > nb
    resize!(qind, nb)
  end

  if dir(qind) != dir(inds(T)[1])
    qind = dag(qind)
  end
  indsQ = setindex(inds(T), qind, 2)
  @show indsQ
  indsL = dag(qind),inds(T)[2]
  @show indsL typeof(indsL)
  nzblocksQ = Vector{Block{4}}(undef, nb)
  nzblocksL = Vector{Block{2}}(undef, nb)
  for n in 1:nb
    blockT=nzblocks(T)[n]
    @show blockT
    blockQ = (blockT[Lpos[1]],UInt(1), blockT[Lpos[2]],blockT[Lpos[3]] )
    @show blockQ
    nzblocksQ[n]=blockQ
    blockL = (blockT[Rpos[1]], UInt(1))
    @show blockL
    nzblocksL[n]=blockL
  end
  Q = BlockSparseTensor(Elt, undef, nzblocksQ, indsQ)
  L = BlockSparseTensor(Elt, undef, nzblocksL, indsL)
  for n in 1:nb
    blockview(Q, nzblocksQ[n]) .= Qs[n]
    blockview(L, nzblocksL[n]) .= Ls[n]
  end
  @show Q
  @show L
  return Q,L
end

# lq decomposition of an order-n tensor according to
# positions Lpos and Rpos
function lq(
  T::NDTensors.DenseTensor{<:Number,N,IndsT},
  Lpos::NTuple{NL,Int},
  Rpos::NTuple{NR,Int};
  kwargs...
) where {N,IndsT,NL,NR}
  M = NDTensors.permute_reshape(T, Lpos, Rpos)
  LM, QM = lq(M; kwargs...)
  l = ind(LM, 2)
  # TODO: simplify this by permuting inds(T) by (Lpos,Rpos)
  # then grab Linds,Rinds
  Rinds = NDTensors.similartype(IndsT, Val{NR})(ntuple(i -> inds(T)[Rpos[i]], Val(NR)))
  Qinds = NDTensors.pushfirst(Rinds, l)
  Q = reshape(QM, Qinds)
  Linds = NDTensors.similartype(IndsT, Val{NL})(ntuple(i -> inds(T)[Lpos[i]], Val(NL)))
  Linds = NDTensors.push(Linds, l)
  L = NDTensors.reshape(LM, Linds)
  return L, Q
end

# lq decomposition of an order-n tensor according to
# positions Lpos and Rpos
function rq(
  T::NDTensors.DenseTensor{<:Number,N,IndsT},
  Lpos::NTuple{NL,Int},
  Rpos::NTuple{NR,Int};
  kwargs...
) where {N,IndsT,NL,NR}
  M = NDTensors.permute_reshape(T, Lpos, Rpos)
  RM, QM = rq(M; kwargs...)
  l = ind(RM, 2)
  # TODO: simplify this by permuting inds(T) by (Lpos,Rpos)
  # then grab Linds,Rinds
  Rinds = NDTensors.similartype(IndsT, Val{NR})(ntuple(i -> inds(T)[Rpos[i]], Val(NR)))
  Qinds = NDTensors.pushfirst(Rinds, l)
  Q = reshape(QM, Qinds)
  Linds = NDTensors.similartype(IndsT, Val{NL})(ntuple(i -> inds(T)[Lpos[i]], Val(NL)))
  Linds = NDTensors.push(Linds, l)
  R = NDTensors.reshape(RM, Linds)
  return R, Q
end

function ql(A::ITensor, Linds...; kwargs...)
    tags::TagSet = get(kwargs, :tags, "Link,ql")
    Lis = commoninds(A, ITensors.indices(Linds))
    Ris = uniqueinds(A, Lis)
    Lpos, Rpos = NDTensors.getperms(inds(A), Lis, Ris)
    QT, LT = ql(ITensors.tensor(A), Lpos, Rpos; kwargs...)
    Q, L = itensor(QT), itensor(LT)
    q::Index = ITensors.commonind(Q, L)
    settags!(Q, tags, q)
    settags!(L, tags, q)
    q = settags(q, tags)
    #
    #  Do row removal for rank revealing LQ
    #
    epsrr::Float64 = get(kwargs, :epsrr , 1e-14)
    if epsrr>0.0 L,Q=trim(L,Q,epsrr) end
    return Q, L, q
end

function lq(A::ITensor, Rinds...; kwargs...)
  tags::TagSet = get(kwargs, :tags, "Link,lq")
  Ris = commoninds(A, ITensors.indices(Rinds))
  Lis = uniqueinds(A, Ris)
  Lpos, Rpos = NDTensors.getperms(inds(A), Lis, Ris)
  LT, QT = lq(ITensors.tensor(A), Lpos, Rpos; kwargs...)
  Q, L = itensor(QT), itensor(LT)
  q::Index = ITensors.commonind(Q, L)
  settags!(Q, tags, q)
  settags!(L, tags, q)
  q = settags(q, tags)
  #
  #  Do row removal for rank revealing LQ
  #
  epsrr::Float64 = get(kwargs, :epsrr , 1e-14)
  if epsrr>0.0 L,Q=trim(L,Q,epsrr) end
  return L, Q, q
end

function rq(A::ITensor, Rinds...; kwargs...)
  tags::TagSet = get(kwargs, :tags, "Link,rq")
  Ris = commoninds(A, ITensors.indices(Rinds))
  Lis = uniqueinds(A, Ris)
  Lpos, Rpos = NDTensors.getperms(inds(A), Lis, Ris)
  RT, QT = rq(ITensors.tensor(A), Lpos, Rpos; kwargs...)
  Q, R = itensor(QT), itensor(RT)
  q::Index = ITensors.commonind(Q, R)
  settags!(Q, tags, q)
  settags!(R, tags, q)
  q = settags(q, tags)
  #
  #  Do row removal for rank revealing RQ
  #
  epsrr::Float64 = get(kwargs, :epsrr , 1e-14)
  if epsrr>0.0 R,Q=trim(R,Q,epsrr) end
  return R, Q, q
end

function qr(A::ITensor, Rinds...; kwargs...)
  Q,R,q=ITensors.qr(A,Rinds...;kwargs...)
  #
  #  Do row removal for rank revealing RQ
  #
  epsrr::Float64 = get(kwargs, :epsrr , 1e-14)
  if epsrr>0.0 R,Q=trim(R,Q,epsrr) end
  return Q, R, q
end

function trim(R::ITensor,Q::ITensor,eps::Float64)
  iq=commonind(R,Q)
  zeros=find_zero_rows(R,iq,eps)
  if sum(zeros)==0
    return R,Q
  end
  #@printf "Rank Reveal removing %4i rows with epsrr=%.1e\n" sum(zeros) eps
  nq=dim(iq)-sum(zeros)
  iRo=noncommoninds(R,iq)
  iQo=noncommoninds(Q,iq)
  iqn=Index(nq,tags(iq))
  T=eltype(R)
  Rn=ITensor(T(0.0),iqn,iRo)
  Qn=ITensor(T(0.0),iQo,iqn)
  ivqn=1
  for ivq in eachindval(iq)
    if zeros[ivq.second]==false
      for ivRo in eachindval(iRo...)
        Rn[iqn=>ivqn,ivRo]=R[ivq,ivRo]
      end #for ivRo
      for ivQo in eachindval(iQo...)
        Qn[iqn=>ivqn,ivQo...]=Q[ivq,ivQo...]
      end #for ivQo
      ivqn+=1
    end #if zero
  end #for ivq
  return Rn,Qn
end

function find_zero_rows(R::ITensor,iq::Index,eps::Float64)::Array{Bool}
  zeros=falses(dim(iq))
  others=noncommoninds(R,iq)
  for iqv in eachindval(iq)
    s=0.0
    for io in eachindval(others...)
      s+=abs(R[iqv,io])
    end
    zeros[iqv.second]= (s<=eps)
  end
  return zeros
end
#=
is=Index(2,"Site,n=1")
ir=Index(4,"Link,l=0")
ic=Index(4,"Link,l=1")
iq=Index(4,"Link,qx")
#AQ=[i*1.0 for i in 1:dim(ic)*dim(iq)*dim(is)*dim(is)]
#Q=ITensor(eltype(AQ),AQ,ic,iq,is,is')
AQ=[i*1.0 for i in 1:dim(ic)*dim(iq)]
Q=ITensor(eltype(AQ),AQ,ic,iq)

A=[1:4;;5:8;;9:12;;13:16]
A[2,:].=0.0
R=ITensor(eltype(A),A,iq,ir)
z=find_zero_rows(R,iq,0.0)
Rn,Qn=trim(R,Q,1e-14)
@show Rn
@show Qn
=#