text:
"""
Create a geometry object of the appropriate type from it's well known
binary (WKB) representation.

### Parameters
* `data`: pointer to the input BLOB data.
* `spatialref`: handle to the spatial reference to be assigned to the created
    geometry object. This may be `NULL` (default).
"""
function unsafe_fromWKB(
        ::Type{G},
        data,
        spatialref::SpatialRef = SpatialRef(C_NULL)
    ) where G <: AbstractGeometry
    geom = Ref{GDALGeometry}()
    result = @gdal(OGR_G_CreateFromWkb::GDAL.OGRErr,
        data::Ptr{Cuchar},
        spatialref.ptr::GDALSpatialRef,
        geom::Ptr{GDALGeometry},
        sizeof(data)::Cint
    )
    @ogrerr result "Failed to create geometry from WKB"
    G(geom[])
end

unsafe_fromWKB(data, args...) = unsafe_fromWKB(Geometry, data, args...)

fromWKB(data, args...) = unsafe_fromWKB(IGeometry, data, args...)

"""
Create a geometry object of the appropriate type from its well known text
(WKT) representation.

### Parameters
* `data`: input zero terminated string containing WKT representation of the
    geometry to be created. The pointer is updated to point just beyond that
    last character consumed.
* `spatialref`: handle to the spatial reference to be assigned to the created
    geometry object. This may be `NULL` (default).
"""
function unsafe_fromWKT(
        ::Type{G},
        data::Vector{String},
        spatialref::SpatialRef = SpatialRef(C_NULL)
    ) where G <: AbstractGeometry
    geom = Ref{GDALGeometry}()
    result = @gdal(OGR_G_CreateFromWkt::GDAL.OGRErr,
        data::StringList,
        spatialref.ptr::GDALSpatialRef,
        geom::Ptr{GDALGeometry}
    )
    @ogrerr result "Failed to create geometry from WKT"
    G(geom[])
end

unsafe_fromWKT(data::Vector{String}, args...) = unsafe_fromWKT(Geometry, data, args...)

fromWKT(data::Vector{String}, args...) = unsafe_fromWKT(IGeometry, data, args...)

unsafe_fromWKT(data::String, args...) = unsafe_fromWKT([data], args...)

fromWKT(data::String, args...) = fromWKT([data], args...)

"""
Destroy geometry object.

Equivalent to invoking delete on a geometry, but it guaranteed to take place
within the context of the GDAL/OGR heap.
"""
function destroy(geom::AbstractGeometry)
    GDAL.destroygeometry(geom.ptr)
    geom.ptr = C_NULL
end

"Returns a copy of the geometry with the original spatial reference system."
function unsafe_clone(::Type{G}, geom::AbstractGeometry) where G <: AbstractGeometry
    G(GDAL.clone(geom.ptr))
end

unsafe_clone(geom::G) where {G <: AbstractGeometry} = unsafe_clone(G, geom)

clone(geom::AbstractGeometry) = unsafe_clone(IGeometry, geom)

"""
Create an empty geometry of desired type.

This is equivalent to allocating the desired geometry with new, but the
allocation is guaranteed to take place in the context of the GDAL/OGR heap.
"""
function unsafe_creategeom(::Type{G}, geomtype::OGRwkbGeometryType ) where G <: AbstractGeometry
    G(GDAL.checknull(@gdal(OGR_G_CreateGeometry::GDALGeometry, geomtype::GDAL.OGRwkbGeometryType)))
end

unsafe_creategeom(geomtype::OGRwkbGeometryType) = unsafe_creategeom(Geometry, geomtype)

creategeom(geomtype::OGRwkbGeometryType) = unsafe_creategeom(IGeometry, geomtype)

"""
Tries to force the provided geometry to the specified geometry type.

### Parameters
* `geom`: the input geometry.
* `targettype`: target output geometry type.
# `options`: (optional) options as a null-terminated vector of strings

It can promote 'single' geometry type to their corresponding collection type
(see OGR_GT_GetCollection()) or the reverse. non-linear geometry type to their
corresponding linear geometry type (see OGR_GT_GetLinear()), by possibly
approximating circular arcs they may contain. Regarding conversion from linear
geometry types to curve geometry types, only "wraping" will be done. No attempt
to retrieve potential circular arcs by de-approximating stroking will be done.
For that, OGRGeometry::getCurveGeometry() can be used.

The passed in geometry is cloned and a new one returned.
"""
function unsafe_forceto(
        ::Type{G},
        geom::AbstractGeometry,
        targettype::OGRwkbGeometryType,
        options = StringList(C_NULL)
    ) where G <: AbstractGeometry
    G(GDAL.checknull(@gdal(OGR_G_ForceTo::GDALGeometry,
        unsafe_clone(geom).ptr::GDALGeometry,
        targettype::GDAL.OGRwkbGeometryType,
        options::StringList
    )))
end
unsafe_forceto(geom::AbstractGeometry, args...) =
    unsafe_forceto(Geometry, geom, args...)
forceto(geom::AbstractGeometry, args...) =
    unsafe_forceto(IGeometry, geom, args...)

"""
Get the dimension of the geometry. 0 for points, 1 for lines and 2 for surfaces.

This function corresponds to the SFCOM IGeometry::GetDimension() method. It
indicates the dimension of the geometry, but does not indicate the dimension of
the underlying space (as indicated by OGR_G_GetCoordinateDimension() function).
"""
getdim(geom::AbstractGeometry) = GDAL.getdimension(geom.ptr)

"""
Get the dimension of the coordinates in this geometry.

### Returns
In practice this will return 2 or 3. It can also return 0 in the case of an
empty point.
"""
getcoorddim(geom::AbstractGeometry) = GDAL.getcoordinatedimension(geom.ptr)

"""
Set the coordinate dimension.

This method sets the explicit coordinate dimension. Setting the coordinate
dimension of a geometry to 2 should zero out any existing Z values. Setting the
dimension of a geometry collection, a compound curve, a polygon, etc. will
affect the children geometries. This will also remove the M dimension if present
before this call.
"""
function setcoorddim!(geom::AbstractGeometry, dim::Integer)
    GDAL.setcoordinatedimension(geom.ptr ,dim)
end

"Computes and returns the bounding envelope for this geometry"
function getenvelope(geom::AbstractGeometry)
    envelope = Ref{GDAL.OGREnvelope}(GDAL.OGREnvelope(0, 0, 0, 0))
    GDAL.getenvelope(geom.ptr, envelope)
    envelope[]
end

"Computes and returns the bounding envelope (3D) for this geometry"
function getenvelope3d(geom::AbstractGeometry)
    envelope = Ref{GDAL.OGREnvelope3D}(GDAL.OGREnvelope3D(0, 0, 0, 0, 0, 0))
    GDAL.getenvelope3d(geom.ptr, envelope)
    envelope[]
end

"""
Convert a geometry well known binary format.

### Parameters
* `geom`: handle on the geometry to convert to a well know binary data from.
* `order`: One of wkbXDR or [wkbNDR] indicating MSB or LSB byte order resp.
"""
function toWKB(geom::AbstractGeometry, order::OGRwkbByteOrder=GDAL.wkbNDR)
    buffer = Array{Cuchar}(wkbsize(geom))
    result = @gdal(OGR_G_ExportToWkb::GDAL.OGRErr,
        geom.ptr::GDALGeometry,
        order::GDAL.OGRwkbByteOrder,
        buffer::Ptr{Cuchar}
    )
    @ogrerr result "Failed to export geometry to WKB"
    buffer
end

"""
Convert a geometry into SFSQL 1.2 / ISO SQL/MM Part 3 well known binary format.

### Parameters
* `geom`: handle on the geometry to convert to a well know binary data from.
* `order`: One of wkbXDR or [wkbNDR] indicating MSB or LSB byte order resp.
"""
function toISOWKB(geom::AbstractGeometry, order::OGRwkbByteOrder=GDAL.wkbNDR)
    buffer = Array{Cuchar}(wkbsize(geom))
    result = @gdal(OGR_G_ExportToIsoWkb::GDAL.OGRErr,
        geom.ptr::GDALGeometry,
        order::GDAL.OGRwkbByteOrder,
        buffer::Ptr{Cuchar}
    )
    @ogrerr result "Failed to export geometry to ISO WKB"
    buffer
end

"Returns size (in bytes) of related binary representation."
wkbsize(geom::AbstractGeometry) = GDAL.wkbsize(geom.ptr)

"Convert a geometry into well known text format."
function toWKT(geom::AbstractGeometry)
    wkt_ptr = Ref{Ptr{UInt8}}()
    result = GDAL.exporttowkt(geom.ptr, wkt_ptr)
    @ogrerr result "OGRErr $result: failed to export geometry to WKT"
    wkt = unsafe_string(wkt_ptr[])
    GDAL.C.OGRFree(Ptr{UInt8}(wkt_ptr[]))
    wkt
end

"Convert a geometry into SFSQL 1.2 / ISO SQL/MM Part 3 well known text format."
function toISOWKT(geom::AbstractGeometry)
    isowkt_ptr = Ref{Ptr{UInt8}}()
    result = GDAL.exporttoisowkt(geom.ptr, isowkt_ptr)
    @ogrerr result "OGRErr $result: failed to export geometry to ISOWKT"
    wkt = unsafe_string(isowkt_ptr[])
    GDAL.C.OGRFree(Ptr{UInt8}(isowkt_ptr[]))
    wkt
end

"Fetch geometry type code"
getgeomtype(geom::AbstractGeometry) =
    OGRwkbGeometryType(GDAL.getgeometrytype(geom.ptr))

"Fetch WKT name for geometry type."
getgeomname(geom::AbstractGeometry) = GDAL.getgeometryname(geom.ptr)

"Convert geometry to strictly 2D."
function flattento2d!(geom::G) where G <: AbstractGeometry
    GDAL.flattento2d(geom.ptr)
    geom
end

"""
Force rings to be closed.

If this geometry, or any contained geometries has polygon rings that are not
closed, they will be closed by adding the starting point at the end.
"""
closerings!(geom::G) where {G <: AbstractGeometry} = (GDAL.closerings(geom.ptr); geom)

"""
Create geometry from GML.

This method translates a fragment of GML containing only the geometry portion
into a corresponding OGRGeometry. There are many limitations on the forms of GML
geometries supported by this parser, but they are too numerous to list here.

The following GML2 elements are parsed : Point, LineString, Polygon, MultiPoint,
MultiLineString, MultiPolygon, MultiGeometry.
"""
unsafe_fromGML(::Type{G}, data) where {G <: AbstractGeometry} =
    G(GDAL.createfromgml(data))

unsafe_fromGML(data) = unsafe_fromGML(Geometry, data)

fromGML(data) = unsafe_fromGML(IGeometry, data)

"Convert a geometry into GML format."
toGML(geom::AbstractGeometry) = GDAL.exporttogml(geom.ptr)

"Convert a geometry into KML format."
# * `altitudemode`: value to write in altitudeMode element, or NULL.
toKML(geom::AbstractGeometry, altitudemode = Ptr{UInt8}(C_NULL)) =
    GDAL.exporttokml(geom.ptr, altitudemode)

"Convert a geometry into GeoJSON format."
toJSON(geom::AbstractGeometry) = GDAL.exporttojson(geom.ptr)

