text:
# Attempt at include dmstag functions
const CDMStag = Ptr{Cvoid}
const CDMStagType = Cstring

mutable struct DMStag{PetscLib} <: AbstractDM{PetscLib}
    ptr::CDMStag
    opts::Options{PetscLib}

    DMStag{PetscLib}(ptr, opts = Options(PetscLib)) where {PetscLib} =
    new{PetscLib}(ptr, opts)
end

"""
    empty(dm::DMStag)

return an uninitialized `DMStag` struct.
"""
Base.empty(::DMStag{PetscLib}) where {PetscLib} = DMStag{PetscLib}(C_NULL)

mutable struct DMSTAGSTENCIL{PetscInt}
    loc::DMStagStencilLocation
    i::PetscInt
    j::PetscInt
    k::PetscInt
    c::PetscInt
end

const DMStagStencil     = DMSTAGSTENCIL

# allows us to pass XXMat objects directly into CMat ccall signatures
Base.cconvert(::Type{CDMStag}, obj::DMStag) = obj.ptr
# allows us to pass XXMat objects directly into Ptr{CMat} ccall signatures
Base.unsafe_convert(::Type{Ptr{CDMStag}}, obj::DMStag) =
    convert(Ptr{CDMStag}, pointer_from_objref(obj))

#Base.eltype(::DMStag{T,PetscInt}) where {T,PetscInt} = T,PetscInt

"""
    dm = DMStagCreate1d(::PetscLib,
        comm::MPI.Comm,
        bndx::DMBoundaryType,
        M,
        dofVertex,
        dofCenter,
        stencilType::DMStagStencilType=DMSTAG_STENCIL_BOX,
        stencilWidth=2,
        lx=C_NULL;
        dmsetfromoptions=true,
        dmsetup=true,
        options...
        )

Creates a 1D DMStag object.
        ::PetscLib      -   PETSc library,
        comm            -   MPI communicator
        bndx            -   boundary type: DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, or DM_BOUNDARY_GHOSTED.
        M               -   global number of grid points
        dofVertex       -   [=1] number of degrees of freedom per vertex/point/node/0-cell
        dofCenter       -   [=1] number of degrees of freedom per element/edge/1-cell
        stencilType     -   ghost/halo region type: DMSTAG_STENCIL_BOX or DMSTAG_STENCIL_NONE
        stencilWidth    -   width, in elements, of halo/ghost region
        lx              -   [Optional] Vector of local sizes, of length equal to the comm size, summing to M
        options...      -   [Optional] keyword arguments (see PETSc webpage), specifiable as stag_grid_x=100, etc.

Creates a 1-D distributed staggered array with the options specified using keyword
arguments.

If keyword argument `dmsetfromoptions == true` then `setfromoptions!` called.
If keyword argument `dmsetup == true` then `setup!` is called.

# External Links
$(_doc_external("DMSTAG/DMStagCreate1d"))

"""
function DMStagCreate1d end

@for_petsc function DMStagCreate1d(
    ::$UnionPetscLib,
    comm::MPI.Comm,
    bndx::DMBoundaryType,
    M,
    dofVertex=1,
    dofCenter=1,
    stencilType::DMStagStencilType=DMSTAG_STENCIL_BOX,
    stencilWidth=2,
    lx=C_NULL;
    dmsetfromoptions = true,
    dmsetup = true,
    options...,
)

    if isempty(lx); lx = C_NULL; end
    opts = Options($petsclib, options...)
    dm   = DMStag{$PetscLib}(C_NULL, opts)   # retrieve options
    with(dm.opts) do

        @chk ccall(
                (:DMStagCreate1d, $petsc_library),
                PetscErrorCode,
                (
                    MPI.MPI_Comm,
                    DMBoundaryType,
                    $PetscInt,
                    $PetscInt,
                    $PetscInt,
                    DMStagStencilType,
                    $PetscInt,
                    Ptr{$PetscInt},
                    Ptr{CDMStag}
                ),
                comm,
                bndx,
                M,
                dofVertex,
                dofCenter,
                stencilType,
                stencilWidth,
                lx,
                dm
                )
    end
    dmsetfromoptions && setfromoptions!(dm)
    dmsetup && setup!(dm)

    if comm == MPI.COMM_SELF
        finalizer(destroy, dm)
    end

    return dm
end

"""
    dm = DMStagCreate2d(
        ::PetscLib,
        comm::MPI.Comm,
        bndx::DMBoundaryType,
        bndy::DMBoundaryType,
        M, N,
        m=C_NULL, n=C_NULL,
        dofVertex=1,
        dofEdge=1,
        dofElement=1,
        stencilType::DMStagStencilType=DMSTAG_STENCIL_BOX,
        stencilWidth=2,
        lx=C_NULL, ly=C_NULL;
        dmsetfromoptions=true,
        dmsetup=true,
        options...
        )

Creates a 2D DMStag object.

If keyword argument `dmsetfromoptions == true` then `setfromoptions!` called.
If keyword argument `dmsetup == true` then `setup!` is called.

# External Links
$(_doc_external("DMSTAG/DMStagCreate2d"))

"""
function DMStagCreate2d end

@for_petsc function DMStagCreate2d(
    ::$UnionPetscLib,
    comm::MPI.Comm,
    bndx::DMBoundaryType,
    bndy::DMBoundaryType,
    M, N,
    m=C_NULL, n=C_NULL,
    dofVertex=1,
    dofEdge=1,
    dofElement=1,
    stencilType::DMStagStencilType=DMSTAG_STENCIL_BOX,
    stencilWidth=2,
    lx=C_NULL, ly=C_NULL;
    dmsetfromoptions=true,
    dmsetup=true,
    options...,
)

    if isempty(lx); lx = C_NULL; end
    if isempty(ly); ly = C_NULL; end
    opts = Options($petsclib, options...)

    dm = DMStag{$PetscLib}(C_NULL, opts)

    with(dm.opts) do
        @chk ccall(
            (:DMStagCreate2d, $petsc_library),
            PetscErrorCode,
                (
                    MPI.MPI_Comm,
                    DMBoundaryType, DMBoundaryType,
                    $PetscInt, $PetscInt,
                    $PetscInt, $PetscInt,
                    $PetscInt,
                    $PetscInt,
                    $PetscInt,
                    DMStagStencilType,
                    $PetscInt,
                    Ptr{$PetscInt},  Ptr{$PetscInt},
                    Ptr{CDMStag}
                ),
                comm,
                bndx, bndy,
                M, N,
                m, n ,
                dofVertex ,
                dofEdge ,
                dofElement ,
                stencilType ,
                stencilWidth ,
                lx ,ly ,
                dm
            )
    end
    dmsetfromoptions && setfromoptions!(dm)
    dmsetup && setup!(dm)

    if comm == MPI.COMM_SELF
        finalizer(destroy, dm)
    end

    return dm
end

"""

    dm = DMStagCreate3d(
        ::PetscLib,
        comm::MPI.Comm,
        bndx::DMBoundaryType, bndy::DMBoundaryType, bndz::DMBoundaryType,
        M, N, P,
        m, n, p,
        dofVertex,
        dofEdge,
        dofFace,
        dofElement,
        stencilType::DMStagStencilType=DMSTAG_STENCIL_BOX,
        stencilWidth,
        lx, ly, lz;
        dmsetfromoptions=true,
        dmsetup=true,
        options...
        )

Creates a 3D DMStag object.

If keyword argument `dmsetfromoptions == true` then `setfromoptions!` called.
If keyword argument `dmsetup == true` then `setup!` is called.

# External Links
$(_doc_external("DMSTAG/DMStagCreate3d"))
"""
function DMStagCreate3d end

@for_petsc function DMStagCreate3d(
    ::$UnionPetscLib,
    comm::MPI.Comm,
    bndx::DMBoundaryType, bndy::DMBoundaryType, bndz::DMBoundaryType,
    M, N, P,
    m=C_NULL, n=C_NULL, p=C_NULL,
    dofVertex=1,
    dofEdge=1,
    dofFace=1,
    dofElement=1,
    stencilType::DMStagStencilType=DMSTAG_STENCIL_BOX,
    stencilWidth=2,
    lx=C_NULL, ly=C_NULL, lz=C_NULL;
    dmsetfromoptions=true,
    dmsetup=true,
    options...,
    )

    if isempty(lx); lx = C_NULL; end
    if isempty(ly); ly = C_NULL; end
    if isempty(lz); lz = C_NULL; end
    opts = Options($petsclib, options...)

    dm = DMStag{$PetscLib}(C_NULL, opts)

    with(dm.opts) do
        @chk ccall((:DMStagCreate3d, $petsc_library), PetscErrorCode,
                (
                    MPI.MPI_Comm,
                    DMBoundaryType, DMBoundaryType, DMBoundaryType,
                    $PetscInt, $PetscInt, $PetscInt,
                    $PetscInt, $PetscInt, $PetscInt,
                    $PetscInt,
                    $PetscInt,
                    $PetscInt,
                    $PetscInt,
                    DMStagStencilType,
                    $PetscInt,
                    Ptr{$PetscInt},  Ptr{$PetscInt},  Ptr{$PetscInt},
                    Ptr{CDMStag}
                ),
                comm,
                bndx, bndy, bndz,
                M, N, P,
                m, n ,p ,
                dofVertex ,
                dofEdge ,
                dofFace ,
                dofElement ,
                stencilType ,
                stencilWidth ,
                lx ,ly ,lz ,
                dm
                )
    end
    dmsetfromoptions && setfromoptions!(dm)
    dmsetup && setup!(dm)

    # We can only let the garbage collect finalize when we do not need to
    # worry about MPI (since garbage collection is asyncronous)
    if comm == MPI.COMM_SELF
        finalizer(destroy, dm)
    end

    return dm
end

"""
    dm = DMStagCreateCompatibleDMStag(
        dm::DMStag,
        dofVertex,
        dofEdge,
        dofFace,
        dofElement;
        kwargs...)

Creates a compatible DMStag with different dof/stratum

        dm              -   the DMStag object
        dofVertex       -   [=0] number of degrees of freedom per vertex/point/node/0-cell
        dofEdge         -   [=0] number of degrees of freedom per edge/1-cell
        dofFace         -   [=0] number of degrees of freedom per face/2-cell
        dofElement      -   [=0] number of degrees of freedom per element/3-cell
        kwargs...       -   [Optional] keyword arguments (see PETSc webpage), specifiable as stag_grid_x=100, etc.

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagCreateCompatibleDMStag.html)
"""
function DMStagCreateCompatibleDMStag end

@for_petsc function DMStagCreateCompatibleDMStag(
    dm::DMStag{$PetscLib},
    dofVertex=0,
    dofEdge=0,
    dofFace=0,
    dofElement=0;
    dmsetfromoptions=true,
    dmsetup=true,
    options...
    )

    opts = Options($petsclib, options...)

    dmnew = DMStag{$PetscLib}(C_NULL, opts)
    comm  = getcomm(dm);

    with(dm.opts) do
        @chk ccall((:DMStagCreateCompatibleDMStag, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            $PetscInt,
            $PetscInt,
            $PetscInt,
            $PetscInt,
            Ptr{CDMStag}
        ),
        dm,
        dofVertex,
        dofEdge,
        dofFace,
        dofElement,
        dmnew
        )
    end

    dmsetfromoptions && setfromoptions!(dm)
    dmsetup && setup!(dm)

    if comm == MPI.COMM_SELF
        finalizer(destroy, dmnew)
    end

    return dmnew

end

"""
    dof0,dof1,dof2,dof3 = DMStagGetDOF(dm::DMStag)

Get number of DOF associated with each stratum of the grid.

    dm      - the DMStag object
    dof0 	- the number of points per 0-cell (vertex/node)
    dof1 	- the number of points per 1-cell (element in 1D, edge in 2D and 3D)
    dof2 	- the number of points per 2-cell (element in 2D, face in 3D)
    dof3 	- the number of points per 3-cell (element in 3D)

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetDOF.html)
"""
function DMStagGetDOF end

@for_petsc function DMStagGetDOF(dm::DMStag{$PetscLib})

    dof0 = Ref{$PetscInt}()
    dof1 = Ref{$PetscInt}()
    dof2 = Ref{$PetscInt}()
    dof3 = Ref{$PetscInt}()

    @chk ccall((:DMStagGetDOF, $petsc_library), PetscErrorCode,
    (
        CDMStag,
        Ptr{$PetscInt},
        Ptr{$PetscInt},
        Ptr{$PetscInt},
        Ptr{$PetscInt}),
    dm,
    dof0,
    dof1,
    dof2,
    dof3
    )

    dim   = getdimension(dm)

    if dim==1
        return  dof0[],dof1[]
    elseif dim==2
        return dof0[],dof1[],dof2[]
    elseif dim==3
        return dof0[],dof1[],dof2[],dof3[]
    end

