text:
@testset "ANOVA Tests" begin
    @testset "Multidimensional Array" begin
        @testset "1-way ANOVA" begin
            expected = [AnovaValue( "Total", 4822.4575, 19),
                        AnovaResult(    "A", 4684.9975,  3, 1561.66583, 181.773995, 1.43627554e-12),
                        AnovaFactor("Error",  137.46,   16,    8.59125)]

            @testset "Replicate Vectors" begin
                observations = Array{Vector{Float64}, 1}(undef, 4)
                observations[1] = [ 60.8,  57.0,  65.0, 58.6,  61.7]
                observations[2] = [ 68.7,  67.7,  74.9, 66.3,  69.8]
                observations[3] = [102.6, 102.1, 100.2, 96.5, 100.4]
                observations[4] = [ 87.9,  84.2,  83.1, 85.7,  90.3]

                results = anova(observations)

                @test all(expected .≈ results.effects)
            end

            @testset "Replicate First Dimension" begin
                observations = [60.8 68.7 102.6 87.9;
                                57.0 67.7 102.1 84.2;
                                65.0 74.9 100.2 83.1;
                                58.6 66.3  96.5 85.7;
                                61.7 69.8 100.4 90.3]

                results = anova(observations)

                @test all(expected .≈ results.effects)
            end
        end

        @testset "2-way ANOVA" begin
            observations1 = Array{Vector{Float64}, 2}(undef, 2, 2)
            observations1[1,1] = [16.5, 18.4, 12.7, 14.0, 12.8]
            observations1[1,2] = [14.5, 11.0, 10.8, 14.3, 10.0]
            observations1[2,1] = [39.1, 26.2, 21.3, 35.8, 40.2]
            observations1[2,2] = [32.0, 23.8, 28.8, 25.0, 29.3]

            observations2 = cat(hcat([16.5, 18.4, 12.7, 14.0, 12.8], [39.1, 26.2, 21.3, 35.8, 40.2]),
                                hcat([14.5, 11.0, 10.8, 14.3, 10.0], [32.0, 23.8, 28.8, 25.0, 29.3]), dims = 3)

            @testset "Fixed-effects (Model I)" begin
                expected = [AnovaValue( "Total", 1827.6975, 19),
                            AnovaResult(    "A",   70.3125,  1,   70.3125,  3.0706495,  0.098856175),
                            AnovaResult(    "B", 1386.1125,  1, 1386.1125, 60.533556,   7.9430782e-7),
                            AnovaResult("A × B",    4.9005,  1,    4.9005,  0.21401199, 0.64987001),
                            AnovaFactor("Error",  366.372,  16,   22.89825)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1)
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2)
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "Random-effects (Model II)" begin
                expected = [AnovaValue( "Total", 1827.6975, 19),
                            AnovaResult(    "A",   70.3125,  1,   70.3125,  14.3480257,  0.16431864),
                            AnovaResult(    "B", 1386.1125,  1, 1386.1125, 282.85124,    0.037808553),
                            AnovaResult("A × B",    4.9005,  1,    4.9005,   0.21401199, 0.64987001),
                            AnovaFactor("Error",  366.372,  16,   22.89825)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random, random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random, random])
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "Mixed-effects (Model III)" begin
                expected = [AnovaValue( "Total", 1827.6975, 19),
                            AnovaResult(    "A",   70.3125,  1,   70.3125, 14.3480257,  0.16431864),
                            AnovaResult(    "B", 1386.1125,  1, 1386.1125, 60.533556,   7.9430782e-7),
                            AnovaResult("A × B",    4.9005,  1,    4.9005,  0.21401199, 0.64987001),
                            AnovaFactor("Error",  366.372,  16,   22.89825)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random])
                    @test all(expected .≈ results.effects)
                end
            end
        end

        @testset "2-way ANOVA without replicates" begin
            observations1 = Array{Vector{Float64}, 2}(undef, 3, 4)
            observations1[1,1] = [123]
            observations1[1,2] = [138]
            observations1[1,3] = [110]
            observations1[1,4] = [151]
            observations1[2,1] = [145]
            observations1[2,2] = [165]
            observations1[2,3] = [140]
            observations1[2,4] = [167]
            observations1[3,1] = [156]
            observations1[3,2] = [176]
            observations1[3,3] = [185]
            observations1[3,4] = [175]

            observations2 = [123 138 110 151; 145 165 140 167; 156 176 185 175]

            expected = [AnovaValue(     "Total", 5594.9167,  11),
                        AnovaResult(        "A", 1116.91667,  3,  372.30556,  2.6316513, 0.144561215),
                        AnovaResult(        "B", 3629.1667,   2, 1814.58333, 12.8264284, 0.0068110648),
                        AnovaFactor("Remainder",  848.83333,  6,  141.472222)]

            @testset "Fixed-effects (Model I)" begin
                @testset "Replicate Vectors" begin
                    results = anova(observations1)
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, hasreplicates = false)
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "Random-effects (Model II)" begin
                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random, random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random, random], hasreplicates = false)
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "Mixed-effects (Model III)" begin
                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random], hasreplicates = false)
                    @test all(expected .≈ results.effects)
                end
            end
        end

        @testset "1-way ANOVA with 1 nested factor" begin
            expected = [AnovaValue( "Total", 71.666667, 11),
                        AnovaResult(    "A", 61.166667,  2, 30.583333, 61.166667,   0.0037032412),
                        AnovaResult(    "B",  1.5,       3,  0.5,       0.33333333, 0.80220227),
                        AnovaFactor("Error",  9.0,       6,  1.5)]

            @testset "Replicate Vectors" begin
                observations = Array{Vector{Float64}, 2}(undef, 2, 3)
                observations[1,1] = [102, 104]
                observations[1,2] = [108, 110]
                observations[1,3] = [104, 106]
                observations[2,1] = [103, 104]
                observations[2,2] = [109, 108]
                observations[2,3] = [105, 107]

                results = anova(observations, [nested])
                @test all(expected .≈ results.effects)
            end

            @testset "Replicate First Dimension" begin
                observations = cat(hcat([102, 104], [103, 104]),
                                   hcat([108, 110], [109, 108]),
                                   hcat([104, 106], [105, 107]), dims = 3)

                results = anova(observations, [nested])
                @test all(expected .≈ results.effects)
            end
        end

        @testset "1-way ANOVA with 2 nested factors" begin
            expected = [AnovaValue( "Total", 101.58753,   47),
                        AnovaResult(    "A",   7.5108875,  2, 3.7554438, 2.00918315, 0.21481422),
                        AnovaResult(    "B",   5.5708812,  3, 1.8569604, 0.99348408, 0.45717102),
                        AnovaResult(    "C",  11.2148375,  6, 1.8691396, 0.87059412, 0.52587983),
                        AnovaFactor("Error",  77.290925,  36, 2.14697014)]

            @testset "Replicate Vectors" begin
                observations = Array{Vector{Float64}, 3}(undef, 2, 2, 3)
                observations[1,1,1] = [3.17, 4.41, 1.81, 1.74]
                observations[2,1,1] = [2.81, 4.98, 2.62, 2.53]
                observations[1,2,1] = [3.0,  3.02, 4.73, 1.77]
                observations[2,2,1] = [4.13, 0.71, 3.18, 3.34]
                observations[1,1,2] = [2.42, 1.28, 1.4,  2.56]
                observations[2,1,2] = [1.26, 1.08, 1.42, 0.85]
                observations[1,2,2] = [0.36, 0.35, 2.64, 3.75]
                observations[2,2,2] = [3.86, 2.53, 3.97, 3.03]
                observations[1,1,3] = [5.24, 2.24, 0.18, 3.06]
                observations[2,1,3] = [4.04, 4.14, 0.33, 4.61]
                observations[1,2,3] = [6.06, 1.61, 2.25, 2.44]
                observations[2,2,3] = [0.02, 3.95, 0.87, 2.0]

                results = anova(observations, [nested, nested])
                @test all(expected .≈ results.effects)
            end

            @testset "Replicate First Dimension" begin
                observations = cat(cat(hcat([3.17, 4.41, 1.81, 1.74], [2.81, 4.98, 2.62, 2.53]),
                                       hcat([3.0,  3.02, 4.73, 1.77], [4.13, 0.71, 3.18, 3.34]), dims = 3),
                                   cat(hcat([2.42, 1.28, 1.4,  2.56], [1.26, 1.08, 1.42, 0.85]),
                                       hcat([0.36, 0.35, 2.64, 3.75], [3.86, 2.53, 3.97, 3.03]), dims = 3),
                                   cat(hcat([5.24, 2.24, 0.18, 3.06], [4.04, 4.14, 0.33, 4.61]),
                                       hcat([6.06, 1.61, 2.25, 2.44], [0.02, 3.95, 0.87, 2.0]), dims = 3), dims = 4)

                results = anova(observations, [nested, nested])
                @test all(expected .≈ results.effects)
            end
        end

        @testset "2-way ANOVA with 1 nested factor, no replicates" begin
            observations1 = Array{Vector{Float64}, 3}(undef, 5, 2, 2)
            observations1[1,1,1] = [16.5]
            observations1[2,1,1] = [18.4]
            observations1[3,1,1] = [12.7]
            observations1[4,1,1] = [14.0]
            observations1[5,1,1] = [12.8]
            observations1[1,2,1] = [39.1]
            observations1[2,2,1] = [26.2]
            observations1[3,2,1] = [21.3]
            observations1[4,2,1] = [35.8]
            observations1[5,2,1] = [40.2]
            observations1[1,1,2] = [14.5]
            observations1[2,1,2] = [11.0]
            observations1[3,1,2] = [10.8]
            observations1[4,1,2] = [14.3]
            observations1[5,1,2] = [10.0]
            observations1[1,2,2] = [32.0]
            observations1[2,2,2] = [23.8]
            observations1[3,2,2] = [28.8]
            observations1[4,2,2] = [25.0]
            observations1[5,2,2] = [29.3]

            observations2 = cat(hcat([16.5, 18.4, 12.7, 14.0, 12.8], [39.1, 26.2, 21.3, 35.8, 40.2]),
                                hcat([14.5, 11.0, 10.8, 14.3, 10.0], [32.0, 23.8, 28.8, 25.0, 29.3]), dims = 3)

            @testset "Fixed-effects (Model I)" begin
               expected = [AnovaValue( "Total", 1827.6975, 19),
                           AnovaResult(    "A",   70.3125,  1,   70.3125,  3.0706495,  0.098856175),
                           AnovaResult(    "B", 1386.1125,  1, 1386.1125, 60.533556,   7.9430782e-7),
                           AnovaResult("A × B",    4.9005,  1,    4.9005,  0.21401199, 0.64987001),
                           AnovaFactor(    "C",  366.372,  16,   22.89825)]

               @testset "Replicate Vectors" begin
                   results = anova(observations1, [nested])
                   @test all(expected .≈ results.effects)
               end

               @testset "Replicate First Dimension" begin
                   results = anova(observations2, [nested], hasreplicates = false)
                   @test all(expected .≈ results.effects)
               end
            end

            @testset "Random-effects (Model II)" begin
               expected = [AnovaValue( "Total", 1827.6975, 19),
                           AnovaResult(    "A",   70.3125,  1,   70.3125,  14.3480257,  0.16431864),
                           AnovaResult(    "B", 1386.1125,  1, 1386.1125, 282.85124,    0.037808553),
                           AnovaResult("A × B",    4.9005,  1,    4.9005,   0.21401199, 0.64987001),
                           AnovaFactor(    "C",  366.372,  16,   22.89825)]

               @testset "Replicate Vectors" begin
                   results = anova(observations1, [nested, random, random])
                   @test all(expected .≈ results.effects)
               end

               @testset "Replicate First Dimension" begin
                   results = anova(observations2, [nested, random, random], hasreplicates = false)
                   @test all(expected .≈ results.effects)
               end
            end

            @testset "Mixed-effects (Model III)" begin
               expected = [AnovaValue( "Total", 1827.6975, 19),
                           AnovaResult(    "A",   70.3125,  1,   70.3125, 14.3480257,  0.16431864),
                           AnovaResult(    "B", 1386.1125,  1, 1386.1125, 60.533556,   7.9430782e-7),
                           AnovaResult("A × B",    4.9005,  1,    4.9005,  0.21401199, 0.64987001),
                           AnovaFactor(    "C",  366.372,  16,   22.89825)]

               @testset "Replicate Vectors" begin
                   results = anova(observations1, [nested, random])
                   @test all(expected .≈ results.effects)
               end

               @testset "Replicate First Dimension" begin
                   results = anova(observations2, [nested, random], hasreplicates = false)
                   @test all(expected .≈ results.effects)
               end
            end
        end

        @testset "3-way ANOVA" begin
            observations1 = Array{Vector{Float64}, 3}(undef, 2, 3, 3)
            observations1[1,1,1] = [1.9, 1.8, 1.6, 1.4]
            observations1[1,1,2] = [2.1, 2.0, 1.8, 2.2]
            observations1[1,1,3] = [1.1, 1.2, 1.0, 1.4]
            observations1[1,2,1] = [2.3, 2.1, 2.0, 2.6]
            observations1[1,2,2] = [2.4, 2.6, 2.7, 2.3]
            observations1[1,2,3] = [2.0, 2.1, 1.9, 2.2]
            observations1[1,3,1] = [2.9, 2.8, 3.4, 3.2]
            observations1[1,3,2] = [3.6, 3.1, 3.4, 3.2]
            observations1[1,3,3] = [2.9, 2.8, 3.0, 3.1]
            observations1[2,1,1] = [1.8, 1.7, 1.4, 1.5]
            observations1[2,1,2] = [2.3, 2.0, 1.9, 1.7]
            observations1[2,1,3] = [1.4, 1.0, 1.3, 1.2]
            observations1[2,2,1] = [2.4, 2.7, 2.4, 2.6]
            observations1[2,2,2] = [2.0, 2.3, 2.1, 2.4]
            observations1[2,2,3] = [2.4, 2.6, 2.3, 2.2]
            observations1[2,3,1] = [3.0, 3.1, 3.0, 2.7]
            observations1[2,3,2] = [3.1, 3.0, 2.8, 3.2]
            observations1[2,3,3] = [3.2, 2.9, 2.8, 2.9]

            observations2 = cat(cat(hcat([1.9, 1.8, 1.6, 1.4], [1.8, 1.7, 1.4, 1.5]),
                                    hcat([2.3, 2.1, 2.0, 2.6], [2.4, 2.7, 2.4, 2.6]),
                                    hcat([2.9, 2.8, 3.4, 3.2], [3.0, 3.1, 3.0, 2.7]), dims = 3),
                                cat(hcat([2.1, 2.0, 1.8, 2.2], [2.3, 2.0, 1.9, 1.7]),
                                    hcat([2.4, 2.6, 2.7, 2.3], [2.0, 2.3, 2.1, 2.4]),
                                    hcat([3.6, 3.1, 3.4, 3.2], [3.1, 3.0, 2.8, 3.2]), dims = 3),
                                cat(hcat([1.1, 1.2, 1.0, 1.4], [1.4, 1.0, 1.3, 1.2]),
                                    hcat([2.0, 2.1, 1.9, 2.2], [2.4, 2.6, 2.3, 2.2]),
                                    hcat([2.9, 2.8, 3.0, 3.1], [3.2, 2.9, 2.8, 2.9]), dims = 3), dims = 4)

            @testset "All Fixed" begin
                expected = [AnovaValue(     "Total", 30.355,       71),
                            AnovaResult(        "A",  1.8175,       2,  0.90875,      24.475062,    2.71459995e-8),
                            AnovaResult(        "B", 24.655833,     2, 12.3279167,   332.02369,     4.5550981e-31),
                            AnovaResult(        "C",  0.008888889,  1,  0.008888889,   0.239401496, 0.62662035),
                            AnovaResult(    "A × B",  1.10166667,   4,  0.27541667,    7.4177057,   7.7516681e-5),
                            AnovaResult(    "A × C",  0.37027778,   2,  0.18513889,    4.9862843,   0.0102990988),
                            AnovaResult(    "B × C",  0.17527778,   2,  0.08763889,    2.3603491,   0.104056404),
                            AnovaResult("A × B × C",  0.220555556,  4,  0.055138889,   1.4850374,   0.21958095),
                            AnovaFactor(    "Error",  2.005,       54,  0.03712963)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1)
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2)
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "All Random" begin
                expected = [AnovaValue(     "Total", 30.355,       71),
                            AnovaResult(        "A",  1.8175,       2,  0.90875,      4.1754946,   0.179725175),
                            AnovaResult(        "B", 24.655833,     2, 12.3279167,   40.036536,    0.0022084069),
                            AnovaResult(        "C",  0.008888889,  1,  0.008888889,  0.021925317, 0.88874495),
                            AnovaResult(    "A × B",  1.10166667,   4,  0.27541667,   4.9949622,   0.074190827),
                            AnovaResult(    "A × C",  0.37027778,   2,  0.18513889,   3.3576826,   0.139349696),
                            AnovaResult(    "B × C",  0.17527778,   2,  0.08763889,   1.58942065,  0.31046402),
                            AnovaResult("A × B × C",  0.220555556,  4,  0.055138889,  1.4850374,   0.21958095),
                            AnovaFactor(    "Error",  2.005,       54,  0.03712963)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random, random, random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random, random, random])
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "1 random, 2 fixed" begin
                expected = [AnovaValue(     "Total", 30.355,       71),
                            AnovaResult(        "A",  1.8175,       2,  0.90875,      4.9084771,   0.16924835),
                            AnovaResult(        "B", 24.655833,     2, 12.3279167,   44.760968,    0.001829334),
                            AnovaResult(        "C",  0.008888889,  1,  0.008888889,  0.239401496, 0.62662035),
                            AnovaResult(    "A × B",  1.10166667,   4,  0.27541667,   4.9949622,   0.074190827),
                            AnovaResult(    "A × C",  0.37027778,   2,  0.18513889,   4.9862843,   0.0102990988),
                            AnovaResult(    "B × C",  0.17527778,   2,  0.08763889,   2.3603491,   0.104056404),
                            AnovaResult("A × B × C",  0.220555556,  4,  0.055138889,  1.4850374,   0.21958095),
                            AnovaFactor(    "Error",  2.005,       54,  0.03712963)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random])
                    @test all(expected .≈ results.effects)
                end
            end

            @testset "2 random, 1 fixed" begin
                expected = [AnovaValue(     "Total", 30.355,       71),
                            AnovaResult(        "A",  1.8175,       2,  0.90875,      4.1754946,    0.179725175),
                            AnovaResult(        "B", 24.655833,     2, 12.3279167,   44.760968,    0.001829334),
                            AnovaResult(        "C",  0.008888889,  1,  0.008888889,  0.032274332, 0.8661604),
                            AnovaResult(    "A × B",  1.10166667,   4,  0.27541667,   4.9949622,   0.074190827),
                            AnovaResult(    "A × C",  0.37027778,   2,  0.18513889,   3.3576826,   0.139349696),
                            AnovaResult(    "B × C",  0.17527778,   2,  0.08763889,   2.3603491,   0.104056404),
                            AnovaResult("A × B × C",  0.220555556,  4,  0.055138889,  1.4850374,   0.21958095),
                            AnovaFactor(    "Error",  2.005,       54,  0.03712963)]

                @testset "Replicate Vectors" begin
                    results = anova(observations1, [random, random])
                    @test all(expected .≈ results.effects)
                end

                @testset "Replicate First Dimension" begin
                    results = anova(observations2, [random, random])
                    @test all(expected .≈ results.effects)
                end
            end
        end
    end

    @testset "Vectors" begin
        observations = [1.9, 1.8, 1.6, 1.4, 2.1, 2.0, 1.8, 2.2, 1.1, 1.2, 1.0, 1.4,
                        2.3, 2.1, 2.0, 2.6, 2.4, 2.6, 2.7, 2.3, 2.0, 2.1, 1.9, 2.2,
                        2.9, 2.8, 3.4, 3.2, 3.6, 3.1, 3.4, 3.2, 2.9, 2.8, 3.0, 3.1,
                        1.8, 1.7, 1.4, 1.5, 2.3, 2.0, 1.9, 1.7, 1.4, 1.0, 1.3, 1.2,
                        2.4, 2.7, 2.4, 2.6, 2.0, 2.3, 2.1, 2.4, 2.4, 2.6, 2.3, 2.2,
                        3.0, 3.1, 3.0, 2.7, 3.1, 3.0, 2.8, 3.2, 3.2, 2.9, 2.8, 2.9]

        factorassignments = [[repeat([1], 36); repeat([2], 36)],
                             repeat([repeat([1], 12); repeat([2], 12); repeat([3], 12)], 2),
                             repeat([repeat([1], 4); repeat([2], 4); repeat([3], 4)], 6)]

        expected = [AnovaValue(     "Total", 30.355,       71),
                    AnovaResult(        "A",  1.8175,       2,  0.90875,      24.475062,    2.71459995e-8),
                    AnovaResult(        "B", 24.655833,     2, 12.3279167,   332.02369,     4.5550981e-31),
                    AnovaResult(        "C",  0.008888889,  1,  0.008888889,   0.239401496, 0.62662035),
                    AnovaResult(    "A × B",  1.10166667,   4,  0.27541667,    7.4177057,   7.7516681e-5),
                    AnovaResult(    "A × C",  0.37027778,   2,  0.18513889,    4.9862843,   0.0102990988),
                    AnovaResult(    "B × C",  0.17527778,   2,  0.08763889,    2.3603491,   0.104056404),
                    AnovaResult("A × B × C",  0.220555556,  4,  0.055138889,   1.4850374,   0.21958095),
                    AnovaFactor(    "Error",  2.005,       54,  0.03712963)]

        @testset "3-way ANOVA ordered" begin
            results = anova(observations, factorassignments)
            @test all(expected .≈ results.effects)
        end

        @testset "3-way ANOVA shuffled" begin
            randorder = [65,  6, 31, 70, 45,  2, 41, 59,  3, 39, 40, 47, 61, 55, 13, 68, 67, 25,
                         60, 66, 44, 33, 42, 52, 56, 15, 49, 43, 69, 10, 51, 58, 19, 22,  8, 37,
                         71, 23, 46, 24, 54, 53, 16, 18, 11, 32, 29, 35, 63, 34, 64, 62, 28, 72,
                         30, 38, 17, 21,  9,  5, 50, 12,  7, 36, 48, 27, 14, 20,  4,  1, 57, 26]
            shuffledobservations = observations[randorder]
            shuffledfactorassignments = [fa[randorder] for fa ∈ factorassignments]

            results = anova(shuffledobservations, shuffledfactorassignments)

            @test all(expected .≈ results.effects)
        end

        @testset "3-way ANOVA tolerates nonconsecutive factor levels" begin
            nonconsecutivefactorassignments = [fa .* 2 for fa ∈ factorassignments]
            results = anova(observations, nonconsecutivefactorassignments)
            @test all(expected .≈ results.effects)
        end

        @testset "3-way ANOVA tolerates noninteger factor levels" begin
            factorassignmentstrings = [repeat(["1"], 36); repeat(["2"], 36)]
            stringsinfactorassignments = [[factorassignmentstrings]; factorassignments[2:3]]
            results = anova(observations, stringsinfactorassignments)
            @test all(expected .≈ results.effects)
        end

        @testset "2-way ANOVA rejects nested factor with exclusive levels" begin
            factorlevels = factorassignments .|> unique .|> sort
            nfactorlevels = length.(factorlevels)
            nlevels = [4; nfactorlevels]
            nestedfactorassignment = (repeat(1:4, Int(length(observations) / 4)) .+
                                      sum([(factorassignments[i] .- 1) .* prod(nlevels[1:i]) for i ∈ 1:length(factorassignments)])) ./ 4 .|> ceil .|> Int
            nestedfactorassignments = [[nestedfactorassignment]; factorassignments[2:3]]

            @test_throws Exception anova(observations, nestedfactorassignments, [nested])
        end
    end

    @testset "DataFrame" begin
        observations = [1.9, 1.8, 1.6, 1.4, 2.1, 2.0, 1.8, 2.2, 1.1, 1.2, 1.0, 1.4,
                        2.3, 2.1, 2.0, 2.6, 2.4, 2.6, 2.7, 2.3, 2.0, 2.1, 1.9, 2.2,
                        2.9, 2.8, 3.4, 3.2, 3.6, 3.1, 3.4, 3.2, 2.9, 2.8, 3.0, 3.1,
                        1.8, 1.7, 1.4, 1.5, 2.3, 2.0, 1.9, 1.7, 1.4, 1.0, 1.3, 1.2,
                        2.4, 2.7, 2.4, 2.6, 2.0, 2.3, 2.1, 2.4, 2.4, 2.6, 2.3, 2.2,
                        3.0, 3.1, 3.0, 2.7, 3.1, 3.0, 2.8, 3.2, 3.2, 2.9, 2.8, 2.9]

        factorassignments = [[repeat([1], 36); repeat([2], 36)],
                             repeat([repeat([1], 12); repeat([2], 12); repeat([3], 12)], 2),
                             repeat([repeat([1], 4); repeat([2], 4); repeat([3], 4)], 6)]

        expected = [AnovaValue(                       "Total", 30.355,       71),
                    AnovaResult(                    "FactorA",  1.8175,       2,  0.90875,      24.475062,    2.71459995e-8),
                    AnovaResult(                    "FactorB", 24.655833,     2, 12.3279167,   332.02369,     4.5550981e-31),
                    AnovaResult(                    "FactorC",  0.008888889,  1,  0.008888889,   0.239401496, 0.62662035),
                    AnovaResult(          "FactorA × FactorB",  1.10166667,   4,  0.27541667,    7.4177057,   7.7516681e-5),
                    AnovaResult(          "FactorA × FactorC",  0.37027778,   2,  0.18513889,    4.9862843,   0.0102990988),
                    AnovaResult(          "FactorB × FactorC",  0.17527778,   2,  0.08763889,    2.3603491,   0.104056404),
                    AnovaResult("FactorA × FactorB × FactorC",  0.220555556,  4,  0.055138889,   1.4850374,   0.21958095),
                    AnovaFactor(                      "Error",  2.005,       54,  0.03712963)]

        @testset "3-way ANOVA shuffled" begin
            df = DataFrame([[observations]; factorassignments], [:observations; :FactorC; :FactorB; :FactorA])
            results = anova(df, :observations, [:FactorC; :FactorB; :FactorA])
            @test all(expected .≈ results.effects)
        end
    end
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/11/2014

