text:
function _readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = _context)
    result = GEOSWKTReader_read_r(context.ptr, wktreader.ptr, pointer(wktstring))
    if result == C_NULL
        error("LibGEOS: Error in GEOSWKTReader_read_r while reading $wktstring")
    end
    result
end
_readgeom(wktstring::String, context::GEOSContext = _context) =
    _readgeom(wktstring, WKTReader(context), context)

function _readgeom(
    wkbbuffer::Vector{Cuchar},
    wkbreader::WKBReader,
    context::GEOSContext = _context,
)
    result = GEOSWKBReader_read_r(
        context.ptr,
        wkbreader.ptr,
        pointer(wkbbuffer),
        length(wkbbuffer),
    )
    if result == C_NULL
        error("LibGEOS: Error in GEOSWKBReader_read_r")
    end
    result
end
_readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = _context) =
    _readgeom(wkbbuffer, WKBReader(context), context)

function _writegeom(geom::GEOSGeom, wktwriter::WKTWriter, context::GEOSContext = _context)
    GEOSWKTWriter_write_r(context.ptr, wktwriter.ptr, geom)
end

function _writegeom(geom::GEOSGeom, wkbwriter::WKBWriter, context::GEOSContext = _context)
    wkbsize = Ref{Csize_t}()
    result = GEOSWKBWriter_write_r(context.ptr, wkbwriter.ptr, geom, wkbsize)
    unsafe_wrap(Array, result, wkbsize[], own = true)
end

_writegeom(geom::GEOSGeom, context::GEOSContext = _context) =
    _writegeom(geom, WKTWriter(context), context)

# -----
# Coordinate Sequence functions
# -----
"""
    createCoordSeq(size::Integer; ndim::Integer=2) -> Ptr{Ptr{Void}}

Create a Coordinate sequence with ``size'' coordinates of ``dims'' dimensions (Return NULL on exception)
"""
function createCoordSeq(size::Integer, context::GEOSContext = _context; ndim::Integer = 2)
    @assert ndim >= 2
    result = GEOSCoordSeq_create_r(context.ptr, size, ndim)
    if result == C_NULL
        error("LibGEOS: Error in GEOSCoordSeq_create")
    end
    result
end

# Clone a Coordinate Sequence (Return NULL on exception)
function cloneCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSCoordSeq_clone_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSCoordSeq_clone")
    end
    result
end

function destroyCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSCoordSeq_destroy_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSCoordSeq_destroy")
    end
    result
end

# Set ordinate values in a Coordinate Sequence (Return 0 on exception)
function setX!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
    result = GEOSCoordSeq_setX_r(context.ptr, ptr, i - 1, value)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_setX")
    end
    result
end

function setY!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
    result = GEOSCoordSeq_setY_r(context.ptr, ptr, i - 1, value)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_setY")
    end
    result
end

function setZ!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
    result = GEOSCoordSeq_setZ_r(context.ptr, ptr, i - 1, value)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_setZ")
    end
    result
end