end

"""
    M,N,P = DMStagGetGlobalSizes(dm::DMStag)

Gets the global size of the DMStag object

    dm      - the DMStag object
    M,N,P   - size in x,y,z

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetGlobalSizes.html)
"""
function DMStagGetGlobalSizes end

@for_petsc function DMStagGetGlobalSizes(dm::DMStag{$PetscLib})

    M = Ref{$PetscInt}()
    N = Ref{$PetscInt}()
    P = Ref{$PetscInt}()

    @chk ccall((:DMStagGetGlobalSizes, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt}
        ),
        dm,
        M, N, P
        )

        dim   = getdimension(dm)

    if dim==1
        return M[]
    elseif dim==2
        return M[], N[]
    elseif dim==3
        return M[], N[], P[]
    end
end

@for_petsc function Base.size(dm::DMStag{$PetscLib})
    size = DMStagGetGlobalSizes(dm)
    return size
end

"""
    M,N,P = DMStagGetLocalSizes(dm::DMStag)

Gets the local size of the DMStag object

    dm      - the DMStag object
    M,N,P   - size in x,y,z

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetLocalSizes.html)
"""
function DMStagGetLocalSizes end

@for_petsc function DMStagGetLocalSizes(dm::DMStag{$PetscLib})

    M = Ref{$PetscInt}()
    N = Ref{$PetscInt}()
    P = Ref{$PetscInt}()

    @chk ccall((:DMStagGetLocalSizes, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt}
        ),
        dm,
        M, N, P
        )

        dim   = getdimension(dm)

    if dim==1
        return M[]
    elseif dim==2
        return M[], N[]
    elseif dim==3
        return M[], N[], P[]
    end
end

"""
    setuniformcoordinatesproduct!(
        dm::DMStag,
        xyzmin::NTuple{N, Real},
        xyzmax::NTuple{N, Real},
        )

Set uniform coordinates for the `dmstag` using the lower and upper corners defined
    by the `NTuple`s `xyzmin` and `xyzmax`. If `N` is less than the dimension of the
    `dmstag` then the value of the trailing coordinates is set to `0`.

# External Links
$(_doc_external("DMSTAG/DMStagSetUniformCoordinatesProduct"))

"""
function setuniformcoordinatesproduct! end

@for_petsc function setuniformcoordinatesproduct!(
    dm::DMStag{$PetscLib},
    xyzmin::NTuple{N, Real},
    xyzmax::NTuple{N, Real},
    ) where {N}

    xmin = $PetscReal(xyzmin[1])
    xmax = $PetscReal(xyzmax[1])

    ymin = (N > 1) ? $PetscReal(xyzmin[2]) : $PetscReal(0)
    ymax = (N > 1) ? $PetscReal(xyzmax[2]) : $PetscReal(0)

    zmin = (N > 2) ? $PetscReal(xyzmin[3]) : $PetscReal(0)
    zmax = (N > 2) ? $PetscReal(xyzmax[3]) : $PetscReal(0)

    @chk ccall((:DMStagSetUniformCoordinatesProduct, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    $PetscReal,
                    $PetscReal,
                    $PetscReal,
                    $PetscReal,
                    $PetscReal,
                    $PetscReal,
                ),
                dm,
                xmin, xmax,
                ymin, ymax,
                zmin, zmax)

    return nothing
end

# NEED TO BE REPAIRED
@for_petsc function DMStagGetProductCoordinateArraysRead(dm::DMStag{$PetscLib})

    Arrx = Ref{$PetscScalar}()
    Arry = Ref{$PetscScalar}()
    Arrz = Ref{$PetscScalar}()

    @chk ccall((:DMStagGetProductCoordinateArraysRead, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$PetscScalar}, Ptr{$PetscScalar}, Ptr{$PetscScalar}
        ),
        dm,
        Arrx, Arry, Arrz
        )

    return Arrx[],Arry[],Arrz[]
end

# TO BE REMOVED, AS setuniformcoordinates! DOES THE SAME?
"""
    DMStagSetUniformCoordinatesExplicit(
        dm::DMStag,
        xmin, xmax,
        ymin=0, ymax=0,
        zmin=0, zmax=0
        )

Set DMStag coordinates to be a uniform grid, storing all values.

    dm                            - the DMStag object
    xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagSetUniformCoordinatesExplicit.html)
"""
function DMStagSetUniformCoordinatesExplicit end

@for_petsc function DMStagSetUniformCoordinatesExplicit(
    dm::DMStag{$PetscLib},
    xmin, xmax,
    ymin=0, ymax=0,
    zmin=0, zmax=0
    )

    @chk ccall((:DMStagSetUniformCoordinatesExplicit, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    $PetscScalar, $PetscScalar,
                    $PetscScalar, $PetscScalar,
                    $PetscScalar, $PetscScalar
                ),
                dm,
                xmin, xmax,
                ymin, ymax,
                zmin, zmax
                )

    return nothing
end

#=
# replaced with createglobalvector in dm.jl
"""
    vec = DMCreateGlobalVector(
        dm::DMStag;
        write_val=true, read_val=true
        )

Creates a global vector from a DM object.

NOTE: for now this is initialized sequentially; MPI should be added

        dm 	- the DM object
        vec - the global vector

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMCreateGlobalVector.html)
"""
function DMCreateGlobalVector end

@for_petsc function DMCreateGlobalVector(
    dm::DMStag{$PetscLib};
    write_val=true, read_val=true
    )

    v = VecSeq(C_NULL, [0.0])  # empty vector

    @chk ccall((:DMCreateGlobalVector, $petsc_library), PetscErrorCode, (CDMStag, Ptr{CVec}), dm, v)

    # Link a julia array to the values from the new vector
    # If we modify values here, it will automatically be changed in the PetcVec as well
    v.array = unsafe_localarray($PetscScalar, v.ptr; write=write_val, read=read_val)

    return v
end
=#

#=
# replaced with createlocalvector?
"""
    vec = DMCreateLocalVector(dm::DMStag;
    write_val=true, read_val=true
    )

Creates a local vector from a DM object.

NOTE: for now this is initialized sequentially; MPI should be added

        dm 	- the DM object
        vec - the local vector

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMCreateLocalVector.html)
"""
function DMCreateLocalVector end

@for_petsc function DMCreateLocalVector(
    dm::DMStag{$PetscLib};
    write_val=true, read_val=true
    )

    v = VecSeq(C_NULL, [0.0])  # empty vector

    @chk ccall((:DMCreateLocalVector, $petsc_library), PetscErrorCode, (CDMStag, Ptr{CVec}), dm, v)

    # Link a julia array to the values from the new vector
    # If we modify values here, it will automatically be changed in the PetcVec as well
    v.array = unsafe_localarray($PetscScalar, v.ptr; write=write_val, read=read_val)

    return v
end
=#

"""
    Array =  DMStagVecGetArray(dm::DMStag, v::AbstractVec)

Get access to local array (including ghost points) of the DMStag.

    dm 	  - the DMStag object
    vec   - the Vec object
    Array - the array

Once you are done with work on the array, you MUST release the memory with

    Base.finalize(Array)

Otherwise the values are not returned correctly to v
"""
function DMStagVecGetArray end

@for_petsc function DMStagVecGetArray(dm::DMStag{$PetscLib}, v::AbstractVec)
    # Note: there is actually no need to call PETSc again, as Julia has the possibility
    # to wrap an existing array into another one. Our vec already has the array wrapper,
    # so we reshape that

    # Extract array from vector. Note: we need to release this by calling
    # Base.finalize on X1!
    array       =   unsafe_localarray($PetscScalar, v.ptr;  write=true, read=true)

    X1          =   DMStagVecGetArray(dm, array)

    return X1
end

"""
    Array =  DMStagVecGetArrayRead(dm::DMStag, v::AbstractVec)

Get read-only access to a local array (including ghost points) of the DMStag

    dm    - the DMStag object
    vec   - the Vec object
    Array - the read-only array
"""
function DMStagVecGetArrayRead end

@for_petsc function DMStagVecGetArrayRead(dm::DMStag{$PetscLib}, v::AbstractVec)
    # Note: there is actually no need to call PETSc again, as Julia has the possibility
    # to wrap an existing array into another one. Our vec already has the array wrapper,
    # so we reshape that

    # Extract array from vector. Note: we need to release this by calling
    # finalize on X1!
    array       =   unsafe_localarray($PetscScalar, v.ptr;  write=false, read=true)

    X1          =   DMStagVecGetArray(dm, array)

    return X1
end

"""
    X1 = DMStagVecGetArray(dm::DMStag, v::Vector)

Returns a julia array from a vector `v`, in the same shape as the DMSTAG, which can be used to set values.
"""
function DMStagVecGetArray end

@for_petsc function DMStagVecGetArray(dm::DMStag{$PetscLib}, v::Vector)

    entriesPerElement   =   DMStagGetEntriesPerElement(dm)
    ghost_corners       =   getghostcorners(dm);
    dim                 =   getdimension(dm);

    # Dimensions of new array (see the PETSc DMStagVecGetArrayRead routine)
    dim_vec             =   [entriesPerElement; collect(ghost_corners.size[1:dim])];

    # Wrap julia vector to new vector.
    X                   =    Base.view(v,:);

    # reshape to correct format
    X                   =   reshape(v, Tuple(dim_vec))
    X1                  =   PermutedDimsArray(X, Tuple([2:dim+1;1]));   # permute to take care of different array ordering in C/Julia

    return X1
end

"""
    Array = DMStagGetGhostArrayLocationSlot(
        dm::DMStag,
        v::AbstractVec,
        loc::DMStagStencilLocation,
        dof::Int
        )

Julia routine that extracts an array related to a certain DOF. Modifying values in the array will change them in the local PetscVec. Use LocalToGlobal to update global vector values.

Usage:

    Input:
        dm           -   the DMStag object
        v,ArrayFull  -   the local vector as obtained with DMCreateLocalVector, can also be a local array
        loc          -   a DMStagStencilLocation
        dof          -   the degree of freedom on loc, which you want to extracts

    Output:
        Array       -   local array that includes the ghost points, that is linked to the vector `v`.
                            Modifying values in Array will update `v`

"""
function DMStagGetGhostArrayLocationSlot end

@for_petsc function DMStagGetGhostArrayLocationSlot(
    dm::DMStag{$PetscLib},
    v::AbstractVec{$PetscScalar},
    loc::DMStagStencilLocation,
    dof::Int
    )

    entriesPerElement   =   DMStagGetEntriesPerElement(dm)
    dim                 =   getdimension(dm);
    slot                =   DMStagGetLocationSlot(dm, loc, dof);
    slot_start          =   mod(slot,entriesPerElement);          # figure out which component we are interested in

    ArrayFull           =   DMStagVecGetArray(dm, v);             # obtain access to full array

    # now extract only the dimension belonging to the current point
    Array               =   selectdim(ArrayFull,dim+1, slot_start+1);

    return Array
end

@for_petsc function DMStagGetGhostArrayLocationSlot(
    dm::DMStag{$PetscLib},
    ArrayFull::PermutedDimsArray,
    loc::DMStagStencilLocation,
    dof::Int
    )

    entriesPerElement   =   DMStagGetEntriesPerElement(dm)
    dim                 =   getdimension(dm);
    slot                =   DMStagGetLocationSlot(dm, loc, dof);
    slot_start          =   mod(slot,entriesPerElement);          # figure out which component we are interested in

    # now extract only the dimension belonging to the current point
    Array               =   selectdim(ArrayFull,dim+1, slot_start+1);

    return Array
end

"""
    slot = DMStagGetProductCoordinateLocationSlot(
        dm::DMStag,
        loc::DMStagStencilLocation
        )

Get slot for use with local product coordinate arrays.

    dm      - the DMStag object
    loc     - the grid location
    slot    - the index to use in local arrays

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetProductCoordinateLocationSlot.html)
"""
function DMStagGetProductCoordinateLocationSlot end

#REPAIR THAT AS WELL (OR GET RID OF...)
@for_petsc function DMStagGetProductCoordinateLocationSlot(
    dm::DMStag{$PetscLib},
    loc::DMStagStencilLocation
    )

    slot = Ref{$PetscInt}()
    @chk ccall((:DMStagGetProductCoordinateLocationSlot, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    DMStagStencilLocation,
                    Ptr{$PetscInt}
                ),
                dm,
                loc,
                slot
                )

    return slot[]