module AbstractCASCoordImpl

export AbstractCASCoord

abstract AbstractCASCoord

end

module AbstractCASCoordInterfaces

export
    setRecord,
    getRecord,
    getAll

function setRecord()
end

function getRecord()
end

function getAll()
end

end

ext: jl


text:
#===============================================================================
                            Plotting callback
===============================================================================#
mutable struct PlottingCallback{T,K} <: eMCMC.Callback
    segm_len::Int64
    num_segm::Int64
    plot_type::K
    size::Tuple{Int64,Int64}
    p1::T
    p2::T
    p::T

    function PlottingCallback(;segm_len, num_segm, size=(1400,1000), first_updt_only=true)
        p1, p2 = plot([1,2],[2,3]), plot()
        p = plot(p1, p2, layout=(1,2), size=size)
        plot_type = ( first_updt_only ? Val(:first) : Val(:all) )
        new{typeof(p1),typeof(plot_type)}(
            segm_len,
            num_segm,
            plot_type,
            size,
            p1,
            p2,
            p,
        )
    end
end

function eMCMC.check_if_execute(callback::PlottingCallback, iter)
    (
        (iter.mcmciter-1) % callback.segm_len == 0 &&
        (iter.pidx == 1) &&
        iter.mcmciter > 1
    )
end

