text:
function destroy(layer::AbstractFeatureLayer)::Nothing
    layer.ptr = C_NULL
    return nothing
end

function destroy(layer::IFeatureLayer)::Nothing
    layer.ptr = C_NULL
    layer.ownedby = Dataset()
    layer.spatialref = SpatialRef()
    return nothing
end

"""
    createlayer(name, dataset, geom, spatialref, options)

This function attempts to create a new layer on the dataset with the indicated
`name`, `spatialref`, and geometry type.

### Keyword Arguments
* `name`: the name for the new layer. This should ideally not match any
    existing layer on the datasource. Defaults to an empty string.
* `dataset`: the dataset. Defaults to creating a new in memory dataset.
* `geom`: the geometry type for the layer. Use `wkbUnknown` (default) if
    there are no constraints on the types geometry to be written.
* `spatialref`: the coordinate system to use for the new layer.
* `options`: a StringList of name=value (driver-specific) options.
"""
function createlayer(;
    name::AbstractString = "",
    dataset::AbstractDataset = create(getdriver("Memory")),
    geom::OGRwkbGeometryType = wkbUnknown,
    spatialref::AbstractSpatialRef = SpatialRef(),
    options = StringList(C_NULL),
)::IFeatureLayer
    return IFeatureLayer(
        GDAL.gdaldatasetcreatelayer(
            dataset.ptr,
            name,
            spatialref.ptr,
            geom,
            options,
        ),
        ownedby = dataset,
        spatialref = spatialref,
    )
end

function unsafe_createlayer(;
    name::AbstractString = "",
    dataset::AbstractDataset = create(getdriver("Memory")),
    geom::OGRwkbGeometryType = wkbUnknown,
    spatialref::AbstractSpatialRef = SpatialRef(),
    options = StringList(C_NULL),
)::FeatureLayer
    return FeatureLayer(
        GDAL.gdaldatasetcreatelayer(
            dataset.ptr,
            name,
            spatialref.ptr,
            geom,
            options,
        ),
        ownedby = dataset,
        spatialref = spatialref,
    )
end

"""
    copy(layer, dataset, name, options)

Copy an existing layer.

This method creates a new layer, duplicate the field definitions of the source
layer, and then duplicates each feature of the source layer.

### Parameters
* `layer`: source layer to be copied.

### Keyword Arguments
* `dataset`: the dataset handle. (Creates a new dataset in memory by default.)
* `name`: the name of the layer to create on the dataset.
* `options`: a StringList of name=value (driver-specific) options.
"""
function copy(
    layer::AbstractFeatureLayer;
    dataset::AbstractDataset = create(getdriver("Memory")),
    name::AbstractString = "copy($(getname(layer)))",
    options = StringList(C_NULL),
)::IFeatureLayer
    return IFeatureLayer(
        GDAL.gdaldatasetcopylayer(dataset.ptr, layer.ptr, name, options),
        ownedby = dataset,
    )
end

function unsafe_copy(
    layer::AbstractFeatureLayer;
    dataset::AbstractDataset = create(getdriver("Memory")),
    name::AbstractString = "copy($(getname(layer)))",
    options = StringList(C_NULL),
)::FeatureLayer
    return FeatureLayer(
        GDAL.gdaldatasetcopylayer(dataset.ptr, layer.ptr, name, options),
    )
end

"""
    getname(layer::AbstractFeatureLayer)

Return the layer name.
"""
getname(layer::AbstractFeatureLayer)::String = GDAL.ogr_l_getname(layer.ptr)

"""
    getgeomtype(layer::AbstractFeatureLayer)

Return the layer geometry type.
"""
getgeomtype(layer::AbstractFeatureLayer)::OGRwkbGeometryType =
    GDAL.ogr_l_getgeomtype(layer.ptr)

"""
    getspatialfilter(layer::AbstractFeatureLayer)

Returns the current spatial filter for this layer.
"""
function getspatialfilter(layer::AbstractFeatureLayer)::IGeometry
    result = GDAL.ogr_l_getspatialfilter(Ptr{Cvoid}(layer.ptr))
    return if result == C_NULL
        IGeometry(result)
    else
        # NOTE(yeesian): we make a clone here so that the geometry does not
        # depend on the FeatureLayer.
        IGeometry(GDAL.ogr_g_clone(result))
    end
end

"""
    getspatialref(layer::AbstractFeatureLayer)

Returns a clone of the spatial reference system for this layer.
"""
function getspatialref(layer::AbstractFeatureLayer)::ISpatialRef
    result = GDAL.ogr_l_getspatialref(layer.ptr)
    return if result == C_NULL
        ISpatialRef()
    else
        # NOTE(yeesian): we make a clone here so that the spatialref does not
        # depend on the FeatureLayer/Dataset.
        ISpatialRef(GDAL.osrclone(result))
    end
end

function unsafe_getspatialref(layer::AbstractFeatureLayer)::SpatialRef
    result = GDAL.ogr_l_getspatialref(layer.ptr)
    return if result == C_NULL
        SpatialRef()
    else
        # NOTE(yeesian): we make a clone here so that the spatialref does not
        # depend on the FeatureLayer/Dataset.
        SpatialRef(GDAL.osrclone(result))
    end
end

"""
    setspatialfilter!(layer::AbstractFeatureLayer, geom::Geometry)

Set a new spatial filter for the layer, using the geom.

This method set the geometry to be used as a spatial filter when fetching
features via the `nextfeature()` method. Only features that geometrically
intersect the filter geometry will be returned.

### Parameters
* `layer`  handle to the layer on which to set the spatial filter.
* `geom`   handle to the geometry to use as a filtering region. NULL may be
           passed indicating that the current spatial filter should be cleared,
           but no new one instituted.

### Remarks
Currently this test may be inaccurately implemented, but it is guaranteed
that all features whose envelope (as returned by OGRGeometry::getEnvelope())
overlaps the envelope of the spatial filter will be returned. This can result
in more shapes being returned that should strictly be the case.

For the time being the passed filter geometry should be in the same SRS as the
geometry field definition it corresponds to (as returned by
`GetLayerDefn()->OGRFeatureDefn::GetGeomFieldDefn(i)->GetSpatialRef()`).
In the future this may be generalized.

Note that only the last spatial filter set is applied, even if several
successive calls are done with different iGeomField values.
"""
function setspatialfilter!(
    layer::L,
    geom::Geometry,
)::L where {L<:AbstractFeatureLayer}
    # This method makes an internal copy of `geom`. The input `geom` remains
    # the responsibility of the caller, and may be safely destroyed.
    GDAL.ogr_l_setspatialfilter(layer.ptr, geom.ptr)
    return layer
end

function clearspatialfilter!(layer::L)::L where {L<:AbstractFeatureLayer}
    GDAL.ogr_l_setspatialfilter(layer.ptr, C_NULL)
    return layer
end

"""
    setspatialfilter!(layer::AbstractFeatureLayer, xmin, ymin, xmax, ymax)

Set a new rectangular spatial filter for the layer.

This method set rectangle to be used as a spatial filter when fetching features
via the `nextfeature()` method. Only features that geometrically intersect the
given rectangle will be returned.

The x/y values should be in the same coordinate system as the layer as a whole
(as returned by OGRLayer::GetSpatialRef()). Internally this method is normally
implemented as creating a 5 vertex closed rectangular polygon and passing it to
OGRLayer::SetSpatialFilter(). It exists as a convenience.

The only way to clear a spatial filter set with this method is to call
`OGRLayer::SetSpatialFilter(NULL)`.
"""
function setspatialfilter!(
    layer::L,
    xmin::Real,
    ymin::Real,
    xmax::Real,
    ymax::Real,
)::L where {L<:AbstractFeatureLayer}
    GDAL.ogr_l_setspatialfilterrect(layer.ptr, xmin, ymin, xmax, ymax)
    return layer
end

"""
    setspatialfilter!(layer::AbstractFeatureLayer, i::Integer,
        geom::AbstractGeometry)

Set a new spatial filter.

This method set the geometry to be used as a spatial filter when fetching
features via the `nextfeature()` method. Only features that geometrically
intersect the filter geometry will be returned.

### Parameters
* `layer`: the layer on which to set the spatial filter.
* `i`: index of the geometry field on which the spatial filter operates.
* `geom`: the geometry to use as a filtering region. NULL may be passed
    indicating that the current spatial filter should be cleared, but
    no new one instituted.

### Remarks
Currently this test is may be inaccurately implemented, but it is guaranteed
that all features who's envelope (as returned by OGRGeometry::getEnvelope())
overlaps the envelope of the spatial filter will be returned. This can result
in more shapes being returned that should strictly be the case.

For the time being the passed filter geometry should be in the same SRS as the
layer (as returned by OGRLayer::GetSpatialRef()). In the future this may be
generalized.
"""
function setspatialfilter!(
    layer::L,
    i::Integer,
    geom::AbstractGeometry,
)::L where {L<:AbstractFeatureLayer}
    # This method makes an internal copy of `geom`. The input `geom` remains
    # the responsibility of the caller, and may be safely destroyed.
    GDAL.ogr_l_setspatialfilterex(layer.ptr, i, geom.ptr)
    return layer
end

function clearspatialfilter!(
    layer::L,
    i::Integer,
)::L where {L<:AbstractFeatureLayer}
    GDAL.ogr_l_setspatialfilterex(layer.ptr, i, C_NULL)
    return layer
end

"""
    setspatialfilter!(layer::AbstractFeatureLayer, i::Integer, xmin, ymin, xmax,
        ymax)

Set a new rectangular spatial filter.

### Parameters
* `layer`: the feature layer on which to set the spatial filter.
* `i`: index of the geometry field on which the spatial filter operates.
* `xmin`: the minimum X coordinate for the rectangular region.
* `ymin`: the minimum Y coordinate for the rectangular region.
* `xmax`: the maximum X coordinate for the rectangular region.
* `ymax`: the maximum Y coordinate for the rectangular region.
"""
function setspatialfilter!(
    layer::L,
    i::Integer,
    xmin::Real,
    ymin::Real,
    xmax::Real,
    ymax::Real,
)::L where {L<:AbstractFeatureLayer}
    GDAL.ogr_l_setspatialfilterrectex(layer.ptr, i, xmin, ymin, xmax, ymax)
    return layer
end

"""
    setattributefilter!(layer::AbstractFeatureLayer, query::AbstractString)

Set a new attribute query.

This method sets the attribute query string to be used when fetching features
via the `nextfeature()` method. Only features for which the query evaluates as
`true` will be returned.

### Parameters
* `layer`: handle to the layer on which attribute query will be executed.
* `query`: query in restricted SQL WHERE format.

### Remarks
The query string should be in the format of an SQL WHERE clause. For instance
`"population > 1000000 and population < 5000000"` where population is an
attribute in the layer. The query format is normally a restricted form of
SQL WHERE clause as described in the "WHERE" section of the OGR SQL tutorial.
In some cases (RDBMS backed drivers) the native capabilities of the database may
be used to interpret the WHERE clause in which case the capabilities will be
broader than those of OGR SQL.

Note that installing a query string will generally result in resetting the
current reading position (ala `resetreading!()`).
"""
function setattributefilter!(
    layer::L,
    query::AbstractString,
)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_setattributefilter(layer.ptr, query)
    @ogrerr result """Failed to set a new attribute query. The query expression
    might be in error."""
    return layer
end

function clearattributefilter!(layer::L)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_setattributefilter(layer.ptr, C_NULL)
    @ogrerr result "OGRErr $result: Failed to clear attribute query."
    return layer
end

