text:
using Intervals: Endpoint, Left, Right, LeftEndpoint, RightEndpoint

@testset "Endpoint" begin
    @testset "constructors" begin
        for D in (Left, Right)
            @test Endpoint{Int, D, Closed}(0).endpoint == 0
            @test Endpoint{Int, D, Open}(0).endpoint == 0
            @test Endpoint{Int, D, Unbounded}(nothing).endpoint isa Int

            @test_throws MethodError Endpoint{Int, D, Closed}(nothing)
            @test_throws MethodError Endpoint{Int, D, Open}(nothing)
            @test_throws MethodError Endpoint{Int, D, Unbounded}(0)

            @test_throws MethodError Endpoint{Nothing, D, Closed}(nothing)
            @test_throws MethodError Endpoint{Nothing, D, Open}(nothing)
            @test Endpoint{Nothing, D, Unbounded}(nothing).endpoint === nothing

            @test Endpoint{Int, D, Closed}(0.0).endpoint == 0
            @test Endpoint{Int, D, Open}(0.0).endpoint == 0
            @test_throws MethodError Endpoint{Int, D, Unbounded}(0.0)
        end
    end

    @testset "bounded" begin
        @testset "LeftEndpoint < LeftEndpoint" begin
            @test LeftEndpoint{Open}(1) < LeftEndpoint{Open}(2.0)
            @test LeftEndpoint{Closed}(1) < LeftEndpoint{Open}(2.0)
            @test LeftEndpoint{Open}(1) < LeftEndpoint{Closed}(2.0)
            @test LeftEndpoint{Closed}(1) < LeftEndpoint{Closed}(2.0)

            @test !(LeftEndpoint{Open}(1) < LeftEndpoint{Open}(1.0))
            @test LeftEndpoint{Closed}(1) < LeftEndpoint{Open}(1.0)
            @test !(LeftEndpoint{Open}(1) < LeftEndpoint{Closed}(1.0))
            @test !(LeftEndpoint{Closed}(1) < LeftEndpoint{Closed}(1.0))

            @test !(LeftEndpoint{Open}(2) < LeftEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Closed}(2) < LeftEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Open}(2) < LeftEndpoint{Closed}(1.0))
            @test !(LeftEndpoint{Closed}(2) < LeftEndpoint{Closed}(1.0))
        end

        @testset "LeftEndpoint <= LeftEndpoint" begin
            @test LeftEndpoint{Open}(1) <= LeftEndpoint{Open}(2.0)
            @test LeftEndpoint{Closed}(1) <= LeftEndpoint{Open}(2.0)
            @test LeftEndpoint{Open}(1) <= LeftEndpoint{Closed}(2.0)
            @test LeftEndpoint{Closed}(1) <= LeftEndpoint{Closed}(2.0)

            @test LeftEndpoint{Open}(1) <= LeftEndpoint{Open}(1.0)
            @test LeftEndpoint{Closed}(1) <= LeftEndpoint{Open}(1.0)
            @test !(LeftEndpoint{Open}(1) <= LeftEndpoint{Closed}(1.0))
            @test LeftEndpoint{Closed}(1) <= LeftEndpoint{Closed}(1.0)

            @test !(LeftEndpoint{Open}(2) <= LeftEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Closed}(2) <= LeftEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Open}(2) <= LeftEndpoint{Closed}(1.0))
            @test !(LeftEndpoint{Closed}(2) <= LeftEndpoint{Closed}(1.0))
        end

        @testset "RightEndpoint < RightEndpoint" begin
            @test RightEndpoint{Open}(1) < RightEndpoint{Open}(2.0)
            @test RightEndpoint{Closed}(1) < RightEndpoint{Open}(2.0)
            @test RightEndpoint{Open}(1) < RightEndpoint{Closed}(2.0)
            @test RightEndpoint{Closed}(1) < RightEndpoint{Closed}(2.0)

            @test !(RightEndpoint{Open}(1) < RightEndpoint{Open}(1.0))
            @test !(RightEndpoint{Closed}(1) < RightEndpoint{Open}(1.0))
            @test RightEndpoint{Open}(1) < RightEndpoint{Closed}(1.0)
            @test !(RightEndpoint{Closed}(1) < RightEndpoint{Closed}(1.0))

            @test !(RightEndpoint{Open}(2) < RightEndpoint{Open}(1.0))
            @test !(RightEndpoint{Closed}(2) < RightEndpoint{Open}(1.0))
            @test !(RightEndpoint{Open}(2) < RightEndpoint{Closed}(1.0))
            @test !(RightEndpoint{Closed}(2) < RightEndpoint{Closed}(1.0))
        end

        @testset "RightEndpoint <= RightEndpoint" begin
            @test RightEndpoint{Open}(1) <= RightEndpoint{Open}(2.0)
            @test RightEndpoint{Closed}(1) <= RightEndpoint{Open}(2.0)
            @test RightEndpoint{Open}(1) <= RightEndpoint{Closed}(2.0)
            @test RightEndpoint{Closed}(1) <= RightEndpoint{Closed}(2.0)

            @test RightEndpoint{Open}(1) <= RightEndpoint{Open}(1.0)
            @test !(RightEndpoint{Closed}(1) <= RightEndpoint{Open}(1.0))
            @test RightEndpoint{Open}(1) <= RightEndpoint{Closed}(1.0)
            @test RightEndpoint{Closed}(1) <= RightEndpoint{Closed}(1.0)

            @test !(RightEndpoint{Open}(2) <= RightEndpoint{Open}(1.0))
            @test !(RightEndpoint{Closed}(2) <= RightEndpoint{Open}(1.0))
            @test !(RightEndpoint{Open}(2) <= RightEndpoint{Closed}(1.0))
            @test !(RightEndpoint{Closed}(2) <= RightEndpoint{Closed}(1.0))
        end

        @testset "LeftEndpoint < RightEndpoint" begin
            @test LeftEndpoint{Open}(1) < RightEndpoint{Open}(2.0)
            @test LeftEndpoint{Closed}(1) < RightEndpoint{Open}(2.0)
            @test LeftEndpoint{Open}(1) < RightEndpoint{Closed}(2.0)
            @test LeftEndpoint{Closed}(1) < RightEndpoint{Closed}(2.0)

            @test !(LeftEndpoint{Open}(1) < RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Closed}(1) < RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Open}(1) < RightEndpoint{Closed}(1.0))
            @test !(LeftEndpoint{Closed}(1) < RightEndpoint{Closed}(1.0))

            @test !(LeftEndpoint{Open}(2) < RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Closed}(2) < RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Open}(2) < RightEndpoint{Closed}(1.0))
            @test !(LeftEndpoint{Closed}(2) < RightEndpoint{Closed}(1.0))
        end

        @testset "LeftEndpoint <= RightEndpoint" begin
            @test LeftEndpoint{Open}(1) <= RightEndpoint{Open}(2.0)
            @test LeftEndpoint{Closed}(1) <= RightEndpoint{Open}(2.0)
            @test LeftEndpoint{Open}(1) <= RightEndpoint{Closed}(2.0)
            @test LeftEndpoint{Closed}(1) <= RightEndpoint{Closed}(2.0)

            @test !(LeftEndpoint{Open}(1) <= RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Closed}(1) <= RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Open}(1) <= RightEndpoint{Closed}(1.0))
            @test LeftEndpoint{Closed}(1) <= RightEndpoint{Closed}(1.0)

            @test !(LeftEndpoint{Open}(2) <= RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Closed}(2) <= RightEndpoint{Open}(1.0))
            @test !(LeftEndpoint{Open}(2) <= RightEndpoint{Closed}(1.0))
            @test !(LeftEndpoint{Closed}(2) <= RightEndpoint{Closed}(1.0))
        end

        @testset "RightEndpoint < LeftEndpoint" begin
            @test RightEndpoint{Open}(1) < LeftEndpoint{Open}(2.0)
            @test RightEndpoint{Closed}(1) < LeftEndpoint{Open}(2.0)
            @test RightEndpoint{Open}(1) < LeftEndpoint{Closed}(2.0)
            @test RightEndpoint{Closed}(1) < LeftEndpoint{Closed}(2.0)

            @test RightEndpoint{Open}(1) < LeftEndpoint{Open}(1.0)
            @test RightEndpoint{Closed}(1) < LeftEndpoint{Open}(1.0)
            @test RightEndpoint{Open}(1) < LeftEndpoint{Closed}(1.0)
            @test !(RightEndpoint{Closed}(1) < LeftEndpoint{Closed}(1.0))

            @test !(RightEndpoint{Open}(2) < LeftEndpoint{Open}(1.0))
            @test !(RightEndpoint{Closed}(2) < LeftEndpoint{Open}(1.0))
            @test !(RightEndpoint{Open}(2) < LeftEndpoint{Closed}(1.0))
            @test !(RightEndpoint{Closed}(2) < LeftEndpoint{Closed}(1.0))
        end

        @testset "RightEndpoint <= LeftEndpoint" begin
            @test RightEndpoint{Open}(1) <= LeftEndpoint{Open}(2.0)
            @test RightEndpoint{Closed}(1) <= LeftEndpoint{Open}(2.0)
            @test RightEndpoint{Open}(1) <= LeftEndpoint{Closed}(2.0)
            @test RightEndpoint{Closed}(1) <= LeftEndpoint{Closed}(2.0)

            @test RightEndpoint{Open}(1) <= LeftEndpoint{Open}(1.0)
            @test RightEndpoint{Closed}(1) <= LeftEndpoint{Open}(1.0)
            @test RightEndpoint{Open}(1) <= LeftEndpoint{Closed}(1.0)
            @test RightEndpoint{Closed}(1) <= LeftEndpoint{Closed}(1.0)

            @test !(RightEndpoint{Open}(2) <= LeftEndpoint{Open}(1.0))
            @test !(RightEndpoint{Closed}(2) <= LeftEndpoint{Open}(1.0))
            @test !(RightEndpoint{Open}(2) <= LeftEndpoint{Closed}(1.0))
            @test !(RightEndpoint{Closed}(2) <= LeftEndpoint{Closed}(1.0))
        end

        @testset "$T < Scalar" for T in (LeftEndpoint, RightEndpoint)
            @test T{Open}(1) < 2.0
            @test T{Closed}(1) < 2.0

            @test (T{Open}(1) < 1.0) == (T === RightEndpoint)
            @test !(T{Closed}(1) < 1.0)

            @test !(T{Open}(1) < 0.0)
            @test !(T{Closed}(1) < 0.0)
        end

        @testset "$T <= Scalar" for T in (LeftEndpoint, RightEndpoint)
            @test T{Open}(1) <= 2.0
            @test T{Closed}(1) <= 2.0

            @test (T{Open}(1) <= 1.0) == (T === RightEndpoint)
            @test T{Closed}(1) <= 1.0

            @test !(T{Open}(1) <= 0.0)
            @test !(T{Closed}(1) <= 0.0)
        end

        @testset "Scalar < $T" for T in (LeftEndpoint, RightEndpoint)
            @test 0 < T{Open}(1.0)
            @test 0 < T{Closed}(1.0)

            @test (1 < T{Open}(1.0)) == (T === LeftEndpoint)
            @test !(1 < T{Closed}(1.0))

            @test !(2 < T{Open}(1.0))
            @test !(2 < T{Closed}(1.0))
        end

        @testset "Scalar <= $T" for T in (LeftEndpoint, RightEndpoint)
            @test 0 <= T{Open}(1.0)
            @test 0 <= T{Closed}(1.0)

            @test (1 <= T{Open}(1.0)) == (T === LeftEndpoint)
            @test 1 <= T{Closed}(1.0)

            @test !(2 <= T{Open}(1.0))
            @test !(2 <= T{Closed}(1.0))
        end

        @testset "LeftEndpoint == LeftEndpoint" begin
            @test LeftEndpoint{Open}(1) != LeftEndpoint{Open}(2.0)
            @test LeftEndpoint{Closed}(1) != LeftEndpoint{Open}(2.0)
            @test LeftEndpoint{Open}(1) != LeftEndpoint{Closed}(2.0)
            @test LeftEndpoint{Closed}(1) != LeftEndpoint{Closed}(2.0)

            @test LeftEndpoint{Open}(1) == LeftEndpoint{Open}(1.0)
            @test LeftEndpoint{Closed}(1) != LeftEndpoint{Open}(1.0)
            @test LeftEndpoint{Open}(1) != LeftEndpoint{Closed}(1.0)
            @test LeftEndpoint{Closed}(1) == LeftEndpoint{Closed}(1.0)
        end

        @testset "RightEndpoint == RightEndpoint" begin
            @test RightEndpoint{Open}(1) != RightEndpoint{Open}(2.0)
            @test RightEndpoint{Closed}(1) != RightEndpoint{Open}(2.0)
            @test RightEndpoint{Open}(1) != RightEndpoint{Closed}(2.0)
            @test RightEndpoint{Closed}(1) != RightEndpoint{Closed}(2.0)

            @test RightEndpoint{Open}(1) == RightEndpoint{Open}(1.0)
            @test RightEndpoint{Closed}(1) != RightEndpoint{Open}(1.0)
            @test RightEndpoint{Open}(1) != RightEndpoint{Closed}(1.0)
            @test RightEndpoint{Closed}(1) == RightEndpoint{Closed}(1.0)
        end

        @testset "LeftEndpoint == RightEndpoint" begin
            @test LeftEndpoint{Open}(1) != RightEndpoint{Open}(2.0)
            @test LeftEndpoint{Closed}(1) != RightEndpoint{Open}(2.0)
            @test LeftEndpoint{Open}(1) != RightEndpoint{Closed}(2.0)
            @test LeftEndpoint{Closed}(1) != RightEndpoint{Closed}(2.0)

            @test LeftEndpoint{Open}(1) != RightEndpoint{Open}(1.0)
            @test LeftEndpoint{Closed}(1) != RightEndpoint{Open}(1.0)
            @test LeftEndpoint{Open}(1) != RightEndpoint{Closed}(1.0)
            @test LeftEndpoint{Closed}(1) == RightEndpoint{Closed}(1.0)
        end

        @testset "RightEndpoint == LeftEndpoint" begin
            @test RightEndpoint{Open}(1) != LeftEndpoint{Open}(2.0)
            @test RightEndpoint{Closed}(1) != LeftEndpoint{Open}(2.0)
            @test RightEndpoint{Open}(1) != LeftEndpoint{Closed}(2.0)
            @test RightEndpoint{Closed}(1) != LeftEndpoint{Closed}(2.0)

            @test RightEndpoint{Open}(1) != LeftEndpoint{Open}(1.0)
            @test RightEndpoint{Closed}(1) != LeftEndpoint{Open}(1.0)
            @test RightEndpoint{Open}(1) != LeftEndpoint{Closed}(1.0)
            @test RightEndpoint{Closed}(1) == LeftEndpoint{Closed}(1.0)
        end

        @testset "$T == Scalar" for T in (LeftEndpoint, RightEndpoint)
            @test T{Open}(0) != 1.0
            @test T{Closed}(0) != 1.0

            @test T{Open}(1) != 1.0
            @test T{Closed}(1) == 1.0

            @test T{Open}(2) != 1.0
            @test T{Closed}(2) != 1.0
        end

        @testset "Scalar == $T" for T in (LeftEndpoint, RightEndpoint)
            @test 1.0 != T{Open}(0)
            @test 1.0 != T{Closed}(0)

            @test 1.0 != T{Open}(1)
            @test 1.0 == T{Closed}(1)

            @test 1.0 != T{Open}(2)
            @test 1.0 != T{Closed}(2)
        end

        @testset "isequal" begin
            @test isequal(LeftEndpoint{Closed}(0.0), LeftEndpoint{Closed}(0.0))
            @test isequal(LeftEndpoint{Open}(0.0), LeftEndpoint{Open}(0.0))
            @test !isequal(LeftEndpoint{Closed}(-0.0), LeftEndpoint{Open}(0.0))
            @test !isequal(LeftEndpoint{Open}(-0.0), LeftEndpoint{Closed}(0.0))
            @test !isequal(LeftEndpoint{Closed}(-0.0), LeftEndpoint{Closed}(0.0))
            @test !isequal(LeftEndpoint{Open}(-0.0), LeftEndpoint{Open}(0.0))

            @test isequal(RightEndpoint{Closed}(0.0), LeftEndpoint{Closed}(0.0))
            @test !isequal(LeftEndpoint{Closed}(-0.0), RightEndpoint{Open}(0.0))
            @test !isequal(RightEndpoint{Open}(-0.0), LeftEndpoint{Closed}(0.0))
            @test !isequal(LeftEndpoint{Closed}(-0.0), RightEndpoint{Closed}(0.0))
        end

        @testset "hash" begin
            # Need a complicated enough element type for this test to possibly fail. Using a
            # ZonedDateTime with a VariableTimeZone should do the trick.
            a = now(tz"Europe/London")
            b = deepcopy(a)
            @test hash(a) == hash(b)  # Double check

            @test hash(LeftEndpoint{Open}(a)) == hash(LeftEndpoint{Open}(b))
            @test hash(LeftEndpoint{Closed}(a)) != hash(LeftEndpoint{Open}(b))
            @test hash(LeftEndpoint{Open}(a)) != hash(LeftEndpoint{Closed}(b))
            @test hash(LeftEndpoint{Closed}(a)) == hash(LeftEndpoint{Closed}(b))

            @test hash(RightEndpoint{Open}(a)) == hash(RightEndpoint{Open}(b))
            @test hash(RightEndpoint{Closed}(a)) != hash(RightEndpoint{Open}(b))
            @test hash(RightEndpoint{Open}(a)) != hash(RightEndpoint{Closed}(b))
            @test hash(RightEndpoint{Closed}(a)) == hash(RightEndpoint{Closed}(b))

            @test hash(LeftEndpoint{Open}(a)) != hash(RightEndpoint{Open}(b))
            @test hash(LeftEndpoint{Closed}(a)) != hash(RightEndpoint{Open}(b))
            @test hash(LeftEndpoint{Open}(a)) != hash(RightEndpoint{Closed}(b))
            @test hash(LeftEndpoint{Closed}(a)) != hash(RightEndpoint{Closed}(b))
        end
    end

    # Note: The value for unbounded endpoints is irrelevant
    @testset "unbounded" begin
        @testset "LeftEndpoint < LeftEndpoint" begin
            @test !(LeftEndpoint{Unbounded}(nothing) < LeftEndpoint{Unbounded}(nothing))
            @test LeftEndpoint{Unbounded}(nothing) < LeftEndpoint{Open}(0.0)
            @test LeftEndpoint{Unbounded}(nothing) < LeftEndpoint{Closed}(0.0)
            @test !(LeftEndpoint{Open}(0) < LeftEndpoint{Unbounded}(nothing))
            @test !(LeftEndpoint{Closed}(0) < LeftEndpoint{Unbounded}(nothing))
        end

        @testset "LeftEndpoint <= LeftEndpoint" begin
            @test LeftEndpoint{Unbounded}(nothing) <= LeftEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Unbounded}(nothing) <= LeftEndpoint{Open}(0.0)
            @test LeftEndpoint{Unbounded}(nothing) <= LeftEndpoint{Closed}(0.0)
            @test !(LeftEndpoint{Open}(0) <= LeftEndpoint{Unbounded}(nothing))
            @test !(LeftEndpoint{Closed}(0) <= LeftEndpoint{Unbounded}(nothing))
        end

        @testset "RightEndpoint < RightEndpoint" begin
            @test !(RightEndpoint{Unbounded}(nothing) < RightEndpoint{Unbounded}(nothing))
            @test !(RightEndpoint{Unbounded}(nothing) < RightEndpoint{Open}(0.0))
            @test !(RightEndpoint{Unbounded}(nothing) < RightEndpoint{Closed}(0.0))
            @test RightEndpoint{Open}(0) < RightEndpoint{Unbounded}(nothing)
            @test RightEndpoint{Closed}(0) < RightEndpoint{Unbounded}(nothing)
        end

        @testset "RightEndpoint <= RightEndpoint" begin
            @test RightEndpoint{Unbounded}(nothing) <= RightEndpoint{Unbounded}(nothing)
            @test !(RightEndpoint{Unbounded}(nothing) <= RightEndpoint{Open}(0.0))
            @test !(RightEndpoint{Unbounded}(nothing) <= RightEndpoint{Closed}(0.0))
            @test RightEndpoint{Open}(0) <= RightEndpoint{Unbounded}(nothing)
            @test RightEndpoint{Closed}(0) <= RightEndpoint{Unbounded}(nothing)
        end

        @testset "LeftEndpoint < RightEndpoint" begin
            @test LeftEndpoint{Unbounded}(nothing) < RightEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Unbounded}(nothing) < RightEndpoint{Open}(0.0)
            @test LeftEndpoint{Unbounded}(nothing) < RightEndpoint{Closed}(0.0)
            @test LeftEndpoint{Open}(0) < RightEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Closed}(0) < RightEndpoint{Unbounded}(nothing)
        end

        @testset "LeftEndpoint <= RightEndpoint" begin
            @test LeftEndpoint{Unbounded}(nothing) <= RightEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Unbounded}(nothing) <= RightEndpoint{Open}(0.0)
            @test LeftEndpoint{Unbounded}(nothing) <= RightEndpoint{Closed}(0.0)
            @test LeftEndpoint{Open}(0) <= RightEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Closed}(0) <= RightEndpoint{Unbounded}(nothing)
        end

        @testset "RightEndpoint < LeftEndpoint" begin
            @test !(RightEndpoint{Unbounded}(nothing) < LeftEndpoint{Unbounded}(nothing))
            @test !(RightEndpoint{Unbounded}(nothing) < LeftEndpoint{Open}(0.0))
            @test !(RightEndpoint{Unbounded}(nothing) < LeftEndpoint{Closed}(0.0))
            @test !(RightEndpoint{Open}(0) < LeftEndpoint{Unbounded}(nothing))
            @test !(RightEndpoint{Closed}(0) < LeftEndpoint{Unbounded}(nothing))
        end

        @testset "RightEndpoint <= LeftEndpoint" begin
            @test !(RightEndpoint{Unbounded}(nothing) < LeftEndpoint{Unbounded}(nothing))
            @test !(RightEndpoint{Unbounded}(nothing) < LeftEndpoint{Open}(0.0))
            @test !(RightEndpoint{Unbounded}(nothing) < LeftEndpoint{Closed}(0.0))
            @test !(RightEndpoint{Open}(0) < LeftEndpoint{Unbounded}(nothing))
            @test !(RightEndpoint{Closed}(0) < LeftEndpoint{Unbounded}(nothing))
        end

        @testset "LeftEndpoint < Scalar" begin
            @test LeftEndpoint{Unbounded}(nothing) < -Inf
            @test LeftEndpoint{Unbounded}(nothing) < Inf
        end

         @testset "LeftEndpoint <= Scalar" begin
            @test LeftEndpoint{Unbounded}(nothing) <= -Inf
            @test LeftEndpoint{Unbounded}(nothing) <= Inf
        end

        @testset "RightEndpoint < Scalar" begin
            @test !(RightEndpoint{Unbounded}(nothing) < -Inf)
            @test !(RightEndpoint{Unbounded}(nothing) < Inf)
        end

        @testset "RightEndpoint <= Scalar" begin
            @test !(RightEndpoint{Unbounded}(nothing) <= -Inf)
            @test !(RightEndpoint{Unbounded}(nothing) <= Inf)
        end

        @testset "Scalar < LeftEndpoint" begin
            @test !(-Inf < LeftEndpoint{Unbounded}(nothing))
            @test !(Inf < LeftEndpoint{Unbounded}(nothing))
        end

        @testset "Scalar <= LeftEndpoint" begin
            @test !(-Inf <= LeftEndpoint{Unbounded}(nothing))
            @test !(Inf <= LeftEndpoint{Unbounded}(nothing))
        end

        @testset "Scalar < RightEndpoint" begin
            @test -Inf < RightEndpoint{Unbounded}(nothing)
            @test Inf < RightEndpoint{Unbounded}(nothing)
        end

        @testset "Scalar < RightEndpoint" begin
            @test -Inf <= RightEndpoint{Unbounded}(nothing)
            @test Inf <= RightEndpoint{Unbounded}(nothing)
        end

        @testset "LeftEndpoint == LeftEndpoint" begin
            @test LeftEndpoint{Unbounded}(nothing) == LeftEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Unbounded}(nothing) != LeftEndpoint{Open}(0.0)
            @test LeftEndpoint{Unbounded}(nothing) != LeftEndpoint{Closed}(0.0)
            @test LeftEndpoint{Open}(0) != LeftEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Closed}(0) != LeftEndpoint{Unbounded}(nothing)
        end

        @testset "RightEndpoint == RightEndpoint" begin
            @test RightEndpoint{Unbounded}(nothing) == RightEndpoint{Unbounded}(nothing)
            @test RightEndpoint{Unbounded}(nothing) != RightEndpoint{Open}(0.0)
            @test RightEndpoint{Unbounded}(nothing) != RightEndpoint{Closed}(0.0)
            @test RightEndpoint{Open}(0) != RightEndpoint{Unbounded}(nothing)
            @test RightEndpoint{Closed}(0) != RightEndpoint{Unbounded}(nothing)
        end

        @testset "LeftEndpoint == RightEndpoint" begin
            @test LeftEndpoint{Unbounded}(nothing) != RightEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Unbounded}(nothing) != RightEndpoint{Open}(0.0)
            @test LeftEndpoint{Unbounded}(nothing) != RightEndpoint{Closed}(0.0)
            @test LeftEndpoint{Open}(0) != RightEndpoint{Unbounded}(nothing)
            @test LeftEndpoint{Closed}(0) != RightEndpoint{Unbounded}(nothing)
        end

        @testset "RightEndpoint == LeftEndpoint" begin
            @test RightEndpoint{Unbounded}(nothing) != LeftEndpoint{Unbounded}(nothing)
            @test RightEndpoint{Unbounded}(nothing) != LeftEndpoint{Open}(0.0)
            @test RightEndpoint{Unbounded}(nothing) != LeftEndpoint{Closed}(0.0)
            @test RightEndpoint{Open}(0) != LeftEndpoint{Unbounded}(nothing)
            @test RightEndpoint{Closed}(0) != LeftEndpoint{Unbounded}(nothing)
        end

        @testset "isequal" begin
            T = Float64

            @test isequal(LeftEndpoint{T,Unbounded}(nothing), LeftEndpoint{T,Unbounded}(nothing))
            @test !isequal(LeftEndpoint{T,Unbounded}(nothing), LeftEndpoint{T,Open}(0.0))
            @test !isequal(LeftEndpoint{T,Unbounded}(nothing), LeftEndpoint{T,Closed}(0.0))
            @test !isequal(LeftEndpoint{T,Open}(-0.0), LeftEndpoint{T,Unbounded}(nothing))
            @test !isequal(LeftEndpoint{T,Closed}(-0.0), LeftEndpoint{T,Unbounded}(nothing))

            @test !isequal(RightEndpoint{Unbounded}(nothing), LeftEndpoint{Unbounded}(nothing))
            @test !isequal(LeftEndpoint{Unbounded}(nothing), RightEndpoint{Unbounded}(nothing))
        end

         @testset "hash" begin
            # Note: Unbounded endpoints should ignore the value
            T = Int

            left_unbounded = LeftEndpoint{T,Unbounded}(nothing)
            @test hash(left_unbounded) == hash(LeftEndpoint{T,Unbounded}(nothing))
            @test hash(left_unbounded) != hash(LeftEndpoint{Unbounded}(nothing))
            @test hash(left_unbounded) != hash(LeftEndpoint{Open}(left_unbounded.endpoint))
            @test hash(left_unbounded) != hash(LeftEndpoint{Closed}(left_unbounded.endpoint))

            right_unbounded = RightEndpoint{T,Unbounded}(nothing)
            @test hash(right_unbounded) == hash(RightEndpoint{T,Unbounded}(nothing))
            @test hash(right_unbounded) != hash(RightEndpoint{Unbounded}(nothing))
            @test hash(right_unbounded) != hash(RightEndpoint{Open}(right_unbounded.endpoint))
            @test hash(right_unbounded) != hash(RightEndpoint{Closed}(right_unbounded.endpoint))

            @test hash(LeftEndpoint{T,Unbounded}(nothing)) != hash(RightEndpoint{T,Unbounded}(nothing))
            @test hash(LeftEndpoint{T,Unbounded}(nothing)) != hash(RightEndpoint{T,Open}(0))
            @test hash(LeftEndpoint{T,Unbounded}(nothing)) != hash(RightEndpoint{T,Closed}(0))
            @test hash(LeftEndpoint{T,Open}(0)) != hash(RightEndpoint{T,Unbounded}(nothing))
            @test hash(LeftEndpoint{T,Closed}(0)) != hash(RightEndpoint{T,Unbounded}(nothing))
        end
    end

    @testset "broadcast" begin
        test = [
            LeftEndpoint{Open}(0),
            LeftEndpoint{Closed}(0),
            LeftEndpoint{Unbounded}(nothing),
            RightEndpoint{Open}(0),
            RightEndpoint{Closed}(0),
            RightEndpoint{Unbounded}(nothing),
        ]

        # Verify that Endpoint is treated as a scalar during broadcast
        result = test .== 0
        @test result == [false, true, false, false, true, false]
    end