end

"""
    entriesPerElement = DMStagGetEntriesPerElement(dm::DMStag)

Get number of entries per element in the local representation.

    dm                - the DMStag objects
    entriesPerElement - number of entries associated with each element in the local representation

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMSTAG/DMStagGetEntriesPerElement.html)
"""
function DMStagGetEntriesPerElement end

@for_petsc function DMStagGetEntriesPerElement(dm::DMStag)
    entriesPerElement = Ref{$PetscInt}()
    @chk ccall((:DMStagGetEntriesPerElement, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    Ptr{$PetscInt}
                ),
                dm,
                entriesPerElement
                )

    return entriesPerElement[]
end

"""
    stencilWidth = DMStagGetStencilWidth(dm::DMStag)

Get elementwise stencil width.

    dm           - the DMStag objects
    stencilWidth - stencil/halo/ghost width in elements

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetStencilWidth.html)
"""
function DMStagGetStencilWidth end

@for_petsc function DMStagGetStencilWidth(dm::DMStag)
    stencilWidth = Ref{$PetscInt}()
    @chk ccall((:DMStagGetStencilWidth, $petsc_library), PetscErrorCode,
                ( CDMStag,  Ptr{$PetscInt}), dm,  stencilWidth)

    return stencilWidth[]
end

"""

    slot = DMStagGetLocationSlot(
        dm::DMStag,
        loc::DMStagStencilLocation,
        c
        )

Get index to use in accessing raw local arrays.

    dm      - the DMStag object
    loc     - location relative to an element
    c       - component ( the degree of freedom)
    slot    - index to use

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetLocationSlot.html)
"""
function DMStagGetLocationSlot end

@for_petsc function DMStagGetLocationSlot(dm::DMStag{$PetscLib},loc::DMStagStencilLocation, c)

    slot = Ref{$PetscInt}()
    @chk ccall((:DMStagGetLocationSlot, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    DMStagStencilLocation,
                    $PetscInt,
                    Ptr{$PetscInt}
                ),
                dm,
                loc,
                c,
                slot
                )

    return slot[]
end

"""
    destroy(dm::DMStag)

Destroys a DMSTAG object and releases the memory

    dm 	- the DM object to destroy

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMDestroy.html)
"""
function destroy end

@for_petsc function destroy(dm::DMStag{$PetscLib})
    finalized($petsclib) ||
        @chk ccall((:DMDestroy, $petsc_library), PetscErrorCode, (Ptr{CDMStag},), dm)
    return nothing
end

#=
# MOVED TO dm.jl
"""
    type = gettype(dm::DMStag)

Gets the DM type name (as a string) from the DM.

    dm  - The DM
    type- The DM type name

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetType.html)
"""
function gettype end

@for_petsc function gettype(dm::DMStag{$PetscLib})
    t_r = Ref{CDMStagType}()
    @chk ccall((:DMGetType, $petsc_library), PetscErrorCode, (CDMStag, Ptr{CDMStagType}), dm, t_r)
    return unsafe_string(t_r[])
end

# REPLACED WITH ROUTINE IN DM.JL
"""
    view(dm::DMStag, viewer::Viewer=ViewerStdout(dm.comm))

Views a DMSTAG object.

    dm      - the DM object to view
    viewer  - the viewer

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMView.html)
"""
#function view end

#@for_petsc function view(dm::DMStag{$PetscLib}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, MPI.MPI_COMM_SELF))
#    @chk ccall((:DMView, $petsc_library), PetscErrorCode,
#                (CDMStag, CPetscViewer),
#            dm, viewer);
#    return nothing
#end
=#

#=
# REPLACED WITH PETSc.getcorners
"""

    x,m,nExtrax = DMStagGetCorners(dm:DMStag)   in 1D
    x[],m[],nExtrax[] = DMStagGetCorners(dm:DMStag)   in 2D or 3D

Returns the global element indices of the local region (excluding ghost points).

    dm 	    - the DMStag object
    x,y,z 	- starting element indices in each direction
    m,n,p 	- element widths in each direction
    nExtrax,nExtray,nExtraz 	- number of extra partial elements in each direction.

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetCorners.html)
"""
function DMStagGetCorners end

@for_petsc function  DMStagGetCorners(dm::DMStag{$PetscLib})

    x = Ref{$PetscInt}()
    y = Ref{$PetscInt}()
    z = Ref{$PetscInt}()
    m = Ref{$PetscInt}()
    n = Ref{$PetscInt}()
    p = Ref{$PetscInt}()
    nExtrax = Ref{$PetscInt}()
    nExtray = Ref{$PetscInt}()
    nExtraz = Ref{$PetscInt}()

    @chk ccall((:DMStagGetCorners, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt},
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt},
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt}
        ),
        dm,
        x,y,z,
        m,n,p,
        nExtrax,nExtray,nExtraz
        )

        dim = getdimension(dm);

        if dim==1
            X = (x[],)
            M = (m[],)
            NEXTRA = (nExtrax[],)
            return X[1], M[1], NEXTRA[1]
        elseif dim==2
            return (x[], y[]), (m[],n[]), (nExtrax[],nExtray[])
        elseif dim==3
            return (x[], y[], z[]), (m[],n[],p[]), (nExtrax[],nExtray[],nExtraz[])
        end
end
=#

#@for_petsc PETSc.getcorners(dm::DMStag{$PetscLib}) = DMStagGetCorners(dm::DMStag{$PetscLib})

#=
# REPLACED WITH PETSc.getghostcorners
"""
    x,m = DMStagGetGhostCorners(dm:DMStag)   in 1D
    x[],m[] = DMStagGetGhostCorners(dm:DMStag)   in 2D or 3D

Return global element indices of the local region (including ghost points).

    dm 	    - the DMStag object
    x[] 	- starting element indices in each direction
    m[] 	- element widths in each direction

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetGhostCorners.html)
"""
function DMStagGetGhostCorners end

@for_petsc function  DMStagGetGhostCorners(dm::DMStag)

    x = Ref{$PetscInt}()
    y = Ref{$PetscInt}()
    z = Ref{$PetscInt}()
    m = Ref{$PetscInt}()
    n = Ref{$PetscInt}()
    p = Ref{$PetscInt}()

    @chk ccall((:DMStagGetGhostCorners, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt},
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt}
        ),
        dm,
        x,y,z,
        m,n,p
        )

        dim =   getdimension(dm);

        if dim==1
            X = (x[],)
            M = (m[],)
            return X[1], M[1]
        elseif dim==2
            return (x[], y[]), (m[],n[])
        elseif dim==3
            return (x[], y[], z[]), (m[],n[],p[])
        end
end
=#

"""
    Indices = DMStagGetCentralNodes(dm::DMStag)

Return indices of start and end of the central nodes of a local array built from the input `dm` (excluding ghost nodes).

    dm 	        - the DMStag object
    Indices 	- indices of lower and upper range of center and vertex nodes
"""
function DMStagGetCentralNodes end

@for_petsc function DMStagGetCentralNodes(dm::DMStag)
    # In Julia, indices in arrays start @ 1, whereas they can go negative in C

    #g_start, g_N    =   getghostcorners(dm);  # MODIFY
    gc              =   getghostcorners(dm);

    #start,N, nExtra =   getcorners(dm);       # MODIFY
    c               =   getcorners(dm);

    # Note that we add the +1 for julia/petsc consistency
    center = (  x= (c.lower[1]+1):(c.upper[1]+1),
                y= (c.lower[2]+1):(c.upper[2]+1),
                z= (c.lower[3]+1):(c.upper[3]+1) )

    vertex = (  x=c.lower[1]+1:c.upper[1]+2 ,
                y=c.lower[2]+1:c.upper[2]+2 ,
                z=c.lower[3]+1:c.upper[3]+2 )

    return (center=center, vertex=vertex)

end

"""
    Bx = DMStagGetBoundaryTypes(dm::DMStag) in 1D
    Bx,By,Bz = DMStagGetBoundaryTypes(dm::DMStag) in 3D

Get boundary types.

        dm 	     - the DMStag object
        Bx,By,Bz - boundary types

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetBoundaryTypes.html)
"""
function DMStagGetBoundaryTypes end

@for_petsc function DMStagGetBoundaryTypes(dm::DMStag)

    Bx = Ref{$DMBoundaryType}()
    By = Ref{$DMBoundaryType}()
    Bz = Ref{$DMBoundaryType}()

    @chk ccall((:DMStagGetBoundaryTypes, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$DMBoundaryType}, Ptr{$DMBoundaryType}, Ptr{$DMBoundaryType}
        ),
        dm,
        Bx,By,Bz
        )

        dim = getdimension(dm);

        if dim==1
            return Bx[]
        elseif dim==2
            return Bx[], By[]
        elseif dim==3
            return Bx[], By[], Bz[]
        end
end

"""
    nRanks0 = DMStagGetNumRanks(dm::DMStag) in 1D
    nRanks0,nRanks1,nRanks2 = DMStagGetNumRanks(dm::DMStag) in 3D

Get number of ranks in each direction in the global grid decomposition.

        dm 	                     - the DMStag object
        nRanks0,nRanks1,nRanks2  - number of ranks in each direction in the grid decomposition

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetNumRanks.html)
"""
function DMStagGetNumRanks end

@for_petsc function  DMStagGetNumRanks(dm::DMStag)

    nRanks0 = Ref{$PetscInt}()
    nRanks1 = Ref{$PetscInt}()
    nRanks2 = Ref{$PetscInt}()

    @chk ccall((:DMStagGetNumRanks, $petsc_library), PetscErrorCode,
        (
            CDMStag,
            Ptr{$PetscInt}, Ptr{$PetscInt}, Ptr{$PetscInt}
        ),
        dm,
        nRanks0,nRanks1,nRanks2
        )

        dim = getdimension(dm);

    if dim==1
        return nRanks0[]
    elseif dim==2
        return nRanks0[], nRanks1[]
    elseif dim==3
        return nRanks0[], nRanks1[], nRanks2[]
    end
end

"""
    DMStagVecSetValuesStencil(
        dm::DMStag,
        vec::Abstractvec,
        pos::DMStagStencil,
        val::Float64,
        insertMode::InsertMode
        )

This puts a value inside a global vector using DMStagStencil.

    dm  - the DMStag object
    vec - the Vec
    n   - the number of values (do not fill if only 1)
    pos - the location of the set values, given by a DMStagStencil struct
    val - the value to be set
    insertMode  - INSERT_VALUES or ADD_VALUES

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagVecSetValuesStencil.html)
"""
function DMStagVecSetValuesStencil end

@for_petsc function  DMStagVecSetValuesStencil(
        dm::DMStag{$PetscLib},
        vec::AbstractVec{$PetscScalar},
        pos::DMStagStencil{$PetscInt},
        val,
        insertMode::InsertMode)

    n=1;
    @chk ccall((:DMStagVecSetValuesStencil, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    CVec,
                    $PetscInt,
                    Ptr{DMStagStencil{$PetscInt}},
                    Ptr{$PetscScalar},
                    InsertMode
                ),
                dm,
                vec.ptr,
                n,
                Ref{DMStagStencil{$PetscInt}}(pos),
                Ref{$PetscScalar}(val),
                insertMode
                )

    return nothing
end

"""

    DMStagVecSetValuesStencil(
        dm::DMStag,
        vec::AbstractVec{PetscScalar},
        n,
        pos::Vector{DMStagStencil},
        values::Vector{PetscScalar},
        insertMode::InsertMode
        )

This puts values inside a global vector using DMStagStencil

    dm  - the DMStag object
    vec - the Vec
    n   - the number of values (do not fill if only 1)
    pos - the location of the set values, given by a DMStagStencil struct
    val - the value to be set
    insertMode  - INSERT_VALUES or ADD_VALUES

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagVecSetValuesStencil.html)
"""
function DMStagVecSetValuesStencil end

@for_petsc function  DMStagVecSetValuesStencil(
    dm::DMStag{$PetscLib},
    vec::AbstractVec{$PetscScalar},
    n::Integer,
    pos::Vector{$DMStagStencil{$PetscInt}},
    values::Vector{$PetscScalar},
    insertMode::InsertMode
    )

    i = 1;
    while i <= n
        pos0 = pos[i];
        val  = values[i];
        m=1;
        @chk ccall((:DMStagVecSetValuesStencil, $petsc_library), PetscErrorCode,
                    (
                        CDMStag,
                        CVec,
                        $PetscInt,
                        Ptr{DMStagStencil{$PetscInt}},
                        Ptr{$PetscScalar},
                        InsertMode
                    ),
                    dm,
                    vec.ptr,
                    m,
                    Ref{DMStagStencil{$PetscInt}}(pos0),
                    Ref{$PetscScalar}(val),
                    insertMode
                    )
    i += 1;
    end
    return nothing