function eMCMC.init!(callback::PlottingCallback, ws::eMCMC.GlobalWorkspace)
    callback.p1 = plot(ws, callback.plot_type; indices=nothing)
    callback.p2 = plot(ws, Val(:ar); indices=nothing)
    callback.p = plot(callback.p1, callback.p2, layout=(1,2), size=callback.size)
    display(callback.p)
end

function reinit!(callback::PlottingCallback, ws::eMCMC.GlobalWorkspace, iter)
    callback.p1 = plot(
        ws,
        callback.plot_type;
        indices = 1:(iter.mcmciter-1)
        #yguide = [
        #    (i == length(ws.state)+1 ? "log-likelihood" : "θ$i")
        #    for i in 1:length(ws.state)+1
        #]
    )
    callback.p2 = plot(
        ws,
        Val(:ar);
        indices = 1:(iter.mcmciter-1),
        #yguide = [
        #    "update $i" for i in 1:length(ws.state_history[1])
        #]
    )
    callback.p = plot(
        callback.p1,
        callback.p2,
        layout = (1,2),
        size = callback.size,
    )
    display(callback.p)
end

function update!(callback::PlottingCallback, ws::eMCMC.GlobalWorkspace, iter)
    plot!(
        callback.p1,
        ws,
        callback.plot_type;
        indices = (iter.mcmciter-callback.segm_len):(iter.mcmciter-1)
    )
    plot!(
        callback.p2,
        ws,
        Val(:ar);
        indices = (iter.mcmciter-callback.segm_len):(iter.mcmciter-1)
    )
    display(callback.p)
end

"""
    execute!(callback::Callback, ws::GlobalWorkspace, iter)

Execute the `callback`.
"""
function eMCMC.execute!(callback::PlottingCallback, lws::eMCMC.LocalWorkspace, ws::eMCMC.GlobalWorkspace, iter)
    (
        div(iter.mcmciter-1, callback.segm_len) %
        callback.num_segm == 0
    ) && return reinit!(callback, ws, iter)
    update!(callback, ws, iter)
end


text:
########################################
## File Name: control_types.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: Control Type Definitions for SACBP
########################################

import Base.vec

abstract type Control end

struct PosControl{T<:Real} <: Control
    t::Float64
    vel::Vector{T}
    dim::Int64
    function PosControl{T}(t,vel,dim) where {T<:Real}
        if dim != length(vel)
            error(ArgumentError("Invalid control vector length."));
        else
            return new(t,vel,dim);
        end
    end
end
PosControl(t::Real,vel::Vector{T},dim::Int64=length(vel)) where {T<:Real} = PosControl{T}(t,vel,dim);
PosControl(t::Real,vel::T) where {T<:Real} = PosControl{T}(t,[vel],1);
vec(u::PosControl) = u.vel;

struct MControl2D{T<:Real} <: Control
    t::Float64
    fx::T  # Force x
    fy::T  # Force y
    tr::T  # Torque
end
function MControl2D(t::Real,params::Vector{T}) where {T<:Real}
    if length(params) != 3
        error(ArgumentError("Invalid parameter vector length."))
    else
        return MControl2D{T}(t,params[1],params[2],params[3]);
    end
end
vec(u::MControl2D) = [u.fx,u.fy,u.tr];


text:
@recipe function f(im::IntrinsicMetric{T}) where T
    2 ≤ ambientdim(im) ≤ 3 || error("Can't plot $(ambientdim(im))d points!")
    flat = ambientdim(im) == 2

    legend := false
    pts = points(im)
    g = adjgraph(im)

    @series begin
        seriestype := :scatter
        markersize --> 0.5
        label := :vertex
        if flat
            getindex.(pts, 1),
            getindex.(pts, 2)
        else
            getindex.(pts, 1),
            getindex.(pts, 2),
            getindex.(pts, 3)
        end
    end

    xs = T[]; ys = T[]; zs = T[]
    for e in edges(g)
        p1 = pts[src(e)]
        p2 = pts[dst(e)]
        append!(xs, [p1[1], p2[1], NaN])
        append!(ys, [p1[2], p2[2], NaN])
        !flat && append!(zs, [p1[3], p2[3], NaN])
    end

    @series begin
        markersize := :none
        seriestype := :path
        linewidth --> 1
        label := :edge
        if flat
            xs, ys
        else
            xs, ys, zs
        end
    end
end


text:
using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq
using Random
using StatsBase: sample

using Distributions: Exponential, Geometric, cdf

using Plots
using LinearAlgebra
using ProgressBars

include("utils.jl")

# age structure ------------------------------------------------------------
@present ThAgeSIR(FreeSchema) begin
    (S,I,R,Agent,Age)::Ob
    s::Hom(S, Agent)
    i::Hom(I, Agent)
    r::Hom(R, Agent)
    age::Hom(Age, Agent)

    AgeValue::AttrType
    agevalue::Attr(Age, AgeValue)
end

@acset_type AgeSIR(ThAgeSIR, index = [:s, :i, :r, :age])

getage(st::AgeSIR, x::Symbol) = st[vcat(incident(st, st[x], :age)...), :agevalue]

# infection rules
I = @acset AgeSIR{Int64} begin Agent=2; I=1; i=[2] end # need 2 agents otherwise have dangling edges
I2 = @acset AgeSIR{Int64} begin Agent=2; I=2; i=[1,2] end
SI = @acset AgeSIR{Int64} begin Agent=2; I=1; S=1; s=[1]; i=[2] end

L_infect = ACSetTransformation(I, SI; Agent = [1,2], I = [1])
R_infect = ACSetTransformation(I, I2; Agent = [1,2], I = [2])

# recovery rules
I1 = @acset AgeSIR{Int64} begin Agent=1; I=1; i=[1] end
R1 = @acset AgeSIR{Int64} begin Agent=1; R=1; r=[1] end
A1 = @acset AgeSIR{Int64} begin Agent=1 end

L_recovery = ACSetTransformation(A1, I1; Agent = [1])
R_recovery = ACSetTransformation(A1, R1; Agent = [1])

# contact matrix and parameters
N = 1000
I0 = 10
S0 = N - I0
Δt = 0.1
tmax = 100
steps = Int(tmax/Δt)
γ = 1/10 # recovery rate
R0 = 2.5

pop_TW_actual = [1011137, 1041749, 975803, 1213008, 1502279, 1618075, 1621625, 1926961, 2032452, 1755391, 1806638, 1852580, 1684376, 1416638, 903130, 1454933]
pop_TW = N * (pop_TW_actual / sum(pop_TW_actual))
pop_TW = Int.(round.(pop_TW))
N_ages = Float64.(pop_TW)

contact_unnorm = [1.13133995414878	0.557101193757009	0.277403588569526	0.166283484217312	0.243637116538954	0.501429856885401	0.773456505042611	0.664734678057495	0.373030323264845	0.207403718009079	0.219768245286902	0.183713318566528	0.100598824431586	0.0723635363224441	0.0451987854057666	0.0264190015627899
0.522116359073187	5.13993582380529	0.999881750990576	0.256872404884024	0.152965888814777	0.440037571740546	0.81718457664922	0.901816339582795	0.776770796127256	0.324978689513501	0.203946781656598	0.170308924697168	0.114410112258284	0.0766043987718478	0.0377257124824009	0.0309979523804595
0.233720556183555	1.92985936809443	11.0692392115508	0.983950046684279	0.331933259964229	0.36113787699594	0.563341231085396	0.903032808933502	1.17790738474796	0.651636521136264	0.342842258726613	0.173034004951767	0.085317189948624	0.0792891389004607	0.0627981424580935	0.0524757101169723
0.117682987086876	0.345910173819922	3.30103123371644	10.0791939893023	1.50916450884044	0.791285884658043	0.640497694236322	0.824953404104752	1.08583985409548	1.09515652926247	0.535792916941561	0.197364288066501	0.0713874830715285	0.0516466911367817	0.0315463617042823	0.0232625811949096
0.184283846973784	0.180619839511097	0.291480297355471	2.42727014100382	4.36737726042213	2.04934997728767	1.39723781506946	1.20032801881538	1.01070904475343	1.23524796447696	0.800832127889001	0.39902019518293	0.0944429030797382	0.0416382822716035	0.0470826453808244	0.0394951197797532
0.56156192817194	0.342992537769556	0.200122810772951	0.916989232128664	2.3980129020946	4.30977975618343	2.22681401915761	1.59261033191491	1.30238635661337	1.05207634128304	1.00949814315475	0.536928086612407	0.168547286689823	0.0624501078266029	0.0317257423796732	0.0269431049949397
0.817794867047499	1.1441520682733	0.841196707142586	0.551555851206977	1.22865834990464	2.17367993548303	3.72979423098863	2.14929466728635	1.55306730764128	1.11703404262428	0.859302034849221	0.601599257124208	0.256942068421412	0.106203837431624	0.0555526815130332	0.0511357782774989
0.665898107120277	1.15314199919716	0.9740043186808	0.697233734859189	0.825555277703876	1.54165452016297	1.99243882433986	3.25956559702894	2.06558139725526	1.22850364764608	0.791670106907652	0.438605150114473	0.244794735730732	0.148777023280774	0.0869686510606626	0.0349673148163165
0.344192143343282	0.807995410469339	1.14415043227526	1.22438695056131	0.988030619038175	1.30424464223099	1.72028688625099	1.87814450829015	2.90470388452872	1.49758830410338	0.934324050254837	0.330600903622231	0.182537686520201	0.114580140863156	0.0813846406012246	0.0400784886603155
0.235791530653026	0.585194161111923	0.776189668181761	1.67596626847397	0.98479644076933	1.01459833535067	1.25482711892082	1.36074838781177	1.44439037412333	2.06915205205562	0.978401494457824	0.39999986541145	0.147110687316517	0.0835852658511084	0.0775170890988546	0.0761516911425131
0.247247471139332	0.643310594857866	1.05811260803268	1.31212557944165	1.09555967723574	1.35073741525286	1.16660463514035	1.06882688561953	1.37118566862	1.49192984016939	1.73933382380217	0.740193371293642	0.243542236335681	0.109266882350377	0.0798121160310869	0.0833033422190585
0.442684358021657	0.702087727411373	0.662177052932204	0.816859526244507	0.738982978496004	1.2358668685183	1.27567127460055	0.896236373088294	0.950326986928672	0.797183567958301	1.02749599715298	1.51681655407295	0.451707898396268	0.195099577202936	0.0961333966552319	0.0769572529007728
0.327638342195265	0.328285016941782	0.238165901349415	0.354821937305308	0.320966945158153	0.529029787102833	0.677350792373831	0.650949037046765	0.501900491233428	0.407240319726415	0.395920768113556	0.541100083728198	0.836769891969799	0.267847588746602	0.145638827229407	0.059691205799689
0.193568331575816	0.313138348784034	0.2640962452464	0.160094097408356	0.204553272242237	0.325505767060564	0.516407578574901	0.473341415992634	0.438186035505819	0.281001067157269	0.296302451590512	0.375946282497057	0.331749859023468	0.679792406234649	0.173000119847818	0.0702778253123571
0.0963863876721027	0.282393513167451	0.295599962015119	0.30895124164141	0.1384465024105	0.257953175229946	0.27963826844241	0.436956466193058	0.509042820432503	0.403187215044375	0.316140357668137	0.277393058839257	0.385660305905714	0.370945733661385	0.606083091791185	0.185944796718132
0.188903590501605	0.259894019899239	0.37803587788717	0.31873977192663	0.125868047138266	0.160125548829692	0.269195296463481	0.305611946209338	0.348244256184042	0.390220623327373	0.393237304756271	0.240097073151863	0.144076910885202	0.19701354194507	0.184345002762545	0.325239463933643];