end


text:
# Autogenerated wrapper script for IVerilog_jll for armv7l-linux-gnueabihf-cxx11
export iverilog, iverilog_vpi

using Readline_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("IVerilog")
JLLWrappers.@declare_executable_product(iverilog)
JLLWrappers.@declare_executable_product(iverilog_vpi)
function __init__()
    JLLWrappers.@generate_init_header(Readline_jll, Zlib_jll)
    JLLWrappers.@init_executable_product(
        iverilog,
        "bin/iverilog",
    )

    JLLWrappers.@init_executable_product(
        iverilog_vpi,
        "bin/iverilog-vpi",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
isfiniterange(::,k) = false
isfiniterange(B::KroneckerOperator,k::Integer) = isfinite(size(B.ops[k],1))
isfiniterange(B::PlusOperator,k::Integer) = isfiniterange(first(B.ops),k)

function findfunctionals(A::Vector,k::Integer)
    T=eltype(eltype(eltype(A)))
    indsBx=find(f->isfiniterange(f,k),A)
    if k==1
        indsBx,Operator{T}[(@assert dekron(Ai,2)==ConstantOperator(Float64,1.0); dekron(Ai,1)) for Ai in A[indsBx]]
    else
        @assert k==2
        indsBx,Operator{T}[(@assert dekron(Ai,1)==ConstantOperator(Float64,1.0); dekron(Ai,2)) for Ai in A[indsBx]]
    end
end

# Operators that are diagonal in a dimension can be solved in O(n) operations
# isdiagop is used to inspect an operator to see if it is diagonal

Base.isdiag{T<:Number}(B::Operator{T}) = bandinds(B)==(0,0)

# multivariate case
isdiagop(K::KroneckerOperator,k) = isdiag(K.ops[k])
isdiagop(S::WrapperOperator,k) = isdiagop(S.op,k)
isdiagop(A::Union{PlusOperator,TimesOperator},k) = all(op->isdiagop(op,k),A.ops)

# diagop gets out the op corresponding to the k-th column
# TODO: diag in x direction

diagop(A::KroneckerOperator,col)=A.ops[2]

diagop(A::PlusOperator,col) = mapreduce(op->diagop(op,col),+,A.ops)
diagop(A::TimesOperator,col) = mapreduce(op->diagop(op,col),*,A.ops)
diagop(A::SpaceOperator,col) = diagop(A.op,col)
# TODO: Space for ConstantOperator
diagop(A::ConstantOperator,col) = ConstantOperator(A.c)
diagop(A::ConstantTimesOperator,col) = A.c*diagop(A.op,col)

diagop(L,k)=error("Override diagop for "*string(typeof(L)))


text:
f = readlines(homedir()*"/data/Facebook100/Facebook_Sets.txt")
graphs = split(f[1])

##  Try something else

C = [0.8 0 0;
0.8 0.8 0;
0 0 .8;
0 .8 .8;
.5 0 .1;
.2 .3 .1;
0 .5 0;
0 .2 .2;
0 0 1;
0.1 0.5 0.6;
0 .75 .75;
.4 .5 .3]

RGBA(.5, .6, 0,1)

RGBA(.1,.5,0,1)
nextcol = 1
Psize = plot(grid = false, legend = false)
Pmax = plot(grid = false, legend = false)
Prun = plot(grid = false)
Pdensity = plot(grid = false,legend = false)
Pmean = plot(grid = false, legend = false)
counter = 0
pathtoFB100 = homedir()*"/data/Facebook100/"

xx = Vector{Float64}()
yy = Vector{Float64}()
N = Vector{Int64}()
TheMeans = zeros(100,2)
for j = 1:100
    global counter, nextcol
    ms = 3
    lw = 1.5
    graph = graphs[j]
    F = matread(pathtoFB100*"$graph.mat")
    A = F["A"]
    n = size(A,1)
    push!(N,n)
    println("$j \t $graph $n")

    M = matread("fb_output/$(graph)_output.mat")
    Ps = M["Ps"]
    Objectives = M["Objectives"]
    Sizes = M["Sizes"]
    EdgeDens = M["EdgeDens"]
    MeanDegree = M["MeanDegree"]
    MaxDegree = M["MaxDegree"]
    Times = M["Times"]

    P = Ps[1:end-1]
    # Save mean degree scores for p = 1 and p = 1.05 solution
    TheMeans[j,:] = [MeanDegree[4], MeanDegree[end]]
    if MeanDegree[end] > MeanDegree[4]
        counter += 1
        plot(P,MeanDegree[1:end-1],grid = false, legend = false)
        # savefig("Figures/$graph)_mean_degree.pdf")
    end

    color = :gray
    s1 = 250
    s2 = 200
    if Sizes[end-1]> 3.4*Sizes[1]
        color = RGBA(C[nextcol,1],C[nextcol,2],C[nextcol,3],1)
        plot!(Psize,P,Sizes[1:end-1],color = color, size = (s1,s2), markershape = :circle, markerstrokecolor = color,markersize = ms, xlabel = "p",ylabel = "|S|", linewidth = lw)
        plot!(Pmax,P,MaxDegree[1:end-1],color = color, size = (s1,s2), markershape = :circle, markerstrokecolor = color, markersize = ms,xlabel = "p",ylabel  = "Maximum Degree", linewidth = lw)
        plot!(Pdensity,P,EdgeDens[1:end-1],color = color, size = (s1,s2), markershape = :circle, markerstrokecolor = color,markersize = ms, xlabel = "p",ylabel  = "Edge density", linewidth = lw)
        plot!(Pmean,P,MeanDegree[1:end-1],color = color, size = (s1,s2), markershape = :circle, markerstrokecolor = color, markersize = ms,xlabel = "p",ylabel = "Mean Degree", linewidth = lw)
        nextcol += 1
    else
        lw =1
        al = 0.1
        plot!(Psize,P,Sizes[1:end-1],color = color,size = (s1,s2),  alpha = al, xlabel = "p",ylabel = "Subgraph Size", linewidth = lw)
        plot!(Pmax,P,MaxDegree[1:end-1],color = color,size = (s1,s2),  alpha = al, xlabel = "p",ylabel  = "Maximum Degree", linewidth = lw)
        plot!(Pdensity,P,EdgeDens[1:end-1],color = color, size = (s1,s2), alpha = al, xlabel = "p",ylabel  = "Edge Density", linewidth = lw)
        plot!(Pmean,P,MeanDegree[1:end-1],color = color, size = (s1,s2), alpha = al, xlabel = "p",ylabel = "Average Degree", linewidth = lw)
    end

    if MeanDegree[5] > MeanDegree[4]
        counter += 1
        push!(xx,P[5])
        push!(yy,MeanDegree[5])
    end

    # Runtime is scatter plot
    scatter!(Prun,vec([n n n]),Times[1:3],size = (300,250),color = :blue, markerstrokecolor = :blue,label = "", xlabel = "n",ylabel = "Runtime (seconds)")
    scatter!(Prun,[n],[Times[4]],color = :red, markerstrokecolor = :red,label = "")
    scatter!(Prun,vec([n n n n]),Times[5:8],color = :gray, markerstrokecolor = :gray,label = "")

end

scatter!(Pmean,xx,yy,markershape = :x,color = :black, markersize = 4)

## Save them

plot!(Pmean)
savefig("Figures/Pmean2.pdf")

plot!(Pmax)
savefig("Figures/Pmax2.pdf")

plot!(Psize)
savefig("Figures/Psize2.pdf")

plot!(Pdensity)
savefig("Figures/Pdensity2.pdf")

plot!(Prun)
scatter!(Prun,[-1],[-1],color = :blue, legend = :topleft, ylim = [0,6.5],markerstrokecolor = :blue, xlabel = "n",label = "p<1",ylabel = "Runtime (seconds)")
scatter!(Prun,[-1],[-1],color = :red, markerstrokecolor = :red,label = "p=1")
scatter!(Prun,[-1],[-1],color = :gray, markerstrokecolor = :gray,label = "p>1")

savefig("Figures/Prun.pdf")

## prun

max_stars_repo_name: nveldt/GenMeanDSG

max_issues_repo_name: nveldt/GenMeanDSG

max_forks_repo_name: nveldt/GenMeanDSG


text:
iterresultsname(λ⁺, ::Val{:ols}) = "ans_$(PARAM[:crspdatalabel])" *
  boundrangestr() *
  (PARAM[:placebovol] ? "_SHUFFLED" : "") *
  "_$(PARAM[:measuretype])" *
  "_$(PARAM[:itersolvetype])" *
  "_$(PARAM[:iterresultssuffix])" *
  "_$(Int(round(λ⁺)))" *
  "_$(Dates.format(now(),"yyyymmdd_HHMM"))"
iterresultsname(::Val{:ols}) = "ans_$(PARAM[:crspdatalabel])" *
  boundrangestr() *
  (PARAM[:placebovol] ? "_SHUFFLED" : "") *
  "_$(PARAM[:measuretype])" *
  "_$(PARAM[:itersolvetype])" *
  "_$(PARAM[:iterresultssuffix])" *
  "_$(Dates.format(now(),"yyyymmdd_HHMM"))"

iterresultsname(λ⁺) = iterresultsname(λ⁺, Val(PARAM[:itermodel]))
iterresultsname() = iterresultsname(Val(PARAM[:itermodel]))

iterresultsname(::Val{:bootstrap}) = iterresultsname(Val(:ols))
iterresultsname(λ⁺, ::Val{:bootstrap}) = iterresultsname(λ⁺, Val(:ols))

#optresultsname() = optresultsname(Val(PARAM[:itermodel]))
#primary function for estimating the measure using optimizations
#note that this still uses the same volume parts construct as flux
function formiterregression(panel::DataFrame, ms::MeasureSpec,
  IterModel::Val = Val(PARAM[:itermodel]);
  returns::DataFrame, funds,
  estimationmindate= PARAM[:globalmindate],
  estimationmaxdate= PARAM[:globalmaxdate])

  #issorted(panel, [:permno, :date]) || error("panel must be sorted")
  sort!(panel, :date)

  #bound the dates
  #the lag of the first date is generally expected to be missing
  #if there is no cutoff then it shouldn't matter, but if the dates are constrained then it will
  measuredatekey = Symbol(PARAM[:crspfrequencysuffix],PARAM[:crsptypesuffix])
  if haskey(estimationmindate, measuredatekey)
    mindate=estimationmindate[measuredatekey]
    @info "overriding mindate with mindate=$mindate"
    panel = panel[panel.date .≥ mindate,:]
    for FRLw ∈ Fweights(ms, :FRLw)
      panel[panel.date .== minimum(panel.date), FRLw] .= missing
    end
  end
  if haskey(estimationmaxdate, measuredatekey)
    panel = panel[panel.date .≤ estimationmaxdate,:]
  end

  #in contrast to other methods, drop any column where we don't have
  #the lag weight (improves performance I think)

  #conditionforturing!(panel::SubDataFrame, ms::MeasureSpec)

  zs::ZSpec = condition!(panel, ms)
  #=
  for p ∈ propertynames(zs)
    (p===:originalms) && continue
    prop = getproperty(zs, p)
    if length(propertynames(prop)) == 0
      println("$p: $prop")
    else
      for subp ∈ propertynames(prop)
        subprop = getproperty(prop, subp)
        println("$subp: $subprop")
      end
    end
  end

  throw("stop")=#

  #testing
  if PARAM[:testmeasure]
    testiter(panel, zs)
  end

  if PARAM[:refreshmeasure]
    λ⁺, results::DataFrame = itermodel(panel, zs, IterModel; funds)
    destandardize!(results, zs, intercept=false)
    results = computemeasurestatistics(panel, results, ms, returns)

    #this will be the name of the results file
    resultname::String = iterresultsname(λ⁺)

    results |> CSV.write("$(PARAM[:outputpath])\\$(resultname).csv")

  end

  return nothing
end

function computemeasurestatistics(panel, results, ms::MeasureSpec, returns;
    regressiontype::Symbol = PARAM[:iterregressiontype],
    volscale::Float64 = PARAM[:volscale])

  Fw::Vector{Symbol} = Fweights(ms, :Fw)
  FRtw::Vector{Symbol} = Fweights(ms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(ms, :FRLw)
  FW::Vector{Symbol} = ms.FW
  FWgroup = Fgroupedcontrols(ms)

  #construct the model parts using cpu only data types and no controls
  p = modelparts(panel, Float64, Vector{Float64}, Matrix{Float64},
    weightdata = (;Fw, FRtw, FRLw), Fvol=ms.Fvol,
    ;FW, FWgroup, ms)

  dims = p.dims
  ws = p.dat[:Fw]
  Rtws = p.dat[:FRtw]
  RLws = p.dat[:FRLw]
  W = p.dat[:FW]
  Wgroup = p.dat[:FWgroup]
  ts::Vector{Int} = p.ts
  v = p.dat[:Fvol]

  #construct simplest possible structures for computing predicted volumes
  @info "Computing raw controls for output file"
  Xvraw = XYIterControl(ws, Rtws, RLws, v, ts) #raw implies no controls or transformations
  @info "Computing actual controls for output file"
  Xv = XYIterControl(ws, Rtws, RLws, v, ts, W=W, Wgroup=Wgroup)
  Θ = VolumePartsIter(dims.T, dims.K, ts, Float64, Matrix{Float64}, Vector{Float64})

  #read in the estimations from the results file
  Fcoefs::Vector{Symbol} = findcols(results, :asset)
  Fzcoefs::Vector{Symbol} = findcols(results, :z)
  Fgcoefs::Vector{Symbol} = findcols(results, :growth)
  Fvcoefs::Vector{Symbol} = (s->Symbol(replace(string(s), r"^A_"=>"V_"))).(Fcoefs)
  Fvscoefs::Vector{Symbol} = (s->Symbol(replace(string(s), r"^A_"=>"Vs_"))).(Fcoefs)
  A = Matrix{Float64}(results[!, Fcoefs])
  Z = Matrix{Float64}(results[!, Fzcoefs])
  G = Matrix{Float64}(results[2:end, Fgcoefs])

  #now follow the sequence for estimating the volume, minus the actual estimation
  someones = ones(Float64, dims.K, 1)
  prodmatG = cumprodbyrow(G' |> Matrix{Float64})
  Ã = hcat(someones, prodmatG) |> Matrix{Float64}

  @assert Xvraw.xsection.xM === nothing
  RHSraw::Matrix{Float64} = projectbyt(Ã, Xvraw.xsection) #no controls or transofmrations
  RHS::Matrix{Float64} = hcat(projectbyt(Ã, Xv.xsection), Xv.W̃)

  #now read in the coefficients from the results file
  #we can also compare to the A in the results file as a check
  local β::Vector{Float64}=Vector{Float64}(undef, dims.K)
  if regressiontype === :none
    noregressionbase::Float64 = PARAM[:iternoregressionbase]
    β .= ones(Float64, dims.K) .* noregressionbase .* A[1,:] ./ Z[1,:] .* 2
  else
    β .= reg(RHS, Xv.ṽ, Val(regressiontype))[1:dims.K] .* 2
  end
  all(β .≈ A[1,:]) || throw(
    "Fcoefs: $Fcoefs
    A[1,:]: $(A[1,:])
    β: $β")
  @assert all(β' .* Ã' .≈ A)

  @info "Computing the volume for each strategy/stock/date"
  #now compute the volume for each strategy x stock x date
  volmat::Matrix{Float64} = RHSraw .* β'

  RHSpurchases = projectpurchasesbyt(Ã, Xvraw.xsection) #no controls or transofmrations
  purchasemat = RHSpurchases .* β'
  @assert abs.(purchasemat) ≈ volmat

  @assert size(volmat) === (length(ts), dims.K)
  @assert issorted(results, :date)
  alldates = (t->results.date[t]).(ts)
  volbystock::DataFrame = hcat(
    DataFrame(date=alldates),
    DataFrame(volmat, Fvcoefs),
    DataFrame(purchasemat, Fvscoefs),
    )

  #aggregate over all stocks
  strategyvol = combine(groupby(volbystock, :date),
    Fvcoefs .=> sum,
    Fvscoefs .=> sum, renamecols=false)

  aggvol = combine(groupby(panel, :date),
    :date => last => :date,
    :Wdvol => ((v)->sum(skipmissing(v))) => :desc_Wvolume,
    :dvolscaled => ((v)->sum(skipmissing(v))) => :desc_volume,
    :dvolscaled => ((v)->std(skipmissing(v))) => :desc_vol_std
    )

  stats = leftjoin(strategyvol, aggvol, on=:date)
  sort!(stats, :date)

  #compute the sum of all market caps in each quantile, differentiating long and short
  @info "Computing market cap statistics"
  Fmcgrosss = (s->Symbol(:mcgross_, s)).(Fw)
  Fmclongs = (s->Symbol(:mclong_, s)).(Fw)
  Fmcshorts = (s->Symbol(:mcshort_, s)).(Fw)
  Fmctotal = :mctotal
  aggmc = combine(groupby(panel,:date),
    :date => last => :date,
    :mc => ((x)->sum(x) |> (x)->x/volscale) => :mctotal,
    zip(Fw, [:mc for i ∈ 1:dims.K]) .|> collect .=>
      ((w, mc) -> (w .> zero(Float64)) .* mc |> sum |> (x)->x/volscale) .=>
      Fmclongs,
    zip(Fw, [:mc for i ∈ 1:dims.K]) .|> collect .=>
      ((w, mc) -> (w .< zero(Float64)) .* mc |> sum |> (x)->x/volscale) .=>
      Fmcshorts,)
  transform!(aggmc, zip(Fmclongs, Fmcshorts) .|> collect .=>
    ((l,s)-> l + s) .=>
    Fmcgrosss)

  stats = outerjoin(stats, aggmc, on=:date, validate=(true,true))
  sort!(stats, :date)

  #started working on showing the controls but no time
  #create the control names
  dtstrings = (d->Dates.format(d, "_yy_m_d")).(alldates|> unique |> sort!)
  FWgroupnames::Vector{Vector{Symbol}} = (dtstr->(Symbol.(FWgroup,dtstr))).(dtstrings)
  #FWgroupnames::Vector{Vector{Symbol}} = (F->(Symbol.(F,dtstrings))).(FWgroup)
  #FWgroupnames::Vector{Vector{Symbol}} = dtstrings .|> (dtstr->(F->(Symbol(F,dtstr)).(FWgroup)))

  #create the expanded full control matrix
  #NOTE: below commented code is useful for verificaiton, but it doesn't check for singular
  #columns so I don't use it anymore
  #=@info "Creating expanded control matrix"
  @time begin
    Wdf = ((xM,Fnames)->DataFrame(xM.W,Fnames)).(Xv.xsection.xM, FWgroupnames)
    Wgroupfull = foldl((df1,df2)->vcat(df1,df2,cols=:union), Wdf)
    for c ∈ eachcol(Wgroupfull)
      c .= c .|> (x)->coalesce(x, 0.0)
    end
    expanded::Matrix{Float64} = Wgroupfull |> Matrix{Float64}
  end

  @info "Verifying expanded control matrix"
  @time @assert (expanded .≈ expandcontrolsbyt(Wgroup, ts)) |> all=#

  @info "Checking groupped regression results match full regression results"
  expanded::Matrix{Float64}, controlsused = expandcontrolsbyt(Wgroup, ts) |> (nt
    )->(nt.expanded, nt.Wgroupcolsused)

  #run the full regression to get the controls
  #while we are at it, we can check that the regression ran as expected
  if regressiontype === :none
    RHSfull = hcat(expanded, Xv.W)
    βfull = svd(Symmetric(RHSfull' * RHSfull))\(RHSfull'*v) |> vec
    βgroup = βfull[1:(end-size(Xv.W,2))]
    fullerr = abs.(sum(RHS[:,1:dims.K],dims=2) .+ RHSfull*βfull .- Xv.v)
    βW̃ = (size(Xv.W̃,2) > 0) ? reg(Xv.W̃, Xv.ṽ, Val{:svd}()) : TV()
    βproj = [ones(dims.K); βW̃]
  else
    RHSfull = hcat(RHSraw,expanded, Xv.W)
    βfull =  cholesky!(Symmetric(RHSfull' * RHSfull))\(RHSfull'*v) |> vec
    βgroup = βfull[(size(RHSraw,2) + 1):(end-size(Xv.W,2))]
    fullerr = abs.(RHSfull*βfull .- Xv.v)
    βproj = reg(RHS, Xv.ṽ, Val(regressiontype))
  end

  #the below check verifies that the residuals are projected as expected
  projerr = abs.(RHS*βproj .- Xv.ṽ)
  parterr = abs.(sum(RHSraw,dims=2))
  @assert (fullerr ≈ projerr) "
    fullerr[1:5]: $(fullerr[1:5])
    projerr[1:5]: $(projerr[1:5])
    parterr[1:5]: $(parterr[1:5])
    Median err: Full: $(median(fullerr)) proj: $(median(projerr))
    Mean err: Full: $(mean(fullerr)) proj: $(mean(projerr))
    Median dev: Full: $(median(fullerr .- projerr))
    mean dev: Full: $(mean(fullerr .- projerr))
    Max abs dif: $(maximum(fullerr .- projerr))
    SSE err: Full: $(sum(fullerr.^2)) Proj: $(sum(projerr.^2))"

  @info "verificaiton of projection matrix passed"
  #reshape the control coefficients into time series with as much care as possible
  #@assert length(βgroup) == length(unique(ts)) * size(Wgroup,2)
  βgroupconformed = Vector{MFloat64}(undef, length(unique(ts)) * size(Wgroup,2))
  βgroupconformed[controlsused] .= βgroup
  βgroupmat = reshape(βgroupconformed, size(Wgroup,2), dims.T-1)'
  @assert (vec(βgroupmat') .≡ βgroupconformed) |> all

  rhsdf::DataFrame = DataFrame(date=alldates|>deepcopy, resid=projerr)
  controldf = combine(groupby(rhsdf,:date),
      :date => last => :date,
      :resid=> ((v)->std(skipmissing(v)))=>:se_resid)
  @assert issorted(controldf, :date)
  for (c,F) ∈ enumerate(FWgroup)
    #controldf[!, F] = [missing; βgroupmat[:,c]]
    controldf[!, F] = βgroupmat[:,c]
  end

  @info "combining controls, extra stats, and returns with the results dataframe"

  stats = outerjoin(stats, controldf, on=:date, validate=(true,true))
  sort!(stats, :date)
  results = leftjoin(results, stats, on=:date, validate=(true,true))

  results = leftjoin(results, returns, on=:date, validate=(true,true))
  sort!(results, :date)

  #patch in data from the return df
  #@assert all(completecases(results[2:end, [:volume; :Wvolume]], [:volume; :Wvolume]))
  if !all(completecases(results[2:end, [:volume; :Wvolume]], [:volume; :Wvolume]))
    results |> CSV.write("$(PARAM[:testpath])\\results_error.csv")
    @assert all(completecases(results[2:end, [:volume; :Wvolume]], [:volume; :Wvolume]))
  end

  resultsdescvolumes = view(results, completecases(results, [:desc_volume, :desc_Wvolume]), :)
  @assert all(resultsdescvolumes[:, :volume] .≈ resultsdescvolumes[:, :desc_volume] )
  @assert all(resultsdescvolumes[:, :Wvolume] .≈ resultsdescvolumes[:, :desc_Wvolume])

  #the data from the returns file are more compelte
  select!(results, Not([:volume, :Wvolume]))

  for ξ ∈ ms.ξs
    @assert all(completecases(results[2:end, [ξ.Fξ]]))
    rename!(results, Dict(ξ.Fξ=>Symbol(:R_, ξ.Fξ)))
  end

  #####now merge in the returns data

  results = ivoltests(panel, results, ms, returns)

  return  results
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module GTD

using DataFrames
using SQLite
using LibCURL
using Base.Markdown
using DataStreams
using Lazy

#include("this.jl")
#include("that.jl")

export add_gtd_db, collect_stuff, show_stuff, process_inbox

# first step: create our database
function add_gtd_db(path::AbstractString, fname::AbstractString)
  cd(path)
  db = SQLite.DB(path * "/" * fname * ".db")
  try
    SQLite.execute!(db,
    """create table stuff (item varchar(140), notes varchar(140), time varchar(140), id varchar(36));""")
  finally
    return db
  end
end

function collect_stuff(db::SQLite.DB, item::AbstractString, notes::AbstractString = "")
  SQLite.execute!(db, """insert into stuff (item, notes, time, id) values ('"""
                  * item * "', '" * notes * "', '" * string(now()) * "', '"
                  * string(Base.Random.uuid4()) * "');")
end

function show_stuff(db::SQLite.DB)
  @>> """select * from stuff""" SQLite.query(db) DataFrame
end

function process_inbox(db::SQLite.DB)
  stuff = @>> """select * from stuff""" SQLite.query(db) DataFrame
  for i in 1:nrow(stuff)
    process_item(stuff[i,:])
  end
end

function process_item(stuff::DataFrame)
  @show stuff
  # is it actionable
    # yes
      # can be done under 2 min && do
      # requires multiple steps && project
      # else make action
        # can it be delegated? delegate
        # can it be postponed? postpone
    # no
      # is it someday/maybe
      # is it reference
end

#=
function process(item)
  if actionable
    process_actionable(item)
  else
    process_other(item)
  end
end

function process_actionable(item)
  if project
    create_project(item)
  else
    create_action(item)
  end
end

function process_other(item)
  if smdmb
    create_someday_maybe(item)
  elseif reference
    create_reference(item) # references must be tieable to projects
  end
end=#

  # end of module
end


text:
# API (Two-stage) #
# ========================== #
"""
    instantiate(stochasticmodel::StochasticModel{2},
                scenarios::Vector{<:AbstractScenario};
                instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                optimizer = nothing;
                defer::Bool = false,
                kw...)

Instantiate a two-stage stochastic program using the model definition stored in the two-stage `stochasticmodel`, and the given collection of `scenarios`. Optionally, supply an `optimizer`. If no explicit `instantiation` is provided, the structure is induced by the optimizer. The structure is `Deterministic` by default.
"""
function instantiate(sm::StochasticModel{2},
                     scenarios::Vector{<:AbstractScenario};
                     instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                     optimizer = nothing,
                     defer::Bool = false,
                     direct_model::Bool = false,
                     kw...)
    sp = StochasticProgram(parameters(sm.parameters[1]; kw...),
                           parameters(sm.parameters[2]; kw...),
                           scenarios,
                           instantiation,
                           optimizer)
    # Generate model recipes
    sm.generator(sp)
    # Return deterministic equivalent as JuMP model if direct model has been requested
    if direct_model
        optimizer(sp) isa MOI.AbstractOptimizer || error("MOI optimizer required to create direct model.")
        model = direct_model(optimizer_constructor(sp))
        model.ext[:decisions] = IgnoreDecisions()
        _generate_deterministic_equivalent!(sp, model)
        return model
    end
    # Generate stochastic program if not deferred
    if !defer
        generate!(sp)
    end
    return sp
end
"""
    instantiate(stochasticmodel::StochasticModel{2};
                instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                optimizer = nothing,
                scenario_type::Type{S} = Scenario,
                defer::Bool = false,
                kw...) where S <: AbstractScenario

Instantiate a deferred two-stage stochastic program using the model definition stored in the two-stage `stochasticmodel` over the scenario type `S`. Optionally, supply an `optimizer`. If no explicit `instantiation` is provided, the structure is induced by the optimizer. The structure is `Deterministic` by default.
"""
function instantiate(sm::StochasticModel{2};
                     instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                     optimizer = nothing,
                     scenario_type::Type{S} = Scenario,
                     kw...) where S <: AbstractScenario
    sp = StochasticProgram(parameters(sm.parameters[1]; kw...),
                           parameters(sm.parameters[2]; kw...),
                           scenario_type,
                           instantiation,
                           optimizer)
    # Generate model recipes
    sm.generator(sp)
    return sp
end
"""
    instantiate(stochasticmodel::StochasticModel,
                scenarios::Vector{<:AbstractScenario};
                instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                optimizer = nothing,
                defer::Bool = false,
                kw...)

Instantiate a stochastic program using the model definition stored in `stochasticmodel`, and the given collection of `scenarios`. Optionally, supply an `optimizer`. If no explicit `instantiation` is provided, the structure is induced by the optimizer. The structure is `Deterministic` by default.
"""
function instantiate(sm::StochasticModel{N},
                     scenarios::NTuple{M,Vector{<:AbstractScenario}};
                     instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                     optimizer = nothing,
                     defer::Bool = false,
                     kw...) where {N,M}
    M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
    stages = ntuple(Val(N)) do i
        Stage(parameters(sm.parameters[i]; kw...))
    end
    sp = StochasticProgram(stages,
                           scenarios,
                           instantiation,
                           optimizer)
    # Generate model recipes
    sm.generator(sp)
    # Generate stochastic program if not deferred
    if !defer
        generate!(sp)
    end
    return sp
end
"""
    instantiate(stochasticmodel::StochasticModel,
                sampler::AbstractSampler,
                n::Integer;
                instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                optimizer = nothing,
                defer::Bool = false,
                kw...)

Generate a sampled instance of size `n` using the model stored in the two-stage `stochasticmodel`, and the provided `sampler`. Optionally, supply an `optimizer`. If no explicit `instantiation` is provided, the structure is induced by the optimizer. The structure is `Deterministic` by default.

"""
function instantiate(stochasticmodel::StochasticModel{2},
                     sampler::AbstractSampler{S},
                     n::Integer;
                     instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                     optimizer = nothing,
                     defer::Bool = false,
                     kw...) where S <: AbstractScenario
    # Create new stochastic program instance
    sp = StochasticProgram(parameters(stochasticmodel.parameters[1]; kw...),
                           parameters(stochasticmodel.parameters[2]; kw...),
                           S,
                           instantiation,
                           optimizer)
    # Generate model recipes
    stochasticmodel.generator(sp)
    # Sample n scenarios
    add_scenarios!(sp, n) do
        return sample(sampler, 1/n)
    end
    # Generate stochastic program if not deferred
    if !defer
        generate!(sp)
    end
    # Return the sample instance
    return sp
end
function instantiate(stochasticmodel::StochasticModel{2},
                     sampler::AbstractSampler{Scenario},
                     n::Integer;
                     instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
                     optimizer = nothing,
                     defer::Bool = false,
                     kw...)
    # Get concrete scenario type
    S = typeof(sampler())
    # Create new stochastic program instance
    sp = StochasticProgram(parameters(stochasticmodel.parameters[1]; kw...),
                           parameters(stochasticmodel.parameters[2]; kw...),
                           S,
                           instantiation,
                           optimizer)
    # Generate model recipes
    stochasticmodel.generator(sp)
    # Sample n scenarios
    add_scenarios!(sp, n) do
        return sample(sampler, 1/n)
    end
    # Generate stochastic program if not deferred
    if !defer
        generate!(sp)
    end
    # Return the sample instance
    return sp
end
"""
    optimize!(stochasticprogram::StochasticProgram; crash::AbstractCrash = Crash.None(); kw...)

Optimize the `stochasticprogram` in expectation. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown. An optional crash procedure can be set to warm-start.

## Examples

The following solves the stochastic program `sp` using the L-shaped algorithm.

```julia
set_optimizer(sp, LShaped.Optimizer)
set_optimizer_attribute(sp, MasterOptimizer(), GLPK.Optimizer)
set_optimizer_attribute(sp, SubproblemOptimizer(), GLPK.Optimizer)
optimize!(sp);

# output

L-Shaped Gap  Time: 0:00:02 (6 iterations)
  Objective:       -855.8333333333358
  Gap:             0.0
  Number of cuts:  8
  Iterations:      6
```

The following solves the stochastic program `sp` using GLPK on the extended form.

```julia
using GLPK

set_optimizer(sp, GLPK.Optimizer)
optimize!(sp)

```

See also: [`VRP`](@ref)
"""
function JuMP.optimize!(stochasticprogram::TwoStageStochasticProgram; crash::AbstractCrash = Crash.None(), cache::Bool = false, kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Ensure stochastic program has been initialized at this point
    if deferred(stochasticprogram)
        generate!(stochasticprogram)
    end
    # Crash initial decision
    x₀ = crash(stochasticprogram)
    # Switch on structure and solver type
    optimize!(structure(stochasticprogram), optimizer(stochasticprogram), x₀; kw...)
    # Cache solution (if requested)
    if cache
        cache_solution!(stochasticprogram, structure(stochasticprogram), optimizer(stochasticprogram))
    end
    return nothing
end
"""
    termination_status(stochasticprogram::StochasticProgram)

Return the reason why the solver of the `stochasticprogram` stopped.
"""
function JuMP.termination_status(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.TerminationStatus())::MOI.TerminationStatusCode
end
"""
    termination_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

Return the reason why the solver of the node at stage `stage` and scenario `scenario_index`
stopped.
"""
function JuMP.termination_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.TerminationStatus())
    return MOI.get(stochasticprogram, attr)::MOI.TerminationStatusCode
end
"""
    termination_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)

Return the reason why the solver of scenario `scenario_index` stopped.
"""
function JuMP.termination_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)
    return terminationstatus(stochasticprogram, 2, scenario_index)
end
"""
    raw_status(stochasticprogram::StochasticProgram)

Return the reason why the solver of the `stochasticprogram` stopped in its own words.
"""
function JuMP.raw_status(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.RawStatusString())
end
"""
    raw_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

Return the reason why the solver of the node at stage `stage` and scenario `scenario_index`
stopped in its own words.
"""
function JuMP.raw_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.RawStatusString())
    return MOI.get(stochasticprogram, attr)
end
"""
    raw_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)

Return the reason why the solver of scenario `scenario_index` stopped in its own words.
"""
function JuMP.raw_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)
    return raw_status(stochasticprogram, 2, scenario_index)
end
"""
    primal_status(stochasticprogram::StochasticProgram; result::Int = 1)

Return the status of the most recent primal solution of the solver of
the `stochasticprogram`.
"""
function JuMP.primal_status(stochasticprogram::StochasticProgram; result::Int = 1)
    return MOI.get(stochasticprogram, MOI.PrimalStatus(result))::MOI.ResultStatusCode
end
"""
    primal_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)

Return the status of the most recent primal solution of the solver of
the node at stage `stage` and scenario `scenario_index`.
"""
function JuMP.primal_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.PrimalStatus(result))
    return MOI.get(stochasticprogram, attr)::MOI.ResultStatusCode
end
"""
    primal_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)

Return the status of the most recent primal solution of the solver of
the scenario `scenario_index`.
"""
function JuMP.primal_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)
    return primal_status(stochasticprogram, 2, scenario_index; result = result)
end
"""
    dual_status(stochasticprogram::StochasticProgram; result::Int = 1)

Return the status of the most recent dual solution of the solver of
the `stochasticprogram`.
"""
function JuMP.dual_status(stochasticprogram::StochasticProgram; result::Int = 1)
    return MOI.get(stochasticprogram, MOI.DualStatus(result))::MOI.ResultStatusCode
end
"""
    dual_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)

Return the status of the most recent dual solution of the solver of
the node at stage `stage` and scenario `scenario_index`.
"""
function JuMP.dual_status(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.DualStatus(result))
    return MOI.get(stochasticprogram, attr)::MOI.ResultStatusCode
end
"""
    dual_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)

Return the status of the most recent dual solution of the solver of
scenario `scenario_index`.
"""
function JuMP.dual_status(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)
    return dual_status(stochasticprogram, 2, scenario_index; result = result)
end
"""
    solve_time(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

If available, returns the solve time reported by the solver of
the `stochasticprogram`.
"""
function JuMP.solve_time(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.SolveTime())
end
"""
    solve_time(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

If available, returns the solve time reported by the solver of
the node at stage `stage` and scenario `scenario_index`.
"""
function JuMP.solve_time(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.SolveTime())
    return MOI.get(stochasticprogram, attr)
end
"""
    solve_time(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)

If available, returns the solve time reported by the solver of
scenario `scenario_index`.
"""
function JuMP.solve_time(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)
    return solve_time(stochasticprogram, 2, scenario_index)
end
"""
    optimal_decision(stochasticprogram::StochasticProgram)

Return the optimal first stage decision of `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::StochasticProgram)
    if termination_status(stochasticprogram) == MOI.OPTIMIZE_NOT_CALLED
        throw(OptimizeNotCalled())
    end
    return JuMP.value.(all_decision_variables(stochasticprogram, 1))
end
"""
    optimal_decision(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

Return the optimal decision in the node at stage `stage` and scenario `scenario_index`
of `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)
    if termination_status(stochasticprogram, stage, scenario_index) == MOI.OPTIMIZE_NOT_CALLED
        throw(OptimizeNotCalled())
    end
    return JuMP.value.(all_decision_variables(stochasticprogram, stage), scenario_index)
end
"""
    optimal_decision(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)

Return the optimal decision in scenario `scenario_index` of the two-stage `stochasticprogram`,
after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)
    return optimal_decision(stochasticprogram, 2, scenario_index)
end
"""
    optimal_recourse_decision(stochasticprogram::StochasticProgram, scenario_index::Integer)

Return the optimal recourse decision in the final stage of `stochasticprogram` in
the scenario `scenario_index`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_recourse_decision(stochasticprogram::StochasticProgram{N}, scenario_index::Integer) where N
    return optimal_decision(stochasticprogram, N, scenario_index)
end
"""
    optimize!(stochasticmodel::StochasticModel, sampler::AbstractSampler; crash::AbstractCrash = Crash.None(), kw...)

Approximately optimize the `stochasticmodel` when the underlying scenario distribution is inferred by `sampler`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function JuMP.optimize!(stochasticmodel::StochasticModel,
                        sampler::AbstractSampler;
                        crash::AbstractCrash = Crash.None(),
                        kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Crash initial decision
    x₀ = crash(stochasticmodel, sampler)
    # Switch on solver type
    return _optimize!(stochasticmodel, sampler, optimizer(stochasticmodel), x₀; kw...)
end
function _optimize!(stochasticmodel::StochasticModel, sampler::AbstractSampler, optimizer::AbstractSampledOptimizer, x₀; kw...)
    # Load optimizer
    load_model!(optimizer, stochasticmodel, sampler, x₀)
    # Run sample-based optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end
function JuMP.termination_status(stochasticmodel::StochasticModel)
    return MOI.get(stochasticmodel, MOI.TerminationStatus())::MOI.TerminationStatusCode
end
function JuMP.objective_sense(stochasticmodel::StochasticModel)
    return MOI.get(stochasticmodel, MOI.ObjectiveSense())::MOI.OptimizationSense
end
"""
    optimal_decision(stochasticmodel::StochasticModel)

Return the optimal first stage decision of `stochasticmodel`, after a call to `optimize!(stochasticmodel)`.
"""
function optimal_decision(stochasticmodel::StochasticModel)
    if termination_status(stochasticmodel) == MOI.OPTIMIZE_NOT_CALLED
        throw(OptimizeNotCalled())
    end
    return optimal_decision(instance(stochasticmodel))
end
"""
    objective_value(stochasticmodel::StochasticModel; result::Int = 1)

Returns the value of the recourse problem after a call to `optimize!(stochasticmodel)`.
"""
function JuMP.objective_value(stochasticmodel::StochasticModel; result::Int = 1)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Switch on termination status
    status = JuMP.termination_status(stochasticmodel)
    if status in AcceptableTermination
        return MOI.get(stochasticmodel, MOI.ObjectiveValue(result))
    else
        if status == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        elseif status == MOI.INFEASIBLE
            @warn "Stochastic program is infeasible"
            return objective_sense(stochasticmodel) == MOI.MAX_SENSE ? -Inf : Inf
        elseif status == MOI.DUAL_INFEASIBLE
            @warn "Stochastic program is unbounded"
            return objective_sense(stochasticmodel) == MOI.MAX_SENSE ? Inf : -Inf
        else
            @warn("Stochastic program could not be solved, returned status: $status")
            return NaN
        end
    end
end
"""
    stage_parameters(stochasticprogram::StochasticProgram, stage::Integer)

Return the parameters at `stage` in `stochasticprogram`.
"""
function stage_parameters(stochasticprogram::StochasticProgram{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return stochasticprogram.stages[stage].parameters
end
"""
    decision(stochasticprogram::StochasticProgram, index::MOI.VariableIndex)

Return the current value of the first-stage decision at `index` of `stochasticprogram`.
"""
function decision(stochasticprogram::StochasticProgram, index::MOI.VariableIndex)
    return decision(structure(stochasticprogram), 1, index)
end
"""
    decision(stochasticprogram::StochasticProgram)

Return the current first-stage decision values of `stochasticprogram`.
"""
function decision(stochasticprogram::StochasticProgram)
    return JuMP.value(all_decision_variables(stochasticprogram))
end
"""
    all_decision_variables(stochasticprogram::StochasticProgram{N}) where N

Returns a stage-wise list of all decisions currently in the `stochasticprogram`. The decisions are
ordered by creation time. Defaults to the first stage.
"""
function all_decision_variables(stochasticprogram::StochasticProgram{N}) where N
    return ntuple(Val{N}()) do s
        return all_decision_variables(stochasticprogram, s)
    end
end
"""
    all_decision_variables(stochasticprogram::StochasticProgram{N}, stage::Integer = 1) where N

Returns a list of all decisions currently in the `stochasticprogram` at `stage`. The decisions are
ordered by creation time.
"""
function all_decision_variables(stochasticprogram::StochasticProgram{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    decisions::NTuple{stage,Decisions} = stochasticprogram.proxy[stage].ext[:decisions]
    return map(all_decisions(decisions[stage])) do index
        DecisionVariable(stochasticprogram, stage, index)
    end
end
"""
    structure(stochasticprogram::StochasticProgram)

Return the underlying structure of the `stochasticprogram`.
"""
function structure(stochasticprogram::StochasticProgram)
    return stochasticprogram.structure
end
"""
    proxy(stochasticprogram::StochasticProgram, stage::Integer)

Return the proxy model of the `stochasticprogram` at `stage`.
"""
function proxy(stochasticprogram::StochasticProgram{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return stochasticprogram.proxy[stage]
end
"""
    Base.getindex(stochasticprogram::StochasticProgram, stage::Integer, name::Symbol)

Returns the decision, or group of decisions, or decision constraint, or group of decision constraints, of the given `name` and `stage` of the `stochasticprogram` which were added to the model. This errors if regular variables or constraints are queried. If so, either annotate the relevant variables with [`@decision`](@ref) or first query the relevant JuMP subproblem and use the regular `[]` syntax.
"""
function Base.getindex(stochasticprogram::StochasticProgram, stage::Integer, name::Symbol)
    obj_dict = object_dictionary(proxy(stochasticprogram, stage))
    if !haskey(obj_dict, name)
        throw(KeyError(name))
    else
        obj = obj_dict[name]
        if obj isa DecisionRef
            return DecisionVariable(stochasticprogram, stage, index(obj))
        elseif obj isa AbstractArray{<:DecisionRef}
            return map(obj) do dvar
                return DecisionVariable(stochasticprogram, stage, index(dvar))
            end
        elseif obj isa ConstraintRef{Model, <:CI{<:DecisionLike}}
            return SPConstraintRef(stochasticprogram, stage, index(obj), obj.shape)
        elseif obj isa AbstractArray{<:ConstraintRef{Model, <:CI{<:DecisionLike}}}
            return map(obj) do cref
                return SPConstraintRef(stochasticprogram, stage, index(cref), cref.shape)
            end
        else
            error("Only decisions and decision constraints can be queried using this syntax. For regular variables and constraints, either annotate the relevant variable with @decision or first query the relevant JuMP subproblem and use the regular `[]` syntax.")
        end
    end
end
"""
    structure_name(stochasticprogram::StochasticProgram)

Return the name of the underlying structure of the `stochasticprogram`.
"""
function structure_name(stochasticprogram::StochasticProgram)
    return structure_name(structure(stochasticprogram))
end
"""
    num_decisions(stochasticprogram::StochasticProgram, stage::Integer = 1)

Return the number of decisions at `stage` in the `stochasticprogram`. Defaults to the first stage.
"""
function num_decisions(stochasticprogram::StochasticProgram{N}, stage::Integer = 1) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return num_decisions(proxy(stochasticprogram, stage), stage)
end
"""
    num_variables(stochasticprogram::StochasticProgram, stage::Integer = 1)

Return the total number of variables at `stage` in the `stochasticprogram`. Defaults to the first stage.
"""
function JuMP.num_variables(stochasticprogram::StochasticProgram{N}, stage::Integer = 1) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return num_variables(proxy(stochasticprogram, stage))
end
"""
    first_stage(stochasticprogram::StochasticProgram)

Return the first stage of `stochasticprogram`.
"""
function first_stage(stochasticprogram::StochasticProgram; optimizer = nothing)
    return first_stage(stochasticprogram, structure(stochasticprogram); optimizer = optimizer)
end
function first_stage(stochasticprogram::StochasticProgram, ::AbstractStochasticStructure; optimizer = nothing)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache, :stage_1)
        stage_one = cache[:stage_1]
        optimizer != nothing && set_optimizer(stage_one, optimizer)
        return stage_one
    end
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    # Generate and cache first stage
    generate_stage_one!(stochasticprogram)
    stage_one = cache[:stage_1]
    optimizer != nothing && set_optimizer(stage_one, optimizer)
    return stage_one
end
"""
    scenario(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

Return the scenario at `scenario_index` of `stochasticprogram` at stage `stage`.
"""
function scenario(stochasticprogram::StochasticProgram{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && error("The first stage does not have scenarios.")
    return scenario(structure(stochasticprogram), stage, scenario_index)
end
"""
    scenario(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)

Return the scenario at `scenario_index` of the two-stage `stochasticprogram`
"""
function scenario(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer) where N
    return scenario(structure(stochasticprogram), 2, scenario_index)
end
"""
    scenarios(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return an array of all scenarios of the `stochasticprogram` at `stage`. Defaults to the second stage.
"""
function scenarios(stochasticprogram::StochasticProgram{N}, stage::Integer = 2) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && error("The first stage does not have scenarios.")
    return scenarios(structure(stochasticprogram), stage)
end
"""
    expected(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return the exected scenario of all scenarios of the `stochasticprogram` at `stage`. Defaults to the second stage.
"""
function expected(stochasticprogram::StochasticProgram{N}, stage::Integer = 2) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && error("The first stage does not have scenarios.")
    return expected(structure(stochasticprogram), stage)
end
"""
    scenario_type(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return the type of the scenario structure associated with `stochasticprogram` at `stage`. Defaults to the second stage.
"""
function scenario_type(stochasticprogram::StochasticProgram{N}, stage::Integer = 2) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && error("The first stage does not have scenarios.")
    return scenario_type(structure(stochasticprogram), stage)
end
"""
    scenario_types(stochasticprogram::StochasticProgram)

Return a stage-wise list of the scenario types of `stochasticprogram`.
"""
function scenario_types(stochasticprogram::StochasticProgram)
    return ntuple(Val(num_stages(stochasticprogram) - 1)) do s
        scenario_type(stochasticprogram, s + 1)
    end
end
"""
    probability(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

Return the probability of the scenario at `scenario_index` in the `stochasticprogram` at `stage` occuring.
"""
function probability(stochasticprogram::StochasticProgram{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    s == 1 && error("The first stage does not have scenarios.")
    return probability(structure(stochasticprogram), stage, scenario_index)
end
"""
    probability(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)

Return the probability of the scenario at `scenario_index` in the two-stage `stochasticprogram` occuring.
"""
function probability(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)
    return probability(structure(stochasticprogram), 2, scenario_index)
end
"""
    stage_probability(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return the probability of any scenario in the `stochasticprogram` at `stage` occuring. A well defined model should return 1. Defaults to the second stage.
"""
function stage_probability(stochasticprogram::StochasticProgram{N}, stage::Integer = 2) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && error("The first stage does not have scenarios.")
    return stage_probability(structure(stochasticprogram), stage)
end
"""
    has_generator(stochasticprogram::StochasticProgram, key::Symbol)

Return true if a problem generator with `key` exists in `stochasticprogram`.
"""
function has_generator(stochasticprogram::StochasticProgram, key::Symbol)
    return haskey(stochasticprogram.generator, key)
end
"""
    generator(stochasticprogram::StochasticProgram, key::Symbol)

Return the problem generator associated with `key` in `stochasticprogram`.
"""
function generator(stochasticprogram::StochasticProgram, key::Symbol)
    return stochasticprogram.generator[key]
end
"""
    subproblem(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)

Return the subproblem at `scenario_index` of the `stochasticprogram` at `stage`.
"""
function subproblem(stochasticprogram::StochasticProgram{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 || error("The first-stage does not have subproblems.")
    return subproblem(structure(stochasticprogram), stage, scenario_index)
end
"""
    subproblem(stochasticprogram::TwoStasgeStochasticProgram, scenario_index::Integer)

Return the subproblem at `scenario_index` of the two-stage `stochasticprogram`.
"""
function subproblem(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer)
    return subproblem(structure(stochasticprogram), 2, scenario_index)
end
"""
    subproblems(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return an array of all subproblems of the `stochasticprogram` at `stage`. Defaults to the second stage.
"""
function subproblems(stochasticprogram::StochasticProgram{N}, stage::Integer = 2) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 || error("The first-stage does not have subproblems.")
    return subproblems(structure(stochasticprogram), stage)
end
"""
    num_subproblems(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return the number of subproblems in the `stochasticprogram` at `stage`. Defaults to the second stage.
"""
function num_subproblems(stochasticprogram::StochasticProgram, stage::Integer = 2)
    stage == 1 && return 0
    return num_subproblems(structure(stochasticprogram), stage)
end
"""
    num_scenarios(stochasticprogram::StochasticProgram, stage::Integer = 2)

Return the number of scenarios in the `stochasticprogram` at `stage`. Defaults to the second stage.
"""
function num_scenarios(stochasticprogram::StochasticProgram, stage::Integer = 2)
    stage == 1 && return 0
    return num_scenarios(structure(stochasticprogram), stage)
end
"""
    num_stages(stochasticprogram::StochasticProgram)

Return the number of stages in `stochasticprogram`.
"""
num_stages(::StochasticProgram{N}) where N = N
"""
    distributed(stochasticprogram::StochasticProgram, s::Integer = 2)

Return true if the `stochasticprogram` is memory distributed at stage `s`. Defaults to the second stage.
"""
distributed(stochasticprogram::StochasticProgram, s::Integer = 2) = distributed(structure(stochasticprogram), s)
"""
    deferred(stochasticprogram::StochasticProgram)

Return true if `stochasticprogram` is not fully generated.
"""
deferred(stochasticprogram::StochasticProgram) = num_scenarios(stochasticprogram) == 0 || deferred(structure(stochasticprogram))
"""
    optimizer_constructor(stochasticmodel::StochasticModel)

Return any optimizer constructor supplied to the `stochasticmodel`.
"""
function optimizer_constructor(stochasticmodel::StochasticModel)
    return stochasticmodel.optimizer.optimizer_constructor
end
"""
    optimizer_constructor(stochasticprogram::StochasticProgram)

Return any optimizer constructor supplied to the `stochasticprogram`.
"""
function optimizer_constructor(stochasticprogram::StochasticProgram)
    return stochasticprogram.optimizer.optimizer_constructor
end
"""
    optimizer_name(stochasticmodel::StochasticModel)

Return the currently provided optimizer type of `stochasticmodel`.
"""
function optimizer_name(stochasticmodel::StochasticModel)
    return optimizer_name(optimizer(stochasticmodel))
end
"""
    optimizer_name(stochasticprogram::StochasticProgram)

Return the currently provided optimizer type of `stochasticprogram`.
"""
function optimizer_name(stochasticprogram::StochasticProgram)
    deferred(stochasticprogram) && return "Uninitialized"
    return optimizer_name(optimizer(stochasticprogram))
end
function optimizer_name(optimizer::MOI.AbstractOptimizer)
    if JuMP.moi_mode(optimizer) != DIRECT && MOIU.state(optimizer) == MOIU.NO_OPTIMIZER
        return "No optimizer attached."
    end
    return JuMP._try_get_solver_name(optimizer)
end
"""
    master_optimizer(stochasticprogram::StochasticProgram)

Return a MOI optimizer using the currently provided optimizer of `stochasticprogram`.
"""
function master_optimizer(stochasticprogram::StochasticProgram)
    return master_optimizer(stochasticprogram.optimizer)
end
"""
    subproblem_optimizer(stochasticprogram::StochasticProgram)

Return a MOI optimizer for solving subproblems using the currently provided optimizer of `stochasticprogram`.
"""
function subproblem_optimizer(stochasticprogram::StochasticProgram)
    return subproblem_optimizer(stochasticprogram.optimizer)
end
"""
    optimizer(stochasticprogram::StochasticProgram)

Return the optimizer attached to `stochasticprogram`.
"""
function optimizer(stochasticprogram::StochasticProgram)
    return stochasticprogram.optimizer.optimizer
end
"""
    optimizer(stochasticmodel::StochasticModel)

Return the optimizer attached to `stochasticmodel`.
"""
function optimizer(stochasticmodel::StochasticModel)
    return stochasticmodel.optimizer.optimizer
end
"""
    get_optimizer_attribute(stochasticprogram::StochasticProgram, name::String)

Return the value associated with the solver-specific attribute named `name` in `stochasticprogram`.

See also: [`set_optimizer_attribute`](@ref), [`set_optimizer_attributes`](@ref).
"""
function get_optimizer_attribute(stochasticprogram::StochasticProgram, name::String)
    return get_optimizer_attribute(stochasticprogram, MOI.RawParameter(name))
end
"""
    get_optimizer_attribute(stochasticprogram::StochasticProgram, attr::MOI.AbstractOptimizerAttribute)

Return the value of the solver-specific attribute `attr` in `stochasticprogram`.

See also: [`set_optimizer_attribute`](@ref), [`set_optimizer_attributes`](@ref).
"""
function get_optimizer_attribute(stochasticprogram::StochasticProgram, attr::MOI.AbstractOptimizerAttribute)
    return MOI.get(stochasticprogram, attr)
end
# ========================== #

# Setters
# ========================== #
"""
    update_known_decisions!(stochasticprogram::Stochasticprogram)

Update all known decision values in the first-stage of `stochasticprogram`.
"""
function update_known_decisions!(stochasticprogram::StochasticProgram)
    update_known_decisions!(structure(stochasticprogram))
end
"""
    update_decisions!(stochasticprogram::Stochasticprogram, stage::Integer, scenario_index::Integer)

Update all known decision values in the node at stage `stage` and scenario `scenario_index` of `stochasticprogram`.
"""
function update_decisions!(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer)
    update_known_decisions!(structure(stochasticprogram), stage, scenario_index)
end
"""
    set_optimizer(stochasticmodel::StochasticModel, optimizer)

Set the optimizer of the `stochasticmodel`.
"""
function JuMP.set_optimizer(stochasticmodel::StochasticModel, optimizer)
    set_optimizer!(stochasticmodel.optimizer, optimizer)
    if stochasticmodel.optimizer.optimizer isa AbstractSampledOptimizer
        return nothing
    end
    # Default to SAA
    set_optimizer(stochasticmodel, SAA.Optimizer)
    set_optimizer_attribute(stochasticmodel, InstanceOptimizer(), optimizer)
    return nothing
end
"""
    set_optimizer(stochasticprogram::StochasticProgram, optimizer)

Set the optimizer of the `stochasticprogram`.
"""
function JuMP.set_optimizer(stochasticprogram::StochasticProgram, optimizer)
    if !supports_structure(optimizer(), structure(stochasticprogram))
        @warn "The provided optimizer does not support the underlying structure of the stochastic program. Consider reinstantiating (see `instantiate`) the stochastic program with the optimizer instead."
        return nothing
    end
    set_optimizer!(stochasticprogram.optimizer, optimizer)
    master_opt = master_optimizer(stochasticprogram)
    if master_opt != nothing
        set_master_optimizer!(structure(stochasticprogram), master_opt)
    end
    sub_opt = subproblem_optimizer(stochasticprogram)
    if sub_opt != nothing
        set_subproblem_optimizer!(structure(stochasticprogram), sub_opt)
    end
    return nothing
end
"""
    set_optimizer_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)

Sets solver-specific attribute identified by `name` to `value` in the `stochasticprogram`.

See also: [`get_optimizer_attribute`](@ref)
"""
function JuMP.set_optimizer_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)
    return set_optimizer_attribute(stochasticprogram, MOI.RawParameter(String(name)), value)
end
"""
    set_optimizer_attribute(stochasticprogram::StochasticProgram, attr::MOI.AbstractOptimizerAttribute, value)

Set the solver-specific attribute `attr` in `stochasticprogram` to `value`.

See also: [`get_optimizer_attribute`](@ref)
"""
function JuMP.set_optimizer_attribute(stochasticprogram::StochasticProgram, attr::MOI.AbstractOptimizerAttribute, value)
    return MOI.set(stochasticprogram, attr, value)
end
"""
    set_optimizer_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)

Given a list of `attribute => value` pairs or a collection of keyword arguments, calls
`set_optimizer_attribute(stochasticprogram, attribute, value)` for each pair.

See also: [`get_optimizer_attribute`](@ref)
"""
function JuMP.set_optimizer_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)
    for (name, value) in pairs
        set_optimizer_attribute(stochasticprogram, name, value)
    end
end
function JuMP.set_optimizer_attributes(stochasticprogram::StochasticProgram; kw...)
    for (name, value) in kw
        set_optimizer_attribute(stochasticprogram, name, value)
    end
end
function JuMP.set_silent(stochasticprogram::StochasticProgram)
    return MOI.set(stochasticprogram, MOI.Silent(), true)
end
function JuMP.unset_silent(stochasticprogram::StochasticProgram)
    return MOI.set(stochasticprogram, MOI.Silent(), false)
end
function JuMP.set_time_limit_sec(stochasticprogram::StochasticProgram, limit)
    return MOI.set(stochasticprogram, MOI.TimeLimitSec(), limit)
end
function JuMP.unset_time_limit_sec(stochasticprogram::StochasticProgram)
    return MOI.set(stochasticprogram, MOI.TimeLimitSec(), nothing)
end
function JuMP.time_limit_sec(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.TimeLimitSec())
end
function JuMP.simplex_iterations(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.SimplexIterations())
end
function JuMP.barrier_iterations(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.BarrierIterations())
end
function JuMP.node_count(stochasticprogram::StochasticProgram)
    return MOI.get(stochasticprogram, MOI.NodeCount())
end
function JuMP.set_optimizer_attribute(stochasticmodel::StochasticModel, name::Union{Symbol, String}, value)
    return set_optimizer_attribute(stochasticmodel, MOI.RawParameter(String(name)), value)
end
function JuMP.set_optimizer_attribute(stochasticmodel::StochasticModel, attr::MOI.AbstractOptimizerAttribute, value)
    return MOI.set(stochasticmodel, attr, value)
end
function JuMP.set_optimizer_attributes(stochasticmodel::StochasticModel, pairs::Pair...)
    for (name, value) in pairs
        set_optimizer_attribute(stochasticmodel, name, value)
    end
end
function JuMP.set_optimizer_attributes(stochasticmodel::StochasticModel; kw...)
    for (name, value) in kw
        set_optimizer_attribute(stochasticmodel, name, value)
    end
end
function JuMP.set_silent(stochasticmodel::StochasticModel)
    return MOI.set(stochasticmodel, MOI.Silent(), true)
end
function JuMP.unset_silent(stochasticmodel::StochasticModel)
    return MOI.set(stochasticmodel, MOI.Silent(), false)
end
"""
    add_scenario!(stochasticprogram::StochasticProgram, stage::Integer, scenario::AbstractScenario)

Store the second stage `scenario` in the `stochasticprogram` at `stage`.

If the `stochasticprogram` is distributed, the scenario will be defined on the node that currently has the fewest scenarios.
"""
function add_scenario!(stochasticprogram::StochasticProgram, stage::Integer, scenario::AbstractScenario)
    add_scenario!(structure(stochasticprogram), stage, scenario)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenario!(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenario)

Store the second stage `scenario` in the two-stage `stochasticprogram`.
"""
function add_scenario!(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenario)
    add_scenario!(structure(stochasticprogram), 2,scenario)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenario!(stochasticprogram::StochasticProgram, stage::Integer, scenario::AbstractScenario, w::Integer)

Store the second stage `scenario` in worker node `w` of the `stochasticprogram` at `stage`.
"""
function add_worker_scenario!(stochasticprogram::StochasticProgram, stage::Integer, scenario::AbstractScenario, w::Integer)
    add_scenario!(structure(stochasticprogram), stage, scenario, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenario!(stochasticprogram::StochasticProgram, stage::Integer, scenario::AbstractScenario, w::Integer)

Store the second stage `scenario` in worker node `w` of the two-stage `stochasticprogram`.
"""
function add_worker_scenario!(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenario, w::Integer)
    add_scenario!(structure(stochasticprogram), 2, scenario, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer)

Store the scenario returned by `scenariogenerator` in the `stage` of the `stochasticprogram`. If the `stochasticprogram` is distributed, the scenario will be defined on the node that currently has the fewest scenarios.
"""
function add_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer)
    add_scenario!(scenariogenerator, structure(stochasticprogram), stage)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram)

Store the scenario returned by `scenariogenerator` in the two-stage `stochasticprogram`. If the `stochasticprogram` is distributed, the scenario will be defined on the node that currently has the fewest scenarios.
"""
function add_scenario!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram)
    add_scenario!(scenariogenerator, structure(stochasticprogram), 2)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer, w::Integer)

Store the scenario returned by `scenariogenerator` in worker node `w` of the `stochasticprogram` at `stage`.
"""
function add_worker_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer, w::Integer)
    add_scenario!(scenariogenerator, structure(stochasticprogram), stage, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenario!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram, w::Integer)

Store the scenario returned by `scenariogenerator` in worker node `w` of the two-stage `stochasticprogram`.
"""
function add_worker_scenario!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram, w::Integer)
    add_scenario!(scenariogenerator, structure(stochasticprogram), 2, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenarios!(stochasticprogram::StochasticProgram, stage::Integer, scenarios::Vector{<:AbstractScenario})

Store the collection of `scenarios` in the `stochasticprogram` at `stage`. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function add_scenarios!(stochasticprogram::StochasticProgram, stage::Integer, scenarios::Vector{<:AbstractScenario})
    add_scenarios!(structure(stochasticprogram), stage, scenarios)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenarios!(stochasticprogram::TwoStageStochasticProgram, scenarios::Vector{<:AbstractScenario})

Store the collection of `scenarios` in the two-stage `stochasticprogram`. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function add_scenarios!(stochasticprogram::TwoStageStochasticProgram, scenarios::Vector{<:AbstractScenario})
    add_scenarios!(structure(stochasticprogram), 2, scenarios)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenarios!(stochasticprogram::StochasticProgram, stage::Integer, scenarios::Vector{<:AbstractScenario}, w::Integer)

Store the collection of `scenarios` in in worker node `w` of the `stochasticprogram` at `stage`.
"""
function add_worker_scenarios!(stochasticprogram::StochasticProgram, stage::Integer, scenarios::Vector{<:AbstractScenario}, w::Integer)
    add_scenarios!(structure(stochasticprogram), stage, scenarios, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenarios!(stochasticprogram::StochasticProgram, stage::Integer, scenarios::Vector{<:AbstractScenario}, w::Integer)

Store the collection of `scenarios` in in worker node `w` of the two-stage `stochasticprogram`.
"""
function add_worker_scenarios!(stochasticprogram::TwoStageStochasticProgram, scenarios::Vector{<:AbstractScenario}, w::Integer)
    add_scenarios!(structure(stochasticprogram), 2, scenarios, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer, n::Integer)

Generate `n` scenarios using `scenariogenerator` and store in the `stochasticprogram` at `stage`. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function add_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer, n::Integer)
    add_scenarios!(scenariogenerator, structure(stochasticprogram), stage, n)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_scenarios!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram, n::Integer)

Generate `n` scenarios using `scenariogenerator` and store in the two-stage `stochasticprogram`. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function add_scenarios!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram, n::Integer)
    add_scenarios!(scenariogenerator, structure(stochasticprogram), 2, n)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer, n::Integer, w::Integer)

Generate `n` scenarios using `scenariogenerator` and store them in worker node `w` of the `stochasticprogram` at `stage`.
"""
function add_worker_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer, n::Integer, w::Integer)
    add_scenarios!(scenariogenerator, structure(stochasticprogram), stage, n, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    add_worker_scenarios!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram, n::Integer, w::Integer)

Generate `n` scenarios using `scenariogenerator` and store them in worker node `w` of the two-stage `stochasticprogram`.
"""
function add_worker_scenarios!(scenariogenerator::Function, stochasticprogram::TwoStageStochasticProgram, n::Integer, w::Integer)
    add_scenarios!(scenariogenerator, structure(stochasticprogram), 2, n, w)
    invalidate_cache!(stochasticprogram)
    return stochasticprogram
end
"""
    sample!(stochasticprogram::StochasticProgram, stage::Integer, sampler::AbstractSampler, n::Integer)

Sample `n` scenarios using `sampler` and add to the `stochasticprogram` at `stage`. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function sample!(stochasticprogram::StochasticProgram, stage::Integer, sampler::AbstractSampler, n::Integer)
    sample!(structure(stochasticprogram), stage, sampler, n)
    return stochasticprogram
end
"""
    sample!(stochasticprogram::TwoStageStochasticProgram, sampler::AbstractSampler, n::Integer)

Sample `n` scenarios using `sampler` and add to the two-stage `stochasticprogram`. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function sample!(stochasticprogram::TwoStageStochasticProgram, sampler::AbstractSampler, n::Integer)
    sample!(structure(stochasticprogram), 2, sampler, n)
    return stochasticprogram
end
# ========================== #


text:
module MCMC_LR_Tests

using ArgCheck
using Distributions
using MCMCDiagnostics

export mean_LR_pvalue, cov_LR_pvalue

"""
    column_ρs(X)

For observations in the rows of `X`, return the estimated `ρ=ESS/n` for each column (variable).

Provides the default value to [`mean_LR_statistic`](@ref), otherwise only useful for debugging.
"""
column_ρs(X::AbstractMatrix) = vec(mapslices(first ∘ ess_factor_estimate, X, 1))

"""
    mean_LR_pvalue(μ, X, [ρ])

P-value for a likelihood ratio statistic for testing that the columnwise mean of `X` are equal to the vector `μ`, if the rows of `X` are normally distributed (with unknown variance, estimated from the sample).

Assume that some form of a CLT holds, so

``n(x̄(n)-x̄) → N(0, Σ)`` in distribution, where `x̄(n)` is the running mean.

The test statistic is

``(x̄(n)-μ)' Σ⁻¹ (x̄(n)-μ) ∼ T²(k, n-1)``

where `T²` is [Hotelling's T²-distribution](https://en.wikipedia.org/wiki/Hotelling%27s_T-squared_distribution), if `x̄ = μ` (Dootika et al 2017).

For IID samples, `Σ` can be estimated from the sample. However, MCMC draws are autocorrelated, so this is more difficult. Let ρᵢ = ESSᵢ/n denote the scaling factor to obtain the effective sample size. For a crude approximation, we rescale the variance by ρᵢ, and covariances by √(ρᵢρⱼ).
"""
function mean_LR_pvalue(μ::AbstractVector, X::AbstractMatrix, ρ = column_ρs(X))
    N, K = size(X)
    @argcheck K == length(μ) == length(ρ)
    R = Diagonal(.√ρ)           # correction by ρ = ESS/n
    S = cov(X, 1, false)        # sample covariance
    δ = vec(mean(X, 1)) - μ     # deviations from μ
    h = (δ' * (R*inv(S)*R) * δ) * N    # has Hotelling T²(K, N-1) distribution
    ccdf(FDist(K, N-K), h*((N-K)/((N-1)*K))) # convert to F(K, N-K) distribution
end

"""
    cov_LR_pvalue(Σ, X)

P-value for likelihood ratio statistic for comparing the covariance of `X` (with observations in rows) to `Σ`.

Note that this is an *asymptotic* statistic, so it has better coverage when `N ≫ K`. Coverage may be very bad when `N` is small, see Bai et al (2009).
"""
function cov_LR_pvalue(Σ::AbstractMatrix, X::AbstractMatrix)
    N, K = size(X)
    @argcheck size(Σ) == (K, K)
    U = chol(Symmetric(Σ))
    S = cov(X / U, 1)
    L = trace(S) - logdet(S) - K
    ccdf(Chisq(K*(K+1)/2), N*L)
end

end # module