text:
using Test, LinearAlgebra

@testset "test/unit/LocalPolytope.jl" begin

using BellScenario: LocalPolytope

@testset "strategies()" begin

    @testset "trivial strategies" begin
        @test LocalPolytope.strategies(1,1) == [fill(1,(1,1))]
        @test LocalPolytope.strategies(1,2) == [[1 0]',[0 1]']
        @test LocalPolytope.strategies(2,1) == [[1 1]]
        @test LocalPolytope.strategies(3,1) == [[1 1 1]]
    end

    @testset "dichotomic strategies" begin
        s = LocalPolytope.strategies(2,2)
        @test size(s) == (4,)

        @test s[1] == [1 1; 0 0]
        @test s[2] == [0 1; 1 0]
        @test s[3] == [1 0; 0 1]
        @test s[4] == [0 0; 1 1]
    end

    @testset "two_input_strategies are unique amongst themselves" begin
        strategies = LocalPolytope.strategies(4,2)
        @test size(strategies) == (16,)
        for test_strategy in strategies
            duplicate_count = 0
            for match_strategy in strategies
                if test_strategy == match_strategy
                    duplicate_count += 1
                end
            end
            # test_strategy should only match with itself
            @test duplicate_count == 1
        end
    end

    @testset "checking size of set" begin
        s = LocalPolytope.strategies(10,3)
        @test size(s) == (59049,)
        @test size(s[1]) == (3,10)

        s = LocalPolytope.strategies(3,10)
        @test size(s) == (1000,)
        @test size(s[1]) == (10,3)

        s = LocalPolytope.strategies(5,5)
        @test size(s) == (3125,)
        @test size(s[1]) == (5,5)
    end
end

@testset "input_set()" begin
    LocalPolytope.input_set(1)
    @testset "num_inputs = $N" for N in 1:10
        input_set = LocalPolytope.input_set(N)
        test_matrix = zeros(N,N)
        for i in 1:N
            test_matrix[:,i] += input_set[i]
        end
        identity_matrix = diagm(0 => fill(1,N))

        # the input set should form the identity matrix
        @test test_matrix == identity_matrix
    end
end

