text:
"Return the layer name."
getname(layer::FeatureLayer) = GDAL.getname(layer)

"Return the layer geometry type."
getgeomtype(layer::FeatureLayer) = GDAL.getgeomtype(layer)

"Returns the current spatial filter for this layer."
getspatialfilter(layer::FeatureLayer) =
    GDAL.getspatialfilter(layer)

# """
#     OGR_L_GetSpatialRef(OGRLayerH) -> OGRSpatialReferenceH
# Fetch the spatial reference system for this layer.
# ### Parameters
# * **hLayer**: handle to the layer to get the spatial reference from.
# ### Returns
# spatial reference, or NULL if there isn't one.
# """
# function getspatialref{T <: OGRLayerH}(arg1::Ptr{T})
#     checknull(ccall((:OGR_L_GetSpatialRef,libgdal),Ptr{OGRSpatialReferenceH},(Ptr{OGRLayerH},),arg1))
# end

"Set a new spatial filter for the layer, using the geom."
setspatialfilter(layer::FeatureLayer, geom::Geometry) =
    GDAL.setspatialfilter(layer, geom)

"Set a new rectangular spatial filter for the layer."
setspatialfilter(layer::FeatureLayer,
                 xmin::Real, ymin::Real, xmax::Real, ymax::Real) =
    GDAL.setspatialfilterrect(layer, xmin, ymin, xmax, ymax)

"""
Set a new spatial filter.

### Parameters
* **layer**: the layer on which to set the spatial filter.
* **index**: 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.
"""
setspatialfilter(layer::FeatureLayer, index::Integer, geom::Geometry) =
    GDAL.setspatialfilterex(layer, index, geom)

clearspatialfilter(layer::FeatureLayer, index::Integer) =
    GDAL.setspatialfilterex(layer, index, C_NULL)

"""
Set a new rectangular spatial filter.
### Parameters
* **layer**: handle to the layer on which to set the spatial filter.
* **index**: 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.
"""
setspatialfilter(layer::FeatureLayer, index::Integer,
                 xmin::Real, ymin::Real, xmax::Real, ymax::Real) =
    GDAL.setspatialfilterrectex(layer, index, xmin, ymin, xmax, ymax)

"""
    OGR_L_SetAttributeFilter(OGRLayerH,
                             const char *) -> OGRErr
Set a new attribute query.
### Parameters
* **hLayer**: handle to the layer on which attribute query will be executed.
* **pszQuery**: query in restricted SQL WHERE format, or NULL to clear the current query.
### Returns
OGRERR_NONE if successfully installed, or an error code if the query expression is in error, or some other failure occurs.
"""
function setattributefilter(layer::FeatureLayer, query::AbstractString)
    result = GDAL.setattributefilter(layer, query)
    (result != GDAL.OGRERR_NONE) && error("Failed to set a new attribute query")
end

"Reset feature reading to start on the first feature."
resetreading(layer::FeatureLayer) = GDAL.resetreading(layer)

"Fetch the next available feature from this layer."
fetchnextfeature(layer::FeatureLayer) = GDAL.getnextfeature(layer)

function fetchnextfeature(f::Function, args...)
    feature = fetchnextfeature(args...)
    try
        f(feature)
    finally
        destroy(feature)
    end
end

"Move read cursor to the nIndex'th feature in the current resultset.
### Parameters
* **hLayer**: handle to the layer
* **nIndex**: the index indicating how many steps into the result set to seek.
### Returns
OGRERR_NONE on success or an error code.
"""
function setnextbyindex(layer::FeatureLayer, i::Integer)
    result = GDAL.setnextbyindex(layer, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to move the cursor to index $i")
    end
end

"Fetch a feature (now owned by the caller) by its identifier."
fetchfeature(layer::FeatureLayer, i::Integer) =
    GDAL.getfeature(layer, i)

function fetchfeature(f::Function, args...)
    feature = fetchfeature(args...)
    try
        f(feature)
    finally
        destroy(feature)
    end
end

"""
Rewrite an existing feature.

### Returns
OGRERR_NONE if the operation works, otherwise an appropriate error code
(e.g OGRERR_NON_EXISTING_FEATURE if the feature does not exist).
"""
function setfeature(layer::FeatureLayer, feature::Feature)
    result = GDAL.setfeature(layer, feature)
    if result != GDAL.OGRERR_NONE
        error("Failed to set feature.")
    end
end

"""
Create and write a new feature within a layer.

### Returns
OGRERR_NONE on success.
"""
function createfeature(layer::FeatureLayer, feature::Feature)
    result = GDAL.createfeature(layer, feature)
    if result != GDAL.OGRERR_NONE
        error("Failed to create and write feature.")
    end
end

createfeature(layer::FeatureLayer) = createfeature(getlayerdefn(layer))

function createfeature(f::Function, layer::FeatureLayer)
    feature = createfeature(layer)
    try
        f(feature)
        createfeature(layer, feature)
    finally
        destroy(feature)
    end
end

function createfeature(f::Function, featuredefn::FeatureDefn)
    feature = createfeature(featuredefn)
    try
        f(feature)
    finally
        destroy(feature)
    end
end

"""
Delete feature with fid i from layer.

### Returns
OGRERR_NONE if the operation works, otherwise an appropriate error code
(e.g OGRERR_NON_EXISTING_FEATURE if the feature does not exist).
"""
function deletefeature(layer::FeatureLayer, i::Integer)
    result = GDAL.deletefeature(layer, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to delete feature $i.")
    end
end

"""
Set feature geometry.

### Returns
OGRERR_NONE if successful, or OGR_UNSUPPORTED_GEOMETRY_TYPE if the geometry
type is illegal for the OGRFeatureDefn (checking not yet implemented).
"""
function setgeomdirectly(feature::Feature, geom::Geometry)
    result = GDAL.setgeometrydirectly(feature, geom)
    if result != GDAL.OGRERR_NONE
        error("Failed to set feature geometry.")
    end
end

"""
Set feature geometry.

### Parameters
* **feature**: the feature on which new geometry is applied to.
* **geometry**: the new geometry to apply to feature.
### Returns
OGRERR_NONE if successful, or OGR_UNSUPPORTED_GEOMETRY_TYPE if the geometry
type is illegal for the OGRFeatureDefn (checking not yet implemented).
"""
function setgeom(feature::Feature, geom::Geometry)
    result = GDAL.setgeometry(feature, geom)
    if result != GDAL.OGRERR_NONE
        error("Failed to set feature geometry.")
    end
end

"""
Fetch an handle to feature geometry.

### Returns
an handle to internal feature geometry. This object should not be modified.
"""
getgeom(feature::Feature) = GDAL.getgeometryref(feature)

"Test if two features are the same."
equal(feature1::Feature, feature2::Feature) =
    Bool(GDAL.equal(feature1, feature2))

"""
Fetch number of fields on this feature.

This will always be the same as the field count for the OGRFeatureDefn.
"""
nfield(feature::Feature) = GDAL.getfieldcount(feature)

"""
Fetch definition for this field.

### Parameters
* **feature**: the feature on which the field is found.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.

### Returns
an handle to the field definition (from the OGRFeatureDefn). This is an
internal reference, and should not be deleted or modified.
"""
fetchfielddefn(feature::Feature, i::Integer) = GDAL.getfielddefnref(feature, i)

fetchfields(feature::Feature) =
    Dict([getname(fetchfielddefn(feature, i-1)) => fetchfield(feature, i-1)
         for i in 1:nfield(feature)])

fetchfields{T <: Integer}(feature::Feature, indices::UnitRange{T}) =
    Dict([getname(fetchfielddefn(feature, i)) => fetchfield(feature, i)
         for i in indices])

fetchfields{T <: Integer}(feature::Feature, indices::Vector{T}) =
    Dict([getname(fetchfielddefn(feature, i)) => fetchfield(feature, i)
         for i in indices])

fetchfields(feature::Feature, names::Vector{ASCIIString}) =
    Dict([name => fetchfield(feature, getfieldindex(feature, name))
         for name in names])

fetchgeomfields(feature::Feature) =
    Dict([getname(fetchgeomfielddefn(feature, i)) =>
          toWKT(fetchgeomfield(feature, i))
          for i in 1:ngeomfield(feature)])

"""
Fetch the field index given field name.

### Parameters
* **feature**: the feature on which the field is found.
* **name**: the name of the field to search for.

### Returns
the field index, or -1 if no matching field is found.
"""
getfieldindex(feature::Feature, name::AbstractString) =
    GDAL.getfieldindex(feature, name)

"""Test if a field has ever been assigned a value or not.

### Parameters
* **feature**: the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
"""
isfieldset(feature::Feature, i::Integer) =
    Bool(GDAL.isfieldset(feature, i))

"""
Clear a field, marking it as unset.

### Parameters
* **feature**: the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
"""
unsetfield(feature::Feature, i::Integer) = GDAL.unsetfield(feature, i)

# """
#     OGR_F_GetRawFieldRef(OGRFeatureH hFeat,
#                          int iField) -> OGRField *
# Fetch an handle to the internal field value given the index.
# ### Parameters
# * **hFeat**: handle to the feature on which field is found.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# ### Returns
# the returned handle is to an internal data structure, and should not be freed, or modified.
# """
# function getrawfieldref(arg1::Ptr{OGRFeatureH},arg2::Integer)
#     ccall((:OGR_F_GetRawFieldRef,libgdal),Ptr{OGRField},(Ptr{OGRFeatureH},Cint),arg1,arg2)
# end

"""Fetch field value as integer.

### Parameters
* **feature**: the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
"""
asint(feature::Feature, i::Integer) = GDAL.getfieldasinteger(feature, i)

"""Fetch field value as integer 64 bit.

### Parameters
* **feature**: the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
"""
asint64(feature::Feature, i::Integer) = GDAL.getfieldasinteger64(feature, i)

"""Fetch field value as a double.

### Parameters
* **feature**: the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
"""
asdouble(feature::Feature, i::Integer) = GDAL.getfieldasdouble(feature, i)

"""
Fetch field value as a string.

### Parameters
* **feature**: the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
"""
asstring(feature::Feature, i::Integer) = GDAL.getfieldasstring(feature, i)

# """
#     OGR_F_GetFieldAsIntegerList(OGRFeatureH hFeat,
#                                 int iField,
#                                 int * pnCount) -> const int *
# Fetch field value as a list of integers.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **pnCount**: an integer to put the list count (number of integers) into.
# ### Returns
# the field value. This list is internal, and should not be modified, or freed. Its lifetime may be very brief. If *pnCount is zero on return the returned pointer may be NULL or non-NULL.
# """
# function getfieldasintegerlist(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_GetFieldAsIntegerList,libgdal),Ptr{Cint},(Ptr{OGRFeatureH},Cint,Ptr{Cint}),arg1,arg2,arg3)
# end

# """
#     OGR_F_GetFieldAsInteger64List(OGRFeatureH hFeat,
#                                   int iField,
#                                   int * pnCount) -> const GIntBig *
# Fetch field value as a list of 64 bit integers.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **pnCount**: an integer to put the list count (number of integers) into.
# ### Returns
# the field value. This list is internal, and should not be modified, or freed. Its lifetime may be very brief. If *pnCount is zero on return the returned pointer may be NULL or non-NULL.
# """
# function getfieldasinteger64list(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_GetFieldAsInteger64List,libgdal),Ptr{GIntBig},(Ptr{OGRFeatureH},Cint,Ptr{Cint}),arg1,arg2,arg3)
# end

# """
#     OGR_F_GetFieldAsDoubleList(OGRFeatureH hFeat,
#                                int iField,
#                                int * pnCount) -> const double *
# Fetch field value as a list of doubles.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **pnCount**: an integer to put the list count (number of doubles) into.
# ### Returns
# the field value. This list is internal, and should not be modified, or freed. Its lifetime may be very brief. If *pnCount is zero on return the returned pointer may be NULL or non-NULL.
# """
# function getfieldasdoublelist(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_GetFieldAsDoubleList,libgdal),Ptr{Cdouble},(Ptr{OGRFeatureH},Cint,Ptr{Cint}),arg1,arg2,arg3)
# end