contact = (contact_unnorm + transpose(contact_unnorm) .* ( pop_TW * transpose(1 ./ pop_TW))) ./ 2

C = contact
for i = 1:16
    for j = 1:16
        C[i,j] = contact[i,j]*N_ages[i]/N_ages[j]
    end
end

β = R0 * γ / max(eigvals(C)...);

# Setup initial conditions
ages = vcat([fill(i, pop_TW[i]) for i in 1:16]...)
shuffle!(ages)

sir_index = vcat(fill.(["S", "I", "R"],[S0, I0, N-S0-I0])...)
shuffle!(sir_index)

state = AgeSIR{Int64}()
add_parts!(state, :S, sum(sir_index .== "S"))
add_parts!(state, :I, sum(sir_index .== "I"))
add_parts!(state, :R, sum(sir_index .== "R"))

add_parts!(state, :Agent, N)
add_parts!(state, :Age, N)

set_subpart!(state, 1:nparts(state, :S), :s, findall(sir_index .== "S"))
set_subpart!(state, 1:nparts(state, :I), :i, findall(sir_index .== "I"))
set_subpart!(state, 1:nparts(state, :R), :r, findall(sir_index .== "R"))

set_subpart!(state, 1:N, :age, 1:N)
set_subpart!(state, 1:N, :agevalue, ages)

# Simulation loop
out = fill(-1, (steps, 3))

for t = ProgressBar(1:steps)

    # possible events
    infections_m = homomorphisms(SI, state)
    recoveries_m = homomorphisms(I1, state)

    # queued updates objects
    infections = queued_updates(infections_m, L_infect, R_infect)
    recoveries = queued_updates(recoveries_m, L_recovery, R_recovery)

    # get infection firing probs
    age_i = [state[only(incident(state, state[only(collect(x[:S])), :s],:age)), :agevalue]  for x in infections.matches]
    age_j = [state[only(incident(state, state[only(collect(x[:I])), :i],:age)), :agevalue]  for x in infections.matches]
    N_j = [N_ages[j] for j in age_j]
    C_ij = [C[i,j] for (i,j) in zip(age_i, age_j)]
    # hazard of each individual S-I possible effective infective contact occuring
    r = β * C_ij .* (1 ./ N_j)

    sample_matches(infections, r, Δt)

    # sample recovery occurances
    sample_matches(recoveries, γ, Δt)

    global state = fire_events(state, [infections, recoveries])

    # write output
    out[t, :] = [nparts(state, x) for x in [:S, :I, :R]]
end

# plots
plot(
    (1:steps) * Δt,
    out,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)

ever_infected_ages = [getage(state, :i); getage(state, :r)]
ever_infected_ages = [sum(ever_infected_ages .== i) for i = 1:16]

bar(pop_TW, label = false, xlabel = "Age bins", ylabel = "Number")
bar!(ever_infected_ages, label = false)


text:
function massflow1(k,m2,T,lambda,Nf=4)
    (1/6)*Nf*lambda*ThresholdFunctions.Fb1(k, m2, T)
end

function massflow2(k,m2,T,lambda,Nf=4)
    (1/6)*(Nf+2)*lambda*ThresholdFunctions.Fb1(k, m2, T)
end

function lambdaflow1(k,m2,T,lambda)
    -(1/3) *lambda^2*ThresholdFunctions.Fb2(k, m2, T)
end

function lambdaflow2(k,m2,T,lambda,Nf=4)
    -(1/6) *lambda^2*(Nf+8)*ThresholdFunctions.Fb2(k, m2, T)
end

function Imlambdaflow1(p0, k, m2,lambda, T, Imlambda,Relambda, ϵ = 0.1)
    -(1 / 3) * (
        -Imlambda^2 * ThresholdFunctions.ImFb2(p0, k, m2, T, ϵ) +
        2 * ThresholdFunctions.ReFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        ThresholdFunctions.ImFb2(p0, k, m2, T, ϵ) * Relambda^2
    )
end

function Relambdaflow1(p0, k, m2, lambda,T, Imlambda,Relambda, ϵ = 0.1)
    -(1 / 3) * (
        -Imlambda^2 * ThresholdFunctions.ReFb2(p0, k, m2, T, ϵ) -
        2 * ThresholdFunctions.ImFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        ThresholdFunctions.ReFb2(p0, k, m2, T, ϵ) * Relambda^2
    )
end

function Imlambdaflow2(p0, k, m2,lambda, T, Imlambda,Relambda,ϵ = 0.1)
    -(2 / 3) * (
        -Imlambda^2 * ThresholdFunctions.ImFb2(p0, k, m2, T, ϵ) +
        2 * ThresholdFunctions.ReFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        ThresholdFunctions.ImFb2(p0, k, m2, T, ϵ) * Relambda^2
    )
end

function Relambdaflow2(p0, k, m2,lambda, T, Imlambda,Relambda,ϵ = 0.1,Nf=4)
    -(2 / 3) * (
        -Imlambda^2 * ThresholdFunctions.ReFb2(p0, k, m2, T, ϵ) -
        2 * ThresholdFunctions.ImFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        ThresholdFunctions.ReFb2(p0, k, m2, T, ϵ) * Relambda^2
    )-(1/6) *lambda^2*(Nf+4)*ThresholdFunctions.Fb2(k, m2, T)
end

function twopointflow(p0, k, m2, T,lambda)
    -1 / 48 * (k^4 * csch(sqrt(k^2 + m2) / (2 * T))^2 * (lambda(k, abs(sqrt(k^2 + m2) - p0)) + lambda(k, abs(sqrt(k^2 + m2) + p0))) * (sqrt(k^2 + m2) + T * sinh(sqrt(k^2 + m2) / T))) / ((k^2 + m2)^(3 / 2) * pi^2 * T)
end

function twopointflow2(p0, k, m2, T, lambda)
    -1 / 24 * (
        k^4 *
        csch(sqrt(k^2 + m2) / (2 * T))^2 *
        lambda(k) *
        (sqrt(k^2 + m2) + T * sinh(sqrt(k^2 + m2) / T))
    ) / ((k^2 + m2)^(3 / 2) * pi^2 * T)
end


text:
module TestSpecialArrays

include("preamble.jl")

function testfiles()
    x = String[]
    for (root, dirs, files) in walkdir(@__DIR__), file in files
        if match(r"^test_.*\.jl$", file) !== nothing
            push!(x, relpath(joinpath(root, file), @__DIR__))
        end
    end
    return x
end

@testset "SpecialArrays" begin
    @testset "$f" for f in testfiles()
        include(f)
    end
end

end # module


text:
import Multroot: convolute

conv(a, b) = convolute(a, b)
conv(a::Poly, b::Poly) = conv(p.a, b.a)


text:
using Base.Test
using QuantumOptics

@testset "timecorrelations" begin

ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

fockbasis = FockBasis(10)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H = Ha + Hc + Hint

Ja = embed(basis, 1, sqrt(γ)*sm)
Ja2 = embed(basis, 1, sqrt(0.5*γ)*sp)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J = [Ja, Ja2, Jc]

Ψ₀ = basisstate(spinbasis, 2) ⊗ fockstate(fockbasis, 5)
ρ₀ = dm(Ψ₀)

tspan = [0.:10:100.;]

op = embed(basis, 1, sqrt(γ)*sz)
exp_values = timecorrelations.correlation(tspan, ρ₀, H, J, dagger(op), op)

ρ₀ = dm(Ψ₀)

tout, exp_values2 = timecorrelations.correlation(ρ₀, H, J, dagger(op), op; tol=1e-5)

@test length(exp_values) == length(tspan)
@test length(exp_values2) == length(tout)
@test norm(exp_values[1]-exp_values2[1]) < 1e-15
@test norm(exp_values[end]-exp_values2[end]) < 1e-4

n = length(tspan)
omega_sample = mod(n, 2) == 0 ? [-n/2:n/2-1;] : [-(n-1)/2:(n-1)/2;]
omega_sample .*= 2pi/tspan[end]
omega, S = timecorrelations.spectrum(omega_sample, H, J, op; rho_ss=ρ₀)

omega2, S2 = timecorrelations.spectrum(H, J, op)
@test length(omega2) == length(S2)

omegaFFT, SFFT = timecorrelations.correlation2spectrum(tspan, exp_values)

@test omega_sample == omegaFFT && S == SFFT

@test_throws ArgumentError timecorrelations.correlation2spectrum(tout, exp_values)

tspan[5] = tspan[4]
@test_throws ArgumentError timecorrelations.correlation2spectrum(tspan, exp_values)

end # testset


text:
include("vpg.jl")
include("A2C.jl")
include("ppo.jl")
include("A2CGAE.jl")
include("MAC.jl")
include("ddpg.jl")
include("td3.jl")
include("sac.jl")


text:
export RobustLoss,
       HuberRho, Huber, AndrewsRho, Andrews,
       BisquareRho, Bisquare, LogisticRho, Logistic,
       FairRho, Fair, TalwarRho, Talwar, QuantileRho, Quantile

abstract type RobustRho end

abstract type RobustRho1P{δ} <: RobustRho end # one parameter

struct RobustLoss{ρ} <: AtomicLoss where ρ <: RobustRho
   rho::ρ
end

(rl::RobustLoss)(x::AVR, y::AVR) = rl(x .- y)
(rl::RobustLoss)(r::AVR) = rl.rho(r)

# ψ(r) = ρ'(r)     (first derivative)
# ω(r) = ψ(r)/r    (weighing function) a threshold can be passed to clip weights
# ϕ(r) = ψ'(r)     (second derivative)

"""
$TYPEDEF

Huber weighing of the residuals corresponding to

``ρ(z) = z²/2``  if `|z|≤δ` and `ρ(z)=δ(|z|-δ/2)` otherwise.
"""
struct HuberRho{δ} <: RobustRho1P{δ}
   HuberRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end
Huber(δ::Real=1.0; delta::Real=δ) = HuberRho(delta)

(::HuberRho{δ})(r::AVR) where δ = begin
   ar = abs.(r)
   w  = ar .<= δ
   return sum( r.^2/2 .* w .+ δ .* (ar .- δ/2) .* .!w )
end

ψ(::Type{HuberRho{δ}}   ) where δ = (r, w) -> r * w + δ * sign(r) * (1.0 - w)
ω(::Type{HuberRho{δ}}, _) where δ = (r, w) -> w + (δ / abs(r)) * (1.0 - w)
ϕ(::Type{HuberRho{δ}}   ) where δ = (r, w) -> w

"""
$TYPEDEF

Andrews weighing of the residuals corresponding to

``ρ(z) = -cos(πz/δ)/(π/δ)²`` if `|z|≤δ` and `ρ(δ)` otherwise.
"""
struct AndrewsRho{δ} <: RobustRho1P{δ}
   AndrewsRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end
Andrews(δ::Real=1.0; delta::Real=δ) = AndrewsRho(delta)

(::AndrewsRho{δ})(r::AVR) where δ = begin
   ar = abs.(r)
   w  = ar .<= δ
   c  = π/δ
   κ  = (δ/π)^2
   return sum( -cos.(c .* r) .* κ .* w .+ κ .* .!w )
