text:
using EvalMetrics.Encodings
import EvalMetrics.Encodings: positive_label, negative_label

@testset "ispositive" begin
    @test ispositive(OneZero(), 1)
    @test ispositive(OneTwo(), 1)
    @test ispositive(OneMinusOne(), 1)
    @test ispositive(OneVsRest(1, [2,3]), 1)
    @test ispositive(RestVsOne([1,2], 3), 1)
    @test ispositive(RestVsOne([1,2], 3), 2)
end

@testset "ispositive broadcasting" begin
    @test all(ispositive.(OneZero(), [1]))
    @test all(ispositive.(OneTwo(), [1]))
    @test all(ispositive.(OneMinusOne(), [1]))
    @test all(ispositive.(OneVsRest(1, [2,3]), [1]))
    @test all(ispositive.(RestVsOne([1,2], 3), [1]))
    @test all(ispositive.(RestVsOne([1,2], 3), [2]))
end

@testset "isnegative" begin
    @test isnegative(OneZero(), 0)
    @test isnegative(OneTwo(), 2)
    @test isnegative(OneMinusOne(), -1)
    @test isnegative(OneVsRest(1, [2,3]), 2)
    @test isnegative(OneVsRest(1, [2,3]), 3)
    @test isnegative(RestVsOne([1,2], 3), 3)
end

@testset "isnegative broadcasting" begin
    @test all(isnegative.(OneZero(), [0]))
    @test all(isnegative.(OneTwo(), [2]))
    @test all(isnegative.(OneMinusOne(), [-1]))
    @test all(isnegative.(OneVsRest(1, [2,3]), [2]))
    @test all(isnegative.(OneVsRest(1, [2,3]), [3]))
    @test all(isnegative.(RestVsOne([1,2], 3), [3]))
end

@testset "positive_label" begin
    @test positive_label(OneZero()) == 1
    @test positive_label(OneTwo()) == 1
    @test positive_label(OneMinusOne()) == 1
    @test positive_label(OneVsRest(1, [2,3])) == 1
    @test positive_label(RestVsOne([1,2], 3)) == 1
    @test positive_label(RestVsOne([2,1], 3)) == 2
end

@testset "negative_label" begin
    @test negative_label(OneZero()) == 0
    @test negative_label(OneTwo()) == 2
    @test negative_label(OneMinusOne()) == -1
    @test negative_label(OneVsRest(1, [2,3])) == 2
    @test negative_label(OneVsRest(1, [3,2])) == 3
    @test negative_label(RestVsOne([1,2], 3)) == 3
end

@testset "check_encoding" begin
    @test check_encoding(OneZero(), rand([0,1], 1000))
    @test !check_encoding(OneZero(), rand([0,1,2], 1000))

    @test check_encoding(OneTwo(), rand([1,2], 1000))
    @test !check_encoding(OneTwo(), rand([1,2,3], 1000))

    @test check_encoding(OneMinusOne(), rand([1,-1], 1000))
    @test !check_encoding(OneMinusOne(), rand([1,-1, 2], 1000))

    @test check_encoding(OneVsOne(3, 4), rand([3,4], 1000))
    @test !check_encoding(OneVsOne(3, 4), rand([1,2], 1000))
    @test check_encoding(OneVsOne(:three, :four), rand([:three, :four], 1000))
    @test !check_encoding(OneVsOne(:three, :four), rand([:one, :two], 1000))
    @test check_encoding(OneVsOne("three", "four"), rand(["three", "four"], 1000))
    @test !check_encoding(OneVsOne("three", "four"), rand(["one", "two"], 1000))

    @test check_encoding(OneVsRest(1, [2,3]), rand([1,2,3], 1000))
    @test !check_encoding(OneVsRest(1, [2,3]), rand([1,2,3,4], 1000))
    @test check_encoding(OneVsRest(:one, [:two, :three]), rand([:one, :two, :three], 1000))
    @test !check_encoding(OneVsRest(:one, [:two, :three]), rand([:four, :five], 1000))
    @test check_encoding(OneVsRest("one", ["two", "three"]), rand(["one", "two", "three"], 1000))
    @test !check_encoding(OneVsRest("one", ["two", "three"]), rand(["four", "five"], 1000))

    @test check_encoding(RestVsOne([1,2], 3), rand([1,2,3], 1000))
    @test !check_encoding(RestVsOne([1,2], 3), rand([1,2,3,4], 1000))
    @test check_encoding(RestVsOne([:one, :two], :three), rand([:one, :two, :three], 1000))
    @test !check_encoding(RestVsOne([:one, :two], :three), rand([:four, :five], 1000))
    @test check_encoding(RestVsOne(["one", "two"], "three"), rand(["one", "two", "three"], 1000))
    @test !check_encoding(RestVsOne(["one", "two"], "three"), rand(["four", "five"], 1000))
end

@testset "recode" begin
    @test recode.(OneZero(), OneMinusOne(), [0,1,0,0,1]) == [-1,1,-1,-1,1]
    @test recode.(OneZero(), OneTwo(), [0,1,0,0,1]) == [2,1,2,2,1]
    @test recode.(OneZero(), OneVsRest(2, [3,4,5]), [0,1,0,0,1]) == [3,2,3,3,2]
    @test recode.(OneZero(), OneVsRest(:two, [:three,:four,:five]), [0,1,0,0,1]) == [:three,:two,:three,:three,:two]
    @test recode.(OneZero(), OneVsRest("two", ["three","four","five"]), [0,1,0,0,1]) == ["three","two","three","three","two"]
    @test recode.(OneZero(), RestVsOne([2,3,4], 5), [0,1,0,0,1]) == [5,2,5,5,2]
    @test recode.(OneZero(), RestVsOne([:two, :three,:four], :five), [0,1,0,0,1]) == [:five,:two,:five,:five,:two]
    @test recode.(OneZero(), RestVsOne(["two","three","four"], "five"), [0,1,0,0,1]) == ["five","two","five","five","two"]

    @test recode.(OneMinusOne(), OneZero(), [-1,1,-1,-1,1]) == [0,1,0,0,1]
    @test recode.(OneMinusOne(), OneTwo(), [-1,1,-1,-1,1]) == [2,1,2,2,1]
    @test recode.(OneMinusOne(), OneVsRest(2, [3,4,5]), [-1,1,-1,-1,1]) == [3,2,3,3,2]
    @test recode.(OneMinusOne(), OneVsRest(:two, [:three,:four,:five]), [-1,1,-1,-1,1]) == [:three,:two,:three,:three,:two]
    @test recode.(OneMinusOne(), OneVsRest("two", ["three","four","five"]), [-1,1,-1,-1,1]) == ["three","two","three","three","two"]
    @test recode.(OneMinusOne(), RestVsOne([2,3,4], 5), [-1,1,-1,-1,1]) == [5,2,5,5,2]
    @test recode.(OneMinusOne(), RestVsOne([:two, :three,:four], :five), [-1,1,-1,-1,1]) == [:five,:two,:five,:five,:two]
    @test recode.(OneMinusOne(), RestVsOne(["two","three","four"], "five"), [-1,1,-1,-1,1]) == ["five","two","five","five","two"]

    @test recode.(OneTwo(), OneZero(), [2,1,2,2,1]) == [0,1,0,0,1]
    @test recode.(OneTwo(), OneMinusOne(), [2,1,2,2,1]) == [-1,1,-1,-1,1]
    @test recode.(OneTwo(), OneVsRest(2, [3,4,5]), [2,1,2,2,1]) == [3,2,3,3,2]
    @test recode.(OneTwo(), OneVsRest(:two, [:three,:four,:five]), [2,1,2,2,1]) == [:three,:two,:three,:three,:two]
    @test recode.(OneTwo(), OneVsRest("two", ["three","four","five"]), [2,1,2,2,1]) == ["three","two","three","three","two"]
    @test recode.(OneTwo(), RestVsOne([2,3,4], 5), [2,1,2,2,1]) == [5,2,5,5,2]
    @test recode.(OneTwo(), RestVsOne([:two, :three,:four], :five), [2,1,2,2,1]) == [:five,:two,:five,:five,:two]
    @test recode.(OneTwo(), RestVsOne(["two","three","four"], "five"), [2,1,2,2,1]) == ["five","two","five","five","two"]

    @test recode.(OneVsOne(3, 4), OneZero(), [4,3,4,4,3]) == [0,1,0,0,1]
    @test recode.(OneVsOne(3, 4), OneMinusOne(), [4,3,4,4,3]) == [-1,1,-1,-1,1]
    @test recode.(OneVsOne(3, 4), OneTwo(), [4,3,4,4,3]) == [2,1,2,2,1]
    @test recode.(OneVsOne(3, 4), OneVsRest(1, [2, 3]), [4,3,4,4,3]) == [2,1,2,2,1]
    @test recode.(OneVsOne(3, 4), RestVsOne([2,3,4], 5), [4,3,4,4,3]) == [5,2,5,5,2]

    @test recode.(OneVsOne(:three, :four), OneZero(), [:four,:three,:four,:four,:three]) == [0,1,0,0,1]
    @test recode.(OneVsOne(:three, :four), OneMinusOne(), [:four,:three,:four,:four,:three]) == [-1,1,-1,-1,1]
    @test recode.(OneVsOne(:three, :four), OneTwo(), [:four,:three,:four,:four,:three]) == [2,1,2,2,1]
    @test recode.(OneVsOne(:three, :four), OneVsRest(1, [2, 3]), [:four,:three,:four,:four,:three]) == [2,1,2,2,1]
    @test recode.(OneVsOne(:three, :four), RestVsOne([2,3,4], 5), [:four,:three,:four,:four,:three]) == [5,2,5,5,2]

    @test recode.(OneVsOne("three", "four"), OneZero(), ["four","three","four","four","three"]) == [0,1,0,0,1]
    @test recode.(OneVsOne("three", "four"), OneMinusOne(), ["four","three","four","four","three"]) == [-1,1,-1,-1,1]
    @test recode.(OneVsOne("three", "four"), OneTwo(), ["four","three","four","four","three"]) == [2,1,2,2,1]
    @test recode.(OneVsOne("three", "four"), OneVsRest(1, [2, 3]), ["four","three","four","four","three"]) == [2,1,2,2,1]
    @test recode.(OneVsOne("three", "four"), RestVsOne([2,3,4], 5), ["four","three","four","four","three"]) == [5,2,5,5,2]

    @test recode.(OneVsRest(2, [3,4,5]), OneZero(), [3,2,4,5,2]) == [0,1,0,0,1]
    @test recode.(OneVsRest(2, [3,4,5]), OneMinusOne(), [3,2,4,5,2]) == [-1,1,-1,-1,1]
    @test recode.(OneVsRest(2, [3,4,5]), OneTwo(), [3,2,4,5,2]) == [2,1,2,2,1]
    @test recode.(OneVsRest(2, [3,4,5]), OneVsRest(1, [2, 3]), [3,2,4,5,2]) == [2,1,2,2,1]
    @test recode.(OneVsRest(2, [3,4,5]), RestVsOne([2,3,4], 5), [3,2,4,5,2]) == [5,2,5,5,2]

    @test recode.(OneVsRest(:two, [:three,:four,:five]), OneZero(), [:three,:two,:four,:five,:two]) == [0,1,0,0,1]
    @test recode.(OneVsRest(:two, [:three,:four,:five]), OneMinusOne(), [:three,:two,:four,:five,:two]) == [-1,1,-1,-1,1]
    @test recode.(OneVsRest(:two, [:three,:four,:five]), OneTwo(), [:three,:two,:four,:five,:two]) == [2,1,2,2,1]
    @test recode.(OneVsRest(:two, [:three,:four,:five]), OneVsRest(1, [2, 3]), [:three,:two,:four,:five,:two]) == [2,1,2,2,1]
    @test recode.(OneVsRest(:two, [:three,:four,:five]), RestVsOne([2,3,4], 5), [:three,:two,:four,:five,:two]) == [5,2,5,5,2]

    @test recode.(OneVsRest("two", ["three","four","five"]), OneZero(), ["three","two","four","five","two"]) == [0,1,0,0,1]
    @test recode.(OneVsRest("two", ["three","four","five"]), OneMinusOne(), ["three","two","four","five","two"]) == [-1,1,-1,-1,1]
    @test recode.(OneVsRest("two", ["three","four","five"]), OneTwo(), ["three","two","four","five","two"]) == [2,1,2,2,1]
    @test recode.(OneVsRest("two", ["three","four","five"]), OneVsRest(1, [2, 3]), ["three","two","four","five","two"]) == [2,1,2,2,1]
    @test recode.(OneVsRest("two", ["three","four","five"]), RestVsOne([2,3,4], 5), ["three","two","four","five","two"]) == [5,2,5,5,2]

    @test recode.(RestVsOne([2,3,4], 5), OneZero(), [5,2,5,5,3]) == [0,1,0,0,1]
    @test recode.(RestVsOne([2,3,4], 5), OneMinusOne(), [5,2,5,5,3]) == [-1,1,-1,-1,1]
    @test recode.(RestVsOne([2,3,4], 5), OneTwo(), [5,2,5,5,3]) == [2,1,2,2,1]
    @test recode.(RestVsOne([2,3,4], 5), OneVsRest(1, [2, 3]), [5,2,5,5,3]) == [2,1,2,2,1]
    @test recode.(RestVsOne([2,3,4], 5), RestVsOne([1,2,3], 4), [5,2,5,5,3]) == [4,1,4,4,1]

    @test recode.(RestVsOne([:two, :three,:four], :five), OneZero(), [:five,:two,:five,:five,:three]) == [0,1,0,0,1]
    @test recode.(RestVsOne([:two, :three,:four], :five), OneMinusOne(), [:five,:two,:five,:five,:three]) == [-1,1,-1,-1,1]
    @test recode.(RestVsOne([:two, :three,:four], :five), OneTwo(), [:five,:two,:five,:five,:three]) == [2,1,2,2,1]
    @test recode.(RestVsOne([:two, :three,:four], :five), OneVsRest(1, [2, 3]), [:five,:two,:five,:five,:three]) == [2,1,2,2,1]
    @test recode.(RestVsOne([:two, :three,:four], :five), RestVsOne([1,2,3], 4), [:five,:two,:five,:five,:three]) == [4,1,4,4,1]

    @test recode.(RestVsOne(["two","three","four"], "five"), OneZero(), ["five","two","five","five","three"]) == [0,1,0,0,1]
    @test recode.(RestVsOne(["two","three","four"], "five"), OneMinusOne(), ["five","two","five","five","three"]) == [-1,1,-1,-1,1]
    @test recode.(RestVsOne(["two","three","four"], "five"), OneTwo(), ["five","two","five","five","three"]) == [2,1,2,2,1]
    @test recode.(RestVsOne(["two","three","four"], "five"), OneVsRest(1, [2, 3]), ["five","two","five","five","three"]) == [2,1,2,2,1]
    @test recode.(RestVsOne(["two","three","four"], "five"), RestVsOne([1,2,3], 4), ["five","two","five","five","three"]) == [4,1,4,4,1]