# """
#     OGR_F_GetFieldAsStringList(OGRFeatureH hFeat,
#                                int iField) -> char **
# Fetch field value as a list of strings.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# ### Returns
# the field value. This list is internal, and should not be modified, or freed. Its lifetime may be very brief.
# """
# function getfieldasstringlist(arg1::Ptr{OGRFeatureH},arg2::Integer)
#     bytestring(unsafe_load(ccall((:OGR_F_GetFieldAsStringList,libgdal),Ptr{Cstring},(Ptr{OGRFeatureH},Cint),arg1,arg2)))
# end

# """
#     OGR_F_GetFieldAsBinary(OGRFeatureH hFeat,
#                            int iField,
#                            int * pnBytes) -> GByte *
# Fetch field value as binary.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **pnBytes**: location to place count of bytes returned.
# ### Returns
# the field value. This list is internal, and should not be modified, or freed. Its lifetime may be very brief.
# """
# function getfieldasbinary(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_GetFieldAsBinary,libgdal),Ptr{GByte},(Ptr{OGRFeatureH},Cint,Ptr{Cint}),arg1,arg2,arg3)
# end

# """
#     OGR_F_GetFieldAsDateTime(OGRFeatureH hFeat,
#                              int iField,
#                              int * pnYear,
#                              int * pnMonth,
#                              int * pnDay,
#                              int * pnHour,
#                              int * pnMinute,
#                              int * pnSecond,
#                              int * pnTZFlag) -> int
# Fetch field value as date and time.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **pnYear**: (including century)
# * **pnMonth**: (1-12)
# * **pnDay**: (1-31)
# * **pnHour**: (0-23)
# * **pnMinute**: (0-59)
# * **pnSecond**: (0-59)
# * **pnTZFlag**: (0=unknown, 1=localtime, 100=GMT, see data model for details)
# ### Returns
# TRUE on success or FALSE on failure.
# """
# function getfieldasdatetime(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3,arg4,arg5,arg6,arg7,arg8,arg9)
#     ccall((:OGR_F_GetFieldAsDateTime,libgdal),Cint,(Ptr{OGRFeatureH},Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9)
# end

# """
#     OGR_F_GetFieldAsDateTimeEx(OGRFeatureH hFeat,
#                                int iField,
#                                int * pnYear,
#                                int * pnMonth,
#                                int * pnDay,
#                                int * pnHour,
#                                int * pnMinute,
#                                float * pfSecond,
#                                int * pnTZFlag) -> int
# Fetch field value as date and time.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **pnYear**: (including century)
# * **pnMonth**: (1-12)
# * **pnDay**: (1-31)
# * **pnHour**: (0-23)
# * **pnMinute**: (0-59)
# * **pfSecond**: (0-59 with millisecond accuracy)
# * **pnTZFlag**: (0=unknown, 1=localtime, 100=GMT, see data model for details)
# ### Returns
# TRUE on success or FALSE on failure.
# """
# function getfieldasdatetimeex(hFeat::Ptr{OGRFeatureH},iField::Integer,pnYear,pnMonth,pnDay,pnHour,pnMinute,pfSecond,pnTZFlag)
#     ccall((:OGR_F_GetFieldAsDateTimeEx,libgdal),Cint,(Ptr{OGRFeatureH},Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cfloat},Ptr{Cint}),hFeat,iField,pnYear,pnMonth,pnDay,pnHour,pnMinute,pfSecond,pnTZFlag)
# end

function fetchfield(feature::Feature, i::Integer)
    const FETCHFIELD = Dict{GDAL.OGRFieldType, Function}(
                        GDAL.OFTInteger => asint,
                        GDAL.OFTInteger64 => asint64,
                        GDAL.OFTReal => asdouble,
                        GDAL.OFTString => asstring)
    asnothing(feature::Feature, i::Integer) = nothing
    if isfieldset(feature, i)
        # _fielddefn =
        # _fieldname = getname(feature, i)
        _fieldtype = gettype(fetchfielddefn(feature, i))
        _fetchfield = get(FETCHFIELD, _fieldtype, asnothing)
        return _fetchfield(feature, i)
    end
end

fetchfield(feature::Feature, name::AbstractString) =
    fetchfield(feature, getfieldindex(feature, name))

"""
Set field to integer value.

### Parameters
* **feature**: handle to the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
* **value**: the value to assign.
"""
setfield(feature::Feature, i::Integer, value::Integer) =
    GDAL.setfieldinteger(feature, i, value)

"""
Set field to 64 bit integer value.

### Parameters
* **feature**: handle to the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
* **value**: the value to assign.
"""
setfield(feature::Feature, i::Integer, value::Int64) =
    GDAL.setfieldinteger64(feature, i, value)

"""
Set field to double value.

### Parameters
* **feature**: handle to the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
* **value**: the value to assign.
"""
setfield(feature::Feature, i::Integer, value::Cdouble) =
    GDAL.setfielddouble(feature, i, value)

"""
Set field to string value.

### Parameters
* **feature**: handle to the feature that owned the field.
* **i**: the field to fetch, from 0 to GetFieldCount()-1.
* **value**: the value to assign.
"""
setfield(feature::Feature, i::Integer, value::AbstractString) =
    GDAL.setfieldstring(feature, i, value)

# """
#     OGR_F_SetFieldIntegerList(OGRFeatureH hFeat,
#                               int iField,
#                               int nCount,
#                               int * panValues) -> void
# Set field to list of integers value.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **nCount**: the number of values in the list being assigned.
# * **panValues**: the values to assign.
# """
# function setfieldintegerlist(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3::Integer,arg4)
#     ccall((:OGR_F_SetFieldIntegerList,libgdal),Void,(Ptr{OGRFeatureH},Cint,Cint,Ptr{Cint}),arg1,arg2,arg3,arg4)
# end

# """
#     OGR_F_SetFieldInteger64List(OGRFeatureH hFeat,
#                                 int iField,
#                                 int nCount,
#                                 const GIntBig * panValues) -> void
# Set field to list of 64 bit integers value.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **nCount**: the number of values in the list being assigned.
# * **panValues**: the values to assign.
# """
# function setfieldinteger64list(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3::Integer,arg4)
#     ccall((:OGR_F_SetFieldInteger64List,libgdal),Void,(Ptr{OGRFeatureH},Cint,Cint,Ptr{GIntBig}),arg1,arg2,arg3,arg4)
# end

# """
#     OGR_F_SetFieldDoubleList(OGRFeatureH hFeat,
#                              int iField,
#                              int nCount,
#                              double * padfValues) -> void
# Set field to list of doubles value.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **nCount**: the number of values in the list being assigned.
# * **padfValues**: the values to assign.
# """
# function setfielddoublelist(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3::Integer,arg4)
#     ccall((:OGR_F_SetFieldDoubleList,libgdal),Void,(Ptr{OGRFeatureH},Cint,Cint,Ptr{Cdouble}),arg1,arg2,arg3,arg4)
# end

# """
#     OGR_F_SetFieldStringList(OGRFeatureH hFeat,
#                              int iField,
#                              char ** papszValues) -> void
# Set field to list of strings value.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **papszValues**: the values to assign.
# """
# function setfieldstringlist(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_SetFieldStringList,libgdal),Void,(Ptr{OGRFeatureH},Cint,Ptr{Cstring}),arg1,arg2,arg3)
# end

# """
#     OGR_F_SetFieldRaw(OGRFeatureH hFeat,
#                       int iField,
#                       OGRField * psValue) -> void
# Set field.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to fetch, from 0 to GetFieldCount()-1.
# * **psValue**: handle on the value to assign.
# """
# function setfieldraw(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_SetFieldRaw,libgdal),Void,(Ptr{OGRFeatureH},Cint,Ptr{OGRField}),arg1,arg2,arg3)
# end

# """
#     OGR_F_SetFieldBinary(OGRFeatureH hFeat,
#                          int iField,
#                          int nBytes,
#                          GByte * pabyData) -> void
# Set field to binary data.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **nBytes**: the number of bytes in pabyData array.
# * **pabyData**: the data to apply.
# """
# function setfieldbinary(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3::Integer,arg4)
#     ccall((:OGR_F_SetFieldBinary,libgdal),Void,(Ptr{OGRFeatureH},Cint,Cint,Ptr{GByte}),arg1,arg2,arg3,arg4)
# end

# """
#     OGR_F_SetFieldDateTime(OGRFeatureH hFeat,
#                            int iField,
#                            int nYear,
#                            int nMonth,
#                            int nDay,
#                            int nHour,
#                            int nMinute,
#                            int nSecond,
#                            int nTZFlag) -> void
# Set field to datetime.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **nYear**: (including century)
# * **nMonth**: (1-12)
# * **nDay**: (1-31)
# * **nHour**: (0-23)
# * **nMinute**: (0-59)
# * **nSecond**: (0-59)
# * **nTZFlag**: (0=unknown, 1=localtime, 100=GMT, see data model for details)
# """
# function setfielddatetime(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3::Integer,arg4::Integer,arg5::Integer,arg6::Integer,arg7::Integer,arg8::Integer,arg9::Integer)
#     ccall((:OGR_F_SetFieldDateTime,libgdal),Void,(Ptr{OGRFeatureH},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9)
# end

# """
#     OGR_F_SetFieldDateTimeEx(OGRFeatureH hFeat,
#                              int iField,
#                              int nYear,
#                              int nMonth,
#                              int nDay,
#                              int nHour,
#                              int nMinute,
#                              float fSecond,
#                              int nTZFlag) -> void
# Set field to datetime.
# ### Parameters
# * **hFeat**: handle to the feature that owned the field.
# * **iField**: the field to set, from 0 to GetFieldCount()-1.
# * **nYear**: (including century)
# * **nMonth**: (1-12)
# * **nDay**: (1-31)
# * **nHour**: (0-23)
# * **nMinute**: (0-59)
# * **fSecond**: (0-59, with millisecond accuracy)
# * **nTZFlag**: (0=unknown, 1=localtime, 100=GMT, see data model for details)
# """
# function setfielddatetimeex(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3::Integer,arg4::Integer,arg5::Integer,arg6::Integer,arg7::Integer,arg8::Cfloat,arg9::Integer)
#     ccall((:OGR_F_SetFieldDateTimeEx,libgdal),Void,(Ptr{OGRFeatureH},Cint,Cint,Cint,Cint,Cint,Cint,Cfloat,Cint),arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9)
# end

"""
Fetch number of geometry fields on this feature.

This will always be the same as the geometry field count for OGRFeatureDefn.
"""
ngeomfield(feature::Feature) = GDAL.getgeomfieldcount(feature)

"""
Fetch definition for this geometry field.

### Parameters
* **feature**: the feature on which the field is found.
* **i**: the field to fetch, from 0 to GetGeomFieldCount()-1.

### Returns
The field definition (from the OGRFeatureDefn). This is an
internal reference, and should not be deleted or modified.
"""
fetchgeomfielddefn(feature::Feature, i::Integer) =
    GDAL.getgeomfielddefnref(feature, i)

"""
Fetch the geometry field index given geometry field name.

### Parameters
* **feature**: the feature on which the geometry field is found.
* **name**: the name of the geometry field to search for.

### Returns
the geometry field index, or -1 if no matching geometry field is found.
"""
getgeomfieldindex(feature::Feature, name::AbstractString) =
    GDAL.getgeomfieldindex(feature, i)

"""
Fetch the feature geometry.

### Parameters
* **feature**: handle to the feature to get geometry from.
* **i**: geometry field to get.

### Returns
an internal feature geometry. This object should not be modified.
"""
fetchgeomfield(feature::Feature, i::Integer) =
    GDAL.getgeomfieldref(feature, i)

"""
Set feature geometry of a specified geometry field.

### Parameters
* **feature**: the feature on which to apply the geometry.
* **i**: geometry field to set.
* **geom**: the new geometry to apply to feature.

### Returns
OGRERR_NONE if successful, or OGRERR_FAILURE if the index is invalid,
or OGR_UNSUPPORTED_GEOMETRY_TYPE if the geometry type is illegal for
the OGRFeatureDefn (checking not yet implemented).
"""
function setgeomfielddirectly(feature::Feature, i::Integer, geom::Geometry)
    result = GDAL.setgeomfielddirectly(feature, i, geom)
    if result != GDAL.OGRERR_NONE
        error("Failed to set feature geometry")
    end
end