end

# Note, sinc(x) = sin(πx)/πx, well defined everywhere

ψ(::Type{AndrewsRho{δ}}   ) where δ = (r, w) -> (c = π/δ; w * sin(c * r) / c)
ω(::Type{AndrewsRho{δ}}, _) where δ = (r, w) -> w * sinc(r / δ)
ϕ(::Type{AndrewsRho{δ}}   ) where δ = (r, w) -> (cr = (π/δ) * r; w * cos(cr))

"""
$TYPEDEF

Bisquare weighing of the residuals corresponding to

``ρ(z) = δ²/6 (1-(1-(z/δ)²)³)`` if `|z|≤δ` and `δ²/6` otherwise.
"""
struct BisquareRho{δ} <: RobustRho1P{δ}
   BisquareRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end
Bisquare(δ::Real=1.0; delta::Real=δ) = BisquareRho(delta)

(::BisquareRho{δ})(r::AVR) where δ = begin
   ar = abs.(r)
   w  = ar .<= δ
   κ  = δ^2/6
   return sum( κ * (1.0 .- (1.0 .- (r ./ δ).^2).^3) .* w + κ .* .!w )
end

ψ(::Type{BisquareRho{δ}}   ) where δ = (r, w) -> w * r * (1.0 - (r / δ)^2)^2
ω(::Type{BisquareRho{δ}}, _) where δ = (r, w) -> w * (1.0 - (r / δ)^2)^2
ϕ(::Type{BisquareRho{δ}}   ) where δ = (r, w) -> (sr = r / δ; w * (1.0 + 5sr^4 - 6sr^2))

"""
$TYPEDEF

Logistic weighing of the residuals corresponding to

``ρ(z) = δ² log(cosh(z/δ))``
"""
struct LogisticRho{δ} <: RobustRho1P{δ}
   LogisticRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end
Logistic(δ::Real=1.0; delta::Real=δ) = LogisticRho(delta)

(::LogisticRho{δ})(r::AVR) where δ = begin
   return sum( δ^2 .* log.(cosh.(r ./ δ)) )
end

# similar to sinc, to avoid NaNs if tanh(0)/0 (lim is 1.0)
tanhc(x::Real) = ifelse(iszero(x), one(x), tanh(x)/x)

ψ(::Type{LogisticRho{δ}}   ) where δ = (r, _) -> δ * tanh(r / δ)
ω(::Type{LogisticRho{δ}}, _) where δ = (r, _) -> tanhc(r / δ)
ϕ(::Type{LogisticRho{δ}}   ) where δ = (r, w) -> sech(r / δ)^2

"""
$TYPEDEF

Fair weighing of the residuals corresponding to

``ρ(z) = δ² (|z|/δ - log(1+|z|/δ))``
"""
struct FairRho{δ} <: RobustRho1P{δ}
   FairRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end
Fair(δ::Real=1.0; delta::Real=δ) = FairRho(delta)

(::FairRho{δ})(r::AVR) where δ = begin
   sr = abs.(r) ./ δ
   return sum( δ^2 .* (sr .- log1p.(sr)) )
end

ψ(::Type{FairRho{δ}}   ) where δ = (r, _) -> δ * r / (abs(r) + δ)
ω(::Type{FairRho{δ}}, _) where δ = (r, _) -> δ / (abs.(r) + δ)
ϕ(::Type{FairRho{δ}}   ) where δ = (r, _) -> δ^2 / (abs.(r) + δ)^2

"""
$TYPEDEF

Talwar weighing of the residuals corresponding to

``ρ(z) = z²/2`` if `|z|≤δ` and `ρ(z)=ρ(δ)` otherwise.
"""
struct TalwarRho{δ} <: RobustRho1P{δ}
   TalwarRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end
Talwar(δ::Real=1.0; delta::Real=δ) = TalwarRho(delta)

(::TalwarRho{δ})(r::AVR) where δ = begin
   w = abs.(r) .<= δ
   return sum( r.^2 ./ 2 .* w .+ δ^2/2 .* .!w)
end

ψ(::Type{TalwarRho{δ}}   ) where δ = (r, w) -> w * r
ω(::Type{TalwarRho{δ}}, _) where δ = (_, w) -> w
ϕ(::Type{TalwarRho{δ}}   ) where δ = (_, w) -> w

"""
$TYPEDEF

Quantile regression weighing of the residuals corresponding to

``ρ(z) = z(δ - 1(z<0))``
"""
struct QuantileRho{δ} <: RobustRho1P{δ}
   QuantileRho(δ::Real=1.0; delta::Real=δ) = new{delta}()
end

Quantile(δ::Real=1.0; delta::Real=δ) = QuantileRho(delta)

(::QuantileRho{δ})(r::AVR) where δ = begin
   return sum( r .* (δ .- (r .<= 0.0)) )
end

ψ(::Type{QuantileRho{δ}}   ) where δ = (r, _) -> (δ - (r <= 0.0))
ω(::Type{QuantileRho{δ}}, τ) where δ = (r, _) -> (δ - (r <= 0.0)) / clip(r, τ)
ϕ(::Type{QuantileRho{δ}}   ) where δ = (_, _) -> error("Newton(CG) not available for Quantile Reg.")


text:
__precompile__()
module MathOptInterfaceXpress

export XpressOptimizer

using Xpress
const XPR = Xpress
using MathOptInterface
const MOI = MathOptInterface
using LinQuadOptInterface
const LQOI = LinQuadOptInterface

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear,
    LQOI.Quad
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.Linear, LQOI.IV),
    (LQOI.Quad, LQOI.EQ),
    (LQOI.Quad, LQOI.LE),
    (LQOI.Quad, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, LQOI.SOS1),
    (LQOI.VecVar, LQOI.SOS2),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct XpressOptimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    params::Dict{Any,Any}
    XpressOptimizer(::Void) = new()
end

LQOI.LinearQuadraticModel(::Type{XpressOptimizer}, env) = XPR.Model()

function XpressOptimizer(; kwargs...)

    env = nothing
    m = XpressOptimizer(nothing)
    m.params = Dict{Any,Any}()
    MOI.empty!(m)
    for (name,value) in kwargs
        m.params[name] = value
        XPR.setparam!(m.inner, XPR.XPRS_CONTROLS_DICT[name], value)
    end
    return m
end

function MOI.empty!(m::XpressOptimizer)
    MOI.empty!(m,nothing)
    for (name,value) in m.params
        XPR.setparam!(m.inner, XPR.XPRS_CONTROLS_DICT[name], value)
    end
end

LQOI.supported_constraints(s::XpressOptimizer) = SUPPORTED_CONSTRAINTS
LQOI.supported_objectives(s::XpressOptimizer) = SUPPORTED_OBJECTIVES

backend_type(m::XpressOptimizer, ::MOI.GreaterThan{T}) where T = Cchar('G')
backend_type(m::XpressOptimizer, ::MOI.LessThan{T}) where T    = Cchar('L')
backend_type(m::XpressOptimizer, ::MOI.EqualTo{T}) where T     = Cchar('E')
# Implemented separately
# backend_type(m::XpressOptimizer, ::MOI.Interval{T}) where T    = Cchar('R')

backend_type(m::XpressOptimizer, ::MOI.Zeros)        = Cchar('E')
backend_type(m::XpressOptimizer, ::MOI.Nonpositives) = Cchar('L')
backend_type(m::XpressOptimizer, ::MOI.Nonnegatives) = Cchar('G')

#=
not in LinQuad
=#

setparam!(instance::XpressOptimizer, name, val) = XPR.setparam!(instance.inner, XPR.XPRS_CONTROLS_DICT[name], val)

setlogfile!(instance::XpressOptimizer, path) = XPR.setlogfile(instance.inner, path::String)

cintvec(v::Vector) = convert(Vector{Int32}, v)

#=
    inner wrapper
=#

#=
Constraints
=#

MOI.canset(::XpressOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.canset(::XpressOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}) = true
MOI.canget(::XpressOptimizer, ::MOI.ConstraintSet, ::Type{LQOI.LCI{LQOI.IV}}) = true

LQOI.change_variable_bounds!(instance::XpressOptimizer, colvec, valvec, sensevec) = XPR.chgbounds!(instance.inner, cintvec(colvec), sensevec, valvec)

LQOI.get_variable_lowerbound(instance::XpressOptimizer, col) = XPR.get_lb(instance.inner, col, col)[1]
LQOI.get_variable_upperbound(instance::XpressOptimizer, col) = XPR.get_ub(instance.inner, col, col)[1]

LQOI.get_number_linear_constraints(instance::XpressOptimizer) = XPR.num_linconstrs(instance.inner)

LQOI.add_linear_constraints!(instance::XpressOptimizer, A::LQOI.CSRMatrix{Float64}, sensevec, rhsvec) = XPR.add_constrs!(instance.inner, A.row_pointers, A.columns, A.coefficients, sensevec, rhsvec)

function LQOI.add_ranged_constraints!(instance::XpressOptimizer, A::LQOI.CSRMatrix{Float64}, lowerbound, upperbound)
    newrows = length(lowerbound)
    rows = XPR.num_linconstrs(instance.inner)
    addedrows = collect((rows+1):(rows+newrows))

    sensevec = fill(Cchar('E'),newrows)
    XPR.add_constrs!(instance.inner, A.row_pointers, A.columns, A.coefficients, sensevec, upperbound)

    XPR.chg_rhsrange!(instance.inner, cintvec(addedrows), +upperbound-lowerbound)
end

function LQOI.modify_ranged_constraints!(instance::XpressOptimizer, rows::Vector{Int}, lowerbound::Vector{Float64}, upperbound::Vector{Float64})
    XPR.set_rhs!(instance.inner, rows, upperbound)
    XPR.chg_rhsrange!(instance.inner, cintvec(rows), +upperbound-lowerbound)
end

LQOI.get_rhs(instance::XpressOptimizer, row) = XPR.get_rhs(instance.inner, row, row)[1]

function LQOI.get_range(instance::XpressOptimizer, row)
    ub = XPR.get_rhs(instance.inner, row, row)[1]
    r = XPR.get_rhsrange(instance.inner, row, row)[1]
    return ub-r, ub
end

# TODO improve
function LQOI.get_linear_constraint(instance::XpressOptimizer, idx)
    A = XPR.get_rows(instance.inner, idx, idx)'
    return A.rowval-1, A.nzval
end

function LQOI.get_quadratic_constraint(instance::XpressOptimizer, idx)
    A = XPR.get_rows(instance.inner, idx, idx)'

    Q = XPR.get_qrowmatrix_triplets(instance.inner, idx)

    return A.rowval-1, A.nzval, Q
end

# notin LQOI
# TODO improve
function getcoef(instance::XpressOptimizer, row, col)
    A = XPR.get_rows(instance.inner, row, row)'
    cols = A.rowval
    vals = A.nzval

    pos = findfirst(cols, col)
    if pos > 0
        return vals[pos]
    else
        return 0.0
    end
end

LQOI.change_matrix_coefficient!(instance::XpressOptimizer, row, col, coef) = XPR.chg_coeffs!(instance.inner, row, col, coef)

LQOI.change_objective_coefficient!(instance::XpressOptimizer, col, coef) = XPR.set_objcoeffs!(instance.inner, Int32[col], Float64[coef])

LQOI.change_rhs_coefficient!(instance::XpressOptimizer, row, coef) = XPR.set_rhs!(instance.inner, Int32[row], Float64[coef])

LQOI.delete_linear_constraints!(instance::XpressOptimizer, rowbeg, rowend) = XPR.del_constrs!(instance.inner, cintvec(collect(rowbeg:rowend)))

LQOI.delete_quadratic_constraints!(instance::XpressOptimizer, rowbeg, rowend) = XPR.del_constrs!(instance.inner, cintvec(collect(rowbeg:rowend)))

LQOI.change_variable_types!(instance::XpressOptimizer, colvec, typevec) = XPR.chgcoltypes!(instance.inner, colvec, typevec)

LQOI.change_linear_constraint_sense!(instance::XpressOptimizer, rowvec, sensevec) = XPR.set_rowtype!(instance.inner, rowvec, sensevec)

LQOI.add_sos_constraint!(instance::XpressOptimizer, colvec, valvec, typ) = XPR.add_sos!(instance.inner, typ, colvec, valvec)

LQOI.delete_sos!(instance::XpressOptimizer, idx1, idx2) = XPR.del_sos!(instance.inner, cintvec(collect(idx1:idx2)))

# TODO improve getting processes
function LQOI.get_sos_constraint(instance::XpressOptimizer, idx)
    A, types = XPR.get_sos_matrix(instance.inner)
    line = A[idx,:] #sparse vec
    cols = line.nzind
    vals = line.nzval
    typ = types[idx] == Cchar('1') ? :SOS1 : :SOS2
    return cols, vals, typ
