text:
function equivalent_to_wkt(geom::GeoInterface.AbstractGeometry, wkt::String)
    test_geom = readgeom(wkt)
    @test LibGEOS.equals(geom, test_geom)
end

function factcheck(f::Function, geom::String, expected::String)
    result = f(readgeom(geom))
    equivalent_to_wkt(result, expected)
end

function factcheck(f::Function, geom::String, expected::Bool)
    @test f(readgeom(geom)) == expected
end

function factcheck(f::Function, g1::String, g2::String, expected::String)
    result = f(readgeom(g1), readgeom(g2))
    equivalent_to_wkt(result, expected)
end

function factcheck(f::Function, g1::String, g2::String, expected::Bool)
    @test f(readgeom(g1), readgeom(g2)) == expected
    @test f(prepareGeom(readgeom(g1)), readgeom(g2)) == expected
end

@testset "GEOS operations" begin
    ls = LineString(Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2]])
    pt = interpolate(ls, 2.5)
    @test GeoInterface.coordinates(pt) ≈ [8.5, 2.0] atol = 1e-5
    for (pt, dist, dest) in [
        (Point(10.0, 1.0), 1.0, Point(9.0, 1.0)),
        (Point(9.0, 1.0), 1.0, Point(9.0, 1.0)),
        (Point(10.0, 0.0), 1.0, Point(9.0, 1.0)),
        (Point(9.0, 2.0), 2.0, Point(9.0, 2.0)),
        (Point(8.7, 1.5), 1.5, Point(9.0, 1.5)),
    ]
        test_dist = project(ls, pt)
        @test test_dist ≈ dist atol = 1e-5
        @test equals(interpolate(ls, test_dist), dest)
    end

    # GEOSConvexHullTest
    input = readgeom(
        "MULTIPOINT (130 240, 130 240, 130 240, 570 240, 570 240, 570 240, 650 240)",
    )
    expected = readgeom("LINESTRING (130 240, 650 240)")
    output = convexhull(input)
    @test !isEmpty(output)
    @test writegeom(output) == writegeom(expected)

    # LibGEOS.delaunayTriangulationTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = delaunayTriangulationEdges(g1)
    @test isEmpty(g1)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :MultiLineString

    g1 = readgeom("POINT(0 0)")
    g2 = delaunayTriangulation(g1)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :GeometryCollection

    g1 = readgeom("MULTIPOINT(0 0, 5 0, 10 0)")
    g2 = delaunayTriangulation(g1, 0.0)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :GeometryCollection
    g2 = delaunayTriangulationEdges(g1, 0.0)
    equivalent_to_wkt(g2, "MULTILINESTRING ((5 0, 10 0), (0 0, 5 0))")

    g1 = readgeom("MULTIPOINT(0 0, 10 0, 10 10, 11 10)")
    g2 = delaunayTriangulationEdges(g1, 2.0)
    equivalent_to_wkt(g2, "MULTILINESTRING ((0 0, 10 10), (0 0, 10 0), (10 0, 10 10))")

    # LibGEOS.constrainedDelaunayTriangulationTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = constrainedDelaunayTriangulation(g1)
    @test isEmpty(g1)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :GeometryCollection

    g1 = readgeom("POINT(0 0)")
    g2 = constrainedDelaunayTriangulation(g1)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :GeometryCollection

    g1 = readgeom("POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))")
    g2 = constrainedDelaunayTriangulation(g1)
    equivalent_to_wkt(g2, "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 90 10, 90 90)))")

    # GEOSDistanceTest
    g1 = readgeom("POINT(10 10)")
    g2 = readgeom("POINT(3 6)")
    @test distance(g1, g2) ≈ 8.06225774829855 atol = 1e-12

    # GEOSGeom_extractUniquePointsTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = uniquePoints(g1)
    @test isEmpty(g2)

    g1 = readgeom("MULTIPOINT(0 0, 0 0, 1 1)")
    g2 = uniquePoints(g1)
    @test equals(g2, readgeom("MULTIPOINT(0 0, 1 1)"))

    g1 = readgeom(
        "GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))",
    )
    @test LibGEOS.equals(uniquePoints(g1), readgeom("MULTIPOINT(0 0, 1 1, 2 2, 5 5, 0 2)"))

    # GEOSGetCentroidTest
    test_centroid(geom::String, expected::String) = factcheck(centroid, geom, expected)
    test_centroid("POINT(10 0)", "POINT (10 0)")
    test_centroid("LINESTRING(0 0, 10 0)", "POINT (5 0)")
    test_centroid("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))", "POINT (5 5)")
    test_centroid("LINESTRING EMPTY", "POINT EMPTY")

    # GEOSIntersectionTest
    test_intersection(g1::String, g2::String, expected::String) =
        factcheck(intersection, g1, g2, expected)
    test_intersection("POLYGON EMPTY", "POLYGON EMPTY", "GEOMETRYCOLLECTION EMPTY")
    test_intersection("POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", "POINT(2 2)")
    test_intersection(
        "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
        "POLYGON((-1 1,-1 2,2 2,2 1,-1 1))",
        "POLYGON ((0 1, 0 2, 2 2, 2 1, 0 1))",
    )
    test_intersection(
        "MULTIPOLYGON(((0 0,5 10,10 0,0 0),(1 1,1 2,2 2,2 1,1 1),(100 100,100 102,102 102,102 100,100 100)))",
        "POLYGON((0 1,0 2,10 2,10 1,0 1))",
        "GEOMETRYCOLLECTION (LINESTRING (1 2, 2 2), LINESTRING (2 1, 1 1), POLYGON ((0.5 1, 1 2, 1 1, 0.5 1)), POLYGON ((9 2, 9.5 1, 2 1, 2 2, 9 2)))",
    )

    # LineString_PointTest
    g1 = readgeom("LINESTRING(0 0, 5 5, 10 10)")
    @test !isClosed(g1)
    @test GeoInterface.geotype(g1) == :LineString
    @test numPoints(g1) == 3
    @test geomLength(g1) ≈ sqrt(100 + 100) atol = 1e-5
    @test GeoInterface.coordinates(startPoint(g1)) ≈ [0, 0] atol = 1e-5
    @test GeoInterface.coordinates(endPoint(g1)) ≈ [10, 10] atol = 1e-5

    # GEOSNearestPointsTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = readgeom("POLYGON EMPTY")
    @test length(nearestPoints(g1, g2)) == 0

    g1 = readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    g2 = readgeom("POLYGON((8 8, 9 9, 9 10, 8 8))")
    points = nearestPoints(g1, g2)
    @test length(points) == 2
    @test GeoInterface.coordinates(points[1])[1:2] == [5.0, 5.0]
    @test GeoInterface.coordinates(points[2])[1:2] == [8.0, 8.0]

    # GEOSNodeTest
    g1 = node(readgeom("LINESTRING(0 0, 10 10, 10 0, 0 10)"))
    LibGEOS.normalize!(g1)
    equivalent_to_wkt(
        g1,
        "MULTILINESTRING ((5 5, 10 0, 10 10, 5 5), (0 10, 5 5), (0 0, 5 5))",
    )

    g1 = node(readgeom("MULTILINESTRING((0 0, 2 0, 4 0),(5 0, 3 0, 1 0))"))
    LibGEOS.normalize!(g1)
    equivalent_to_wkt(
        g1,
        "MULTILINESTRING ((4 0, 5 0), (3 0, 4 0), (2 0, 3 0), (1 0, 2 0), (0 0, 1 0))",
    )

    g1 = node(readgeom("MULTILINESTRING((0 0, 2 0, 4 0),(0 0, 2 0, 4 0))"))
    LibGEOS.normalize!(g1)
    equivalent_to_wkt(g1, "MULTILINESTRING ((2 0, 4 0), (0 0, 2 0))")

    # GEOSPointOnSurfaceTest
    test_pointonsurface(geom::String, expected::String) =
        factcheck(pointOnSurface, geom, expected)
    test_pointonsurface("POINT(10 0)", "POINT (10 0)")
    test_pointonsurface("LINESTRING(0 0, 5 0, 10 0)", "POINT (5 0)")
    test_pointonsurface("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))", "POINT (5 5)")
    test_pointonsurface("LINESTRING EMPTY", "POINT EMPTY")
    test_pointonsurface("LINESTRING(0 0, 0 0)", "POINT (0 0)")

    g1 = readgeom("""POLYGON((
                  56.528666666700 25.2101666667,
                  56.529000000000 25.2105000000,
                  56.528833333300 25.2103333333,
                  56.528666666700 25.2101666667))""")
    @test GeoInterface.coordinates(pointOnSurface(g1)) ≈ [56.5286666667, 25.2101666667] atol =
        1e-5

    # GEOSSharedPathsTest
    factcheck(
        sharedPaths,
        "LINESTRING (-30 -20, 50 60, 50 70, 50 0)",
        "LINESTRING (-29 -20, 50 60, 50 70, 51 0)",
        "GEOMETRYCOLLECTION (MULTILINESTRING ((50 60, 50 70)), MULTILINESTRING EMPTY)",
    )

    # GEOSSimplifyTest
    g1 = readgeom("POLYGON EMPTY")
    @test isEmpty(g1)
    g2 = simplify(g1, 43.2)
    @test isEmpty(g2)
    g1 = readgeom("""POLYGON((
                  56.528666666700 25.2101666667,
                  56.529000000000 25.2105000000,
                  56.528833333300 25.2103333333,
                  56.528666666700 25.2101666667))""")
    equivalent_to_wkt(simplify(g1, 0.0), "POLYGON EMPTY")
    @test equals(g1, topologyPreserveSimplify(g1, 43.2))

    # GEOSSnapTest
    function test_snap(g1::String, g2::String, expected::String, tol::Float64 = 0.0)
        equivalent_to_wkt(snap(readgeom(g1), readgeom(g2), tol), expected)
    end
    test_snap(
        "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))",
        "POINT(0.5 0)",
        "POLYGON ((0.5 0, 10 0, 10 10, 0 10, 0.5 0))",
        1.0,
    )
    test_snap(
        "LINESTRING (-30 -20, 50 60, 50 0)",
        "LINESTRING (-29 -20, 40 60, 51 0)",
        "LINESTRING (-29 -20, 50 60, 51 0)",
        2.0,
    )
    test_snap(
        "LINESTRING (-20 -20, 50 50, 100 100)",
        "LINESTRING (-10 -9, 40 20, 80 79)",
        "LINESTRING (-20 -20, -10 -9, 50 50, 80 79, 100 100)",
        2.0,
    )
    test_snap("LINESTRING(0 0, 10 0)", "LINESTRING(0 0, 9 0)", "LINESTRING(0 0, 9 0)", 2.0)
    test_snap(
        "LINESTRING(0 0, 10 0)",
        "LINESTRING(0 0, 9 0, 10 0, 11 0)",
        "LINESTRING(0 0, 9 0, 10 0, 11 0)",
        2.0,
    )
    test_snap(
        "LINESTRING(0 3,4 1,0 1)",
        "MULTIPOINT(5 0,4 1)",
        "LINESTRING (0 3, 4 1, 5 0, 0 1)",
        2.0,
    )
    test_snap(
        "LINESTRING(0 3,4 1,0 1)",
        "MULTIPOINT(4 1,5 0)",
        "LINESTRING (0 3, 4 1, 5 0, 0 1)",
        2.0,
    )
    test_snap(
        "LINESTRING(0 0,10 0,10 10,0 10,0 0)",
        "MULTIPOINT(0 0,-1 0)",
        "LINESTRING (-1 0, 0 0, 10 0, 10 10, 0 10, -1 0)",
        3.0,
    )
    test_snap(
        "LINESTRING(0 2,5 2,9 2,5 0)",
        "POINT(5 0)",
        "LINESTRING (0 2, 5 2, 9 2, 5 0)",
        3.0,
    )
    test_snap(
        "LINESTRING(-71.1317 42.2511,-71.1317 42.2509)",
        "MULTIPOINT(-71.1261 42.2703,-71.1257 42.2703,-71.1261 42.2702)",
        "LINESTRING (-71.1257 42.2703, -71.1261 42.2703, -71.1261 42.2702, -71.1317 42.2509)",
        0.5,
    )

    # GEOSUnaryUnionTest
    test_unaryunion(geom::String, expected::String) = factcheck(unaryUnion, geom, expected)
    test_unaryunion("POINT EMPTY", "POINT EMPTY")
    test_unaryunion("POINT (6 3)", "POINT (6 3)")
    test_unaryunion("POINT (4 5 6)", "POINT Z (4 5 6)")
    test_unaryunion("MULTIPOINT (4 5, 6 7, 4 5, 6 5, 6 7)", "MULTIPOINT (4 5, 6 5, 6 7)")
    test_unaryunion(
        "GEOMETRYCOLLECTION (POINT(4 5), MULTIPOINT(6 7, 6 5, 6 7), LINESTRING(0 5, 10 5), LINESTRING(4 -10, 4 10))",
        "GEOMETRYCOLLECTION (POINT (6 7), LINESTRING (0 5, 4 5), LINESTRING (4 5, 10 5), LINESTRING (4 -10, 4 5), LINESTRING (4 5, 4 10))",
    )
    test_unaryunion(
        "GEOMETRYCOLLECTION (POINT(4 5), MULTIPOINT(6 7, 6 5, 6 7), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)))",
        "GEOMETRYCOLLECTION (POINT (6 7), POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (5 6, 7 6, 7 8, 5 8, 5 6)))",
    )
    test_unaryunion(
        "GEOMETRYCOLLECTION (MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5)), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)))",
        "GEOMETRYCOLLECTION (LINESTRING (5 7, 7 7), LINESTRING (10 7, 12 7), LINESTRING (5.5 7.5, 6.5 7.5), POLYGON ((10 7, 10 0, 0 0, 0 10, 10 10, 10 7), (5 6, 7 6, 7 7, 7 8, 5 8, 5 7, 5 6)))",
    )
    test_unaryunion(
        "GEOMETRYCOLLECTION (MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5)), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)), MULTIPOINT(6 6.5, 6 1, 12 2, 6 1))",
        "GEOMETRYCOLLECTION (POINT (6 6.5), POINT (12 2), LINESTRING (5 7, 7 7), LINESTRING (10 7, 12 7), LINESTRING (5.5 7.5, 6.5 7.5), POLYGON ((10 7, 10 0, 0 0, 0 10, 10 10, 10 7), (5 6, 7 6, 7 7, 7 8, 5 8, 5 7, 5 6)))",
    )

    for (g1, g2, testvalue) in (
        ("POLYGON EMPTY", "POLYGON EMPTY", false),
        ("POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", false),
        ("POINT(2 2)", "POLYGON((1 1,1 5,5 5,5 1,1 1))", true),
        (
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            false,
        ),
        (
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            true,
        ),
        ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 201, 300 300)", false),
        ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 200, 300 300)", false),
        ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 300 100)", false),
    )
        for f in (within, coveredby)
            factcheck(f, g1, g2, testvalue)
        end
        for f in (contains, covers)
            factcheck(f, g2, g1, testvalue)
        end
    end

    for (g1, g2, testvalue) in (
        ("POLYGON EMPTY", "POLYGON EMPTY", false),
        ("POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", true),
        ("POINT(2 2)", "POLYGON((1 1,1 5,5 5,5 1,1 1))", true),
        (
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            true,
        ),
        (
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            true,
        ),
        ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 201, 300 300)", false),
        ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 200, 300 300)", true),
        ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 300 100)", true),
    )
        factcheck(intersects, g1, g2, testvalue)
        factcheck(disjoint, g1, g2, !testvalue)
    end

    # GEOSisClosedTest
    @test !isClosed(readgeom("LINESTRING(0 0, 1 0, 1 1)"))
    @test isClosed(readgeom("LINESTRING(0 0, 0 1, 1 1, 0 0)"))

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

    @testset "numGeometries" begin
        @test numGeometries(readgeom("POINT(2 2)")) == 1
        @test numGeometries(readgeom("MULTIPOINT(0 0, 5 0, 10 0)")) == 3
        @test numGeometries(readgeom("LINESTRING(0 0, 0 1, 1 1, 0 0)")) == 1
        @test numGeometries(
            readgeom("MULTILINESTRING ((5 5, 10 0, 10 10, 5 5), (0 10, 5 5), (0 0, 5 5))"),
        ) == 3
        @test numGeometries(readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")) == 1
        # Polygon with a hole
        @test numGeometries(
            readgeom(
                "POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10),(20 30, 35 35, 30 20, 20 30))",
            ),
        ) == 1
        @test numGeometries(
            readgeom(
                "MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)), ((15 5, 40 10, 10 20, 5 10, 15 5)))",
            ),
        ) == 2
        # MultiPolygon with holes
        @test numGeometries(
            readgeom(
                "MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),((20 35, 10 30, 10 10, 30 5, 45 20, 20 35),(30 20, 20 15, 20 25, 30 20)))",
            ),
        ) == 2
        @test numGeometries(
            readgeom(
                "GEOMETRYCOLLECTION (POINT (40 10), LINESTRING (10 10, 20 20, 10 40), POLYGON ((40 40, 20 45, 45 30, 40 40)))",
            ),
        ) == 3
    end

    # Buffer should return Polygon or MultiPolygon
    @test buffer(MultiPoint([[1.0, 1.0], [2.0, 2.0], [2.0, 0.0]]), 0.1) isa
          LibGEOS.MultiPolygon
    @test buffer(MultiPoint([[1.0, 1.0], [2.0, 2.0], [2.0, 0.0]]), 10) isa LibGEOS.Polygon

    # bufferWithStyle
    g1 = bufferWithStyle(
        readgeom("LINESTRING(0 0,0 1,1 1)"),
        0.1,
        endCapStyle = LibGEOS.GEOSBUF_CAP_FLAT,
        joinStyle = LibGEOS.GEOSBUF_JOIN_BEVEL,
    )
    g2 = readgeom(
        "POLYGON((-0.1 0.0,-0.1 1.0,0.0 1.1,1.0 1.1,1.0 0.9,0.1 0.9,0.1 0.0,-0.1 0.0))",
    )
    @test equals(g1, g2)

    g1 = bufferWithStyle(
        readgeom("LINESTRING(0 0,0 1,1 1)"),
        0.1,
        endCapStyle = LibGEOS.GEOSBUF_CAP_SQUARE,
        joinStyle = LibGEOS.GEOSBUF_JOIN_MITRE,
    )
    g2 =
        readgeom("POLYGON((-0.1 -0.1,-0.1 1.1,1.1 1.1,1.1 0.9,0.1 0.9,0.1 -0.1,-0.1 -0.1))")
    @test equals(g1, g2)

    g1 = bufferWithStyle(
        readgeom("POLYGON((-1 -1,1 -1,1 1,-1 1,-1 -1))"),
        0.2,
        joinStyle = LibGEOS.GEOSBUF_JOIN_MITRE,
    )
    g2 = readgeom("POLYGON((-1.2 1.2,1.2 1.2,1.2 -1.2,-1.2 -1.2,-1.2 1.2))")
    @test equals(g1, g2)

    @testset "getXMin et al." begin
        # taken from https://github.com/libgeos/geos/blob/main/tests/unit/capi/GEOSGeom_extentTest.cpp
        g = readgeom(("LINESTRING (3 8, -12 -4)"))
        @test LibGEOS.getXMin(g) == -12
        @test LibGEOS.getXMax(g) == 3
        @test LibGEOS.getYMin(g) == -4
        @test LibGEOS.getYMax(g) == 8
        g = readgeom("POLYGON EMPTY")
        @test (@test_throws ErrorException LibGEOS.getXMin(g) == 0).value.msg ==
              "LibGEOS: Error in GEOSGeom_getXMin_r"
        @test (@test_throws ErrorException LibGEOS.getXMax(g) == 0).value.msg ==
              "LibGEOS: Error in GEOSGeom_getXMax_r"
        @test (@test_throws ErrorException LibGEOS.getYMin(g) == 0).value.msg ==
              "LibGEOS: Error in GEOSGeom_getYMin_r"
        @test (@test_throws ErrorException LibGEOS.getYMax(g) == 0).value.msg ==
              "LibGEOS: Error in GEOSGeom_getYMax_r"
    end

    @testset "getGeometry/getGeometries" begin
        @test writegeom(getGeometry(readgeom("POLYGON EMPTY"), 1)) == "POLYGON EMPTY"
        @test writegeom.(getGeometries(readgeom("POLYGON EMPTY"))) == ["POLYGON EMPTY"]

        multipoint = readgeom("MULTIPOINT (4 5, 6 7, 4 5, 6 5, 9 10)")
        @test writegeom(getGeometry(multipoint, 1)) == "POINT (4 5)"
        @test writegeom(getGeometry(multipoint, 2)) == "POINT (6 7)"
        @test writegeom(getGeometry(multipoint, 5)) == "POINT (9 10)"
        @test (@test_throws ErrorException getGeometry(multipoint, 7)).value.msg ==
              "GEOSGetGeometryN: 5 sub-geometries in geom, therefore n should be in 1:5"
        @test writegeom.(getGeometries(multipoint)) ==
              ["POINT (4 5)", "POINT (6 7)", "POINT (4 5)", "POINT (6 5)", "POINT (9 10)"]

        multilinestring =
            readgeom("MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5))")
        @test writegeom(getGeometry(multilinestring, 1)) == "LINESTRING (5 7, 12 7)"
        @test writegeom(getGeometry(multilinestring, 2)) == "LINESTRING (4 5, 6 5)"
        @test writegeom(getGeometry(multilinestring, 3)) == "LINESTRING (5.5 7.5, 6.5 7.5)"
        @test writegeom.(getGeometries(multilinestring)) == [
            "LINESTRING (5 7, 12 7)",
            "LINESTRING (4 5, 6 5)",
            "LINESTRING (5.5 7.5, 6.5 7.5)",
        ]

        multipolygon = readgeom(
            "MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)), ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))",
        )
        @test writegeom(getGeometry(multipolygon, 1)) ==
              "POLYGON ((40 40, 20 45, 45 30, 40 40))"
        @test writegeom(getGeometry(multipolygon, 2)) ==
              "POLYGON ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20))"
        @test writegeom.(getGeometries(multipolygon)) == [
            "POLYGON ((40 40, 20 45, 45 30, 40 40))",
            "POLYGON ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20))",
        ]

        geomcollection = readgeom(
            "GEOMETRYCOLLECTION (POINT (40 10), LINESTRING (10 10, 20 20, 10 40), POLYGON ((40 40, 20 45, 45 30, 40 40)))",
        )
        @test writegeom(getGeometry(geomcollection, 1)) == "POINT (40 10)"
        @test writegeom(getGeometry(geomcollection, 2)) ==
              "LINESTRING (10 10, 20 20, 10 40)"
        @test writegeom(getGeometry(geomcollection, 3)) ==
              "POLYGON ((40 40, 20 45, 45 30, 40 40))"
        @test writegeom.(getGeometries(geomcollection)) == [
            "POINT (40 10)",
            "LINESTRING (10 10, 20 20, 10 40)",
            "POLYGON ((40 40, 20 45, 45 30, 40 40))",
        ]
    end
end


text:
# Various graph layouts.

function layout_random(g::GenericGraph, dim=2)
    # Generate random locations in unit-square.
    #
    # Each node is associated with a dim-dimensional random vector in
    # [0,1]^{dim}.
    #
    # g : graph to embed
    # dim : number of dimensions to use
    #
    # Returns: Dict keyed by vertices with positions as values.
    #
    # Ported from networkx.

    n = num_vertices(g)
    pos = rand(n, dim)

    return Dict(g.vertices, pos)
end

function layout_circular(g::GenericGraph, dim=2)
    error("Circular layout not implemented yet!")
end

function layout_shell(g::GenericGraph, dim=2)
    error("Shell layout not implemented yet!")
end

function layout_fruchterman_reingold(g::GenericGraph, dim=2)
    error("Fruchterman-Reingold layout not implemented yet!")
end

function layout_kamada_kawai(g::GenericGraph, dim=2)
    error("Kamada-Kawai layout not implemented yet!")
end

function _spectral_dense{T<:Number}(L::Matrix{T}, dim=2)
    # Helper function for spectral embedding
end

function _spectral_sparse{T<:Real}(L::CSCMatrix{T}, dim=2)
    # Helper function for spectral embedding
end

function layout_spectral(g::GenericGraph, dim=2, scale=1)
    # Use eigenvectors of graph Laplacian as coordinates of
    # nodes.
    #
    # g : graph to embed
    # dim : number of dimensions to use
    #
    # Returns: Dict keyed by vertices with positions as values.
    #
    # Notes:
    #   * Sparse eigen-solver used for large graphs.
    #
    # Based on networkx version.

    error("Spectral not done yet")

end

function layout_spectral{W}(g::GenericGraph, eweights::AbstractVector{W},
                         dim=2, scale=1)
    # Use eigenvectors of weighted graph Laplacian as coordinates of
    # nodes.
    #
    # g : graph to embed
    # dim : number of dimensions to use
    #
    # Returns: Dict keyed by vertices with positions as values.
    #
    # Notes:
    #   * Sparse eigen-solver used for large graphs.
    #
    # Based on networkx version.

    error("Weighted spectral not done yet")

end

ext: jl


text:
#<----Bayes Theorem---->
#P = prob, A = prior, B = Evidence,
"""
    ## Bayes Theorem
    ### Description
      Bayes theorem takes probability, prior observations, and evidence
      as percentages and will return a posterior.\n
      --------------------\n
    ### Input
      bay_ther(p, a, b)\n
      --------------------\n
      #### Positional Arguments
      Float64 - p:: A percentage representing probability.\n
      Float64 - a:: A percentage representing prior. <- Float64\n
      Float64 - b:: A percentage representing evidence. <- Float64\n
      --------------------\n
     ### Output
     P:: Posterior value.
       """
function bay_ther(p,a,b)
    psterior = (p*(b|a) * p*(a)) / (p*b)
    return(psterior)
end
"""
    ## Conditional Probability
    ### Description
      Performs Bayesian Conditional Probability with probability, prior, and
      evidence.\n
      --------------------\n
      ### Input
        cond_prob(p, a, b)\n
        --------------------\n
        #### Positional Arguments
        Float64 - p:: A percentage representing probability.\n
        Float64 - a:: A percentage representing prior. <- Float64\n
        Float64 - b:: A percentage representing evidence. <- Float64\n
        --------------------\n
       ### Output
       P:: Posterior value.
       """
function cond_prob(p,a,b)
    psterior = bay_ther(p,a,b)
    cond = p*(a|b)
    return(cond)
end


text:
function Taylor_Chebyshev_I(n::Integer, R0::Number)::Function
    # Antenna Theory 7-26
    A_(R0) = 1 / π * acosh(R0)
    # Antenna Theory 7-28
    σ_(n, A) = n / √(A^2 + (n - 0.5)^2)
    # Antenna Theory 7-29
    μ_(n, σ, A) = π * σ * √(A^2 + (n - 0.5)^2)

    A = A_(R0)
    σ = σ_(n, A)
    μ(n) = μ_(n, σ, A)

    # Antenna Theory 7-30a
    SF(p) = begin
        item1 = factorial(n - 1)^2 / (factorial(n - 1 + p)factorial(n - 1 - p))
        item2 = prod([1 - (p * π / μ(i))^2 for i in 1:n-1])
        item1 * item2
    end
    # Antenna Theory 7-30
    I(l, z) = begin
        item1 = 1 / l
        item2 = sum(SF(p)cos(2π * p * z / l) for p in 1:n-1)
        item1 * (1 + 2 * item2)
    end
end


text:
function destructure(m)
  xs = []
  Flux.mapleaves(m) do x
    x isa AbstractArray && push!(xs, x)
    return x
  end
  return vcat(vec.(xs)...)
end

function restructure(m, xs)
  i = 0
  Flux.mapleaves(m) do x
    x isa AbstractArray || return x
    x = reshape(xs[i.+(1:length(x))], size(x))
    i += length(x)
    return x
  end
end


text:
# ---
# title: 775. Global and Local Inversions
# id: problem775
# author: Indigo
# date: 2021-01-28
# difficulty: Medium
# categories: Array, Math
# link: <https://leetcode.com/problems/global-and-local-inversions/description/>
# hidden: true
# ---
#
# We have some permutation `A` of `[0, 1, ..., N - 1]`, where `N` is the length
# of `A`.
#
# The number of (global) inversions is the number of `i < j` with `0 <= i < j <
# N` and `A[i] > A[j]`.
#
# The number of local inversions is the number of `i` with `0 <= i < N` and
# `A[i] > A[i+1]`.
#
# Return `true` if and only if the number of global inversions is equal to the
# number of local inversions.
#
# **Example 1:**
#
#
#
#     Input: A = [1,0,2]
#     Output: true
#     Explanation: There is 1 global inversion, and 1 local inversion.
#
#
# **Example 2:**
#
#
#
#     Input: A = [1,2,0]
#     Output: false
#     Explanation: There are 2 global inversions, and 1 local inversion.
#
#
# **Note:**
#
#   * `A` will be a permutation of `[0, 1, ..., A.length - 1]`.
#   * `A` will have length in range `[1, 5000]`.
#   * The time limit for this problem has been reduced.
#
#
## @lc code=start
using LeetCode

function is_ideal_permutation(A::Vector{Int})
    m = A[end]
    for i in length(A):-1:3
        m = min(m, A[i])
        if A[i - 2] > m
            return false
        end
    end
    return true
end
## @lc code=end


text:
module AdventOfCode21

include("./day1.jl")
include("./day2.jl")
include("./day3.jl")
include("./day4.jl")
include("./day5.jl")
include("./day6.jl")
include("./day7.jl")
include("./day8.jl")

function run(data_dir::String)
    run_bit("Day 1 (part 1)", Day1.run_part1, data_dir)
    run_bit("Day 1 (part 2)", Day1.run_part2, data_dir)
    run_bit("Day 2 (part 1)", Day2.run_part1, data_dir)
    run_bit("Day 2 (part 2)", Day2.run_part2, data_dir)
    run_bit("Day 3 (part 1)", Day3.run_part1, data_dir)
    run_bit("Day 3 (part 2)", Day3.run_part2, data_dir)
    run_bit("Day 4 (part 1)", Day4.run_part1, data_dir)
    run_bit("Day 4 (part 2)", Day4.run_part2, data_dir)
    run_bit("Day 5 (part 1)", Day5.run_part1, data_dir)
    run_bit("Day 5 (part 2)", Day5.run_part2, data_dir)
    run_bit("Day 6 (part 1)", Day6.run_part1, data_dir)
    run_bit("Day 6 (part 2)", Day6.run_part2, data_dir)
    run_bit("Day 7 (part 1)", Day7.run_part1, data_dir)
    run_bit("Day 7 (part 2)", Day7.run_part2, data_dir)
    run_bit("Day 8 (part 1)", Day8.run_part1, data_dir)
    run_bit("Day 8 (part 2)", Day8.run_part2, data_dir)
end

function run_bit(name::String, func::Function, data_dir::String)
    printstyled("$name:", bold=true, color=:green)
    println("\t\t$(func(data_dir))")
end

end # module


text:
using BenchmarkTools

using KernelAbstractions: @kernel, @index

using Oceananigans
using Oceananigans.Operators: ∇²ᶜᶜᶜ
using Oceananigans.Architectures: device
using Oceananigans.BoundaryConditions: fill_halo_regions!

@info "Using $(Base.Threads.nthreads()) threads"

@kernel function _∇²_KA!(∇²f, grid, f)
    i, j, k = @index(Global, NTuple)
    @inbounds ∇²f[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, f)
end

function _∇²_base_threads!(∇²f, grid, f)

    Nx, Ny, Nz = size(grid)

    Base.Threads.@threads for k = 1:Nz
        for j = 1:Ny
            for i = 1:Nx
                @inbounds ∇²f[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, f)
            end
        end
    end

    return nothing
end

function ∇²_KA!(∇²ϕ, ϕ)
    arch = ϕ.architecture
    grid = ϕ.grid

    Nx, Ny, Nz = worksize = size(grid)
    workgroup = (16, 16)

    fill_halo_regions!(ϕ)
    loop! = _∇²_KA!(device(arch), workgroup, worksize)
    event = loop!(∇²ϕ, grid, ϕ)
    wait(device(arch), event)

    return nothing
end

function ∇²_base_threads!(∇²ϕ, ϕ)
    arch = ϕ.architecture
    fill_halo_regions!(ϕ)
    _∇²_base_threads!(∇²ϕ, grid, ϕ)
    return nothing
end

grid = RectilinearGrid(size=(64, 64, 64), extent=(2π, 2π, 2π))

ϕ = CenterField(CPU(), grid)
∇²ϕ = CenterField(CPU(), grid)

set!(ϕ, (x, y, z) -> randn())

@btime ∇²_KA!(∇²ϕ, ϕ)

@btime ∇²_base_threads!(∇²ϕ, ϕ)


text:
module TestPipeline

using Test
using AutoMLPipeline
using AutoMLPipeline.Pipelines
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.SKPreprocessors
using AutoMLPipeline.DecisionTreeLearners
using AutoMLPipeline.Utils
using AutoMLPipeline.FeatureSelectors

function test_pipeline()
  data = getiris()
  X=data[:,1:5]
  Y=data[:,5] |> Vector
  X[!,5]= X[!,5] .|> string
  ohe = OneHotEncoder()
  ohe1 = OneHotEncoder()
  linear1 = LinearPipeline(Dict(:name=>"lp",:machines => [ohe]))
  linear2 = LinearPipeline(Dict(:name=>"lp",:machines => [ohe]))
  combo1 = ComboPipeline(Dict(:machines=>[ohe,ohe]))
  combo2 = ComboPipeline(Dict(:machines=>[linear1,linear2]))
  linear1 = LinearPipeline([ohe])
  linear2 = LinearPipeline([ohe])
  combo1 = ComboPipeline([ohe,ohe])
  combo2 = ComboPipeline([linear1,linear2])
  fit!(combo1,X)
  res1=transform!(combo1,X)
  res2=fit_transform!(combo1,X)
  @test (res1 .== res2) |> Matrix |> sum == 2100
  fit!(combo2,X)
  res3=transform!(combo2,X)
  res4=fit_transform!(combo2,X)
  @test (res3 .== res3) |> Matrix |> sum == 2100
  pcombo1 = @pipeline ohe1 * ohe1
  pres1 = fit_transform!(pcombo1,X)
  @test (pres1 .== res1) |> Matrix |> sum == 2100
  features = data[:,1:4]
  pca = SKPreprocessor("PCA")
  ica = SKPreprocessor("FastICA")
  fa = SKPreprocessor("FactorAnalysis")
  pcombo2 = @pipeline (pca + ica)*ica*pca
  @test fit_transform!(pcombo2,features) |> Matrix |> size |> collect |> sum == 162
  pcombo2 = @pipeline pca+ica+fa
  @test fit_transform!(pcombo2,features) |> Matrix |> size |> collect |> sum == 154
  disc = CatNumDiscriminator()
  catf = CatFeatureSelector()
  numf = NumFeatureSelector()
  rf = RandomForest()
  pcombo3 = @pipeline disc + ((catf*numf) * (numf+pca) * (numf+ica) * (catf+ohe)) + rf
  (fit_transform!(pcombo3,X,Y)  .== Y) |> sum == 150
end
@testset "Pipelines" begin
    test_pipeline()
end

end


text:
function call(a)
    for i in 1:a
        println("Hello")
    end
end


text:
export FGraph, legmap
export dangling_legs, occupied_legs, rem_leg!
export eliminate
export line_graph

mutable struct FGraph{T, TT<:AbstractArray{T}, ET, VTT<:Union{NTuple{<:Any, TT}, Vector{TT}}} <: IncidenceMatrix
    tensors::VTT
    legmap::Matrix{ET}
end

legmap(fg::FGraph) = fg.legmap
Base.eltype(x::FGraph{T}) where T = T
Base.copy(fg::FGraph) = FGraph(fg.tensors[:], fg.legmap)

################# IncidenceMatrix interfaces #################
imap(fg::FGraph) = legmap(fg)

# remove edge i
function rem_edge!(fg::FGraph, ie)
    fg.legmap = deleteat(fg.legmap, 2, ie)
    fg
end

function rem_vertex!(fg::FGraph, it)
    deleteat!(fg.tensors, it)
    fg.legmap = deleteat(fg.legmap, 1, it)
    fg
end

function add_vertex!(fg::FGraph{T, TT}, ts::TT, edges) where {T, TT<:AbstractArray{T}}
    lm = similar(fg |> legmap, nv(fg)+1, ne(fg))
    @inbounds copy!(view(lm, 1:nv(fg), :), fg |> legmap)
    lm[end, :] .= 0
    for (i, ie) in enumerate(edges)
        @inbounds ie isa Nothing || (lm[end, ie] = i)
    end
    fg.legmap = lm
    push!(fg.tensors, ts)
    fg
end

### NEW INTERFACE
"""remove a leg"""
function rem_leg!(fg::FGraph, iv::Int, idim::Int)
    @inbounds for ie in 1:ne(fg)
        if fg.legmap[iv, ie] == idim
            fg.legmap[iv, ie] = 0
        elseif fg.legmap[iv,ie] > idim
            fg.legmap[iv, ie] -= 1
        end
    end
    fg
end

function occupied_legs(fg::FGraph, it::Int)
    legmap(fg)[it, edges(fg, it)]
end

function dangling_legs(fg::FGraph, it::Int)
    setdiff(1:ndims(fg.tensors[it]), occupied_legs(fg, it))
end

for FUNC in [:dangling_legs, :occupied_legs]
    @eval $FUNC(fg::FGraph) = $FUNC.(Ref(fg), fg |> vertices)
end

"""
    eliminate([fg::FGraph], ie::Int) -> FGraph

eliminate a variable.
"""
function eliminate(fg::FGraph, ie::Int)
    Ne = ne(fg)
    vs = ie isa Tuple ? union(vertices.(Ref(fg), ie)...) : vertices(fg, ie)
    Nvs = length(vs)
    Nvs == 0 && return fg |> rem_edge(ie)
    Nvs == 1 && return _eliminate_loop(fg, ie, vs[])
    #Nvs != 2 && throw(NotSimpleGraphError("Star contraction is not implemented yet: $fg"))
    IVS = edges.(Ref(fg), vs)
    #ICON = intersect(IV...)
    IALL = union(IVS...)
    IC = setdiff(IALL, ie)
    ID = []
    #IDS = dangling_legs.(Ref(fg), vs)
    k = Ne
    for i = 1:Nvs
        IV, v = IVS[i], vs[i]
        Nd = ndims(fg.tensors[v])
        _code = zeros(Int, Nd)
        for ie in IV
            @inbounds _code[legmap(fg)[v,ie]] = ie
        end
         @inbounds for idim = 1:Nd
            if _code[idim] == 0
                k += 1
               _code[idim] = k
               push!(ID, k)
            end
        end
        IVS[i] = _code
    end

    #new_tensor = eincontract((i%2==1 ? fg.tensors[vs[(i+1)÷2]] : Tuple(IVS[(i+1)÷2]) for i=1:2*length(IVS))..., union(IC, ID)|>Tuple)
    xs = Tuple(fg.tensors[vs])
    ixs = Tuple(Tuple.(IVS))
    iy = Tuple(union(IC, ID))
    new_tensor = einsum(EinCode{ixs,iy}(), xs)
    NIC = Tuple(i > ie ? i-1 : i for i in IC)
    fg |> rem_edge(ie) |> rem_vertex(vs) |> add_vertex(new_tensor, NIC)
end
eliminate(args...) = fg::FGraph -> eliminate(fg, args...)

function _eliminate_loop(fg::FGraph, ie::Int, iv::Int)
    il = legmap(fg)[iv, ie]
    fg = fg |> rem_edge(ie)
    fg.tensors[iv] = dropdims(sum(fg.tensors[iv], dims=il), dims=il)
    rem_leg!(fg, iv, il)
end

# check tensor shapes, leg indices
function check_tensors(fg::FGraph)
end

# TODO line graph
function line_graph(fg::FGraph, edge_tensors)
    FGraph(edge_tensors, fg.legmap')
end


text:
# Aiyagari model
# Endogenous Grid Points with IID Income
# Greg Kaplan 2017
# Translated by Tom Sweeney Jan 2021

using Random, Interpolations, NLsolve, Plots
include("discrete_normal.jl")
include("lininterp1.jl")

# PARAMETERS

## preferences
risk_aver = 2
beta = 0.95

## production
deprec = 0.10
capshare = 0.4

## income risk: discretized N(mu,sigma^2)
mu_y = 1
sd_y = 0.2
ny = 5

## asset grids
na = 40
amax = 50
borrow_lim = 0
agrid_par = 0.5 # 1 for linear, 0 for L-shaped

## computation
max_iter = 1000
tol_iter = 1.0e-6
Nsim = 50000
Tsim = 500

maxiter_KL = 70
tol_KL = 1.0e-5
step_KL = 0.005
rguess = 1/beta-1-0.001 # a bit lower than inverse of discount rate
KLratioguess = ((rguess + deprec)/capshare)^(1/(capshare-1))

# OPTIONS
Display = 1
MakePlots = 1

## which function to interpolation
InterpCon = 0
InterpEMUC = 1

## tolerance for non-linear solver
TolX=1.0e-6

# UTILITY FUNCTION

if risk_aver==1
    u(c) = log.(c)
else
    u(c) = (c.^(1-risk_aver).-1)./(1-risk_aver)
end

u1(c) = c.^(-risk_aver)
u1inv(u) = u.^(-1 ./risk_aver)

# SET UP GRIDS

## assets
agrid = range(0,1,length=na)
agrid = agrid.^(1 ./ agrid_par)
agrid = borrow_lim .+ (amax.-borrow_lim).*agrid

## income: disretize normal distribution
width = nlsolve(x -> discrete_normal(ny,mu_y,sd_y,x...)[1],[2.0]).zero
temp, ygrid, ydist = discrete_normal(ny,mu_y,sd_y,width...)
ycumdist = cumsum(ydist)

# DRAW RANDOM NUMBERS
Random.seed!(2021)
yrand = rand(Nsim,Tsim)

# SIMULATE LABOR EFFICIENCY REALIZATIONS
if Display>=1
    println("Simulating labor efficiency realizations in advance")
end
yindsim = zeros(Int,Nsim,Tsim)

for it = 1:Tsim

    # income realization: note we vectorize simulations at once because
    # of matlab, in other languages we would loop over individuals
    yindsim[yrand[:,it].<=ycumdist[1],it] .= 1
    for iy = 2:ny
        yindsim[(yrand[:,it].>ycumdist[iy-1]) .& (yrand[:,it].<=ycumdist[iy]),it] .= iy;
    end
end

ysim = ygrid[yindsim]

# ITERATE OVER KL RATIO
KLratio = KLratioguess

iterKL = 0
KLdiff = 1

while iterKL<=maxiter_KL && abs(KLdiff)>tol_KL
    iterKL = iterKL + 1

    r = capshare.*KLratio^(capshare-1) - deprec
    R = 1+r
    wage = (1-capshare).* KLratio^capshare

    ## rescale efficiency units of labor so that output = 1
    yscale = (KLratio^(-capshare))./(ygrid'*ydist)

    ## initialize consumption function in first iteration only
    if iterKL==1
        conguess = zeros(na,ny)
        for iy = 1:ny
            conguess[:,iy] = r.*agrid .+ wage.* yscale.*ygrid[iy]
        end
        global con = copy(conguess)
     end

    ## solve for policy functions with EGP

    iter = 0
    cdiff = 1000
    while iter <= max_iter && cdiff>tol_iter
        iter = iter + 1
        global sav = zeros(na,ny)

        conlast = copy(con)

        emuc = u1(conlast)*ydist
        muc1 = beta.*R.*emuc
        con1 = u1inv(muc1)

        ## loop over income
        ass1 = zeros(na,ny)
        for iy = 1:ny

            ass1[:,iy] = (con1 .+ agrid .-wage.* yscale.*ygrid[iy])./R

            ## loop over current period ssets
            for ia  = 1:na
                if agrid[ia]<ass1[1,iy] # borrowing constraint binds
                    sav[ia,iy] = borrow_lim
                else # borrowing constraint does not bind;
                    sav[ia,iy] = lininterp1(ass1[:,iy],agrid,agrid[ia])
                end
            end
            con[:,iy] = R.*agrid .+ wage.*yscale.*ygrid[iy] - sav[:,iy]
        end

        cdiff = maximum(abs.(con - conlast))
        if Display>=2
            println("Iteration no. " * string(iter), " max con fn diff is " * string(cdiff))
        end
    end

    ## simulate: start at assets from last interation
    if iterKL==1
        global asim = zeros(Nsim,Tsim)
    elseif iterKL>1
        # asim[:,1] = Ea.*ones(Nsim,1)
        asim[:,1] = asim[:,Tsim]
    end

    ## create interpolating function
    savinterp = Array{Any}(undef,ny)
    for iy = 1:ny
        savinterp[iy] = interpolate((agrid,), sav[:,iy], Gridded(Linear()))
    end

    ## loop over time periods
    for it = 1:Tsim
        if Display>=2 && mod(it,100)==0
            println("Simulating, time period " * string(it))
        end

        ## asset choice
        if it<Tsim
            for iy = 1:ny
                asim[yindsim[:,it].==iy,it+1] = savinterp[iy](asim[yindsim[:,it].==iy,it])
            end
        end
    end

    ## assign actual labor income values
    labincsim = wage.*yscale.*ysim

    ## mean assets and efficiency units
    global Ea = mean(asim[:,Tsim])
    L = yscale.*mean(ysim[:,Tsim])

    KLrationew = Ea./ L

    KLdiff = KLrationew./KLratio - 1
    if Display>=1
        println("Equm iter " * string(iterKL) * ", r = " * string(r) * ", KL ratio: " * string(KLrationew) * " KL diff: " * string(KLdiff*100) * "%")
    end

    KLratio = (1-step_KL)*KLratio + step_KL*KLrationew
end

# MAKE PLOTS
if MakePlots==1

    ## consumption policy function
    p1 = plot(agrid, [con[:,1] con[:,ny]], xlims=(0,amax), title="Consumption Policy Function", color=[:blue :red], label=["Lowest income state" "Highest income state"])
    display(p1)

    ## savings policy function
    p2 = plot(agrid, [sav[:,1].-agrid[:,1] sav[:,ny].-agrid[:,1]], xlims=(0,amax), title="Savings Policy Function", color=[:blue :red], legend=false)
    plot!(agrid, zeros(na,1), color=:black, lw=0.5)
    display(p2)

    ## nice zoom
    xlimits = (0,1)
    xlimind = trues(na)
    if minimum(agrid) < xlimits[1]
        xlimind = xlimind .& (agrid.>=maximum(agrid[agrid<xlimits[1]]))
    elseif minimum(agrid) > xlimits[2]
        xlimind .= 0
    end
    if maximum(agrid) > xlimits[2]
        xlimind = xlimind .& (agrid.<=minimum(agrid[agrid.>xlimits[2]]))
    elseif maximum(agrid) < xlimits[1]
        xlimind .= 0
    end

    ## consumption policy function: zoomed in
    p3 = plot(agrid[xlimind], [con[xlimind,1] con[xlimind,ny]], xlims=xlimits, title="Consumption: Zoomed", marker=:circle, color=[:blue :red], linewidth=2, legend=false)
    plot!(show=true)
    display(p3)

    ## savings policy function: zoomed in
    p4 = plot(agrid[xlimind], [sav[xlimind,1].-agrid[xlimind] sav[xlimind,ny].-agrid[xlimind]], xlims=xlimits, title="Savings: Zoomed (a'-a)", marker=:circle, color=[:blue :red], linewidth=2, legend=false)
    plot!(agrid, zeros(na,1), color=:black, lw=0.5)
    display(p4)

    ## income distribution
    p5 = histogram(ysim[:,Tsim], bins=[2*ygrid[1]-ygrid[2];ygrid].+(ygrid[2]-ygrid[1])/2, title="Income distribution", color=RGB(0,0.5,0.5), linecolor=:blue, legend=false)
    display(p5)

    ## asset distribution
    p6 = histogram(asim[:,Tsim], nbins=100, title="Asset distribution", color=RGB(.7,.7,.7), linecolor=:black, legend=false)
    display(p6)

    ## convergence check
    p7 = plot(1:Tsim, mean(asim,dims=1)', ylims=(0,2*Ea), title="Mean Asset Convergence", xlabel="Time Period", color=:black, lw=1.5, legend=false)
    display(p7)

    ## asset distribution statistics
    aysim = asim[:,Tsim]./mean(ysim[:,Tsim])
    println("Mean assets: " * string(mean(aysim)))
    println("Fraction borrowing constrained: " * string(sum(aysim.==borrow_lim)/Nsim * 100) * '%')
    println("10th Percentile: " * string(quantile(aysim,.1)))
    println("50th Percentile: " * string(quantile(aysim,.5)))
    println("90th Percentile: " * string(quantile(aysim,.9)))
    println("99th Percentile: " * string(quantile(aysim,.99)))
end


text:
using Clang

find_llvm_config() = begin
    if haskey(ENV, "LLVM_CONFIG")
        return readchomp(`which $(LLVM_CONFIG)`)
    end
    tmp = joinpath(JULIA_HOME, "llvm-config")
    isfile(tmp) && return tmp
    tmp = readchomp(`which llvm-config`)
    isfile(tmp) && return tmp
end

if !haskey(ENV, "LLVM_CONFIG") ENV["LLVM_CONFIG"] = find_llvm_config() end

clang_include_path = readchomp(`$(ENV["LLVM_CONFIG"]) --includedir`)
ncurses_include_path = joinpath(Pkg.dir("Ncurses"), "deps", "usr", "include")

indexh = joinpath(clang_include_path, "clang-c/Index.h")

_clang_includes = ASCIIString[]
push!(_clang_includes, "/usr/local/include")
push!(_clang_includes, clang_include_path)
push!(_clang_includes, dirname(indexh))
push!(_clang_includes, ncurses_include_path)
# joinpath("/usr/local", "include"), joinpath(clang_include_path), dirname(indexh), ncurses_include_path]

headers = ["ncurses.h"]
# headers = ["ncurses.h", "form.h", "menu.h", "panel.h"]
headers = [joinpath(ncurses_include_path, "ncursestw", h) for h in headers]

out_path = joinpath(Pkg.dir("Ncurses"), "test")
out_file = joinpath(out_path, "gen_ncurses.jl")
out_common = joinpath(out_path, "gen_ncurses_h.jl")

# clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS", "-D", "NCURSES_WIDECHAR", "-D", "NCURSES_TPARM_VARARGS", "-D", "NCURSES_INTEROP_FUNCS", "-D", "NCURSES_REENTRANT", "-D", "NCURSES_OPAQUE"]

context=wrap_c.init(
                    headers = headers,
                    output_file="gen_ncurses.jl",
                    common_file="gen_ncurses_h.jl",
                    header_library="ncurses",
                    clang_includes=_clang_includes,
                    clang_diagnostics=true
                    )

context.options = wrap_c.InternalOptions(true,true)  # wrap structs, too
context.options.wrap_structs = true

# Execute the wrap
run(context)


text:
function ec_scatter_init()
    opt = ec_init_cartesian()
    opt["tooltip"]["trigger"] = "item"
    opt["toolbox"]["feature"]["dataZoom"] = true
    opt["yAxis"] = [
        Dict("type" => "value",
        "splitArea" => Dict("show" => true),
         "scale" => true)
        ]
    opt["xAxis"] = [
        Dict("type" => "value",
        "scale" => true)
        ]

    return opt
end

function ec_scatter{S<:String}(legendlist::Vector{S}, datas::Array...)
    opt = ec_scatter_init()
    opt["legend"]["data"] = legendlist

    m = length(datas)
    resize!(opt["series"], m)

    for i = 1:m
        opt["series"][i] = Dict(
            "name" => legendlist[i],
            "type" => "scatter",
            "data" => datas[i].'
            )
    end

    return opt
end

function ec_scatter_large{S<:String}(legendlist::Vector{S}, datas::Array...)
    opt = ec_scatter(legendlist, datas...)

    for i = 1:length(opt["series"])
    	opt["series"][i]["large"] = true
    end

    return opt
end


text:
using Fleetdm
using Test

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


text:
# Wrapper of htsFile with customized backend
# ------------------------------------------

## ? It seems that htsindex does not have to associate with htsFile?
## ? May consider removing sethtsindex! and use a new struct to manage index's lifetime.

export SamHeader
export HTSReadWriter

function remote_io(source::AbstractString; write=false, gcsclient = nothing)::IO
    if startswith(source, "gs://")
        if isnothing(gcsclient)
            gcsclient = HttpIO.GCSClient()
        end
        write == true && throw(ArgumentError("GCSFileIO does not support write"))
        return GCSFileIO(source, gcsclient)
    end
    if startswith(source, "http://") || startswith(source, "https://")
        write == true && throw(ArgumentError("HttpFileIO does not support write"))
        return HttpFileIO(source)
    end
    # assume local file
    mode = (write == true) ? "r+" : "r"
    return open(source, mode)
end

## SamHeader should be treated like an immutable object
mutable struct SamHeader
    sam_hdr::Ptr{htslib.sam_hdr_t}
    # below are cached fields
    # TODO: remove these cached fields?
    text::String
    lines::Vector{SubString}
    function SamHeader(ptr::Ptr)
        ptr == C_NULL && error("pointer is invalid")
        ans = new(ptr)
        Base.finalizer(ans) do ans
            ans.sam_hdr == C_NULL && error("pointer is invalid")
            htslib.sam_hdr_destroy(ans.sam_hdr)
        end
        ans.text = unsafe_string(htslib.sam_hdr_str(ans.sam_hdr))
        ans.lines = collect(filter(s -> s != "", split(ans.text, "\n")))
        ans
    end
end

function Base.pointer(x::SamHeader)
    x.sam_hdr
end

function SamHeader(x::SamHeader)
    dup = htslib.sam_hdr_dup(x.sam_hdr)
    dup == C_NULL && error("failed to duplicate sam_hdr: $(x.sam_hdr)")
    SamHeader(dup)
end

function BioGenerics.seqname(hdr::SamHeader, tid::Integer)
    GC.@preserve hdr begin
        ctid = Cint(tid - Cint(1)) # one-based to zero-based
        ptr = htslib.sam_hdr_tid2name(pointer(header(hf)), ctid)
        ptr == C_NULL && throw(KeyError(tid))
        chr_str = unsafe_string(ptr)
        chr_str
    end
end

function refname(hdr::SamHeader, tid::Integer)
    seqname(hdr, tid)
end

function refid(hdr::SamHeader, chrname::AbstractString)::Int
    GC.@preserve hdr begin
        chrname = convert(String, chrname)
        ret::Int = htslib.sam_hdr_name2tid(pointer(hdr), pointer(chrname))
        ret == -1 && throw(KeyError(chrname))
        ret == -2 && error("header could not be parsed")
        ans = ret + 1 # zero-based to one-based
        ans
    end
end

### TODO: wrap sam_hdr_nref and sam_hdr_tid2len

mutable struct HTSReadWriter{T<:IO}
    # hold the reference to objects below to avoid gc
    io::T
    hfile_backend::Ref{htslib.hFILE_backend}
    # htsFile pointer
    htsFile::Ptr{htslib.htsFile}
    # hFILE pointer
    hfile::Ptr{htslib.hFILE_julia}
    # samheader struct
    samheader::Union{SamHeader,Nothing}
    function HTSReadWriter{T}() where {T}
        new{T}()
    end
end

function Base.pointer(x::HTSReadWriter)
    x.htsFile
end

function BioGenerics.header(hf::HTSReadWriter)
    hf.samheader
end

function BioGenerics.IO.stream(hf::HTSReadWriter)
    hf.io
end

function SamHeader(hf::HTSReadWriter)
    SamHeader(header(hf)) # copy
end

function SamHeader(x::String)
    ptr = htslib.sam_hdr_parse(sizeof(x), pointer(x))
    ptr == C_NULL && error("error parsing samheader")
    SamHeader(ptr)
end

function Base.string(x::SamHeader)
    x.text
end

function Base.show(io::IO, x::SamHeader)
    s = string(x)
    lines = split(s, "\n")
    lines = filter(l -> l != "", lines)
    SQlines = filter(l -> l[1:3] == "@SQ", lines)
    first_seq = split(split(SQlines[1], "\t")[2], ":")[2]
    last_seq  = split(split(SQlines[end], "\t")[2], ":")[2]
    print(io, typeof(x), ": ", length(lines), " lines; ", length(SQlines), " sequences '$(first_seq)...$(last_seq)';")
    x
end

"""
    HTSReadWriter(source, mode="r"; index=nothing, gcsclient=nothing, kwargs...)

Create a reader/writer for SAM/BAM/CRAM files.

`mode` may be specified to indicate the file format when writing to a file. For example,
"wbz" means writing to bam file compressed with bgzf. See below for a list of available options.

- 'r'/'w'/'a': read or write or append to a file.
- 'b'/'c'/'f'/'F': binary format, cram file, fastq format, fasta format
- 'z'/'g'/'u': bgzf compressed, gzip compressed, no compression
"""
function HTSReadWriter(source::AbstractString, mode="r"; index=nothing, gcsclient=nothing, kwargs...)
    # reuse gcsclient to download bam index and construct the reader
    if startswith(source, "gs://") && isnothing(gcsclient)
        gcsclient = GCSClient()
    end
    if isa(index, AbstractString) && startswith(index, "gs://")
        index = read(remote_io(index; gcsclient=gcsclient))::Vector{UInt8}
    end

    ## TODO
    write = startswith(mode, "r") ? false : true
    bamio = remote_io(source; write=write, gcsclient=gcsclient)
    HTSReadWriter(bamio, mode; index=index, kwargs...)
end

function HTSReadWriter(source::IO, mode="r"; index=nothing, gcsclient=nothing)
    # reusing gcsclient to download index data if gcsclient is specified
    if isa(index, AbstractString) && startswith(index, "gs://") && (!isnothing(gcsclient))
        index = read(remote_io(index; gcsclient=gcsclient))
    end
    hfileptr = htslib.hfile_init(sizeof(htslib.hFILE_julia), pointer(mode), 0)
    hfileptr = convert(Ptr{htslib.hFILE_julia}, hfileptr)
    hfile_backend = Ref(htslib.hfile_julia_backend())
    ## Add argument's pointer to the struct and assign julia's backend implementation
    let
        hfile = unsafe_load(hfileptr)
        hfile_base2 = htslib.hFILE_base(
            hfile.base.buffer   ,
            hfile.base.beginz   ,
            hfile.base.endz     ,
            hfile.base.limit    ,
            pointer_from_objref(hfile_backend),
            hfile.base.offset   ,
            hfile.base.remaining,
        )
        hfile2 = htslib.hFILE_julia(hfile_base2, pointer_from_objref(source))
        unsafe_store!(hfileptr, hfile2)
    end
    ## htsFile pointer
    htsFileptr = htslib.hts_hopen(hfileptr, pointer("/dummy"), pointer(mode))
    ## wrapper type holding both julia object and pointers managed in C
    #h = HTSReadWriter{typeof(source)}(htsFileptr, hfileptr, source, hfile_backend)
    h = HTSReadWriter{typeof(source)}()
    let
        h.io = source
        h.hfile_backend = hfile_backend
        h.htsFile = htsFileptr
        h.hfile = hfileptr

        Base.finalizer(h) do h
            #h.hfile == C_NULL && error("invalid pointer")
            #h.htsFile == C_NULL && error("invalid pointer")
            if h.htsFile != C_NULL
                ret = htslib.hts_close(h.htsFile)
                ret < 0 && Base.systemerror("error running hts_close with exit code $ret")
            end
            # Looks like that hts_close will also free the underlying hfile. So we don't need it here.
            # htslib.hfile_destroy(hfileptr)
        end
    end
    # TODO
    # let
    #     # set number of threads to use
    #     hts_set_threads!(h, nthreads)
    # end
    let
        # read sam header
        sam_hdr = htslib.sam_hdr_read(htsFileptr)
        if sam_hdr == C_NULL
            h.samheader = nothing
        else
            h.samheader = SamHeader(sam_hdr)
        end
    end
    # load index.
    htssetindex!(h, index)
    h
end

function format_description(hf::HTSReadWriter)
    GC.@preserve hf begin
        unsafe_string(htslib.hts_format_description(htslib.hts_get_format(pointer(hf))))
    end
end

# TODO: also implement it for HTSRegionsIterator

function Base.open(::Type{HTSReadWriter}, args...; kwargs...)
    HTSReadWriter(args...; kwargs...)
end

function Base.open(f::Function, ::Type{HTSReadWriter}, args...; kwargs...)
    io = HTSReadWriter(args...; kwargs...)
    try
        f(io)
    finally
        close(io)
    end
end

"""
    htssetindex!(htsio::HTSReadWriter, index::AbstractString)

Load index of hts file. Supports local file, http url and gs url.
"""
#function htssetindex!(htsio::HTSReadWriter, index)
#    htssetindex!(htsio, HTSIndexData(index))
#end

function htssetindex!(htsio::HTSReadWriter, index::AbstractString)
    # short circuit for local file (avoid current roundtrip of disk->memory->disk->memory)
    if index isa AbstractString && (~occursin(r"^.*://.*$", index))
        ptr = htslib.hts_idx_load2(pointer("/dummy"), pointer(index))
        ptr == C_NULL && error("error loading index file $(index)")
        let
            tmp = unsafe_load(pointer(htsio))
            tmp = @set tmp.idx = ptr
            unsafe_store!(pointer(htsio), tmp)
        end
        return htsio
    end

    htssetindex!(htsio, remote_io(index))
end

# remove index
function htssetindex!(htsio::HTSReadWriter, index::Nothing)
    # hts_close will also free the index, otherwise we need to call hts_idx_destroy
    oldidx = unsafe_load(pointer(htsio)).idx
    let
        tmp = unsafe_load(pointer(htsio))
        tmp = @set tmp.idx = C_NULL
        unsafe_store!(pointer(htsio), tmp)
    end
    oldidx != C_NULL && htslib.hts_idx_destroy(oldidx)
    return htsio
end

function htssetindex!(htsio::HTSReadWriter, index::IO)
    data = read(index)::Vector{UInt8}
    htssetindex!(htsio, data)
end

function htssetindex!(htsio::HTSReadWriter, index::Vector{UInt8})
    # this is a dirty solution
    ptr = mktemp() do tmppath, tmpio
        write(tmpio, index)
        flush(tmpio)
        ptr = GC.@preserve tmppath htslib.hts_idx_load2(pointer("/dummy"), pointer(tmppath))
        ptr == C_NULL && throw("Failed to load hts index file")
        ptr
    end

    # if index already exists, call hts_idx_destroy on the old index
    oldidx = unsafe_load(pointer(htsio)).idx

    # assign new index
    let
        tmp = unsafe_load(pointer(htsio))
        tmp = @set tmp.idx = ptr
        unsafe_store!(pointer(htsio), tmp)
    end

    oldidx != C_NULL && htslib.hts_idx_destroy(oldidx)
    return htsio
end

function Base.show(io::IO, hf::HTSReadWriter)
    if pointer(hf) == C_NULL
        print(io, typeof(hf), " Invalid HTSReadWriter")
        return nothing
    end
    if isnothing(header(hf))
        print(io, typeof(hf), " ", format_description(hf), " with no header")
    else
        nseq = length(filter(startswith("@SQ"), split(string(header(hf)), "\n")))
        print(io, typeof(hf), " ", format_description(hf), " with ", nseq, " reference sequences")
    end
end

# Methods
# -------

# During initialization

function hts_set_threads!(hf::HTSReadWriter, n)
    ret = htslib.hts_set_threads(hf.htsFile, convert(Int32, n))
    ret != 0 && error("error running hts_set_threads")
    hf
end

function hts_set_cache_size!(hf::HTSReadWriter, n)
    htslib.hts_set_cache_size(hf.htsFile, convert(Int32, n))
end

# Read

function BioGenerics.IO.tryread!(hf::HTSReadWriter, record::HTSRecord)::Union{HTSRecord,Nothing}
    ret = htslib.sam_read1(pointer(hf), pointer(hf.samheader), record.ptr)
    ret < -1 && error("error reading record")
    ret == -1 && return nothing
    return record
end

function Base.read!(hf::HTSReadWriter, record::HTSRecord)
    isnothing(tryread!(hf, record)) && throw(EOFError())
    record
end

function Base.read(hf::HTSReadWriter)
    read!(hf, HTSRecord())
end

# eof, flush, close

# TODO: eof, see https://github.com/samtools/htslib/issues/720
# TODO: flush, use htslib.hflush

function Base.close(hf::HTSReadWriter)
    # finalizer will run hts_close, which will close the underlying IO from julia
    Base.finalize(hf)
    hf.htsFile = C_NULL
    return nothing
end

# TODO: iswriteable, isreadable (false after hts_close) and isreadonly

# write

function Base.write(hf::HTSReadWriter, samheader::SamHeader)
    ## TODO: check writable
    GC.@preserve hf samheader begin
        ret = htslib.sam_hdr_write(pointer(hf), pointer(samheader))
        ret != 0 && error("error writing samheader to file ($ret)")
        # if header has been successfully written, we update the cached header associated with the struct
        hf.samheader = samheader
    end
    return nothing # base api returns the number of bytes written
end

function Base.write(hf::HTSReadWriter, samheader::SamHeader, record::HTSRecord)::Int
    GC.@preserve hf samheader record begin
        ret = htslib.sam_write1(pointer(hf), pointer(samheader), pointer(record))
        ret < 0 && error("error writing record to file ($ret)")
        return ret
    end
end

function Base.write(hf::HTSReadWriter, record::HTSRecord)
    samheader = header(hf)
    write(hf, samheader, record)
end


text:
using BenchmarkTools

cd(dirname(@__FILE__))
const lines = open("one.txt", "r") do f
    readlines(f)
end

ints = map(lines) do l     # don't make this const, otherwise @btime does not give valid results
    parse(Int64, l)
end
const len = length(ints)   # but this can be const and makes a big difference for the benchmarks

x = y = 0
@btime for i in 1:len
    for j in i+1:len
        if ints[i]+ints[j] == 2020
            global x = ints[i]
            global y = ints[j]
            break
        end
    end
end
println(x," ", y, " ", x+y, " ", x*y)

# More functional-style version that uses list comprehensions as generators for the indices
pair_product(a,b) = ((i,j) for i in a:b for j in i+1:b)

x = y = 0
@btime begin
    (i,j) = first((i,j) for (i,j) in pair_product(1, len) if ints[i]+ints[j] == 2020)
    global x = ints[i]
    global y = ints[j]
end
println(x," ", y, " ", x+y, " ", x*y)


text:
module Sound
	const SPEED_OF_SOUND = 343 # m/s @ 20C
	const AIR_DENSITY = 1.225 # kg/m^3 @ 20C

	struct SoundPlaneOpts
		width::Int # m
		height::Int # m
		pointsPerUnit::Int # How many data points per meter
	end

	"""
		Calculates amplitude at a point
	"""
	function calcAmplitude(f::Float64, d::Float64)
		return sin((d * f * 2 * pi) / SPEED_OF_SOUND)
	end

	"""
		Calculates distance falloff coefficient at a point
	"""
	function calcDistFalloff(d::Float64)
		return 1 / d^2
	end

	"""
		Generates a base sound plane (A (w * ppu) x (h * ppu) array of 1)
	"""
	function genSoundPlane(opts::SoundPlaneOpts)
		return ones(Float64, (opts.height * opts.pointsPerUnit, opts.width * opts.pointsPerUnit))
	end

	"""
		Precalculate the distances to audio source
	"""
	function genDistPlane(opts::SoundPlaneOpts, origin::Tuple{Float64, Float64})
		arrDims = (opts.width * opts.pointsPerUnit, opts.height * opts.pointsPerUnit)
		arr = Array{Float64}(undef, arrDims)

		for c in range(1, arrDims[1])
			for r in range(1, arrDims[2])
				x = c / opts.pointsPerUnit - origin[1]
				y = r / opts.pointsPerUnit - origin[2]

				arr[c,r] = sqrt(x^2 + y^2)
			end
		end

		return arr
	end

	"""
		Generates an amplitude plane from a single audio source
	"""
	function genAmplitudePlane(distPlane::Array{Float64}, f::Float64)
		return calcAmplitude.(f, distPlane)
	end

	"""
		Generates a distance falloff plane
	"""
	function genFalloffPlane(distPlane::Array{Float64})
		return calcDistFalloff.(distPlane)
	end

	function genDbPlane(ampPlane::Array{Float64}, falloffPlane::Array{Float64})
		return 10 .* log10.(ampPlane .^ 2 .* falloffPlane .^ 2)
	end

	export SoundPlaneOpts
	export genSoundPlane, genAmplitudePlane, genFalloffPlane, genDbPlane, genDistPlane
end


text:
srand(1)

let
@testset "orthogonality $((m, n))" for (m, n) in ((5, 5), (10, 4), (3, 7))
    X = randorth(m, n)
    @test isorth(X)
end

@testset "singular values $((m, n))" for (m, n) in ((5, 5), (10, 4), (3, 7))
    sv = [0.01; 1.0; 0.1]
    k = min(m, n)
    X = randrealsv(sv, m, n)
    U, S, V = svd(X)
    @test S ≈ sort([sv; zeros(k-length(sv))], rev=true)
    @test size(U) == (m, k)
    @test size(V) == (n, k)
end

@testset "sparse matrices $((m, n, k))" for (m, n, k) in ((10, 5, 3), (5, 10, 4))
    X = randsparse(m, n, k, 0.3)
    @test count(diag(qrfact(X).R) .!= 0) == k
end

end # let


text:
include("../../../opti/src/opti_module.jl")

module ntest
using ..opti

using ..noisydiff
using PyPlot
PyPlot.pygui(true)

function test()

    t = collect( 0:0.01:7 )
    y_theo = sin.(t) # abs.(t .- 3.5)# abs.(t .- 3.5)
    noise = 0.01 * sin.(100 .* t)
    y_noisy = y_theo .+ noise # 0.01 .* randn(size(y_theo))
    y_est = noisydiff.pt1_smooth(t, y_noisy; T = 0.1)

    # @time (x_sol, y_sol, iter, stopping_crit) = noisydiff.tvd(t, y_noisy, lamda = 1 * 0.2,
    #         tol_options = opti.ToleranceOptions(max_iter = 500, step_size_tol = 1e-20, f_abs_tol = 1e-15))
    # dy_tvd_fw = stat.x_sol
    # dy_tvd_fw = x_sol
    # show(stat)
    # y_tvd = y_noisy[1] .+ noisydiff.cumtrapezoidal(t, dy_tvd_fw)

    t_mids = 0.5 * (t[1:end-1] + t[2:end])
    dy_theo = cos.(t_mids)
    dy_straightfw = noisydiff.diff_fw(t, y_noisy)
    dy_pt1_fw = noisydiff.diff_fw(t, y_est)

    w1 = 50
    t_lsq_d1 = t[(w1+1):end]
    dy_lsq_d1 = noisydiff.diff_sliding_least_squares_d1(t, y_noisy, w = w1)

    w3 = w1# 3*w1
    t_lsq_d3 = t[(w3+1):end]
    dy_lsq_d3 = noisydiff.diff_sliding_least_squares_d3(t, y_noisy, w = w3)

    # dy_tvd_fw = noisydiff.diff_fw(t, y_tvd)

    # PyPlot.figure()
    # PyPlot.grid()
    # PyPlot.plot(t, y_noisy)
    # # PyPlot.plot(t, y_theo, linewidth = 2)
    # PyPlot.plot(t, y_est, linewidth = 2)
    # PyPlot.plot(t, y_tvd, linewidth = 2)

    # PyPlot.figure()

    # PyPlot.subplot(2,1,1)
    # PyPlot.grid()
    # PyPlot.plot(t, y_noisy)

    # PyPlot.subplot(2,1,2)
    # PyPlot.grid()
    # # PyPlot.plot(t_mids, dy_theo, linewidth = 2)
    # # PyPlot.plot(t_mids, dy_straightfw)
    # PyPlot.plot(t_mids, dy_pt1_fw, linewidth = 2)
    # PyPlot.plot(t, dy_tvd_fw, linewidth = 4)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(t_mids, dy_pt1_fw, linewidth = 2)
    PyPlot.plot(t_lsq_d1, dy_lsq_d1, linewidth = 2)
    PyPlot.plot(t_lsq_d3, dy_lsq_d3, linewidth = 2)

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
type Errors
    a::Float64 # parameter distribution smoothing
    h::Float64

    odev::Vector{Float64} # σ for observation, param distributions
    pdev::Vector{Float64}

    lb::Vector{Float64} # lower/upper bounds for truncated normal param distrs.
    ub::Vector{Float64}

    # Amin::Vector{Float64} # minimum values for A0 for non-measured arteries
    # bmax::Vector{Float64}

    function Errors(nparams=11)
        this = new();
        δ = 0.99; # ∃ (0, 1, typically within 0.95-0.99), lower values = less parameter dispersion
        this.a = (3*δ-1)/(2*δ);
        this.h = sqrt.(1-this.a^2);
        this.odev = [25*mmTom,25*mmTom,25*mmTom];
        # this.pdev = 1e-12*ones(nparams);
        this.pdev = [1e9,2e9,1e10,4e-11,2e-11,2e-12,0.2,0.2,0.2,0.2,0.2,0.2,8e6,8e6,1e7];
        # this.pdev = [1e9,2.3e9,1e10,4e-11,2e-11,4.8e-12,0.1,0.3,0.05,0.05,0.05,4e-7,4.3e-7,2e-7];
        # this.pdev = [1e9,2.3e9,1e10,4e-11,2e-11,4.8e-12,0.1,0.3,7e5,7e5,7e5];
        # this.pdev = [1e9,2.3e9,1e10,4e-11,2e-11,4.8e-12,0.1,0.1,4e-7,4.3e-7,2e-7];
        this.lb = -Inf*ones(nparams);
        this.ub = Inf*ones(nparams);
        this.lb[1] = 1e10; # R1
        this.lb[2] = 1e10; # R2
        this.lb[3] = 1e10; # R3
        this.lb[4] = 1e-10; # C1
        this.lb[5] = 1e-10; # C2
        this.lb[6] = 7e-12; # C3
        this.lb[7] = 0.1; # Cfac
        this.lb[8] = 0.1; # Rfac
        this.lb[9] = 0.1; # Cfac
        this.lb[10] = 0.1; # Rfac
        this.lb[11] = 0.1; # Cfac
        this.lb[12] = 0.1; # Rfac
        # this.lb[9] = 2e-6; # A1
        # this.lb[10] = 2.15e-6; # A2
        # this.lb[11] = 9.5e-7; # A3
        this.lb[13] = 2e7; # β1
        this.lb[14] = 2e7; # β2
        this.lb[15] = 2e7; # β3
        this.ub[1] = 1e11;
        this.ub[2] = 1e11;
        this.ub[3] = 1e11;
        this.ub[4] = 1e-9;
        this.ub[5] = 1e-9; #2.5e-10
        this.ub[6] = 2e-11;
        # this.ub[7] = 1.5;
        # this.ub[8] = 1.2;
        # this.ub[9] = 1.5;
        # this.ub[10] = 1.2;
        # this.ub[11] = 1.5;
        # this.ub[12] = 1.2;
        this.ub[13] = 2e8; # β1 8.5e7
        this.ub[14] = 2e8; # β2 1.5e8
        this.ub[15] = 2e8; # β3 5.5e8
        # this.ub[9] = 6e7; # E1
        # this.ub[10] = 6e7; # E2
        # this.ub[11] = 6e7; # E3

        # pull in artery data from text file
        temp = LegHemodynamics.loadtexttree("arterylist_2.txt");
        # this.Amin = 0.8*[temp[1,:A0_m2]];
        # for i = 2:length(temp[:A0_m2])
        #     push!(this.Amin,0.8*temp[i,:A0_m2])
        # end
        # this.bmax = 3*[temp[1,:beta_Pam]];
        # for i = 2:length(temp[:beta_Pam])
        #     push!(this.bmax,3*temp[i,:beta_Pam])
        # end
        return this
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using DrWatson
@quickactivate "ChaosThroughBilliards"
include(srcdir("style.jl"))

x = 1.0
y = 0.54

si = billiard_sinai(y/4, x, y)

ps = [Particle(4x/5, y/5, π/9), Particle(4x/5, y/5 + 1e-5, π/9)]

# interactive_billiard(re, ps;
# backgroundcolor = :black, res = MAXRES, particle_size = 2.0, add_controls = false)

billiard_video(
    videodir("two_initial.mp4"), si, deepcopy(ps);
    frames = 1200, backgroundcolor = :black, colors = COLORS[1:2],
    tailwidth = 3.5, particle_size = 2, res = (1200*x, 1200*y),
    plot_particles = true,
    dt = 0.0001, speed = 150, tail = 10000, # this makes ultra fine temporal resolution
)

# %% Same video, with timeseries now
f(p) = p.pos[2] # the function that obtains the data from the particle
ylabel = "y"

billiard_video_timeseries(
    videodir("two_initial_timeseries.mp4"), si, deepcopy(ps), f;
    frames = 1200, backgroundcolor = :black, colors = COLORS[1:2],
    tailwidth = 3.5, particle_size = 2, res = (1200*x, 1200*y*2),
    plot_particles = true,
    dt = 0.0001, speed = 150, tail = 10000, # this makes ultra fine temporal resolution
    ylabel
)

# %% Also plot of timeseries and hopefully show exponential divergence
x1, y1, vx1, vy1, t1 = DynamicalBilliards.timeseries(ps[1], si, 20.0; dt = 0.1)
x2, y2, vx2, vy2, t2 = DynamicalBilliards.timeseries(ps[2], si, 20.0; dt = 0.1)
using LinearAlgebra
δ = Float64[]
for i in 1:min(length(t1), length(t2))
    push!(δ, sqrt((x1[i]-x2[i])^2 + (y1[i]-y2[i])^2 + (vx1[i]-vx2[i])^2 + (vy1[i]-vy2[i])^2))
end

fig = Figure()
ax = fig[1,1] = Axis(fig, xlabel = L"t", ylabel = L"\log(\delta)", xlabelsize= 32, ylabelsize= 32)
lines!(ax, log.(δ); color = COLORS[3])
# ax.xlabel = "time"
# ax.ylabel = "log(δ)"
display(fig)
# δ = abs.(y2[1:length(y1)] .- y1)
GLMakie.save(videodir("divergence.png"), fig)

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
## Test plot functions
function testModifyY()
    p0, xArray, tArray = getInitialDist()

    pMax = maximum(p0)
    p = copy(p0)

    D = 1E-5
    dx = 1/(xArray.len-1);
    dt = 0.5*dx*dx/D;
    s = D*dt/(dx*dx);
    function stepTime(p, s, nx)
        p[2:end-1] = s.*p[3:nx]-(2*s-1).*p[2:nx-1]+s.*p[1:nx-2];
    end

    YplotObject = @makeLivePlot myPlot(xArray, p0)
    for tt in tArray
        stepTime(p, s, xArray.len)
        modifyPlotObject!(YplotObject, arg2 = p)
    end
end
function testModifyX()
    p0, xArray, tArray = getInitialDist()
    pMax = maximum(p0)

    XplotObject = @makeLivePlot myPlot(xArray, p0)
    for tt in tArray
        modifyPlotObject!(XplotObject, arg1 = XplotObject[][1]*0.99)
    end
end
function testModifyXY()
    p0, xArray, tArray = getInitialDist()

    pMax = maximum(p0)

    XYplotObject = @makeLivePlot myPlot(xArray, p0)
    for tt in tArray
        modifyPlotObject!(XYplotObject, arg2 = XYplotObject[][2]*0.99, arg1 = XYplotObject[][1]*1.01)
    end
end
function testModifyZ()
    x = 1:0.05:20
    y = 1:0.05:10
    f(x, y, t) = begin
        (3x + y^2) * abs(sin(x) + (1-t)*cos(y))
    end
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    Z = map((x, y) -> f(x, y, 0.0), X, Y)

    ZplotObject = @makeLivePlot myPlotZ(x, y, Z)

    ttt = 0.0:0.1:1.0
    for tt in ttt
        Z = map((x, y) -> f(x, y, tt), X, Y)
        modifyPlotObject!(ZplotObject, arg3 = Z)
    end
end
function testModifyXText()
    p0, xArray, tArray = getInitialDist()

    pMax = maximum(p0)
    p = copy(p0)
    titleText = "Title, step: "

    XtextPlotObject = @makeLivePlot myPlotTitle(xArray, p0, titleText)
    for tt in 1:50
        modifyPlotObject!(XtextPlotObject, arg3 = string(titleText, tt), arg1 = XtextPlotObject[][1]*0.99)
    end
end
function testModifyXColour()
    p0, xArray, tArray = getInitialDist()

    pMax = maximum(p0)
    p = copy(p0)
    scatterColour = range(0, 1, length=50)

    XColourPlotObject = @makeLivePlot myPlotColour(xArray, p0, scatterColour[1])

    for tt in 1:50
        modifyPlotObject!(XColourPlotObject, arg3 = scatterColour[tt], arg1 = XColourPlotObject[][1]*0.99)
    end
end
function testAddXY()
    p0, xArray, tArray = getInitialDist()

    pMax = maximum(p0)
    p = copy(p0)

    addXYPlotObject = @makeLivePlot myPlot(xArray, p0)

    for tt in 1:length(xArray)
        modifyPlotObject!(addXYPlotObject, arg1 = xArray[1:tt], arg2 = p0[1:tt])
    end
end
function testAddXYZ()
    attractor = Lorenz()
    plt = plot3d(
        1,
        xlim = (-30, 30),
        ylim = (-30, 30),
        zlim = (0, 60),
        title = "Lorenz Attractor",
        label = "",
        marker = 2,
    )

    addXYZPlotObject = @makeLivePlot myPlotXYZ(plt)

    for tt in 1:500
        step!(attractor)
        push!(plt, attractor.x, attractor.y, attractor.z)
        modifyPlotObject!(addXYZPlotObject, arg1 = plt)
    end
end
function getInitialDist()
    tStart = 0.0
    tEnd = 2.0
    ntSteps = 100
    xStart = -5.0
    xEnd = 5.0
    nxSteps = 20

    tArray = range(tStart, tEnd, length = ntSteps)
    xArray = range(xStart, xEnd, length = nxSteps)

    x0 = 0.0
    sig = 2.0
    f(x) = exp(-((x - x0)/sig)^2)
    p = f.(xArray)

    return p, xArray, tArray
end
Base.@kwdef mutable struct Lorenz
    dt::Float32 = 0.02
    σ::Float32 = 10
    ρ::Float32 = 28
    β::Float32 = 8/3
    x::Float32 = 1
    y::Float32 = 1
    z::Float32 = 1
end
function step!(l::Lorenz)
    dx = l.σ * (l.y - l.x);         l.x += l.dt * dx
    dy = l.x * (l.ρ - l.z) - l.y;   l.y += l.dt * dy
    dz = l.x * l.y - l.β * l.z;     l.z += l.dt * dz
end
function myPlot(xx, yy)
    sleep(0.0001)
    plot(
        xx, yy,
        label = "",
        color = :red,
        xlim = (-5, 5),
        ylim = (0, 1),
        title = "Title",
        xlabel = "X label",
        ylabel = "Y label"
    )
    scatter!(xx, yy, label = "")
end
function myPlotTitle(xx, yy, titleText)
    sleep(0.0001)
    plot(
        xx, yy,
        label = "",
        color = :red,
        xlim = (-5, 5),
        ylim = (0, 1),
        title = titleText,
        xlabel = "X label",
        ylabel = "Y label"
    )
    scatter!(xx, yy, label = "")
end
function myPlotColour(xx, yy, cc)
    sleep(0.0001)
    plot(
        xx, yy,
        label = "",
        color = :red,
        xlim = (-5, 5),
        ylim = (0, 1),
        title = "Title text",
        xlabel = "X label",
        ylabel = "Y label"
        )
    scatter!(xx, yy, markersize = 20, color = RGBA(cc, 0.5, 0, 0), label = "")
end
function myPlotZ(xx, yy, ZZ)
    sleep(0.0001)
    p1 = contour(
        xx, yy, ZZ,
        fill = true, clims=(0,300),
        title = "Title text",
        xlabel = "X label",
        ylabel = "Y label"
    );
end
function myPlotXYZ(pltObj)
    sleep(0.0001)
    xx = getindex(pltObj.series_list[1].plotattributes, :x)
    yy = getindex(pltObj.series_list[1].plotattributes, :y)
    zz = getindex(pltObj.series_list[1].plotattributes, :z)
    plot3d(
        xx, yy, zz,
        xlim = (-30, 30),
        ylim = (-30, 30),
        zlim = (0, 60),
        title = "Lorenz Attractor",
        label = "",
        marker = 2,
    )
end


text:
using DSGE
using Nullables, DataFrames, OrderedCollections, Dates

##############
# Model Setup
##############
# Instantiate the FRBNY DSGE model object
m = Model1002()

# estimate as of 2015-Q3 using the default data vintage from 2015 Nov 27
m <= Setting(:data_vintage, "181115")
m <= Setting(:date_forecast_start, quartertodate("2018-Q4"))

#############
# Estimation
#############
# reoptimize parameter vector, compute Hessian at mode, and full posterior
# parameter sampling
estimate(m)

# produce LaTeX tables of parameter moments
moment_tables(m)

############
# Forecasts
############
# (1) forecast_one produces forecasts of the variables specified in output_vars
# and stores the forecasts in model units (log deviations from steady state).
# (2) compute_meansbands stores the forecast in a MeansBands type, defined in DSGE.jl,
# which transforms the forecast from model units to "human readable" units, i.e.
# the GDP forecast is transformed from annualized quarterly log per-capita growth rates
# to annualized quarterly aggregate percent change.

# The variables we want to forecast. In this case, all of the model observables
output_vars = [:histobs, :forecastobs]

m <= Setting(:use_population_forecast, false)

# Modal forecast (point forecast)
forecast_one(m, :mode, :none, output_vars)
compute_meansbands(m, :mode, :none, output_vars)

# Full-distribution forecast (point forecast (mean) and uncertainty bands)
# Optionally add 10 processes to run the forecast in parallel (uncomment the 3 lines below).
# Alternatively, you can load the ClusterManagers package and add processes
# using one of the schedulers such as SGE or Slurm.
addprocs(10)
@everywhere using DSGE
m <= Setting(:use_parallel_workers, true)

forecast_one(m, :full, :none, output_vars)
compute_meansbands(m, :full, :none, output_vars)

rmprocs(procs())


text:
module Stat31020
using JuMP
using MathProgBase
"""
...
# Arguments
- `x::Array{Float64,1}`: current iterate
- `n::Int64`: number of derivatives to query
- `kkt::Bool`: T/F to return KKT matrix
- `μ::Float64`: objective function multiplier (default = 1.0)
- `λ::Array{Float64,1}`: constraint multiplier (default = [1.0, ..., 1.0])
# Output
- `f::SparseMatrixCSC`: function evaulation
- `g::SparseMatrixCSC`: gradient evaulation
- `H::SparseMatrixCSC`: sparse, symmetric hessian evaulation
- `J::SparseMatrixCSC`: sparse constraint jacobian evaluation
- `KKT::SparseMatrixCSC`: sparse KKT system evaulation
...
"""
function constrained_wrap(x::Array{Float64,1}, nderivs::Int64;
                          kkt=true, μ=nothing, λ=nothing)
    ## setup
    n = length(x)                          ## problem dimension
    m = 0                                  ## number of constraints (TBD below)
    I = 1:n-1                              ## index set

    ## -------------------------------------------------------------------------
    ## TODO: fill in function evaluation
    ## -------------------------------------------------------------------------
    f = sum(x[i]^4 - 2.0*x[i]^2*x[i+1] + x[i]^2 + x[i+1]^2 - 2*x[i] for i=I)
    ## -------------------------------------------------------------------------

    if nderivs == 0
        return f

    elseif nderivs >= 1
        ## initialize JuMP model and MathProgBase interface
        NLmodel = Model()                   ## JuMP object
        @variable(NLmodel, y[1:n])
        ## ---------------------------------------------------------------------
        ## TODO: fill in objective function and constraints
        ## ---------------------------------------------------------------------
        @NLobjective(NLmodel, Min, sum(y[i]^4 -
                                       2.0*y[i]^2*y[i+1] +
                                       y[i]^2 +
                                       y[i+1]^2 -
                                       2*y[i] for i=I))
        @NLconstraint(NLmodel, -(y[1] + 0.25)^2 + 0.75*y[2] >= 0)
        @NLconstraint(NLmodel, -(y[1] + 0.25) + 0.75*y[2]^2 >= 0)
        ## ---------------------------------------------------------------------
        d = JuMP.NLPEvaluator(NLmodel)      ## MPB object
        m = MathProgBase.numconstr(NLmodel) ## number of contraints
        ## initialize dual variable values (if not specified)
        if μ == nothing
            μ = 1.
        end
        if λ == nothing
            λ = ones(m)
        end

        ## evaluate derivatives
        if nderivs == 1
            MathProgBase.initialize(d, [:Grad])
            ## evaluate gradient
            g = zeros(n)
            MathProgBase.eval_grad_f(d, g, x)
            return (f, g)

        elseif nderivs == 2
            MathProgBase.initialize(d, [:Grad, :Hess, :Jac])
            ## evaluate Hessian
            ij_H = MathProgBase.hesslag_structure(d)
            h = ones(length(ij_H[1]))
            MathProgBase.eval_hesslag(d, h, x, μ, λ)
            H = populate_hess_sparse(ij_H[1], ij_H[2], h, n)
            ## evalute gradient (partially pre-computed from above)
            g = zeros(n)
            MathProgBase.eval_grad_f(d, g, x)
            ## evaluate Jacobian (constraints)
            if kkt == true
                ij_J = MathProgBase.jac_structure(d)
                j = ones(length(ij_J[1]))
                MathProgBase.eval_jac_g(d, j, x)
                J = sparse(ij_J[1], ij_J[2], j)
                KKT = [H J'; J zeros(m, m)]
                println(KKT)
                return (f, g, H, J, KKT)
            else
               return (f, g, H)
            end
        end
    end
end

"""
...
# Arguments
- `i::Array{Int64,1}`: row indices
- `j::Array{Int64,1}`: col indices
- `h::Array{Int64,1}`: values s.t. h[1] = H[i[1], j[1]]
- `n::Int64`: dimension of hessian
# Output
- `H::SparseMatrixCSC`: sparse, symmetric hessian
...
"""
function populate_hess_sparse(
    i::Array{Int64,1},
    j::Array{Int64,1},
    h::Array{Float64,1},
    n::Int64)
    ## build matrix; probably a better way to do this...
    H = sparse([i;j], [j;i], [h;h], n, n, +)  ## NOTE: symmetrizing doubles the diagonal
    H = H - Diagonal(H)/2.                    ## NOTE: required removal of 1/2 digonal
    return H
end
end # module


text:
"""
    indent(str::AbstractString, n::Int)

Indent each non-blank line by `n` spaces.

# Examples
```jldoctest; setup = :(using MultilineStrings)
julia> indent("a\\nb", 4)
"    a\\n    b"

julia> indent("  a\\n  \\n  b", 2)
"    a\\n  \\n    b"
```

See also `Base.unindent` and `Base.indentation`.
"""
function indent(str::AbstractString, n::Int)
    n == 0 && return str
    # Note: this loses the type of the original string
    buf = IOBuffer(sizehint=sizeof(str))
    indent_str = ' ' ^ n

    line_start = firstindex(str)
    blank_line = true
    for (i, ch) in enumerate(str)
        if ch == '\n'
            !blank_line && print(buf, indent_str)
            print(buf, SubString(str, line_start, i))
            line_start = nextind(str, i)
            blank_line = true
        elseif blank_line && !isspace(ch)
            blank_line = false
        end
    end

    # Last line of string that doesn't contain a newline
    !blank_line && print(buf, indent_str)
    print(buf, SubString(str, line_start, lastindex(str)))

    String(take!(buf))
end


text:
using Zygote
@testset "ChainRulesCore" begin
    # 1D example
    x = 1:10
    y = sin.(x)

    # Simple test that rrules give equivalent results to internal gradient function

    itp = interpolate(y,BSpline(Linear()))

    if VERSION ≥ v"1.3"
        @test Zygote.gradient(itp, 1)[1] == Interpolations.gradient(itp, 1)
    else
        @test_skip Zygote.gradient(itp, 1)[1] == Interpolations.gradient(itp, 1)
    end

    # 2D example
    x2 = [i*j for i=1:5, j=1:5]
    y2 = sin.(x2)

    itp2 = interpolate(y2, BSpline(Cubic(Line(OnGrid()))))

    if VERSION ≥ v"1.3"
        @test Zygote.gradient(itp2,1,2)[1] == Interpolations.gradient(itp2,1,2)
    else
        @test_skip Zygote.gradient(itp2,1,2)[1] == Interpolations.gradient(itp2,1,2)
    end
end


text:
using CImGui
using Test

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


text:
using AmplNLWriter, Test

@testset "MOI" begin
    include("MOI_wrapper.jl")
end

@testset "Base" begin
    include("nl_convert.jl")
    include("nl_linearity.jl")
    include("nl_write.jl")
    include("sol_file_parser.jl")
end


text:
module irtfun

using StatsFuns
using FastGaussQuadrature
import Distributions: Normal, MvNormal, Beta, Gamma, pdf, cdf, quantile
using DataFrames: DataFrameRow, DataFrame, All, groupby, nrow, combine, Not
using ForwardDiff: gradient, hessian, Dual
using Printf: @sprintf
using DataFrames: AbstractDataFrame, mapcols
using StaticArrays
using LinearAlgebra: LowerTriangular, diagm, Cholesky, cholesky, det, Hermitian, tril!
using Statistics: mean
using Optim: optimize, BFGS, Newton, minimizer

include("Models.jl")
include("IRTmodels/graded.jl")
include("IRTmodels/guessing.jl")
include("IRTmodels/lpe.jl")
include("GaussQuadrature.jl")
include("Estep.jl")
include("Mstep.jl")
include("EMcycle.jl")
include("StartingValues.jl")
include("index.jl")
include("Rescore.jl")
include("Qmatrix.jl")
include("Estimate.jl")
include("UpdateLatentDistributions.jl")
include("CheckParameters.jl")

export IRTmodel
export irf
export distribute
export distribute!
export take
export rescale!
export graded
export guessing
export LPE

export DNq
export initializenode

export BAEstep!
export initializematrix

export ell
export Newton
export expectedcount

export BAEMcycle
export checkconvergence

export init!
export init

export takeitemindex
export takepersonindex
export setgroupindex!

export rescore

export explanatoryMIRT

export Estep
export Mstep
export estimate

export Moments
export updateAGH!
export updateDN!

export checitemparameters

end


text:
##########################################################################################
######################### Este archivo contiene los tipos usados y    ####################
######################### las funciones que los crean                 ####################
##########################################################################################

type Gato                     ## Tipo gato con el tablero, el cual será un vector
    tablero::Vector{Int8}     ## 0 = vacía, 1 = X, -1 = 0
end

Gato() = Gato(zeros(Int8, 9)) ## Se inicia el tablero vacío

type Jugador        ## Tipo jugador, con su número y marcador
    num::Int8       ## Judador 1 o Jugador 2
    marcador::Int8
end

Jugador(n::Int64) = Jugador(Int8(n), Int8(0))  # Función que inicia un jugador nuevo con su número
                                               # marcador cero


text:
# Test the extensions of Base.convert
@testset "Convert Extensions" begin
    m = Model()
    @variable(m, x[1:2])
    # Array -> SparseAxisArray
    @testset "Array->Sparse" begin
        @test convert(JuMPC.SparseAxisArray,
                      x) isa JuMPC.SparseAxisArray
        @test length(convert(JuMPC.SparseAxisArray, x)) == length(x)
        inds = CartesianIndices(x)
        @test convert(JuMPC.SparseAxisArray,
                      x)[inds[1][1]] == x[inds[1]]
        @test convert(JuMPC.SparseAxisArray,
                      x)[inds[2][1]] == x[inds[2]]
    end
    # DenseAxisArray -> SparseAxisArray
    @testset "DenseArray->Sparse" begin
        @variable(m, y[2:3])
        @test convert(JuMPC.SparseAxisArray,
                      y) isa JuMPC.SparseAxisArray
        @test length(convert(JuMPC.SparseAxisArray, y)) == length(y)
        inds = collect(keys(y))
        @test convert(JuMPC.SparseAxisArray,
                      y)[inds[1][1]] == y[inds[1]]
        @test convert(JuMPC.SparseAxisArray,
                      y)[inds[2][1]] == y[inds[2]]
    end
    # Number -> AffExpr
    @testset "Number->AffExpr" begin
        expr = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) + 42
        @test convert(JuMP.AbstractJuMPScalar, 42) == expr
    end
end

# Test extension of keys
@testset "_keys" begin
    m = Model()
    @variable(m, x[1:2], container = SparseAxisArray)
    @test InfiniteOpt._keys(x) == keys(x.data)
    @test InfiniteOpt._keys(ones(2)) == keys(ones(2))
end

# Test allequal
@testset "_allequal" begin
    @test InfiniteOpt._allequal([1])
    @test !InfiniteOpt._allequal([1, 2])
    @test InfiniteOpt._allequal([1, 1])
    @test InfiniteOpt._allequal([1 1; 1 1])
end

# Test _make_float_info
@testset "_make_float_info" begin
    info1 = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false)
    num = Float64(0)
    info2 = VariableInfo(true, num, true, num, true, num, true, num, true, true)
    @test InfiniteOpt._make_float_info(info1) isa VariableInfo{Float64, Float64, Float64, Float64}
    @test InfiniteOpt._make_float_info(info2) == info2
end

# Test _kwargs2string
@testset "_kwargs2string" begin
    nt1 = (d = 1,)
    nt2 = (d = 1, e = 2)
    @test InfiniteOpt._kwargs2string(nt1) == "d = 1"
    @test InfiniteOpt._kwargs2string(nt2) == "d = 1, e = 2"
end


text:
mutable struct EmpiricalMean{T<:Real,V<:AbstractVector{<:Real}} <: PriorMean{T}
    C::V
    opt::Optimizer
end

"""
**EmpiricalMean**
```julia`
function EmpiricalMean(c::V=1.0;opt::Optimizer=Adam(α=0.01)) where {V<:AbstractVector{<:Real}}
```
Construct a constant mean with values `c`
Optionally give an optimizer `opt` (`Adam(α=0.01)` by default)
"""
function EmpiricalMean(c::V=1.0;opt::Optimizer=Adam(α=0.01)) where {V<:AbstractVector{<:Real}}
    EmpiricalMean{eltype(c),V}(c,opt)
end

function update!(μ::EmpiricalMean{T},grad::AbstractVector{T}) where {T<:Real}
    μ.C .+= update!(μ.opt,grad)
end

Base.:+(x::Real,y::EmpiricalMean{<:Real}) = x.+y.C
Base.:+(x::AbstractVector{<:Real},y::EmpiricalMean{<:Real}) = x+y.C
Base.:+(x::EmpiricalMean{<:Real},y::Real) = y.+x.C
Base.:+(x::EmpiricalMean{<:Real},y::AbstractVector{<:Real}) = y+x.C
Base.:+(x::EmpiricalMean{<:Real},y::EmpiricalMean{<:Real}) = EmpiricalMean(x.C+y.C)
Base.:-(x::Real,y::EmpiricalMean) = x .- y.C
Base.:-(x::AbstractVector{<:Real},y::EmpiricalMean) = x - y.C
Base.:-(x::EmpiricalMean{<:Real},y::Real) = x.C .- y
Base.:-(x::EmpiricalMean{<:Real},y::AbstractVector{<:Real}) = x.C - y
Base.:-(x::EmpiricalMean{<:Real},y::EmpiricalMean{<:Real}) = EmpiricalMean(x.C-y.C)
Base.:*(A::AbstractMatrix{<:Real},y::EmpiricalMean{T}) where {T<:Real} = A*y.C
Base.:*(y::EmpiricalMean{T},A::AbstractMatrix{<:Real}) where {T<:Real} = transpose(y)*A
Base.:convert(::T1,x::EmpiricalMean{T2}) where {T1<:Real,T2<:Real} = T1(x.C)


text:
using MeasureTheory
using Test

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


text:
# For optimal performance, run this script with
#   julia --threads=1 --check-bounds=no detect_crashes_2d.jl
# to let Julia use half of your `Sys.CPU_THREADS` cores. You can also specify
# the number of processes to use via
#   julia --threads=1 --procs=N --check-bounds=no detect_crashes_2d.jl
# where `N` is the number of processes you want to use.

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.precompile()

using Distributed

# Add worker processes if none are present
if nprocs() == 1
  addprocs(Sys.CPU_THREADS ÷ 2)
end

# Package setup required on each processor
@everywhere begin
  using Pkg
  Pkg.activate(@__DIR__)
end
@everywhere using Trixi
@everywhere using OrdinaryDiffEq

# Setup on the driver process
using CSV, DataFrames

@everywhere begin

   Base.@kwdef struct ProblemSetup{NameT, InitialConditionT, TspanT, CoordinateT, SourceT, PeriodicityT, BCT}
        problem_name::NameT
        initial_condition::InitialConditionT
        tspan::TspanT
        coordinates_min::CoordinateT
        coordinates_max::CoordinateT
        source_terms::SourceT
        periodicity::PeriodicityT
        boundary_conditions::BCT
    end

    ## Create problem setups

    function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D; slope = 15)
        # change discontinuity to tanh, domain size is [-1,+1]^2
        # discontinuous function with value 2 for -.5 <= x <= .5 and 0 elsewhere
        B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
        rho = 0.5 + 0.75 * B # rho ∈ [.5, 2]
        v1 = 0.5 * (B - 1)
        v2 = 0.1 * sin(2 * pi * x[1])
        p = 1.0
        return prim2cons(SVector(rho, v1, v2, p), equations)
    end

    @inline function initial_condition_rayleigh_taylor_instability(x, t,
                                                                equations::CompressibleEulerEquations2D,
                                                                slope=15)

        rho1 = 1.0
        rho2 = 2.0
        if x[2] < 0.5
            p = 2 * x[2] + 1
        else
            p = x[2] + 3 / 2
        end

        # smooth the discontinuity to avoid ambiguity at element interfaces
        smoothed_heaviside(x, left, right; slope=slope) = left + 0.5 * (1 + tanh(slope * x)) * (right - left)
        rho = smoothed_heaviside(x[2] - 0.5, rho2, rho1)

        c = sqrt(equations.gamma * p / rho)
        # the velocity is multiplied by sin(pi*y)^6 as in Remacle et al. 2003 to ensure that the
        # initial condition satisfies reflective boundary conditions at the top/bottom boundaries.
        k = 1
        v = -0.025 * c * cos(k * 8 * pi * x[1]) * sin(pi * x[2])^6
        u = 0.0

        return prim2cons(SVector(rho, u, v, p), equations)
    end

    @inline function source_terms_rayleigh_taylor_instability(u, x, t,
                                                            equations::CompressibleEulerEquations2D)
        rho, rho_v1, rho_v2, rho_e = u
        g = 1.0
        return SVector(0.0, 0.0, g * rho, g * rho_v2)
    end

    @inline function initial_condition_richtmeyer_meshkov_instability(x, t,
                                                                    equations::CompressibleEulerEquations2D;
                                                                    slope = 2.0)

        #   Setup used for the Richtmeyer-Meshkov instability. Initial condition adapted from the description of
        #   https://www.youtube.com/watch?v=8K_oe_GhKzM.
        #   The domain is [0, 40/3] x [0, 40]. Boundary conditions are all reflective walls.

        # smooth the discontinuity to avoid ambiguity at element interfaces
        smoothed_heaviside(x, left, right; slope=slope) = left + 0.5 * (1 + tanh(slope * x)) * (right - left)

        L = 40 # domain size
        rho = smoothed_heaviside(x[2] - (18 + 2 * cos(2 * pi * 3 / L * x[1])), 1.0, .25)
        rho = rho + smoothed_heaviside(abs(x[2] - 4) - 2, 3.22, 0.0) # 2 < x < 6
        p = smoothed_heaviside(abs(x[2] - 4) - 2, 4.9, 1.0)
        u = 0.0
        v = 0.0

        return prim2cons(SVector(rho, u, v, p), equations)
    end

    function blow_up_time(element_type, approximation_type,
                          polydeg, cells_per_dimension, problem_setup::ProblemSetup;
                          abstol=1.0e-7, reltol=1.0e-7)

        @unpack tspan, periodicity, coordinates_min, coordinates_max = problem_setup
        @unpack initial_condition, source_terms, boundary_conditions = problem_setup

	@show problem_setup.problem_name, polydeg, element_type, approximation_type, cells_per_dimension

        # 3D compressible Euler equations with ratio of specific heats 7/5
        equations = CompressibleEulerEquations2D(1.4)
        surface_flux = flux_lax_friedrichs
        volume_flux  = flux_ranocha

        dg = DGMulti(polydeg=polydeg, element_type=element_type, approximation_type=approximation_type,
                     surface_integral=SurfaceIntegralWeakForm(surface_flux),
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

        mesh = DGMultiMesh(dg; coordinates_min, coordinates_max,
                           cells_per_dimension=cells_per_dimension,
                           periodicity=periodicity)

        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
                                            source_terms = source_terms,
                                            boundary_conditions = boundary_conditions)

        ode = semidiscretize(semi, tspan)

        alive_callback = AliveCallback(alive_interval = 10000)

        integrator = init(ode, RDPK3SpFSAL49(), abstol=abstol, reltol=reltol,
                          dt = .1 * estimate_dt(mesh, dg), save_everystep=false,
                          callback=alive_callback)

        try
          solve!(integrator)
        catch error
          @info "Blow-up" integrator.t error
        end

        return integrator.t
    end
end

khi_setup = ProblemSetup(problem_name="KHI", initial_condition = initial_condition_kelvin_helmholtz_instability,
                         tspan=(0, 15.0), coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0),
                         source_terms = nothing, periodicity = true,
                         boundary_conditions = Trixi.BoundaryConditionPeriodic())

wall_bcs = (; :entire_boundary => boundary_condition_slip_wall)

rti_setup = ProblemSetup(problem_name="RTI", initial_condition = initial_condition_rayleigh_taylor_instability,
                         tspan=(0, 15.0), coordinates_min = (0.0, 0.0), coordinates_max = (0.25, 1.0),
                         source_terms = source_terms_rayleigh_taylor_instability, periodicity = (false, false),
                         boundary_conditions = wall_bcs)

rmi_setup = ProblemSetup(problem_name="RMI", initial_condition = initial_condition_richtmeyer_meshkov_instability,
                         tspan=(0, 30.0), coordinates_min = (0.0, 0.0), coordinates_max = (40.0 / 3.0, 40.0),
                         source_terms = nothing, periodicity = (false, false),
                         boundary_conditions = wall_bcs)

# Prepare parameters
## DG
mesh_resolutions = [16, 32]

element_types = (Quad(), Tri())
polydegs(elem::Quad) = 1:7
polydegs(elem::Tri) = 1:6 # SBP operators only available up to polydeg = 6
approximation_types(elem::Quad) = (SBP(), GaussSBP())
approximation_types(elem::Tri) = (SBP{StartUpDG.Kubatko{StartUpDG.LegendreFaceNodes}}(), Polynomial())

# TODO: remove later. This is to use Chen-Shu nodes in the SBP tri operators
element_types = (Tri(),)
approximation_types(elem::Tri) = (SBP{StartUpDG.Kubatko{StartUpDG.LegendreFaceNodes}}(),)

khi_parameters = [(element_type, approximation_type, polydeg, (cells_per_dimension, cells_per_dimension), khi_setup)
                  for element_type in element_types for approximation_type in approximation_types(element_type)
                  for cells_per_dimension in mesh_resolutions for polydeg in polydegs(element_type)]
rti_parameters = [(element_type, approximation_type, polydeg, (cells_per_dimension, 4 * cells_per_dimension), rti_setup)
                  for element_type in element_types for approximation_type in approximation_types(element_type)
                  for cells_per_dimension in mesh_resolutions for polydeg in polydegs(element_type)]
rmi_parameters = [(element_type, approximation_type, polydeg, (cells_per_dimension, 3 * cells_per_dimension), rmi_setup)
                  for element_type in element_types for approximation_type in approximation_types(element_type)
                  for cells_per_dimension in mesh_resolutions for polydeg in polydegs(element_type)]

parameters = vcat(khi_parameters, rti_parameters, rmi_parameters)
println("Total $(length(parameters)) parameters runs. ")

blow_up_times = pmap(args -> blow_up_time(args...), parameters)

data = DataFrame(element_type = map(x -> x[1], parameters),
                 approximation_type = map(x -> x[2], parameters),
                 polydeg = map(x -> x[3], parameters),
                 cells_per_dimension = map(x -> x[4], parameters),
                 problem_name = map(x -> x[5].problem_name, parameters),
                 blow_up_time = blow_up_times)

CSV.write(joinpath(@__DIR__, "crash_detection_results_2d.txt"), data)