text:
"""

    findCOM(pos0, startatom=1, endatom=size(pos0, 2))

Return the unweighted center of mass of a block `pos0[:, startatom:endatom]` of coordinates,
where each column `pos0[:, p]` has 3 elements representing the Cartesian coordinates of a
particle labeled ``p``.

"""
function findCOM(pos0::AbstractArray{<:AbstractFloat, 2}, startatom::Integer=1,
                 endatom::Integer=size(pos0, 2))

    if (size(pos0, 1) != 3)
        error("findCOM(): pos0 must by a 3-by-N real matrix");
    end
    if (startatom < 1)
        error("findCOM(): startatom must be >= 1");
    end
    if (endatom > size(pos0, 2))
        error("findCOM(): endatom must be <= N = size(pos0, 2)");
    end
    if (startatom > endatom)
        error("findCOM(): startatom must be <= endatom");
    end
    currCOM = zeros(eltype(pos0), 3);
    @inbounds for pp=startatom:endatom
        @inbounds for jj=1:3
            currCOM[jj] += pos0[jj, pp];
        end
    end

    return currCOM./(endatom - startatom + 1);

end

"""

    translateCOM!(pos0, dispCOM, startatom=1, endatom=size(pos0, 2))

Translate the block `pos0[:, startatom:endatom]` of coordinates in-place (overwriting it),
where each column `pos0[:, p]` has 3 elements representing the Cartesian coordinates of a
particle labeled ``p`` and `dispCOM` is a 3-element vector representing the displacement,
such that for each `p` in the range, `pos0[:, p] .+= dispCOM`.

See also: [`translateCOM`](@ref)

"""
function translateCOM!(pos0::AbstractArray{FT, 2}, dispCOM::Array{FT, 1},
                       startatom::Integer=1,
                       endatom::Integer=size(pos0, 2)) where FT<:AbstractFloat

    if (length(dispCOM) != 3)
        error("translateCOM!(): dispCOM must be a 3-element real vector");
    end
    if (size(pos0, 1) != 3)
        error("translateCOM!(): pos0 must by a 3-by-N real matrix");
    end
    if (startatom < 1)
        error("translateCOM!(): startatom must be >= 1");
    end
    if (endatom > size(pos0, 2))
        error("translateCOM!(): endatom must be <= N = size(pos0, 2)");
    end
    if (startatom > endatom)
        error("translateCOM!(): startatom must be <= endatom");
    end

    if (dispCOM != zero(dispCOM))
        @inbounds for pp=startatom:endatom
            @inbounds for jj=1:3
                pos0[jj, pp] += dispCOM[jj];
            end
        end
    end

end

"""

    rotatecoords!(pos0, rotmat, startatom=1, endatom=size(pos0, 2))

Rotate the block `pos0[:, startatom:endatom]` of coordinates in-place (overwriting it),
where each column `pos0[:, p]` has 3 elements representing the Cartesian coordinates of a
particle labeled ``p`` and `rotmat` is a 3-by-3 matrix representing the (alibi) rotation,
such that for each `p` in the range, `pos0[:, p] .= rotmat * pos0[:, p]`.

See also: [`rotatecoords`](@ref)

!!! warning

    This function will not check whether `rotmat` is real and orthogonal. It is the
    responsibility of the caller to ensure that `rotmat` represents a physical rotation;
    otherwise, unexpected behavior could emerge.

"""
function rotatecoords!(pos0::AbstractArray{FT, 2}, rotmat::AbstractArray{FT, 2},
                       startatom::Integer=1,
                       endatom::Integer=size(pos0, 2)) where FT<:AbstractFloat
    if (size(rotmat, 1) != 3 || size(rotmat, 2) != 3)
        error("rotatecoords!(): rotmat must be a 3-by-3 orthogonal real matrix");
    end
    if (size(pos0, 1) != 3)
        error("rotatecoords!(): pos0 must by a 3-by-N real matrix");
    end
    if (startatom < 1)
        error("rotatecoords!(): startatom must be >= 1");
    end
    if (endatom > size(pos0, 2))
        error("rotatecoords!(): endatom must be <= N = size(pos0, 2)");
    end
    if (startatom > endatom)
        error("rotatecoords!(): startatom must be <= endatom");
    end

    if (rotmat != Matrix{FT}(I, 3, 3))
        newvec = zero(pos0[:, 1]);
        @inbounds for pp=startatom:endatom
            setindex!(newvec, rotmat * pos0[:, pp], 1:3);
            @inbounds for jj=1:3
                pos0[jj, pp] = newvec[jj];
            end
        end
    end

end

"""

    translateCOM(pos0, dispCOM, startatom=1, endatom=size(pos0, 2))

Translate the block `pos0[:, startatom:endatom]` of coordinates not in-place, returning the
result, where each column `pos0[:, p]` has 3 elements representing the Cartesian coordinates
of a particle labeled ``p`` and `dispCOM` is a 3-element vector representing the
displacement, such that for each `p` in the range, `pos0[:, p] .+= dispCOM`.

See also: [`translateCOM!`](@ref)

"""
function translateCOM(pos0::AbstractArray{FT, 2}, dispCOM::Array{FT, 1},
                      startatom::Integer=1,
                      endatom::Integer=size(pos0, 2)) where FT<:AbstractFloat
    B = deepcopy(pos0);
    translateCOM!(B, dispCOM, startatom, endatom);
    return B;
end

"""

    rotatecoords(pos0, rotmat, startatom=1, endatom=size(pos0, 2))

Rotate the block `pos0[:, startatom:endatom]` of coordinates not in-place, returning the
result, where each column `pos0[:, p]` has 3 elements representing the Cartesian coordinates
of a particle labeled ``p`` and `rotmat` is a 3-by-3 matrix representing the rotation,
such that for each `p` in the range, `pos0[:, p] .= rotmat * pos0[:, p]`.

!!! warning

    This function will not check whether `rotmat` is real and orthogonal. It is the
    responsibility of the caller to ensure that `rotmat` represents a physical rotation;
    otherwise, unexpected behavior could emerge.

See also: [`rotatecoords!`](@ref)

"""
function rotatecoords(pos0::AbstractArray{FT, 2}, rotmat::AbstractArray{FT, 2},
                      startatom::Integer=1,
                      endatom::Integer=size(pos0, 2)) where FT<:AbstractFloat
    B = deepcopy(pos0);
    rotatecoords!(B, rotmat, startatom, endatom);
    return B;
end

"""

    rotate3Nmat!(A, rotmat, startatom=1, endatom=div(size(A, 2), 3))

Rotate the contents in the diagonal block of the complex matrix
`A[3*startatom-2:3*endatom, 3*startatom-2:3*endatom]` by rotating each 3-by-3 block therein
by the 3-by-3 real orthogonal (alibi) rotation matrix `rotmat`; this is done in-place,
overwriting `A`.

See also: [`rotate3Nmat`](@ref)

!!! warning

    This function will not check whether `rotmat` is real and orthogonal. It is the
    responsibility of the caller to ensure that `rotmat` represents a physical rotation;
    otherwise, unexpected behavior could emerge.

"""
function rotate3Nmat!(A::AbstractArray{Complex{FT}, 2}, rotmat::AbstractArray{FT, 2},
                      startatom::Integer=1,
                      endatom::Integer=div(size(A, 2), 3)) where FT<:AbstractFloat

    if (size(rotmat, 1) != 3 || size(rotmat, 2) != 3)
        error("rotate3Nmat!(): rotmat must be a 3-by-3 orthogonal real matrix");
    end
    if (size(A, 1) != size(A, 2) || mod(size(A, 1), 3) != 0)
        error("rotate3Nmat!(): A must by a 3N-by-3N real or complex matrix");
    end
    if (startatom < 1)
        error("rotate3Nmat!(): startatom must be >= 1");
    end
    if (endatom > div(size(A, 2), 3))
        error("rotate3Nmat!(): endatom must be <= N = size(A, 2)/3");
    end
    if (startatom > endatom)
        error("rotate3Nmat!(): startatom must be <= endatom");
    end

    if (rotmat != Matrix{FT}(I, 3, 3))
        @inbounds for qq=startatom:endatom
            @inbounds for pp=startatom:endatom
                tempmat = view(A, 3*pp-2:3*pp, 3*qq-2:3*qq);
                tempmat .= rotmat * tempmat * transpose(rotmat);
            end
        end
    end

end

"""

    rotate3Nmat!(A, rotmat, startatom=1, endatom=div(size(A, 2), 3))

Rotate the contents in the diagonal block of the complex matrix
`A[3*startatom-2:3*endatom, 3*startatom-2:3*endatom]` by rotating each 3-by-3 block therein
by the 3-by-3 real orthogonal (alibi) rotation matrix `rotmat`; this is done not in-place,
returning the result.

See also: [`rotate3Nmat!`](@ref)

!!! warning

    This function will not check whether `rotmat` is real and orthogonal. It is the
    responsibility of the caller to ensure that `rotmat` represents a physical rotation;
    otherwise, unexpected behavior could emerge.

"""
function rotate3Nmat(A::AbstractArray{Complex{FT}, 2}, rotmat::AbstractArray{FT, 2},
                     startatom::Integer=1,
                     endatom::Integer=div(size(A, 2), 3)) where FT<:AbstractFloat
    B = deepcopy(A);
    rotate3Nmat!(B, rotmat, startatom, endatom);
    return B;
end

"""

    makeCOM0!(pos0, startatom=1, endatom=size(pos0, 2))

Translate the block `pos0[:, startatom:endatom]` of coordinates in-place (overwriting it) to
have a center of mass of 0, where each column `pos0[:, p]` has 3 elements representing the
Cartesian coordinates of a particle labeled ``p``.

See also: [`makeCOM0`](@ref)

"""
function makeCOM0!(pos0::AbstractArray{<:AbstractFloat, 2}, startatom::Integer=1,
                   endatom::Integer=size(pos0, 2))

    currCOM = findCOM(pos0, startatom, endatom);
    lmul!(-1, currCOM);
    translateCOM!(pos0, currCOM, startatom, endatom);

end

"""

    makeCOM0(pos0, startatom=1, endatom=size(pos0, 2))

Translate the block `pos0[:, startatom:endatom]` of coordinates not in-place, returning the
result, to have a center of mass of 0, where each column `pos0[:, p]` has 3 elements
representing the Cartesian coordinates of a particle labeled ``p``.

See also: [`makeCOM0!`](@ref)

"""
function makeCOM0(pos0::AbstractArray{<:AbstractFloat, 2}, startatom::Integer=1,
                  endatom::Integer=size(pos0, 2))

    currCOM = findCOM(pos0, startatom, endatom);
    lmul!(-1.0, currCOM);
    return translateCOM(pos0, currCOM, startatom, endatom);

end


text:
function add_two_digit(d1::Float64, d2::Float64)
    return d1 + d2
end

ext: jl


text:
export MtlBuffer, device, content, alloc, free, handle

const MTLBuffer = Ptr{MtBuffer}

# From docs: "MSL implements a buffer as a pointer to a built-in or user defined data type described in the
# device, constant, or threadgroup address space.
struct MtlBuffer{T} <: MtlResource
    handle::MTLBuffer
end

Base.unsafe_convert(::Type{MTLBuffer}, buf::MtlBuffer) = buf.handle

Base.:(==)(a::MtlBuffer, b::MtlBuffer) = a.handle == b.handle
Base.hash(buf::MtlBuffer, h::UInt) = hash(buf.handle, h)

## properties

Base.propertynames(o::MtlBuffer) = (
    :length,
    invoke(propertynames, Tuple{MtlResource}, o)...
)

function Base.getproperty(o::MtlBuffer, f::Symbol)
    if f === :length
        mtBufferLength(o)
    else
        invoke(getproperty, Tuple{MtlResource, Symbol}, o, f)
    end
end

Base.sizeof(buf::MtlBuffer)          = Int(buf.length)
Base.length(d::MtlBuffer{T}) where T = sizeof(d) ÷ sizeof(T)

