text:
"""
    fromWKB(data)

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

### Parameters
* `data`: pointer to the input BLOB data.
"""
function fromWKB(data)
    geom = Ref{GDALGeometry}()
    result = @gdal(OGR_G_CreateFromWkb::GDAL.OGRErr,
        data::Ptr{Cuchar},
        C_NULL::GDALSpatialRef,
        geom::Ptr{GDALGeometry},
        sizeof(data)::Cint
    )
    @ogrerr result "Failed to create geometry from WKB"
    return IGeometry(geom[])
end

function unsafe_fromWKB(data)
    geom = Ref{GDALGeometry}()
    result = @gdal(OGR_G_CreateFromWkb::GDAL.OGRErr,
        data::Ptr{Cuchar},
        C_NULL::GDALSpatialRef,
        geom::Ptr{GDALGeometry},
        sizeof(data)::Cint
    )
    @ogrerr result "Failed to create geometry from WKB"
    return Geometry(geom[])
end

"""
    fromWKT(data::Vector{String})

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.
"""
function fromWKT(data::Vector{String})
    geom = Ref{GDALGeometry}()
    result = @gdal(OGR_G_CreateFromWkt::GDAL.OGRErr,
        data::StringList,
        C_NULL::GDALSpatialRef,
        geom::Ptr{GDALGeometry}
    )
    @ogrerr result "Failed to create geometry from WKT"
    return IGeometry(geom[])
end

function unsafe_fromWKT(data::Vector{String})
    geom = Ref{GDALGeometry}()
    result = @gdal(OGR_G_CreateFromWkt::GDAL.OGRErr,
        data::StringList,
        C_NULL::GDALSpatialRef,
        geom::Ptr{GDALGeometry}
    )
    @ogrerr result "Failed to create geometry from WKT"
    return Geometry(geom[])
end

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

unsafe_fromWKT(data::String, args...) = unsafe_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.ogr_g_destroygeometry(geom.ptr)
    geom.ptr = C_NULL
end

"""
    clone(geom::AbstractGeometry)

Returns a copy of the geometry with the original spatial reference system.
"""
function clone(geom::AbstractGeometry)
    if geom.ptr == C_NULL
        return IGeometry()
    else
        return IGeometry(GDAL.ogr_g_clone(geom.ptr))
    end
end

function unsafe_clone(geom::AbstractGeometry)
    if geom.ptr == C_NULL
        return Geometry()
    else
        return Geometry(GDAL.ogr_g_clone(geom.ptr))
    end
end

"""
    creategeom(geomtype::OGRwkbGeometryType)

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.
"""
creategeom(geomtype::OGRwkbGeometryType) =
    IGeometry(GDAL.ogr_g_creategeometry(geomtype))

unsafe_creategeom(geomtype::OGRwkbGeometryType) =
    Geometry(GDAL.ogr_g_creategeometry(geomtype))

"""
    forceto(geom::AbstractGeometry, targettype::OGRwkbGeometryType, [options])

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 forceto(
        geom::AbstractGeometry,
        targettype::OGRwkbGeometryType,
        options = StringList(C_NULL)
    )
    return IGeometry(GDAL.ogr_g_forceto(unsafe_clone(geom).ptr, targettype,
        options))
end

function unsafe_forceto(
        geom::AbstractGeometry,
        targettype::OGRwkbGeometryType,
        options = StringList(C_NULL)
    )
    return Geometry(GDAL.ogr_g_forceto(unsafe_clone(geom).ptr, targettype,
        options))
end

"""
    geomdim(geom::AbstractGeometry)

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).
"""
geomdim(geom::AbstractGeometry) = GDAL.ogr_g_getdimension(geom.ptr)

"""
    getcoorddim(geom::AbstractGeometry)

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.ogr_g_getcoordinatedimension(geom.ptr)

"""
    setcoorddim!(geom::AbstractGeometry, dim::Integer)

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.ogr_g_setcoordinatedimension(geom.ptr ,dim)
    return geom
end

"""
    envelope(geom::AbstractGeometry)

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

"""
    envelope3d(geom::AbstractGeometry)

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

"""
    boundingbox(geom::AbstractGeometry)

Returns a bounding box polygon (CW) built from envelope coordinates
"""
function boundingbox(geom::AbstractGeometry)
    coordinates = envelope(geom)
    MinX, MaxX, MinY, MaxY = coordinates.MinX, coordinates.MaxX, coordinates.MinY, coordinates.MaxY
    # creates a CW closed ring polygon
    return createpolygon([[MaxY, MinX], [MaxY, MaxX], [MinY, MaxX], [MinY, MinX], [MaxY, MinX]])
end

"""
    toWKB(geom::AbstractGeometry, order::OGRwkbByteOrder = GDAL.wkbNDR)

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}(undef, wkbsize(geom))
    result = GDAL.ogr_g_exporttowkb(geom.ptr, order, buffer)
    @ogrerr result "Failed to export geometry to WKB"
    return buffer
end

"""
    toISOWKB(geom::AbstractGeometry, order::OGRwkbByteOrder = GDAL.wkbNDR)

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}(undef, wkbsize(geom))
    result = GDAL.ogr_g_exporttoisowkb(geom.ptr, order, buffer)
    @ogrerr result "Failed to export geometry to ISO WKB"
    return buffer
end

"""
    wkbsize(geom::AbstractGeometry)

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

"""
    toWKT(geom::AbstractGeometry)

Convert a geometry into well known text format.
"""
function toWKT(geom::AbstractGeometry)
    wkt_ptr = Ref(Cstring(C_NULL))
    result = GDAL.ogr_g_exporttowkt(geom.ptr, wkt_ptr)
    @ogrerr result "OGRErr $result: failed to export geometry to WKT"
    wkt = unsafe_string(wkt_ptr[])
    GDAL.vsifree(pointer(wkt_ptr[]))
    return wkt
end

"""
    toISOWKT(geom::AbstractGeometry)

Convert a geometry into SFSQL 1.2 / ISO SQL/MM Part 3 well known text format.
"""
function toISOWKT(geom::AbstractGeometry)
    isowkt_ptr = Ref(Cstring(C_NULL))
    result = GDAL.ogr_g_exporttoisowkt(geom.ptr, isowkt_ptr)
    @ogrerr result "OGRErr $result: failed to export geometry to ISOWKT"
    wkt = unsafe_string(isowkt_ptr[])
    GDAL.vsifree(pointer(isowkt_ptr[]))
    return wkt
end

"""
    getgeomtype(geom::AbstractGeometry)

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

"""
    geomname(geom::AbstractGeometry)

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

"""
    flattento2d!(geom::AbstractGeometry)

Convert geometry to strictly 2D.
"""
function flattento2d!(geom::AbstractGeometry)
    GDAL.ogr_g_flattento2d(geom.ptr)
    return geom
end

"""
    closerings!(geom::AbstractGeometry)

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.
"""
function closerings!(geom::AbstractGeometry)
    GDAL.ogr_g_closerings(geom.ptr)
    return geom
end

"""
    fromGML(data)

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.
"""
fromGML(data) = IGeometry(GDAL.ogr_g_createfromgml(data))

unsafe_fromGML(data) = Geometry(GDAL.ogr_g_createfromgml(data))

"""
    toGML(geom::AbstractGeometry)

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

"""
    toKML(geom::AbstractGeometry, altitudemode = C_NULL)

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

"""
    toJSON(geom::AbstractGeometry)

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

"""
    toJSON(geom::AbstractGeometry, options)

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) =
    GDAL.ogr_g_exporttojsonex(geom.ptr, options)

"""
    fromJSON(data::String)

Create a geometry object from its GeoJSON representation.
"""
fromJSON(data::String) = IGeometry(GDAL.ogr_g_creategeometryfromjson(data))

unsafe_fromJSON(data::String) =
    Geometry(GDAL.ogr_g_creategeometryfromjson(data))

# """
# Assign spatial reference to this object.

# Any existing spatial reference is replaced, but under no circumstances
# does this result in the object being reprojected. It is just changing
# the interpretation of the existing geometry. Note that assigning a
# spatial reference increments the reference count on the
# OGRSpatialReference, but does not copy it.

# Starting with GDAL 2.3, this will also assign the spatial reference to
# potential sub-geometries of the geometry (OGRGeometryCollection,
# OGRCurvePolygon/OGRPolygon, OGRCompoundCurve, OGRPolyhedralSurface and
# their derived classes).
# """
# function setspatialref!(geom::Geometry, spatialref::AbstractSpatialRef)
#     GDAL.assignspatialreference(geom.ptr, spatialref.ptr)
#     return geom
# end

"""
    getspatialref(geom::AbstractGeometry)

Returns a clone of the spatial reference system for the geometry.

(The original SRS may be shared with many objects, and should not be modified.)
"""
function getspatialref(geom::AbstractGeometry)
    if geom.ptr == C_NULL
        return ISpatialRef()
    end
    result = GDAL.ogr_g_getspatialreference(geom.ptr)
    if result == C_NULL
        return ISpatialRef()
    else
        return ISpatialRef(GDAL.osrclone(result))
    end
end

function unsafe_getspatialref(geom::AbstractGeometry)
    if geom.ptr == C_NULL
        return SpatialRef()
    end
    result = GDAL.ogr_g_getspatialreference(geom.ptr)
    if result == C_NULL
        return SpatialRef()
    else
        return SpatialRef(GDAL.osrclone(result))
    end
end

"""
    transform!(geom::AbstractGeometry, coordtransform::CoordTransform)

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::AbstractGeometry, coordtransform::CoordTransform)
    result = GDAL.ogr_g_transform(geom.ptr, coordtransform.ptr)
    @ogrerr result "Failed to transform geometry"
    return geom
end

# """
# Transform geometry to new spatial reference system.

# This function will transform the coordinates of a geometry from their
# current spatial reference system to a new target spatial reference
# system. Normally this means reprojecting the vectors, but it could
# include datum shifts, and changes of units.

# This function will only work if the geometry already has an assigned
# spatial reference system, and if it is transformable to the target
# coordinate system.

# Because this function requires internal creation and initialization of
# an OGRCoordinateTransformation object it is significantly more
# expensive to use this function to transform many geometries than it is
# to create the OGRCoordinateTransformation in advance, and call
# transform() with that transformation. This function exists primarily
# for convenience when only transforming a single geometry.

# ### Parameters
# * `geom`: handle on the geometry to apply the transformation.
# * `spatialref`: Target spatial reference system.
# """
# function transform!(geom::AbstractGeometry, spatialref::AbstractSpatialRef)
#     result = GDAL.ogr_g_transformto(geom.ptr, spatialref.ptr)
#     @ogrerr result "Failed to transform geometry to the new SRS"
#     return geom
# end

"""
    simplify(geom::AbstractGeometry, tol::Real)

Compute a simplified geometry.

### Parameters
* `geom`: the geometry.
* `tol`: the distance tolerance for the simplification.
"""
simplify(geom::AbstractGeometry, tol::Real) =
    IGeometry(GDAL.ogr_g_simplify(geom.ptr, tol))

unsafe_simplify(geom::AbstractGeometry, tol::Real) =
    Geometry(GDAL.ogr_g_simplify(geom.ptr, tol))