"""
    resetreading!(layer::AbstractFeatureLayer)

Reset feature reading to start on the first feature.

This affects `nextfeature()`.
"""
function resetreading!(layer::L)::L where {L<:AbstractFeatureLayer}
    GDAL.ogr_l_resetreading(layer.ptr)
    return layer
end

"""
    unsafe_nextfeature(layer::AbstractFeatureLayer)

Fetch the next available feature from this layer.

### Parameters
* `layer`: the feature layer to be read from.

### Remarks
This method implements sequential access to the features of a layer. The
`resetreading!()` method can be used to start at the beginning again. Only
features matching the current spatial filter (set with `setspatialfilter!()`)
will be returned.

The returned feature becomes the responsibility of the caller to delete with
`destroy()`. It is critical that all features associated with a `FeatureLayer`
(more specifically a `FeatureDefn`) be destroyed before that layer is destroyed.

Features returned by `nextfeature()` may or may not be affected by concurrent
modifications depending on drivers. A guaranteed way of seeing modifications in
effect is to call `resetreading!()` on layers where `nextfeature()` has been
called, before reading again. Structural changes in layers (field addition,
deletion, ...) when a read is in progress may or may not be possible depending
on drivers. If a transaction is committed/aborted, the current sequential
reading may or may not be valid after that operation and a call to
`resetreading!()` might be needed.
"""
function unsafe_nextfeature(layer::AbstractFeatureLayer)::Feature
    return Feature(GDAL.ogr_l_getnextfeature(layer.ptr))
end

"""
    setnextbyindex!(layer::AbstractFeatureLayer, i::Integer)

Move read cursor to the `i`-th feature in the current resultset.

This method allows positioning of a layer such that the `nextfeature()` call
will read the requested feature, where `i` is an absolute index into the
current result set. So, setting it to `3` would mean the next feature read with
`nextfeature()` would have been the fourth feature to have been read if
sequential reading took place from the beginning of the layer, including
accounting for spatial and attribute filters.

### Parameters
* `layer`: handle to the layer
* `i`: the index indicating how many steps into the result set to seek.

### Remarks
Only in rare circumstances is `setnextbyindex!()` efficiently implemented. In
all other cases the default implementation which calls `resetreading!()` and
then calls `nextfeature()` `i` times is used. To determine if fast seeking is
available on the layer, use the `testcapability()` method with a value of
`OLCFastSetNextByIndex`.
"""
function setnextbyindex!(
    layer::L,
    i::Integer,
)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_setnextbyindex(layer.ptr, i)
    @ogrerr result "Failed to move the cursor to index $i"
    return layer
end

"""
    unsafe_getfeature(layer::AbstractFeatureLayer, i::Integer)

Return a feature (now owned by the caller) by its identifier or NULL on failure.

### Parameters
* `layer`: the feature layer to be read from.
* `i`: the index of the feature to be returned.

### Remarks
This function will attempt to read the identified feature. The nFID value cannot
be OGRNullFID. Success or failure of this operation is unaffected by the spatial
or attribute filters (and specialized implementations in drivers should make
sure that they do not take into account spatial or attribute filters).

If this function returns a non-NULL feature, it is guaranteed that its feature
id (OGR_F_GetFID()) will be the same as nFID.

Use OGR_L_TestCapability(OLCRandomRead) to establish if this layer supports
efficient random access reading via OGR_L_GetFeature(); however, the call should
always work if the feature exists as a fallback implementation just scans all
the features in the layer looking for the desired feature.

Sequential reads (with OGR_L_GetNextFeature()) are generally considered
interrupted by a OGR_L_GetFeature() call.

The returned feature is now owned by the caller, and should be freed with
`destroy()`.
"""
unsafe_getfeature(layer::AbstractFeatureLayer, i::Integer)::Feature =
    Feature(GDAL.ogr_l_getfeature(layer.ptr, i))

"""
    setfeature!(layer::AbstractFeatureLayer, feature::Feature)

Rewrite an existing feature.

This function will write a feature to the layer, based on the feature id within
the OGRFeature.

### Remarks
Use OGR_L_TestCapability(OLCRandomWrite) to establish if this layer supports
random access writing via OGR_L_SetFeature().
"""
function setfeature!(layer::AbstractFeatureLayer, feature::Feature)
    result = GDAL.ogr_l_setfeature(layer.ptr, feature.ptr)
    # OGRERR_NONE if the operation works, otherwise an appropriate error code
    # (e.g OGRERR_NON_EXISTING_FEATURE if the feature does not exist).
    @ogrerr result "Failed to set feature."
    return layer
end

"""
    addfeature!(layer::AbstractFeatureLayer, feature::Feature)

Write a new feature within a layer.

### Remarks
The passed feature is written to the layer as a new feature, rather than
overwriting an existing one. If the feature has a feature id other than
OGRNullFID, then the native implementation may use that as the feature id of
the new feature, but not necessarily. Upon successful return the passed feature
will have been updated with the new feature id.
"""
function addfeature!(
    layer::L,
    feature::Feature,
)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_createfeature(layer.ptr, feature.ptr)
    @ogrerr result "Failed to create and write feature in layer."
    return layer
end

function addfeature(f::Function, layer::L)::L where {L<:AbstractFeatureLayer}
    feature = unsafe_createfeature(layer)
    return try
        f(feature)
        addfeature!(layer, feature)
    finally
        destroy(feature)
    end
end

"""
    unsafe_createfeature(layer::AbstractFeatureLayer)

Create and returns a new feature based on the layer definition.

The newly feature is owned by the layer (it will increase the number of features
the layer by one), but the feature has not been written to the layer yet.
"""
unsafe_createfeature(layer::AbstractFeatureLayer)::Feature =
    unsafe_createfeature(layerdefn(layer))

function createfeature(f::Function, layer::L)::L where {L<:AbstractFeatureLayer}
    feature = unsafe_createfeature(layer)
    return try
        f(feature)
        setfeature!(layer, feature)
    finally
        destroy(feature)
    end
end

"""
    deletefeature!(layer::AbstractFeatureLayer, i::Integer)

Delete feature with fid `i` from layer.

### Remarks
The feature with the indicated feature id is deleted from the layer if supported
by the driver. Most drivers do not support feature deletion, and will return
OGRERR_UNSUPPORTED_OPERATION. The OGR_L_TestCapability() function may be called
with OLCDeleteFeature to check if the driver supports feature deletion.
"""
function deletefeature!(layer::L, i::Integer)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_deletefeature(layer.ptr, i)
    # OGRERR_NONE if the operation works, otherwise an appropriate error code
    # (e.g OGRERR_NON_EXISTING_FEATURE if the feature does not exist).
    @ogrerr result "OGRErr $result: Failed to delete feature $i"
    return layer
end

"""
    layerdefn(layer::AbstractFeatureLayer)

Returns a view of the schema information for this layer.

### Remarks
The `featuredefn` is owned by the `layer` and should not be modified.
"""
layerdefn(layer::AbstractFeatureLayer)::IFeatureDefnView =
    IFeatureDefnView(GDAL.ogr_l_getlayerdefn(layer.ptr))

"""
    findfieldindex(layer::AbstractFeatureLayer,
        field::Union{AbstractString, Symbol}, exactmatch::Bool)

Find the index of the field in a layer, or -1 if the field doesn't exist.

If `exactmatch` is set to `false` and the field doesn't exists in the given form
the driver might apply some changes to make it match, like those it might do if
the layer was created (eg. like `LAUNDER` in the OCI driver).
"""
function findfieldindex(
    layer::AbstractFeatureLayer,
    field::Union{AbstractString,Symbol},
    exactmatch::Bool,
)::Integer
    return GDAL.ogr_l_findfieldindex(layer.ptr, field, exactmatch)
end

"""
    nfeature(layer::AbstractFeatureLayer, force::Bool = false)

Fetch the feature count in this layer, or `-1` if the count is not known.

### Parameters
* `layer`: handle to the layer that owned the features.
* `force`: flag indicating whether the count should be computed even if it is
    expensive. (`false` by default.)
"""
nfeature(layer::AbstractFeatureLayer, force::Bool = false)::Integer =
    GDAL.ogr_l_getfeaturecount(layer.ptr, force)

"""
    ngeom(layer::AbstractFeatureLayer)

Fetch number of geometry fields on the feature layer.
"""
ngeom(layer::AbstractFeatureLayer)::Integer = ngeom(layerdefn(layer))

"""
    nfield(layer::AbstractFeatureLayer)

Fetch number of fields on the feature layer.
"""
nfield(layer::AbstractFeatureLayer)::Integer = nfield(layerdefn(layer))

"""
    envelope(layer::AbstractFeatureLayer, force::Bool = false)
    envelope(layer::AbstractFeatureLayer, i::Integer, force::Bool = false)

Fetch the extent of this layer.

Returns the extent (MBR) of the data in the layer. If `force` is `false`, and it
would be expensive to establish the extent then OGRERR_FAILURE will be returned
indicating that the extent isn't know. If `force` is `true` then some
implementations will actually scan the entire layer once to compute the MBR of
all the features in the layer.

### Parameters
* `layer`: handle to the layer from which to get extent.
* `i`:     (optional) the index of the geometry field to compute the extent.
* `force`: Flag indicating whether the extent should be computed even if it is
            expensive.

### Additional Remarks
Depending on the drivers, the returned extent may or may not take the spatial
filter into account. So it is safer to call GetExtent() without setting a
spatial filter.

Layers without any geometry may return OGRERR_FAILURE just indicating that no
meaningful extents could be collected.

Note that some implementations of this method may alter the read cursor of the
layer.
"""
function envelope(
    layer::AbstractFeatureLayer,
    i::Integer,
    force::Bool = false,
)::GDAL.OGREnvelope
    envelope = Ref{GDAL.OGREnvelope}(GDAL.OGREnvelope(0, 0, 0, 0))
    result = GDAL.ogr_l_getextentex(layer.ptr, i, envelope, force)
    @ogrerr result "Extent not known"
    return envelope[]
end

function envelope(
    layer::AbstractFeatureLayer,
    force::Bool = false,
)::GDAL.OGREnvelope
    envelope = Ref{GDAL.OGREnvelope}(GDAL.OGREnvelope(0, 0, 0, 0))
    result = GDAL.ogr_l_getextent(layer.ptr, envelope, force)
    @ogrerr result "Extent not known"
    return envelope[]
end

