text:
using Test
import GeoInterface
import ArchGDAL as AG
import GeoFormatTypes as GFT

@testset "test_geometry.jl" begin
    @testset "Incomplete GeoInterface geometries" begin
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbCircularString),
        )
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbCompoundCurve),
        )
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbCurvePolygon),
        )
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbMultiSurface),
        )
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbPolyhedralSurface),
        )
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbTIN),
        )
        @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(
            AG.creategeom(AG.wkbTriangle),
        )
    end

    @testset "Create a Point" begin
        # Method 1
        AG.createpoint(100, 70) do point
            @test GeoInterface.geotype(point) == :Point
            @test isapprox(
                GeoInterface.coordinates(point),
                [100, 70],
                atol = 1e-6,
            )
            @test AG.geomdim(point) == 0
            @test AG.getcoorddim(point) == 2
            AG.setcoorddim!(point, 3)
            @test AG.getcoorddim(point) == 3
            @test AG.isvalid(point) == true
            @test AG.issimple(point) == true
            @test AG.isring(point) == false
            @test AG.getz(point, 0) == 0

            @test sprint(print, AG.envelope(point)) ==
                  "GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)"
            @test sprint(print, AG.envelope3d(point)) ==
                  "GDAL.OGREnvelope3D(100.0, 100.0, 70.0, 70.0, 0.0, 0.0)"
            @test AG.toISOWKB(point, AG.wkbNDR) == UInt8[
                0x01,
                0xe9,
                0x03,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x59,
                0x40,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x80,
                0x51,
                0x40,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
            ]
            @test AG.toISOWKB(point, AG.wkbXDR) == UInt8[
                0x00,
                0x00,
                0x00,
                0x03,
                0xe9,
                0x40,
                0x59,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x40,
                0x51,
                0x80,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
                0x00,
            ]
            @test AG.toKML(point, "relativeToGround") ==
                  "<Point><altitudeMode>relativeToGround</altitudeMode>" *
                  "<coordinates>100,70,0</coordinates></Point>"
            @test AG.toKML(point, "clampToGround") ==
                  "<Point><altitudeMode>clampToGround</altitudeMode>" *
                  "<coordinates>100,70,0</coordinates></Point>"
            @test AG.toKML(point) ==
                  "<Point><coordinates>100,70,0</coordinates></Point>"
            @test AG.toJSON(point) ==
                  "{ \"type\": \"Point\", \"coordinates\": " *
                  "[ 100.0, 70.0, 0.0 ] }"
            @test startswith(
                AG.toJSON(point, SIGNIFICANT_FIGURES = 1),
                "{ \"type\": \"Point\", \"coordinates\": [",
            )
            @test startswith(
                AG.toJSON(point, ["SIGNIFICANT_FIGURES=1"]),
                "{ \"type\": \"Point\", \"coordinates\": [",
            )
            AG.createpoint(100, 70, 0) do point2
                @test isapprox(
                    GeoInterface.coordinates(point2),
                    [100, 70, 0],
                    atol = 1e-6,
                )
                @test AG.equals(point, point2) == true
            end
            AG.createpoint((100, 70, 0)) do point3
                @test AG.equals(point, point3) == true
            end
            AG.createpoint([100, 70, 0]) do point4
                @test AG.equals(point, point4) == true
            end
            point5 = AG.createpoint([100, 70, 0])
            @test AG.equals(point, point5) == true
            AG.flattento2d!(point)
            @test AG.getcoorddim(point) == 2
            @test AG.getnonlineargeomflag() == true
            AG.setnonlineargeomflag!(false)
            @test AG.getnonlineargeomflag() == false
            AG.setnonlineargeomflag!(true)
            @test AG.getnonlineargeomflag() == true
            AG.closerings!(point)
            @test AG.toJSON(point) ==
                  "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0 ] }"
        end

        # Method 2
        point = AG.createpoint(100, 70)
        @test AG.geomdim(point) == 0
        @test AG.getcoorddim(point) == 2
        AG.setcoorddim!(point, 3)
        @test AG.getcoorddim(point) == 3
        @test AG.isvalid(point) == true
        @test AG.issimple(point) == true
        @test AG.isring(point) == false
        @test AG.getz(point, 0) == 0
        @test typeof(point) == AG.IGeometry{AG.wkbPoint}
        @test sprint(print, AG.envelope(point)) ==
              "GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)"
        @test sprint(print, AG.envelope3d(point)) ==
              "GDAL.OGREnvelope3D(100.0, 100.0, 70.0, 70.0, 0.0, 0.0)"
        @test AG.toISOWKB(point, AG.wkbNDR) == UInt8[
            0x01,
            0xe9,
            0x03,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x59,
            0x40,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x80,
            0x51,
            0x40,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
        ]
        @test AG.toISOWKB(point, AG.wkbXDR) == UInt8[
            0x00,
            0x00,
            0x00,
            0x03,
            0xe9,
            0x40,
            0x59,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x40,
            0x51,
            0x80,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
            0x00,
        ]
        @test AG.toKML(point, "relativeToGround") ==
              "<Point><altitudeMode>relativeToGround</altitudeMode>" *
              "<coordinates>100,70,0</coordinates></Point>"
        @test AG.toKML(point, "clampToGround") ==
              "<Point><altitudeMode>clampToGround</altitudeMode>" *
              "<coordinates>100,70,0</coordinates></Point>"
        @test AG.toKML(point) ==
              "<Point><coordinates>100,70,0</coordinates></Point>"
        @test AG.toJSON(point) ==
              "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0, 0.0 ] }"
        @test AG.equals(point, AG.createpoint(100, 70, 0)) == true
        @test AG.equals(point, AG.createpoint((100, 70, 0))) == true
        AG.flattento2d!(point)
        @test AG.getcoorddim(point) == 2
        @test AG.getnonlineargeomflag() == true
        AG.setnonlineargeomflag!(false)
        @test AG.getnonlineargeomflag() == false
        AG.setnonlineargeomflag!(true)
        @test AG.getnonlineargeomflag() == true
        AG.closerings!(point)
        @test AG.toJSON(point) ==
              "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0 ] }"
    end

    @testset "Testing construction of complex geometries" begin
        @test AG.toWKT(AG.createlinestring([1.0, 2.0, 3.0], [4.0, 5.0, 6.0])) ==
              "LINESTRING (1 4,2 5,3 6)"
        AG.createlinestring([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) do geom
            @test GeoInterface.geotype(geom) == :LineString
            @test isapprox(
                GeoInterface.coordinates(geom),
                [[1, 4], [2, 5], [3, 6]],
                atol = 1e-6,
            )
            @test AG.toWKT(geom) == "LINESTRING (1 4,2 5,3 6)"
            AG.closerings!(geom)
            @test AG.toWKT(geom) == "LINESTRING (1 4,2 5,3 6)"
            AG.setpoint!(geom, 1, 10, 10)
            @test AG.toWKT(geom) == "LINESTRING (1 4,10 10,3 6)"
            @test GFT.val(convert(GFT.WellKnownText, geom)) == AG.toWKT(geom)
            @test typeof(geom) == AG.Geometry{AG.wkbLineString}
        end
        AG.createlinestring(
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        ) do geom
            @test AG.toWKT(geom) == "LINESTRING (1 4 7,2 5 8,3 6 9)"
            AG.setpoint!(geom, 1, 10, 10, 10)
            @test AG.toWKT(geom) == "LINESTRING (1 4 7,10 10 10,3 6 9)"
            AG.addpoint!(geom, 11, 11, 11)
            @test AG.toWKT(geom) == "LINESTRING (1 4 7,10 10 10,3 6 9,11 11 11)"
        end

        @test AG.toWKT(AG.createlinearring([1.0, 2.0, 3.0], [4.0, 5.0, 6.0])) ==
              "LINEARRING (1 4,2 5,3 6)"
        AG.createlinearring([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) do geom
            @test GeoInterface.geotype(geom) == :LineString
            @test isapprox(
                GeoInterface.coordinates(geom),
                [[1, 4], [2, 5], [3, 6]],
                atol = 1e-6,
            )
            @test AG.toWKT(geom) == "LINEARRING (1 4,2 5,3 6)"
            AG.setpointcount!(geom, 5)
            @test AG.toWKT(geom) == "LINEARRING (1 4,2 5,3 6,0 0,0 0)"
            AG.empty!(geom)
            @test AG.toWKT(geom) == "LINEARRING EMPTY"
            @test typeof(geom) == AG.Geometry{AG.wkbLineString} # this seems odd
        end
        AG.createlinearring(
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        ) do geom
            @test AG.toWKT(geom) == "LINEARRING (1 4 7,2 5 8,3 6 9)"
            AG.closerings!(geom)
            @test AG.toWKT(geom) == "LINEARRING (1 4 7,2 5 8,3 6 9,1 4 7)"
        end

        @test AG.toWKT(AG.createpolygon([1.0, 2.0, 3.0], [4.0, 5.0, 6.0])) ==
              "POLYGON ((1 4,2 5,3 6))"
        AG.createpolygon([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) do geom
            @test GeoInterface.geotype(geom) == :Polygon
            @test isapprox(
                GeoInterface.coordinates(geom),
                [[[1, 4], [2, 5], [3, 6]]],
                atol = 1e-6,
            )
            @test AG.toWKT(geom) == "POLYGON ((1 4,2 5,3 6))"
            @test typeof(geom) == AG.Geometry{AG.wkbPolygon}
        end
        AG.createpolygon(
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        ) do geom
            @test AG.toWKT(geom) == "POLYGON ((1 4 7,2 5 8,3 6 9))"
            AG.closerings!(geom)
            @test AG.toWKT(geom) == "POLYGON ((1 4 7,2 5 8,3 6 9,1 4 7))"
        end

        @test AG.toWKT(AG.createmultipoint([1.0, 2.0, 3.0], [4.0, 5.0, 6.0])) ==
              "MULTIPOINT (1 4,2 5,3 6)"
        AG.createmultipoint([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) do geom
            @test GeoInterface.geotype(geom) == :MultiPoint
            @test isapprox(
                GeoInterface.coordinates(geom),
                [[1, 4], [2, 5], [3, 6]],
                atol = 1e-6,
            )
            @test AG.toWKT(geom) == "MULTIPOINT (1 4,2 5,3 6)"
            @test typeof(geom) == AG.Geometry{AG.wkbMultiPoint}
        end
        AG.createmultipoint(
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        ) do geom
            @test AG.toWKT(geom) == "MULTIPOINT (1 4 7,2 5 8,3 6 9)"
        end

        @test AG.toWKT(
            AG.createmultipolygon(
                Vector{Vector{Tuple{Cdouble,Cdouble}}}[
                    Vector{Tuple{Cdouble,Cdouble}}[
                        [(0, 0), (0, 4), (4, 4), (4, 0)],
                        [(1, 1), (1, 3), (3, 3), (3, 1)],
                    ],
                    Vector{Tuple{Cdouble,Cdouble}}[
                        [(10, 0), (10, 4), (14, 4), (14, 0)],
                        [(11, 1), (11, 3), (13, 3), (13, 1)],
                    ],
                ],
            ),
        ) ==
              "MULTIPOLYGON (" *
              "((0 0,0 4,4 4,4 0),(1 1,1 3,3 3,3 1))," *
              "((10 0,10 4,14 4,14 0),(11 1,11 3,13 3,13 1)))"
        AG.createmultipolygon(
            Vector{Vector{Tuple{Cdouble,Cdouble}}}[
                Vector{Tuple{Cdouble,Cdouble}}[
                    [(0, 0), (0, 4), (4, 4), (4, 0)],
                    [(1, 1), (1, 3), (3, 3), (3, 1)],
                ],
                Vector{Tuple{Cdouble,Cdouble}}[
                    [(10, 0), (10, 4), (14, 4), (14, 0)],
                    [(11, 1), (11, 3), (13, 3), (13, 1)],
                ],
            ],
        ) do geom
            @test GeoInterface.geotype(geom) == :MultiPolygon
            @test isapprox(
                GeoInterface.coordinates(geom),
                [
                    [
                        [[0, 0], [0, 4], [4, 4], [4, 0]],
                        [[1, 1], [1, 3], [3, 3], [3, 1]],
                    ],
                    [
                        [[10, 0], [10, 4], [14, 4], [14, 0]],
                        [[11, 1], [11, 3], [13, 3], [13, 1]],
                    ],
                ],
                atol = 1e-6,
            )
            @test AG.toWKT(geom) ==
                  "MULTIPOLYGON (" *
                  "((0 0,0 4,4 4,4 0),(1 1,1 3,3 3,3 1))," *
                  "((10 0,10 4,14 4,14 0),(11 1,11 3,13 3,13 1)))"
            @test typeof(geom) == AG.Geometry{AG.wkbMultiPolygon}
        end

        AG.fromWKT(
            "CURVEPOLYGON (" *
            "CIRCULARSTRING (-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0)," *
            "(-1 0,0 0.5,1 0,0 1,-1 0))",
        ) do geom
            @test typeof(geom) == AG.Geometry{AG.wkbCurvePolygon}
            @test AG.toWKT(AG.curvegeom(AG.lineargeom(geom, 0.5))) ==
                  "CURVEPOLYGON (" *
                  "CIRCULARSTRING (-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0)," *
                  "(-1 0,0.0 0.5,1 0,0 1,-1 0))"
            AG.lineargeom(geom, 0.5) do lgeom
                @test typeof(lgeom) == AG.Geometry{AG.wkbPolygon}
                AG.curvegeom(lgeom) do clgeom
                    @test AG.toWKT(clgeom) ==
                          "CURVEPOLYGON (" *
                          "CIRCULARSTRING (-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0)," *
                          "(-1 0,0.0 0.5,1 0,0 1,-1 0))"
                    @test typeof(clgeom) == AG.Geometry{AG.wkbCurvePolygon}
                end
                @test AG.ngeom(
                    AG.polygonize(AG.forceto(lgeom, AG.wkbMultiLineString)),
                ) == 2
                AG.forceto(lgeom, AG.wkbMultiLineString) do mlsgeom
                    @test typeof(mlsgeom) == AG.Geometry{AG.wkbMultiLineString}
                    AG.polygonize(mlsgeom) do plgeom
                        @test AG.ngeom(plgeom) == 2
                        @test typeof(plgeom) ==
                              AG.Geometry{AG.wkbGeometryCollection}
                    end
                end
            end

            @test startswith(
                AG.toWKT(
                    AG.curvegeom(
                        AG.lineargeom(geom, 0.5, ADD_INTERMEDIATE_POINT = "NO"),
                    ),
                ),
                "CURVEPOLYGON (CIRCULARSTRING (",
            )
            AG.lineargeom(geom, 0.5, ADD_INTERMEDIATE_POINT = "NO") do lgeom
                AG.curvegeom(lgeom) do clgeom
                    @test startswith(
                        AG.toWKT(clgeom),
                        "CURVEPOLYGON (CIRCULARSTRING (",
                    )
                end
                @test AG.ngeom(
                    AG.polygonize(AG.forceto(lgeom, AG.wkbMultiLineString)),
                ) == 2
                AG.forceto(lgeom, AG.wkbMultiLineString) do mlsgeom
                    AG.polygonize(mlsgeom) do plgeom
                        @test AG.ngeom(plgeom) == 2
                    end
                end
            end

            @test startswith(
                AG.toWKT(
                    AG.curvegeom(
                        AG.lineargeom(geom, ["ADD_INTERMEDIATE_POINT=NO"], 0.5),
                    ),
                ),
                "CURVEPOLYGON (CIRCULARSTRING (",
            )
            AG.lineargeom(geom, ["ADD_INTERMEDIATE_POINT=NO"], 0.5) do lgeom
                AG.curvegeom(lgeom) do clgeom
                    @test startswith(
                        AG.toWKT(clgeom),
                        "CURVEPOLYGON (CIRCULARSTRING (",
                    )
                end
                @test AG.ngeom(
                    AG.polygonize(AG.forceto(lgeom, AG.wkbMultiLineString)),
                ) == 2
                AG.forceto(lgeom, AG.wkbMultiLineString) do mlsgeom
                    AG.polygonize(mlsgeom) do plgeom
                        @test AG.ngeom(plgeom) == 2
                    end
                end
            end
        end
    end

    @testset "Testing remaining methods for geometries" begin
        geom1 = AG.createmultipolygon(
            Vector{Vector{Tuple{Cdouble,Cdouble}}}[
                Vector{Tuple{Cdouble,Cdouble}}[
                    [(0, 0), (0, 4), (4, 4), (4, 0)],
                    [(1, 1), (1, 3), (3, 3), (3, 1)],
                ],
                Vector{Tuple{Cdouble,Cdouble}}[
                    [(10, 0), (10, 4), (14, 4), (14, 0)],
                    [(11, 1), (11, 3), (13, 3), (13, 1)],
                ],
            ],
        )
        geom2 = AG.createmultipoint(
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        )

        AG.closerings!(geom1)
        @test AG.disjoint(geom1, geom2) == false
        @test AG.touches(geom1, geom2) == true
        @test AG.crosses(geom1, geom2) == false
        @test AG.overlaps(geom1, geom2) == false

        @test AG.toWKT(AG.boundary(geom2)) == "GEOMETRYCOLLECTION EMPTY"
        AG.boundary(geom2) do result
            @test AG.toWKT(result) == "GEOMETRYCOLLECTION EMPTY"
        end

        @test AG.toWKT(AG.union(geom1, geom2)) ==
              "GEOMETRYCOLLECTION (" *
              "POLYGON (" *
              "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
              "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
              "POLYGON (" *
              "(10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
              "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8))," *
              "POINT (2 5 8),POINT (3 6 9))"
        AG.union(geom1, geom2) do result
            @test AG.toWKT(result) ==
                  "GEOMETRYCOLLECTION (" *
                  "POLYGON (" *
                  "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
                  "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
                  "POLYGON (" *
                  "(10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
                  "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8))," *
                  "POINT (2 5 8),POINT (3 6 9))"
            @test AG.hascurvegeom(result, true) == false
            @test AG.hascurvegeom(result, false) == false
        end

        @test AG.toWKT(AG.difference(geom1, geom2)) ==
              "MULTIPOLYGON (" *
              "((0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
              "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
              "((10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
              "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8)))"
        AG.difference(geom1, geom2) do result
            @test AG.toWKT(result) ==
                  "MULTIPOLYGON (" *
                  "((0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
                  "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
                  "((10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
                  "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8)))"
            AG.segmentize!(result, 20)
            @test AG.toWKT(result) ==
                  "MULTIPOLYGON (" *
                  "((0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
                  "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
                  "((10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
                  "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8)))"
            AG.segmentize!(result, 2)
            @test AG.toWKT(result) ==
                  "MULTIPOLYGON (" *
                  "(" *
                  "(" *
                  "0 4 8," *
                  "2 4 8," *
                  "4 4 8," *
                  "4 2 8," *
                  "4 0 8," *
                  "2 0 8," *
                  "0 0 8," *
                  "0 2 8," *
                  "0 4 8)," *
                  "(" *
                  "3 1 8," *
                  "3 3 8," *
                  "1 3 8," *
                  "1 1 8," *
                  "3 1 8))," *
                  "(" *
                  "(" *
                  "10 4 8," *
                  "12 4 8," *
                  "14 4 8," *
                  "14 2 8," *
                  "14 0 8," *
                  "12 0 8," *
                  "10 0 8," *
                  "10 2 8," *
                  "10 4 8)," *
                  "(" *
                  "13 1 8," *
                  "13 3 8," *
                  "11 3 8," *
                  "11 1 8," *
                  "13 1 8)))"
            @test typeof(result) == AG.Geometry{AG.wkbMultiPolygon25D}
        end

        @test AG.toWKT(AG.symdifference(geom1, geom2)) ==
              "GEOMETRYCOLLECTION (" *
              "POLYGON (" *
              "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
              "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
              "POLYGON (" *
              "(10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
              "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8))," *
              "POINT (2 5 8),POINT (3 6 9))"
        AG.symdifference(geom1, geom2) do result
            @test GeoInterface.geotype(result) == :GeometryCollection
            @test AG.toWKT(result) ==
                  "GEOMETRYCOLLECTION (" *
                  "POLYGON (" *
                  "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
                  "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
                  "POLYGON (" *
                  "(10 4 8,14 4 8,14 0 8,10 0 8,10 4 8)," *
                  "(13 1 8,13 3 8,11 3 8,11 1 8,13 1 8))," *
                  "POINT (2 5 8)," *
                  "POINT (3 6 9))"
            AG.removegeom!(result, 1)
            @test AG.toWKT(result) ==
                  "GEOMETRYCOLLECTION (" *
                  "POLYGON (" *
                  "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
                  "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
                  "POINT (2 5 8)," *
                  "POINT (3 6 9))"
            AG.removeallgeoms!(result)
            @test AG.toWKT(result) == "GEOMETRYCOLLECTION EMPTY"
            @test typeof(result) == AG.Geometry{AG.wkbGeometryCollection25D}
        end

        geom3 = AG.fromWKT(
            "GEOMETRYCOLLECTION (" *
            "POINT (2 5 8)," *
            "POLYGON (" *
            "(0 0 8,0 4 8,4 4 8,4 0 8,0 0 8)," *
            "(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8))," *
            "POLYGON (" *
            "(10 0 8,10 4 8,14 4 8,14 0 8,10 0 8)," *
            "(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8))," *
            "POINT EMPTY)",
        )
        AG.clone(geom3) do geom4
            @test sprint(print, AG.clone(geom3)) ==
                  "Geometry: GEOMETRYCOLLECTION (" *
                  "POINT (2 5 8)," *
                  "POLYGON ((0 0 8," *
                  " ... MPTY)"
            @test sprint(print, AG.clone(geom4)) ==
                  "Geometry: GEOMETRYCOLLECTION (" *
                  "POINT (2 5 8)," *
                  "POLYGON ((0 0 8," *
                  " ... MPTY)"
            @test typeof(geom4) == AG.Geometry{AG.wkbGeometryCollection25D}
        end
        AG.clone(AG.getgeom(geom3, 3)) do geom4
            @test sprint(print, geom4) == "Geometry: POINT EMPTY"
        end

        @test AG.toISOWKT(geom3) ==
              "GEOMETRYCOLLECTION Z (" *
              "POINT Z (2 5 8)," *
              "POLYGON Z (" *
              "(0 0 8,0 4 8,4 4 8,4 0 8,0 0 8)," *
              "(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8))," *
              "POLYGON Z (" *
              "(10 0 8,10 4 8,14 4 8,14 0 8,10 0 8)," *
              "(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8))," *
              "POINT Z EMPTY)"
        # the JSON driver in GDAL 3.0 does not handle null geometries well yet
        AG.removegeom!(geom3, AG.ngeom(geom3) - 1)
        @test AG.toJSON(geom3) ==
              """{ "type": "GeometryCollection", "geometries": [ """ *
              """{ "type": "Point", "coordinates": [ 2.0, 5.0, 8.0 ] }, """ *
              """{ "type": "Polygon", "coordinates": [ """ *
              "[ " *
              "[ 0.0, 0.0, 8.0 ], " *
              "[ 0.0, 4.0, 8.0 ], " *
              "[ 4.0, 4.0, 8.0 ], " *
              "[ 4.0, 0.0, 8.0 ], " *
              "[ 0.0, 0.0, 8.0 ] ], " *
              "[ " *
              "[ 1.0, 1.0, 8.0 ], " *
              "[ 3.0, 1.0, 8.0 ], " *
              "[ 3.0, 3.0, 8.0 ], " *
              "[ 1.0, 3.0, 8.0 ], " *
              "[ 1.0, 1.0, 8.0 ] ] ] }, " *
              """{ "type": "Polygon", "coordinates": [ """ *
              "[ " *
              "[ 10.0, 0.0, 8.0 ], " *
              "[ 10.0, 4.0, 8.0 ], " *
              "[ 14.0, 4.0, 8.0 ], " *
              "[ 14.0, 0.0, 8.0 ], " *
              "[ 10.0, 0.0, 8.0 ] ], " *
              "[ " *
              "[ 11.0, 1.0, 8.0 ], " *
              "[ 13.0, 1.0, 8.0 ], " *
              "[ 13.0, 3.0, 8.0 ], " *
              "[ 11.0, 3.0, 8.0 ], " *
              "[ 11.0, 1.0, 8.0 ] ] ] } ] }"

        AG.createmultilinestring([[
            [1.0, 4.0],
            [2.0, 5.0],
            [3.0, 6.0],
            [1.0, 4.0],
        ]]) do geom4
            @test AG.toWKT(geom4) == "MULTILINESTRING ((1 4,2 5,3 6,1 4))"
            @test AG.toWKT(AG.polygonfromedges(geom4, 0.1)) ==
                  "POLYGON ((1 4,2 5,3 6,1 4))"
            AG.polygonfromedges(geom4, 0.1) do geom5
                @test AG.toWKT(geom5) == "POLYGON ((1 4,2 5,3 6,1 4))"
            end
        end

        @test AG.getgeomtype(AG.getgeom(geom3, 0)) == AG.wkbPoint25D
        @test AG.getgeomtype(AG.getgeom(geom3, 1)) == AG.wkbPolygon25D
        @test AG.getgeomtype(AG.getgeom(geom3, 2)) == AG.wkbPolygon25D
        @test sprint(print, AG.getgeom(geom3, 3)) == "NULL Geometry"
        @test sprint(print, AG.getgeom(AG.IGeometry(), 3)) == "NULL Geometry"
        AG.getgeom(geom3, 0) do geom4
            @test AG.getgeomtype(geom4) == AG.wkbPoint25D
        end
        AG.getgeom(geom3, 1) do geom4
            @test AG.getgeomtype(geom4) == AG.wkbPolygon25D
        end
        AG.getgeom(geom3, 2) do geom4
            @test AG.getgeomtype(geom4) == AG.wkbPolygon25D
        end
        AG.getgeom(geom3, 3) do geom4
            @test sprint(print, geom4) == "NULL Geometry"
        end
        AG.getgeom(AG.IGeometry(), 0) do geom
            @test sprint(print, geom) == "NULL Geometry"
        end
    end

    @testset "Spatial Reference Systems" begin
        @test sprint(print, AG.getspatialref(AG.IGeometry())) ==
              "NULL Spatial Reference System"
        AG.getspatialref(AG.IGeometry()) do spatialref
            @test sprint(print, spatialref) == "NULL Spatial Reference System"
        end

        AG.createpoint(100, 70, 0) do geom
            @test sprint(print, AG.getspatialref(geom)) ==
                  "NULL Spatial Reference System"
            AG.getspatialref(geom) do spatialref
                @test sprint(print, spatialref) ==
                      "NULL Spatial Reference System"
            end
        end

        AG.read("data/point.geojson") do dataset
            layer = AG.getlayer(dataset, 0)
            AG.nextfeature(layer) do feature
                geom = AG.getgeom(feature)
                @test AG.toPROJ4(AG.getspatialref(geom)) ==
                      "+proj=longlat +datum=WGS84 +no_defs"
                AG.getspatialref(geom) do spatialref
                    @test AG.toPROJ4(spatialref) ==
                          "+proj=longlat +datum=WGS84 +no_defs"
                end
            end
            AG.createpoint(1, 2) do point
                @test sprint(print, AG.getspatialref(point)) ==
                      "NULL Spatial Reference System"
                AG.getspatialref(point) do spatialref
                    @test sprint(print, spatialref) ==
                          "NULL Spatial Reference System"
                end
            end
        end

        AG.importEPSG(2927) do source
            AG.importEPSG(4326) do target
                AG.createcoordtrans(source, target) do transform
                    AG.fromWKT("POINT (1120351.57 741921.42)") do point
                        @test AG.toWKT(point) == "POINT (1120351.57 741921.42)"
                        AG.transform!(point, transform)
                        @test GeoInterface.coordinates(point) ≈
                              [47.3488070138318, -122.5981499431438]
                    end
                end
            end
        end
    end

    @testset "Cloning NULL geometries" begin
        geom = AG.IGeometry()
        @test AG.geomname(geom) === missing
        @test sprint(print, AG.clone(geom)) == "NULL Geometry"
        AG.clone(geom) do g
            @test sprint(print, g) == "NULL Geometry"
        end
    end
end


text:
using Test
using WordTokenizers

@testset "Targetted" begin
    @testset "Initials" begin
        @test 1 == length(rulebased_split_sentences("It is by A. Adamson, the famous author."))
        @test 1 == length(rulebased_split_sentences("It is by Z. Zeckerson, the famous author."))
        @test 1 == length(rulebased_split_sentences("It is by Bill R. Emerson, the famous author."))
        @test 1 == length(rulebased_split_sentences("It is by A. A. Adamson, the famous author."))

        @test_broken 1 == (rulebased_split_sentences("It is by Ian I. Irving, the famous author."))

        @test 2 == length(rulebased_split_sentences("He doesn't and nor will I. It is best this way."))
    end

    @testset "period only" begin
        @test length(rulebased_split_sentences("a good day . . Yes"))==2
    end

    @testset "Acronyms" begin
        @test length(rulebased_split_sentences("Adamson is not from USA. They are from Europe."))==2
    end
end

@testset "Simple" begin
    split_sentences("Never going to give you up. Never going to let you down.") ==
        ["Never going to give you up.", "Never going to let you down."]
end

function test_sentence_splitting(split_sentences, raw)
    sents = split(raw, "\n") # Originally 1 sentence per line
    task = join(sents, " ") # Now all one lines
    output = strip.(split_sentences(task)) #strip off the extra spaces added

    @test output == sents
end

@testset "turtles" begin
    #  https://simple.wikipedia.org/wiki/Turtle ,
    doc1 = """Turtles are the reptile order Testudines.
    They have a special bony or cartilaginous shell developed from their ribs that acts as a shield.
    The order Testudines includes both living and extinct species.
    The earliest fossil turtles date from about 220 million years ago.
    So turtles are one of the oldest surviving reptile groups and a more ancient group than lizards, snakes and crocodiles.
    Turtle have been very successful, and have almost world-wide distribution.
    But, of the many species alive today, some are highly endangered."""

    #  https://simple.wikipedia.org/wiki/Sea_turtle ,
    doc2 = """Sea turtles (Chelonioidea) are turtles found in all the world's oceans except the Arctic Ocean, and some species travel between oceans.
    The term is US English.
    In British English they are simply called "turtles"; fresh-water chelonians are called "terrapins" and land chelonians are called tortoises.
    There are seven types of sea turtles: Kemp's Ridley, Flatback, Green, Olive Ridley, Loggerhead, Hawksbill and the leatherback.
    All but the leatherback are in the family Chelonioidea.
    The leatherback belongs to the family Dermochelyidae and is its only member.
    The leatherback sea turtle is the largest, measuring six or seven feet (2 m) in length at maturity, and three to five feet (1 to 1.5 m) in width, weighing up to 2000 pounds (about 900 kg).
    Most other species are smaller, being two to four feet in length (0.5 to 1 m) and proportionally less wide.
    The Flatback turtle is found solely on the northern coast of Australia."""

    #  https://simple.wikipedia.org/wiki/Green_turtle
    doc3 = """Chelonia mydas, commonly known as the green turtle, is a large sea turtle belonging to the family Cheloniidae.
    It is the only species in its genus.
    It is one of the seven marine turtles, which are all endangered.
    Although it might have some green on its carapace (shell), the green turtle is not green.
    It gets its name from the fact that its body fat is green.
    It can grow up to 1 m (3 ft) long and weigh up to 160 kg (353 lb).
    They are an endangered species, especially in Florida and the Pacific coast of Mexico.
    They can also be found in warm waters around the world and are found along the coast of 140 countries.
    The female turtle lays eggs in nests she builds in the sand on the beaches.
    She uses the same beach that she was born on.
    During the nesting season in summer she can make up to five nests.
    She can lay as many as 135 eggs in a nest.
    The eggs take about two months to hatch.
    The baby turtles are about 50 mm (2 in) in length."""

    test_sentence_splitting(rulebased_split_sentences, doc1)
    test_sentence_splitting(rulebased_split_sentences, doc2)
    test_sentence_splitting(rulebased_split_sentences, doc3)
end

@testset "Difficult" begin
    test_sentence_splitting(rulebased_split_sentences,
        """Punkt knows that the periods in Mr. Smith and Johann S. Bach do not mark sentence boundaries.
        And sometimes sentences can start with non-capitalized words.
        i is a good variable name.""")
end

ext: jl


text:
"""
$(SIGNATURES)

Shakhov non-equilibrium part

* @arg: particle velocity quadrature points
* @arg: discrete Maxwellian
* @arg: primitive variables, Prandtl number, heat flux, inner degree of freedom
"""
function shakhov(u::AV{X}, M::AV{Y}, q, prim::AV{Z}, Pr) where {X<:FN,Y<:FN,Z<:Real} # 1F1V
    M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 - 5.0) *
       M

    return M_plus
end

#--- 2F1V ---#
function shakhov(
    u::AV{X},
    H::Y,
    B::Y,
    q,
    prim::AV{Z},
    Pr,
    K,
) where {X<:FN,Y<:AV{<:FN},Z<:Real}

    H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 + K - 5.0) *
       H
    B_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 + K - 3.0) *
       B

    return H_plus, B_plus

end

#--- 1F2V ---#
function shakhov(
    u::T,
    v::T,
    M::AM{X},
    q::AV{Y},
    prim::AV{Z},
    Pr,
) where {T<:AA{<:FN,2},X<:FN,Y<:FN,Z<:Real}

    M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) - 5.0) *
       M

    return M_plus