end


text:
using NNlib, Test

maxpool_answer_dict = Dict(
    1 => Dict(
        "y"          => [2, 4.],
        "y_nostride" => [2, 3, 4, 5.],
        "y_pad"      => [1, 3, 5.],

        "dx"          => [0, 2, 0, 4, 0.],
        "dx_nostride" => [0, 2, 3, 4, 5.],
        "dx_pad"      => [1, 0, 3, 0, 5.],
    ),
    2 => Dict(
        "y" => [
            7 17.;
            9 19.
        ],
        "y_nostride" => [
            7  12 17;
            8  13 18;
            9  14 19;
            10 15 20.
        ],
        "y_pad" => [
            1  11 16;
            3  13 18;
            5  15 20.
        ],

        "dx" => [
            0 0 0  0;
            0 7 0 17;
            0 0 0  0;
            0 9 0 19;
            0 0 0  0.
        ],
        "dx_nostride" => [
            0  0  0  0;
            0  7 12 17;
            0  8 13 18;
            0  9 14 19;
            0 10 15 20.
        ],
        "dx_pad"      => [
            1 0 11 16;
            0 0  0  0;
            3 0 13 18;
            0 0  0  0;
            5 0 15 20.
        ],
    ),
    3 => Dict(
        "y" => reshape([
            27, 29,
            37, 39.
        ], (2, 2, 1)),
        "y_nostride" => reshape([
            27, 28, 29, 30,
            32, 33, 34, 35,
            37, 38, 39, 40,

            47, 48, 49, 50,
            52, 53, 54, 55,
            57, 58, 59, 60.
        ], (4, 3, 2)),
        "y_pad" => reshape([
             1,  3, 5,
            11, 13, 15,
            16, 18, 20,

            41, 43, 45,
            51, 53, 55,
            56, 58, 60.
        ], (3, 3, 2)),

        "dx" => reshape([
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,

            0,  0, 0,  0, 0,
            0, 27, 0, 29, 0,
            0,  0, 0,  0, 0,
            0, 37, 0, 39, 0,

            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0.
        ], (5, 4, 3)),
        "dx_nostride" => reshape([
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,

            0,  0,  0,  0,  0,
            0, 27, 28, 29, 30,
            0, 32, 33, 34, 35,
            0, 37, 38, 39, 40,

            0,  0,  0,  0,  0,
            0, 47, 48, 49, 50,
            0, 52, 53, 54, 55,
            0, 57, 58, 59, 60.
        ], (5, 4, 3)),
        "dx_pad" => reshape([
             1, 0,  3, 0,  5,
             0, 0,  0, 0,  0,
            11, 0, 13, 0, 15,
            16, 0, 18, 0, 20,

            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,

            41, 0, 43, 0, 45,
             0, 0,  0, 0,  0,
            51, 0, 53, 0, 55,
            56, 0, 58, 0, 60.
        ], (5, 4, 3)),
    )
)

meanpool_answer_dict = Dict(
    1 => Dict(
        "y"          => [1.5, 3.5],
        "y_nostride" => [1.5, 2.5, 3.5, 4.5],
        "y_pad"      => [0.5, 2.5, 4.5],

        "dx"          => [0.75, 0.75, 1.75, 1.75,  0.0],
        "dx_nostride" => [0.75,  2.0,  3.0,  4.0, 2.25],
        "dx_pad"      => [0.25, 1.25, 1.25, 2.25, 2.25],
    ),
    2 => Dict(
        "y" => [
            4.0 14.0;
            6.0 16.0
        ],
        "y_nostride" => [
            4.0  9.0 14.0
            5.0 10.0 15.0
            6.0 11.0 16.0
            7.0 12.0 17.0
        ],
        "y_pad" => [
            0.25  4.25 4.0
            1.25 10.0  8.75
            2.25 12.0  9.75
        ],

        "dx" => [
            1.0 1.0 3.5 3.5;
            1.0 1.0 3.5 3.5;
            1.5 1.5 4.0 4.0;
            1.5 1.5 4.0 4.0;
            0.0 0.0 0.0 0.0
        ],
        "dx_nostride" => [
            1.0  3.25  5.75 3.5;
            2.25 7.0  12.0  7.25;
            2.75 8.0  13.0  7.75;
            3.25 9.0  14.0  8.25;
            1.75 4.75  7.25 4.25
        ],
        "dx_pad"      => [
            0.0625 1.0625 1.0625 1.0;
            0.3125 2.5    2.5    2.1875;
            0.3125 2.5    2.5    2.1875;
            0.5625 3.0    3.0    2.4375;
            0.5625 3.0    3.0    2.4375
        ],
    ),
    3 => Dict(
        "y" => reshape([
            14.0, 16.0,
            24.0, 26.0
        ], (2, 2, 1)),
        "y_nostride" => reshape([
            14.0, 15.0, 16.0, 17.0,
            19.0, 20.0, 21.0, 22.0,
            24.0, 25.0, 26.0, 27.0,

            34.0, 35.0, 36.0, 37.0,
            39.0, 40.0, 41.0, 42.0,
            44.0, 45.0, 46.0, 47.0
        ], (4, 3, 2)),
        "y_pad" => reshape([
            0.125, 0.625, 1.125,
            2.125, 5.0,   6.0,
            2.0,   4.375, 4.875,

             7.75, 16.25, 17.25,
            19.25, 40.0,  42.0,
            11.5,  23.75, 24.75,
        ], (3, 3, 2)),

        "dx" => reshape([
            1.75, 1.75, 2.0, 2.0, 0.0,
            1.75, 1.75, 2.0, 2.0, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,

            1.75, 1.75, 2.0, 2.0, 0.0,
            1.75, 1.75, 2.0, 2.0, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,

            0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0,
        ], (5, 4, 3)),
        "dx_nostride" => reshape([
            1.75,   3.625,  3.875,  4.125, 2.125,
            4.125,  8.5,    9.0,    9.5,   4.875,
            5.375, 11.0,   11.5,   12.0,   6.125,
            3.0,    6.125,  6.375,  6.625, 3.375,

             6.0,  12.25, 12.75, 13.25,  6.75,
            13.25, 27.0,  28.0,  29.0,  14.75,
            15.75, 32.0,  33.0,  34.0,  17.25,
             8.5,  17.25, 17.75, 18.25,  9.25,

             4.25,   8.625,  8.875,  9.125,  4.625,
             9.125, 18.5,   19.0,   19.5,    9.875,
            10.375, 21.0,   21.5,   22.0,   11.125,
             5.5,   11.125, 11.375, 11.625,  5.875
        ], (5, 4, 3)),
        "dx_pad" => reshape([
            0.015625, 0.078125, 0.078125, 0.140625, 0.140625,
            0.265625, 0.625,    0.625,    0.75,     0.75,
            0.265625, 0.625,    0.625,    0.75,     0.75,
            0.25,     0.546875, 0.546875, 0.609375, 0.609375,

            0.96875, 2.03125, 2.03125, 2.15625, 2.15625,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            1.4375,  2.96875, 2.96875, 3.09375, 3.09375,

            0.96875, 2.03125, 2.03125, 2.15625, 2.15625,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            1.4375,  2.96875, 2.96875, 3.09375, 3.09375
        ], (5, 4, 3)),
    )
)

for rank in (1, 2, 3)
    @testset "pool$(rank)d" begin
        for (pool, ∇pool, answer_dict) in (
                # Main API name
                (maxpool, ∇maxpool, maxpool_answer_dict),
                (meanpool, ∇meanpool, meanpool_answer_dict),

                # _direct name
                (NNlib.maxpool_direct, NNlib.∇maxpool_direct, maxpool_answer_dict),
                (NNlib.meanpool_direct, NNlib.∇meanpool_direct, meanpool_answer_dict),
            )

            @testset "$(pool)$(rank)d" begin
                y = answer_dict[rank]["y"]
                y_nostride = answer_dict[rank]["y_nostride"]
                y_pad = answer_dict[rank]["y_pad"]
                dx = answer_dict[rank]["dx"]
                dx_nostride = answer_dict[rank]["dx_nostride"]
                dx_pad = answer_dict[rank]["dx_pad"]

                x = reshape(Float64[1:prod(size(dx));], size(dx)..., 1, 1)

                # A "drop channels and batch dimension" helper
                ddims(x) = dropdims(x, dims=(rank+1, rank+2))

                # Let's ensure that a 1x1x1 pooling kernel always just returns `x`
                @test pool(x, PoolDims(x, 1)) == x

                # Test vanilla pooling
                pdims = PoolDims(x, 2)
                y_hat = pool(x, pdims)
                @test ddims(y_hat) == y
                @test ddims(∇pool(y_hat, y_hat, x, pdims)) == dx

                # Strided pooling
                pdims = PoolDims(x, 2; stride=1)
                y_hat = pool(x, pdims)
                @test ddims(y_hat) == y_nostride
                @test ddims(∇pool(y_hat, y_hat, x, pdims)) == dx_nostride

                # Padded pooling
                pdims = PoolDims(x, 2; padding=1)
                y_hat = pool(x, pdims)
                @test ddims(y_hat) == y_pad
                @test ddims(∇pool(y_hat, y_hat, x, pdims)) == dx_pad
            end
        end
    end
end

x = rand(10, 10, 3, 10)
@test size(maxpool(x, (2, 2))) == (5, 5, 3, 10)
@test size(maxpool(x, (2, 2); pad = (2, 2), stride = (2, 2))) == (7, 7, 3, 10)
@test size(meanpool(x, (2, 2))) == (5, 5, 3, 10)
@test size(meanpool(x, (2, 2); pad = (2, 2), stride = (2, 2))) == (7, 7, 3, 10)