"""
Set feature geometry of a specified geometry field.

### Parameters
* **feature**: the feature on which to apply the geometry.
* **i**: geometry field to set.
* **geom**: the new geometry to apply to feature.

### Returns
OGRERR_NONE if successful, or OGRERR_FAILURE if the index is invalid,
or OGR_UNSUPPORTED_GEOMETRY_TYPE if the geometry type is illegal for
the OGRFeatureDefn (checking not yet implemented).
"""
function setgeomfield(feature::Feature, i::Integer, geom::Geometry)
    GDAL.setgeomfield(feature, i, geom)
end

"""
Get feature identifier.

### Returns
feature id or OGRNullFID if none has been assigned.
"""
getfid(feature::Feature) = GDAL.getfid(feature)

"""
Set the feature identifier.

### Parameters
* **feature**: handle to the feature to set the feature id to.
* **i**: the new feature identifier value to assign.

### Returns
On success OGRERR_NONE, or on failure some other value.
"""
function setfid(feature::Feature, i::Integer)
    result = GDAL.setfid(feature, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to set FID")
    end
end

# """
#     OGR_F_DumpReadable(OGRFeatureH hFeat,
#                        FILE * fpOut) -> void
# Dump this feature in a human readable form.
# ### Parameters
# * **hFeat**: handle to the feature to dump.
# * **fpOut**: the stream to write to, such as strout.
# """
# function dumpreadable(arg1::Ptr{OGRFeatureH},arg2)
#     ccall((:OGR_F_DumpReadable,libgdal),Void,(Ptr{OGRFeatureH},Ptr{FILE}),arg1,arg2)
# end

"""
Set one feature from another.

### Parameters
* **feature1**: handle to the feature to set to.
* **feature2**: handle to the feature from which geometry, and field values
    will be copied.
* **forgiving**: TRUE if the operation should continue despite lacking output
    fields matching some of the source fields.

### Returns
OGRERR_NONE if the operation succeeds, even if some values are not transferred,
otherwise an error code.
"""
function setfrom(feature1::Feature, feature2::Feature, forgiving::Bool=false)
    result = GDAL.setfrom(feature1, feature2, forgiving)
    if result != GDAL.OGRERR_NONE
        error("Failed to set feature")
    end
end

"""
Set one feature from another.

### Parameters
* **feature1**: the feature to set to.
* **feature2**: the feature from which geometry and field values will be copied
* **indices**: indices of the destination feature's fields stored at the
    corresponding index of the source feature's fields. A value of -1 should be
    used to ignore the source's field. The array should not be NULL and be as
    long as the number of fields in the source feature.
* **bForgiving**: TRUE if the operation should continue despite lacking output
    fields matching some of the source fields.

### Returns
OGRERR_NONE if the operation succeeds, even if some values are not transferred,
otherwise an error code.
"""
function setfrom(feature1::Feature, feature2::Feature, indices::Vector{Cint},
                 forgiving::Bool=false)
    result = GDAL.setfromwithmap(feature1, feature2, forgiving, pointer(indices))
    if result != GDAL.OGRERR_NONE
        error("Failed to set feature with map")
    end
end

"Fetch style string for this feature."
getstylestring(feature::Feature) = GDAL.getstylestring(feature)

"Set feature style string."
setstylestring(feature::Feature, style::AbstractString) =
    GDAL.setstylestring(feature, style)

"Set feature style string."
setstylestringdirectly(feature::Feature, style::AbstractString) =
    GDAL.setstylestringdirectly(feature, style)

# """
#     OGR_F_GetStyleTable(OGRFeatureH hFeat) -> OGRStyleTableH
# """
# function getstyletable(arg1::Ptr{OGRFeatureH})
#     checknull(ccall((:OGR_F_GetStyleTable,libgdal),Ptr{OGRStyleTableH},(Ptr{OGRFeatureH},),arg1))
# end

# """
#     OGR_F_SetStyleTableDirectly(OGRFeatureH hFeat,
#                                 OGRStyleTableH hStyleTable) -> void
# """
# function setstyletabledirectly(arg1::Ptr{OGRFeatureH},arg2::Ptr{OGRStyleTableH})
#     ccall((:OGR_F_SetStyleTableDirectly,libgdal),Void,(Ptr{OGRFeatureH},Ptr{OGRStyleTableH}),arg1,arg2)
# end

# """
#     OGR_F_SetStyleTable(OGRFeatureH hFeat,
#                         OGRStyleTableH hStyleTable) -> void
# """
# function setstyletable(arg1::Ptr{OGRFeatureH},arg2::Ptr{OGRStyleTableH})
#     ccall((:OGR_F_SetStyleTable,libgdal),Void,(Ptr{OGRFeatureH},Ptr{OGRStyleTableH}),arg1,arg2)
# end

"Returns the native data for the feature."
getnativedata(feature::Feature) = GDAL.getnativedata(feature)

"Sets the native data for the feature."
setnativedata(feature::Feature, data::AbstractString) =
    GDAL.setnativedata(feature, data)

"Returns the native media type for the feature."
getmediatype(feature::Feature) = GDAL.getnativemediatype(feature)

"Sets the native media type for the feature."
setmediatype(feature::Feature, mediatype::AbstractString) =
    GDAL.setnativemediatype(feature, mediatype)

"""
Fill unset fields with default values that might be defined.

### Parameters
* **feature**: handle to the feature.
* **notnull**: if we should fill only unset fields with a not-null constraint.
* **papszOptions**: unused currently. Must be set to NULL.
"""
fillunset(feature::Feature, notnull::Bool=true) =
    GDAL.fillunsetwithdefault(feature, notnull)

"""
Validate that a feature meets constraints of its schema.

### Parameters
* **feature**: handle to the feature to validate.
* **flags**: OGR_F_VAL_ALL or combination of OGR_F_VAL_NULL,
    OGR_F_VAL_GEOM_TYPE, OGR_F_VAL_WIDTH and OGR_F_VAL_ALLOW_NULL_WHEN_DEFAULT
    with '|' operator
* **emiterror**: TRUE if a CPLError() must be emitted when a check fails

### Returns
TRUE if all enabled validation tests pass.
"""
validate(feature::Feature, flags::Integer, emiterror::Bool) =
    Bool(GDAL.validate(feature, flags, emiterror))

"Fetch the schema information for this layer."
getlayerdefn(layer::FeatureLayer) = GDAL.getlayerdefn(layer)

"""
Find the index of field in a layer.

### Returns
field index, or -1 if the field doesn't exist
"""
findfieldindex(layer::FeatureLayer, field::AbstractString, exactmatch::Bool) =
    GDAL.findfieldindex(layer, field, exactmatch)

"""
Fetch the feature count in this layer.

### Parameters
* **layer**: handle to the layer that owned the features.
* **force**: Flag indicating whether the count should be computed even if it
    is expensive.

### Returns
feature count, -1 if count not known.
"""
nfeature(layer::FeatureLayer, force::Bool=false) =
    GDAL.getfeaturecount(layer, force)

# """
#     OGR_L_GetExtent(OGRLayerH,
#                     OGREnvelope *,
#                     int) -> OGRErr
# Fetch the extent of this layer.
# ### Parameters
# * **hLayer**: handle to the layer from which to get extent.
# * **psExtent**: the structure in which the extent value will be returned.
# * **bForce**: Flag indicating whether the extent should be computed even if it is expensive.
# ### Returns
# OGRERR_NONE on success, OGRERR_FAILURE if extent not known.
# """
# function getextent{T <: OGRLayerH}(arg1::Ptr{T},arg2,arg3::Integer)
#     ccall((:OGR_L_GetExtent,libgdal),OGRErr,(Ptr{OGRLayerH},Ptr{OGREnvelope},Cint),arg1,arg2,arg3)
# end

# """
#     OGR_L_GetExtentEx(OGRLayerH,
#                       int iGeomField,
#                       OGREnvelope * psExtent,
#                       int bForce) -> OGRErr
# Fetch the extent of this layer, on the specified geometry field.
# ### Parameters
# * **hLayer**: handle to the layer from which to get extent.
# * **iGeomField**: the index of the geometry field on which to compute the extent.
# * **psExtent**: the structure in which the extent value will be returned.
# * **bForce**: Flag indicating whether the extent should be computed even if it is expensive.
# ### Returns
# OGRERR_NONE on success, OGRERR_FAILURE if extent not known.
# """
# function getextentex{T <: OGRLayerH}(arg1::Ptr{T},iGeomField::Integer,psExtent,bForce::Integer)
#     ccall((:OGR_L_GetExtentEx,libgdal),OGRErr,(Ptr{OGRLayerH},Cint,Ptr{OGREnvelope},Cint),arg1,iGeomField,psExtent,bForce)
# end

# """
#     OGR_L_TestCapability(OGRLayerH,
#                          const char *) -> int
# Test if this layer supported the named capability.
# ### Parameters
# * **hLayer**: handle to the layer to get the capability from.
# * **pszCap**: the name of the capability to test.
# ### Returns
# TRUE if the layer has the requested capability, or FALSE otherwise. OGRLayers will return FALSE for any unrecognized capabilities.
# """
# function testcapability{T <: OGRLayerH}(arg1::Ptr{T},arg2)
#     ccall((:OGR_L_TestCapability,libgdal),Cint,(Ptr{OGRLayerH},Cstring),arg1,arg2)
# end

"""
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.

### Returns
OGRERR_NONE on success.
"""
function createfield(layer::FeatureLayer, field::FieldDefn,
                     approx::Bool = false)
    result = GDAL.createfield(layer, field, approx)
    if result != GDAL.OGRERR_NONE
        error("Failed to create new field")
    end
end

"""
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.

### Returns
OGRERR_NONE on success.
"""
function creategeomfield(layer::FeatureLayer, field::GeomFieldDefn,
                     approx::Bool = false)
    result = GDAL.creategeomfield(layer, field, approx)
    if result != GDAL.OGRERR_NONE
        error("Failed to create new field")
    end
end

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

### Parameters
* **layer**: handle to the layer.
* **i**: index of the field to delete.

### Returns
OGRERR_NONE on success.
"""
function deletefield(layer::FeatureLayer, i::Integer)
    result = GDAL.deletefield(layer, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to delete field")
    end
end

# """
#     OGR_L_ReorderFields(OGRLayerH,
#                         int * panMap) -> OGRErr
# Reorder all the fields of a layer.
# ### Parameters
# * **hLayer**: handle to the layer.
# * **panMap**: an array of GetLayerDefn()->OGRFeatureDefn::GetFieldCount() elements which is a permutation of [0, GetLayerDefn()->OGRFeatureDefn::GetFieldCount()-1].
# ### Returns
# OGRERR_NONE on success.
# """
# function reorderfields{T <: OGRLayerH}(arg1::Ptr{T},panMap)
#     ccall((:OGR_L_ReorderFields,libgdal),OGRErr,(Ptr{OGRLayerH},Ptr{Cint}),arg1,panMap)
# end

# """
#     OGR_L_ReorderField(OGRLayerH,
#                        int iOldFieldPos,
#                        int iNewFieldPos) -> OGRErr
# Reorder an existing field on a layer.
# ### Parameters
# * **hLayer**: handle to the layer.
# * **iOldFieldPos**: previous position of the field to move. Must be in the range [0,GetFieldCount()-1].
# * **iNewFieldPos**: new position of the field to move. Must be in the range [0,GetFieldCount()-1].
# ### Returns
# OGRERR_NONE on success.
# """
# function reorderfield{T <: OGRLayerH}(arg1::Ptr{T},iOldFieldPos::Integer,iNewFieldPos::Integer)
#     ccall((:OGR_L_ReorderField,libgdal),OGRErr,(Ptr{OGRLayerH},Cint,Cint),arg1,iOldFieldPos,iNewFieldPos)
# end

# """
#     OGR_L_AlterFieldDefn(OGRLayerH,
#                          int iField,
#                          OGRFieldDefnH hNewFieldDefn,
#                          int nFlags) -> OGRErr
# Alter the definition of an existing field on a layer.
# ### Parameters
# * **hLayer**: handle to the layer.
# * **iField**: index of the field whose definition must be altered.
# * **hNewFieldDefn**: new field definition
# * **nFlags**: 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.
# ### Returns
# OGRERR_NONE on success.
# """
# function alterfielddefn{T <: OGRLayerH}(arg1::Ptr{T},iField::Integer,hNewFieldDefn::Ptr{OGRFieldDefnH},nFlags::Integer)
#     ccall((:OGR_L_AlterFieldDefn,libgdal),OGRErr,(Ptr{OGRLayerH},Cint,Ptr{OGRFieldDefnH},Cint),arg1,iField,hNewFieldDefn,nFlags)
# end

# """
#     OGR_L_StartTransaction(OGRLayerH) -> OGRErr
# For datasources which support transactions, StartTransaction creates a transaction.
# ### Parameters
# * **hLayer**: handle to the layer
# ### Returns
# OGRERR_NONE on success.
# """
# function starttransaction{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_StartTransaction,libgdal),OGRErr,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_CommitTransaction(OGRLayerH) -> OGRErr
# For datasources which support transactions, CommitTransaction commits a transaction.
# ### Parameters
# * **hLayer**: handle to the layer
# ### Returns
# OGRERR_NONE on success.
# """
# function committransaction{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_CommitTransaction,libgdal),OGRErr,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_RollbackTransaction(OGRLayerH) -> OGRErr
# For datasources which support transactions, RollbackTransaction will roll back a datasource to its state before the start of the current transaction.
# ### Parameters
# * **hLayer**: handle to the layer
# ### Returns
# OGRERR_NONE on success.
# """
# function rollbacktransaction{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_RollbackTransaction,libgdal),OGRErr,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_Reference(OGRLayerH hLayer) -> int
# """
# function reference{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_Reference,libgdal),Cint,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_Dereference(OGRLayerH hLayer) -> int
# """
# function dereference{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_Dereference,libgdal),Cint,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_GetRefCount(OGRLayerH hLayer) -> int
# """
# function getrefcount{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_GetRefCount,libgdal),Cint,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_SyncToDisk(OGRLayerH) -> OGRErr
# Flush pending changes to disk.
# ### Parameters
# * **hLayer**: handle to the layer
# ### Returns
# OGRERR_NONE if no error occurs (even if nothing is done) or an error code.
# """
# function synctodisk{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_SyncToDisk,libgdal),OGRErr,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_GetFeaturesRead(OGRLayerH hLayer) -> GIntBig
# """
# function getfeaturesread{T <: OGRLayerH}(arg1::Ptr{T})
#     ccall((:OGR_L_GetFeaturesRead,libgdal),GIntBig,(Ptr{OGRLayerH},),arg1)
# end

# """
#     OGR_L_GetFIDColumn(OGRLayerH) -> const char *
# This method returns the name of the underlying database column being used as the FID column, or "" if not supported.
# ### Parameters
# * **hLayer**: handle to the layer
# ### Returns
# fid column name.
# """
# function getfidcolumn{T <: OGRLayerH}(arg1::Ptr{T})
#     bytestring(ccall((:OGR_L_GetFIDColumn,libgdal),Cstring,(Ptr{OGRLayerH},),arg1))
# end

"""
This method returns the name of the underlying database column being used as
the geometry column, or "" if not supported.
"""
getgeomcolumn(layer::FeatureLayer) = GDAL.getgeometrycolumn(layer)

"""
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.

### Returns
OGRERR_NONE if all field names have been resolved (even if the driver does not
support this method)
"""
function setignoredfields(layer::FeatureLayer, fieldnames)
    GDAL.setignoredfields(layer, Ptr{Ptr{UInt8}}(pointer(fieldnames)))
end

"""
Intersection 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 above 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.
"""
function intersection(input::FeatureLayer, method::FeatureLayer,
                      result::FeatureLayer)
    result = GDAL.intersection(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to compute the intersection of the two layers")
    end
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 above 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.
"""
function union(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer)
    result = GDAL.union(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to compute the union of the two layers")
    end
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 above 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.
"""
function symdifference(input::FeatureLayer, method::FeatureLayer,
                       result::FeatureLayer)
    result = GDAL.symdifference(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to compute the sym difference of the two layers")
    end
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 above 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.
"""
function identity(input::FeatureLayer, method::FeatureLayer,
                  result::FeatureLayer)
    result = GDAL.identity(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to compute the identity of the two layers")
    end
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 above 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.
"""
function update(input::FeatureLayer, method::FeatureLayer,
                result::FeatureLayer)
    result = GDAL.update(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to update the layer")
    end
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 above 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.
"""
function clip(input::FeatureLayer, method::FeatureLayer,
              result::FeatureLayer)
    result = GDAL.clip(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to clip the input layer")
    end
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 above 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.
"""
function erase(input::FeatureLayer, method::FeatureLayer,
               result::FeatureLayer)
    result = GDAL.erase(input, method, result, C_NULL, C_NULL, C_NULL)
    if result != GDAL.OGRERR_NONE
        error("Failed to remove areas that are covered by the method layer.")
    end
end


text:
using Convex

type DiffVariable{T,D}
    v::Array{T,D}
    dv::AbstractExpr
end

DiffVariable{T,D}(v::Array{T,D}) = DiffVariable{T,D}(v,Variable(size(v)...))

import Base.transpose
transpose(v::DiffVariable) = DiffVariable(v.v',v.dv')

import Base.*
*(a::DiffVariable, b::DiffVariable) = (a.v + a.dv)*b.v + a.v*b.dv
*(a::DiffVariable, b::Array) = (a.v + a.dv)*b
*(a::Array, b::DiffVariable) = a*(b.v + b.dv)

export DiffVariable, transpose, *

ext: jl


text:
export TrivilGate, Wait

abstract type TrivilGate{N} <: PrimitiveBlock{N} end

mat(d::TrivilGate{N}) where N = IMatrix{1<<N}()
apply!(reg::DefaultRegister, d::TrivilGate) = reg
Base.adjoint(g::TrivilGate) = g

"""
    Wait{N, T} <: TrivilGate{N}
    Wait{N}(t)

Wait the experimental signals for time `t` (empty run).
"""
struct Wait{N, T} <: TrivilGate{N}
    t::T
    Wait{N}(t::T) where {N,T} = new{N, T}(t)
end
YaoBlocks.print_block(io::IO, d::Wait) = print(io, "Wait → $(d.t)")

export EchoBlock
struct EchoBlock{N,OT<:IO} <: PrimitiveBlock{N}
    sym::Symbol
    io::OT
end

EchoBlock(nbits::Int, sym::Symbol) = EchoBlock{nbits, typeof(stdout)}(sym, stdout)
EchoBlock(nbits::Int) = EchoBlock(nbits, :ECHO)
EchoBlock() = nbits->EchoBlock(nbits)

Yao.apply!(reg::AbstractRegister, ec::EchoBlock) = (println(ec.io, "apply!(::$(typeof(reg)), $(ec.sym))"); reg)
Yao.mat(::Type{T}, ec::EchoBlock{N}) where {T,N} = (println(ec.io, "mat(::Type{$T}, $(ec.sym))"); IMatrix{1<<N}())
Yao.ishermitian(ec::EchoBlock{N}) where N = (println(ec.io, "ishermitian($(ec.sym))"); true)
Yao.isunitary(ec::EchoBlock{N}) where N = (println(ec.io, "isunitary($(ec.sym))"); true)
Yao.isreflexive(ec::EchoBlock{N}) where N = (println(ec.io, "isreflexive($(ec.sym))"); true)
Yao.getiparams(ec::EchoBlock{N}) where N = (println(ec.io, "getiparams($(ec.sym))");())
Yao.print_block(io::IO, ec::EchoBlock) = print(io, "EchoBlock($(ec.sym))")


text:
module TimeManager
export now
export programStartTime
export currentTime
export programTime
export programTimeStr
export OnTime

now() = Dates.time()
programStartTime = now()

currentTime(startTime::Real) = (now() - startTime)
programTime() = currentTime(programStartTime)
programTimeStr() = @sprintf("%.3f", programTime())

function OnTime(milisec::Number, prevTime::Base.RefValue{Float64})
	time=now()
	r=(time - Base.getindex(prevTime)) >= milisec
	if r Base.setindex!(prevTime, time) end
	r
end

end

module LoggerManager

import Base.print, Base.println, Base.info, Base.warn, Base.error

export @printf # cannot replace Base.@printf

using Suppressor
using TimeManager

export LOGGER_OUT
export LOGGER_ERROR

export print
export printf
export println
export info
export warn
export error

LOGGER_OUT = "out.log"
LOGGER_ERROR = "error.log"

function __init__()
	open(LOGGER_OUT, "w+")
	open(LOGGER_ERROR, "w+")
end

@suppress begin print(xs...) = open(f -> (print(f, xs...); print(STDOUT, xs...)), LOGGER_OUT, "a+") end
macro printf(xs...) open(f -> (:(Base.@printf(f, $xs...)); :(Base.@printf(:STDOUT, $xs...))),LOGGER_OUT, "a+") end
@suppress begin println(xs...) = open(f -> (println(f, programTimeStr(), " ", xs...); println(STDOUT, xs...)), LOGGER_OUT, "a+") end

@suppress begin info(xs...) = open(f -> (info(f, programTimeStr(), " ", xs...); info(STDOUT, xs...)), LOGGER_OUT, "a+") end
@suppress begin error(xs...) = open(f -> (error(f, programTimeStr()," ", xs...); error(STDERR, xs...)), LOGGER_ERROR, "a+") end
@suppress begin warn(xs...) = open(f -> (warn(f, programTimeStr()," ", xs...); warn(STDERR, xs...)), LOGGER_ERROR, "a+") end

function logException(ex::Exception, title="")
	time=programTimeStr()
  m=msg(:Exception, title, "", "See '$LOGGER_ERROR' for more info.")
	open(f -> println(f, m), LOGGER_OUT, "a+")
	open(function(f)
			print(f, msg(:Exception, title, "", ":"))
			Base.showerror(f, ex, catch_backtrace())
			println(f, "\n----------")
	end, LOGGER_ERROR, "a+")

  println(STDERR, "\x1b[91m", m, "\x1b[0m")
  print("Program will close in 3 sec")
  for i=1:3 Libc.systemsleep(1); print(".") end
  println()
  quit()
end

log(title, f::Function, args...) = try; f(args...); catch err; logException(err,title); end

msg(msg::String,lb=true) = string(programTimeStr(), " ",  msg, lb?"\n":"")

function msg(mode::Union{Symbol,String}, title::Union{Symbol,String}, name::Union{Symbol,String}, msg::String, lb=true)
  title=title!=""?string("(",title,") "):""
  mode=mode!=""?string("[",mode,"] "):""
  name=name!=""?string("'",name,"' "):""
  lb=lb?"\n":""
  string(programTimeStr(), " $mode$title$name$msg$lb")
end

end #LoggerManager

module ThreadManager

export ThreadID
export tcall

Mutex = Threads.Mutex()

ThreadID() = Threads.threadid()

function tcall(m::Threads.Mutex, f::Function, args...)
	Threads.lock(m)
	f(args...)
	Threads.unlock(m)
end

tcall(f::Function, args...) = tcall(Mutex, f, args...)

function run(a::Array{Function,1})
	global Mutex = Threads.Mutex()
	max=Threads.nthreads()
	i=0; Threads.@threads for start in a
		start()
		if i >= max break end; i+=1
  end
end

end # ThreadManager

tinfo(args...) = ThreadManager.tcall(() -> info(args...))
twarn(args...) = ThreadManager.tcall(() -> warn(args...))
terror(args...) = ThreadManager.tcall(() -> error(args...))

tprint(args...) = ThreadManager.tcall(() -> print(args...))
tprintln(args...) = ThreadManager.tcall(() -> println(args...))

tpush(args...) = ThreadManager.tcall(() -> push!(args...))

tsleep(sec) = Libc.systemsleep(sec)

tcheck(f::Function, args...) = ThreadManager.tcall(()->LoggerManager.log(string("T",ThreadManager.ThreadID()),f, args...))

tinit(p::Tuple{String,Function}) = function()
	tprint(LoggerManager.msg(:Debug, string("T",ThreadManager.ThreadID()), p[1], "start"))
	tcheck(p[2])
end

Messages = String[] #shared object among threads
clearMessages() = (global Messages = Array{String,1}())
getMessages() = Messages

showMessages() = ThreadManager.tcall(function()
					if length(Messages) > 0
						println("---[Messages]---")
						for msg in Messages print(msg) end
						clearMessages()
						println("----------------")
					end
				end)

pushMsg(mode, name, msg) = tpush(Messages, LoggerManager.msg(mode , string("T",ThreadManager.ThreadID()), name, msg))
pushMsg(name, msg)   = pushMsg("", name, msg)
pushInfo(name, msg)  = pushMsg(:Info, name, msg)
pushDebug(name, msg) = pushMsg(:Debug, name, msg)
pushError(name, msg) = pushMsg(:Error, name, msg)
pushWarn(name, msg)  = pushMsg(:Warning, name, msg)

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
#///////////////////////////////////////
#// File Name: runtests.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: Test script for RiskSensitiveSAC package
#///////////////////////////////////////

using Test
using RiskSensitiveSAC
using PyCall

py"""
import torch
"""

@testset "RiskSensitiveSAC Unit Tests" begin
@info "Executing Type Definition Test"
include("type_definition_test.jl");
@info "Executing State Transition Test"
include("state_transition_test.jl");
@info "Executing Cost Test"
include("cost_test.jl")
@info "Executing Cost CUDA Test"
include("cost_gpu_test.jl")
@info "Executing Scene Loader Test"
include("scene_loader_test.jl")
@info "Executing Predictor Test"
include("predictor_test.jl")
@info "Executing Forward-Backward Simulation Test"
include("forward_backward_simulation_test.jl")
@info "Executing RiskSensitiveSAC Controller Test"
include("rs_sac_controller_test.jl")
@info "Executing BIC Test"
include("bic_test.jl")
@info "Executing CrowdNav Test"
include("crowd_nav_controller_test.jl")
@info "Executing Utils Test"
include("utils_test.jl")
end


text:
using JLD2, Test

for ioty in [JLD2.MmapIO, IOStream]
    fn = joinpath(tempdir(), "test.jld")

    # Ensure we can handle cases where the same file is open multiple times
    f = jldopen(fn, true, true, true, ioty)
    write(f, "hello", "world")
    @test_throws ArgumentError jldopen(fn, true, true, true, ioty)
    f2 = jldopen(fn, true, false, false, ioty)
    write(f2, "this", "works")
    close(f2)
    write(f, "still", "works")
    close(f)

    f = jldopen(fn, false, false, false, ioty)
    @test read(f, "hello") == "world"
    @test read(f, "this") == "works"
    @test read(f, "still") == "works"
    @test_throws ArgumentError jldopen(fn, true, false, false, ioty)
    @test_throws ArgumentError jldopen(fn, false, false, false, ioty; compress=true)
    @test_throws ArgumentError jldopen(fn, false, false, false, ioty; mmaparrays=true)
    close(f)

    # Make sure that the file gets written properly regardless of the order in which the
    # IO and JLDFile are finalized
    for order in [[1, 2], [2, 1]]
        f = jldopen(fn, true, true, true, ioty)
        write(f, "x", 1:10)

        objs = isa(f, JLD2.JLDFile{IOStream}) ? [f, f.io] : [f, f.io.f]
        for i in order
            finalize(objs[i])
        end

        f = jldopen(fn, false, false, false, ioty)
        @test read(f, "x") == 1:10
        close(f)
    end

    cmd = """
    using JLD2
    f = jldopen($(repr(fn)), true, true, true, $(ioty))
    write(f, "x", 1:10)
    """
    run(`$(Base.julia_cmd()) -e $cmd`)

    f = jldopen(fn, false, false, false, ioty)
    @test read(f, "x") == 1:10
    close(f)
end


text:
# ==========================================================================
# General low-level interface to fit Kernel-SVMs

@inline function fit(spec::SVMSpec, X::AbstractMatrix, y⃗::AbstractVector;
                     solver::SvmSolver = DualCD(),
                     loss::MarginBasedLoss = L2HingeLoss(),
                     nargs...)
  encoding = SignedClassEncoding(y⃗)
  t = labelencode(encoding, y⃗)
  model = fit(spec, X, t; solver = solver, loss = loss, nargs...)
  SVC(model, encoding)
end

# --------------------------------------------------------------------------

@inline function fit{T<:Real}(spec::SVMSpec, X::AbstractMatrix, y⃗::AbstractVector{T};
                              solver::SvmSolver = DualCD(), # TODO: This must not be DualCD for Kernel SVMs
                              nargs...)
  fit(spec, X, y⃗, solver; nargs...)
end

# ==========================================================================
# General low-level interface to fit Linear SVMs
# For them it is possible to request the primal or dual solution

@inline function fit{T<:Real}(spec::LinearSVMSpec, X::AbstractMatrix, y⃗::AbstractVector{T};
                              solver::SvmSolver = DualCD(),
                              dual::Bool = false,
                              nargs...)
  native_model = fit(spec, X, y⃗, solver; dual = dual, nargs...)
  # return the native_fit if it is the desired one
  dual == (typeof(native_model) <: DualSVM) && return native_model
  try # convert to desired format
    if dual == true
      info("Solver was unable to provide the dual solution. Attempting automatic conversion...")
      return convert(DualSVM, native_model)
    elseif dual == false
      info("Solver was unable to provide the primal solution. Attempting automatic conversion...")
      return convert(PrimalSVM, native_model)
    end
  catch
    warn("Automatic conversion failed! Returning solution natively provided by the solver.")
    return native_model
  end
end

# ==========================================================================
# This in-between function will compute the appropriate prediction function

function fit{T<:Real}(spec::SVMSpec, X::AbstractMatrix, y⃗::AbstractVector{T}, solver::SvmSolver;
                      bias::Real = 1.,
                      dual::Bool = true,
                      nargs...)
  k, l = size(X)
  l == length(y⃗) || throw(DimensionMismatch("X and y⃗ have to have the same number of observations (columns)"))
  predtype = prediction_type(spec.loss, y⃗, bias)
  problem = dual ? DualProblem() : PrimalProblem()
  solution = fit(spec, X, y⃗, solver, predtype, problem; nargs...)
  svmModel(spec, solution, LinearPredictor(bias), X, y⃗)
end

# ==========================================================================
# Special method overload to support the "do" notation for the callback
#   model = fit(spec, X, Y) do t, alpha, v, g
#     println("Iteration $t: $v")
#   end

@inline function fit(cb::Function, spec::SVMSpec, X::AbstractMatrix, y⃗::AbstractVector;
                     nargs...)
  fit(spec, X, y⃗; callback = cb, nargs...)
end

# ==========================================================================
# Fallback for svm specific solver that are not implemented (for the given arguments)

function fit{T<:Real}(
    spec::SVMSpec,
    X::AbstractMatrix, y⃗::AbstractVector{T},
    ::SvmSolver,
    ::MvPredicton,
    ::SvmProblem;
    nargs...)
  throw(MvNotNativelyHandled())
end

# --------------------------------------------------------------------------

function fit{T<:Real}(
    spec::SVMSpec,
    X::AbstractMatrix, y⃗::AbstractVector{T},
    solver::SvmSolver,
    ::PredictionType,
    ::SvmProblem;
    nargs...)
  throw(ArgumentError("The types of the given arguments are not compatible with $(typeof(solver))"))
end


text:
using Random
using LinearAlgebra

"""
Generate an N×N matrix with uniformly random entries in [0,1].
"""
function uniform_matrix(N)
	return rand(N, N)
end

"""
Generate a random N×N matrix with entries satisfying the triangle inequality.
"""
function metric_matrix(N; dim=5, p=2)
	points = rand(N, dim)
	M = [norm(points[i,:] - points[j,:], p) for i in 1:N, j in 1:N]
	return M
end

"""
Generate a random 0-1 N×N matrix with edge probability p.
"""
function zeroone_matrix(N; p=0.25, selfedges=true, symmetric=false)
	M = rand(N, N) .≤ p
	if symmetric
		M = triu(M) .| triu(M)'
	end
	if !selfedges
		for i in 1:N
			M[i,i] = 0
		end
	end
	return M
end

"""
Generate a pair of N×N matricies with known optimal permutation using
Algorithm 4 from "Generating Quadratic Assignment Test Problems with
Known Optimal Permutations" by Li and Pardalos.
"""
function generate_qap(N)
	D = Matrix(I, N, N)
	O = .~D

	ΔA = rand(1:100)
	ΔB = rand(1:100)

	A = fill(ΔA, N, N)
	A[D] .= 0

	B = Array{Int,2}(undef, N, N)
	for i in 1:N, j in 1:N
		if i == j
			B[i,j] = 0
		else
			B[i,j] = rand(0:ΔB)
		end
	end

	R = findall(O)
	sort!(R, by=(ij -> B[ij]))

	X = rand(1:ΔA, N*(N-1))
	sort!(X, rev=true)

	for (k, ij) in enumerate(R)
		A[ij] = A[ij] - X[k]
	end

	σ = randperm(N)
	σinv = invperm(σ)
	B = B[σinv,σinv]

	return A, B, σ
end


text:
using Woody
using Base.Test

# write your own tests here
function test_workers()
    ctr = Controller(5543)
    r = Reporter(5543, ctr.key)
    report(r, "hello")
    report(r, "bonjour")
    report(r, "goodbye")
    timetable = finalize(ctr)
end

function test_collector()
    c = Collector(5543)
    run_collector_loop(c)
end

@test 1 == 1


text:
################################################################################
# Penalty
################################################################################

mutable struct Penalty{T}
	ρ::SVector{1,T}
	ρ_trial::SVector{1,T}
end

function Penalty(ρ::T) where {T}
	return Penalty(SVector{1,T}(ρ), SVector{1,T}(ρ))
end

################################################################################
# GameProblem
################################################################################

mutable struct GameProblem{KN,n,m,T,SVd,SVx}
    probsize::ProblemSize
	model::AbstractGameModel
    core::NewtonCore
	x0::SVx
    game_obj::GameObjective
    game_con::GameConstraintValues
	# colmult::CollisionMultiplier
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}
	pdtraj_trial::PrimalDualTraj{KN,n,m,T,SVd}
    Δpdtraj::PrimalDualTraj{KN,n,m,T,SVd}
	pen::Penalty{T}
    opts::Options
	stats::Statistics
end

function GameProblem(N::Int, dt::T, x0::SVx, model::AbstractGameModel, opts::Options,
	game_obj::GameObjective, game_con::GameConstraintValues
	) where {T,SVx}

	probsize = ProblemSize(N,model)
	pdtraj = PrimalDualTraj(probsize, dt)
	pdtraj_trial = PrimalDualTraj(probsize, dt)
	Δpdtraj = PrimalDualTraj(probsize, dt)
	core = NewtonCore(probsize)
	# colmult = CollisionMultiplier(probsize)
	stats = Statistics()
	pen = Penalty(opts.ρ_0)
	# Set the constraint parameters according to the options selected for the problem.
	set_constraint_params!(game_con, opts)

	TYPE = (eltype(pdtraj.pr), model.n, model.m, T, eltype(pdtraj.du), typeof(x0))
	return GameProblem{TYPE...}(probsize, model, core, x0, game_obj, game_con, #colmult,
		pdtraj, pdtraj_trial, Δpdtraj, pen, opts, stats)
end


text:
# ============================================================================ #
# CORE TYPES
# ============================================================================ #
"""
    AbstractNode{L<:Integer, T}

Can be a `SXTree` or a `SXTNode`. `L` is the type of label and `T` is the type
of data.
"""
abstract type AbstractNode{L<:Integer, T} end

"""
    SXTNode{L, T} <: AbstractNode{L, T}

A node in a simplex tree. The `data::T` field is nullable and empty by default.

Fields in constructor:

* `label::L`
* `depth::L`
* `parent::AbstractNode{L, T}`
* `data::Nullable{T}`
"""
struct SXTNode{L, T} <: AbstractNode{L, T}
    label::L
    depth::L

    data::Nullable{T}

    parent::AbstractNode{L, T}
    children::Dict{L, SXTNode{L, T}}

    next::Base.RefValue{SXTNode{L, T}}
    prev::Base.RefValue{SXTNode{L, T}}
end

function SXTNode(l::L, d::L,
              p::AbstractNode{L, T}, data=Nullable{T}()) where {L<:Integer, T}

    v = SXTNode(l, d, data, p, Dict{L, SXTNode{L, T}}(),
             Base.RefValue{SXTNode{L, T}}(), Base.RefValue{SXTNode{L, T}}())
    v.next[] = v
    v
end

Base.show(io::IO, nd::SXTNode) =
    print(io, "SXTNode{label=$(label(nd)), depth=$(depth(nd))",
          !isnull(nd.data) ? ", data=$(get(nd.data))" : "", "}")

"""
SXTree{L, T} <: AbstractNode{L, T}

The simplex tree.

TODO
"""
struct SXTree{L, T} <: AbstractNode{L, T}
    children::Vector{SXTNode{L, T}}
    cousin_list::Vector{Dict{L, SXTNode{L, T}}}
end

function SXTree(n_nodes::L, ::Type{T} = Float64) where {L<:Integer, T}
    sxtree = SXTree(Vector{SXTNode{L, T}}(n_nodes), Dict{L, SXTNode{L, T}}[])

    for l in 1:n_nodes
        sxtree[l] = SXTNode(l, 0, sxtree)
    end

    sxtree
end

Base.show(io::IO, sxt::SXTree{L, T}) where {L, T} =
    print(io, "SXTree{$L, $T} with $(length(sxt.children)) vertices.")

# ============================================================================ #
# ITERATORS
# ============================================================================ #
mutable struct CousinIter{L, T}
    start::SXTNode{L, T}
    first::Bool
end

Base.start(ci::CousinIter) = ci.start
Base.done(ci::CousinIter, nd::SXTNode) = !ci.first && nd === ci.start
Base.next(ci::CousinIter, nd::SXTNode) =
    if ci.first
        ci.first = false
        nd, nd.next[]
    else
        nd, nd.next[]
    end

Base.iteratorsize(::Type{CousinIter}) = Base.SizeUnknown()
Base.iteratoreltype(::Type{CousinIter}) = Base.HasEltype()
Base.eltype(::Type{CousinIter{L, T}}) where {L, T} = SXTNode{L, T}

Base.show(io::IO, ci::CousinIter) =
    show(io, "CousinIter($(ci.start))")

"""
Iterate over all nodes with the same label on the same level.
"""
cousins(nd::SXTNode) = CousinIter(nd, true)

# ---------------------------------------------------------------------------- #
struct UpIter{L, T}
    start::SXTNode{L, T}
end
Base.start(ui::UpIter) = ui.start
Base.done(ui::UpIter, nd::SXTNode) = false
Base.done(ui::UpIter, sxt::SXTree) = true
Base.next(ui::UpIter, nd::SXTNode) = nd, nd.parent

Base.iteratorsize(::Type{UpIter}) = Base.HasLength()
Base.iteratoreltype(::Type{UpIter}) = Base.HasEltype()
Base.eltype(::Type{UpIter{L, T}}) where {L, T} = SXTNode{L, T}
Base.length(ui::UpIter) = ui.start.depth + 1

Base.show(io::IO, ci::UpIter) =
    show(io, "UpIter($(ci.start))")

"""
Iterate ancestors upwards until reaching the root node.
"""
ancestors(nd::SXTNode) = UpIter(nd)

# ---------------------------------------------------------------------------- #
struct BFIter{L, T}
    start::AbstractNode{L, T}
end

function Base.start(bi::BFIter{L, T}) where {L, T}
    q = Queue(SXTNode{L, T})
    foreach(n -> enqueue!(q, n), children(bi.start))
    q
end
Base.done(::BFIter, q) = isempty(q)
function Base.next(::BFIter, q)
    nd = dequeue!(q)
    foreach(n -> enqueue!(q, n), children(nd))
    nd, q
end
Base.show(bi::BFIter) = println("BFIter($(bi.start))")

Base.iteratorsize(::Type{BFIter}) = Base.SizeUnknown()
Base.iteratoreltype(::Type{BFIter}) = Base.HasEltype()
Base.eltype(::Type{BFIter{L, T}}) where {L, T} = SXTNode{L, T}

struct DFIter{L, T}
    start::AbstractNode{L, T}
end

function Base.start(di::DFIter{L, T}) where {L, T}
    s = Stack(SXTNode{L, T})
    foreach(n -> push!(s, n), (reverse ∘ collect ∘ children)(di.start))
    s
end
Base.done(::DFIter, s) = isempty(s)
function Base.next(::DFIter, s)
    nd = pop!(s)
    foreach(n -> push!(s, n), (reverse ∘ collect ∘ children)(nd))
    nd, s
end
Base.show(bi::DFIter) = println("DFIter($(bi.start))")

Base.iteratorsize(::Type{DFIter}) = Base.SizeUnknown()
Base.iteratoreltype(::Type{DFIter}) = Base.HasEltype()
Base.eltype(::Type{DFIter{L, T}}) where {L, T} = SXTNode{L, T}

"""
Traverse successors of node. Us BFIter for breadth first search (default) and
DFIter for depth first search.
"""
successors(an::AbstractNode, ::Type{BFIter}) = BFIter(an)
successors(an::AbstractNode) = BFIter(an)
successors(an::AbstractNode, ::Type{DFIter}) = DFIter(an)

"""
Iterate over children of a node.
"""
children(nd::SXTNode) = values(nd.children)
children(sxt::SXTree) = sxt.children

# ============================================================================ #
# GETTERS, SETTERS, UTILS
# ============================================================================ #
Base.getindex(an::AbstractNode, i) = an.children[i]
function Base.getindex(an::AbstractNode, is...)
    nd = an.children[is[1]]
    for i in is[2:end]
        nd = nd.children[i]
    end
    nd
end
Base.setindex!(nd::AbstractNode, m, i) = nd.children[i] = m

label(::SXTree{L, T}) where {L, T} = zero(L)
label(nd::SXTNode) = nd.label

data(nd::SXTNode) = get(nd.data)
n_vertices(sxt::SXTree) = length(sxt.children)

Base.vec(nd::SXTNode) = reverse!(collect(label(n) for n in ancestors(nd)))

dim(sxt::SXTree) = length(sxt.cousin_list)

function n_simplices(sxt::SXTree)
    n = 0
    for _ in successors(sxt)
        n += 1
    end
    n
end

# ============================================================================ #
# PRIVATE HELPERS
# ============================================================================ #
parent(nd::SXTNode) = nd.parent

has_child(sxt::SXTree, l) = l <= length(sxt.children)
has_child(nd::SXTNode, l) = haskey(nd.children, l)
has_children(an::AbstractNode, l) = all(has_child.(an, l))
num_children(an::AbstractNode) = length(an.children)
has_children(an::AbstractNode) = num_children(an) > 0

depth(nd::SXTNode) = nd.depth

function connect_cousins!(n1::SXTNode, n2::SXTNode)
    n1_next = n1.next[]
    n1.next[] = n2
    n1_next.prev[] = n2
    n2.prev[] = n1
    n2.next[] = n1_next
end

function add_child!(sxt::SXTree, nd::SXTNode, l)
    d = depth(nd) + 1
    new = SXTNode(l, d, nd)
    nd[l] = new

    # Connect to cousins or add a new entry in cousin list.
    if haskey(sxt.cousin_list[d], l)
        connect_cousins!(sxt.cousin_list[d][l], new)
    else
        sxt.cousin_list[d][l] = new
    end

    new
end

# ============================================================================ #
# OPERATIONS
# ============================================================================ #
"""
    insert_simplex!(sxtree::SXTree, simplex)

Insert a simplex into sxtree.
The simplex is given as an array of labels<:Integer.
"""
function insert_simplex!(sxtree::SXTree{L, T},
                         simplex::AbstractArray{L, 1}) where {L, T}

    simplex = sort(simplex)
    check_simplex(sxtree, simplex)

    while dim(sxtree) < length(simplex)-1
        push!(sxtree.cousin_list, Dict())
    end

    q = Queue(Tuple{SXTNode{L, T}, typeof(simplex)})
    for i in 2:length(simplex)
        enqueue!(q, (sxtree[simplex[i-1]], simplex[i:end]))
    end

    while !isempty(q)
        nd, sx = dequeue!(q)
        for i in 1:length(sx)
            if !has_child(nd, sx[i])
                add_child!(sxtree, nd, sx[i])
            end
            enqueue!(q, (nd[sx[i]], sx[i+1:end]))
        end
    end
end

"""
    is_face(sx1, sx2)

Check whether sx1 is a face of sx2.
"""
function is_face(sx::AbstractArray, nd::SXTNode)
    sx = sort(sx)
    check_simplex(sx)

    i = endof(sx)
    for p in ancestors(nd)
        i == 0 && break
        i > p.depth + 1 && break

        if label(p) == sx[i]
            i -= 1
        end
    end

    i == 0
end

is_face(nd1::SXTNode, nd2::SXTNode) = is_face(vec(nd1), nd2)
is_face(sx1::AbstractArray, sx2::AbstractArray) = issubset(sx1, sx2)
is_face(nd::SXTNode, sx::AbstractArray) = issubset(vec(nd), sx)

"""
    cofaces(sxtree, sx)

Find all cofaces of simplex `sx`, including itself.
"""
function cofaces(sxtree::SXTree{L, T}, sx::AbstractArray{L}) where {L, T}
    sx = sort(sx)
    check_simplex(sxtree, sx)

    d = length(sx) - 1
    cofaces = SXTNode{L, T}[]

    dopt = d
    for list in sxtree.cousin_list[d:end]
        if haskey(list, sx[end])
            for nd in cousins(list[sx[end]])
                if is_face(sx, nd)
                    push!(cofaces, nd)
                    foreach(n -> push!(cofaces, n), successors(nd))
                end
            end
        end
    end

    cofaces
end

"""
Remove node WITHOUT fixing cofaces or SXTree dimension.
Do not use this, use rem_simplex! instead.
"""
function rem_node!(sxtree::SXTree{L, T}, nd::SXTNode{L, T}) where {L, T}
    depth(nd) == 0 && throw(ArgumentError("Can't remove top nodes!"))

    cousins = sxtree.cousin_list[depth(nd)]
    if nd.next[] === nd
        # Only cousin.
        delete!(sxtree.cousin_list[depth(nd)], label(nd))
    else
        # Remove from list.
        cousins[label(nd)] === nd && (cousins[label(nd)] = nd.next[])
        nd.next[].prev[] = nd.prev[]
        nd.prev[].next[] = nd.next[]
    end

    delete!(parent(nd).children, label(nd))
end

"""
    rem_simplex!(sxtree, sx)

Remove simplex `sx` from `sxtree`. `sx` can be given as a node or as an array.
"""
function rem_simplex!(sxtree::SXTree{L, T}, sx::AbstractArray{L}) where {L, T}
    sort(sx)
    check_simplex(sxtree, sx)

    foreach(n -> rem_node!(sxtree, n), cofaces(sxtree, sx))

    d = dim(sxtree)
    while isempty(sxtree.cousin_list[d])
        pop!(sxtree.cousin_list)
        d -= 1
    end
end

function rem_simplex!(sxtree::SXTree{L, T}, nd::SXTNode{L, T}) where {L, T}
    rem_simplex!(sxtree, vec(nd))
end


text:
#=
Julia implementation of the SPIKE distance, from the description of the Scholarpedia
http://www.scholarpedia.org/article/SPIKE-distance and the paper
Kreuz et al. (2013), [Monitoring spike train synchrony](https://doi.org/10.1152/jn.00873.2012)

Implementation by Geoge Datseris.
=#
export SPIKE_distance_profile, SPIKE_distance

"""
    SPIKE_distance_profile(y1, y2) -> t, S(t)
Calculate the SPIKE distance profile function ``S(t)``, Eq. (19) of [1], given input
spike trains `y1, y2`. Return `t, S(t)`.

The keywords `t0 = min(y1[1], y2[1]) - 1` and
`tf = max(y1[end], y2[end]) + 1` add auxilary spikes to both trains
at these times, as demanded by the algorithm. (Auxilary spikes are only added
in case `y1, y2` do not start or end with the given `t0, tf`)

[1] : Kreuz et al. (2013), [Monitoring spike train synchrony](https://doi.org/10.1152/jn.00873.2012)
"""
function SPIKE_distance_profile(y1, y2;
        t0 = min(y1[1], y2[1]) - 1,
        tf = max(y1[end], y2[end]) + 1,
    )
    @assert issorted(y1) && issorted(y2)
    @assert y1 == unique(y1) && y2 == unique(y2)
    y1 = copy(y1); y2 = copy(y2)
    # Add auxilary spikes:
    y1[1] != t0 && pushfirst!(y1, t0)
    y2[1] != t0 && pushfirst!(y2, t0)
    y1[end] != tf && push!(y1, tf)
    y2[end] != tf && push!(y2, tf)
    tvec, indices = _initialize_SPIKE_data(y1, y2)
    # S are the values of S(t), for t ∈ tvec.
    # Its values are only stored exactly before and exactly on a spike
    # (as only these values are need for the full S(t)) function)
    S⁻ = zeros(Float64, length(tvec))
    S⁺ = zeros(Float64, length(tvec))
    i1, i2 = 1, 1 # indices of the previous spike for each train
    @inbounds for k in 2:length(tvec)-1
        t = tvec[k]
        # First compute S given at the time of the next spike, but before
        # being exactly ON the spike
        S⁻[k] = _compute_S(t, y1, y2, i1, i2)
        # Then find out which indices to increment (where is current spike from)
        if indices[k] == 3 # both trains share current spike
            i1 += 1; i2 += 1
        elseif indices[k] == 1 # current spike from train 1
            i1 += 1
        else # current spike from train 2
            i2 += 1
        end
        # Then compute S exactly ON the current spike
        S⁺[k] = _compute_S(t, y1, y2, i1, i2)
    end
    tvec = sort!(vcat(tvec, tvec)) # duplicate time vector
    S = zeros(length(tvec))
    S[1:2:end] .= S⁻
    S[2:2:end] .= S⁺
    return tvec, S
end

"""
Core computatation for  ``S(t)`` at time ``t``. Whether we are exactly
before or exactly on top of a spike is taken care of by `i1, i2`.
"""
function _compute_S(t, y1, y2, i1, i2)
    # Because of the way `_initialize_SPIKE_data` works,
    # this is guaranteed to be correct from the data preperation:
    tP1 = y1[i1]; tP2 = y2[i2]
    tF1 = y1[i1+1]; tF2 = y2[i2+1]
    xISI1 = tF1 - tP1;  xISI2 = tF2 - tP2
    xP1 = t - tP1; xP2 = t - tP2
    xF1 = tF1 - t; xF2 = tF2 - t
    ΔtP1 = minimum_Δt(tP1, y2, i2); ΔtF1 = minimum_Δt(tF1, y2, i2);
    ΔtP2 = minimum_Δt(tP2, y1, i1); ΔtF2 = minimum_Δt(tF2, y1, i1)
    S1 = (ΔtP1*xF1 + ΔtF1*xP1)/xISI1; S2 = (ΔtP2*xF2 + ΔtF2*xP2)/xISI2
    S = 2(S1*xISI2 + S2*xISI1)/(xISI1 + xISI2)^2
    return S
end

"""
Create smart containers for the time vector (that has ALL spikes sorted) as
well as the indices of where each spike comes from. This allows for exceptionally
simplified source code for the SPIKE distance profile ``S(t)``.
"""
function _initialize_SPIKE_data(y1::AbstractVector{T1}, y2::AbstractVector{T2}) where {T1, T2}
    T = promote_type(T1, T2)
    tvec = Vector{T}()
    indices = Vector{Int8}()
    i1 = i2 = 1 # indices of the spike trains
    while i1 ≤ length(y1) && i2 ≤ length(y2)
        v1, v2 = y1[i1], y2[i2]
        if v1 == v2
            push!(tvec, v1); push!(indices, 3)
            i1 += 1; i2 += 1;
        elseif v1 < v2
            push!(tvec, v1); push!(indices, 1)
            i1 += 1
        else
            push!(tvec, v2); push!(indices, 2)
            i2 += 1
        end
    end
    @assert tvec == sort!(unique(vcat(y1, y2)))
    return tvec, indices
end

"""
    minimum_Δt(t, y, i = 1) -> Δt
Return the minimum distance `Δt` between time point `t` and the spikes in `y`.
Start searching from index `i` of `y` and efficiently scan indices starting from
`i` and incrementally increasing and decreasing.
"""
function minimum_Δt(t::T, y, i = 1) where {T}
    dleft, dright = typemax(t), typemax(t)
    for j in i:-1:1 # scan from i until 1
        z = abs(t - y[j])
        z > dleft ? break : (dleft = z)
    end
    for j in i+1:1:length(y) # scan from i+1 until the end
        z = abs(t - y[j])
        z > dright ? break : (dright = z)
    end
    return min(dright, dleft)
end

"""
    SPIKE_distance(y1, y2) -> D_S
Calculate the SPIKE distance ``D_S``, which is the average of
[`SPIKE_distance_profile`](@ref).
"""
function SPIKE_distance(y1, y2; kwargs...)
    t, S = SPIKE_distance_profile(y1, y2; kwargs...)
    return trapezoid_integral(t, S)/(t[end] - t[1])
end

function trapezoid_integral(t, S::AbstractVector{T}) where {T}
    I = zero(T)
    for i in 2:length(t)
        @inbounds I += (S[i] + S[i-1])*(t[i] - t[i-1])
    end
    return T(0.5) * I
end


text:
# define atom and module

struct Atom
    name::String
    Z::Int
    Ne::Int
    orbitals::Vector{STONG}
    R::Vec3{Float64}
    Atom(name::AbstractString, Z::Integer, Ne::Integer, orbitals::Vector{STONG{N}}, R::Vec3) where N = new(
        name, Z%Int, Ne%Int, orbitals, float(R)
    )
end

Atom(name::AbstractString, R::Vec3) = begin
    name = titlecase(name)
    number = atomic_number[name]
    Rf = float(R)
    number <= 2 && return Atom(name, number, number, [
        STONG(coefficient_data[name]["1s"], Rf)
    ], Rf)
    number <= 10 && return Atom(name, number, number, [
        STONG(coefficient_data[name]["1s"], Rf),
        STONG(coefficient_data[name]["2s"], Rf),
        STONG(coefficient_data[name]["2p"], Rf; type=:px),
        STONG(coefficient_data[name]["2p"], Rf; type=:py),
        STONG(coefficient_data[name]["2p"], Rf; type=:pz)
    ], Rf)
    error("atomic_number > 10 atoms are not support yet")
end

const atomic_number = Dict{String, Int}(
    "H"  => 1,
    "He" => 2,
    "Li" => 3,
    "Be" => 4,
    "B"  => 5,
    "C"  => 6,
    "N"  => 7,
    "O"  => 8,
    "F"  => 9,
    "Ne" => 10,
    "Na" => 11,
    "Mg" => 12,
    "Al" => 13,
    "Si" => 14,
    "P"  => 15,
    "S"  => 16,
    "Cl" => 17,
    "Ar" => 18,
    "K"  => 19,
    "Ca" => 20
)

hydrogen_atom(R::Vec3) = Atom("H", R)
helium_atom(R::Vec3)   = Atom("He", R)
carbon_atom(R::Vec3)   = Atom("C", R)
oxygen_atom(R::Vec3)   = Atom("O", R)
fluorine_atom(R::Vec3) = Atom("F", R)

struct Molecule
    name::String
    charge::Int
    atoms::Vector{Atom}
    Molecule(name::AbstractString, charge::Integer, atoms::Vector{Atom}) = new(name, charge%Int, atoms)
end

Molecule(name::AbstractString, atoms::Vector{Atom}) = Molecule(name, 0, atoms)
Molecule(atoms::Vector{Atom}) = Molecule("", 0, atoms)


text:
_hbox(args...) = Widgets.div(args...; style = Dict("display" => "flex", "flex-direction"=>"row"))
_vbox(args...) = Widgets.div(args...; style = Dict("display" => "flex", "flex-direction"=>"column"))

get_kwargs(; kwargs...) = kwargs
string2kwargs(s::AbstractString) = eval(Meta.parse("get_kwargs($s)"))

const analysis_options = OrderedDict(
    "" => nothing,
    "Cumulative" => Recombinase.cumulative,
    "Density" => Recombinase.density,
    "Hazard" => Recombinase.hazard,
    "Prediction" => Recombinase.prediction,
)

function is_small(col, n = 100)
    for (ind, _) in enumerate(IterTools.distinct(col))
        ind > n && return false
    end
    return true
end

function take_small(t, vec::AbstractVector{Symbol}, n = 100)
    filter(vec) do sym
        col = get(columns(t), sym, ())
        return is_small(col, n)
    end
end

"""
`gui(data, plotters; postprocess = NamedTuple())`

Create a gui around `data::IndexedTable` given a list of plotting
functions plotters.

## Examples

```julia
using StatsPlots, Recombinase, JuliaDB, Interact
school = loadtable(joinpath(Recombinase.datafolder, "school.csv"))
plotters = [plot, scatter, groupedbar]
Recombinase.gui(school, plotters)
```
"""
function gui(data′, plotters; postprocess = NamedTuple())
    (data′ isa AbstractObservable) || (data′ = Observable{Any}(data′))
    data = @map table(&data′, copy = false, presorted = true)
    ns = @map collect(colnames(&data))
    maybens = @map vcat(Symbol(), &ns)
    xaxis = dropdown(ns,label = "X")
    yaxis = dropdown(maybens,label = "Y")
    an_opt = dropdown(analysis_options, label = "Analysis")
    axis_type = dropdown([:automatic, :continuous, :discrete, :vectorial], label = "Axis type")
    error = dropdown(@map(vcat(automatic, &ns)), label="Error")
    styles = collect(keys(style_dict))
    sort!(styles)
    smallns = map(take_small, data, maybens)
    splitters = [dropdown(smallns, label = string(style)) for style in styles]
    plotter = dropdown(plotters, label = "Plotter")
    ribbon = toggle("Ribbon", value = false)
    btn = button("Plot")
    output = Observable{Any}("Set the dropdown menus and press plot to get started.")
    plot_kwargs = Widgets.textbox("Insert optional plot attributes")
    @map! output begin
        &btn
        select = yaxis[] == Symbol() ? xaxis[] : (xaxis[], yaxis[])
        grps = Dict(key => val[] for (key, val) in zip(styles, splitters) if val[] != Symbol())
        an = an_opt[]
        an_inf = isnothing(an) ? nothing : Analysis{axis_type[]}(an)
        args, kwargs = series2D(
                                an_inf,
                                &data,
                                Group(; grps...);
                                postprocess = postprocess,
                                select = select,
                                error = error[],
                                ribbon = ribbon[]
                               )
        plotter[](args...; kwargs..., string2kwargs(plot_kwargs[])...)
    end
    ui = Widget(
        OrderedDict(
            :xaxis => xaxis,
            :yaxis => yaxis,
            :analysis => an_opt,
            :axis_type => axis_type,
            :error => error,
            :plotter => plotter,
            :plot_button => btn,
            :plot_kwargs => plot_kwargs,
            :ribbon => ribbon,
            :splitters => splitters,
        ),
        output = output
    )
    Widgets.@layout! ui Widgets.div(
                                    _hbox(
                                          :xaxis,
                                          :yaxis,
                                          :analysis,
                                          :axis_type,
                                          :error,
                                          :plotter
                                         ),
                                    :ribbon,
                                    :plot_button,
                                    _hbox(
                                          _vbox(:splitters...),
                                          _vbox(output, :plot_kwargs)
                                         )
                                   )
end


text:
using Glob, DelimitedFiles, PyPlot#, Plots

files = glob("boundary[0-9]*.txt", "data")

xy = [zeros(Float32, 0,2) for _ in 1:length(files)]

for i in 1:length(files)
   bc = readdlm(files[i], ',', Float32, header=true)
   xy[i] = bc[1][:,1:2]
end

figure()

for i ∈ 1:length(files)
   plot(xy[i][:,1],xy[i][:,2], label="boundary")
   xlim(-2.1,4.2)
   ylim(-3.0,3.0)
   title("t = $i")
   savefig("figure/PIC_test"*lpad(i,4,"0")*".png")
   clf()
end

#=
anim = @animate for i ∈ 1:length(files)
   plot(xy[i][:,1],xy[i][:,2], xlims=(-2.1,4.2), ylims=(-3.0,3.0), label="boundary")
   title!("t = $i")
end

gif(anim, "test.mp4", fps = 15)
=#


text:
"""
A base abstract type for `OptimizationProblem` search space.
A concrete `SearchSpace` subtype specifies the valid
candidate points that could be considered in a search/optimization.
"""
abstract type SearchSpace end

"""
A base abstract type for search spaces with a fixed finite number of dimensions.
Applicable to the vast majority of problems.
"""
abstract type FixedDimensionSearchSpace <: SearchSpace end

"""
A `SearchSpace` with `N`-dimensional rectangle as a set of valid points.
I.e. `dimmin(ss, i)` ≤ `x[i]` ≤ `dimmax(ss, i)` for each dimension `i`.
"""
abstract type RectSearchSpace <: FixedDimensionSearchSpace end

Base.@deprecate_binding ContinuousSearchSpace RectSearchSpace

"""
The point of the `SearchSpace`.

The abstract type. It allows different actual implementations to be used,
e.g `Vector` or `SubArray`.
"""
const AbstractIndividual = AbstractVector{Float64}

"""
The point of the `SearchSpace`.

The concrete type that could be used for storage.
"""
const Individual = Vector{Float64}

"""
The valid range of values for a specific dimension of `SearchSpace`.
"""
const ParamBounds = Tuple{Float64,Float64}

numdims(ss::RectSearchSpace) = length(dimmin(ss))

"""
    dimmin(ss::SearchSpace, [i::Integer])

A minimal valid value for `i`-th dimension of `ss`, or a vector of
minimal valid values for each dimension of `ss` if no `i` was given.
"""
dimmin(ss::RectSearchSpace) = ss.dimmin

dimmin(ss::RectSearchSpace, i::Integer) = ss.dimmin[i]

"""
    dimmax(ss::SearchSpace, [i::Integer])

A maximal valid value for `i`-th dimension of `ss`, or a vector of
maximal valid values for each dimension of `ss` if no `i` was given.
"""
dimmax(ss::RectSearchSpace) = ss.dimmax

dimmax(ss::RectSearchSpace, i::Integer) = ss.dimmax[i]

"""
    dimdelta(ss::SearchSpace, [i::Integer])

A delta of maximal and minimal valid values for `i`-th dimension of `ss`
(*diameter* of `i`-th dimension), or a vector of deltas for each dimension
if no `i` given.
"""
dimdelta(ss::RectSearchSpace) = ss.dimdelta

dimdelta(ss::RectSearchSpace, i::Integer) = ss.dimdelta[i]

@deprecate mins(ss) dimmin(ss)
@deprecate maxs(ss) dimmax(ss)
@deprecate deltas(ss) dimdelta(ss)
@deprecate diameters(ss) dimdelta(ss)

"""
    dimrange(ss::SearchSpace, [i::Integer])

Gets a `ParamsRange` tuple of minimal and maximal valid values for
`i`-th dimension of `ss`, or a vector of `ParamsRange` tuples
for each dimension if no `i` given.
"""
dimrange(ss::RectSearchSpace, i::Integer) = (dimmin(ss, i), dimmax(ss, i))
@deprecate range_for_dim(ss, i) dimrange(ss, i)

dimrange(ss::RectSearchSpace) = tuple.(dimmin(ss), dimmax(ss))
@deprecate ranges(ss) dimrange(ss)

"""
    in(ind::AbstractIndividual, ss::SearchSpace)

Check if given individual lies in the given search space.
"""
function Base.in(ind::AbstractIndividual, ss::RectSearchSpace)
    # FIXME checks for the discrete space?
    @assert length(ind) == numdims(ss)
    @inbounds for i in eachindex(ind)
        (dimmin(ss, i) <= ind[i] <= dimmax(ss, i)) || return false
    end
    return true
end

"""
`SearchSpace` defined by a continuous range of valid values for each dimension.
"""
struct ContinuousRectSearchSpace <: RectSearchSpace
    dimmin::Vector{Float64}    # minimal valid value per dimension
    dimmax::Vector{Float64}    # maximal valid value per dimension
    dimdelta::Vector{Float64}  # delta/diameter (dimmax-dimmin) per dimension

    function ContinuousRectSearchSpace(
        dimmin::AbstractVector{<:Real},
        dimmax::AbstractVector{<:Real}
    )
        length(dimmin) == length(dimmax) ||
            throw(DimensionMismatch("dimmin and dimmax should have the same length"))
        all(xy -> xy[1] <= xy[2], zip(dimmin, dimmax)) ||
            throw(ArgumentError("dimmin should not exceed dimmax"))
        new(dimmin, dimmax, dimmax .- dimmin)
    end
end

ContinuousRectSearchSpace(ranges) =
    ContinuousRectSearchSpace(getindex.(ranges, 1), getindex.(ranges, 2))

Base.:(==)(a::ContinuousRectSearchSpace,
           b::ContinuousRectSearchSpace) =
    (numdims(a) == numdims(b)) &&
    (dimmin(a) == dimmin(b)) &&
    (dimmax(a) == dimmax(b))

"""
Projects a given point onto the search space coordinate-wise.
"""
feasible(v::AbstractIndividual, ss::RectSearchSpace) =
    clamp.(v, dimmin(ss), dimmax(ss))

# concatenates two range-based search spaces
Base.vcat(ss1::ContinuousRectSearchSpace,
          ss2::ContinuousRectSearchSpace) =
    ContinuousRectSearchSpace(vcat(dimmin(ss1), dimmin(ss2)),
                              vcat(dimmax(ss1), dimmax(ss2)))

# all dimensions are continuous
dimdigits(ss::ContinuousRectSearchSpace, i::Integer) = -1
dimdigits(ss::ContinuousRectSearchSpace) = fill(-1, numdims(ss))

"""
0-dimensional search space.
Could be used as a placeholder for optional `SearchSpace` parameters.
"""
const ZERO_SEARCH_SPACE = ContinuousRectSearchSpace(Vector{Float64}(), Vector{Float64}())

"""
`RectSearchSpace` that allows both continuous and discrete dimensions.
Discrete dimensions are defined by the number of digits (`dimdigits(ss)`) after decimal point (i.e. precision)

If `dimdigits(ss, i)` is negative, `i`-th dimension is considered continuous.
"""
struct MixedPrecisionRectSearchSpace <: RectSearchSpace
    dimmin::Vector{Float64}
    dimmax::Vector{Float64}
    dimdelta::Vector{Float64}
    dimdigits::Vector{Int}

    function MixedPrecisionRectSearchSpace(dimmin::AbstractVector,
                                           dimmax::AbstractVector,
                                           dimdigits::AbstractVector{<:Integer})
        length(dimmin) == length(dimmax) == length(dimdigits) ||
            throw(DimensionMismatch("dimmin, dimmax and dimdigits should have the same length"))
        dmin = Float64[dimdigits[i] >= 0 ? round(dimmin[i], digits=dimdigits[i]) : dimmin[i] for i in eachindex(dimdigits)]
        dmax = Float64[dimdigits[i] >= 0 ? round(dimmax[i], digits=dimdigits[i]) : dimmax[i] for i in eachindex(dimdigits)]
        all(xy -> xy[1] <= xy[2], zip(dmin, dmax)) ||
            throw(ArgumentError("discretized dimmin should not exceed dimmax"))
        new(dmin, dmax, dmax .- dmin,
            copyto!(Vector{Int}(undef, length(dimdigits)), dimdigits))
    end
end

MixedPrecisionRectSearchSpace(dimranges::AbstractVector{<:Tuple{<:Number, <:Number}},
                              dimdigits::AbstractVector{<:Integer}) =
    MixedPrecisionRectSearchSpace(getindex.(dimranges, 1), getindex.(dimranges, 2), dimdigits)

dimdigits(ss::MixedPrecisionRectSearchSpace, i::Integer) = ss.dimdigits[i]
dimdigits(ss::MixedPrecisionRectSearchSpace) = ss.dimdigits

feasible(x::Real, u::Real, v::Real, dimdigits::Integer) =
    clamp(dimdigits >= 0 ? round(x, digits=dimdigits) : x, u, v)

"""
Projects a given point onto the search space coordinate-wise.
"""
feasible(v::AbstractIndividual, ss::MixedPrecisionRectSearchSpace) =
    feasible.(v, dimmin(ss), dimmax(ss), dimdigits(ss))

# concatenates two rectangular search spaces
# one of the should be MixedPrecisionRectSearchSpace, otherwise
# ContinuousRectSearchSpace-only method would be called
Base.vcat(ss1::RectSearchSpace, ss2::RectSearchSpace) =
    MixedPrecisionRectSearchSpace(vcat(dimmin(ss1), dimmin(ss2)),
                                  vcat(dimmax(ss1), dimmax(ss2)),
                                  vcat(dimdigits(ss1), dimdigits(ss2)))

"""
    RectSearchSpace(dimranges::AbstractVector; [dimdigits = nothing])

Create `RectSearchSpace` with given range of valid values for each dimension and,
if specified, `dimdigits` precision for each dimension.
Returns `MixedPrecisionRectSearchSpace` if there's at least one dimension
with specified precision (dimdigits[i] ≥ 0), otherwise `ContinuousRectSearchSpace`.
"""
RectSearchSpace(dimranges::AbstractVector;
                dimdigits::Union{AbstractVector{<:Integer}, Nothing} = nothing) =
    dimdigits === nothing || all(ddigits -> ddigits < 0, dimdigits) ?
        ContinuousRectSearchSpace(dimranges) :
        MixedPrecisionRectSearchSpace(dimranges, dimdigits)

"""
Create `RectSearchSpace` with given number of dimensions
and given range of valid values for each dimension.
"""
RectSearchSpace(numdims::Integer, range=(0.0, 1.0); dimdigits::Union{Integer, Nothing} = nothing) =
    RectSearchSpace(fill(range, numdims);
                    dimdigits = dimdigits !== nothing && dimdigits >= 0 ? fill(dimdigits, numdims) : nothing)

@deprecate symmetric_search_space(numdims, range=(0.0, 1.0); dimdigits=nothing) RectSearchSpace(numdims, range, dimdigits=dimdigits)

# round x to fit ss requirements
# by default it does nothing
_round!(x::Union{AbstractVector, AbstractMatrix}, ss::RectSearchSpace) = x

function _round!(x::AbstractVector, ss::MixedPrecisionRectSearchSpace)
    length(x) == numdims(ss) ||
        throw(DimensionMismatch("Vector length should match search space dimensions"))
    ddigits = dimdigits(ss)
    @inbounds for j in eachindex(x)
        if ddigits[j] >= 0
            x[j] = round(x[j], digits=ddigits[j])
        end
    end
    return x
end

function _round!(x::AbstractMatrix, ss::MixedPrecisionRectSearchSpace)
    size(x, 1) == numdims(ss) ||
        throw(DimensionMismatch("Matrix rows should match search space dimensions"))
    ddigits = dimdigits(ss)
    @inbounds for i in axes(x, 2)
        for j in axes(x, 1)
            if ddigits[j] >= 0
                x[j, i] = round(x[j, i], digits=ddigits[j])
            end
        end
    end
    return x
end

"""
Generate one random candidate.
"""
rand_individual(ss::RectSearchSpace) =
    _round!(rand(numdims(ss)) .* dimdelta(ss) .+ dimmin(ss), ss)

"""
    rand_individuals(ss, n; [method=:latin_hypercube])

Generate `n` individuals by randomly sampling in the `ss` search space using
specified sampling `method`. The supported methods are:
 * `uniform`: uniform independent sampling for each dimension
 * `latin_hypercube` (the default): use latin hypercube sampling method
"""
function rand_individuals(ss::RectSearchSpace, n::Integer; method::Symbol=:latin_hypercube)
    if method == :uniform
        return _round!(rand(numdims(ss), n) .* dimdelta(ss) .+ dimmin(ss), ss)
    elseif method == :latin_hypercube
        return _round!(Utils.latin_hypercube_sampling(dimmin(ss), dimmax(ss), n), ss)
    else
        throw(ArgumentError("Unknown sampling method \"$method\""))
    end
end

@deprecate rand_individuals_lhs(ss::RectSearchSpace, n::Integer) rand_individuals(ss, n, method=:latin_hypercube)


text:
import Base.hash
import StatsBase.counts

hash{T<:Integer}(X::Array{T,1}) = X

function hash{T<:Integer}(X::Array{T,2})
	ndims = size(X,1)
	base = maximum(X)+1
	x = base.^[0:ndims-1]
	return X'*x
end

function counts{T<:Integer}(X::Array{T,2})
	y = hash(X)
	n = StatsBase.counts(y,0:maximum(y))
	return n
end

function counts{T<:Integer}(X::Array{T,1})
	n = StatsBase.counts(X,0:maximum(X))
	return n
end

function stratify(n::Array{Int64,1})
	k = StatsBase.counts(n,1:maximum(n))
	kx = k[k.>0]'
	nx = sort(unique(n))
	nx = nx[nx.>0]'
	return nx,kx
end