"""
    testcapability(layer::AbstractFeatureLayer, capability::AbstractString)

Test if this layer supported the named capability.

### Parameters
* `capability`  the name of the capability to test.

### Returns
`true` if the layer has the requested capability, `false` otherwise. It will
return `false` for any unrecognized capabilities.

### Additional Remarks
The capability codes that can be tested are represented as strings, but
`#defined` constants exists to ensure correct spelling. Specific layer types may
implement class specific capabilities, but this can't generally be discovered by
the caller.

* `OLCRandomRead` / \"RandomRead\": `true` if the GetFeature() method is
    implemented in an optimized way for this layer, as opposed to the default
    implementation using `resetreading!()` and `nextfeature()` to find the
    requested feature id.

* `OLCSequentialWrite` / \"SequentialWrite\": `true` if the CreateFeature()
    method works for this layer. Note this means that this particular layer is
    writable. The same OGRLayer class may returned `false` for other layer
    instances that are effectively read-only.

* `OLCRandomWrite` / \"RandomWrite\": `true` if the SetFeature() method is
    operational on this layer. Note this means that this particular layer is
    writable. The same OGRLayer class may returned `false` for other layer
    instances that are effectively read-only.

* `OLCFastSpatialFilter` / \"FastSpatialFilter\": `true` if this layer
    implements spatial filtering efficiently. Layers that effectively read all
    features, and test them with the OGRFeature intersection methods should
    return `false`. This can be used as a clue by the application whether it
    should build and maintain its own spatial index for features in this layer.

* `OLCFastFeatureCount` / \"FastFeatureCount\": `true` if this layer can return
    a feature count (via GetFeatureCount()) efficiently. i.e. without counting
    the features. In some cases this will return `true` until a spatial filter
    is installed after which it will return `false`.

* `OLCFastGetExtent` / \"FastGetExtent\": `true` if this layer can return its
    data extent (via GetExtent()) efficiently, i.e. without scanning all the
    features. In some cases this will return `true` until a spatial filter is
    installed after which it will return `false`.

* `OLCFastSetNextByIndex` / \"FastSetNextByIndex\": `true` if this layer can
    perform the SetNextByIndex() call efficiently, otherwise `false`.

* `OLCCreateField` / \"CreateField\": `true` if this layer can create new fields
    on the current layer using CreateField(), otherwise `false`.

* `OLCCreateGeomField` / \"CreateGeomField\": (GDAL >= 1.11) `true` if this
    layer can create new geometry fields on the current layer using
    CreateGeomField(), otherwise `false`.

* `OLCDeleteField` / \"DeleteField\": `true` if this layer can delete existing
    fields on the current layer using DeleteField(), otherwise `false`.

* `OLCReorderFields` / \"ReorderFields\": `true` if this layer can reorder
    existing fields on the current layer using ReorderField() or
    ReorderFields(), otherwise `false`.

* `OLCAlterFieldDefn` / \"AlterFieldDefn\": `true` if this layer can alter the
    definition of an existing field on the current layer using AlterFieldDefn(),
    otherwise `false`.

* `OLCDeleteFeature` / \"DeleteFeature\": `true` if the DeleteFeature() method
    is supported on this layer, otherwise `false`.

* `OLCStringsAsUTF8` / \"StringsAsUTF8\": `true` if values of OFTString fields
    are assured to be in UTF-8 format. If `false` the encoding of fields is
    uncertain, though it might still be UTF-8.

* `OLCTransactions` / \"Transactions\": `true` if the StartTransaction(),
    CommitTransaction() and RollbackTransaction() methods work in a meaningful
    way, otherwise `false`.

* `OLCIgnoreFields` / \"IgnoreFields\": `true` if fields, geometry and style
    will be omitted when fetching features as set by SetIgnoredFields() method.

* `OLCCurveGeometries` / \"CurveGeometries\": `true` if this layer supports
    writing curve geometries or may return such geometries. (GDAL 2.0).
"""
testcapability(layer::AbstractFeatureLayer, capability::AbstractString)::Bool =
    Bool(GDAL.ogr_l_testcapability(layer.ptr, capability))

function listcapability(
    layer::AbstractFeatureLayer,
    capabilities = (
        GDAL.OLCRandomRead,
        GDAL.OLCSequentialWrite,
        GDAL.OLCRandomWrite,
        GDAL.OLCFastSpatialFilter,
        GDAL.OLCFastFeatureCount,
        GDAL.OLCFastGetExtent,
        GDAL.OLCCreateField,
        GDAL.OLCDeleteField,
        GDAL.OLCReorderFields,
        GDAL.OLCAlterFieldDefn,
        GDAL.OLCTransactions,
        GDAL.OLCDeleteFeature,
        GDAL.OLCFastSetNextByIndex,
        GDAL.OLCStringsAsUTF8,
        GDAL.OLCIgnoreFields,
        GDAL.OLCCreateGeomField,
        GDAL.OLCCurveGeometries,
        GDAL.OLCMeasuredGeometries,
    ),
)::Dict{String,Bool}
    return Dict{String,Bool}([
        c => testcapability(layer, c) for c in capabilities
    ],)
end

# TODO use syntax below once v0.4 support is dropped (not in Compat.jl)
# listcapability(layer::AbstractFeatureLayer) = Dict(
#     c => testcapability(layer,c) for c in
#     (GDAL.OLCRandomRead,        GDAL.OLCSequentialWrite, GDAL.OLCRandomWrite,
#      GDAL.OLCFastSpatialFilter, GDAL.OLCFastFeatureCount,GDAL.OLCFastGetExtent,
#      GDAL.OLCCreateField,       GDAL.OLCDeleteField,     GDAL.OLCReorderFields,
#      GDAL.OLCAlterFieldDefn,    GDAL.OLCTransactions,    GDAL.OLCDeleteFeature,
#      GDAL.OLCFastSetNextByIndex,GDAL.OLCStringsAsUTF8,   GDAL.OLCIgnoreFields,
#      GDAL.OLCCreateGeomField,   GDAL.OLCCurveGeometries,
#      GDAL.OLCMeasuredGeometries)
# )

"""
    addfielddefn!(layer::AbstractFeatureLayer, field::AbstractFieldDefn,
        approx = false)

Create a new field on a layer.

### Parameters
* `layer`:  the layer to write the field definition.
* `field`:  the field definition to write to disk.
* `approx`: If `true`, the field may be created in a slightly different form
            depending on the limitations of the format driver.

### Remarks
You must use this to create new fields on a real layer. Internally the
OGRFeatureDefn for the layer will be updated to reflect the new field.
Applications should never modify the OGRFeatureDefn used by a layer directly.

This function should not be called while there are feature objects in existence
that were obtained or created with the previous layer definition.

Not all drivers support this function. You can query a layer to check if it
supports it with the `GDAL.OLCCreateField` capability. Some drivers may only
support this method while there are still no features in the layer. When it is
supported, the existing features of the backing file/database should be updated
accordingly.

Drivers may or may not support not-null constraints. If they support creating
fields with not-null constraints, this is generally before creating any feature
to the layer.
"""
function addfielddefn!(
    layer::L,
    field::AbstractFieldDefn,
    approx::Bool = false,
)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_createfield(layer.ptr, field.ptr, approx)
    @ogrerr result "Failed to create new field"
    return layer
end

"""
    addgeomdefn!(layer::AbstractFeatureLayer, field::AbstractGeomFieldDefn,
        approx = false)

Create a new geometry field on a layer.

### Parameters
* `layer`:  the layer to write the field definition.
* `field`:  the geometry field definition to write to disk.
* `approx`: If `true`, the field may be created in a slightly different form
            depending on the limitations of the format driver.

### Remarks
You must use this to create new geometry fields on a real layer. Internally the
OGRFeatureDefn for the layer will be updated to reflect the new field.
Applications should never modify the OGRFeatureDefn used by a layer directly.

This function should not be called while there are feature objects in existence
that were obtained or created with the previous layer definition.

Not all drivers support this function. You can query a layer to check if it
supports it with the `GDAL.OLCCreateGeomField` capability. Some drivers may only
support this method while there are still no features in the layer. When it is
supported, the existing features of the backing file/database should be updated
accordingly.

Drivers may or may not support not-null constraints. If they support creating
fields with not-null constraints, this is generally before creating any feature
to the layer.
"""
function addgeomdefn!(
    layer::L,
    field::AbstractGeomFieldDefn,
    approx::Bool = false,
)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_creategeomfield(layer.ptr, field.ptr, approx)
    # OGRERR_NONE on success.
    @ogrerr result "Failed to create new geometry field"
    return layer
end

# """
# Delete the field at index `i` on a layer.

# You must use this to delete existing fields on a real layer. Internally the
# OGRFeatureDefn for the layer will be updated to reflect the deleted field.
# Applications should never modify the OGRFeatureDefn used by a layer directly.

# This function should not be called while there are feature objects in existence
# that were obtained or created with the previous layer definition.

# Not all drivers support this function. You can query a layer to check if it
# supports it with the OLCDeleteField capability. Some drivers may only support
# this method while there are still no features in the layer. When it is
# supported, the existing features of the backing file/database should be updated
# accordingly.

# ### Parameters
# * `layer`: handle to the layer.
# * `i`: index of the field to delete.
# """
# function deletefield!(layer::AbstractFeatureLayer, i::Integer)
#     result = GDAL.ogr_l_deletefield(layer.ptr, i)
#     @ogrerr result "Failed to delete field $i"
#     layer
# end

# # function deletegeomdefn!(layer::AbstractFeatureLayer, i::Integer)
# #     error("OGR_L_DeleteGeomField() is not implemented in GDAL yet.")
# # end

# """
# Reorder all the fields of a layer.

# You must use this to reorder existing fields on a real layer. Internally the
# OGRFeatureDefn for the layer will be updated to reflect the reordering of the
# fields. Applications should never modify the OGRFeatureDefn used by a layer
# directly.

# This method should not be called while there are feature objects in existence
# that were obtained or created with the previous layer definition.

# panMap is such that,for each field definition at position i after reordering,
# its position before reordering was panMap[i].

# For example, let suppose the fields were "0","1","2","3","4" initially.
# ReorderFields([0,2,3,1,4]) will reorder them as "0","2","3","1","4".

# Not all drivers support this method. You can query a layer to check if it
# supports it with the OLCReorderFields capability. Some drivers may only support
# this method while there are still no features in the layer. When it is
# supported, the existing features of the backing file/database should be updated
# accordingly.

# ### Parameters
# * `layer`: handle to the layer.
# * `indices`: an array of GetLayerDefn()->OGRFeatureDefn::GetFieldCount()
#             elements which is a permutation of
#                 `[0, GetLayerDefn()->OGRFeatureDefn::GetFieldCount()-1]`.
# """
# function reorderfields!(layer::AbstractFeatureLayer, indices::Vector{Cint})
#     result = GDAL.ogr_l_reorderfields(layer.ptr, indices)
#     @ogrerr result "Failed to reorder the fields of layer according to $indices"
#     layer
# end

# """
# Reorder an existing field on a layer.

# This method is a convenience wrapper of ReorderFields() dedicated to move a
# single field. It is a non-virtual method, so drivers should implement
# ReorderFields() instead.

# You must use this to reorder existing fields on a real layer. Internally the
# OGRFeatureDefn for the layer will be updated to reflect the reordering of the
# fields. Applications should never modify the OGRFeatureDefn used by a layer
# directly.

# This method should not be called while there are feature objects in existence
# that were obtained or created with the previous layer definition.

# The field definition that was at initial position iOldFieldPos will be moved at
# position iNewFieldPos, and elements between will be shuffled accordingly.

# For example, let suppose the fields were "0","1","2","3","4" initially.
# ReorderField(1, 3) will reorder them as "0","2","3","1","4".

# Not all drivers support this method. You can query a layer to check if it
# supports it with the OLCReorderFields capability. Some drivers may only support
# this method while there are still no features in the layer. When it is
# supported, the existing features of the backing file/database should be updated
# accordingly.

# ### Parameters
# * `layer`: handle to the layer.
# * `oldpos`: previous position of the field to move. Must be in the range
#             [0,GetFieldCount()-1].
# * `newpos`: new position of the field to move. Must be in the range
#             [0,GetFieldCount()-1].
# """
# function reorderfield!(
#         layer::AbstractFeatureLayer,
#         oldpos::Integer,
#         newpos::Integer
#     )
#     result = GDAL.ogr_l_reorderfield(layer.ptr, oldpos, newpos)
#     @ogrerr result "Failed to reorder field from $oldpos to $newpos."
#     layer
# end

# """
# Alter the definition of an existing field on a layer.