end

"""
    val = DMStagVecGetValuesStencil(
        dm::DMStag,
        vec::AbstractVec,
        pos::DMStagStencil
        )

Get vector values using grid indexing

    dm  - the DMStag object
    vec - the vector object
    n   - the number of values to obtain (do not fill if only one)
    pos - locations to obtain values from (as an array of DMStagStencil values)
    val - value at the point

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagVecGetValuesStencil.html)
"""
function DMStagVecGetValuesStencil end

@for_petsc function  DMStagVecGetValuesStencil(
    dm::DMStag{$PetscLib},
    vec::AbstractVec{$PetscScalar},
    pos::DMStagStencil
    )

    n=1;
    val = Ref{$PetscScalar}()
    @chk ccall((:DMStagVecGetValuesStencil, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    CVec,
                    $PetscInt,
                    Ptr{DMStagStencil{$PetscInt}},
                    Ptr{$PetscScalar}
                ),
                dm,
                vec.ptr,
                n,
                Ref{DMStagStencil{$PetscInt}}(pos),
                val
                )

    return val[]
end

"""
    val = DMStagVecGetValuesStencil(
        dm::DMStag,
        vec::AbstractVec,
        n,
        pos::Vector{DMStagStencil}
        )

Get vector values using grid indexing.

    dm  - the DMStag object
    vec - the vector object
    n   - the number of values to obtain (do not fill if only one)
    pos - locations to obtain values from (as an array of DMStagStencil values)
    val - value at the point

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagVecGetValuesStencil.html)
"""
function DMStagVecGetValuesStencil end

@for_petsc function  DMStagVecGetValuesStencil(
    dm::DMStag{$PetscLib},
    vec::AbstractVec{$PetscScalar},
    n,
    pos::Vector{$DMStagStencil{$PetscInt}}
    )

    i = 1;
    values = zeros(n);
    while i <= n
        pos0 = pos[i];
        m=1;
        val = Ref{$PetscScalar}()
        @chk ccall((:DMStagVecGetValuesStencil, $petsc_library), PetscErrorCode,
                    (
                        CDMStag,
                        CVec,
                        $PetscInt,
                        Ptr{DMStagStencil{$PetscInt}},
                        Ptr{$PetscScalar}
                    ),
                    dm,
                    vec.ptr,
                    m,
                    Ref{DMStagStencil{$PetscInt}}(pos0),
                    val
                    )

        values[i] = val[];
        i += 1;
    end
    return values
end

"""
    val =  DMStagMatGetValuesStencil(
        dm::DMStag,
        mat::AbstractMat,
        posRow::DMStagStencil,
        posCol::DMStagStencil
        )

This reads a single value from a matrix DMStagStencil

    dm      - the DMStag object
    mat     - the Mat
    posRow  - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    posCol  - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    val     - the value

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagMatGetValuesStencil.html)
"""
function DMStagMatGetValuesStencil end

@for_petsc function  DMStagMatGetValuesStencil(
    dm::DMStag{$PetscLib},
    mat::AbstractMat{$PetscScalar},
    posRow::DMStagStencil,
    posCol::DMStagStencil
    )

    nRow= 1;
    nCol= 1;
    val = Ref{$PetscScalar}()
    @chk  ccall((:DMStagMatGetValuesStencil, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    CMat,
                    $PetscInt,
                    Ptr{DMStagStencil{$PetscInt}},
                    $PetscInt,
                    Ptr{DMStagStencil{$PetscInt}},
                    Ptr{$PetscScalar}
                ),
                dm,
                mat.ptr,
                nRow,
                Ref{DMStagStencil{$PetscInt}}(posRow),
                nCol,
                Ref{DMStagStencil{$PetscInt}}(posCol),
                val
                )

    return val[]
end

"""
    val =  DMStagMatGetValuesStencil(
        dm::DMStag,
        mat::AbstractMat{PetscScalar},
        nRow,
        posRow::Vector{DMStagStencil},
        nCol,
        posCol::Vector{DMStagStencil}
        )

This reads a single value from a matrix DMStagStencil.

    dm      - the DMStag object
    mat     - the Mat
    posRow  - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    posCol  - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    val     - the value

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagMatGetValuesStencil.html)
"""
function DMStagMatGetValuesStencil end

@for_petsc function  DMStagMatGetValuesStencil(
    dm::DMStag{$PetscLib},
    mat::AbstractMat{$PetscScalar},
    nRow,
    posRow::Vector{$DMStagStencil{$PetscInt}},
    nCol,
    posCol::Vector{$DMStagStencil{$PetscInt}}
    )

    i = 1;
    j = 1;
    values = zeros(nRow*nCol);
    while i <= nRow
        while j <= nCol
            posr = posRow[i];
            posc = posCol[j];
            n_Row= 1;
            n_Col= 1;
            val = Ref{$PetscScalar}()
            @chk  ccall((:DMStagMatGetValuesStencil, $petsc_library), PetscErrorCode,
                        (
                            CDMStag,
                            CMat,
                            $PetscInt,
                            Ptr{DMStagStencil{$PetscInt}},
                            $PetscInt,
                            Ptr{DMStagStencil{$PetscInt}},
                            Ptr{$PetscScalar}
                        ),
                        dm,
                        mat.ptr,
                        n_Row,
                        Ref{DMStagStencil{$PetscInt}}(posr),
                        n_Col,
                        Ref{DMStagStencil{$PetscInt}}(posc),
                        val
                        )
            values[i*j] = val[]
            j += 1;
        end
        i += 1;
    end

    return values
end

"""
    indices = LocalInGlobalIndices(dm::DMStag)

Give the non-ghosted indices in the local vector that contribute to the global vector.

    dm      - the DMStag object
    indices - local indices
"""
function LocalInGlobalIndices end

function LocalInGlobalIndices(dm::DMStag)
    # note: this can likely be done more efficiently and will have to be modified in parallel
    ind_g   =   createglobalvector(dm)
    v_ind_l =   createlocalvector(dm)

    # Set indices in local vector
    for i=1:length(v_ind_l)
        v_ind_l[i] = i
    end
    @show v_ind_l
    update!(ind_g, v_ind_l, INSERT_VALUES); # update global vector

    #@show v_ind_l ind_g

    return Int64.(ind_g)
end

"""

    DMStagMatSetValuesStencil(
        dm::DMStag,
        mat::AbstractMat,
        nRow,
        posRow::Vector{DMStagStencil},
        nCol,
        posCol::Vector{DMStagStencil},
        values::Vector{PetscScalar},
        insertMode::InsertMode
        )

This puts values inside a matrix using DMStagStencil position.

    dm	        - the DMStag object
    mat	        - the Mat
    posRow	    - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    posCol	    - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    val	        - the value to be set
    insertMode	- INSERT_VALUES or ADD_VALUES

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagMatSetValuesStencil.html)
"""
function DMStagMatSetValuesStencil end

@for_petsc function  DMStagMatSetValuesStencil(
    dm::DMStag{$PetscLib},
    mat::AbstractMat{$PetscScalar},
    posRow::DMStagStencil,
    posCol::DMStagStencil,
    val,
    insertMode::InsertMode
    )

    nRow= 1;
    nCol= 1;
    @chk ccall((:DMStagMatSetValuesStencil, $petsc_library), PetscErrorCode,
                (
                    CDMStag,
                    CMat,
                    $PetscInt,
                    Ptr{DMStagStencil{$PetscInt}},
                    $PetscInt,
                    Ptr{DMStagStencil{$PetscInt}},
                    Ptr{$PetscScalar},
                    InsertMode
                ),
                dm,
                mat.ptr,
                nRow,
                Ref{DMStagStencil{$PetscInt}}(posRow),
                nCol,
                Ref{$DMStagStencil{$PetscInt}}(posCol),
                Ref{$PetscScalar}(val),
                insertMode
                )

    return nothing
end

"""
    DMStagMatSetValuesStencil(
        dm::DMStag,
        mat::AbstractMat,
        nRow,
        posRow::Vector{DMStagStencil},
        nCol,
        posCol::Vector{DMStagStencil},
        values::Vector{PetscScalar},
        insertMode::InsertMode
        )

This puts values inside a matrix using DMStagStencil position

    dm	        - the DMStag object
    mat	        - the Mat
    posRow	    - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    posCol	    - the location of the row of the set value, given by a DMStagStencil struct (as a vector)
    val	        - the value to be set
    insertMode	- INSERT_VALUES or ADD_VALUES

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagMatSetValuesStencil.html)
"""
function DMStagMatSetValuesStencil end

@for_petsc function  DMStagMatSetValuesStencil(
    dm::DMStag{$PetscLib},
    mat::AbstractMat{$PetscScalar},
    nRow,
    posRow::Vector{$DMStagStencil{$PetscInt}},
    nCol,
    posCol::Vector{$DMStagStencil{$PetscInt}},
    values::Vector{$PetscScalar},
    insertMode::InsertMode
    )

    i = 1;
    j = 1;
    while i <= nRow
        while j <= nCol
            posr = posRow[i];
            posc = posCol[j];
            val  = values[i*j];
            n_Row= 1;
            n_Col= 1;
            @chk ccall((:DMStagMatSetValuesStencil, $petsc_library), PetscErrorCode,
                        (
                            CDMStag,
                            CMat,
                            $PetscInt,
                            Ptr{DMStagStencil{$PetscInt}},
                            $PetscInt,
                            Ptr{DMStagStencil{$PetscInt}},
                            Ptr{$PetscScalar},
                            InsertMode
                        ),
                        dm,
                        mat.ptr,
                        n_Row,
                        Ref{DMStagStencil{$PetscInt}}(posr),
                        n_Col,
                        Ref{$DMStagStencil{$PetscInt}}(posc),
                        Ref{$PetscScalar}(val),
                        insertMode
                        )
        j += 1;
        end
    i += 1;
    end
    return nothing
end

#=
# REPLACED WITH   PETSc.getdimension
"""
    dim = DMGetDimension(dm::DMStag)

Return the topological dimension of the DM.

    dm 	- The DM
    dim - dimensions

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetDimension.html)
"""
function DMGetDimension end

@for_petsc function DMGetDimension(dm::DMStag)
    dim = Ref{$PetscInt}()

    @chk ccall((:DMGetDimension, $petsc_library), PetscErrorCode, (CDMStag,Ptr{$PetscInt}), dm, dim )

    return dim[]
end
=#

"""
    DMLocalToGlobal(
        dm::DMStag,
        l::AbstractVec,
        mode::InsertMode,
        g::AbstractVec
        )

Updates global vectors from local vectors.

    dm 	 - the DM object
    l 	 - the local vector
    mode - if INSERT_VALUES then no parallel communication is used, if ADD_VALUES then all ghost points from the same base point accumulate into that base point.
    g 	 - the global vector

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMLocalToGlobal.html)
"""
function DMLocalToGlobal end

@for_petsc function DMLocalToGlobal(
    l::AbstractVec{$PetscScalar},
    mode::InsertMode,
    g::AbstractVec{$PetscScalar}
    )

    #DMLocalToGlobal(dm,l.ptr, mode,g.ptr)
    update!(l,g,mode);

    return nothing
end

@for_petsc function DMLocalToGlobal(
    l::AbstractVec{$PetscScalar},
    mode::InsertMode,
    g::CVec
    )

    #DMLocalToGlobal(dm,l.ptr, mode,g)
    update!(l,g,mode);

    return nothing
end

#=
# replaced with update!
@for_petsc function DMLocalToGlobal(
    dm::DMStag{$PetscLib},
    l::CVec,
    mode::InsertMode,
    g::CVec
    )

    @chk ccall((:DMLocalToGlobal, $petsc_library), PetscErrorCode,
    (
        CDMStag,
        CVec,
        InsertMode,
        CVec
    ),
    dm,
    l,
    mode,
    g)

    return nothing
end
=#

@for_petsc function DMGlobalToLocal(
    g::AbstractVec{$PetscScalar},
    mode::InsertMode,
    l::AbstractVec{$PetscScalar}
    )

    #DMGlobalToLocal(dm,g.ptr, mode::InsertMode,l.ptr)
    update!(g,l,mode)

    return nothing
end

@for_petsc function DMGlobalToLocal(
    g::CVec,
    mode::InsertMode,
    l::AbstractVec{$PetscScalar}
    )

   # DMGlobalToLocal(dm,g, mode::InsertMode,l.ptr)
    update!(g,l,mode)

    return nothing