# Add another test for 2d maxpool that uses an odd-length size:
@testset "Issue #133" begin
    x = reshape([(1.:9.)...], 3, 3, 1, 1)
    pdims = PoolDims(size(x), (2,2), padding = (1,1), stride = (2,2))
    y = maxpool(x, pdims)

    dy = y .* 0 .+ 1
    dx = ∇maxpool(dy, y, x, pdims)
    @test dx[:,:,1,1] == [1.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 1.0]
end

ext: jl


text:
"""
    finetune!(learner, nepochs[, base_lr = 0.002; kwargs...])

Behaves like the fastai implementation
[`fastai.Learner.fine_tune`](https://github.com/fastai/fastai/blob/master/fastai/callback/schedule.py#L151).

## Keyword arguments

- `freezeepochs = 1`: Number of epochs to train with the backbone completely frozen.
- `grouper = FastAI.defaultgrouper(learner.model)`: [`ParamGrouper`](#) which assigns
    groups `1` (backbone) or `2` (head) for every parameter in `learner.model`. The
    default expects `learner.model` to be a `Chain(backbone, head)`.
- `backbone_factor = 0.1`: Factor by which updates to backbone model are discounted
    during the second phase of training.

Any additional keyword arguments are passed to [`fitonecycle!`](#).
"""
function finetune!(
        learner,
        nepochs,
        base_lr=0.002;
        freezeepochs=1,
        grouper=defaultgrouper(learner.model),
        backbone_factor=0.1,
        div=5,
        kwargs...)

    # Freeze backbone and train head
    foptim = frozen_optimizer(learner.optimizer, grouper, learner.model)
    withfields(learner, optimizer=foptim) do
        fitonecycle!(learner, freezeepochs, base_lr, pct_start=0.99; kwargs...)
    end

    # Use discriminative learning rates on backbone and train some more
    doptim = discrlr_optimizer(learner.optimizer, grouper, learner.model, backbone_factor)
    withfields(learner, optimizer=doptim) do
        fitonecycle!(
            learner, nepochs, base_lr / 2;
            div=div, kwargs...)
    end

    return learner
end

"""
    frozen_optimizer(optim, grouper, model)

Create an optimizer that only updates parameters which [`ParamGrouper`](#)
puts into group `2`.
"""
frozen_optimizer(optim, grouper, model) = discrlr_optimizer(optim, grouper, model, 0.)

"""
    frozen_optimizer(optim, grouper, model, factor)

Create an optimizer that discounts updates parameters which [`ParamGrouper`](#)
puts into group `1` by `factor`.
"""
function discrlr_optimizer(optim, grouper, model, factor)
    paramgroups = ParamGroups(grouper, model)
    return Optimiser(
        DiscriminativeLRs(paramgroups, Dict(1 => factor, 2 => 1.)),
        optim,
    )
end

