text:
@testset "WKTWriter" begin
    # default writing options
    p = readgeom("POINT(0.12345 2.000 0.1)")
    @test writegeom(p) == "POINT Z (0.12345 2 0.1)"

    p = readgeom("POINT(0.12345 2.000)")
    @test writegeom(p) == "POINT (0.12345 2)"

    # round to 2 decimals
    writer = LibGEOS.WKTWriter(LibGEOS._context, trim=true, outputdim=3, roundingprecision=2)
    @test writegeom(p, writer) == "POINT (0.12 2)"

    # round to 2 decimals and don't trim trailing zeros
    writer = LibGEOS.WKTWriter(LibGEOS._context, trim=false, outputdim=3, roundingprecision=2)
    @test writegeom(p, writer) == "POINT (0.12 2.00)"

    # don't output the Z dimension
    p = readgeom("POINT(0.12345 2.000 0.1)")
    writer = LibGEOS.WKTWriter(LibGEOS._context, trim=false, outputdim=2, roundingprecision=2)
    @test writegeom(p, writer) == "POINT (0.12 2.00)"
end

@testset "GEOS functions" begin
    a = LibGEOS.createCoordSeq(Vector{Float64}[[1,2,3],[4,5,6]])
    b = LibGEOS.cloneCoordSeq(a)
    @test LibGEOS.getCoordinates(b) == LibGEOS.getCoordinates(a)
    a = LibGEOS.createCoordSeq(Vector{Float64}[[1,2,3],[4,5,6],[7,8,9],[10,11,12]])
    b = LibGEOS.cloneCoordSeq(a)
    @test LibGEOS.getCoordinates(b) == LibGEOS.getCoordinates(a)
    LibGEOS.setCoordSeq!(b, 2, [3.0, 3.0, 3.0])
    @test LibGEOS.getCoordinates(a) == Vector{Float64}[[1,2,3],[4,5,6],[7,8,9],[10,11,12]]
    @test LibGEOS.getCoordinates(b) == Vector{Float64}[[1,2,3],[3,3,3],[7,8,9],[10,11,12]]
    c = LibGEOS.createPoint(LibGEOS.createCoordSeq(Vector{Float64}[[1,2]]))
    @test LibGEOS.getCoordinates(LibGEOS.getCoordSeq(c))[1] ≈ [1,2] atol=1e-5

    # Polygons and Holes
    shell = LibGEOS.createLinearRing(Vector{Float64}[[0,0],[10,0],[10,10],[0,10],[0,0]])
    hole1 = LibGEOS.createLinearRing(Vector{Float64}[[1,8],[2,8],[2,9],[1,9],[1,8]])
    hole2 = LibGEOS.createLinearRing(Vector{Float64}[[8,1],[9,1],[9,2],[8,2],[8,1]])
    polygon = LibGEOS.createPolygon(shell,LibGEOS.GEOSGeom[hole1,hole2])
    @test LibGEOS.getGeomDimensions(polygon) == 2
    @test LibGEOS.geomTypeId(polygon) == LibGEOS.GEOS_POLYGON
    @test LibGEOS.geomArea(polygon) ≈ 98.0 atol=1e-5
    exterior = LibGEOS.exteriorRing(polygon)
    @test LibGEOS.getCoordinates(LibGEOS.getCoordSeq(exterior)) == Vector{Float64}[[0,0],[10,0],[10,10],[0,10],[0,0]]
    interiors = LibGEOS.interiorRings(polygon)
    @test LibGEOS.getCoordinates(LibGEOS.getCoordSeq(interiors[1])) == Vector{Float64}[[1,8],[2,8],[2,9],[1,9],[1,8]]
    @test LibGEOS.getCoordinates(LibGEOS.getCoordSeq(interiors[2])) == Vector{Float64}[[8,1],[9,1],[9,2],[8,2],[8,1]]

    # Interpolation and Projection
    ls = LibGEOS.createLineString(Vector{Float64}[[8,1],[9,1],[9,2],[8,2]])
    pt = LibGEOS.interpolate(ls, 2.5)
    coords = LibGEOS.getCoordinates(LibGEOS.getCoordSeq(pt))
    @test length(coords) == 1
    @test coords[1] ≈ [8.5, 2.0] atol=1e-5
    p1 = LibGEOS.createPoint(Float64[10,1])
    p2 = LibGEOS.createPoint(Float64[9,1])
    p3 = LibGEOS.createPoint(Float64[10,0])
    p4 = LibGEOS.createPoint(Float64[9,2])
    p5 = LibGEOS.createPoint(Float64[8.7,1.5])
    dist = LibGEOS.project(ls, p1)
    @test dist ≈ 1 atol=1e-5
    @test LibGEOS.equals(LibGEOS.interpolate(ls, dist), p2)
    @test LibGEOS.project(ls, p2) ≈ 1 atol=1e-5
    @test LibGEOS.project(ls, p3) ≈ 1 atol=1e-5
    @test LibGEOS.project(ls, p4) ≈ 2 atol=1e-5
    @test LibGEOS.project(ls, p5) ≈ 1.5 atol=1e-5

    # Taken from https://svn.osgeo.org/geos/trunk/tests/unit/capi/

    # GEOSContainsTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom2_ = LibGEOS._readgeom("POLYGON EMPTY")
    @test !LibGEOS.contains(geom1_, geom2_)
    @test !LibGEOS.contains(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    geom2_ = LibGEOS._readgeom("POINT(2 2)")
    @test LibGEOS.contains(geom1_, geom2_)
    @test !LibGEOS.contains(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))")
    geom2_ = LibGEOS._readgeom("POLYGON((1 1,1 2,2 2,2 1,1 1))")
    @test LibGEOS.contains(geom1_, geom2_)
    @test !LibGEOS.contains(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    # GEOSConvexHullTest
    input_ = LibGEOS._readgeom("MULTIPOINT (130 240, 130 240, 130 240, 570 240, 570 240, 570 240, 650 240)")
    expected_ = LibGEOS._readgeom("LINESTRING (130 240, 650 240)")
    output_ = LibGEOS.convexhull(input_)
    @test !LibGEOS.isEmpty(output_)
    @test LibGEOS.equals(output_, expected_)
    LibGEOS.destroyGeom(input_)
    LibGEOS.destroyGeom(expected_)
    LibGEOS.destroyGeom(output_)

    # GEOSCoordSeqTest
    cs_ = LibGEOS.createCoordSeq(5, 3)
    @test LibGEOS.getSize(cs_) == 1
    @test LibGEOS.getDimensions(cs_) == 2
    @test LibGEOS.getCoordinates(cs_) == [[5.0, 3.0]]

    cs_2 = LibGEOS.createCoordSeq([5.0, 3.0])
    @test LibGEOS.getSize(cs_2) == 1
    @test LibGEOS.getDimensions(cs_2) == 2
    @test LibGEOS.getCoordinates(cs_2) == [[5.0, 3.0]]

    cs_3 = LibGEOS.createCoordSeq(Vector{Float64}[[5.0, 3.0], [1.0, 2.0], [1.0, 3.0]])
    @test LibGEOS.getSize(cs_3) == 3
    @test LibGEOS.getDimensions(cs_3) == 2
    @test LibGEOS.getCoordinates(cs_3)[1] == [5.0, 3.0]
    @test LibGEOS.getCoordinates(cs_3)[2] == [1.0, 2.0]
    @test LibGEOS.getCoordinates(cs_3)[3] == [1.0, 3.0]
    @test LibGEOS.getCoordinates(cs_3, 1) == [5.0, 3.0]
    @test LibGEOS.getCoordinates(cs_3, 2) == [1.0, 2.0]
    @test LibGEOS.getCoordinates(cs_3, 3) == [1.0, 3.0]

    for i=1:5
        x = i*10.0
        y = i*10.0+1.0
        z = i*10.0+2.0

        LibGEOS.setX!(cs_, 1, x)
        LibGEOS.setY!(cs_, 1, y)
        LibGEOS.setZ!(cs_, 1, z)
        @test LibGEOS.getX(cs_, 1) ≈ x atol=1e-5
        @test LibGEOS.getY(cs_, 1) ≈ y atol=1e-5
        @test LibGEOS.getZ(cs_, 1) ≈ z atol=1e-5
    end

    cs_ = LibGEOS.createCoordSeq(1, ndim=3)
    @test LibGEOS.getSize(cs_) == 1
    @test LibGEOS.getDimensions(cs_) == 3
    x,y,z = 10.0, 11.0, 12.0

    LibGEOS.setX!(cs_, 1, x)
    LibGEOS.setY!(cs_, 1, y)
    LibGEOS.setZ!(cs_, 1, z)
    @test LibGEOS.getX(cs_, 1) ≈ x atol=1e-5
    @test LibGEOS.getY(cs_, 1) ≈ y atol=1e-5
    @test LibGEOS.getZ(cs_, 1) ≈ z atol=1e-5

    cs_ = LibGEOS.createCoordSeq(1, ndim=3)
    @test LibGEOS.getSize(cs_) == 1
    @test LibGEOS.getDimensions(cs_) == 3
    x,y,z = 10.0, 11.0, 12.0

    LibGEOS.setX!(cs_, 1, x)
    LibGEOS.setY!(cs_, 1, y)
    LibGEOS.setZ!(cs_, 1, z)
    @test LibGEOS.getX(cs_, 1) ≈ x atol=1e-5
    @test LibGEOS.getY(cs_, 1) ≈ y atol=1e-5
    @test LibGEOS.getZ(cs_, 1) ≈ z atol=1e-5

    # LibGEOS.delaunayTriangulationTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    @test LibGEOS.isEmpty(geom1_)
    geom2_ = LibGEOS.delaunayTriangulation(geom1_,0.0,true)
    @test LibGEOS.isEmpty(geom2_)
    @test LibGEOS.geomTypeId(geom2_) == LibGEOS.GEOS_MULTILINESTRING
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("POINT(0 0)")
    geom2_ = LibGEOS.delaunayTriangulation(geom1_, 0.0, false)
    @test LibGEOS.isEmpty(geom2_)
    @test LibGEOS.geomTypeId(geom2_) == LibGEOS.GEOS_GEOMETRYCOLLECTION
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("MULTIPOINT(0 0, 5 0, 10 0)")
    geom2_ = LibGEOS.delaunayTriangulation(geom1_, 0.0, false)
    @test LibGEOS.isEmpty(geom2_)
    @test LibGEOS.geomTypeId(geom2_) == LibGEOS.GEOS_GEOMETRYCOLLECTION
    LibGEOS.destroyGeom(geom2_)
    geom2_ = LibGEOS.delaunayTriangulation(geom1_, 0.0, true)
    geom3_ = LibGEOS._readgeom("MULTILINESTRING ((5 0, 10 0), (0 0, 5 0))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("MULTIPOINT(0 0, 10 0, 10 10, 11 10)")
    geom2_ = LibGEOS.delaunayTriangulation(geom1_, 2.0, true)
    geom3_ = LibGEOS._readgeom("MULTILINESTRING ((0 0, 10 10), (0 0, 10 0), (10 0, 10 10))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    # GEOSDistanceTest
    geom1_ = LibGEOS._readgeom("POINT(10 10)")
    geom2_ = LibGEOS._readgeom("POINT(3 6)")
    @test LibGEOS.geomDistance(geom1_, geom2_) ≈ 8.06225774829855 atol=1e-12
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    # GEOSGeom_extractUniquePointsTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom2_ = LibGEOS.uniquePoints(geom1_)
    @test LibGEOS.isEmpty(geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("MULTIPOINT(0 0, 0 0, 1 1)")
    geom2_ = LibGEOS._readgeom("MULTIPOINT(0 0, 1 1)")
    geom3_ = LibGEOS.uniquePoints(geom1_)
    @test LibGEOS.equals(geom3_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))")
    geom2_ = LibGEOS._readgeom("MULTIPOINT(0 0, 1 1, 2 2, 5 5, 0 2)")
    geom3_ = LibGEOS.uniquePoints(geom1_)
    @test LibGEOS.equals(geom3_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    # GEOSGetCentroidTest
    geom1_ = LibGEOS._readgeom("POINT(10 0)")
    geom2_ = LibGEOS.centroid(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (10 0)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 10 0)")
    geom2_ = LibGEOS.centroid(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (5 0)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
    geom2_ = LibGEOS.centroid(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (5 5)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("LINESTRING EMPTY")
    geom2_ = LibGEOS.centroid(geom1_)
    geom3_ = LibGEOS._readgeom("POINT EMPTY")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    # GEOSIntersectionTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom2_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom3_ = LibGEOS.intersection(geom1_, geom2_)
    geom4_ = LibGEOS._readgeom("GEOMETRYCOLLECTION EMPTY")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    geom2_ = LibGEOS._readgeom("POINT(2 2)")
    geom3_ = LibGEOS.intersection(geom1_, geom2_)
    @test LibGEOS.equals(geom3_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))")
    geom2_ = LibGEOS._readgeom("POLYGON((-1 1,-1 2,2 2,2 1,-1 1))")
    geom3_ = LibGEOS.intersection(geom1_, geom2_)
    geom4_ = LibGEOS._readgeom("POLYGON ((0 1, 0 2, 2 2, 2 1, 0 1))")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("MULTIPOLYGON(((0 0,5 10,10 0,0 0),(1 1,1 2,2 2,2 1,1 1),(100 100,100 102,102 102,102 100,100 100)))")
    geom2_ = LibGEOS._readgeom("POLYGON((0 1,0 2,10 2,10 1,0 1))")
    geom3_ = LibGEOS.intersection(geom1_, geom2_)
    geom4_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (LINESTRING (1 2, 2 2), LINESTRING (2 1, 1 1), POLYGON ((0.5 1, 1 2, 1 1, 0.5 1)), POLYGON ((9 2, 9.5 1, 2 1, 2 2, 9 2)))")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    # GEOSIntersectsTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom2_ = LibGEOS._readgeom("POLYGON EMPTY")
    @test !LibGEOS.intersects(geom1_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    geom2_ = LibGEOS._readgeom("POINT(2 2)")
    @test LibGEOS.intersects(geom1_, geom2_)
    @test LibGEOS.intersects(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))")
    geom2_ = LibGEOS._readgeom("POLYGON((1 1,1 2,2 2,2 1,1 1))")
    @test LibGEOS.intersects(geom1_, geom2_)
    @test LibGEOS.intersects(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    # LineString_PointTest
    geom1 = LibGEOS._readgeom("LINESTRING(0 0, 5 5, 10 10)")
    @test !LibGEOS.isClosed(geom1)
    @test LibGEOS.geomTypeId(geom1) == LibGEOS.GEOS_LINESTRING
    @test LibGEOS.numPoints(geom1) == 3
    @test LibGEOS.geomLength(geom1) ≈ sqrt(100 + 100) atol=1e-5
    geom2 = LibGEOS.getPoint(geom1, 1)
    @test LibGEOS.getGeomX(geom2) ≈ 0.0 atol=1e-5
    @test LibGEOS.getGeomY(geom2) ≈ 0.0 atol=1e-5
    geom2 = LibGEOS.getPoint(geom1, 2)
    @test LibGEOS.getGeomX(geom2) ≈ 5.0 atol=1e-5
    @test LibGEOS.getGeomY(geom2) ≈ 5.0 atol=1e-5
    LibGEOS.destroyGeom(geom2)
    geom2 = LibGEOS.getPoint(geom1, 3)
    @test LibGEOS.getGeomX(geom2) ≈ 10.0 atol=1e-5
    @test LibGEOS.getGeomY(geom2) ≈ 10.0 atol=1e-5
    geom2 = LibGEOS.startPoint(geom1)
    @test LibGEOS.getGeomX(geom2) ≈ 0.0 atol=1e-5
    @test LibGEOS.getGeomY(geom2) ≈ 0.0 atol=1e-5
    LibGEOS.destroyGeom(geom2)
    geom2 = LibGEOS.endPoint(geom1)
    @test LibGEOS.getGeomX(geom2) ≈ 10.0 atol=1e-5
    @test LibGEOS.getGeomY(geom2) ≈ 10.0 atol=1e-5
    LibGEOS.destroyGeom(geom2)

    # GEOSNearestPointsTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom2_ = LibGEOS._readgeom("POLYGON EMPTY")
    @test LibGEOS.nearestPoints(geom1_, geom2_) == C_NULL
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    geom2_ = LibGEOS._readgeom("POLYGON((8 8, 9 9, 9 10, 8 8))")
    coords_ = LibGEOS.nearestPoints(geom1_, geom2_)
    @test LibGEOS.getSize(coords_) == 2
    @test LibGEOS.getCoordinates(coords_,1)[1:2] ≈ [5.0,5.0] atol=1e-5
    @test LibGEOS.getCoordinates(coords_,2)[1:2] ≈ [8.0,8.0] atol=1e-5
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyCoordSeq(coords_)

    # GEOSNodeTest
    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 10 10, 10 0, 0 10)")
    geom2_ = LibGEOS.node(geom1_)
    LibGEOS.normalize!(geom2_)
    geom3_ = LibGEOS._readgeom("MULTILINESTRING ((5 5, 10 0, 10 10, 5 5), (0 10, 5 5), (0 0, 5 5))")
    @test equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("MULTILINESTRING((0 0, 2 0, 4 0),(5 0, 3 0, 1 0))")
    geom2_ = LibGEOS.node(geom1_)
    LibGEOS.normalize!(geom2_)
    geom3_ = LibGEOS._readgeom("MULTILINESTRING ((4 0, 5 0), (3 0, 4 0), (2 0, 3 0), (1 0, 2 0), (0 0, 1 0))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("MULTILINESTRING((0 0, 2 0, 4 0),(0 0, 2 0, 4 0))")
    geom2_ = LibGEOS.node(geom1_)
    LibGEOS.normalize!(geom2_)
    geom3_ = LibGEOS._readgeom("MULTILINESTRING ((2 0, 4 0), (0 0, 2 0))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    # GEOSPointOnSurfaceTest
    geom1_ = LibGEOS._readgeom("POINT(10 0)")
    geom2_ = LibGEOS.pointOnSurface(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (10 0)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 5 0, 10 0)")
    geom2_ = LibGEOS.pointOnSurface(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (5 0)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
    geom2_ = LibGEOS.pointOnSurface(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (5 5)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom(
    """POLYGON((
    56.528666666700 25.2101666667,
    56.529000000000 25.2105000000,
    56.528833333300 25.2103333333,
    56.528666666700 25.2101666667))""")
    geom2_ = LibGEOS.pointOnSurface(geom1_)
    @test LibGEOS.getCoordinates(LibGEOS.getCoordSeq(geom2_))[1] ≈ [56.5286666667, 25.2101666667] atol=1e-5
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("LINESTRING EMPTY")
    geom2_ = LibGEOS.pointOnSurface(geom1_)
    geom3_ = LibGEOS._readgeom("POINT EMPTY")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 0 0)")
    geom2_ = LibGEOS.pointOnSurface(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (0 0)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    # GEOSPreparedGeometryTest

    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyPreparedGeom(prepGeom1_)

    geom1_ = LibGEOS._readgeom("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")
    geom2_ = LibGEOS._readgeom("POLYGON((2 2, 2 3, 3 3, 3 2, 2 2))")
    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
    @test LibGEOS.prepcontainsproperly(prepGeom1_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyPreparedGeom(prepGeom1_)

    geom1_ = LibGEOS._readgeom("POLYGON((2 2, 2 3, 3 3, 3 2, 2 2))")
    geom2_ = LibGEOS._readgeom("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")
    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
    @test !LibGEOS.prepcontainsproperly(prepGeom1_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyPreparedGeom(prepGeom1_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 10 10)")
    geom2_ = LibGEOS._readgeom("LINESTRING(0 10, 10 0)")
    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
    @test LibGEOS.prepintersects(prepGeom1_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyPreparedGeom(prepGeom1_)

    geom1_ = LibGEOS._readgeom("POLYGON((0 0, 0 10, 10 11, 10 0, 0 0))")
    geom2_ = LibGEOS._readgeom("POLYGON((0 0, 2 0, 2 2, 0 2, 0 0))")
    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
    @test LibGEOS.prepcovers(prepGeom1_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyPreparedGeom(prepGeom1_)

    geom1_ = LibGEOS._readgeom("POLYGON((0 0, 0 10, 10 11, 10 0, 0 0))")
    geom2_ = LibGEOS._readgeom("POLYGON((0 0, 2 0, 2 2, 0 2, 0 0))")
    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
    @test LibGEOS.prepcontains(prepGeom1_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyPreparedGeom(prepGeom1_)

    # GEOSSharedPathsTest
    geom1_ = LibGEOS._readgeom("LINESTRING (-30 -20, 50 60, 50 70, 50 0)")
    geom2_ = LibGEOS._readgeom("LINESTRING (-29 -20, 50 60, 50 70, 51 0)")
    geom3_ = LibGEOS.sharedPaths(geom1_, geom2_)
    geom4_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (MULTILINESTRING ((50 60, 50 70)), MULTILINESTRING EMPTY)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    # /// Illegal case (point-poly)
    # void object::test<1>()
    # {
    #     geom1_ = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")
    #     geom2_ = GEOSGeomFromWKT("POINT(0.5 0)")
    #     geom3_ = GEOSSharedPaths(geom1_, geom2_)

    #     ensure(!geom3_)
    # }

    # http://trac.osgeo.org/postgis/ticket/670#comment:3
    # geom1_ = LibGEOS._readgeom("POINT(-11.1111111 40)")
    # geom2_ = LibGEOS._readgeom("POLYGON((-8.1111111 60,-8.16875525879031 59.4147290339516,-8.33947250246614 58.8519497029047,-8.61670226309236 58.3332893009412,-8.98979075644036 57.8786796564404,-9.44440040094119 57.5055911630924,-9.96306080290473 57.2283614024661,-10.5258401339516 57.0576441587903,-11.1111111 57,-11.6963820660484 57.0576441587903,-12.2591613970953 57.2283614024661,-12.7778217990588 57.5055911630924,-13.2324314435596 57.8786796564404,-13.6055199369076 58.3332893009412,-13.8827496975339 58.8519497029047,-14.0534669412097 59.4147290339516,-14.1111111 60,-14.0534669412097 60.5852709660484,-13.8827496975339 61.1480502970953,-13.6055199369076 61.6667106990588,-13.2324314435597 62.1213203435596,-12.7778217990588 62.4944088369076,-12.2591613970953 62.7716385975339,-11.6963820660484 62.9423558412097,-11.1111111 63,-10.5258401339516 62.9423558412097,-9.96306080290474 62.7716385975339,-9.4444004009412 62.4944088369076,-8.98979075644036 62.1213203435596,-8.61670226309237 61.6667106990588,-8.33947250246614 61.1480502970953,-8.16875525879031 60.5852709660484,-8.1111111 60))")
    # geom3_ = LibGEOS.sharedPaths(geom1_, geom2_)
    # ensure(!geom3_)

    # GEOSSimplifyTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    @test LibGEOS.isEmpty(geom1_)
    geom2_ = LibGEOS.simplify(geom1_, 43.2)
    @test LibGEOS.isEmpty(geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    # GEOSSnapTest

    geom1_ = LibGEOS._readgeom("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")
    geom2_ = LibGEOS._readgeom("POINT(0.5 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 1.0)
    geom4_ = LibGEOS._readgeom("POLYGON ((0.5 0, 10 0, 10 10, 0 10, 0.5 0))")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING (-30 -20, 50 60, 50 0)")
    geom2_ = LibGEOS._readgeom("LINESTRING (-29 -20, 40 60, 51 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 2.0)
    geom4_ = LibGEOS._readgeom("LINESTRING (-29 -20, 50 60, 51 0)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING (-20 -20, 50 50, 100 100)")
    geom2_ = LibGEOS._readgeom("LINESTRING (-10 -9, 40 20, 80 79)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 2.0)
    geom4_ = LibGEOS._readgeom("LINESTRING (-20 -20, -10 -9, 50 50, 80 79, 100 100)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 10 0)")
    geom2_ = LibGEOS._readgeom("LINESTRING(0 0, 9 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 2.0)
    @test LibGEOS.equals(geom3_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0, 10 0)")
    geom2_ = LibGEOS._readgeom("LINESTRING(0 0, 9 0, 10 0, 11 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 2.0)
    @test LibGEOS.equals(geom3_, geom2_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 3,4 1,0 1)")
    geom2_ = LibGEOS._readgeom("MULTIPOINT(5 0,4 1)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 2.0)
    geom4_ = LibGEOS._readgeom("LINESTRING (0 3, 4 1, 5 0, 0 1)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 3,4 1,0 1)")
    geom2_ = LibGEOS._readgeom("MULTIPOINT(4 1,5 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 2.0)
    geom4_ = LibGEOS._readgeom("LINESTRING (0 3, 4 1, 5 0, 0 1)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 0,10 0,10 10,0 10,0 0)")
    geom2_ = LibGEOS._readgeom("MULTIPOINT(0 0,-1 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 3.0)
    geom4_ = LibGEOS._readgeom("LINESTRING (-1 0, 0 0, 10 0, 10 10, 0 10, -1 0)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING(0 2,5 2,9 2,5 0)")
    geom2_ = LibGEOS._readgeom("POINT(5 0)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 3.0)
    geom4_ = LibGEOS._readgeom("LINESTRING (0 2, 5 2, 9 2, 5 0)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    geom1_ = LibGEOS._readgeom("LINESTRING(-71.1317 42.2511,-71.1317 42.2509)")
    geom2_ = LibGEOS._readgeom("MULTIPOINT(-71.1261 42.2703,-71.1257 42.2703,-71.1261 42.2702)")
    geom3_ = LibGEOS.snap(geom1_, geom2_, 0.5)
    geom4_ = LibGEOS._readgeom("LINESTRING (-71.1257 42.2703, -71.1261 42.2703, -71.1261 42.2702, -71.1317 42.2509)")
    @test LibGEOS.equals(geom3_, geom4_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(geom4_)

    # GEOSUnaryUnionTest
    geom1_ = LibGEOS._readgeom("POINT EMPTY")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("POINT EMPTY")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("POINT (6 3)")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("POINT (6 3)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("POINT (4 5 6)")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("POINT Z (4 5 6)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("MULTIPOINT (4 5, 6 7, 4 5, 6 5, 6 7)")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("MULTIPOINT (4 5, 6 5, 6 7)")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (POINT(4 5), MULTIPOINT(6 7, 6 5, 6 7), LINESTRING(0 5, 10 5), LINESTRING(4 -10, 4 10))")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (POINT (6 7), LINESTRING (0 5, 4 5), LINESTRING (4 5, 10 5), LINESTRING (4 -10, 4 5), LINESTRING (4 5, 4 10))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (POINT(4 5), MULTIPOINT(6 7, 6 5, 6 7), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)))")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (POINT (6 7), POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (5 6, 7 6, 7 8, 5 8, 5 6)))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5)), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)))")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (LINESTRING (5 7, 7 7), LINESTRING (10 7, 12 7), LINESTRING (5.5 7.5, 6.5 7.5), POLYGON ((10 7, 10 0, 0 0, 0 10, 10 10, 10 7), (5 6, 7 6, 7 7, 7 8, 5 8, 5 7, 5 6)))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    geom1_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5)), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)), MULTIPOINT(6 6.5, 6 1, 12 2, 6 1))")
    geom2_ = LibGEOS.unaryUnion(geom1_)
    geom3_ = LibGEOS._readgeom("GEOMETRYCOLLECTION (POINT (6 6.5), POINT (12 2), LINESTRING (5 7, 7 7), LINESTRING (10 7, 12 7), LINESTRING (5.5 7.5, 6.5 7.5), POLYGON ((10 7, 10 0, 0 0, 0 10, 10 10, 10 7), (5 6, 7 6, 7 7, 7 8, 5 8, 5 7, 5 6)))")
    @test LibGEOS.equals(geom2_, geom3_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)

    # GEOSWithinTest
    geom1_ = LibGEOS._readgeom("POLYGON EMPTY")
    geom2_ = LibGEOS._readgeom("POLYGON EMPTY")
    @test !LibGEOS.within(geom1_, geom2_)
    @test !LibGEOS.within(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    geom2_ = LibGEOS._readgeom("POINT(2 2)")
    @test !LibGEOS.within(geom1_, geom2_)
    @test LibGEOS.within(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    geom1_ = LibGEOS._readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))")
    geom2_ = LibGEOS._readgeom("POLYGON((1 1,1 2,2 2,2 1,1 1))")
    @test !LibGEOS.within(geom1_, geom2_)
    @test LibGEOS.within(geom2_, geom1_)
    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)

    # GEOSisClosedTest
    geom_ = LibGEOS._readgeom("LINESTRING(0 0, 1 0, 1 1)")
    @test !LibGEOS.isClosed(geom_)
    LibGEOS.destroyGeom(geom_)
    geom_ = LibGEOS._readgeom("LINESTRING(0 0, 0 1, 1 1, 0 0)")
    @test LibGEOS.isClosed(geom_)
    LibGEOS.destroyGeom(geom_)

    # setPrecision, getPrecision
    # Taken from https://git.osgeo.org/gitea/geos/geos/src/branch/master/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp

    geom1_ = readgeom("LINESTRING(2 10, 4 30)")
    geom2_ = readgeom("LINESTRING(4 10, 2 30)")
    geom3_ = intersection(geom1_, geom2_)
    @test equals(geom3_, readgeom("POINT(3 20)"))

    g = setPrecision(geom1_, 2.0)
    LibGEOS.destroyGeom(geom1_)
    geom1_ = g
    g = setPrecision(geom2_, 2.0)
    LibGEOS.destroyGeom(geom2_)
    geom2_ = g
    geom3_ = intersection(geom1_, geom2_)
    @test equals(geom3_, readgeom("POINT(4 20)"))

    g = setPrecision(geom1_, 0.5)
    LibGEOS.destroyGeom(geom1_)
    geom1_ = g
    geom3_ = intersection(geom1_, geom2_)
    @test equals(geom3_, readgeom("POINT(3 20)"))
    @test getPrecision(geom1_) == 0.5
    @test getPrecision(geom2_) == 2.0

    geom1_ = readgeom("POLYGON((10 10,20 10,16 15,20 20, 10 20, 14 15, 10 10))")
    geom2_ = setPrecision(geom1_, 5.0)
    @test equals(geom2_, readgeom(
        "MULTIPOLYGON (((10 10, 15 15, 20 10, 10 10)), ((15 15, 10 20, 20 20, 15 15)))"))
    geom3_ = setPrecision(geom1_, 5.0, flags = LibGEOS.GEOS_PREC_NO_TOPO)
    @test equals(geom3_, readgeom(
        "POLYGON ((10 10, 20 10, 15 15, 20 20, 10 20, 15 15, 10 10))"))

    geom1_ = readgeom("LINESTRING(1 0, 2 0)")
    geom2_ = setPrecision(geom1_, 5.0)
    @test equals(geom2_, readgeom("LINESTRING EMPTY"))
    geom3_ = setPrecision(geom1_, 5.0, flags = LibGEOS.GEOS_PREC_KEEP_COLLAPSED)
    # @test equals(geom3_, readgeom("LINESTRING (0 0, 0 0)")) # false ??
    @test writegeom(geom3_) == "LINESTRING (0 0, 0 0)"

    LibGEOS.destroyGeom(geom1_)
    LibGEOS.destroyGeom(geom2_)
    LibGEOS.destroyGeom(geom3_)
    LibGEOS.destroyGeom(g)

    # GEOSMinimumRotatedRectangleTest
    input_ = readgeom("POLYGON ((1 6, 6 11, 11 6, 6 1, 1 6))")
    output_ = minimumRotatedRectangle(input_)
    expected_ = readgeom("POLYGON ((6 1, 11 6, 6 11, 1 6, 6 1))")
    @test LibGEOS.equals(output_, expected_)
    LibGEOS.destroyGeom(input_)
    LibGEOS.destroyGeom(output_)
    LibGEOS.destroyGeom(expected_)
end


text:
mutable struct SGD{F} <: Optimizer
    target_layers :: F
    lr
end

SGD(target_layers; lr=1e-3) = SGD(target_layers, lr)

function do_optimize!(model::T, optimizer::SGD) where T <: Model
    for layer in optimizer.target_layers(model)
        for p in parameters(layer)
            p.values -= optimizer.lr * p.grad.values
        end
    end
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 736856ce-f490-11eb-3349-057c86edfe7e
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using MLJ, MLJLinearModels, Plots, DataFrames, Random
end

# ╔═╡ c693088f-7f80-4cdd-b9b5-65a50da732ac
begin
    using MLCourse
    import MLCourse: poly
    MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ e1077092-f72a-42af-b6e0-a616f938cba8
using PlutoUI; PlutoUI.TableOfContents()

# ╔═╡ 7eb6a060-f948-4d85-881a-4909e74c15bd
md"# Validation Set Approach

In the following cell we define a `data_generator` that creates data sets with a polynomial relationship of degree 3 between the input x and the average output f(x). We split these data sets into a training set, a validation set and a test set. The training set will be used to find the parameters for a machine learning method with given hyper-parameters, the validation set will be used to find the hyper-parameters and the test set will be used to evaluate all our models.
"

# ╔═╡ 7dd7e9a7-9245-4c64-af0c-8f7d2f62b2bf
begin
    f(x) = 0.3 + 2x - 0.8x^2 - 0.4x^3
    function data_generator(; n = 500, seed = 12)
        rng = MersenneTwister(seed)
        x = randn(rng, n)
        DataFrame(x = x, y = f.(x) .+ randn(rng, n))
    end
    function data_split(data;
                        shuffle = false,
                        idx_train = 1:50,
                        idx_valid = 51:100,
                        idx_test = 101:500)
        idxs = if shuffle
                randperm(size(data, 1))
            else
                1:size(data, 1)
            end
        (train = data[idxs[idx_train], :],
         valid = data[idxs[idx_valid], :],
         test = data[idxs[idx_test], :])
    end
    function fit_and_evaluate(model, data)
        mach = fit!(machine(model, select(data.train, :x), data.train.y),
                    verbosity = 0)
        (train = rmse(predict(mach, select(data.train, :x)), data.train.y),
         valid = rmse(predict(mach, select(data.valid, :x)), data.valid.y),
         test = rmse(predict(mach, select(data.test, :x)), data.test.y))
    end
end;

# ╔═╡ fa9b4b9e-4e97-4e5c-865d-ad3ef288e4cf
data1 = data_split(data_generator(seed = 1))

# ╔═╡ 751880ec-1a82-4142-b875-177d436bbc72
md"To simplify our training and test procedure we define a new type of a machine called `PolynomialRegressor` that takes as keyword argument the `degree` of the polynomial. You do not need to understand the code in the following hidden cell. But if you are interested in how one can define custom machines in the `MLJ` framework you are welcome to look at the code and the [MLJ Documentation](https://alan-turing-institute.github.io/MLJ.jl/dev/simple_user_defined_models/); it is actually quite simple to write a custom machine."

# ╔═╡ da2bbb4d-d45d-4d2a-9c8a-21395bcbb851
begin
    # Use the Live Docs to learn more about Base.@kwdef
    Base.@kwdef mutable struct PolynomialRegressor{T} <: Deterministic
        degree::Int = 3
        regressor::T = LinearRegressor()
    end
    function MLJ.MLJBase.fit(model::PolynomialRegressor, verbosity, X, y)
        Xpoly = poly(X, model.degree) # here we transform the input
        MLJ.MLJBase.fit(model.regressor, verbosity, Xpoly, y)
    end
    function MLJ.MLJBase.predict(model::PolynomialRegressor, fitresult, Xnew)
        Xpoly = poly(Xnew, model.degree) # here we transform the input
        MLJ.MLJBase.predict(model.regressor, fitresult, Xpoly)
    end
end

# ╔═╡ 2e02ac0d-c4d0-47ba-be57-445adeb6ab8b
losses1 = [fit_and_evaluate(PolynomialRegressor(; degree), data1) for degree in 1:10]

# ╔═╡ 91eecd2b-af18-4a63-9684-28950e604d1a
let validlosses = getproperty.(losses1, :valid), i = argmin(validlosses)
    plot(1:10, validlosses, label = nothing)
    scatter!([i], [validlosses[i]], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ b45b32b1-8a65-4823-b3bb-f0b7cc57604b
md"For this seed of the random number generator the optimal degree (x coordinate of the red point in the figure above) found by the validation set approach is close to the actual degree of the data generating process.
Let us now look at other seeds."

# ╔═╡ ac3c7e84-6c47-4dc1-b862-de4cfb05dad9
losses = [fit_and_evaluate(PolynomialRegressor(; degree),
                           data_split(data_generator(seed = seed)))
          for degree in 1:10, seed in 1:20]

# ╔═╡ 26362233-b006-423d-8fb5-7cd9150405b4
let validlosses = getproperty.(losses, :valid), i = argmin(validlosses, dims = 1)
    plot(validlosses, label = nothing, ylims = (0, 10))
    scatter!((x -> x[1]).(i), validlosses[i], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 5cff19c7-3426-4d91-b555-f059f5f41886
md"Here we see a large variability of the optimal degree (x coordinates of all the points in the figure above).

In the following cell we compare the validation set estimate of the test error to the test set estimation of the test error. We can confirm that the validation set estimate is much lower than the test set estimate."

# ╔═╡ df16313f-76be-4d8b-88e7-7cc618ff49d4
let validlosses = getproperty.(losses, :valid),
	testlosses = getproperty.(losses, :test),
	winners = argmin(validlosses, dims = 1)
	(mean_valid_winner = mean(validlosses[winners]),
	 mean_test_winner = mean(testlosses[winners]))
end

# ╔═╡ 46bffe01-1a7c-405b-b3ec-bfe7570b8a3c
md"# Cross-Validation

In the following cell we implement cross-validation ourselves."

# ╔═╡ b3aed705-5072-4d9b-bb8f-865ac1561bf6
begin
    function cross_validation_sets(idx, K)
        n = length(idx)
        r = n ÷ K
        [let idx_valid = idx[(i-1)*r+1:(i == K ? n : i*r)]
             (idx_valid = idx_valid, idx_train = setdiff(idx, idx_valid))
         end
         for i in 1:K]
    end
    function cross_validation(model, data; K = 5)
        losses = [fit_and_evaluate(model,
                                   data_split(data; idxs...))
                  for idxs in cross_validation_sets(1:100, K)]
        (train = mean(getproperty.(losses, :train)),
         valid = mean(getproperty.(losses, :valid)),
         test = mean(getproperty.(losses, :test)))
    end
end

# ╔═╡ 6916273a-d116-4173-acaa-2bcac1d1753b
cross_validation_sets(1:100, 4)

# ╔═╡ c14c9e49-a993-456e-9115-97da86f8e498
losses1_cv10 = [cross_validation(PolynomialRegressor(; degree),
		                         data_generator(seed = 1),
		                         K = 10) for degree in 1:10]

# ╔═╡ 4b435518-2f12-4921-bb1f-fdd049ddfaed
let validlosses = getproperty.(losses1_cv10, :valid), i = argmin(validlosses)
    plot(1:10, validlosses, label = nothing)
    scatter!([i], [validlosses[i]], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 1e584a38-2fef-4877-87f6-92237d71c4b3
losses_cv10 = [cross_validation(PolynomialRegressor(; degree),
                                data_generator(seed = seed),
                                K = 10)
               for degree in 1:10, seed in 1:20]

# ╔═╡ a7c88b3f-92cb-4253-a889-c78683722c1d
let validlosses = getproperty.(losses_cv10, :valid), i = argmin(validlosses, dims = 1)
    plot(validlosses, label = nothing, ylims = (0, 10))
    scatter!((x -> x[1]).(i), validlosses[i], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 0ab1238e-8279-43d5-9803-848cb00156c2
let validlosses = getproperty.(losses_cv10, :valid),
	testlosses = getproperty.(losses_cv10, :test),
	winners = argmin(validlosses, dims = 1)
	(mean_valid_winner = mean(validlosses[winners]),
	 mean_test_winner = mean(testlosses[winners]))
end

# ╔═╡ eaacf529-1727-4743-941b-360c53088b1d
md"Instead of our own cross-validation function we can also use the builtin functions of MLJ. The `train_test_pairs` function does basically the same as our `cross_validation_sets` function."

# ╔═╡ 967a6e08-f0e6-45b7-988b-d0df237f3ddf
MLJ.MLJBase.train_test_pairs(CV(nfolds = 10), 1:100)

# ╔═╡ 535a11d9-833e-4237-949e-9a3e289c600b
md"`MLJ` has the very useful function `evaluate!`. Have a look at the Live docs to learn more about this function. In the following you see an example where we perform 10-fold cross-validation on the `rmse` measure to estimate the expected test error for a polynomial regressor of degree 4."

# ╔═╡ 8ae790b3-3987-4f41-8e21-adbb71081eb9
let data = data_generator(seed = 1, n = 100)
    evaluate!(machine(PolynomialRegressor(degree = 4), select(data, :x), data.y),
              resampling = CV(nfolds = 10), measure = rmse)
end

# ╔═╡ 134df37f-b737-4c80-a5f9-1149aeec970c
md"We can use this now to find the best degree with 10-fold cross-validation for 100 different seeds. You will see that - even though not perfect - the hyper-parameters found with 10-fold cross-validation have lower variance and are usually closer to the true value than with the validation set approach."

# ╔═╡ abb71af3-8aae-4806-9d5a-d144c15d22ef
losses_mlj_cv10 = [let data = data_generator(seed = seed)[1:100, :]
                       evaluate!(machine(PolynomialRegressor(; degree),
                                         select(data, :x),
                                         data.y),
                                resampling = CV(nfolds = 10),
                                measure = rmse,
                                verbosity = 0).measurement[]
                   end
                   for degree in 1:10, seed in 1:20]

# ╔═╡ 6599d2b4-68f4-4c22-8e40-bf3722597692
let validlosses  = losses_mlj_cv10, i = argmin(validlosses, dims = 1)
    plot(validlosses, label = nothing, ylims = (0, 10))
    scatter!((x -> x[1]).(i), validlosses[i], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 29683c99-6a6a-4f65-bea2-d592895d887e
md"# Model Tuning

Finding good hyper-parameters (tuning) is such an important step in the process of finding good machine learning models that there exist some nice utility functions to tune the hyper-parameters ([tuning section in the MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/)). In the cell below you see an example where a `PolynomialRegressor` is automatically tuned on a given dataset by performing 10-fold cross-validation on all degrees on a \"grid\", i.e. of all degrees from 1 to 17 are tested."

# ╔═╡ f93f20db-4fed-481f-b085-ca744b68fa8f
begin
    model = PolynomialRegressor()
    data2 = data_generator(seed = 2, n = 100)
    self_tuning_model = TunedModel(model = model,
                                   resampling = CV(nfolds = 10),
                                   tuning = Grid(),
                                   range = range(model, :degree, values = 1:17),
                                   measure = rmse)
    self_tuning_mach = machine(self_tuning_model, select(data2, :x), data2.y) |> fit!
end

# ╔═╡ d72609a1-f93f-4ca6-8759-727662233e97
md"The self-tuned model actually found the true degree 3 for this data set, as we can see from the report of the tuned machine. The estimated test error can also be found in the report under `best_history_entry.measurement`."

# ╔═╡ 59245b3e-ddfc-46c4-ba44-86ce191672ae
report(self_tuning_mach)

# ╔═╡ f8b48f50-09cb-498a-89d1-9ac9b5722d0c
md"The result of the tuned machine looks quite good here: the reducible error (average distance between red and green curve) looks small."

# ╔═╡ 47a2c14b-a02d-43d7-ac8a-4d95a0d91fa8
let x = -3:.1:3
    scatter(data2.x, data2.y, label = "data", legend = :bottomleft)
    plot!(f, label = "generator", w = 2, xlims = (-3, 3))
    plot!(x, predict(self_tuning_mach, DataFrame(x = x)), label = "self tuning fit", w = 2)
end

# ╔═╡ ee89c448-1e69-4fd1-a4b8-7297a09f2685
md"# Exercises

## Conceptual
1. We review k-fold cross-validation.
    - Explain how k-fold cross-validation is implemented.
    - What are the advantages and disadvantages of k-fold cross-validation relative to:
        - The validation set approach?
        - LOOCV?
2. Suppose you receive the following email of a colleague of yours. Write an answer to this email.
```
    Hi

    In my internship I am doing this research project with a company that
    manufactures constituents for medical devices. We need to assure high
    quality and detect the pieces that are not fully functional. To do so we
    perform 23 measurements on each piece. My goal is to find a machine
    learning method that detects defective pieces automatically. Both, my
    training set and my test set consist of 4000 pieces that were fine and
    30 defective ones.  I ran logistic regression on that data and found a
    training error of 0.1% and a test error of 0.3%, which is already pretty
    good, I think, no? But with kNN classification it is even better. For
    k = 7 I found a test error of 0.05% which is by far lower than all test
    errors I obtained with other values of k. I was really impressed. Now
    we have a method that predicts with 99.95% accuracy whether a piece is
    defective or not!

    Because you are taking this machine learning class now, I wanted to ask you
    for advice. Does it all make sense to you what I described? If not, do you
    have any suggestion to get even better results?

    Ã bientôt
    Jamie
```

## Applied
1. Take the `classification_data` in our notebook on
   \"flexibility and bias-variance-decomposition notebook\" and find with 10-fold
   cross-validation the optimal number ``k`` of neighbors of kNN
   classification, using the AUC measure. Hint: `MLJ` has the builtin function `auc`.
   Plot the validation AUC for ``k = 1, \ldots, 50``.
2. In this exercise you apply our \"recipe for supervised learning\" (see slides). The goal is to predict the miles a car can drive per gallon fuel (mpg) as a function of its horsepower. You can download a dataset with `using OpenML; cars = DataFrame(OpenML.load(455))`. In the cleaning step we will remove all rows that contain missing values (you can use the function `dropmissing`). We select the machine learning methods `PolynomialRegressor` and `KNNRegressor` and we take as measure the `rmse`. Make sure to go trough the steps 2, 5, 9 of the recipe. Plot the predictions of the best method you found.
"

# 1. Perform k-nearest neighbors regression on data generated with our `data_generator`
#    defined in the first cell of this notebook and find the optimal number k of neighbors
#    with k-fold cross validation.

# ╔═╡ 0651292e-3f4e-4263-8235-4caa563403ec
MLCourse.footer()

# ╔═╡ Cell order:
# ╠═736856ce-f490-11eb-3349-057c86edfe7e
# ╟─7eb6a060-f948-4d85-881a-4909e74c15bd
# ╠═7dd7e9a7-9245-4c64-af0c-8f7d2f62b2bf
# ╠═fa9b4b9e-4e97-4e5c-865d-ad3ef288e4cf
# ╟─751880ec-1a82-4142-b875-177d436bbc72
# ╟─da2bbb4d-d45d-4d2a-9c8a-21395bcbb851
# ╠═2e02ac0d-c4d0-47ba-be57-445adeb6ab8b
# ╟─91eecd2b-af18-4a63-9684-28950e604d1a
# ╟─b45b32b1-8a65-4823-b3bb-f0b7cc57604b
# ╠═ac3c7e84-6c47-4dc1-b862-de4cfb05dad9
# ╟─26362233-b006-423d-8fb5-7cd9150405b4
# ╟─5cff19c7-3426-4d91-b555-f059f5f41886
# ╠═df16313f-76be-4d8b-88e7-7cc618ff49d4
# ╟─46bffe01-1a7c-405b-b3ec-bfe7570b8a3c
# ╠═b3aed705-5072-4d9b-bb8f-865ac1561bf6
# ╠═6916273a-d116-4173-acaa-2bcac1d1753b
# ╠═c14c9e49-a993-456e-9115-97da86f8e498
# ╟─4b435518-2f12-4921-bb1f-fdd049ddfaed
# ╠═1e584a38-2fef-4877-87f6-92237d71c4b3
# ╟─a7c88b3f-92cb-4253-a889-c78683722c1d
# ╠═0ab1238e-8279-43d5-9803-848cb00156c2
# ╟─eaacf529-1727-4743-941b-360c53088b1d
# ╠═967a6e08-f0e6-45b7-988b-d0df237f3ddf
# ╟─535a11d9-833e-4237-949e-9a3e289c600b
# ╠═8ae790b3-3987-4f41-8e21-adbb71081eb9
# ╟─134df37f-b737-4c80-a5f9-1149aeec970c
# ╠═abb71af3-8aae-4806-9d5a-d144c15d22ef
# ╟─6599d2b4-68f4-4c22-8e40-bf3722597692
# ╟─29683c99-6a6a-4f65-bea2-d592895d887e
# ╠═f93f20db-4fed-481f-b085-ca744b68fa8f
# ╟─d72609a1-f93f-4ca6-8759-727662233e97
# ╠═59245b3e-ddfc-46c4-ba44-86ce191672ae
# ╟─f8b48f50-09cb-498a-89d1-9ac9b5722d0c
# ╟─47a2c14b-a02d-43d7-ac8a-4d95a0d91fa8
# ╟─ee89c448-1e69-4fd1-a4b8-7297a09f2685
# ╟─c693088f-7f80-4cdd-b9b5-65a50da732ac
# ╟─e1077092-f72a-42af-b6e0-a616f938cba8
# ╟─0651292e-3f4e-4263-8235-4caa563403ec

max_stars_repo_name: jbrea/MLCourse

max_issues_repo_name: davibarreira/MLCourse

max_forks_repo_name: davibarreira/MLCourse


text:
import Random

"""
    Simulation

A container to hold simulations.

All arguments are arrays of the same ength + a model from which the simulations came

# Arguments:
- `t::Array{Float64, 1}`:
- `φ::Array{Int, 1}`:
- `X::Array{Float64, 1}`:
- `n::Array{Int, 1}`:
- `model::Model`:
"""
struct Simulation
    # t, times at which SFM is observed (realisations of stopping times)
    # X, array of X(t) values
    # φ, array of φ(t) values
    # n, number of phase changes by t, n(t)
    # model, the model from which the sims came
    t::Array{Float64,1}
    φ::Array{Int,1}
    X::Array{Float64,1}
    n::Array{Int,1}
    model::Model
    function Simulation(t::Array{Float64,1},φ::Array{Int,1},
        X::Array{Float64,1},n::Array{Int,1},model::Model)
        !(length(t)==length(φ)==length(X)==length(n))&&throw(DomainError("lenght must all be the same"))
        new(t,φ,X,n,model)
    end
end

"""
Simulates a SFM defined by `Model` until the `StoppingTime` has occured,
given the `InitialCondition` on (φ(0),X(0)).

    simulate(
        model::Model,
        StoppingTime::Function,
        InitCondition::NamedTuple{(:φ, :X)},
    )

# Arguments
- `model`: A Model object
- `StoppingTime`: A function which takes the value of the process at the current
    time and at the time of the last jump of the phase process, as well as the
    `model` object.
    i.e. `StoppingTime(;model,SFM,SFM0,lwer,upr)` where `SFM` and `SFM0` are tuples with
    keys `(:t::Float64, :φ::Int, :X::Float64, :n::Int)` which are the value of
    the SFM at the current time, and time of the previous jump of the phase
    process, repsectively. The `StoppingTime` must return a
    `NamedTuple{(:Ind, :SFM)}` type where `:Ind` is a `:Bool` value stating
    whether the stopping time has occured or not and `:SFM` is a tuple in the
    same form as the input `SFM` but which contains the value of the `SFM` at
    the stopping time.
- `InitCondition`: `NamedTuple` with keys `(:φ, :X)`, `InitCondition.φ` is a
    vector of length `M` of initial states for the phase, and `InitCondition.X`
    is a vector of length `M` of initial states for the level. `M` is the number
    of simulations to be done.

# Output
- A `Simulation` object.
"""
function simulate(
    model::BoundedFluidQueue,
    StoppingTime::Function,
    InitCondition::NamedTuple{(:X,:φ)},
    rng::Random.AbstractRNG=Random.default_rng(),
)
    for n in 1:length(InitCondition.X)
        !(0.0<=InitCondition.X[n]<=model.b)&&throw(DomainError("inital X values must be within the model bounds"))
        !(InitCondition.φ[n]∈(1:n_phases(model)))&&throw(DomainError("inital φ values must be within the model phases"))
    end
    # the transition matrix of the jump chain
    d = LinearAlgebra.diag(model.T)
    P = (model.T - LinearAlgebra.diagm(0 => d)) ./ -d
    CumP = cumsum(P, dims = 2)

    idx₋ = findall(rates(model).<0.0)
    CumP_lwr = zeros(size(P))
    CumP_lwr[idx₋,:] = cumsum(model.P_lwr, dims = 2)
    CumP_upr = zeros(size(P))
    idx₊ = findall(rates(model).>0.0)
    CumP_upr[idx₊,:] = cumsum(model.P_upr, dims = 2)

    Λ = LinearAlgebra.diag(model.T)

    M = length(InitCondition.φ)
    tSims = zeros(Float64, M)
    φSims = zeros(Int, M)
    XSims = zeros(Float64, M)
    nSims = zeros(Int, M)

    for m = 1:M
        SFM0 = (t = 0.0, φ = InitCondition.φ[m], X = InitCondition.X[m], n = 0)
        while true
            S = log(rand(rng)) / Λ[SFM0.φ] # generate exp random variable
            X, φ, t = Update(model, SFM0, S, CumP, CumP_lwr, CumP_upr, rng)
            n = SFM0.n + 1
            SFM = (t = t, φ = φ, X = X, n = n)
            τ = StoppingTime(model, SFM, SFM0)
            if τ.Ind # if the stopping time occurs
                (tSims[m], φSims[m], XSims[m], nSims[m]) = τ.SFM
                break
            end
            SFM0 = SFM
        end
    end
    return Simulation(tSims, φSims, XSims, nSims, model)
end

function pdf(s::Simulation)
    throw(DomainError("cannot reliably construct pdf from simulation,
        either do it yourself, or construct the cdf, see cdf() function"))
end

"""
    cdf(s::Simulation)

Return a empirical cdf of the fluid queue given a simulation.

Output is a function of two variables (x,i) and gives the empirical distribution
    function ``P(X(τ)≤x,φ(τ)=i)``
"""
function cdf(s::Simulation)
    n_sims = length(s.t)
    function F(x::Real,i::Int)
        !(i∈phases(s.model)) && throw(DomainError("phase i must be in the support of the model"))
        i_idx = s.φ.==i
        Fxi = sum(s.X[i_idx].<=x)/n_sims
        return Fxi
    end
    return F
end

function cell_probs(s::Simulation,nodes)
    n_sims = length(s.t)
    function p(x::Real,i::Int)
        !(i∈phases(s.model)) && throw(DomainError("phase i must be in the support of the model"))
        _x_in_bounds = (x>nodes[1])&&(x<nodes[end])
        if _x_in_bounds
            idx = findfirst(x.<nodes)
            left_edge = nodes[idx-1]
            right_edge = nodes[idx]
            i_idx = s.φ.==i
            return (membership(s.model.S,i)>0.0) ? (sum(left_edge.<=s.X[i_idx].<right_edge)/n_sims) : (sum(left_edge.<s.X[i_idx].<=right_edge)/n_sims)
        else
            return 0.0
        end
    end
    return p
end

"""
Returns ``X(t+S) = min(max(X(t) + cᵢS,0),model.b)`` where ``model.b`` is the upper boundary of the model
on the process.

    UpdateX(
        model::Model,
        SFM0::NamedTuple,
        S::Real,
    )

# Arguments
- `model`: a Model object
- `SFM0::NamedTuple` containing at least the keys `:X` giving the value of
    ``X(t)`` at the current time, and `:φ` giving the value of
    ``φ(t)`` at the current time.
- `S::Real`: an elapsed amount of time to evaluate ``X`` at, i.e. ``X(t+S)``.
"""
function UpdateX(model::BoundedFluidQueue,SFM0::NamedTuple,S::Real)
    # given the last position of a SFM, SFM0, a time step of size s, find the
    # position of X at time t

    # X moves freely unless a boundary is hit
    X = min( max(SFM0.X+rates(model,SFM0.φ)*S, 0.0), model.b)
    return X
end

function Update(model::BoundedFluidQueue,SFM0::NamedTuple,S::Real,CumP::Matrix{Float64},
    CumP_lwr::Matrix{Float64},CumP_upr::Matrix{Float64},rng::Random.AbstractRNG=Random.default_rng())

    # evolve freely until a boundary is hit
    X =  UpdateX(model,SFM0,S)

    # if boundary is hit, change phase, update time to time when boundary is hit
    if (X==0.0)&&(SFM0.X!=0.0)#&&(rates(model,SFM0.φ)<0.0)
        φ = findfirst(rand(rng).<=CumP_lwr[SFM0.φ,:])
        t = SFM0.t + (X-SFM0.X)/rates(model,SFM0.φ)
    elseif (X==model.b)&&(SFM0.X!=model.b)#&&(rates(model,SFM0.φ)>0.0)
        φ = findfirst(rand(rng).<=CumP_upr[SFM0.φ,:])
        t = SFM0.t + (X-SFM0.X)/rates(model,SFM0.φ)
    else # update phase
        φ = findfirst(rand(rng).<=CumP[SFM0.φ, :])
        t = SFM0.t + S
    end
    return X, φ, t
end

"""
Constructs the `StoppingTime` ``1(t>T)``

    fixed_time( T::Real)

# Arguments
- `T`: a time at which to stop the process

# Output
- `fixed_timeFun`: a function with one methods
    - `fixed_timeFun(
        model::Model,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
    )`: a stopping time for a SFM.
"""
function fixed_time(T::Float64)
    # Defines a simple stopping time, 1(t>T).
    # SFM method
    function fixed_timeFun(
        model::Model,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
    )
        Ind = SFM.t > T
        if Ind
            s = T - SFM0.t
            X = UpdateX(model, SFM0, s)
            SFM = (t = T, φ = SFM0.φ, X = X, n = SFM0.n)
        end
        return (Ind = Ind, SFM = SFM)
    end
    return fixed_timeFun
end

"""
Constructs the `StoppingTime` ``1(N(t)>n)`` where ``N(t)`` is the number of
jumps of ``φ`` by time ``t``.

    n_jumps( N::Int)

# Arguments
- `N`: a desired number of jumps

# Output
- `n_jumpsFun`: a function with one methods
    - `n_jumpsFun(
        model::Model,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
    )`: a stopping time for a SFM.
"""
function n_jumps( N::Int)
    # Defines a simple stopping time, 1(n>N), where n is the number of jumps of φ.
    # SFM method
    function n_jumpsFun(
        model::Model,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
    )
        Ind = SFM.n >= N
        return (Ind = Ind, SFM = SFM)
    end
    return n_jumpsFun
end

"""
Constructs the `StoppingTime` which is the first exit of the process ``X(t)``
from the interval ``[u,v]``.

    first_exit_x( u::Real, v::Real)

# Arguments
- `u`: a lower boundary
- `v`: an upper boundary

# Output
- `first_exit_xFun`: a function with one methods
    - `first_exit_xFun(
        model::Model,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
    )`: a stopping time for a SFM.
"""
function first_exit_x( u::Real, v::Real)
    # SFM Method
    function first_exit_xFun(
        model::Model,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
    )
        Ind = ((SFM.X > v) || (SFM.X < u))
        if Ind
            if SFM.X > v
                X = v
            else
                X = u
            end
            s = (X - SFM0.X) / rates(model,SFM0.φ) # can't exit with c = 0
            t = SFM0.t + s
            SFM = (t = t, φ = SFM0.φ, X = X, n = SFM0.n)
        end
        return (Ind = Ind, SFM = SFM)
    end
    return first_exit_xFun
end


text:
# Energy calculation

export total_energy,
    kinetic_energy,
    temperature,
    potential_energy

"""
    total_energy(s, neighbors=nothing)

Compute the total energy of the system.
"""
total_energy(s, neighbors=nothing) = kinetic_energy(s) + potential_energy(s, neighbors)

kinetic_energy_noconvert(s) = sum(mass.(s.atoms) .* sum.(abs2, s.velocities)) / 2

"""
    kinetic_energy(s)

Compute the kinetic energy of the system.
"""
function kinetic_energy(s)
    ke = kinetic_energy_noconvert(s)
    # Convert energy to per mol if required
    if dimension(s.energy_units) == u"𝐋^2 * 𝐌 * 𝐍^-1 * 𝐓^-2"
        T = typeof(ustrip(ke))
        return uconvert(s.energy_units, ke * T(Unitful.Na))
    else
        return uconvert(s.energy_units, ke)
    end
end

"""
    temperature(system)

Calculate the temperature of a system from the kinetic energy of the atoms.
"""
function temperature(s)
    ke = kinetic_energy_noconvert(s)
    df = 3 * length(s) - 3
    T = typeof(ustrip(ke))
    k = unit(ke) == NoUnits ? one(T) : uconvert(u"K^-1" * unit(ke), T(Unitful.k))
    return 2 * ke / (df * k)
end

"""
    potential_energy(s, neighbors=nothing)

Compute the potential energy of the system.
If the interactions use neighbor lists, the neighbors should be computed
first and passed to the function.
"""
function potential_energy(s, neighbors=nothing)
    n_atoms = length(s)
    potential = zero(ustrip(s.box_size[1])) * s.energy_units

    for inter in values(s.general_inters)
        if inter.nl_only
            @inbounds for ni in 1:neighbors.n
                i, j, weight_14 = neighbors.list[ni]
                if weight_14
                    potential += potential_energy(inter, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size, true)
                else
                    potential += potential_energy(inter, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size)
                end
            end
        else
            for i in 1:n_atoms
                for j in (i + 1):n_atoms
                    potential += potential_energy(inter, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size)
                end
            end
        end
    end

    for inter_list in values(s.specific_inter_lists)
        potential += potential_energy(inter_list, s.coords, s.box_size)
    end

    return uconvert(s.energy_units, potential)
end

@views function potential_energy(inter_list::InteractionList2Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    (box_size,)))
end

@views function potential_energy(inter_list::InteractionList3Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], (box_size,)))
end

@views function potential_energy(inter_list::InteractionList4Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], coords[inter_list.ls], (box_size,)))
end


text:
function Spin1()

    Spin{Float64}(1, 1000, 100, 0)
    Spin{Float64}(1, Int32(1000), 100, 0)
    Spin{Float64}(1f0, 1000, 100, 0, Position(0, 1, 2))
    Spin{Float64}(1, Int32(1000), 100, 0, Position(0.0, -1.0, -1))
    Spin{Float64}(1f0, 1000f0, 100f0, 0f0)
    Spin{Float64}(Magnetization(1, 1, Int32(1)), 1, 1000, 100, 0)
    Spin(1, 1000, 100, 0)
    Spin(1, Int32(1000), 100, 0)
    Spin(1f0, 1000, 100, 0, Position(0, 1, 2))
    Spin(1, Int32(1000), 100, 0, Position(0.0, -1.0, -1))
    Spin(1f0, 1000f0, 100f0, 0f0)
    s = Spin(Magnetization(1, 1, Int32(1)), 1, 1000, 100, 0)

    show(devnull, "text/plain", s)
    return s.N == 1

end

function SpinMC1()

    SpinMC{Float64}(1, (0.8, 0.2), (1000, 100), (100, 10), (0, 0), (100, 25))
    SpinMC{Float64}(1, (0.8, 0.2), (1000f0, 100), (100, Int16(10)), (0, 0), (100, 25))
    SpinMC{Float64}((Magnetization(0, 0, 1.0), Magnetization(0, 0, 1)), 1, (1, 1),
           (100, 100), (10.0, 10), (0, 0.0), (1f0, 1.0))
    SpinMC{Float64}(((0, 0, 1.0), (0, 0, 1)), 1, (1, 1),
           (100, 100), (10.0, 10), (0, 0.0), (1f0, 1.0))
    SpinMC(1, (0.8, 0.2), (1000, 100), (100, 10), (0, 0), (100, 25))
    SpinMC(1, [0.8, 0.2], (1000f0, 100), (100, Int16(10)), (0, 0), (100, 25))
    SpinMC((Magnetization(0, 0, 1.0), Magnetization(0, 0, 1)), 1, (1, 1),
           [100, 100], (10.0, 10), [0, 0.0], (1f0, 1.0))
    s = SpinMC(((0, 0, 1.0), (0, 0, 1)), 1, (1, 1),
           [100, 100], (10.0, 10), [0, 0.0], (1f0, 1.0))

    show(devnull, s)
    show(devnull, "text/plain", s)
    return s.N == 2

end

function SpinMCError1()

    SpinMC((Magnetization(0,0,1),), 1, (1,), (1000,), (100,), (0,), ())

end

function SpinMCError2()

    SpinMC(1, 1, 1000, 100, 0, ())

end

function SpinMCError3()

    SpinMC(1, (1, 1), 1000, 100, 0, (2, 3))

end

function Position1()

    P1 = Position(1, 2, 3)
    P2 = Position(1.0, 2.0, 3.0)

    show(devnull, P1)
    show(devnull, "text/plain", P2)
    return P1 == P2

end

@testset "AbstractSpin" begin

    @testset "Spin" begin

        @test Spin1()

    end

    @testset "SpinMC" begin

        @test SpinMC1()
        @test_throws ErrorException SpinMCError1()
        @test_throws ErrorException SpinMCError2()
        @test_throws MethodError SpinMCError3()

    end

end

@testset "Position" begin

    @test Position1()

end


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

@doc raw"""ServicePort contains information on service's port.

    IoK8sApiCoreV1ServicePort(;
        name=nothing,
        nodePort=nothing,
        port=nothing,
        protocol=nothing,
        targetPort=nothing,
    )

    - name::String : The name of this port within the service. This must be a DNS_LABEL. All ports within a ServiceSpec must have unique names. When considering the endpoints for a Service, this must match the 'name' field in the EndpointPort. Optional if only one ServicePort is defined on this service.
    - nodePort::Int32 : The port on each node on which this service is exposed when type=NodePort or LoadBalancer. Usually assigned by the system. If specified, it will be allocated to the service if unused or else creation of the service will fail. Default is to auto-allocate a port if the ServiceType of this Service requires one. More info: https://kubernetes.io/docs/concepts/services-networking/service/#type-nodeport
    - port::Int32 : The port that will be exposed by this service.
    - protocol::String : The IP protocol for this port. Supports \"TCP\", \"UDP\", and \"SCTP\". Default is TCP.
    - targetPort::IoK8sApimachineryPkgUtilIntstrIntOrString : Number or name of the port to access on the pods targeted by the service. Number must be in the range 1 to 65535. Name must be an IANA_SVC_NAME. If this is a string, it will be looked up as a named port in the target Pod's container ports. If this is not specified, the value of the 'port' field is used (an identity map). This field is ignored for services with clusterIP=None, and should be omitted or set equal to the 'port' field. More info: https://kubernetes.io/docs/concepts/services-networking/service/#defining-a-service
"""
mutable struct IoK8sApiCoreV1ServicePort <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    nodePort::Any # spec type: Union{ Nothing, Int32 } # spec name: nodePort
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    targetPort::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: targetPort

    function IoK8sApiCoreV1ServicePort(;name=nothing, nodePort=nothing, port=nothing, protocol=nothing, targetPort=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("nodePort"), nodePort)
        setfield!(o, Symbol("nodePort"), nodePort)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("targetPort"), targetPort)
        setfield!(o, Symbol("targetPort"), targetPort)
        o
    end
end # type IoK8sApiCoreV1ServicePort

const _property_map_IoK8sApiCoreV1ServicePort = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("nodePort")=>Symbol("nodePort"), Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"), Symbol("targetPort")=>Symbol("targetPort"))
const _property_types_IoK8sApiCoreV1ServicePort = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("nodePort")=>"Int32", Symbol("port")=>"Int32", Symbol("protocol")=>"String", Symbol("targetPort")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiCoreV1ServicePort }) = collect(keys(_property_map_IoK8sApiCoreV1ServicePort))
Swagger.property_type(::Type{ IoK8sApiCoreV1ServicePort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ServicePort[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ServicePort }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ServicePort[property_name]

function check_required(o::IoK8sApiCoreV1ServicePort)
    (getproperty(o, Symbol("port")) === nothing) && (return false)
    true
end

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


text:
abstract type WiringT <: Function end

struct Wiring <: WiringT
    n_in::Int
    out::Int
    n_sensory::Int
    n_inter::Int
    n_command::Int
    n_motor::Int
    n_total::Int

    sensory_out::Int
    inter_out::Int
    rec_command_out::Int
    motor_in::Int

    sens_mask::Matrix{Float32}
    syn_mask::Matrix{Float32}

    sens_pol::Matrix{Float32}
    syn_pol::Matrix{Float32}

    # function Wiring(n_in,out,n_sensory,n_inter,n_command,n_motor,n_total,sensory_out,inter_out,rec_command_out,motor_in,sens_mask,syn_mask,sens_pol,syn_pol)
    #     new{typeof(sens_mask),typeof(syn_mask),typeof(sens_pol),typeof(syn_pol)}(
    #                 n_in,out,n_sensory,n_inter,n_command,n_motor,n_total,sensory_out,inter_out,rec_command_out,motor_in,sens_mask,syn_mask,sens_pol,syn_pol)
    # end
end

get_n_in(m::Wiring) = m.n_in
get_n_total(m::Wiring) = m.n_total

function Wiring(in::Int, out::Int;
                n_sensory=2, n_inter=5, n_command=0, n_motor=1,
                sensory_in=-1, rec_sensory=-1, sensory_out=-1,
                rec_inter=-1, inter_out=-1,                       # inter_in = sensory_out
                rec_command=-1, command_out=-1,                   # command_in = inter_out
                rec_motor=-1)                                     # motor_in = command_out, motor_out = out
  sensory_s = 1
  inter_s   = n_sensory + 1
  command_s = n_sensory + n_inter + 1
  motor_s   = n_sensory + n_inter + n_command + 1
  n_total   = n_sensory + n_inter + n_command + n_motor

  sens_mask = ones(Float32, in, n_total)
  syn_mask  = ones(Float32, n_total, n_total)

  sens_pol = ones(Float32, in, n_total)
  syn_pol  = ones(Float32, n_total, n_total)
  for i in eachindex(sens_pol)
    sens_pol[i] = [-1,1,1][rand(1:3)]
  end
  for i in eachindex(syn_pol)
    syn_pol[i] = [-1,1,1][rand(1:3)]
  end

  Wiring(in,out,n_sensory, n_inter,n_command,n_motor,n_total,sensory_out,inter_out,rec_command,command_out,sens_mask,syn_mask,sens_pol,syn_pol)
end

function add_synapses!(srcs,dsts,n_out,bitmask,polmask)
    n_out == 0 && return
    n_out = n_out == -1 ? length(dsts) : n_out
    n_out = min(n_out,length(dsts))
    for src in srcs
        for _ in n_out
            dst = dsts[rand(1:length(dsts))]
            bitmask[src,dst] = 1
            polmask[src,dst] = [-1,1,1][rand(1:3)]
        end
    end
end

function add_missing_synapses!(srcs, dsts, bitmask, polmask; n=1)
    for dst in dsts
        bm = @view bitmask[:,dst]
        pm = @view polmask[:,dst]

        if 1 ∉ filter(x -> x != dst, bm[srcs])
        #if 1 ∉ bm[srcs][bm[srcs] .!= dst]
            for i in 1:n
                random_src = srcs[rand(1:length(srcs))]
                bm[random_src] = 1
                pm[random_src] = [-1,1,1][rand(1:3)]
            end
        end
    end
end

function NCPWiring(in::Int, out::Int;
                n_sensory=2, n_inter=3, n_command=5, n_motor=1,
                sensory_in=-1, rec_sensory=0, sensory_inter=2, sensory_command=0, sensory_motor=0,
                inter_in=0, rec_inter=2, inter_command=2, inter_motor=0,                       # inter_in = sensory_out
                command_in=0, rec_command=1, command_motor=3,                   # command_in = inter_out
                motor_in=0, rec_motor=1)                                     # motor_in = command_out, motor_out = out
  sensory_s = 1
  inter_s   = n_sensory + 1
  command_s = n_sensory + n_inter + 1
  motor_s   = n_sensory + n_inter + n_command + 1
  n_total   = n_sensory + n_inter + n_command + n_motor

  out = min(n_total, out)

  sens_mask = zeros(Float32, in, n_total)
  syn_mask  = zeros(Float32, n_total, n_total)
  sens_pol = ones(Float32, in, n_total)
  syn_pol  = ones(Float32, n_total, n_total)

  input_range = range(1, length=in)
  sens_range = range(sensory_s, length=n_sensory)
  inter_range = range(inter_s, length=n_inter)
  command_range = range(command_s, length=n_command)
  motor_range = range(motor_s, length=n_motor)

  add_synapses!(input_range,    sens_range,    sensory_in, sens_mask, sens_pol)
  add_synapses!(input_range,    inter_range,    inter_in, sens_mask, sens_pol)
  add_synapses!(input_range,    command_range,    command_in, sens_mask, sens_pol)
  add_synapses!(input_range,    motor_range,    motor_in, sens_mask, sens_pol)

  add_synapses!(sens_range,    sens_range,    rec_sensory, syn_mask, syn_pol)
  add_synapses!(sens_range,    inter_range,   sensory_inter, syn_mask, syn_pol)
  add_synapses!(sens_range,    command_range, sensory_command, syn_mask, syn_pol)
  add_synapses!(sens_range,    motor_range,   sensory_motor, syn_mask, syn_pol)
  add_synapses!(inter_range,   inter_range,   rec_inter, syn_mask, syn_pol)
  add_synapses!(inter_range,   command_range, inter_command, syn_mask, syn_pol)
  add_synapses!(inter_range,   motor_range,   inter_motor, syn_mask, syn_pol)
  add_synapses!(command_range, command_range, rec_command, syn_mask, syn_pol)
  add_synapses!(command_range, motor_range,   command_motor, syn_mask, syn_pol)
  add_synapses!(motor_range,   motor_range,   rec_motor, syn_mask, syn_pol)

  add_missing_synapses!(sens_range, inter_range, syn_mask, syn_pol)
  add_missing_synapses!(inter_range, command_range, syn_mask, syn_pol)
  add_missing_synapses!(command_range, motor_range, syn_mask, syn_pol)

  Wiring(in,out,n_sensory, n_inter,n_command,n_motor,n_total,sensory_inter,inter_command,rec_command,command_motor,sens_mask,syn_mask,sens_pol,syn_pol)
end

Flux.trainable(m::Wiring) = ()


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

"""
    rotation_x(ϑ::Float64) :: Transformation

Encoding a rotation around the x-axis of an angle `ϑ` _**in radiant**_.

The positive sign is given by the right-hand rule, therefore clockwise
rotation for entering x-axis corresponds to a `ϑ>0` rotation angle.

See also: [`Transformation`](@ref)
"""
function rotation_x(ϑ::Float64)
    Transformation(
        [1.0    0.0     0.0     0.0 ;
         0.0    cos(ϑ)  -sin(ϑ) 0.0 ;
         0.0    sin(ϑ)  cos(ϑ)  0.0 ;
         0.0    0.0     0.0     1.0]
         ,
        [1.0    0.0     0.0     0.0 ;
         0.0    cos(ϑ)  sin(ϑ)  0.0 ;
         0.0    -sin(ϑ) cos(ϑ)  0.0 ;
         0.0    0.0     0.0     1.0]
    )
end

"""
    rotation_y(ϑ::Float64) :: Transformation

Encoding a rotation around the y-axis of an angle `ϑ` _**in radiant**_.

The positive sign is given by the right-hand rule, therefore clockwise
rotation for entering y-axis corresponds to a `ϑ>0` rotation angle.

See also: [`Transformation`](@ref)
"""
function rotation_y(ϑ::Float64)
    Transformation(
        [cos(ϑ)     0.0     sin(ϑ)  0.0 ;
         0.0        1.      0.0     0.0 ;
         -sin(ϑ)    0.0     cos(ϑ)  0.0 ;
         0.0        0.0     0.0     1.  ]
         ,
        [cos(ϑ)     0.0     -sin(ϑ) 0.0 ;
         0.0        1.      0.0     0.0 ;
         sin(ϑ)     0.0     cos(ϑ)  0.0 ;
         0.0        0.0     0.0     1.  ]
    )
end

"""
    rotation_z(ϑ::Float64) :: Transformation

Encoding a rotation around the z-axis of an angle `ϑ` _**in radiant**_.

The positive sign is given by the right-hand rule, therefore clockwise
rotation for entering z-axis corresponds to a `ϑ>0` rotation angle.

See also: [`Transformation`](@ref)
"""
function rotation_z(ϑ::Float64)
    Transformation(
        [cos(ϑ) -sin(ϑ) 0.0     0.0 ;
         sin(ϑ) cos(ϑ)  0.0     0.0 ;
         0.0    0.0     1.0     0.0 ;
         0.0    0.0     0.0     1.0]
         ,
        [cos(ϑ)     sin(ϑ)  0.0     0.0 ;
         -sin(ϑ)    cos(ϑ)  0.0     0.0 ;
         0.0        0.0     1.0     0.0 ;
         0.0        0.0     0.0     1.0]
    )
end

##########################################################################################92

function scaling(v::Vec)
    Transformation(
        [v.x    0.0     0.0     0.0 ;
         0.0    v.y     0.0     0.0 ;
         0.0    0.0     v.z     0.0 ;
         0.0    0.0     0.0     1.0]
         ,
        [1/v.x  0.0     0.0     0.0 ;
         0.0    1/v.y   0.0     0.0 ;
         0.0    0.0     1/v.z   0.0 ;
         0.0    0.0     0.0     1.0]
    )
end

scaling(x::Float64, y::Float64, z::Float64) = scaling(Vec(x,y,z))

"""
    scaling(v::Vec) :: Transformation
    scaling(x::Float64, y::Float64, z::Float64) = scaling(Vec(x,y,z))

Encoding a scaling of the 3 spatial coordinates according to the
vector `v` (negative values codify spatial reflections).

Each component of  `v` must be different from zero.

See also: [`Transformation`](@ref)
"""
scaling

##########################################################################################92

function translation(v::Vec)
   Transformation(
        [1.0    0.0     0.0     v.x ;
         0.0    1.0     0.0     v.y ;
         0.0    0.0     1.0     v.z ;
         0.0    0.0     0.0     1.0]
         ,
        [1.0    0.0     0.0     -v.x ;
         0.0    1.0     0.0     -v.y ;
         0.0    0.0     1.0     -v.z ;
         0.0    0.0     0.0     1.0]
    )
end

translation(x::Float64, y::Float64, z::Float64) = translation(Vec(x,y,z))

"""
    translation(v::Vec) :: Transformation
    translation(x::Float64, y::Float64, z::Float64) = translation(Vec(x,y,z))

Encoding a rigid translation of the 3 spatial coordinates according to the
vector `v`, which specifies the amount of shift to be applied along the three axes.

See also: [`Transformation`](@ref)
"""
translation

##########################################################################################92

"""
    inverse(T::Transformation) :: Transformation

Return the inverse affine transformation of `T`.
This method is very cheap to call.

See also: [`Transformation`](@ref)
"""
function inverse(T::Transformation)
    return Transformation(T.invM, T.M)
end

"""
    is_consistent(T::Transformation) :: Bool

Check the internal consistency of the  input transformation,
returning a bool variable indicating whether `T.M==T.invM`.
This method is useful when writing tests.

See also: [`Transformation`](@ref)
"""
function is_consistent(T::Transformation)
    p = T.M * T.invM
    I = SMatrix{4,4}( Diagonal(ones(4)) )
    return p ≈ I
end


text:
using LinearAlgebra: Eigen, eigen
using Tensorial: stress_invariants, deviatoric_stress_invariants, vol, dev

export principal_values,
    principal_axes, principal_invariants, main_invariants, hydrostatic, deviatoric

principal_values(x::Union{TensorStress,TensorStrain}) = _eigen(x).values

principal_axes(x::Union{TensorStress,TensorStrain}) = _eigen(x).vectors

principal_invariants(x::Union{TensorStress,TensorStrain}) = stress_invariants(x.data)

main_invariants(x::Union{TensorStress,TensorStrain}) = deviatoric_stress_invariants(x.data)

for T in (:TensorStress, :TensorStrain)
    @eval begin
        hydrostatic(x::$T) = $T(vol(SymmetricSecondOrderTensor{3}(float(x))))
        deviatoric(x::$T) = $T(dev(SymmetricSecondOrderTensor{3}(float(x))))
    end
end

function _eigen(x)
    unitless = @. x / oneunit(x)
    eg = eigen(unitless)
    return Eigen(eg.values * oneunit(x[1]), eg.vectors)
end


text:
using DataFrames
using CSV
using Plots
using Statistics
using GLM

using Latexify

df_e = CSV.read("zombiecar2_sweep_res5_endtrue.csv");
df_e.pc_median_step = df_e.median_step./df_e.steps
df = aggregate(df_e, [:n_agents, :discretize_m, :jump, :g_nodes, :g_edges], mean)
df.assymetry = df.left_tail_sum_mean./df.right_tail_sum_mean
df.pc_steps_mean = df.median_step_mean ./ df.steps_mean

df.e_sim_time = 2df.g_nodes.*log.(df.n_agents)./df.n_agents
df.log_n_x_n = log.(df.n_agents)./df.n_agents

df.log_agents = log.(df.n_agents)

r = DataFrame(discretization=Int[],jump=String[], R2=Float64[])
for d in [25., 50., 75.]
    for jj in [false, true]
        dat = df[(df.jump.==jj) .& (df.discretize_m.==d),:]
        println("COR log_n_x_n $(cor(dat.log_n_x_n, dat.steps_mean))")
        println("COR n_agents $(cor(dat.n_agents, dat.steps_mean))")
        println("COR e_sim_time $(cor(dat.e_sim_time, dat.steps_mean))")
        ols = lm(@formula(steps_mean~e_sim_time), dat)
        #display(ols) #coef(ols)
        println("R2=$(r2(ols))")
        push!(r,(d,(jj ? "Y" : "N"),round(cor(dat.log_n_x_n, dat.steps_mean),digits=6) ))
    end
end
println(latexify(r, env=:table))

r = DataFrame(jump=String[], R2=Float64[])
jj = false
for jj in [false, true]
    dat = df[(df.jump.==1) ,:]
    println("COR log_n_x_n $(cor(dat.log_n_x_n, dat.steps_mean))")
    println("COR n_agents $(cor(dat.n_agents, dat.steps_mean))")
    println("COR e_sim_time $(cor(dat.e_sim_time, dat.steps_mean))")
    ols = lm(@formula(steps_mean~e_sim_time), dat)
    #display(ols) #coef(ols)
    println("R2=$(r2(ols))")
    push!(r,((jj ? "Y" : "N"),round(cor(dat.e_sim_time, dat.steps_mean),digits=6) ))
end
println(latexify(r, env=:table))

function savefig2(plot,filename; twidth=raw"\columnwidth", replace_y_lab_pos=nothing)
    savefig(plot,filename)
    dat = read(filename, String)
    open(filename, "w") do f
        println(f,raw"\resizebox{"*twidth*"}{!}{%")
        if replace_y_lab_pos != nothing
            dat = replace(dat, "ylabel style={"=>"ylabel style={at={(axis description cs:$(replace_y_lab_pos))}")
        end

        println(f, dat)
        println(f,raw"}")
    end
end

function scatterplot(df, jumps=[true,false])
    p = Plots.scatter(xlabel = "Steps requred to deliver message to all participants",
        lab="", legend=:bottomright, size=(320,310))
    ylabel!(p, "Theoretical time: \$n \\ln k/k\$")
    shps=Dict(true=>:cross, false=>:star5)
    colrs=Dict(true=>:blue, false=>:red)
    for jj in jumps
        dd = df[(df.jump .== jj), :]
        p = Plots.scatter!(p,dd.steps_mean,dd.e_sim_time,
            markershape=shps[jj],
            markerstrokecolor=colrs[jj],
            lab="Jump-over = $(jj ? "Yes" : "No" )")
    end
    p
end

Plots.pgfplotsx()

fig = scatterplot(df[ df.n_agents .> 1, :])
savefig(fig, "../5d63f3b5c3791641ba87e19f/jump_yes_vs_no.pdf")

function scatterplot2(df)
    p = Plots.scatter(
        xlabel="Number of agents --- \$k\$", lab="", legend=:topright,
        size=(500,400))
    ylabel!(p, "Steps requred to deliver message to all participants")
    lcls=Dict(25=>:green, 50=>:blue, 75=>:red)
    shps=Dict(true=>:cross, false=>:star5)
    for discretize_m in [25, 50, 75]
        for jj in [true,false]
            dd = df[ (df.discretize_m .== discretize_m) .& (df.jump .== jj), :]
            p = Plots.scatter!(p,dd.n_agents, dd.steps_mean,
                markershape=shps[jj],
                markerstrokecolor=lcls[discretize_m],
                lab="Distance=$(discretize_m)m \\ \\ \\ \\ \\ Jump-over=$(jj ? "Yes" : "No" )")
        end
    end
    p
end

Plots.pgfplotsx()

fig = scatterplot2(df[ df.n_agents .> 40 , :])
savefig(fig, "../5d63f3b5c3791641ba87e19f/number_of_agents_vs_steps_min_agent70.pdf")

function scatterplot2logk(df)
    p = Plots.scatter(xlabel = "Steps requred to deliver message to all participants",
        ylabel=raw"$\ln k/k$", lab="", legend=:bottomright, size=(400,300))
    lcls=Dict(25=>:green, 50=>:blue, 75=>:red)
    shps=Dict(true=>:cross, false=>:star5)
    for discretize_m in [25, 50, 75]
        for jj in [true,false]
            dd = df[ (df.discretize_m .== discretize_m) .& (df.jump .== jj), :]
            p = Plots.scatter!(p,dd.steps_mean,log.(dd.n_agents)./dd.n_agents,
                markershape=shps[jj],
                markerstrokecolor=lcls[discretize_m],
                lab="Discr.=$(discretize_m), \\textit{Jump-over=$(jj ? "Yes" : "No" )}")
        end
    end
    p
end

Plots.pgfplotsx();
fig = scatterplot2logk(df[ df.n_agents .> 1, :]);
savefig(fig, "../5d63f3b5c3791641ba87e19f/scatterplotlogk_over_k.pdf")

function scatterplot_assy_median(df, jumps=[true,false])
    p = Plots.scatter(
        xlabel=raw"Number of agents --- $k$", lab="", legend=:bottomright,
        size=(400,400))
    ylabel!(p, "Percentage of simulation steps where half agents have received the message")
    lcls=Dict(25=>:green, 50=>:blue, 75=>:red)
    shps=Dict(true=>:cross, false=>:star5)
    for discretize_m in [25,50,75]
        for jj in jumps
            dd = df[ (df.discretize_m .== discretize_m) .& (df.jump .== jj), :]
            Plots.scatter!(p,dd.n_agents,dd.pc_steps_mean,
                markershape=shps[jj],
                markerstrokecolor=lcls[discretize_m],
                lab="Distance=$(discretize_m)m, Jump-over: $(jj ? "Yes" : "No" )")
        end
    end
    p
end

Plots.pgfplotsx()
fig = scatterplot_assy_median(df[ df.n_agents .> 1, :])
savefig(fig, "../5d63f3b5c3791641ba87e19f/scatterplot_assy_median.pdf")

function scatter_assymetry_left_right(df, jumps=[true, false])
    p = Plots.scatter(
        xlabel="Number of agents --- \$k\$", lab="", legend=:bottomright,
        size=(400,300))
    ylabel!(p, raw"Assymetry = $ left\_tail / upper\_right\_tail $")
    lcls=Dict(25=>:green, 50=>:blue, 75=>:red)
    shps=Dict(true=>:cross, false=>:star5)
    for discretize_m in [25,50,75]
        for jj in jumps
            dd = df[ (df.discretize_m .== discretize_m) .& (df.jump .== jj), :]
            Plots.scatter!(p,dd.n_agents,dd.assymetry,
                markershape=shps[jj],
                markerstrokecolor=lcls[discretize_m],
                lab="Distance=$(discretize_m)m, Jump-over: $(jj ? "Yes" : "No" )")
        end
    end
    p
end

fig = scatter_assymetry_left_right(df[ df.n_agents .> 1, :])
savefig(fig, "../5d63f3b5c3791641ba87e19f/scatter_assymetry_left_right.pdf")

function scatter_assymetry_both_median_left_right(df)
    p = Plots.scatter(
        xlabel="Percentage of simulation steps where half agents have received the message", lab="", legend=:bottomright)
    ylabel!(p, raw"Assymetry = $ left\_tail / upper\_right\_tail $")
    shps=Dict(25=>:circle, 50=>:rect, 75=>:diamond)

    for discretize_m in [25,50,75]
        jj = true
        dd = df[ (df.discretize_m .== discretize_m) .& (df.jump .== jj), :]
        Plots.scatter!(p, dd.pc_steps_mean[1:1],dd.assymetry[1:1],
            lab="Distance=$(discretize_m)m, Jump-over: $(jj ? "Yes" : "No" )",
            markershape=shps[discretize_m], color=:black )
    end

    for discretize_m in [25,50,75]
        for jj in [true]
            dd = df[ (df.discretize_m .== discretize_m) .& (df.jump .== jj), :]

            Plots.scatter!(p,dd.pc_steps_mean,dd.assymetry,
                markershape=shps[discretize_m],
                zcolor=log.(dd.n_agents), palette=:balance,
                markerstrokecolor=nothing,
                lab="")

        end
    end
    p
end

fig = scatter_assymetry_both_median_left_right(df[ df.n_agents .> 1, :])
savefig2(fig, "../5d63f3b5c3791641ba87e19f/scatter_assymetry_both_median_left_right_pdf.pdf", replace_y_lab_pos="-0.028,0.5")

function plot_sim_res(df::DataFrame,
        kw_list::Vector{NamedTuple{(:n_agents, :discretize_m, :jump),Tuple{Int64,Float64,Int64}}};
        filename::Union{AbstractString,Nothing}=nothing)

    ds = [get_slice(df; kw...) for kw in kw_list]
    for i in 1:length(kw_list)
        @assert nrow(ds[i]) > 0 "ERRROR no rows found for $(string(kw_list[i]))"
    end

    labs = ["n=$(kw.n_agents), j=$(kw.jump==1 ? "Y" : "N"), d=$(round(Int, kw.discretize_m))m" for kw in kw_list]
    ordering = sortperm(ds; lt = (x,y) -> x.step[end] < y.step[end] )
    ds = ds[ordering]
    labs = labs[ordering]

    ncolor = max(length(kw_list),2)
    cols = reshape( range(colorant"blue", stop=colorant"red",length=ncolor), 1, ncolor );
    p = Plots.plot(;xlabel = "Simulation steps", ylabel="Share of agent who received the message", ylim=(0.0,1.1), legend=:bottomright)
    for i in 1:length(kw_list)
        plot_res!(p, ds[i], cols[i], labs[i])
    end
    filename!=nothing && savefig(p, filename)
    p
end


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/misc_utils/convert_koi_cat_csv_to_jld.jl
## (c) 2018 Danley C. Hsu
# Converts KOI catalog data file from csv to jld format

using ExoplanetsSysSim
using HDF5, JLD

include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2018", "christiansen_func.jl"))

koi_catalog_file_in = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q16_koi_cand.csv")
koi_catalog_file_out = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q16_koi_cand.jld")

koi_catalog, koi_usable = read_koi_catalog(ascii(koi_catalog_file_in))

save(koi_catalog_file_out,"koi_catalog", koi_catalog, "koi_catalog_usable", koi_usable)


text:
module RocketTakeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: take()" begin

    println("Testing: operator take()")

    run_proxyshowcheck("Take", take(1))

    run_testset([
        (
            source = from(1:5) |> take(3),
            values = @ts([ 1:3, c ])
        ),
        (
            source = from(1:5) |> take(0),
            values = @ts(c)
        ),
        (
            source = from(1:5) |> async(0) |> take(3),
            values = @ts([ 1 ] ~ [ 2 ] ~ [ 3, c ])
        ),
        (
            source = timer(100, 30) |> take(3),
            values = @ts(100 ~ [ 0 ] ~ 30 ~ [ 1 ] ~ 30 ~ [ 2, c ])
        ),
        (
            source = completed() |> take(10),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> take(10),
            values = @ts(e("e"))
        ),
        (
            source = never() |> take(10),
            values = @ts()
        )
    ])

    @testset "Infinite reaction" begin
        source = Subject(Int)
        events = []

        subscription = subscribe!(source |> take(3), lambda(
            on_next     = (state) -> begin
                push!(events, state)
                next!(source, state + 1)
            end,
            on_complete = () -> push!(events, "c")
        ))

        next!(source, 1)

        @test events == [ 1, 2, 3, "c" ]
    end

end

end