content(buf::MtlBuffer{T}) where T   = Base.bitcast(Ptr{T}, mtBufferContents(buf))

## allocation

alloc_buffer(dev::MtlDevice, bytesize, opts::MtlResourceOptions) =
    mtDeviceNewBufferWithLength(dev, bytesize, opts)
alloc_buffer(dev::MtlHeap, bytesize, opts::MtlResourceOptions) =
    mtHeapNewBufferWithLength(heap, bytesize, opts)
alloc_buffer(dev::MtlDevice, bytesize, opts::MtlResourceOptions, ptr::Ptr) =
    mtDeviceNewBufferWithBytes(dev, ptr, bytesize, opts)
alloc_buffer(dev::MtlHeap, bytesize, opts::MtlResourceOptions, ptr::Ptr) =
    mtHeapNewBufferWithBytes(heap, ptr, bytesize, opts)

alloc_buffer(dev, bytesize, opts::Integer) =
    alloc_buffer(dev, bytesize, Base.bitcast(MtlResourceOptions, UInt32(opts)))
alloc_buffer(dev, bytesize, opts::Integer, ptr) =
    alloc_buffer(dev, bytesize, Base.bitcast(MtlResourceOptions, UInt32(opts)), ptr)

# TODO: buffer should be untyped
function MtlBuffer{T}(dev::Union{MtlDevice,MtlHeap},
                      length::Integer;
                      storage = Private,
                      hazard_tracking = DefaultTracking,
                      cache_mode = DefaultCPUCache) where {T}
    opts = storage | hazard_tracking | cache_mode

    bytesize = length * sizeof(T)
    ptr = alloc_buffer(dev, bytesize, opts)

    return MtlBuffer{T}(ptr)
end

function MtlBuffer{T}(dev::Union{MtlDevice,MtlHeap},
                      length::Integer,
                      ptr::Ptr;
                      storage = Managed,
                      hazard_tracking = DefaultTracking,
                      cache_mode = DefaultCPUCache) where {T}
    storage == Private && error("Can't create a Private copy-allocated buffer.")
    opts =  storage | hazard_tracking | cache_mode

    bytesize = length * sizeof(T)
    ptr = alloc_buffer(dev, bytesize, opts, ptr)

    dev = dev isa MtlDevice ? dev : dev.device
    return MtlBuffer{T}(ptr)
end

MtlBuffer(T::Type, dev::Union{MtlDevice,MtlHeap}, args...; kwargs...) =
    MtlBuffer{T}(dev, args...;kwargs...)
MtlBuffer(dev::Union{MtlDevice,MtlHeap}, args...; kwargs...) =
    MtlBuffer(Cvoid, dev, args...;kwargs...)

"""
    alloc(T, device, length, [ptr=nothing]; storage=Default, hazard_tracking=Default, chache_mode=Default)
    MtlBuffer{T}(device, length...)

Allocates a Metal Buffer on `device` of bytes equal to `length * sizeof(T)`. If
a CPU-pointer is passed as last argument, then the Metal Buffer is initialized
with the content of the memory starting at `ptr`, otherwise it's
zero-initialized.

! Note: You are responsible for freeing the returned buffer

The storage kwarg controls where the buffer is stored. Possible values are:
 - Private : Residing on the device
 - Shared  : Residing on the host
 - Managed : Keeps two copies of the buffer, on device and on host. Explicit
   calls must be given to syncronize the two
 - Memoryless : an iOs specific thing that won't work on Mac.

Note that `Private` buffers can't be directly accessed from the CPU, therefore
you cannot use this option if you pass a ptr to initialize the memory.
"""
alloc(args...; kwargs...) = MtlBuffer(args...; kwargs...)

"""
    free(buffer::MtlBuffer)

Frees the buffer if the handle is valid.
This does not protect against double-freeing of the same buffer!
"""
free(buf::MtlBuffer) = mtRelease(buf.handle)

"""
    DidModifyRange!(buf::MtlBuffer{T}, range::UnitRange)

Notifies the GPU that the range of elements `range`, corresponding to the bytes
`sizeof(T)*first(range):sizeof(T)*last(range)` have been modified on the CPU,
and that they should be transferred to the device before executing any following
command.

Only valid for `Managed` buffers.
"""
function DidModifyRange!(buf::MtlBuffer{T}, range::UnitRange) where {T}
    mtBufferDidModifyRange(buf, range*sizeof(T))
end

# Views on different device
NewBuffer(buf::MtlBuffer, d::MtlDevice) =
    mtBufferNewRemoteBufferViewForDevice(buf, d);

function ParentBuffer(buf::MtlBuffer)
    orig = mtBufferRemoteStorageBuffer(buf);
    if orig == C_NULL
        return nothing
    else
        return MtlBuffer(orig)
    end
end

## interop with CPU array

function Base.unsafe_wrap(t::Type{<:Array}, buf::MtlBuffer{T}, dims; own=false) where {T}
    ptr = content(buf)
    ptr == C_NULL && error("Can't unsafe_wrap a GPU Private array.")
    return unsafe_wrap(t, ptr, dims; own=own)
end

handle_array(vec::Vector{<:MtlBuffer}) = [buf.handle for buf in vec]


text:
using egads

# This is the tutorial case provided with EGADS

function periodicSeam(eedge::ego,   # Edge associated with seam
                      sense::Cint ) # sense of PCurve

    status = EGADS_SUCCESS
    range = fill(0.0,2)
    data  = fill(0.0,4)
    per   = Ref{Cint}(0)

    (context,status) = EG_getContext(eedge)
    if (status != EGADS_SUCCESS) return status end

    status = EG_getRange(eedge, range, per)
    if (status != EGADS_SUCCESS) return status end

    # set up u and v at ends
    data = [2*pi, range[1], 0.0, sense]
    if (sense == -1) data[2] = range[2] end

    # make (linear) PCurve
    (pcurve,status) =  EG_makeGeometry(context, PCURVE, LINE, NULL_E, C_NULL, data)

    return (pcurve, status)

end

# set the parameter values
width  =  5.0
minrad =  8.0
maxrad = 12.0
fillet =  2.0
thick  =  0.5
bolts  =  5
crad   =  5.0
brad   =  1.0

# initialize arrays
enodes   = fill( ego(0), 8 )
ecurve   = fill( ego(0), 16 )
eedges   = fill( ego(0), 16 )
elist    = fill( ego(0), 20 )
esurface = fill( ego(0), 4 )
efaces   = fill( ego(0), 8 )
epcurve  = fill( ego(0), 4 )
trange   = fill( 0.0, 2 )

# initialize pointers
NULL_E = ego(0)
NULL_I = Ref{Cint}(0)
NULL_D = Ref{Cdouble}(0.0)

(context, status) = EG_open( )
if (status < 0) error("Can't open, failure code: %i", status) end

# node locations
node1 = [-minrad, 0.0, -width/2.0]
node2 = [-minrad, 0.0,  width/2.0]
node3 = [-maxrad, 0.0,  width/2.0]
node4 = [-maxrad, 0.0, -width/2.0]
node5 = [ minrad, 0.0, -width/2.0]
node6 = [ maxrad, 0.0, -width/2.0]
node7 = [ maxrad, 0.0,  width/2.0]
node8 = [ minrad, 0.0,  width/2.0]

# make nodes