function defaultgrouper(model)
    if !((model isa Chain) && length(model) == 2)
        error(
            "Cannot freeze `learner.model` automatically since it is not a `Chain`.
            Please provide a `ParamGrouper` with the `grouper` keyword argument.
            The `grouper` should assign groups `1` (backbone) and `2` (head).
            ")
    else
        return IndexGrouper([1:length(model) - 1, length(model)])
    end
end


text:
# # Little Group in 3D

using TightBindingLattice

unitcell = makeunitcell([1.0 0.0 0.0; 0.0 1.0 0.0; 0 0 1]; SiteType=String)
addsite!(unitcell, "Ox", FractCoord([0,0,0], [0.5, 0.0, 0.0]))
addsite!(unitcell, "Oy", FractCoord([0,0,0], [0.0, 0.5, 0.0]))

# Currently, makelattice in three-dimension is not well supported

lattice = makelattice(unitcell, [4 0 0; 0 4 0; 0 0 3], [1 0 0; 0 1 0; 0 0 1])
tsym = TranslationSymmetry(lattice.orthocube, [1 0 0; 0 1 0; 0 0 1])
psym = project(PointSymmetryDatabase.get(15), [1 0 0; 0 1 0; 0 0 1])

for idx in 1:num_irreps(tsym)
    kf = tsym.fractional_momenta[idx]
    k = lattice.unitcell.reducedreciprocallatticevectors * kf
    psym_little1 = little_symmetry(tsym, idx, psym)
    psym_little2 = TightBindingLattice.little_symmetry_iso(tsym, idx, psym)
    println("- irrep_index: $(idx)")
    println("  momentum: $(k)")
    println("  little_point_group1: { name: \"$(psym_little1.hermann_mauguin)\", order: $(group_order(psym_little1)) }")
    println("  little_point_group2: { name: \"$(psym_little2.hermann_mauguin)\", order: $(group_order(psym_little2)) }")
    println("  is_psym_compatible: $(iscompatible(tsym, idx, psym))")
    println("  is_psym_little1_compatible: $(iscompatible(tsym, idx, psym_little1))")
    println("  is_psym_little2_compatible: $(iscompatible(tsym, idx, psym_little2))")
end


text:
## --- Load (and install if necessary) the Chron.jl package

    try
        using Chron
    catch
        using Pkg
        Pkg.add("Chron")
        using Chron
    end

    using Statistics, StatsBase, DelimitedFiles, SpecialFunctions
    using KernelDensity: kde
    using Plots; gr(); default(fmt = :png)

## --- Test Bayesian eruption age estimation with a synthetic dataset

    dt_sigma = 30; # Timescale relative to analytical uncertainty
    N = 42; # Number of zircons

    # Draw set of pseudorandom ages from MELTS volcanic zircon distribution,
    # with true minimum == 0 and analytical sigma == 1
    ages = draw_from_distribution(MeltsVolcanicZirconDistribution,N).*dt_sigma + randn(N);
    uncert = ones(N)

    # Calculate the weighted mean age
    (wx, wsigma, mswd) = awmean(ages, uncert)

    h = plot(ylabel="Time before eruption (sigma)", xlabel="N", fg_color_legend=:white, legend=:topleft)
    plot!(h,1:length(ages),sort(ages),yerror=uncert*2,seriestype=:scatter, label="Synthetic zircon ages")
    plot!(h,collect(xlims()), [0,0], label="Eruption age")
    plot!(h,collect(xlims()), [wx,wx], label="Weighted mean, MSWD $(round(mswd,sigdigits=2))")

## --- Calculate bootstrapped $\ \mathcal{\vec{f}}_{xtal}(t_r)$

    # Bootstrap the crystallization distribution,
    # accounting for any expected analytical "tail" beyond eruption/deposition
    dist = BootstrapCrystDistributionKDE(ages, uncert)
    dist ./= mean(dist) # Normalize

    # Plot bootstrapped distribution
    plot(range(0,1.3,length=length(dist)),dist, label="bootstrapped", ylabel="Probability Density", xlabel="Time before eruption (scaled)", legend=:bottomleft, fg_color_legend=:white)
    plot!(range(0,1,length=100),MeltsVolcanicZirconDistribution,label="original")

## ---  Run MCMC to estimate eruption/deposition age distribution of synthetic dataset
    # Configure model
    nsteps = 200000 # Length of Markov chain
    burnin = 100000 # Number of steps to discard at beginning of Markov chain

    # Run MCMC
    tminDist = metropolis_min(nsteps,dist,ages,uncert)

    # Print results
    AgeEst = mean(tminDist[burnin:end]);
    AgeEst_sigma = std(tminDist[burnin:end]);
    print("\nEstimated eruption age of synthetic dataset:\n $AgeEst +/- $(2*AgeEst_sigma) Ma (2σ)\n (True synthetic age 0 Ma)")

    # Plot results
    h = histogram(tminDist[burnin:end],nbins=50,label="Posterior distribution",xlabel="Eruption Age (Ma)",ylabel="N")
    plot!(h,[0,0],collect(ylims()),line=(3),label="True (synthetic) age",fg_color_legend=:white)
    plot!(h,[wx,wx],collect(ylims()),line=(3),label="Weighted mean age",legend=:topright)
    display(h)

## --- Estimate eruption age for real zircon data
    # The example dataset here is from [Wotzlaw et al., 2013](https://doi.org/10.1130/G34366.1) FCT+MLX
    # Age and one-sigma uncertainty.
    ages = [28.196, 28.206, 28.215, 28.224, 28.232, 28.241, 28.246, 28.289, 28.308, 28.332, 28.341, 28.359, 28.379, 28.383, 28.395, 28.4, 28.405, 28.413, 28.415, 28.418, 28.42, 28.422, 28.428, 28.452, 28.454, 28.454, 28.458, 28.468, 28.471, 28.475, 28.482, 28.485, 28.502, 28.52, 28.551, 28.561, 28.565, 28.582, 28.584, 28.586, 28.611, 28.638, 28.655]
    uncert = [0.019, 0.0155, 0.019, 0.0215, 0.018, 0.023, 0.013, 0.029, 0.0175, 0.0315, 0.0095, 0.0245, 0.0255, 0.0175, 0.0235, 0.014, 0.021, 0.022, 0.0125, 0.0135, 0.016, 0.0195, 0.0175, 0.0125, 0.01, 0.014, 0.015, 0.0205, 0.0155, 0.011, 0.0115, 0.0185, 0.0255, 0.014, 0.0125, 0.013, 0.015, 0.014, 0.012, 0.016, 0.0215, 0.0125, 0.0215]

    # Sort by age (just to make rank-order plots prettier)
    t = sortperm(ages)
    ages = ages[t];
    uncert = uncert[t];

## --- Calculate bootstrapped $\ \mathcal{\vec{f}}_{xtal}(t_r)$
    # Bootstrap the crystallization distribution,
    # accounting for any expected analytical "tail" beyond eruption/deposition
    dist = BootstrapCrystDistributionKDE(ages, uncert)
    dist ./= mean(dist) # Normalize

    # Plot bootstrapped distribution
    plot(range(0,1,length=length(dist)),dist, label="bootstrapped f_xtal", ylabel="Probability Density", xlabel="Time before eruption (unitless)", fg_color_legend=:white)

## --  Run MCMC to estimate eruption age
    # Configure model
    nsteps = 400000; # Length of Markov chain
    burnin = 150000; # Number of steps to discard at beginning of Markov chain

    # Run MCMC
    tminDist = metropolis_min(nsteps,dist,ages,uncert);
    # Consider also:
    # tminDist = metropolis_min(nsteps,UniformDistribution,ages,uncert);
    # tminDist = metropolis_min(nsteps,TriangularDistribution,ages,uncert);
    # tminDist = metropolis_min(nsteps,HalfNormalDistribution,ages,uncert);
    # tminDist = metropolis_min(nsteps,MeltsVolcanicZirconDistribution,ages,uncert);

    # Print results
    AgeEst = mean(tminDist[burnin:end]);
    AgeEst_sigma = std(tminDist[burnin:end]);
    print("\nEstimated eruption age:\n $AgeEst +/- $(2*AgeEst_sigma) Ma (2σ)\n")

    # Plot results
    h = histogram(tminDist[burnin:end],nbins=100,label="Posterior distribution",xlabel="Eruption Age (Ma)",ylabel="N",legend=:topleft,fg_color_legend=:white)
    # plot!(h,[wx,wx],collect(ylims()),line=(3),label="Weighted mean age",legend=:topleft)
    display(h)

## ---  Plot eruption age estimate relative to rank-order plot of raw zircon ages
    h = plot(ylabel="Age (Ma)", xlabel="N", legend=:topleft, fg_color_legend=:white)
    plot!(h,1:length(ages),ages,yerror=uncert*2,seriestype=:scatter, markerstrokecolor=:auto, label="Observed ages")
    plot!(h,[length(ages)],[AgeEst],yerror=2*AgeEst_sigma, markerstrokecolor=:auto, label="Bayesian eruption age estimate",color=:red)
    plot!(h,collect(xlims()),[AgeEst,AgeEst],color=:red, label="")

## ---


text:
using Pkg
#=Executing this file under windows, creates an executable:
DemoLoopVectorization_compiled_for_windows\bin\DemoLoopVectorization.exe
or under linux creates:
DemoLoopVectorization_compiled_for_linux/bin/DemoLoopVectorization
=#

packagename = "DemoLoopVectorization"

Pkg.activate()#Activate default enviroment
Pkg.add("PackageCompiler")

using PackageCompiler

if Sys.iswindows()
    os = "windows"
elseif Sys.islinux()
    os = "linux"
else
    throw("Operating system not supported")
end

srcfolder = @__DIR__
packagefolder = normpath(joinpath(srcfolder, ".."))
startfolder = normpath(joinpath(packagefolder, ".."))
compiled_app = joinpath(startfolder, "$(packagename)_compiled_for_$os")
precompile_execution_file = joinpath(srcfolder,"precompile_execution_file.jl")

cd(startfolder)
Pkg.activate(packagename)
Pkg.instantiate()
create_app(packagename, compiled_app, force=true,
           precompile_execution_file=precompile_execution_file,
           include_transitive_dependencies=false)


text:
# This file is a part of SimilaritySearch.jl

export LocalSearchAlgorithm, SearchGraph, SearchGraphPools, SearchGraphCallbacks, VisitedVertices, NeighborhoodReduction, index!, push_item!
export BeamSearch, Callback
export KDisjointHints, DisjointHints, RandomHints
export IdentityNeighborhood, DistalSatNeighborhood, SatNeighborhood, find_neighborhood, push_neighborhood!, NeighborhoodSize

include("graph.jl")
include("rebuild.jl")
include("insertions.jl")

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
"""
    Return an offset vector `p`, where p[i] is the prior
    probability of any hypothesis such that sum_{j=1}^m γj = i,
    according to Womack (2015)'s proposal
"""
function ph0(m, ζ)
    p = OffsetArray([zeros(m); 1.0], 0:m)
    for l = m-1:-1:0
        for j = 1:m-l
            p[l] += ζ * p[l + j] * binomial(l + j, l)
        end
    end
    p /= sum(p)
    for l = 0:m-1
        p[l] /= binomial(m, l)
    end
    return p
end

function γcode(γ)
    code = 1
    @inbounds for j = 2:length(γ)
        code += γ[j] * 2^(j - 2)
    end
    code
end

function γvector(J, code)
    γ = zeros(Int, J)
    code -= 1
    for j = J:-1:2
        if (code ÷ 2^(j-2)) > 0
            code -= 2^(j-2)
            γ[j] = 1
        end
    end
    return γ
end

"Return the vector γ associated with code `code`, given `J`"
γvec(J, code) = γvector(J, code)[2:J]

"Return the vector γ associated with code `code`, given `J`, as a string"
γstr(J, code) = string(γvec(J, code)...)


text:
# ---
# title: 748. Shortest Completing Word
# id: problem748
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/shortest-completing-word/description/>
# hidden: true
# ---
#
# Given a string `licensePlate` and an array of strings `words`, find the
# **shortest completing** word in `words`.
#
# A **completing** word is a word that **contains all the letters** in
# `licensePlate`. **Ignore numbers and spaces** in `licensePlate`, and treat
# letters as **case insensitive**. If a letter appears more than once in
# `licensePlate`, then it must appear in the word the same number of times or
# more.
#
# For example, if `licensePlate`` = "aBc 12c"`, then it contains letters `'a'`,
# `'b'` (ignoring case), and `'c'` twice. Possible **completing** words are
# `"abccdef"`, `"caaacab"`, and `"cbca"`.
#
# Return _the shortest **completing** word in _`words` _._ It is guaranteed an
# answer exists. If there are multiple shortest **completing** words, return the
# **first** one that occurs in `words`.
#
#
#
# **Example 1:**
#
#
#
#     Input: licensePlate = "1s3 PSt", words = ["step","steps","stripe","stepple"]
#     Output: "steps"
#     Explanation: licensePlate contains letters 's', 'p', 's' (ignoring case), and 't'.
#     "step" contains 't' and 'p', but only contains 1 's'.
#     "steps" contains 't', 'p', and both 's' characters.
#     "stripe" is missing an 's'.
#     "stepple" is missing an 's'.
#     Since "steps" is the only word containing all the letters, that is the answer.
#
#
# **Example 2:**
#
#
#
#     Input: licensePlate = "1s3 456", words = ["looks","pest","stew","show"]
#     Output: "pest"
#     Explanation: licensePlate only contains the letter 's'. All the words contain 's', but among these "pest", "stew", and "show" are shortest. The answer is "pest" because it is the word that appears earliest of the 3.
#
#
# **Example 3:**
#
#
#
#     Input: licensePlate = "Ah71752", words = ["suggest","letter","of","husband","easy","education","drug","prevent","writer","old"]
#     Output: "husband"
#
#
# **Example 4:**
#
#
#
#     Input: licensePlate = "OgEu755", words = ["enough","these","play","wide","wonder","box","arrive","money","tax","thus"]
#     Output: "enough"
#
#
# **Example 5:**
#
#
#
#     Input: licensePlate = "iMSlpe4", words = ["claim","consumer","student","camera","public","never","wonder","simple","thought","use"]
#     Output: "simple"
#
#
#
#
# **Constraints:**
#
#   * `1 <= licensePlate.length <= 7`
#   * `licensePlate` contains digits, letters (uppercase or lowercase), or space `' '`.
#   * `1 <= words.length <= 1000`
#   * `1 <= words[i].length <= 15`
#   * `words[i]` consists of lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module SPHtoGrid

    using Distributed
    using SPHKernels
    using Unitful
    using Printf
    using ProgressMeter

    function output_time(t1, t2)
        return @sprintf("%0.3e", Float64((t2-t1))*1.e-9)
    end

    # shared functionality
    include("shared/parameters.jl")
    include("shared/reconstruct_grid.jl")
    include("shared/periodic_mapping.jl")
    include("shared/rotate_particles.jl")
    include("shared/io.jl")

    # multi-core functionality
    include("parallel/domain_decomp.jl")

    # functions for smac
    include("smac/smac1_utility.jl")
    include("smac/smac2_utility.jl")

    # tsc interpolation
    include("tsc_interpolation/tsc_interpolation.jl")

    # cic interpolation
    include("cic_interpolation/filter_shift.jl")
    include("cic_interpolation/reduce_image.jl")
    include("cic_interpolation/cic_shared.jl")
    include("cic_interpolation/cic_2D.jl")
    include("cic_interpolation/cic_3D.jl")
    include("cic_interpolation/weight_functions.jl")
    include("cic_interpolation/cic_interpolation.jl")

    # effect functions
    include("effects/constants.jl")
    include("effects/density.jl")
    include("effects/synchrotron_analytic.jl")
    include("effects/synchrotron_spectrum.jl")
    include("effects/sz_effect.jl")
    include("effects/x_ray.jl")

    export mappingParameters,                         # parameters for SPH mapping
           sphMapping,                                # main function for mapping
           filter_particles_in_image,                 # helper function to preselect particles
           get_map_grid_2D,
           get_map_grid_3D,
           read_smac1_binary_info,
           read_smac1_binary_image,
           read_smac2_image,
           read_smac2_info,
           write_smac1_par,
           write_smac2_par,
           # rotate particles
           rotate_3D, rotate_3D!,
           project_along_axis,
           project_along_axis!,
           # Weight functions
           part_weight_one,
           part_weight_physical,
           part_weight_emission,
           part_weight_XrayBand,
           part_weight_spectroscopic,
           # effect functions
           density_2D,
           comptonY,
           kinetic_SZ,
           thermal_SZ,
           x_ray_emission,
           analytic_synchrotron_emission,
           spectral_synchrotron_emission,
           # IO
           read_fits_image,
           write_fits_image

end # module


text:
let
    doc = open("$testdir/example.html") do example
        example |> readstring |> parsehtml
    end
    io = IOBuffer()
    print(io, doc)
    seek(io, 0)
    newdoc = io |> readstring |> parsehtml
    @test newdoc == doc
end


text:
# We wish to extend operations to identically named methods dispatched
# on `Machine`s and `NodalMachine`s. For example, we have from the model API
#
# `predict(model::M, fitresult, X) where M<:Supervised`
#
# but want also want
#
# `predict(machine::Machine, X)` where `X` is data
#
# and "networks.jl" requires us to define
#
# `predict(machine::NodalMachine, X)` where `X` is data
#
# and we would like the syntactic sugar (for `X` a node):
#
# `predict(machine::NodalMachine, X::Node)=node(predict, machine, X)`
#
# (If an operation has zero arguments, we cannot achieve the last
# desire because of ambiguity with the preceding one.)

## TODO: need to add checks on the arguments of
## predict(::AbstractMachine, ) and transform(::AbstractMachine, )

for operation in (:predict, :predict_mean, :predict_mode, :predict_median,
                  :transform, :inverse_transform)
    ex = quote
        function $(operation)(machine::AbstractMachine, args...)
            if isdefined(machine, :fitresult)
                return $(operation)(machine.model, machine.fitresult, args...)
            else
                throw(error("$machine has not been trained."))
            end
        end
        $(operation)(machine::Machine; rows=:) =
            $(operation)(machine, selectrows(machine.args[1], rows))
        $(operation)(machine::NodalMachine, args::AbstractNode...) =
            node($(operation), machine, args...)
    end
    eval(ex)
end

# the zero argument special cases:
function fitted_params(machine::AbstractMachine)
    if isdefined(machine, :fitresult)
        return fitted_params(machine.model, machine.fitresult)
    else
        throw(error("$machine has not been trained."))
    end
end


text:
# returns Magnus and Neudecker's commutation matrix of dimensions n by m
function commutation(n::Int64, m::Int64)
  k = reshape(kron(vec(eye(n)), eye(m)), n*m, n*m)
  return k
end


text:
@shorthands meanplot
@shorthands autocorplot
@shorthands mixeddensity
@shorthands pooleddensity
@shorthands traceplot
@shorthands corner
@userplot RidgelinePlot
@userplot ForestPlot

struct _TracePlot; c; val; end
struct _MeanPlot; c; val;  end
struct _DensityPlot; c; val;  end
struct _HistogramPlot; c; val;  end
struct _AutocorPlot; lags; val;  end

# define alias functions for old syntax
const translationdict = Dict(
                        :traceplot => _TracePlot,
                        :meanplot => _MeanPlot,
                        :density => _DensityPlot,
                        :histogram => _HistogramPlot,
                        :autocorplot => _AutocorPlot,
                        :pooleddensity => _DensityPlot
                      )

const supportedplots = push!(collect(keys(translationdict)), :mixeddensity, :corner)

@recipe f(c::Chains, s::Symbol) = c, [s]

@recipe function f(
    chains::Chains, i::Int;
    colordim = :chain,
    barbounds = (-Inf, Inf),
    maxlag = nothing,
    append_chains = false
)
    st = get(plotattributes, :seriestype, :traceplot)
    c = append_chains || st == :pooleddensity ? pool_chain(chains) : chains

    if colordim == :parameter
        title --> "Chain $(MCMCChains.chains(c)[i])"
        label --> string.(names(c))
        val = c.value[:, :, i]
    elseif colordim == :chain
        title --> string(names(c)[i])
        label --> map(x -> "Chain $x", MCMCChains.chains(c))
        val = c.value[:, i, :]
    else
        throw(ArgumentError("`colordim` must be one of `:chain` or `:parameter`"))
    end

    if st == :mixeddensity || st == :pooleddensity
        discrete = indiscretesupport(c, barbounds)
        st = if colordim == :chain
            discrete[i] ? :histogram : :density
        else
            # NOTE: It might make sense to overlay histograms and density plots here.
            :density
        end
        seriestype := st
    end

    if st == :autocorplot
        lags = 0:(maxlag === nothing ? round(Int, 10 * log10(length(range(c)))) : maxlag)
        ac = autocor(c; sections = nothing, lags = lags)
        ac_mat = convert(Array, ac)
        val = colordim == :parameter ? ac_mat[:, :, i]' : ac_mat[i, :, :]
        _AutocorPlot(lags, val)
    elseif st ∈ supportedplots
        translationdict[st](c, val)
    else
        range(c), val
    end
end

@recipe function f(p::_DensityPlot)
    xaxis --> "Sample value"
    yaxis --> "Density"
    trim --> true
    [collect(skipmissing(p.val[:,k])) for k in 1:size(p.val, 2)]
end

@recipe function f(p::_HistogramPlot)
    xaxis --> "Sample value"
    yaxis --> "Frequency"
    fillalpha --> 0.7
    bins --> 25
    trim --> true
    [collect(skipmissing(p.val[:,k])) for k in 1:size(p.val, 2)]
end

@recipe function f(p::_MeanPlot)
    seriestype := :path
    xaxis --> "Iteration"
    yaxis --> "Mean"
    range(p.c), cummean(p.val)
end

@recipe function f(p::_AutocorPlot)
    seriestype := :path
    xaxis --> "Lag"
    yaxis --> "Autocorrelation"
    p.lags, p.val
end

@recipe function f(p::_TracePlot)
    seriestype := :path
    xaxis --> "Iteration"
    yaxis --> "Sample value"
    range(p.c), p.val
end

@recipe function f(
    chains::Chains,
    parameters::AbstractVector{Symbol};
    colordim = :chain
)
    colordim != :chain &&
        error("Symbol names are interpreted as parameter names, only compatible with ",
              "`colordim = :chain`")

    ret = indexin(parameters, names(chains))
    any(y === nothing for y in ret) && error("Parameter not found")

    return chains, Int.(ret)
end

@recipe function f(
    chains::Chains,
    parameters::AbstractVector{<:Integer} = Int[];
    sections = _default_sections(chains),
    width = 500,
    height = 250,
    colordim = :chain,
    append_chains = false
)
    _chains = isempty(parameters) ? Chains(chains, _clean_sections(chains, sections)) : chains
    c = append_chains ? pool_chain(_chains) : _chains
    ptypes = get(plotattributes, :seriestype, (:traceplot, :mixeddensity))
    ptypes = ptypes isa Symbol ? (ptypes,) : ptypes
    @assert all(ptype -> ptype ∈ supportedplots, ptypes)
    ntypes = length(ptypes)
    nrows, nvars, nchains = size(c)
    isempty(parameters) && (parameters = colordim == :chain ? (1:nvars) : (1:nchains))
    N = length(parameters)

    if :corner ∉ ptypes
        size --> (ntypes*width, N*height)
        legend --> false

        multiple_plots = N * ntypes > 1
        if multiple_plots
            layout := (N, ntypes)
        end

        i = 0
        for par in parameters
            for ptype in ptypes
                i += 1

                @series begin
                    if multiple_plots
                        subplot := i
                    end
                    colordim := colordim
                    seriestype := ptype
                    c, par
                end
            end
        end
    else
        ntypes > 1 && error(":corner is not compatible with multiple seriestypes")
        Corner(c, names(c)[parameters])
    end
end

struct Corner
    c
    parameters
end

@recipe function f(corner::Corner)
    label --> permutedims(corner.parameters)
    compact --> true
    size --> (600, 600)
    ar = collect(Array(corner.c.value[:, corner.parameters,i]) for i in chains(corner.c))
    RecipesBase.recipetype(:cornerplot, vcat(ar...))
end

function _compute_plot_data(
    i::Integer,
    chains::Chains,
    par_names::AbstractVector{Symbol};
    hpd_val = [0.05, 0.2],
    q = [0.1, 0.9],
    spacer = 0.4,
    _riser = 0.2,
    barbounds = (-Inf, Inf),
    show_mean = true,
    show_median = true,
    show_qi = false,
    show_hpdi = true,
    fill_q = true,
    fill_hpd = false,
    ordered = false
)

    chain_dic = Dict(zip(quantile(chains)[:,1], quantile(chains)[:,4]))
    sorted_chain = sort(collect(zip(values(chain_dic), keys(chain_dic))))
    sorted_par = [sorted_chain[i][2] for i in 1:length(par_names)]
    par = (ordered ? sorted_par : par_names)
    hpdi = sort(hpd_val)

    chain_sections = MCMCChains.group(chains, Symbol(par[i]))
    chain_vec = vec(chain_sections.value.data)
    lower_hpd = [MCMCChains.hpd(chain_sections, alpha = hpdi[j]).nt.lower
        for j in 1:length(hpdi)]
    upper_hpd = [MCMCChains.hpd(chain_sections, alpha = hpdi[j]).nt.upper
        for j in 1:length(hpdi)]
    h = _riser + spacer*(i-1)
    qs = quantile(chain_vec, q)
    k_density = kde(chain_vec)
    if fill_hpd
        x_int = filter(x -> lower_hpd[1][1] <= x <= upper_hpd[1][1], k_density.x)
        val = pdf(k_density, x_int) .+ h
    elseif fill_q
        x_int = filter(x -> qs[1] <= x <= qs[2], k_density.x)
        val = pdf(k_density, x_int) .+ h
    else
        x_int = k_density.x
        val = k_density.density .+ h
    end
    chain_med = median(chain_vec)
    chain_mean = mean(chain_vec)
    min = minimum(k_density.density .+ h)
    q_int = (show_qi ? [qs[1], chain_med, qs[2]] : [chain_med])

    return par, hpdi, lower_hpd, upper_hpd, h, qs, k_density, x_int, val, chain_med,
        chain_mean, min, q_int
end

@recipe function f(
    p::RidgelinePlot;
    hpd_val = [0.05, 0.2],
    q = [0.1, 0.9],
    spacer = 0.5,
    _riser = 0.2,
    show_mean = true,
    show_median = true,
    show_qi = false,
    show_hpdi = true,
    fill_q = true,
    fill_hpd = false,
    ordered = false
)

    chn = p.args[1]
    par_names = p.args[2]

    for i in 1:length(par_names)
        par, hpdi, lower_hpd, upper_hpd, h, qs, k_density, x_int, val, chain_med, chain_mean,
            min, q_int = _compute_plot_data(i, chn, par_names; hpd_val = hpd_val, q = q,
            spacer = spacer, _riser = _riser, show_mean = show_mean, show_median = show_median,
            show_qi = show_qi, show_hpdi = show_hpdi, fill_q = fill_q, fill_hpd = fill_hpd,
            ordered = ordered)

        yticks --> (length(par_names) > 1 ?
            (_riser .+ ((1:length(par_names)) .- 1) .* spacer, string.(par)) : :default)
        yaxis --> (length(par_names) > 1 ? "Parameters" : "Density" )
        @series begin
            seriestype := :hline
            label := nothing
            linecolor := "#BBBBBB"
            linewidth --> 1.2
            [h]
        end
        @series begin
            seriestype := :path
            label := nothing
            fillrange --> min
            fillalpha --> 0.8
            x_int, val
        end
        @series begin
            seriestype := :path
            label := nothing
            linecolor --> "#000000"
            k_density.x, k_density.density .+ h
        end
        @series begin
            seriestype := :path
            label --> (show_mean ? (i == 1 ? "Mean" : nothing) : nothing)
            linecolor --> "dark red"
            linewidth --> (show_mean ? 1.2 : 0)
            [chain_mean, chain_mean], [min, min + pdf(k_density, chain_mean)]
        end
        @series begin
            seriestype := :path
            label --> (show_median ? (i == 1 ? "Median" : nothing) : nothing)
            linecolor --> "#000000"
            linewidth --> (show_median ? 1.2 : 0)
            [chain_med, chain_med], [min, min + pdf(k_density, chain_med)]
        end
        @series begin
            seriestype := :scatter
            label := (show_qi ? (i == 1 ? "Q$(q[1]), Q$(q[2])" : nothing) : nothing)
            markershape --> (show_qi ? :diamond : :circle)
            markercolor --> "#000000"
            markersize --> (show_qi ? 2 : 0)
            q_int, [h]
        end
        @series begin
            seriestype := :path
            label := nothing
            linecolor := "#000000"
            linewidth --> (show_qi ? 1.2 : 0)
            [qs[1], qs[2]], [h, h]
        end
        @series begin
            seriestype := :path
            label := (show_hpdi ? (i == 1 ? "$(Integer((1-hpdi[1])*100))% HPDI" : nothing)
                : nothing)
            linewidth --> (show_hpdi ? 2 : 0)
            seriesalpha --> 0.80
            linecolor --> :darkblue
            [lower_hpd[1][1], upper_hpd[1][1]], [h, h]
        end
    end
end

@recipe function f(
    p::ForestPlot;
    hpd_val = [0.05, 0.2],
    q = [0.1, 0.9],
    spacer = 0.5,
    _riser = 0.2,
    show_mean = true,
    show_median = true,
    show_qi = false,
    show_hpdi = true,
    fill_q = true,
    fill_hpd = false,
    ordered = false
)

    chn = p.args[1]
    par_names = p.args[2]

    for i in 1:length(par_names)
        par, hpdi, lower_hpd, upper_hpd, h, qs, k_density, x_int, val, chain_med, chain_mean,
            min, q_int = _compute_plot_data(i, chn, par_names; hpd_val = hpd_val, q = q,
            spacer = spacer, _riser = _riser, show_mean = show_mean, show_median = show_median,
            show_qi = show_qi, show_hpdi = show_hpdi, fill_q = fill_q, fill_hpd = fill_hpd,
            ordered = ordered)

        yticks --> (length(par_names) > 1 ?
            (_riser .+ ((1:length(par_names)) .- 1) .* spacer, string.(par)) : :default)
        yaxis --> (length(par_names) > 1 ? "Parameters" : "Density" )

        for j in 1:length(hpdi)
            @series begin
                seriestype := :path
                label := (show_hpdi ?
                    (i == 1 ? "$(Integer((1-hpdi[j])*100))% HPDI" : nothing) : nothing)
                linecolor --> j
                linewidth --> (show_hpdi ? 1.5*j : 0)
                seriesalpha --> 0.80
                [lower_hpd[j][1], upper_hpd[j][1]], [h, h]
            end
        end
        @series begin
            seriestype := :scatter
            label := (show_median ? (i == 1 ? "Median" : nothing) : nothing)
            markershape --> :diamond
            markercolor --> "#000000"
            markersize --> (show_median ? length(hpdi) : 0)
            [chain_med], [h]
        end
        @series begin
            seriestype := :scatter
            label := (show_mean ? (i == 1 ? "Mean" : nothing) : nothing)
            markershape --> :circle
            markercolor --> :gray
            markersize --> (show_mean ? length(hpdi) : 0)
            [chain_mean], [h]
        end
        @series begin
            seriestype := :scatter
            label := (show_qi ? (i == 1 ? "Q1 = $(q[1]), Q3 = $(q[2])" : nothing) : nothing)
            markershape --> (show_qi ? :diamond : :circle)
            markercolor --> "#000000"
            markersize --> (show_qi ? 2 : 0)
            q_int, [h]
        end
        @series begin
            seriestype := :path
            label := nothing
            linecolor := "#000000"
            linewidth --> (show_qi ? 1.2 : 0.0)
            [qs[1], qs[2]], [h, h]
        end
    end
end


text:
# Functions for working with missings in RLE values

# Using implementations matching those in Missings.jl for AbstractArray
Missings.allowmissing(x::RLEVector{T1,T2}) where {T1,T2} =
    convert(RLEVector{Union{T1,Missing},T2}, x)
Missings.disallowmissing(x::RLEVector{T1,T2}) where {T1,T2} =
    convert(RLEVector{nonmissingtype(T1),T2}, x)

# RLEVector methods on functions for arrays that may have missings
# see https://github.com/JuliaLang/julia/blob/c6da87ff4bc7a855e217856757ad3413cf6d1f79/base/missing.jl


text:
using Video4Linux
using ImageView
using Colors, ColorTypes

# Warning! no memory and device protection is implemented yet, therefore doing things out of order will cause julia to crash!

## NOTE: if your device does not support read try using Video4Linux.IO_METHOD_MMAP
set_io_method(Video4Linux.IO_METHOD_READ)

## open device
fid = open_device("/dev/video0")

## init_device(fd, force_format);
init_device(fid)

## start_capturing(fd);
start_capturing(fid)

## mainloop(fd, frame_count);
mainloop( fid, 1 )

## copy_buffer_bytes, copy the image buffer bytes to uint8 vector, the lenght will depend on the pixel format
imbuff = copy_buffer_bytes(640*480*2)
## copy frame data: raw + width + height + pixelformat
frame = Video4Linux.copy_buffer_frame()
#do some convertions on frame
img = Video4Linux.convert(Gray, frame)
imc = Video4Linux.convert(YCbCr, frame)
imshow(img)

## stop_capturing(fd);
stop_capturing(fid)

## uninit_device();
uninit_device(fid)

## close device
close_device(fid)

## kinect color image; manually set the driver to 640x480 YCbCr 4:2:2 for this test (can be done with qv4l2)

imYCbCr = convertUYVYtoYCbCr(imbuff, 640, 480)
#convert to rgb and display
imshow(convert.(RGB,imYCbCr))

## kinect 1 depth image, set kernel to depth with:
    #sudo rmmod gspca_kinect
    #sudo modprobe gspca_kinect depth_mode = 1

depthvec = convertY10BtoU16(imbuff[1:384000])

depthim = reshape(depthvec,(640,480))'
imshow(depthim)


text:
using Random

@testset "surrogate" begin
    N, M = 1024, 2
    X = zeros(Float64, N, M)
    for i in 1:M
        X[:,i] = map(n -> sin(2π*n/200), 1:N)
    end

    # Test 1D surrogate norm conservation
    Y = surrogate(X[:,1])
    normX = sum(map(t -> abs(t)^2, X[:,1]))
    normY = sum(map(t -> abs(t)^2, Y))
    @test normX ≈ normY rtol=1/N

    # Test 1D in-place surrogate norm conservation
    Y = copy(X[:,1]); surrogate!(Y)
    normX = sum(map(t -> abs(t)^2, X[:,1]))
    normY = sum(map(t -> abs(t)^2, Y))
    @test normX ≈ normY rtol=1/N

    # Test 2D surrogate norm conservation
    Y = surrogate(X)
    normX = sum(map(t -> abs(t)^2, X))
    normY = sum(map(t -> abs(t)^2, Y))
    @test normX ≈ normY rtol=1/N

    # Test 2D in-place surrogate norm conservation
    Y = copy(X); surrogate!(Y)
    normX = sum(map(t -> abs(t)^2, X))
    normY = sum(map(t -> abs(t)^2, Y))
    @test normX ≈ normY rtol=1/N
end


text:
# This file was generated, do not modify it. # hide
ẑ = 0.3ẑ₁ + 0.7ẑ₂;

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
module TwoPeriodAltruism

include("structures.jl")

end # module


text:
export kdownload

struct Auth
    username::String
    key::String
end

StructTypes.StructType(::Type{Auth}) = StructTypes.Struct()

function gen_auth_key()
    auth_file = joinpath(homedir(), ".kaggle", "kaggle.json")

    if isfile(auth_file)
        f = open(auth_file, "r")
        auth = JSON3.read(f, Auth)
        close(f)
    else
        auth = Auth(ENV["KAGGLE_USERNAME"], ENV["KAGGLE_KEY"])
    end

    auth_str = Base64.base64encode("$(auth.username):$(auth.key)")

    return "Basic $(auth_str)"
end

function gen_kaggle_url(dataset)
    return "https://www.kaggle.com/api/v1/datasets/download/$dataset"
end

function kdownload(dataset, localdir)
    url = gen_kaggle_url(dataset)
    filepath = joinpath(localdir, "$(replace(dataset, '/'=>'_')).zip")

    HTTP.open("GET", url, ["Authorization"=>gen_auth_key()]) do stream
        HTTP.startread(stream)
        eof(stream) && return

        total_bytes = tryparse(Int64, HTTP.header(stream, "Content-Length"))
        (total_bytes === nothing) && (total_bytes = NaN)
        println("Total: $total_bytes bytes")

        downloaded_bytes = progress = 0
        print("Downloaded:\e[s")
        Base.open(filepath, "w") do f
            while !eof(stream)
                downloaded_bytes += write(f, readavailable(stream))
                new_progress = 100downloaded_bytes ÷ total_bytes
                (new_progress > progress) && print("\e[u $downloaded_bytes bytes ($new_progress%)")
                progress = new_progress
            end
        end
    end

    return filepath
end


text:
@inline function Transducers.__foldl__(rf, val, A::UnionVector)
    @assert size(A.data) == size(A.typeid)
    for i in axes(A, 1)
        val = @return_if_reduced let i = i, val = val
            foldlargs(1, A.views...) do id, v
                Base.@_inline_meta
                if @inbounds A.typeid[i] == id
                    # input = @inbounds v[i]
                    input = GC.@preserve v unsafe_load(pointer(v, i))
                    Reduced(next(rf, val, input))
                else
                    id + 1
                end
            end :: Reduced |> unreduced
        end
    end
    return complete(rf, val)
end
# NOTE: Using `Reduced` and not `reduced` to exit `foldltupletype` but
# not `__foldl__` unless `next` returns a `Reduced` as well.  In this
# case (i.e., `next` returned a `Reduced`), the returned value of
# `foldltupletype` is doubly wrapped in `Reduced`.


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory
using MUMPSjInv
using Test

include("getDivGrad.jl")
A = getDivGrad(24, 23, 25)
A2 = getDivGrad(34, 32, 36)
n = size(A, 1)
n2 = size(A2, 1)
A = A + im * spdiagm(0 => rand(n))

nrhs = 10
rhs = randn(n, nrhs) + im*randn(n, nrhs)
rhs2 = randn(n2, nrhs)

println("Factorize complex matrix")
@time F1 = factorMUMPS(A,1)

println("Factorize real matrix")
@time F2 = factorMUMPS(A2,1)

println("Solve complex system")
@time x = applyMUMPS(F1,rhs)
err = zeros(nrhs)
for i=1:nrhs
    err[i] = norm(A * x[:, i] - rhs[:, i]) / norm(rhs[:, i])
end
@test maximum(err) < 1e-14

println("Solve real system")
@time x2 = applyMUMPS(F2, rhs2)
err = zeros(nrhs)
for i = 1:nrhs
    err[i] = norm(A2 * x2[:, i] - rhs2[:, i]) / norm(rhs2[:, i])
end
@test maximum(err) < 1e-14

println("Free memory")
destroyMUMPS(F1)
destroyMUMPS(F2)


text:
# Test for infinite Bi-Lanczos

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "Infbilanczos σ=0" begin
    nep=nep_gallery("qdep0");
    nept=SPMF_NEP([copy(nep.A[1]'), copy(nep.A[2]'), copy(nep.A[3]')], nep.fi)
    n=size(nep,1);

    m=40;
    λ,V,T = infbilanczos(Float64,nep,nept,maxit=m,neigs=3,σ=0,logger=displaylevel,
                         v=ones(Float64,n),u=ones(Float64,n),check_error_every=3,
                         tol=1e-7, errmeasure=ResidualErrmeasure(nep));

    # Produced with a different implementation
    Tstar=[  -1.665117675679600   5.780562035399026                   0                   0
             5.780562035399026  11.562308485001218 -18.839546184493731                   0
             0  18.839546184493734 -15.213756300995186   9.788512505128466
             0                   0   9.788512505128464  -0.120825360586847                                ]
    n0=min(size(Tstar,1),size(λ,1));
    @test opnorm(Tstar[1:n0,1:n0]-T[1:n0,1:n0])<1e-10

    thiserr=ones(m)*NaN
    for i=1:length(λ)
        thiserr[i]=norm(compute_Mlincomb(nep,λ[i],V[:,i]));
    end
    @test length(findall(thiserr .< 1e-7)) == 3

    @testset "Compute as many eigenpairs as possible (neigs=Inf)" begin
        m=30;
        λ,V,T = infbilanczos(Float64,nep,nept,maxit=m,neigs=Inf,σ=0,logger=displaylevel,
                             v=ones(Float64,n),u=ones(Float64,n),check_error_every=3,
                             tol=1e-7, errmeasure=ResidualErrmeasure(nep));
        verify_lambdas(3, nep, λ, V, 1e-6)
    end

    @testset "Errors thrown" begin
        @test_throws NEPCore.NoConvergenceException λ,V,T = infbilanczos(Float64,nep,nept,maxit=9,neigs=8,σ=0,logger=displaylevel,
                             v=ones(Float64,n),u=ones(Float64,n),check_error_every=3,
                             tol=1e-7, errmeasure=ResidualErrmeasure(nep));
    end
end

# Disabled to improve unit test performance
#@testset "Infbilanczos σ=$x" for x in (1.0, 1.0+0.1im)
#    m=30;
#    λ,V,T = infbilanczos(nep,nept,maxit=m,neigs=2,σ=x,tol=1e-7,
#                         v=ones(n),u=ones(n));
#    thiserr=ones(m)*NaN
#    for i=1:length(λ)
#        thiserr[i]=norm(compute_Mlincomb(nep,λ[i],V[:,i]));
#    end
#    @test length(find(thiserr .< 1e-7)) >= 2
#end

##for i=1:m
## semilogy(1:m, err[1:m,i], color="red", linewidth=2.0, linestyle="--")
##end
##
#errormeasure=default_errmeasure(nep);
#for i=1:length(λ)
# @info "Eigenvalue = $(λ[i]), residual = $(errormeasure(λ[i],Q[:,i]))"
#end
#


text:
function cnf_and(exprs::BoolExprVec)::And
    cnf_exprs = map(cnf, exprs)
    And(flatmap(children, cnf_exprs))
end

function cnf_or(exprs::BoolExprVec)::And
    cnf_exprs = map(cnf, exprs)
    # Each item in cnf_exprs is in CNF form, that is
    # * A variable
    # * A conjunction (And) of disjunctions (Or) of variables
    prod = cartprod(map(children, cnf_exprs))
    And(map(a -> Or(children(a)), prod))
end

cnf_not(parent::Not, expr::BoolExpr) = parent
cnf_not(parent::Not, expr::And) = cnf(Or(negate(expr.exprs))) # DeMorgan
cnf_not(parent::Not, expr::Or)  = cnf(And(negate(expr.exprs))) # DeMorgan
cnf_not(parent::Not, expr::Not) = cnf(expr.expr) # Double negation
cnf_not(parent::Not, expr::True) = f
cnf_not(parent::Not, expr::False) = t

cnf(expr::Not) = cnf_not(expr, expr.expr)
cnf(expr::And) = cnf_and(expr.exprs)
cnf(expr::Or)  = cnf_or(expr.exprs)
cnf(expr::BoolExpr) = expr


text:
using DataAPI

using CategoricalArrays: CategoricalVector, CategoricalArray, CategoricalPool

compress_then_write(b::CategoricalVector{T,IntType}, io) where {T, IntType<:Integer} = begin
    #println("abc")
    compress_refs = compress_then_write(b.refs, io)
    compress_poolindex = compress_then_write(DataAPI.levels(b), io)

    (
        type = CategoricalVector,
        refs = compress_refs,
        poolindex = compress_poolindex,
        ordered = b.pool.ordered,
    )
end

column_loader(b::Type{CategoricalVector}, io, metadata) = begin
    refs_meta = metadata.refs
    pi_meta = metadata.poolindex
    ref = column_loader(refs_meta.type, io, refs_meta)
    poolindex = column_loader(pi_meta.type, io, pi_meta)
    CategoricalArray{pi_meta.type,1}(
        ref,
        CategoricalPool{eltype(poolindex),eltype(ref)}(Array(poolindex), metadata.ordered),
    )
end

if false
    # # compress_then_write-cate.jl
    # using Revise
    # using JDF
    # using DataFrames, Debugger
    # b = categorical(["a", "b", "a", "c"])
    # io = open("c:/data/test.io", "w")
    # metadata = compress_then_write(a, io)
    # close(io)
    #
    # io = open("c:/data/test.io", "r")
    # aa = column_loader(CategoricalVector, io, metadata)
    # close(io)
    #
    # df = DataFrame(a = a)
    #
    # savejdf("c:/data/pls_del.jdf", df)
    #
    # loadjdf("c:/data/pls_del.jdf")
end


text:
module CompressedIndices

import Base: eltype, length, endof, size, checkbounds, eachindex, getindex, start, next, done, push!, append!, intersect, intersect!, +, -
export Indices

AbstractIndices = Union{Int, UnitRange{Int}, Vector{Int}}

immutable Indices <: AbstractVector{Int}
    vectors::Vector{Int}
    ranges::Vector{Int}
    length::Vector{Int}
    isrange::BitVector
end

eltype(I::Indices) = Int
length(I::Indices) = sum(I.length)
endof(I::Indices) = length(I)
size(I::Indices) = (length(I),)
eachindex(I::Indices) = 1:length(I)

start(I::Indices) = (1, 0, 0, 1)
function next(I::Indices, state::NTuple{4, Int})
    @inbounds l = state[1]
    @inbounds vi = state[2]
    @inbounds ri = state[3]
    @inbounds li = state[4]
    @inbounds isrange = I.isrange[li]

    item = 0
    if isrange
        @inbounds item = I.ranges[ri + 1] + l - 1
    else
        @inbounds item = I.vectors[vi + l]
    end

    l += 1
    @inbounds l_li = I.length[li]
    if l > l_li
        if isrange
            ri += 1
        else
            vi += l - 1
        end

        l = 1
        li += 1
    end

    return item, (l, vi, ri, li)
end

done(I::Indices, state::NTuple{4, Int}) = (@inbounds s = state[4] ; s > length(I.length))

function push!(I::Indices, index::Int)
    if length(I.length) == 0
        push!(I.ranges, index)
        push!(I.length, 1)
        push!(I.isrange, true)

        return I
    end

    li = endof(I.length)
    @inbounds isrange = I.isrange[li]

    if isrange
        @inbounds rs = I.ranges[end]
        @inbounds rl = I.length[li]
        if index == rs + rl
            @inbounds I.length[li] += 1
        else
            push!(I.vectors, index)
            push!(I.length, 1)
            push!(I.isrange, false)
        end
    else
        @inbounds vl = I.vectors[end]
        if index == vl + 1
            pop!(I.vectors)
            @inbounds l_li = I.length[li]
            if l_li > 1
                @inbounds I.length[li] -= 1
                push!(I.length, 2)
                push!(I.isrange, true)
            else
                @inbounds I.length[li] += 1
                @inbounds I.isrange[li] = true
            end
            push!(I.ranges, vl)
        else
            push!(I.vectors, index)
            @inbounds I.length[li] += 1
        end
    end

    return I
end

function getindex(I::Indices, index::Int)
    n = 0
    ri = 0
    vi = 0

    for (i, l) in enumerate(I.length)
        @inbounds isrange = I.isrange[i]
        if n + l >= index
            if isrange
                return I.ranges[ri + 1] + (index - n - 1)
            else
                return I.vectors[vi + index - n]
            end
        end

        n += l
        if isrange
            ri += 1
        else
            vi += l
        end
    end

    return 0
end

@generated function Indices(ind::AbstractIndices...)
    length(ind) == 0 && return :(Indices(Vector{Int}(), Vector{Int}(), Vector{Int}(), BitVector()))
    num_ranges = count(x -> x == UnitRange{Int}, ind)

    isrange = BitVector([ind[1] == UnitRange{Int}])

    for i in 2:endof(ind)
        b = ind[i] == UnitRange{Int}
        if last(isrange) || b
            push!(isrange, b)
        end
    end

    ex = Expr[
        :(vectors = Vector{Int}())
        :(ranges = Vector{Int}($num_ranges))
        :(len = Vector{Int}($(length(isrange))))
        :(isrange = BitVector($isrange))
    ]

    ri = 0
    oi = 0
    last_isrange = true
    for (i, i_type) in enumerate(ind)
        if i_type == UnitRange{Int}
            ri += 1
            push!(ex, :(ranges[$ri] = first(ind[$i])))

            oi += 1
            push!(ex, :(len[$oi] = length(ind[$i])))

            last_isrange = true
        elseif i_type == Int
            if last_isrange
                oi += 1
                push!(ex, :(len[$oi] = 1))
                last_isrange = false
            else
                push!(ex, :(len[$oi] += 1))
            end

            push!(ex, :(push!(vectors, ind[$i])))
        elseif i_type == Vector{Int}
            if last_isrange
                oi += 1
                push!(ex, :(len[$oi] = length(ind[$i])))
                last_isrange = false
            else
                push!(ex, :(len[$oi] += length(ind[$i])))
            end

            push!(ex, :(append!(vectors, ind[$i])))
        else
            throw(ArgumentError("Invalid types in constructor Indices()"))
        end
    end

    push!(ex, :(Indices(vectors, ranges, len, isrange)))

    return Expr(:block, ex...)
end

end


text:
"""
Algorithm for solving Physics-Informed Neural Networks problems.

Arguments:
* `dx` is the discretization of the grid
* `chain` is a Flux.jl chain with a d-dimensional input and a 1-dimensional output,
* `init_params` is the initial parameter of the neural network,
* `phi` is a trial solution,
* `autodiff` is a boolean variable that determines whether to use automatic, differentiation (not supported while) or numerical,
* `derivative` is method that calculates the derivative,
* `strategy` determines which training strategy will be used.
"""

struct PhysicsInformedNN{D,C,P,PH,DER,T,K}
  dx::D
  chain::C
  initθ::P
  phi::PH
  autodiff::Bool
  derivative::DER
  strategy::T
  kwargs::K
end

function PhysicsInformedNN(dx,
                           chain,
                           init_params = nothing;
                           _phi = nothing,
                           autodiff=false,
                           _derivative = nothing,
                           strategy = GridTraining(),
                           kwargs...)
    if init_params === nothing
        initθ = DiffEqFlux.initial_params(chain)
    else
        initθ = init_params
    end

    if _phi == nothing
        phi = get_phi(chain)
    else
        phi = _phi
    end

    if _derivative == nothing
        derivative = get_derivative(autodiff)
    else
        derivative = _derivative
    end

    PhysicsInformedNN(dx,chain,initθ,phi,autodiff,derivative, strategy, kwargs)
end

abstract type TrainingStrategies  end

struct GridTraining <: TrainingStrategies end

struct StochasticTraining <:TrainingStrategies
    include_frac::Float64
end
function StochasticTraining(;include_frac=0.75)
    StochasticTraining(include_frac)
end

"""
Create dictionary: variable => unique number for variable

# Example 1

Dict{Symbol,Int64} with 3 entries:
  :y => 2
  :t => 3
  :x => 1

# Example 2

 Dict{Symbol,Int64} with 2 entries:
  :u1 => 1
  :u2 => 2
"""
get_dict_vars(vars) = Dict( [Symbol(v) .=> i for (i,v) in enumerate(vars)])

# Wrapper for _transform_derivative
function transform_derivative(ex,dict_indvars,dict_depvars)
    if ex isa Expr
        ex.args = _transform_derivative(ex.args,dict_indvars,dict_depvars)
    end
    return ex
end

function get_ε(dim, der_num)
    epsilon = cbrt(eps(Float32))
    ε = zeros(Float32, dim)
    ε[der_num] = epsilon
    ε
end

"""
Transform the derivative expression to inner representation

# Examples

1. First compute the derivative of function 'u(x,y)' with respect to x.

Take expressions in the form: `derivative(u(x,y,θ), x)` to `derivative(u, [x, y], εs, order, θ)`,

where
 u_d - derived function
 x,y - coordinates of point
 εs - epsilon mask
 order - order of derivative
 θ - parameters of neural network
"""
function _transform_derivative(_args,dict_indvars,dict_depvars)
    dim = length(dict_indvars)
    εs = [get_ε(dim,d) for d in 1:dim]
    for (i,e) in enumerate(_args)
        if !(e isa Expr)
            if e == :derivative && _args[end] != :θ
                derivative_variables = Symbol[]
                order = 0
                while (_args[1] == :derivative)
                    order += 1
                    push!(derivative_variables, _args[end])
                    _args = _args[2].args
                end
                depvar = _args[1]
                indvars = _args[2:end-1]
                undv = [dict_indvars[d_p] for d_p  in derivative_variables]
                εs_dnv = [εs[d] for d in undv]
                _args = [:derivative, Symbol(:($depvar),:_d), :([$(indvars...)]), εs_dnv, order, :θ]
                break
            end
        else
            _args[i].args = _transform_derivative(_args[i].args,dict_indvars,dict_depvars)
        end
    end
    return _args
end

"""
Parse ModelingToolkit equation form to the inner representation.

Example:

1)  1-D ODE: Dt(u(t,θ)) ~ t +1

    Take expressions in the form: 'Equation(derivative(u(t, θ), t), t + 1)' to 'derivative(u_d, [t], [[ε]], 1, θ) - (t + 1)'

2)  2-D PDE: Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)

    Take expressions in the form:
     Equation(derivative(derivative(u(x, y, θ), x), x) + derivative(derivative(u(x, y, θ), y), y), -(sin(πx)) * sin(πy))
    to
     (derivative(u_d, [x, y], [[ε,0],[ε,0]], 2, θ) + derivative(u_d, [x, y], [[0,ε],[0,ε]], 2, θ)) - -(sin(πx)) * sin(πy)

3)  System of PDEs: [Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0,
                    Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0]

    Take expressions in the form:
    2-element Array{Equation,1}:
        Equation(derivative(u1(x, y, θ), x) + 4 * derivative(u2(x, y, θ), y), ModelingToolkit.Constant(0))
        Equation(derivative(u2(x, y, θ), x) + 9 * derivative(u1(x, y, θ), y), ModelingToolkit.Constant(0))
    to
      [(derivative(u1_d, [x, y], [[ε,0]], 1, θ) + 4 * derivative(u2_d, [x, y], [[0,ε]], 1, θ)) - 0,
       (derivative(u2_d, [x, y], [[ε,0]], 1, θ) + 9 * derivative(u1_d, [x, y], [[0,ε]], 1, θ)) - 0]
"""
function parse_equation(eq,dict_indvars,dict_depvars)
    left_expr = transform_derivative(toexpr(eq.lhs),
                                     dict_indvars,dict_depvars)
    right_expr = transform_derivative(toexpr(eq.rhs),
                                     dict_indvars,dict_depvars)
    loss_func = :($left_expr - $right_expr)
end

"""
Build a loss function for a PDE or a boundary condition

# Examples: System of PDEs:

Take expressions in the form:

[Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0,
 Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0]

to

:((phi, cord, θ, derivative)->begin
          #= none:35 =#
          #= none:35 =#
          begin
              begin
                  u1 = ((x, y, θ)->begin
                              #= none:16 =#
                              (phi([x, y], θ))[1]
                          end)
                  u2 = ((x, y, θ)->begin
                              #= none:16 =#
                              (phi([x, y], θ))[2]
                          end)
                  u1_d = ((cord, θ)->begin
                              #= none:20 =#
                              (phi(cord, θ))[1]
                          end)
                  u2_d = ((cord, θ)->begin
                              #= none:20 =#
                              (phi(cord, θ))[2]
                          end)
              end
              let (x, y) = (cord[1], cord[2])
                  [(derivative(u1_d, [x, y], [[ε, 0]], 1, θ) + 4 * derivative(u2_d, [x, y], [[0, ε]], 1, θ)) - 0,
                   (derivative(u2_d, [x, y], [[ε, 0]], 1, θ) + 9 * derivative(u1_d, [x, y], [[0, ε]], 1, θ)) - 0]
              end
          end
      end)
"""
function build_loss_function(eqs,_indvars,_depvars)
    # dictionaries: variable -> unique number
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return build_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars)
end

function build_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars)
    if !(eqs isa Array)
        eqs = [eqs]
    end
    loss_functions= []
    for eq in eqs
        loss_function = parse_equation(eq,dict_indvars,dict_depvars)
        push!(loss_functions,loss_function)
    end

    vars = :(phi, cord, θ, derivative)
    ex = Expr(:block)
    us = []
    for v in depvars
        var_num = dict_depvars[v]
        push!(us,:($v = ($(indvars...), θ) -> phi([$(indvars...)],θ)[$var_num]))
    end
    for v in depvars
        var_num = dict_depvars[v]
        push!(us,:($(Symbol(:($v),:_d)) = (cord, θ) -> phi(cord,θ)[$var_num]))
    end

    u_ex = build_expr(:block, us)
    push!(ex.args,  u_ex)

    indvars_ex = [:($:cord[$i]) for (i, u) ∈ enumerate(indvars)]
    arg_pairs_indvars = indvars,indvars_ex

    left_arg_pairs, right_arg_pairs = arg_pairs_indvars
    vars_eq = Expr(:(=), build_expr(:tuple, left_arg_pairs), build_expr(:tuple, right_arg_pairs))
    let_ex = Expr(:let, vars_eq, build_expr(:vect, loss_functions))
    push!(ex.args,  let_ex)

    return :(($vars) -> begin $ex end)
end

# Get arguments from boundary condition functions
function get_bc_argument(bcs,indvars,depvars,dict_indvars,dict_depvars)
    bc_args = []
    for _bc in bcs
        _bc isa Array && error("boundary conditions must be represented as a one-dimensional array")
        left_expr = transform_derivative(toexpr(_bc.lhs),
                                         dict_indvars,dict_depvars)
        bc_arg = if (left_expr.args[1] == :derivative)
            left_expr.args[3].args
        else
            left_expr.args[2:end-1]
        end
        push!(bc_args, bc_arg)
    end
    return bc_args
end

function generate_training_sets(domains,dx,bcs,_indvars,_depvars)
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return generate_training_sets(domains,dx,bcs,indvars,depvars,dict_indvars,dict_depvars)
end
# Generate training set in the domain and on the boundary
function generate_training_sets(domains,dx,bcs,indvars,depvars,dict_indvars,dict_depvars)
    if dx isa Array
        dxs = dx
    else
        dxs = fill(dx,length(domains))
    end

    bound_args = get_bc_argument(bcs,indvars,depvars,dict_indvars,dict_depvars)

    spans = [d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)]
    dict_var_span = Dict([Symbol(d.variables) => d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)])

    train_set = map(points -> collect(points), Iterators.product(spans...))

    train_bound_set = map(bound_args) do bt
        span = map(b -> get(dict_var_span, b, [b]), bt)
        _set = map(points -> collect(points), Iterators.product(span...))
    end

    flat_train_bound_set = map(x -> x, Iterators.flatten(train_bound_set))
    train_domain_set =  setdiff(train_set, flat_train_bound_set)

    [train_domain_set,train_bound_set,train_set]