end

LQOI.get_number_quadratic_constraints(instance::XpressOptimizer) = XPR.num_qconstrs(instance.inner)

function scalediagonal!(V, I, J, scale)
    #  LQOI assumes 0.5 x' Q x, but Gurobi requires the list of terms, e.g.,
    #  2x^2 + xy + y^2, so we multiply the diagonal of V by 0.5. We don't
    #  multiply the off-diagonal terms since we assume they are symmetric and we
    #  only need to give one.
    #
    #  We also need to make sure that after adding the constraint we un-scale
    #  the vector because we can't modify user-data.
    for i in 1:length(I)
        if I[i] == J[i]
            V[i] *= scale
        end
    end
end

function LQOI.add_quadratic_constraint!(instance::XpressOptimizer, cols, coefs, rhs, sense, I, J, V)
    @assert length(I) == length(J) == length(V)
    scalediagonal!(V, I, J, 0.5)
    XPR.add_qconstr!(instance.inner, cols, coefs, I, J, V, sense, rhs)
    scalediagonal!(V, I, J, 2.0)
end
#=
    Objective
=#

function LQOI.set_quadratic_objective!(instance::XpressOptimizer, I, J, V)
    @assert length(I) == length(J) == length(V)
    XPR.delq!(instance.inner)
    # scalediagonal!(V, I, J, 0.5)
    XPR.add_qpterms!(instance.inner, I, J, V)
    # scalediagonal!(V, I, J, 2.0)
    return nothing
end

function LQOI.set_linear_objective!(instance::XpressOptimizer, colvec, coefvec)
    nvars = XPR.num_vars(instance.inner)
    obj = zeros(Float64, nvars)

    for i in eachindex(colvec)
        obj[colvec[i]] = coefvec[i]
    end

    XPR.set_obj!(instance.inner, obj)
    nothing
end

function LQOI.change_objective_sense!(instance::XpressOptimizer, symbol)
    if symbol == :min
        XPR.set_sense!(instance.inner, :minimize)
    else
        XPR.set_sense!(instance.inner, :maximize)
    end
end

function LQOI.get_linear_objective!(instance::XpressOptimizer, x::Vector{Float64})
    obj = XPR.get_obj(instance.inner)
    @assert length(x) == length(obj)
    for i in 1:length(obj)
        x[i] = obj[i]
    end
end

function LQOI.get_quadratic_terms_objective(instance::XpressOptimizer)
    I, J, V = XPR.getq(instance.inner)
    return I, J, V
end

function LQOI.get_objectivesense(instance::XpressOptimizer)
    s = XPR.model_sense(instance.inner)
    if s == :maximize
        return MOI.MaxSense
    else
        return MOI.MinSense
    end
end

#=
    Variables
=#

LQOI.get_number_variables(instance::XpressOptimizer) = XPR.num_vars(instance.inner)

LQOI.add_variables!(instance::XpressOptimizer, int) = XPR.add_cvars!(instance.inner, zeros(int))

LQOI.delete_variables!(instance::XpressOptimizer, col, col2) = XPR.del_vars!(instance.inner, col)

# LQOI.lqs_addmipstarts(m, colvec, valvec)
function LQOI.add_mip_starts!(instance::XpressOptimizer, colvec, valvec)
    x = zeros(XPR.num_vars(instance.inner))
    for i in eachindex(colvec)
        x[colvec[i]] = valvec[i]
    end
    XPR.loadbasis(instance.inner, x)
end

#=
    Solve
=#

LQOI.solve_mip_problem!(instance::XpressOptimizer) = XPR.mipoptimize(instance.inner)

LQOI.solve_quadratic_problem!(instance::XpressOptimizer) = ( writeproblem(instance, "db", "l");LQOI.solve_linear_problem!(instance) )

LQOI.solve_linear_problem!(instance::XpressOptimizer) = XPR.lpoptimize(instance.inner)

function LQOI.get_termination_status(instance::XpressOptimizer)
    stat_lp = XPR.get_lp_status2(instance.inner)
    if XPR.is_mip(instance.inner)
        stat_mip = XPR.get_mip_status2(instance.inner)
        if stat_mip == XPR.MIP_NotLoaded
            return MOI.OtherError
        elseif stat_mip == XPR.MIP_LPNotOptimal
            # MIP search incomplete but there is no linear sol
            # return MOI.OtherError
            return MOI.InfeasibleOrUnbounded
        elseif stat_mip == XPR.MIP_NoSolFound
            # MIP search incomplete but there is no integer sol
            other = xprsmoi_stopstatus(instance.inner)
            if other == MOI.OtherError
                return MOI.SlowProgress#OtherLimit
            else
                return other
            end

        elseif stat_mip == XPR.MIP_Solution
            # MIP search incomplete but there is a solution
            other = xprsmoi_stopstatus(instance.inner)
            if other == MOI.OtherError
                return MOI.OtherLimit
            else
                return other
            end

        elseif stat_mip == XPR.MIP_Infeasible
            if XPR.hasdualray(instance.inner)
                return MOI.Success
            else
                return MOI.InfeasibleNoResult
            end
        elseif stat_mip == XPR.MIP_Optimal
            return MOI.Success
        elseif stat_mip == XPR.MIP_Unbounded
            if XPR.hasprimalray(instance.inner)
                return MOI.Success
            else
                return MOI.UnboundedNoResult
            end
        end
        return MOI.OtherError
    else
        if stat_lp == XPR.LP_Unstarted
            return MOI.OtherError
        elseif stat_lp == XPR.LP_Optimal
            return MOI.Success
        elseif stat_lp == XPR.LP_Infeasible
            if XPR.hasdualray(instance.inner)
                return MOI.Success
            else
                return MOI.InfeasibleNoResult
            end
        elseif stat_lp == XPR.LP_CutOff
            return MOI.ObjectiveLimit
        elseif stat_lp == XPR.LP_Unfinished
            return xprsmoi_stopstatus(instance.inner)
        elseif stat_lp == XPR.LP_Unbounded
            if XPR.hasprimalray(instance.inner)
                return MOI.Success
            else
                return MOI.UnboundedNoResult
            end
        elseif stat_lp == XPR.LP_CutOffInDual
            return MOI.ObjectiveLimit
        elseif stat_lp == XPR.LP_Unsolved
            return MOI.OtherError
        elseif stat_lp == XPR.LP_NonConvex
            return MOI.InvalidModel
        end
        return MOI.OtherError
    end
end

function xprsmoi_stopstatus(instance::XpressOptimizer)
    ss = XPR.get_stopstatus(instance.inner)
    if ss == XPR.StopTimeLimit
        return MOI.TimeLimit
    elseif ss == XPR.StopControlC
        return MOI.Interrupted
    elseif ss == XPR.StopNodeLimit
        # should not be here
        warn("should not be here")
        return MOI.NodeLimit
    elseif ss == XPR.StopIterLimit
        return MOI.IterationLimit
    elseif ss == XPR.StopMIPGap
        return MOI.ObjectiveLimit
    elseif ss == XPR.StopSolLimit
        return MOI.SolutionLimit
    elseif ss == XPR.StopUser
        return MOI.Interrupted
    end
    return MOI.OtherError
end

function LQOI.get_primal_status(instance::XpressOptimizer)
    if XPR.is_mip(instance.inner)
        stat_mip = XPR.get_mip_status2(instance.inner)
        if stat_mip in [XPR.MIP_Solution, XPR.MIP_Optimal]
            return MOI.FeasiblePoint
        elseif XPR.MIP_Infeasible && XPR.hasdualray(instance.inner)
            return MOI.InfeasibilityCertificate
        elseif XPR.MIP_Unbounded && XPR.hasprimalray(instance.inner)
            return MOI.InfeasibilityCertificate
        elseif stat_mip in [XPR.MIP_LPOptimal, XPR.MIP_NoSolFound]
            return MOI.InfeasiblePoint
        end
        return MOI.UnknownResultStatus
    else
        stat_lp = XPR.get_lp_status2(instance.inner)
        if stat_lp == XPR.LP_Optimal
            return MOI.FeasiblePoint
        elseif stat_lp == XPR.LP_Unbounded && XPR.hasprimalray(instance.inner)
            return MOI.InfeasibilityCertificate
        # elseif stat_lp == LP_Infeasible
        #     return MOI.InfeasiblePoint - xpress wont return
        # elseif cutoff//cutoffindual ???
        else
            return MOI.UnknownResultStatus
        end
    end
end

function LQOI.get_dual_status(instance::XpressOptimizer)
    if XPR.is_mip(instance.inner)
        return MOI.UnknownResultStatus
    else
        stat_lp = XPR.get_lp_status2(instance.inner)
        if stat_lp == XPR.LP_Optimal
            return MOI.FeasiblePoint
        elseif stat_lp == XPR.LP_Infeasible && XPR.hasdualray(instance.inner)
            return MOI.InfeasibilityCertificate
        # elseif stat_lp == LP_Unbounded
        #     return MOI.InfeasiblePoint - xpress wont return
        # elseif cutoff//cutoffindual ???
        else
            return MOI.UnknownResultStatus
        end
    end
end

LQOI.get_variable_primal_solution!(instance::XpressOptimizer, place) = XPR.get_solution!(instance.inner, place)

function LQOI.get_linear_primal_solution!(instance::XpressOptimizer, place)
    if num_qconstrs(instance.inner) == 0
        XPR.get_slack_lin!(instance.inner, place)
        rhs = XPR.get_rhs(instance.inner)
        for i in eachindex(place)
            place[i] = -place[i]+rhs[i]
        end
    else
        XPR.get_slack_lin!(instance.inner, place)
        rhs = XPR.get_rhs(instance.inner)
        lrows = XPR.get_lrows(instance.inner)
        for (i,v) in enumerate(lrows)
            place[i] = -place[i]+rhs[v]
        end
    end
    nothing
end

function moi_lrows(m::Xpress.Model)
    tt_rows = collect(1:num_constrs(m))
    if num_qconstrs(m) == 0
        return tt_rows
    else
        return setdiff(tt_rows, XPR.get_qrows(m))
    end
end

function LQOI.get_quadratic_primal_solution!(instance::XpressOptimizer, place)
    if num_qconstrs(instance.inner) == 0
        return nothing
    else
        qrows = XPR.get_qrows(instance.inner)
        newplace = zeros(num_constrs(instance.inner))
        XPR.get_slack!(instance.inner, newplace)
        rhs = XPR.get_rhs(instance.inner)
        for (i,v) in enumerate(qrows)
            place[i] = -newplace[v]+rhs[v]
        end
    end
    nothing
end

LQOI.get_variable_dual_solution!(instance::XpressOptimizer, place) = XPR.get_reducedcost!(instance.inner, place)

LQOI.get_linear_dual_solution!(instance::XpressOptimizer, place) = XPR.get_dual_lin!(instance.inner, place)

function LQOI.get_quadratic_dual_solution!(instance::XpressOptimizer, place)
    if num_qconstrs(instance.inner) == 0
        return nothing
    else
        qrows = XPR.get_qrows(instance.inner)
        newplace = zeros(num_constrs(instance.inner))
        XPR.get_dual!(instance.inner, newplace)
        for (i,v) in enumerate(qrows)
            place[i] = newplace[v]
        end
    end
end

LQOI.get_objective_value(instance::XpressOptimizer) = XPR.get_objval(instance.inner)

LQOI.get_objective_bound(instance::XpressOptimizer) = XPR.get_mip_objval(instance.inner)

function LQOI.get_relative_mip_gap(instance::XpressOptimizer)
    L = XPR.get_mip_objval(instance.inner)
    U = XPR.get_bestbound(instance.inner)
    return abs(U-L)/U
end

LQOI.get_iteration_count(instance::XpressOptimizer)  = XPR.get_simplex_iter_count(instance.inner)

LQOI.get_barrier_iterations(instance::XpressOptimizer) = XPR.get_barrier_iter_count(instance.inner)

LQOI.get_node_count(instance::XpressOptimizer) = XPR.get_node_count(instance.inner)

LQOI.get_farkas_dual!(instance::XpressOptimizer, place) = XPR.getdualray!(instance.inner, place)

LQOI.get_unbounded_ray!(instance::XpressOptimizer, place) = XPR.getprimalray!(instance.inner, place)

MOI.free!(m::XpressOptimizer) = XPR.free_model(m.onner)

"""
    writeproblem(m: :MOI.AbstractOptimizer, filename::String)
Writes the current problem data to the given file.
Supported file types are solver-dependent.
"""
writeproblem(instance::XpressOptimizer, filename::String, flags::String="") = XPR.write_model(instance.inner, filename, flags)