# You must use this to alter the definition of an existing field of a real layer.
# Internally the OGRFeatureDefn for the layer will be updated to reflect the
# altered field. Applications should never modify the OGRFeatureDefn used by a
# layer directly.

# This method should not be called while there are feature objects in existence
# that were obtained or created with the previous layer definition.

# Not all drivers support this method. You can query a layer to check if it
# supports it with the OLCAlterFieldDefn capability. Some drivers may only
# support this method while there are still no features in the layer. When it is
# supported, the existing features of the backing file/database should be updated
# accordingly. Some drivers might also not support all update flags.

# ### Parameters
# * `layer`:        handle to the layer.
# * `i`:            index of the field whose definition must be altered.
# * `newfielddefn`: new field definition
# * `flags`: combination of ALTER_NAME_FLAG, ALTER_TYPE_FLAG,
#             ALTER_WIDTH_PRECISION_FLAG, ALTER_NULLABLE_FLAG and
#             ALTER_DEFAULT_FLAG to indicate which of the name and/or type and/or
#             width and precision fields and/or nullability from the new field
#             definition must be taken into account.
# """
# function alterfielddefn!(
#         layer::AbstractFeatureLayer,
#         i::Integer,
#         newfielddefn::FieldDefn,
#         flags::UInt8
#     )
#     result = OGR.alterfielddefn(layer.ptr, i, newfielddefn.ptr, flags)
#     @ogrerr result "Failed to alter fielddefn of field $i."
#     layer
# end

# "For datasources which support transactions, creates a transaction."
# function starttransaction(layer::AbstractFeatureLayer)
#     result = GDAL.ogr_l_starttransaction(layer.ptr)
#     @ogrerr result "Failed to start transaction."
#     layer
# end

# "For datasources which support transactions, commits a transaction."
# function committransaction(layer::AbstractFeatureLayer)
#     result = GDAL.ogr_l_committransaction(layer.ptr)
#     @ogrerr result "Failed to commit transaction."
#     layer
# end

# """
# For datasources which support transactions, RollbackTransaction will roll back
# a datasource to its state before the start of the current transaction.
# """
# function rollbacktransaction(layer::AbstractFeatureLayer)
#     result = GDAL.ogr_l_rollbacktransaction(layer.ptr)
#     @ogrerr result "Failed to rollback transaction."
#     layer
# end

"""
    reference(layer::AbstractFeatureLayer)

Increment layer reference count.

### Returns
The reference count after incrementing.
"""
reference(layer::AbstractFeatureLayer)::Integer =
    GDAL.ogr_l_reference(layer.ptr)

"""
    dereference(layer::AbstractFeatureLayer)

Decrement layer reference count.

### Returns
The reference count after decrementing.
"""
dereference(layer::AbstractFeatureLayer)::Integer =
    GDAL.ogr_l_dereference(layer.ptr)

"""
    nreference(layer::AbstractFeatureLayer)

The current reference count for the layer object itself.
"""
nreference(layer::AbstractFeatureLayer)::Integer =
    GDAL.ogr_l_getrefcount(layer.ptr)

# """
# Flush pending changes to disk.

# This call is intended to force the layer to flush any pending writes to disk,
# and leave the disk file in a consistent state. It would not normally have any
# effect on read-only datasources.

# Some layers do not implement this method, and will still return OGRERR_NONE.
# The default implementation just returns OGRERR_NONE. An error is only returned
# if an error occurs while attempting to flush to disk.

# In any event, you should always close any opened datasource with
# DestroyDataSource() that will ensure all data is correctly flushed.

# ### Returns
# OGRERR_NONE if no error occurs (even if nothing is done) or an error code.
# """
# function synctodisk!(layer::AbstractFeatureLayer)
#     result = GDAL.ogr_l_synctodisk(layer.ptr)
#     @ogrerr result "Failed to flush pending changes to disk"
#     layer.ptr = C_NULL
# end

# """
# Return the total number of features read.

# Warning: not all drivers seem to update this count properly.
# """
# getfeaturesread(layer::AbstractFeatureLayer) =
#     GDAL.ogr_l_getfeaturesread(layer.ptr)

"""
    fidcolumnname(layer::AbstractFeatureLayer)

The name of the FID column in the database, or \"\" if not supported.
"""
fidcolumnname(layer::AbstractFeatureLayer)::String =
    GDAL.ogr_l_getfidcolumn(layer.ptr)

"""
    geomcolumnname(layer::AbstractFeatureLayer)

The name of the geometry column in the database, or \"\" if not supported.
"""
geomcolumnname(layer::AbstractFeatureLayer)::String =
    GDAL.ogr_l_getgeometrycolumn(layer.ptr)

"""
    setignoredfields!(layer::AbstractFeatureLayer, fieldnames)

Set which fields can be omitted when retrieving features from the layer.

### Parameters
* `fieldnames`: an array of field names terminated by NULL item. If NULL is
passed, the ignored list is cleared.

### Remarks
If the driver supports this functionality (testable using `OLCIgnoreFields`
capability), it will not fetch the specified fields in subsequent calls to
`GetFeature()`/`nextfeature()` and thus save some processing time and/or
bandwidth.

Besides field names of the layers, the following special fields can be passed:
`"OGR_GEOMETRY"` to ignore geometry and `"OGR_STYLE"` to ignore layer style.

By default, no fields are ignored.
"""
function setignoredfields!(
    layer::L,
    fieldnames,
)::L where {L<:AbstractFeatureLayer}
    result = GDAL.ogr_l_setignoredfields(layer.ptr, fieldnames)
    # OGRERR_NONE if all field names have been resolved (even if the driver
    # does not support this method)
    @ogrerr result "Failed to set ignored fields $fieldnames."
    return layer
end

# """
# Intersection of two layers.

# The result layer contains features whose geometries represent areas that are
# common between features in the input layer and in the method layer. The features
# in the result layer have attributes from both input and method layers. The
# schema of the result layer can be set by the user or, if it is empty, is
# initialized to contain all fields in the input and method layers.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `output`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.

# ### Keyword Arguments
# * `options`: NULL terminated list of options (may be NULL).
# * `progressfunc`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `progressdata`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is:

# * `SKIP_FAILURES=YES/NO`. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO`. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# * `USE_PREPARED_GEOMETRIES=YES/NO`. Set to NO to not use prepared geometries to
#     pretest intersection of features of method layer with features of this layer
# * `PRETEST_CONTAINMENT=YES/NO`. Set to YES to pretest the containment of
#     features of method layer within the features of this layer. This will speed
#     up the method significantly in some cases. Requires that the prepared
#     geometries are in effect.
# * `KEEP_LOWER_DIMENSION_GEOMETRIES=YES/NO`. Set to NO to skip result features
#     with lower dimension geometry that would otherwise be added to the result
#     layer. The default is to add but only if the result layer has an unknown
#     geometry type.
# """
# function intersection(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options                 = StringList(C_NULL),
#         progressfunc::Function  = GDAL.gdaldummyprogress,
#         progressdata            = C_NULL
#     )
#     result = GDAL.ogr_l_intersection(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to compute the intersection of the two layers"
#     output
# end

# """
# Union of two layers.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `result`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.
# * `papszOptions`: NULL terminated list of options (may be NULL).
# * `pfnProgress`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `pProgressArg`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is:

# * `SKIP_FAILURES=YES/NO. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# * `USE_PREPARED_GEOMETRIES=YES/NO. Set to NO to not use prepared geometries to
#     pretest intersection of features of method layer with features of this layer
# * `KEEP_LOWER_DIMENSION_GEOMETRIES=YES/NO. Set to NO to skip result features
#     with lower dimension geometry that would otherwise be added to the result
#     layer. The default is to add but only if the result layer has an unknown
#     geometry type.
# """
# function union(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options                 = StringList(C_NULL),
#         progressfunc::Function  = GDAL.gdaldummyprogress,
#         progressdata            = C_NULL
#     )
#     result = GDAL.ogr_l_union(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to compute the union of the two layers"
#     output
# end

# """
# Symmetrical difference of two layers.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `result`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.
# * `papszOptions`: NULL terminated list of options (may be NULL).
# * `pfnProgress`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `pProgressArg`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is :

# * `SKIP_FAILURES=YES/NO`. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO`. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# """
# function symdifference(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options                 = StringList(C_NULL),
#         progressfunc::Function  = GDAL.gdaldummyprogress,
#         progressdata            = C_NULL
#     )
#     result = GDAL.ogr_l_symdifference(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to compute the sym difference of the two layers"
#     output
# end

# """
# Identify the features of this layer with the ones from the identity layer.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `result`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.
# * `papszOptions`: NULL terminated list of options (may be NULL).
# * `pfnProgress`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `pProgressArg`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is :

# * `SKIP_FAILURES=YES/NO`. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO`. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# * `USE_PREPARED_GEOMETRIES=YES/NO`. Set to NO to not use prepared geometries to
#     pretest intersection of features of method layer with features of this layer
# * `KEEP_LOWER_DIMENSION_GEOMETRIES=YES/NO`. Set to NO to skip result features
#     with lower dimension geometry that would otherwise be added to the result
#     layer. The default is to add but only if the result layer has an unknown
#     geometry type.
# """
# function identity(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options                 = StringList(C_NULL),
#         progressfunc::Function  = GDAL.gdaldummyprogress,
#         progressdata            = C_NULL
#     )
#     result = GDAL.ogr_l_identity(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to compute the identity of the two layers"
#     output
# end

# """
# Update this layer with features from the update layer.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `result`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.
# * `papszOptions`: NULL terminated list of options (may be NULL).
# * `pfnProgress`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `pProgressArg`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is :

# * `SKIP_FAILURES=YES/NO`. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO`. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# """
# function update(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options                 = StringList(C_NULL),
#         progressfunc::Function  = GDAL.gdaldummyprogress,
#         progressdata            = C_NULL
#     )
#     result = GDAL.ogr_l_update(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to update the layer"
#     output
# end

# """
# Clip off areas that are not covered by the method layer.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `result`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.
# * `papszOptions`: NULL terminated list of options (may be NULL).
# * `pfnProgress`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `pProgressArg`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is :

# * `SKIP_FAILURES=YES/NO`. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO`. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# """
# function clip(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options                 = StringList(C_NULL),
#         progressfunc::Function  = GDAL.gdaldummyprogress,
#         progressdata            = C_NULL
#     )
#     result = GDAL.ogr_l_clip(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to clip the input layer"
#     output
# end

# """
# Remove areas that are covered by the method layer.

# ### Parameters
# * `input`: the input layer. Should not be NULL.
# * `method`: the method layer. Should not be NULL.
# * `result`: the layer where the features resulting from the operation
#     are inserted. Should not be NULL. See the note about the schema.
# * `papszOptions`: NULL terminated list of options (may be NULL).
# * `pfnProgress`: a GDALProgressFunc() compatible callback function for
#     reporting progress or NULL.
# * `pProgressArg`: argument to be passed to pfnProgress. May be NULL.

# ### Returns
# an error code if there was an error or the execution was interrupted,
# OGRERR_NONE otherwise.

# ### Additional Remarks
# If the schema of the result is set by user and contains fields that have the
# same name as a field in input and in method layer, then the attribute in the
# result feature will get the value from the feature of the method layer.

# For best performance use the minimum amount of features in the method layer and
# copy it into a memory layer. This method relies on GEOS support. Do not use
# unless the GEOS support is compiled in. The recognized list of options is :