end

#=
# replaced with update!
"""
    DMGlobalToLocal(
        dm::DMStag,
        g::CVec,
        mode::InsertMode,
        l::CVec
        )

Update local vectors from global vector.

    dm   - the DM object
    g    - the global vector
    mode - INSERT_VALUES or ADD_VALUES
    l    - the local vector

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGlobalToLocal.html)
"""
function DMGlobalToLocal end

@for_petsc function DMGlobalToLocal(
    dm::DMStag{$PetscLib},
    g::CVec,
    mode::InsertMode,
    l::CVec
    )

    @chk ccall((:DMGlobalToLocal, $petsc_library), PetscErrorCode,
    (
        CDMStag,
        CVec,
        InsertMode,
        CVec
    ),
    dm,
    g,
    mode,
    l
    )

    return nothing
end
=#

#=
# REPLACED WITH PETSc.creatematrix?
"""
    mat = DMCreateMatrix(dm::DMStag)

Generates a matrix from a DMStag object. The type is a MatSeqAIJ is we are on 1 core.

    dm  - the DMStag object
    mat - the matrix of type MatSeqAIJ (on 1 core)

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMCreateMatrix.html)
"""
function DMCreateMatrix end

@for_petsc function DMCreateMatrix(dm::DMStag)
    # Note: the matrix cannot be viewed yet, as it remains unassembled
    #  ideally, we should modify the viewer to take care of this case

    if dm.comm==MPI.COMM_SELF
        mat = MatSeqAIJ{$PetscScalar}(C_NULL)
    elseif dm.comm==MPI.COMM_WORLD
        error("MatMPIAIJ still to be implemented")
    end

    @chk ccall((:DMCreateMatrix, $petsc_library), PetscErrorCode,
                (CDMStag, Ptr{CMat}), dm, mat)

    return mat
end
=#

"""
    stencilType = DMStagGetStencilType(dm::DMStag)

Get elementwise ghost/halo stencil type.

    dm          - the DMStag object
    stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetStencilType.html)
"""
function DMStagGetStencilType end

@for_petsc function DMStagGetStencilType(dm::DMStag)
    stencilType =  Ref{DMStagStencilType}()

    @chk ccall((:DMStagGetStencilType, $petsc_library), PetscErrorCode,
    (
        CDMStag,
        Ptr{DMStagStencilType}
    ),
    dm,
    stencilType
    )

    return stencilType[]
end

"""
    fr_X,fr_Y,fr_Z = DMStagGetIsFirstRank(dm::DMStag)

Returns boolean value to indicate whether this rank is first in each direction in the rank grid. Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

    dm             - the DMStag object
    fr_X,fr_Y,fr_Z - whether this rank is first in each direction

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetIsFirstRank.html)
"""
function DMStagGetIsFirstRank end

@for_petsc function DMStagGetIsFirstRank(dm::DMStag)
    fr_X = Ref{PetscBool}()
    fr_Y = Ref{PetscBool}()
    fr_Z = Ref{PetscBool}()

    @chk ccall((:DMStagGetIsFirstRank, $petsc_library), PetscErrorCode,
    (
        CDMStag,
        Ptr{PetscBool}, Ptr{PetscBool}, Ptr{PetscBool}
    ),
    dm,
    fr_X, fr_Y, fr_Z
    )

    return fr_X[]== PETSC_TRUE, fr_Y[]== PETSC_TRUE, fr_Z[]== PETSC_TRUE
end

"""
    fr_X,fr_Y,fr_Z = DMStagGetIsLastRank(dm::DMStag)

Returns boolean value to indicate whether this rank is last in each direction in the rank grid.

    dm             - the DMStag object
    fr_X,fr_Y,fr_Z - whether this rank is last in each direction

From [PETSc Manual])(https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagGetIsLastRank.html)
"""
function DMStagGetIsLastRank end

@for_petsc function DMStagGetIsLastRank(dm::DMStag)
    fr_X = Ref{PetscBool}()
    fr_Y = Ref{PetscBool}()
    fr_Z = Ref{PetscBool}()

    @chk ccall((:DMStagGetIsLastRank, $petsc_library), PetscErrorCode,
    (
        CDMStag,
        Ptr{PetscBool}, Ptr{PetscBool}, Ptr{PetscBool}
    ),
    dm,
    fr_X, fr_Y, fr_Z
    )

    return fr_X[]== PETSC_TRUE, fr_Y[]== PETSC_TRUE, fr_Z[]== PETSC_TRUE
end

#=
# REPLACED WITH getcoordinateDM
"""
    dmnew = DMGetCoordinateDM(dm::DMStag; kwargs...)

Gets the DM that prescribes coordinate layout and scatters between global and local coordinates.

    dm          -    the DMStag object
    kwargs...   -   [Optional] keyword arguments (see PETSc webpage), specifiable as stag_grid_x=100, etc.
    dmnew       -    Coordinate DM

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetCoordinateDM.html)
"""
function DMGetCoordinateDM end

@for_petsc function DMGetCoordinateDM(
    dm::DMStag{$PetscLib}; kwargs...)

    opts = Options($PetscLib, kwargs...)

    dmnew = DMStag{$PetscLib}(C_NULL, opts)   # retrieve options

    @chk ccall((:DMGetCoordinateDM, $petsc_library), PetscErrorCode,
    (CDMStag, Ptr{CDMStag}),
     dm, dmnew)

    return dmnew
end
=#

#=
# replaced with getcoordinateslocal
"""
    v = DMGetCoordinatesLocal(dm::DMStag; write_val=true, read_val=true)

Gets a local vector with the coordinates associated with the DM.

    dm 	- the DMStag object
    v   - coordinate local vector

From [PETSc Manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetCoordinatesLocal.html)
"""
function DMGetCoordinatesLocal end

@for_petsc function DMGetCoordinatesLocal(
    ::$UnionPetscLib,
    dm::DMStag; write_val=true, read_val=true)

    # Should probably be a parallel vector
    v = VecSeq(C_NULL, [0.0])  # empty vector

    @chk ccall((:DMGetCoordinatesLocal, $petsc_library), PetscErrorCode, (CDMStag, Ptr{CVec}), dm, v)

    # Link a julia array to the values from the new vector
    # If we modify values here, it will automatically be changed in the PetcVec as well
    v.array = unsafe_localarray($PetscScalar, v.ptr; write=write_val, read=read_val)

    return v
end
=#

"""
    setuniformcoordinates!(
        dm::DMStag,
        xyzmin::NTuple{N, Real},
        xyzmax::NTuple{N, Real},
    ) where {N}

Set uniform coordinates for the `dmstag` using the lower and upper corners defined
by the `NTuples` `xyzmin` and `xyzmax`. If `N` is less than the dimension of the
`dm` then the value of the trailing coordinates is set to `0`.

# External Links
$(_doc_external("DMSTAG/DMStagSetUniformCoordinatesExplicit"))
"""
function setuniformcoordinates!(dm::DMStag,xyzmin,xyzmax) end

@for_petsc function setuniformcoordinates!(
    dm::DMStag{$PetscLib},
    xyzmin::NTuple{N, Real},
    xyzmax::NTuple{N, Real},
) where {N}
    xmin = $PetscReal(xyzmin[1])
    xmax = $PetscReal(xyzmax[1])

    ymin = (N > 1) ? $PetscReal(xyzmin[2]) : $PetscReal(0)
    ymax = (N > 1) ? $PetscReal(xyzmax[2]) : $PetscReal(0)

    zmin = (N > 2) ? $PetscReal(xyzmin[3]) : $PetscReal(0)
    zmax = (N > 2) ? $PetscReal(xyzmax[3]) : $PetscReal(0)

    @chk ccall(
        (:DMStagSetUniformCoordinatesExplicit, $petsc_library),
        PetscErrorCode,
        (
            CDMStag,
            $PetscReal,
            $PetscReal,
            $PetscReal,
            $PetscReal,
            $PetscReal,
            $PetscReal,
        ),
        dm,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
    )
    return nothing
end

"""
    getcorners(dm::DMSTAG)

Returns a `NamedTuple` with the global indices (excluding ghost points) of the
`lower` and `upper` corners as well as the `size`.

# External Links
$(_doc_external("DMSTAG/DMStagGetCorners"))
"""
function getcorners(dm::DMStag) end

@for_petsc function getcorners(dm::DMStag{$PetscLib})
  #  info = DMDALocalInfo{$PetscInt}()
    corners     = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    local_size  = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    nExtra      = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    @chk ccall(
        (:DMStagGetCorners, $petsc_library),
        PetscErrorCode,
        (
            CDMStag,
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
        ),
        dm,
        Ref(corners, 1),
        Ref(corners, 2),
        Ref(corners, 3),
        Ref(local_size, 1),
        Ref(local_size, 2),
        Ref(local_size, 3),
        Ref(nExtra,     1),
        Ref(nExtra,     2),
        Ref(nExtra,     3),
    )
    corners .+= 1
    return (
        lower = corners,
        upper = corners .+ local_size .- $PetscInt(1),
        size  = local_size,
        extra = nExtra
    )
end

"""
    getghostcorners(dm::DMStag)

Returns a `NamedTuple` with the global indices (including ghost points) of the
`lower` and `upper` corners as well as the `size`.

# External Links
$(_doc_external("DMSTAG/DMStagGetGhostCorners"))
"""
function getghostcorners(dm::DMStag) end

@for_petsc function getghostcorners(dm::DMStag)
    corners     = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    local_size  = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    nExtra      = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    @chk ccall(
        (:DMStagGetGhostCorners, $petsc_library),
        PetscErrorCode,
        (
            CDMStag,
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
        ),
        dm,
        Ref(corners, 1),
        Ref(corners, 2),
        Ref(corners, 3),
        Ref(local_size, 1),
        Ref(local_size, 2),
        Ref(local_size, 3),
        Ref(nExtra,     1),
        Ref(nExtra,     2),
        Ref(nExtra,     3),
    )
    corners .+= 1
    return (
        lower = corners,
        upper = corners .+ local_size .- 1,
        size  = local_size,
        extra = nExtra
    )
end

Base.show(io::IO, dm::DMStag) = _show(io, dm)


text:
abstract type AbstractNodeInput end

struct Linear{TA<:AbstractActivation} <: AbstractNodeInput
    activation::TA
end

Linear(::Type{TA}) where TA = Linear(TA())

# special case the event when both x and Wt are vectors
input_to_node(::Type{Linear}, x::AbstractVector, Wt::AbstractVector, d) =
    dot(x, Wt) + d'

# otherwise let dispatch work its magic
input_to_node(::Type{T}, x, Wt, d) where {T<:Linear} = x*Wt .+ d'
input_to_node(::Linear, x, Wt, d) = input_to_node(Linear, x, Wt, d)
activate(l::Linear, h) = activate(l.activation, h)
hidden_out(l::Linear, x, Wt, d) = activate(l, input_to_node(l, x, Wt, d))

# Radial basis functions
abstract type AbstractRBFFamily end

struct RBF{TF<:AbstractRBFFamily,TD<:Union{PreMetric,SemiMetric}} <: AbstractNodeInput
    dist::TD
end

function RBF(::Union{TF,Type{TF}}, d::TD=SqEuclidean()) where TF<:AbstractRBFFamily where TD<:Union{PreMetric,SemiMetric}
    RBF{TF,TD}(d)
end

# Gaussian RBF
struct Gaussian <: AbstractRBFFamily end

input_to_node(r::RBF{Gaussian}, x::AbstractVector, c::AbstractVector, σ::Real) =
    exp(-evaluate(r.dist, x, c)/σ)