end

#--- 2F2V ---#
function shakhov(
    u::T,
    v::T,
    H::X,
    B::X,
    q::AV{Y},
    prim::AV{Z},
    Pr,
    K,
) where {T<:AA{<:FN,2},X<:AA{<:FN,2},Y<:Real,Z<:Real}

    H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 5.0) *
       H
    B_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 3.0) *
       B

    return H_plus, B_plus

end

#--- 1F3V ---#
function shakhov(
    u::T,
    v::T,
    w::T,
    M::AA{X,3},
    q::AV{Y},
    prim::AV{Z},
    Pr,
) where {T<:AA{<:FN,3},X<:FN,Y<:Real,Z<:Real}

    M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2] + (w - prim[4]) * q[3]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) - 5.0) *
       M

    return M_plus

end

"""
    shakhov!(S, u, M, q, prim, Pr)
    shakhov!(SH, SB, u, H, B, q, prim, Pr, K)
    shakhov!(S, u, v, M, q, prim, Pr)
    shakhov!(SH, SB, u, v, H, B, q, prim, Pr, K)
    shakhov!(S, u, v, w, M, q, prim, Pr)

In-place Shakhov non-equilibrium part

* @arg: particle velocity quadrature points
* @arg: discrete Maxwellian
* @arg: primitive variables, Prandtl number, heat flux, inner degree of freedom

"""
function shakhov!(S::T1, u::AV{T2}, M::T1, q, prim, Pr) where {T1<:AA{<:FN,1},T2<:FN} # 1F1V
    @. S = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 - 5.0) *
       M

    return nothing