@testset "communication_protocols()" begin

    @testset "fewer inputs than outputs" begin
        @testset "simple cases" begin
            @test LocalPolytope.communication_protocols(1,2) == [[1 0]']
            @test LocalPolytope.communication_protocols(2,3) == [[1 0;0 1;0 0]]
            @test LocalPolytope.communication_protocols(2,4) == [[1 0;0 1;0 0;0 0]]
            @test LocalPolytope.communication_protocols(3,4) == [[1 0 0;0 1 0;0 0 1;0 0 0]]
        end

        @testset "scales like binomial coefficient: outputs = $num_out " for num_out in 2:7
            @testset "inputs = $num_in" for num_in in 1:num_out
                @test length(LocalPolytope.communication_protocols(num_in,num_out)) == 1
            end
        end
    end

    @testset "trivial edge cases" begin
        @test LocalPolytope.communication_protocols(2,2) == [[1 0;0 1]]
        @test LocalPolytope.communication_protocols(3,2) == [[1 1 0;0 0 1],[1 0 1;0 1 0],[1 0 0;0 1 1]]
    end

    @testset "there exists a single protocol for $N:$N and $N:1" for N in 1:10
        @test 1 == size(LocalPolytope.communication_protocols(N,N))[1]
        @test 1 == size(LocalPolytope.communication_protocols(N,1))[1]
    end

    stirling_base_cases = [
        (3,2,3),
        (4,2,7),
        (5,2,15),
        (6,2,31),
        (7,2,63),
        (8,2,127),
        (4,3,6),
        (5,3,25),
        (6,3,90),
        (5,4,10),
        (6,4,65),
        (7,4,350),
        (6,5,15),
        (7,6,21),
        (8,7,28),
    ]

    @testset "verify stirling number base case: n=$(case[1]), k=$(case[2])" for case in stirling_base_cases
        protocols = LocalPolytope.communication_protocols(case[1],case[2])
        # protocols match stirling's number of the second kind
        @test case[3] == size(protocols)[1]
        # protocols are unique
        @test size(protocols) == size(unique(protocols))

        # protocols map a single input to a single output
        for protocol in protocols
            row_sum = zeros(case[1])
            for row_id in 1:case[2]
                row_sum += protocol[row_id,:]
            end
            @test row_sum == ones(case[1])
        end
    end
end

@testset "is_no_signaling_duplicate()" begin
    @test LocalPolytope.is_no_signaling_duplicate([1 1 1 1;0 0 0 0], 2, 2)
    @test LocalPolytope.is_no_signaling_duplicate([1 1 0 0;0 0 1 1], 2, 2)
    @test !( LocalPolytope.is_no_signaling_duplicate([1 0 1 1;0 1 0 0], 2, 2) )
    @test LocalPolytope.is_no_signaling_duplicate([1 1 1 0 0 0;0 0 0 1 1 1], 2, 3)
end

@testset "vertices()" begin
    @testset "trivial no-signaling cases" begin
        @test LocalPolytope.vertices((1,1),(1,1)) == [fill(1,(1,1))]

        @test LocalPolytope.vertices((2,1),(1,1)) == [[1 1]']
        @test LocalPolytope.vertices((1,2),(1,1)) == [fill(1,(1,1)),fill(0,(1,1))]
        @test LocalPolytope.vertices((1,1),(2,1)) == [[1 1]']
        @test LocalPolytope.vertices((1,1),(1,2)) == [fill(1,(1,1)),fill(0,(1,1))]

        @test LocalPolytope.vertices((2,2),(1,1)) == [[1 1]',[0 1]',[1 0]',[0 0]']
        @test LocalPolytope.vertices((2,1),(2,1)) == [[1 1 1 1]']
        @test LocalPolytope.vertices((2,1),(1,2)) == [[1 1]',[0 0]']
        @test LocalPolytope.vertices((1,2),(2,1)) == [[1 1]',[0 0]']
        @test LocalPolytope.vertices((1,2),(1,2)) == [
            [1 1 1]',[1 0 0]',[0 1 0]',[0 0 0]'
        ]
        @test LocalPolytope.vertices((1,1),(2,2)) == [[1 1]',[0 1]',[1 0]',[0 0]']

        @test LocalPolytope.vertices((2,2),(2,1)) == [
            [1 1 1 1]',[0 0 1 1]',[1 1 0 0]',[0 0 0 0]'
        ]
        @test LocalPolytope.vertices((2,2),(1,2)) == [
            [1 1 1 1 1]',[1 1 0 0 0]',[0 1 1 0 1]',[0 1 0 0 0]',
            [1 0 1 1 0]',[1 0 0 0 0]',[0 0 1 0 0]',[0 0 0 0 0]'
        ]
        @test LocalPolytope.vertices((2,1),(2,2)) == [
            [1 1 1 1]',[0 1 0 1]',[1 0 1 0]',[0 0 0 0]'
        ]
        @test LocalPolytope.vertices((1,2),(2,2)) == [
            [1 1 1 1 1]',[1 0 1 0 1]',[1 1 0 1 0]',[1 0 0 0 0]',
            [0 1 1 0 0]',[0 0 1 0 0]',[0 1 0 0 0]',[0 0 0 0 0]'
        ]

        @test LocalPolytope.vertices((3,1),(1,2)) == [[1 1 1]',[0 0 0]']
        @test LocalPolytope.vertices((3,1),(1,3)) == [[1 1 1 0 0 0]',[0 0 0 1 1 1]',[0 0 0 0 0 0]']
    end

    @testset "trivial fixed-direction signaling cases" begin
        @test LocalPolytope.vertices((1,1),(1,1),bits=1) == [fill(1,(1,1))]

        @test LocalPolytope.vertices((2,1),(1,1),bits=1) == [[1 1]']
        @test LocalPolytope.vertices((1,2),(1,1),bits=1) == [fill(1,(1,1)),fill(0,(1,1))]
        @test LocalPolytope.vertices((1,1),(2,1),bits=1) == [[1 1]']
        @test LocalPolytope.vertices((1,1),(1,2),bits=1) == [fill(1,(1,1)),fill(0,(1,1))]

        @test LocalPolytope.vertices((2,2),(1,1),bits=1) == [[1 1]',[0 1]',[1 0]',[0 0]']
        @test LocalPolytope.vertices((2,1),(2,1),bits=1) == [[1 1 1 1]']
        @test LocalPolytope.vertices((2,1),(1,2),bits=1) == [[1 1]',[0 0]',[0 1]',[1 0]']
        @test LocalPolytope.vertices((1,1),(2,2),bits=1) == [[1 1]',[0 1]',[1 0]',[0 0]']

        @test LocalPolytope.vertices((2,2),(2,1),bits=1) == [
            [1 1 1 1]',[0 0 1 1]',[1 1 0 0]',[0 0 0 0]'
        ]
        @test issetequal(LocalPolytope.vertices((2,2),(1,2),bits=1), [
            [0 0 0 0 0 0]',[0 0 0 0 1 1]',[0 0 0 0 1 0]',[0 0 0 0 0 1]',
            [1 1 0 0 0 0]',[1 1 1 1 0 0]',[1 1 1 0 0 0]',[1 1 0 1 0 0]',
            [0 1 0 0 0 0]',[0 1 0 1 1 0]',[0 1 0 1 0 0]',[0 1 0 0 1 0]',
            [1 0 0 0 0 0]',[1 0 1 0 0 1]',[1 0 1 0 0 0]',[1 0 0 0 0 1]'
        ])
        @test issetequal(LocalPolytope.vertices((2,1),(2,2),bits=1), [
                [0 0 0 0]',[0 0 0 1]',[0 0 1 0]',[0 0 1 1]',
                [0 1 0 0]',[0 1 0 1]',[0 1 1 0]',[0 1 1 1]',
                [1 0 0 0]',[1 0 0 1]',[1 0 1 0]',[1 0 1 1]',
                [1 1 0 0]',[1 1 0 1]',[1 1 1 0]',[1 1 1 1]'
        ])

        @test issetequal(LocalPolytope.vertices((3,1),(1,2),bits=1), [
            [1 1 1]',[0 0 0]',[1 0 0]',[1 1 0]',[0 1 0]',[0 1 1]',[0 0 1]',[1 0 1]'
        ])
        # there are no strategies which allow for 3 distinct outputs simultaneously
        @test issetequal(LocalPolytope.vertices((3,1),(1,3),bits=1), [
            [1 0 0 0 0 0]',[1 1 0 0 0 0]',[1 0 1 0 0 0]',[1 1 1 0 0 0]',[1 1 0 0 0 1]',[1 0 1 0 1 0]',[1 0 0 0 1 1]',
            [0 0 0 1 0 0]',[0 0 0 1 1 0]',[0 0 0 1 0 1]',[0 1 1 1 0 0]',[0 1 0 1 0 1]',[0 0 1 1 1 0]',[0 0 0 1 1 1]',
            [0 0 0 0 0 0]',[0 1 0 0 0 0]',[0 0 1 0 0 0]',[0 1 1 0 0 0]',[0 0 0 0 1 1]',[0 0 0 0 1 0]',[0 0 0 0 0 1]'
        ])
    end

    @testset "trivial one-way signaling cases" begin
        @test LocalPolytope.vertices((1,1),(1,1),bits=1,fixed_direction=false) == [fill(1,(1,1))]

        @test LocalPolytope.vertices((2,1),(1,1),bits=1,fixed_direction=false) == [[1 1]']
        @test LocalPolytope.vertices((1,2),(1,1),bits=1,fixed_direction=false) == [fill(1,(1,1)),fill(0,(1,1))]
        @test LocalPolytope.vertices((1,1),(2,1),bits=1,fixed_direction=false) == [[1 1]']
        @test LocalPolytope.vertices((1,1),(1,2),bits=1,fixed_direction=false) == [fill(1,(1,1)),fill(0,(1,1))]

        @test LocalPolytope.vertices((2,2),(1,1),bits=1,fixed_direction=false) == [[1 1]',[0 1]',[1 0]',[0 0]']
        @test LocalPolytope.vertices((2,1),(2,1),bits=1,fixed_direction=false) == [[1 1 1 1]']
        @test LocalPolytope.vertices((2,1),(1,2),bits=1,fixed_direction=false) == [[1 1]',[0 0]',[0 1]',[1 0]']
        @test LocalPolytope.vertices((1,2),(2,1),bits=1,fixed_direction=false) == [[1 1]',[0 0]',[0 1]',[1 0]']
        @test LocalPolytope.vertices((1,2),(1,2),bits=1,fixed_direction=false) == [[1 0 0]',[0 1 0]',[0 0 1]',[0 0 0]']
        @test LocalPolytope.vertices((1,1),(2,2),bits=1,fixed_direction=false) == [[1 1]',[0 1]',[1 0]',[0 0]']

        @test issetequal(LocalPolytope.vertices((2,2),(2,1),bits=1,fixed_direction=false), [
            [0 0 0 0]',[0 0 0 1]',[0 0 1 0]',[0 0 1 1]',
            [0 1 0 0]',[0 1 0 1]',[0 1 1 0]',[0 1 1 1]',
            [1 0 0 0]',[1 0 0 1]',[1 0 1 0]',[1 0 1 1]',
            [1 1 0 0]',[1 1 0 1]',[1 1 1 0]',[1 1 1 1]'
        ])
        @test issetequal(LocalPolytope.vertices((2,2),(1,2),bits=1,fixed_direction=false), [
            [1 1 0 0 0 0]',[1 0 0 1 0 0]',[1 0 0 0 0 1]',[1 0 0 0 0 0]',
            [0 1 1 0 0 0]',[0 0 1 1 0 0]',[0 0 1 0 0 1]',[0 0 1 0 0 0]',
            [0 1 0 0 1 0]',[0 0 0 1 1 0]',[0 0 0 0 1 1]',[0 0 0 0 1 0]',
            [0 1 0 0 0 0]',[0 0 0 1 0 0]',[0 0 0 0 0 1]',[0 0 0 0 0 0]'
        ])
        @test issetequal(LocalPolytope.vertices((2,1),(2,2),bits=1,fixed_direction=false), [
            [0 0 0 0]',[0 0 0 1]',[0 0 1 0]',[0 0 1 1]',
            [0 1 0 0]',[0 1 0 1]',[0 1 1 0]',[0 1 1 1]',
            [1 0 0 0]',[1 0 0 1]',[1 0 1 0]',[1 0 1 1]',
            [1 1 0 0]',[1 1 0 1]',[1 1 1 0]',[1 1 1 1]'
        ])
        @test issetequal(LocalPolytope.vertices((1,2),(2,2),bits=1,fixed_direction=false), [
            [1 1 0 0 0 0]',[1 0 0 1 0 0]',[1 0 0 0 0 1]',[1 0 0 0 0 0]',
            [0 1 1 0 0 0]',[0 0 1 1 0 0]',[0 0 1 0 0 1]',[0 0 1 0 0 0]',
            [0 1 0 0 1 0]',[0 0 0 1 1 0]',[0 0 0 0 1 1]',[0 0 0 0 1 0]',
            [0 1 0 0 0 0]',[0 0 0 1 0 0]',[0 0 0 0 0 1]',[0 0 0 0 0 0]'
        ])
    end

    @testset "verify (22)(22) no-signaling setup" begin
        vertices = LocalPolytope.vertices((2,2), (2,2))
        @test size(vertices) == (16,)

        @testset "(22)(22) vertices are well-formed and unique" begin
            for vertex in vertices
                @test size(vertex) == (8,1)
            end

            @test size(vertices) == size(unique(vertices))
        end

        # check results against specialized (22)(22) function
        vertices_match_p8 = map(
            x -> begin
                catv = cat([1],x,dims=1)
                v8 = (Behavior.gen_to_ns_proj((2,2),(2,2))*catv)[2:end]
                reshape(v8,(8,1))
            end,
            DichotomicLocalPolytope.deterministic_no_signaling_behaviors()
        )

        @test issetequal(vertices, vertices_match_p8)
    end

    @testset "verify (22)-1->(22) fixed direction signaling case" begin
        alice_expt = (2,2)
        bob_expt = (2,2)

        vertices = LocalPolytope.vertices(alice_expt, bob_expt, bits=1)
        @test size(vertices) == (64,)
        @test size(vertices) == size(unique(vertices))

        # check against existing results
        vertices_p10 = DichotomicLocalPolytope.deterministic_fixed_direction_signaling_behaviors(subspace="p10")

        # need to recast vertices as matrix and not array
        adjoint_match = []
        for vertex in vertices_p10
            push!(adjoint_match, reshape(vertex, length(vertex),1))
        end
        adjoint_match

        @test issetequal(vertices, adjoint_match)
    end

    @testset "verify (22)<-1->(22) one-way signaling case" begin
        vertices = LocalPolytope.vertices((2,2), (2,2), bits=1, fixed_direction=false)

        @test size(vertices) == (112,)
        @test size(vertices) == size(unique(vertices))

        match_vertices = map(
            (x) -> reshape(x,(12,1)),
            DichotomicLocalPolytope.deterministic_one_way_signaling_behaviors(subspace="p12")
        )

        @test issetequal(vertices, match_vertices)
    end

    test_cases = [
        ((2,2),(2,2),1),
        ((3,2),(3,2),1),
    ]
    @testset "equivalence of dits and bits" for case in test_cases

        vbits = LocalPolytope.vertices(case[1], case[2], bits=case[3])
        vdits = LocalPolytope.vertices(case[1], case[2], dits=2^case[3])

        @test vbits == vdits
    end

    @testset "one dit is the same as no-signaling" for case in test_cases
        vdits = LocalPolytope.vertices(case[1],case[2],dits=1)
        vns = LocalPolytope.vertices(case[1],case[2],rep="fixed-direction")

        @test vdits == vns
    end

    @testset "odd numbered dits" begin
        vdits = LocalPolytope.vertices((3,1),(1,3),dits=3)
        v_match = [
            [0 0 0 0 0 0]',[0 0 1 0 0 0]',[0 1 0 0 0 0]',[0 1 1 0 0 0]',[1 0 0 0 0 0]',
            [1 0 1 0 0 0]',[1 1 0 0 0 0]',[1 1 1 0 0 0]',[0 0 0 0 0 1]',[0 0 0 0 1 0]',
            [0 0 0 0 1 1]',[0 0 0 1 0 0]',[0 0 0 1 0 1]',[0 0 0 1 1 0]',[0 0 0 1 1 1]',
            [0 0 1 0 1 0]',[0 0 1 1 0 0]',[0 0 1 1 1 0]',[0 1 0 0 0 1]',[0 1 0 1 0 0]',
            [0 1 0 1 0 1]',[1 0 0 0 0 1]',[1 0 0 0 1 0]',[1 0 0 0 1 1]',[0 1 1 1 0 0]',
            [1 0 1 0 1 0]',[1 1 0 0 0 1]'
        ]

        @test issetequal(vdits, v_match)
    end
end

@testset "compute_no_signaling_behavior()" begin
    @testset "trivial cases" begin
        @test LocalPolytope.compute_no_signaling_behavior(
            ([1], [fill(1,(1,1))], 1, 1),
            ([1], [fill(1,(1,1))], 1, 1),
            subspace="generalized"
        ) == fill(1,(1,1))
        @test LocalPolytope.compute_no_signaling_behavior(
            ([1 0], [[1 0]',[0 1]'], 2, 1),
            ([1], [fill(1,(1,1))], 1, 1),
            subspace="generalized"
        ) == [1 0]'
    end
end

@testset "strategy_to_behavior()" begin
    @test LocalPolytope.strategy_to_behavior([1 2 3;4 5 6], constant = false) == [1 2 3 4 5 6]'
    @test LocalPolytope.strategy_to_behavior([1 2;3 4;5 6], constant = false) == [1 2 3 4 5 6]'
    @test LocalPolytope.strategy_to_behavior([1 2 3;4 5 6;7 8 9], constant = false) == [1 2 3 4 5 6 7 8 9]'

    @test LocalPolytope.strategy_to_behavior([1 2 3;4 5 6]) == [1 1 2 3 4 5 6]'
    @test LocalPolytope.strategy_to_behavior([1 2;3 4;5 6]) == [1 1 2 3 4 5 6]'
    @test LocalPolytope.strategy_to_behavior([1 2 3;4 5 6;7 8 9]) == [1 1 2 3 4 5 6 7 8 9]'
end

@testset "behavior_to_strategy()" begin
    @test LocalPolytope.behavior_to_strategy(3, 3, [1 1 0 0 0 1 0 0 0 1]') == [1 0 0;0 1 0;0 0 1]
    @test LocalPolytope.behavior_to_strategy(3, 3, [1 0 0 0 1 0 0 0 1]') == [1 0 0;0 1 0;0 0 1]
    @test LocalPolytope.behavior_to_strategy(3, 3, [1 0 0 0 1 1 0 0 0 1]') == [0 0 0;1 1 0;0 0 1]
    @test LocalPolytope.behavior_to_strategy(4, 2, [1 0 1 0 0 1 0 1]') == [1 0 1 0;0 1 0 1]

    @test_throws DomainError LocalPolytope.behavior_to_strategy(3, 3, [1 1 0 0 0 1 0]')
end

@testset "num_prepare_and_measure_vertices()" begin
    for N in 1:5
        for d in 1:5
            α = (N,1)
            β = (1,N)
            vertices = unique(LocalPolytope.vertices(α, β, dits=d))
            @test length(vertices) == LocalPolytope.num_prepare_and_measure_vertices(N,d)
        end
    end
end

@testset "num_inhomogeneous_prepare_and_measure_vertices()" begin
    for X in 1:5
        for B in 1:5
            for d in 1:5
                α = (X,1)
                β = (1,B)
                vertices = unique(LocalPolytope.vertices(α, β, dits=d))
                @test length(vertices) == LocalPolytope.num_inhomogeneous_prepare_and_measure_vertices(X,B,d)
            end
        end
    end
end

@testset "facet_to_matrix()" begin
    @test LocalPolytope.facet_to_matrix(3,3, [-2 1 0 0 0 1 0 0 0 1]) == (2, [1 0 0;0 1 0;0 0 1])
    @test LocalPolytope.facet_to_matrix(3,3, [-1 1 0 -1 0 1 -1 0 0 0]) == (2, [1 0 0;0 1 0;0 0 1])

    @test LocalPolytope.facet_to_matrix(5,5,
        [-1 0 0 -1 1 -1 0 0 0 0 0 0 0 1 -1 -1 1 0 -1 -1 0 0 0 0 0 0]
    ) == (4, [0 0 0 2 0;0 0 1 1 1;0 0 2 0 0;1 0 0 0 1;0 0 1 1 1])
end

end


text:
using Documenter, Juno, AsyPlots

makedocs(
    modules = [AsyPlots],
    format = Documenter.HTML(),
    sitename = "AsyPlots.jl",
    doctest = false,
    strict = false,
    checkdocs = :none,
    pages = Any[
        "Home" => "index.md",
        "Manual" => [
            "man/overview.md",
            "man/asy2d.md",
            "man/asy3d.md",
            "man/pens.md",
	    "man/heatmaps.md",
	    "man/isolines.md",
            "man/graphing.md",
            "man/output.md",
            "man/settings.md",
            "man/defaults.md",
            "man/animate.md"
            ],
        "Examples" => [
            "examples/examples2d.md",
            "examples/examples3d.md"
        ]
    ]
)

deploydocs(
    repo="github.com/sswatson/AsyPlots.jl.git",
    target="build",
    deps=nothing,
    make=nothing
)

ext: jl


text:
abstract type AbstractComplexTerm{T} <: Symbolic{Complex{T}}
end

struct ComplexTerm{T} <: AbstractComplexTerm{T}
    re
    im
end

Base.imag(c::Symbolic{Complex{T}}) where {T} = term(imag, c)
SymbolicUtils.promote_type(::typeof(imag), ::Type{Complex{T}}) where {T} = T
Base.promote_rule(::Type{Complex{T}}, ::Type{S}) where {T<:Real, S<:Num} =  Complex{S} # 283

has_symwrapper(::Type{<:Complex{T}}) where {T<:Real} = true
wraps_type(::Type{Complex{Num}}) = Complex{Real}
iswrapped(::Complex{Num}) = true
function wrapper_type(::Type{Complex{T}}) where T
    Symbolics.has_symwrapper(T) ? Complex{wrapper_type(T)} : Complex{T}
end

TermInterface.symtype(a::ComplexTerm{T}) where T = Complex{T}
TermInterface.istree(a::ComplexTerm) = true
TermInterface.operation(a::ComplexTerm{T}) where T = Complex{T}
TermInterface.arguments(a::ComplexTerm) = [a.re, a.im]

function TermInterface.similarterm(t::ComplexTerm, f, args, symtype; metadata=nothing, exprhead=exprhead(t))
    if f <: Complex
        ComplexTerm{real(f)}(args...)
    else
        similarterm(first(args), f, args, symtype; metadata=metadata, exprhead=exprhead)
    end
end

function Base.show(io::IO, mime::MIME"text/plain", a::ComplexTerm)
    print(io, "ComplexTerm(")
    show(io, mime, wrap(a))
    print(io, ")")
end

function Base.show(io::IO, a::Complex{Num})
    rr = unwrap(real(a))
    ii = unwrap(imag(a))

    if istree(rr) && (operation(rr) === real) &&
        istree(ii) && (operation(ii) === imag) &&
        isequal(arguments(rr)[1], arguments(ii)[1])

        return print(io, arguments(rr)[1])
    end

    i = Sym{Real}(:im)
    show(io, real(a) + i * imag(a))
end

function unwrap(a::Complex{<:Num})
    re, img = unwrap(real(a)), unwrap(imag(a))
    if re isa Real && img isa Real
        return re + im * img
    end
    T = promote_type(symtype(re), symtype(img))
    ComplexTerm{T}(re, img)
end
wrap(a::ComplexTerm) = Complex(wrap.(arguments(a))...)
wrap(a::Symbolic{<:Complex}) = Complex(wrap(real(a)), wrap(imag(a)))

SymbolicUtils.@number_methods(
                              ComplexTerm,
                              unwrap(f(wrap(a))),
                              unwrap(f(wrap(a), wrap(b))),
                             )

function Base.isequal(a::ComplexTerm{T}, b::ComplexTerm{S}) where {T,S}
    T === S && isequal(a.re, b.re) && isequal(a.im, b.im)
end

function Base.hash(a::ComplexTerm{T}, h::UInt) where T
    hash(hash(a.im, hash(a.re, hash(T, hash(h ⊻ 0x1af5d7582250ac4a)))))
end

Base.iszero(x::Complex{<:Num}) = iszero(real(x)) && iszero(imag(x))
Base.isone(x::Complex{<:Num}) = isone(real(x)) && iszero(imag(x))
_iszero(x::Complex{<:Num}) = _iszero(unwrap(x))
_isone(x::Complex{<:Num}) = _isone(unwrap(x))


text:
module SessionActions

import ..Pluto: ServerSession, Notebook, Cell, emptynotebook, tamepath, new_notebooks_directory, without_pluto_file_extension, numbered_until_new, readwrite, update_save_run!, update_from_file, wait_until_file_unchanged, putnotebookupdates!, putplutoupdates!, load_notebook, clientupdate_notebook_list, WorkspaceManager, @asynclog
using FileWatching

struct NotebookIsRunningException <: Exception
    notebook::Notebook
end

abstract type AbstractUserError <: Exception end
struct UserError <: AbstractUserError
    msg::String
end
function Base.showerror(io::IO, e::UserError)
    print(io, e.msg)
end

function open_url(session::ServerSession, url::AbstractString; kwargs...)
    path = download(url, emptynotebook().path)
    open(session, path; kwargs...)
end

function open(session::ServerSession, path::AbstractString; run_async=true, compiler_options=nothing, as_sample=false)
    if as_sample
        new_filename = "sample " * without_pluto_file_extension(basename(path))
        new_path = numbered_until_new(joinpath(new_notebooks_directory(), new_filename); suffix=".jl")

        readwrite(path, new_path)
        path = new_path
    end

    for nb in values(session.notebooks)
        if realpath(nb.path) == realpath(tamepath(path))
            throw(NotebookIsRunningException(nb))
        end
    end

    nb = load_notebook(tamepath(path); disable_writing_notebook_files=session.options.server.disable_writing_notebook_files)

    # overwrites the notebook environment if specified
    if compiler_options !== nothing
        nb.compiler_options = compiler_options
    end

    session.notebooks[nb.notebook_id] = nb
    if session.options.evaluation.run_notebook_on_load
        for c in nb.cells
            c.queued = true
        end
        update_save_run!(session, nb, nb.cells; run_async=run_async, prerender_text=true)
    end

    add(session, nb; run_async=run_async)

    nb
end

function add(session::ServerSession, nb::Notebook; run_async::Bool=true)
    session.notebooks[nb.notebook_id] = nb

    if run_async
        @asynclog putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    else
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    end

    running = Ref(false)
    function update_from_file_throttled()
        if !running[]
            running[] = true

            @info "Updating from file..."

		    sleep(0.1) ## There seems to be a synchronization issue if your OS is VERYFAST
            wait_until_file_unchanged(nb.path, .3)
            update_from_file(session, nb)

            @info "Updating from file done!"

            running[] = false
        end
    end

    session.options.server.auto_reload_from_file && @asynclog while true
        if !isfile(nb.path)
            # notebook file deleted... let's ignore this, changing the notebook will cause it to save again. Fine for now
            sleep(2)
        else
            watch_file(nb.path)
            # the above call is blocking until the file changes

            current_time = time()
            modified_time = mtime(nb.path)
            @info "File changed" (current_time - nb.last_save_time) (modified_time - nb.last_save_time) (current_time - modified_time)

            # if current_time - nb.last_save_time < 2.0
            if modified_time - nb.last_save_time < session.options.server.auto_reload_from_file_cooldown
                # @info "Notebook was saved by me very recently, not reloading from file."
                @info "Modified time is very close to my last save time, not reloading from file."
            else
                update_from_file_throttled()
            end
        end
    end

    nb
end

function save_upload(content::Vector{UInt8})
    save_path = emptynotebook().path
    Base.open(save_path, "w") do io
        write(io, content)
    end

    save_path
end

function new(session::ServerSession; run_async=true)
    nb = if session.options.server.init_with_file_viewer

        @warn "DEPRECATED: init_with_file_viewer will be removed soon."

        file_viewer_code = """html\"\"\"

        <script>

        const nbfile_url = window.location.href.replace("edit", "notebookfile")

        const pre = html`<pre style="font-size: .6rem;">Loading...</pre>`

        const handle = setInterval(async () => {

        pre.innerText = await (await fetch(nbfile_url)).text()
        }, 500)

        invalidation.then(() => {
        clearInterval(handle)
        })

        return pre

        </script>

        \"\"\"
        """
        Notebook([Cell(), Cell(code=file_viewer_code, code_folded=true)])

    else
        emptynotebook()
    end
    update_save_run!(session, nb, nb.cells; run_async=run_async, prerender_text=true)

    add(session, nb; run_async=run_async)

    nb
end

function shutdown(session::ServerSession, notebook::Notebook; keep_in_session=false, async=false)
    notebook.nbpkg_restart_recommended_msg = nothing
    notebook.nbpkg_restart_required_msg = nothing

    if !keep_in_session
        listeners = putnotebookupdates!(session, notebook) # TODO: shutdown message
        delete!(session.notebooks, notebook.notebook_id)
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
        for client in listeners
            @async close(client.stream)
        end
    end
    WorkspaceManager.unmake_workspace((session, notebook); async=async)
end

end


text:
mutable struct ParamTable
    m :: Matrix{String}
end

"""
    $(SIGNATURES)

Holds symbol, name, description, value.
All as formatted text.
Used for constructing formatted parameter tables.
"""
ParamTable() = ParamTable(Matrix{String}(undef, 0, 4));
ParamTable(n :: Integer) = ParamTable(fill("", n, 4));

Base.length(pt :: ParamTable) = size(pt.m, 1);

get_symbol(pt :: ParamTable, ir :: Integer) = pt.m[ir, 1];
get_description(pt :: ParamTable, ir :: Integer) = pt.m[ir, 2];
get_value(pt :: ParamTable, ir :: Integer) = pt.m[ir, 3];
get_name(pt :: ParamTable, ir :: Integer) = pt.m[ir, 4];

get_symbols(pt :: ParamTable) = pt.m[:, 1];
get_descriptions(pt :: ParamTable) = pt.m[:, 2];
get_values(pt :: ParamTable) = pt.m[:, 3];
get_names(pt :: ParamTable) = pt.m[:, 4];

"""
	$(SIGNATURES)

Add a row to a `ParamTable`.
"""
function add_row!(pt :: ParamTable, name, lsymbol, descr, value)
    append!(pt.m, [lsymbol descr value name]);
end

"""
	$(SIGNATURES)

Replace an existing row in a `ParamTable`.
"""
function set_row!(pt :: ParamTable, ir :: Integer, name, lsymbol, descr, value)
    pt.m[ir,:] .= [lsymbol, descr, value, name];
end

# Return a vector of latex table rows with
#   symbol / description / value
# Contains only the "bodies" of the rows. Not newlines etc.
# First row is multicolumn description, if provided.
function latex_param_table(tbM :: ParamTable, descr)
    n = length(tbM);
    if isempty(descr)
        lOffset = 0;
    else
        lOffset = 1;
    end
    sV = fill("", n + lOffset);
    if !isempty(descr)
        sV[1] = "\\multicolumn{3}{l}{$descr} \\\\";
    end
    for j = 1 : n
        sV[j + lOffset] = "\$$(get_symbol(tbM, j))\$ & $(get_description(tbM, j)) & $(get_value(tbM, j))"
    end
    return sV
end

# ------------


text:
module SpringCollab2020AttachedIceWall

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/AttachedIceWall" AttachedIceWall(x::Integer, y::Integer, height::Integer=8, left::Bool=false)

const placements = Ahorn.PlacementDict(
    "Attached Ice Wall (Right) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        AttachedIceWall,
        "rectangle",
        Dict{String, Any}(
            "left" => true
        )
    ),
    "Attached Ice Wall (Left) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        AttachedIceWall,
        "rectangle",
        Dict{String, Any}(
            "left" => false
        )
    )
)

Ahorn.minimumSize(entity::AttachedIceWall) = 0, 8
Ahorn.resizable(entity::AttachedIceWall) = false, true

function Ahorn.selection(entity::AttachedIceWall)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, 8, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::AttachedIceWall, room::Maple.Room)
    left = get(entity.data, "left", false)

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 8))
    tileHeight = div(height, 8)

    if left
        for i in 2:tileHeight - 1
            Ahorn.drawImage(ctx, "objects/wallBooster/iceMid00", 0, (i - 1) * 8)
        end

        Ahorn.drawImage(ctx, "objects/wallBooster/iceTop00", 0, 0)
        Ahorn.drawImage(ctx, "objects/wallBooster/iceBottom00", 0, (tileHeight - 1) * 8)

    else
        Ahorn.Cairo.save(ctx)
        Ahorn.scale(ctx, -1, 1)

        for i in 2:tileHeight - 1
            Ahorn.drawImage(ctx, "objects/wallBooster/iceMid00", -8, (i - 1) * 8)
        end

        Ahorn.drawImage(ctx, "objects/wallBooster/iceTop00", -8, 0)
        Ahorn.drawImage(ctx, "objects/wallBooster/iceBottom00", -8, (tileHeight - 1) * 8)

        Ahorn.restore(ctx)
    end
end

# Offset X position so it flips in place
function Ahorn.flipped(entity::AttachedIceWall, horizontal::Bool)
    if horizontal
        entity.left = !entity.left
        entity.x += entity.left ? 8 : -8

        return entity
    end
end

end


text:
using ReachabilityAnalysis, SparseArrays
using ReachabilityAnalysis: add_dimension

LazySets.set_ztol(Float64, 1e-14)

using LazySets: HalfSpace

# physical units of the state variables
# [x₁] = rad
# [x₂] = Nm
# [x₃] = rad
# [x₄] = rad/s
# [x₅] = rad
# [x₆] = rad/s
# [x₇] = rad/s
# [x₈] = rad
# [x₉] = rad/s
# ...
# [x_(2θ+6)] = rad
# [x_(2θ+7)] = rad/s

# physical constants
# indices 'm' (ₘ) and 'l' (ₗ) refer to the motor and the load
# indices 'i' (ᵢ) refer to the numbering of additional rotating masses

const t_init = 0.2  # time to stay in the initial location
const α = 0.03  # backlash size (half of the gap width)
const τ_eng = 0.1  # engine time constant
const γ = 12.0 # gearbox ratio (dimensionless)
const u = 5.0 # requested engine torque

# moments of inertia [kg m²]
const Jₗ = 140.
const Jₘ = 0.3
const Jᵢ = 100.0  # TODO the paper says 0.01, the SpaceEx model uses 100

# viscous friction constants [Nm s/rad]
const bₗ = 5.6
const bₘ = 0.0
const bᵢ = 1.0

# shaft stiffness [Nm/rad]
const kᵢ = 1e5
const kₛ = 1e4

# PID parameters
const k_P = 0.5  # [Nms/rad]
const k_I = 0.5  # [Nm/rad]
const k_D = 0.5  # [Nms²/rad]

function print_dynamics(A, b, location_name)
    println("dynamics of location $location_name:")
    for i in 1:size(A, 1)-1  # ignore the last dimension (time)
        print("x_$i' = ")
        for j in 1:size(A, 2)
            if !iszero(A[i,j])
                print("$(A[i,j]) x_$j + ")
            end
        end
        println("$(b[i])\n")
    end
end

function get_dynamics(kₛ, α, u, n, θ)
    # linear dynamics
    A = spzeros(n, n)

    A[1, 7] = 1.0 / γ
    A[1, 9] = -1.0

    A[2, 1] = (-k_I * γ + k_D * kₛ / (γ * Jₘ)) / τ_eng
    A[2, 2] = (-k_D / Jₘ - 1.0) / τ_eng
    A[2, 3] = k_I * γ / τ_eng
    A[2, 4] = k_P * γ / τ_eng
    A[2, 7] = (-k_P + k_D * bₘ /Jₘ) / τ_eng
    A[2, 8] = -k_I * γ / τ_eng

    A[3, 4] = 1.0

    A[5, 6] = 1.0

    A[6, 5] = -kᵢ / Jₗ
    A[6, 6] = -bₗ / Jₗ
    A[6, 2*θ+6] = kᵢ / Jₗ

    A[7, 1] = -kₛ / (Jₘ * γ)
    A[7, 2] = 1.0 / Jₘ
    A[7, 7] = -bₘ / Jₘ

    i = 8
    while i < n-1
        A[i, i+1] = 1.0

        if i == 8
            # x9 has special dynamics
            A[i+1, 1] = kₛ / Jᵢ
            A[i+1, i] = -kᵢ / Jᵢ
        else
            A[i+1, i-2] = kᵢ / Jᵢ
            A[i+1, i] = -2. * kᵢ / Jᵢ
        end
        A[i+1, i+1] = -bᵢ / Jᵢ

        # wrap-around to x5 in the last step
        j = (i == n-2) ? 5 : i+2
        A[i+1, j] = kᵢ / Jᵢ

        i += 2
    end

    # affine vector
    b = spzeros(n)
    b[2] = k_D * (γ * u - kₛ * α / (Jₘ * γ)) / τ_eng
    b[4] = u
    b[7] = kₛ * α / (Jₘ * γ)
    b[9] = -kₛ * α / Jᵢ
    b[n] = 1.0  # time

    return A, b
end

function get_initial_condition(n, X0_scale)
    c = Vector{Float64}(undef, n)
    g = Vector{Float64}(undef, n)
    c[1:7] = [-0.0432, -11., 0., 30., 0., 30., 360.]
    g[1:7] = [0.0056, 4.67, 0., 10., 0., 10., 120.]
    i = 8
    while i < n
        c[i] = -0.0013
        g[i] = 0.0006
        i += 1
        c[i] = 30.
        g[i] = 10.
        i += 1
    end
    c[n] = 0.0
    g[n] = 0.0
    if X0_scale < 1.0
        g = X0_scale * g
    end

    X0 = Zonotope(c, hcat(g))
    return X0
end

function powertrain(; θ::Int=1, X0_scale::Float64=1.0)
    @assert θ > 0 "θ must be positive, but was $θ"
    @assert (X0_scale > 0.0 && X0_scale <= 1.0) "scale $X0_scale ∉ (0, 1]"

    # dimension of state space (last dimension is time)
    n = 2 * θ + 7 + 1

    # hybrid automaton
    automaton = LightAutomaton(4)

    # negAngle
    A, b = get_dynamics(kₛ, -α, u, n, θ)
    X = HalfSpace(sparsevec([1], [1.], n), -α)  # x1 <= -α
    m_negAngle = @system(x' = A * x + b, x ∈ X)

    # deadzone
    A, b = get_dynamics(0., -α, u, n, θ)
    X = HPolyhedron([HalfSpace(sparsevec([1], [-1.], n), α),  # x1 >= -α
                     HalfSpace(sparsevec([1], [1.], n), α)])  # x1 <= α
    m_deadzone = @system(x' = A * x + b, x ∈ X)

    # posAngle
    A, b = get_dynamics(kₛ, α, u, n, θ)
    X = HalfSpace(sparsevec([1], [-1.], n), -α)  # x1 >= α
    m_posAngle = @system(x' = A * x + b, x ∈ X)

    # negAngleInit
    A, b = get_dynamics(kₛ, -α, -u, n, θ)
    X = HalfSpace(sparsevec([n], [1.], n), t_init)  # t <= t_init
    m_negAngleInit = @system(x' = A * x + b, x ∈ X)

    # modes
    modes = [m_negAngle, m_deadzone, m_posAngle, m_negAngleInit]

    # transition negAngleInit -> negAngle
    add_transition!(automaton, 4, 1, 1)
    guard = HalfSpace(sparsevec([n], [-1.], n), -t_init)  # t >= t_init
    r_41 = ConstrainedIdentityMap(n, guard)

    # transition negAngle -> deadzone
    add_transition!(automaton, 1, 2, 2)
    guard = HalfSpace(sparsevec([1], [-1.], n), α)  # x1 >= -α
    r_12 = ConstrainedIdentityMap(n, guard)

    # transition deadzone -> posAngle
    add_transition!(automaton, 2, 3, 3)
    guard = HalfSpace(sparsevec([1], [-1.], n), -α)  # x1 >= α
    r_23 = ConstrainedIdentityMap(n, guard)

    # TODO the SpaceEx model does not contain the following transitions
#     # transition deadzone -> negAngle
#     add_transition!(automaton, 2, 1, 4)
#     guard = HalfSpace(sparsevec([1], [1.], n), -α)  # x1 <= -α
#     r_21 = ConstrainedIdentityMap(n, guard)
#     # transition posAngle -> deadzone
#     add_transition!(automaton, 3, 2, 5)
#     guard = HalfSpace(sparsevec([1], [1.], n), α)  # x1 <= α
#     r_32 = ConstrainedIdentityMap(n, guard)

    # transition annotations
    resetmaps = [r_41, r_12, r_23]

    # switching
    switchings = [HybridSystems.AutonomousSwitching()]

    H = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition TODO check initial mode
    X0 = get_initial_condition(n, X0_scale)
    initial_condition = [(4, X0)]

    return InitialValueProblem(H, initial_condition)
end

function powertrain_homog(; θ::Int=1, X0_scale::Float64=1.0)
    @assert θ > 0 "θ must be positive, but was $θ"
    @assert (X0_scale > 0.0 && X0_scale <= 1.0) "scale $X0_scale ∉ (0, 1]"

    # dimension of state space (last dimension is time)
    n = 2 * θ + 7 + 1

    # hybrid automaton
    automaton = LightAutomaton(4)

    # negAngle
    A, b = get_dynamics(kₛ, -α, u, n, θ)
    X = HalfSpace(sparsevec([1], [1.], n+1), -α)  # x1 <= -α
    Aext = add_dimension(A, 1)
    Aext[1:n, n+1] .= b
    m_negAngle = @system(x' = Aext * x, x ∈ X)

    # deadzone
    A, b = get_dynamics(0., -α, u, n, θ)
    X = HPolyhedron([HalfSpace(sparsevec([1], [-1.], n+1), α),  # x1 >= -α
                     HalfSpace(sparsevec([1], [1.], n+1), α)])  # x1 <= α
    Aext = add_dimension(A, 1)
    Aext[1:n, n+1] .= b
    m_deadzone = @system(x' = Aext * x, x ∈ X)

    # posAngle
    A, b = get_dynamics(kₛ, α, u, n, θ)
    X = HalfSpace(sparsevec([1], [-1.], n+1), -α)  # x1 >= α
    Aext = add_dimension(A, 1)
    Aext[1:n, n+1] .= b
    m_posAngle = @system(x' = Aext * x, x ∈ X)

    # negAngleInit
    A, b = get_dynamics(kₛ, -α, -u, n, θ)
    X = HalfSpace(sparsevec([n], [1.], n+1), t_init)  # t <= t_init
    Aext = add_dimension(A, 1)
    Aext[1:n, n+1] .= b
    m_negAngleInit = @system(x' = Aext * x, x ∈ X)

    # modes
    modes = [m_negAngle, m_deadzone, m_posAngle, m_negAngleInit]

    # transition negAngleInit -> negAngle
    add_transition!(automaton, 4, 1, 1)
    guard = HalfSpace(sparsevec([n], [-1.], n+1), -t_init)  # t >= t_init
    r_41 = ConstrainedIdentityMap(n, guard)

    # transition negAngle -> deadzone
    add_transition!(automaton, 1, 2, 2)
    guard = HalfSpace(sparsevec([1], [-1.], n+1), α)  # x1 >= -α
    r_12 = ConstrainedIdentityMap(n, guard)

    # transition deadzone -> posAngle
    add_transition!(automaton, 2, 3, 3)
    guard = HalfSpace(sparsevec([1], [-1.], n+1), -α)  # x1 >= α
    r_23 = ConstrainedIdentityMap(n, guard)

    # transition annotations
    resetmaps = [r_41, r_12, r_23]

    # switching
    switchings = [HybridSystems.AutonomousSwitching()]

    H = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition TODO check initial mode
    X0 = get_initial_condition(n, X0_scale)
    @assert isa(X0, Zonotope)
    X0 = Zonotope(vcat(X0.center, 1.0), vcat(X0.generators, 0.0))
    initial_condition = [(4, X0)]

    return InitialValueProblem(H, initial_condition)
end

const IA = ReachabilityAnalysis.IntervalArithmetic
using ReachabilityAnalysis: post

# requirement
function req(sol)
    idx = findall(x -> 2.0 ∈ tspan(x), sol.Fk.u)
    !isempty(idx) && all(x -> location(x) == 3, sol.Fk.u[idx])
end

function _solve_powertrain(prob)

    # solve mode 4
    p1 = IVP(mode(prob.s, 4), prob.x0[1][2])
    alg = GLGM06(δ=0.001, approx_model=Forward(setops=:zono))
    sol1F = post(alg, p1, IA.Interval(0, 2.0))
    sol1F.ext[:loc_id] = 4
    aux = sol1F(t_init) # clustered sets
    aux2 = overapproximate(ConvexHullArray([set(Xi) for Xi in aux]), Zonotope)

    # solve mode 1
    p2 = IVP(mode(prob.s, 1), aux2)
    Δt0 = tspan(aux)
    tprev = tstart(aux)
    sol2F = post(alg, p2, IA.Interval(0, 2.0 - tprev); Δt0=Δt0)
    sol2F.ext[:loc_id] = 1

    auxvec = zeros(dim(sol2F)) # = n + 1
    auxvec[1] = 1.0
    auxh = Hyperplane(auxvec, -α) # x₁ = -0.03
    jump_rset_idx = findall(R -> !isdisjoint(R, auxh), sol2F)
    sol2F_rset = view(sol2F, jump_rset_idx)
    aux3 = overapproximate(ConvexHullArray([set(Xi) for Xi in sol2F_rset]), Zonotope)

    # solve mode 2
    p3 = IVP(mode(prob.s, 2), aux3)
    alg2 = GLGM06(δ=0.001, approx_model=Forward(setops=:zono), disjointness_method=BoxEnclosure())
    Δt0 = tspan(sol2F_rset)
    tprev = tstart(sol2F_rset)
    sol3F = post(alg2, p3, IA.Interval(0, 2.0 - tprev); Δt0=Δt0)
    sol3F.ext[:loc_id] = 2
    auxh = Hyperplane(auxvec, α) # x₁ = 0.03
    jump_rset_idx = findall(R -> !isdisjoint(R, auxh), sol3F)
    sol3F_rset = view(sol3F, jump_rset_idx)
    aux4 = overapproximate(ConvexHullArray([set(Xi) for Xi in sol3F_rset]), Zonotope)

    # solve mode 3
    p4 = [IVP(mode(prob.s, 3), set(Xi)) for Xi in sol3F_rset]
    Δt0 = tspan(sol3F_rset)
    tprev = tstart(sol3F_rset)
    sol4F = [post(alg, p4i, IA.Interval(0, 2.0 - tprev); Δt0=Δt0) for p4i in p4]
    for F in sol4F
        F.ext[:loc_id] = 3
    end

    return HybridFlowpipe(reduce(vcat, (sol1F, sol2F, sol3F, sol4F)))
end


text:
#==============================================================================#
# IAM/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSIAM
using Base.Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config(region="ap-southeast-2")

#-------------------------------------------------------------------------------
# IAM tests
#-------------------------------------------------------------------------------

@test ismatch(r"^arn:aws:iam::[0-9]+:[^:]+$", iam_whoami(aws))

#==============================================================================#
# End of file.
#==============================================================================#


text:
using Unitful, Gadfly, DataFrames

set_default_plot_size(6inch, 6inch)

a = -9.81u"m/s^2"
t = (1:0.5:10)u"s"
v = a .* t
v2 = v .+ 2u"m/s"
h = 0.5 * a .* t.^2
df = DataFrame(time=t, velocity=v, position=h,
               unitlesst=ustrip.(t), unitlessv=ustrip.(v), unitlessh=ustrip.(h),
               position2=reverse(h))

# test basics of point/line plots with Unitful
gridstack(reshape([plot(df, x=:time, y=:velocity),
                   plot(df, x=:position, y=:time),
                   plot(df, x=:unitlesst, y=:velocity),
                   plot(df, x=:time, y=:unitlessh),
                   plot(df, x=:time, y=:position, Geom.line),
                   plot(df, layer(x=:time, y=:position, Geom.line),
                        layer(x=:time, y=:position2))], (3,2)))


text:
function findλmax(IT::iterator, xdy, d::dataCD)::Vector{Float64}

    λmax = Vector{Float64}(TOTALCOMPONENTS)
    temp = 0.0
    maxPos = 0

    for i in IT.components

        λmax[i] = maximum(abs(xdy[i]))
        # temp, k = findmax(abs(xdy[i]))#.*d.σ[i]))
        # temp = temp #/ d.σ[i][k]

        # if temp > λmax
        #   λmax = temp
        # end
    end

    # print(λmax)

    return λmax / IT.obs
end

function compute_λ_path(IT::iterator, xdy, numλ::Int, d::dataCD; logarit::Bool=true)::Tuple{Vector{Float64}, Vector{Float64}}

    λ_max = findλmax(IT, xdy, d)

    if logarit
        path = collect(logspace(1, 0.001, numλ)-1.0)/9.0
    else
        path = collect(linspace(1, 0.001, numλ))
    end
    return path, λ_max
end

function compute_γ_path(IT::iterator, xdy, numγ::Int, d::dataCD; logarit=true)::Vector{Float64}

    γ_max = 2

    if logarit
        vec = γ_max * collect(logspace(1, 0.001, numγ) - 1.0) / 9.0
    else
        vec = γ_max * collect(linspace(1, 0.001, numγ))
    end

    return vec
end

function compute_BIC(y_hat::Vector{Float64}, y::Vector{Float64}, β, IT::iterator; ɛ = 1e-5::Float64)::Float64

    BIC = 0.0
    err = zeros(y)

    N = IT.obs

    @inbounds for i in 1:IT.obs
        err[i] = y[i] - y_hat[i]
    end

    k = 0::Int64
    p = 0::Int64

    @inbounds for i in IT.components
        for j in IT.elements[i]
            p += 1
            if abs(β[i][j]) != 0.0
                k += 1
            end
        end
    end

    BIC = N * log(var(err)) + k * log(N) + convert(Float64, log(binomial(BigInt(p), BigInt(k))))

    return BIC
end

function compute_BIC(
    y::Vector{Float64}, β, activeSet, IT::iterator, d, xdy, ɛ::Float64 = 1e-5, std::Int = 1,
    lower_bounds::Vector{Float64} = -10e+7*ones(TOTALCOMPONENTS), upper_bounds::Vector{Float64} = 10e+7*ones(TOTALCOMPONENTS)
    )

    β_ols = compute_OLS(β, activeSet, IT, xdy, d, lower_bounds, upper_bounds)
    y_hat, β_new = compute_estimate(zeros(y), IT, β, d)

    if std == 1
        BIC = compute_BIC(y_hat::Vector{Float64}, y::Vector{Float64}, β_ols, IT)
    else
        BIC = compute_BIC(y_hat::Vector{Float64}, y::Vector{Float64}, β_ols, IT)
    end

    return BIC, y_hat
end

function compute_OLS(β, activeSet, IT::iterator, xdy, d::dataCD,
  lower_bounds::Vector{Float64}, upper_bounds::Vector{Float64})

    β_ols = deepcopy(β)

    const N = IT.obs

    @inbounds for c1 in IT.components
        @inbounds for j in IT.elements[c1]
            if β[c1][j] != 0.0
                partial_fit = 0.0
                @inbounds for c2 in IT.components
                    @inbounds for l in IT.elements[c2]
                        if β[c2][l] != 0.0 && (c1, j) != (c2, l)
                          @inbounds partial_fit += GM2(c1, c2, j, l, d, N) * β[c2][l]
                        end
                    end
                end

                # univariate ordinary leasts squares coefficient
                β_ols[c1][j] = (xdy[c1][j] - partial_fit) / (IT.obs * (d.σ[c1][j]^2))
            else
                β_ols[c1][j] = 0.0
            end

            β_ols[c1][j] = max(β_ols[c1][j], lower_bounds[c1])
            β_ols[c1][j] = min(β_ols[c1][j], upper_bounds[c1])

        end
    end

    return β_ols
end

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
"""
    x = axb(a,b)

Solve the linear equation ``ax=b`` with scalar univariate polynomials.

In the equation `a` and `b` are given polynomials in some variable,say, `s`, and `x` is a polynomial to be determined. The equation can be interpretted as dividing the the polynomial `b` by `a` to obtain a polynomial `x`. Obviously this will not always be possible, in which case the function will return `nothing`.

# Examples

```julia
julia> a = Polynomial([1,2,3],:s);
julia> b = a*Polynomial([1,2],:s);
julia> x = axb(a,b)
Polynomial(0.9999999999999997 + 2.0000000000000004*s)

julia> a*x≈b
true

julia> a = Polynomial([1,2],:s);
julia> b = Polynomial([1,2,3],:s);
julia> x = axb(a,b);
julia> isnothing(x)
true
```
"""
function axb(a::Polynomial,b::Polynomial)
    da = degree(a)
    db = degree(b)
    dx = db-da
    A = ltbtmatrix(a,dx+1)              # The matrix A will be (generally) tall rectangular.
    bc = coeffs(b)
    F = qr(A)                           # Instead of just solving the least squares problem,
    R = F.R                             # we want to check if an exact solution exists
    qb = F.Q'*bc                        # and do it efficiently (not by substituting into Ax=b).
    if norm(qb[dx+2:end],Inf) < 1e-8    # Here should be perhaps something related to eps.
        xc = R\(qb[1:dx+1])             # will Julia's backslash operator take advantage of triangularity of R? Perhaps yes, R has the "triangular" flag.
        x = Polynomial(xc,a.var)
    else
        x = nothing
    end
    return x
end

"""
    x,y = axby0(a,b)

Solve the homogeneous linear Diophantine equation ``ax+by=0`` with scalar univariate polynomials.

The problem corresponds to finding the minimum-order representation ``x/y`` of the fraction ``b/a`` (up to the minus sign), that is, it solves ``b/a = -x/y`` for a given fraction ``b/a`` such that `y` is of minimum possible degree.

# Examples

```julia
julia> a = Polynomial([1,2,3],:s); b = Polynomial([4,5],:s); c = Polynomial([1,1],:s);
julia> a = a*c; b = b*c;
julia> x, y = axby0(a,b)
(Polynomial(-0.5393598899705949 - 0.6741998624632413*s), Polynomial(0.1348399724926485 + 0.2696799449852973*s + 0.4045199174779448*s^2))
```
"""
function axby0(a::Polynomial,b::Polynomial)
    da = degree(a)
    db = degree(b)
    dx = db-1             # no need to start with deg(x) = deg(b) because such solution is guaranteed to exist:
    x = -b
    y = a
    dx = degree(x)
    dy = degree(y)
    S = sylvestermatrix(a,b,(dx+1,dy+1))
    xy = nullspace(S)
    ds = size(xy,2)
    #println("Degree of x: $dx  Dimension of nullspace of S: $ds")
    while !isempty(xy)
        x = Polynomial(xy[1:dx+1],a.var)
        y = Polynomial(xy[dx+2:end],a.var)
        dx = dx-1         # simply reducing the anticipated degrees of the solutions
        dy = dy-1
        S = sylvestermatrix(a,b,(dx+1,dy+1))
        xy = nullspace(S)
        ds = size(xy,2)
        #println("Degree of x: $dx  Dimension of nullspace of S: $ds")
    end
    return (x,y)
end

"""
    x,y = axbyc(a,b,c[;deg=:miny])

Solve the linear Diophantine equation ``ax+by=c`` with univariate polynomials.

The extra input parameter `deg` gives some specification of degrees of the two polynomials `x` and `y`. By default, `deg=:miny` specifies that a solution minimizing the degree of `y` is sought. If `deg=:minx`, the degree of `x` is minimized.

The problem corresponds to finding a (negative) feedback controller given by the transfer function ``y/x`` for a system modelled by the transfer function ``b/a`` such the that closed-loop characteristic polynomial is `c`. Choosing a solution that minimizes the degree of `y`, causal controller `y/x` is obtained.

# Examples

```julia
julia> a = Polynomial([1,2,3],:s);
julia> b = Polynomial([4,5],:s);
julia> c = Polynomial([5,6,7,8,9],:s);

julia> x,y = axbyc(a,b,c,deg=:miny)
(Polynomial(1.8484848484848486 + 0.6666666666666666*s + 3.0*s^2), Polynomial(0.7878787878787878 - 0.5757575757575759*s))

julia> a*x+b*y≈c
true

julia> x,y = axbyc(a,b,c,deg=:minx)
(Polynomial(0.6909090909090905), Polynomial(1.3272727272727274 - 0.25454545454545435*s + 1.8*s^2))

julia> a*x+b*y≈c
true
```
"""
function axbyc(a::Polynomial,b::Polynomial,c::Polynomial;deg=:miny)
    da = degree(a)
    db = degree(b)
    dc = degree(c)
    if deg == :miny
        dy = da-1                       # degree of y < degree of a
        m2 = db+1+dy                    # minimum number of rows of the b-block in the Sylvester
        m = max(da+1,m2,dc+1)           # number of rows of the whole Sylvester matrix
        n1 = m-da                       # number of columns of the a-block in the Sylvester matrix
        dx = n1-1                       # degree of x
        w = (n1,dy+1)
    elseif deg == :minx
        dx = db-1
        m1 = da+1+dx
        m = max(db+1,m1,dc+1)
        n2 = m-db
        dy = n2-1
        w = (dx+1,n2)
    else
        throw(ArgumentError("deg argument with unacceptable value"))
    end
    cc = zeros(Float64,m)
    cc[1:dc+1] = coeffs(c)
    S = sylvestermatrix(a,b,w)          # The matrix S will be square but may be singular.
    xy = S\cc                           # Could be perhaps done more efficiently, see axb().
    x = xy[1:dx+1]
    y = xy[dx+2:end]
    x = Polynomial(x,a.var)
    y = Polynomial(y,a.var)
    if a*x+b*y≈c
        return (x,y)
    else
        return (nothing,nothing)
    end
end

"""
    x = axaxbb(a,b)

Solve the symmetric linear Diophantine equation ``ãx+ax̃=b+b̃`` with univariate polynomials.

The conjugation is understood with respect to the imaginary axis, that is, ã(s)=a(-s) for polynomials with real coefficients. For complex polynomials, the coefficients are additionally complex conjugated.

# Examples

```julia
julia> a = Polynomial([-0.12, -0.29, 1],:s);
julia> b = Polynomial([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:s);
julia> x = axaxbb(a,b)
Polynomial(-15.50000000000003 + 50.0096551724139*s + 1.19*s^2)
```
"""
function axaxbb(a::Polynomial,b::Polynomial)
    da = degree(a)
    db = degree(b)
    bb = b+cconj(b)
    dbb = degree(bb)
    dx = dbb - da
    ã = cconj(a)
    T̃ = ltbtmatrix(ã,dx+1)
    i = fill(1,dx+1)
    i[2:2:end] .= -1
    T = ltbtmatrix(a,dx+1)
    T =  T*Diagonal(i)
    A = T̃+T
    x = A\coeffs(bb)
    x = Polynomial(x,a.var)
    return x
end

"""
    x,y = axbycd(a,b,c)

Solve the symmetric linear Diophantine equation ``ãx+b̃y=c+d̃`` with univariate polynomials.

The conjugation is understood with respect to the imaginary axis, that is, ã(s)=a(-s) for polynomials with real coefficients. For complex polynomials, the coefficients are additionally complex conjugated.

# Examples

```julia
julia> a = Polynomial([-0.12, -0.29, 1],:s);
julia> b = Polynomial([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:s);
julia> x = axxabb(a,b)
Polynomial(-15.50000000000003 + 50.0096551724139*s + 1.19*s^2)
```
"""
function axbycd(a::Polynomial,b::Polynomial,c::Polynomial,d::Polynomial)
    da = degree(a)
    db = degree(b)
    dc = degree(c)
    dd = degree(d)
    ccd = c+cconj(d)
    dccd = degree(ccd)
    dx = db-1
    dy = da-1
    Sa = ltbtmatrix(a,dx+1)
    i = fill(1.0,dx+1)
    i[2:2:end] .= -1.0
    E = Diagonal(i)
    Sa =  Sa*E
    b̃ = cconj(b)
    Sb = ltbtmatrix(b̃,dy+1)
    S = [Sa Sb]
    cxy = S\coeffs(ccd)
    cx = cxy[1:dx+1]
    cy = cxy[dx+2:end]
    x = Polynomial(cx,a.var)
    y = Polynomial(cy,a.var)
    return x,y
end

"""
    axycminl1(a,c[,dymax=100,emax=1e-4])

Solve the linear Diophantine equation ``ax+y=c`` for given univariate polynomials `a` and `c`, that is, find two univariate polynomials `x` and `y` satisfying the equation and, moreover, with the 1-norm of the vector of coefficients of the polynomial `y` minimized.

The degree of the minimizing polynomial `y` is not known apriori. It is, however, known that such minimizing polynomials of a finite degree exists. The solver then starts with smalest possible degree of `y` and goes up to `dmax` while checking if the improvement with respect to the previous iteration is better than `emax`.

# Examples

```julia
julia> a = Poly([-0.12, -0.29, 0.71, 2.79, 2.92, 1],:d)
Poly(-0.12 - 0.29*d + 0.71*d^2 + 2.79*d^3 + 2.92*d^4 + 1.0*d^5)

julia> c = Poly([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:d)
Poly(1.86 - 0.34*d - 1.14*d^2 - 0.21*d^3 + 1.19*d^4 - 1.12*d^5)

julia> x, y, val = axycminl1(a,c);
```
Verify that the solutions are truly satisfying the original equation
```julia
julia> a*x+y≈c
true
```
Note that the degree of the resulting polynomial `y` might look high but it is only an artefact of the numerical solution.
```julia
julia> degree(chop(y,atol=1e-3))
```

# References

Hurak, Z., A. Böttcher, and M. Sebek, "Minimum Distance to the Range of a Banded Lower Triangular Toeplitz Operator in l1 and Application in l1-Optimal Control", SIAM Journal on Control and Optimization, vol. 45, issue 1: SIAM, pp. 107-122, 2006. [DOI: 10.1137/S0363012903437940](http://dx.doi.org/10.1137/S0363012903437940).
"""
function axycminl1(a::Polynomial,c::Polynomial;dymax=100,emax=1e-4)
    da = degree(a)                      # degree of a
    dc = degree(c)                      # degree of c
    dymin = max(da,dc)                  # initial degree of y
    b = append!(copy(coeffs(c)),zeros(dymin-dc)) # initial length of the right hand column
    l1y_pre = Inf                       # initial optimal value (l1 norm of the coefficients of y)
    for dy = dymin:dymax
        dx = dy-da
        T = ltbtmatrix(a,dx+1)          # lower triangular banded Toeplitz matrix
        model = Model(Clp.Optimizer)
        @variable(model, x[1:dx+1])     # coefficients of the x polynomial
        @variable(model, y⁺[1:dy+1]>=0) # coefficients of the y polynomial, the plus part
        @variable(model, y⁻[1:dy+1]>=0) # coefficients of the y polynomial, the minus part
        @constraint(model, con, T*x + y⁺ - y⁻ .== b)
        @objective(model, Min, sum(y⁺+y⁻))
        set_silent(model::Model)
        optimize!(model)
        l1y = objective_value(model)
        xₛ = value.(x)
        y⁺ₛ = value.(y⁺)
        y⁻ₛ = value.(y⁻)
        yₛ = y⁺ₛ-y⁻ₛ
        @printf("|y|₁=%f \t\t Deg(x)=%d \t Deg(y)=%d\n",l1y,dx,dy)   # print the value just for diagnostics
        b = append!(b,0)                # extend the b for the next iteration
        global xₛ
        global yₛ
        global l1y
        e = l1y_pre-l1y                 # diff. between the value of the current and the previous
        abs(e)<emax && break            # if not improved enough, break
        l1y_pre = l1y
    end
    x = Polynomial(xₛ,a.var)
    y = Polynomial(yₛ,a.var)
    return (x,y,l1y)
end


text:
# const RealOrMissing = Union{Missing,Real}
"""
  AbstractNaArray = AbstractArray{<:Union{T, Missing}} where T <: Real
  AbstractNanArray = AbstractArray{Union{T, Missing}} where T <: Real

- `AbstractNanArray`: must have missing value
- `AbstractNaArray` : may have missing value

# Examples
```
f(x::AbstractNanArray) = x
# f([1, 2]) # not work
f([1, 2, missing])
f([missing])

f2(x::AbstractNaArray) = x
f2([1, 2])
f2([1, 2, missing])
f2([missing])
```
"""
# 可以含有missing、也可以不含有
AbstractNaArray = AbstractArray{<:Union{T,Missing}} where {T<:Real}
# 必须含有missing
AbstractNanArray = AbstractArray{Union{T,Missing}} where {T<:Real}

to_missing(x::AbstractNanArray) = x

function to_missing(x::AbstractArray{<:Real}, replacement=0)
  x2 = zeros(Union{eltype(x),Missing}, size(x)...)
  x2 .= x
  for i = eachindex(x2)
    if x[i] == replacement
      x2[i] = missing
    end
  end
  x2
end

function replace_value!(x::AbstractNaArray, org, new)
  for i in eachindex(x)
    if !ismissing(x[i]) && x[i] == org
      x[i] = new
    end
  end
  x
end

"""
  to_missing(x::AbstractNanArray)
  to_missing(x::AbstractArray{<:Real}, replacement = 0)

  drop_missing(x::AbstractArray{<:Real}, replacement = 0)
  drop_missing(x::AbstractNanArray, replacement=0)
"""
drop_missing(x::AbstractArray{<:Real}, replacement=0) = x;

function drop_missing(x::AbstractNanArray, replacement=0)
  x2 = ones(eltype(x).b, size(x)...) .* replacement
  for i = eachindex(x)
    if !ismissing(x[i])
      x2[i] = x[i]
    end
  end
  x2
end

## for `AbstractNaArray`: should use `===` to compare
# Base.:-, Base.isequal
function Base.isequal(a::AbstractNaArray, b::AbstractNaArray)
  size(a) != size(b) && return (false)
  for i = eachindex(a)
    if (a[i] !== b[i])
      return (false)
    end
  end
  return (true)
end

export AbstractNaArray, AbstractNanArray,
  isequal,
  replace_value!,
  drop_missing, to_missing;


text:
using MetEP
using MetEP.MetLP
using MetEP.MetNets
using MetEP.ProgressMeter
using Test

@testset "MetEP.jl" begin
    include("maxent_ep_test.jl")
end


text:
using Compat

import Compat: String

depsfile = joinpath(dirname(@__FILE__),"deps.jl")
if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    open(depsfile,"w") do f
        print(f,"const libcplex = ")
        show(f, path) # print with backslashes excaped on windows
        println(f)
    end
end

if is_apple()
    Libdl.dlopen("libstdc++",Libdl.RTLD_GLOBAL)
end

cpxvers = ["125","1251","1260","1261","1262","1263"]

libnames = String["cplex"]
for v in reverse(cpxvers)
    if is_apple()
        push!(libnames, "libcplex$v.dylib")
    elseif is_unix()
        push!(libnames, "libcplex$v.so")
    end
end

wincpxvers = ["126","1261","1262","1263"]
if is_windows()
    for v in reverse(wincpxvers)
        env = "CPLEX_STUDIO_BINARIES$v"
        if haskey(ENV,env)
            for d in split(ENV[env],';')
                contains(d,"cplex") || continue
                if v == "126" # annoying inconsistency
                    push!(libnames,joinpath(d,"cplex1260"))
                else
                    push!(libnames,joinpath(d,"cplex$(v)"))
                end
            end
        end
    end
end

found = false

for l in libnames
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        found = true
        write_depsfile(Libdl.dlpath(d))
        break
    end
end

if !found
    error("Unable to locate CPLEX installation. Note this must be downloaded separately. See the CPLEX.jl README for further instructions.")
end


text:
module HeatEquationTests

using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

import Gridap: ∇
import GridapODEs.TransientFETools: ∂t

θ = 0.2

u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
u(t::Real) = x -> u(x,t)
f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 2

reffe = ReferenceFE(lagrangian,Float64,order)
V0 = FESpace(
  model,
  reffe,
  conformity=:H1,
  dirichlet_tags="boundary"
)
U = TransientTrialFESpace(V0,u)

Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)

#
a(u,v) = ∫(∇(v)⋅∇(u))dΩ
b(v,t) = ∫(v*f(t))dΩ

res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t)
jac(t,u,du,v) = a(du,v)
jac_t(t,u,dut,v) = ∫(dut*v)dΩ

op = TransientFEOperator(res,jac,jac_t,U,V0)

t0 = 0.0
tF = 1.0
dt = 0.1

U0 = U(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)

ls = LUSolver()
using Gridap.Algebra: NewtonRaphsonSolver
nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
ode_solver = ThetaMethod(ls,dt,θ)

sol_t = solve(ode_solver,op,uh0,t0,tF)

# Juno.@enter Base.iterate(sol_t)

l2(w) = w*w

tol = 1.0e-6
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
end

end #module


text:
module SWResample
  include(joinpath(Pkg.dir("VideoIO"),"src","init.jl"))
  w(f) = joinpath(swresample_dir, f)

  using AVUtil
  include(w("LIBSWRESAMPLE.jl"))
end


text:
module RGFPrecompTest
    using RuntimeGeneratedFunctions
    RuntimeGeneratedFunctions.init(@__MODULE__)

    f = @RuntimeGeneratedFunction(:((x,y)->x+y))
end


text:
module Utils

export chain, flatmap, snakecase_from_camelcase, @verify

function chain(itrs...)
    Base.Iterators.flatten(itrs)
end

function flatmap(f, c)
    map(f, c) |> Base.Iterators.flatten
end

function snakecase_from_camelcase(str::String)
    stringbuilder = IOBuffer()
    inuppercase = false
    for c in str
        if !inuppercase && isuppercase(c)
            print(stringbuilder, '_')
            inuppercase = true
        elseif inuppercase && islowercase(c)
            inuppercase = false
        end
        print(stringbuilder, lowercase(c))
    end
    String(take!(stringbuilder))
end

macro verify(expr, msgs...)
    msg = isempty(msgs) ? string(expr) : :(join($(esc(msgs)), " "))
    :($(esc(expr)) ? nothing : error($msg))
end

end


text:
# To ensure that tests can run on travis we have to do a little
# hackadoodle here. The tests require a license file. We include
# a license file that is only valid for one day (the day when
# change is submitted).
# If there is no valid license file, we default to that file.

using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const FALLBACK_URL = "mosek://solve.mosek.com:30080"

using MosekTools
const optimizer = Mosek.Optimizer(fallback = FALLBACK_URL)
MOI.set(optimizer, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "Mosek"
end

@testset "Parameters" begin
    optimizer = Mosek.Optimizer(fallback = FALLBACK_URL)
    @testset "Double Parameter" begin
        MOI.set(optimizer, MOI.RawParameter("INTPNT_CO_TOL_DFEAS"), 1e-7)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-7
        MOI.set(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1e-8)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-8
        @testset "with integer value" begin
            MOI.set(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1)
            @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1
        end
    end
    @testset "Integer Parameter" begin
        MOI.set(optimizer, MOI.RawParameter("MSK_IPAR_INTPNT_MAX_ITERATIONS"), 100)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 100
        MOI.set(optimizer, MOI.RawParameter("INTPNT_MAX_ITERATIONS"), 200)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 200
        @testset "with enum value" begin
            MOI.set(optimizer, MOI.RawParameter("MSK_IPAR_OPTIMIZER"), MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
            @test MOI.get(optimizer, MOI.RawParameter("MSK_IPAR_OPTIMIZER")) == convert(Int32, MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
        end
    end
    @testset "String Parameter" begin
        MOI.set(optimizer, MOI.RawParameter("PARAM_WRITE_FILE_NAME"), "foo.txt")
        # Needs https://github.com/JuliaOpt/Mosek.jl/pull/174
        #@test MOI.get(optimizer, MOI.RawParameter("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "foo.txt"
        MOI.set(optimizer, MOI.RawParameter("MSK_SPAR_PARAM_WRITE_FILE_NAME"), "bar.txt")
        #@test MOI.get(optimizer, MOI.RawParameter("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "bar.txt"
    end
    @testset "TimeLimitSec" begin
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
        MOI.set(optimizer, MOI.TimeLimitSec(), 1.0)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_OPTIMIZER_MAX_TIME")) == 1.0
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === 1.0
        MOI.set(optimizer, MOI.TimeLimitSec(), nothing)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
    end
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(optimizer, false)
    @test MOIU.supports_default_copy_to(optimizer, true)
end

const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3)

@testset "Basic" begin
    @testset "Linear" begin
        MOIT.basic_constraint_tests(
            optimizer, config,
            include=[
                (MOI.SingleVariable, MOI.EqualTo{Float64}),
                (MOI.SingleVariable, MOI.LessThan{Float64}),
                (MOI.SingleVariable, MOI.GreaterThan{Float64}),
                (MOI.SingleVariable, MOI.Interval{Float64}),
                (MOI.SingleVariable, MOI.ZeroOne),
                (MOI.SingleVariable, MOI.Integer),
                (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}),
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}),
                (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
        ])
    end
    @testset "Conic" begin
        MOIT.basic_constraint_tests(
            optimizer, config,
            include=[
                (MOI.VectorOfVariables, MOI.ExponentialCone),
                (MOI.VectorOfVariables, MOI.DualExponentialCone),
                (MOI.VectorOfVariables, MOI.PowerCone{Float64}),
                (MOI.VectorOfVariables, MOI.DualPowerCone{Float64}),
                (MOI.VectorOfVariables, MOI.SecondOrderCone),
                (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
        ])
    end
end

const bridged = MOIB.full_bridge_optimizer(optimizer, Float64)

# Mosek errors during `MOI.set` instead of `MOI.get` when there are duplicates.
#@testset "Name" begin
#    MOIT.nametest(bridged)
#end

@testset "Copy" begin
    model = MOIB.full_bridge_optimizer(Mosek.Optimizer(), Float64)
    MOIT.copytest(bridged, model)
end

@testset "Start" begin
    # TODO this should be checked somewhere in MOI
    @test MOI.supports(bridged, MOI.VariablePrimalStart(), MOI.VariableIndex)
end

@testset "Unit" begin
    # Mosek does not support names
    MOIT.unittest(bridged, config, [
        # TODO
        "number_threads",
        # Find objective bound of 0.0 which is lower than 4.0
        "solve_objbound_edge_cases",
        # Cannot put multiple bound sets of the same type on a variable
        "solve_integer_edge_cases",
        # Cannot mix `ZeroOne` with `GreaterThan`/`LessThan`
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3"])
end

@testset "Continuous Linear" begin
    MOIT.contlineartest(bridged, config)
end

@testset "Continuous Quadratic" begin
    MOIT.contquadratictest(bridged, config, [
        # Non-convex
        "ncqcp",
        # QuadtoSOC does not work as the matrix is not SDP
        "socp"
    ])
end

@testset "Continuous Conic" begin
    MOIT.contconictest(bridged, config, ["rootdets", "logdet"])
end

@testset "Integer Linear" begin
    MOIT.intlineartest(bridged, config, [
        # SOS constraints not supported:
        # int2 uses SOS and indicator can be bridged to SOS
        "int2", "indicator1", "indicator2", "indicator3", "indicator4"
    ])
end


text:
ti = time()

include("discrete_systems.jl")
include("discrete_lyapunov.jl")
include("continuous_systems.jl")
include("continuous_lyapunov.jl")

ti = time() - ti
println("Test took total time of:")
println(round(ti, 3), " seconds or ", round(ti/60, 3), " minutes")


text:
module Cthulhu

Base.Experimental.@compiler_options compile=min optimize=1

using CodeTracking: definition, whereis
using InteractiveUtils
using UUIDs
using REPL: REPL, AbstractTerminal

using Core: MethodInstance
const Compiler = Core.Compiler
import Core.Compiler: MethodMatch, LimitedAccuracy, ignorelimited, specialize_method
import Base: unwrapva, isvarargtype, unwrap_unionall, rewrap_unionall
const mapany = Base.mapany

# branch on https://github.com/JuliaLang/julia/pull/42125
# TODO remove me once v1.7 is released
@static if isdefined(Base, Symbol("@constprop"))
    import Base: @constprop
else
    macro constprop(_, ex); esc(ex); end
end

Base.@kwdef mutable struct CthulhuConfig
    enable_highlighter::Bool = false
    highlighter::Cmd = `pygmentize -l`
    asm_syntax::Symbol = :att
    dead_code_elimination::Bool = true
    pretty_ast::Bool = false
end

"""
    Cthulhu.CONFIG

# Options
- `enable_highlighter::Bool`: Use command line `highlighter` to syntax highlight
  Julia, LLVM and native code.
- `highlighter::Cmd`: A command line program that receives "julia" as an argument and julia
   code as stdin. Defaults to `$(CthulhuConfig().highlighter)`.
- `asm_syntax::Symbol`: Set the syntax of assembly code being used.
  Defaults to `$(CthulhuConfig().asm_syntax)`.
- `dead_code_elimination::Bool`: Enable dead-code elimination for high-level Julia IR.
  Defaults to `true`. DCE is known to be buggy and you may want to disable it if you
  encounter errors. Please report such bugs with a MWE to Julia or Cthulhu.
- `pretty_ast::Bool`: Use a pretty printer for the ast dump. Defaults to false.
"""
const CONFIG = CthulhuConfig()

module DInfo
    @enum DebugInfo none compact source
end
using .DInfo: DebugInfo
const AnyDebugInfo = Union{DebugInfo,Symbol}

include("interpreter.jl")
include("callsite.jl")
include("reflection.jl")
include("ui.jl")
include("codeview.jl")
include("backedges.jl")

export descend, @descend, descend_code_typed, descend_code_warntype, @descend_code_typed, @descend_code_warntype
export ascend

"""
    @descend_code_typed

Evaluates the arguments to the function or macro call, determines their
types, and calls `code_typed` on the resulting expression.
"""
macro descend_code_typed(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_typed, ex0)
end

"""
    @interp

For debugging. Returns a CthulhuInterpreter from the appropriate entrypoint.
"""
macro interp(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :mkinterp, ex0)
end

"""
    @descend_code_warntype

Evaluates the arguments to the function or macro call, determines their
types, and calls `code_warntype` on the resulting expression.
"""
macro descend_code_warntype(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_warntype, ex0)
end

"""
    @descend

Shortcut for [`@descend_code_typed`](@ref).
"""
const var"@descend" = var"@descend_code_typed"

"""
    descend_code_typed(f, argtypes=Tuple{}; kwargs...)
    descend_code_typed(tt::Type{<:Tuple}; kwargs...)
    descend_code_typed(Cthulhu.BOOKMARKS[i]; kwargs...)
    descend_code_typed(mi::MethodInstance; kwargs...)

Given a function and a tuple-type, interactively explore the output of
`code_typed` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
julia> descend_code_typed(sin, (Int,))
[...]

julia> descend_code_typed(sin, Tuple{Int})
[...]

julia> descend_code_typed(Tuple{typeof(sin), Int})
[...]

julia> descend_code_typed() do
           T = rand() > 0.5 ? Int64 : Float64
           sin(rand(T)
       end
[...]
```
"""
descend_code_typed(@nospecialize(args...); kwargs...) =
    _descend_with_error_handling(args...; iswarn=false, kwargs...)

"""
    descend_code_warntype(f, argtypes=Tuple{}; kwargs...)
    descend_code_warntype(tt::Type{<:Tuple}; kwargs...)
    descend_code_warntype(Cthulhu.BOOKMARKS[i])
    descend_code_warntype(mi::MethodInstance; kwargs...)

Given a function and a tuple-type, interactively explore the output of
`code_warntype` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
julia> descend_code_warntype(sin, (Int,))
[...]

julia> descend_code_warntype(sin, Tuple{Int})
[...]

julia> descend_code_warntype(Tuple{typeof(sin), Int})
[...]

julia> descend_code_warntype() do
           T = rand() > 0.5 ? Int64 : Float64
           sin(rand(T)
       end
[...]
```
"""
descend_code_warntype(@nospecialize(args...); kwargs...) =
    _descend_with_error_handling(args...; iswarn=true, kwargs...)

function _descend_with_error_handling(@nospecialize(f), @nospecialize(argtypes = Tuple{}); kwargs...)
    ft = Core.Typeof(f)
    if isa(argtypes, Type)
        u = unwrap_unionall(argtypes)
        tt = rewrap_unionall(Tuple{ft, u.parameters...}, argtypes)
    else
        tt = Tuple{ft, argtypes...}
    end
    __descend_with_error_handling(tt; kwargs...)
end
_descend_with_error_handling(@nospecialize(tt::Type{<:Tuple}); kwargs...) =
    __descend_with_error_handling(tt; kwargs...)
_descend_with_error_handling(interp::AbstractInterpreter, mi::MethodInstance; kwargs...) =
    __descend_with_error_handling(interp, mi; kwargs...)
function __descend_with_error_handling(args...; terminal=default_terminal(), kwargs...)
    @nospecialize
    try
        _descend(terminal, args...; kwargs...)
    catch x
        if x isa InterruptException
            return nothing
        else
            rethrow(x)
        end
    end
    return nothing
end

default_terminal() = REPL.LineEdit.terminal(Base.active_repl)

"""
    descend

Shortcut for [`descend_code_typed`](@ref).
"""
const descend = descend_code_typed

function descend_code_typed(mi::MethodInstance; terminal=default_terminal(), kwargs...)
    interp = Cthulhu.CthulhuInterpreter()
    Cthulhu.do_typeinf!(interp, mi)
    _descend(terminal, interp, mi; iswarn=false, interruptexc=false, kwargs...)
end
function descend_code_warntype(mi::MethodInstance; terminal=default_terminal(), kwargs...)
    interp = Cthulhu.CthulhuInterpreter()
    Cthulhu.do_typeinf!(interp, mi)
    _descend(terminal, interp, mi; iswarn=true, interruptexc=false, kwargs...)
end

descend(interp::CthulhuInterpreter, mi::MethodInstance; kwargs...) = _descend(interp, mi; iswarn=false, interruptexc=false, kwargs...)

function codeinst_rt(code::CodeInstance)
    rettype = code.rettype
    if isdefined(code, :rettype_const)
        rettype_const = code.rettype_const
        if isa(rettype_const, Vector{Any}) && !(Vector{Any} <: rettype)
            return Core.PartialStruct(rettype, rettype_const)
        elseif isa(rettype_const, Core.PartialOpaque) && rettype <: Core.OpaqueClosure
            return rettype_const
        elseif isa(rettype_const, Core.InterConditional) && !(Core.InterConditional <: rettype)
            return rettype_const
        else
            return Const(rettype_const)
        end
    else
        return rettype
    end
end

# `@constprop :aggressive` here in order to make sure the constant propagation of `allow_no_src`
@constprop :aggressive function lookup(interp::CthulhuInterpreter, mi::MethodInstance, optimize::Bool; allow_no_src::Bool=false)
    if !optimize
        codeinf = src = copy(interp.unopt[mi].src)
        rt = interp.unopt[mi].rt
        infos = interp.unopt[mi].stmt_info
        slottypes = src.slottypes
        if isnothing(slottypes)
            slottypes = Any[ Any for i = 1:length(src.slotflags) ]
        end
    else
        codeinst = interp.opt[mi]
        rt = codeinst_rt(codeinst)
        opt = codeinst.inferred
        if opt !== nothing
            opt = opt::OptimizedSource
            src = Core.Compiler.copy(opt.ir)
            codeinf = opt.src
            infos = src.stmts.info
            slottypes = src.argtypes
        elseif allow_no_src
            # This doesn't showed up as covered, but it is (see the CI test with `coverage=false`).
            # But with coverage on, the empty function body isn't empty due to :code_coverage_effect expressions.
            codeinf = src = nothing
            infos = []
            slottypes = Any[Base.unwrap_unionall(mi.specTypes).parameters...]
        else
            Core.eval(Main, quote
                interp = $interp
                mi = $mi
                optimize = $optimize
            end)
            error("couldn't find the source; inspect `Main.interp` and `Main.mi`")
        end
    end
    # NOTE return `codeinf::CodeInfo` in any case since it can provide additional information on slot names
    (; src, rt, infos, slottypes, codeinf)
end

##
# _descend is the main driver function.
# src/reflection.jl has the tools to discover methods
# src/ui.jl provides the user facing interface to which _descend responds
##
function _descend(term::AbstractTerminal, interp::CthulhuInterpreter, mi::MethodInstance;
    override::Union{Nothing,InferenceResult}=nothing, debuginfo::Union{Symbol,DebugInfo}=DInfo.compact, # default is compact debuginfo
    optimize::Bool=true, interruptexc::Bool=true,
    iswarn::Bool=false, hide_type_stable::Union{Nothing,Bool}=nothing, verbose::Union{Nothing,Bool}=nothing,
    remarks::Bool=false, inline_cost::Bool=false, type_annotations::Bool=true)
    if isnothing(hide_type_stable)
        hide_type_stable = something(verbose, false)
    end
    isnothing(verbose) || Base.depwarn("The `verbose` keyword argument to `Cthulhu.descend` is deprecated. Use `hide_type_stable` instead.")
    if isa(debuginfo, Symbol)
        debuginfo = getfield(DInfo, debuginfo)::DebugInfo
    end

    is_cached(key, opt::Bool) = haskey(opt ? interp.opt : interp.unopt, key)

    menu_options = (cursor = '•', scroll_wrap = true)
    display_CI = true
    view_cmd = cthulhu_typed
    iostream = term.out_stream::IO
    while true
        if override === nothing && optimize && interp.opt[mi].inferred === nothing
            # This was inferred to a pure constant - we have no code to show,
            # but make something up that looks plausible.
            if display_CI
                println(iostream)
                println(iostream, "│ ─ $(string(Callsite(-1, MICallInfo(mi, interp.opt[mi].rettype), :invoke)))")
                println(iostream, "│  return ", Const(interp.opt[mi].rettype_const))
                println(iostream)
            end
            callsites = Callsite[]
        else
            if override !== nothing
                if optimize
                    opt = override.src
                    rt = override.result
                    if isa(opt, Compiler.OptimizationState)
                        src = Compiler.copy(src.ir::IRCode)
                        codeinf = src.src
                        infos = opt.stmts.info
                        slottypes = opt.argtypes
                    else # source might be unavailable at this point, when a result is fully constant-folded etc.
                        (; src, rt, infos, slottypes, codeinf) = lookup(interp, mi, optimize)
                    end
                else
                    codeinf = src = copy(interp.unopt[override].src)
                    rt = interp.unopt[override].rt
                    infos = interp.unopt[override].stmt_info
                    slottypes = src.slottypes
                    if isnothing(slottypes)
                        slottypes = Any[ Any for i = 1:length(src.slotflags) ]
                    end
                end
            else
                (; src, rt, infos, slottypes, codeinf) = lookup(interp, mi, optimize)
            end
            src = preprocess_ci!(src, mi, optimize, CONFIG)
            if optimize # optimization might have deleted some statements
                infos = src.stmts.info
            else
                @assert length(src.code) == length(infos)
            end
            callsites = find_callsites(interp, src, infos, mi, slottypes, optimize)

            if display_CI
                _remarks = remarks ? get(interp.remarks, mi, nothing) : nothing
                printstyled(IOContext(iostream, :limit=>true), mi.def, '\n'; bold=true)
                if debuginfo == DInfo.compact
                    str = let debuginfo=debuginfo, src=src, codeinf=codeinf, rt=rt, mi=mi,
                              iswarn=iswarn, hide_type_stable=hide_type_stable,
                              remarks=_remarks, inline_cost=inline_cost, type_annotations=type_annotations
                        stringify() do io # eliminate trailing indentation (see first item in bullet list in PR #189)
                            lambda_io = IOContext(io, :SOURCE_SLOTNAMES => Base.sourceinfo_slotnames(codeinf))
                            cthulhu_typed(lambda_io, debuginfo, src, rt, mi;
                                          iswarn, hide_type_stable,
                                          remarks, inline_cost, type_annotations,
                                          interp)
                        end
                    end
                    rmatch = findfirst(r"\u001B\[90m\u001B\[(\d+)G( *)\u001B\[1G\u001B\[39m\u001B\[90m( *)\u001B\[39m$", str)
                    if rmatch !== nothing
                        str = str[begin:prevind(str, first(rmatch))]
                    end
                    print(iostream, str)
                else
                    lambda_io = IOContext(iostream, :SOURCE_SLOTNAMES => Base.sourceinfo_slotnames(codeinf))
                    cthulhu_typed(lambda_io, debuginfo, src, rt, mi;
                                  iswarn, hide_type_stable,
                                  remarks=_remarks, inline_cost, type_annotations,
                                  interp)
                end
                view_cmd = cthulhu_typed
            end
            display_CI = true
        end

        menu = CthulhuMenu(callsites, optimize, iswarn&get(iostream, :color, false)::Bool; menu_options...)
        usg = usage(view_cmd, optimize, iswarn, hide_type_stable, debuginfo, remarks, inline_cost, type_annotations, CONFIG.enable_highlighter)
        cid = request(term, usg, menu)
        toggle = menu.toggle

        if toggle === nothing
            if cid == length(callsites) + 1
                break
            end
            if cid == -1
                interruptexc ? throw(InterruptException()) : break
            end
            callsite = callsites[cid]

            info = callsite.info
            if info isa MultiCallInfo
                sub_callsites = let callsite=callsite
                    map(ci->Callsite(callsite.id, ci, callsite.head), info.callinfos)
                end
                if isempty(sub_callsites)
                    @warn "Expected multiple callsites, but found none. Please fill an issue with a reproducing example" info
                    continue
                end
                menu = CthulhuMenu(sub_callsites, optimize, false; sub_menu=true, menu_options...)
                cid = request(term, "", menu)
                if cid == length(sub_callsites) + 1
                    continue
                end
                if cid == -1
                    interruptexc ? throw(InterruptException()) : break
                end

                callsite = sub_callsites[cid]
                info = callsite.info
            end

            # forcibly enter and inspect the frame, although the native interpreter gave up
            if info isa UncachedCallInfo || info isa TaskCallInfo
                # XXX: this may use a different native interpreter
                next_interp = CthulhuInterpreter()
                next_mi = get_mi(info)::MethodInstance
                do_typeinf!(next_interp, next_mi)
                _descend(term, next_interp, next_mi;
                         debuginfo,
                         optimize, interruptexc,
                         iswarn, hide_type_stable,
                         remarks, inline_cost, type_annotations)
                continue
            end

            if info isa GeneratedCallInfo || callsite.info isa FailedCallInfo
                @error "Callsite %$(callsite.id) failed to be extracted" callsite
            end

            # recurse
            next_mi = get_mi(callsite)::Union{MethodInstance,Nothing}
            if next_mi === nothing
                continue
            end

            _descend(term, interp, next_mi;
                     override = isa(info, ConstPropCallInfo) ? info.result : nothing, debuginfo,
                     optimize, interruptexc,
                     iswarn, hide_type_stable,
                     remarks, inline_cost, type_annotations)

        elseif toggle === :warn
            iswarn ⊻= true
        elseif toggle === :hide_type_stable
            hide_type_stable ⊻= true
        elseif toggle === :optimize
            optimize ⊻= true
            if !is_cached(mi, optimize)
                @warn "can't switch to post-optimization state, since this inference frame isn't cached"
                optimize ⊻= true
            end
        elseif toggle === :debuginfo
            debuginfo = DebugInfo((Int(debuginfo) + 1) % 3)
        elseif toggle === :remarks
            remarks ⊻= true
            if remarks && optimize
                @warn "disable optimization to see the inference remarks"
            end
        elseif toggle === :inline_cost
            inline_cost ⊻= true
            if inline_cost && !optimize
                @warn "enable optimization to see the inlining costs"
            end
        elseif toggle === :type_annotations
            type_annotations ⊻= true
        elseif toggle === :highlighter
            CONFIG.enable_highlighter ⊻= true
            if CONFIG.enable_highlighter
                @info "Using syntax highlighter $(CONFIG.highlighter)"
            else
                @info "Turned off syntax highlighter for Julia, LLVM and native code."
            end
            display_CI = false
        elseif toggle === :dump_params
            @info "Dumping inference cache"
            Core.show(mapany(((i, x),) -> (i, x.result, x.linfo), enumerate(get_inference_cache(interp))))
            Core.println()
            display_CI = false
        elseif toggle === :bookmark
            push!(BOOKMARKS, Bookmark(mi, interp))
            @info "The method is pushed at the end of `Cthulhu.BOOKMARKS`."
            display_CI = false
        elseif toggle === :revise
            # Call Revise.revise() without introducing a dependency on Revise
            id = Base.PkgId(UUID("295af30f-e4ad-537b-8983-00126c2a3abe"), "Revise")
            mod = get(Base.loaded_modules, id, nothing)
            if mod !== nothing
                revise = getfield(mod, :revise)::Function
                revise()
                mi = get_specialization(mi.specTypes)::MethodInstance
                do_typeinf!(interp, mi)
            end
        elseif toggle === :edit
            edit(whereis(mi.def::Method)...)
            display_CI = false
        else
            #Handle Standard alternative view, e.g. :native, :llvm
            if toggle === :typed
                view_cmd = cthulhu_typed
                display_CI = true
            else
                view_cmd = get(CODEVIEWS, toggle, nothing)
                @assert !isnothing(view_cmd) "invalid option $toggle"
                println(iostream)
                view_cmd(iostream, mi, optimize, debuginfo, interp, CONFIG)
                display_CI = false
            end
        end
        println(iostream)
    end
end
_descend(interp::CthulhuInterpreter, mi::MethodInstance; kwargs...) =
    _descend(default_terminal(), interp::CthulhuInterpreter, mi::MethodInstance; kwargs...)

function do_typeinf!(interp::CthulhuInterpreter, mi::MethodInstance)
    result = InferenceResult(mi)
    # we may want to handle the case when `InferenceState(...)` returns `nothing`,
    # which indicates code generation of a `@generated` has been failed,
    # and show it in the UI in some way ?
    # branch on https://github.com/JuliaLang/julia/pull/42082
    frame = @static hasmethod(InferenceState, (InferenceResult,Symbol,AbstractInterpreter)) ?
            InferenceState(result, #=cache=# :global, interp)::InferenceState :
            InferenceState(result, #=cached=# true, interp)::InferenceState
    Core.Compiler.typeinf(interp, frame)
    return nothing
end

function get_specialization(@nospecialize(TT))
    match = Base._which(TT)
    mi = specialize_method(match)
    return mi
end

function get_specialization(@nospecialize(F), @nospecialize(TT))
    return get_specialization(Base.signature_type(F, TT))
end

function mkinterp(interp::AbstractInterpreter, @nospecialize(args...))
    interp′ = CthulhuInterpreter(interp)
    mi = get_specialization(args...)
    do_typeinf!(interp′, mi)
    (interp′, mi)
end

mkinterp(@nospecialize(args...); interp::AbstractInterpreter=NativeInterpreter()) = mkinterp(interp, args...)

function _descend(@nospecialize(args...);
                  interp::AbstractInterpreter=NativeInterpreter(), kwargs...)
    (interp′, mi) = mkinterp(interp, args...)
    _descend(interp′, mi; kwargs...)
end
function _descend(term::AbstractTerminal, @nospecialize(args...);
                  interp::AbstractInterpreter=NativeInterpreter(), kwargs...)
    (interp′, mi) = mkinterp(interp, args...)
    _descend(term, interp′, mi; kwargs...)
end

descend_code_typed(b::Bookmark; kw...) =
    _descend_with_error_handling(b.interp, b.mi; iswarn=false, kw...)

descend_code_warntype(b::Bookmark; kw...) =
    _descend_with_error_handling(b.interp, b.mi; iswarn=true, kw...)

FoldingTrees.writeoption(buf::IO, data::Data, charsused::Int) = FoldingTrees.writeoption(buf, data.callstr, charsused)

function ascend(term, mi; interp::AbstractInterpreter=NativeInterpreter(), kwargs...)
    root = treelist(mi)
    menu = TreeMenu(root)
    choice = menu.current
    while choice !== nothing
        menu.chosen = false
        choice = TerminalMenus.request(term, "Choose a call for analysis (q to quit):", menu; cursor=menu.currentidx)
        browsecodetyped = true
        if choice !== nothing
            node = menu.current
            mi = instance(node.data.nd)
            if !isroot(node)
                # Help user find the sites calling the parent
                miparent = instance(node.parent.data.nd)
                ulocs = find_caller_of(interp, miparent, mi)
                if !isempty(ulocs)
                    strlocs = [string(" "^k[3] * '"', k[2], "\", ", k[1], ": lines ", v) for (k, v) in ulocs]
                    push!(strlocs, "Browse typed code")
                    linemenu = TerminalMenus.RadioMenu(strlocs; charset=:ascii)
                    browsecodetyped = false
                    choice2 = 1
                    while choice2 != -1
                        choice2 = TerminalMenus.request(term, "\nChoose caller of $miparent or proceed to typed code:", linemenu; cursor=choice2)
                        if 0 < choice2 < length(strlocs)
                            loc = ulocs[choice2]
                            edit(String(loc[1][2]), first(loc[2]))
                        elseif choice2 == length(strlocs)
                            browsecodetyped = true
                            break
                        end
                    end
                end
            end
            # The main application of `ascend` is finding cases of non-inferrability, so the
            # warn highlighting is useful.
            interp′ = CthulhuInterpreter(interp)
            do_typeinf!(interp′, mi)
            browsecodetyped && _descend(term, interp′, mi; iswarn=true, optimize=false, interruptexc=false, kwargs...)
        end
    end
end
ascend(mi; kwargs...) = ascend(default_terminal(), mi; kwargs...)

"""
    ascend(mi::MethodInstance; kwargs...)
    ascend(bt; kwargs...)

Follow a chain of calls (either through a backtrace `bt` or the backedges of a `MethodInstance` `mi`),
with the option to `descend` into intermediate calls. `kwargs` are passed to [`descend`](@ref).
"""
ascend

if ccall(:jl_generating_output, Cint, ()) == 1
    input = Pipe()
    Base.link_pipe!(input, reader_supports_async=true, writer_supports_async=true)
    term = REPL.Terminals.TTYTerminal("dumb", input.out, IOBuffer(), IOBuffer())
    write(input.in, 'q')
    descend(gcd, (Int, Int); terminal=term)
    nothing
end

end


text:
# Preprocessing script, written in julia6

# Some useful constants
const unkchar = '\x11'
# Train data:
const traindir = "/ai/data/nlp/semeval2018/t12/train-full.txt"

# Dev data:
const devdir = "/ai/data/nlp/semeval2018/t12/dev-full.txt"

type Instance
    id
    w0
    w1
    correctLabel
    reason
    claim
    title
    inform
end

function Instance(x, chvoc)
    @assert(length(x) == 8, "The datapoint seems broken $x")
    # The change in struct initialization happens here
    f(t) = word2Indx(t, chvoc)
    id = x[1]
    w0 = map(f, split(x[2]))
    w1 = map(f, split(x[3]))
    correctLabel = x[4]
    reason = map(f, split(x[5]))
    claim = map(f, split(x[6]))
    title = map(f, split(x[7]))
    inform = map(f, split(x[8]))
    return Instance(id, w0, w1, correctLabel, reason, claim, title, inform)
end

function readfile(tfile, chvoc)
    [ Instance(split(line, "\t"), chvoc) for line in eachline(tfile) ]
end

word2Indx(word, chvoc) = [ get(chvoc, item, chvoc[unkchar]) for item in word ]

# Task description:
# Argument == claim + reason
# debateTitle == Topic
# debateInfo  == Additional Information
# w0, w1

# The initially created language model can be found in:
# ../deeparse/data/english_chmodel.jld

max_stars_repo_name: kirnap/ArgReasoning

max_issues_repo_name: kirnap/ArgReasoning

max_forks_repo_name: kirnap/ArgReasoning


text:
module GLFWBackend

using Libdl
using GLFW
using ..LibCImGui
using CImGui: FLT_MAX, AddInputCharacter
using CImGui: GetIO, Set_KeyMap, Set_KeysDown, Get_KeysDown, Set_MouseDown, GetMouseCursor

include("callbacks.jl")
export SetCustomMouseButtonCallback
export SetCustomScrollCallback
export SetCustomKeyCallback
export SetCustomCharCallback

include("impl.jl")
export GlfwClientApi, GlfwClientApi_Unknown, GlfwClientApi_OpenGL, GlfwClientApi_Vulkan
export ImGui_ImplGlfw_Init
export ImGui_ImplGlfw_InitForOpenGL, ImGui_ImplGlfw_InitForVulkan
export ImGui_ImplGlfw_Shutdown
export ImGui_ImplGlfw_NewFrame

const g_Window = Ref{GLFW.Window}()
const g_ClientApi = Ref(GlfwClientApi_Unknown)
const g_Time = Ref(0.0)
const g_MouseJustPressed = [false, false, false, false, false]
const g_MouseCursors = [GLFW.Cursor(C_NULL) for i = 1:ImGuiMouseCursor_COUNT]
const g_ImplGlfw_GetClipboardText = Ref(C_NULL)
const g_ImplGlfw_SetClipboardText = Ref(C_NULL)
const g_CustomCallbackMousebutton = Ref{Any}(C_NULL)
const g_CustomCallbackScroll = Ref{Any}(C_NULL)
const g_CustomCallbackKey = Ref{Any}(C_NULL)
const g_CustomCallbackChar = Ref{Any}(C_NULL)

function __init__()
    g_ImplGlfw_GetClipboardText[] = dlsym(dlopen(GLFW.libglfw), :glfwGetClipboardString)
    g_ImplGlfw_SetClipboardText[] = dlsym(dlopen(GLFW.libglfw), :glfwSetClipboardString)
end

end # module