"""
Convert a geometry into GeoJSON format.
### Parameters
* `geom`: handle to the geometry.
* `options`: a list of options.
### Returns
A GeoJSON fragment or NULL in case of error.
"""
toJSON(geom::AbstractGeometry, options) =
    unsafe_string(@gdal(OGR_G_ExportToJsonEx::Cstring,
        geom.ptr::GDALGeometry,
        options::StringList
    ))

"Create a geometry object from its GeoJSON representation"
unsafe_fromJSON(::Type{G}, data::String) where {G <: AbstractGeometry} =
    G(GDAL.creategeometryfromjson(data))

unsafe_fromJSON(data::String) = unsafe_fromJSON(Geometry, data)

fromJSON(data::String) = unsafe_fromJSON(IGeometry, data)

"Assign spatial reference to this object."
setspatialref!(geom::G, spatialref::SpatialRef) where {G <: AbstractGeometry} =
    (GDAL.assignspatialreference(geom.ptr, spatialref.ptr); geom)

"""
Returns spatial reference system for geometry.

The object may be shared with many geometry objects, and should not be modified.
"""
getspatialref(geom::AbstractGeometry) =
    SpatialRef(GDAL.getspatialreference(geom.ptr))

"""
Apply arbitrary coordinate transformation to geometry.

### Parameters
* `geom`: handle on the geometry to apply the transform to.
* `coordtransform`: handle on the transformation to apply.
"""
function transform!(geom::G,
        coordtransform::CoordTransform
    ) where G <: AbstractGeometry
    result = GDAL.transform(geom.ptr, coordtransform.ptr)
    @ogrerr result "Failed to transform geometry"
    geom
end

"Transform geometry to new spatial reference system."
function transform!(geom::G, spatialref::SpatialRef) where G <: AbstractGeometry
    result = GDAL.transformto(geom.ptr, spatialref.ptr)
    @ogrerr result "Failed to transform geometry to the new SRS"
    geom
end

"""
Compute a simplified geometry.

### Parameters
* `geom`: the geometry.
* `tol`: the distance tolerance for the simplification.
"""
function unsafe_simplify(::Type{G},
        geom::AbstractGeometry, tol::Real
    ) where G <: AbstractGeometry
    G(GDAL.simplify(geom.ptr, tol))
end

unsafe_simplify(geom::AbstractGeometry, args...) =
    unsafe_simplify(Geometry, geom, args...)

simplify(geom::AbstractGeometry, args...) =
    unsafe_simplify(IGeometry, geom, args...)

"""
Simplify the geometry while preserving topology.

### Parameters
* `geom`: the geometry.
* `tol`: the distance tolerance for the simplification.
"""
function unsafe_simplifypreservetopology(::Type{G},
        geom::AbstractGeometry, tol::Real
    ) where G <: AbstractGeometry
    G(GDAL.simplifypreservetopology(geom.ptr, tol))
end

unsafe_simplifypreservetopology(geom::AbstractGeometry, args...) =
    unsafe_simplifypreservetopology(Geometry, geom, args...)

simplifypreservetopology(geom::AbstractGeometry, args...) =
    unsafe_simplifypreservetopology(IGeometry, geom, args...)

"""
Return a Delaunay triangulation of the vertices of the geometry.

### Parameters
* `geom`: the geometry.
* `tol`: optional snapping tolerance to use for improved robustness
* `onlyedges`: if TRUE, will return a MULTILINESTRING, otherwise it
    will return a GEOMETRYCOLLECTION containing triangular POLYGONs.
"""
function unsafe_delaunaytriangulation(::Type{G},
        geom::AbstractGeometry, tol::Real, onlyedges::Bool
    ) where G <: AbstractGeometry
    G(GDAL.delaunaytriangulation(geom.ptr, tol, onlyedges))
end

unsafe_delaunaytriangulation(geom::AbstractGeometry, args...) =
    unsafe_delaunaytriangulation(Geometry, geom, args...)

delaunaytriangulation(geom::AbstractGeometry, args...) =
    unsafe_delaunaytriangulation(IGeometry, geom, args...)

"""
Modify the geometry such it has no segment longer than the given distance.

Interpolated points will have Z and M values (if needed) set to 0. Distance
computation is performed in 2d only

### Parameters
* `geom`: the geometry to segmentize
* `maxlength`: the maximum distance between 2 points after segmentization
"""
function segmentize!(geom::G, maxlength::Real) where G <: AbstractGeometry
    GDAL.segmentize(geom.ptr, maxlength)
    geom
end

"""
Returns whether the geometries intersect

Determines whether two geometries intersect. If GEOS is enabled, then this is
done in rigorous fashion otherwise TRUE is returned if the envelopes (bounding
boxes) of the two geometries overlap.
"""
function intersects(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.intersects(g1.ptr, g2.ptr))
end

"Returns TRUE if the geometries are equivalent."
function equals(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.equals(g1.ptr, g2.ptr))
end

"Returns TRUE if the geometries are disjoint."
function disjoint(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.disjoint(g1.ptr, g2.ptr))
end

"Returns TRUE if the geometries are touching."
function touches(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.touches(g1.ptr, g2.ptr))
end

"Returns TRUE if the geometries are crossing."
function crosses(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.crosses(g1.ptr, g2.ptr))
end

"Returns TRUE if g1 is contained within g2."
function within(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.within(g1.ptr, g2.ptr))
end

"Returns TRUE if g1 contains g2."
function contains(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.contains(g1.ptr, g2.ptr))
end

"Returns TRUE if the geometries overlap."
function overlaps(g1::AbstractGeometry, g2::AbstractGeometry)
    Bool(GDAL.overlaps(g1.ptr, g2.ptr))
end

"""
Returns the boundary of the geometry.

A new geometry object is created and returned containing the boundary of the
geometry on which the method is invoked.
"""
unsafe_boundary(::Type{G}, geom::AbstractGeometry) where {G <: AbstractGeometry} =
    G(GDAL.boundary(geom.ptr))

unsafe_boundary(geom::AbstractGeometry, args...) =
    unsafe_boundary(Geometry, geom, args...)

boundary(geom::AbstractGeometry, args...) =
    unsafe_boundary(IGeometry, geom, args...)

"""
Returns the convex hull of the geometry.

A new geometry object is created and returned containing the convex hull of the
geometry on which the method is invoked.
"""
unsafe_convexhull(::Type{G}, geom::AbstractGeometry) where {G <: AbstractGeometry} =
    G(GDAL.convexhull(geom.ptr))

unsafe_convexhull(geom::AbstractGeometry, args...) =
    unsafe_convexhull(Geometry, geom, args...)

convexhull(geom::AbstractGeometry, args...) =
    unsafe_convexhull(IGeometry, geom, args...)

"""
Compute buffer of geometry.

Builds a new geometry containing the buffer region around the geometry on which
it is invoked. The buffer is a polygon containing the region within the buffer
distance of the original geometry.

Some buffer sections are properly described as curves, but are converted to
approximate polygons. The nQuadSegs parameter can be used to control how many
segments should be used to define a 90 degree curve - a quadrant of a circle.
A value of 30 is a reasonable default. Large values result in large numbers of
vertices in the resulting buffer geometry while small numbers reduce the
accuracy of the result.

### Parameters
* `geom`: the geometry.
* `dist`: the buffer distance to be applied. Should be expressed into the
    same unit as the coordinates of the geometry.
* `quadsegs`: the number of segments used to approximate a 90 degree
    (quadrant) of curvature.
"""
function unsafe_buffer(::Type{G},
        geom::AbstractGeometry, dist::Real, quadsegs::Integer = 30
    ) where G <: AbstractGeometry
    G(GDAL.buffer(geom.ptr, dist, quadsegs))
end

unsafe_buffer(geom::AbstractGeometry, args...) =
    unsafe_buffer(Geometry, geom, args...)

buffer(geom::AbstractGeometry, args...) =
    unsafe_buffer(IGeometry, geom, args...)

"""
Returns a new geometry representing the intersection of the geometries, or NULL
if there is no intersection or an error occurs.

Generates a new geometry which is the region of intersection of the two
geometries operated on. The OGR_G_Intersects() function can be used to test if
two geometries intersect.
"""
function unsafe_intersection(::Type{G},
        g1::AbstractGeometry, g2::AbstractGeometry
    ) where G <: AbstractGeometry
    G(GDAL.intersection(g1.ptr, g2.ptr))
end

unsafe_intersection(g1::AbstractGeometry, g2::AbstractGeometry, args...) =
    unsafe_intersection(Geometry, g1, g2, args...)

intersection(g1::AbstractGeometry, g2::AbstractGeometry, args...) =
    unsafe_intersection(IGeometry, g1, g2, args...)

"""
Returns a new geometry representing the union of the geometries.

Generates a new geometry which is the region of union of the two geometries
operated on.
"""
function unsafe_union(::Type{G},
        g1::AbstractGeometry, g2::AbstractGeometry
    ) where G <: AbstractGeometry
    G(GDAL.union(g1.ptr, g2.ptr))
end

"Compute the union of the geometry using cascading."
unsafe_union(::Type{G}, geom::AbstractGeometry) where {G <: AbstractGeometry} =
    G(GDAL.unioncascaded(geom.ptr))

unsafe_union(geom::AbstractGeometry, args...) =
    unsafe_union(Geometry, geom, args...)

union(geom::AbstractGeometry, args...) =
    unsafe_union(IGeometry, geom, args...)

"""
Returns a point guaranteed to lie on the surface.

This method relates to the SFCOM ISurface::get_PointOnSurface() method however
the current implementation based on GEOS can operate on other geometry types
than the types that are supported by SQL/MM-Part 3 : surfaces (polygons) and
multisurfaces (multipolygons).
"""
function unsafe_pointonsurface(::Type{G},
        geom::AbstractGeometry
    ) where G <: AbstractGeometry
    G(GDAL.pointonsurface(geom.ptr))
end

unsafe_pointonsurface(geom::AbstractGeometry, args...) =
    unsafe_pointonsurface(Geometry, geom, args...)

pointonsurface(geom::AbstractGeometry, args...) =
    unsafe_pointonsurface(IGeometry, geom, args...)

"""
Generates a new geometry which is the region of this geometry with the region of
the other geometry removed.

### Returns
A new geometry representing the difference of the geometries, or NULL
if the difference is empty.
"""
function unsafe_difference(::Type{G},
        g1::AbstractGeometry, g2::AbstractGeometry
    ) where G <: AbstractGeometry
    G(GDAL.difference(g1.ptr, g2.ptr))
end

unsafe_difference(geom::AbstractGeometry, args...) =
    unsafe_difference(Geometry, geom, args...)

difference(geom::AbstractGeometry, args...) =
    unsafe_difference(IGeometry, geom, args...)

"""
Returns a new geometry representing the symmetric difference of the geometries
or NULL if the difference is empty or an error occurs.
"""
function unsafe_symdifference(::Type{G},
        g1::AbstractGeometry, g2::AbstractGeometry
    ) where G <: AbstractGeometry
    G(GDAL.symdifference(g1.ptr, g2.ptr))
end

unsafe_symdifference(geom::AbstractGeometry, args...) =
    unsafe_symdifference(Geometry, geom, args...)

symdifference(geom::AbstractGeometry, args...) =
    unsafe_symdifference(IGeometry, geom, args...)

"Returns the distance between the geometries or -1 if an error occurs."
distance(g1::AbstractGeometry, g2::AbstractGeometry) = GDAL.distance(g1.ptr, g2.ptr)