end

#--- 2F1V ---#
function shakhov!(
    SH::T1,
    SB::T1,
    u::AV{T2},
    H::T1,
    B::T1,
    q,
    prim,
    Pr,
    K,
) where {T1<:AA{<:FN,1},T2<:FN}

    @. SH =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        (u - prim[2]) *
        q *
        (2.0 * prim[end] * (u - prim[2])^2 + K - 5.0) *
        H
    @. SB =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        (u - prim[2]) *
        q *
        (2.0 * prim[end] * (u - prim[2])^2 + K - 3.0) *
        B

    return nothing

end

#--- 1F2V ---#
function shakhov!(
    S::T1,
    u::T2,
    v::T2,
    M::T1,
    q::AV{T3},
    prim,
    Pr,
) where {T1<:AA{<:FN,2},T2<:AA{<:FN,2},T3<:Real}

    @. S =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) - 5.0) *
        M

    return nothing

end

#--- 2F2V ---#
function shakhov!(
    SH::T1,
    SB::T1,
    u::T2,
    v::T2,
    H::T1,
    B::T1,
    q::AV{T3},
    prim,
    Pr,
    K,
) where {T1<:AA{<:FN,2},T2<:AA{<:FN,2},T3<:Real}

    @. SH =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 5.0) *
        H
    @. SB =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 3.0) *
        B

    return nothing

end

#--- 1F3V ---#
function shakhov!(
    S::T1,
    u::T2,
    v::T2,
    w::T2,
    M::T1,
    q::AV{T3},
    prim,
    Pr,
) where {T1<:AA{<:FN,3},T2<:AA{<:FN,3},T3<:Real}

    @. S =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2] + (w - prim[4]) * q[3]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) - 5.0) *
        M

    return nothing