(enodes[1], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node1, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[2], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node2, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[3], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node3, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[4], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node4, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[5], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node5, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[6], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node6, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[7], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node7, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(enodes[8], status) = EG_makeTopology(context, NULL_E, NODE, Cint(0), node8, Cint(0), NULL_E, NULL_I)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 1
data = [node1[1], node1[2], node1[3], node2[1] - node1[1], node2[2] - node1[2], node2[3] - node1[3] ]
(ecurve[1], status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1], status) = EG_invEvaluate(ecurve[1], node1, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[1], node2, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[1]
elist[2] = enodes[2]
(eedges[1],status) = EG_makeTopology(context, ecurve[1], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 2
data = [ node2[1], node2[2], node2[3], node3[1] - node2[1], node3[2] - node2[2], node3[3] - node2[3] ]
(ecurve[2], status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[2], node2, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[2], node3, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[2]
elist[2] = enodes[3]
(eedges[2],status) = EG_makeTopology(context, ecurve[2], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 3
data = [ node3[1], node3[2], node3[3], node4[1] - node3[1], node4[2] - node3[2], node4[3] - node3[3] ]
(ecurve[3], status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[3], node3, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[3], node4, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[3]
elist[2] = enodes[4]
(eedges[3],status) = EG_makeTopology(context, ecurve[3], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 4
data = [ node4[1], node4[2], node4[3], node1[1] - node4[1], node1[2] - node4[2], node1[3] - node4[3] ]
(ecurve[4],status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[4], node4, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[4], node1, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[4]
elist[2] = enodes[1]
(eedges[4],status) = EG_makeTopology(context, ecurve[4], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 5
data = [node5[1], node5[2], node5[3], node6[1] - node5[1], node6[2] - node5[2], node6[3] - node5[3] ]
(ecurve[5],status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[5], node5, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[5], node6, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[5]
elist[2] = enodes[6]
(eedges[5],status) = EG_makeTopology(context, ecurve[5], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 6
data = [node6[1], node6[2], node6[3], node7[1] - node6[1], node7[2] - node6[2], node7[3] - node6[3] ]
(ecurve[6],status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[6], node6, data);
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[6], node7, data);
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[6]
elist[2] = enodes[7]
(eedges[6],status) = EG_makeTopology(context, ecurve[6], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 7
data = [node7[1], node7[2], node7[3], node8[1] - node7[1], node8[2] - node7[2], node8[3] - node7[3] ]
(ecurve[7],status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[7], node7, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[7], node8, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[7]
elist[2] = enodes[8]
(eedges[7],status) = EG_makeTopology(context, ecurve[7], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (linear) Edge 8
data = [node8[1], node8[2], node8[3], node5[1] - node8[1], node5[2] - node8[2], node5[3] - node8[3] ]
(ecurve[8],status) = EG_makeGeometry(context, CURVE, LINE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[8], node8, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[8], node5, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

elist[1] = enodes[8]
elist[2] = enodes[5]
(eedges[8],status) = EG_makeTopology(context, ecurve[8], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# data used in creating the arcs
axis1 = [1.0, 0.0, 0.0]
axis2 = [0.0, 1.0, 0.0]
axis3 = [0.0, 0.0, 1.0]

cent1 = [0.0, 0.0, -width / 2.0]
cent2 = [0.0, 0.0,  width / 2.0]

# make (circular) Edge 9
data = [cent1[1], cent1[2], cent1[3],
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3],
        minrad]
(ecurve[9],status) = EG_makeGeometry(context, CURVE, CIRCLE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[9], node5, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[9], node1, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end # ensure trange[2] > trange[1]

elist[1] = enodes[5];
elist[2] = enodes[1];
(eedges[9],status) = EG_makeTopology(context, ecurve[9], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (circular) Edge 10
data = [cent2[1], cent2[2], cent2[3],
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3],
        minrad]
(ecurve[10],status) = EG_makeGeometry(context, CURVE, CIRCLE, NULL_E, NULL_I, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[10], node8, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[10], node2, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end # ensure trange[2] > trange[1]

elist[1] = enodes[8]
elist[2] = enodes[2]
(eedges[10],status) = EG_makeTopology(context, ecurve[9], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (circular) Edge 11
data = [cent1[1], cent1[2], cent1[3],
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3], maxrad]
(ecurve[11],status) = EG_makeGeometry(context, CURVE, CIRCLE, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[11], node6, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[11], node4, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end  # ensure trange[2] > trange[1]

elist[1] = enodes[6]
elist[2] = enodes[4]
(eedges[11],status) = EG_makeTopology(context, ecurve[11], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (circular) Edge 12
data = [cent2[1], cent2[2], cent2[3],
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3], maxrad]
(ecurve[12], status) = EG_makeGeometry(context, CURVE, CIRCLE, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[1],status) = EG_invEvaluate(ecurve[12], node7, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[12], node3, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end   # ensure trange[2] > trange[1]

elist[1] = enodes[7]
elist[2] = enodes[3]
(eedges[12],status) = EG_makeTopology(context, ecurve[12], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make the outer cylindrical surface
data = [cent1[1], cent1[2], cent1[3],
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3],
        axis3[1], axis3[2], axis3[3], maxrad]
(esurface[1],status) = EG_makeGeometry(context, SURFACE, CYLINDRICAL, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(esurface[2],status) = EG_makeGeometry(context, SURFACE, CYLINDRICAL, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make the inner cylindrical surface
data = [cent1[1], cent1[2], cent1[3],
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3],
        axis3[1], axis3[2], axis3[3], minrad]
(esurface[3],status) = EG_makeGeometry(context, SURFACE, CYLINDRICAL, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(esurface[4],status) = EG_makeGeometry(context, SURFACE, CYLINDRICAL, NULL_E, C_NULL, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (planar) Face 1
sense = [SFORWARD, SREVERSE, SFORWARD, SFORWARD]
elist[1] = eedges[4]
elist[2] = eedges[9]
elist[3] = eedges[5]
elist[4] = eedges[11]
(eloop,status) = EG_makeTopology(context, NULL_E, LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[1],status) = EG_makeFace(eloop, SFORWARD, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (planar) Face 2
sense = [SFORWARD, SREVERSE, SFORWARD, SFORWARD]
elist[1] = eedges[2]
elist[2] = eedges[12]
elist[3] = eedges[7]
elist[4] = eedges[10]
(eloop,status) = EG_makeTopology(context, NULL_E, LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[2],status) = EG_makeFace(eloop, SFORWARD, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (cylindrical) Face 3
(epcurve[1],status) = EG_otherCurve(esurface[1], ecurve[3], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[2],status) = EG_otherCurve(esurface[1], ecurve[11], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[3],status) = EG_otherCurve(esurface[1], ecurve[6], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[4],status) = EG_otherCurve(esurface[1], ecurve[12], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

sense = [SFORWARD, SREVERSE, SFORWARD, SFORWARD]
elist = [eedges[3],  eedges[11], eedges[6],  eedges[12],
         epcurve[1], epcurve[2], epcurve[3], epcurve[4]]
(eloop,status) = EG_makeTopology(context, esurface[1], LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[3],status) = EG_makeTopology(context, esurface[1], FACE, SREVERSE, C_NULL, Cint(1), Ref{ego}(eloop), sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (cylindrical) Face 4
(epcurve[1],status) = EG_otherCurve(esurface[2], ecurve[1], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[2],status) = EG_otherCurve(esurface[2], ecurve[10], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[3],status) = EG_otherCurve(esurface[2], ecurve[8], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[4],status) = EG_otherCurve(esurface[2], ecurve[9], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

sense = [SFORWARD, SREVERSE, SFORWARD, SFORWARD]
elist = [eedges[1],  eedges[10], eedges[8], eedges[9],
         epcurve[1], epcurve[2], epcurve[3], epcurve[4]]
(eloop,status) = EG_makeTopology(context, esurface[2], LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[4],status) = EG_makeTopology(context, esurface[2], FACE, SFORWARD, C_NULL, Cint(1), Ref{ego}(eloop), sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("edge 12 done")
# make (circular) Edge 13
ecurve[13] = ecurve[9] # reuse

(trange[1],status) = EG_invEvaluate(ecurve[13], node1, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[13], node5, data);
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end   # ensure trange[2] > trange[1]

elist[1] = enodes[1]
elist[2] = enodes[5]
(eedges[13],status) = EG_makeTopology(context, ecurve[13], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("edge 13 done")
# make (circular) Edge 14
ecurve[14] = ecurve[10] # reuse

(trange[1],status) = EG_invEvaluate(ecurve[14], node2, data);
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[14], node8, data);
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end   # ensure trange[2] > trange[1]

elist[1] = enodes[2]
elist[2] = enodes[8]
(eedges[14],status) = EG_makeTopology(context, ecurve[14], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("edge 14 done")
# make (circular) Edge 15
ecurve[15] = ecurve[11]

(trange[1],status) = EG_invEvaluate(ecurve[14], node4, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[14], node6, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end   # ensure trange[2] > trange[1]

elist[1] = enodes[4]
elist[2] = enodes[6]
(eedges[15],status) = EG_makeTopology(context, ecurve[15], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("edge 15 done")
# make (circular) Edge 16
ecurve[16] = ecurve[12]

(trange[1],status) = EG_invEvaluate(ecurve[16], node3, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(trange[2],status) = EG_invEvaluate(ecurve[16], node7, data)
if (status != EGADS_SUCCESS) cleanup(status, context) end

if (trange[1] > trange[2]) trange[2] += 2*pi end   # ensure trange[2] > trange[1]

elist[1] = enodes[3]
elist[2] = enodes[7]
(eedges[16],status) = EG_makeTopology(context, ecurve[16], EDGE, TWONODE, trange, Cint(2), elist, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("edge 16 done")
# make (planar) Face 5
sense = [SFORWARD, SREVERSE, SREVERSE, SREVERSE]
elist = [eedges[15], eedges[5], eedges[13], eedges[4]]
(eloop,status) = EG_makeTopology(context, NULL_E, LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[5],status) = EG_makeFace(eloop, SFORWARD, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("face 5 done")
# make (planar) Face 6
sense = [SFORWARD, SREVERSE, SFORWARD, SFORWARD]
elist = [eedges[7], eedges[14], eedges[2], eedges[16]]
(eloop,status) = EG_makeTopology(context, NULL_E, LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[6],status) = EG_makeFace(eloop, SFORWARD, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("face 6 done")
# make (cylindrical) Face 7
(epcurve[1],status) = EG_otherCurve(esurface[3], ecurve[3], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[2],status) = EG_otherCurve(esurface[3], ecurve[15], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end
println("face 7 done")
# EG_otherCurve cannot be used along periodic seam
(epcurve[3],status) = periodicSeam(eedges[6], SFORWARD)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[4],status) = EG_otherCurve(esurface[3], ecurve[16], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

sense = [SFORWARD, SFORWARD, SFORWARD, SREVERSE]
elist = [eedges[3],  eedges[15], eedges[6],  eedges[16],
         epcurve[1], epcurve[2], epcurve[3], epcurve[4]]
(eloop,status) = EG_makeTopology(context, esurface[3], LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[7],status) = EG_makeTopology(context, esurface[3], FACE, SFORWARD, C_NULL, Cint(1), Ref{ego}(eloop), sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make (cylindrical) Face 8
(epcurve[1],status) = EG_otherCurve(esurface[4], ecurve[1], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[2],status) = EG_otherCurve(esurface[4], ecurve[14], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# EG_otherCurve cannot be used along periodic seam
(epcurve[3],status) = periodicSeam(eedges[8], SREVERSE)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(epcurve[4],status) = EG_otherCurve(esurface[4], ecurve[13], 0.0)
if (status != EGADS_SUCCESS) cleanup(status, context) end

sense = [SFORWARD, SFORWARD, SFORWARD, SREVERSE]
elist = [eedges[1],  eedges[14], eedges[8],  eedges[13],
         epcurve[1], epcurve[2], epcurve[3], epcurve[4]]
(eloop,status) = EG_makeTopology(context, esurface[4], LOOP, CLOSED, C_NULL, Cint(4), elist, sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(efaces[8],status) = EG_makeTopology(context, esurface[4], FACE, SREVERSE, C_NULL, Cint(1), Ref{ego}(eloop), sense)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# make the shell and initial Body
(eshell,status) = EG_makeTopology(context, NULL_E, SHELL, CLOSED, C_NULL, Cint(8), efaces, C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

(ebody1,status) = EG_makeTopology(context, NULL_E, BODY, SOLIDBODY, C_NULL, Cint(1), Ref{ego}(eshell), C_NULL)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# add fillets if desired (result is ebody2)
if (fillet > 0)
    elist = [ eedges[11], eedges[12], eedges[15], eedges[16] ]
    (ebody,status) = EG_filletBody(ebody1, Cint(4), elist, fillet, C_NULL)
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    status = EG_deleteObject(ebody1);
    if (status != EGADS_SUCCESS)
        cleanup(status, context)
    else
        ebody2 = ebody1
    end
end
# add wheel if desired (result is source)
if (thick > 0)
    data = [ 0.0, 0.0,  thick / 2.0,
             0.0, 0.0, -thick / 2.0,
             (minrad + maxrad) / 2.0 ]
    (ebody3,status) = EG_makeSolidBody(context, CYLINDER, data)
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    (emodel,status) = EG_solidBoolean(ebody2, ebody3, FUSION)
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    status = EG_deleteObject(ebody2)
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    status = EG_deleteObject(ebody3)
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    status = EG_getTopology(emodel, &eref, &oclass, &mtype, data, &nchild, &echilds, &senses);
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    if (oclass != MODEL || nchild != 1)
      @printf("You didn't input a model or are returning more than one body ochild = %d, nchild = %d/n", oclass, nchild);
      status = -999
      cleanup(status, context)
    end

    (source, status) = EG_copyObject(echilds[1], C_NULL)
    if (status != EGADS_SUCCESS) cleanup(status, context) end

    status = EG_deleteObject(emodel);
    cleanup(status, context)

    # add bolt holes
    for i in 1:bolts
        data[1] =  crad * cos(i * (2*pi / bolts))
        data[2] =  crad * sin(i * (2*pi / bolts))
        data[3] =  thick / 2.0
        data[4] =  crad * cos(i * (2*pi / bolts))
        data[5] =  crad * sin(i * (2*pi / bolts))
        data[6] = -thick / 2.0
        data[7] =  brad

        (ebody4,status) = EG_makeSolidBody(context, CYLINDER, data)
        if (status != EGADS_SUCCESS) cleanup(status, context) end

        (emodel,status) = EG_solidBoolean(source, ebody4, SUBTRACTION)
        if (status != EGADS_SUCCESS) cleanup(status, context) end

        status = EG_deleteObject(source)
        if (status != EGADS_SUCCESS) cleanup(status, context) end

        status = EG_deleteObject(ebody4)
        if (status != EGADS_SUCCESS) cleanup(status, context) end

        status = EG_getTopology(emodel, &eref, &oclass, &mtype, data, &nchild, &echilds2, &senses);
        if (status != EGADS_SUCCESS) cleanup(status, context) end

        if (oclass != MODEL || nchild != 1)
            @printf("You didn't input a model or are returning more than one body ochild = %d, nchild = %d/n", oclass, nchild);
            status = -999;
            cleanup(status, context)
        end

        (source,status) = EG_copyObject(echilds2[0], C_NULL)
        if (status != EGADS_SUCCESS) cleanup(status, context) end

        status = EG_deleteObject(emodel)
        if (status != EGADS_SUCCESS) cleanup(status, context) end
    end
    ebody = source
else
    ebody = ebody2
end

# make and dump the model
(emodel, status) = EG_makeTopology(context, C_NULL, MODEL, 0, C_NULL, 1, Ref{ego}(ebody), C_NULL)
printf("EG_makeTopology -> status=%d\n", status);
if (status < 0) exit(EXIT_FAILURE) end

status = EG_saveModel(emodel, "tire.egads");
printf("EG_saveModel -> status=%d\n", status);
if (status < 0) exit(EXIT_FAILURE) end

# cleanup
status = EG_deleteObject(emodel);
printf("EG_close -> status=%d\n", status);


text:
#= src/prange.jl

References:

  B. Engquist, L. Ying. A fast directional algorithm for high frequency
    acoustic scattering in two dimensions. Commun. Math. Sci. 7 (2): 327-345,
    2009.

  N. Halko, P.G. Martinsson, J.A. Tropp. Finding structure with randomness:
    Probabilistic algorithms for constructing approximate matrix
    decompositions. SIAM Rev. 53 (2): 217-288, 2011.
=#

for sfx in ("", "!")
  f = Symbol("prange", sfx)
  g = Symbol("pqrfact", sfx)
  @eval begin
    function $f(
        trans::Symbol, A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T);
        args...) where T
      prange_chktrans(trans)
      opts = copy(opts; args...)
      chkopts!(opts, A)
      if trans == :b
        checksquare(A)
        ishermitian(A) && return $f(:n, A, opts)
        opts.pqrfact_retval = "qr"
        if opts.sketch == :none
          Fr = pqrfact!(A', opts)
          Fc =       $g(A , opts)
        else
          Fr = sketchfact(:right, :c, A, opts)
          Fc = sketchfact(:right, :n, A, opts)
        end
        kr = Fr[:k]
        kc = Fc[:k]
        B = Array{T}(undef, size(A,1), kr+kc)
        B[:,   1:kr   ] = Fr[:Q]
        B[:,kr+1:kr+kc] = Fc[:Q]
        Rr = view(Fr.R, 1:kr, 1:kr)
        Rc = view(Fc.R, 1:kc, 1:kc)
        rmul!(view(B,:,   1:kr   ), UpperTriangular(Rr))
        rmul!(view(B,:,kr+1:kr+kc), UpperTriangular(Rc))
        opts.pqrfact_retval="q"
        return pqrfact_backend!(B, opts)[:Q]
      else
        opts.pqrfact_retval = "q"
        if opts.sketch == :none
          if trans == :n  Q =       $g(A , opts)[:Q]
          else            Q = pqrfact!(A', opts)[:Q]
          end
        elseif opts.sketch == :sub  Q = prange_sub(trans, A, opts)
        else                        Q = sketchfact(:right, trans, A, opts)[:Q]
        end
        Q
      end
    end
    $f(trans::Symbol, A, args...; kwargs...) =
      $f(trans, LinOp(A), args...; kwargs...)
    $f(A, args...; kwargs...) = $f(:n, A, args...; kwargs...)
  end
end

function prange_sub(trans::Symbol, A::AbstractMatrix{T}, opts::LRAOptions) where T
  F = sketchfact(:left, trans, A, opts)
  k = F[:k]
  if trans == :n
    B = A[:,F[:p][1:k]]
  else
    n = size(A, 2)
    B = Array{T}(undef, n, k)
    @inbounds for j = 1:k, i = 1:n
      B[i,j] = conj(A[F[:p][j],i])
    end
  end
  orthcols!(B)
end

prange_chktrans(trans::Symbol) =
  trans in (:n, :c, :b) || throw(ArgumentError("trans"))


text:
# Batch posterior samples
# Enable multiple threads and sysimage by starting julia using:
# JULIA_NUM_THREADS=10 julia --sysimage=/home/student/.julia/environments/v1.4/JuliaSysimage.so

using EstimateSDE
using Serialization
using Base.Threads

## Batch settings
f = joinpath(@__DIR__,"batch-settings.toml")
setttings_batch = batchsettings(f)

## Maximum likelihood results
dd = EstimateSDE.resultsdir("ICMLA",setttings_batch.name,string(setttings_batch.proclabel))
datafilename = joinpath(dd,"maximum-likelihood.jld")
println(datafilename)
output_mle = deserialize(datafilename)
@info "Read maximum likelihood results from $datafilename"

## Settings
settings = (
  tor_kwargs=(N=200,Nwarmup=200,nchains=1),
  Nest=1 # Nothing: do all
)
printstyled("\nSettings posterior samples\n",color=:blue)
println(settings)

## Batch posterior samples
ENV["JULIA_DEBUG"] = EstimateSDE
if nthreads()==1
  @warn "nthreads==1. Increase number of threads for faster processing"
end
output_posterior_samples = @time batch_posterior_samples(settings,output_mle)
nothing


text:
# # Analog forecasting

using Plots, DifferentialEquations, AnalogDataAssimilation

σ = 10.0
ρ = 28.0
β = 8.0 / 3

dt_integration = 0.01
dt_states = 1
dt_obs = 8
var_obs = [1]
nb_loop_train = 100
nb_loop_test = 10
sigma2_catalog = 0.0
sigma2_obs = 2.0

ssm = StateSpaceModel(
    lorenz63,
    dt_integration,
    dt_states,
    dt_obs,
    [σ, ρ, β],
    var_obs,
    nb_loop_train,
    nb_loop_test,
    sigma2_catalog,
    sigma2_obs,
)

# # compute u0 to be in the attractor space

u0 = [8.0; 0.0; 30.0]
tspan = (0.0, 5.0)
prob = ODEProblem(ssm.model, u0, tspan, ssm.params)
u0 = last(solve(prob, reltol = 1e-6, save_everystep = false))

# # Generate the data

xt, yo, catalog = generate_data(ssm, u0)

# # Create the forecasting function

af = AnalogForecasting(50, xt, catalog; regression = :local_linear, sampling = :multinomial)

# # Data assimilation

np = 100
DA = DataAssimilation(af, xt, ssm.sigma2_obs)
x̂ = forecast(DA, yo, AnEnKS(np), progress = false);
println(RMSE(xt, x̂))

# # Plot results

plot(xt.t, xt[1], label = "true")
plot!(xt.t, x̂[1], label = "forecasted")
scatter!(yo.t, yo[1], markersize = 2, label = "observed")


text:
"""
    pplus, eplus, pminus, eminus = dssplit2d(u, ABCD, p)

Split a convex polygon `p` into "plus" and "minus" polygons.
C(`pplus`) + D₁`u` ≥ 0, C(`pminus`) + D₁`u` ≤ 0.

`p` is given as a sequential list of vertices with the first vertex replicated
at the end of the list. `ABCD` describes the modulator structure, and `u` is the
modulator input.
"""
function dssplit2d(u, ABCD, p)
    n = size(ABCD, 1) - 1
    C = ABCD[n+1, 1:n]
    D1 = ABCD[n+1, n+1]	# D2=ABCD(n+1,n+2) must be zero
    N = size(p, 2)
    D1u = D1 * u
    y = p'*C + D1u*ones(N)

    sign1 = sign(y[1])
    i = findall(x -> sign(x) != sign1, y)
    i1 = i[1]       # First change of sign.
    pa = dscut(p[:,i1-1], y[i1-1], p[:,i1], y[i1])
    i2 = i[length(i)]   # Second change of sign.
    pb = dscut(p[:,i2], y[i2], p[:,i2+1], y[i2+1])
    if sign1 > 0
        pminus = [pa p[:,i] pb pa]
        pplus = [p[:,1:i1-1] pa pb p[:,i2+1:N]]
    else
        pplus = [pa p[:,i] pb pa ]
        pminus = [p[:,1:i1-1] pa pb p[:,i2+1:N]]
    end
    ne = size(pplus, 2)
    eplus = [(1:ne)'; (2:ne)' 1]
    ne = size(pminus, 2)
    eminus = [(1:ne)'; (2:ne)' 1]

    return pplus, eplus, pminus, eminus
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{R, D},
                           ::Type{S};
                           kwargs...) where {R<:PSY.RenewableGen,
                                             D<:AbstractRenewableDispatchFormulation,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(R, sys)

    if validate_available_devices(devices, R)
        return
    end

    #Variables
    activepower_variables!(canonical, devices);

    reactivepower_variables!(canonical, devices);

    #Constraints
    activepower_constraints!(canonical, devices, D, S)

    reactivepower_constraints!(canonical, devices, D, S)

    feedforward!(canonical, R, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{R, D},
                           ::Type{S};
                           kwargs...) where {R<:PSY.RenewableGen,
                                             D<:AbstractRenewableDispatchFormulation,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(R, sys)

    if validate_available_devices(devices, R)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    #Constraints
    activepower_constraints!(canonical, devices, D, S)

    feedforward!(canonical, R, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{R, RenewableFixed},
                           system_formulation::Type{S};
                           kwargs...) where {R<:PSY.RenewableGen,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(R, sys)

    if validate_available_devices(devices, R)
        return
    end

    nodal_expression!(canonical, devices, system_formulation)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{PSY.RenewableFix, D},
                           system_formulation::Type{S};
                           kwargs...) where {D<:AbstractRenewableDispatchFormulation,
                                             S<:PM.AbstractPowerModel}

    @warn("The Formulation $(D) only applies to FormulationControllable Renewable Resources, \n Consider Changing the Device Formulation to RenewableFixed")

    construct_device!(canonical,
                      sys,
                      DeviceModel(PSY.RenewableFix,RenewableFixed),
                      system_formulation;
                      kwargs...)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{PSY.RenewableFix, RenewableFixed},
                           system_formulation::Type{S};
                           kwargs...) where {S<:PM.AbstractPowerModel}

    devices = PSY.get_components(PSY.RenewableFix, sys)

    if validate_available_devices(devices, PSY.RenewableFix)
        return
    end

    nodal_expression!(canonical, devices, system_formulation)

    return

end


text:
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport GLS_Apr_weekday_AM

xi_list = GLS_Apr_weekday_AM.xi_list
P = GLS_Apr_weekday_AM.P
L = GLS_Apr_weekday_AM.L  # dimension of xi

size(P, 1), size(P, 2)

using JuMP

mGLSJulia = Model()

@defVar(mGLSJulia, lam[1:size(P,1)] >= 0)

@defVar(mGLSJulia, p[1:size(P,1), 1:size(P,2)] >= 0)

for i = 1:size(P,1)
    for j = 1:size(P,2)
        if P[i,j] == 0
            @addConstraint(mGLSJulia, p[i,j] == 0)
        end
    end
end

for i = 1:size(P,1)
    @addNLConstraint(mGLSJulia, sum{p[i,j], j = 1:size(P,2)} == 1)
end

for l = 1:L
    @addNLConstraint(mGLSJulia, sum{p[i,l] * lam[i], i = 1:size(P,1)} == xi_list[l])
end

@setNLObjective(mGLSJulia, Min, sum{p[1,j], j = 1:size(P,2)})  # play no actual role, but could not use zero objective

solve(mGLSJulia)

getValue(lam)

getObjectiveValue(mGLSJulia)

GLS_Apr_weekday_AM.saveDemandVec(getValue(lam))

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
const plotly_version = "1.53.0"

function __init__()
  is_ijulia = isdefined(Main, :IJulia) && Main.IJulia.inited
  is_juno = isdefined(Main, :Juno) && Main.Juno.isactive()
  is_vscode = isdefined(Main, :_vscodeserver)

  is_ide = is_ijulia || is_juno || is_vscode
  global is_repl = !is_ide

  if is_repl
    work_size = (62, 16)

    work_palette = [
      "green", "blue", "red",
      "yellow", "cyan", "magenta"
    ]
  else
    work_size = (600, 400)

    work_palette = [
      "#636EFA", "#EF553B", "#00CC96", "#AB63FA", "#FFA15A",
      "#19D3F3", "#FF6692", "#B6E880", "#FF97FF", "#FECB52"
    ]
  end

  global _palette = work_palette
  global default_plot_size = work_size

  global interact_comms = Dict{AbstractString, Observable}()
  empty!(widget_observers)

  SimplePlot() # sets global _plots variable

  if is_repl
    include(joinpath(@__DIR__, "repl.jl"))
  else
    _init_notebook()
  end
end

function _init_notebook()
  assets_url = "https://cdn.jsdelivr.net/gh/djsegal/SimplePlots.jl/dist/"

  plotly_url = "https://cdn.plot.ly/plotly-"
  plotly_url *= "$(plotly_version).min"

  plotly_javascript = """
    <script type="text/javascript" class="js-plotly-script">
      globalComms = {};
      demoData = {};

      \$(".js-plotly-script").parent().css('padding', 0);

      function customPlotlyReact(curPlot, curData, curLayout, curConfig) {
        for (var i = 0; i < curData.length; i++) {
          curData[i]["visible"] = true;
          for (var j = 0; j < curPlot.data.length; j++) {
            if ( curData[i].name !== curPlot.data[j].name ) { continue; }
            if ( curPlot.data[j].visible !== "legendonly" ) { continue; }

            curData[i]["visible"] = "legendonly";
          }
        }

        Plotly.react(curPlot, curData, curLayout, curConfig);
      }

      function customPlotLoader(curCallback) {
        if ( \$(".js-nouislider-css").length == 0 ) {
          \$("head").append(
            '<link class="js-nouislider-css" href="$( assets_url * "nouislider.min.css" )" rel="stylesheet">'
          );
        }

        if ( \$(".js-custom-css").length == 0 ) {
          \$("head").append(
            '<link class="js-custom-css" href="$( assets_url * "custom.min.css" )" rel="stylesheet">'
          );
        }

        if ( typeof Plotly !== "undefined" ) {
          if ( typeof curCallback !== "undefined" ) {
            curCallback();
          }
          return;
        }

        var plotlyScripts = document.getElementsByClassName("js-plotly-script");

        for (var i = 0; i < plotlyScripts.length; i++) {
          var scriptParent = plotlyScripts[i].parentElement;
          scriptParent.style.margin = "0";
          scriptParent.style.padding = "0";
        }

        require.config({
          paths: {
            Plotly: "$( plotly_url )",
            noUiSlider: "$( assets_url * "nouislider.min" )",
            wNumb: "$( assets_url * "wnumb.min" )"
          }
        });

        require(["Plotly", "noUiSlider", "wNumb"], function(Plotly, noUiSlider, wNumb){
          window.Plotly = Plotly;
          window.noUiSlider = noUiSlider;
          window.wNumb = wNumb;

          if ( typeof curCallback !== "undefined" ) {
            curCallback();
          }

          \$("head").append(\$(".js-plotly-script"));
          \$("head .js-plotly-script:not(:first)").remove();
          \$("body .js-plotly-script").remove();
        });
      }

      \$(document).ready(customPlotLoader);
      customPlotLoader();
    </script>
  """

  display(HTML(plotly_javascript))

end


text:
date = Dates.DateTime(2020, 2, 20, 11, 11, 0)
lon, lat = [FT(80.0), FT(20.0)]

sza, azi, d = instantaneous_zenith_angle(
    date,
    lon,
    lat,
    param_set,
    eot_correction=false,
    milankovitch=false)

F = insolation(sza, d, param_set)
@test typeof(sza) == FT
@test typeof(azi) == FT
@test typeof(d) == FT
@test typeof(F) == FT

sza, azi, d = instantaneous_zenith_angle(
    date,
    lon,
    lat,
    param_set,
    eot_correction=false,
    milankovitch=true)

F = insolation(sza, d, param_set)
@test typeof(sza) == FT
@test typeof(azi) == FT
@test typeof(d) == FT
@test typeof(F) == FT

sza, azi, d = instantaneous_zenith_angle(
    date,
    lon,
    lat,
    param_set,
    eot_correction=true,
    milankovitch=false)

F = insolation(sza, d, param_set)
@test typeof(sza) == FT
@test typeof(azi) == FT
@test typeof(d) == FT
@test typeof(F) == FT

sza, azi, d = instantaneous_zenith_angle(
    date,
    lon,
    lat,
    param_set,
    eot_correction=true,
    milankovitch=true)

F = insolation(sza, d, param_set)
@test typeof(sza) == FT
@test typeof(azi) == FT
@test typeof(d) == FT
@test typeof(F) == FT


text:
#push!(LOAD_PATH,joinpath("..","src"))

using Multilane
using MCTS
using POMDPs
using POMDPSimulators
using Random
using Test
using Cairo

#Set up problem configuration
nb_lanes = 4
pp = PhysicalParam(nb_lanes,lane_length=100.) #2.=>col_length=8
_discount = 1.
nb_cars=10

rmodel = NoCrashRewardModel()

dmodel = NoCrashIDMMOBILModel(nb_cars, pp)

mdp = NoCrashMDP{typeof(rmodel), typeof(dmodel.behaviors)}(dmodel, rmodel, _discount, true);

rng = MersenneTwister(5)

s = rand(rng, initialstate(mdp))
# @show s.cars[1]
# visualize(mdp,s,MLAction(0,0))

policy = Multilane.BehaviorPolicy(mdp, Multilane.NORMAL, false, rng)

sim = HistoryRecorder(rng=rng, max_steps=100) # initialize a random number generator

hist = POMDPs.simulate(sim, mdp, policy, s)

# check for crashes
for i in 1:length(state_hist(hist))-1
    if is_crash(mdp, state_hist(hist)[i], state_hist(hist)[i+1])
        println("Crash:")
        println("mdp = $mdp\n")
        println("s = $(state_hist(hist)[i])\n")
        println("a = $(action_hist(hist)[i])\n")
        println("Saving gif...")
        f = write_tmp_gif(mdp, hist)
        println("gif written to $f")
    end
    @test !is_crash(mdp, state_hist(hist)[i], state_hist(hist)[i+1])
end

a_hist = collect(action_hist(hist))
s_hist = collect(state_hist(hist))
for i in 1:length(s_hist)-1
    if is_crash(mdp, s_hist[i], s_hist[i+1])
        visualize(mdp, s_hist[i], a_hist[i], s_hist[i+1], two_frame_crash=true)
        # println(repr(mdp))
        # println(repr(s_hist[i]))
        println("Crash after step $i")
        println("Chosen Action: $(a_hist[i])")
        println("Available actions:")
        for a in actions(mdp, s_hist[i], actions(mdp))
            println(a)
        end
        println("Press Enter to continue.")
        readline(STDIN)
    end
end

behaviors = standard_uniform(correlation=0.75)
dmodel = NoCrashIDMMOBILModel(dmodel, behaviors)
pomdp = NoCrashPOMDP{typeof(rmodel), typeof(dmodel.behaviors)}(dmodel, rmodel, _discount, true);

rng = MersenneTwister(5)

s = rand(rng, initialstate(mdp))

policy = Multilane.BehaviorPolicy(pomdp, Multilane.NORMAL, false, rng)

sim = HistoryRecorder(rng=rng, max_steps=100) # initialize a random number generator

wup = WeightUpdateParams(smoothing=0.0, wrong_lane_factor=0.5)
up = BehaviorParticleUpdater(pomdp, 1000, 0.1, 0.1, wup, MersenneTwister(50000))

@time hist = POMDPs.simulate(sim, pomdp, policy, up, MLPhysicalState(s), s)
@show n_steps(hist)

# check for crashes
a_hist = collect(action_hist(hist))
s_hist = collect(state_hist(hist))
for i in 1:length(state_hist(hist))-1
    if is_crash(mdp, s_hist[i], s_hist[i+1])
        println("Crash:")
        println("mdp = $mdp\n")
        println("s = $(s_hist[i])\n")
        println("a = $(a_hist[i])\n")
        println("Saving gif...")
        f = write_tmp_gif(mdp, hist)
        println("gif written to $f")
    end
    @test !is_crash(mdp, s_hist[i], s_hist[i+1])
end


text:
#=
Test suite across dimensions and element types.
An aggregation of individual test cases (e.g. see test/2D/tiny.jl)

See below.
Loosely speaking:
    data types: color, grey, char
    dimensions: 2, 3

>cd <top>
>julia --project=test test/dimensions.jl
=#

# The tested code
include("../src/apps/inpaint.jl")

# the Julia testing framework
using ReferenceTests
using Test
using Profile
using TestImages    # for "testimage" of packaged, remote  test files ?
using FileIO    # For "load" local  files
using ColorTypes

#=
This is executable, not a definition
Each let;test_reference;end is a test of the testset
Use let;end for separate scope in each test (no spillover of bindings)
=#
@testset "Dimensions" begin

    # 1D text
    # Not working, index is Int64 instead of Cartesian{1}

    # 2D text
    # exactly 72 bytes long
    let
        inPath = "/home/bootch/git/resynthesizerJulia/test/data/in/tiny.txt"
        # !!! load will not load a text, use read()
        file = open(inPath)
        textStream = read(file, 72)
        # reshape to an array, a pseudo-poem
        data = reshape(textStream, 6, 12)

        mask = falses(size(data))
        mask[3:4, 2:3] .= true

        outPath = "data/out/tiny.txt"
        @test_reference outPath inpaint(data, mask)
    end

    # 2D color image
    let
        inPath = "/home/bootch/git/resynthesizerJulia/test/data/in/tinylighthouse.png"
        data = load(inPath)
        mask = falses(size(data))
        mask[3:4, 2:3] .= true

        outPath = "data/out/tinylighthouse.png"
        @test_reference outPath inpaint(data, mask)
    end

    # 3D gray (slices of an mri)
    let
        using TestImages
        fullImage = testimage("mri");
        print(summary(fullImage))

        # See JuliaImages. The image has named axes. We only want to convert the data.
        data = fullImage.data

        # mask
        mask = falses(size(data))
        # synthesize a small cube in the horizontal slices 5 and 6
        mask[3:4, 2:3, 5:6] .= true

        out = inpaint(data, mask)

        outPath = "data/out/mri"
        #=
        TODO this fails because test_reference fails to handle 3D???

        @test_reference outPath inpaint(data, mask)

        Consequently, we iteratively compare some slices of the image.
        =#

        # The first horizontal (superior axis) plane is at z=5
        slice = out[:,:,5]
        print(summary(slice))

        # compare one slice to the reference
        @test_reference outPath slice
    end

    # 3D color video (time lapsed 2D images)
    # TODO

    # 4D e.g. time lapsed 3D images, or time lapses 2D images with explicit channels?
    # copied from test/4D/multi.jl
    # "multi" is short name
    let
        fullImage = testimage("multi");

        println("fullImage is:")
        println(summary(fullImage))
        # 167×439×3×7

        #=
        See JuliaImages.
        The image has named axes and metadata??
        We only want to convert the data.
        =#
        using ImageMetadata

        # Extract: data(fullImage), not fullImage.data?? data() is a method of ImageMetadata
        imageWithAxes = data(fullImage)

        imageView = imageWithAxes.data
        #=
        imageView is still reshape(reinterpret(...)) which throws type errors.
        Its a view. So copy it.
        =#
        image = copy(imageView)

        # Create mask
        mask = falses(size(image))
        # ??? The third dimension is 3 channels.  Arbitrarily take them all.
        mask[3:4, 2:3, 1:3, 5:6] .= true

        # test
        out = inpaint(image, mask)

        outPath = "data/out/multi"
        #=
        test_reference fails to handle 3D, 4D ???
        Consequently, we iteratively compare some 2D slices of the image.
        =#

        #=
        I'm not sure what this slice represents.
        z=1 is the first of three channels?
        u=5 is the 5th of 7 time slots?
        Anyway, it must be 2D for test_reference.
        And the slice is through the mask, the synth region.
        =#
        slice = out[:,:,1,5]
        print(summary(slice))

        # compare one slice to the reference
        @test_reference outPath slice
    end
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using Weave

function weaveit(name)
    @info "Building $name"
    weave(name, out_path=joinpath("..","docs","src","$(splitext(name)[1]).md"), doctype="github")
end

let start_dir = pwd()
    cd(@__DIR__)
    outpath = normpath(joinpath(@__DIR__, "..", "docs", "src"))
    @show outpath
    if !isdirpath(outpath)
        mkpath(outpath)
    end

    weaveit.(("gettingstarted.jmd", "resistors.jmd"))
    if !isfile(joinpath(outpath,"resistor.png"))
        cp("resistor.png",joinpath(outpath,"resistor.png"))
    end
    cd(start_dir)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Statistics
X = rand(Float32, 10^7)

x = copy(X); y = quantile!(x, [0.02, 0.98])
x = copy(X); z = [quantile!(x, 0.02), quantile!(x, 0.98)]
@show y == z;

# %%
x = copy(X); @time quantile!(x, [0.02, 0.98])
x = copy(X); @time quantile!(x, [0.02, 0.98])
x = copy(X); @time quantile!(x, [0.02, 0.98]);

# %%
x = copy(X); @time [quantile!(x, 0.02), quantile!(x, 0.98)]
x = copy(X); @time [quantile!(x, 0.02), quantile!(x, 0.98)]
x = copy(X); @time [quantile!(x, 0.02), quantile!(x, 0.98)];

# %%
# `@time`'s added version of https://github.com/JuliaLang/Statistics.jl/blob/master/src/Statistics.jl#L949-L959

function myquantile!(v::AbstractVector, p::Union{AbstractArray, Tuple{Vararg{Real}}};
                   sorted::Bool=false, alpha::Real=1., beta::Real=alpha)
    if !isempty(p)
        minp, maxp = extrema(p)
        print("Statistics._quantilesort!(v, sorted, minp, maxp):             ")
        @time Statistics._quantilesort!(v, sorted, minp, maxp)
    end
    print("map(x->Statistics._quantile(v, x, alpha=alpha, beta=beta), p):")
    return @time map(x->Statistics._quantile(v, x, alpha=alpha, beta=beta), p)
end

myquantile!(v::AbstractVector, p::Real; sorted::Bool=false, alpha::Real=1., beta::Real=alpha) = begin
    print("Statistics._quantilesort!(v, sorted, p, p):        ")
    @time Statistics._quantilesort!(v, sorted, p, p)
    print("Statistics._quantile(v, p, alpha=alpha, beta=beta):")
    @time Statistics._quantile(v, p, alpha=alpha, beta=beta)
end

println("myquantile!(x, [0.02, 0.98]):")
x = copy(X); u = myquantile!(x, [0.02, 0.98])
println("\n[myquantile!(x, 0.02), myquantile!(x, 0.98)]:")
x = copy(X); v = [myquantile!(x, 0.02), myquantile!(x, 0.98)]
println(); @show u == v;

# %%
x = copy(X); @time y = quantile!(x, [0.499, 0.501])
x = copy(X); @time y = quantile!(x, [0.499, 0.501])
x = copy(X); @time y = quantile!(x, [0.499, 0.501]);

# %%
x = copy(X); @time z = [quantile!(x, 0.499), quantile!(x, 0.501)]
x = copy(X); @time z = [quantile!(x, 0.499), quantile!(x, 0.501)]
x = copy(X); @time z = [quantile!(x, 0.499), quantile!(x, 0.501)];

# %%
@which quantile!(x, [0.02, 0.98])

# %%
@which quantile!(x, 0.02)

# %%
@which sort!(x, 1, length(x), Base.Sort.PartialQuickSort(2*10^5:98*10^5), Base.Sort.Forward)

# %%

max_stars_repo_path: 0013/quantile! problem.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/quantile! problem.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/quantile! problem.jl

max_forks_repo_name: genkuroki/public


text:
using DelimitedFiles
using NCDatasets
using Interpolations

const fvars_name = ("δ", "V", "τ̂", "ψ")
const svars_name = ("δ", "V", "τ̂", "ψ")

"""
    init_fault_data(filename::String, var::String, nn::Integer)

creates a NetCDF file called `filename` in which fault time series of variable `var` is stored at along a total of `nn` nodes.

"""
function init_fault_data(filename::String, nn::Integer, depth::Array{Float64, 1})

    ds = NCDataset(filename, "c")

    defDim(ds, "time index", Inf)
    defDim(ds, "depth index", nn)

    defVar(ds, "time", Float64, ("time index",))
    defVar(ds, "depth", Float64, ("depth index",))
    defVar(ds, "maximum V", Float64, ("time index",))
    defVar(ds, "maximum v", Float64, ("time index",))

    for var in fvars_name
        defVar(ds, var, Float64, ("depth index", "time index"))
    end

    ds["depth"][:] .= depth

    close(ds)

end

"""
    init_station_data(filename::String, lendepths::Integer)

 creates a NetCDF file called `filename` in which station time series data is stored in `lendepths` total stations.

"""
function init_station_data(filename::String, stations::AbstractVector)

    ds = NCDataset(filename, "c")

    defDim(ds, "station index", length(stations))
    defDim(ds, "time index", Inf)

    defVar(ds, "time", Float64, ("time index",))
    defVar(ds, "maximum V", Float64, ("time index",))
    defVar(ds, "maxR", Float64, ("time index",))
    defVar(ds, "stations", Float64, ("station index",))
    defVar(ds, "δ", Float64, ("time index", "station index"))
    defVar(ds, "V", Float64, ("time index", "station index"))
    defVar(ds, "τ̂", Float64, ("time index", "station index"))
    defVar(ds, "ψ", Float64, ("time index", "station index"))

    ds["stations"][:] .= stations

    close(ds)

end

function init_volume_data(filename::String, x::Array{Float64, 1}, y::Array{Float64, 1})

    ds = NCDataset(filename, "c")

    defDim(ds, "time index", Inf)
    defDim(ds, "x index", length(x))
    defDim(ds, "y index", length(y))

    defVar(ds, "time", Float64, ("time index",))
    defVar(ds, "x", Float64, ("x index",))
    defVar(ds, "y", Float64, ("y index",))
    defVar(ds, "maximum V", Float64, ("time index",))
    defVar(ds, "u", Float64, ("x index", "y index", "time index"))
    defVar(ds, "v", Float64, ("x index", "y index", "time index"))
    defVar(ds, "σ", Float64, ("x index", "y index", "time index"))

    ds["x"][:] .= x
    ds["y"][:] .= y

    close(ds)

end

"""
    new_dir(new_dir::String, stations::AbstractVector, nn::Integer)

creates a new directory called `new_dir` to store data from `stations`, volume, and fault variables, for solutions with `nn` nodes per dimension.

"""
function new_dir(new_dir::String, input_file::String, stations::Array{Float64, 1}, depth::Array{Float64,1}, x::Array{Float64,1}, y::Array{Float64,1})

    if !isdir(new_dir)
        mkdir(new_dir)
    else
        error("new directory already exists.")
    end

    nn = length(depth)

    fault_name = string(new_dir, "fault.nc")
    stations_name = string(new_dir, "stations.nc")
    remote_name = string(new_dir, "remote.nc")
    volume_name = string(new_dir, "volume.nc")

    init_fault_data(fault_name, nn, depth)
    init_fault_data(remote_name, nn, depth)
    init_station_data(stations_name, stations)

    init_volume_data(volume_name, x, y)

    cp(input_file, string(new_dir, "input_file.dat"))
    write_depth_grid(string(new_dir, "depth_grid.dat"), depth)

    return fault_name, stations_name, remote_name, volume_name

end

"""
    write_out_fault_data(filenames::Tuple, vars::Tuple, t::Float64)

writes out `vars` fault varibles at time `t` to NetCDF `filenames`.

"""
function write_out_fault_data(filename::String, vars::Tuple, maxv::Float64, t::Float64)

    file = NCDataset(filename, "a")
    max_V = maximum(vars[2])

    t_ind = size(file["time"])[1] + 1
    file["time"][t_ind] = t
    file["maximum V"][t_ind] = max_V
    file["maximum v"][t_ind] = maxv
    for i in 1:length(vars)
        file[fvars_name[i]][:, t_ind] .= vars[i]
    end

    close(file)

end

"""
    write_out_stations(station_file::String, stations::AbstractVector, depth:: Array{Float64,1}, vars::Tuple)

writes out interpolated station data `vars` at `stations` using grid spacing `depth` at time `t` to netCDF file `station_file`.

"""
function write_out_stations(station_file::String, stations::Array{Float64,1}, depth::Array{Float64,1}, vars::Tuple, maxR::Float64, t::Float64)

    file = NCDataset(station_file, "a")
    t_ind = size(file["time"])[1] + 1
    file["time"][t_ind] = t
    file["maximum V"][t_ind] = maximum(vars[2])
    file["maxR"][t_ind] = maxR
    #file["maximum v"][t_ind] = maximum(maxv)
    for (i, var) in enumerate(vars)
        interp = interpolate((depth,), var, Gridded(Linear()))
        var_stations = interp(stations)
        file[svars_name[i]][t_ind, :] .= var_stations
    end

    close(file)

end

function write_out_volume(volume_file::String, volume_vars::Tuple, V::Array{Float64,1}, nn::Integer, t::Float64)

    # this is probably pretty inefficient.....
    u = reshape(volume_vars[1], (nn,nn))
    v = reshape(volume_vars[2], (nn,nn))

    file = NCDataset(volume_file, "a")
    t_ind = size(file["time"])[1] + 1
    file["time"][t_ind] = t
    file["maximum V"][t_ind] = maximum(V)
    file["u"][:, :, t_ind] .= u[1:2:end, 1:2:end]
    file["v"][:, :, t_ind] .= v[1:2:end, 1:2:end]

end

function write_depth_grid(filename, fc)

    io = open(filename, "w")

    for i in 2:length(fc)
        write(io, string(fc[i], ", ", fc[i] - fc[i-1], "\n"))
    end

    close(io)

end


text:
# struct HybridProblem{uType,tType,iip,P<:AbstractODEProblem{uType,tType,iip},C} <: AbstractODEProblem{uType,tType,iip}
#     problem::P
#     callbacks::C
#     function HybridProblem(prob::AbstractODEProblem{uType,tType,iip}, cb) where {uType,tType,iip}
#         new{uType, tType, iip, typeof(prob), typeof(cb)}(prob, cb)
#     end
# end

# function HybridProblem(cont_prob::AbstractODEProblem{Tu,Tt,iip}, d_f, d_u0, Δt) where {Tu,Tt,iip}
#     saved_values = SavedValues(typeof(Δt), typeof(d_u0))

#     cb_affect! = function ()
#     end

#     p_cb = PeriodicCallback(d_f, Δt; initial_affect=false)
#     s_cb = SavingCallback(d_f, saved_values, saveat=Δt)
#     return HybridProblem(cont_prob, cb)
# end
# function HybridProblem{iip}(cont_f, disc_f, cont_u0, disc_u0, Δt, args...; kwargs...) where {iip}
#     prob = ODEProblem{iip}(cont_f, cont_u0, args...; kwargs...)
#     return HybridProblem(prob, disc_f, disc_u0, Δt)
# end
# HybridProblem(c_f, args...; kwargs...) = HybridProblem{isinplace(c_f, 5)}(c_f, args...; kwargs...)

# struct HybridFunction{iip, recompile, FC<:AbstractODEFunction{iip}, FD} <: AbstractODEFunction{iip}
#     f_cont::FC
#     f_disc::FD
# end

# hybrid_call(f!, du, cont_u, (disc_u, p), t) = f!(du, cont_u, disc_u, p, t)
# hybrid_call(f, cont_u, (disc_u, p), t) = f(cont_u, disc_u, p, t)

function DiscreteVariableCallback(f, saved_values, Δt; initial_affect=false, kwargs...)
    callback_fun! = function (integrator)
        @unpack p, t = integrator
        u_discrete = saved_values.u[end]
        push!(saved_values.t, t)
        push!(saved_values.u, f(integrator.u, u_discrete, p, t))
        return nothing
    end

    return PeriodicCallback(
        callback_fun!,
        Δt;
        initial_affect=initial_affect,
        save_positions=(false, true),
        kwargs...
    )
end

"""
    ContinuousDiscreteProblem
"""
struct ContinuousDiscreteProblem

end

"""
    ContinuousDiscreteSolution(cont_sol::AbstractTimeseriesSolution, disc_sol::DiffEqArray)
"""
struct ContinuousDiscreteSolution{C<:AbstractTimeseriesSolution, D<:DiffEqArray}
    continuous::C
    discrete::D
end

(sol::ContinuousDiscreteSolution)(t) = (continuous=sol.continuous(t), discrete=zoh_interp(sol.discrete, t))

function Base.getindex(sol::ContinuousDiscreteSolution, k_cont)
    cont = sol.continuous
    return (continuous=cont[k_cont], discrete=zoh_interp(sol.discrete, cont.t[k_cont]))
end

zoh_interp(dea::DiffEqArray, t) = dea.u[max(1,searchsortedlast(dea.t, t))]

# function myfun!(dx, x, x_disc, p, t)

# end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SubnetPropertiesFormat <: SwaggerModel
    addressPrefix::Any # spec type: Union{ Nothing, String } # spec name: addressPrefix
    networkSecurityGroup::Any # spec type: Union{ Nothing, NetworkSecurityGroup } # spec name: networkSecurityGroup
    routeTable::Any # spec type: Union{ Nothing, RouteTable } # spec name: routeTable
    ipConfigurations::Any # spec type: Union{ Nothing, Vector{IPConfiguration} } # spec name: ipConfigurations
    resourceNavigationLinks::Any # spec type: Union{ Nothing, Vector{ResourceNavigationLink} } # spec name: resourceNavigationLinks
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function SubnetPropertiesFormat(;addressPrefix=nothing, networkSecurityGroup=nothing, routeTable=nothing, ipConfigurations=nothing, resourceNavigationLinks=nothing, provisioningState=nothing)
        o = new()
        validate_property(SubnetPropertiesFormat, Symbol("addressPrefix"), addressPrefix)
        setfield!(o, Symbol("addressPrefix"), addressPrefix)
        validate_property(SubnetPropertiesFormat, Symbol("networkSecurityGroup"), networkSecurityGroup)
        setfield!(o, Symbol("networkSecurityGroup"), networkSecurityGroup)
        validate_property(SubnetPropertiesFormat, Symbol("routeTable"), routeTable)
        setfield!(o, Symbol("routeTable"), routeTable)
        validate_property(SubnetPropertiesFormat, Symbol("ipConfigurations"), ipConfigurations)
        setfield!(o, Symbol("ipConfigurations"), ipConfigurations)
        validate_property(SubnetPropertiesFormat, Symbol("resourceNavigationLinks"), resourceNavigationLinks)
        setfield!(o, Symbol("resourceNavigationLinks"), resourceNavigationLinks)
        validate_property(SubnetPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type SubnetPropertiesFormat

const _property_map_SubnetPropertiesFormat = Dict{Symbol,Symbol}(Symbol("addressPrefix")=>Symbol("addressPrefix"), Symbol("networkSecurityGroup")=>Symbol("networkSecurityGroup"), Symbol("routeTable")=>Symbol("routeTable"), Symbol("ipConfigurations")=>Symbol("ipConfigurations"), Symbol("resourceNavigationLinks")=>Symbol("resourceNavigationLinks"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_SubnetPropertiesFormat = Dict{Symbol,String}(Symbol("addressPrefix")=>"String", Symbol("networkSecurityGroup")=>"NetworkSecurityGroup", Symbol("routeTable")=>"RouteTable", Symbol("ipConfigurations")=>"Vector{IPConfiguration}", Symbol("resourceNavigationLinks")=>"Vector{ResourceNavigationLink}", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ SubnetPropertiesFormat }) = collect(keys(_property_map_SubnetPropertiesFormat))
Swagger.property_type(::Type{ SubnetPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_SubnetPropertiesFormat[name]))}
Swagger.field_name(::Type{ SubnetPropertiesFormat }, property_name::Symbol) =  _property_map_SubnetPropertiesFormat[property_name]

function check_required(o::SubnetPropertiesFormat)
    true
end

function validate_property(::Type{ SubnetPropertiesFormat }, name::Symbol, val)
end


text:
# series of deconvolution tools

export selectFactorType
export approxDeconv, deconvSolveKey
export approxDeconvBelief

## Initial version of selecting the dimension of a factor -- will be consolidated with existing infrastructure later

"""
    $SIGNATURES

Inverse solve of predicted noise value and returns tuple of (newly calculated-predicted, and known measurements) values.

Notes
- Only works for first value in measurement::Tuple at this stage.
- "measured" is used as starting point for the "calculated-predicted" values solve.
- Not all factor evaluation cases are support yet.
- NOTE only works on `.threadid()==1` at present, see #1094
- This function is still part of the initial implementation and needs a lot of generalization improvements.

DevNotes
- TODO Test for various cases with multiple variables.
- TODO make multithread-safe, and able, see #1094
- TODO Test for cases with `nullhypo`
- FIXME FactorMetadata object for all use-cases, not just empty object.
- TODO resolve #1096 (multihypo)
  - TODO Test cases for `multihypo`.
- TODO figure out if there is a way to consolidate with evalFactor and approxConv?
  - basically how to do deconv for just one sample with unique values (wrt TAF)
- TODO N should not be hardcoded to 100

Related

[`approxDeconv`](@ref), [`_solveCCWNumeric!`](@ref)
"""
function approxDeconv(fcto::DFGFactor,
                      ccw::CommonConvWrapper = _getCCW(fcto);
                      N::Int=100,
                      measurement::AbstractVector=sampleFactor(ccw, N),
                      retries::Int=3  )
  #
  # but what if this is a partial factor -- is that important for general cases in deconv?
  _setCCWDecisionDimsConv!(ccw)

  # FIXME This does not incorporate multihypo??
  varsyms = getVariableOrder(fcto)
  # vars = getPoints.(getBelief.(dfg, varsyms, solveKey) )
  fcttype = getFactorType(fcto)

  # get measurement dimension
  zDim = _getZDim(fcto)
  # TODO consider using ccw.cpt[thrid].res # likely needs resizing
  res_ = zeros(zDim)
  # TODO, consolidate fmd with getSample/sampleFactor and _buildLambda
  fctSmpls = deepcopy(measurement)
  fmd = _getFMdThread(ccw)

  # TODO assuming vector on only first container in measurement::Tuple
  makeTarget = (i) -> measurement[i] # TODO does not support copy-primitive types like Float64, only Ref()
  # makeTarget = (i) -> view(measurement[1][i],:)
  # makeTarget = (i) -> view(measurement[1], :, i)

  # NOTE
  # build a lambda that incorporates the multihypo selections
  # set these first
  # ccw.cpt[].activehypo / .p / .params  # params should already be set from construction
  hyporecipe = _prepareHypoRecipe!(nothing, N, 0, length(varsyms))
  # Juila 1.7 allows destructure assign `(;a,b) = namedtype`
  # certainidx, allelements, activehypo, mhidx =
  # only doing the current active hypo
  @assert hyporecipe.activehypo[2][1] == 1 "deconv was expecting hypothesis nr == (1, 1:d)"

  islen1 = zDim == 1

  for idx in 1:N
    # towards each particle in their own thread (not 100% ready yet, factors should be separate memory)
    thrid = Threads.threadid()
    cpt_ = ccw.cpt[thrid]
    targeti_ = makeTarget(idx)

    # TODO must first resolve hypothesis selection before unrolling them -- deferred #1096
    cpt_.activehypo = hyporecipe.activehypo[2][2]

    onehypo!, _ = _buildCalcFactorLambdaSample( ccw,
                                                idx,
                                                cpt_,
                                                targeti_,
                                                measurement  )
    #

    # lambda with which to find best measurement values
    hypoObj = (tgt) -> (targeti_.=tgt; onehypo!() )

    # find solution via SubArray view pointing to original memory location
    if fcttype isa AbstractManifoldMinimize
      sfidx = ccw.varidx
      targeti_ .= _solveLambdaNumericMeas(fcttype, hypoObj, res_, measurement[idx], ccw.vartypes[sfidx](), islen1)
    else
      targeti_ .= _solveLambdaNumeric(fcttype, hypoObj, res_, measurement[idx], islen1)
    end

  end

  # return (deconv-prediction-result, independent-measurement)
  # r_meas = map(m->m[1], measurement)
  # r_fctSmpls = map(m->m[1], fctSmpls)
  return measurement, fctSmpls
end

"""
    $SIGNATURES

Generalized deconvolution to find the predicted measurement values of the factor `fctsym` in `dfg`.
Inverse solve of predicted noise value and returns tuple of (newly predicted, and known "measured" noise) values.

Notes
- Opposite operation contained in `approxConv`.
- For more notes see [`solveFactorMeasurements`](@ref).

Related

[`approxConv`](@ref), `deconvSolveKey`
"""
function approxDeconv(dfg::AbstractDFG,
                      fctsym::Symbol,
                      solveKey::Symbol=:default;
                      retries::Int=3  )
  #

  # which factor
  fct = getFactor(dfg, fctsym)
  pts = getPoints(getBelief(dfg,getVariableOrder(fct)[1], solveKey))
  N = length(pts)
  pts = approxDeconv(fct, N=N, retries=retries )
  return pts
end

function approxDeconv(dfg::AbstractDFG,
                      fctlbl::Symbol,
                      factorType::AbstractRelative,
                      solveKey::Symbol=:default;
                      tfg::AbstractDFG=initfg(),
                      retries::Int=3 )
  #

  # build a local temporary graph copy containing the same values but user requested factor type.
  fct = getFactor(dfg, fctlbl)
  fctT = getFactorType(fct)
  lbls = getVariableOrder(fct)
  for lb in lbls
    exists(tfg, lb) ? nothing : addVariable!(tfg, lb, getVariableType(dfg, lb))
    initManual!(tfg, lb, getBelief(dfg, lb, solveKey))
  end

  # add factor type requested by user
  f_ = addFactor!(tfg, lbls, factorType, graphinit=false)

  # peform the deconvolution operation on the temporary graph with user desired factor instead.
  approxDeconv(tfg, getLabel(f_); retries=retries)
end

# try default constructor
approxDeconv( dfg::AbstractDFG,
              fctlbl::Symbol,
              factorType::Type{<:AbstractRelative},
              w...;
              kw... ) = approxDeconv(dfg, fctlbl, factorType(), w...; kw...)
#

approxDeconvBelief(dfg::AbstractDFG, lb::Symbol, w...;kw...) = manikde!(approxDeconv(dfg, lb, w...;kw...)[1], getManifold(getFactorType(dfg, lb)))

"""
    $SIGNATURES

Calculate the relative difference between two variables and across solveKeys.

Example

```julia
fg = generateCanonicalFG_Kaess()
solveTree!(fg, storeOld=true)
# calculate the relative motion induced by the solver from init to solve.
pts = deconvSolveKey(fg, :x1, :default, :x1, :graphinit)
```

Notes
- Can pass user `tfg::AbstractDFG` for better in-place operations.

DevNotes
- TODO use dfg, rather than building new tfg internally.

Related

[`approxDeconv`](@ref), [`mmd`](@ref)
"""
function deconvSolveKey(dfg::AbstractDFG,
                        refSym::Symbol,
                        refKey::Symbol,
                        tstSym::Symbol,
                        tstKey::Symbol;
                        tfg = initfg()  )
  #
  # create a new temporary factor graph for calculations

  # add the first "reference" variable
  Xref = getBelief(dfg, refSym, refKey)
  refSym_ = Symbol(refSym, "_ref")
  refVarType = getVariableType(dfg, refSym)
  if !exists(tfg, refSym_)
    addVariable!(tfg, refSym_, refVarType)
  end
  initManual!(tfg, refSym_, Xref)

  # add the second "test" variable
  tstVarType = getVariableType(dfg, tstSym)
  Xtst = getBelief(dfg, tstSym, tstKey)
  tstSym_ = Symbol(tstSym, "_tst")
  if !exists(tfg, tstSym_)
    addVariable!(tfg, tstSym_, tstVarType)
  end
  initManual!(tfg, tstSym_, Xtst)

  # add the new dummy factor with default manifold for computations
  fctType = selectFactorType(refVarType, tstVarType)
  nf = addFactor!(tfg, [refSym_; tstSym_], fctType())

  # TODO connect from dfg all other data that might form part of FactorMetadata in tfg
  pts = approxDeconv(tfg, nf.label)

  # assuming tfg was passed in by the user
  deleteFactor!(tfg, nf.label)

  # return result
  return pts, fctType
end

#


text:
"""
yeeDielectric = cellToYeeDielectric(cell_σ, cell_ϵr, x_axis, y_axis, z_axis).\n

Input the conductivity and permittivity on a cell grid (3D Array).\n
And input the node (corners of the cells) coordinates that make of the grid (x,y,z 1D Vectors).\n

The function calculates what the dielectric properties are on the edges of each cell.\n

returns yeeDielectric = (;σˣ, ϵˣ, σʸ, ϵʸ, σᶻ, ϵᶻ).\n
"""
function cell_to_yee_dielectric(cell_σ, cell_ϵr, x, y, z)
    #get voxel lengths
    xLengths = diff(vec(x))
    yLengths = diff(vec(y))
    zLengths = diff(vec(z))

    #get surface areas to average the dielectric over
    x_surface = yLengths*zLengths'
    y_surface = xLengths*zLengths'
    z_surface = xLengths*yLengths'

    #translate the voxel dielectric to the yee cell grid
    σˣ = _x_dielectric(cell_σ, x_surface)
    ϵˣ = _x_dielectric(cell_ϵr, x_surface)
    σʸ = _y_dielectric(cell_σ, y_surface)
    ϵʸ = _y_dielectric(cell_ϵr, y_surface)
    σᶻ = _z_dielectric(cell_σ, z_surface)
    ϵᶻ = _z_dielectric(cell_ϵr, z_surface)

    return (;σˣ, ϵˣ, σʸ, ϵʸ, σᶻ, ϵᶻ)
end

function _x_dielectric(A, x_surface)
    typeToUse = eltype(A)
    dimensions = size(A)

    Aˣ = zeros(typeToUse, dimensions[1], dimensions[2]+1, dimensions[3]+1)
    #average the center of the grid
    Aˣ[:,2:end-1,2:end-1]   = permutedims((permutedims(A[:,1:end-1,1:end-1], [2,3,1]).*x_surface[1:end-1,1:end-1]+
                                            permutedims(A[:,1:end-1,2:end], [2,3,1]).*x_surface[1:end-1,2:end]+
                                            permutedims(A[:,2:end,1:end-1], [2,3,1]).*x_surface[2:end,1:end-1]+
                                            permutedims(A[:,2:end,2:end], [2,3,1]).*x_surface[2:end,2:end])./
                                            (x_surface[1:end-1,1:end-1]+x_surface[1:end-1,2:end]+
                                            x_surface[2:end,1:end-1]+x_surface[2:end,2:end]), [3,1,2])
    #average the sides
    Aˣ[:,1,2:end-1]      = permutedims((permutedims(A[:,1,1:end-1], [2,1]).*x_surface[1,1:end-1]+
                                        permutedims(A[:,1,2:end], [2,1]).*x_surface[1,2:end])./
                                        (x_surface[1,1:end-1]+x_surface[1,2:end]), [2,1])

    Aˣ[:,end,2:end-1]    = permutedims((permutedims(A[:,end,1:end-1], [2,1]).*x_surface[end,1:end-1]+
                                        permutedims(A[:,end,2:end], [2,1]).*x_surface[end,2:end])./
                                        (x_surface[end,1:end-1]+x_surface[end,2:end]), [2,1])

    Aˣ[:,2:end-1,1]      = permutedims((permutedims(A[:,1:end-1,1], [2,1]).*x_surface[1:end-1,1]+
                                        permutedims(A[:,2:end,1], [2,1]).*x_surface[2:end,1])./
                                        (x_surface[1:end-1,1]+x_surface[2:end,1]), [2,1])

    Aˣ[:,2:end-1,end]    = permutedims((permutedims(A[:,1:end-1,end], [2,1]).*x_surface[1:end-1,end]+
                                        permutedims(A[:,2:end,end], [2,1]).*x_surface[2:end,end])./
                                        (x_surface[1:end-1,end]+x_surface[2:end,end]), [2,1])
    #set the edges of the grid
    Aˣ[:,1,1]     = A[:,1,1]
    Aˣ[:,end,1]   = A[:,end,1]
    Aˣ[:,1,end]   = A[:,1,end]
    Aˣ[:,end,end] = A[:,end,end]
    return Aˣ
end

function _y_dielectric(A, y_surface)
    typeToUse = eltype(A)
    dimensions = size(A)

    Aʸ = zeros(typeToUse, dimensions[1]+1, dimensions[2], dimensions[3]+1)
    #average the center of the grid
    Aʸ[2:end-1,:,2:end-1] = permutedims((permutedims(A[1:end-1,:,1:end-1], [1,3,2]).*y_surface[1:end-1,1:end-1]+
                                         permutedims(A[1:end-1,:,2:end], [1,3,2]).*y_surface[1:end-1,2:end]+
                                         permutedims(A[2:end,:,1:end-1], [1,3,2]).*y_surface[2:end,1:end-1]+
                                         permutedims(A[2:end,:,2:end], [1,3,2]).*y_surface[2:end,2:end])./
                                         (y_surface[1:end-1,1:end-1]+y_surface[1:end-1,2:end]+
                                         y_surface[2:end,1:end-1]+y_surface[2:end,2:end]), [1,3,2])
    #average the sides
    Aʸ[1,:,2:end-1] = permutedims((permutedims(A[1,:,1:end-1], [2,1]).*y_surface[1,1:end-1]+
                                   permutedims(A[1,:,2:end], [2,1]).*y_surface[1,2:end])./
                                   (y_surface[1,1:end-1]+y_surface[1,2:end]), [2,1])

    Aʸ[end,:,2:end-1] = permutedims((permutedims(A[end,:,1:end-1], [2,1]).*y_surface[end,1:end-1]+
                                     permutedims(A[end,:,2:end], [2,1]).*y_surface[end,2:end])./
                                     (y_surface[end,1:end-1]+y_surface[end,2:end]), [2,1])

    Aʸ[2:end-1,:,1] = (A[1:end-1,:,1].*y_surface[1:end-1,1]+
                       A[2:end,:,1].*y_surface[2:end,1])./
                       (y_surface[1:end-1,1]+y_surface[2:end,1])

    Aʸ[2:end-1,:,end] = (A[1:end-1,:,end].*y_surface[1:end-1,end]+
                         A[2:end,:,end].*y_surface[2:end,end])./
                         (y_surface[1:end-1,end]+y_surface[2:end,end])
    #set the edges of the grid
    Aʸ[1,:,1]     = A[1,:,1]
    Aʸ[end,:,1]   = A[end,:,1]
    Aʸ[1,:,end]   = A[1,:,end]
    Aʸ[end,:,end] = A[end,:,end]

    return Aʸ
end

function _z_dielectric(A, z_surface)
    typeToUse = eltype(A)
    dimensions = size(A)

    Aᶻ = zeros(typeToUse, dimensions[1]+1, dimensions[2]+1, dimensions[3])
    #average the center of the grid
    Aᶻ[2:end-1,2:end-1,:] = (A[1:end-1,1:end-1,:].*z_surface[1:end-1,1:end-1]+
                             A[1:end-1,2:end,:].*z_surface[1:end-1,2:end]+
                             A[2:end,1:end-1,:].*z_surface[2:end,1:end-1]+
                             A[2:end,2:end,:].*z_surface[2:end,2:end])./
                             (z_surface[1:end-1,1:end-1]+z_surface[2:end,2:end]+
                              z_surface[1:end-1,2:end]+z_surface[2:end,1:end-1])
    #side squares
    Aᶻ[1,2:end-1,:] = (A[1,1:end-1,:].*z_surface[1,1:end-1]+
                       A[1,2:end,:].*z_surface[1,2:end])./
                       (z_surface[1,1:end-1] + z_surface[1,2:end])

    Aᶻ[end,2:end-1,:] = (A[end,1:end-1,:].*z_surface[end,1:end-1]+
                         A[end,2:end,:].*z_surface[end,2:end])./
                         (z_surface[end,1:end-1] + z_surface[end,2:end])

    Aᶻ[2:end-1,1,:] = (A[1:end-1,1,:].*z_surface[1:end-1,1]+
                       A[2:end,1,:].*z_surface[2:end,1])./
                       (z_surface[1:end-1,1] + z_surface[2:end,1])

    Aᶻ[2:end-1,end,:] = (A[1:end-1,end,:].*z_surface[1:end-1,end]+
                         A[2:end,end,:].*z_surface[2:end,end])./
                         (z_surface[1:end-1,end] + z_surface[2:end,end])

    #edges
    Aᶻ[1,1,:]     = A[1,1,:]
    Aᶻ[end,1,:]   = A[end,1,:]
    Aᶻ[1,end,:]   = A[1,end,:]
    Aᶻ[end,end,:] = A[end,end,:]

    return Aᶻ
end


text:
register(DataDep(
    "pudl-release-v1.0.0",
    """
    $VERSION_DESCRIPTION_V1_0_0

    # Included Data Packages
    This release consists of three tabular data packages, conforming to the standards published by Frictionless Data and the Open Knowledge Foundation. The data are stored in CSV files (some of which are compressed using gzip), and the associated metadata is stored as JSON. These tabular data can be used to populate a relational database.

    ## pudl-eia860-eia923
    $(DATASET_DESCRIPTION_V1_0_0["pudl-eia860-eia923"])
    ## pudl-eia860-eia923-epacems
    $(DATASET_DESCRIPTION_V1_0_0["pudl-eia860-eia923-epacems"])
    ## pudl-ferc1
    $(DATASET_DESCRIPTION_V1_0_0["pudl-ferc1"])

    # Using the Data
    The data packages are just CSVs (data) and JSON (metadata) files. They can be used with a variety of tools on many platforms.
    However, the data is organized primarily with the idea that it will be loaded into a relational database, and the PUDL Python package that was used to generate this data release can facilitate that process.
    Once the data is loaded into a database, you can access that DB however you like.
    """,
    [
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/archived-environment.yml",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/data-release.log",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/data-release-settings.yml",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/data-release.sh",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/load-pudl.sh",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/pudl-eia860-eia923-epacems.tgz",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/pudl-eia860-eia923.tgz",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/pudl-ferc1.tgz",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/pudl-input-data.tgz",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/README.md",
        "https://zenodo.org/api/files/ff3b621f-6d31-40a8-b578-e87d7041f27a/reproduce-data-release.sh",
    ],
))


text:
using StrainRecon
using Base.Test

 @testset "findw" begin

n    = 2
R     = StrainRecon.sdiag(5*ones(2))
M     = BitArray(eye(2))
Data  = ones(2)

paramGauss = getStrainReconParam(Data,n,R=R,pwr=2.0)
#paramLaplace = getStrainReconParam(Data,n,R=R,pwr=1.0)
#param = paramLaplace
for param=[paramGauss]
	@testset "findw pwr=$(param.pwr)" begin
		a     = findw(param,M)

		@test norm(a-[.5;.5])<1e-10

		param.d = ones(2)
		a    = findw(param,M)
		@test norm(a-[.5;.5])<1e-10

		param.d = [190;0]
		a    = findw(param,M)
		@test norm(a-[1;0])<1e-10

		param.d = [0;20]
		a    = findw(param,M)
		@test norm(a-[.5;.5])<1e-10

		M    = BitArray([1.0 0; 0 0])
		param.d = ones(2)
		a    = findw(param,M)
		@test norm(a-[1;0])<1e-10
	end
end

    n    = 2
    R     = StrainRecon.sdiag(5*ones(2))
    M     = BitArray([1 0; 1 0])
    Data  = ones(2)
    paramLaplace = getStrainReconParam(Data,n,R=R,pwr=1.0)

for param=[paramLaplace]
        @testset "findw pwr=$(param.pwr)" begin
            a     = findw(param,M)
            # println("a = $a")
            @test norm(a-[1;0])<1e-10

            param.d = zeros(2)
            a    = findw(param,M)
            # println("a = $a")
            @test norm(a-[0.5;0.5])<1e-10

            M    = BitArray([1.0 0; 0 0])
            param.d = ones(2)
            a    = findw(param,M)
            # println("a = $a")
            @test norm(a-[1;0])<1e-10

            m = 24
            for n = 2:5
                param.R = StrainRecon.sdiag(5*ones(m))
                param.n = n
                # srand(s)
                M = rand([0,1], m, n)
                atrue = randomw(n)
                param.d = M*atrue
                M = BitArray(M)
                a    = findw(param,M)
                # println("a = $a , atrue = $atrue")
                @test norm(a-atrue)<1e-10
            end

        end
    end
end