# * `SKIP_FAILURES=YES/NO`. Set it to YES to go on, even when a feature could not
#     be inserted or a GEOS call failed.
# * `PROMOTE_TO_MULTI=YES/NO`. Set it to YES to convert Polygons into
#     MultiPolygons, or LineStrings to MultiLineStrings.
# * `INPUT_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the input layer.
# * `METHOD_PREFIX=string`. Set a prefix for the field names that will be created
#     from the fields of the method layer.
# """
# function erase(
#         input::AbstractFeatureLayer,
#         method::AbstractFeatureLayer,
#         output::AbstractFeatureLayer;
#         options = StringList(C_NULL),
#         progressfunc::Function = GDAL.gdaldummyprogress,
#         progressdata = C_NULL
#     )
#     result = GDAL.ogr_l_erase(
#         input.ptr,
#         method.ptr,
#         output.ptr,
#         options,
#         @cplprogress(progressfunc),
#         progressdata
#     )
#     @ogrerr result "Failed to remove areas covered by the method layer."
#     output
# end


text:
module TensorLab
	import MATLAB: mxarray
	using MATLAB
	include("matlabTypes.jl")
	export Options, MatFcn

	# Create a table of functions and the number of outputs they return.
	# TODO: These functions have a variable number of outputs
	# - cpd (probably variants as well) (1 or 2)
	# - similar mechanism for lmlra
	# - mlsvd and variants (2 or 3)
	functions = [
		(:cpd, 1),
		(:cpd_als, 1),
		(:cpd_els, 1),
		(:cpd_gevd, 1),
		(:cpd_minf, 1),
		(:cpd_rbs, 1),
		(:cpd_rnd, 1),
		(:cpd3_sd, 1),
		(:cpd3_sgsd, 1),
		(:cpderr, 1),
		(:cpdgen, 1),
		(:cpdres, 1),
		(:fmt, 1),
		(:frob, 1),
		(:frobcpdres, 1),
		(:frobmlrares, 1),
		(:fromll1res, 1),
		(:ful, 1),
		(:hankelize, 1),
		(:isvalidtensor, 1),
		(:kr, 1),
		(:kron, 1),
		(:ll1, 1),
		(:ll1_gevd, 1),
		(:ll1_minf, 1),
		(:ll1_rnd, 1),
		(:ll1convert, 1),
		(:ll1gen, 1),
		(:ll1res, 1),
		(:lmlra, 2),
		(:lmlra3_dgn, 2),
		(:lmlra3_rtr, 2),
		(:lmlra_aca, 2),
		(:lmlra_hooi, 2),
		(:lmlra_minf, 2),
		(:lmlra_nls, 2),
		(:lmlra_rnd, 2),
		(:lmlraerr, 1),
		(:lmlragen, 1),
		(:lmlrares, 1),
		(:mat2tens, 1),
		(:mlrank, 1),
		(:mlrankest, 0),
		(:mlsvd, 3),
		(:mlsvd_rsi, 3),
		(:mlsvds, 3),
		(:rankest, 0),
		(:surf3,0),
		(:slice3,0),
		(:tens2mat, 1),
		(:tmprod, 1),
		(:visualize, 0),
		(:voxel3,0),
	]

	for (fcn, numout) in functions
		@eval $fcn(args...) = mxcall(Symbol($fcn), $numout, args...)
		@eval export $fcn
	end
end

ext: jl


text:
"""
    convert_to_option(::Type{OptionType}, ::Type{ValueType}, x) where {OptionType, ValueType}

Convert `x` to type `ValueType` for option type `OptionType`. This is similar to `Base.convert`,
when creating an instance of the option type, but one can use this to avoid type piracy.
"""
convert_to_option(::Type, ::Type{A}, x) where {A} = nothing

convert_to_option(::Type, ::Type{VersionNumber}, x::String) = VersionNumber(x)

function convert_union_to_option(::Type{T}, ::Type{A}, x) where {T, A}
    if !(A isa Union)
        v = convert_to_option(T, A, x)
        v === nothing || return v
        return
    end

    v = convert_union_to_option(T, A.a, x)
    v === nothing || return v
    return convert_union_to_option(T, A.b, x)
end


text:
"""
    homodyne(mag, phase)

    homodyne(mag, phase; dims, σ)

    homodyne(I; kwargs...)

Performs homodyne filtering via division of complex complex smoothing.
"""
homodyne, homodyne!

homodyne(mag, phase; kwargs...) = homodyne(mag .* exp.(1im .* phase); kwargs...)
homodyne(I; kwargs...) = homodyne!(copy(I); kwargs...)
homodyne!(I; dims=1:2, σ=8 .* ones(length(dims))) = I ./= gaussiansmooth3d(I, σ; dims)

"""
    NumART2star(image::AbstractArray{T,4}, TEs) where T

Performs T2* calculation on 4D-multi-echo magnitude data.
https://doi.org/10.1002/mrm.10283
"""
function NumART2star(image::AbstractArray{T,4}, TEs) where T
    neco = length(TEs)
    t2star(m) = (TEs[end]-TEs[1]) / 2(neco-1) * (m[1]+m[end]+2sum(m[2:end-1])) / (m[1]-m[end])
    return [t2star(image[I,:]) for I in CartesianIndices(size(image)[1:3])]
end

"""
    r2s_from_t2s(t2s) = 1000 ./ t2s

Converts from T2* [ms] to R2* [1/s] values.
"""
r2s_from_t2s(t2s) = 1000 ./ t2s


text:
using ArgParse
using BenchmarkTools
using Cambrian
using CartesianGeneticProgramming
using IICGP
using LinearAlgebra
import Cambrian.mutate  # mutate function scope definition

# rom_names = setdiff(getROMList(), ["pacman", "surround"])
rom_names = [
    "boxing",
    "centipede",
    "demon_attack",
    "enduro",
    "freeway",
    "kung_fu_master",
    "space_invaders",
    "riverraid",
    "pong"
]

function generate_io(rom_name::String="freeway", frame_number::Int64=30)
    inps = Array{Array{Array{UInt8,2},1},1}()
    outs = Array{Array{Float64,1},1}()
    for rom in rom_names
        r, g, b = IICGP.load_rgb(rom, 30)
        i = IICGP.CGPFunctions.f_binary(r, g, [0.5])
        j = IICGP.CGPFunctions.f_erode(i, i, [0.5])
        k = IICGP.CGPFunctions.f_subtract(j, g, [0.5])
        feature = IICGP.ReducingFunctions.max_pool_reduction(k, 5)
        t1 = 0.5 * (feature[1, 1] + feature[2, 2]) * feature[1, 1]
        t2 = feature[3, 3]
        push!(inps, [r, g, b])
        push!(outs, [t1, t2])
    end
    return inps, outs
end

function fitness(encoder::CGPInd, controller::CGPInd,
                 inps::Array{Array{Array{UInt8,2},1},1},
                 outs::Array{Array{Float64,1},1})
    score = 0.0
    for i in eachindex(inps)
        pred = IICGP.process(encoder, controller, inps[i], encoder_cfg.features_size)
        score -= LinearAlgebra.norm(pred - outs[i])
    end
    [score]
end

# Read configuration file
s = ArgParseSettings()
@add_arg_table! s begin
    "--encoder_cfg"
    help = "configuration script"
    default = "cfg/encoder.yaml"
    "--controller_cfg"
    help = "configuration script"
    default = "cfg/controller.yaml"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end
args = parse_args(ARGS, s)
n_in = 3  # RGB images
n_out = 2  # Two scalar values
inps, outs = generate_io()
img_size = size(inps[1][1])

encoder_cfg = get_config(
    args["encoder_cfg"];
    function_module=IICGP.CGPFunctions,
    n_in=n_in,
    img_size=img_size
)

n_in_controller = encoder_cfg.n_out * encoder_cfg.features_size^2
controller_cfg = get_config(
    args["controller_cfg"];
    n_in=n_in_controller,
    n_out=n_out
)

#=
# Test I/O without evolution
foo = IPCGPInd(encoder_cfg)
bar = CGPInd(controller_cfg)
out = IICGP.process(foo, bar, inps[1], encoder_cfg.features_size)
=#

if length(args["ind"]) > 0
    ind = CGPInd(cfg, read(args["ind"], String))
    ftn = fitness(ind, inps, outs)
    println("Fitness: ", ftn)
else
    # Define mutate and fit functions
    function mutate(ind::CGPInd, ind_type::String)
        if ind_type == "encoder"
            return goldman_mutate(encoder_cfg, ind, init_function=IPCGPInd)
        elseif ind_type == "controller"
            return goldman_mutate(controller_cfg, ind)
        end
    end
    fit(encoder::CGPInd, controller::CGPInd) = fitness(encoder, controller, inps, outs)
    # Create an evolution framework
    e = IICGP.DualCGPEvolution(encoder_cfg, controller_cfg, fit,
                               encoder_init_function=IPCGPInd)
    # Run evolution
    run!(e)
end


text:
# Tests to simulate traits

## Get a network and make it ultrametric
net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");
#plot(net, useEdgeLength = true,  showEdgeNumber=true)

@testset "Simulate function against fixed values" begin

## Simulate a BM
Random.seed!(17920921); # fix the seed
pars = ParamsBM(1, 0.1); # params of a BM
@test_logs show(devnull, pars)

sim = simulate(net, pars); # simulate according to a BM
@test_logs show(devnull, sim)
@test_throws ErrorException sim[:Tips, :Broken]

# Extract simulated values
traitsTips = sim[:Tips];
traitsNodes = sim[:InternalNodes];

# Expected values
traitsTipsExp = [0.6455995230091043,-0.22588106270381064,0.05703904710270408,-0.692650796714688,1.578622599565194,1.4106438068675058,1.9166557600811194,1.0579005662214953,1.2340762902144904,1.4130757789427886,0.7115737497673081,2.201943319276716];
traitsNodesExp = [-0.3481603206484607,-0.6698437934551933,-0.018135478212541654,-0.33844527112230455,-0.0717742134084467,0.19417331380691694,1.3919535151447147,1.5106942025265466,1.2526948727806593,1.1552248152172964,1.224823113083187,1.0617270280846993,1.0436547766241817,1.0];

@test traitsTips ≈ traitsTipsExp
@test traitsNodes ≈ traitsNodesExp

end

###############################################################################
## Test of distibution
###############################################################################
@testset "Simulate test distribution" begin

## Generate some values
Random.seed!(18480224); # fix the seed
pars = ParamsBM(1, 0.1); # params of a BM
N = 50000
S = length(tipLabels(net));
values = zeros(Float64, (S, N));
for i = 1:N
    values[:,i] = simulate(net, pars)[:Tips]
end

## Check that each tip has same mean (1)
for s in 1:S
    @test mean(values[s, :]) ≈ pars.mu atol=0.01
end

## Check for variances
V = sharedPathMatrix(net);
Sig = V[:Tips] * pars.sigma2;
for s in 1:S
    for t in s:S
        @test cov(values[s, :], values[t, :]) ≈ Sig[s, t] atol=0.01
    end
end

end
###############################################################################
## With Shifts
###############################################################################
@testset "Simulate with Shifts" begin
global net
net = readTopology("(A:2.5,((B:1,#H1:0.5::0.4):1,(C:1,(D:0.5)#H1:0.5::0.6):1):0.5);")

## Test construction function
@test_throws ErrorException ShiftNet(net.edge[7], 3.0,  net) # can't put a shift on hybrid branch
@test_throws ErrorException ShiftNet(net.node[6], 3.0,  net) # can't put a shift on hybrid branch
@test ShiftNet(net).shift ≈ zeros(length(net.node))
@test ShiftNet(net.edge[8], 3.0,  net).shift ≈ ShiftNet([net.edge[8]], [3.0],  net).shift
@test ShiftNet(net.edge[8], 3.0,  net).shift ≈ ShiftNet(net.node[7], 3.0,  net).shift
@test ShiftNet(net.node[7], 3.0,  net).shift ≈ ShiftNet([net.node[7]], [3.0],  net).shift