end

function get_phi(chain)
    # The phi trial solution
    if chain isa FastChain
        phi = (x,θ) -> chain(adapt(typeof(θ),x),θ)
    else
        _,re  = Flux.destructure(chain)
        phi = (x,θ) -> re(θ)(adapt(typeof(θ),x))
    end
    phi
end

# the method to calculate the derivative
function get_derivative(autodiff)
    epsilon = 2*cbrt(eps(Float32))
    if autodiff # automatic differentiation (not implemented yet)
        error("automatic differentiation is not implemented yet)")
    else # numerical differentiation
        derivative = (u,x,εs,order,θ) ->
        begin
            ε = εs[order]
            if order > 1
                return (derivative(u,x+ε,εs,order-1,θ)
                      - derivative(u,x-ε,εs,order-1,θ))/epsilon
            else
                return (u(x+ε,θ) - u(x-ε,θ))/epsilon
            end
        end
    end
    derivative
end

function get_loss_function(loss_functions, train_sets, phi, derivative, strategy)

    # norm coefficient for loss function
    τ = sum(length(train_set) for train_set in train_sets)

    function inner_loss(loss_functions,phi,x,θ,derivative)
        sum(loss_functions(phi, x, θ, derivative))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        train_sets = [train_sets]
    end

    if strategy isa StochasticTraining
        include_frac = strategy.include_frac
        count_elements = []
        sets_size = []
        for j in 1:length(train_sets)
            size_set = size(train_sets[j])[1]
            count_element = convert(Int64,round(include_frac*size_set, digits=0))
            if count_element <= 2
                count_element = size_set
            end
            push!(sets_size,size_set)
            push!(count_elements,count_element)
        end
        loss = (θ) -> begin
            total = 0.
            for (j,l) in enumerate(loss_functions)
                size_set = sets_size[j]
                for i in 1:count_elements[j]
                    index = rand(1:size_set)
                    total += inner_loss(l,phi,train_sets[j][index],θ,derivative)^2
                end
            end
            return (1.0f0/τ) * total
        end

    elseif strategy isa GridTraining
        f = (loss,train_set,θ) -> sum(abs2,[inner_loss(loss,phi,x,θ,derivative) for x in train_set])
        loss = (θ) ->  (1.0f0/τ) * sum(f(loss_function,train_set,θ) for (loss_function,train_set) in zip(loss_functions,train_sets))
    end
    return loss