"Get size info from a Coordinate Sequence (Return 0 on exception)"
function getSize(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    out = Ref{UInt32}()
    result = GEOSCoordSeq_getSize_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getSize")
    end
    Int(out[])
end

"Get dimensions info from a Coordinate Sequence (Return 0 on exception)"
function getDimensions(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    out = Ref{UInt32}()
    result = GEOSCoordSeq_getDimensions_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getDimensions")
    end
    Int(out[])
end

# convenience functions
function setCoordSeq!(
    ptr::GEOSCoordSeq,
    i::Integer,
    coords::Vector{Float64},
    context::GEOSContext = _context,
)
    ndim = length(coords)
    @assert ndim >= 2
    setX!(ptr, i, coords[1], context)
    setY!(ptr, i, coords[2], context)
    ndim >= 3 && setZ!(ptr, i, coords[3], context)
    ptr
end

"""
    createCoordSeq(x::Real, y::Real) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a single 2D coordinate
"""
function createCoordSeq(x::Real, y::Real, context::GEOSContext = _context)
    coordinates = createCoordSeq(1, context, ndim = 2)
    setX!(coordinates, 1, x, context)
    setY!(coordinates, 1, y, context)
    coordinates
end

"""
    createCoordSeq(x::Real, y::Real, z::Real) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a single 3D coordinate
"""
function createCoordSeq(x::Real, y::Real, z::Real, context::GEOSContext = _context)
    coordinates = createCoordSeq(1, context, ndim = 3)
    setX!(coordinates, 1, x, context)
    setY!(coordinates, 1, y, context)
    setZ!(coordinates, 1, z, context)
    coordinates
end

"""
    createCoordSeq(coords::Vector{Float64}) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a single N dimensional coordinate
"""
function createCoordSeq(coords::Vector{Float64}, context::GEOSContext = _context)
    ndim = length(coords)
    @assert ndim >= 2
    coordinates = createCoordSeq(1, context, ndim = ndim)
    setCoordSeq!(coordinates, 1, coords, context)
end

"""
    createCoordSeq(coords::Vector{Float64}) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a multiple N dimensional coordinate
"""
function createCoordSeq(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
    ncoords = length(coords)
    @assert ncoords > 0
    ndim = length(coords[1])
    coordinates = createCoordSeq(ncoords, context, ndim = ndim)
    for (i, coord) in enumerate(coords)
        setCoordSeq!(coordinates, i, coord, context)
    end
    coordinates
end

function getX(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSCoordSeq_getX_r(context.ptr, ptr, i - 1, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getX")
    end
    out[]
end

function getX(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ncoords = getSize(ptr, context)
    xcoords = Array{Float64}(ncoords)
    start = pointer(xcoords)
    floatsize = sizeof(Float64)
    for i = 0:ncoords-1
        GEOSCoordSeq_getX_r(context.ptr, ptr, i, start + i * floatsize)
    end
    xcoords
end

function getY(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSCoordSeq_getY_r(context.ptr, ptr, i - 1, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getY")
    end
    out[]
end

function getY(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ncoords = getSize(ptr, context)
    ycoords = Array{Float64}(ncoords)
    start = pointer(ycoords)
    floatsize = sizeof(Float64)
    for i = 0:ncoords-1
        GEOSCoordSeq_getY_r(context.ptr, ptr, i, start + i * floatsize)
    end
    ycoords
end

function getZ(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSCoordSeq_getZ_r(context.ptr, ptr, i - 1, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getZ")
    end
    out[]
end

function getZ(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ncoords = getSize(ptr, context)
    zcoords = Array{Float64}(ncoords)
    start = pointer(zcoords)
    floatsize = sizeof(Float64)
    for i = 0:ncoords-1
        GEOSCoordSeq_getZ_r(context.ptr, ptr, i, start + i * floatsize)
    end
    zcoords
end

function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    ndim = getDimensions(ptr, context)
    coord = Array{Float64}(undef, ndim)
    start = pointer(coord)
    floatsize = sizeof(Float64)
    GEOSCoordSeq_getX_r(context.ptr, ptr, i - 1, start)
    GEOSCoordSeq_getY_r(context.ptr, ptr, i - 1, start + floatsize)
    if ndim == 3
        GEOSCoordSeq_getZ_r(context.ptr, ptr, i - 1, start + 2 * floatsize)
    end
    coord
end

function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ndim = getDimensions(ptr, context)
    ncoords = getSize(ptr, context)
    coordseq = Vector{Float64}[]
    sizehint!(coordseq, ncoords)
    for i = 1:ncoords
        push!(coordseq, getCoordinates(ptr, i, context))
    end
    coordseq
end

function isCCW(ptr::GEOSCoordSeq, context::GEOSContext=_context)::Bool
    d = UInt8[1]
    GC.@preserve result = GEOSCoordSeq_isCCW_r(context.ptr, ptr, pointer(d))
    if result == C_NULL
        error("LibGEOS: Error in GEOSInterpolateNormalized")
    end
    Bool(d[1])
end

# -----
# Linear referencing functions -- there are more, but these are probably sufficient for most purposes
# -----
# (GEOSGeometry ownership is retained by caller)

# Return distance of point 'p' projected on 'g' from origin of 'g'. Geometry 'g' must be a lineal geometry
project(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = _context) =
    GEOSProject_r(context.ptr, g, p)

# Return closest point to given distance within geometry (Geometry must be a LineString)
function interpolate(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
    result = GEOSInterpolate_r(context.ptr, ptr, d)
    if result == C_NULL
        error("LibGEOS: Error in GEOSInterpolate")
    end
    result
end

projectNormalized(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = _context) =
    GEOSProjectNormalized_r(context.ptr, g, p)

function interpolateNormalized(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
    result = GEOSInterpolateNormalized_r(context.ptr, ptr, d)
    if result == C_NULL
        error("LibGEOS: Error in GEOSInterpolateNormalized")
    end
    result
end

# -----
# Buffer related functions
# -----
# Computes the buffer of a geometry, for both positive and negative buffer distances.
# Since true buffer curves may contain circular arcs, computed buffer polygons can only be approximations to the true geometry.
# The user can control the accuracy of the curve approximation by specifying the number of linear segments with which to approximate a curve.

# Always returns a polygon. The negative or zero-distance buffer of lines and points is always an empty Polygon.
buffer(ptr::GEOSGeom, width::Real, quadsegs::Integer = 8, context::GEOSContext = _context) =
    GEOSBuffer_r(context.ptr, ptr, width, Int32(quadsegs))

bufferWithStyle(
    ptr::GEOSGeom,
    width::Real,
    quadsegs::Integer = 8,
    endCapStyle::GEOSBufCapStyles = GEOSBUF_CAP_ROUND,
    joinStyle::GEOSBufJoinStyles = GEOSBUF_JOIN_ROUND,
    mitreLimit::Real = 5.0,
    context::GEOSContext = _context,
) = GEOSBufferWithStyle_r(
    context.ptr,
    ptr,
    width,
    Int32(quadsegs),
    Int32(endCapStyle),
    Int32(joinStyle),
    mitreLimit,
)

# GEOSBufferParams_create
# GEOSBufferParams_destroy
# GEOSBufferParams_setEndCapStyle
# GEOSBufferParams_setJoinStyle
# GEOSBufferParams_setMitreLimit
# GEOSBufferParams_setQuadrantSegments
# GEOSBufferParams_setSingleSided
# GEOSBufferWithParams
# GEOSOffsetCurve

# -----
# Geometry Constructors -- All functions return NULL on exception
# -----
# (GEOSCoordSequence* arguments will become ownership of the returned object.)

function createPoint(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSGeom_createPoint_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createPoint")
    end
    result
end
createPoint(x::Real, y::Real, context::GEOSContext = _context) =
    createPoint(createCoordSeq(x, y, context), context)
createPoint(x::Real, y::Real, z::Real, context::GEOSContext = _context) =
    createPoint(createCoordSeq(x, y, z, context), context)
createPoint(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
    createPoint(createCoordSeq(coords, context), context)
createPoint(coords::Vector{Float64}, context::GEOSContext = _context) =
    createPoint(createCoordSeq(Vector{Float64}[coords], context), context)

function createLinearRing(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSGeom_createLinearRing_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createLinearRing")
    end
    result
end
createLinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
    GEOSGeom_createLinearRing_r(context.ptr, createCoordSeq(coords, context))

function createLineString(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSGeom_createLineString_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createLineString")
    end
    result
end
createLineString(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
    GEOSGeom_createLineString_r(context.ptr, createCoordSeq(coords, context))

# Second argument is an array of GEOSGeometry* objects.
# The caller remains owner of the array, but pointed-to
# objects become ownership of the returned GEOSGeometry.
function createPolygon(
    shell::GEOSGeom,
    holes::Vector{GEOSGeom},
    context::GEOSContext = _context,
)
    result = GEOSGeom_createPolygon_r(context.ptr, shell, pointer(holes), length(holes))
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createPolygon")
    end
    result
end

function createCollection(
    geomtype::GEOSGeomTypes,
    geoms::Vector{GEOSGeom},
    context::GEOSContext = _context,
)
    result =
        GEOSGeom_createCollection_r(context.ptr, geomtype, pointer(geoms), length(geoms))
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createCollection")
    end
    result
end

function createEmptyCollection(geomtype::GEOSGeomTypes, context::GEOSContext = _context)
    result = GEOSGeom_createEmptyCollection_r(context.ptr, geomtype)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createEmptyCollection")
    end
    result
end

# Memory management
function cloneGeom(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeom_clone_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_clone")
    end
    result
end

destroyGeom(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom_destroy_r(context.ptr, ptr)

# -----
# Topology operations - return NULL on exception.
# -----

function envelope(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSEnvelope_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSEnvelope")
    end
    result
end

"""
    minimumRotatedRectangle(geom)

Returns the minimum rotated rectangular POLYGON which encloses the input geometry. The rectangle
has width equal to the minimum diameter, and a longer length. If the convex hill of the input is
degenerate (a line or point) a LINESTRING or POINT is returned. The minimum rotated rectangle can
be used as an extremely generalized representation for the given geometry.
"""
function minimumRotatedRectangle(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSMinimumRotatedRectangle_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSMinimumRotatedRectangle")
    end
    result
end

function intersection(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSIntersection_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSIntersection")
    end
    result
end

# Returns a Geometry that represents the convex hull of the input geometry. The returned geometry contains the minimal number of points needed to represent the convex hull. In particular, no more than two consecutive points will be collinear.

# Returns:
# if the convex hull contains 3 or more points, a Polygon; 2 points, a LineString; 1 point, a Point; 0 points, an empty GeometryCollection.
function convexhull(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSConvexHull_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSConvexHull")
    end
    result
end

function difference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSDifference_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSDifference")
    end
    result
end

function symmetricDifference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSSymDifference_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSymDifference")
    end
    result
end

function boundary(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSBoundary_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSBoundary")
    end
    result
end

function union(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSUnion_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSUnion")
    end
    result
end

function unaryUnion(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSUnaryUnion_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSUnaryUnion")
    end
    result
end

function pointOnSurface(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPointOnSurface_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSPointOnSurface")
    end
    result
end

function centroid(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetCentroid_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetCentroid")
    end
    result
end

function node(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSNode_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSNode")
    end
    result
end

# all arguments remain ownership of the caller (both Geometries and pointers)
function polygonize(geoms::Vector{GEOSGeom}, context::GEOSContext = _context)
    result = GEOSPolygonize_r(context.ptr, pointer(geoms), length(geoms))
    if result == C_NULL
        error("LibGEOS: Error in GEOSPolygonize")
    end
    result
end
# GEOSPolygonizer_getCutEdges
# GEOSPolygonize_full

function lineMerge(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSLineMerge_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSLineMerge")
    end
    result
end

function simplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSSimplify_r(context.ptr, ptr, tol)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSimplify")
    end
    result
end

function topologyPreserveSimplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSTopologyPreserveSimplify_r(context.ptr, ptr, tol)
    if result == C_NULL
        error("LibGEOS: Error in GEOSTopologyPreserveSimplify")
    end
    result
end

# Return all distinct vertices of input geometry as a MULTIPOINT.
# (Note that only 2 dimensions of the vertices are considered when testing for equality)
function uniquePoints(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeom_extractUniquePoints_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_extractUniquePoints")
    end
    result
end

# Find paths shared between the two given lineal geometries.
# Returns a GEOMETRYCOLLECTION having two elements:
# - first element is a MULTILINESTRING containing shared paths
#    having the _same_ direction on both inputs
#  - second element is a MULTILINESTRING containing shared paths
#    having the _opposite_ direction on the two inputs
# (Returns NULL on exception)
function sharedPaths(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSSharedPaths_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSharedPaths")
    end
    result
end

# Snap first geometry on to second with given tolerance
# (Returns a newly allocated geometry, or NULL on exception)
function snap(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSSnap_r(context.ptr, g1, g2, tol)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSnap")
    end
    result
end

# Return a Delaunay triangulation of the vertex of the given geometry
#
# @param g the input geometry whose vertex will be used as "sites"
# @param tolerance optional snapping tolerance to use for improved robustness
# @param onlyEdges if non-zero will return a MULTILINESTRING, otherwise it will
#                  return a GEOMETRYCOLLECTION containing triangular POLYGONs.
#
# @return  a newly allocated geometry, or NULL on exception
function delaunayTriangulation(
    ptr::GEOSGeom,
    tol::Real = 0.0,
    onlyEdges::Bool = false,
    context::GEOSContext = _context,
)
    result = GEOSDelaunayTriangulation_r(context.ptr, ptr, tol, Int32(onlyEdges))
    if result == C_NULL
        error("LibGEOS: Error in GEOSDelaunayTriangulation")
    end
    result
end

function constrainedDelaunayTriangulation(
    ptr::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSConstrainedDelaunayTriangulation_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSConstrainedDelaunayTriangulation")
    end
    result
end

# -----
# Binary predicates - return 2 on exception, 1 on true, 0 on false
# -----
function disjoint(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSDisjoint_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSDisjoint")
    end
    result != 0x00
end

function touches(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSTouches_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSTouches")
    end
    result != 0x00
end

function intersects(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSIntersects_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSIntersects")
    end
    result != 0x00
end

function crosses(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSCrosses_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSCrosses")
    end
    result != 0x00
end

function within(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSWithin_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSWithin")
    end
    result != 0x00
end

function Base.contains(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSContains_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSContains")
    end
    result != 0x00
end

function overlaps(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSOverlaps_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSOverlaps")
    end
    result != 0x00
end

function equals(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSEquals_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSEquals")
    end
    result != 0x00
end

function equalsexact(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSEqualsExact_r(context.ptr, g1, g2, tol)
    if result == 0x02
        error("LibGEOS: Error in GEOSEqualsExact")
    end
    result != 0x00
end

function covers(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSCovers_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSCovers")
    end
    result != 0x00
end

function coveredby(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSCoveredBy_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSCoveredBy")
    end
    result != 0x00
end

# -----
# Prepared Geometry Binary predicates - return 2 on exception, 1 on true, 0 on false
# -----

# GEOSGeometry ownership is retained by caller
function prepareGeom(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPrepare_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSPrepare")
    end
    result
end

destroyPreparedGeom(ptr::Ptr{GEOSPreparedGeometry}, context::GEOSContext = _context) =
    GEOSPreparedGeom_destroy_r(context.ptr, ptr)

function prepcontains(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedContains_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedContains")
    end
    result != 0x00
end

function prepcontainsproperly(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedContainsProperly_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedContainsProperly")
    end
    result != 0x00
end

function prepcoveredby(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedCoveredBy_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedCoveredBy")
    end
    result != 0x00
end

function prepcovers(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedCovers_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedCovers")
    end
    result != 0x00
end

function prepcrosses(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedCrosses_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedCrosses")
    end
    result != 0x00
end

function prepdisjoint(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedDisjoint_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedDisjoint")
    end
    result != 0x00
end

function prepintersects(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedIntersects_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedIntersects")
    end
    result != 0x00
end

function prepoverlaps(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedOverlaps_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedOverlaps")
    end
    result != 0x00
end

function preptouches(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedTouches_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedTouches")
    end
    result != 0x00
end

function prepwithin(
    g1::Ptr{GEOSPreparedGeometry},
    g2::GEOSGeom,
    context::GEOSContext = _context,
)
    result = GEOSPreparedWithin_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedWithin")
    end
    result != 0x00
end

# -----
# STRtree functions
# -----
# GEOSSTRtree_create
# GEOSSTRtree_insert
# GEOSSTRtree_query
# GEOSSTRtree_iterate
# GEOSSTRtree_remove
# GEOSSTRtree_destroy

# -----
# Unary predicate - return 2 on exception, 1 on true, 0 on false
# -----
function isEmpty(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisEmpty_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisEmpty")
    end
    result != 0x00
end

function isSimple(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisSimple_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisSimple")
    end
    result != 0x00
end

function isRing(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisRing_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisRing")
    end
    result != 0x00
end

function hasZ(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSHasZ_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSHasZ")
    end
    result != 0x00
end

# Call only on LINESTRING (return 2 on exception, 1 on true, 0 on false)
function isClosed(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisClosed_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisClosed")
    end
    result != 0x00
end

# -----
# Dimensionally Extended 9 Intersection Model related
# -----

# GEOSRelatePattern (return 2 on exception, 1 on true, 0 on false)
# GEOSRelate (return NULL on exception, a string to GEOSFree otherwise)
# GEOSRelatePatternMatch (return 2 on exception, 1 on true, 0 on false)
# GEOSRelateBoundaryNodeRule (return NULL on exception, a string to GEOSFree otherwise)

# -----
# Validity checking -- return 2 on exception, 1 on true, 0 on false
# -----

# /* These are for use with GEOSisValidDetail (flags param) */
# enum GEOSValidFlags {
#     GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE=1
# };

function isValid(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisValid_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisValid")
    end
    result != 0x00
end

# * return NULL on exception, a string to GEOSFree otherwise
# GEOSisValidReason

# * Caller has the responsibility to destroy 'reason' (GEOSFree)
# * and 'location' (GEOSGeom_destroy) params
# * return 2 on exception, 1 when valid, 0 when invalid
# GEOSisValidDetail

# -----
# Geometry info
# -----

# Return NULL on exception, result must be freed by caller
# function geomType(ptr::GEOSGeom)
#     result = GEOSGeomType(ptr)
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSGeomType")
#     end
#     result
# end

# Return -1 on exception
function geomTypeId(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomTypeId_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGeomTypeId")
    end
    result
end

# Return 0 on exception
function getSRID(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetSRID_r(context.ptr, ptr)
    if result == 0
        error("LibGEOS: Error in GEOSGeomTypeId")
    end
    result
end

setSRID(ptr::GEOSGeom, context::GEOSContext = _context) = GEOSSetSRID_r(context.ptr, ptr)

# May be called on all geometries in GEOS 3.x, returns -1 on error and 1
# for non-multi geometries. Older GEOS versions only accept
# GeometryCollections or Multi* geometries here, and are likely to crash
# when fed simple geometries, so beware if you need compatibility with
# old GEOS versions.
"""
    numGeometries(geom, context=_context)

Returns the number of sub-geometries immediately under a
multi-geometry or collection or 1 for a simple geometry.
For nested collections, remember to check if returned
sub-geometries are **themselves** also collections.
"""
function numGeometries(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetNumGeometries_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGeomTypeId")
    end
    result
end

# Call only on GEOMETRYCOLLECTION or MULTI*
# (Return a pointer to the internal Geometry. Return NULL on exception.)
# Returned object is a pointer to internal storage:
# it must NOT be destroyed directly.
# Up to GEOS 3.2.0 the input geometry must be a Collection, in
# later version it doesn't matter (i.e. getGeometryN(0) for a single will return the input).
"""
    getGeometry(geom, n)

Returns a copy of the specified sub-geometry of a collection.
Numbering in one-based
For a simple geometry, returns a copy of the input.
"""
function getGeometry(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
    n in 1:numGeometries(ptr, context) || error(
        "GEOSGetGeometryN: $(numGeometries(ptr, context)) sub-geometries in geom, therefore n should be in 1:$(numGeometries(ptr, context))",
    )
    result = GEOSGetGeometryN_r(context.ptr, ptr, n - 1)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetGeometryN")
    end
    cloneGeom(result, context)
end

"""
    getGeometries(geom)

Returns a vector of copy of the sub-geometries of a collection.
For a simple geometry, returns the geometry in a vector of length one.
"""
getGeometries(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom[getGeometry(ptr, i, context) for i = 1:numGeometries(ptr, context)]

# Converts Geometry to normal form (or canonical form).
# Return -1 on exception, 0 otherwise.
function normalize!(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSNormalize_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSNormalize")
    end
    result
end

# Return -1 on exception
function numInteriorRings(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetNumInteriorRings_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGetNumInteriorRings")
    end
    result
end

# Call only on LINESTRING (returns -1 on exception)
function numPoints(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomGetNumPoints_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGeomGetNumPoints")
    end
    result
end

# Return -1 on exception, Geometry must be a Point.
function getGeomX(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeomGetX_r(context.ptr, ptr, out)
    if result == -1
        error("LibGEOS: Error in GEOSGeomGetX")
    end
    out[]
end

function getGeomY(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeomGetY_r(context.ptr, ptr, out)
    if result == -1
        error("LibGEOS: Error in GEOSGeomGetY")
    end
    out[]
end

# Return NULL on exception, Geometry must be a Polygon.
# Returned object is a pointer to internal storage: it must NOT be destroyed directly.
function interiorRing(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
    result = GEOSGetInteriorRingN_r(context.ptr, ptr, n - 1)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetInteriorRingN")
    end
    cloneGeom(result, context)
end

function interiorRings(ptr::GEOSGeom, context::GEOSContext = _context)
    n = numInteriorRings(ptr, context)
    if n == 0
        return GEOSGeom[]
    else
        return GEOSGeom[interiorRing(ptr, i, context) for i = 1:n]
    end
end

# Return NULL on exception, Geometry must be a Polygon.
# Returned object is a pointer to internal storage: it must NOT be destroyed directly.
function exteriorRing(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetExteriorRing_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetExteriorRing")
    end
    cloneGeom(result, context)
end

# Return -1 on exception
function numCoordinates(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetNumCoordinates_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGetNumCoordinates")
    end
    result
end

# Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
function getCoordSeq(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeom_getCoordSeq_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_getCoordSeq")
    end
    result
end
# getGeomCoordinates(ptr::GEOSGeom) = getCoordinates(getCoordSeq(ptr))

# Return 0 on exception (or empty geometry)
getGeomDimensions(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom_getDimensions_r(context.ptr, ptr)

# Return 2 or 3.
getCoordinateDimension(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom_getCoordinateDimension_r(context.ptr, ptr)

"""
    getXMin(geom, context=_context)

Finds the minimum X value in the geometry
"""
function getXMin(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeom_getXMin_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSGeom_getXMin_r")
    end
    return out[]
end

"""
    getYMin(geom, context=_context)

Finds the minimum Y value in the geometry
"""
function getYMin(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeom_getYMin_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSGeom_getYMin_r")
    end
    return out[]
end

"""
    getXMax(geom, context=_context)

Finds the maximum X value in the geometry
"""
function getXMax(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeom_getXMax_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSGeom_getXMax_r")
    end
    return out[]
end

"""
    getYMax(geom, context=_context)

Finds the maximum Y value in the geometry
"""
function getYMax(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeom_getYMax_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSGeom_getYMax_r")
    end
    return out[]
end

# TODO 02/2022: wait for libgeos release beyond 3.10.2 which will in include GEOSGeom_getExtent_r
# # Finds the extent (minimum and maximum X and Y value) of the geometry
# function getExtent(ptr::GEOSGeom, context::GEOSContext = _context)
#     out = Vector{Float64}(undef, 4)
#     result = GEOSGeom_getExtent_r(context.ptr, ptr, Ref(out, 1), Ref(out, 2), Ref(out, 3), Ref(out, 4))
#     if result == 0
#         error("LibGEOS: Error in GEOSGeom_getExtent_r")
#     end
#     return out
# end

# Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
function getPoint(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
    result = GEOSGeomGetPointN_r(context.ptr, ptr, n - 1)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeomGetPointN")
    end
    result
end

# Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
function startPoint(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomGetStartPoint_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeomGetStartPoint")
    end
    result
end

# Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
function endPoint(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomGetEndPoint_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeomGetEndPoint")
    end
    result
end

# -----
# Misc functions
# -----
function geomArea(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSArea_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSArea")
    end
    out[]
end

function geomLength(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSLength_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSLength")
    end
    out[]
end

function geomDistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSDistance_r(context.ptr, g1, g2, out)
    if result == 0
        error("LibGEOS: Error in GEOSDistance")
    end
    out[]
end

function hausdorffdistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSHausdorffDistance_r(context.ptr, g1, g2, out)
    if result == 0
        error("LibGEOS: Error in GEOSHausdorffDistance")
    end
    out[]
end

function hausdorffdistance(
    g1::GEOSGeom,
    g2::GEOSGeom,
    densifyFrac::Real,
    context::GEOSContext = _context,
)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSHausdorffDistanceDensify_r(context.ptr, g1, g2, densifyFrac, out)
    if result == 0
        error("LibGEOS: Error in GEOSHausdorffDistanceDensify")
    end
    out[]
end

# Return 0 on exception, the closest points of the two geometries otherwise.
# The first point comes from g1 geometry and the second point comes from g2.
nearestPoints(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context) =
    GEOSNearestPoints_r(context.ptr, g1, g2)

# -----
# Precision functions
# -----

"""
    getPrecision(geom)

Return the size of the geometry's precision grid, 0 for FLOATING precision.
"""
function getPrecision(geom::GEOSGeom, context::GEOSContext = _context)
    GEOSGeom_getPrecision_r(context.ptr, geom)
end

"""
    setPrecision(geom, gridSize; flags = GEOS_PREC_VALID_OUTPUT)

Change the rounding precision on a geometry. This will affect the precision of the existing geometry as well as any geometries derived from this geometry using overlay functions. The output will be a valid GEOSGeometry.

Note that operations will always be performed in the precision of the geometry with higher precision (smaller "gridSize"). That same precision will be attached to the operation outputs.

In the Default and `GEOS_PREC_KEEP_COLLAPSED` modes invalid input may cause an error to occur, unless the invalidity is below the scale of the requested precision

There are only 3 modes. The `GEOS_PREC_NO_TOPO` mode takes precedence over `GEOS_PREC_KEEP_COLLAPSED`. So the combination `GEOS_PREC_NO_TOPO || GEOS_PREC_KEEP_COLLAPSED` has the same semantics as `GEOS_PREC_NO_TOPO`

# Parameters
* `g`	Input geometry
* `gridSize`	cell size of grid to round coordinates to, or 0 for FLOATING precision
* `flags`	The bitwise OR of members of the `GEOSPrecisionRules` enum

# Returns
The precision reduced result. Caller must free with `GEOSGeom_destroy()` NULL on exception.

"""
function setPrecision(
    geom::GEOSGeom,
    gridSize::Real,
    flags::GEOSPrecisionRules,
    context::GEOSContext = _context,
)
    result = geomFromGEOS(GEOSGeom_setPrecision_r(context.ptr, geom, gridSize, flags))
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_setPrecision_r")
    end
    result
end
@eval function setPrecision(
    geom::GEOSGeom,
    gridSize::Real,
    flags::Int,
    context::GEOSContext = _context,
)
    !(flags in $(Int.(instances(GEOSPrecisionRules)))) && error(
        "flags value should be in $(Int.(instances(GEOSPrecisionRules))) or use GEOSPrecisionRules enum members",
    )
    return setPrecision(geom, gridSize, GEOSPrecisionRules(flags), context)
end


text:
include("web_render.jl")
include("repl_render.jl")

ext: jl


text:
using LaTeXStrings
using Printf

# * User flags

ABBREVIATE_POWER_OF_NEGATIVE_UNITY=false
ADVANCE_DELIMITERS=true

"""
    abbreviate_power_of_negative_unity(flag)

If `flag` is `true`, powers of `(-1)` will be printed as `(-)ᵖ`, if
`false` as `(-1)ᵖ`.
"""
function abbreviate_power_of_negative_unity(flag::Bool)
    global ABBREVIATE_POWER_OF_NEGATIVE_UNITY
    ABBREVIATE_POWER_OF_NEGATIVE_UNITY=flag
end

"""
    advance_delimiters(flag)

If `flag` is `true`, nested delimiters will alternate as `{x + [y + (z
+ w)]}` to highlight the level of nesting; if `false`, the delimiters
will not alternate: `(x + [y + (z + w)])`.
"""
function advance_delimiters(flag::Bool)
    global ADVANCE_DELIMITERS
    ADVANCE_DELIMITERS=flag
end

# * LaTeX conversion

# ** Numbers

function latex(sym::Symbolic)
    sym_str = replace(replace("$(sym)", "{"=>"\\{"), "}"=>"\\}")
    "\\textrm{$(sym_str)}",0
end
latex(s::Real) = "$(s)",0

function base_exp(v::T) where {T<:AbstractFloat}
    v == zero(T) && return (zero(T),0)
    r = log10(abs(v))
    e = floor(Int,r)
    b = sign(v)*T(10)^(r-e)
    b,e
end

function latex(s::AbstractFloat)
    b,e = base_exp(s)
    if abs(e) < 3
        @sprintf("%0.10g", s)
    else
        "$b\\times10^{$e}"
    end,0
end

function latex(z::Complex{T}) where T
    a,b = real(z),imag(z)
    la,lb = latex(a),latex(b)
    sla,slb = if T == Bool
        "",""
    else
        la[1],lb[1]
    end
    if a == zero(T)
        slb * "\\mathrm{i}"
    elseif b == zero(T)
        sla
    else
        sla * (b < 0 ? "" : "+") * slb * "\\mathrm{i}"
    end,max(la[2],lb[2])
end
latex(s::Sym) = "$(s)",0

# ** Operators

LaTeX_operators = Dict(:(+) => "+",
                       :(-) => "-",
                       :(*) => "",
                       :(/) => "/",
                       :(%) => "\\mod", # rem(::Symbolic) not implemented
                       :log => "\\ln",
                       :log10 => "\\log_{10}",
                       :log2 => "\\log_2",
                       :asin => "\\arcsin",
                       :asinh => "\\operatorname{arcsinh}",
                       :acos => "\\arccos",
                       :acosh => "\\operatorname{arccosh}",
                       :atan => "\\arctan",
                       :atanh => "\\operatorname{arctanh}")

# These functions have the same names in LaTeX
LaTeX_functions = [:sin, :sinh,
                   :cos, :cosh,
                   :tan, :tanh,
                   :exp,
                   :min, :max,
                   :mod,
                   :sqrt]

function latex_op(op::Sym)
    if op.name in keys(LaTeX_operators)
        LaTeX_operators[op.name]
    elseif op.name in LaTeX_functions
        "\\$(op.name)"
    else
        "$(op)"
    end
end

latex_op(op) = latex(op)

# ** Delimiters

get_next_delim(max_delim) =
    [("(",")"), ("[","]"), ("\\{", "\\}")][mod(max_delim,3)+1]

function delimit(val, max_delim, doit::Bool=true)
    !doit && return (val,max_delim)
    a,b = get_next_delim(max_delim)
    "\\left$(a)$(val)\\right$(b)", max_delim + (ADVANCE_DELIMITERS ? 1 : 0)
end

should_be_delimited(::Sym) = false
should_be_delimited(::Number) = true
should_be_delimited(r::Real) = isnegated(r)
should_be_delimited(v::Vector) = length(v) > 1 || should_be_delimited(v[1])

brace(s::String) = "{$(s)}"

# ** Arguments

function latex_arg(arg,should_delimit::Bool)
    larg = latex(arg)
    if should_delimit && (arg isa SymExpr && arg.op.name ∈ [:(+)] ||
                         arg isa Complex ||
                         arg isa Real && arg < 0)
        delimit(larg...)
    else
        larg
    end
end

latex_args(args::Vector,should_delimit::Bool) =
    map(arg -> latex_arg(arg, should_delimit), args)

latex_args(expr::SymExpr,should_delimit::Bool) =
    latex_args(expr.args, should_delimit)

latex_args(s::Number,::Bool) = [latex_arg(s, false)]

# ** Expressions

function latex(expr::SymExpr)
    op = expr.op
    lop = latex_op(op)
    # The arguments need only be wrapped in delimiters for products
    # and powers.
    largs = latex_args(expr.args, op.name ∈ [:(*), :(^)])
    max_delim = maximum(last.(largs))

    if op.name == :(+)
        S = first(largs[1])
        for a in first.(largs[2:end])
            S *= first(lstrip(a, ['{'])) == '-' ? "" : "+"
            S *= a
        end
        S,max_delim
    elseif op.name == :(^)
        arg1 = expr.args[1]
        earg2 = latex_arg(expr.args[2], false)
        earg2s = "^{$(earg2[1])}"
        # If the user has requested thus, (-1)^z will be printed as (-)^z
        if ABBREVIATE_POWER_OF_NEGATIVE_UNITY &&
            arg1 isa Number && !(arg1 isa SymExpr) && arg1 isa Real && arg1 == -1
            "(-)$(earg2s)",max((ADVANCE_DELIMITERS ? 1 : 0),earg2[2])
        elseif arg1 isa SymExpr && arg1.op.name ∉ [:(+), :(*), :(^), :sqrt]
            larg1 = latex_args(arg1.args,false)
            max_delim = max(maximum(last.(larg1)), earg2[2])
            slarg1 = join(first.(larg1), ",")
            if should_be_delimited(arg1.args)
                slarg1,max_delim = delimit(slarg1,max_delim)
            end
            "$(latex_op(arg1.op))$(earg2s)$(slarg1)",max_delim
        else
            "$(largs[1][1])$(earg2s)",max(largs[1][2],earg2[2])
        end
    elseif op.name == :(*)
        i = findfirst(isequal(-1), expr.args)
        isneg = i !== nothing
        expr = !isneg ? expr : prod(expr.args[vcat(1:i-1,i+1:length(expr.args))])

        num = numerator(expr)
        den = denominator(expr)

        max_delim = 0
        # If there is a factor -1, represent that as a prefactor -.
        (isneg ? "-" : "")*if den == 1
            # Only numerator
            if num isa SymExpr && num.op.name == :(*)
                nlargs = latex_args(num, true)
                max_delim = maximum(last.(nlargs))
                join(first.(nlargs),"")
            else
                # Dropping a possible -1 factor may have yielded an
                # expression that is not a multiplication anymore.
                nlarg,max_delim = delimit(latex(num)...,
                                          isneg && should_be_delimited(num))
                nlarg
            end
        else
            # Fraction
            lnum = latex(num)
            lden = latex(den)
            max_delim = max(lnum[2],lden[2])
            "\\frac{$(lnum[1])}{$(lden[1])}"
        end,max_delim
    else
        # Any other function
        slargs = join(first.(largs), ",")

        dargs,max_delim = if op.name != :sqrt && (length(largs) > 1 ||
                                                 should_be_delimited(first(expr.args)))
            delimit(slargs, max_delim)
        else
            brace(slargs),max_delim
        end

        "$(lop)$(dargs)",max_delim
    end
end

# * Convert to LaTeXString, display

Base.convert(::Type{LaTeXString}, s::S) where {S<:Symbolic} =
    latexstring("\$$(latex(s)[1])\$")

Base.show(io::IO, ::MIME"text/latex", s::S) where {S<:Symbolic} =
    show(io, "text/latex", convert(LaTeXString, s))


text:
module Transformation

export Vec, Quat, Transf, dot, cross, magnitude, norm, ×, slerp, unit

"""
    Vec - three dimensional vector.

# Examples
```julia-repl
julia> Vec(1, 2, 3)
Vec(1, 2, 3)
```
"""
struct Vec
    x::Real
    y::Real
    z::Real
    function Vec(x::Real=0, y::Real=0, z::Real=0)::Vec
        new(x, y, z)
    end
end

function Base.:+(a::Vec, b::Vec)::Vec
    Vec(
        a.x + b.x,
        a.y + b.y,
        a.z + b.z
    )
end

function Base.:-(a::Vec, b::Vec)::Vec
    Vec(
        a.x - b.x,
        a.y - b.y,
        a.z - b.z
    )
end

function Base.:-(a::Vec)::Vec
    Vec(
        -a.x,
        -a.y,
        -a.z
    )
end

"""
    dot(a::Vec, b::Vec)::Vec

Compute the dot product of two vectors.
"""
function dot(a::Vec, b::Vec)
    a.x * b.x + a.y * b.y + a.z * b.z
end

"""
    cross(a::Vec, b::Vec)::Vec

Compute the cross product of two vectors.
"""
function cross(a::Vec, b::Vec)::Vec
    Vec(
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x
    )
end

"""
    ×(a::Vec, b::Vec)::Vec

Compute the cross product of two vectors.
"""
function ×(a::Vec, b::Vec)::Vec
    cross(a, b)
end

function Base.:*(a::Real, b::Vec)::Vec
    Vec(
        a * b.x,
        a * b.y,
        a * b.z
    )
end

function Base.:*(a::Vec, b::Real)::Vec
    Vec(
        a.x * b,
        a.y * b,
        a.z * b
    )
end

function Base.:/(a::Vec, b::Real)::Vec
    Vec(
        a.x / b,
        a.y / b,
        a.z / b
    )
end

function Base.isequal(a::Vec, b::Vec)::Bool
    a.x == b.x && a.y == b.y && a.z == b.z
end

"""
    norm(a::Vec)::Real

Compute the Euclidean norm of the given vector.
"""
function norm(a::Vec)::Real
    √(a.x^2 + a.y^2 + a.z^2)
end

"""
    unit(a::Vec)::Vec

Compute the unit vector colinear to the given one.
"""
function unit(a::Vec)::Vec
    a / norm(a)
end

"""
    round(a::Vec)::Vec

Compute the vector with it's elements rounder to nearest value.
"""
function Base.round(a::Vec)::Vec
    Vec(
        round(a.x),
        round(a.y),
        round(a.z)
    )
end

"""
    Quat - hypercomplex number.

Unit quaternions are used for rotation representation.

# Examples
```julia-repl
julia> # 120° counterclockwise rotation around the (1, 1, 1) axis
julia> Quat(2π/3, Vec(1, 1, 1))
Quat(0.5, 0.5, 0.5, 0.5)

julia> # 45° clockwise rotation around the (1, 0, 0) axis
julia> Quat(-π/4, Vec(1, 0, 0))
Quat(0.9238795325112867, -0.3826834323650898, -0.0, -0.0)
```
"""
struct Quat
    w::Real
    x::Real
    y::Real
    z::Real
    function Quat(w::Real=1, x::Real=0, y::Real=0, z::Real=0)::Quat
        new(w, x, y, z)
    end

    function Quat(angle::Real, axis::Vec)::Quat
        w = cos(angle / 2)
        xyz = unit(axis) * sin(angle / 2)
        Quat(w, xyz.x, xyz.y, xyz.z)
    end
end

function Base.:+(a::Quat, b::Quat)::Quat
    Quat(
        a.w + b.w,
        a.x + b.x,
        a.y + b.y,
        a.z + b.z
    )
end

function Base.:-(a::Quat, b::Quat)::Quat
    Quat(
        a.w - b.w,
        a.x - b.x,
        a.y - b.y,
        a.z - b.z
    )
end

function Base.:*(a::Quat, b::Quat)::Quat
    Quat(
        a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
        a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
        a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
        a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w
    )
end

function Base.:*(q::Quat, v::Vec)::Vec
    r = q * Quat(0, v.x, v.y, v.z) * conj(q)
    Vec(r.x, r.y, r.z)
end

function Base.:*(a::Quat, b::Real)::Quat
    Quat(
        a.w * b,
        a.x * b,
        a.y * b,
        a.z * b
    )
end

function Base.:*(a::Real, b::Quat)::Quat
    Quat(
        a * b.w,
        a * b.x,
        a * b.y,
        a * b.z
    )
end

function dot(a::Quat, b::Quat)::Real
    a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z
end

function Base.:/(a::Quat, b::Real)::Quat
    Quat(
        a.w / b,
        a.x / b,
        a.y / b,
        a.z / b
    )
end

"""
    conj(a::Quat)::Quat

Compute the quaternion conjugate.
If the quaternion is unit, it's conjugate represents the reverse rotation.

# Examples
```julia-repl
julia> q = Quat(π, Vec(1, 1, 1))
julia> q * conj(q)
Quat(1.0, 0.0, 0.0, 0.0)
```
"""
function Base.conj(a::Quat)::Quat
    Quat(
        a.w,
        -a.x,
        -a.y,
        -a.z
    )
end

"""
    norm(a::Quat)::Real

Compute the quaternion norm.
"""
function norm(a::Quat)::Real
    √(a.w^2 + a.x^2 + a.y^2 + a.z^2)
end

"""
    unit(a::Quat)::Quat

Compute the unit quaternion colinear to the given one.
"""
function unit(a::Quat)::Quat
    a / norm(a)
end

function slerp(a::Quat, b::Quat, t)::Quat
    d = dot(a, b)
    if d < 0
        b = -b
        d = -d
    end
    if d > 0.9995
        result = a + (a - b) * t
        return norm(result)
    end
    theta0 = acos(d)
    theta = theta0 * t
    sinTheta = sin(theta)
    sinTheta0 = sin(theta0)

    s1 = sinTheta / sinTheta0
    s0 = cos(theta) - d * s1
    return a * s0 + b * s1
end

"""
    Transf - spatial transformation representation.

The transformation consists of consecutive translation on vector `v` and rotation on `q` quaternion.
"""
struct Transf
    v::Vec
    q::Quat
    function Transf(v::Vec=Vec(), q::Quat=Quat())::Transf
        new(v, q)
    end
end

function Base.:+(a::Transf, b::Transf)::Transf
    Transf(
        a.v + a.q * b.v,
        a.q * b.q
    )
end

function Base.:+(a::Transf, b::Vec)::Vec
    a.v + a.q * b
end

"""
    conj(a::Transf)::Transf

Compute such transformation, that `a + conj(a)` equals to zero transformation.

# Examples
```julia-repl
julia> t = Transf(Vec(1, 2, 3), Quat(π/3, Vec(1, 1, 0)))
julia> t + conj(t)
Transf(Vec(0.0, 0.0, 0.0), Quat(1.0, 0.0, 0.0, 0.0))
```
"""
function Base.conj(a::Transf)::Transf
    q = conj(a.q)
    Transf(
        q * (-a.v),
        q
    )
end

end

max_stars_repo_name: red-hara/Transformation

max_issues_repo_name: red-hara/Transformation

max_forks_repo_name: red-hara/Transformation


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VariableTemplates: flattenednames
using CLIMA.HydrostaticBoussinesq
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK
using CLIMA.PlanetParameters: grav
using CLIMA.HydrostaticBoussinesq: AbstractHydrostaticBoussinesqProblem
import CLIMA.HydrostaticBoussinesq: ocean_init_aux!, ocean_init_state!,
                                    ocean_boundary_state!,
                                    CoastlineFreeSlip, CoastlineNoSlip,
                                    OceanFloorFreeSlip, OceanFloorNoSlip,
                                    OceanSurfaceNoStressNoForcing,
                                    OceanSurfaceStressNoForcing,
                                    OceanSurfaceNoStressForcing,
                                    OceanSurfaceStressForcing
import CLIMA.DGmethods: update_aux!, update_aux_diffusive!,
                        vars_state, vars_aux
using GPUifyLoops

const ArrayType = CLIMA.array_type()

HBModel   = HydrostaticBoussinesqModel

struct OceanGyre{T} <: AbstractHydrostaticBoussinesqProblem
  Lˣ::T
  Lʸ::T
  H::T
  τₒ::T
  fₒ::T
  β::T
  λʳ::T
  θᴱ::T
end

@inline function ocean_boundary_state!(m::HBModel, p::OceanGyre, bctype, x...)
  if bctype == 1
    ocean_boundary_state!(m, CoastlineNoSlip(), x...)
  elseif bctype == 2
    ocean_boundary_state!(m, OceanFloorNoSlip(), x...)
  elseif bctype == 3
    ocean_boundary_state!(m, OceanSurfaceStressForcing(), x...)
  end
end

# A is Filled afer the state
function ocean_init_aux!(m::HBModel, P::OceanGyre, A, geom)
  FT = eltype(A)
  @inbounds y = geom.coord[2]

  Lʸ = P.Lʸ
  τₒ = P.τₒ
  fₒ = P.fₒ
  β  = P.β
  θᴱ = P.θᴱ

  A.τ  = -τₒ * cos(y * π / Lʸ)
  A.f  =  fₒ + β * y
  A.θʳ =  θᴱ * (1 - y / Lʸ)

  A.ν = @SVector [m.νʰ, m.νʰ, m.νᶻ]
  A.κ = @SVector [m.κʰ, m.κʰ, m.κᶻ]
end

function ocean_init_state!(P::OceanGyre, Q, A, coords, t)
  @inbounds z = coords[3]
  @inbounds H = P.H

  Q.u = @SVector [0,0]
  Q.η = 0
  Q.θ = 9 + 8z/H
end

###################
# PARAM SELECTION #
###################
FT = Float64
vtkpath = "vtk_ekman_spiral"

const timeend = 3 * 30 * 86400   # s
const tout    = 6 * 24 * 60 * 60 # s

const N  = 4
const Nˣ = 20
const Nʸ = 20
const Nᶻ = 20
const Lˣ = 4e6  # m
const Lʸ = 4e6  # m
const H  = 1000 # m

xrange = range(FT(0);  length=Nˣ+1, stop=Lˣ)
yrange = range(FT(0);  length=Nʸ+1, stop=Lʸ)
zrange = range(FT(-H); length=Nᶻ+1, stop=0)

const cʰ = sqrt(grav * H)
const cᶻ = 0

const τₒ = 1e-1  # (m/s)^2
const fₒ = 1e-4  # Hz
const β  = 1e-11 # Hz / m
const θᴱ = 25    # K

const αᵀ = 2e-4  # (m/s)^2 / K
const νʰ = 5e3   # m^2 / s
const νᶻ = 5e-3  # m^2 / s
const κʰ = 1e3   # m^2 / s
const κᶻ = 1e-4  # m^2 / s
const λʳ = 4 // 86400 # m / s
let
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD

  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
             ll == "WARN"  ? Logging.Warn  :
             ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  brickrange = (xrange, yrange, zrange)
  topl = StackedBrickTopology(mpicomm, brickrange;
                              periodicity = (false, false, false),
                              boundary = ((1, 1), (1, 1), (2, 3)))

  minΔx = Lˣ / Nˣ / (N + 1)
  minΔz = H  / Nᶻ / (N + 1)
  CFL_gravity = minΔx / cʰ
  CFL_diffusive = minΔz^2 / (1000 * κᶻ)
  CFL_viscous = minΔz^2 / νᶻ
  dt = 1//2 * minimum([CFL_gravity, CFL_diffusive, CFL_viscous])
  nout = ceil(Int64, tout / dt)
  dt = tout / nout

  @info @sprintf("""Update
                    Gravity CFL   = %.1f
                    Diffusive CFL = %.1f
                    Viscous CFL   = %.1f
                    Timestep      = %.1f""",
                 CFL_gravity, CFL_diffusive, CFL_viscous, dt)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  prob = OceanGyre{FT}(Lˣ, Lʸ, H, τₒ, fₒ, β, λʳ, θᴱ)

  model = HBModel{FT}(prob, cʰ = cʰ)

  dg = OceanDGModel(model,
                    grid,
                    Rusanov(),
                    CentralNumericalFluxDiffusive(),
                    CentralNumericalFluxGradient())

  Q = init_ode_state(dg, FT(0); forcecpu=true)
  update_aux!(dg, model, Q, FT(0))
  update_aux_diffusive!(dg, model, Q, FT(0))

  if isdir(vtkpath)
    rm(vtkpath, recursive=true)
  end
  mkpath(vtkpath)
  mkpath(vtkpath*"/weekly")
  mkpath(vtkpath*"/monthly")

  step = [0, 0]
  function do_output(span, step)
    outprefix = @sprintf("%s/%s/mpirank%04d_step%04d",vtkpath, span,
                         MPI.Comm_rank(mpicomm), step)
    @info "doing VTK output" outprefix
    statenames = flattenednames(vars_state(model, eltype(Q)))
    auxnames = flattenednames(vars_aux(model, eltype(Q)))
    writevtk(outprefix, Q, dg, statenames, dg.auxstate, auxnames)
  end

  do_output("weekly", step[1])
  cbvtkw = GenericCallbacks.EveryXSimulationSteps(nout)  do (init=false)
    do_output("weekly", step[1])
    step[1] += 1
    nothing
  end

  do_output("monthly", step[2])
  cbvtkm = GenericCallbacks.EveryXSimulationSteps(5*nout)  do (init=false)
    do_output("monthly", step[2])
    step[2] += 1
    nothing
  end

  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  lsrk = LSRK144NiegemannDiehlBusch(dg, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e
  ArrayType = %s""" eng0 ArrayType

  solve!(Q, lsrk, nothing; timeend=timeend, callbacks=(cbinfo,cbvtkw,cbvtkm))

  return nothing
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
"""
    RotatingTimeValue{R, T}

A wrapper around a time value for distinguishing potentially different
rotation group it could belong to in a rotating sampling design.
See also [`rotatingtime`](@ref) and [`settime`](@ref).

# Fields
- `rotation::R`: a rotation group in a rotating sampling design.
- `time::T`: a time value belonged to the rotation group.
"""
struct RotatingTimeValue{R, T}
    rotation::R
    time::T
end

RotatingTimeValue(::Type{RotatingTimeValue{R,T}}, rotation, time) where {R,T} =
    RotatingTimeValue(convert(R, rotation), convert(T, time))

"""
    rotatingtime(rotation, time)

Construct [`RotatingTimeValue`](@ref)s from `rotation` and `time`.
This method simply broadcasts the default constructor over the arguments.
"""
rotatingtime(rotation, time) = RotatingTimeValue.(rotation, time)

# Compare time first
isless(x::RotatingTimeValue, y::RotatingTimeValue) =
    isequal(x.time, y.time) ? isless(x.rotation, y.rotation) : isless(x.time, y.time)

isless(x::RotatingTimeValue, y) = isless(x.time, y)
isless(x, y::RotatingTimeValue) = isless(x, y.time)
isless(::RotatingTimeValue, ::Missing) = true
isless(::Missing, ::RotatingTimeValue) = false

Base.isequal(x::RotatingTimeValue, y::RotatingTimeValue) =
    isequal(x.rotation, y.rotation) && isequal(x.time, y.time)

function ==(x::RotatingTimeValue, y::RotatingTimeValue)
    req = x.rotation == y.rotation
    isequal(req, missing) && return missing
    teq = x.time == y.time
    isequal(teq, missing) && return missing
    return req && teq
end

==(x::RotatingTimeValue, y) = x.time == y
==(x, y::RotatingTimeValue) = x == y.time
==(::RotatingTimeValue, ::Missing) = missing
==(::Missing, ::RotatingTimeValue) = missing

Base.zero(::Type{RotatingTimeValue{R,T}}) where {R,T} = RotatingTimeValue(zero(R), zero(T))
Base.iszero(x::RotatingTimeValue) = iszero(x.time)
Base.one(::Type{RotatingTimeValue{R,T}}) where {R,T} = RotatingTimeValue(one(R), one(T))
Base.isone(x::RotatingTimeValue) = isone(x.time)

Base.convert(::Type{RotatingTimeValue{R,T}}, x::RotatingTimeValue) where {R,T} =
    RotatingTimeValue(convert(R, x.rotation), convert(T, x.time))

Base.nonmissingtype(::Type{RotatingTimeValue{R,T}}) where {R,T} =
    RotatingTimeValue{nonmissingtype(R), nonmissingtype(T)}

Base.iterate(x::RotatingTimeValue) = (x, nothing)
Base.iterate(x::RotatingTimeValue, ::Any) = nothing

Base.length(x::RotatingTimeValue) = 1

show(io::IO, x::RotatingTimeValue) = print(io, x.rotation, '_', x.time)
function show(io::IO, ::MIME"text/plain", x::RotatingTimeValue)
    println(io, typeof(x), ':')
    println(io, "  rotation: ", x.rotation)
    print(io, "  time:     ", x.time)
end

"""
    RotatingTimeArray{T<:RotatingTimeValue,N,C,I} <: AbstractArray{T,N}

Array type for [`RotatingTimeValue`](@ref)s that stores
the field values `rotation` and `time` in two arrays for efficiency.
The two arrays that hold the field values for all elements can be accessed as properties.
"""
struct RotatingTimeArray{T<:RotatingTimeValue,N,C,I} <: AbstractArray{T,N}
    a::StructArray{T,N,C,I}
    RotatingTimeArray(a::StructArray{T,N,C,I}) where {T,N,C,I} = new{T,N,C,I}(a)
end

"""
    RotatingTimeArray(rotation::AbstractArray, time::AbstractArray)

Construct a [`RotatingTimeValue`](@ref) from arrays of `rotation` and `time`.
"""
function RotatingTimeArray(rotation::AbstractArray, time::AbstractArray)
    a = StructArray{RotatingTimeValue{eltype(rotation), eltype(time)}}((rotation, time))
    return RotatingTimeArray(a)
end

_getarray(a::RotatingTimeArray) = getfield(a, :a)

Base.size(a::RotatingTimeArray) = size(_getarray(a))
Base.IndexStyle(::Type{<:RotatingTimeArray{<:Any,<:Any,<:Any,I}}) where I =
    I === Int ? IndexLinear() : IndexCartesian()

Base.@propagate_inbounds Base.getindex(a::RotatingTimeArray, i::Int) = _getarray(a)[i]
Base.@propagate_inbounds Base.getindex(a::RotatingTimeArray, I) =
    RotatingTimeArray(_getarray(a).rotation[I], _getarray(a).time[I])

Base.@propagate_inbounds Base.setindex!(a::RotatingTimeArray, v, i::Int) =
    setindex!(_getarray(a), v, i)
Base.@propagate_inbounds Base.setindex!(a::RotatingTimeArray, v, I) =
    setindex!(_getarray(a), v, I)

@inline Base.view(a::RotatingTimeArray, I...) = RotatingTimeArray(view(_getarray(a), I...))

Base.similar(a::RotatingTimeArray, dims::Dims=size(a)) =
    RotatingTimeArray(similar(_getarray(a), dims))

Base.similar(a::RotatingTimeArray, dims::Int...) = similar(a, dims)

Base.getproperty(a::RotatingTimeArray, n::Symbol) = getproperty(_getarray(a), n)

DataAPI.refarray(a::RotatingTimeArray) =
    RotatingTimeArray(DataAPI.refarray(a.rotation), DataAPI.refarray(a.time))

const ValidTimeType = Union{Signed, TimeType, Period, RotatingTimeValue}


text:
module WaterReservoir

using Interpolations

include("reservoir-geometry.jl")
include("inflow.jl")
include("outlet.jl")

Base.@kwdef struct ReservoirModel{G<:ReservoirGeometry,I<:Inflow,O<:Outlet}
    geometry::G
    inflow::I
    outlet::O
    z₀::Float64 = 0.0
end

Base.@kwdef struct Simulation{M <: ReservoirModel}
    model::M
    Δt::Float64
    n::Int
end

function run(sim)
    # Result storage
    t = Vector{Float64}(undef, sim.n + 1)
    z = similar(t)
    Qout = similar(t)
    # Initial conditions
    t[1] = 0.0
    z[1] = sim.model.z₀
    Qout[1] = discharge(sim.model.outlet, z = z[1])
    # Iterations
    for i = 1:sim.n
        Qin = discharge(sim.model.inflow, t[i])
        S = surface_area(sim.model.geometry, z[i])
        t[i+1] = i * sim.Δt
        z[i+1] = z[i] + sim.Δt * (Qin - Qout[i]) / S
        Qout[i+1] = discharge(sim.model.outlet, z = z[i+1])
    end
    return t, z, Qout
end

end # module


text:
using StatsModels
import Splines2.bs
import Base.show
const bsPLINE_CONTEXT = Any
spl(x,df) = Splines2.bs(x,df=df,intercept=true) # assumes intercept
mutable struct uf_bsplineTerm{T,D} <: AbstractTerm
    term::T
    df::D
    fun # function handle
end
function uf_bsplineTerm(term,df)
    uf_bsplineTerm(term,df,nothing)
end

Base.show(io::IO, p::uf_bsplineTerm) = print(io, "spl($(p.term), $(p.df))")

function StatsModels.apply_schema(t::FunctionTerm{typeof(spl)},
                                  sch::StatsModels.Schema,
                                  Mod::Type{<:bsPLINE_CONTEXT})
    apply_schema(uf_bsplineTerm(t.args_parsed...), sch, Mod)
end
function StatsModels.apply_schema(t::uf_bsplineTerm,
                                  sch::StatsModels.Schema,
                                  Mod::Type{<:bsPLINE_CONTEXT})
    term = apply_schema(t.term, sch, Mod)
    isa(term, ContinuousTerm) ||
        throw(ArgumentError("uf_bsplineTerm only works with continuous terms (got $term)"))

    isa(t.df, ConstantTerm) ||
            throw(ArgumentError("uf_bsplineTerm df must be a number (got $(t.df))"))
    uf_bsplineTerm(term, t.df.n)
end
function StatsModels.modelcols(p::uf_bsplineTerm, d::NamedTuple)
    col = modelcols(p.term, d)
    if isnothing(p.fun)
        p.fun = Splines2.bs_(col,df=p.df+1,intercept=true)
    end
        #X = Splines2.bs(col, df=p.df+1,intercept=true)
    X = p.fun(col)
    # remove middle X to negate intercept = true, generating a pseudo effect code
    X[:,Not(Int(ceil(end/2)))]
end
StatsModels.terms(p::uf_bsplineTerm) = terms(p.term)
StatsModels.termvars(p::uf_bsplineTerm) = StatsModels.termvars(p.term)
StatsModels.width(p::uf_bsplineTerm) = 1
StatsModels.coefnames(p::uf_bsplineTerm) = "spl(" .* coefnames(p.term) .* "," .* string.(1:p.df) .* ")"


text:
#--------------------------------------------------------------------------------------------------
# Additional functions tests
#--------------------------------------------------------------------------------------------------
using Test
using Hensel

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

@test Hensel.logistic(Hensel.logistic(Hensel.logistic(0.3, 4.0), 4.0), 4.0) ==  Hensel.logistic(0.3,4.0,3)
@test Hensel.dyadic(Hensel.dyadic(Hensel.dyadic(0.3))) ==  Hensel.dyadic(0.3,3)
@test Hensel.nyadic(Hensel.nyadic(Hensel.nyadic(0.3, 5), 5), 5) ==  Hensel.nyadic(0.3,5,3)
@test Hensel.cfshift(Hensel.cfshift(Hensel.cfshift(0.3))) ==  Hensel.cfshift(0.3,3)

@test ((y) -> Hensel.FunEnv((x) -> Hensel.dyadic(0.35,y), (0, 1), 2, 8)(y, return_type="integer")).(0:10) == [
    154,
    205,
    102,
     51,
    153,
    204,
    102,
     51,
    153,
    204,
    102
]

@test ((y) -> Hensel.FunEnv((x) -> Hensel.nyadic(0.35,2,y), (0, 1), 2, 8)(y, return_type="integer")).(0:10) == [
    154,
    205,
    102,
     51,
    153,
    204,
    102,
     51,
    153,
    204,
    102
]

@test Hensel.FunEnv(Hensel.dyadic, (0, 1), 2, 8)(1,2,return_type="integer") == 64
@test Hensel.FunEnv(Hensel.dyadic, (0, 1), 2, 8)((1,2,"integer")) == 64
@test Hensel.FunEnv(Hensel.dyadic, (0, 1), 2, 8)((1,2)) == 64
@test Hensel.FunEnv(Hensel.dyadic, (0, 1), 2, 8).([(1,1),(1,2),(1,3)]) == [128, 64, 32]
@test Hensel.FunEnv(Hensel.dyadic, (0, 1), 2, 8).([(1,1,"hensel"),(1,2,"hensel"),(1,3,"hensel")]) == [
    [0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 1, 0],
    [0, 0, 0, 0, 0, 1, 0, 0]
]

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

max_stars_repo_name: ayarodionov/Hensel

max_issues_repo_name: ayarodionov/Hensel

max_forks_repo_name: ayarodionov/Hensel


text:
type ExponentialRungeKuttaEquidistantTimeStepperIterator
    psi::WaveFunction
    t0::Real
    tend::Real
    dt::Real
    steps::Int
    scheme::Int
    s::Int
    G::WaveFunction
    U::Vector{WaveFunction}
end

function exponential_runge_kutta_equidistant_time_stepper(psi::WaveFunction,
         t0::Real, tend::Real, dt::Real; steps::Int=-1, scheme::Symbol=:krogstad)
    m = psi.m
    set_time!(psi, t0)
    G = wave_function(m)
    scheme1 = 42
    s = 4
    if scheme==:krogstad
        scheme1=42
        s = 4
#    elseif scheme==:cox_mathews
#        scheme1=41
#        s = 4
    elseif scheme==:strehmel_weiner
        scheme1=43
        s = 4
    else
        warn("scheme not known, I use :krogstad")
    end
    U = [wave_function(m) for j=1:s]
    ExponentialRungeKuttaEquidistantTimeStepperIterator(psi, t0, tend, dt, steps, scheme1, s, G, U)
end

Base.start(tsi::ExponentialRungeKuttaEquidistantTimeStepperIterator) = (tsi.steps<0 ? tsi.t0 : 0)

Base.done(tsi::ExponentialRungeKuttaEquidistantTimeStepperIterator, t) = (tsi.steps<0 ? t >= tsi.tend : t>=tsi.steps )

function ERK4_step_krogstad!(psi::WaveFunction, dt::Number;
    #Krogstad, eq. (2.42) in Hochbruck/Ostermann
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])

    s = 4
    c = [0.0, 0.5, 0.5, 1.0]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[2], c[2], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 2, -1.0*dt)
    add_phi_A!(G, U[4], c[4], 1, 1.0*dt)
    add_phi_A!(G, U[4], c[4], 2, -2.0*dt)
    add_phi_A!(G, psi, dt, 1, 1.0*dt)
    add_phi_A!(G, psi, dt, 2, -3.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[3], c[3], 2, 1.0*dt)
    add_phi_A!(G, psi, dt, 2, 2.0*dt)
    add_phi_A!(G, psi, dt, 3, -4.0*dt)

    gen_rhs!(G, U[3])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[4], c[4], 2, 2.0*dt)
    add_phi_A!(G, psi, dt, 2, 2.0*dt)
    add_phi_A!(G, psi, dt, 3, -4.0*dt)

    gen_rhs!(G, U[4])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, psi, dt, 2, -1.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)
    psi
end

function ERK4_step_strehmel_weiner!(psi::WaveFunction, dt::Number;
    #Strehmel/Weiner, eq. (2.43) in Hochbruck/Ostermann
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])

    s = 4
    c = [0.0, 0.5, 0.5, 1.0]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[2], c[2], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 2, -0.5*dt)
    add_phi_A!(G, U[4], c[4], 1, 1.0*dt)
    add_phi_A!(G, U[4], c[4], 2, -2.0*dt)
    add_phi_A!(G, psi, dt, 1, 1.0*dt)
    add_phi_A!(G, psi, dt, 2, -3.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[3], c[3], 2, 0.5*dt)
    add_phi_A!(G, U[4], c[4], 2, -2.0*dt)

    gen_rhs!(G, U[3])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[4], c[4], 2, 4.0*dt)
    add_phi_A!(G, psi, dt, 2, 4.0*dt)
    add_phi_A!(G, psi, dt, 3, -8.0*dt)

    gen_rhs!(G, U[4])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, psi, dt, 2, -1.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)
    psi
end

function Base.next(tsi::ExponentialRungeKuttaEquidistantTimeStepperIterator, t)
    if tsi.scheme==42
        ERK4_step_krogstad!(tsi.psi, tsi.dt, G=tsi.G, U=tsi.U)
    elseif tsi.scheme==43
        ERK4_step_strehmel_weiner!(tsi.psi, tsi.dt, G=tsi.G, U=tsi.U)
    elseif tsi.scheme==41
        ERK4_step_cox_mathews!(tsi.psi, tsi.dt, G=tsi.G, U=tsi.U)
    end

    if tsi.steps<0
        t1 = t + tsi.dt < tsi.tend ? t + tsi.dt : tsi.tend
    else
        t1 = t+1
    end
    return t1, t1
 end

 function global_orders(psi::WaveFunction, reference_solution::WaveFunction,
                       t0::Real, tend::Real, dt::Real; scheme::Symbol=:krogstad, rows=8)
    @assert psi.m==reference_solution.m
    tab = Array(Float64, rows, 3)

    wf_save_initial_value = clone(psi)
    copy!(wf_save_initial_value, psi)

    steps = Int(floor((tend-t0)/dt))
    dt1 = dt
    err_old = 0.0
    println("             dt         err      p")
    println("-----------------------------------")
    for row=1:rows
        for t in exponential_runge_kutta_equidistant_time_stepper(psi, t0, tend, dt1, steps=steps, scheme=scheme)
        end
        err = distance(psi, reference_solution)
        if (row==1) then
            @printf("%3i%12.3e%12.3e\n", row, Float64(dt1), Float64(err))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = 0
        else
            p = log(err_old/err)/log(2.0);
            @printf("%3i%12.3e%12.3e%7.2f\n", row, Float64(dt1), Float64(err), Float64(p))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = p
        end
        err_old = err
        dt1 = 0.5*dt1
        steps = 2*steps
        copy!(psi,wf_save_initial_value)
    end
end


text:
import Base: start, next, done, iteratorsize, iteratoreltype, eltype, length

"""
    KalmanFilter(y, M)

Kalman filter as iterator, iterating over `Gaussian`s representing
the filtered distribution of `x`. Arguments `y` iterates over signal values.

# Example

```
kf = KalmanFilter(Y, M) #
est = collect(kf)
```
"""
struct KalmanFilter{Tit,TM}
    it::Tit
    M::TM
end

Base.iteratorsize(::KalmanFilter{Tit}) where {Tit} = Base.iteratorsize(Tit) == Base.HasShape() ? Base.HasLength() : Base.iteratorsize(Tit)
Base.iteratoreltype(::KalmanFilter) = Base.HasEltype()
Base.eltype(::Type{KalmanFilter{Tit, T}}) where {Tit, T<:LinearHomogSystem{Tx, TP}} where {Tx, TP} = Gaussian{Tx,TP}

start(kf::KalmanFilter) = start(kf.it), prior(kf.M)
function next(kf::KalmanFilter, state)
    st, N = state
    y, st2 = next(kf.it, st)
    _, x, P = kalman_kernel(0.0, mean(N), cov(N), 1.0, y, kf.M)
    N2 = Gaussian(x, P)
    N2, (st2, N2)
end

done(kf::KalmanFilter, state) = done(kf.it, state[1])

length(kf::KalmanFilter) = length(kf.it)

"""
    MappedKalmanFilter(ty, M, f)

Kalman filter with model `M` as iterator, iterating over the result `ret` of
`ret = f(t, x, P, Ppred, ll, K)`. `ty` iterates over pairs `(t, y)`
of signal time and signal value.

Eltype is determined by calling `mappedreturntype(::typeof(f))`
with fallback any.

# Example
```
f(t, x, P, Ppred, ll, K) = ll
kf = MappedKalmanFilter(zip(Base.Iterators.countfrom(1), Y), M, f)
ll = sum(kf)
```

For performance, define:
```
import Kalman: mappedreturntype
Kalman.mappedreturntype(_, ::Type{typeof(f)}) = Float64
```
"""
struct MappedKalmanFilter{Tit,TM,F}
    it::Tit
    M::TM
    f::F
end

Base.iteratorsize(::MappedKalmanFilter{Tit}) where {Tit} = Base.iteratorsize(Tit) == Base.HasShape() ? Base.HasLength() : Base.iteratorsize(Tit)
Base.iteratoreltype(::MappedKalmanFilter) = Base.HasEltype()
Base.eltype(::Type{MappedKalmanFilter{Tit, TM, F}}) where {Tit, TM, F} = mappedreturntype(TM, F)
mappedreturntype(TM, F) = Any

start(kf::MappedKalmanFilter) = start(kf.it), -Inf, prior(kf.M)
function next(kf::MappedKalmanFilter, state)
    st, s, N = state

    ty, st2 = next(kf.it, st)
    t, y = ty
    t, x, P, Ppred, ll, K = kalman_kernel(s, mean(N), cov(N), t, y, kf.M)
    ret = kf.f(t, x, P, Ppred, ll, K)
    ret, (st2, t, Gaussian(x, P))
end

done(kf::MappedKalmanFilter, state) = done(kf.it, state[1])

length(kf::MappedKalmanFilter) = length(kf.it)


text:
export ARC_step_computation

function ARC_step_computation(h::Float64, g::Float64, Δ::Float64; kwargs...)

    discr1 = (h) ^ 2 - 4.0 * (g / Δ)
    discr2 = (h) ^ 2 + 4.0 * (g / Δ)

    if h > 0.0
      if g > 0.0
        d=(-2.0 * g) / (h + sqrt(discr2))
      else
        d=(-2.0 * g) / (h + sqrt(discr1))
      end
    else
      if g < 0.0
        d=(-h + sqrt(discr1)) / (2.0 / Δ)
      else
        d=(-h + sqrt(discr2)) / (-2.0 / Δ)
      end
    end

    return d
  end


text:
for filename in ["ta_func_api_gen_level0.jl", "ta_func_api_gen_level1.jl",
        "ta_func_api_gen_level2_dataframe.jl", "ta_func_api_gen_level2_timearray.jl"]
    include(filename)
end

function generate_header()
    "# Auto generated file"
end

function generate_footer()
    "# end of auto generated file"
end

function generate_code(d::OrderedDict{Symbol,Any}, _generate_header::Function, _generate_code::Function, _generate_footer::Function)
    s_code = _generate_header()
    s_code *= "\n"
    for symb in keys(d)
        #println("## Generate code for $symb")
        s_code *= _generate_code(d, symb)
    end
    s_code *= "\n"
    s_code *= _generate_footer()
    s_code
end

function generate_all_code(d::OrderedDict{Symbol,Any})
    path = joinpath(basepath(), "generated")
    println("# Generate level 0 code")
    filename = joinpath(path, "ta_func_api_code_level0.jl")
    s_code = generate_code(d, generate_header, generate_ta_func_raw, generate_footer)
    f = open(filename, "w")
    write(f, s_code)
    close(f)

    println("# Generate level 1 code")
    filename = joinpath(path, "ta_func_api_code_level1.jl")
    s_code = generate_code(d, generate_header, generate_ta_func_with_arrays, generate_footer)
    f = open(filename, "w")
    write(f, s_code)
    close(f)

    println("# Generate level 2 code - DataFrames")
    filename = joinpath(path, "ta_func_api_code_level2_dataframes.jl")
    s_code = generate_code(d, generate_header_ta_func_with_dataframes, generate_ta_func_with_dataframes, generate_footer)
    f = open(filename, "w")
    write(f, s_code)
    close(f)

    println("# Generate level 2 code - TimeArrays")
    filename = joinpath(path, "ta_func_api_code_level2_timearrays.jl")
    s_code = generate_code(d, generate_header_ta_func_with_timearray, generate_ta_func_with_timearray, generate_footer)
    f = open(filename, "w")
    write(f, s_code)
    close(f)
end

generate_all_code(D_INDICATORS)