end


text:
export Input

mutable struct Input end

function handleMessage!(i::Input, m::TICK)

  Log.@debug "Input tick"
  evt_ref = Ref{SDL_Event}()
  while Bool(SDL_PollEvent(evt_ref))
    evt = evt_ref[]
    ty = evt.type
    if ty ∈ [SDL_AUDIODEVICEADDED, SDL_AUDIODEVICEREMOVED]
      sendMessage(evt.adevice)
    elseif ty == SDL_CONTROLLERAXISMOTION
      sendMessage(evt.caxis)
    elseif ty ∈ [SDL_CONTROLLERBUTTONDOWN, SDL_CONTROLLERBUTTONUP]
      sendMessage(evt.cbutton)
    elseif ty ∈ [SDL_CONTROLLERDEVICEADDED, SDL_CONTROLLERDEVICEREMOVED, SDL_CONTROLLERDEVICEREMAPPED]
      sendMessage(evt.cdevice)
    elseif ty ∈ [SDL_DOLLARGESTURE, SDL_DOLLARRECORD]
      sendMessage(evt.dgesture)
    elseif ty ∈ [SDL_DROPFILE, SDL_DROPTEXT, SDL_DROPBEGIN, SDL_DROPCOMPLETE]
      sendMessage(evt.drop)
    elseif ty ∈ [SDL_FINGERMOTION, SDL_FINGERDOWN, SDL_FINGERUP]
      sendMessage(evt.tfinger)
    elseif ty ∈ [SDL_KEYDOWN, SDL_KEYUP]
      sendMessage(evt.key)
    elseif ty == SDL_JOYAXISMOTION
      sendMessage(evt.jaxis)
    elseif ty == SDL_JOYBALLMOTION
      sendMessage(evt.jball)
    elseif ty == SDL_JOYHATMOTION
      sendMessage(evt.jhat)
    elseif ty ∈ [SDL_JOYBUTTONDOWN, SDL_JOYBUTTONUP]
      sendMessage(evt.jbutton)
    elseif ty ∈ [SDL_JOYDEVICEADDED, SDL_JOYDEVICEREMOVED]
      sendMessage(evt.jdevice)
    elseif ty == SDL_MOUSEMOTION
      sendMessage(evt.motion)
    elseif ty ∈ [SDL_MOUSEBUTTONDOWN, SDL_MOUSEBUTTONUP]
      sendMessage(evt.button)
    elseif ty == SDL_MOUSEWHEEL
      sendMessage(evt.wheel)
    elseif ty == SDL_MULTIGESTURE
      sendMessage(evt.mgesture)
    elseif ty == SDL_QUIT
      sendMessage(evt.quit)
    elseif ty == SDL_SYSWMEVENT
      sendMessage(evt.syswm)
    elseif ty == SDL_TEXTEDITING
      sendMessage(evt.edit)
    elseif ty == SDL_TEXTINPUT
      sendMessage(evt.text)
    elseif ty == SDL_USEREVENT
      sendMessage(evt.user)
    elseif ty == SDL_WINDOWEVENT
      sendMessage(evt.window)
    else
      sendMessage(evt)
    end
  end
end

function awake!(i::Input)
  Log.@debug "Input awake!"
  listenFor(i, TICK)
end

function shutdown!(i::Input)
  Log.@debug "Input shutdown!"
  unlistenFrom(i, TICK)
end


text:
"""
    TypeT()

Lorem ipsum dolor sit amet.
"""
struct TypeT

end


text:
"""
    type Cardinality

A trait describing the cardinality of a [`Transform`]. Available cardinalities are:
[`OneToOne`](@ref), [`ManyToOne`](@ref), [`OneToMany`](@ref), and [`ManyToMany`](@ref).
"""
abstract type Cardinality end

"""
    OneToOne <: Cardinality

Transforms that map each input to exactly one output: `x → y`.
Examples: [`Power`](@ref), [`Periodic`](@ref).
"""
struct OneToOne <: Cardinality end

"""
    ManyToOne <: Cardinality

Transforms that map many inputs to one output: `(x_1, x_2, ..., x_n) → y`.
These are typically reduction operations.
Examples: [`LinearCombination`](@ref).
"""
struct ManyToOne <: Cardinality end

"""
    OneToMany <: Cardinality

Transforms that map one input to many outputs: `x → (y_1, y_2, ..., y_n)`.
Examples: [`OneHotEncoding`](@ref).
"""
struct OneToMany <: Cardinality end

"""
    ManyToMany <: Cardinality

Transforms that map many inputs to many outputs: `(x_1, x_2, ..., x_m) → (y_1, y_2, ..., y_n)`.
Examples: Principle Component Analysis (not implemented).
"""
struct ManyToMany <: Cardinality end

"""
    cardinality(transform) -> Cardinality

Returns the [`Cardinality`](@ref) of the `transform`.
"""
function cardinality end


text:
#include("../model0.jl")
include("../computation/buffers0.jl")
include("../kernels/kernel.jl")
using Cubature
using TimerOutputs

if !@isdefined(to)
    const to = TimerOutput()
end