end

# Convert a PDE problem into an OptimizationProblem
function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eq
    bcs = pde_system.bcs

    domains = pde_system.domain
    # dimensionality of equation
    dim = length(domains)

    depvars = [nameof(value(d)) for d in pde_system.depvars]
    indvars = [nameof(value(i)) for i in pde_system.indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)

    dx = discretization.dx
    chain = discretization.chain
    initθ = discretization.initθ
    phi = discretization.phi
    autodiff = discretization.autodiff
    derivative = discretization.derivative
    autodiff == true && error("Automatic differentiation is not support yet")
    strategy = discretization.strategy

    # length(domains) != chain.layers[1].in && error("the input of chain should equal the length of domains according to the dimensionality of the task")
    # length(depvars) != chain.layers[end].out && error("the output of chain should equal the number of variables")

    train_sets = generate_training_sets(domains,dx,bcs,
                                        indvars,depvars,
                                        dict_indvars,dict_depvars)

    # the points in the domain and on the boundary
    train_domain_set,train_bound_set,train_set = train_sets

    expr_pde_loss_function = build_loss_function(eqs,indvars,depvars,
                                                 dict_indvars,dict_depvars)
    expr_bc_loss_functions = [build_loss_function(bc,indvars,depvars,
                                                  dict_indvars,dict_depvars) for bc in bcs]

    pde_loss_function = get_loss_function(eval(expr_pde_loss_function),
                                          train_domain_set,
                                          phi,
                                          derivative,
                                          strategy)
    bc_loss_function = get_loss_function(eval.(expr_bc_loss_functions),
                                         train_bound_set,
                                         phi,
                                         derivative,
                                         strategy)

    function loss_function(θ,p)
        return pde_loss_function(θ) + bc_loss_function(θ)
    end

    f = OptimizationFunction(loss_function, GalacticOptim.AutoZygote())
    GalacticOptim.OptimizationProblem(f, initθ)