end # module


text:
import brainflow

brainflow.enable_dev_board_logger()

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(Integer(brainflow.SYNTHETIC_BOARD))
sampling_rate = brainflow.get_sampling_rate(Integer(brainflow.SYNTHETIC_BOARD))

data_first_channel = data[eeg_channels[1], :]
println("Original Data First Channel")
println(data_first_channel)
brainflow.perform_rolling_filter(data_first_channel, 3, Integer(brainflow.MEAN))
println("After Rolling Filter")
println(data_first_channel)

data_second_channel = data[eeg_channels[2], :]
println("Original Data Second Channel")
println(data_second_channel)
brainflow.perform_wavelet_denoising(data_second_channel, "db4", 2)
println("After Wavelet Denoising")
println(data_second_channel)

max_stars_repo_name: kusumikakd/Brainflow

max_issues_repo_name: kusumikakd/Brainflow

max_forks_repo_name: kusumikakd/Brainflow


text:
# Metadata containing whether to use the feats-col or not

const use_feats = Dict{String, Bool}(
    "ar_padt"      =>false,
    "ru_syntagrus" =>true,
    "cs_pdt"       =>true,
    "en_ewt"       =>true,
    "tr_imst"      =>false,
    "sv_lines"     =>false,
    "en_gum"       =>false,
    "ca_ancora"    =>true,
    "hu_szeged"    =>false,
    "ro_rrt"       =>true,
    "sv_talbanken" =>true,
    "hi_hdtb"      =>true,
    "ko_gsd"       =>false,
    "it_isdt"      =>true,
    "pt_bosque"    =>true,
    "bg_btb"       =>true,
    "got_proiel"   =>true,
    "vi_vtb"       =>false,
    "da_ddt"       =>false,
    "eu_bdt"       =>true,
    "id_gsd"       =>true,
    "fi_ftb"       =>true,
    "fa_seraji"    =>true,
    "fi_tdt"       =>true,
    "cs_fictree"   =>true,
    "fro_srcmf"    =>true,
    "he_htb"       =>true,
    "ja_gsd"       =>true,
    "el_gdt"       =>false,
    "nl_alpino"    =>true,
    "fr_sequoia"   =>false,
    "en_lines"     =>true,
    "lv_lvtb"      =>true,
    "la_proiel"    =>true,
    "nl_lassysmall"=>true,
    "et_edt"       =>true,
    "cu_proiel"    =>true,
    "it_postwita"  =>false,
    "ko_kaist"     =>true,
    "ur_udtb"      =>true,
    "zh_gsd"       =>true,
    "pl_sz"        =>true,
    "fr_gsd"       =>true,
    "cs_cac"       =>true,
    "sk_snk"       =>true,
    "no_nynorsk"   =>true,
    "sl_ssj"       =>true,
    "fr_spoken"    =>false,
    "grc_perseus"  =>true,
    "af_afribooms" =>false,
    "uk_iu"        =>true,
    "ug_udt"       =>false,
    "pl_lfg"       =>true,
    "gl_ctg"       =>true,
    "bxr_bdt"      =>false,
    "hsb_ufal"     =>false,
    "hy_armtdp"    =>false,
    "kk_ktb"       =>false,
    "kmr_mg"       =>false,
    "es_ancora"    =>true,
    "ga_idt"       =>false,
    "la_ittb"      =>true,
    "sme_giella"   =>false,
    "no_nynorsklia"=>true,
    "ru_taiga"     =>true,
    "sl_sst"       =>true,
    "gl_treegal"   =>true,
    "de_gsd"       =>true,
    "grc_proiel"   =>true,
    "hr_set"       =>true,
    "no_bokmaal"   =>true,
    "sr_set"       =>false,
    "la_perseus"   =>true,
)


text:
__precompile__()

module Primes

import Base

const offwheel = Array{Bool}(210)
fill!(offwheel, false)
for p in [2, 3, 5, 7], i in 1:210÷p
    offwheel[i*p] = true
end

const incs = Int[]

function gen_incs()
    c=1
    for i=2:210
        if offwheel[i]
            c+=1
        else
            push!(incs, c)
            c=1
        end
    end
    push!(incs,c)
end

gen_incs()

struct Iter{P <: Integer}
    small_primes :: Array{P}
    sieve :: Dict{P, P}
    offwheel :: Array{Bool}

end

function Iter{P}() where P <: Integer
    first_primes = P[   2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31
                    ,  37,  41,  43,  47,  53,  59,  61,  67,  71 , 73,  79
                    ,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137
                    , 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193
                    , 197, 199, 211]

    sieve = Dict{P, P}()

    res = Iter{P}(first_primes, sieve, offwheel)
end

const IterState{P} = Tuple{P, Int}

function Base.start(it::Iter{P})::IterState{P} where P <: Integer
    (1, length(incs)-1)
end

function is_wheeled_off(n :: P) :: Bool where P <: Integer
    offwheel[(n-1) % 210 + 1]
end

function Base.next(it::Iter{P}, state::IterState{P}) where P <: Integer
    previous, inc_idx = state

    if isempty(it.small_primes)
        candidate = previous + incs[inc_idx]

        inc_idx = inc_idx == length(incs) ? 1 : inc_idx+1

        period = pop!(it.sieve, candidate, 0)

        while period != 0
            scratched = candidate + period
            while is_wheeled_off(scratched) || haskey(it.sieve, scratched)
                scratched += period
            end
            it.sieve[scratched] = period

            candidate += incs[inc_idx]
            inc_idx = inc_idx == length(incs) ? 1 : inc_idx+1

            period = get(it.sieve, candidate, 0)
        end

        it.sieve[candidate^2] = candidate
        (candidate, (candidate, inc_idx))
    else
        prime = shift!(it.small_primes)

        if prime > 7
            scratched = prime^2
            while scratched < 210 || is_wheeled_off(scratched)
                scratched += prime
            end
            it.sieve[scratched] = prime
        end

        (prime, (prime, 1))
    end
end

Base.done(::Iter{P}, ::IterState{P}) where P <: Integer = false

function test(limit)
    it=Iter{Int}()
    i=1

    for p in it
        # println("$(i): $(p)")
        i += 1
        if i > limit
            return p
        end
    end
end

function iter(T)
    Iter{T}()
end

end # Module

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
struct Empty end

struct TypedEmpty{T} end

@testset "Struct with no fields" begin
    @test begin
        b = JlrsReflect.reflect([Empty])
        sb = JlrsReflect.StringWrappers(b)

        sb[Empty] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.Empty", zero_sized_type)]
        pub struct Empty {
        }"""
    end
end

@testset "Struct with type parameter but no fields" begin
    @test begin
        b = JlrsReflect.reflect([TypedEmpty])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(TypedEmpty)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.TypedEmpty")]
        pub struct TypedEmpty {
        }"""
    end
end


text:
#
# Autogenerated by Thrift Compiler (0.11.0)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

# struct _enum_TDeviceType
#   CPU::Int32
#   GPU::Int32
# end
# const TDeviceType = _enum_TDeviceType(Int32(0), Int32(1))

# struct _enum_TDatumType
#   SMALLINT::Int32
#   INT::Int32
#   BIGINT::Int32
#   FLOAT::Int32
#   DECIMAL::Int32
#   DOUBLE::Int32
#   STR::Int32
#   TIME::Int32
#   TIMESTAMP::Int32
#   DATE::Int32
#   BOOL::Int32
#   INTERVAL_DAY_TIME::Int32
#   INTERVAL_YEAR_MONTH::Int32
#   POINT::Int32
#   LINESTRING::Int32
#   POLYGON::Int32
#   MULTIPOLYGON::Int32
#   TINYINT::Int32
#   GEOMETRY::Int32
#   GEOGRAPHY::Int32
# end
# const TDatumType = _enum_TDatumType(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7), Int32(8), Int32(9), Int32(10), Int32(11), Int32(12), Int32(13), Int32(14), Int32(15), Int32(16), Int32(17), Int32(18), Int32(19))

# struct _enum_TEncodingType
#   NONE::Int32
#   FIXED::Int32
#   RL::Int32
#   DIFF::Int32
#   DICT::Int32
#   SPARSE::Int32
#   GEOINT::Int32
#   DATE_IN_DAYS::Int32
# end
# const TEncodingType = _enum_TEncodingType(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

mutable struct TTypeInfo <: Thrift.TMsg
  _type::Int32
  encoding::Int32
  nullable::Bool
  is_array::Bool
  precision::Int32
  scale::Int32
  comp_param::Int32
  size::Int32
  TTypeInfo() = (o=new(); fillunset(o); o)
end # mutable struct TTypeInfo
meta(t::Type{TTypeInfo}) = meta(t, Symbol[:size], Int[1,4,2,3,5,6,7,8], Dict{Symbol,Any}(:size => Int32(-1)))


text:
using Pkg

Pkg.activate(".")
include("scrape_utils.jl")

docs = extract_docs()

build_reference(docs)
build_methods(docs)


text:
module BamReader

using DataStructures
using Libz
using BufferedStreams

import Base: show, print, seek, +, isless, isempty, copy, iterate, eltype, size, length,
             getindex, checkbounds, string, view, position, eof, read, readbytes!

export
	BamFile,
	BamRead,
	CIGAR,
	#CIGARState,
	CIGARElement,
	Seq,
	#SeqState,
	SeqElement,
	Quality,
	QualityState,
	QualityElement,
	filename,
	header,
	sequences,
	sequence_name,

	refID,
	pos,
	bin_mq_nl,
	flag_nc,
	l_seq,
	next_refID,
	next_pos,
	tlen,

	# bin, # TODO: find a name that doesn't clash with Base.bin
	mapq,
	l_read_name,
	flag,
	n_cigar_op,
	flagset,

	ref_name,
	next_ref_name,

	read_name,
	cigar,
	seq,
	qual,

	op,
	op_char,
	len,
	cigar_str,
	alignedlength,
	endpos,
	base_mask,
	base_index,
	base_char,
	seq_str,
	# quality_array,
	@tag_str,
	hastag,
	tagint,
	tagchar,
	tagfloat,
	tagstring,

	sam,

	BAM_FPAIRED,
	BAM_FPROPER_PAIR,
	BAM_FUNMAP,
	BAM_FMUNMAP,
	BAM_FREVERSE,
	BAM_FMREVERSE,
	BAM_FREAD1,
	BAM_FREAD2,
	BAM_FSECONDARY,
	BAM_FQCFAIL,
	BAM_FDUP,
	BAM_FSUPPLEMENTARY,
	CIGAR_M,
	CIGAR_I,
	CIGAR_D,
	CIGAR_N,
	CIGAR_S,
	CIGAR_H,
	CIGAR_P,
	CIGAR_Eq,
	CIGAR_X,

	readpos,
	refpos,

	BamIndex,
	reg2bin,
	reg2bins,
	reg2binranges,

	region

include("constants.jl")
include("BamIndex.jl")
include("BamFile.jl")
include("BamRead.jl")
include("BamState.jl")
include("BGZFBlockStream.jl")
include("BamRegion.jl")
include("readiterators.jl")

end


text:
findMinSVSV(a::SVector{3,Float64}, b::SVector{3,Float64}) = min.(a, b)
findMinSVSV(sv::SVector{3, SVector{3,Float64}}) = min.(min.(sv[1], sv[2]),      sv[3])
findMinSVSV(sv::SVector{4, SVector{3,Float64}}) = min.(min.(sv[1], sv[2]), min.(sv[3], sv[4]))

findMaxSVSV(a::SVector{3,Float64}, b::SVector{3,Float64}) = max.(a, b)
findMaxSVSV(sv::SVector{3, SVector{3,Float64}}) = max.(max.(sv[1], sv[2]),      sv[3])
findMaxSVSV(sv::SVector{4, SVector{3,Float64}}) = max.(max.(sv[1], sv[2]), max.(sv[3], sv[4]))

function minMaxToCenterExtent(box_min::SVector{3,Float64}, box_max::SVector{3,Float64})
    center = (box_max + box_min) * 0.5
    extent = (box_max - box_min) * 0.5
    return center, extent
end

### calc_min_max
function calc_min_max(a::OBB)
    δ = abs.(a.R) * a.e
    return a.c - δ, a.c + δ
end
function calc_min_max(point::Vector{SVector{3,T}}) where {T}
    min_val = SVector{3,Float64}(+Inf, +Inf, +Inf)
    max_val = SVector{3,Float64}(-Inf, -Inf, -Inf)
    for k = 1:length(point)
        point_k = point[k]
        min_val = min.(min_val, point_k)
        max_val = max.(max_val, point_k)
    end
    return min_val, max_val