"Returns the length of the geometry, or 0.0 for unsupported geometry types."
geomlength(geom::AbstractGeometry) = GDAL.length(geom.ptr)

"Returns the area of the geometry or 0.0 for unsupported geometry types."
geomarea(geom::AbstractGeometry) = GDAL.area(geom.ptr)

"""
Compute the geometry centroid.

The centroid location is applied to the passed in OGRPoint object. The centroid
is not necessarily within the geometry.

This method relates to the SFCOM ISurface::get_Centroid() method however the
current implementation based on GEOS can operate on other geometry types such as
multipoint, linestring, geometrycollection such as multipolygons. OGC SF SQL 1.1
defines the operation for surfaces (polygons). SQL/MM-Part 3 defines the
operation for surfaces and multisurfaces (multipolygons).
"""
function centroid!(
        geom::AbstractGeometry, centroid::G
    ) where G <: AbstractGeometry
    result = GDAL.centroid(geom.ptr, centroid.ptr)
    @ogrerr result "Failed to compute the geometry centroid"
    centroid
end

unsafe_centroid(geom::AbstractGeometry) = centroid!(geom, unsafe_createpoint())

centroid(geom::AbstractGeometry) = centroid!(geom, createpoint())

"""
Fetch point at given distance along curve.

### Parameters
* `geom`: curve geometry.
* `distance`: distance along the curve at which to sample position. This
    distance should be between zero and geomlength() for this curve.

### Returns
a point or NULL.
"""
function unsafe_pointalongline(::Type{G},
        geom::AbstractGeometry, distance::Real
    ) where G <: AbstractGeometry
    G(GDAL.value(geom.ptr, distance))
end

unsafe_pointalongline(geom::AbstractGeometry, args...) =
    unsafe_pointalongline(Geometry, geom, args...)

pointalongline(geom::AbstractGeometry, args...) =
    unsafe_pointalongline(IGeometry, geom, args...)

"""
Clear geometry information.

This restores the geometry to its initial state after construction, and before
assignment of actual geometry.
"""
empty!(geom::G) where {G <: AbstractGeometry} = (GDAL.empty(geom.ptr); geom)

"Returns TRUE if the geometry has no points, otherwise FALSE."
isempty(geom::AbstractGeometry) = Bool(GDAL.isempty(geom.ptr))

"Returns TRUE if the geometry is valid, otherwise FALSE."
isvalid(geom::AbstractGeometry) = Bool(GDAL.isvalid(geom.ptr))

"Returns TRUE if the geometry is simple, otherwise FALSE."
issimple(geom::AbstractGeometry) = Bool(GDAL.issimple(geom.ptr))

"Returns TRUE if the geometry is a ring, otherwise FALSE."
isring(geom::AbstractGeometry) = Bool(GDAL.isring(geom.ptr))

"""
Polygonizes a set of sparse edges.

A new geometry object is created and returned containing a collection of
reassembled Polygons: NULL will be returned if the input collection doesn't
correspond to a MultiLinestring, or when reassembling Edges into Polygons is
impossible due to topological inconsistencies.

### Returns
a handle to a newly allocated geometry now owned by the caller, or NULL on
failure.
"""
unsafe_polygonize(::Type{G}, geom::AbstractGeometry) where {G <: AbstractGeometry} =
    G(GDAL.polygonize(geom.ptr))

unsafe_polygonize(geom::AbstractGeometry, args...) =
    unsafe_polygonize(Geometry, geom, args...)

polygonize(geom::AbstractGeometry, args...) =
    unsafe_polygonize(IGeometry, geom, args...)

# """
#     OGR_G_GetPoints(OGRGeometryH hGeom,
#                     void * pabyX,
#                     int nXStride,
#                     void * pabyY,
#                     int nYStride,
#                     void * pabyZ,
#                     int nZStride) -> int
# Returns all points of line string.
# ### Parameters
# * `hGeom`: handle to the geometry from which to get the coordinates.
# * `pabyX`: a buffer of at least (sizeof(double) * nXStride * nPointCount)
# bytes, may be NULL.
# * `nXStride`: the number of bytes between 2 elements of pabyX.
# * `pabyY`: a buffer of at least (sizeof(double) * nYStride * nPointCount)
# bytes, may be NULL.
# * `nYStride`: the number of bytes between 2 elements of pabyY.
# * `pabyZ`: a buffer of at last size (sizeof(double) * nZStride * nPointCount)
# bytes, may be NULL.
# * `nZStride`: the number of bytes between 2 elements of pabyZ.
# ### Returns
# the number of points
# """
# function getpoints(hGeom::Ptr{OGRGeometryH},pabyX,nXStride::Integer,pabyY,
# nYStride::Integer,pabyZ,nZStride::Integer)
#     ccall((:OGR_G_GetPoints,libgdal),Cint,(Ptr{OGRGeometryH},Ptr{Void},Cint,
# Ptr{Void},Cint,Ptr{Void},Cint),hGeom,pabyX,nXStride,pabyY,nYStride,pabyZ,
# nZStride)
# end

"Fetch the x coordinate of a point from a geometry, at index i."
getx(geom::AbstractGeometry, i::Integer) = GDAL.getx(geom.ptr, i)

"Fetch the y coordinate of a point from a geometry, at index i."
gety(geom::AbstractGeometry, i::Integer) = GDAL.gety(geom.ptr, i)

"Fetch the z coordinate of a point from a geometry, at index i."
getz(geom::AbstractGeometry, i::Integer) = GDAL.getz(geom.ptr, i)

"""
Fetch a point in line string or a point geometry, at index i.

### Parameters
* `i`: the vertex to fetch, from 0 to getNumPoints()-1, zero for a point.
"""
getpoint!(geom::AbstractGeometry, i::Integer, x, y, z) =
    (GDAL.getpoint(geom.ptr, i, x, y, z); (x[], y[], z[]))

getpoint(geom::AbstractGeometry, i::Integer) =
    getpoint!(geom, i, Ref{Cdouble}(), Ref{Cdouble}(), Ref{Cdouble}())

"""
Set number of points in a geometry.

### Parameters
* `geom`: the geometry.
* `n`: the new number of points for geometry.
"""
setpointcount!(geom::AbstractGeometry, n::Integer) =
    (GDAL.setpointcount(geom.ptr, n); geom)

"""
Set the location of a vertex in a point or linestring geometry.

### Parameters
* `geom`: handle to the geometry to add a vertex to.
* `i`: the index of the vertex to assign (zero based) or zero for a point.
* `x`: input X coordinate to assign.
* `y`: input Y coordinate to assign.
* `z`: input Z coordinate to assign (defaults to zero).
"""
setpoint!(geom::AbstractGeometry, i::Integer, x::Real, y::Real, z::Real) =
    (GDAL.setpoint(geom.ptr, i, x, y, z); geom)

"""
Set the location of a vertex in a point or linestring geometry.

### Parameters
* `geom`: handle to the geometry to add a vertex to.
* `i`: the index of the vertex to assign (zero based) or zero for a point.
* `x`: input X coordinate to assign.
* `y`: input Y coordinate to assign.
"""
setpoint!(geom::AbstractGeometry, i::Integer, x::Real, y::Real) =
    (GDAL.setpoint_2d(geom.ptr, i, x, y); geom)

"""
Add a point to a geometry (line string or point).

### Parameters
* `geom`: the geometry to add a point to.
* `x`: x coordinate of point to add.
* `y`: y coordinate of point to add.
* `z`: z coordinate of point to add.
"""
addpoint!(geom::AbstractGeometry, x::Real, y::Real, z::Real) =
    (GDAL.addpoint(geom.ptr, x, y, z); geom)

"""
Add a point to a geometry (line string or point).

### Parameters
* `geom`: the geometry to add a point to.
* `x`: x coordinate of point to add.
* `y`: y coordinate of point to add.
"""
addpoint!(geom::AbstractGeometry, x::Real, y::Real) =
    (GDAL.addpoint_2d(geom.ptr, x, y); geom)

# """
#     OGR_G_SetPoints(OGRGeometryH hGeom,
#                     int nPointsIn,
#                     void * pabyX,
#                     int nXStride,
#                     void * pabyY,
#                     int nYStride,
#                     void * pabyZ,
#                     int nZStride) -> void
# Assign all points in a point or a line string geometry.
# ### Parameters
# * `hGeom`: handle to the geometry to set the coordinates.
# * `nPointsIn`: number of points being passed in padfX and padfY.
# * `pabyX`: list of X coordinates (double values) of points being assigned.
# * `nXStride`: the number of bytes between 2 elements of pabyX.
# * `pabyY`: list of Y coordinates (double values) of points being assigned.
# * `nYStride`: the number of bytes between 2 elements of pabyY.
# * `pabyZ`: list of Z coordinates (double values) of points being assigned
# (defaults to NULL for 2D objects).
# * `nZStride`: the number of bytes between 2 elements of pabyZ.
# """
# function setpoints(hGeom::Ptr{OGRGeometryH},nPointsIn::Integer,pabyX,
# nXStride::Integer,pabyY,nYStride::Integer,pabyZ,nZStride::Integer)
#     ccall((:OGR_G_SetPoints,libgdal),Void,(Ptr{OGRGeometryH},Cint,Ptr{Void},
# Cint,Ptr{Void},Cint,Ptr{Void},Cint),hGeom,nPointsIn,pabyX,nXStride,pabyY,
# nYStride,pabyZ,nZStride)
# end

"""
The number of elements in a geometry or number of geometries in container.

This corresponds to

* `OGR_G_GetPointCount` for wkbPoint[25D] or wkbLineString[25D],
* `OGR_G_GetGeometryCount` for geometries of type wkbPolygon[25D], wkbMultiPoint[25D], wkbMultiLineString[25D], wkbMultiPolygon[25D] or wkbGeometryCollection[25D], and
* `0` for other geometry types.
"""
function ngeom(geom::AbstractGeometry)
    n = GDAL.getpointcount(geom.ptr)
    n == 0 ? GDAL.getgeometrycount(geom.ptr) : n
end

"""
Fetch geometry from a geometry container.

This function returns an handle to a geometry within the container. The returned
geometry remains owned by the container, and should not be modified. The handle
is only valid until the next change to the geometry container. Use OGR_G_Clone()
to make a copy.

For a polygon, `getgeom(polygon,i)` returns the exterior ring if
`i == 0`, and the interior rings for `i > 0`.

### Parameters
* `geom`: the geometry container from which to get a geometry from.
* `i`: index of the geometry to fetch, between 0 and getNumGeometries() - 1.
"""
function unsafe_getgeom(::Type{G},
        geom::AbstractGeometry, i::Integer
    ) where G <: AbstractGeometry
    G(GDAL.getgeometryref(geom.ptr, i))
end

unsafe_getgeom(geom::AbstractGeometry, args...) =
    unsafe_getgeom(Geometry, geom, args...)

getgeom(geom::AbstractGeometry, args...) =
    unsafe_getgeom(IGeometry, unsafe_clone(geom), args...)

"""
Add a geometry to a geometry container.

Some subclasses of OGRGeometryCollection restrict the types of geometry that can
be added, and may return an error. The passed geometry is cloned to make an
internal copy.

For a polygon, `subgeom` must be a linearring. If the polygon is empty, the
first added subgeometry will be the exterior ring. The next ones will be the
interior rings.

### Parameters
* `geomcontainer`: existing geometry.
* `subgeom`: geometry to add to the existing geometry.
"""
function addgeom!(geomcontainer::G, subgeom::AbstractGeometry) where G <: AbstractGeometry
    result = GDAL.addgeometry(geomcontainer.ptr, subgeom.ptr)
    @ogrerr result "Failed to add geometry. The geometry type could be illegal"
    geomcontainer