end


text:
using Test

using StaticArrays
using iLQGames:
    ControlSystem,
    Lorenz3D,
    dx,
    integrate,
    linearize

l3 = Lorenz3D{0.1}(10., 28., 8//3)

# test some linearization
@testset "control_system.jl" begin
    @test @inferred_with_info integrate(l3, @SVector([0., 0., 0.]), @SVector([0., 1.]), 0., 0.1)
    @test @inferred_with_info dx(l3, @SVector([2., 1., 7.]), @SVector([0., 1.]), 0.)
    @test @inferred_with_info linearize(l3, @SVector([2., 1., 7.]), @SVector([0., 1.]), 0.)
end;


text:
include("../src/RationalFunctionFields.jl")

using .RationalFunctionFields
using .RationalFunctionFields: load_generators, new_generating_set,
             field_generators, compute_groebner!, simplify_generators!,
             RationalFunctionField, field_generators

Singular = RationalFunctionFields.Singular
Logging = RationalFunctionFields.Logging

logger = Logging.SimpleLogger(stderr, Logging.Debug)
Logging.global_logger(logger)

set = load_generators("../data/LV_gens.txt")

println(set)
FF = RationalFunctionField(set)
compute_groebner!(FF)
simplify_generators!(FF)

println( "\ngenerators ", field_generators(FF) )


text:
using PDMats
using Test
using LinearAlgebra: LinearAlgebra

_randPDMat(T, n) = (X = randn(T, n, n); PDMat(X * X' + LinearAlgebra.I))
_randPDiagMat(T, n) = PDiagMat(rand(T, n))
_randScalMat(T, n) = ScalMat(n, rand(T))

function _pd_compare(A::AbstractPDMat, B::AbstractPDMat)
    @test dim(A) == dim(B)
    @test Matrix(A) ≈ Matrix(B)
    @test cholesky(A).L ≈ cholesky(B).L
    @test cholesky(A).U ≈ cholesky(B).U
end

function _pd_kron_compare(A::AbstractPDMat, B::AbstractPDMat)
    PDAkB_kron = kron(A, B)
    PDAkB_dense = PDMat( kron( Matrix(A), Matrix(B) ) )
    _pd_compare(PDAkB_kron, PDAkB_dense)
end

n = 4
m = 7

@testset "Kronecker product" begin
    for T in [Float32, Float64]
        _pd_kron_compare( _randPDMat(T, n),    _randPDMat(T, m) )
        _pd_kron_compare( _randPDiagMat(T, n), _randPDiagMat(T, m) )
        _pd_kron_compare( _randScalMat(T, n),  _randScalMat(T, m) )
        _pd_kron_compare( _randPDMat(T, n),    _randPDiagMat(T, m) )
        _pd_kron_compare( _randPDiagMat(T, m), _randPDMat(T, n) )
        _pd_kron_compare( _randPDMat(T, n),    _randScalMat(T, m) )
        _pd_kron_compare( _randScalMat(T, m),  _randPDMat(T, n) )
        _pd_kron_compare( _randPDiagMat(T, n), _randScalMat(T, m) )
        _pd_kron_compare( _randScalMat(T, m),  _randPDiagMat(T, n) )
    end
end


text:
using SpikingNetworks
using Base.Test

parseParameters"""
τₘ = 20
τₑ = 5
τᵢ = 10
Vₜ = -50
Vᵣ = -60
Eₗ = -49
N=4000;M=4*N÷5
wₑ=spzeros(N,N)
wᵢ=spzeros(N,N)
wₑ[1:N,1:M]=sprand(N,M,0.02)*(60*0.27/10)
wᵢ[1:N,M+1:N]=sprand(N,N-M,0.02)*(-20*4.5/10)"""

parseEquations"""
dv/dt=(gₑ+gᵢ-(v-Eₗ))/τₘ
dgₑ/dt=-gₑ/τₑ
dgᵢ/dt=-gᵢ/τᵢ"""

net=init"""
N=N
v=Vᵣ+rand(N)*(Vₜ-Vᵣ)
gₑ=zeros(N)
gᵢ=zeros(N)
wₑ=wₑ
wᵢ=wᵢ"""

spike"v.>Vₜ"
reset"""
v[i]=Vᵣ"""

synapse"""
gₑ+=wₑ
gᵢ+=wᵢ"""

net.dt=0.1
@time run!(net,100/net.dt)


text:
# Autogenerated wrapper script for lrslib_jll for armv6l-linux-gnueabihf
export liblrs, liblrsnash, lrs, lrsnash

using GMP_jll
JLLWrappers.@generate_wrapper_header("lrslib")
JLLWrappers.@declare_library_product(liblrs, "liblrs.so.1")
JLLWrappers.@declare_library_product(liblrsnash, "liblrsnash.so")
JLLWrappers.@declare_executable_product(lrs)
JLLWrappers.@declare_executable_product(lrsnash)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        liblrs,
        "lib/liblrs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblrsnash,
        "lib/liblrsnash.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        lrs,
        "bin/lrs",
    )

    JLLWrappers.@init_executable_product(
        lrsnash,
        "bin/lrsnash",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module KhepriPOVRay
using KhepriBase

# functions that need specialization
include(khepribase_interface_file())
include("POVRay.jl")

function __init__()
  set_material(povray, material_basic, povray_neutral)
  set_material(povray, material_metal, povray_metal)
  set_material(povray, material_glass, povray_glass)
  set_material(povray, material_grass, povray_grass)
  set_material(povray, material_wood, povray_wood)
  set_material(povray, material_concrete, povray_material("Concrete", gray=0.3))
  set_material(povray, material_plaster, povray_material("InteriorWall70", gray=0.7))

  add_current_backend(povray)
end
end


text:
# This mechanism allows us to import a package list.
# Source: https://discourse.julialang.org/t/building-a-dockerfile-with-packages/37272/2

using Pkg
pkg"add Colors"
pkg"add Conda"
pkg"add CSV"
pkg"add CSVFiles"
pkg"add DataFrames"
pkg"add DelimitedFiles"
pkg"add Distributions"
pkg"add FileIO"
pkg"add GMT"
pkg"add GR"
pkg"add GtkReactive"
pkg"add IJulia"
pkg"add ImageFiltering"
pkg"add ImageInTerminal"
pkg"add ImageFeatures"
pkg"add Images"
pkg"add ImageView"
pkg"add Interact"
pkg"add LaTeXStrings"
pkg"add LinearAlgebra"
pkg"add OSQP"
pkg"add Plotly"
pkg"add Plots"
pkg"add ProgressBars"
pkg"add PyCall"
pkg"add PyPlot"
pkg"add Random"
pkg"add Roots"
pkg"add SparseArrays"
pkg"add Statistics"
pkg"add SymEngine"
pkg"add Symbolics"
pkg"add SymbolicUtils"
pkg"add TestImages"
pkg"add WebIO"
pkg"add https://github.com/VMLS-book/VMLS.jl"
pkg"add https://github.com/korsbo/Latexify.jl"

max_stars_repo_name: IllumiDesk/umich-stacks

max_issues_repo_name: IllumiDesk/umich-stacks

max_forks_repo_name: IllumiDesk/umich-stacks