## Concatenate function
sh1 = ShiftNet(net.node[7], 3.0,  net)*ShiftNet(net.node[9], -2.1,  net)
@test sh1.shift ≈ ShiftNet([net.node[7], net.node[9]], [3.0, -2.1],  net).shift
@test_throws ErrorException sh1*ShiftNet(net.edge[7], 2.0,  net) # can't concatenate if the two affect the same edges
@test sh1.shift ≈ (sh1*ShiftNet([net.node[7]], [3.0],  net)).shift

## Values and edge numbers functions
sh = ShiftNet(net.node[7], 3.0,  net)
@test getShiftEdgeNumber(sh) == [8]
@test all(getShiftValue(sh) .== [3.0])

## Hybrid shifts
@test shiftHybrid([2.0], net).shift ≈ ShiftNet(net.edge[6], 2.0, net).shift
@test shiftHybrid(2.0, net).shift ≈ shiftHybrid([2.0], net).shift

## Test simulate

# No shift on root
@test_throws ErrorException simulate(net, ParamsBM(1.0, 0.1, ShiftNet(net.node[9], 3.0,  net)))

@test ParamsBM(1.0, 1.0, net).shift.shift ≈ ParamsBM(1.0, 1.0, ShiftNet(net)).shift.shift

pars = ParamsBM(1, 0.1, ShiftNet(net.edge[8], 3.0,  net)); # params of a BM
@test_logs show(devnull, pars)
@test_logs show(devnull, pars.shift)

Random.seed!(17920921); # fix the seed
sim = simulate(net, pars); # simulate according to a BM
@test_logs show(devnull, sim)

traitsTips = sim[:Tips];
traitsNodes = sim[:InternalNodes];

meansTips = sim[:Tips, :Exp];
meansNodes = sim[:InternalNodes, :Exp];

# Expected values
meansTipsExp = [1.0 1.0 1.0+3.0 1.0+3.0*0.6]';
traitsTipsExp = [0.9230584254019785 1.4363450675116494 4.180396447825185 3.299820483104897]';

meansNodesExp = [1.0 1.0+3.0*0.6 1.0+3.0 1.0 1.0]';
traitsNodesExp = [1.50594336537754 3.296894371572107 4.346436961253621 1.3212328642182367 1.0]';

@test traitsTips ≈ traitsTipsExp
@test traitsNodes ≈ traitsNodesExp
@test meansTips ≈ meansTipsExp
@test meansNodes ≈ meansNodesExp

# Test same as MultiBM
pars = ParamsMultiBM([1.0], 0.1*ones(1,1), ShiftNet(net.edge[8], 3.0,  net));
simMulti = simulate(net, pars);
@test simMulti[:Tips, :Exp] ≈ sim[:Tips, :Exp]'
@test simMulti[:InternalNodes, :Exp] ≈ sim[:InternalNodes, :Exp]'

###############################################################################
## Test of distibution - with shifts
###############################################################################

## Generate some values
Random.seed!(18480224); # fix the seed
@test_throws ErrorException ParamsBM(1, 0.1, ShiftNet(net.edge[8], [3.0, 1.0],  net))
pars = ParamsBM(1, 0.1, ShiftNet(net.edge[8], 3.0,  net)); # params of a BM
N = 50000
S = length(tipLabels(net));
values = zeros(Float64, (S, N));
for i = 1:N
    values[:,i] = simulate(net, pars)[:Tips]
end

## Check that each tip has same mean (1)
expectations = simulate(net, pars)[:Tips,:Exp]
for s in 1:S
    @test mean(values[s, :]) ≈ expectations[s] atol=1e-2
end

## Check for variances
V = sharedPathMatrix(net);
Sig = V[:Tips] * pars.sigma2;
for s in 1:S
    for t in s:S
        @test cov(values[s, :], values[t,:]) ≈ Sig[s, t] atol=1e-2
    end
end

end


text:
############################################################################
#
#  Interpolation.jl
#
#  Purpose:    This file contains functions to perform linear interpolation
#              of multi-dimensional tables
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################

function interp1( x, X, Y )
    ## Initialize Indicies
    ii = 1;
    VAL = 0.0;
    x1 = 0.0;
    x2 = 0.0;

    n1 = length(X);
    ifound = 0;
    xfac = 0;

    ## Search for X
    while (ifound == 0)
        tmp1 = x - X[ii];
        tmp2 = x - X[ii + 1];
        tmp3 = tmp1 * tmp2;

        if (tmp3 <= 0)
            x1 = X[ii];
            x2 = X[ii + 1];
            ifound = 1;
        else
            ii = ii + 1;
        end
    end

    f1 = Y[ii];
    f2 = Y[ii + 1];

    VAL = (f1 * (x2 - x) + f2 * (x - x1)) / (x2 - x1);
    return VAL;
end

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

function interp2(x, y, X, Y, Z)

    n1 = length(X);
    n2 = length(Y);
    ifound = 1;
    jfound = 1;
    xfac = 0;
    yfac = 0;

    # Search for X
    for i = 2:n1
        if ( (x >= X[i-1]) && ( x <= X[i] ) )
            ifound = i-1;
            xfac = (x - X[i-1])/( X[i] -  X[i-1] );
            break
        end

        if (i == n1) && (ifound == 0)
            xfac = 1.0;
            ifound = i-1;
        end

    end

    # Search for Y
    for j = 2:n2
        if ( (y >= Y[j-1]) && ( y <= Y[j] ) )
            jfound = j-1;
            yfac = (y - Y[j-1])/( Y[j] -  Y[j-1] );
            break;
        end

        if (j > n2) && (jfound == 0)
            yfac = 1.0;
            jfound = j-1;
        end

    end

    ZZ = (Z[ifound+1,:] - Z[ifound,:] )*xfac + Z[ifound,:];
    zout = (ZZ[jfound+1] - ZZ[jfound] )*yfac + ZZ[jfound];

    return zout
end

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

function interp3( x, y, z, X, Y, Z, VV )

    ##Initialize Indicies
    ii = 1;
    jj = 1;
    kk = 1;
    VAL = 0.0;
    x1 = 0.0;
    x2 = 0.0;
    y1 = 0.0;
    y2 = 0.0;
    z1 = 0.0;
    z2 = 0.0;

    ##Search for X
    ifound = 0;
    while (ifound == 0)
        tmp1 = x - X[ii];
        tmp2 = x - X[ii + 1];
        tmp3 = tmp1 * tmp2;

        if (tmp3 <= 0)
            x1 = X[ii];
            x2 = X[ii + 1];
            ifound = 1;

        else
            ii = ii + 1;
        end
    end

    ##Search for Y
    jfound = 0;
    while (jfound == 0)
        tmp1 = y - Y[jj];
        tmp2 = y - Y[jj + 1];
        tmp3 = tmp1 * tmp2;

        if (tmp3 <= 0)
            y1 = Y[jj];
            y2 = Y[jj + 1];
            ifound = 1;
        else
            jj = jj + 1;
        end
    end

    ##Search for Z
    kfound = 0;
    while (ifound == 0)
        tmp1 = z - Z[kk];
        tmp2 = z - Z[kk + 1];
        tmp3 = tmp1 * tmp2;

        if (tmp3 <= 0)
            z1 = Z[kk];
            z2 = Z[kk + 1];
            ifound = 1;
        else
            kk = kk + 1;
        end
    end

    f111 = VV[ii][jj][kk];
    f211 = VV[ii + 1][jj][kk];
    f121 = VV[ii][jj + 1][kk];
    f112 = VV[ii][jj][kk + 1];
    f221 = VV[ii + 1][jj + 1][kk];
    f212 = VV[ii + 1][jj][kk + 1];
    f122 = VV[ii][jj + 1][kk + 1];
    f222 = VV[ii + 1][jj + 1][kk + 1];

    VAL = (-f112 * y * z * x2 - f112 * x * z * y2 - f112 * x * y * z1 + f112 * z * x2 * y2 + f112 * y * x2 * z1 + f112 * x * y2 * z1 - f221 * y * x1 * z2 - f221 * x * y1 * z2 + f221 * x1 * y1 * z2 + f112 * x * y * z - f221 * z * x1 * y1 - f112 * x2 * y2 * z1 + f121 * x * y * z - f121 * y * z * x2 - f121 * x * z * y1 - f121 * x * y * z2 + f121 * z * x2 * y1 + f121 * y * x2 * z2 + f121 * x * y1 * z2 - f121 * x2 * y1 * z2 + f211 * x * y * z - f211 * y * z * x1 - f212 * z * x1 * y2 - f212 * y * x1 * z1 + f212 * x * z * y2 + f212 * x * y * z1 - f111 * x * y * z + f111 * y * z * x2 + f111 * x * z * y2 + f111 * x * y * z2 - f111 * z * x2 * y2 - f111 * y * x2 * z2 - f111 * x * y2 * z2 + f111 * x2 * y2 * z2 + f211 * x * y2 * z2 - f211 * x1 * y2 * z2 - f211 * x * z * y2 - f211 * x * y * z2 + f211 * z * x1 * y2 + f211 * y * x1 * z2 + f212 * y * z * x1 + f221 * x * y * z2 + f221 * y * z * x1 + f221 * x * z * y1 - f221 * x * y * z - f212 * x * y2 * z1 + f212 * x1 * y2 * z1 - f122 * x * y1 * z1 + f122 * x2 * y1 * z1 - f212 * x * y * z - f122 * y * x2 * z1 + f122 * y * z * x2 + f122 * x * z * y1 + f122 * x * y * z1 - f122 * z * x2 * y1 - f222 * y * z * x1 - f222 * x * z * y1 - f222 * x * y * z1 + f222 * x * y * z - f122 * x * y * z - f222 * x1 * y1 * z1 + f222 * x * y1 * z1 + f222 * z * x1 * y1 + f222 * y * x1 * z1) / (-z1 + z2) / (-y2 + y1) / (-x2 + x1);

    return VAL;
end

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


text:
using Documenter, PairwiseListMatrices

makedocs(
    doctest = true,
    format = :html,
    sitename = "PairwiseListMatrices",
    modules = [PairwiseListMatrices],
    pages = [
        "index.md",
        "api.md"
    ]
)