"""
Compute cdf, pdf, and pdf_deriv of T-distribution
"""
function comp_tdist(btg::btg, θ::Union{Array{T, 1}, T} where T<:Real, λ::Real; validate = 0)
    g = getfield(Main, Symbol(options.transform_type))() #nonlinear transform, e.g. BoxCox
    invg = (x, λ) -> inverse(g, x, λ)
    (x, Fx, y, _, n, p) = unpack(btg.trainingData) #unpack training data
    @assert validate >= 0 && n >= validate
    n = validate == 0 ? n : n-1 #downsize by 1, if we will be deleting 1 data point
    θλpair = (θ, λ)
    choleskyΣθ = btg.train_buffer_dict[θ].choleskyΣθ
    if validate == 0
        #retrieve qtilde, gλy, βhat, Σθ_inv_y
        βhat = btg.θλbuffer_dict[θλpair].βhat
        qtilde = btg.θλbuffer_dict[θλpair].qtilde
        @assert qtilde > 0
        Σθ_inv_y = btg.θλbuffer_dict[θλpair].Σθ_inv_y
        #(_, _, βhat, qtilde, _, Σθ_inv_y) = unpack(btg.θλbuffer_dict[θλpair]) #unpack can be dangerous, because if the order of elements change and you don't know, you're screwed. safer tp access directly.
        #(_, _, _, choleskyΣθ, _, _, _) = unpack(btg.train_buffer_dict[θ])
        #(_, gλy, _) = unpack(btg.λbuffer_dict[λ])
    else
        @warn "doing cross validation in tdist"
        #cur_val_θλ_buf = btg.validation_θλ_buffer_dict[θλpair]
        #qtilde = cur_val_θλ_buf.
        (_, _, _, βhat, qtilde, Σθ_inv_y)  = unpack(btg.validation_θλ_buffer_dict[θλpair])  #depends on theta and lambda
        #retrieve qtilde_minus_i, etc.
    end
    dg = (y, λ) -> partialx(g, y, λ) #derivative w.r.t z
    dg2 = (y, λ) -> partialxx(g, y, λ) #second derivative w.r.t z
    jac = x -> abs(reduce(*, map(z -> dg(z, λ), x))) #Jacobian function
    #println("unpack in tdist: ", unpack)
    function compute_qmC(x0, Fx0) #assume q and C are 1x1
        x0 = reshape(x0, 1, length(x0)) #reshape into 1 x d vector
        x1, y1 = getDimension(btg.trainingData), getDimension(btg.testingData)
        x2, y2 = getCovDimension(btg.trainingData), getCovDimension(btg.testingData)
        #@info "validate" validate
        if validate == 0
            #@warn "not in validation mode"
            @timeit to "update testingData" update!(btg.testingData, x0, Fx0)# update testing data with x0, Fx0 with USER INPUT
            @timeit to "update Dθ, Hθ, Cθ,... " update!(btg.train_buffer_dict[θ], btg.test_buffer_dict[θ], btg.trainingData, btg.testingData) #update θ-testing buffer (which holds Dθ, Hθ, Cθ,...) with testingData
            #print((btg.test_buffer_dict[θ]))
            #######       unpack throws an odd error here, still need to figure out why      ##############################
            @timeit to "load Bθ, Hθ, Cθ, Dθ, Eθ" begin
                Bθ = btg.test_buffer_dict[θ].Bθ
                Hθ = btg.test_buffer_dict[θ].Hθ
                Cθ = btg.test_buffer_dict[θ].Cθ
                Dθ = btg.test_buffer_dict[θ].Dθ
                Eθ = btg.test_buffer_dict[θ].Eθ
                ΣθinvBθ = btg.test_buffer_dict[θ].ΣθinvBθ
            end

            @timeit to "BθΣθinvBθ, m, qC" begin
                BθΣθinvBθ = Bθ*ΣθinvBθ
                #(_, Bθ, _, _, Hθ, Cθ) = anotherone(btg.test_buffer_dict[θ])
                #(_, Bθ, _, _, Hθ, Cθ) = unpack(btg.test_buffer_dict[θ])
                ###############################################################
                #println("past the breakpoint")
                #@info "θ", θ
                m = Bθ*Σθ_inv_y + Hθ*βhat #recompute mean
                #if m[1] < 0
                #    @info "in tdist computeqMC..."
                #    @info "Bθ", Bθ
                #    @info "Σθ_inv_y", Σθ_inv_y
                #    @info "Hθ", Hθ
                #    @info "βhat", βhat
                #end
                qC = qtilde[1]*Cθ[1] #both q and C are 1x1 for single-point prediction
                # sigma_m = qC/(n-p-2) + m[1]^2 # E[T_i^2] for quantile estimation
                # temporary definition for testing
                expr2 = [1]
            end

        else # want to do validation instead
            @info "validate", validate
            @warn "doing cross validation in tdist"
            Fx_i = btg.trainingData.Fx[validate:validate, :] #ith row of trainingData
            Fx_minus_i = btg.trainingData.Fx[[1:validate-1; validate+1:end], :] #all but ith row of trainingData
            x_i = btg.trainingData.x[validate:validate, :]
            x_minus_i = btg.trainingData.x[[1:validate-1; validate+1:end], :]
            update!(btg.testingData, x_i, Fx_i) #testingData is subset of trainingData. N.B: USER INPUT FOR x0, Fx0 are IGNORED
            update!(btg.train_buffer_dict[θ], btg.test_buffer_dict[θ], btg.trainingData, btg.testingData) #update θ-testing buffer (which holds Dθ, Hθ, Cθ,...) with testingData
            #update!(btg.validation_test_buffer_dict[θ], btg.train_buffer_dict[θ]::train_buffer, btg.test_buffer_dict[θ]::test_buffer, validate)
            Σθ_inv_X_minus_i = btg.validation_train_buffer_dict[θ].Σθ_inv_X_minus_i
            ΣθinvBθ = btg.test_buffer_dict[θ].ΣθinvBθ
            ΣθinvBθ_minus_i = lin_sys_loocv_IC(ΣθinvBθ, choleskyΣθ, validate)

            (Eθ, Bθ, ΣθinvBθ, Dθ, Hθ, Cθ) = unpack(btg.test_buffer_dict[θ]) #a critical assumption is that the covariates Fx0 remain constant throughout cross-validation

             Bθ_minus_i = @view Bθ[:, [1:validate-1; validate+1:end]] #discard ith entry
             Dθ_minus_i = Eθ - Bθ_minus_i * ΣθinvBθ_minus_i #Eθ is going to be 1 for

             Hθ_minus_i = Fx_i - Bθ_minus_i * Σθ_inv_X_minus_i
             expr2 = Hθ_minus_i*((Fx_minus_i'*Σθ_inv_X_minus_i)\Hθ_minus_i')
             Cθ = Dθ_minus_i + expr2 #O(p^3) operation, not bad if p is small/ How to update choleskyXΣX to get choleskyXΣX_minus_i?

             #we refrain from calling the last expression Cθ_minus_i so it can be returned by this func as usual
             m = Bθ_minus_i * Σθ_inv_y + Hθ_minus_i * βhat

            #temporary definition for testing
            Dθ = Dθ_minus_i
            #@info "ΣθinvBθ_minus_i", ΣθinvBθ_minus_i
            #@info "Eθ", Eθ
            #@info Bθ
            #@info Σθ_inv_y
            #@info Hθ
            #@info βhat
        end
        #gλz = btg.λbuffer_dict[λ].λ
        sigma = sqrt(qtilde[1]*Cθ[1]/(n-p))
        return m[1], qtilde[1], Cθ[1], sigma, βhat, Dθ, expr2, validate,λ, BθΣθinvBθ, Eθ
    end

    function compute(f, x0, Fx0, y0)#updates testingData and test_buffer, but leaves train_buffer and trainingData alone
        @timeit to "compute q, m, C, etc." m, q, C, sigma, βhat, Dθ, expr2, validate, λ, BθΣθinvBθ, Eθ  = compute_qmC(x0, Fx0)
        @assert typeof(validate) <: Int64
        #@warn "pushing to debug log"
        #push!(btg.debug_log, (m, C, q))
        qC = q*C
        @timeit to "try-catch statement" begin
        try
            @assert Dθ[1]>0
            @assert expr2[1]>0
            @assert q[1] > 0
            @assert C[1] > 0
            @assert n-p > 0
        catch e
            @info "θ", θ
            @info "D_theta", Dθ[1]
            @info "BθΣθinvBθ" BθΣθinvBθ
            @info "Eθ" Eθ
            @info "H(XsigmaX)'H", expr2[1]
            @info "q", q
            @info "C", C
            @info "validate", validate
            @info "logdet choleskyΣθ", logdet(choleskyΣθ)
        end
        end
        #@assert m>=0
        t = LocationScale(m, sqrt(qC/(n-p)), TDist(n-p)) #avail ourselves of built-in tdist
        #@info "n", n; @info "p", p; @info "g(y0, λ)", g(y0, λ);
        #return (g(y0, λ) .- m)/sqrt(qC/(n-p))
        return f(y0, t, m, qC) # return location parameter to utilize T mixture structure
        #return C
        #return m
        #(_, _, _, _, n, p) = unpack(btg.trainingData)
        #return f
        #return  sqrt(qC/(n-p))
        #return m
        #return m
        #return g(y0, λ)
    end

    function main_pdf(y0, t, m, qC)
        #@assert y0 >= 0
        # condition tdistpdf on assumption that data is positive
        temp = Distributions.pdf.(t, g(y0, λ)) * jac(y0)
        if g == BoxCox() || g == ShiftedBoxCox()
            if λ > 0 # g(0) finite, truncated from left
                a = Distributions.cdf.(t, -1/λ)
                temp = temp/(1 - a)
            elseif λ < 0 # g(+inf) finite, truncated from right
                a = Distributions.cdf.(t, -1/λ)
                temp = temp/a
            end
        elseif g == IdentityTransform()
            a = Distributions.cdf.(t, 0)
            temp = temp/(1 - a)
        else # can add other transform here
            throw(ArgumentError("Type of transform not supported."))
        end
        return temp
        # @timeit to "atomic pdf eval" Distributions.pdf.(t, g(y0, λ)) * jac(y0)
    end
    function main_cdf(y0, t, m, qC)
        #@assert y0 >= 0
        temp = Distributions.cdf.(t, g(y0, λ))
        if g == BoxCox() || g == ShiftedBoxCox()
            if λ > 0 # g(0) finite, truncated from left
                a = Distributions.cdf.(t, -1/λ)
                temp = (temp-a)/(1 - a)
            elseif λ < 0 # g(+inf) finite, truncated from right
                a = Distributions.cdf.(t, -1/λ)
                temp = temp/a
            end
        elseif g == IdentityTransform()
            a = Distributions.cdf.(t, 0)
            temp = (temp-a)/(1 - a)
        else # can add other transform here
            throw(ArgumentError("Type of transform not supported."))
        end
        return temp
        #@timeit to "atomic cdf eval" return Distributions.cdf.(t, g(y0, λ))
    end
    function main_pdf_deriv_helper(y0, t, m, qC)
        #@assert y0 >= 0
        (k = size(qC, 1); gλy0 = g(y0, λ); Ty0 = Distributions.pdf.(t, gλy0);
                            Ty0 .* (-(n-p+k)) .* ( gλy0 .- m) ./ (qC .+ (gλy0 .- m) .^2) .* dg(y0, λ)) #this is a stable computation of the derivative of the tdist
    end

    function main_pdf_deriv(y0, t, m, qC)
        #@assert y0>=0
        (gλy0 = g(y0, λ); Ty0 = Distributions.pdf.(t, gλy0); (dg2(y0, λ) .* Ty0 .+ abs.(dg(y0, λ)) .* main_pdf_deriv_helper(y0, t, m, qC))[1])
    end

    function pdf_deriv(x0, Fx0, y0)
        #@assert y0 >= 0
        @timeit to "atomic pdf deriv" compute(main_pdf_deriv, x0, Fx0, y0)
    end
    function pdf(x0, Fx0, y0)
        #@assert y0>=0
        return compute(main_pdf, x0, Fx0, y0)
    end
    function cdf(x0, Fx0, y0)
        #@assert y0 >= 0
        return compute(main_cdf, x0, Fx0, y0)
    end

    """
    Gradient of CDF w.r.t augmented (u, s) vector (u is for value, and y is for location)
    """
    function cdf_grad_us(x0, Fx0, y0)
        @timeit to "m, q, c" m, q, C, βhat = compute_qmC(x0, Fx0)
        @timeit to "cdf jacobian" begin
        qC = q*C
        k = size(qC, 1)
        gλy0 = g(y0, λ)
        #dgλy0 = partialx(g, y0, λ)
        jac = x -> abs(reduce(*, map(z -> dg(z, λ), x))) #Jacobian function
        jacy0 = jac(y0)
        ## N.B: please note the difference between t and vanilla t. By default, the build-in julia tdist pdf
        ## will divide by the covariance, which isn't what we want if we are differentiating with respect to, say
        ## the location s. Since we still want to take advantage of the built-in tdist, we create a tdist with
        ## unit covariance and put in the argument arg ourselves, which corrects for mean and covariance.
        arg = ((gλy0 .- m)/sqrt(qC/(n-p)))[1] #1 x 1
        t = LocationScale(m, sqrt(qC/(n-p)), TDist(n-p))
        vanillat = LocationScale(0, 1, TDist(n-p))
        cdf_eval =  Distributions.cdf.(vanillat, arg)
        cdf_deriv = Distributions.pdf.(vanillat, arg) #chain rule term is computed later
        cdf_deriv_alternate = Distributions.pdf.(t, gλy0)
        #standard unpacking
        x0 = reshape(x0, 1, length(x0))
        #unpack pertinent quantities
        Σθ_inv_X  = btg.train_buffer_dict[θ].Σθ_inv_X
        choleskyΣθ = btg.train_buffer_dict[θ].choleskyΣθ
        choleskyXΣX = btg.train_buffer_dict[θ].choleskyXΣX
        (x, Fx, y, d, n, p) = unpack(btg.trainingData)

        ## update block
        @timeit to "update testingData" update!(btg.testingData, x0, Fx0)# update testing data with x0, Fx0 with USER INPUT
        @timeit to "update test buffer" update!(btg.train_buffer_dict[θ], btg.test_buffer_dict[θ], btg.trainingData, btg.testingData)  #attempt to update test buffer
        (Eθ, Bθ, ΣθinvBθ, Dθ, Hθ, Cθ) = unpack(btg.test_buffer_dict[θ])
        #@timeit to "intermediates_old" (jacC, jacm, jacB, jacH, jacD) = intermediates_old(x0, Bθ, Hθ, Σθ_inv_X, choleskyXΣX)
        @timeit to "intermediates_new" (jacC, jacm, jacB, jacH, jacD) = intermediates(btg.jac_buffer_dict[θ], btg.testingData, btg.trainingData, btg.test_buffer_dict[θ],
        btg.train_buffer_dict[θ], btg.θλbuffer_dict[θλpair])

        #we need to define covariance (sigma here ...)
        #sigma is covariance sqrt(qC/(n-p))
        qC = (q .* Cθ[1])
        sigma = sqrt(qC/(n-p)) #covariance
        dsigma = sqrt(q/(n-p)) .* 0.5 * Cθ[1]^(-1/2) .* jacC
        function Y(i, sigma, dsigma)
            return - jacm[i]/sigma - (gλy0 - m)/sigma^2 * dsigma[i]
        end
        jacG = zeros(1, d) #1 x d
        for i = 1:d
            jacG[i] = (cdf_deriv .* Y(i, sigma, dsigma))[1]
        end
        if false #check derivative of arg
            check_arg = zeros(1, d)
            for i = 1:d
                check_arg[i] = Y(i, sigma, dsigma)
            end
        end
        #@info "n", n; @info "p", p; @info "gλy0", gλy0
        #return check_arg
        #return jacG
        #return hcat(cdf_deriv * dgλy0, jacG)
        #@info jacG
        #return dsigma
        #return jacD
        #return jacG'
        #return jacm
    end
        return hcat(cdf_deriv_alternate*jacy0, jacG)
        #return cdf_deriv_alternate*jacy0
        #return jacC'
        #return jacm
        #return [dgλy0]
        #return cdf_deriv
    end

    """
    Gradient of CDF w.r.t augmented (u, s) vector (u is for value, and y is for location)
    """
    function cdf_grad_us!(x0, Fx0, y0, store)
        m, q, C, βhat = compute_qmC(x0, Fx0)
        qC = q*C
        k = size(qC, 1)
        gλy0 = g(y0, λ)
        #dgλy0 = partialx(g, y0, λ)
        jac = x -> abs(reduce(*, map(z -> dg(z, λ), x))) #Jacobian function
        jacy0 = jac(y0)
        ## N.B: please note the difference between t and vanilla t. By default, the build-in julia tdist pdf
        ## will divide by the covariance, which isn't what we want if we are differentiating with respect to, say
        ## the location s. Since we still want to take advantage of the built-in tdist, we create a tdist with
        ## unit covariance and put in the argument arg ourselves, which corrects for mean and covariance.
        arg = ((gλy0 .- m)/sqrt(qC/(n-p)))[1] #1 x 1
        t = LocationScale(m, sqrt(qC/(n-p)), TDist(n-p))
        vanillat = LocationScale(0, 1, TDist(n-p))
        cdf_eval =  Distributions.cdf.(vanillat, arg)
        cdf_deriv = Distributions.pdf.(vanillat, arg) #chain rule term is computed later
        cdf_deriv_alternate = Distributions.pdf.(t, gλy0)
        #standard unpacking
        x0 = reshape(x0, 1, length(x0))

        #unpack pertinent quantities
        Σθ_inv_X  = btg.train_buffer_dict[θ].Σθ_inv_X
        choleskyΣθ = btg.train_buffer_dict[θ].choleskyΣθ
        choleskyXΣX = btg.train_buffer_dict[θ].choleskyXΣX
        (x, Fx, y, d, n, p) = unpack(btg.trainingData)

         ## update block
         @timeit to "update testingData" update!(btg.testingData, x0, Fx0)# update testing data with x0, Fx0 with USER INPUT
         @timeit to "update test buffer" update!(btg.train_buffer_dict[θ], btg.test_buffer_dict[θ], btg.trainingData, btg.testingData)  #attempt to update test buffer
         (Eθ, Bθ, ΣθinvBθ, Dθ, Hθ, Cθ) = unpack(btg.test_buffer_dict[θ])
         #@timeit to "intermediates_old" (jacC, jacm, jacB, jacH, jacD) = intermediates_old(x0, Bθ, Hθ, Σθ_inv_X, choleskyXΣX)
         @timeit to "intermediates_new" (jacC, jacm, jacB, jacH, jacD) = intermediates(btg.jac_buffer_dict[θ], btg.testingData, btg.trainingData, btg.test_buffer_dict[θ],
         btg.train_buffer_dict[θ], btg.θλbuffer_dict[θλpair])

       #we need to define covariance (sigma here ...)
        #sigma is covariance sqrt(qC/(n-p))
        qC = (q .* Cθ[1])
        sigma = sqrt(qC/(n-p)) #covariance
        dsigma = sqrt(q/(n-p)) .* 0.5 * Cθ[1]^(-1/2) .* jacC
        function Y(i, sigma, dsigma)
            return - jacm[i]/sigma - (gλy0 - m)/sigma^2 * dsigma[i]
        end
        #jacG = zeros(1, d) #1 x d
        @assert length(store) == d+1
        store[1] = cdf_deriv_alternate*jacy0
        for i = 1:d
            store[i+1] = (cdf_deriv .* Y(i, sigma, dsigma))[1]
        end
        check_arg = zeros(1, d)
        for i = 1:d
            check_arg[i] = Y(i, sigma, dsigma)
        end
        #return check_arg
        #return jacG
        return hcat(cdf_deriv * dgλy0, jacG)
        #@info jacG
        #return dsigma
        #return jacG'
        #return hcat(cdf_deriv_alternate*jacy0, jacG)
        #return nothing
        #return jacD
        #return cdf_deriv_alternate*jacy0
        #return jacC'
        #return jacm
        #return [dgλy0]
        #return cdf_deriv
    end

    function cdf_hess_us(x0, Fx0, y0)
        @timeit to "cdf hessian eval" begin
        m, q, C, βhat = compute_qmC(x0, Fx0)
        qC = q*C
        k = size(qC, 1)
        gλy0 = g(y0, λ)
        dgλy0 = dg(y0, λ)
        d2gλy0 = dg2(y0, λ)
        x0 = reshape(x0, 1, length(x0))
        #unpack pertinent quantities
        Σθ_inv_X  = btg.train_buffer_dict[θ].Σθ_inv_X
        choleskyΣθ = btg.train_buffer_dict[θ].choleskyΣθ
        choleskyXΣX = btg.train_buffer_dict[θ].choleskyXΣX
        (x, Fx, y, d, n, p) = unpack(btg.trainingData)
        @timeit to "update testingData" update!(btg.testingData, x0, Fx0)
        @timeit to "update Dθ, Hθ, Cθ,... " update!(btg.train_buffer_dict[θ], btg.test_buffer_dict[θ], btg.trainingData, btg.testingData)#attempt to update
        (Eθ, Bθ, ΣθinvBθ, Dθ, Hθ, Cθ) = unpack(btg.test_buffer_dict[θ])
        #(jacC, jacm, jacB, jacH, jacD) = intermediates(x0, Bθ, Hθ, Σθ_inv_X, choleskyXΣX, ΣθinvBθ)
        (jacC, jacm, jacB, jacH, jacD) = intermediates(btg.jac_buffer_dict[θ], btg.testingData, btg.trainingData, btg.test_buffer_dict[θ],
        btg.train_buffer_dict[θ], btg.θλbuffer_dict[θλpair])
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(compute hessian of D)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ΣinvjacB  = choleskyΣθ\jacB
        W = zeros(d, d)
        S = repeat(x0, n, 1) .- x # n x d
        flattenedθ = typeof(θ)<:Real ? [θ for i=1:d] : θ[:]
        Stheta = S * diagm(- flattenedθ) # n x d
        for i = 1:d
            for j = 1:d
                arg = Bθ' .* Stheta[:, i:i] .* Stheta[:, j:j]
                if i==j
                    eq = Bθ * (choleskyΣθ \ (-flattenedθ[i] .* Bθ'))
                else
                    eq = 0
                end
                W[i, j] = (Bθ * (choleskyΣθ\arg) .+ eq)[1]
            end
        end
        hessD = -2*jacB' * ΣinvjacB - 2*W #d x d
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(compute hessian of C)~~~~~~~~~~~~~~~~~~~~~~~~~~~
        V = zeros(d, d)
        for i = 1:d
            for j = 1:d
                arg = Bθ' .* Stheta[:, i:i] .* Stheta[:, j:j]
                if i==j
                    eq = (-flattenedθ[i] .* Bθ')
                else
                    eq = 0
                end
                #println("size arg: ", size(arg))
                #println("size eq: ", size(eq))
                # println(arg .+ eq)
                # println((Fx' * (choleskyΣθ \ (arg .+ eq))))
                # println((choleskyXΣX \ (Fx' * (choleskyΣθ \ (arg .+ eq)))))
                # println(Hθ)
                V[i, j] =  (Hθ * (choleskyXΣX \ ( - Fx' * (choleskyΣθ \ (arg .+ eq)))))[1]
            end
        end
        hessC = hessD + 2*jacH*(choleskyXΣX\jacH') + 2*V #d x d
        #~~~~~~~~~~~~~~~~~~~~~~~~~(compute derivative and hessian of m)~~~~~~~~~~~~~~~~~~~~~~~~~~
        #@info "size(jacB')", size(jacB')
        #@info "size(Σθ_inv_y)",  size(Σθ_inv_y)
        #@info "size(jacH)", size(jacH)
        #@info " size(βhat)", size(βhat)
        jacm = jacB' * Σθ_inv_y + jacH * βhat  # d x 1
        hess_m = zeros(d, d)
        for i = 1:d
            for j = 1:d
                arg = Bθ' .* Stheta[:, i:i] .* Stheta[:, j:j]
                if i==j
                    eq = (-flattenedθ[i] .* Bθ')
                else
                    eq = 0
                end
                #println()
                #println("sigma inv y: ", size(Σθ_inv_y))
                #println("size arg .+ eq: ", size(arg .+ eq))
                expr1 = dot(Σθ_inv_y,  (arg .+ eq))
                expr2 = dot(( - Fx' * (choleskyΣθ \ (arg .+ eq))), βhat)
                #println("size expr1: ", size(expr1))
                #println("size expr2: ", size(expr2))
                hess_m[i, j] = (expr1 + expr2)[1] #d_i d_j H
            end
        end
        #we need to define covariance (sigma here ...)
        #sigma is covariance sqrt(qC/(n-p))
        sigma = sqrt(qC/(n-p)) #covariance
        dsigma = sqrt(q/(n-p)) .* 0.5 * Cθ[1]^(-1/2) .* jacC
        hess_sigma = zeros(d, d)
        for i = 1:d
            for j = 1:d
                hess_sigma[i, j] =  sqrt(q/(n-p)) * ( -0.25 * Cθ[1]^(-3/2) .* jacC[i]jacC[j] + 0.5 * Cθ[1]^(-1/2) * hessC[i, j])
            end
        end
        function Y(i, sigma, dsigma)
            return - jacm[i]/sigma .- (gλy0 .- m)/sigma^2 * dsigma[i]
        end
        function R(i, j, sigma, dsigma, hess_sigma) #assume jacm  and hess_m defined already
            e1 = -hess_m[i, j]/sigma
            e2 = jacm[j] * dsigma[i] / sigma^2
            e3 = jacm[i] * dsigma[j] / sigma^2
            e4 = 2*(gλy0 .- m) * dsigma[j]dsigma[i]/sigma^3
            e5 = -(gλy0 .- m)*hess_sigma[i, j] / sigma^2
            return (e1 .+ e2 .+ e3 .+ e4 .+ e5)[1]
        end
        check_arg = zeros(d, d)
        for i = 1:d
            for j = 1:d
                check_arg[i, j] = R(i, j, sigma, dsigma, hess_sigma)
            end
        end
        arg = ((gλy0 .- m)/sqrt(qC/(n-p)))[1] #1 x 1
        vanillat = LocationScale(0, 1, TDist(n-p))
        cdf_eval =  Distributions.cdf.(vanillat, arg)
        cdf_deriv = Distributions.pdf.(vanillat, arg) #chain rule term is computed later
        cdf_second_deriv = -(n-p+k) * arg/( n - p  + arg^2) * cdf_deriv #chain rule term computed later
        gλy0 = g(y0, λ)
        dgλy0 = dg(y0, λ)
        d2gλy0 = dg2(y0, λ)
        hess_G = zeros(1+d, 1+d)
        for i = 1:d  #can make this loop 2x faster by using symmetry
            for j = 1:i
                expr1 = cdf_second_deriv * Y(i, sigma, dsigma) * Y(j, sigma, dsigma)
                expr2 = cdf_deriv * R(i, j, sigma, dsigma, hess_sigma)
                hess_G[i+1, j+1] = hess_G[j+1, i+1] = (expr1 .+ expr2)[1]
            end
        end
        for i = 1:d
            expr1 = cdf_second_deriv * Y(i, sigma, dsigma) * dgλy0 / sigma
            expr2 = cdf_deriv * (-dgλy0) / sigma^2 * dsigma[i]
            hess_G[1, i+1] = hess_G[i+1, 1] = (expr1 + expr2)[1]
        end
        hess_G[1, 1] = cdf_second_deriv * (dgλy0/sigma)^2 + cdf_deriv * d2gλy0 / sigma
        #return hessC
        #return hess_sigma
        #return hessD
        #return hess_m
        #@info "n", n; @info "p", p; @info "gλy0", gλy0
        #return check_arg
        end
        return hess_G
        #jac = x -> abs(reduce(*, map(z -> dg(z, λ), x))) #Jacobian function
        #jacy0 = jac(y0)
        #return (cdf_eval, cdf_deriv* dg(y0,  λ) /sqrt(qC/(n-p)), cdf_second_deriv * (dg(y0, λ)/ sqrt(qC/(n-p)))^2 + cdf_deriv * dg2(y0,λ) / sqrt(qC/(n-p)))
        #return (gλy0, dg(y0,  λ), dg2(y0,λ))
    end

    #function intermediates(x0, Bθ, Hθ, Σθ_inv_X, choleskyXΣX, ΣθinvBθ)
    function intermediates(jac_buffer::jac_buffer, test_data::testingData, train_data::AbstractTrainingData, test_buffer::test_buffer, train_buffer::train_buffer, θλbuffer::θλbuffer)
        update!(jac_buffer, test_data, train_data, test_buffer, train_buffer) #attemot to update jac_buffer
        (jacC, jacB, jacH, jacD) = unpack(jac_buffer)
        Σθ_inv_y = θλbuffer.Σθ_inv_y
        βhat = θλbuffer.βhat
        @timeit to "jacm" jacm = jacB' * Σθ_inv_y + jacH * βhat  # d x 1
        return (jacC, jacm, jacB, jacH, jacD)
    end
    function intermediates_old(x0, Bθ, Hθ, Σθ_inv_X, choleskyXΣX)
        #@assert typeof(θ) <: Array{T, 1} where T
        #@assert length(θ) ==d
        S = repeat(x0, n, 1) .- x # n x d
        d = size(x0, 2)
        #println("size of S: ", size(S))
        flattenedθ = typeof(θ)<:Real ? [θ for i=1:d] : θ[:]
        jacB =   diagm(Bθ[:]) * S * diagm(- flattenedθ) #n x d
        #println("size of jacB: ", size(jacB))
        #println("size of Bθ: ", size(Bθ))
        #println("size of choleskyΣθ: ", size(choleskyΣθ))
        #println("choleskyΣθ inv Bθ': ", size(choleskyΣθ \ Bθ'))
        @timeit to "choleskyΣθ inv Bθ" jacD = -2* jacB' * (choleskyΣθ \ Bθ') #d x 1
        #println("size of jacD: ", size(jacD))
        #assuming linear polynomial basis
        jacFx0 = vcat(zeros(1, d), diagm(ones(d))) #p x d
        #println("size of Fx0: ", size(Fx0))
        jacH = jacFx0' - jacB' * Σθ_inv_X #d x p
        #println("size of jacH: ", size(jacH))
        jacC = jacD + 2 * jacH * (choleskyXΣX \ Hθ') #d x 1
        jacm = jacB' * Σθ_inv_y + jacH * βhat  # d x 1
        #return (jacB, jacD, jacFx0, jacC, jacm)
        return (jacC, jacm, jacB, jacH, jacD)
    end

    # compute quantile q for each component
    function q_fun(x0, Fx0, quant)
        m, q, C = compute_qmC(x0, Fx0)
        t = LocationScale(m, sqrt(q*C/(n-p)), TDist(n-p))
        quant_tdist = Distributions.quantile(t, quant)
        return invg(quant_tdist, λ)
    end
    #return (pdf_deriv, pdf, cdf, cdf_prime_loc, m, Ex2, q_fun)

    return (pdf_deriv, pdf, cdf, cdf_grad_us, cdf_hess_us, q_fun, compute_qmC)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Test
using DataFrames

import Flatten

using Agtor, Dates, CSV

function setup_zone(data_dir::String="test/data/")
    zone_spec_dir::String = "$(data_dir)zones"
    zone_params::Dict{Symbol, Dict} = load_spec(zone_spec_dir)

    return create(zone_params[:Zone_1])
end

function test_in_season()
    dt = Date("2010-03-01")
    s_start = Date("2010-01-01")
    s_end = Date("2010-05-01")

    @test Agtor.in_season(dt, s_start, s_end) == true
end

function test_out_of_season()
    dt = Date("2009-03-01")
    s_start = Date("2010-01-01")
    s_end = Date("2010-05-01")

    @test Agtor.in_season(dt, s_start, s_end) == false

    # Should return false because in_season is start and end exclusive
    dt = Date("2010-01-01")
    @test Agtor.in_season(dt, s_start, s_end) == false
end

function test_crop_plantings()
    z1 = setup_zone()
    field = z1.fields[1]

    for _ in field.crop_rotation
        crop = field.crop
        @test (crop.plant_date + crop.harvest_offset) == crop.harvest_date

        Agtor.set_next_crop!(field)
    end
end

function test_crop_update_stages()
    z1 = setup_zone()
    crop = z1.fields[1].crop

    Agtor.update_stages!(crop, Date("2010-01-01"))

    for (k,v) in pairs(crop.growth_stages)
        @test (v[:start] + Dates.Day(v[:stage_length])) == v[:end]
    end
end

function test_pump_creation(data_dir::String="test/data/")
    pump_spec_dir::String = "$(data_dir)pumps/"
    pump_specs::Dict{Symbol, Dict} = load_spec(pump_spec_dir)

    pumps = Pump[create(spec) for (pn, spec) in pump_specs]

    for created in pumps
        @test created isa Pump
        Flatten.flatten(created)
    end
end

@testset "Crop in-season" begin
    test_in_season()
    test_out_of_season()
end

@testset "Crop plantings and stages" begin
    test_crop_plantings()
    test_crop_update_stages()
end

@testset "Pumping" begin
    test_pump_creation()
end

include("test_parameters.jl")
include("test_scenario_run.jl")

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
abstract TruncatedContinuousUnivariateDistribution <: ContinuousUnivariateDistribution
abstract TruncatedDiscreteUnivariateDistribution <: DiscreteUnivariateDistribution
typealias TruncatedUnivariateDistribution Union(TruncatedContinuousUnivariateDistribution, TruncatedDiscreteUnivariateDistribution)

macro truncate(dname::Any)
    new_dname = esc(symbol(string("Truncated", string(dname))))
    # TODO: Are we not supposed to run eval() in a macro?
    if eval(dname) <: ContinuousUnivariateDistribution
        dtype = esc(TruncatedContinuousUnivariateDistribution)
    else
        dtype = esc(TruncatedDiscreteUnivariateDistribution)
    end
    dname = esc(dname)
    quote
        immutable $new_dname <: $dtype
            untruncated::$dname
            lower::Float64
            upper::Float64
            nc::Float64 # Normalization constant
            function ($new_dname)(d::$dname, l::Real, u::Real, nc::Real)
                if l >= u
                    error("upper must be > lower")
                end
                new(d, float64(l), float64(u), float64(nc))
            end
        end
        function ($new_dname)(d::$dname, l::Real, u::Real)
            return ($new_dname)(d, l, u, cdf(d, u) - cdf(d, l))
        end
    end
end

function insupport(d::TruncatedUnivariateDistribution, x::Number)
    return x >= d.lower && x <= d.upper && insupport(d.untruncated, x)
end

function pdf(d::TruncatedUnivariateDistribution, x::Real)
    if !insupport(d, x)
        return 0.0
    else
        return pdf(d.untruncated, x) / d.nc
    end
end

function logpdf(d::TruncatedUnivariateDistribution, x::Real)
    if !insupport(d, x)
        return -Inf
    else
        return logpdf(d.untruncated, x) - log(d.nc)
    end
end

function cdf(d::TruncatedUnivariateDistribution, x::Real)
    if x < d.lower
        return 0.0
    elseif x > d.upper
        return 1.0
    else
        return (cdf(d.untruncated, x) - cdf(d.untruncated, d.lower)) / d.nc
    end
end

function quantile(d::TruncatedUnivariateDistribution, p::Real)
    top = cdf(d.untruncated, d.upper)
    bottom = cdf(d.untruncated, d.lower)
    return quantile(d.untruncated, bottom + p * (top - bottom))
end

function rand(d::TruncatedUnivariateDistribution)
    while true
        r = rand(d.untruncated)
        if d.lower <= r <= d.upper
            return r
        end
    end
end


text:
import Base: UInt64, convert
using Dates
using CategoricalArrays

# ============================================================================ #
# [Constants]                                                                  #
# ============================================================================ #

const CLIENT_NAME = "ClickHouseJL"
const DBMS_VER_MAJOR = 19
const DBMS_VER_MINOR = 11
const DBMS_VER_REV = 54423

# ============================================================================ #
# [Context structs]                                                            #
# ============================================================================ #

"ClickHouse client socket. Created using `connect`."
mutable struct ClickHouseSock
    io::IO
    server_tz::Union{String, Nothing}
    dirty::Bool
    client_info

    function ClickHouseSock(io::IO)::ClickHouseSock
        new(
            io,
            nothing,
            false,
            ClientInfo(
                0x01,
                "",
                "",
                "0.0.0.0:0",
                0x01,
                "",
                "",
                CLIENT_NAME,
                DBMS_VER_MAJOR,
                DBMS_VER_MINOR,
                DBMS_VER_REV,
                "",
                2,
            ),
        )
    end
end

# ============================================================================ #
# [Variable length integer]                                                    #
# ============================================================================ #

primitive type VarUInt <: Unsigned 64 end

VarUInt(x::Number) = reinterpret(VarUInt, UInt64(x))
UInt64(x::VarUInt) = reinterpret(UInt64, x)
Base.show(io::IO, x::VarUInt) = print(io, UInt64(x))

function chwrite(sock::ClickHouseSock, x::VarUInt)
    mx::UInt64 = x
    while mx >= 0x80
        write(sock.io, UInt8(mx & 0xFF) | 0x80)
        mx >>= 7
    end
    write(sock.io, UInt8(mx & 0xFF))
end

function chread(sock::ClickHouseSock, ::Type{VarUInt})::VarUInt
    x::UInt64 = 0
    s::UInt32 = 0
    i::UInt64 = 0
    while true
        b = read(sock.io, UInt8)
        if b < 0x80
            if i > 9 || (i == 9 && b > 1)
                throw(OverflowError("varint would overflow"))
            end
            return x | UInt64(b) << s
        end

        x |= UInt64(b & 0x7F) << s
        s += 7
        i += 1
    end
end

# ============================================================================ #
# [chread impls for primitive types]                                           #
# ============================================================================ #

# Scalar reads
chread(sock::ClickHouseSock, ::Type{T}) where T <: Number =
    read(sock.io, T)

chread(sock::ClickHouseSock, x::UInt64)::Vector{UInt8} =
    read(sock.io, x)

function chread(sock::ClickHouseSock, ::Type{String})::String
    len = chread(sock, VarUInt) |> UInt64
    chread(sock, len) |> String
end

# Vector reads
function chread(
    sock::ClickHouseSock,
    ::Type{Vector{T}},
    count::VarUInt,
)::Vector{T} where T <: Number
    data = Vector{T}(undef, UInt64(count))
    read!(sock.io, data)
    data
end

chread(
    sock::ClickHouseSock,
    ::Type{Vector{String}},
    count::VarUInt,
)::Vector{String} = [chread(sock, String) for _ ∈ 1:UInt64(count)]

# ============================================================================ #
# [chwrite impls for primitive types]                                          #
# ============================================================================ #

# Scalar writes
chwrite(sock::ClickHouseSock, x::Number) =
    write(sock.io, x)

function chwrite(sock::ClickHouseSock, x::String)
    chwrite(sock, x |> length |> VarUInt)
    chwrite(sock, x |> Array{UInt8})
end

# Vector writes
chwrite(sock::ClickHouseSock, x::AbstractVector{T}) where T <: Number =
    write(sock.io, x)

chwrite(sock::ClickHouseSock, x::AbstractVector{String}) =
    foreach(x -> chwrite(sock, x), x)

# ============================================================================ #
# [Parse helpers]                                                              #
# ============================================================================ #

function chread(sock::ClickHouseSock, ::Type{T}) where {T}
    if @generated
        arg_exprs = [:(chread(sock, $t)) for t ∈ fieldtypes(T)]
        return :($T($(arg_exprs...)))
    else
        return T((chread(sock, t) for t in fieldtypes(T)))
    end
end

function chwrite(sock::ClickHouseSock, x::T) where {T}
    if @generated
        block = Expr(:block)
        for i = 1:fieldcount(T)
            push!(block.args, :(chwrite(sock, getfield(x, $i))))
        end
        return block
    else
        for i = 1:fieldcount(T)
            chwrite(sock, getfield(x, i))
        end
    end
end

# ============================================================================ #
# [Shared messages (wire format)]                                              #
# ============================================================================ #

const BLOCK_INFO_FIELD_STOP = UInt64(0)
const BLOCK_INFO_FIELD_OVERFLOWS = UInt64(1)
const BLOCK_INFO_FIELD_BUCKET_NUM = UInt64(2)

struct BlockInfo
    is_overflows::Bool
    bucket_num::Int32

    BlockInfo() = new(false, -1)
    BlockInfo(is_overflows, bucket_num) = new(is_overflows, bucket_num)
end

function chread(sock::ClickHouseSock, ::Type{BlockInfo})::BlockInfo
    is_overflows = false
    bucket_num = -1

    while (field_num = UInt64(chread(sock, VarUInt))) != BLOCK_INFO_FIELD_STOP
        if field_num == BLOCK_INFO_FIELD_OVERFLOWS
            is_overflows = chread(sock, Bool)
        elseif field_num == BLOCK_INFO_FIELD_BUCKET_NUM
            bucket_num = chread(sock, Int32)
        else
            throw("Unknown block info field")
        end
    end

    BlockInfo(is_overflows, bucket_num)
end

function chwrite(sock::ClickHouseSock, x::BlockInfo)
    # This mirrors what the C++ client does.
    chwrite(sock, VarUInt(BLOCK_INFO_FIELD_OVERFLOWS))
    chwrite(sock, x.is_overflows)
    chwrite(sock, VarUInt(BLOCK_INFO_FIELD_BUCKET_NUM))
    chwrite(sock, x.bucket_num)
    chwrite(sock, VarUInt(BLOCK_INFO_FIELD_STOP))
end

struct Column
    name::String
    type::String
    data::Any
end

const COL_TYPE_MAP = Dict(
    # Unsigned
    "UInt8"    => UInt8,
    "UInt16"   => UInt16,
    "UInt32"   => UInt32,
    "UInt64"   => UInt64,

    # Signed
    "Int8"     => Int8,
    "Int16"    => Int16,
    "Int32"    => Int32,
    "Int64"    => Int64,

    # Floats
    "Float32"  => Float32,
    "Float64"  => Float64,

    "String"   => String,

    "DateTime" => Int32,
    "Date"     => Int16,
)

const COL_TYPE_REV_MAP = Dict(v => k for (k, v) ∈ COL_TYPE_MAP)
const SECS_IN_DAY = 24 * 60 * 60

const ENUM_RE_OUTER = r"Enum(\d{1,2})\(\s*(.*)\)$"
const ENUM_RE_INNER = r"""
    (?:
    '((?:(?:[^'])|(?:\\'))+)'
    \s*=\s*
    (-?\d+)
    \s*,?\s*
    )+?
"""x

function parse_enum_def(str::String)
    def = match(ENUM_RE_OUTER, str)
    matches = eachmatch(ENUM_RE_INNER, def[2])
    map = Dict(x[1] => parse(Int64, x[2]) for x ∈ matches)
    type = "Int" * def[1]
    type, map
end

# We can't just use chread here because we need the size to be passed
# in from the `Block` decoder that holds the row count.
function read_col(sock::ClickHouseSock, num_rows::VarUInt)::Column
    name = chread(sock, String)
    type_name = chread(sock, String)

    decode_type_name, enum_def = if startswith(type_name, "Enum")
        parse_enum_def(type_name)
    else
        type_name, nothing
    end

    decode_type = try
        COL_TYPE_MAP[decode_type_name]
    catch exc
        if exc isa KeyError
            error("Unsupported data type: $(decode_type_name)")
        end
        rethrow()
    end

    data = chread(sock, Vector{decode_type}, num_rows)

    if type_name == "DateTime"
        data = unix2datetime.(data)
    elseif type_name == "Date"
        data = convert(Array{Int64}, data)
        data .*= SECS_IN_DAY
        data = unix2datetime.(data)
    elseif startswith(type_name, "Enum")
        inv_map = Dict(v => k for (k, v) ∈ enum_def)
        data = recode(data, inv_map...)
    end

    Column(name, type_name, data)
end

function chwrite(sock::ClickHouseSock, x::Column)
    chwrite(sock, x.name)
    chwrite(sock, x.type)

    data = if x.type == "DateTime"
        d = datetime2unix.(x.data)
        d = round.(d)
        d = convert(Array{Int32}, d)
    elseif x.type == "Date"
        d = datetime2unix.(x.data)
        d ./= convert(Float64, SECS_IN_DAY)
        d = round.(d)
        d = convert(Array{Int16}, d)
    elseif startswith(x.type, "Enum")
        ty, map = parse_enum_def(x.type)
        # TODO: this could be faster for categorical arrays.
        try
            d = [map[x] for x ∈ x.data]
        catch exc
            if exc isa KeyError
                error("Value is not a valid enum variant: $(exc.key)")
            end
            rethrow()
        end
        ty = COL_TYPE_MAP[ty]
        convert(Array{ty}, d)
    else
        x.data
    end

    chwrite(sock, data)
end

struct Block
    temp_table::String
    block_info::BlockInfo
    num_columns::VarUInt
    num_rows::VarUInt
    columns::Array{Column}
end

function chread(sock::ClickHouseSock, ::Type{Block})::Block
    temp_table = chread(sock, String)
    block_info = chread(sock, BlockInfo)
    num_columns = chread(sock, VarUInt)
    num_rows = chread(sock, VarUInt)
    columns = [read_col(sock, num_rows) for _ ∈ 1:UInt64(num_columns)]
    Block(temp_table, block_info, num_columns, num_rows, columns)
end

function chwrite(sock::ClickHouseSock, x::Block)
    chwrite(sock, x.temp_table)
    chwrite(sock, x.block_info)
    chwrite(sock, x.num_columns)
    chwrite(sock, x.num_rows)
    for x ∈ x.columns
        chwrite(sock, x)
    end
end

# ============================================================================ #
# [Server messages (wire format)]                                              #
# ============================================================================ #

struct ServerInfo
    server_name::String
    server_major_ver::VarUInt
    server_minor_ver::VarUInt
    server_rev::VarUInt

    # DBMS_MIN_REVISION_WITH_SERVER_TIMEZONE
    server_timezone::String

    # DBMS_MIN_REVISION_WITH_SERVER_DISPLAY_NAME
    server_display_name::String

    # DBMS_MIN_REVISION_WITH_VERSION_PATCH
    server_version_patch::VarUInt
end

struct ServerPong
end

struct ServerProgress
    rows::VarUInt
    bytes::VarUInt
    total_rows::VarUInt

    # DBMS_MIN_REVISION_WITH_CLIENT_WRITE_INFO
    written_rows::VarUInt
    written_bytes::VarUInt
end

struct ServerProfileInfo
    rows::VarUInt
    blocks::VarUInt
    bytes::VarUInt
    applied_limit::Bool
    rows_before_limit::VarUInt
    calc_rows_before_limit::Bool
end

struct ServerException
    code::UInt32
    name::String
    message::String
    stack_trace::String
    nested::Union{Nothing, ServerException}
end

function chread(sock::ClickHouseSock, ::Type{ServerException})::ServerException
    code = chread(sock, UInt32)
    name = chread(sock, String)
    message = chread(sock, String)
    stack_trace = chread(sock, String)
    has_nested = chread(sock, Bool)
    nested = has_nested ? chread(sock, ServerException) : nothing
    ServerException(code, name, message, stack_trace, nested)
end

struct ServerEndOfStream
end

struct ServerTableColumns
    external_table_name::String
    columns::String
    sample_block::Block
end

# ============================================================================ #
# [Client messages (wire format)]                                              #
# ============================================================================ #

struct ClientHello
    client_name::String
    client_dbms_ver_major::VarUInt
    client_dbms_ver_minor::VarUInt
    client_dbms_ver_rev::VarUInt
    database::String
    username::String
    password::String
end

struct ClientPing
end

struct ClientInfo
    query_kind::UInt8
    initial_user::String
    initial_query_id::String
    initial_address_string::String
    read_interface::UInt8
    os_user::String
    client_hostname::String
    client_name::String
    client_ver_major::VarUInt
    client_ver_minor::VarUInt
    client_rev::VarUInt
    quota_key::String # DBMS_MIN_REVISION_WITH_QUOTA_KEY_IN_CLIENT_INFO
    client_ver_patch::VarUInt # DBMS_MIN_REVISION_WITH_VERSION_PATCH
end

struct ClientQuery
    query_id::String
    client_info::ClientInfo
    settings::String
    query_stage::VarUInt
    compression::VarUInt
    query::String
end

# ============================================================================ #
# [Opcodes]                                                                    #
# ============================================================================ #

const CLIENT_HELLO = UInt64(0)
const CLIENT_QUERY = UInt64(1)
const CLIENT_DATA = UInt64(2)
const CLIENT_CANCEL = UInt64(3)
const CLIENT_PING = UInt64(4)
const CLIENT_TABLE_STATUS_REQ = UInt64(5)
const CLIENT_KEEP_ALIVE = UInt64(6)

const SERVER_HELLO = UInt64(0)
const SERVER_DATA = UInt64(1)
const SERVER_EXCEPTION = UInt64(2)
const SERVER_PROGRESS = UInt64(3)
const SERVER_PONG = UInt64(4)
const SERVER_END_OF_STREAM = UInt64(5)
const SERVER_PROFILE_INFO = UInt64(6)
const SERVER_TOTALS = UInt64(7)
const SERVER_EXTREMES = UInt64(8)
const SERVER_TABLES_STATUS_REPORT = UInt64(9)
const SERVER_TABLES_LOG = UInt64(10)
const SERVER_TABLE_COLUMNS = UInt64(11)

# ============================================================================ #
# [Message decoding]                                                           #
# ============================================================================ #

const SERVER_OPCODE_TY_MAP = Dict(
    SERVER_HELLO => ServerInfo,
    SERVER_PONG => ServerPong,
    SERVER_PROGRESS => ServerProgress,
    SERVER_PROFILE_INFO => ServerProfileInfo,
    SERVER_EXCEPTION => ServerException,
    SERVER_DATA => Block,
    SERVER_TOTALS => Block,
    SERVER_EXTREMES => Block,
    SERVER_TABLE_COLUMNS => ServerTableColumns,
    SERVER_END_OF_STREAM => ServerEndOfStream,
)

const CLIENT_OPCODE_TY_MAP = Dict(
    CLIENT_HELLO => ClientHello,
    CLIENT_QUERY => ClientQuery,
    CLIENT_DATA => Block,
    CLIENT_PING => ClientPing,
)

function read_packet(
    sock::ClickHouseSock,
    opcode_map::Dict{UInt64, DataType},
)::Any
    opcode = chread(sock, VarUInt)
    ty = opcode_map[UInt64(opcode)]
    chread(sock, ty)
end

"ClickHouse server-side exception."
struct ClickHouseServerException <: Exception
    exc::ServerException
end

read_client_packet(sock::ClickHouseSock)::Any =
    read_packet(sock, CLIENT_OPCODE_TY_MAP)

function read_server_packet(sock::ClickHouseSock)::Any
    packet = read_packet(sock, SERVER_OPCODE_TY_MAP)

    if typeof(packet) == ServerException
        throw(ClickHouseServerException(packet))
    end

    packet
end

# ============================================================================ #
# [Message encoding]                                                           #
# ============================================================================ #

const CLIENT_TY_OPCODE_MAP = Dict(v => k for (k, v) ∈ CLIENT_OPCODE_TY_MAP)

function write_packet(sock::ClickHouseSock, packet::Any)
    opcode = CLIENT_TY_OPCODE_MAP[typeof(packet)]
    chwrite(sock, VarUInt(opcode))
    chwrite(sock, packet)
end

# ============================================================================ #


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "openbabel"
version = v"3.1.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openbabel/openbabel/archive/refs/tags/openbabel-3-1-1.tar.gz", "c97023ac6300d26176c97d4ef39957f06e68848d64f1a04b0b284ccff2744f02")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/openbabel-*
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("obfit", :obfit),
    ExecutableProduct("obspectrophore", :obspectrophore),
    ExecutableProduct("roundtrip", :roundtrip),
    ExecutableProduct("obprop", :obprop),
    ExecutableProduct("obgrep", :obgrep),
    ExecutableProduct("obsym", :obsym),
    ExecutableProduct("obabel", :obabel),
    ExecutableProduct("obthermo", :obthermo),
    LibraryProduct("libopenbabel", :libopenbabel),
    ExecutableProduct("obfitall", :obfitall),
    ExecutableProduct("obprobe", :obprobe),
    ExecutableProduct("obminimize", :obminimize),
    ExecutableProduct("obrotamer", :obrotamer),
    ExecutableProduct("obdistgen", :obdistgen),
    ExecutableProduct("obgen", :obgen),
    ExecutableProduct("obenergy", :obenergy),
    ExecutableProduct("obrms", :obrms),
    ExecutableProduct("obrotate", :obrotate),
    ExecutableProduct("obconformer", :obconformer),
    ExecutableProduct("obmm", :obmm),
    ExecutableProduct("obtautomer", :obtautomer),
    LibraryProduct("libinchi", :libinchi)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="Eigen_jll", uuid="bc6bbf8a-a594-5541-9c57-10b0d0312c70"))
    Dependency(PackageSpec(name="XML2_jll", uuid="02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
_cycle(λ0, λ1, x) = abs(λ0 - λ1) * x + min(λ0, λ1)
_tri(t, period) = (2 / π) * abs(asin(sin(π * (t - 1) / period)))
_sin(t, period) = abs(sin(π * (t - 1) / period))
_cos(t, period) = (1 + cos(π * mod(t - 1, period) / period)) / 2

@testset "Triangle" begin
    λ0 = rand()
    λ1 = rand()
    period = rand(1:10)
    s = Triangle(λ0 = λ0, λ1 = λ1, period = period)
    @test s == Triangle(λ0, λ1, period)
    @test all(_cycle(λ0, λ1, _tri(t, period)) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "TriangleDecay2" begin
    λ0 = rand()
    λ1 = rand()
    period = rand(1:10)
    s = TriangleDecay2(λ0 = λ0, λ1 = λ1, period = period)
    @test s == TriangleDecay2(λ0, λ1, period)
    @test all(_cycle(λ0, λ1, _tri(t, period) / (2^fld(t - 1, period))) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "TriangleExp" begin
    λ0 = rand()
    λ1 = rand()
    γ = rand()
    period = rand(1:10)
    s = TriangleExp(λ0 = λ0, λ1 = λ1, period = period, γ = γ)
    @test s == TriangleExp(λ0, λ1, period, γ)
    @test all(_cycle(λ0, λ1, _tri(t, period) * γ^(t - 1)) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
end

@testset "Sin" begin
    λ0 = rand()
    λ1 = rand()
    period = rand(1:10)
    s = Sin(λ0 = λ0, λ1 = λ1, period = period)
    @test s == Sin(λ0, λ1, period)
    @test all(_cycle(λ0, λ1, _sin(t, period)) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "SinDecay2" begin
    λ0 = rand()
    λ1 = rand()
    period = rand(1:10)
    s = SinDecay2(λ0 = λ0, λ1 = λ1, period = period)
    @test s == SinDecay2(λ0, λ1, period)
    @test all(_cycle(λ0, λ1, _sin(t, period) / (2^fld(t - 1, period))) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "SinExp" begin
    λ0 = rand()
    λ1 = rand()
    γ = rand()
    period = rand(1:10)
    s = SinExp(λ0 = λ0, λ1 = λ1, period = period, γ = γ)
    @test s == SinExp(λ0, λ1, period, γ)
    @test all(_cycle(λ0, λ1, _sin(t, period) * γ^(t - 1)) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "Cos" begin
    λ0 = rand()
    λ1 = rand()
    period = rand(1:10)
    s = Cos(λ0 = λ0, λ1 = λ1, period = period)
    @test s == Cos(λ0, λ1, period)
    @test all(_cycle(λ0, λ1, _cos(t, period)) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ0)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end