input_to_node(r::RBF{Gaussian}, x::AbstractMatrix, c::AbstractVector, σ::Real) =
    exp(-colwise(r.dist, x', c) ./ σ)

function input_to_node(r::RBF{Gaussian}, x::AbstractMatrix, c::AbstractMatrix,
                       σ::AbstractVector)
    _out = Array{eltype(x)}(undef, size(c, 2), size(x, 1))
    xp = x'
    for i in 1:size(x, 1)
        _out[:, i] = input_to_node(r, view(xp, :, i), c, σ)
    end
    _out'
end

function input_to_node(r::RBF{Gaussian}, x::AbstractVector, c::AbstractMatrix,
                       σ::AbstractVector)
    if size(c, 1) == 1
        # q = 1, s > 1; many points
        exp(-pairwise(r.dist, x', c) ./ σ')
    else
        exp(-colwise(r.dist, x, c') ./ σ)
    end
end

# Special case scalar x version
input_to_node(r::RBF{Gaussian,SqEuclidean}, x::Real, c, σ) =
    exp(-(x-c).^2 ./ σ')

input_to_node(r::RBF{Gaussian,Euclidean}, x::Real, c, σ) =
    exp(-sqrt((x-c).^2) ./ σ')

hidden_out(rbf::RBF, x, c, σ) = input_to_node(rbf, x, c, σ)

ext: jl


text:
using LinearAlgebra, NLPModels, Stopping

include("NewtonSolver.jl")
include("Armijo.jl")
include("NewtonStop.jl")
include("prepare_LS.jl")
include("NewtonStopLS.jl")


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	investment_energy(EP::Model, inputs::Dict)

This function defines the expressions and constraints keeping track of total available storage charge capacity across all resources as well as constraints on capacity retirements. The function also adds investment and fixed O\&M related costs related to charge capacity to the objective function.

The total capacity of each resource is defined as the sum of the existing capacity plus the newly invested capacity minus any retired capacity.

```math
\begin{aligned}
& \Delta^{total,energy}_{y,z} =(\overline{\Delta^{energy}_{y,z}}+\Omega^{energy}_{y,z}-\Delta^{energy}_{y,z}) \forall y \in \mathcal{O}, z \in \mathcal{Z}
\end{aligned}
```

One cannot retire more capacity than existing capacity.

```math
\begin{aligned}
&\Delta^{energy}_{y,z} \leq \overline{\Delta^{energy}_{y,z}}
		\hspace{4 cm}  \forall y \in \mathcal{O}, z \in \mathcal{Z}
\end{aligned}
```

For resources where $\overline{\Omega_{y,z}^{energy}}$ and $\underline{\Omega_{y,z}^{energy}}$ is defined, then we impose constraints on minimum and maximum power capacity.

```math
\begin{aligned}
& \Delta^{total,energy}_{y,z} \leq \overline{\Omega}^{energy}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{O}, z \in \mathcal{Z} \\
& \Delta^{total,energy}_{y,z}  \geq \underline{\Omega}^{energy}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{O}, z \in \mathcal{Z}
\end{aligned}
```

In addition, this function adds investment and fixed O\&M related costs related to charge capacity to the objective function:

```math
\begin{aligned}
& 	\sum_{y \in \mathcal{O} } \sum_{z \in \mathcal{Z}}
	\left( (\pi^{INVEST,energy}_{y,z} \times    \Omega^{energy}_{y,z})
	+ (\pi^{FOM,energy}_{y,z} \times  \Delta^{total,energy}_{y,z})\right)
\end{aligned}
```
"""
function investment_energy(EP::Model, inputs::Dict)

	println("Storage Investment Module")

	dfGen = inputs["dfGen"]

	G = inputs["G"] # Number of resources (generators, storage, DR, and DERs)

	STOR_ALL = inputs["STOR_ALL"] # Set of all storage resources
	NEW_CAP_ENERGY = inputs["NEW_CAP_ENERGY"] # Set of all storage resources eligible for new energy capacity
	RET_CAP_ENERGY = inputs["RET_CAP_ENERGY"] # Set of all storage resources eligible for energy capacity retirements

	### Variables ###

	## Energy storage reservoir capacity (MWh capacity) built/retired for storage with variable power to energy ratio (STOR=1 or STOR=2)

	# New installed energy capacity of resource "y"
	@variable(EP, vCAPENERGY[y in NEW_CAP_ENERGY] >= 0)

	# Retired energy capacity of resource "y" from existing capacity
	@variable(EP, vRETCAPENERGY[y in RET_CAP_ENERGY] >= 0)

	### Expressions ###

	@expression(EP, eTotalCapEnergy[y in STOR_ALL],
		if (y in intersect(NEW_CAP_ENERGY, RET_CAP_ENERGY))
			dfGen[!,:Existing_Cap_MWh][y] + EP[:vCAPENERGY][y] - EP[:vRETCAPENERGY][y]
		elseif (y in setdiff(NEW_CAP_ENERGY, RET_CAP_ENERGY))
			dfGen[!,:Existing_Cap_MWh][y] + EP[:vCAPENERGY][y]
		elseif (y in setdiff(RET_CAP_ENERGY, NEW_CAP_ENERGY))
			dfGen[!,:Existing_Cap_MWh][y] - EP[:vRETCAPENERGY][y]
		else
			dfGen[!,:Existing_Cap_MWh][y] + EP[:vZERO]
		end
	)

	## Objective Function Expressions ##

	# Fixed costs for resource "y" = annuitized investment cost plus fixed O&M costs
	# If resource is not eligible for new energy capacity, fixed costs are only O&M costs
	@expression(EP, eCFixEnergy[y in STOR_ALL],
		if y in NEW_CAP_ENERGY # Resources eligible for new capacity
			dfGen[!,:Inv_Cost_per_MWhyr][y]*vCAPENERGY[y] + dfGen[!,:Fixed_OM_Cost_per_MWhyr][y]*eTotalCapEnergy[y]
		else
			dfGen[!,:Fixed_OM_Cost_per_MWhyr][y]*eTotalCapEnergy[y]
		end
	)

	# Sum individual resource contributions to fixed costs to get total fixed costs
	@expression(EP, eTotalCFixEnergy, sum(EP[:eCFixEnergy][y] for y in STOR_ALL))

	# Add term to objective function expression
	EP[:eObj] += eTotalCFixEnergy

	### Constratints ###

	## Constraints on retirements and capacity additions
	# Cannot retire more energy capacity than existing energy capacity
	@constraint(EP, cMaxRetEnergy[y in RET_CAP_ENERGY], vRETCAPENERGY[y] <= dfGen[!,:Existing_Cap_MWh][y])

	## Constraints on new built energy capacity
	# Constraint on maximum energy capacity (if applicable) [set input to -1 if no constraint on maximum energy capacity]
	# DEV NOTE: This constraint may be violated in some cases where Existing_Cap_MWh is >= Max_Cap_MWh and lead to infeasabilty
	@constraint(EP, cMaxCapEnergy[y in intersect(dfGen[dfGen.Max_Cap_MWh.>0,:R_ID], STOR_ALL)], eTotalCapEnergy[y] <= dfGen[!,:Max_Cap_MWh][y])

	# Constraint on minimum energy capacity (if applicable) [set input to -1 if no constraint on minimum energy apacity]
	# DEV NOTE: This constraint may be violated in some cases where Existing_Cap_MWh is <= Min_Cap_MWh and lead to infeasabilty
	@constraint(EP, cMinCapEnergy[y in intersect(dfGen[dfGen.Min_Cap_MWh.>0,:R_ID], STOR_ALL)], eTotalCapEnergy[y] >= dfGen[!,:Min_Cap_MWh][y])

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
function get_beidou_orbital_elements(obj::BeidouPropagator)
    return jcall(obj, "getBeidouOrbitalElements", BeidouOrbitalElements, ())
end


text:
import CairoMakie
using Pollen
using FastAI
using FluxTraining
using DLPipelines
import DataAugmentation
using FilePathsBase
using Colors

function serve(lazy=true; kwargs...)
    refmodules = [FastAI, FluxTraining, DLPipelines, DataAugmentation, DataLoaders, FastAI.Datasets]
    project = Pollen.documentationproject(FastAI; refmodules, watchpackage=true, kwargs...)
    Pollen.serve(project, lazy=lazy)
end
serve()

##

#=
project = Pollen.documentationproject(FastAI; refmodules, inlineincludes = false, )
=#


text:
using GenomicRegion

path = "/home/gilesc/Data/transcripts/hg19/eg.bed.gz"
println (grload (path))

path = "/home/gilesc/Data/RNAseq/bg/DRR001622.bg.gz"
for item in gropen (path)
    println (item)
end

#function write_regions_by_name(outdir :: String, regions)
#    mkdir_p(outdir)
#    for region in regions
#        file = replace(region.name, [' ', '/'], '_')".bed"
#        file = uppercase(file[1:1])file[2:end]
#        path = joinpath(outdir, file)
#        open(path, "a+") do io
#            write(io, convert(String, region)"\n")
#        end
#    end
#end


text:
using DuneIstlSolvers
using Test

@testset "DuneIstlSolvers.jl" begin
    # Write your tests here.
end


text:
"""
    FEMMAcoustModule

Module for operations on interiors of domains to construct system matrices and
system vectors for linear acoustics.
"""
module FEMMAcoustModule

import Base.Complex

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, gradN!, nodesperelem, manifdim
import FinEtools.MatAcoustFluidModule: MatAcoustFluid
import FinEtools.IntegDataModule: IntegData, integrationdata, Jacobianvolume
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, SysvecAssembler, makevector!
import FinEtools.FEMMBaseModule: FEMMAbstractBase
import FinEtools.MatrixUtilityModule: add_mggt_ut_only!, add_nnt_ut_only!, complete_lt!, locjac!
import LinearAlgebra: norm

"""
    FEMMAcoust{S<:FESet}

Type for linear acoustics finite element modeling machine.
"""
mutable struct FEMMAcoust{S<:FESet, F<:Function, M} <: FEMMAbstractBase
    integdata::IntegData{S, F} # geometry data
    material::M # material object
end

function  buffers(self::FEMMAcoust, geom::NodalField{FFlt}, P::NodalField{F}) where {F}
    fes = self.integdata.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(P); # number of degrees of freedom per node
    nne =  nodesperelem(fes); # number of nodes per element
    sdim =  ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    Cedim = ndn*nne;          # dimension of the element matrix
    # Prepare assembler and temporaries
    dofnums = fill(zero(FInt), Cedim); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    gradN = fill(zero(FFlt), nne, mdim); # intermediate result -- used as a buffer
    elmat = fill(zero(FFlt), Cedim, Cedim);       # element matrix -- used as a buffer
    elvec = fill(zero(F), Cedim); # buffer
    elvecfix = fill(zero(F), Cedim); # buffer
    return dofnums, loc, J, gradN, elmat, elvec, elvecfix
end

"""
    acousticmass(self::FEMMAcoust,
      assembler::A, geom::NodalField,
      P::NodalField{T}) where {T<:Number, A<:SysmatAssemblerBase}

Compute the acoustic mass matrix.

Return K as a matrix.
Arguments
self   =  acoustics model
assembler  =  matrix assembler
geom = geometry field
P = acoustic (perturbation) pressure field
"""
function acousticmass(self::FEMMAcoust, assembler::A, geom::NodalField, P::NodalField{T}) where {T<:Number, A<:SysmatAssemblerBase}
    fes = self.integdata.fes
    dofnums, loc, J, gradN, elmat, elvec, elvecfix =   buffers(self, geom, P)
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdata);
    Jac = 0.0;
    afactor = 0.0;
    startassembly!(assembler, size(elmat,1), size(elmat,2), count(fes),
        P.nfreedofs, P.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        fill!(elmat, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
            # gradient WRT global Cartesian coordinates
            gradN!(fes, gradN, gradNparams[j], J);
            afactor = (Jac*w[j]);
            add_mggt_ut_only!(elmat, gradN, afactor)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(P, dofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function acousticmass(self::FEMMAcoust, geom::NodalField, P::NodalField{T}) where {T<:Number}
    # Make the default assembler object.
    assembler  =  SysmatAssemblerSparseSymm();
    return acousticmass(self, assembler, geom, P);
end

"""
    nzebcloadsacousticmass(self::FEMMAcoust, assembler::A,
      geom::NodalField, P::NodalField{T}) where {T<:Number,
      A<:SysvecAssemblerBase}

Compute load vector for nonzero EBC for fixed pressure..
"""
function nzebcloadsacousticmass(self::FEMMAcoust, assembler::A, geom::NodalField, P::NodalField{T}) where {T<:Number, A<:SysvecAssemblerBase}
    fes = self.integdata.fes
    dofnums, loc, J, gradN, elmat, elvec, elvecfix = buffers(self, geom, P)
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdata);
    startassembly!(assembler, P.nfreedofs);
    # Now loop over all finite elements in the set
    for i = 1:count(fes) # Loop over elements
        gatherfixedvalues_asvec!(P, elvecfix, fes.conn[i]);# retrieve element coordinates
        if norm(elvecfix, Inf) !=  0.0     # Is the load nonzero?
            fill!(elmat, 0.0);
            for j = 1:npts # Loop over quadrature points
                locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
                Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
                # gradient WRT global Cartesian coordinates
                gradN!(fes, gradN, gradNparams[j], J);
                afactor = (Jac*w[j]);
                add_mggt_ut_only!(elmat, gradN, afactor)
            end # Loop over quadrature points
            complete_lt!(elmat)
            gatherdofnums!(P, dofnums, fes.conn[i]);# retrieve degrees of freedom
            assemble!(assembler, -elmat*elvecfix, dofnums); # assemble element load vector
        end
    end
    return  makevector!(assembler);
end

function nzebcloadsacousticmass(self::FEMMAcoust, geom::NodalField, P::NodalField{T}) where {T<:Number}
    assembler  =  SysvecAssembler(P.values[1])
    return nzebcloadsacousticmass(self, assembler, geom, P);
end

"""
    acousticstiffness(self::FEMMAcoust, assembler::A,
      geom::NodalField,
      Pddot::NodalField{T}) where {T<:Number,
      A<:SysmatAssemblerBase}

Compute the acoustic stiffness matrix.
"""
function acousticstiffness(self::FEMMAcoust, assembler::A, geom::NodalField, Pddot::NodalField{T}) where {T<:Number, A<:SysmatAssemblerBase}
    fes = self.integdata.fes
    dofnums, loc, J, gradN, elmat, elvec, elvecfix = buffers(self, geom, Pddot)
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdata);
    # Material
    bulk_modulus  =   self.material.bulk_modulus;
    mass_density  =   self.material.mass_density;
    c  =  sqrt(bulk_modulus/mass_density); # sound speed
    oc2 = 1.0/c^2;
    startassembly!(assembler, size(elmat,1), size(elmat,2), count(fes),
        Pddot.nfreedofs, Pddot.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        fill!(elmat, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
            ffactor = Jac*oc2*w[j]
            add_nnt_ut_only!(elmat, Ns[j], ffactor)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(Pddot, dofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function acousticstiffness(self::FEMMAcoust, geom::NodalField, Pddot::NodalField{T}) where {T<:Number}
    # Make the default assembler object.
    assembler  =  SysmatAssemblerSparseSymm();
    return acousticstiffness(self, assembler, geom, Pddot);
end

"""
    nzebcloadsacousticstiffness(self::FEMMAcoust, assembler::A,
      geom::NodalField,
      Pddot::NodalField{T}) where {T<:Number,
      A<:SysvecAssemblerBase}

Compute load vector for nonzero EBC for fixed second-order pressure rate.
"""
function nzebcloadsacousticstiffness(self::FEMMAcoust, assembler::A, geom::NodalField, Pddot::NodalField{T}) where {T<:Number, A<:SysvecAssemblerBase}
    fes = self.integdata.fes
    dofnums, loc, J, gradN, elmat, elvec, elvecfix = buffers(self, geom, Pddot)
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdata);
    # Material
    bulk_modulus  =   self.material.bulk_modulus;
    mass_density  =   self.material.mass_density;
    c  =  sqrt(bulk_modulus/mass_density); # sound speed
    oc2 = 1.0/c^2;
    startassembly!(assembler, Pddot.nfreedofs);
    # Now loop over all finite elements in the set
    for i = 1:count(fes) # Loop over elements
        gatherfixedvalues_asvec!(Pddot, elvecfix, fes.conn[i]);# retrieve element coordinates
        if norm(elvecfix, Inf) !=  0.0  # Is the load nonzero?
            fill!(elmat, 0.0);
            for j = 1:npts # Loop over quadrature points
                locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
                Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
                ffactor = Jac*oc2*w[j]
                add_nnt_ut_only!(elmat, Ns[j], ffactor)
            end # Loop over quadrature points
            complete_lt!(elmat)
            gatherdofnums!(Pddot, dofnums, fes.conn[i]); # retrieve degrees of freedom
            assemble!(assembler, -elmat*elvecfix, dofnums); # assemble element load vector
        end
    end
    return makevector!(assembler);
end

function nzebcloadsacousticstiffness(self::FEMMAcoust, geom::NodalField, Pddot::NodalField{T}) where {T<:Number}
    assembler  =  SysvecAssembler(Pddot.values[1])#T(0.0)
    return nzebcloadsacousticstiffness(self, assembler, geom, Pddot)
end

end


text:
################################################################################
# --------------- Base type for model:
################################################################################
"""
AbstractFittedModel

Base type for fitted models.
"""
abstract type AbstractFittedModel end

################################################################################
# --------------- Base type for differentiable model:
################################################################################
abstract type AbstractDifferentiableModel <: AbstractFittedModel end

include("differentiable/julia.jl") # Julia models
include("differentiable/R.jl") # R models
include("differentiable/Python.jl") # Python models


text:
function to_rna(dna)
    transcription_dict = Dict(
        'A' => 'U', 'G' => 'C',
        'C' => 'G', 'T' => 'A'
    )
    rna = ""

    for node in dna
        try
            rna *= transcription_dict[node]
        catch e
            error("Invalid DNA node: " * node)
        end
    end

return rna

end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
"""
    rand_tangent([rng::AbstractRNG,] x)

Returns a arbitary tangent vector _appropriate_ for the primal value `x`.
Note that despite the name, no promises on the statistical randomness are made.
Rather it is an arbitary value, that is generated using the `rng`.
"""
rand_tangent(x) = rand_tangent(Random.GLOBAL_RNG, x)

rand_tangent(rng::AbstractRNG, x::Symbol) = NoTangent()
rand_tangent(rng::AbstractRNG, x::AbstractChar) = NoTangent()
rand_tangent(rng::AbstractRNG, x::AbstractString) = NoTangent()

rand_tangent(rng::AbstractRNG, x::Integer) = NoTangent()

# Try and make nice numbers with short decimal representations for good error messages
# while also not biasing the sample space too much
function rand_tangent(rng::AbstractRNG, x::T) where {T<:Number}
    # multiply by 9 to give a bigger range of values tested: no so tightly clustered around 0.
    return round(9 * randn(rng, T), digits=5, base=2)
end
rand_tangent(rng::AbstractRNG, x::Float64) = rand(rng, -9:0.01:9)
function rand_tangent(rng::AbstractRNG, x::ComplexF64)
    return ComplexF64(rand(rng, -9:0.1:9), rand(rng, -9:0.1:9))
end

#BigFloat/MPFR is finicky about short numbers, this doesn't always work as well as it should

# multiply by 9 to give a bigger range of values tested: no so tightly clustered around 0.
rand_tangent(rng::AbstractRNG, ::BigFloat) = round(big(9 * randn(rng)), digits=5, base=2)

rand_tangent(rng::AbstractRNG, x::Array{<:Any, 0}) = _compress_notangent(fill(rand_tangent(rng, x[])))
rand_tangent(rng::AbstractRNG, x::Array) = _compress_notangent(rand_tangent.(Ref(rng), x))

# All other AbstractArray's can be handled using the ProjectTo mechanics.
# and follow the same requirements
function rand_tangent(rng::AbstractRNG, x::AbstractArray)
    return _compress_notangent(ProjectTo(x)(rand_tangent(rng, collect(x))))
end

# TODO: arguably ProjectTo should handle this for us for AbstactArrays
# https://github.com/JuliaDiff/ChainRulesCore.jl/issues/410
_compress_notangent(::AbstractArray{NoTangent}) = NoTangent()
_compress_notangent(x) = x

function rand_tangent(rng::AbstractRNG, x::T) where {T}
    if !isstructtype(T)
        throw(ArgumentError("Non-struct types are not supported by this fallback."))
    end

    field_names = fieldnames(T)
    tangents = map(field_names) do field_name
        rand_tangent(rng, getfield(x, field_name))
    end
    if all(tangent isa NoTangent for tangent in tangents)
        # if none of my fields can be perturbed then I can't be perturbed
        return NoTangent()
    end

    if T <: Tuple
        return Tangent{T}(tangents...)
    else
        return Tangent{T}(; NamedTuple{field_names}(tangents)...)
    end
end

rand_tangent(rng::AbstractRNG, ::Type) = NoTangent()
rand_tangent(rng::AbstractRNG, ::Module) = NoTangent()


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using LinearAlgebra
using Random
using Distributions
using StatsPlots

# %%
A = [
    2 -1
    -1 2
]

# %%
dist = MvNormal(A)

# %%
A[1,2]/√(A[1,1]*A[2,2])

# %%
n = 2^10
XY = rand(dist, n)
X, Y = XY[1,:], XY[2,:]
@show cor(X, Y)

scatter(X, Y)

# %%
function sim(A, n, L=10^5)
    dist = MvNormal(A)
    R = zeros(L)
    XY = zeros(2, n)
    for i in eachindex(R)
        rand!(dist, XY)
        R[i] = @views cor(XY[1,:], XY[2,:])
    end
    R
end

# %%
n = 2^6
R = sim(A, n)
Z = atanh.(R)
r0 = A[1,2]/√(A[1,1]*A[2,2])
z0 = atanh(r0)
histogram(Z; norm=true, alpha=0.3, bin=200)
plot!(Normal(z0, 1/√(n - 3)))
plot!(fit_mle(Normal, Z))

# %%
n = 2^8
x, y = rand(Uniform(-1, 1), n), rand(Uniform(-1, 1), n)
X, Y = x, x + y
@show cor(X, Y)
scatter(X, Y; label="")

# %%
function sim2(n, L=10^5)
    R = zeros(L)
    d = Uniform(-1, 1)
    x = zeros(n)
    y = similar(x)
    Y = similar(x)
    for i in eachindex(R)
        rand!(d, x)
        rand!(d, y)
        @. Y = x + y
        R[i] = cor(x, Y)
    end
    R
end

# %%
n = 2^6
R = sim2(n)
Z = atanh.(R)
r0 = 1/√2
z0 = atanh(r0)
histogram(Z; norm=true, alpha=0.3, bin=200)
plot!(Normal(z0, 1/√(n - 3)))
plot!(fit_mle(Normal, Z))

# %%
n = 2^8
x, y = rand(Exponential(), n), rand(Exponential(), n)
X, Y = x, x + y
@show cor(X, Y)
scatter(X, Y; label="")

# %%
function sim3(n, L=10^5)
    R = zeros(L)
    d = Exponential()
    x = zeros(n)
    y = similar(x)
    Y = similar(x)
    for i in eachindex(R)
        rand!(d, x)
        rand!(d, y)
        @. Y = x + y
        R[i] = cor(x, Y)
    end
    R
end

# %%
n = 2^6
R = sim3(n)
Z = atanh.(R)
r0 = 1/√2
z0 = atanh(r0)
histogram(Z; norm=true, alpha=0.3, bin=200)
plot!(Normal(z0, 1/√(n - 3)))
plot!(fit_mle(Normal, Z))

# %%
function simplot(dist, n, L=10^5)
    Z = zeros(L)
    x = zeros(n)
    y = similar(x)
    Y = similar(x)
    for i in eachindex(R)
        rand!(dist, x)
        rand!(dist, y)
        @. Y = x + y
        Z[i] = atanh(cor(x, Y))
    end
    r0 = 1/√2
    z0 = atanh(r0)
    histogram(Z; norm=true, alpha=0.3, bin=200)
    plot!(Normal(z0, 1/√(n - 3)))
    plot!(fit_mle(Normal, Z))
end

# %%
simplot(Normal(), 2^6)

# %%
simplot(Uniform(), 2^6)

# %%
simplot(Exponential(), 2^6)

# %%
simplot(MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05]), 2^4)

# %%
simplot(MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05]), 2^6)

# %%

max_stars_repo_path: 0023/Fisher's Z transform.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0023/Fisher's Z transform.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0023/Fisher's Z transform.jl

max_forks_repo_name: genkuroki/public


text:
n_k = [80; 100; 120]
k = length(n_k)
n = sum(n_k)
p = 4  # make greater than or equal to k

y, X = sampledata(n_k, p)
M = MOD.classmeans(Val{:row}, X, y)
H = MOD.centerclasses!(Val{:row}, copy(X), M, y)
Σ = H'H/(n-1)
priors = ones(k)./k

info("Testing ", MOD.lda!)
for order in (:row, :col)
    for T in FloatingPointTypes
        isrowmajor = order == :row
        X_tst = isrowmajor ? copy(convert(Matrix{T}, X)) : transpose(convert(Matrix{T}, X))
        M_tst = isrowmajor ? copy(convert(Matrix{T}, M)) : transpose(convert(Matrix{T}, M))
        H_tst = isrowmajor ? copy(convert(Matrix{T}, H)) : transpose(convert(Matrix{T}, H))
        Σ_tst = isrowmajor ? copy(convert(Matrix{T}, Σ)) : transpose(convert(Matrix{T}, Σ))

        for γ in (Nullable{T}(), Nullable(zero(T)), Nullable(convert(T,0.5)), Nullable(one(T)))
            W_tmp = MOD.lda!(Val{order}, copy(X_tst), M_tst, y, γ)
            S = isnull(γ) ? Σ_tst : (1-get(γ))*Σ_tst + (get(γ)/p)*trace(Σ_tst)*I
            if isrowmajor
                @test_approx_eq (W_tmp')*S*W_tmp eye(T,p)  # Validate whitening of covariance
            else
                @test_approx_eq W_tmp*S*(W_tmp') eye(T,p)
            end
        end
        γ = Nullable{T}()

        # Check that dimension problems are caught
        X_tst1 = isrowmajor ? rand(T,1,p) : rand(T,p,1)
        X_tst2 = isrowmajor ? rand(T,n,0) : rand(T,0,n)
        @test_throws DimensionMismatch MOD.lda!(Val{order}, X_tst1, M_tst, y, γ)
        @test_throws DimensionMismatch MOD.lda!(Val{order}, X_tst2, M_tst, y, γ)
    end
end

info("Testing ", MOD.cda!)
for order in (:row, :col)
    for T in FloatingPointTypes
        isrowmajor = order == :row
        X_tst = isrowmajor ? copy(convert(Matrix{T}, X)) : transpose(convert(Matrix{T}, X))
        M_tst = isrowmajor ? copy(convert(Matrix{T}, M)) : transpose(convert(Matrix{T}, M))
        H_tst = isrowmajor ? copy(convert(Matrix{T}, H)) : transpose(convert(Matrix{T}, H))
        Σ_tst = isrowmajor ? copy(convert(Matrix{T}, Σ)) : transpose(convert(Matrix{T}, Σ))
        priors_tst = copy(convert(Vector{T}, priors))

        # Regular model tests -> Validated that covariance matrices are whitened
        for γ in (Nullable{T}(), Nullable(zero(T)), Nullable(convert(T,0.5)), Nullable(one(T)))
            S = isnull(γ) ? Σ_tst : (1-get(γ))*Σ_tst + (get(γ)/p)*trace(Σ_tst)*I
            W_tmp = MOD.cda!(Val{order}, copy(X_tst), copy(M_tst), y, γ, priors_tst)
            if isrowmajor
                @test_approx_eq (W_tmp')*S*W_tmp eye(T,k-1)
            else
                @test_approx_eq W_tmp*S*(W_tmp') eye(T,k-1)
            end
        end
        γ = Nullable{T}()

        # Check that LDA solution returned when p < k
        X_tst = isrowmajor ? copy(convert(Matrix{T}, X[:,1:k-1])) :
                             transpose(convert(Matrix{T}, X[:,1:k-1]))
        M_tst = isrowmajor ? copy(convert(Matrix{T}, M[:,1:k-1])) :
                             transpose(convert(Matrix{T}, M[:,1:k-1]))
        W_tmp = MOD.cda!(Val{order}, copy(X_tst), copy(M_tst), y, γ, priors_tst)
        W_tst = MOD.lda!(Val{order}, copy(X_tst), copy(M_tst), y, γ)
        @test_approx_eq W_tmp W_tst

        # Check that dimension problems are caught
        X_tst1 = isrowmajor ? rand(T,1,p) : rand(T,p,1)
        X_tst2 = isrowmajor ? rand(T,n,0) : rand(T,0,n)
        @test_throws DimensionMismatch MOD.cda!(Val{order}, X_tst1, M_tst, y, γ, priors_tst)
        @test_throws DimensionMismatch MOD.cda!(Val{order}, X_tst2, M_tst, y, γ, priors_tst)
        @test_throws DimensionMismatch MOD.cda!(Val{order}, X_tst,  M_tst, y, γ, rand(T,k+1))
    end
end

info("Testing ", MOD.lda)
for order in (:row, :col)
    for T in FloatingPointTypes
        isrowmajor = order == :row
        X_tst = isrowmajor ? copy(convert(Matrix{T}, X)) : transpose(convert(Matrix{T}, X))
        M_tst = isrowmajor ? copy(convert(Matrix{T}, M)) : transpose(convert(Matrix{T}, M))

        for γ in (Nullable{T}(), zero(T), Nullable(convert(T, 0.5)), convert(T, 0.5), one(T))
            model = MOD.lda(X_tst, y, order=Val{order}, gamma=γ)
            γ_tst = isa(γ, Nullable) ? γ : Nullable(γ)
            W_tst = MOD.lda!(Val{order}, copy(X_tst), copy(M_tst), y, γ_tst)
            @test model.is_cda == false
            @test_approx_eq model.W W_tst
            @test_approx_eq model.M M_tst
        end
    end
end

info("Testing ", MOD.cda)
for order in (:row, :col)
    for T in FloatingPointTypes
        isrowmajor = order == :row
        X_tst = isrowmajor ? copy(convert(Matrix{T}, X)) : transpose(convert(Matrix{T}, X))
        M_tst = isrowmajor ? copy(convert(Matrix{T}, M)) : transpose(convert(Matrix{T}, M))
        priors_tst = copy(convert(Vector{T}, priors))

        for γ in (Nullable{T}(), zero(T), Nullable(convert(T, 0.5)), convert(T, 0.5), one(T))
            model = MOD.cda(X_tst, y, order=Val{order}, gamma=γ)
            γ_tst = isa(γ, Nullable) ? γ : Nullable(γ)
            W_tst = MOD.cda!(Val{order}, copy(X_tst), copy(M_tst), y, γ_tst, priors_tst)
            @test model.is_cda == true
            @test_approx_eq model.W W_tst
            @test_approx_eq model.M M_tst
        end
    end
end

info("Testing ", MOD.discriminants_lda)
for T in FloatingPointTypes
    X_tmp = copy(convert(Matrix{T}, X))
    priors_tmp = convert(Vector{T}, priors)

    model1 = lda(X_tmp,  y, order=Val{:row})
    model2 = lda(X_tmp', y, order=Val{:col})
    D = hcat(-[MOD.dotvectors(Val{:row}, (X_tmp .- M[i,:])*model1.W)/2
                for i in eachindex(priors_tmp)]...) .+ log(priors_tmp)'

    @test_approx_eq D  MOD.discriminants(model1, X_tmp)
    @test_approx_eq D' MOD.discriminants(model2, X_tmp')
end

info("Testing ", MOD.classify)
for order in (:row, :col)
    for T in FloatingPointTypes
        isrowmajor = order == :row
        X_tst = isrowmajor ? copy(convert(Matrix{T}, X)) : transpose(convert(Matrix{T}, X))
        priors_tst = convert(Vector{T}, priors)

        model1 = lda(X_tst, y, order=Val{order})
        model2 = cda(X_tst, y, order=Val{order})

        D1 = vec(mapslices(indmax, MOD.discriminants(model1, X_tst), isrowmajor ? 2 : 1))
        D2 = vec(mapslices(indmax, MOD.discriminants(model2, X_tst), isrowmajor ? 2 : 1))

        @test D1 == MOD.classify(model1, X_tst)
        @test D2 == MOD.classify(model2, X_tst)
    end
end

info("Testing ", MOD.ModelLDA)
show(DevNull, lda(X, y, order=Val{:row}))


text:
"""
    area(x, y)

Creates an `EChart` where region below plotted line filled with color.

## Methods
```julia
area(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}})
area(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2})
area(df::AbstractDataFrame, x::Symbol, y::Symbol)
area(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol)
area(k::KernelDensity.UnivariateKDE)
```

## Arguments
* `mark::Union{String, AbstractVector} = "line"` : how to display plotted points
* `fill::Union{Bool, AbstractVector} = true` : fill area below marks with color?
* `stack::Union{Bool, AbstractVector, Void} = nothing` : stack (add together) when multple series present?
* `step::Union{String, Void} = nothing` : one of {"start", "end", "middle", nothing}
* `legend::Bool` : display legend?
* `scale::Bool = false` : show full Y-axis or truncated
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `area`, such as displaying a legend when two or more series present.

## Examples
```julia
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y)
```
"""
function area(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}};
			mark::Union{String, AbstractVector} = "line",
			fill::Union{Bool, AbstractVector} = true,
			stack::Union{Bool, AbstractVector, Nothing} = nothing,
			step::Union{String, Nothing} = nothing,
			legend::Bool = false,
			scale::Bool = false,
			kwargs...)

	ec = xy_plot(x, y; stack = stack, mark = mark, step = step, legend = legend, scale = scale, kwargs...)
	ec.xAxis[1].boundaryGap = false

	# Fill area if requested
	fill!(ec, 1, fill)

	return ec

end

#Basically same function as above, except automatically puts a legend when multiple series present
function area(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2};
			mark::Union{String, AbstractVector} = "line",
			fill::Union{Bool, AbstractVector} = true,
			stack::Union{Bool, AbstractVector, Nothing} = true,
			step::Union{String, Nothing} = nothing,
			legend::Bool = true,
			scale::Bool = false,
			kwargs...)

	ec = xy_plot(x, y; stack = stack, mark = mark, step = step, legend = legend, scale = scale, kwargs...)
	ec.xAxis[1].boundaryGap = false

	# Fill area if requested
	fill!(ec, size(y)[2], fill)

	return ec

end

#dataframe, single series
function area(df::AbstractDataFrame, x::Symbol, y::Symbol;
			mark::Union{String, AbstractVector} = "line",
			fill::Union{Bool, AbstractVector} = true,
			stack::Union{Bool, AbstractVector, Nothing} = nothing,
			step::Union{String, Nothing} = nothing,
			legend::Bool = false,
			scale::Bool = false,
			kwargs...)

			ec = xy_plot(df, x, y, mark = mark, stack = stack, step = step, legend = legend, scale = scale, kwargs...)

			ec.xAxis[1].boundaryGap = false

			fill!(ec, 1, fill)

			return ec

end

#dataframe, group argument
function area(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol;
			mark::Union{String, AbstractVector} = "line",
			fill::Union{Bool, AbstractVector} = true,
			stack::Union{Bool, AbstractVector, Nothing} = true,
			step::Union{String, Nothing} = nothing,
			legend::Bool = true,
			scale::Bool = false,
			kwargs...)

			ec = xy_plot(df, x, y, group, mark = mark, stack = stack, step = step, legend = legend, scale = scale, kwargs...)

			ec.xAxis[1].boundaryGap = false

			fill!(ec, length(ec.series), fill)

			return ec

end

function area(k::KernelDensity.UnivariateKDE;
				mark::Union{String, AbstractVector} = "line",
				fill::Union{Bool, AbstractVector} = true,
				stack::Union{Bool, AbstractVector, Nothing} = nothing,
				step::Union{String, Nothing} = nothing,
				legend::Bool = false,
				scale::Bool = false,
				kwargs...)

	ec = area(collect(k.x), k.density, mark = mark, fill = fill, stack = stack, step = step, legend = legend, scale = scale, kwargs...)

	return ec

end


text:
# Demo test functions
using Dates: Date
using Lazy

# Account (child object)
function test_account()
    acct = Account("1234", 100.00, Date(2019, 1, 1))
    @show acct
    @show deposit!(acct, 25)
    println("-------------------------------------")

    dest = Account("4321", 500.00, Date(2019, 2, 1))
    @show dest
    @show withdraw!(dest, 50.00)
    println("-------------------------------------")

    @show transfer!(acct, dest, 10.00)
    @show acct
    @show dest

    return nothing
end

# SavingsAccount (parent object - manual forwarding)
function test_savings_account()
    acct = SavingsAccount("1234", 3000.00, Date(2019, 3, 1), 0.03)
    acct2 = SavingsAccount("5678", 1000.00, Date(2018, 3, 1), 0.03)
    test_savings_account_logic(acct, acct2)
end

# SavingsAccount2 (parent object - using @forward)
function test_savings_account2()
    acct = SavingsAccount2("1234", 3000.00, Date(2019, 3, 1), 0.03)
    acct2 = SavingsAccount2("5678", 1000.00, Date(2018, 3, 1), 0.03)
    test_savings_account_logic(acct, acct2)
end

# common test logic
function test_savings_account_logic(acct, acct2)
    @show acct acct2
    @show deposit!(acct, 200)
    @show withdraw!(acct2, 30)
    println("-------------------------------------")
    @show acct acct2
    @show transfer!(acct, acct2, 500)
    @show acct acct2
    println("-------------------------------------")
    @show accrue_daily_interest!(acct)
    @show accrue_daily_interest!(acct)
    return nothing
end

# Check how @forward expands code
function check_forward_macro()
    @macroexpand(@forward SavingsAccount.acct balance,deposit!)
end