deploydocs(
    repo   = "github.com/diegozea/PairwiseListMatrices.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
using Test, AutoARIMA, Statistics

@testset "Time Series Analysis Forecasting and Control" begin
    @testset "Table6.5" begin
        arma = fit(ARMAParams(true, 1,1), seriesA)
        @test arma.μ ≈ 2.45 rtol = 0.2
        @test arma.ϕ ≈ [0.87] rtol = 0.05
        @test arma.θ ≈ [0.48] rtol = 0.05
        @test arma.σ2 ≈ 0.098 rtol = 0.05

        θ, σ2 = innovations(diff(seriesA), 1, n=10)
        @test θ ≈ [0.53] rtol = 0.01
        @test σ2 ≈ 0.107 rtol = 0.01

        θ, σ2 = innovations(diff(seriesB), 1)
        @test θ ≈ [-0.09] rtol = 0.05
        @test σ2 ≈ 52.2 rtol = 0.05

        ar = fit(ARParams(false, 1), diff(seriesC))
        @test ar.μ ≈ 0.0 atol = 0.01
        @test ar.ϕ ≈ [0.81] rtol = 0.05
        @test ar.σ2 ≈ 0.019 rtol = 0.05

        θ, σ2 = innovations(diff(diff(seriesC)), 2, n=30)
        @test θ ≈ [0.09,0.07] rtol = 0.05
        @test σ2 ≈ 0.020 rtol = 0.05

        ar = fit(ARParams(1), seriesD)
        @test ar.μ ≈ 1.32 rtol = 0.05
        @test ar.ϕ ≈ [0.86] rtol = 0.05
        @test ar.σ2 ≈ 0.093 rtol = 0.05

        θ, σ2 = innovations(diff(seriesD), 1)
        @test θ ≈ [0.05] rtol = 0.05
        @test σ2 ≈ 0.096 rtol = 0.01

        ar = fit(ARParams(2), seriesE)
        @test ar.μ ≈ 14.9 rtol = 0.05
        @test ar.ϕ ≈ [1.32,-0.63] rtol = 0.05
        @test ar.σ2 ≈ 289.0 rtol = 0.05

        ar = fit(ARParams(3), seriesE)
        @test ar.μ ≈ 13.7 rtol = 0.05
        @test ar.ϕ ≈ [1.37,-0.74,0.08] rtol = 0.05
        @test ar.σ2 ≈ 287.0 rtol = 0.05

        ar = fit(ARParams(2), seriesF)
        @test ar.μ ≈ 58.3 rtol = 0.05
        @test ar.ϕ ≈ [-0.32,0.18] rtol = 0.05
        @test ar.σ2 ≈ 115.0 rtol = 0.05
    end

    @testset "Table3.1" begin
        ϕ = map(i -> partial_autocorrelation(seriesF, i), 1:15)
        @test ϕ ≈ [-0.39,0.18,0.00,-0.04,-0.07,-0.12,0.02,0.00,-0.06,0.00,0.14,-0.01,0.09,0.17,0.00] rtol = 0.05
    end

    @testset "Table2.1" begin
        ϕ = map(i -> autocorrelation(seriesF, i), 1:15)
        @test ϕ ≈ [-0.39,0.30,-0.17,0.07,-0.10,-0.05,0.04,-0.04,0.00,0.01,0.11,-0.07,0.15,0.04,0.01] rtol = 0.05
    end

    @testset "Section9.1" begin
        sarima = fit(MSARIMAParams([0,0], [1,1], [1,1], [1,12]), log.(seriesG))
        @test sarima.ϕ ≈ [1,0,0,0,0,0,0,0,0,0,0,1,-1] rtol = 0.05
        # @test sarima.θ[1] ≈ 0.4 rtol = 0.05
        # @test sarima.θ[12] ≈ 0.6 rtol = 0.05
        @test sarima.θ[13] ≈ -sarima.θ[1] * sarima.θ[12] rtol = 0.05
    end

    # @testset "Section10.1" begin
    #     arma = MA∞(fit(ARMAParams(1, 1), seriesA), 30)
    #     @test arma.θ ≈ [0.34,0.32,0.29,0.265,0.24,0.22,0.21,0.19,0.17,0.16,0.15,0.14,0.12,0.11,0.10,0.10,0.09,0.082,0.075,0.07,0.065,0.06,0.054,0.05,0.046,0.043,0.039,0.36,0.033,0.030] rtol = 0.05
    #     arima = MA∞(fit(ARIMAParams(0, 1, 1), seriesA), 30)
    #     @test arima.θ ≈ fill(0.29, 30)
    # end

    @testset "Table4.2" begin
        @test AR∞(toarma(ARIMAParams(1, 1, 1), ARMAModel{1,1}(1.0, [-0.3], [0.5], 1.0)), 7).ϕ ≈ [0.2,0.4,0.2,0.1,0.05,0.025,0.0125]
    end

end


text:
function collect{T}(enumerable::Enumerable{T})
    ret = Array(T,0)
    for i in enumerable
        push!(ret, i)
    end
    return ret
end

function collect{TS,Provider}(source::Queryable{TS,Provider})
    collect(Provider, source)
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ a96058a5-e563-4fa1-b014-feec6c4380ee
#hide
begin
#for sharing environments between notebooks
import Pkg
Pkg.activate(Base.current_project())
end

# ╔═╡ 68ed9c30-c6dc-11ec-191c-832624c25a17
#export
using TOML

# ╔═╡ 2c72c7dc-4bdb-4809-b33c-a34ae555120e
#hide
using PlutoUI

# ╔═╡ c02ea829-0a9b-4185-9784-daa42b290da2
#hide
using PlutoTest

# ╔═╡ b795e345-d5d4-4966-a77c-168e21156f46
#hide
include("../src/Export.jl")

# ╔═╡ cd7a8ece-670e-4db4-a689-bf1fe9a51c42
#hide
TableOfContents()

# ╔═╡ 7d2ff452-0cec-4307-950c-228adf3ad8d5
#hide
a = TOML.parsefile(joinpath("..", "Project.toml"))

# ╔═╡ 15629b65-7f23-47cd-967d-4b65a4cbc250
#hide
a["settings"]["github_url"]

# ╔═╡ 14591bed-affc-4fff-b37c-f8869028e803
md"## getsettings"

# ╔═╡ 1df8ad2c-612c-4c6e-a86f-ef03e9574c75
#export
begin
function getsetting(f:: AbstractString, k:: AbstractString)
	a = TOML.parsefile(f)
	return a[k]
end

function getsettings(f:: AbstractString, ks)
	res = TOML.parsefile(f)
	for k in ks
		#"$res"
		res = res[k]
	end
	return res
end
end

# ╔═╡ 871faecd-7481-4398-8d1c-2db69466eb6d
begin
f = joinpath("..", "Project.toml")
path = "https://github.com/sapal6/Nbdev.jl"
res = Dict("github_url" => "https://github.com/sapal6/Nbdev.jl")
end

# ╔═╡ 6e0d170e-126e-4999-b6ce-8f063871ffe8
#noop
@test getsetting(f, "settings") == res

# ╔═╡ c6ad861d-c0d6-4edf-88e2-267be1184bce
#noop
@test getsettings(f, ["settings", "github_url"]) == res["github_url"]

# ╔═╡ 7386505d-240c-46cf-9655-8ac30f41f148
#hide
Export.notebook2script(joinpath("..", "nbs"), joinpath("..", "src"))

# ╔═╡ Cell order:
# ╠═a96058a5-e563-4fa1-b014-feec6c4380ee
# ╠═68ed9c30-c6dc-11ec-191c-832624c25a17
# ╠═2c72c7dc-4bdb-4809-b33c-a34ae555120e
# ╠═c02ea829-0a9b-4185-9784-daa42b290da2
# ╠═cd7a8ece-670e-4db4-a689-bf1fe9a51c42
# ╠═b795e345-d5d4-4966-a77c-168e21156f46
# ╠═7d2ff452-0cec-4307-950c-228adf3ad8d5
# ╠═15629b65-7f23-47cd-967d-4b65a4cbc250
# ╠═14591bed-affc-4fff-b37c-f8869028e803
# ╠═1df8ad2c-612c-4c6e-a86f-ef03e9574c75
# ╠═871faecd-7481-4398-8d1c-2db69466eb6d
# ╠═6e0d170e-126e-4999-b6ce-8f063871ffe8
# ╠═c6ad861d-c0d6-4edf-88e2-267be1184bce
# ╠═7386505d-240c-46cf-9655-8ac30f41f148


text:
# ---
# title: 519. Random Flip Matrix
# id: problem519
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Random
# link: <https://leetcode.com/problems/random-flip-matrix/description/>
# hidden: true
# ---
#
# You are given the number of rows `n_rows` and number of columns `n_cols` of a
# 2D binary matrix where all values are initially 0. Write a function `flip`
# which chooses a 0 value [uniformly at
# random](https://en.wikipedia.org/wiki/Discrete_uniform_distribution), changes
# it to 1, and then returns the position `[row.id, col.id]` of that value. Also,
# write a function `reset` which sets all values back to 0.  **Try to minimize
# the number of calls to system 's Math.random()** and optimize the time and
# space complexity.
#
# Note:
#
#   1. `1 <= n_rows, n_cols <= 10000`
#   2. `0 <= row.id < n_rows` and `0 <= col.id < n_cols`
#   3. `flip` will not be called when the matrix has no 0 values left.
#   4. the total number of calls to `flip` and `reset` will not exceed 1000.
#
# **Example 1:**
#
#
#
#     Input: ["Solution","flip","flip","flip","flip"]
#     [[2,3],[],[],[],[]]
#     Output: [null,[0,1],[1,2],[1,0],[1,1]]
#
#
# **Example 2:**
#
#
#
#     Input: ["Solution","flip","flip","reset","flip"]
#     [[1,2],[],[],[],[]]
#     Output: [null,[0,0],[0,1],null,[0,0]]
#
# **Explanation of Input Syntax:**
#
# The input is two lists: the subroutines called and their arguments.
# `Solution`'s constructor has two arguments, `n_rows` and `n_cols`. `flip` and
# `reset` have no arguments. Arguments are always wrapped with a list, even if
# there aren't any.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using LinearAlgebra
using LMCLUS
using Statistics
using Distributions
using Test
import Random

@testset "MDL" begin

    @testset "Quantization" begin
        @test LMCLUS.MDL.univar([1.0]) ≈ [1 ./ 12.0]
        @test LMCLUS.MDL.optbins([1.], 1.) == [one(UInt64)]
        bins, ɛ, c, itr = LMCLUS.MDL.optquant([1.], 1e-2)
        @test bins[1] == 29
        @test ɛ < 1e-2
        bins, ɛ, c, itr = LMCLUS.MDL.optquant([Inf], 1e-2)
        @test bins[1] == 29
        @test ɛ < 1e-2
        bins, ɛ, c, itr = LMCLUS.MDL.optquant([NaN], 1e-2)
        @test bins[1] == 1
        @test ɛ < 1e-2
    end

    @testset "Calculations" begin
        # Generate manifold data
        function generate_lm(N::Int, M::Int, C::Int,
                        B::Matrix{Float64},
                        bounds::Matrix{Float64}, θ::Float64,
                        D::Symbol = :Uniform;  σs::Vector{Float64} = ones(N))
            @assert size(bounds) == (N,2) "Define bounds for every dimension"
            @assert size(B) == (N,M) "Define bounds for every dimension"
            manifold = Manifold(M, zeros(N), B, round.(Int, range(1, stop=C, length=C)), θ, 0.0)

            c = 1
            X = zeros(N,C)
            while c <= C
                if D == :Gaussian
                    for i in 1:N
                        R = abs(bounds[i,2]-bounds[i,1])
                        X[i,c] = rand(Normal(0.,σs[i]))
                    end
                else
                    for i in 1:N
                        R = abs(bounds[i,2]-bounds[i,1])
                        X[i,c] = rand()*R - R/2.
                    end
                end
                if distance_to_manifold(X[:,c], B)[1] < θ
                    c += 1
                end
            end

            return X, manifold
        end

        Pm = 64      # Precision encoding constant for model
        Pd = 32      # Precision encoding constant for data
        N = 2        # Space dimension
        M = 1        # Linear manifold dimension
        C = 100      # Size of a LM cluster
        B = Matrix(I,N,M) # Basis vectors
        bounds = hcat(fill(-1.,N), fill(1., N)) # LM cluster bounds
        θ = 0.8      # distance threshold
        σs = [1.0, 0.25] # diag covariances
        ɛ = 1e-2
        tot = 1000
        tol = 1e-8

        Random.seed!(923487298)
        B *= rand()

        Random.seed!(923487298)
        Xg, Mg = generate_lm(N, M, C, B, bounds, θ, :Gausian; σs = σs)
        LMCLUS.adjustbasis!(Mg, Xg)
        @test LMCLUS.hasfullbasis(Mg)

        @test LMCLUS.MDL.calculate(LMCLUS.MDL.Raw, Mg, Xg, Pm, Pd) == Pm*N*C
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.ZeroDim, Mg, Xg, Pm, Pd) == 221

        @test LMCLUS.MDL.calculate(LMCLUS.MDL.Uniform, Mg, Xg, Pm, Pd)  == 5623
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.Gaussian, Mg, Xg, Pm, Pd) == 3534

        # Empirical entropy from optimal quantization
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.Empirical, Mg, Xg, Pm, Pd, ɛ = 1e-2) == 3889
        # Empirical entropy from fixed bin size quantization
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.Empirical, Mg, Xg, Pm, Pd, ɛ = 20.0) == 3808

        # Optimal quantizing
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.OptimalQuant, Mg, Xg, Pm, Pd, ɛ = 1e-2) == 4093
        @test LMCLUS.mdl(Mg, Xg, Pm, Pd, ɛ = 1e-2) == 4093

        Mg.d = 0
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.SizeIndependent, Mg, Xg, Pm, Pd, ɛ = 1e-2) == 9152
        Mg.d = 1
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.SizeIndependent, Mg, Xg, Pm, Pd, ɛ = 1e-2) == 3008

        # Test size dependence
        Random.seed!(923487298)
        Xg, Mg = generate_lm(N, M, 10*C, B, bounds, θ, :Gausian; σs = σs)
        LMCLUS.adjustbasis!(Mg, Xg)
        @test LMCLUS.hasfullbasis(Mg)
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.OptimalQuant, Mg, Xg, Pm, Pd, ɛ = 1e-2) == 39306
        @test LMCLUS.MDL.calculate(LMCLUS.MDL.SizeIndependent, Mg, Xg, Pm, Pd, ɛ = 1e-2) == 3136

    end