end

"""
Add a geometry directly to an existing geometry container.

Some subclasses of OGRGeometryCollection restrict the types of geometry that can
be added, and may return an error. Ownership of the passed geometry is taken by
the container rather than cloning as addGeometry() does.

For a polygon, hNewSubGeom must be a linearring. If the polygon is empty, the
first added subgeometry will be the exterior ring. The next ones will be the
interior rings.

### Parameters
* `geomcontainer`: existing geometry.
* `subgeom`: geometry to add to the existing geometry.
"""
function addgeomdirectly!(geomcontainer::G, subgeom::AbstractGeometry) where G <: AbstractGeometry
    result = GDAL.addgeometrydirectly(geomcontainer.ptr, subgeom.ptr)
    @ogrerr result "Failed to add geometry. The geometry type could be illegal"
    geomcontainer
end

"""
Remove a geometry from an exiting geometry container.

### Parameters
* `geom`: the existing geometry to delete from.
* `i`: the index of the geometry to delete. A value of -1 is a special flag
    meaning that all geometries should be removed.
* `todelete`: if TRUE the geometry will be destroyed, otherwise it will not.
    The default is TRUE as the existing geometry is considered to own the
    geometries in it.
"""
function removegeom!(geom::AbstractGeometry, i::Integer, todelete::Bool = true)
    result = GDAL.removegeometry(geom.ptr, i, todelete)
    @ogrerr result "Failed to remove geometry. The index could be out of range."
    geom
end

"""
Remove all geometries from an exiting geometry container.

### Parameters
* `geom`: the existing geometry to delete from.
* `todelete`: if TRUE the geometry will be destroyed, otherwise it will not.
    The default is TRUE as the existing geometry is considered to own the
    geometries in it.
"""
function removeallgeoms!(geom::AbstractGeometry, todelete::Bool = true)
    result = GDAL.removegeometry(geom.ptr, -1, todelete)
    @ogrerr result "Failed to remove all geometries."
    geom
end

"""
Returns if this geometry is or has curve geometry.

### Parameters
* `geom`: the geometry to operate on.
* `nonlinear`: set it to TRUE to check if the geometry is or contains a
    CIRCULARSTRING.
"""
hascurvegeom(geom::AbstractGeometry, nonlinear::Bool) =
    Bool(GDAL.hascurvegeometry(geom.ptr, nonlinear))

"""
Return, possibly approximate, linear version of this geometry.

Returns a geometry that has no CIRCULARSTRING, COMPOUNDCURVE, CURVEPOLYGON,
MULTICURVE or MULTISURFACE in it, by approximating curve geometries.

The ownership of the returned geometry belongs to the caller.

### Parameters
* `geom`: the geometry to operate on.
* `stepsize`: the largest step in degrees along the arc, zero to use the
    default setting.
* `options`: options as a null-terminated list of strings or NULL.
    See OGRGeometryFactory::curveToLineString() for valid options.
"""
function unsafe_getlineargeom(::Type{G},
        geom::AbstractGeometry, stepsize::Real = 0
    ) where G <: AbstractGeometry
    G(GDAL.getlineargeometry(geom.ptr, stepsize, C_NULL))
end

function unsafe_getlineargeom(::Type{G},
        geom::AbstractGeometry, options::Vector
    ) where G <: AbstractGeometry
    G(GDAL.getlineargeometry(geom.ptr, 0, options))
end

function unsafe_getlineargeom(::Type{G},
        geom::AbstractGeometry, stepsize::Real, options::Vector
    ) where G <: AbstractGeometry
    G(GDAL.getlineargeometry(geom.ptr, stepsize, options))
end

unsafe_getlineargeom(geom::AbstractGeometry, args...) =
    unsafe_getlineargeom(Geometry, geom, args...)

getlineargeom(geom::AbstractGeometry, args...) =
    unsafe_getlineargeom(IGeometry, geom, args...)

"""
Return curve version of this geometry.

Returns a geometry that has possibly CIRCULARSTRING, COMPOUNDCURVE,
CURVEPOLYGON, MULTICURVE or MULTISURFACE in it, by de-approximating linear into
curve geometries.

If the geometry has no curve portion, the returned geometry will be a clone.

The ownership of the returned geometry belongs to the caller.

The reverse function is OGR_G_GetLinearGeometry().
"""
unsafe_getcurvegeom(::Type{G}, geom::AbstractGeometry) where {G <: AbstractGeometry} =
    G(GDAL.getcurvegeometry(geom.ptr, C_NULL))

unsafe_getcurvegeom(geom::AbstractGeometry, args...) =
    unsafe_getcurvegeom(Geometry, geom, args...)

getcurvegeom(geom::AbstractGeometry, args...) =
    unsafe_getcurvegeom(IGeometry, geom, args...)

"""
Build a ring from a bunch of arcs.

### Parameters
* `lines`: handle to an OGRGeometryCollection (or OGRMultiLineString)
    containing the line string geometries to be built into rings.
* `besteffort`: not yet implemented???.
* `autoclose`: indicates if the ring should be close when first and last
    points of the ring are the same.
* `tol`: whether two arcs are considered close enough to be joined.
"""
function unsafe_polygonfromedges(
        ::Type{G},
        lines::AbstractGeometry,
        besteffort::Bool,
        autoclose::Bool,
        tol::Real
    ) where G <: AbstractGeometry
    perr = Ref{GDAL.OGRErr}()
    result = GDAL.buildpolygonfromedges(lines, besteffort, autoclose, tol, perr)
    @ogrerr perr[] "Failed to build polygon from edges."
    G(result)
end

unsafe_polygonfromedges(geom::AbstractGeometry, args...) =
    unsafe_polygonfromedges(Geometry, geom, args...)

polygonfromedges(geom::AbstractGeometry, args...) =
    unsafe_polygonfromedges(IGeometry, geom, args...)

"""
Set flag to enable/disable returning non-linear geometries in the C API.

This flag has only an effect on the OGR_F_GetGeometryRef(),
OGR_F_GetGeomFieldRef(), OGR_L_GetGeomType(), OGR_GFld_GetType() and
OGR_FD_GetGeomType() C API methods. It is meant as making it simple for
applications using the OGR C API not to have to deal with non-linear geometries,
even if such geometries might be returned by drivers. In which case, they will
be transformed into their closest linear geometry, by doing linear
approximation, with OGR_G_ForceTo().

Libraries should generally not use that method, since that could interfere with
other libraries or applications.

### Parameters
* `flag`: TRUE if non-linear geometries might be returned (default value).
          FALSE to ask for non-linear geometries to be approximated as linear
          geometries.

### Returns
a point or NULL.
"""
setnonlineargeomflag!(flag::Bool) = GDAL.setnonlineargeometriesenabledflag(flag)

"Get flag to enable/disable returning non-linear geometries in the C API."
getnonlineargeomflag() = Bool(GDAL.getnonlineargeometriesenabledflag())

for (geom, wkbgeom) in ((:geomcollection,   GDAL.wkbGeometryCollection),
                        (:linestring,       GDAL.wkbLineString),
                        (:linearring,       GDAL.wkbLinearRing),
                        (:multilinestring,  GDAL.wkbMultiLineString),
                        (:multipoint,       GDAL.wkbMultiPoint),
                        (:multipolygon,     GDAL.wkbMultiPolygon),
                        (:point,            GDAL.wkbPoint),
                        (:polygon,          GDAL.wkbPolygon))
    eval(quote
        function $(Symbol("unsafe_create$geom"))(::Type{G}) where G <: AbstractGeometry
            unsafe_creategeom(G, $wkbgeom)
        end

        $(Symbol("unsafe_create$geom"))() = $(Symbol("unsafe_create$geom"))(Geometry)

        $(Symbol("create$geom"))() = $(Symbol("unsafe_create$geom"))(IGeometry)
    end)
end

function unsafe_createpoint(::Type{G}, coords::Vector{<:Real}) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbPoint)
    addpoint!(geom, coords...)
    geom
end

unsafe_createpoint(coords::Vector{<:Real}) = unsafe_createpoint(Geometry, coords)

createpoint(coords::Vector{<:Real}) = unsafe_createpoint(IGeometry, coords)

function unsafe_createpoint(::Type{G}, x::Real, y::Real) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbPoint)
    addpoint!(geom, x, y)
    geom
end

unsafe_createpoint(x::Real, y::Real) = unsafe_createpoint(Geometry, x, y)

createpoint(x::Real, y::Real) = unsafe_createpoint(IGeometry, x, y)

function unsafe_createpoint(::Type{G}, x::Real, y::Real, z::Real) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbPoint)
    addpoint!(geom, x, y, z)
    geom
end

unsafe_createpoint(x::Real, y::Real, z::Real) =
    unsafe_createpoint(Geometry, x, y, z)

createpoint(x::Real, y::Real, z::Real) =
    unsafe_createpoint(IGeometry, x, y, z)

function unsafe_createpoint(::Type{G}, xy::Tuple{T,U}) where {T <: Real, U <: Real, G <: AbstractGeometry}
    geom = unsafe_creategeom(G, GDAL.wkbPoint)
    addpoint!(geom, xy...)
    geom
end

unsafe_createpoint(xy::Tuple{T,U}) where {T <: Real, U <: Real} =
    unsafe_createpoint(Geometry, xy)

createpoint(xy::Tuple{T,U}) where {T <: Real, U <: Real} =
    unsafe_createpoint(IGeometry, xy)

function unsafe_createpoint(::Type{G}, xyz::Tuple{T,U,V}) where {T <: Real, U <: Real, V <: Real, G <: AbstractGeometry}
    geom = unsafe_creategeom(G, GDAL.wkbPoint)
    addpoint!(geom, xyz...)
    geom
end

unsafe_createpoint(xyz::Tuple{T,U,V}) where {T <: Real, U <: Real, V <: Real} =
    unsafe_createpoint(Geometry, xyz)

createpoint(args...) = unsafe_createpoint(IGeometry, args...)

# Tuples of Vectors
for (geom, wkbgeom) in ((:linestring, GDAL.wkbLineString),
                        (:linearring, GDAL.wkbLinearRing))
    eval(quote
        function $(Symbol("unsafe_create$geom"))(
                ::Type{G},
                xs::Vector{Cdouble},
                ys::Vector{Cdouble}
            ) where G <: AbstractGeometry
            geom = unsafe_creategeom(G, $wkbgeom)
            for (x,y) in zip(xs, ys)
                addpoint!(geom, x, y)
            end
            geom
        end

        function $(Symbol("unsafe_create$geom"))(
                ::Type{G},
                xs::Vector{Cdouble},
                ys::Vector{Cdouble},
                zs::Vector{Cdouble}
            ) where G <: AbstractGeometry
            geom = unsafe_creategeom(G, $wkbgeom)
            for (x,y,z) in zip(xs, ys, zs)
                addpoint!(geom, x, y, z)
            end
            geom
        end

        $(Symbol("unsafe_create$geom"))(args...) = $(Symbol("unsafe_create$geom"))(Geometry, args...)

        $(Symbol("create$geom"))(args...) = $(Symbol("unsafe_create$geom"))(IGeometry, args...)
    end)
