text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Intersection" begin
    function samplepoints(numsamples, lowu, highu, lowv, highv)
        samples = Array{Tuple{Float64,Float64},1}(undef, 0)

        for i in 0:numsamples
            u = lowu * i / numsamples + (1 - i / numsamples) * highu
            for j in 0:numsamples
                v = lowv * j / numsamples + (1 - j / numsamples) * highv
                push!(samples, (u, v))
            end
        end
        return samples
    end

    @testset "Cylinder" begin
        # Random.seed!(SEED)

        # # random point intersections
        # ur, vr = uvrange(Cylinder)
        # for i in 1:10000
        #     cyl = Cylinder(rand() * 3.0, rand() * 50.0)
        #     u, v = rand() * (ur[2] - ur[1]) + ur[1], rand() * (vr[2] - vr[1]) + vr[1]
        #     pointon = point(cyl, u, v)
        #     origin = SVector(4.0, 4.0, 4.0)
        #     r = Ray(origin, pointon .- origin)
        #     halfspace = surfaceintersection(cyl, r)
        #     @test halfspace !== nothing && !((lower(halfspace) isa RayOrigin) || (upper(halfspace) isa Infinity))
        #     @test isapprox(pointon, point(lower(halfspace)), rtol = RTOLERANCE, atol = ATOLERANCE) || isapprox(pointon, point(upper(halfspace)), rtol = RTOLERANCE, atol = ATOLERANCE)
        # end

        # on xy
        cyl = Cylinder(0.5)
        r = Ray([1.0, 1.0, 0.0], [-1.0, -1.0, 0.0])
        intscts = surfaceintersection(cyl, r)
        xy = sqrt(2) / 4.0
        pt1 = [xy, xy, 0.0]
        pt2 = [-xy, -xy, 0.0]
        cpt1 = point(lower(intscts))
        cpt2 = point(upper(intscts))
        @test isapprox(pt1, cpt1, rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(pt2, cpt2, rtol = RTOLERANCE, atol = ATOLERANCE)

        # starting inside
        r = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
        intscts = surfaceintersection(cyl, r)
        @test lower(intscts) isa RayOrigin && isapprox(pt1, point(upper(intscts)), rtol = RTOLERANCE, atol = ATOLERANCE)

        # inside infinite
        r = Ray([0.0, 0.0, 0.0], [0.0, 0.0, 1.0])
        intscts = surfaceintersection(cyl, r)
        @test lower(intscts) isa RayOrigin && upper(intscts) isa Infinity

        # on surface
        r = Ray([0.5, 0.0, 0.0], [0.0, 0.0, 1.0])
        intscts = surfaceintersection(cyl, r)
        @test lower(intscts) isa RayOrigin && upper(intscts) isa Infinity

        # on diagonal
        r = Ray([1.0, 1.0, 1.0], [-1.0, -1.0, -1.0])
        intscts = surfaceintersection(cyl, r)
        cpt1 = point(lower(intscts))
        cpt2 = point(upper(intscts))
        pt1 = [xy, xy, xy]
        pt2 = [-xy, -xy, -xy]
        @test isapprox(pt1, cpt1, rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(pt2, cpt2, rtol = RTOLERANCE, atol = ATOLERANCE)

        # missing
        r = Ray([5.0, 5.0, 5.0], [0.0, 0.0, -1.0])
        intscts = surfaceintersection(cyl, r)
        @test intscts isa EmptyInterval
    end # testset cylinder

    @testset "Sphere" begin
        # Random.seed!(SEED)

        # # random point intersections
        # ur, vr = uvrange(Sphere)
        # for i in 1:10000
        #     sph = Sphere(rand() * 3)
        #     u, v = rand() * (ur[2] - ur[1]) + ur[1], rand() * (vr[2] - vr[1]) + vr[1]
        #     pointon = point(sph, u, v)
        #     origin = SVector(4.0, 4.0, 4.0)
        #     r = Ray(origin, pointon .- origin)
        #     halfspace = surfaceintersection(sph, r)
        #     @test halfspace !== nothing && !((lower(halfspace) isa RayOrigin) || (upper(halfspace) isa Infinity))
        #     @test isapprox(pointon, point(lower(halfspace)), rtol = RTOLERANCE, atol = ATOLERANCE) || isapprox(pointon, point(upper(halfspace)), rtol = RTOLERANCE, atol = ATOLERANCE)
        # end

        # on xy plane
        sph = Sphere(0.5)
        r = Ray([1.0, 1.0, 0.0], [-1.0, -1.0, 0.0])
        intscts = surfaceintersection(sph, r)
        xy = sqrt(2) / 4.0
        pt1 = [xy, xy, 0.0]
        pt2 = [-xy, -xy, 0.0]
        cpt1 = point(lower(intscts))
        cpt2 = point(upper(intscts))
        @test isapprox(pt1, cpt1, rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(pt2, cpt2, rtol = RTOLERANCE, atol = ATOLERANCE)

        # starting inside
        r = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
        intscts = surfaceintersection(sph, r)
        @test lower(intscts) isa RayOrigin && isapprox(pt1, point(upper(intscts)), rtol = RTOLERANCE, atol = ATOLERANCE)

        # on diagonal
        r = Ray([1.0, 1.0, 1.0], [-1.0, -1.0, -1.0])
        intscts = surfaceintersection(sph, r)
        cpt1 = point(lower(intscts))
        cpt2 = point(upper(intscts))
        xyz = sqrt(3) / 6.0
        pt1 = [xyz, xyz, xyz]
        pt2 = [-xyz, -xyz, -xyz]
        @test isapprox(pt1, cpt1, rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(pt2, cpt2, rtol = RTOLERANCE, atol = ATOLERANCE)

        # missing
        r = Ray([5.0, 5.0, 5.0], [0.0, 0.0, -1.0])
        intscts = surfaceintersection(sph, r)
        @test intscts isa EmptyInterval

        # tangent
        r = Ray([0.5, 0.5, 0.0], [0.0, -1.0, 0.0])
        intscts = surfaceintersection(sph, r)
        @test intscts isa EmptyInterval
    end # testset sphere

    @testset "Spherical Cap" begin
        @test_throws AssertionError SphericalCap(0.0, 1.0)
        @test_throws AssertionError SphericalCap(1.0, 0.0)

        sph = SphericalCap(0.5, π / 3, SVector(1.0, 1.0, 0.0), SVector(1.0, 1.0, 1.0))

        r = Ray([10.0, 10.0, 1.0], [-1.0, -1.0, 0.0])
        int = surfaceintersection(sph, r)
        @test isapprox(point(lower(int)), [1.0, 1.0, 1.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(int) isa Infinity

        r = Ray([0.8, 1.4, 1.2], [1.0, -1.0, 0.0])
        int = surfaceintersection(sph, r)
        @test isapprox(point(lower(int)), [0.898231126982741, 1.301768873017259, 1.2], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [1.301768873017259, 0.8982311269827411, 1.2], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([2.0, 2.0, 1.5], [-1.0, -1.0, 0.0])
        @test surfaceintersection(sph, r) isa EmptyInterval
    end # testset spherical cap

    @testset "Triangle" begin
        Random.seed!(SEED)

        tri = Triangle(SVector{3}(2.0, 1.0, 1.0), SVector{3}(1.0, 2.0, 1.0), SVector{3}(1.0, 1.0, 2.0), SVector{2}(0.0, 0.0), SVector{2}(1.0, 0.0), SVector{2}(0.5, 0.5))

        # random point intersections
        for i in 1:1000
            barycentriccoords = rand(3)
            barycentriccoords /= sum(barycentriccoords)
            pointontri = point(tri, barycentriccoords...)
            origin = SVector(3.0, 3.0, 3.0)
            r = Ray(origin, pointontri .- origin)
            halfspace = surfaceintersection(tri, r)
            @test !(halfspace isa EmptyInterval) && upper(halfspace) isa Infinity
            t = α(lower(halfspace))
            @test isapprox(point(r, t), point(lower(halfspace)), rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(pointontri, point(lower(halfspace)), rtol = RTOLERANCE, atol = ATOLERANCE)
        end

        # front and back face intersections
        tri = Triangle(SVector{3}(0.0, 0.0, 0.0), SVector{3}(1.0, 0.0, 0.0), SVector{3}(0.0, 1.0, 0.0), SVector{2}(0.0, 0.0), SVector{2}(1.0, 0.0), SVector{2}(0.0, 1.0))
        r1 = Ray([0.1, 0.1, 1.0], [0.0, 0.0, -1.0])
        r2 = Ray([0.1, 0.1, -1.0], [0.0, 0.0, 1.0])

        intsct1 = halfspaceintersection(surfaceintersection(tri, r1))
        intsct2 = halfspaceintersection(surfaceintersection(tri, r2))

        @test isapprox(point(intsct1), point(intsct2), rtol = RTOLERANCE, atol = ATOLERANCE)

        # ray should miss
        r = Ray([1.0, 1.0, 1.0], [0.0, 0.0, -1.0])
        intsct = surfaceintersection(tri, r)
        @test intsct isa EmptyInterval
    end # testset triangle

    @testset "Plane" begin
        rinplane = Ray([0.0, 0.0, 1.0], [1.0, 1.0, 0.0])
        routside = Ray([0.0, 0.0, 2.0], [1.0, 1.0, 1.0])
        rinside = Ray([0.0, 0.0, -1.0], [1.0, 1.0, -1.0])
        routintersects = Ray([0.0, 0.0, 2.0], [0.0, 0.0, -1.0])
        rinintersects = Ray([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
        pln = Plane(0.0, 0.0, 1.0, 0.0, 0.0, 1.0)

        # parallel to plane and outside
        @test surfaceintersection(pln, routside) isa EmptyInterval

        # NOTE for coplanar faces to work with visualization, rays in the plane must count as being 'inside'
        # parallel and on plane
        res = surfaceintersection(pln, rinplane)
        @test lower(res) isa RayOrigin && upper(res) isa Infinity

        # inside but not hitting
        res = surfaceintersection(pln, rinside)
        @test lower(res) isa RayOrigin && upper(res) isa Infinity

        # starts outside and hits
        res = surfaceintersection(pln, routintersects)
        @test isapprox(point(lower(res)), [0.0, 0.0, 1.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits
        res = surfaceintersection(pln, rinintersects)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.0, 0.0, 1.0], rtol = RTOLERANCE, atol = ATOLERANCE)
    end # testset plane

    @testset "Rectangle" begin
        @test_throws AssertionError Rectangle(0.0, 1.0)
        @test_throws AssertionError Rectangle(1.0, 0.0)

        rinplane = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
        routside = Ray([0.0, 0.0, 1.0], [1.0, 1.0, 1.0])
        rinside = Ray([0.0, 0.0, -1.0], [1.0, 1.0, -1.0])
        routintersects = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])
        rinintersects = Ray([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
        routmiss = Ray([0.0, 2.0, 1.0], [0.0, 0.0, -1.0])
        rinmiss = Ray([0.0, 2.0, -1.0], [0.0, 0.0, 1.0])
        rinbounds = Ray([0.5, 0.5, -1.0], [0.0, 0.0, 1.0])
        routbounds = Ray([0.5, 0.5, 1.0], [0.0, 0.0, -1.0])

        rect = Rectangle(0.5, 0.5)

        # parallel to plane and outside
        @test surfaceintersection(rect, routside) isa EmptyInterval

        # parallel and on plane
        @test surfaceintersection(rect, rinplane) isa EmptyInterval

        # inside but not hitting
        @test surfaceintersection(rect, rinside) isa EmptyInterval

        # starts outside and hits
        res = surfaceintersection(rect, routintersects)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits
        res = surfaceintersection(rect, rinintersects)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # starts inside and misses bounds
        @test surfaceintersection(rect, rinmiss) isa EmptyInterval

        # starts outside and misses bounds
        @test surfaceintersection(rect, routmiss) isa EmptyInterval

        # starts outside and hits bounds
        res = surfaceintersection(rect, routbounds)
        @test isapprox(point(lower(res)), [0.5, 0.5, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits bounds
        res = surfaceintersection(rect, rinbounds)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.5, 0.5, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # test a rectangle with translation and rotation
        rect2 = Rectangle(0.3, 0.5, SVector(3.0, 1.0, 4.0), SVector(0.2, 0.3, 0.4))
        rayhit = Ray([0.6, 0.6, 0.6], [-0.3, -0.1, -0.3])
        raymiss = Ray([0.7, 0.4, 0.4], [-0.3, -0.1, -0.3])
        res = surfaceintersection(rect2, rayhit)
        @test isapprox(point(lower(res)), [0.2863636363636363, 0.4954545454545454, 0.2863636363636363], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        @test surfaceintersection(rect2, raymiss) isa EmptyInterval
    end # testset Rectangle

    @testset "Ellipse" begin
        @test_throws AssertionError Ellipse(0.0, 1.0)
        @test_nowarn Circle(0.5)

        rinplane = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
        routside = Ray([0.0, 0.0, 1.0], [1.0, 1.0, 1.0])
        rinside = Ray([0.0, 0.0, -1.0], [1.0, 1.0, -1.0])
        routintersects = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])
        rinintersects = Ray([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
        routmiss = Ray([0.0, 2.0, 1.0], [0.0, 0.0, -1.0])
        rinmiss = Ray([0.0, 2.0, -1.0], [0.0, 0.0, 1.0])
        rinbounds = Ray([0.5, 0.0, -1.0], [0.0, 0.0, 1.0])
        routbounds = Ray([0.5, 0.0, 1.0], [0.0, 0.0, -1.0])
        rasymhit = Ray([0.0, 0.9, 1.0], [0.0, 0.0, -1.0])
        rasymmiss = Ray([0.9, 0.0, 1.0], [0.0, 0.0, -1.0])

        ell = Ellipse(0.5, 1.0)

        # parallel to plane and outside
        @test surfaceintersection(ell, routside) isa EmptyInterval

        # parallel and on plane
        @test surfaceintersection(ell, rinplane) isa EmptyInterval

        # inside but not hitting
        @test surfaceintersection(ell, rinside) isa EmptyInterval

        # starts outside and hits
        res = surfaceintersection(ell, routintersects)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits
        res = surfaceintersection(ell, rinintersects)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # starts inside and misses bounds
        @test surfaceintersection(ell, rinmiss) isa EmptyInterval

        # starts outside and misses bounds
        @test surfaceintersection(ell, routmiss) isa EmptyInterval

        # starts outside and hits bounds
        res = surfaceintersection(ell, routbounds)
        @test isapprox(point(lower(res)), [0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits bounds
        res = surfaceintersection(ell, rinbounds)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # asymmetric hit
        res = surfaceintersection(ell, rasymhit)
        @test isapprox(point(lower(res)), [0.0, 0.9, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # asymmetric miss
        @test surfaceintersection(ell, rasymmiss) isa EmptyInterval

        # test an ellipse with translation and rotation
        ell2 = Ellipse(0.3, 0.5, SVector(3.0, 1.0, 4.0), SVector(0.2, 0.3, 0.4))
        rayhit = Ray([0.6, 0.6, 0.6], [-0.3, -0.1, -0.3])
        raymiss = Ray([0.7, 0.4, 0.4], [-0.3, -0.1, -0.3])
        res = surfaceintersection(ell2, rayhit)
        @test isapprox(point(lower(res)), [0.2863636363636363, 0.4954545454545454, 0.2863636363636363], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        @test surfaceintersection(ell2, raymiss) isa EmptyInterval
    end # testset Ellipse

    @testset "Hexagon" begin
        @test_nowarn Hexagon(0.5)

        rinplane = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
        routside = Ray([0.0, 0.0, 1.0], [1.0, 1.0, 1.0])
        rinside = Ray([0.0, 0.0, -1.0], [1.0, 1.0, -1.0])
        routintersects = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])
        rinintersects = Ray([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
        routmiss = Ray([0.0, 2.0, 1.0], [0.0, 0.0, -1.0])
        rinmiss = Ray([0.0, 2.0, -1.0], [0.0, 0.0, 1.0])
        rinbounds = Ray([sqrt(3) / 2, 0.0, -1.0], [0.0, 0.0, 1.0])
        routbounds = Ray([sqrt(3) / 2, 0.0, 1.0], [0.0, 0.0, -1.0])
        rasymhit = Ray([0.0, 1.0, 1.0], [0.0, 0.0, -1.0])
        rasymmiss = Ray([1.0, 0.0, 1.0], [0.0, 0.0, -1.0])

        hex = Hexagon(1.0)

        # parallel to plane and outside
        @test surfaceintersection(hex, routside) isa EmptyInterval

        # parallel and on plane
        @test surfaceintersection(hex, rinplane) isa EmptyInterval

        # inside but not hitting
        @test surfaceintersection(hex, rinside) isa EmptyInterval

        # starts outside and hits
        res = surfaceintersection(hex, routintersects)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits
        res = surfaceintersection(hex, rinintersects)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # starts inside and misses bounds
        @test surfaceintersection(hex, rinmiss) isa EmptyInterval

        # starts outside and misses bounds
        @test surfaceintersection(hex, routmiss) isa EmptyInterval

        # starts outside and hits bounds
        res = surfaceintersection(hex, routbounds)
        @test isapprox(point(lower(res)), [sqrt(3) / 2, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # starts inside and hits bounds
        res = surfaceintersection(hex, rinbounds)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [sqrt(3) / 2, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # asymmetric hit
        res = surfaceintersection(hex, rasymhit)
        @test isapprox(point(lower(res)), [0.0, 1.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # asymmetric miss
        @test surfaceintersection(hex, rasymmiss) isa EmptyInterval

        # test an ellipse with translation and rotation
        hex2 = Hexagon(0.4, SVector(3.0, 1.0, 4.0), SVector(0.2, 0.3, 0.4))
        rayhit = Ray([0.6, 0.6, 0.6], [-0.3, -0.1, -0.3])
        raymiss = Ray([0.7, 0.4, 0.4], [-0.3, -0.1, -0.3])
        res = surfaceintersection(hex2, rayhit)
        @test isapprox(point(lower(res)), [0.2863636363636363, 0.4954545454545454, 0.2863636363636363], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        @test surfaceintersection(hex2, raymiss) isa EmptyInterval
    end # testset Hexagon

    @testset "Stops" begin
        infiniterect = RectangularAperture(0.4, 0.8, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))
        finiterect = RectangularAperture(0.4, 0.8, 1.0, 2.0, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))
        # through hole
        r = Ray([0.0, 1.0, 1.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(infiniterect, r) isa EmptyInterval
        @test surfaceintersection(finiterect, r) isa EmptyInterval
        r = Ray([0.0, -1.0, 1.0], [0.0, 1.0, 0.0])
        @test surfaceintersection(infiniterect, r) isa EmptyInterval
        @test surfaceintersection(finiterect, r) isa EmptyInterval
        # on edge
        r = Ray([0.0, 1.0, 0.6], [0.0, -1.0, 0.0])
        res = surfaceintersection(infiniterect, r)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.6], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        res = surfaceintersection(finiterect, r)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.6], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        # through hole asym
        r = Ray([0.7, 1.0, 1.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(infiniterect, r) isa EmptyInterval
        @test surfaceintersection(finiterect, r) isa EmptyInterval
        # on finite edge
        r = Ray([1.0, -1.0, 2.0], [0.0, 1.0, 0.0])
        res = surfaceintersection(infiniterect, r)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [1.0, 0.0, 2.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        res = surfaceintersection(finiterect, r)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [1.0, 0.0, 2.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        # outside finite
        r = Ray([2.0, 1.0, 3.0], [0.0, -1.0, 0.0])
        @test isapprox(point(surfaceintersection(infiniterect, r)), [2.0, 0.0, 3.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test surfaceintersection(finiterect, r) isa EmptyInterval

        infinitecirc = CircularAperture(0.4, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))
        finitecirc = CircularAperture(0.4, 1.0, 2.0, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))
        # through hole
        r = Ray([0.0, 1.0, 1.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(infinitecirc, r) isa EmptyInterval
        @test surfaceintersection(finitecirc, r) isa EmptyInterval
        r = Ray([0.0, -1.0, 1.0], [0.0, 1.0, 0.0])
        @test surfaceintersection(infinitecirc, r) isa EmptyInterval
        @test surfaceintersection(finitecirc, r) isa EmptyInterval
        # on edge
        r = Ray([0.0, 1.0, 0.6], [0.0, -1.0, 0.0])
        res = surfaceintersection(infinitecirc, r)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.6], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        res = surfaceintersection(finitecirc, r)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.6], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        # through hole 2
        r = Ray([0.3, 1.0, 1.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(infinitecirc, r) isa EmptyInterval
        @test surfaceintersection(finitecirc, r) isa EmptyInterval
        # on finite edge
        r = Ray([1.0, -1.0, 2.0], [0.0, 1.0, 0.0])
        res = surfaceintersection(infinitecirc, r)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [1.0, 0.0, 2.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        res = surfaceintersection(finitecirc, r)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [1.0, 0.0, 2.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        # outside finite
        r = Ray([2.0, 1.0, 3.0], [0.0, -1.0, 0.0])
        @test isapprox(point(surfaceintersection(infinitecirc, r)), [2.0, 0.0, 3.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test surfaceintersection(finitecirc, r) isa EmptyInterval

        annulus = Annulus(0.5, 1.0, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))
        # through hole
        r = Ray([0.0, 1.0, 1.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(annulus, r) isa EmptyInterval
        r = Ray([0.0, -1.0, 1.0], [0.0, 1.0, 0.0])
        @test surfaceintersection(annulus, r) isa EmptyInterval
        # on edge
        r = Ray([0.0, 1.0, 0.5], [0.0, -1.0, 0.0])
        res = surfaceintersection(annulus, r)
        @test isapprox(point(lower(res)), [0.0, 0.0, 0.5], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
        # through hole 2
        r = Ray([0.3, 1.0, 1.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(infiniterect, r) isa EmptyInterval
        @test surfaceintersection(annulus, r) isa EmptyInterval
        # on finite edge
        r = Ray([1.0, -1.0, 1.0], [0.0, 1.0, 0.0])
        res = surfaceintersection(annulus, r)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [1.0, 0.0, 1.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        # outside finite
        r = Ray([0.9, 1.0, 1.9], [0.0, -1.0, 0.0])
        @test surfaceintersection(annulus, r) isa EmptyInterval
    end # testset Stops

    @testset "Bezier" begin
        surf = TestData.beziersurface()
        accelsurf = AcceleratedParametricSurface(surf)
        numsamples = 100

        samples = samplepoints(numsamples, 0.05, 0.95, 0.05, 0.95)
        missedintersections = 0

        for uv in samples
            pt = point(surf, uv[1], uv[2])
            origin = [0.5, 0.5, 5.0]
            r = Ray(origin, pt .- origin)
            allintersections = surfaceintersection(accelsurf, r)
            if allintersections isa EmptyInterval
                # sampled triangle acceleration structure doesn't guarantee all intersecting rays will be detected.
                # Have to use patch subdivision and convex hull polyhedron to do this.
                missedintersections += 1
            else
                if isa(allintersections, Interval)
                    allintersections = (allintersections,)
                end

                for intsct in allintersections
                    @test isapprox(point(halfspaceintersection(intsct)), pt, rtol = RTOLERANCE, atol = ATOLERANCE)
                    @test upper(intsct) isa Infinity
                end
            end
        end

        if missedintersections > 0
            @warn "$missedintersections out of total $(length(samples)) bezier suface intersections were missed"
        end

        # miss from outside
        r = Ray([0.5, 0.5, 5.0], [0.0, 0.0, 1.0])
        @test surfaceintersection(accelsurf, r) isa EmptyInterval
        r = Ray([5.0, 0.5, -5.0], [0.0, 0.0, -1.0])
        @test surfaceintersection(accelsurf, r) isa EmptyInterval

        # miss from inside
        r = Ray([0.5, 0.5, -5.0], [0.0, 0.0, -1.0])
        res = surfaceintersection(accelsurf, r)
        # TODO!! Fix bezier surface to create halfspace
        # @test lower(res) isa RayOrigin && upper(res) isa Infinity

        # hit from inside
        r = Ray([0.5, 0.5, 0.2], [0.0, 1.0, 0.0])
        res = surfaceintersection(accelsurf, r)
        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.5, 0.8996900152986361, 0.2], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit from outside
        r = Ray([0.0, 0.5, 0.2], [1.0, 0.0, -1.0])
        res = surfaceintersection(accelsurf, r)
        @test isapprox(point(lower(res)), [0.06398711204047353, 0.5, 0.13601288795952649], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity

        # two hits from outside
        r = Ray([0.0, 0.5, 0.25], [1.0, 0.0, 0.0])
        res = surfaceintersection(accelsurf, r)
        @test isapprox(point(lower(res)), [0.12607777053030267, 0.5, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res)), [0.8705887077060419, 0.5, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)

        surf = TestData.wavybeziersurface()
        accelsurf = AcceleratedParametricSurface(surf)

        # 3 hit starting outside
        r = Ray([0.5, 0.0, 0.0], [0.0, 1.0, 0.0])
        res = surfaceintersection(accelsurf, r)
        @test isa(res, DisjointUnion) && length(res) == 2 && isapprox(point(lower(res[1])), [0.5, 0.10013694786182059, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[1])), [0.5, 0.49625, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(lower(res[2])), [0.5, 0.8971357794109067, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && (upper(res[2]) isa Infinity)

        # 3 hit starting inside
        r = Ray([0.5, 1.0, 0.0], [0.0, -1.0, 0.0])
        res = surfaceintersection(accelsurf, r)
        @test isa(res, DisjointUnion) && length(res) == 2 && (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.5, 0.8971357794109067, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(lower(res[2])), [0.5, 0.49625, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.5, 0.10013694786182059, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        surf = TestData.verywavybeziersurface()
        accelsurf = AcceleratedParametricSurface(surf, 20)

        # five hits starting outside
        r = Ray([0.9, 0.0, -0.3], [0.0, 1.0, 0.7])
        res = surfaceintersection(accelsurf, r)
        a = isapprox(point(lower(res[1])), [0.9, 0.03172286522032046, -0.2777939943457758], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[1])), [0.9, 0.1733979947040411, -0.17862140370717122], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(res[2])), [0.9, 0.5335760974594397, 0.07350326822160776], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.7767707607392784, 0.24373953251749486], rtol = RTOLERANCE, atol = ATOLERANCE)
        c = isapprox(point(lower(res[3])), [0.9, 0.9830891958374246, 0.3881624370861975], rtol = RTOLERANCE, atol = ATOLERANCE) && (upper(res[3]) isa Infinity)
        @test isa(res, DisjointUnion) && length(res) == 3 && a && b && c

        # five hits starting inside
        r = Ray([0.9, 1.0, 0.4], [0.0, -1.0, -0.7])
        res = surfaceintersection(accelsurf, r)
        a = (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.9, 0.9830891958374246, 0.3881624370861975], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(res[2])), [0.9, 0.7767707607392784, 0.24373953251749486], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.5335760974594397, 0.07350326822160776], rtol = RTOLERANCE, atol = ATOLERANCE)
        c = isapprox(point(lower(res[3])), [0.9, 0.17339799470404108, -0.1786214037071712], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[3])), [0.9, 0.03172286522032046, -0.27779399434577573], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(res, DisjointUnion) && length(res) == 3 && a && b && c

        # 4 hits starting inside
        r = Ray([0.1, 0.0, -0.3], [0.0, 1.0, 0.7])
        res = surfaceintersection(accelsurf, r)
        a = (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.1, 0.2851860296285551, -0.10036977926001144], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(res[2])), [0.1, 0.5166793625025807, 0.06167555375180668], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.1, 0.7770862508789854, 0.24396037561528983], rtol = RTOLERANCE, atol = ATOLERANCE)
        c = isapprox(point(lower(res[3])), [0.1, 0.98308919558696, 0.3881624369108719], rtol = RTOLERANCE, atol = ATOLERANCE) && (upper(res[3]) isa Infinity)
        @test isa(res, DisjointUnion) && length(res) == 3 && a && b && c

        # 4 hits starting outside
        r = Ray([0.9, 0.9, 0.4], [0.0, -0.9, -0.7])
        res = surfaceintersection(accelsurf, r)
        a = isapprox(point(lower(res[1])), [0.9, 0.736072142615238, 0.2725005553674076], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[1])), [0.9, 0.567439326091764, 0.141341698071372], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(res[2])), [0.9, 0.16601081959179267, -0.1708804736508277], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.032434058775915924, -0.274773509840954], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(res, DisjointUnion) && length(res) == 2 && a && b
    end # testset Bezier

    @testset "Zernike" begin
        @test_nowarn ZernikeSurface(1.5)
        @test_throws AssertionError ZernikeSurface(0.0)
        @test_throws AssertionError ZernikeSurface(1.5, aspherics = [(1, 0.1)])

        z1 = TestData.zernikesurface1()
        az1 = AcceleratedParametricSurface(z1, 20)
        numsamples = 100

        # test that an on axis ray doesn't miss
        @test !(surfaceintersection(AcceleratedParametricSurface(z1), Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])) isa EmptyInterval)

        samples = samplepoints(numsamples, 0.01, 0.99, 0.01, 2π - 0.01)
        missedintersections = 0

        for uv in samples
            pt = point(z1, uv[1], uv[2])
            origin = [0.0, 0.0, 5.0]
            r = Ray(origin, pt .- origin)
            allintersections = surfaceintersection(az1, r)
            if (allintersections isa EmptyInterval)
                # sampled triangle acceleration structure doesn't guarantee all intersecting rays will be detected.
                # Have to use patch subdivision and convex hull polyhedron to do this.
                missedintersections += 1
            else
                if isa(allintersections, Interval)
                    allintersections = (allintersections,)
                end

                for intsct in allintersections
                    if lower(intsct) isa RayOrigin
                        missedintersections += 1
                    else
                        # closest point should be on the surface, furthest should be on bounding prism
                        @test isapprox(point(lower(intsct)), pt, rtol = RTOLERANCE, atol = ATOLERANCE)
                        @test isa(upper(intsct), Intersection{Float64,3}) || (OpticSim.direction(r)[3] == -1 && isa(upper(intsct), Infinity{Float64}))
                    end
                end
            end
        end

        if missedintersections > 0
            @warn "$missedintersections out of total $(length(samples)) zernike suface intersections were missed"
        end

        # hit from inside
        r = Ray([0.0, 0.0, -0.5], [0.1, 0.2, 0.5])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.12363711269619936, 0.24727422539239863, 0.11818556348099664], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit from outside
        r = Ray([0.0, 0.0, 0.5], [0.1, 0.2, -0.5])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [0.07118311042701463, 0.1423662208540292, 0.144084447864927], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [0.6708203932499369, 1.3416407864998738, -2.8541019662496843], rtol = RTOLERANCE, atol = ATOLERANCE)

        # miss from inside
        r = Ray([0.2, 0.2, -0.5], [0.0, 0.0, -1.0])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && (upper(intvl) isa Infinity)

        # miss from outside
        r = Ray([0.2, 2.0, -0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 2.0, 0.5], [0.0, 0.0, -1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 0.2, 0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 2.0, 2.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(az1, r) isa EmptyInterval

        # two hits from outside
        r = Ray([2.0, 0.0, 0.25], [-1.0, 0.0, 0.0])
        hit = surfaceintersection(az1, r)
        a = isapprox(point(lower(hit[1])), [1.5, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[1])), [1.3571758210851095, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(hit[2])), [-1.2665828947521165, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[2])), [-1.5, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(hit, DisjointUnion) && length(hit) == 2 && a && b

        # hit cyl from outside
        r = Ray([0.0, 2.0, -0.5], [0.0, -1.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [0.0, 1.5, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [0.0, -1.5, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit cyl from inside
        r = Ray([0.0, 0.0, -0.5], [0.0, -1.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.0, -1.5, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE)

        z2 = TestData.zernikesurface2()
        az2 = AcceleratedParametricSurface(z2, 20)

        # two hits
        r = Ray([2.0, -1.0, 0.2], [-1.0, 0.0, 0.0])
        intvl = surfaceintersection(az2, r)
        @test isapprox(point(lower(intvl)), [0.238496383738369, -1.0, 0.2], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [-0.8790518227874484, -1.0, 0.2], rtol = RTOLERANCE, atol = ATOLERANCE)

        # three hits
        r = Ray([-0.7, 1.0, 0.2], [0.0, -1.0, 0.0])
        hit = surfaceintersection(az2, r)
        a = (lower(hit[1]) isa RayOrigin) && isapprox(point(upper(hit[1])), [-0.7, 0.8732489598020176, 0.2], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(hit[2])), [-0.7, -0.34532502965048606, 0.2], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[2])), [-0.7, -1.1615928003236047, 0.2], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(hit, DisjointUnion) && length(hit) == 2 && a && b

        # four hits
        r = Ray([1.5, 0.1, 0.35], [-1.0, -0.5, 0.0])
        hit = surfaceintersection(az2, r)
        a = isapprox(point(lower(hit[1])), [1.4371467631969683, 0.06857338159848421, 0.35], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[1])), [1.1428338578611368, -0.07858307106943167, 0.35], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(hit[2])), [0.12916355683491865, -0.5854182215825409, 0.35], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[2])), [-0.7290713614371003, -1.0145356807185504, 0.35], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(hit, DisjointUnion) && length(hit) == 2 && a && b

        # failure case, had a bug where rounding error would cause o2 to fail
        o1 = leaf(AcceleratedParametricSurface(ZernikeSurface(1.4 * 1.15)), translation(0.0, 0.0, 3.0))()
        o2 = leaf(AcceleratedParametricSurface(ZernikeSurface(1.4 * 1.15)), translation(2.0, 0.0, 3.0))()
        r = Ray([-5.0, 0.0, 0.0], [1.0, 0.0, 0.0])
        @test !(surfaceintersection(o1, r) isa EmptyInterval)
        @test !(surfaceintersection(o2, r) isa EmptyInterval)
    end # testset Zernike

    @testset "QType" begin
        @test_nowarn QTypeSurface(1.5)
        @test_throws AssertionError QTypeSurface(0.0)

        q1 = TestData.qtypesurface1()
        aq1 = AcceleratedParametricSurface(q1, 20)
        numsamples = 100

        # test that an on axis ray doesn't miss
        @test !(surfaceintersection(AcceleratedParametricSurface(q1), Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])) isa EmptyInterval)

        samples = samplepoints(numsamples, 0.01, 0.99, 0.01, 2π - 0.01)
        missedintersections = 0

        for uv in samples
            pt = point(q1, uv[1], uv[2])
            origin = [0.0, 0.0, 5.0]
            r = Ray(origin, pt .- origin)
            allintersections = surfaceintersection(aq1, r)
            if (allintersections isa EmptyInterval)
                # sampled triangle acceleration structure doesn't guarantee all intersecting rays will be detected.
                # Have to use patch subdivision and convex hull polyhedron to do this.
                missedintersections += 1
            else
                if isa(allintersections, Interval)
                    allintersections = (allintersections,)
                end

                for intsct in allintersections
                    if lower(intsct) isa RayOrigin
                        missedintersections += 1
                    else
                        # closest point should be on the surface, furthest should be on bounding prism
                        @test isapprox(point(lower(intsct)), pt, rtol = RTOLERANCE, atol = ATOLERANCE)
                        @test isa(upper(intsct), Intersection{Float64,3}) || (OpticSim.direction(r)[3] == -1 && isa(upper(intsct), Infinity{Float64}))
                    end
                end
            end
        end

        if missedintersections > 0
            @warn "$missedintersections out of total $(length(samples)) qtype suface intersections were missed"
        end

        # hit from inside
        r = Ray([0.0, 0.0, -0.5], [0.1, 0.2, 0.5])
        intvl = surfaceintersection(aq1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.09758258750208074, 0.19516517500416142, -0.01208706248959643], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit from outside
        r = Ray([0.0, 0.0, 0.5], [0.1, 0.2, -0.5])
        intvl = surfaceintersection(aq1, r)
        @test isapprox(point(lower(intvl)), [0.10265857811957124, 0.20531715623914257, -0.013292890597856405], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [0.6708203932499369, 1.3416407864998738, -2.8541019662496843], rtol = RTOLERANCE, atol = ATOLERANCE)

        # miss from inside
        r = Ray([0.2, 0.2, -0.5], [0.0, 0.0, -1.0])
        intvl = surfaceintersection(aq1, r)
        @test (lower(intvl) isa RayOrigin) && (upper(intvl) isa Infinity)

        # miss from outside
        r = Ray([0.2, 2.0, -0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(aq1, r) isa EmptyInterval
        r = Ray([0.2, 2.0, 0.5], [0.0, 0.0, -1.0])
        @test surfaceintersection(aq1, r) isa EmptyInterval
        r = Ray([0.2, 0.2, 0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(aq1, r) isa EmptyInterval
        r = Ray([0.2, 2.0, 2.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(aq1, r) isa EmptyInterval
    end # testset QType

    @testset "BoundingBox" begin
        function boundingval(a::BoundingBox{T}, axis::Int, plane::Bool) where {T<:Real}
            if axis === 1
                return plane ? a.xmax : a.xmin
            elseif axis === 2
                return plane ? a.ymax : a.ymin
            elseif axis == 3
                return plane ? a.zmax : a.zmin
            else
                throw(ErrorException("Invalid axis: $axis"))
            end
        end

        Random.seed!(SEED)

        axis(x) = (x + 1) ÷ 2
        face(x) = mod(x, 2)
        facevalue(a::BoundingBox, x) = boundingval(a, axis(x), Bool(face(x)))
        boundingindices = ((2, 3), (1, 3), (1, 2))

        function facepoint(a::BoundingBox, facenumber)
            axisindex = axis(facenumber)
            indices = boundingindices[axisindex]
            b1min, b1max = boundingval(a, indices[1], false), boundingval(a, indices[1], true)
            b2min, b2max = boundingval(a, indices[2], false), boundingval(a, indices[2], true)

            step = 0.00001
            pt1val = rand((b1min + step):step:(b1max - step))
            pt2val = rand((b2min + step):step:(b2max - step))

            result = Array{Float64,1}(undef, 3)
            result[indices[1]] = pt1val
            result[indices[2]] = pt2val
            result[axisindex] = facevalue(a, facenumber)
            return result
        end

        bbox = BoundingBox(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)

        for i in 1:10000
            # randomly pick two planes
            face1 = rand(1:6)
            face2 = rand(1:6)
            while face2 == face1
                face2 = rand(1:6)
            end
            # pick a point on each face
            pt1 = facepoint(bbox, face1)
            pt2 = facepoint(bbox, face2)

            direction = normalize(pt1 - pt2)
            origin = pt2 - 3 * direction

            r = Ray(origin, direction)

            @test doesintersect(bbox, r)
        end

        # should miss
        @test !doesintersect(bbox, Ray([-2.0, -2.0, 0.0], [1.0, 0.0, 0.0]))
        @test !doesintersect(bbox, Ray([-2.0, -2.0, 0.0], [0.0, 1.0, 0.0]))
        @test !doesintersect(bbox, Ray([-2.0, -2.0, 0.0], [0.0, 0.0, 1.0]))
        @test !doesintersect(bbox, Ray([-2.0, 0.0, 0.0], [-1.0, 0.0, 0.0]))
        @test !doesintersect(bbox, Ray([0.0, -2.0, 0.0], [0.0, -1.0, 0.0]))
        @test !doesintersect(bbox, Ray([0.0, 0.0, -2.0], [0.0, 0.0, -1.0]))

        bbox = BoundingBox(-0.5, 0.5, -0.75, 0.75, typemin(Float64), typemax(Float64))

        # starting outside
        r = Ray([-1.0, -1.1, 0.0], [1.0, 1.0, 0.0])
        intscts = surfaceintersection(bbox, r)
        @test isapprox(point(lower(intscts)), [-0.5, -0.6, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intscts)), [0.5, 0.4, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # starting inside
        r = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
        intscts = surfaceintersection(bbox, r)
        @test lower(intscts) isa RayOrigin && isapprox(point(upper(intscts)), [0.5, 0.5, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # inside infinite
        r = Ray([0.0, 0.0, 0.0], [0.0, 0.0, 1.0])
        intscts = surfaceintersection(bbox, r)
        @test lower(intscts) isa RayOrigin && upper(intscts) isa Infinity

        # on surface
        r = Ray([0.5, 0.0, 0.0], [0.0, 0.0, 1.0])
        intscts = surfaceintersection(bbox, r)
        @test lower(intscts) isa RayOrigin && upper(intscts) isa Infinity

        # missing
        r = Ray([5.0, 5.0, 5.0], [0.0, 0.0, -1.0])
        intscts = surfaceintersection(bbox, r)
        @test intscts isa EmptyInterval

        r = Ray([5.0, 5.0, 5.0], [0.0, -1.0, 0.0])
        intscts = surfaceintersection(bbox, r)
        @test intscts isa EmptyInterval
    end # testset BoundingBox

    @testset "CSG" begin
        pln1 = Plane([0.0, 0.0, -1.0], [0.0, 0.0, -1.0])
        pln2 = Plane([0.0, 0.0, 1.0], [0.0, 0.0, 1.0])
        r = Ray([0.0, 0, 2.0], [0.0, 0.0, -1.0])
        gen1 = pln1 ∩ pln2
        gen2 = pln1 ∪ pln2

        csg = gen1(identitytransform())
        intsct = evalcsg(csg, r)

        intvlpt1 = lower(intsct)
        intvlpt2 = upper(intsct)

        @test isapprox(point(intvlpt1), SVector{3}(0.0, 0.0, 1.0)) && isapprox(point(intvlpt2), SVector{3}(0.0, 0.0, -1.0))

        csg = gen2(identitytransform())
        intsct = evalcsg(csg, r)
        @test (lower(intsct) isa RayOrigin) && (upper(intsct) isa Infinity)

        # test simple rays for intersection, union and difference operations
        # INTERSECTION
        intersection_obj = (leaf(Cylinder(0.5, 3.0), OpticSim.rotationd(90.0, 0.0, 0.0)) ∩ Sphere(1.0))()
        r = Ray([0.7, 0.0, 0.0], [-1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(intersection_obj, r)
        @test isapprox(point(lower(int)), [0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [-0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([5.0, 0.0, 0.0], [-1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(intersection_obj, r)
        @test isapprox(point(lower(int)), [0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [-0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([0.7, 0.0, 0.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(intersection_obj, r)
        @test (int isa EmptyInterval)

        r = Ray([0.0, 0.0, 0.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(intersection_obj, r)
        @test (lower(int) isa RayOrigin) && isapprox(point(upper(int)), [0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([0.0, 0.0, 0.0], [0.0, 1.0, 0.0])
        int = OpticSim.evalcsg(intersection_obj, r)
        @test (lower(int) isa RayOrigin) && isapprox(point(upper(int)), [0.0, 1.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([0.0, 2.0, 0.0], [0.0, -1.0, 0.0])
        int = OpticSim.evalcsg(intersection_obj, r)
        @test isapprox(point(lower(int)), [0.0, 1.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [0.0, -1.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # UNION
        union_obj = (Cylinder(0.5, 3.0) ∪ Sphere(1.0))(OpticSim.translation(1.0, 0.0, 0.0))
        r = Ray([1.0, 0.0, 0.0], [0.0, 0.0, 1.0])
        int = OpticSim.evalcsg(union_obj, r)
        @test (lower(int) isa RayOrigin) && (upper(int) isa Infinity)

        r = Ray([1.0, 0.0, 0.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(union_obj, r)
        @test (lower(int) isa RayOrigin) && isapprox(point(upper(int)), [2.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([1.6, 0.0, 0.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(union_obj, r)
        @test (lower(int) isa RayOrigin) && isapprox(point(upper(int)), [2.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([1.6, 0.0, 0.0], [-1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(union_obj, r)
        @test (lower(int) isa RayOrigin) && isapprox(point(upper(int)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([-1.0, 0.0, 0.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(union_obj, r)
        @test isapprox(point(lower(int)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [2.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([-1.0, 0.0, 4.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(union_obj, r)
        @test isapprox(point(lower(int)), [0.5, 0.0, 4.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [1.5, 0.0, 4.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        # DIFFERENCE
        difference_obj = (Cylinder(0.5, 3.0) - leaf(Sphere(1.0), OpticSim.translation(0.75, 0.0, 0.2)))()
        r = Ray([0.25, 0.0, 0.0], [-1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(difference_obj, r)
        @test isapprox(point(lower(int)), [-0.2297958971132712, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(int)), [-0.5, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([0.25, 0.0, 0.0], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(difference_obj, r)
        @test int isa EmptyInterval

        r = Ray([0.25, 0.0, 0.0], [0.0, 0.0, 1.0])
        int = OpticSim.evalcsg(difference_obj, r)
        @test isapprox(point(lower(int)), [0.25, 0.0, 1.0660254037844386], rtol = RTOLERANCE, atol = ATOLERANCE) && (upper(int) isa Infinity)

        r = Ray([0.25, 0.0, 1.5], [0.0, 0.0, -1.0])
        int = OpticSim.evalcsg(difference_obj, r)
        @test (lower(int[1]) isa RayOrigin) && isapprox(point(upper(int[1])), [0.25, 0.0, 1.0660254037844386], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(lower(int[2])), [0.25, 0.0, -0.6660254037844386], rtol = RTOLERANCE, atol = ATOLERANCE) && (upper(int[2]) isa Infinity)

        r = Ray([0.25, 0.0, 1.5], [1.0, 0.0, 0.0])
        int = OpticSim.evalcsg(difference_obj, r)
        @test (lower(int) isa RayOrigin) && isapprox(point(upper(int)), [0.5, 0.0, 1.5], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([0.25, 0.0, 1.5], [0.0, 0.0, 1.0])
        int = OpticSim.evalcsg(difference_obj, r)
        @test (lower(int) isa RayOrigin) && (upper(int) isa Infinity)

        # DisjointUnion result on CSG
        surf = TestData.verywavybeziersurface()
        accelsurf = leaf(AcceleratedParametricSurface(surf, 20))()

        # five hits starting outside
        r = Ray([0.9, 0.0, -0.3], [0.0, 1.0, 0.7])
        res = surfaceintersection(accelsurf, r)
        a = isapprox(point(lower(res[1])), [0.9, 0.03172286522032046, -0.2777939943457758], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[1])), [0.9, 0.1733979947040411, -0.17862140370717122], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(res[2])), [0.9, 0.5335760974594397, 0.07350326822160776], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.7767707607392784, 0.24373953251749486], rtol = RTOLERANCE, atol = ATOLERANCE)
        c = isapprox(point(lower(res[3])), [0.9, 0.9830891958374246, 0.3881624370861975], rtol = RTOLERANCE, atol = ATOLERANCE) && (upper(res[3]) isa Infinity)
        @test isa(res, DisjointUnion) && length(res) == 3 && a && b && c

        # five hits starting inside
        r = Ray([0.9, 1.0, 0.4], [0.0, -1.0, -0.7])
        res = surfaceintersection(accelsurf, r)
        a = (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.9, 0.9830891958374246, 0.3881624370861975], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(res[2])), [0.9, 0.7767707607392784, 0.24373953251749486], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.5335760974594397, 0.07350326822160776], rtol = RTOLERANCE, atol = ATOLERANCE)
        c = isapprox(point(lower(res[3])), [0.9, 0.17339799470404108, -0.1786214037071712], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(res[3])), [0.9, 0.03172286522032046, -0.27779399434577573], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(res, DisjointUnion) && length(res) == 3 && a && b && c
    end # testset CSG

    @testset "ThinGrating" begin
        angle_from_ray(raydirection) = 90 + atand(raydirection[3], raydirection[2])
        true_diff(order, λ, period, θi) = asind((order * λ / period + sind(θi)))
        period = 3.0
        int = TestData.transmissivethingrating(period, 2)
        for k in 1:50
            for θi in [-5.0, 0.0, 5.0, 10.0]
                for λ in [0.35, 0.55, 1.0]
                    ray = OpticalRay(SVector(0.0, 0.0, 2.0), SVector(0.0, sind(θi), -cosd(θi)), 1.0, λ)
                    raydir, _, _ = OpticSim.processintersection(int, SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 1.0), ray, 20.0, 1.0, true, true)
                    # order is random so just check that it is correct for one of them
                    @test any([isapprox(x, angle_from_ray(raydir), rtol = RTOLERANCE, atol = ATOLERANCE) for x in [true_diff(m, λ, period, θi) for m in -2:2]])
                end
            end
        end
    end # testset ThinGrating

    # TODO Hologram intersection tests

    @testset "Chebyshev" begin
        @test_throws AssertionError ChebyshevSurface(0.0, 1.0, [(1, 2, 1.0)])
        @test_throws AssertionError ChebyshevSurface(1.0, 0.0, [(1, 2, 1.0)])

        z1 = TestData.chebyshevsurface2()
        az1 = AcceleratedParametricSurface(z1, 20)
        numsamples = 100

        # test that an on axis ray doesn't miss
        @test !(surfaceintersection(AcceleratedParametricSurface(z1), Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])) isa EmptyInterval)

        samples = samplepoints(numsamples, -0.95, 0.95, -0.95, 0.95)
        missedintersections = 0

        for uv in samples
            pt = point(z1, uv[1], uv[2])
            origin = [0.0, 0.0, 5.0]
            r = Ray(origin, pt .- origin)
            allintersections = surfaceintersection(az1, r)
            if (allintersections isa EmptyInterval)
                # sampled triangle acceleration structure doesn't guarantee all intersecting rays will be detected.
                # Have to use patch subdivision and convex hull polyhedron to do this.
                missedintersections += 1
            else
                if isa(allintersections, Interval)
                    allintersections = (allintersections,)
                end

                for intsct in allintersections
                    if lower(intsct) isa RayOrigin
                        missedintersections += 1
                    else
                        # closest point should be on the surface, furthest should be on bounding prism
                        @test isapprox(point(lower(intsct)), pt, rtol = RTOLERANCE, atol = ATOLERANCE)
                        @test isa(upper(intsct), Intersection{Float64,3}) || (OpticSim.direction(r)[3] == -1 && isa(upper(intsct), Infinity{Float64}))
                    end
                end
            end
        end

        if missedintersections > 0
            @warn "$missedintersections out of total $(length(samples)) chebychev suface intersections were missed"
        end

        z1 = TestData.chebyshevsurface1()
        az1 = AcceleratedParametricSurface(z1, 20)

        # hit from inside
        r = Ray([0.0, 0.0, -0.5], [0.1, 0.2, 0.5])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.07870079995991423, 0.15740159991982847, -0.10649600020042879], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit from outside
        r = Ray([0.0, 0.0, 0.5], [0.1, 0.2, -0.5])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [0.13584702907541094, 0.2716940581508219, -0.1792351453770547], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [1.0, 2.0, -4.5], rtol = RTOLERANCE, atol = ATOLERANCE)

        # miss from inside
        r = Ray([0.2, 0.2, -0.5], [0.0, 0.0, -1.0])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && (upper(intvl) isa Infinity)

        # miss from outside
        r = Ray([0.2, 3.0, -0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 3.0, 0.5], [0.0, 0.0, -1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 0.2, 0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 3.0, 2.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(az1, r) isa EmptyInterval

        # two hits from outside
        r = Ray([0.0, -1.5, 0.25], [-1.0, 0.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [-1.030882920068228, -1.5, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [-1.786071230727613, -1.5, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)

        r = Ray([1.5, -0.8, 0.25], [-1.0, 0.0, 0.0])
        hit = surfaceintersection(az1, r)
        a = isapprox(point(lower(hit[1])), [0.21526832427065165, -0.8, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[1])), [-0.25523748548950076, -0.8, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(hit[2])), [-1.9273057166986296, -0.8, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[2])), [-2.0, -0.8, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(hit, DisjointUnion) && length(hit) == 2 && a && b

        # hit cyl from outside
        r = Ray([0.0, 3.0, -0.5], [0.0, -1.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [0.0, 2.0, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [0.0, -2.0, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit cyl from inside
        r = Ray([0.0, 0.0, -0.5], [0.0, -1.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.0, -2.0, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE)

    end # testset Chebyshev

    @testset "GridSag" begin
        z1 = TestData.gridsagsurfacebicubic()
        az1 = AcceleratedParametricSurface(z1, 20)

        numsamples = 100

        # test that an on axis ray doesn't miss
        @test !(surfaceintersection(AcceleratedParametricSurface(z1), Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])) isa EmptyInterval)

        samples1 = samplepoints(numsamples, 0.01, 0.99, 0.01, 2π - 0.01)
        missedintersections = 0

        for uv in samples1
            pt = point(z1, uv[1], uv[2])
            origin = [0.0, 0.0, 5.0]
            r = Ray(origin, pt .- origin)
            allintersections = surfaceintersection(az1, r)
            if (allintersections isa EmptyInterval)
                # sampled triangle acceleration structure doesn't guarantee all intersecting rays will be detected.
                # Have to use patch subdivision and convex hull polyhedron to do this.
                missedintersections += 1
            else
                if isa(allintersections, Interval)
                    allintersections = (allintersections,)
                end

                for intsct in allintersections
                    if lower(intsct) isa RayOrigin
                        missedintersections += 1
                    else
                        # closest point should be on the surface, furthest should be on bounding prism
                        @test isapprox(point(lower(intsct)), pt, rtol = RTOLERANCE, atol = ATOLERANCE)
                        @test isa(upper(intsct), Intersection{Float64,3}) || (OpticSim.direction(r)[3] == -1 && isa(upper(intsct), Infinity{Float64}))
                    end
                end
            end
        end

        z2 = TestData.gridsagsurfacebicubiccheby()
        az2 = AcceleratedParametricSurface(z2, 20)

        samples2 = samplepoints(numsamples, -0.95, 0.95, -0.95, 0.95)
        for uv in samples2
            pt = point(z2, uv[1], uv[2])
            origin = [0.0, 0.0, 5.0]
            r = Ray(origin, pt .- origin)
            allintersections = surfaceintersection(az2, r)
            if (allintersections isa EmptyInterval)
                # sampled triangle acceleration structure doesn't guarantee all intersecting rays will be detected.
                # Have to use patch subdivision and convex hull polyhedron to do this.
                missedintersections += 1
            else
                if isa(allintersections, Interval)
                    allintersections = (allintersections,)
                end

                for intsct in allintersections
                    if lower(intsct) isa RayOrigin
                        missedintersections += 1
                    else
                        # closest point should be on the surface, furthest should be on bounding prism
                        @test isapprox(point(lower(intsct)), pt, rtol = RTOLERANCE, atol = ATOLERANCE)
                        @test isa(upper(intsct), Intersection{Float64,3}) || (OpticSim.direction(r)[3] == -1 && isa(upper(intsct), Infinity{Float64}))
                    end
                end
            end
        end

        if missedintersections > 0
            @warn "$missedintersections out of total $(length(samples1) + length(samples2)) gridsag suface intersections were missed"
        end

        # hit from inside
        r = Ray([0.0, 0.0, -0.5], [0.1, 0.2, 0.5])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.13038780645120746, 0.26077561290241486, 0.15193903225603717], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit from outside
        r = Ray([0.0, 0.0, 0.5], [0.1, 0.2, -0.5])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [0.046677740288643785, 0.0933554805772876, 0.26661129855678095], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [0.6708203932499369, 1.3416407864998738, -2.8541019662496843], rtol = RTOLERANCE, atol = ATOLERANCE)

        # miss from inside
        r = Ray([0.2, 0.2, -0.5], [0.0, 0.0, -1.0])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && (upper(intvl) isa Infinity)

        # miss from outside
        r = Ray([0.2, 2.0, -0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 2.0, 0.5], [0.0, 0.0, -1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 0.2, 0.5], [0.0, 0.0, 1.0])
        @test surfaceintersection(az1, r) isa EmptyInterval
        r = Ray([0.2, 2.0, 2.0], [0.0, -1.0, 0.0])
        @test surfaceintersection(az1, r) isa EmptyInterval

        # two hits from outside
        r = Ray([2.0, 0.0, 0.25], [-1.0, 0.0, 0.0])
        hit = surfaceintersection(az1, r)
        a = isapprox(point(lower(hit[1])), [1.5, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[1])), [1.394132887629247, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)
        b = isapprox(point(lower(hit[2])), [0.30184444700168467, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(hit[2])), [-0.6437764440680096, 0.0, 0.25], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isa(hit, DisjointUnion) && length(hit) == 2 && a && b

        # hit cyl from outside
        r = Ray([0.0, 2.0, -0.5], [0.0, -1.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test isapprox(point(lower(intvl)), [0.0, 1.5, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(point(upper(intvl)), [0.0, -1.5, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE)

        # hit cyl from inside
        r = Ray([0.0, 0.0, -0.5], [0.0, -1.0, 0.0])
        intvl = surfaceintersection(az1, r)
        @test (lower(intvl) isa RayOrigin) && isapprox(point(upper(intvl)), [0.0, -1.5, -0.5], rtol = RTOLERANCE, atol = ATOLERANCE)
    end # testset GridSag
end # testset intersection


text:
"""
_check_add!!_behavour(acc, val)

This checks that `acc + val` is the same as `add!!(acc, val)`.
It matters primarily for types that overload `add!!` such as `InplaceableThunk`s.

`acc` is the value that has been accumulated so far.
`val` is a deriviative, being accumulated into `acc`.

`kwargs` are all passed on to isapprox
"""
function _check_add!!_behavour(acc, val; kwargs...)
    # Note, we don't test that `acc` is actually mutated because it doesn't have to be
    # e.g. if it is immutable. We do test the `add!!` return value.
    # That is what people should rely on. The mutation is just to save allocations.
    acc_mutated = deepcopy(acc)  # prevent this test changing others
    @test isapprox(add!!(acc_mutated, val), acc + val; kwargs...)
end

ext: jl


text:
println("Hello World!")
#> Hello World!


text:
using Test, LinearAlgebra, Random, Statistics

dir = "I:\\My Documents\\00 AFIT\\Research\\Julia Projects\\StrategyGames\\src"
include(joinpath(dir, "cfrkuhn.jl"))

@testset "Kuhn CFR" begin
    Random.seed!(23408)
    @test actions("JQp") == ["p", "b"]
    @test infoset("KQpb") == 13
    @test player("") == 3
    @test player("JQp") == 2

    h = "KJpp"
    i = player(h)
    πi = 1.0
    πo = 1.0
    leafutility(h, i, πi, πo)

    u, r, s, σ = cfr(200_000)
    α = σ[2][2] # α = p(b | J)
    optstrategy_p = [1 - α, 1.0, 1 - 3α, 2/3, 1.0, 0.0, 1.0, 2/3, 0.0, 1.0, 2/3 - α, 0.0] # optimal p(pass) for each info set
    cfrstrategy_p = [σ[i][1] for i = 2:13] # calculated p(pass) for each infoset
    @test isapprox(cfrstrategy_p, optstrategy_p, atol = .05)
    @test mean(u) ≈ -0.055 atol = 0.01
end
#
# using Test, LinearAlgebra, Random, Statistics
#
# dir = "I:\\My Documents\\00 AFIT\\Research\\Julia Projects\\StrategyGames\\src"
# include(joinpath(dir, "cfrkuhn.jl"))
#
Random.seed!(23408)
T = 100_000
u, r, s, σ = cfr(T)
um = [sum(u[1:i]) / i for i in 1:T]
ustar = -1/18
rerr = abs(mean(u[(end - Int(T * 0.1)):end]) - ustar) / abs(ustar)

using Plots; gr()
plot(1:T, um)
plot!(1:T, fill(ustar, T))

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
"""
This type conforms to the mesh interface but is specialised for the case of a segment of the real axis subdived in equally sized intervals. Typical use is a discretisation of the time axis.
"""
mutable struct SegmentedAxis{T}
    timestep::T
    numsteps::Int
end

numcells(ax::SegmentedAxis) = ax.numsteps
#cellvertices(ax::SegmentedAxis, i) = SVector(point((i-1)*ax.timestep), point(i*ax.timestep))

#Base.size(a::SegmentedAxis) = (a.numsteps,)
#Base.linearindexing(a::SegmentedAxis) = Base.LinearFast()
#Base.getindex(a::SegmentedAxis, i) = ((i-1)*a.timestep, i*a.timestep)

function minmaxdist(τ, σ)
	T = eltype(τ[1])
	m = norm(τ[1]-σ[1])
	M = zero(T)
	for i in 1:length(τ)
		p = τ[i]
		for j in 1:length(σ)
			q = σ[j]
			d = norm(p-q)
			d < m && (m=d)
			d > M && (M=d)
		end
	end
	return m, M
end

function rings(τ, σ, ΔR)
	m, M = minmaxdist(τ, σ)
	r0 = floor(Int, m/ΔR) + 1
	r1 = ceil(Int, M/ΔR+1)
	r0 : r1
end

ring(r, ΔR) = ((r-1)*ΔR, r*ΔR)


text:
import Base.ht_keyindex

# Cf. https://github.com/JuliaStats/StatsBase.jl/issues/135
struct UnitWeights <: AbstractVector{Int} end
Base.getindex(w::UnitWeights, ::Integer...) = 1
Base.getindex(w::UnitWeights, ::AbstractVector) = w

# About the type inference limitation which prompts this workaround, see
# https://github.com/JuliaLang/julia/issues/10880
Base.@pure eltypes(T) = Tuple{map(eltype, T.parameters)...}
Base.@pure vectypes(T) = Tuple{map(U -> Vector{U}, T.parameters)...}

# Internal function needed for now so that n is inferred
function _freqtable(x::Tuple,
                    skipmissing::Bool = false,
                    weights::AbstractVector{<:Real} = UnitWeights(),
                    subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing)
    n = length(x)
    n == 0 && throw(ArgumentError("at least one argument must be provided"))

    if !isa(subset, Nothing)
        x = map(y -> y[subset], x)
        weights = weights[subset]
    end

    l = map(length, x)
    vtypes = eltypes(typeof(x))

    for i in 1:n
        if l[1] != l[i]
            error("arguments are not of the same length: $l")
        end
    end

    if !isa(weights, UnitWeights) && length(weights) != l[1]
        error("'weights' (length $(length(weights))) must be of the same length as vectors (length $(l[1]))")
    end

    d = Dict{vtypes, eltype(weights)}()

    for (i, el) in enumerate(zip(x...))
        index = ht_keyindex(d, el)

        if index > 0
            @inbounds d.vals[index] += weights[i]
        else
            @inbounds d[el] = weights[i]
        end
    end

    if skipmissing
        filter!(p -> !any(ismissing, p[1]), d)
    end

    keyvec = collect(keys(d))

    dimnames = Vector{Vector}(undef, n)
    for i in 1:n
        s = Set{vtypes.parameters[i]}()
        for j in 1:length(keyvec)
            push!(s, keyvec[j][i])
        end

        # convert() is needed for Union{T, Missing}, which currently gives a Vector{Any}
        # which breaks inference of the return type
        dimnames[i] = convert(Vector{vtypes.parameters[i]}, unique(s))
        try
            sort!(dimnames[i])
        catch err
            err isa MethodError || rethrow(err)
        end
    end

    a = zeros(eltype(weights), map(length, dimnames)...)::Array{eltype(weights), n}
    na = NamedArray(a, tuple(dimnames...)::vectypes(vtypes), ntuple(i -> "Dim$i", n))

    for (k, v) in d
        na[Name.(k)...] = v
    end

    na
end

"""
    freqtable(x::AbstractVector...;
              skipmissing::Bool = false,
              weights::AbstractVector{<:Real} = UnitWeights(),
              subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing])

    freqtable(t, cols::Symbol...;
              skipmissing::Bool = false,
              weights::AbstractVector{<:Real} = UnitWeights(),
              subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing])

Create frequency table from vectors or table columns.

`t` can be any type of table supported by the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface.

**Examples**

```jldoctest
julia> freqtable([1, 2, 2, 3, 4, 3])
4-element Named Array{Int64,1}
Dim1  │
──────┼──
1     │ 1
2     │ 2
3     │ 2
4     │ 1

julia> df = DataFrame(x=[1, 2, 2, 2], y=[1, 2, 1, 2]);

julia> freqtable(df, :x, :y)
2×2 Named Array{Int64,2}
x ╲ y │ 1  2
──────┼─────
1     │ 1  0
2     │ 1  2

julia> freqtable(df, :x, :y, subset=df.x .> 1)
1×2 Named Array{Int64,2}
x ╲ y │ 1  2
──────┼─────
2     │ 1  2

```
"""

freqtable(x::AbstractVector...;
          skipmissing::Bool = false,
          weights::AbstractVector{<:Real} = UnitWeights(),
          subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing) =
    _freqtable(x, skipmissing, weights, subset)

# Internal function needed for now so that n is inferred
function _freqtable(x::NTuple{n, AbstractCategoricalVector}, skipmissing::Bool = false,
                    weights::AbstractVector{<:Real} = UnitWeights(),
                    subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing) where n
    n == 0 && throw(ArgumentError("at least one argument must be provided"))

    if !isa(subset, Nothing)
        x = map(y -> y[subset], x)
        weights = weights[subset]
    end

    len = map(length, x)
    miss = map(v -> eltype(v) >: Missing, x)
    lev = map(x) do v
        eltype(v) >: Missing && !skipmissing ? [levels(v); missing] : allowmissing(levels(v))
    end
    dims = map(length, lev)
    # First entry is for missing values (only correct and used if present)
    ord = map((v, d) -> Int[d; CategoricalArrays.order(v.pool)], x, dims)

	for i in 1:n
	    if len[1] != len[i]
	        error(string("arguments are not of the same length: ", tuple(len...)))
	    end
	end

    if !isa(weights, UnitWeights) && length(weights) != len[1]
        error("'weights' (length $(length(weights))) must be of the same length as vectors (length $(len[1]))")
    end

    sizes = cumprod([dims...])
    a = zeros(eltype(weights), dims)
    missingpossible = any(miss)

    @inbounds for i in 1:len[1]
        ref = x[1].refs[i]
        el = ord[1][ref + 1]
        anymiss = missingpossible & (ref <= 0)

        for j in 2:n
            ref = x[j].refs[i]
            anymiss |= missingpossible & (ref <= 0)
            el += (ord[j][ref + 1] - 1) * sizes[j - 1]
        end

        if !(missingpossible && skipmissing && anymiss)
            a[el] += weights[i]
        end
    end

    NamedArray(a, lev, ntuple(i -> "Dim$i", n))
end

freqtable(x::AbstractCategoricalVector...; skipmissing::Bool = false,
          weights::AbstractVector{<:Real} = UnitWeights(),
          subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing) =
    _freqtable(x, skipmissing, weights, subset)

function freqtable(t, cols::Symbol...; args...)
    all_cols = Tables.columns(t)
    a = freqtable((getproperty(all_cols, y) for y in cols)...; args...)
    setdimnames!(a, cols)
    a
end

"""
    prop(tbl::AbstractArray{<:Number};
         margins = nothing)

Create a table of proportions from an array `tbl`.

If `margins` is `nothing` (the default), proportions over the whole `tbl` are computed.
If `margins` is an `Integer`, or an iterable of `Integer`s, proportions
sum to `1` over dimensions specified by `margins`.
In particular for a two-dimensional array, when `margins` is `1` row proportions are
calculated, and when `margins` is `2` column proportions are calculated.

This function does not check that `tbl` contains only non-negative values.

Calculating `sum(proptable(..., margins=margins), dims=dims)` with `dims` equal to the complement
of `margins` produces an array containing only `1.0` (see last example below).

**Examples**

```jldoctest
julia> prop([1 2; 3 4])
2×2 Array{Float64,2}:
 0.1  0.2
 0.3  0.4

julia> prop([1 2; 3 4], margins=1)
2×2 Array{Float64,2}:
 0.333333  0.666667
 0.428571  0.571429

julia> prop([1 2; 3 4], margins=2)
2×2 Array{Float64,2}:
 0.25  0.333333
 0.75  0.666667

julia> prop([1 2; 3 4], margins=(1, 2))
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> pt = prop(reshape(1:12, (2, 2, 3)), margins=3)
2×2×3 Array{Float64,3}:
[:, :, 1] =
 0.1  0.3
 0.2  0.4

[:, :, 2] =
 0.192308  0.269231
 0.230769  0.307692

[:, :, 3] =
 0.214286  0.261905
 0.238095  0.285714

julia> sum(pt, dims=(1, 2))
1×1×3 Array{Float64,3}:
[:, :, 1] =
 1.0

[:, :, 2] =
 1.0

[:, :, 3] =
 1.0

```
"""

function prop(tbl::AbstractArray{<:Number,N}; margins=nothing) where N
    if margins === nothing
        return tbl / sum(tbl)
    else
        lo, hi = extrema(margins)
        (lo < 1 || hi > N) && throw(ArgumentError("margins must be a valid dimension"))
        return tbl ./ sum(tbl, dims=tuple(setdiff(1:N, margins)...)::NTuple{N-length(margins),Int})
    end
end

prop(tbl::NamedArray{<:Number}; margins=nothing) =
    NamedArray(prop(convert(Array, tbl); margins=margins), tbl.dicts, tbl.dimnames)

"""
    proptable(x::AbstractVector...;
              margins = nothing,
              skipmissing::Bool = false,
              weights::AbstractVector{<:Real} = UnitWeights(),
              subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing])

    proptable(t, cols::Symbol...;
              margins = nothing,
              skipmissing::Bool = false,
              weights::AbstractVector{<:Real} = UnitWeights(),
              subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing])

Create a frequency table of proportions from vectors or table columns.
This is equivalent to calling `prop(freqtable(...), margins=margins)`.

`t` can be any type of table supported by the [Tables.jl](https://github.com/JuliaData/Tables.jl)
interface.

If `margins` is `nothing` (the default), proportions over the whole table are computed.
If `margins` is an `Integer`, or an iterable of `Integer`s, proportions
sum to `1` over dimensions specified by `margins`.
In particular for a two-dimensional array, when `margins` is `1` row proportions are
calculated, and when `margins` is `2` column proportions are calculated.

Calculating `sum(proptable(..., margins=margins), dims=dims)` with `dims` equal to the complement
of `margins` produces an array containing only `1.0` (see last example below).

**Examples**

```jldoctest
julia> proptable([1, 2, 2, 3, 4, 3])
4-element Named Array{Float64,1}
Dim1  │
──────┼─────────
1     │ 0.166667
2     │ 0.333333
3     │ 0.333333
4     │ 0.166667

julia> df = DataFrame(x=[1, 2, 2, 2, 1, 1], y=[1, 2, 1, 2, 2, 2], z=[1, 1, 1, 2, 2, 1]);

julia> proptable(df, :x, :y)
2×2 Named Array{Float64,2}
x ╲ y │        1         2
──────┼───────────────────
1     │ 0.166667  0.333333
2     │ 0.166667  0.333333

julia> proptable(df, :x, :y, subset=df.x .> 1)
1×2 Named Array{Float64,2}
x ╲ y │        1         2
──────┼───────────────────
2     │ 0.333333  0.666667

julia> proptable([1, 2, 2, 2], [1, 1, 1, 2], margins=1)
2×2 Named Array{Float64,2}
Dim1 ╲ Dim2 │        1         2
────────────┼───────────────────
1           │      1.0       0.0
2           │ 0.666667  0.333333

julia> proptable([1, 2, 2, 2], [1, 1, 1, 2], margins=2)
2×2 Named Array{Float64,2}
Dim1 ╲ Dim2 │        1         2
────────────┼───────────────────
1           │ 0.333333       0.0
2           │ 0.666667       1.0

julia> proptable([1, 2, 2, 2], [1, 1, 1, 2], margins=(1,2))
2×2 Named Array{Float64,2}
Dim1 ╲ Dim2 │   1    2
────────────┼─────────
1           │ 1.0  NaN
2           │ 1.0  1.0

julia> proptable(df.x, df.y, df.z)
2×2×2 Named Array{Float64,3}

[:, :, Dim3=1] =
Dim1 ╲ Dim2 │        1         2
────────────┼───────────────────
1           │ 0.166667  0.166667
2           │ 0.166667  0.166667

[:, :, Dim3=2] =
Dim1 ╲ Dim2 │        1         2
────────────┼───────────────────
1           │      0.0  0.166667
2           │      0.0  0.166667

julia> pt = proptable(df.x, df.y, df.z, margins=(1,2))
2×2×2 Named Array{Float64,3}

[:, :, Dim3=1] =
Dim1 ╲ Dim2 │   1    2
────────────┼─────────
1           │ 1.0  0.5
2           │ 1.0  0.5

[:, :, Dim3=2] =
Dim1 ╲ Dim2 │   1    2
────────────┼─────────
1           │ 0.0  0.5
2           │ 0.0  0.5

julia> sum(pt, dims=3)
2×2×1 Named Array{Float64,3}

[:, :, Dim3=sum(Dim3)] =
Dim1 ╲ Dim2 │   1    2
────────────┼─────────
1           │ 1.0  1.0
2           │ 1.0  1.0

```
"""
proptable(x::AbstractVector...;
          margins = nothing,
          skipmissing::Bool = false,
          weights::AbstractVector{<:Real} = UnitWeights(),
          subset::Union{Nothing, AbstractVector{Int}, AbstractVector{Bool}} = nothing) =
    prop(freqtable(x...,
                   skipmissing=skipmissing, weights=weights, subset=subset), margins=margins)

proptable(t, cols::Symbol...; margins=nothing, kwargs...) =
    prop(freqtable(t, cols...; kwargs...), margins=margins)


text:
using AlgebraicDynamics.DWDDynam
using Catlab.WiringDiagrams, Catlab.Present, Catlab.Theories, Catlab.Programs
using LabelledArrays
using DelayDiffEq
using Test

# Define all of the wiring diagrams that will be used
@present C(FreeBiproductCategory) begin
    A::Ob
    f::Hom(A, A)
    g::Hom(A, A⊗A)
    h::Hom(A⊗A, A)
end

d_id = @program C (x::A -> f(x))
d12 = @program(C, (x::A) -> f(f(x)))
d_copymerge = @program(C, (x::A) -> [f(x), f(x)])

d_trace = copy(d_id); b, = box_ids(d_trace)
add_wire!(d_trace, (b, 1) => (b, 1))

d_tot = ocompose(d_copymerge, [d12, d_id])

d_big = WiringDiagram([:A, :A], [:A, :A, :A])
b1 = add_box!(d_big, Box(:f, [:A, :A], [:A]))
b2 = add_box!(d_big, Box(:g, [:A], [:A, :A]))
b3 = add_box!(d_big, Box(:h, [:A], [:A]))

bin = input_id(d_big); bout = output_id(d_big)

add_wires!(d_big, Pair[
    (bin, 1) => (b1, 2),
    (bin, 1) => (b2, 1),
    (bin, 2) => (b3, 1),
    (b1, 1) => (bout, 1),
    (b2, 1) => (bout, 1),
    (b2, 1) => (bout, 2),
    (b2, 2) => (b2, 1),
    (b3, 1) => (b2, 1),
    (b3, 1) => (bout, 3)
])

# Identity
uf(u, x, p, t) = [x[1] - u[1]]
rf(u, args...) = u
mf = ContinuousMachine{Float64}(1,1,1, uf, rf)

m_id = oapply(d_id, [mf])

x0 = 1
p0 = 200
@test eval_dynamics(m_id, [x0], [p0]) == [p0 - x0]
@test readout(m_id, [x0]) == [x0]

# unfed parameter
m12 = oapply(d12, [mf, mf])

x0 = -1
y0 = 29
p0 = 200
@test eval_dynamics(m12, [x0, y0], [p0]) == [p0 - x0, x0 - y0]
@test readout(m12, [x0,y0]) == [y0]

# break and back together
m = oapply(d_copymerge, Dict(:f => mf))

x0 = -1
y0 = 29
p0 = 200
@test eval_dynamics(m, [x0, y0], [p0]) == [p0 - x0, p0 - y0]
@test readout(m, [x0,y0]) == [x0 + y0]

@test ninputs(m) == 1
@test nstates(m) == 2
@test noutputs(m) == 1

# trace
m_trace = oapply(d_trace, mf)
@test nstates(m_trace) == 1
@test eval_dynamics(m_trace, [x0], [p0]) == [p0]
@test readout(m_trace, [x0]) == [x0]

# oapply and ocompose
m_tot1 = oapply(d_tot, mf)
m_tot2 = oapply(d_copymerge, [m12, m_id])
x0 = [-1, 5.5, 20]
@test nstates(m_tot1) == 3
@test nstates(m_tot2) == 3

@test eval_dynamics(m_tot1, x0, [p0]) == eval_dynamics(m_tot2, x0, [p0])
@test eval_dynamics(m_tot1, x0, [p0]) == [p0 - x0[1], x0[1] - x0[2], p0 - x0[3]]

@test readout(m_tot1, x0) == [x0[2] + x0[3]]
@test readout(m_tot2, x0) == [x0[2] + x0[3]]

# big diagram
m1 = ContinuousMachine{Float64}(2,1,1,
        (u, x, p, t) -> [x[1] * x[2]  - u[1]],
        (u,p,t) -> 2*u)
m2 = ContinuousMachine{Float64}(1, 2, 2,
        (u, x, p, t) -> [u[1]*u[2], x[1]^2*u[2]],
        (u,p,t) -> u)

m3 = ContinuousMachine{Float64}(1,2,1,
        (u, x, p, t) -> [u[1]^2 - x[1], u[2] - u[1]],
        (u,p,t) -> [u[1] + u[2]])

xs = Dict(:f => m1, :g => m2, :h => m3, :j => mf)
m = oapply(d_big, xs)
@test ninputs(m) == 2
@test nstates(m) == 5
@test noutputs(m) == 3
u = [2.0, 3.0, 5.0, -7.0, -0.5]
x = [0.1, 11.0]
@test eval_dynamics(m, u, x) ≈
    [-u[1], u[2]*u[3], (x[1]+u[3]+u[4]+u[5])^2*u[3], u[4]^2 - x[2], u[5] - u[4]]
@test readout(m, u) == [2*u[1] + u[2], u[2], u[4] + u[5]]
@test readout(oapply(d_id, m3), [u[4], u[5]]) == [u[4] + u[5]]

# eulers
h = 0.15
euler_m = oapply(d_big, euler_approx([m1, m2, m3], h))
@test eval_dynamics(euler_m, u, x) == u + h*eval_dynamics(m, u, x)
@test eval_dynamics(euler_m, u, x) == eval_dynamics(euler_approx(m, h), u, x)
euler_m2 = oapply(d_big, euler_approx([m1, m2, m3]))
@test eval_dynamics(euler_m, u, x) == eval_dynamics(euler_m2, u, x, [h])

h = 0.02
xs = Dict(:h => m3, :f => m1, :g => m2)
euler_m = oapply(d_big, euler_approx(xs, h))
@test eval_dynamics(euler_m, u, x) == u + h*eval_dynamics(m, u, x)
@test eval_dynamics(euler_m, u, x) == eval_dynamics(euler_approx(m, h), u, x)
euler_m2 = oapply(d_big, euler_approx(xs))
@test eval_dynamics(euler_m, u, x) == eval_dynamics(euler_m2, u, x, [h])

# delay differential equations
df(u,x,h,p,t) = h(p, t - p.τ)
rf(u,h,p,t) = u
delay = DelayMachine{Float64}(1,1,1,df,rf)

delay_copy = oapply(d_id, delay)

hist(p,t) = [0.0]

u0 = [2.7]
x0 = [10.0]
τ = 10.0
@test eval_dynamics(delay_copy, u0, x0, hist, LVector(τ = τ)) == [0.0]

prob = DDEProblem(delay, u0, x0, hist, (0.0, τ - 0.01), LVector(τ = τ))
alg = MethodOfSteps(Tsit5())
@test last(solve(prob, alg)) == u0

prob = DDEProblem(delay, u0, x0, hist, (0.0, 2*τ), LVector(τ = τ))
prob_copy = DDEProblem(delay_copy, u0, x0, hist, (0.0, 2*τ), LVector(τ = τ))

@test last(solve(prob, alg)) == last(solve(prob_copy, alg))

# delay differential equation with analytic solution via method of steps
# dx/dt = -x(t-1) with history x(t) = 10 for t <= 0
df(u,x,h,p,t) = -h(p, t - 1.0) # dynamics

delay_machine = DelayMachine{Float64}(0,1,0,df,(u,h,p,t) -> u)
u0 = 10.0
p = nothing

hist(p, t) = [u0] # history

prob = DDEProblem(delay_machine, [u0], [], hist, (0.0, 3.0), p)
alg = MethodOfSteps(Tsit5())

sol = solve(prob, alg,abstol=1e-12,reltol=1e-12)

# solve over 0<t<1; we have x(t-1)=10 over this interval therefore
# dx/dt = -10; dx = -10dt; \int{dx} = \int{-10dt}; x + c1 = -10t + c2
# x(t) = c - 10t; x(t) = 10 - 10t
@test sol(0.5)[1] ≈ 10.0 - 10.0*0.5
@test sol(0.9)[1] ≈ 10.0 - 10.0*0.9

# solve over 1<t<2; we have x(t-1)=10-10(t-1) over this interval therefore
# dx/dt = -(10 - 10(t-1)); \int{dx} = \int{-(10 - 10(t-1))dt};
# x + c1 = -10(t-1) + 10\int{(t-1) dt}; x + c1 = -10t + 10((t-1)^2/2) + c2;
# x(t) = -10(t-1) + 5(t-1)^2
@test sol(1.5)[1] ≈ -10.0*(1.5-1.0) + 5.0*(1.5-1.0)^2
@test sol(1.9)[1] ≈ -10.0*(1.9-1.0) + 5.0*(1.9-1.0)^2

# Test composite without readout delay
d = ocompose(d_trace, [d12])
df(u,x,h,p,t) = x[1]*h(p, t - p.τ)
ef(u,x,h,p,t) = x + h(p, t - p.τ)
delay_rf(u,h,p,t) = h(p, t - p.τ)

mult_delay = DelayMachine{Float64}(1,1,1, df, rf)
add_delay = DelayMachine{Float64}(1,1,1, ef, rf)

hist(p,t) = [0.0, 0.0]
u0 = [1.0, 1.0]
x0 = [1.0]
prob = DDEProblem(oapply(d, [mult_delay, add_delay]), u0, x0, hist, (0, 4.0), LVector(τ = 2.0))
sol1 = solve(prob, alg; dtmax = 0.1)
f = (u,h,p,t) -> [ (u[2] + x0[1]) * h(p,t - p.τ)[1], u[1] + h(p,t - p.τ)[2] ]
prob = DDEProblem(f, u0, hist, (0, 4.0), LVector(τ = 2.0))
sol2 = solve(prob, alg; dtmax = 0.1)

@test last(sol1) == last(sol2)

# test composite with readout delay
delay_rf(u,h,p,t) = h(p, t - p.τ)

mult_delay = DelayMachine{Float64}(1,1,1, df, delay_rf)
add_delay = DelayMachine{Float64}(1,1,1, ef, delay_rf)

hist(p,t) = [0.0, 0.0]
u0 = [1.0, 1.0]
x0 = [1.0]
prob = DDEProblem(oapply(d, [mult_delay, add_delay]), u0, x0, hist, (0, 4.0), LVector(τ = 2.0))
sol1 = solve(prob, alg; dtmax = 0.1)
f = (u,h,p,t) -> [ (h(p,t - p.τ)[2] + x0[1]) * h(p,t - p.τ)[1], h(p, t - p.τ)[1] + h(p,t - p.τ)[2] ]
prob = DDEProblem(f, u0, hist, (0, 4.0), LVector(τ = 2.0))
sol2 = solve(prob, alg; dtmax = 0.1)

@test last(sol1) == last(sol2)

# VectorInterface
m1 = ContinuousMachine{Float64,3}(2, 3, 1,
        (u, x, p, t) -> x[1] + x[2],
        (u,p,t) -> [u])
m2 = ContinuousMachine{Float64, 3}(1, 3, 2,
        (u, x, p, t) -> x[1] + u,
        (u,p,t) -> [u, 2*u])

m3 = ContinuousMachine{Float64, 3}(1, 3, 1,
        (u, x, p, t) -> x[1],
        (u,p,t) -> [u])

m = oapply(d_big, [m1, m2, m3])

u1 = 1:3; u2 = 4:6; u3 = 7:9;
x1 = ones(3); x2 = fill(0.5, 3)

@test readout(m, vcat(u1, u2, u3), nothing, 0) == [u1 + u2, u2, u3]
@test eval_dynamics(m, vcat(u1, u2, u3), [x1, x2], nothing, 0) == vcat(x1, u2 + x1 + 2*u2 + u3, x2)


text:
# Tests for distances
function t_ppca()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

for tr in [A01, A02, A03, vectordata]
	for ts in [A01, A02, A03, vectordata]
		Test.@test try
			Wppca = tr |> j4pr.ppca(maxoutdim=1)
			Wppcar = j4pr.ppcar(Wppca)
			ts |> Wppca
			(ts |> Wppca) |> Wppcar
			true
		catch
			false
		end
	end
end

# Short test for ppcar on vector data
Test.@test try
	Wppca = A01 |> j4pr.ppca(maxoutdim=1)
	Wppcar = j4pr.ppcar(Wppca)
	(A01 |> Wppca) |> Wppcar
	true
catch
	false
end

for tr in [A04, A05, A06, matrixdata]
	for ts in [A04, A05, A06, matrixdata]
		Test.@test try
			Wppca = tr |> j4pr.ppca()
			Wppcar = j4pr.ppcar(Wppca)
			ts |> Wppca
			(ts |> Wppca) |> Wppcar
			true
		catch
			false
		end
	end
end

end


text:
using DrWatson
@quickactivate "ExercisesRepo"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

# %% showing embedding
lo = Systems.lorenz([0, 10, 0.0])

tr = trajectory(lo, 1000; Ttr=100)
x, y, z = columns(tr)
w = @. x
τ = estimate_delay(w, "mi_min")
ψ = @. z - y

# figure for Cao's method
fig = figure(figsize=(figx/2.5, figx/2.5));
for (s, l) in zip((w, ψ), ("\$w=x\$", "\$w=z-y\$"))
    τ = estimate_delay(s, "mi_min")
    Ds = estimate_dimension(s, τ, 1:6, "afnn")
    plot(2:7, Ds, marker = "o", label = l)
end
xticks(2:7)
xlabel("\$d\$")
ylabel("\$E_{d}\$")
legend(title="measurement")
fig.tight_layout(pad = 0.1)
# fsave(joinpath(figdir, "caodemonstration"), fig)

# boxes = 10 .^ (-3:0.2:1)
# D_A = generalized_dim(0.0, tr, boxes)
# D_R = generalized_dim(0.0, R, boxes)
# @show (D_A, D_R)

# %%
γ = 2
R = reconstruct(w, γ, τ)

close("all")
N = 5000
fig = figure()
ax1 = subplot(131, projection = "3d")
ax1.plot3D(x[1:N], y[1:N], z[1:N], lw = 1.0, color = "C1")
ax1.set_title("original set \$A\$")
ax = ax1
for a in (ax.xaxis, ax.yaxis, ax.zaxis); a.set_ticklabels([]); end
# ax1.text(-20, -20, 0, "\$D_0(A)\$ = $(round(D_A;digits=3))", color = "C1")

ax2 = subplot(132)
wx = 400:1000
ax2.plot(wx, w[wx])
tau = 700:(700+τ)
ax2.plot(tau, fill(10, length(tau)), color = "C3")
ax2.text(700, 7.5, "\$\\tau\$", color = "C3", size = 32)
ax2.set_title("measurement \$w=x\$")
# make fancy axis
ax2.spines["bottom"].set_position("center")
ax2.spines["right"].set_color("none")
ax2.spines["top"].set_color("none")
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_xlabel("\$t\$")
ax2.set_ylabel("\$w\$", rotation = 0)
ax2.xaxis.set_label_coords(1.0, 0.5)
ax2.yaxis.set_label_coords(-0.05, 0.9)

ax3 = subplot(133, projection = "3d")
ax3.plot(R[1:N, 1], R[1:N, 2], R[1:N, 3], color = "C2", lw = 1.0)
ax = ax3
for a in (ax.xaxis, ax.yaxis, ax.zaxis); a.set_ticklabels([]); end
ax3.set_title("reconstruction \$R\$")
# ax3.set_title("reconstruction \$R\$\n\$(\\gamma=$(γ),\\tau=$(τ))\$")
# ax3.text(-5,-25,-10, "\$D_0(R)\$ = $(round(D_R;digits=3))", color = "C2")

fig.tight_layout(pad = 0.1)
# fsave(joinpath(figdir, "delayembedding"), fig)

# %% value of τ demonstration
γ = 3
τ1 = 2
τ2 = τ
τ3 = 2τ2
θ=-30
e=10

fig = figure()
for (i, τ) in enumerate((0, τ1, τ2, τ3))
    R = i == 1 ? tr : reconstruct(w, γ, τ)
    ax = subplot(1, 4, i, projection = "3d")
    ax.plot(R[1:N, 1], R[1:N, 2], R[1:N, 3],
    color = i == 1 ? "C1" : "C2", lw = 1,
    zorder = 1)
    ax.view_init(elev=e, azim=θ)
    ax.set_title(i == 1 ? "original" : "\$R:\\,\\tau=$(τ)\$", pad = -10)
    for a in (ax.xaxis, ax.yaxis, ax.zaxis); a.set_ticklabels([]); end
    dj = 65
    for j in (456, )
        χ, ψ, ζ = R[j:dj:j+dj, 1], R[j:dj:j+dj, 2], R[j:dj:j+dj, 3]
        ax.scatter3D(χ, ψ, ζ,depthshade = false,c = "C0", zorder = 99, s = 50)
    end
end
# ax = subplot(1, 4, 1, projection = "3d")
# ax.plot3D(x[1:N], y[1:N], z[1:N], lw = 1.0, color = "C1")
# ax.view_init(elev=e, azim=θ)
# for a in (ax.xaxis, ax.yaxis, ax.zaxis); a.set_ticklabels([]); end

fig.tight_layout()
fig.subplots_adjust(wspace = 0.0001, bottom = -0.1, top = 1, left = 0.0, right = 1)
# fsave(joinpath(figdir, "taudemonstration"), fig)

# %% broomhead king
using DynamicalSystems, PyPlot
ds = Systems.lorenz()
tr = trajectory(ds, 100.0; Ttr=10)
w = tr[:, 1]
w .+= rand(length(w)) #add noise
U, S = broomhead_king(w, 42)
R = reconstruct(w, 2, 17)

close("all")
fig = figure(figsize = (1.2figx/2, figx/2))
subplot(121)
plot(U[:, 1], U[:, 2], lw = 1.0)
title("Broomhead-King")
xticks([]); yticks([])
subplot(122)
plot(R[:, 1], R[:, 2], lw = 1.0, color = "C2")
xticks([]); yticks([])
title("delay embedding")
tight_layout()
subplots_adjust(wspace = 0.05, left = 0.025, right = 0.975)
fsave(joinpath(figdir, "broomhead"), fig)

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
using Dance.Router

function hello(headers::Dict{String, String}) :: String
    return "Hello World"
end

route("/", hello; method=GET, endpoint=EP_HTML)


text:
mutable struct Particle{T <: Real}
    x::AbstractArray{T}
    v::AbstractArray{T}
    x_best::AbstractArray{T}
    min_dim::T
    max_dim::T

    function Particle{T}(
        x::AbstractArray{T},
        v::AbstractArray{T},
        x_best::AbstractArray{T},
        a::T,
        b::T
    ) where {T <: AbstractFloat}

        @assert length(x) == length(v) == length(x_best) "Dimension must be unique"

        return new(x, v, x_best, a, b)
    end
end

"""
    Particle(
        x::AbstractArray{T},
        v::AbstractArray{T},
        x_best::AbstractArray{T},
        a::T,
        b::T
    ) where {T <: Real}

A type that can hold information about current position, current velocity,
the _best_ candidate to a solution, as well as defining the bounds.
The dimensions of the `Particle` are inferred from the length of the arrays.

# Arguments
- `x`: Array that holds the **positions** of possible solutions.
- `v`: Array that holds **velocities** related to `x`.
- `x_best`: An element of `x` that determines the best position for the particle.
- `a`: lower bound for `x`
- `b`: upper bound for `v`

# Example
```julia
p = Particle(zeros(3), rand(3), zeros(3), 0.0, 1.0)
```
"""
Particle(
    x::AbstractArray{T},
    v::AbstractArray{T},
    x_best::AbstractArray{T},
    a::T,
    b::T
) where {T <: Real} =
    Particle{T}(x, v, x_best, a, b)

"""
    Particle(a, b, n::Int)

`Particle` that can be created randomly using the bounds and the dimension needed.

# Arguments
- `a`: _lower_ bound for `x` and `v`.
- `b`: _upper_ bound for `x` and `v`.
- `n`: dimension for `x`, `v`, and `x_best`.

# Example
```julia
p = Particle(-1.0, 1.0, 3)
```
"""
function Particle(a, b, n::Int, rng)
    @assert n > 0 "Dimension is always positive"

    a, b = promote(a, b)
    x = a .+ (rand(rng, n) * (b - a))
    v = a .+ (rand(rng, n) * (b - a))
    x_best = rand(rng, n)

    return Particle(x, v, x_best, a, b)
end

mutable struct Population end

"""
    Population(num_particles::T, dim::T, a::V, b::V)
        where {T<:Int, V<:AbstractFloat} -> Vector{Particle}(undef, num_particles)

An array of `Particle`s where each of them are bounded and are given a dimension. This is
essentially a multi-dimensional array. It makes handling `Particle`s much easier.

# Arguments
- `num_particles`: Number of particles in the `Population`.
- `dim`: Dimension for every `Particle`.
- `a`: Lower bound for every `Particle`, this is shared across every instance.
- `b`: Upper bound for every `Particle`, this is shared across every instance.

# Example
```julia
pop = Population(35, 4, -1.0, 1.0)
```
"""
function Population(num_particles::T, dim::T, a, b, rng) where {T <: Int}
    @assert dim > 0 "Dimension is always positive"
    @assert num_particles > 0 "There must be at least 1 Particle in the Population"

    # if isnothing(rng)
    #     rng = Xorshifts.Xoroshiro128Plus()
    # end

    container = Vector{Particle}(undef, num_particles)
    for idx in 1:num_particles
        container[idx] = Particle(a, b, dim, rng)
    end

    return container
end

"""
    Population(num_particles::Int, ranges; seed=nothing)
        -> Vector{Particle}(undef, num_particles)

An array of `Particle`'s where each of them are bounded and are given a dimension.
`x` is a collection of ranges for each *dimension* for the `Particle`'s specified.

# Arguments
- `num_particles`: Number of particles in the `Population`.
- `dim`: Dimension for every `Particle`.
- `x`: Should be a collection of ranges, for instance a tuple or a vector.

# Example
```julia
# Two ranges, one for each dimension
range_a = SVector(-10.0, 10.0)
range_b = SVector(-2.5, 2.0)
pops = Population(2, 20, [ranges_a, range_b])
```
"""
function Population(num_particles::Int, ranges, rng)
    dim = length(ranges)
    @assert dim > 0 "Dimension is always positive"
    @assert num_particles > 0 "There must be at least 1 Particle in the Population"

    container = Vector{Particle}(undef, num_particles)
    element_type_range = eltype(ranges[1]) # Use the type from the first range's elements
    random_values = Matrix{element_type_range}(undef, dim, num_particles)
    rand!(rng, random_values)

    @assert length(ranges) == length(eachrow(random_values))

    @inbounds for (r, idx) in zip(ranges, eachrow(random_values))
        a, b = r # Upper and lower bounds
        # Apply the boundaries for each dimension
        @. idx *= b - a
        @. idx += a
    end

    # Now that each row has random values for each boundary range, we loop over
    # each particle and each row and create the corresponding array
    @inbounds for i in 1:num_particles
        a, b = extrema(random_values[:, i])
        container[i] = Particle(random_values[:, i],
            random_values[:, i], randn(rng, dim), a, b)
    end

    return container
end


text:
using CovidNYT
using Test


text:
using Colors
using Plots
using NoisySignalIntegration
using Random: seed!

gausspeaks(x, p) = sum([@. A * 1/√(2π * σ^2) * exp(-(x - μ)^2 / (2σ^2)) for (A, μ, σ) in p])

c1, c2, c3 = let
    seed!(1)
    x = 0:0.1:100
    y = gausspeaks(x, [(3.3, 30.0, 4.0), (4.1, 70.0, 4.0), (1, 50.0, 20.0)])
    c = Curve(x, y)
    uc = add_noise(c, MvGaussianNoiseModel(0.1, 0.01, 1.0), 3)
    c1, c2, c3 = [NoisySignalIntegration.get_draw(i, uc) for i in 1:3]
    c2 += 0.5
    c3 += 1
    c1, c2, c3
end

p = let
    r = [203, 60, 51]  ./ 255
    v = [149, 88, 178] ./ 255
    b = [64, 99, 216]  ./ 255
    g = [56, 152, 38]  ./ 255
    plot(; background_color = :transparent, color=:black, ticks=:none, legend=false, linewidth=3, axis=false);
    plot!(c1, 18., 40., fillcolor=RGB(r...), fillalpha=1);
    plot!(c1, 60., 84., fillcolor=RGB(v...), fillalpha=1);
    plot!(c2, 19., 42., fillcolor=RGB(b...), fillalpha=1);
    plot!(c2, 59., 84., fillcolor=RGB(g...), fillalpha=1);
    plot!(c1; color=:black, linewidth=8);
    plot!(c2; color=:black, linewidth=8)
    plot!(annotations=(70, 0.13, Plots.text("?", 50, RGBA(1, 1, 1, 0.25))))
end

savefig(p, joinpath(@__DIR__, "src", "assets", "logo.png"))


text:
using EconPDEs

# PDE with 1 state variable
include("./AssetPricing/CampbellCochrane.jl")
@time pdesolve(m, stategrid, yend)
#   0.011112 seconds (286.45 k allocations: 12.268 MiB)

# PDE with 2 state variables
include("./AssetPricing/BansalYaron.jl")
@time pdesolve(m, stategrid, yend)
# 0.081095 seconds (485.84 k allocations: 78.910 MiB, 11.41% gc time)

# System of 4 PDEs with 1 state variable
include("./AssetPricing/GarleanuPanageas.jl")
@time pdesolve(m, stategrid, yend)
# 0.047346 seconds (167.95 k allocations: 65.049 MiB, 12.90% gc time)

# System of 3 PDEs with 2 state variables
include("./AssetPricing/DiTella.jl")
@time pdesolve(m, stategrid, yend)
# 0.912719 seconds (2.55 M allocations: 1.303 GiB, 9.01% gc time)


text:
module TestFactorizationMachines

using FactorizationMachines
using Base.Test

matrix_dim = 4
X = eye(matrix_dim)
y = [.1, .2, .3, .4]

info("Solving system of equations with identity matrix")
fm = train(sparse(X), y, model_params = Models.gauss(k0 = false, k1 = true, num_factors = 0))
@test fm.model.w0 == 0
@test size(fm.model.V) == (0, matrix_dim)
@test_approx_eq_eps 0.0 evaluate(Evaluators.rmse(), fm, sparse(X), y) 1e-1

info("Solving system of equations with identity matrix")
fm = train(sparse(X), y, model_params = Models.gauss(k0 = false, k1 = false))
@test fm.model.w0 == 0
@test fm.model.w == zeros(matrix_dim)
@test_approx_eq_eps 0.0 evaluate(Evaluators.rmse(), fm, sparse(X), y) 1e-1

end