"""
    simplifypreservetopology(geom::AbstractGeometry, tol::Real)

Simplify the geometry while preserving topology.

### Parameters
* `geom`: the geometry.
* `tol`: the distance tolerance for the simplification.
"""
simplifypreservetopology(geom::AbstractGeometry, tol::Real) =
    IGeometry(GDAL.ogr_g_simplifypreservetopology(geom.ptr, tol))

unsafe_simplifypreservetopology(geom::AbstractGeometry, tol::Real) =
    Geometry(GDAL.ogr_g_simplifypreservetopology(geom.ptr, tol))

"""
    delaunaytriangulation(geom::AbstractGeometry, tol::Real, onlyedges::Bool)

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.
"""
delaunaytriangulation(geom::AbstractGeometry, tol::Real, onlyedges::Bool) =
    IGeometry(GDAL.ogr_g_delaunaytriangulation(geom.ptr, tol, onlyedges))

function unsafe_delaunaytriangulation(
        geom::AbstractGeometry,
        tol::Real,
        onlyedges::Bool
    )
    return Geometry(GDAL.ogr_g_delaunaytriangulation(geom.ptr, tol, onlyedges))
end

"""
    segmentize!(geom::AbstractGeometry, maxlength::Real)

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::AbstractGeometry, maxlength::Real)
    GDAL.ogr_g_segmentize(geom.ptr, maxlength)
    return geom
end

"""
    intersects(g1::AbstractGeometry, g2::AbstractGeometry)

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.
"""
intersects(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_intersects(g1.ptr, g2.ptr))

"""
    equals(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if the geometries are equivalent.
"""
equals(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_equals(g1.ptr, g2.ptr))

"""
    disjoint(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if the geometries are disjoint.
"""
disjoint(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_disjoint(g1.ptr, g2.ptr))

"""
    touches(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if the geometries are touching.
"""
touches(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_touches(g1.ptr, g2.ptr))

"""
    crosses(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if the geometries are crossing.
"""
crosses(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_crosses(g1.ptr, g2.ptr))

"""
    within(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if g1 is contained within g2.
"""
within(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_within(g1.ptr, g2.ptr))

"""
    contains(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if g1 contains g2.
"""
contains(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_contains(g1.ptr, g2.ptr))

"""
    overlaps(g1::AbstractGeometry, g2::AbstractGeometry)

Returns `true` if the geometries overlap.
"""
overlaps(g1::AbstractGeometry, g2::AbstractGeometry) =
    Bool(GDAL.ogr_g_overlaps(g1.ptr, g2.ptr))

"""
    boundary(geom::AbstractGeometry)

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.
"""
boundary(geom::AbstractGeometry) = IGeometry(GDAL.ogr_g_boundary(geom.ptr))

unsafe_boundary(geom::AbstractGeometry) =
    Geometry(GDAL.ogr_g_boundary(geom.ptr))

"""
    convexhull(geom::AbstractGeometry)

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.
"""
convexhull(geom::AbstractGeometry) = IGeometry(GDAL.ogr_g_convexhull(geom.ptr))

unsafe_convexhull(geom::AbstractGeometry) =
    Geometry(GDAL.ogr_g_convexhull(geom.ptr))

"""
    buffer(geom::AbstractGeometry, dist::Real, quadsegs::Integer = 30)

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.
"""
buffer(geom::AbstractGeometry, dist::Real, quadsegs::Integer = 30) =
    IGeometry(GDAL.ogr_g_buffer(geom.ptr, dist, quadsegs))

unsafe_buffer(geom::AbstractGeometry, dist::Real, quadsegs::Integer = 30) =
    Geometry(GDAL.ogr_g_buffer(geom.ptr, dist, quadsegs))

"""
    intersection(g1::AbstractGeometry, g2::AbstractGeometry)

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.
"""
intersection(g1::AbstractGeometry, g2::AbstractGeometry) =
    IGeometry(GDAL.ogr_g_intersection(g1.ptr, g2.ptr))

unsafe_intersection(g1::AbstractGeometry, g2::AbstractGeometry) =
    Geometry(GDAL.ogr_g_intersection(g1.ptr, g2.ptr))

"""
    union(g1::AbstractGeometry, g2::AbstractGeometry)

Returns a new geometry representing the union of the geometries.
"""
union(g1::AbstractGeometry, g2::AbstractGeometry) =
    IGeometry(GDAL.ogr_g_union(g1.ptr, g2.ptr))

unsafe_union(g1::AbstractGeometry, g2::AbstractGeometry) =
    Geometry(GDAL.ogr_g_union(g1.ptr, g2.ptr))

"""
    pointonsurface(geom::AbstractGeometry)

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).
"""
pointonsurface(geom::AbstractGeometry) =
    IGeometry(GDAL.ogr_g_pointonsurface(geom.ptr))

unsafe_pointonsurface(geom::AbstractGeometry) =
    Geometry(GDAL.ogr_g_pointonsurface(geom.ptr))

"""
    difference(g1::AbstractGeometry, g2::AbstractGeometry)

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.
"""
difference(g1::AbstractGeometry, g2::AbstractGeometry) =
    IGeometry(GDAL.ogr_g_difference(g1.ptr, g2.ptr))

unsafe_difference(g1::AbstractGeometry, g2::AbstractGeometry) =
    Geometry(GDAL.ogr_g_difference(g1.ptr, g2.ptr))

"""
    symdifference(g1::AbstractGeometry, g2::AbstractGeometry)

Returns a new geometry representing the symmetric difference of the geometries
or NULL if the difference is empty or an error occurs.
"""
symdifference(g1::AbstractGeometry, g2::AbstractGeometry) =
    IGeometry(GDAL.ogr_g_symdifference(g1.ptr, g2.ptr))

unsafe_symdifference(g1::AbstractGeometry, g2::AbstractGeometry) =
    Geometry(GDAL.ogr_g_symdifference(g1.ptr, g2.ptr))

"""
    distance(g1::AbstractGeometry, g2::AbstractGeometry)

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

"""
    geomlength(geom::AbstractGeometry)

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

"""
    geomarea(geom::AbstractGeometry)

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

"""
    centroid!(geom::AbstractGeometry, centroid::AbstractGeometry)

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::AbstractGeometry)
    result = GDAL.ogr_g_centroid(geom.ptr, centroid.ptr)
    @ogrerr result "Failed to compute the geometry centroid"
    return centroid
end

"""
    centroid(geom::AbstractGeometry)

Compute the geometry centroid.

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)
    point = createpoint()
    centroid!(geom, point)
    return point
end

function unsafe_centroid(geom::AbstractGeometry)
    point = unsafe_createpoint()
    centroid!(geom, point)
    return point
end

"""
    pointalongline(geom::AbstractGeometry, distance::Real)

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.
"""
pointalongline(geom::AbstractGeometry, distance::Real) =
    IGeometry(GDAL.ogr_g_value(geom.ptr, distance))

unsafe_pointalongline(geom::AbstractGeometry, distance::Real) =
    Geometry(GDAL.ogr_g_value(geom.ptr, distance))

"""
    empty!(geom::AbstractGeometry)

Clear geometry information.

This restores the geometry to its initial state after construction, and before
assignment of actual geometry.
"""
function empty!(geom::AbstractGeometry)
    GDAL.ogr_g_empty(geom.ptr)
    return geom
end

"""
    isempty(geom::AbstractGeometry)

Returns `true` if the geometry has no points, otherwise `false`.
"""
isempty(geom::AbstractGeometry) = Bool(GDAL.ogr_g_isempty(geom.ptr))

"""
    isvalid(geom::AbstractGeometry)

Returns `true` if the geometry is valid, otherwise `false`.
"""
isvalid(geom::AbstractGeometry) = Bool(GDAL.ogr_g_isvalid(geom.ptr))

"""
    issimple(geom::AbstractGeometry)

Returns `true` if the geometry is simple, otherwise `false`.
"""
issimple(geom::AbstractGeometry) = Bool(GDAL.ogr_g_issimple(geom.ptr))

"""
    isring(geom::AbstractGeometry)

Returns `true` if the geometry is a ring, otherwise `false`.
"""
isring(geom::AbstractGeometry) = Bool(GDAL.ogr_g_isring(geom.ptr))

"""
    polygonize(geom::AbstractGeometry)

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.
"""
polygonize(geom::AbstractGeometry) = IGeometry(GDAL.ogr_g_polygonize(geom.ptr))

unsafe_polygonize(geom::AbstractGeometry) =
    Geometry(GDAL.ogr_g_polygonize(geom.ptr))

# """
#     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{Cvoid},Cint,
# Ptr{Cvoid},Cint,Ptr{Cvoid},Cint),hGeom,pabyX,nXStride,pabyY,nYStride,pabyZ,
# nZStride)
# end

"""
    getx(geom::AbstractGeometry, i::Integer)

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

"""
    gety(geom::AbstractGeometry, i::Integer)

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

"""
    getz(geom::AbstractGeometry, i::Integer)

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

"""
    getpoint(geom::AbstractGeometry, i::Integer)

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) =
    getpoint!(geom, i, Ref{Cdouble}(), Ref{Cdouble}(), Ref{Cdouble}())

function getpoint!(geom::AbstractGeometry, i::Integer, x, y, z)
    GDAL.ogr_g_getpoint(geom.ptr, i, x, y, z)
    return (x[], y[], z[])
end

"""
    setpointcount!(geom::AbstractGeometry, n::Integer)

Set number of points in a geometry.

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

"""
    setpoint!(geom::AbstractGeometry, i::Integer, x, y)
    setpoint!(geom::AbstractGeometry, i::Integer, x, y, z)

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).
"""
function setpoint! end

function setpoint!(
        geom::AbstractGeometry,
        i::Integer,
        x::Real,
        y::Real,
        z::Real
    )
    GDAL.ogr_g_setpoint(geom.ptr, i, x, y, z)
    return geom
end

function setpoint!(geom::AbstractGeometry, i::Integer, x::Real, y::Real)
    GDAL.ogr_g_setpoint_2d(geom.ptr, i, x, y)
    return geom
end

"""
    addpoint!(geom::AbstractGeometry, x, y)
    addpoint!(geom::AbstractGeometry, x, y, z)

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.
"""
function addpoint! end

function addpoint!(geom::AbstractGeometry, x::Real, y::Real, z::Real)
    GDAL.ogr_g_addpoint(geom.ptr, x, y, z)
    return geom
end

function addpoint!(geom::AbstractGeometry, x::Real, y::Real)
    GDAL.ogr_g_addpoint_2d(geom.ptr, x, y)
    return geom
end

# """
#     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{Cvoid},
# Cint,Ptr{Cvoid},Cint,Ptr{Cvoid},Cint),hGeom,nPointsIn,pabyX,nXStride,pabyY,
# nYStride,pabyZ,nZStride)
# end

"""
    ngeom(geom::AbstractGeometry)

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.ogr_g_getpointcount(geom.ptr)
    n == 0 ? GDAL.ogr_g_getgeometrycount(geom.ptr) : n
end

"""
    getgeom(geom::AbstractGeometry, i::Integer)

Fetch geometry from a geometry container.

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 getgeom(geom::AbstractGeometry, i::Integer)
    # NOTE(yeesian): GDAL.ogr_g_getgeometryref(geom, i) 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.
    if geom.ptr == C_NULL
        return Geometry()
    end
    result = GDAL.ogr_g_getgeometryref(geom.ptr, i)
    if result == C_NULL
        return IGeometry()
    else
        return IGeometry(GDAL.ogr_g_clone(result))
    end
end

function unsafe_getgeom(geom::AbstractGeometry, i::Integer)
    # NOTE(yeesian): GDAL.ogr_g_getgeometryref(geom, i) 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.
    if geom.ptr == C_NULL
        return Geometry()
    end
    result = GDAL.ogr_g_getgeometryref(geom.ptr, i)
    if result == C_NULL
        return Geometry()
    else
        return Geometry(GDAL.ogr_g_clone(result))
    end
end

"""
    addgeom!(geomcontainer::AbstractGeometry, subgeom::AbstractGeometry)

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::AbstractGeometry, subgeom::AbstractGeometry)
    result = GDAL.ogr_g_addgeometry(geomcontainer.ptr, subgeom.ptr)
    @ogrerr result "Failed to add geometry. The geometry type could be illegal"
    return 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::AbstractGeometry,
#         subgeom::AbstractGeometry
#     )
#     result = GDAL.ogr_g_addgeometrydirectly(geomcontainer.ptr, subgeom.ptr)
#     @ogrerr result "Failed to add geometry. The geometry type could be illegal"
#     return geomcontainer
# end

"""
    removegeom!(geom::AbstractGeometry, i::Integer, todelete::Bool = true)

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.ogr_g_removegeometry(geom.ptr, i, todelete)
    @ogrerr result "Failed to remove geometry. The index could be out of range."
    return geom
end

"""
    removeallgeoms!(geom::AbstractGeometry, todelete::Bool = true)

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.ogr_g_removegeometry(geom.ptr, -1, todelete)
    @ogrerr result "Failed to remove all geometries."
    return geom
end

"""
    hascurvegeom(geom::AbstractGeometry, nonlinear::Bool)

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.ogr_g_hascurvegeometry(geom.ptr, nonlinear))

"""
    lineargeom(geom::AbstractGeometry, stepsize::Real = 0)

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.

### 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.
"""
lineargeom(geom::AbstractGeometry, stepsize::Real = 0) =
    IGeometry(GDAL.ogr_g_getlineargeometry(geom.ptr, stepsize, C_NULL))

unsafe_lineargeom(geom::AbstractGeometry, stepsize::Real = 0) =
    Geometry(GDAL.ogr_g_getlineargeometry(geom.ptr, stepsize, C_NULL))

lineargeom(geom::AbstractGeometry, options::Vector, stepsize::Real = 0) =
    IGeometry(GDAL.ogr_g_getlineargeometry(geom.ptr, stepsize, options))

function unsafe_lineargeom(
        geom::AbstractGeometry,
        options::Vector,
        stepsize::Real = 0
    )
    return Geometry(GDAL.ogr_g_getlineargeometry(geom.ptr, stepsize, options))
end

"""
    curvegeom(geom::AbstractGeometry)

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 reverse function is OGR_G_GetLinearGeometry().
"""
curvegeom(geom::AbstractGeometry) =
    IGeometry(GDAL.ogr_g_getcurvegeometry(geom.ptr, C_NULL))

unsafe_curvegeom(geom::AbstractGeometry) =
    Geometry(GDAL.ogr_g_getcurvegeometry(geom.ptr, C_NULL))

"""
    polygonfromedges(lines::AbstractGeometry, tol::Real; besteffort = false, autoclose = false)

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.
* `tol`: whether two arcs are considered close enough to be joined.

### Keyword Arguments
* `besteffort`: (defaults to `false`) not yet implemented???.
* `autoclose`: indicates if the ring should be close when first and last
    points of the ring are the same. (defaults to `false`)
"""
function polygonfromedges(
        lines::AbstractGeometry,
        tol::Real;
        besteffort::Bool = false,
        autoclose::Bool = false
    )
    perr = Ref{GDAL.OGRErr}()
    result = GDAL.ogrbuildpolygonfromedges(lines.ptr, besteffort, autoclose,
        tol, perr)
    @ogrerr perr[] "Failed to build polygon from edges."
    return IGeometry(result)
end

function unsafe_polygonfromedges(
        lines::AbstractGeometry,
        tol::Real;
        besteffort::Bool = false,
        autoclose::Bool = false
    )
    perr = Ref{GDAL.OGRErr}()
    result = GDAL.ogrbuildpolygonfromedges(lines.ptr, besteffort, autoclose,
        tol, perr)
    @ogrerr perr[] "Failed to build polygon from edges."
    return Geometry(result)
end

"""
    setnonlineargeomflag!(flag::Bool)

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.ogrsetnonlineargeometriesenabledflag(flag)

"""
    getnonlineargeomflag()

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

for (geom, wkbgeom) in ((:geomcollection,       GDAL.wkbGeometryCollection),
                        (:linestring,           GDAL.wkbLineString),
                        (:linearring,           GDAL.wkbLinearRing),
                        (:multilinestring,      GDAL.wkbMultiLineString),
                        (:multipoint,           GDAL.wkbMultiPoint),
                        (:multipolygon,         GDAL.wkbMultiPolygon),
                        (:multipolygon_noholes, GDAL.wkbMultiPolygon),
                        (:point,                GDAL.wkbPoint),
                        (:polygon,              GDAL.wkbPolygon))
    eval(quote
        $(Symbol("create$geom"))() = creategeom($wkbgeom)
        $(Symbol("unsafe_create$geom"))() = unsafe_creategeom($wkbgeom)
    end)
end

for f in (:create, :unsafe_create)
    for (args, typedargs) in (
            ((:x,:y), (:(x::Real),:(y::Real))),
            ((:x,:y,:z), (:(x::Real),:(y::Real),:(z::Real)))
        )
        eval(quote
            function $(Symbol("$(f)point"))($(typedargs...))
                geom = $(Symbol("$(f)point"))()
                addpoint!(geom, $(args...))
                return geom
            end
        end)
    end

    for (args, typedargs) in (
            ((:xs,:ys),     (:(xs::Vector{Cdouble}),
                             :(ys::Vector{Cdouble}))),
            ((:xs,:ys,:zs), (:(xs::Vector{Cdouble}),
                             :(ys::Vector{Cdouble}),
                             :(zs::Vector{Cdouble})))
        )
        for geom in (:linestring, :linearring)
            eval(quote
                function $(Symbol("$f$geom"))($(typedargs...))
                    geom = $(Symbol("$f$geom"))()
                    for pt in zip($(args...))
                        addpoint!(geom, pt...)
                    end
                    return geom
                end
            end)
        end

        for (geom,component) in ((:polygon, :linearring),)
            eval(quote
                function $(Symbol("$f$geom"))($(typedargs...))
                    geom = $(Symbol("$f$geom"))()
                    subgeom = $(Symbol("unsafe_create$component"))($(args...))
                    result = GDAL.ogr_g_addgeometrydirectly(
                        geom.ptr,
                        subgeom.ptr
                    )
                    @ogrerr result "Failed to add $component."
                    return geom
                end
            end)
        end

        for (geom,component) in ((:multipoint, :point),)
            eval(quote
                function $(Symbol("$f$geom"))($(typedargs...))
                    geom = $(Symbol("$f$geom"))()
                    for pt in zip($(args...))
                        subgeom = $(Symbol("unsafe_create$component"))(pt)
                        result = GDAL.ogr_g_addgeometrydirectly(
                            geom.ptr,
                            subgeom.ptr
                        )
                        @ogrerr result "Failed to add point."
                    end
                    return geom
                end
            end)
        end
    end

    for typeargs in (Vector{<:Real},
                     Tuple{<:Real,<:Real},
                     Tuple{<:Real,<:Real,<:Real})
        eval(quote
            function $(Symbol("$(f)point"))(coords::$typeargs)
                geom = $(Symbol("$(f)point"))()
                addpoint!(geom, coords...)
                return geom
            end
        end)
    end

    for typeargs in (Vector{Tuple{Cdouble,Cdouble}},
                     Vector{Tuple{Cdouble,Cdouble,Cdouble}},
                     Vector{Vector{Cdouble}})
        for geom in (:linestring, :linearring)
            eval(quote
                function $(Symbol("$f$geom"))(coords::$typeargs)
                    geom = $(Symbol("$f$geom"))()
                    for coord in coords
                        addpoint!(geom, coord...)
                    end
                    return geom
                end
            end)
        end

        for (geom,component) in ((:polygon, :linearring),)
            eval(quote
                function $(Symbol("$f$geom"))(coords::$typeargs)
                    geom = $(Symbol("$f$geom"))()
                    subgeom = $(Symbol("unsafe_create$component"))(coords)
                    result = GDAL.ogr_g_addgeometrydirectly(
                        geom.ptr,
                        subgeom.ptr
                    )
                    @ogrerr result "Failed to add $component."
                    return geom
                end
            end)
        end
    end

    for (variants,typeargs) in (
            (((:multipoint, :point),),
             (Vector{Tuple{Cdouble,Cdouble}},
              Vector{Tuple{Cdouble,Cdouble,Cdouble}},
              Vector{Vector{Cdouble}})),

            (((:polygon, :linearring),
              (:multilinestring, :linestring),
              (:multipolygon_noholes, :polygon)),
             (Vector{Vector{Tuple{Cdouble,Cdouble}}},
              Vector{Vector{Tuple{Cdouble,Cdouble,Cdouble}}},
              Vector{Vector{Vector{Cdouble}}})),

            (((:multipolygon, :polygon),),
             (Vector{Vector{Vector{Tuple{Cdouble,Cdouble}}}},
              Vector{Vector{Vector{Tuple{Cdouble,Cdouble,Cdouble}}}},
              Vector{Vector{Vector{Vector{Cdouble}}}}))
        )
        for typearg in typeargs, (geom, component) in variants
            eval(quote
                function $(Symbol("$f$geom"))(coords::$typearg)
                    geom = $(Symbol("$f$geom"))()
                    for coord in coords
                        subgeom = $(Symbol("unsafe_create$component"))(coord)
                        result = GDAL.ogr_g_addgeometrydirectly(
                            geom.ptr,
                            subgeom.ptr
                        )
                        @ogrerr result "Failed to add $component."
                    end
                    return geom
                end
            end)
        end
    end
end


text:
using Downloads, Tar, p7zip_jll, SimpleBufferStream, SHA

export unpack, list_tarball_files, verify, download_verify

function detect_compressor(header::Vector)
    compressor_magic_bytes = Dict(
        "gzip" => [0x1f, 0x8b],
        "xz" => [0xfd, 0x37, 0x7a, 0x58, 0x5A, 0x00],
        "zstd" => [0x28, 0xB5, 0x2F, 0xFD],
        "bzip2" => [0x42, 0x5a, 0x68],
    )
    for (compressor, magic) in compressor_magic_bytes
        lm = length(magic)
        if length(header) >= lm && header[1:lm] == magic
            return compressor
        end
    end
    return nothing
end

function decompress(path::AbstractString)
    # Read the first few bytes of data to classify it:
    compressor = open(path) do io
        detect_compressor(read(io, 6))
    end
    if compressor === nothing
        error("Called decompress() on uncompressed file")
    end

    p7zip() do p7z
        # Launch p7zip as our decompressor engine, clueing it in to the compressor
        p = open(`$p7z e $(path) -so -t$(compressor)`; read=true)
        return p.out
    end
end

# Many functions don't like `PipeEndpoint`, so we interface with a BufferStream
function buff_wrap(io::IO; blocksize = 2*1024*1024)
    buff = BufferStream(blocksize)
    @async begin
        while !eof(io)
            write(buff, readavailable(io))
        end
        close(buff)
    end
    return buff
end

# Reimplement Pkg.PlatformEngines verbs using Tar and Downloader
function unpack(tarball_path::AbstractString, dest::AbstractString)
    mkpath(dest)

    @sync begin
        Tar.extract(decompress(tarball_path), dest)
    end
end

function list_tarball_files(tarball_path::AbstractString)
    @sync begin
        # Decompress to a file, since `Tar.list()` can't feed off of a BufferStream yet
        mktemp() do path, _
            function remove_dotslash(path)
                if startswith(path, "./")
                    return string(path[3:end])
                end
                return path
            end
            return String[remove_dotslash(h.path) for h in Tar.list(buff_wrap(decompress(tarball_path)))]
        end
    end
end

function verify(path::AbstractString, hash::AbstractString; hash_path::AbstractString="$(path).sha256")
    # Check hash string format
    if !occursin(r"^[0-9a-f]{64}$"i, hash)
        msg = "Hash value must be 64 hexadecimal characters (256 bits), "
        if !isascii(hash)
            msg *= "given hash value is non-ASCII"
        elseif occursin(r"^[0-9a-f]*$"i, hash)
            msg *= "given hash value has the wrong length ($(length(hash)))"
        else
            msg *= "given hash value contains non-hexadecimal characters"
        end
        msg *= ": $(repr(hash))"
        error(msg)
    end
    hash = lowercase(hash)

    # First, it must exist
    if isfile(hash_path)
        # Next, it must contain the same hash as what we're verifying against
        if read(hash_path, String) == hash
            # Next, it must be no older than the actual path
            if stat(hash_path).mtime >= stat(path).mtime
                return true
            end
        end
    end

    calc_hash = open(path) do file
        bytes2hex(sha256(file))
    end
    @assert occursin(r"^[0-9a-f]{64}$", calc_hash)

    if calc_hash != hash
        msg  = "Hash Mismatch!\n"
        msg *= "  Expected sha256:   $hash\n"
        msg *= "  Calculated sha256: $calc_hash"
        @error(msg)
        return false
    end

    # Try to save a hash cache if everything worked out fine
    try
        open(hash_path, "w") do file
            write(file, hash)
        end
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
    end

    return true
end

function download_verify(url, hash, path)
    if isfile(path) && verify(path, hash)
        @info "Cached file found in $(path)"
        return true
    else
        mkpath(dirname(path))
        # `curl` doesn't automatically overwrite an existing file, delete it
        # before attempting the download.
        rm(path; force=true)
        @info "Downloading $(url) to $(path)..."
        # Temporarily shell out to `curl` to download, until `Downloads` bugs are squashed
        try
            run(`curl -C - -s -\# -f -L $(url) -o $(path)`)
        catch e
            # Downloads throws an ErrorException, so we'll do the same
            error("download failed: $(e)")
        end
        verify(path, hash) || error("Verification failed")
    end
end

ext: jl


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

main_paths = ["."]
for main_path in main_paths
    format(
        main_path;
        whitespace_ops_in_indices=true,
        remove_extra_newlines=true,
        verbose=true,
        always_for_in=true,
        whitespace_typedefs=true,
        whitespace_in_kwargs=false,
        format_docstrings=true,
        always_use_return=false, # removed since it has false positives.
    )
end


text:
if basename(pwd()) == "aoc"
    cd("2015/4")
end
using Pkg
Pkg.activate(".")
using MD5

function hashsearch(needle, key)
    for i in Iterators.countfrom(1)
        if startswith(bytes2hex(md5(key * string(i))), needle)
            return i
        end
    end
end

@assert hashsearch("00000", "abcdef") == 609043
@assert hashsearch("00000", "pqrstuv") == 1048970
hashsearch("00000", "yzbqklnj")
hashsearch("000000", "yzbqklnj")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random
using ArraysOfArrays, IntervalSets

@testset "SpatialVolume" begin
    lo = [-1., -0.1]
    hi = [2.,1]
    hyperRectVolume = BAT.HyperRectVolume(lo, hi)
    @testset "BAT.SpatialVolume" begin
        @test eltype(@inferred BAT.HyperRectVolume([-1., 0.5], [2.,1])) == Float64
        @test size(rand(MersenneTwister(7002), hyperRectVolume)) == (2,)
        @test size(rand(MersenneTwister(7002), hyperRectVolume, 3)) == (3,)
    end
    @testset "BAT.HyperRectVolume" begin
        @test typeof(@inferred BAT.HyperRectVolume([-1., 0.5], [2.,1])) <: BAT.SpatialVolume{Float64}
        @test_throws ArgumentError BAT.HyperRectVolume([-1.], [2.,1])

        @test [0.0, 0.0] in hyperRectVolume
        @test ([0.5, 2] in hyperRectVolume) == false

        @test ndims(hyperRectVolume) == 2

        res = @inferred rand!(MersenneTwister(7002), hyperRectVolume, zeros(2))
        @test typeof(res) <: AbstractArray{Float64, 1}
        @test size(res) == (2,)
        @test res in hyperRectVolume
        res = @inferred rand!(MersenneTwister(7002), hyperRectVolume, VectorOfSimilarVectors(zeros(2,3)))
        @test typeof(res) <: VectorOfSimilarVectors{Float64,Array{Float64,2}}
        @test size(res) == (3,)

        res = @inferred similar(hyperRectVolume)
        @test typeof(res) <: BAT.HyperRectVolume{Float64}

        a = @inferred BAT.HyperRectVolume([0.0,0.0],[1.0,1.0])
        b = @inferred BAT.HyperRectVolume([0.5,0.9],[0.8,1.5])
        res = @inferred intersect(a, b)
        @test res.lo ≈ [0.5, 0.9]
        @test res.hi ≈ [0.8, 1.0]
        @test isempty(res) == false

        a = @inferred BAT.HyperRectVolume([0.0,0.0],[0.3,0.4])
        res = @inferred intersect(a, b)
        @test res.lo ≈ [0.5, 0.9]
        @test res.hi ≈ [0.3, 0.4]
        @test isempty(res)
    end

    @testset "log_volume" begin
        @test BAT.log_volume(@inferred BAT.HyperRectVolume([-1.,0.],[1.,3.])) ≈
            1.79175946
        @test BAT.log_volume(@inferred BAT.HyperRectVolume([-0.0001],[0.0000])) ≈
            -9.210340371
    end

    @testset "fromuhc" begin
        @test inv(BAT.fromuhc!) == BAT.inv_fromuhc!
        @test inv(BAT.inv_fromuhc!) == BAT.fromuhc!
        @test inv(BAT.fromuhc) == BAT.inv_fromuhc
        @test inv(BAT.inv_fromuhc) == BAT.fromuhc

        y = zeros(Float64,2)
        x = similar(y)
        BAT.fromuhc!(x, y, hyperRectVolume)
        @test x ≈ hyperRectVolume.lo
        u = similar(y)
        inv(BAT.fromuhc!)(u, x, hyperRectVolume)
        @test u ≈ y
        y = ones(Float64, 2)
        res = @inferred BAT.fromuhc(y, hyperRectVolume)
        @test res ≈ hyperRectVolume.hi
        @test inv(BAT.fromuhc)(res, hyperRectVolume) ≈ y
        y = VectorOfSimilarVectors([0.2 0.5 0.5; 0.9 0.7 0.5])
        res = @inferred BAT.fromuhc(y, hyperRectVolume)
        @test res ≈ VectorOfSimilarVectors([-0.4 0.5 0.5; 0.89 0.67 0.45])
    end
end


text:
@testset "Variables basics" begin

    let s = """
            function f()
            x = 0
            x
            end
            """
        @test lintstr(s) == []
    end
    let s = """
            function f()
            x = zeros(10, 10)
            x
            end
            """
        @test lintstr(s) == []
    end
end


text:
@testset "agent_count" begin
    grass1 = Grass(1,1,5)
    grass2 = Grass(2,2,5)
    sheep  = Sheep(3,1,1,1,1)
    wolf   = Wolf(5,2,2,2,2)
    world  = World([sheep,grass1,grass2,wolf])

    @test agent_count(grass1) ≈ 0.2
    @test agent_count(sheep) == 1
    @test agent_count([grass2,grass1]) ≈ 0.6
    res = agent_count(world)
    tst = Dict(:Sheep=>1,:Wolf=>1,:Grass=>0.6)
    for (k,_) in res
        @test res[k] ≈ tst[k]
    end
end

@testset "Mushroom" begin
    sheep = Sheep(5,10,2,1,1)
    mushr = Mushroom(2,5,5)
    world = World([sheep,mushr])

    EcosystemCore.eat!(sheep,mushr,world)
    @test size(mushr) == 0
    @test energy(sheep) == 0
end

@testset "every_nth" begin
    i = 0
    f() = i+=1
    cb = every_nth(f,3)

    cb(); cb();
    @test i == 0
    cb()
    @test i == 1
end


text:
module InfraFlow

include("delim.jl")

using JuMP, GLPK, Test, YAML, Logging
const MOI = JuMP.MathOptInterface

export get_data, formulate_gmcnf!, run
"""
Author: Will Usher
Date: 3rd January 2018

A linear programming new_capacity expansion generalized multi-commodity network
flow implemented using JuMP/Julia based on generalized multi-commodity network
flow as found in:

Ishimatsu, Takuto. “Generalized Multi-Commodity Network Flows : Case Studies
in Space Logistics and Complex Infrastructure Systems.”
Masssachusett Institute of Technology, 2013.
"""

"""
    get_data(file_path)

Load in model data from a YAML file
"""
function get_data(file_path::AbstractString)

    data = YAML.load(open(file_path))

    model_data = Dict{String,Any}()

    model_data["years"] = data["years"]
    years = Dict{Int,Int}([data["years"][i]=>i for i in eachindex(data["years"])])
    model_data["commodities"] = data["commodities"]
    commodities = Dict{String,Int}([data["commodities"][i]=>i for i in eachindex(data["commodities"])])

    model_data["discount_rate"] = data["discount_rate"]

    nodes = Dict{String,Int}()
    node_names = []
    self_loops = []
    for index in eachindex(data["nodes"])
        node = data["nodes"][index]
        nodes[node["name"]] = index
        push!(node_names, node["name"])
        if haskey(node, "requirements")
            push!(self_loops, index)
        end

    end
    model_data["nodes"] = node_names

    num_nodes = length(node_names)
    num_comm = length(commodities)
    num_years = length(years)

    cap2act = zeros(Float64, (num_nodes, num_nodes, num_comm))
    outflow_cost = zeros(Float64, (num_nodes, num_nodes, num_comm))
    inflow_cost = zeros(Float64,(num_nodes, num_nodes, num_comm))
    residual_capacity = zeros(Float64, (num_nodes, num_nodes, num_comm, num_years))
    maximum_new_capacity = fill(99999999999999999999.0, (num_nodes, num_nodes, num_comm, num_years))
    transformation = zeros(Float64, (num_nodes, num_nodes, num_comm, num_comm))

    requirements_outflow = zeros(Float64, (num_nodes, num_nodes, num_comm, num_comm))
    requirements_inflow = zeros(Float64, (num_nodes, num_nodes, num_comm, num_comm))

    edges = []
    for edge in data["edges"]
        source = nodes[edge["source"]]
        sink = nodes[edge["sink"]]
        push!(edges, source, sink)

        for flow in edge["flow"]
            comm_idx = commodities[flow["name"]]
            requirements_inflow[source, sink, comm_idx, comm_idx] = 1
            requirements_outflow[source, sink, comm_idx, comm_idx] = 1
        end

        if haskey(edge, "cap2act")
            for (comm, value) in edge["cap2act"]
                comm_idx = commodities[comm]
                cap2act[source, sink, comm_idx] = value
            end
        end

        if haskey(edge, "operational_cost")
            for (comm, value) in edge["operational_cost"]
                comm_idx = commodities[comm]
                outflow_cost[source, sink, comm_idx] = value
            end
        end

        if haskey(edge, "losses")
            for (comm, value) in edge["losses"]
                comm_idx = commodities[comm]
                transformation[source, sink, comm_idx, comm_idx] = (1.0 - value)
            end
        end

    end

    for loop in self_loops
        push!(edges, loop, loop)
    end

    edges = transpose(reshape(edges, 2, :))
    model_data["edges"] = edges

    demand = zeros((num_nodes, num_comm, num_years))
    capacity_cost = zeros((num_nodes, num_nodes, num_comm))

    for index in eachindex(data["nodes"])
        node = data["nodes"][index]
        source = node["name"]
        node_idx = nodes[source]
        if haskey(node, "demand")
            for (commodity, yearly_data) in node["demand"]
                for (year, value) in yearly_data
                    node_idx = nodes[source]
                    comm_index = commodities[commodity]
                    year_idx = years[year]
                    demand[node_idx, comm_index, year_idx] = value
                end
            end
        end

        if haskey(node, "investment_cost")
            output_commodity = node["output"]
            output_comm_idx = commodities[output_commodity]
            capacity_cost[node_idx, node_idx, output_comm_idx] = node["investment_cost"]
        end

        if haskey(node, "cap2act")
            output_commodity = node["output"]
            output_comm_idx = commodities[output_commodity]
            cap2act[node_idx, node_idx, output_comm_idx] = node["cap2act"]
        end

        if haskey(node, "requirements")
            for requirement in node["requirements"]
                comm = node["output"]
                comm_idx = commodities[comm]
                r_comm = requirement["name"]
                r_comm_idx = commodities[r_comm]
                cap2act[node_idx, node_idx, r_comm_idx] = 1
                transformation[node_idx, node_idx, comm_idx, r_comm_idx] = 1
                requirements_outflow[node_idx, node_idx, r_comm_idx, comm_idx] = requirement["value"]
                requirements_inflow[node_idx, node_idx, comm_idx, comm_idx] = 1
                requirements_outflow[node_idx, node_idx, comm_idx, comm_idx] = 1
                if haskey(node, "residual_capacity")
                    for (year, value) in node["residual_capacity"]
                        year_idx = years[year]
                        residual_capacity[node_idx, node_idx, comm_idx, year_idx] = value
                    end
                end
                if haskey(node, "maximum_new_capacity")
                    for (year, value) in node["maximum_new_capacity"]
                        year_idx = years[year]
                        maximum_new_capacity[node_idx, node_idx, comm_idx, year_idx] = value
                    end
                end
            end

        end

    end

    model_data["demand"] = demand
    model_data["capacity_cost"] = capacity_cost
    model_data["cap2act"] = cap2act
    model_data["outflow_cost"] = outflow_cost
    model_data["inflow_cost"] = inflow_cost
    model_data["residual_capacity"] = residual_capacity
    model_data["maximum_new_capacity"] = maximum_new_capacity
    model_data["transformation"] = transformation
    model_data["requirements_inflow"] = requirements_inflow
    model_data["requirements_outflow"] = requirements_outflow

    source_nodes = view(edges, :, 1)
    sink_nodes = view(edges, :, 2)

    outflow_edges = make_edge_dict(source_nodes, sink_nodes)
    inflow_edges = make_edge_dict(sink_nodes, source_nodes)

    model_data["outflow_edges"] = outflow_edges
    model_data["inflow_edges"] = inflow_edges

    return model_data
end

"""
    make_edge_dict(edge_nodes, other_nodes)

"""
function make_edge_dict(edge_nodes::AbstractArray, other_nodes::AbstractArray)
    edges = Dict{Int8,Array{Int8}}()
    for node in Set(edge_nodes)
        for i in eachindex(edge_nodes)
            if edge_nodes[i] == node
                if haskey(edges, node)
                    push!(edges[node], other_nodes[i])
                else
                    push!(edges, node=>[other_nodes[i]])
                end
            end
        end
    end
    return edges
end

"""
    formulate_gmcnf()

``\\sum_{ijky} c_{ijky}x_{ijky}``

"""
function formulate_gmcnf!(model_data::Dict{String,Any}, model; verbose = true)

    nodes = model_data["nodes"]
    edges = model_data["edges"]
    commodities = model_data["commodities"]
    discount_rate = model_data["discount_rate"]
    years = model_data["years"]

    outflow_edges = model_data["outflow_edges"]
    inflow_edges = model_data["inflow_edges"]

    num_nodes = length(nodes)
    num_edges = length(edges)
    num_comm = length(commodities)
    num_years = length(years)

    # demand at node by commodity
    demand = model_data["demand"]

    capacity_cost = model_data["capacity_cost"]
    cap2act = model_data["cap2act"]

    inflow_cost = model_data["inflow_cost"]
    outflow_cost = model_data["outflow_cost"]

    # describe the commodity requirements for an in- or out-flow
    requirements_outflow = model_data["requirements_outflow"]
    requirements_inflow = model_data["requirements_inflow"]

    # describe the flow gain/loss or transformation between commodities
    transformation = model_data["transformation"]

    # upper bound on operational decision variables (new_capacity)
    residual_capacity = model_data["residual_capacity"]

    # upper bound on new capacity
    maximum_new_capacity = model_data["maximum_new_capacity"]

    new_capacity = @variable(model,
                         new_capacity[i=1:num_nodes, j=1:num_nodes, k=1:num_comm, y=1:num_years],
                         lower_bound = 0)

    outflow = @variable(model,
                        outflow[i=1:num_nodes, j=1:num_nodes, k=1:num_comm, y=1:num_years],
                        lower_bound = 0)
    inflow = @variable(model,
                       inflow[i=1:num_nodes, j=1:num_nodes, k=1:num_comm, y=1:num_years],
                       lower_bound = 0)

    @constraint(
        model,
        max_new_capacity[i=1:num_nodes, j=1:num_nodes, k=1:num_comm, y=1:num_years],
        new_capacity[i, j, k, y] <= maximum_new_capacity[i, j, k, y]
    )

    @variable(
        model,
        total_annual_capacity[i in keys(outflow_edges), j in outflow_edges[i], k=1:num_comm, y=1:num_years],
        lower_bound = 0
    )

    @constraint(
        model,
        accumulate_capacity[i in keys(outflow_edges), j in outflow_edges[i], k=1:num_comm, y=1:num_years],
        total_annual_capacity[i, j, k, y] == residual_capacity[i, j, k, y] + sum(new_capacity[i, j, k, z] for z in 1:y)
    )

    @constraint(
        model,
        capacity_exp_outflow[i in keys(outflow_edges), j in outflow_edges[i], k=1:num_comm, y=1:num_years],
        outflow[i, j, k, y] <= total_annual_capacity[i, j, k, y] * cap2act[i, j, k])

    @constraint(
        model,
        capacity_exp_inflow[i in keys(outflow_edges), j in outflow_edges[i], k=1:num_comm, y=1:num_years],
        inflow[i, j, k, y] <= total_annual_capacity[i, j, k, y] * cap2act[i, j, k])

    function discount_factor(year::Int)
        return ((1 + discount_rate) ^ (years[year] - years[1]))
    end

    discounted_capital_cost = @expression(
        model,
        [i in keys(outflow_edges), j in outflow_edges[i], k=1:num_comm, y=1:num_years],
        new_capacity[i, j, k, y] * capacity_cost[i, j, k] / discount_factor(y)
    )

    discounted_operational_cost = @expression(
        model,
        [i in keys(outflow_edges), j in outflow_edges[i], k=1:num_comm, y=1:num_years],
        outflow_cost[i, j, k] * outflow[i, j, k, y]
        + inflow_cost[i, j, k] * inflow[i, j, k, y]
        / discount_factor(y)
    )

    @objective(
        model,
        Min,
        sum(discounted_operational_cost[i, j, k, y]
            + discounted_capital_cost[i, j, k, y]
            for i in keys(outflow_edges), j in outflow_edges[i], k in 1:num_comm, y=1:num_years)
        )

    requirements_outflow_const = @expression(
        model,
        requirements_outflow_const[i=1:num_nodes, k=1:num_comm, y=1:num_years],
        if haskey(outflow_edges, i)
            sum(
                sum(requirements_outflow[i, j, k, l] for l in 1:num_comm) * outflow[i, j, k, y]
                for j in outflow_edges[i]
            )

        else
            @info "No outflow edge found for ($(nodes[i]), $(commodities[k]), $(years[y]))"
            0
        end
    )

    requirements_inflow_const = @expression(
        model,
        requirements_inflow_const[i=1:num_nodes, k=1:num_comm, y=1:num_years],
        if haskey(inflow_edges, i)
            sum(
                sum(requirements_inflow[h, i, k, m] for m in 1:num_comm) * inflow[h, i, k, y]
                for h in inflow_edges[i]
            )

        else
            @info "No inflow edge found for ($(nodes[i]), $(commodities[k]), $(years[y]))"
            0
        end
    )

    mass_balance = @constraint(
        model,
        mass_balance[i=1:num_nodes, k=1:num_comm, y=1:num_years],
        requirements_outflow_const[i, k, y] - requirements_inflow_const[i, k, y]
        <= demand[i, k, y])

    flow_transformation = @constraint(
        model,
        flow_transformation[i in keys(outflow_edges), j in outflow_edges[i], k in 1:num_comm, l in 1:num_comm, y=1:num_years; transformation[i, j, k, l] != 0],
        transformation[i, j, k, l] * outflow[i, j, l, y] == inflow[i, j, k, y]
        )

    return model
end

function print_vars(var_object)
    for var in var_object
        if JuMP.value(var) != 0.0
            println("$(var): $(JuMP.value(var))")
        end
    end
end

function print_duals(con_object)
    for con in con_object
        if JuMP.shadow_price(con) != 0.0
            println("$(con): $(JuMP.shadow_price(con))")
        end
    end
end

function print_constraint(con_object)
    for con in con_object
        println("$(con)")
    end
end

function run(file_path::String)

    model_data = get_data(file_path)

    model = Model(with_optimizer(GLPK.Optimizer))

    @time model = formulate_gmcnf!(model_data, model, verbose = true)

    println("Compiled model, now running")
    @time JuMP.optimize!(model)
    println("Finished running, objective: £$(JuMP.objective_value(model))")

    status = JuMP.termination_status(model) == MOI.OPTIMAL
    feasible = JuMP.primal_status(model) == MOI.FEASIBLE_POINT

    if status
        println("\nSolution optimal\n")
        outflow = model[:outflow]
        inflow = model[:inflow]
        total_capacity = model[:total_annual_capacity]
        new_capacity = model[:new_capacity]

        write_results(outflow, "outflow", model_data)
        write_results(inflow, "inflow", model_data)
        write_results(total_capacity, "total_capacity", model_data)
        write_results(new_capacity, "new_capacity", model_data)

    else
        println("\nSolution not optimal\n")
        println("i,j=> $(model_data["nodes"])")
        println("k=> $(model_data["commodities"])")
        println("y=> $(model_data["years"])")
        println("mass_balance[i, k, y]")
        print_constraint(model[:mass_balance])
        println("flow_transformation[i, j, k, y]")
        print_constraint(model[:flow_transformation])
        print_constraint(model[:max_new_capacity])
        print_constraint(model[:accumulate_capacity])
    end

end

if length(ARGS) > 0
    run(ARGS[1])
end

end


text:
    __precompile__()
"""
A package that provides a new REPL that has syntax highlighting,
bracket matching and other nifty features.
"""
module OhMyREPL

using Tokenize
using Crayons

export colorscheme!, colorschemes, enable_autocomplete_brackets, test_colorscheme

include("repl_pass.jl")
include("repl.jl")
include(joinpath("passes", "Passes.jl"))

include("BracketInserter.jl")
include("prompt.jl")

# Some backward compatability
module ANSICodes
    using Crayons
    const ANSIToken = Crayon
end

export ANSITokens

function colorscheme!(name::String)
    Passes.SyntaxHighlighter.activate!(Passes.SyntaxHighlighter.SYNTAX_HIGHLIGHTER_SETTINGS,
                                       name)
end

function colorschemes()
    show(Passes.SyntaxHighlighter.SYNTAX_HIGHLIGHTER_SETTINGS)
end

const TEST_STR = """

function funcdef(x::Float64, y::Int64)
    y = 100_000
    x = :foo
    s = "I am a happy string"
    @time 1+1
    #= Comments look like this =#
    z = funccall(x, y)
    5 * 3 + 2 - 1
end
"""

function test_colorscheme(name::String, str::String = TEST_STR)
    syntaxpass = get_pass(PASS_HANDLER, "SyntaxHighlighter")
    active = syntaxpass.active
    try
        colorscheme!(name)
        test_pass(syntaxpass, str)
    finally
        syntaxpass.active = active
    end
    return
end

function test_colorscheme(cs::Passes.SyntaxHighlighter.ColorScheme, str::String = TEST_STR)
    syntaxpass = get_pass(PASS_HANDLER, "SyntaxHighlighter")
    active = syntaxpass.active
    name = "plzdontnameyourcolorschemethis"
    try
        Passes.SyntaxHighlighter.add!(syntaxpass, name, cs)
        colorscheme!(name)
        test_pass(syntaxpass, str)
    finally
        syntaxpass.active = active
        if haskey(syntaxpass.schemes, name)
            delete!(syntaxpass.schemes, name)
        end
    end
    return
end

showpasses(io::IO = STDOUT) = Base.show(io, PASS_HANDLER)

function __init__()
    # Only load in interactive mode, see #68
    if !isempty(ARGS) && !isempty(first(ARGS))
        return
    end
    options = Base.JLOptions()
    # command-line
    if (options.eval != C_NULL) || (options.print != C_NULL)
        return
    end

    if isdefined(Base, :active_repl)
        Prompt.insert_keybindings()
    else
        atreplinit() do repl
            repl.interface = Base.REPL.setup_interface(repl)
            Prompt.insert_keybindings()
            update_interface(repl.interface)
            main_mode = repl.interface.modes[1]
            p = repl.interface.modes[5]
            # These are inserted here because we only want to insert them for the Julia mode
            d = Dict(
            # Up Arrow
            "\e[A" => (s,o...)-> begin
                Base.LineEdit.edit_move_up(s) || Base.LineEdit.enter_prefix_search(s, p, true)
                Prompt.rewrite_with_ANSI(s)
            end,
            # Down Arrow
            "\e[B" => (s,o...)-> begin
                 Base.LineEdit.edit_move_down(s) || Base.LineEdit.enter_prefix_search(s, p, false)
                 Prompt.rewrite_with_ANSI(s)
            end
            )
            main_mode.keymap_dict = Base.LineEdit.keymap([d, main_mode.keymap_dict])
        end
    end

    mktemp() do _, f
        old_stderr = STDERR
        redirect_stderr(f)
        Base.LineEdit.refresh_line(s) = (Base.LineEdit.refresh_multi_line(s); OhMyREPL.Prompt.rewrite_with_ANSI(s))
        include(joinpath(dirname(@__FILE__), "output_prompt_overwrite.jl"))
        redirect_stderr(old_stderr)
    end
end

end # module


text:
@testset "KDTree" begin

tree1 = searchstructure(KDTree, data, Euclidean(); reorder = true)
tree2 = searchstructure(KDTree, data, Euclidean(); reorder = false)

@test datatype(tree1) === datatype(typeof(tree1)) === eltype(data)
@test getmetric(tree1) === tree1.metric
@test datatype(tree2) === datatype(typeof(tree2)) === eltype(data)
@test getmetric(tree2) === tree2.metric

idxs, ds = knn(tree1, query, 5)
@test issorted(ds)
@test isearch(tree1, query, NeighborNumber(5)) == idxs
@test search(tree1, query, NeighborNumber(5)) == (idxs, ds)

ridxs, rds = inrange(tree1, query, maximum(ds))
@test issorted(rds)
@test ridxs == idxs
@test rds == ds
ridxs_srt, rds_srt = inrange(tree2, query,  maximum(ds))
@test ridxs_srt == idxs
@test rds_srt == ds

__idxs, = inrange(tree1, queries[1], 0.005)
@test sort!(__idxs) == 48:52
__idxs, = inrange(tree1, queries[1], 0.005, theiler1(nidxs[1]))
@test sort!(__idxs) == 50:52
__idxs, = inrange(tree1, queries[2], 0.005, theiler2(nidxs[2]))
@test isempty(__idxs)

vec_of_idxs, vec_of_ds = bulksearch(tree1, queries, NeighborNumber(k))
@test length(vec_of_idxs) == length(nidxs)
@test vec_of_idxs[1] == 48:52
@test vec_of_idxs[3] == 52:-1:48
@test sort(vec_of_idxs[1]) == sort(vec_of_idxs[2])
@test vec_of_ds[1] == vec_of_ds[3] == [(i-1)*0.001 for i in 1:5] # also tests sorting

vec_of_idxs, vec_of_ds = bulksearch(tree1, queries, NeighborNumber(k), theiler1)

@test 48 ∉ vec_of_idxs[1]
@test 49 ∉ vec_of_idxs[1]
@test 50 ∈ vec_of_idxs[1]
@test 51 ∉ vec_of_idxs[3]
@test 52 ∉ vec_of_idxs[3]
@test 50 ∈ vec_of_idxs[3]
@test 48 ∈ vec_of_idxs[2]
@test 52 ∈ vec_of_idxs[2]
@test 50 ∉ vec_of_idxs[2]

vec_of_idxs, vec_of_ds = bulksearch(tree1, queries, NeighborNumber(k), theiler2)
@test vec_of_idxs[1] == isearch(tree1, queries[1], NeighborNumber(k), theiler2(48))
for j in 48:52
    @test j ∉ vec_of_idxs[2]
end

_vec_of_idxs = bulkisearch(tree1, queries, WithinRange(0.002))
@test length.(_vec_of_idxs) == [3, 5, 3]

vec_of_idxs, vec_of_ds = bulksearch(tree1, queries, WithinRange(0.002))
for ds in vec_of_ds
    @test issorted(ds)
end
@test vec_of_idxs[1] == [48, 49, 50]
@test sort(vec_of_idxs[2]) == [48, 49, 50, 51, 52]

# final test, theiler window

vec_of_idxs, vec_of_ds = bulksearch(tree1, queries, WithinRange(0.002), theiler1)
@test vec_of_idxs[1] == vec_of_idxs[3] == [50]
@test vec_of_ds[1] == vec_of_ds[3] == [0.002]
@test vec_of_idxs[2] == [48, 52]
@test vec_of_ds[2] == [0.002, 0.002]

vec_of_idxs, vec_of_ds = bulksearch(tree1, queries, WithinRange(0.002), theiler2)
for (x, y) in zip(vec_of_idxs, vec_of_ds)
    @test isempty(x)
    @test isempty(y)
end

end


text:
using TuringModels, StatsFuns

Turing.setadbackend(:reverse_diff);
#nbTuring.turnprogress(false);

μ = 1.4
σ = 1.5
nponds = 60
ni = repeat([5,10,25,35], inner=15);

a_pond = rand(Normal(μ, σ), nponds);

dsim = DataFrame(pond = 1:nponds, ni = ni, true_a = a_pond);

prob = logistic.(Vector{Real}(dsim[:true_a]));

dsim[!, :si] = [rand(Binomial(ni[i], prob[i])) for i = 1:nponds];

# Used only in the continuation of this example
dsim[!, :p_nopool] = dsim[:, :si] ./ dsim[:, :ni];

# Turing model

@model m12_3(pond, si, ni) = begin

    # Separate priors on μ and σ for each pond
    σ ~ Truncated(Cauchy(0, 1), 0, Inf)
    μ ~ Normal(0, 1)

    # Number of ponds in the data set
    N_ponds = length(pond)

    # vector for the priors for each pond
    a_pond = Vector{Real}(undef, N_ponds)

    # For each pond set a prior. Note the [] around Normal(), i.e.,
    a_pond ~ [Normal(μ, σ)]

    # Observation
    logitp = [a_pond[pond[i]] for i = 1:N_ponds]
    si ~ VecBinomialLogit(ni, logitp)

end

# Sample

chns = sample(m12_3(Vector{Int64}(dsim[:, :pond]), Vector{Int64}(dsim[:, :si]),
    Vector{Int64}(dsim[:, :ni])), Turing.NUTS(0.8), 1000);

# Results from rethinking

m123rethinking = "
                  mean   sd  5.5% 94.5% n_eff Rhat
a                 1.30 0.23  0.94  1.67  8064    1
sigma         1.55 0.21  1.24  1.92  3839    1
a_pond[1]   2.57 1.17  0.85  4.57  9688    1
a_pond[2]   2.58 1.19  0.83  4.56  9902    1
a_pond[3]   2.56 1.16  0.84  4.57 12841    1
a_pond[4]   1.49 0.92  0.12  3.03 15532    1
a_pond[5]   1.51 0.95  0.07  3.09 14539    1
a_pond[6]   0.72 0.84 -0.59  2.08 13607    1
a_pond[7]   2.56 1.16  0.86  4.51 12204    1
a_pond[8]   1.50 0.93  0.07  3.05 19903    1
a_pond[9]   2.56 1.15  0.86  4.51 11054    1
a_pond[10]  1.49 0.95  0.05  3.09 14134    1
a_pond[11] -0.64 0.86 -2.06  0.70 15408    1
a_pond[12]  2.56 1.16  0.86  4.53 11512    1
a_pond[13]  1.49 0.95  0.05  3.10 16270    1
a_pond[14]  0.71 0.84 -0.59  2.07 17077    1
a_pond[15]  1.50 0.93  0.10  3.05 16996    1
a_pond[16]  2.98 1.07  1.45  4.84  9033    1
a_pond[17]  2.09 0.84  0.85  3.54 14636    1
a_pond[18]  1.01 0.66  0.00  2.10 12971    1
a_pond[19]  1.01 0.68 -0.03  2.13 12598    1
a_pond[20]  1.48 0.72  0.38  2.67 15500    1
a_pond[21]  2.96 1.09  1.42  4.87 11204    1
a_pond[22] -2.04 0.87 -3.53 -0.75  9065    1
a_pond[23]  0.99 0.67 -0.04  2.11 15365    1
a_pond[24]  1.48 0.72  0.41  2.67 14879    1
a_pond[25]  2.10 0.85  0.85  3.53 13298    1
a_pond[26]  1.00 0.65  0.01  2.06 18583    1
a_pond[27]  3.00 1.08  1.44  4.86  9312    1
a_pond[28]  0.98 0.66 -0.03  2.09 14703    1
a_pond[29]  0.21 0.61 -0.76  1.19 15554    1
a_pond[30]  2.95 1.05  1.45  4.73  9816    1
a_pond[31]  1.70 0.53  0.89  2.59 19148    1
a_pond[32]  0.82 0.42  0.17  1.51 13556    1
a_pond[33]  0.32 0.40 -0.33  0.96 19388    1
a_pond[34] -0.15 0.40 -0.79  0.48 18684    1
a_pond[35]  3.57 0.98  2.19  5.26  8769    1
a_pond[36]  0.16 0.40 -0.46  0.80 17595    1
a_pond[37]  2.00 0.58  1.13  2.99 14669    1
a_pond[38] -1.41 0.49 -2.22 -0.65 12957    1
a_pond[39]  1.21 0.46  0.49  1.97 14185    1
a_pond[40] -1.18 0.46 -1.95 -0.48 16142    1
a_pond[41]  2.86 0.78  1.73  4.18 10508    1
a_pond[42]  0.00 0.39 -0.61  0.63 16138    1
a_pond[43]  1.43 0.48  0.70  2.24 17100    1
a_pond[44]  2.86 0.77  1.75  4.15 12002    1
a_pond[45] -1.40 0.49 -2.21 -0.66 14292    1
a_pond[46]  0.12 0.33 -0.40  0.66 20425    1
a_pond[47] -0.56 0.36 -1.14  0.00 18981    1
a_pond[48]  1.11 0.38  0.52  1.73 14176    1
a_pond[49]  3.81 0.95  2.47  5.45  8841    1
a_pond[50]  2.05 0.50  1.31  2.88 15898    1
a_pond[51] -1.40 0.41 -2.08 -0.76 17188    1
a_pond[52] -0.11 0.34 -0.65  0.43 17158    1
a_pond[53]  1.61 0.44  0.94  2.36 15132    1
a_pond[54]  2.05 0.50  1.30  2.89 15799    1
a_pond[55]  3.14 0.75  2.08  4.40 12702    1
a_pond[56]  3.13 0.74  2.07  4.41 11143    1
a_pond[57]  1.26 0.40  0.65  1.92 14587    1
a_pond[58]  1.11 0.38  0.51  1.74 21740    1
a_pond[59]  2.33 0.56  1.50  3.25 13116    1
a_pond[60]  1.27 0.40  0.66  1.91 15611    1
";

# Draw summary

describe(chns)

# End of `12/m12.3t.jl`


text:
function ylmKCoefficient(l::Integer, m::Integer)
    k = BigInt(4) # BigInt avoids overflow
    for i = l-m+1:l+m
        k *= i
    end

    return sqrt((2.0 * l + 1) / k / π)
end

function ylmCosSinPolynomial(m::Integer, x::Variable, y::Variable)
    sum = 0.0 * (x + y)
    for j = 0:m÷2
        sum += (-1)^j * binomial(m, 2 * j) * (y^(2 * j)) * (x^(m - 2 * j))
    end

    return sum
end

function ylmSinSinPolynomial(m::Integer, x::Variable, y::Variable)
    sum = 0.0 * (x + y)
    for j = 0:(m-1)÷2
        sum += ((-1)^j) * binomial(m, 2 * j + 1) * (y^(2 * j + 1)) * (x^(m - 2 * j - 1))
    end

    return sum
end

"""
    ylm(l, m, x, y, z)

Calculation of the spherical harmonic for a given order (l,m) in Cartesian coordinates

- @arg l: degree of spherical harmonics
- @arg m: order of spherical harmonics
- @arg x, y, z: Cartesian coordinates
- @return: spherical harmonic polynomial

"""
function ylm(l::Integer, m::Integer, x::Variable, y::Variable, z::Variable)
    if abs(m) > l
        throw("-l <= m <= l expected, but m = $m and l = $l.")
    end

    p = 1.0 * (z^2 - 1)^l + 0.0 * (x + y)

    for i = 1:l+abs(m)
        c = i <= l ? 1 / (2 * i) : 1.0
        p = c * differentiate(p, z)
    end

    if m > 0
        return sqrt(2) * ylmKCoefficient(l, m) * ylmCosSinPolynomial(m, x, y) * p
    elseif m < 0
        return sqrt(2) * ylmKCoefficient(l, abs(m)) * ylmSinSinPolynomial(abs(m), x, y) * p
    else
        return ylmKCoefficient(l, 0) * p
    end
end

"""
    rlylm(l, m, x, y, z)

Calculate `r^l * Ylm(x,y,z)``
"""
function rlylm(l::Integer, m::Integer, x::Variable, y::Variable, z::Variable)
    p = ylm(l, m, x, y, z)
    tout = []

    for t in terms(p)
        deg = degree(monomial(t))
        degR = l - deg
        push!(tout, (x^2 + y^2 + z^2)^Int(degR / 2) * t)
    end

    return polynomial(tout)
end

"""
    eval_spherharmonic(points, L)

Evaluate spherical harmonics basis at given quadrature points
"""
function eval_spherharmonic(points::AbstractMatrix{T}, L) where {T<:Real}
    ne = (L + 1)^2
    nq = size(points, 1)
    m = zeros(ne, nq)

    @polyvar x y z
    if L == 0
        spe = [ylm(0, 0, x, y, z)]
    elseif L == 1
        spe = [
            ylm(0, 0, x, y, z),
            ylm(1, -1, x, y, z),
            ylm(1, 0, x, y, z),
            ylm(1, 1, x, y, z),
        ]
    elseif L == 2
        spe = [
            ylm(0, 0, x, y, z),
            ylm(1, -1, x, y, z),
            ylm(1, 0, x, y, z),
            ylm(1, 1, x, y, z),
            ylm(2, -2, x, y, z),
            ylm(2, -1, x, y, z),
            ylm(2, 0, x, y, z),
            ylm(2, 1, x, y, z),
            ylm(2, 2, x, y, z),
        ]
    end

    for j in axes(m, 2), i in axes(m, 1)
        m[i, j] = spe[i](x => points[j, 1], y => points[j, 2], z => points[j, 3])
    end

    return m
end


text:
using StanSample
path = @__DIR__
println(path)
gp_stan = open(path * "/gp_stan.txt") do file
    read(file, String)
end

sm = SampleModel("gaussian_process", gp_stan)

#=
function run_stan(sm, df, strain, tc)
    sub_df = df[(df.tc .== tc) .& (df.strain .== strain) , :]
    println(sub_df)
    analyzed_df = DataFrame(strain=String[], well=String[], tc=Float64[], λ_max=Float64[])

    for (well, date) in collect(Iterators.product(sub_df.well |> unique, sub_df.date |> unique))
        for_df = sub_df[(sub_df.tc .== tc) .& (sub_df.strain .== strain) , :]
        sort!(for_df, :time_min)

        x = for_df[(for_df.well .== well) , "time_min"]
        y = for_df[(for_df.well .== well), "OD600"]

        data = Dict(
            "N" => length(x),
            "t"=> x,
            "y"=> y,
            "N_predict"=> length(x),
            "t_predict"=> x,
            "alpha_param"=> [0, 1],
            "sigma_param"=> [0, 1],
            "rho_param"=> [1000, 50],
        )

        rc = stan_sample(sm; data, num_chains=2);
        st = read_samples(sm)
        dy = filter(x -> occursin("dy",String(x)), names(st))
        mean_dy = [mean(st[x]) for x in dy]
        _df = DataFrame(
            strain=String[strain],
            well=String[well],
            tc=Float64[tc],
            λ_max=Float64[maximum(mean_dy)],
        )
        append!(analyzed_df, _df)
    end
    analyzed_df
end

#analyzed_df = DataFrame(strain=String[], well=String[], tc=Float64[], λ_max=Float64[])#, dy=[])

sm = SampleModel("gaussian_process", gp_stan)
#run_stan(sm, df, "MG1655", 0)=#


text:
doc"""
    BetaPrime(α,β)

The *Beta prime distribution* has probability density function

$f(x; \alpha, \beta) = \frac{1}{B(\alpha, \beta)}
x^{\alpha - 1} (1 + x)^{- (\alpha + \beta)}, \quad x > 0$

The Beta prime distribution is related to the [`Beta`](:func:`Beta`) distribution via the
relation ship that if $X \sim \operatorname{Beta}(\alpha, \beta)$ then $\frac{X}{1 - X}
\sim \operatorname{BetaPrime}(\alpha, \beta)$

```julia
BetaPrime()        # equivalent to BetaPrime(0.0, 1.0)
BetaPrime(a)       # equivalent to BetaPrime(a, a)
BetaPrime(a, b)    # Beta prime distribution with shape parameters a and b

params(d)          # Get the parameters, i.e. (a, b)
```

External links

* [Beta prime distribution on Wikipedia](http://en.wikipedia.org/wiki/Beta_prime_distribution)

    """
immutable BetaPrime <: ContinuousUnivariateDistribution
    α::Float64
    β::Float64

    function BetaPrime(α::Real, β::Real)
        @check_args(BetaPrime, α > zero(α) && β > zero(β))
        new(α, β)
    end
    BetaPrime(α::Real) = BetaPrime(α,α)
    BetaPrime() = new(1.0, 1.0)
end

@distr_support BetaPrime 0.0 Inf

#### Parameters

params(d::BetaPrime) = (d.α, d.β)

#### Statistics

mean(d::BetaPrime) = ((α, β) = params(d); β > 1.0 ? α / (β - 1.0) : NaN)

mode(d::BetaPrime) = ((α, β) = params(d); α > 1.0 ? (α - 1.0) / (β + 1.0) : 0.0)

function var(d::BetaPrime)
    (α, β) = params(d)
    β > 2.0 ? α * (α + β - 1.0) / ((β - 2.0) * (β - 1.0)^2) : NaN
end

function skewness(d::BetaPrime)
    (α, β) = params(d)
    if β > 3.0
        s = α + β - 1.0
        2.0 * (α + s) / (β - 3.0) * sqrt((β - 2.0) / (α * s))
    else
        return NaN
    end
end

#### Evaluation

function logpdf(d::BetaPrime, x::Float64)
    (α, β) = params(d)
    (α - 1.0) * log(x) - (α + β) * log1p(x) - lbeta(α, β)
end

pdf(d::BetaPrime, x::Float64) = exp(logpdf(d, x))

cdf(d::BetaPrime, x::Float64) = betacdf(d.α, d.β, x / (1.0 + x))
ccdf(d::BetaPrime, x::Float64) = betaccdf(d.α, d.β, x / (1.0 + x))
logcdf(d::BetaPrime, x::Float64) = betalogcdf(d.α, d.β, x / (1.0 + x))
logccdf(d::BetaPrime, x::Float64) = betalogccdf(d.α, d.β, x / (1.0 + x))

quantile(d::BetaPrime, p::Float64) = (x = betainvcdf(d.α, d.β, p); x / (1.0 - x))
cquantile(d::BetaPrime, p::Float64) = (x = betainvccdf(d.α, d.β, p); x / (1.0 - x))
invlogcdf(d::BetaPrime, p::Float64) = (x = betainvlogcdf(d.α, d.β, p); x / (1.0 - x))
invlogccdf(d::BetaPrime, p::Float64) = (x = betainvlogccdf(d.α, d.β, p); x / (1.0 - x))

#### Sampling

function rand(d::BetaPrime)
    (α, β) = params(d)
    rand(Gamma(α)) / rand(Gamma(β))
end


text:
"""
    function check_par(par)
---
## Description
Check if enough and right.
"""
function check_par(par)
    @info "Check if the simulation parameters are valid"
    must_keys = Dict(
        :nSire => "Number of sires",
        :nDam  => "Number of dams",
        :nSib  => "Number of sibs",
        :nC7e  => "Number of sibs for challenge(c7e)",
        :nG8n  => "Number of generation(g8n) to breed",
        :h²    => "Heritability for traits",
        :w4t   => "Weight(w4t) for the binary trait",
        :cv    => "Covariance of pleiotropic QTL effects for 2 traits",
        :p17h  => "Percentage-of-death(p17h) after challenge",
        :e19e  => "edit-successful-rate(e19e), [0, 1]",
        :nK3n  => "Number of known(k3n) QTL of the binary trait",
        :dd    => "An ϵ to be added to diag(LHS matrix) to avoid singular",
        :edit  => "Whether to edit known QTL",
        :fix   => "Whether to emphasize known QTL",
        :b4y   => "Whether the binary(b4y) trait to be recorded as binary",
        :u38y  => "use-genetypes-of-current-generation-only(u38y)"
    )
    missingkey = []
    for key in keys(must_keys)
        haskey(par, key) || push!(missingkey, key)
    end
    isempty(missingkey) || begin
        for key in missingkey
            println("Key $key \"$(must_keys[key])\" is missing")
        end
        error("Some must-have key(s) is/are missing")
    end
    haskey(par, :b4y) && !par.b4y &&
        @warn "Continuous disease phenotype was for debugging" par.b4y
    # More to be added
end


text:
module JsonDB

export Database
import Base.Threads.@spawn
import JSON3

"""
`JsonDB.Database`

This type holds instance of database
"""
mutable struct Database
	db::AbstractDict
	thread
	path::String
	auto_dump::Bool
	Database(path::String, auto_dump::Bool) = new(Dict(), nothing, path, auto_dump)
end

"""
`JsonDB.load(path::String, auto_dump::Bool)`

Returns a Database object which is stored at `path`. Flag `auto_dump` determines whether to dump the databse to disk after every set(write), update and delete operation.
"""
function load(path::String, auto_dump::Bool, sig=true)
	db = Database(path, auto_dump)
	_load(db, path, auto_dump)
	if sig
		set_sigterm_handler(db::Database)
	end
	return db
end

function set_sigterm_handler(db::Database)
	item = db
	function sigterm_handler()
	    if item.thread !== nothing
		     fetch(item.thread)
		exit(Base.SIGTERM)
	    end
	end
	atexit(sigterm_handler)
end

function _load(db::Database, path, auto_dump)
	path = expanduser(path)
	db.path = path
	db.auto_dump = auto_dump
	if ispath(path)
		_loaddb(db)
	else
		db.db = Dict()
	end
	return true
end

function _dump(db::Database)
	json = JSON3.write(db.db)
	open(db.path, "w") do file
		write(file, json)
	end
end

"""
`JsonDB.dump(db::Database)`

Force dump databse to memory
"""
function dump(db::Database)
	db.thread = @spawn _dump(db)
	fetch(db.thread)
	return true
end

function _loaddb(db::Database)
		file = open(db.path, "r")
	try
		db.db = JSON3.read(file, Dict)
	catch err
		if filesize(file) == 0
			db.db = Dict()
		elseif isa(err, JSON3.Error)
			throw(err)
		else
			error("Failed to read JSON")
		end
	end
end

function _autodump(db::Database)
	if db.auto_dump
		dump(db)
	end
end

"""
`JsonDB.set(db::Database, key, value)`

Insert key-value pair into the database
"""
function set(db, key, value)
	if isa(key, String)
		db.db[key] = value
		_autodump(db)
		return true
	else
		throw(error("key/name must be string"))
	end
end

"""
`JsonDB.get(db::Database, key)`

Get value of the `key`
"""
function get(db::Database, key)
	Base.get(db.db, key, false)
end

"""
`JsonDB.getall(db::Database)`

Returns iterator over all keys in the database
"""
function getall(db::Database)
	return Base.keys(db.db)
end

"""
`Base.haskey(db::Database, key)`

Returns `true` if `key` is present in the database else returns `false`
"""
function Base.haskey(db::Database, key)
	return Base.haskey(db.db, key)
end

"""
`JsonDB.delete(db::Database, key)`

Deletes `key` from the database
"""
function delete(db::Database, key)
	if !Base.haskey(db.db, key)
		return false
	end
	Base.delete!(db.db, key)
	_autodump(db)
	return true
end

"""
`Base.length(db::Database, name=nothing)`

Returns total number of key/value pairs in the database if `name` is `nothing` else returns length of the value with key `name`
"""
function Base.length(db::Database, name=nothing)
	if isnothing(name)
		total = length(db.db)
		return total
	else
		total = length(db.db[name])
		return total
	end
end

"""
`JsonDB.append(db::Database, key, more)`

Adds `more` to `keys`'s value
"""
function append(db::Database, key, more)
	tmp = db.db[key]
	db.db[key] = tmp * more
	_autodump(db)
	return true
end

"""
`JsonDB.create_array(db::Database, name::String)`

Create an array as key with label `name`
"""
function create_array(db::Database, name::String)
	if isa(name, String)
		db.db[name] = []
		_autodump(db)
		return true
	else
		throw(err("key/name must be string"))
	end
end

"""
`JsonDB.set_array(db::Database, name::String, value)`

Adds `value` to array with label `name`
"""
function set_array(db::Database, name::String, value)
	db.db[name] = vcat(db.db[name], value)
	_autodump(db)
	return true
end

"""
`JsonDB.append_array(db::Database, name::String, seq)`

Appends `more` to the array with label `name`
"""
function append_array(db::Database, name::String, seq)
	Base.push!(db.db[name], seq)
	_autodump(db)
	return true
end

"""
`JsonDB.get(db::Database, name::String, pos::Integer)`

Returns value at `pos` from an array with label `name`
"""
function get(db::Database, name::String, pos::Integer)
	return db.db[name][pos]
end

"""
JsonDB.range(db::Database, name::String, _start, _end)

Returns slice from array with label `name` with starting index `_start` and ending index `_end`
"""
function range(db::Database, name::String, _start, _end)
	return db.db[name][_start:_end]
end

"""
`JsonDB.delete_array(db::Database, name::String)`

Deletes array with key/label `name` from the database
"""
function delete_array(db, name::String)
	len = Base.length(db.db[name])
	Base.delete!(db.db , name)
	_autodump(db)
	return len
end

"""
`JsonDB.delete_value(db::Database, name::String, value)`

Deletes `value` from array with label `name`
"""
function delete_value(db, name::String, value)
	Base.delete!(db.db[name], value)
	_autodump(db)
	return true
end

"""
`JsonDB.pop(db::Database, name::String, pos)`

In case of array deletes the value at index `pos`from array with label `name`
In case of dictionary deletes the  value with key `pos` from dictionary with label `name`
"""
JsonDB.pop(db::Database, name, pos)
function pop(db, name, pos)
	value = pop(db.db[name], pos)
	_autodump(db)
	return value
end
"""
`JsonDB.append(db::Database, name::String, pos::Integer, more)`

Appends `more` to value at index `pos` in array with label `name`
"""
function append(db::Database, name::String, pos::Integer, more)
	tmp = db.db[name][pos]
	db.db[name][pos] = tmp + more
	_autodump(db)
	return true
end

"""
`JsonDB.hasvalue(db::Database, name::String, value)`

Returns `true` if `value` is present in array with label `name` else returns `false`
"""
function hasvalue(db::Database, name::String, value)
	return value in db.db[name]
end

"""
`JsonDB.create_dict(db::Database, name::String)`

Creates a dictionary with label `name`
"""
function create_dict(db::Database, name::String)
	if isa(name, String)
		db.db[name] = Dict()
		_autodump(db)
		return true;
	else
		throw(err("key/name must be string"))
	end
end

"""
`JsonDB.set_dict(db::Database, name::String, pair)`

Adds a key/value pair to dictionary with label `name`. `pair` can be an array or tuple of length 2
"""
function set_dict(db::Database, name::String, pair)
	db.db[name][pair[1]] = pair[2]
	JsonDB._autodump(db)
	return true
end

"""
`JsonDB.delete_dict(db::Database, name::String)`

Deletes dictionary with label `name` and all of its keys
"""
function delete_dict(db::Database, name::String)
	Base.delete!(db.db, name)
	_autodump(db)
	return true
end

"""
`JsonDB.haskey(db::Database, name::String, key)`

Checks whether `key` is in dictionary with label `name`
"""
function haskey(db::Database, name::String, key)
	return key in Base.collect(Base.keys(db.db[name]))
end

"""
`JsonDB.keys(db::Database, name::String)`

Returns iterator over all the keys in dictionary with label `name`
"""
function keys(db::Database, name::String)
	return Base.keys(self.db[name])
end

"""
`JsonDB.values(db::Database, name::String)

Returns iterator over all the values in dictionary with label `name`
"""
function values(db::Database, name::String)
	return Base.values(self.db[name])
end

"""
`JsonDB.merge(db::Database, name1::String, name2::String)`

Merge the dictionary with label `name2` into the dictionary with label `name1`.
See the documentation for `[Base.merge]`(https://docs.julialang.org/en/v1/base/collections/#Base.merge)
"""
function merge(db::Database, name1::String, name2::String)
	merge!(db.db[name1], db.db[name2])
	_autodump(db)
	return true
end

"""
`JsonDB.deletedb(db::Database)`

Deletes everyhing from the database `db`
"""
function deletedb(db::Database)
	db.db = Dict()
	_autodump(db)
	return true
end

end