end

function unsafe_createpolygon(
        ::Type{G},
        xs::Vector{Cdouble},
        ys::Vector{Cdouble}
    ) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbPolygon)
    addgeomdirectly!(geom, unsafe_createlinearring(xs, ys))
    geom
end

function unsafe_createpolygon(
        ::Type{G},
        xs::Vector{Cdouble},
        ys::Vector{Cdouble},
        zs::Vector{Cdouble}
    ) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbPolygon)
    addgeomdirectly!(geom, unsafe_createlinearring(xs, ys, zs))
    geom
end

unsafe_createpolygon(args...) = unsafe_createpolygon(Geometry, args...)

createpolygon(args...) = unsafe_createpolygon(IGeometry, args...)

function unsafe_createmultipoint(
        ::Type{G},
        xs::Vector{Cdouble},
        ys::Vector{Cdouble}
    ) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbMultiPoint)
    for (x, y) in zip(xs, ys)
        addgeomdirectly!(geom, unsafe_createpoint(x, y))
    end
    geom
end

function unsafe_createmultipoint(
        ::Type{G},
        xs::Vector{Cdouble},
        ys::Vector{Cdouble},
        zs::Vector{Cdouble}
    ) where G <: AbstractGeometry
    geom = unsafe_creategeom(G, GDAL.wkbMultiPoint)
    for (x, y, z) in zip(xs, ys, zs)
        addgeomdirectly!(geom, unsafe_createpoint(x, y, z))
    end
    geom
end

unsafe_createmultipoint(args...) = unsafe_createmultipoint(Geometry, args...)

createmultipoint(args...) = unsafe_createmultipoint(IGeometry, args...)

# Vectors of Tuples
for typeargs in (Vector{Tuple{Cdouble,Cdouble}},
                 Vector{Tuple{Cdouble,Cdouble,Cdouble}},
                 Vector{Vector{Cdouble}})
    for (geom, wkbgeom) in ((:linestring, GDAL.wkbLineString),
                            (:linearring, GDAL.wkbLinearRing))
        eval(quote
            function $(Symbol("unsafe_create$geom"))(
                    ::Type{G},
                    coords::$typeargs
                ) where G <: AbstractGeometry
                geom = unsafe_creategeom(G, $wkbgeom)
                for coord in coords; addpoint!(geom, coord...) end
                geom
            end

            $(Symbol("unsafe_create$geom"))(coords::$typeargs) =
                $(Symbol("unsafe_create$geom"))(Geometry, coords)

            $(Symbol("create$geom"))(coords::$typeargs) =
                $(Symbol("unsafe_create$geom"))(IGeometry, coords)
        end)
    end

    eval(quote
        function unsafe_createpolygon(::Type{G}, coords::$typeargs) where G <: AbstractGeometry
            geom = unsafe_creategeom(G, GDAL.wkbPolygon)
            addgeomdirectly!(geom, unsafe_createlinearring(coords))
            geom
        end

        unsafe_createpolygon(coords::$typeargs) = unsafe_createpolygon(Geometry, coords)

        createpolygon(coords::$typeargs) = unsafe_createpolygon(IGeometry, coords)

        function unsafe_createmultipoint(::Type{G}, coords::$typeargs) where G <: AbstractGeometry
            geom = unsafe_creategeom(G, GDAL.wkbMultiPoint)
            for point in coords
                addgeomdirectly!(geom, unsafe_createpoint(point))
            end
            geom
        end

        unsafe_createmultipoint(coords::$typeargs) = unsafe_createmultipoint(Geometry, coords)

        createmultipoint(coords::$typeargs) = unsafe_createmultipoint(IGeometry, coords)
    end)
end

for typeargs in (Vector{Vector{Tuple{Cdouble,Cdouble}}},
                 Vector{Vector{Tuple{Cdouble,Cdouble,Cdouble}}},
                 Vector{Vector{Vector{Cdouble}}})
    for (geom, wkbgeom, f) in ((:polygon, GDAL.wkbPolygon, :unsafe_createlinearring),
                               (:multilinestring, GDAL.wkbMultiLineString, :unsafe_createlinestring),
                               (:multipolygon_noholes, GDAL.wkbMultiPolygon, :unsafe_createpolygon))
        eval(quote
            function $(Symbol("unsafe_create$geom"))(
                    ::Type{G},
                    coords::$typeargs
                ) where G <: AbstractGeometry
                geom = unsafe_creategeom(G, $wkbgeom)
                for coord in coords; addgeomdirectly!(geom, $(f)(coord)) end
                geom
            end

            $(Symbol("unsafe_create$geom"))(coords::$typeargs) =
                $(Symbol("unsafe_create$geom"))(Geometry, coords::$typeargs)

            $(Symbol("create$geom"))(coords::$typeargs) =
                $(Symbol("unsafe_create$geom"))(IGeometry, coords::$typeargs)
        end)
    end
end

for typeargs in (Vector{Vector{Vector{Tuple{Cdouble,Cdouble}}}},
                 Vector{Vector{Vector{Tuple{Cdouble,Cdouble,Cdouble}}}},
                 Vector{Vector{Vector{Vector{Cdouble}}}})
    eval(quote
        function unsafe_createmultipolygon(
                ::Type{G},
                coords::$typeargs
            ) where G <: AbstractGeometry
            geom = unsafe_creategeom(G, GDAL.wkbMultiPolygon)
            for coord in coords
                addgeomdirectly!(geom, unsafe_createpolygon(coord))
            end
            geom
        end

        unsafe_createmultipolygon(coords::$typeargs) = unsafe_createmultipolygon(Geometry, coords)

        createmultipolygon(coords::$typeargs) = unsafe_createmultipolygon(IGeometry, coords)
    end)
end


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

mutable struct Provider <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace
    registrationState::Any # spec type: Union{ Nothing, String } # spec name: registrationState
    registrationPolicy::Any # spec type: Union{ Nothing, String } # spec name: registrationPolicy
    resourceTypes::Any # spec type: Union{ Nothing, Vector{ProviderResourceType} } # spec name: resourceTypes

    function Provider(;id=nothing, namespace=nothing, registrationState=nothing, registrationPolicy=nothing, resourceTypes=nothing)
        o = new()
        validate_property(Provider, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(Provider, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        validate_property(Provider, Symbol("registrationState"), registrationState)
        setfield!(o, Symbol("registrationState"), registrationState)
        validate_property(Provider, Symbol("registrationPolicy"), registrationPolicy)
        setfield!(o, Symbol("registrationPolicy"), registrationPolicy)
        validate_property(Provider, Symbol("resourceTypes"), resourceTypes)
        setfield!(o, Symbol("resourceTypes"), resourceTypes)
        o
    end
end # type Provider

const _property_map_Provider = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("namespace")=>Symbol("namespace"), Symbol("registrationState")=>Symbol("registrationState"), Symbol("registrationPolicy")=>Symbol("registrationPolicy"), Symbol("resourceTypes")=>Symbol("resourceTypes"))
const _property_types_Provider = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("namespace")=>"String", Symbol("registrationState")=>"String", Symbol("registrationPolicy")=>"String", Symbol("resourceTypes")=>"Vector{ProviderResourceType}")
Base.propertynames(::Type{ Provider }) = collect(keys(_property_map_Provider))
Swagger.property_type(::Type{ Provider }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Provider[name]))}
Swagger.field_name(::Type{ Provider }, property_name::Symbol) =  _property_map_Provider[property_name]

function check_required(o::Provider)
    true
end

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

ext: jl


text:
const ARB_CALC_SUCCESS = 0
const ARB_CALC_IMPRECISE_INPUT = 1
const ARB_CALC_NO_CONVERGENCE = 2

const BLOCK_NO_ZERO = 0
const BLOCK_ISOLATED_ZERO = 1
const BLOCK_UNKNOWN = 2

################################################################################
#
#  https://github.com/fredrik-johansson/arb/blob/master/arb_calc/refine_root_bisect.c
#
################################################################################

_arb_sign(t::arb_t) = begin
    arb_is_positive(t) && return 1
    ifelse(arb_is_negative(t), -1, 0)
end

_arf_interval_signs(func, t::arf_interval_t, arb::ArbField) = begin
    x = arb()

    arb_set_arf(x, t.a)
    y = func(x, 1)[1]
    asign = _arb_sign(y)

    arb_set_arf(x, t.b)
    y = func(x, 1)[1]
    bsign = _arb_sign(y)

    asign, bsign
end

arb_calc_partition(L::arf_interval_t, R::arf_interval_t, func, block::arf_interval_t, arb::ArbField) = begin
    # Compute the midpoint (TODO: try other points)
    u = arf()
    arf_add(u, block.a, block.b, ARF_PREC_EXACT, ARF_RND_DOWN)
    arf_mul_2exp_si(u, u, -1)

    # Evaluate and get sign at midpoint
    m = arb()
    arb_set_arf(m, u)
    t = func(m, 1)[1]
    msign = _arb_sign(t)

    # L, R = block, split at midpoint
    arf_set(L.a, block.a)
    arf_set(R.b, block.b)
    arf_set(L.b, u)
    arf_set(R.a, u)

    msign
end

arb_calc_refine_root_bisect(func, start::arf_interval_t, iter, arb::ArbField) = begin
    r = arf_interval_t()

    # must have proper sign changes
    asign, bsign = _arf_interval_signs(func, start, arb)
    (asign == 0 || bsign == 0 || asign == bsign) && return ARB_CALC_IMPRECISE_INPUT, r

    arf_interval_set(r, start)
    t, u = arf_interval(), arf_interval()

    for i in 1:iter
        msign = arb_calc_partition(t, u, func, r, arb)

        # the algorithm fails if the value at the midpoint cannot be distinguished from zero
        msign == 0 && return ARB_CALC_NO_CONVERGENCE, r

        arf_interval_swap(r, ifelse(msign == asign, u, t))
    end

    ARB_CALC_SUCCESS, r
end

arb_calc_refine_root_bisect{T}(func, start::T, iter, arb::ArbField) = begin
    interval = convert(arf_interval_t, start)
    s, r = arb_calc_refine_root_bisect(func, interval, iter, arb)
    s, convert(T, r)
end

################################################################################
#
#  https://github.com/fredrik-johansson/arb/blob/master/arb_calc/isolate_roots.c
#
################################################################################

check_block(func, block::arf_interval_t, asign, bsign, arb::ArbField) = begin
    x = arb()
    arf_interval_get_arb(x, block, prec(arb))

    t = func(x, 1)[1]

    !arb_is_positive(t) || return BLOCK_NO_ZERO
    !arb_is_negative(t) || return BLOCK_NO_ZERO
    asign < 0 < bsign || asign > 0 > bsign || return BLOCK_UNKNOWN

    t = func(x, 2)[2]

    arb_is_finite(t) || return BLOCK_UNKNOWN
    !arb_contains_zero(t) || return BLOCK_UNKNOWN

    BLOCK_ISOLATED_ZERO
end