end
function calc_min_max(a::SVector{3,Float64}, b::SVector{3,Float64})
    min_ = findMinSVSV(a, b)
    max_ = findMaxSVSV(a, b)
    return min_, max_
end

### calc_obb
function calc_obb(vert::SVector{N, SVector{3,Float64}}) where {N}
    box_min = findMinSVSV(vert)
    box_max = findMaxSVSV(vert)
    return calc_obb(box_min, box_max)
end
function calc_obb(arg_in)
    min_val, max_val = calc_min_max(arg_in)
    return calc_obb(min_val, max_val)
end
function calc_obb(a::SVector{3,Float64}, b::SVector{3,Float64})
    min_val, max_val = calc_min_max(a, b)
    center, extent = minMaxToCenterExtent(min_val, max_val)
    return OBB(center, extent, one(SMatrix{3,3,Float64,9}))
end

function sort_so_big_ϵ_last(ϵ::SVector{4,Float64}, thing::SVector{4,T}) where {T}
	ϵ = abs.(ϵ)
	v, i = findmax(ϵ)
	return thing[tet_perm_by_num(i)]
end


text:
"""
Generic implementation of an unconstrained single-field FE space
Private fields and type parameters
"""
struct UnconstrainedFESpace{V} <: SingleFieldFESpace
  vector_type::Type{V}
  nfree::Int
  ndirichlet::Int
  cell_dofs_ids::AbstractArray
  cell_shapefuns::CellField
  cell_dof_basis::CellDof
  dirichlet_dof_tag::Vector{Int8}
  dirichlet_cells::Vector{Int32}
  ntags::Int
end

# FESpace interface

ConstraintStyle(::Type{<:UnconstrainedFESpace}) = UnConstrained()
num_free_dofs(f::UnconstrainedFESpace) = f.nfree
zero_free_values(f::UnconstrainedFESpace) = allocate_vector(f.vector_type,num_free_dofs(f))
get_cell_shapefuns(f::UnconstrainedFESpace) = f.cell_shapefuns
get_cell_dof_basis(f::UnconstrainedFESpace) = f.cell_dof_basis
get_cell_dof_ids(f::UnconstrainedFESpace) = f.cell_dofs_ids
get_triangulation(f::UnconstrainedFESpace) = get_triangulation(f.cell_shapefuns)
get_dof_value_type(f::UnconstrainedFESpace{V}) where V = eltype(V)
get_vector_type(f::UnconstrainedFESpace{V}) where V = V

# SingleFieldFESpace interface

num_dirichlet_dofs(f::UnconstrainedFESpace) = f.ndirichlet
num_dirichlet_tags(f::UnconstrainedFESpace) = f.ntags
zero_dirichlet_values(f::UnconstrainedFESpace) = allocate_vector(f.vector_type,num_dirichlet_dofs(f))
get_dirichlet_dof_tag(f::UnconstrainedFESpace) = f.dirichlet_dof_tag

function scatter_free_and_dirichlet_values(f::UnconstrainedFESpace,free_values,dirichlet_values)
  @check eltype(free_values) == eltype(dirichlet_values) """\n
  The entries stored in free_values and dirichlet_values should be of the same type.

  This error shows up e.g. when trying to build a FEFunction from a vector of integers
  if the Dirichlet values of the underlying space are of type Float64, or when the
  given free values are Float64 and the Dirichlet values ComplexF64.
  """
  cell_dof_ids = get_cell_dof_ids(f)
  lazy_map(Broadcasting(PosNegReindex(free_values,dirichlet_values)),cell_dof_ids)
end

function gather_free_and_dirichlet_values!(free_vals,dirichlet_vals,f::UnconstrainedFESpace,cell_vals)

  cell_dofs = get_cell_dof_ids(f)
  cache_vals = array_cache(cell_vals)
  cache_dofs = array_cache(cell_dofs)
  cells = 1:length(cell_vals)

  _free_and_dirichlet_values_fill!(
    free_vals,
    dirichlet_vals,
    cache_vals,
    cache_dofs,
    cell_vals,
    cell_dofs,
    cells)

  (free_vals,dirichlet_vals)
end

function gather_dirichlet_values!(dirichlet_vals,f::UnconstrainedFESpace,cell_vals)

  cell_dofs = get_cell_dof_ids(f)
  cache_vals = array_cache(cell_vals)
  cache_dofs = array_cache(cell_dofs)
  free_vals = zero_free_values(f)
  cells = f.dirichlet_cells

  _free_and_dirichlet_values_fill!(
    free_vals,
    dirichlet_vals,
    cache_vals,
    cache_dofs,
    cell_vals,
    cell_dofs,
    cells)

  dirichlet_vals
end

function  _free_and_dirichlet_values_fill!(
  free_vals,
  dirichlet_vals,
  cache_vals,
  cache_dofs,
  cell_vals,
  cell_dofs,
  cells)

  for cell in cells
    vals = getindex!(cache_vals,cell_vals,cell)
    dofs = getindex!(cache_dofs,cell_dofs,cell)
    for (i,dof) in enumerate(dofs)
      val = vals[i]
      if dof > 0
        free_vals[dof] = val
      elseif dof < 0
        dirichlet_vals[-dof] = val
      else
        @unreachable "dof ids either positibe or negative, not zero"
      end
    end
  end

end


text:
# # Custom ReLU layer

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/custom-relu.jl)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/custom-relu.ipynb)

# We demonstrate how DiffOpt can be used to generate a simple neural network
# unit - the ReLU layer. A neural network is created using Flux.jl which is
# trained on the MNIST dataset.

# This tutorial uses the following packages

using JuMP
import DiffOpt
import OSQP
import ChainRulesCore
import Flux
import Statistics
import Base.Iterators: repeated

# ## The ReLU and its derivative

# Define a relu through an optimization problem solved by a quadratic solver.
# Return the solution of the problem.
function matrix_relu(
    y::AbstractArray{T};
    model = Model(() -> DiffOpt.diff_optimizer(OSQP.Optimizer))
) where T
    _x = zeros(size(y))
    N = length(y[:, 1])
    empty!(model)
    set_silent(model)
    @variable(model, x[1:N] >= 0)
    for i in 1:size(y)[2]
        @objective(
            model,
            Min,
            x'x -2x'y[:, i]
        )
        optimize!(model)
        _x[:, i] = value.(x)
    end
    return _x
end

# Define the backward differentiation rule, for the function we defined above.
function ChainRulesCore.rrule(
    ::typeof(matrix_relu),
    y::AbstractArray{T};
    model = Model(() -> DiffOpt.diff_optimizer(OSQP.Optimizer))
) where T
    pv = matrix_relu(y, model = model)
    function pullback_matrix_relu(dx)
        x = model[:x]
        dy = zeros(T, size(dx))
        for i in 1:size(y)[2]
            MOI.set.(
                model,
                DiffOpt.BackwardInVariablePrimal(),
                x,
                dx[:, i]
            ) # set sensitivities
            DiffOpt.backward(model) # compute grad
            obj_exp = MOI.get(
                model,
                DiffOpt.BackwardOutObjective()
            ) # return grdiente wrt objective function parameters
            dy[:, i] = JuMP.coefficient.(obj_exp, x) # coeff of `x` in -2x'y
            dy[:, i] = -2 * dy[:, i]
        end
        return (ChainRulesCore.NoTangent(), dy,)
    end
    return pv, pullback_matrix_relu
end

# ## prepare data
imgs = Flux.Data.MNIST.images()
labels = Flux.Data.MNIST.labels();

# Preprocessing
X = hcat(float.(reshape.(imgs, :))...) #stack all the images
Y = Flux.onehotbatch(labels, 0:9); # just a common way to encode categorical variables

N = 1000

train_X = X[:, 1:N]
train_Y = Y[:, 1:N]

test_X = hcat(float.(reshape.(Flux.Data.MNIST.images(:test), :))...)
test_Y = Flux.onehotbatch(Flux.Data.MNIST.labels(:test), 0:9);

# ## Define the Network

# Network structure

inner = 15

m = Flux.Chain(
    Flux.Dense(784, inner), #784 being image linear dimension (28 x 28)
    matrix_relu,
    Flux.Dense(inner, 10), # 10 beinf the number of outcomes (0 to 9)
    Flux.softmax,
)

# Define input data
# The original data is repeated `epochs` times because `Flux.train!` only
# loops through the data set once

epochs = 5

dataset = repeated((train_X, train_Y), epochs);

# Parameters for the network training

custom_loss(x, y) = Flux.crossentropy(m(x), y) # training loss function
opt = Flux.ADAM(); # stochastic gradient descent variant to optimize weights of the neral network
evalcb = () -> @show(custom_loss(X, Y)); # callback to show loss

# Train to optimize network parameters

@time Flux.train!(custom_loss, Flux.params(m), dataset, opt, cb = Flux.throttle(evalcb, 5));

# Although our custom implementation takes time, it is able to reach similar
# accuracy as the usual ReLU function implementation.

accuracy(x, y) = Statistics.mean(Flux.onecold(m(x)) .== Flux.onecold(y)); # average of correct guesses

# Train accuracy

accuracy(train_X, train_Y)

# Test accuracy

accuracy(test_X, test_Y)

# Note that the accuracy is low due to simplified training.
# It is possible to increase the number of samples `N`,
# the number of epochs `epoch` and teh conectivity `inner`.


text:
# ---
# title: 199. Binary Tree Right Side View
# id: problem199
# author: Indigo
# date: 2021-01-20
# difficulty: Medium
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/binary-tree-right-side-view/description/>
# hidden: true
# ---
#
# Given a binary tree, imagine yourself standing on the _right_ side of it,
# return the values of the nodes you can see ordered from top to bottom.
#
# **Example:**
#
#
#
#     Input:  [1,2,3,null,5,null,4]
#     Output:  [1, 3, 4]
#     Explanation:
#        1            <---
#      /   \
#     2     3         <---
#      \     \
#       5     4       <---
#
#
#
## @lc code=start
using LeetCode

function right_side_view(root::TreeNode{Int})::Vector{Int}
    q = Queue{Pair{TreeNode{Int}, Int}}()
    res = Int[]
    enqueue!(q, Pair(root, 1))
    while !isempty(q)
        nd, layer = dequeue!(q)
        (nd.left !== nothing) && enqueue!(q, Pair(nd.left, layer + 1))
        (nd.right !== nothing) && enqueue!(q, Pair(nd.right, layer + 1))
        (isempty(q) || layer < first(q).second) && (push!(res, nd.val))
    end
    res
end
## @lc code=end


text:
for N in [Float64, Rational{Int}, Float32]
    # =======
    #  Ball1
    # =======

    # approximation of a centered unit Ball1
    b = Ball1(N[0, 0, 0], N(1))
    @test norm(b) ≈ N(1)  # in the infinity norm (default)
    @test radius(b) ≈ N(1)  # in the infinity norm (default)
    @test diameter(b) ≈ N(2)  # in the infinity norm (default)

    # approximations of a non-centered unit Ball1
    b = Ball1(N[1, 2, 0], N(1))
    @test norm(b) ≈ N(3)  # in the infinity norm (default)
    @test radius(b) ≈ N(1)  # in the infinity norm (default)
    @test diameter(b) ≈ N(2)  # in the infinity norm (default)

    # ================
    #  Hyperrectangle
    # ================

    # metrics in the infinity norm (default)
    h = Hyperrectangle(N[-1, 2], N[0.2, 0.5])
    @test norm(h, Inf) ≈ N(2.5)
    @test radius(h, Inf) ≈ N(0.5)
    @test diameter(h, Inf) ≈ N(1)

    # metrics in the 2-norm
    @test norm(h, 2) ≈ N(sqrt(2.5^2 + 1.2^2))
    @test radius(h, 2) ≈ N(sqrt(0.5^2 + 0.2^2))
    @test diameter(h, 2) ≈ N(2*sqrt(0.5^2 + 0.2^2))

    # =========
    #  BallInf
    # =========

    # metrics in the infinity norm (default)
    b = BallInf(N[-1, -2, -4], N(0.2))
    @test norm(b, Inf) ≈ N(4.2)
    @test radius(b, Inf) ≈ N(0.2)
    @test diameter(b, Inf) ≈ N(0.4)

    # metrics in the 2-norm
    @test norm(b, 2) ≈ N(sqrt(1.2^2 + 2.2^2 + 4.2^2))
    @test radius(b, 2) ≈ N(sqrt(0.2^2 * 3))
    @test diameter(b, 2) ≈ N(2*sqrt(0.2^2 * 3))

    # ====================================
    #  failing case (not implemented yet)
    # ====================================

    s = MinkowskiSum(b, b)
    @test_throws ErrorException norm(s, 2)
    @test_throws ErrorException radius(s, 2)
    @test_throws ErrorException diameter(s, 2)
end