end


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_tables.jl - https://github.com/strickek/ReadWriteDlm2.jl

"""
    ReadWriteDlm2.matrix2(table)
Materialize any table source input as a `Matrix{Any}`.
Column names - in Symbol type - are written in first row.
"""
function matrix2(table)
    cols = Tables.Columns(table)
    cnames = Tables.columnnames(table)
    nr = Tables.rowcount(cols) + 1
    nc = length(cnames)
    matrix = Matrix{Any}(undef, nr, nc)
    for (i, col) in enumerate(cols)
        matrix[1, i] = cnames[i]
        matrix[2:end, i] = col
    end
    return matrix
end

"""
    vecofvec(a::Matrix, ct::Vector{Type})
Take the columns of Matrix `a`, convert each column to `Array{T,1} where T`
(`ct` provides Type), and return the columns as `Vector{Array{T,1} where T}`.
"""
function vecofvec(a::Matrix, ct::Vector{<:Type})
    cols = size(a, 2)
    cols == length(ct) || throw(ArgumentError(
        "`ct` length ($(length(ct))) must match number of columns in matrix ($(size(a, 2)))"
        ))
    vov = Vector{Array{T,1} where T}(undef, cols)
    for i = 1:cols
        vov[i] = convert(Vector{ct[i]}, view(a, :,i))
    end
    return vov
end

# MatrixTable - ReadWriteDlm2 outputformat for Tables interface
struct MatrixTable <: Tables.AbstractColumns
    names::Vector{Symbol}
    lookup::Dict{Symbol, Int}
    matrix::Vector{Array{T,1} where T}
end

# Overload Table.table for m::Vector{Vector}
"""
    Tables.table(m::Vector{Array{T,1} where T; [header::Vector{Symbol}])
Wrap an vector containing columns as vectors (`Array{T,1} where T`) in a
`MatrixTable`, which satisfies the Tables.jl interface. This allows accesing
the matrix via `Tables.rows` and `Tables.columns`. An optional keyword argument
`header` can be passed as a `Vector{Symbol}` to be used as the column names.
"""
function Tables.table(m::Vector{Array{T,1} where T};
    header::Vector{Symbol} = [Symbol("Column$i") for i = 1:length(m)]
    )
    length(header) == length(m) || throw(ArgumentError(
    "`header` length must match number of columns in matrix ($(length(m)))"))
    lookup = Dict(nm=>i for (i, nm) in enumerate(header))
    return MatrixTable(header, lookup, m)
end

"""
    mttoarray(m::ReadWriteDlm2.MatrixTable; elt::Type=Union{})
Take Vectors of field matrix in MatrixTable `m` and return data as Array{T,2}.
With given `elt` T=elt, otherwise T is the union of element types in matrix.
"""
function mttoarray(m::ReadWriteDlm2.MatrixTable; elt::Type=Union{})
    ma = ReadWriteDlm2.matrix(m)
    nc = length(ma)
    nr = length(ma[1])
    elt == Union{} ? T1 = Any : T1 = elt
    a = Array{T1}(undef, (nr,nc))
    T = elt
    for i in eachindex(ma)
        T = Union{T, eltype(m[i])}
        a[:,i] = m[i]
    end
    elt != Union{} && return a
    return convert(Array{T,2}, a)
end

# declare that MatrixTable is a table
Tables.istable(::Type{<:MatrixTable}) = true
# getter methods to avoid getproperty clash
names(m::MatrixTable) = getfield(m, :names)
matrix(m::MatrixTable) = getfield(m, :matrix)
lookup(m::MatrixTable) = getfield(m, :lookup)
# schema is column names and types
Tables.schema(m::MatrixTable) = Tables.Schema(names(m), eltype.(matrix(m)))

# column interface
Tables.columnaccess(::Type{<:MatrixTable}) = true
Tables.columns(m::MatrixTable) = m
# required Tables.AbstractColumns object methods
Tables.getcolumn(m::MatrixTable, ::Type{T}, col::Int, nm::Symbol) where {T} = matrix(m)[col]
Tables.getcolumn(m::MatrixTable, nm::Symbol) = matrix(m)[lookup(m)[nm]]
Tables.getcolumn(m::MatrixTable, i::Int) = matrix(m)[i]
Tables.columnnames(m::MatrixTable) = names(m)

# declare that any MatrixTable defines its own `Tables.rows` method
Tables.rowaccess(::Type{<:MatrixTable}) = true
# just return itself, which means MatrixTable must iterate `Tables.AbstractRow`-compatible objects
Tables.rows(m::MatrixTable) = m
# the iteration interface, at a minimum, requires `eltype`, `length`, and `iterate`
# for `MatrixTable` `eltype`, we're going to provide a custom row type
Base.eltype(m::MatrixTable) = MatrixRow
Base.length(m::MatrixTable) = length(matrix(m)[1])
Base.iterate(m::MatrixTable, st=1) = st > length(m) ? nothing : (MatrixRow(st, m), st + 1)

# a custom row type; acts as a "view" into a row of an AbstractMatrix
struct MatrixRow <: Tables.AbstractRow
    row::Int
    source::MatrixTable
end
# required `Tables.AbstractRow` interface methods (same as for `Tables.AbstractColumns` object before)
# but this time, on our custom row type
Tables.getcolumn(m::MatrixRow, ::Type, col::Int, nm::Symbol) =
    getfield(getfield(m, :source), :matrix)[col][getfield(m, :row)]
Tables.getcolumn(m::MatrixRow, i::Int) =
    getfield(getfield(m, :source), :matrix)[i][getfield(m, :row)]
Tables.getcolumn(m::MatrixRow, nm::Symbol) =
    getfield(getfield(m, :source), :matrix)[getfield(getfield(m, :source), :lookup)[nm]][getfield(m, :row)]
Tables.columnnames(m::MatrixRow) = names(getfield(m, :source))


text:
"""
    RandomLM(;ny=100, p=200, sy=.1, nbet=5, mbet=2, sdbet=0, Xbr=0.05)

Randomly generates a linear model.

| Regressors  |   Coefficients      |
|-------------|---------------------|
| X[0:nbet]   |   Normal(mbet,sbet) |
| X[nbet:p]   |   0    nbet         |

If dupl=True, the second half of regressors has collinearity with the first half:

```
X[p/2:p]= X[0:p/2]+ Normal(0,Xbr)
```

## Parameters
  - ny: Number of observations.
  - p: Number of variables.
  - nbet: Number of non-zero coefficients
  - mbet: Mean value of these
  - sdbet: Std of these
  - Xbr: Std for collinearity
  - sy: Variance of u

## Returns
  - X : array, shape(n_samples,n_features)
    Regressor matrix. iid N(0,1) but columns are centered.
    If Xbr>0 each entry in the second half of the matrix
    is the one in the first half plus a noise of amplitude Xbr:
```
X[j,k+p/2]= X[j,k]+ Normal(0,Xbr)         k<p/2.
```
  - betas : array, shape(n_features)
    Coefficients. iid N(mbet,sbet). betas[i]=0 for i>=nbet.
  - y : array, shape(n_samples)
    Responses
"""
function RandomLM(; ny = 100, p = 200, sy = 0.1, nbet = 5, mbet = 2, sdbet = 0, Xbr = 0.05)

    betas = randn(nbet) .* sdbet .+ mbet
    length = nbet÷2:nbet
    betas[length] = -betas[length]
    betas = vcat(betas, zeros(Float64, p - nbet))
    X = randn(ny, p)
    if Xbr > 0
        p0 = p ÷ 2
        X[:, p-p0:p] .= X[:, 0:p0] + Xbr .* randn(ny * p0)
    end
    X .= X .- mean(X, dims = 1)
    y = X * betas + randn(ny) .* sy
    return X, y .- mean(y), betas

end


text:
module SeisNoise

using Dates, DataFrames, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO
using Statistics, StatsBase, Plots, Distributed
using CUDA, Adapt, GPUArrays

# check use of cuda
const use_cuda = Ref(false)
if !CUDA.functional()
  else
    use_cuda[] = true
end

# import types first
include("Types/NoiseData.jl")
include("Types/show.jl")

# import pre and post processing tools
include("distance.jl")
include("io.jl")
include("ArrayFuncs.jl")
include("tools.jl")
include("slicing.jl")
include("filter.jl")
include("phase_shift.jl")
include("stacking.jl")

# import  routines for doin' stuff
include("compute_fft.jl")
include("correlation.jl")
include("rotation.jl")
include("Plotting/plotting.jl")

# out with the old
include("deprecated.jl")

end # module


text:
"""
    This example shows how to create Dendogram visualizations.

According to wikipedia: "A dendrogram is a diagram representing a tree",
but really you just need to see to understand:
"""

import Term: Dendogram

mydendo = Dendogram("awesome", "this", :is, "a", "dendogram!")
print("\n"^2)
print(mydendo)

"""
The first argument si the "head" (or trunk) and all the others are leaves
"""

otherdendo = Dendogram("head", "these", "are", "colorful", "leaves")
print("\n"^2)
print(otherdendo)

"""
Let's say you want to create a hierarchical structure in your dendogram,
just link individual elements together.
"""

import Term.Dendograms: link
print("\n"^2)
print(link(mydendo, otherdendo; title = "[red]superdendo[/red]"))

"""
and so on...
"""
smalldendo = Dendogram("head", [1, 2])

print("\n"^2)
print(
    link(
        smalldendo,
        link(
            smalldendo,
            link(smalldendo, smalldendo; title = "a level");
            title = "another level",
        );
        title = "first level",
    ),
)