isolate_roots_recursive(blocks::arf_interval_ptr, flags::Vector{Int}, func, block::arf_interval_t, asign, bsign, depth, eval_count, found_count, arb::ArbField, verbose) = begin

    add_block(block::arf_interval_t, status) = begin
        b = arf_interval()
        arf_interval_set(b, block)
        push!(blocks, b)
        push!(flags, status)
    end

    if eval_count <= 0 || found_count <= 0
        add_block(block, BLOCK_UNKNOWN)
        return eval_count, found_count
    end

    eval_count -= 1
    status = check_block(func, block, asign, bsign, arb)

    status != BLOCK_NO_ZERO || return eval_count, found_count

    if status == BLOCK_ISOLATED_ZERO || depth <= 0
        if status == BLOCK_ISOLATED_ZERO
            if verbose
                print("found isolated root in: ")
                arf_interval_printd(block, 15)
                println()
            end
            found_count -= 1
        end

        add_block(block, status)
    else
        L, R = arf_interval(), arf_interval()
        msign = arb_calc_partition(L, R, func, block, arb)

        if msign == 0 && verbose
            print("possible zero at midpoint: ")
            arf_interval_printd(block, 15)
            println()
        end

        eval_count, found_count = isolate_roots_recursive(blocks, flags, func, L, asign, msign, depth - 1, eval_count, found_count, arb, verbose)
        eval_count, found_count = isolate_roots_recursive(blocks, flags, func, R, msign, bsign, depth - 1, eval_count, found_count, arb, verbose)
    end

    eval_count, found_count
end

arb_calc_isolate_roots(func, block::arf_interval_t, maxdepth, maxeval, maxfound, arb::ArbField, verbose = false) = begin
    blocks = arf_interval_ptr()
    flags = Vector{Int}()

    asign, bsign = _arf_interval_signs(func, block, arb)

    isolate_roots_recursive(blocks, flags, func, block, asign, bsign, maxdepth, maxeval, maxfound, arb, verbose)

    blocks, flags
end

arb_calc_isolate_roots(func, block::T, maxdepth, maxeval, maxfound, arb::ArbField, verbose = false) where {T} = begin
    interval = convert(arf_interval_t, block)
    blocks, flags = arb_calc_isolate_roots(func, interval, maxdepth, maxeval, maxfound, arb, verbose)
    [convert(T, b) for b in blocks], flags
end


text:
using helloIBox
using Sockets
using HTTP
using Base64
using JSON
using Test

localIP =  Sockets.localhost
port = 8080
pathToIBox = "C:/Temp/IBox/IBoxLauncher.exe"
# pathToIBox = "C:/Yuly/IBox/IBoxLauncher.exe"

start_server(""; localIP = localIP, port = port)
r = HTTP.request("GET", "http://$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=8888&IBox_path=$pathToIBox&IBox_host=$localIP")
r = HTTP.request("GET", "http://$localIP:$port/api/getDataTree")
tree = JSON.parse(String(r.body))[1]

#поиск данных в дереве
info1 = helloIBox.findnode(tree, "EcgRecalcChannals") #запрос группы
@test info1["name"]=="EcgRecalcChannals" && length(info1["nodes"])==12
info2 = helloIBox.findnode(tree, "Noise_i_QRS") #запрос внутри группы более вложенное
@test info2["name"]=="Noise_i_QRS" && length(info2["nodes"])==0
info3 = helloIBox.findnode(tree, "Ecg_2")  #запрос внутри группы просто
@test info3["name"]=="Ecg_2" && info3["attrs"]["tag"] == "C1"

Freq = helloIBox.getData(localIP,8888,"Freq")[1] #читаем общую частоту-пока так
TimeSt = helloIBox.getData(localIP,8888,"StartTime")[1] #читаем общую частоту-пока так

attr = info3["attrs"]
attr["Freq"] = Freq
attr["TimeStart"] = TimeSt
#собираем временную сетку и лсб
TG, PhInfo = helloIBox.parseAttr(attr)
@test TG.fs == 257 && PhInfo.lsb.lsb == 1.063

#проверяем парсер масок
attr["Type.X"] = 1
attr["Type.S"] = 2
attr["Type.Z"] = 8
attr["datatype"] = eval(Symbol(info3["type"])) #добавляем в инфорамцию тип данных
mask = getMask(parseType(attr))
#PhysicalDataSet
DS1, dsreader,dsreader2,TG,PI,Mask = ds_new(localIP,8888,"EcgChannals","Ecg",attr)
@test isa(DS1, PhysicalDataSet) && DS1.data[1]("0","","2")==DS1.data[2]("0","","2")
@test DS1.data[1]("10","","2") == [-133, -144] && DS1.data[1]("10","12","") == [-133, -144]

#IndexDataSet
infoQPoint = helloIBox.findnode(tree, "QPoint")
infoQPoint["attrs"]["datatype"] = eval(Symbol(infoQPoint["type"]))
infoQPoint["attrs"]["Freq"] = Freq
infoQPoint["attrs"]["TimeStart"] = TimeSt

DSQ = helloIBox.ds_new(localIP,8888,"QRS","QPoint",infoQPoint["attrs"])
@test isa(DSQ, helloIBox.IndexDataSet) && DSQ.data[1]("0","","2")!=DSQ.data[2]("0","","2")
@test DSQ.data[1]("1","","2") == [221, 403] && DSQ.data[2]("400","500","") == [403] #сейчас не пройдет

#IntervalDataSet
infoWidthQRS = findnode(tree, "WidthQRS")
infoWidthQRS["attrs"]["datatype"] = eval(Symbol(infoWidthQRS["type"]))
DSW,dsreader,dsreader2,TG,PI,Mask = ds_new(localIP,8888,"QRS","WidthQRS",infoWidthQRS["attrs"])
@test isa(DSW, IntervalDataSet) && DSW.data[1]("0","","2")!=DSW.data[2]("0","","2")
@test DSW.data[1]("1","","2") == [0x0015, 0x0016] && DSW.data[2]("400","500","") == [0x0016] #но тут пока не ясно

# 1. правок нет - просто сквозные запросы без декодирования и создания объектов
# 2. пришла правка - делаем объект - вычитываем данные/декодируем/сводим/сохраняем - теперь всё чтение из сведенных результатов
infoQRS = helloIBox.findnode(tree, "QRS") #запрос группы

#данные в интервале
# r = HTTP.request("GET", "http://127.0.0.1:8888/api/getData?dataName=QPoint&res=oxy115829.dat&count=20&from=0&index=0")
# QPoint = reinterpret(Int32, base64decode(r.body)) |> collect
#
# r = HTTP.request("GET", "http://127.0.0.1:8888/api/getStructData?dataName=QRS&fields=QPoint&res=oxy115829.dat&count=20&from=400&index=0")
# QPoint = reinterpret(Int32, base64decode(r.body)) |> collect
#
# r = HTTP.request("GET", "http://127.0.0.1:8888/apibox/getCountDataInInterval?dataName=QPoint&index=0&count=20&from=400")
# cou = reinterpret(Int32, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getStructData?dataName=QRS&fields=QPoint&res=oxy115829.dat&count=20&from=400&index=0")
QPoint = reinterpret(Int32, base64decode(r.body)) |> collect
r = HTTP.request("GET", "http://127.0.0.1:8888/api/getStructData?dataName=QRS&fields=WidthQRS&res=oxy115829.dat&count=500&from=0&index=0")
QPoint = reinterpret(UInt16, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getData?dataName=WidthQRS&res=oxy115829.dat&count=6&from=0&index=0")
QPoint = reinterpret(UInt16, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getStructData?dataName=QRS&fields=NoiseQRS&res=oxy115829.dat&count=500&from=0&index=0")
QPoint = reinterpret(UInt8, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getData?dataName=NoiseQRS&res=oxy115829.dat&count=6&from=0&index=0")
QPoint = reinterpret(UInt8, base64decode(r.body)) |> collect

d = Dict{String,Any}()
d["ss"]=1


text:
function sum_rows(data, n, m)
    s = zero(eltype(data))
    for i in 1:n
        for j in 1:m
            s += data[i,j]
        end
    end
    return s
end

function sum_columns(data, n, m)
    s = zero(eltype(data))
    for j in 1:m
        for i in 1:n
            s += data[i,j]
        end
    end
    return s
end

using Random

Random.seed!(0)
n, m = 500, 10^6
data = rand(n,m)

#dry run for precompilation
sum_columns(data, n, m)
sum_rows(data, n, m)

println("Summing with the column iteration moving fastest")
@time sum_columns(data, n, m)
println("Summing with the row iteration moving fastest")
@time sum_rows(data, n, m);


text:
using Documenter, BenchlingIntegration

makedocs(;
    modules=[BenchlingIntegration],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/AcubeSAT/BenchlingIntegration.jl/blob/{commit}{path}#L{line}",
    sitename="BenchlingIntegration.jl",
    authors="Orestis Ousoultzoglou & Stavros Filosidis",
)

deploydocs(;
    repo="github.com/AcubeSAT/BenchlingIntegration.jl",
)


text:
include("./produce_data.jl")
using NPZ
using HDF5

function main(;
    )
    # create data
    data = produce_data_sampling_bias_illustration()
    npzwrite(@sprintf("%s/data/data_sampling_bias_illustration.npz", path_project),data)
end

# call function
main()


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    promote_rule(Binion{Int}, Float64) == Binion{Float64}
end

@test begin
    promote_rule(Binion{Int}, Binion{Float64}) == Binion{Float64}
end

@test begin
    -(1.0 - Binion(2, 3)) == Binion(1.0, 3.0) == Binion(2, 3) - 1.0
end

@test begin
    1.0 + Binion(2, 3) == Binion(3.0, 3.0) == Binion(2, 3) + 1.0
end

@test begin
    convert(Binion{Float64}, 1) == Binion(1.0, 0)
end

@test begin
    convert(Binion{Float64}, Binion(1,2)) == Binion(1.0, 2.0)
end

@test begin
    convert(Binion, 2) == Binion(2, 0)
end

@test begin
    z = Binion{Int}(1, 2)
    convert(Elliptic{Int}, z) == z
end

@test begin
    convert(Int, Binion(3)) == 3
end

@test_throws InexactError begin
    convert(Int, Binion(3,4))
end

@test begin
    length(unreal(1)) == 0
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Binion(a))
end

@test begin
    a = 1
    !isreal(Binion(a, a))
end

@test begin
    z = random(Binion{Int})
    z == +(z)
end

@test begin
    a = CayleyDickson.randomBigFloat()
    real(Binion(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, Binion(1, 2))
    l = String(take!(io))
    r = "[1: 1, i: 2]"
    l == r
end

@test begin
    x = random(Binion{BigFloat})
    y = random(Binion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Binion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Binion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    z = random(Binion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    z = random(Binion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Binion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Binion{Rational{BigInt}})
    y = random(Binion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Binion{Rational{BigInt}})
    y = random(Binion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Binion{Rational{BigInt}})
    y = random(Binion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    z = random(Binion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Binion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    z = random(Binion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    z = random(Binion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    z = random(Binion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Binion{BigInt})
    quadrance(x) > big(0)
end

@test begin
    x = random(Binion{BigInt})
    y = random(Binion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Binion(0, 0.0))
end

@test_throws ErrorException begin
    random(Binion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Binion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Binion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Binion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Binion{BigInt})
    z == z + zero(Binion{BigInt})
end

@test begin
    z = random(Binion{BigInt})
    z == z * one(Binion{BigInt})
end

@test begin
    z = random(Binion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Binion{BigInt})
    z == z * one(z)
end

@test begin
    x = random(Binion{Rational{BigInt}})
    selfstar(selfstar(x)) == selfstar(x)
end

@test begin
    x = random(Binion{Rational{BigInt}})
    anti_selfstar(anti_selfstar(x)) == anti_selfstar(x)
end


text:
using PackageCompiler

modulepth  = "compile/mkcell.jl"
testpth    = "test/runtests.jl"
builddir   = "_build"

build_executable(modulepth, snoopfile=testpth, builddir=builddir)


text:
using RecyclePermanent
using RNGPool
using JLD2

include("../models/gkModel.jl")
include("compare.jl")

@load "data/gkData.jld" ys θ
model = makeGKModel(θ, ys)

function findN(f::F, N0::Int64) where F<:Function
  c = 2
  N = c*N0
  output = f(N)
  while isnan(output) || output > 1.5
    c *= 2
    N = c*N0
    output = f(N)
  end
  d = 1.0
  prevOutput = output
  prevN = N
  while !isnan(output) && output < 1.5
    d *= 1.15
    prevOutput = output
    prevN = N
    N = ceil(Int64, c*N0/d)
    output = f(N)
  end
  return prevN, prevOutput
end

function makeRecycleN(model, n::Int64, trials::Int64)
  function recyclef(N::Int64)
    nt = Threads.nthreads()
    rpios = Vector{RPIO{model.Tx}}(undef, nt)
    println(N)
    for i in 1:nt
      rpios[i] = RPIO{model.Tx}(N, n, 1)
    end
    vals::Vector{Float64} = Vector{Float64}(undef, trials)
    Threads.@threads for trial in 1:trials
      vals[trial] = rpEstimate(model, rpios[Threads.threadid()])
    end
    logMean::Float64 = logmeanexp(vals)
    relativeVariance::Float64 = var(exp.(vals .- logMean))
    return relativeVariance
  end
  return recyclef
end

function makeSimpleN(model, n::Int64, trials::Int64)
  function simplef(N::Int64)
    nt = Threads.nthreads()
    rpios = Vector{RPIO{model.Tx}}(undef, nt)
    println(N)
    for i in 1:nt
      rpios[i] = RPIO{model.Tx}(N, n, 1)
    end
    vals::Vector{Float64} = Vector{Float64}(undef, trials)
    Threads.@threads for trial in 1:trials
      vals[trial] = simpleEstimate(model, rpios[Threads.threadid()])
    end
    logMean::Float64 = logmeanexp(vals)
    relativeVariance::Float64 = var(exp.(vals .- logMean))
    return relativeVariance
  end
  return simplef
end

function makeKMN(model, n::Int64, ϵ::Float64, trials::Int64)
  function kmf(N::Int64)
    nt = Threads.nthreads()
    kmios = Vector{KMIO{model.Tx}}(undef, nt)
    println(N)
    for i in 1:nt
      kmios[i] = KMIO{model.Tx}(N, n)
    end
    vals::Vector{Float64} = Vector{Float64}(undef, trials)
    Threads.@threads for trial in 1:trials
      vals[trial] = kmEstimate(model, kmios[Threads.threadid()], ϵ)
    end
    logMean::Float64 = logmeanexp(vals)
    relativeVariance::Float64 = var(exp.(vals .- logMean))
    return relativeVariance
  end
  return kmf
end

trials = 1000

setRNGs(12345)

ns = 50:50:250
outR = Vector{Tuple{Int64, Float64}}(undef, length(ns))
@time for i in 1:length(ns)
  output = findN(makeRecycleN(model, ns[i], trials), ns[i])
  outR[i] = output
end

outS = Vector{Tuple{Int64, Float64}}(undef, length(ns))
@time for i in 1:length(ns)
  output = findN(makeSimpleN(model, ns[i], trials), ns[i])
  outS[i] = output
end

outKM1 = Vector{Tuple{Int64, Float64}}(undef, length(ns))
@time for i in 1:length(ns)
  output = findN(makeKMN(model, ns[i], 1.0, trials), ns[i])
  outKM1[i] = output
end

outKM2 = Vector{Tuple{Int64, Float64}}(undef, length(ns))
@time for i in 1:length(ns)
  output = findN(makeKMN(model, ns[i], 0.001, trials), ns[i])
  outKM2[i] = output
end

outKM3 = Vector{Tuple{Int64, Float64}}(undef, length(ns))
@time for i in 1:length(ns)
  output = findN(makeKMN(model, ns[i], 0.0001, trials), ns[i])
  outKM3[i] = output
end

@save "compare/gktau.jld" ns trials outR outS outKM1 outKM2 outKM3


text:
# Microscopic water models.

"""
    WaterModel

Potentials and forces describing water molecules and their interactions.
"""
abstract type WaterModel end

"""
    TermContainer

Container for potential and force terms.
"""
struct TermContainer
    "Intramolecular terms."
    terms_intra::Vector{Function}
    "Intermolecular terms."
    terms_inter::Vector{Function}
end

function build_model(model_T::Type{<:WaterModel},
                     potential_intra::Function, potential_inter::Function,
                     force_intra::Function, force_inter::Function,
                     N::Int)
    potential_intra_terms = []
    potential_inter_terms = []
    force_intra_terms = []
    force_inter_terms = []

    for n1 in 1:N
        push!(potential_intra_terms, potential_intra(n1))
        push!(force_intra_terms, force_intra(n1))

        for n2 in (n1+1):N
            push!(potential_inter_terms, potential_inter(n1, n2))
            push!(force_inter_terms, force_inter(n1, n2))
        end
    end

    potential = TermContainer(potential_intra_terms, potential_inter_terms)
    force = TermContainer(force_intra_terms, force_inter_terms)

    model_T(potential, force)
end

"""
    ForceContainer

Container for evaluated forces. Each value is the total force on that degree of
freedom from all terms.
"""
mutable struct ForceContainer
    "Forces (kJ / mol nm)."
    f::Array{Float64,4}
    "Whether the forces have been set."
    initialized::Bool
end

function ForceContainer(P::Int, D::Int, A::Int, N::Int)
    f = Array{Float64}(undef, P, D, A, N)
    initialized = false

    ForceContainer(f, initialized)
end

"""
    eval_pot(model::WaterModel, s::State)

Evaluate the `model` potential for `s`.
"""
function eval_pot(model::WaterModel, s::State{P,D,A,N}) where {P,D,A,N}
    result = 0.0 # kJ / mol

    for term_f in model.potential.terms_intra
        for j in 1:P
            result += term_f(j, s)
        end
    end

    for term_f in model.potential.terms_inter
        for j in 1:P
            result += term_f(j, s)
        end
    end

    result
end

"""
    eval_force!(result::ForceContainer, model::WaterModel, s::State)

Evaluate the `model` force for `s`, storing the result in `result`.
"""
function eval_force!(result::ForceContainer, model::WaterModel,
                     s::State{P,D,A,N}) where {P,D,A,N}
    for term_f! in model.force.terms_intra
        for j in 1:P
            term_f!(result, j, s)
        end
    end

    for term_f! in model.force.terms_inter
        for j in 1:P
            term_f!(result, j, s)
        end
    end

    nothing
end

const registered_models = Dict{String,Type{<:WaterModel}}()

list_models() = sort(collect(keys(registered_models)))

get_model(name::String) = registered_models[name]

include("model_harmonic.jl")
include("model_qspcfw.jl")


text:
# deprecations should also include a comment stating the version
# for which the deprecation should be removed.

# Deprecated to fix spelling. Can be removed for version 2.0.
@deprecate simplecycles_hadwick_james simplecycles_hawick_james

# Deprecated for more explicit function name. Can be removed for version 2.0.
@deprecate saw self_avoiding_walk


text:
using LoopVectorization

function oneHot(input::Array{Int64}, depth::Int64, dict::Dict; on_value::Float32=1.0f0, off_value::Float32=0.0f0)
    output_data = fill(off_value, (depth, size(input, 1)))
    @avxt for i in 1:length(input)
        output_data[dict[input[i]],i] = on_value
    end
    return output_data
end


text:
#=
The sum of the squares of the first ten natural numbers is,

12 + 22 + ... + 102 = 385
The square of the sum of the first ten natural numbers is,

(1 + 2 + ... + 10)2 = 552 = 3025
Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 − 385 = 2640.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.
=#

function sumsquaredifferences()
  collection = 1:100
  x = sum(map((x)->x^2, collection))
  y = sum(collection)^2
  y - x
end
@time println(sumsquaredifferences())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module WhiteBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "White Block" => Ahorn.EntityPlacement(
        Maple.WhiteBlock
    )
)

sprite = "objects/whiteblock.png"

function Ahorn.selection(entity::Maple.WhiteBlock)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.0, jy=0.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.WhiteBlock, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.0, jy=0.0)

end


text:
using vfld
using Test

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


text:
#!/usr/local/bin/julia -q --color=yes

dataFile = Pkg.dir(homedir(), ".jonepw.jls")

save(data) = begin
  open(f->serialize(f, data), dataFile, "w")
end

load() = begin
  open(deserialize, dataFile)
end

type Onepw
  id::AbstractString
  tag::AbstractString
  pwd::AbstractString
end

if !isfile(dataFile)
  onepw = Onepw[]
  save(onepw)
end

onepw = load()

align(str, len) = begin
  spaces = repeat(" ", len - length(str))
  return "$str$spaces"
end

ls(args) = begin
  for item in onepw
    println("$(align(item.id, 16))$(align(item.tag, 16))$(align(item.pwd, 32))")
  end
end

rm(args) = begin
  for (i, item) in enumerate(onepw)
    if args[1] == item.id
      splice!(onepw, i)
    end
  end
  ls([])
  save(onepw)
end

add(args) = begin
  id = num2hex(time())[9:end]
  push!(onepw, Onepw(id, args[1], args[2]))
  ls([])
  save(onepw)
end

gen(args) = begin
  println(randstring(parse(Int, args[1])))
end

usage(args) = begin
  str = """
    ls  => list
    rm  => remove
    add => add
    gen => generate
  """

  println(str)
end

handler = Dict(
  "ls"  => ls,
  "rm"  => rm,
  "add" => add,
  "gen" => gen
)

main() = begin
  print("password: ")
  pw = readline(STDIN)

  if length(ARGS) == 0
    ls([])
  else
    get(handler, ARGS[1], usage)(ARGS[2:end])
  end
end

main()

max_stars_repo_name: moldray/juliaOnePW

max_issues_repo_name: moldray/juliaOnePW

max_forks_repo_name: moldray/juliaOnePW


text:
dimname{N}(a::Axis{N}) = N

timestamps(A::TimeAxisArray) = A.axes[1].val
interval(A::RegularTimeAxisArray) = A |> timestamps |> step

function isregularspaced(x::AbstractTimeVector)
    spacing = x[2:end] - x[1:end-1]
    i = 1
    regularspaced = true
    while regularspaced && i < length(spacing)
        regularspaced = spacing[i] == spacing[i+1]
        i += 1
    end #while
    return regularspaced
end #isregularspaced

stringfallback(x, fallback) = x
stringfallback(x::AbstractString, fallback) = fallback

date(format::AbstractString="yyyy-mm-dd"; locale="english") = x -> Date(x, format, locale=locale)
datetime(format::AbstractString="yyyy-mm-ddTHH:MM:SS.sss"; locale="english") = x -> DateTime(x, format, locale=locale)

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
using Main.AffineM: affinebrute

while true

    # Get info in
    println("\nText [leave blank to break]:")
    human_input = readline()
    human_input == "" ? break :
    println("\nThe brute only works if the text is English so we assume that's what you're using.")
    human_alphabet = "abcdefghijklmnopqrstuvwxyz"
    println()

    println("Here's the best we got: \n", affinebrute(human_input, human_alphabet))

end

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
# ------------------------------------------------------------------------------------------
# Models
# ------------------------------------------------------------------------------------------
struct TopPushK{S<:Surrogate, T<:Real} <: AbstractTopPush{S}
    K::Int
    C::T
    l::S
    state::State{T}

    function TopPushK(
        K::Int;
        C::Real = 1,
        ϑ::Real = 1,
        surrogate::Type{<:Surrogate} = Hinge,
        T = Float32,
    )
        l = surrogate(T(ϑ))
        return new{typeof(l), T}(K, T(C), l, State(T))
    end
end

function parameters(m::TopPushK{S}) where {S<:Surrogate}
    return (K = m.K, C = m.C, surrogate = S.name.name, ϑ = m.l.ϑ)
end
model_K(model::TopPushK, ::KernelMatrix) = model.K

struct τFPL{S<:Surrogate, T<:Real} <: AbstractTopPush{S}
    τ::T
    C::T
    l::S
    state::State{T}

    function τFPL(
        τ::Real;
        C::Real = 1,
        ϑ::Real = 1,
        surrogate::Type{<:Surrogate} = Hinge,
        T = Float32,
    )
        l = surrogate(T(ϑ))
        return new{typeof(l), T}(T(τ), T(C), l, State(T))
    end
end

function parameters(m::τFPL{S}) where {S<:Surrogate}
    return (τ = m.τ, C = m.C, surrogate = S.name.name, ϑ = m.l.ϑ)
end
model_K(model::τFPL, K::KernelMatrix) = max(1, round(Int, model.τ * K.nβ))

function threshold(model::AbstractTopPush, K::KernelMatrix)
    return - mean(partialsort(model.state.s[inds_β(K)], 1:model_K(model, K)))
end

# TopPush
struct TopPush{S<:Surrogate, T<:Real} <: AbstractTopPush{S}
    C::T
    l::S
    state::State{T}

    function TopPush(;
        C::Real = 1,
        ϑ::Real = 1,
        surrogate::Type{<:Surrogate} = Hinge,
        T = Float32,
    )
        l = surrogate(T(ϑ))
        return new{typeof(l), T}(T(C), l, State(T))
    end
end

function parameters(m::TopPush{S}) where {S<:Surrogate}
    return (C = m.C, surrogate = S.name.name, ϑ = m.l.ϑ)
end
threshold(model::TopPush, K::KernelMatrix) = - minimum(model.state.s[inds_β(K)])

# ------------------------------------------------------------------------------------------
# Auxilliary functions
# ------------------------------------------------------------------------------------------
function initialization!(model::TopPush, K::KernelMatrix)
    T = eltype(K)

    model.state.s = zeros(T, K.n)
    model.state.αβ = zeros(T, K.n)
    return
end

function initialization!(model::AbstractTopPush{S}, K::KernelMatrix) where {S <: Surrogate}
    T = eltype(K)
    α, β, Km = rand(T, K.nα), zeros(T, K.nβ), model_K(model, K)
    if S <: Hinge
        α, β = projection(α, β, model.l.ϑ*model.C, Km)
    elseif S <: Quadratic
        α, β = projection(α, β, Km)
    else
        throw(ArgumentError("$(S.name.name) surrogate function is not supported"))
    end
    αβ = vcat(α, β)

    model.state.s = K * αβ
    model.state.αβ = αβ
    model.state.αsum = sum(αβ[inds_α(K)])
    model.state.βsort = sort(αβ[inds_β(K)], rev = true)
    return
end

function permutation(::AbstractTopPush, y::BitVector)
    perm_α = findall(y)
    perm_β = findall(.~y)

    return length(perm_α), length(perm_β), vcat(perm_α, perm_β)
end

function extract_scores(model::AbstractTopPush, K::KernelMatrix)
    s = copy(model.state.s)
    s[inds_β(K)] .*= -1
    return s[invperm(K.perm)]
end

function extract_state(model::AbstractTopPush, K::KernelMatrix)
    αβ = copy(model.state.αβ)
    return  (α = αβ[inds_α(K)], β = αβ[inds_β(K)])
end

function objective(model::AbstractTopPush{S}, K::KernelMatrix) where {S<:Surrogate}
    @unpack s, αβ = model.state
    ϑ, C = model.l.ϑ, model.C

    α = αβ[inds_α(K)]
    sα = s[inds_α(K)]
    w_norm = s'*αβ/2
    t = threshold(model, K)

    # objectives
    L_primal = w_norm + C*sum(value.(model.l, t .- sα))
    if S <: Hinge
        L_dual = - w_norm + sum(α)/ϑ
    else
        L_dual = - w_norm + sum(α)/ϑ - sum(abs2, α)/(4*ϑ*C)
    end
    return L_primal, L_dual, L_primal - L_dual
end

# ------------------------------------------------------------------------------------------
# Update rules - Hinge loss
# ------------------------------------------------------------------------------------------
function rule_αα(model::AbstractTopPush{<:Hinge}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ = model.state
    C, ϑ = model.C, model.l.ϑ

    num = s[k] - s[l]
    den = K[k, k] - 2*K[k, l] + K[l, l]
    lb = max(- αβ[k], αβ[l] - ϑ*C)
    ub = min(ϑ*C - αβ[k], αβ[l])

    return UpdateRule(num, den, lb, ub, k, l)
end

function rule_αβ(model::AbstractTopPush{<:Hinge}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, αsum = model.state
    C, ϑ = model.C, model.l.ϑ

    num = s[k] + s[l] - 1/ϑ
    den = K[k, k] + 2*K[k, l] + K[l, l]
    if isa(model, TopPush)
        lb = max(- αβ[k], - αβ[l])
        ub = ϑ*C - αβ[k]
    else
        βmax = find_βmax(model, αβ[l])
        Km = model_K(model, K)
        lb = max(- αβ[k], - αβ[l], Km*βmax - αsum)
        ub = min(ϑ*C - αβ[k], (αsum - Km*αβ[l])/(Km - 1))
    end
    return UpdateRule(num, den, lb, ub, k, l)
end

function rule_ββ(model::AbstractTopPush{<:Hinge}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, αsum = model.state

    num = s[k] - s[l]
    den = K[k, k] - 2*K[k, l] + K[l, l]
    if isa(model, TopPush)
        lb = - αβ[k]
        ub = αβ[l]
    else
        Km = model_K(model, K)
        lb = max(- αβ[k], αβ[l] - αsum/Km)
        ub = min(αsum/Km - αβ[k], αβ[l])
    end
    return UpdateRule(num, den, lb, ub, k, l)
end

# ------------------------------------------------------------------------------------------
# Update rules - Quadratic loss
# ------------------------------------------------------------------------------------------
function rule_αα(model::AbstractTopPush{<:Quadratic}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ = model.state
    C, ϑ = model.C, model.l.ϑ

    num = s[k] - s[l] + (αβ[k] - αβ[l])/(2*C*ϑ^2)
    den = K[k, k] - 2*K[k, l] + K[l, l] + 1/(C*ϑ^2)
    lb = - αβ[k]
    ub = αβ[l]

    return UpdateRule(num, den, lb, ub, k, l)
end

function rule_αβ(model::AbstractTopPush{<:Quadratic}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, αsum = model.state
    C, ϑ = model.C, model.l.ϑ

    num = s[k] + s[l] - 1/ϑ + αβ[k]/(2*C*ϑ^2)
    den = K[k, k] + 2*K[k, l] + K[l, l] + 1/(2*C*ϑ^2)
    if isa(model, TopPush)
        lb = max(- αβ[k], - αβ[l])
        ub = eltype(s)(Inf)
    else
        βmax = find_βmax(model, αβ[l])
        Km = model_K(model, K)
        lb = max(- αβ[k], - αβ[l], Km*βmax - αsum)
        ub = (αsum - Km*αβ[l])/(Km - 1)
    end
    return UpdateRule(num, den, lb, ub, k, l)
end

function rule_ββ(model::AbstractTopPush{<:Quadratic}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, αsum = model.state

    num = s[k] - s[l]
    den = K[k, k] - 2*K[k, l] + K[l, l]
    if isa(model, TopPush)
        lb = - αβ[k]
        ub = αβ[l]
    else
        Km = model_K(model, K)
        lb = max(- αβ[k], αβ[l] - αsum/Km)
        ub = min(αsum/Km - αβ[k], αβ[l])
    end
    return UpdateRule(num, den, lb, ub, k, l)
end


text:
function diffeq(du,u,p,t)
    Rb    = u[V.uRb]  - u[V.tE2f]  + u[V.E2f]
    pRb   = p[C.tRb]  - u[V.uRb]
    RbE2f = u[V.tE2f] - u[V.E2f]

    uERK = p[C.tERK] - u[V.pERK] - u[V.pERK_DUSP]

    Ce  = u[V.tCe]  - u[V.CeP21]

    P21 = p[C.tP21]-u[V.CeP21]
    ###
    rSyE2f = p[C.kSyE2f] + p[C.kSyE2fE2f]*u[V.E2f]/(p[C.jSyE2f] + u[V.E2f])
    rPhRb  = p[C.kPhRbCd]*u[V.Cd] + p[C.kPhRbCe]*Ce
    rDeP21 = p[C.kDeP21] + p[C.kDeP21Cy]*p[C.Skp2]*Ce + p[C.pAkt]

    ###
    du[V.pERK] = (p[C.kf]*p[C.mitogen]*uERK - p[C.kr]*u[V.pERK]) - (p[C.kon]*p[C.DUSP]*u[V.pERK] - p[C.koff]*u[V.pERK_DUSP])
    du[V.pERK_DUSP] = p[C.kon]*p[C.DUSP]*u[V.pERK] - p[C.koff]*u[V.pERK_DUSP]
    du[V.pRSK] = p[C.kfRSK]*u[V.pERK]/(p[C.KmRSK]+u[V.pERK]) - p[C.krRSK]*u[V.pRSK]
    du[V.cFos] = p[C.kfFOS]*u[V.pERK]*u[V.pRSK]/(p[C.KmFOS]+u[V.pRSK]) - p[C.krFOS]*u[V.cFos]
    du[V.Cd] = p[C.kSyCd]*u[V.cFos] - p[C.kDeCd]*u[V.Cd]
    du[V.uRb]  = - rPhRb*u[V.uRb] + p[C.kDpRb]*pRb
    du[V.tE2f] = rSyE2f - p[C.kDeE2f]*u[V.tE2f]
    du[V.E2f]  = rSyE2f - p[C.kDeE2f]*u[V.E2f] - p[C.kAsRbE2f]*Rb*u[V.E2f] + (p[C.kDsRbE2f]+rPhRb)*RbE2f
    du[V.tCe]  = p[C.kSyCe] *u[V.E2f] - p[C.kDeCe]*u[V.tCe]
    du[V.CeP21] = p[C.kAsCyP21]*P21*Ce - (p[C.kDsCyP21]+p[C.kDeCe]+rDeP21)*u[V.CeP21]

end

function get_steady_state(p::Vector{Float64})
    tspan_ss::Tuple{Float64,Float64} = (0.0,Inf)
    u0::Vector{Float64} = zeros(V.NUM)

    prob = ODEProblem(diffeq,u0,tspan_ss,p)
    prob = SteadyStateProblem(prob)
    sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)

    return sol.u
end