text:
using Test
using st16SseqJuliaTools
using Suppressor
using DataStructures

@testset "uc2jc" begin
    @testset "Case: C" begin
        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n959\t959\n959\t959\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959\t*\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;AAAAA;AAAAA\t*\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;AAAAA;AAAAA\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n959\t959\n959\t959\n"
    end

    @testset "Case: H" begin
        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n975\t838\n975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=6629\t838\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975\t838\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;cacac;acaac\t838;acac;c;a;ca\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;cacac;acaac\t838;acac;c;a;ca\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975\t838\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=6629\t838\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n975\t838\n975\t838\n975\t838\n975\t838\n"
    end

    @testset "Case: MIXED" begin

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "975\t838\n959\t959\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n975\t838\n959\t959\n975\t838\n959\t959\n975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n959\t959\n959\t959\n959\t959\n975\t838\n975\t838\n975\t838\n975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "S\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "S\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "S\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "S\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "S\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n959\t959\n959\t959\n959\t959\n975\t838\n975\t838\n975\t838\n975\t838\n"

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "X\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "A\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "C\t0\t10115\t*\t*\t*\t*\t*\t959;size=10115\t*\n")
        print(tmp_fo, "G\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "S\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "A\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        print(tmp_fo, "A\t6\t240\t*\t*\t*\t*\t*\t952;size=5124\t*\n")
        print(tmp_fo, "H\t2\t240\t99.6\t+\t0\t0\t240M\t975;size=2856\t838;size=6629\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == "959\t959\n959\t959\n959\t959\n959\t959\n975\t838\n975\t838\n975\t838\n"
    end

    @testset "Case: ANOMALY" begin
        tmp_file, tmp_fo = mktemp()
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == ""

        tmp_file, tmp_fo = mktemp()
        print(tmp_fo, "aCACACACACACA\n")
        close(tmp_fo)
        output = @capture_out begin
            uc2jc(tmp_file)
        end
        @test output == ""

        tmp_file = "AAAAA"
        @test_throws SystemError uc2jc(tmp_file)
    end

end

@testset "mergejc" begin
    @testset "Case: MERGE" begin
        output = ""
        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "4956\t14932\n14980\t14932\n13426\t14932\n")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "14932\t14932\n")
        close(tmp_fo2)
        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '0'
        @test output == "4956\t14932\n14980\t14932\n13426\t14932\n"

        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "14931\t14932\n")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '2'
        @test output == "4956\t14932\n"

        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "14931\t14932\n14932\t14933\n")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '0'
        @test output == "4956\t14932\n14980\t14933\n13426\t14933\n"

        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "14931\t14932\n14932\t14933\n14930\t14938\n")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '0'
        @test output == "4956\t14932\n14980\t14933\n13426\t14933\n"

        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n13450\t14955\n")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "14931\t14932\n14932\t14933\n14930\t14938\n")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '1'
        @test output == "4956\t14932\n14980\t14933\n13426\t14933\n"
    end

    @testset "Case: ANOMALY" begin
        output = ""
        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '0'
        @test output == ""

        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n13450\t14955\n")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '4'
        @test output == ""

        tmp_file1, tmp_fo1 = mktemp()
        print(tmp_fo1, "")
        close(tmp_fo1)
        tmp_file2, tmp_fo2 = mktemp()
        print(tmp_fo2, "14931\t14932\n")
        close(tmp_fo2)

        err = @capture_err begin
            output = @capture_out begin
                mergejc(tmp_file1, tmp_file2)
            end
        end
        @test split(err, "\n")[3][1] == '0'
        @test output == ""

        tmp_file1 = "NO_file"
        tmp_file2 = "NO_file2"
        @suppress_err begin
            @test_throws SystemError mergejc(tmp_file1, tmp_file2)
        end
    end
end

@testset "mark_proper_size_clusters" begin
    test_dict = OrderedDict{String,Array{String,1}}()
    test_dict["Clus1"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus2"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus3"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus4"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus5"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    true_dict = OrderedDict{String,Bool}("Clus1" => 0,"Clus2" => 0,"Clus3" => 0,"Clus4" => 0,"Clus5" => 0)
    @test mark_proper_size_clusters(test_dict, 0.5) == true_dict

    test_dict = OrderedDict{String,Array{String,1}}()
    test_dict["Clus1"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus2"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus3"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus4"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus5"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    true_dict = OrderedDict{String,Bool}("Clus1" => 1,"Clus2" => 1,"Clus3" => 1,"Clus4" => 1,"Clus5" => 1)
    @test mark_proper_size_clusters(test_dict, 0.1) == true_dict

    test_dict = OrderedDict{String,Array{String,1}}()
    test_dict["Clus1"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus2"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus3"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus4"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus5"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    true_dict = OrderedDict{String,Bool}("Clus1" => 1,"Clus2" => 1,"Clus3" => 1,"Clus4" => 1,"Clus5" => 1)
    @test mark_proper_size_clusters(test_dict, 0.2) == true_dict

    test_dict = OrderedDict{String,Array{String,1}}()
    test_dict["Clus1"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus2"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus3"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus4"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    test_dict["Clus5"] = ["ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB", "ABCDABCDAB"]
    true_dict = OrderedDict{String,Bool}("Clus1" => 1,"Clus2" => 0,"Clus3" => 1,"Clus4" => 0,"Clus5" => 1)
    @test mark_proper_size_clusters(test_dict, 0.2) == true_dict
end

@testset "read_jc_clusters" begin
    tmp_file1, tmp_fo1 = mktemp()
    print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n13450\t14955\n")
    close(tmp_fo1)
    true_dict = OrderedDict("14932" => ["14980", "13426"],"14931" => ["4956"],"14955" => ["13450"])
    @test read_jc_clusters(tmp_file1) == true_dict

    tmp_file1, tmp_fo1 = mktemp()
    print(tmp_fo1, "4956\t14931\n14980\t14932\n13426\t14932\n13450\t14955\n")
    close(tmp_fo1)
    true_dict = OrderedDict("14932" => ["14980", "13426"],"14931" => ["4956"],"14955" => ["13450"])
    @test read_jc_clusters(tmp_file1) == true_dict

    tmp_file1, tmp_fo1 = mktemp()
    print(tmp_fo1, "4956\t14931\n14980\t14932\n13450\t14955\n")
    close(tmp_fo1)
    true_dict = OrderedDict("14932" => ["14980"],"14931" => ["4956"],"14955" => ["13450"])
    @test read_jc_clusters(tmp_file1) == true_dict

    tmp_file1, tmp_fo1 = mktemp()
    print(tmp_fo1, "4956\t14931\n4957\t14931\n4958\t14931\n4959\t14931\n")
    close(tmp_fo1)
    true_dict = OrderedDict("14931" => ["4956", "4957", "4958", "4959"])
    @test read_jc_clusters(tmp_file1) == true_dict

    tmp_file1, tmp_fo1 = mktemp()
    print(tmp_fo1, "4956\t14931\n4957\t14931\n4958\t14931\n4959\t14931\n4956\t14932\n4957\t14932\n4958\t14932\n")
    close(tmp_fo1)
    true_dict = OrderedDict("14931" => ["4956", "4957", "4958", "4959"], "14932" => ["4956", "4957", "4958"])
    @test read_jc_clusters(tmp_file1) == true_dict

    tmp_file1, tmp_fo1 = mktemp()
    print(tmp_fo1, "4956\t14932\n4957\t14932\n4958\t14932\n4956\t14931\n4957\t14931\n4958\t14931\n4959\t14931\n")
    close(tmp_fo1)
    true_dict = OrderedDict("14931" => ["4956", "4957", "4958", "4959"], "14932" => ["4956", "4957", "4958"])
    @test read_jc_clusters(tmp_file1) == true_dict

end


text:
using LinearAlgebra

A = [0.5 0.7 0;
       1 1 0;
       0 0 0]

B = [1 0 0;
     0 1 1;
     -1.5 -0.5/4 1]

mat = inv(B)*A

tt = eigen(mat)

abs.(tt.values)

uu = [1;0;0]

dynamics = zeros(11,3)

dynamics[1,:] = uu

for t = 1:10
    dynamics[t+1,:] = mat*dynamics[t,:]

end

using Plots

plot(dynamics[:,1], label = "π")
plot(dynamics[:,2],label = "y")
plot(dynamics[:,3],label = "i")

ext: jl

max_stars_repo_name: danmrc/Julia-Para-Economistas

max_issues_repo_name: danmrc/Julia-Para-Economistas

max_forks_repo_name: danmrc/Julia-Para-Economistas


text:
using Base
using Test
using LinearAlgebra
using Surrogates

#1D
lb = 0.0
ub = 4.0
x = [1.0,2.0,3.0]
y = [4.0,5.0,6.0]
linear = x -> norm(x)
cubic = x -> x^3
λ = 2.3
multiquadr = x -> sqrt(x^2+λ^2)
q = 1
my_rad = RadialBasis(x,y,lb,ub,linear,q)
est = my_rad(3.0)
@test est ≈ 7.875
add_point!(my_rad,4.0,10.0)
est = my_rad(3.0)
@test est ≈ 6.499999999999991
add_point!(my_rad,[3.2,3.3,3.4],[8.0,9.0,10.0])
est_rad2 = my_rad(3.0)
@test est_rad2 ≈ 7.79959017

my_rad = RadialBasis(x,y,lb,ub,cubic,q)
q = 2
my_rad = RadialBasis(x,y,lb,ub,multiquadr,q)

#ND
x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
y = [4.0,5.0,6.0]
bounds = [[0.0,3.0,6.0],[4.0,7.0,10.0]]
my_rad = RadialBasis(x,y,bounds,z->norm(z),1)
est = my_rad((1.0,2.0,3.0))
@test est ≈ 4.0

#WITH ADD_POINT, adding singleton
x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
y = [4.0,5.0,6.0]
bounds = [[0.0,3.0,6.0],[4.0,7.0,10.0]]
my_rad = RadialBasis(x,y,bounds,z->norm(z),1)
add_point!(my_rad,(9.0,10.0,11.0),10.0)
est = my_rad((1.0,2.0,3.0))
@test est ≈ 4.0

#WITH ADD_POINT, adding more
x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
y = [4.0,5.0,6.0]
bounds = [[0.0,3.0,6.0],[4.0,7.0,10.0]]
my_rad = RadialBasis(x,y,bounds,z->norm(z),1)
add_point!(my_rad,[(9.0,10.0,11.0),(12.0,13.0,14.0)],[10.0,11.0])
est = my_rad((1.0,2.0,3.0))
@test est ≈ 4.0

lb = [0.0,0.0,0.0]
ub = [10.0,10.0,10.0]
bounds = [lb,ub]
x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
y = [4.0,5.0,6.0]
linear = x -> norm(x)
cubic = x -> norm(x)^3
λ = 2.3
multiquadr = x -> sqrt(norm(x)^2+λ^2)
q = 1
my_rad_ND = RadialBasis(x,y,bounds,linear,q)
add_point!(my_rad_ND,(3.5,4.5,1.2),18.9)
add_point!(my_rad_ND,[(3.2,1.2,6.7),(3.4,9.5,7.4)],[25.72,239.0])
my_rad_ND = RadialBasis(x,y,bounds,cubic,q)
q = 2
my_rad_ND = RadialBasis(x,y,bounds,multiquadr,q)
prediction = my_rad_ND((1.0,1.0,1.0))


text:
#! /usr/bin/env julia

using Documenter

using GeometryBasics

# Copy the README to serve as the homepage
cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"))

makedocs(
    format = Documenter.HTML(),
    sitename = "GeometryBasics.jl",
    pages = [
        "index.md",
        "primitives.md",
        "rectangles.md",
        "polygons.md",
        "meshes.md",
        "decomposition.md",
        "distancefields.md",
        "metadata.md",
        "api.md",
    ],
    modules = [GeometryBasics]
)

deploydocs(
    repo = "github.com/JuliaGeometry/GeometryBasics.jl.git",
    push_preview = true
)


text:
export get_quantified_variables_form,
    get_unbound_variables

"""
    get_quantified_variables_form(formula)

Return `formula` with all unbound variables bound to existential quantifiers.

The returned formula is equisatisfiable with `formula`.
"""
function get_quantified_variables_form(formula)
    variables = get_unbound_variables(formula)
    quantified_variables_form = formula
    for variable in variables
        quantified_variables_form = EFormula(variable, quantified_variables_form)
    end
    quantified_variables_form
end

function get_unbound_variables(formula, boundvars=Set{Variable}())
    @match formula begin
    f::Conjunction || f::Disjunction => begin
        unboundvars = get_unbound_variables.([f.formula1, f.formula2], boundvars)
        union(unboundvars...)
    end
    f::Negation => get_unbound_variables(f.formula, boundvars)
    f::AFormula || f::EFormula => begin
        push!(boundvars, f.variable)
        unboundvars = get_unbound_variables(f.formula, boundvars)
        pop!(boundvars, f.variable)
        unboundvars
    end
    f::AtomicFormula || f::Function => union(get_unbound_variables.(f.arguments, boundvars)...)
    f::Variable => f in boundvars ? Set() : Set([f])
    end
end


text:
module Move_FourBar

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
massProperties = "Aluminium"
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")
cmat    = Modia3D.defaultContactMaterial()

@assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=Ly) begin
   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vmat1; contactMaterial = cmat)) #
   frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
   frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])
   #cyl    = Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2)
   #cyl1   = Modia3D.Object3D(frame1, cyl)
   #cyl2   = Modia3D.Object3D(frame2, cyl)
end

@assembly Fourbar(;Lx = 0.1, Ly=Lx/5, Lz=Ly, groundWidth=Lx, groundHeight=0.1*Lx) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
   frame1 = Modia3D.Object3D(world ; r=[0.5*Lx, 0.0, groundWidth/2])
   # frame2 = Modia3D.Object3D(frame1; r=[Lx    , 0.0, 0.0])
   ground = Modia3D.Object3D(world , Modia3D.Box(3*Lx, groundWidth, groundHeight; material=groundMaterial);
                                              r = [1.5*Lx, -groundHeight/2, groundWidth/2],
                                              R = ModiaMath.rot_nxy([-1,0,0], [0,0,1]))
   bar1 = Bar(Lx=Lx, Ly=Ly, Lz=Lz)
   bar2 = Bar(Lx=Lx, Ly=Ly, Lz=Lz)
   bar3 = Bar(Lx=Lx, Ly=Ly, Lz=Lz)
   bar4 = Bar(Lx=2*Lx, Ly=Ly, Lz=Lz)
   bar5 = Bar(Lx=Lx, Ly=Ly, Lz=Lz)
   bar6 = Bar(Lx=Lx, Ly=Ly, Lz=Lz)

   Modia3D.connect(frame1, bar4.frame0)

   rev1 = Modia3D.Revolute(bar4.frame1, bar1.frame1; phi_start =  pi/2)
   rev2 = Modia3D.Revolute(bar1.frame2, bar2.frame1; phi_start = -pi/2)
   rev3 = Modia3D.Revolute(bar2.frame2, bar3.frame2; phi_start =  pi/2)
   rev4 = Modia3D.Revolute(bar3.frame1, bar4.frame0; phi_start =  NaN)   # phi_start = -pi/2

   rev5 = Modia3D.Revolute(bar3.frame2, bar5.frame1; phi_start = -pi/2)
   rev6 = Modia3D.Revolute(bar5.frame2, bar6.frame2; phi_start =  pi/2)
   rev7 = Modia3D.Revolute(bar6.frame1, bar4.frame2; phi_start =  NaN)
end

#@signal LinearMovement(;yStart = 0.0, yDelta = 1.0) begin
#   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
#end
#function Modia3D.computeSignal(signal::LinearMovement, sim::ModiaMath.SimulationState)
#    signal.y.value = signal.yStart + signal.yDelta*(sim.time-sim.startTime)/(sim.stopTime-sim.startTime)
#end

@signal Sine(;yStart = 0.0, A = 1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Sine, sim::ModiaMath.SimulationState)
    signal.y.value = signal.yStart + signal.A*sin(2*pi*signal.freqHz*sim.time)
end

Lx = 1.0
@assembly Move begin
   fourbar = Fourbar(Lx=Lx)

   sine = Sine(yStart=fourbar.rev1.phi.start, A=-pi/3, freqHz=0.5)
   sig  = Modia3D.SignalToFlangeAngle(sine.y)
   Modia3D.connect(sig, fourbar.rev1)
end

model  = Modia3D.SimulationModel( Move(), analysis=ModiaMath.KinematicAnalysis )
result = ModiaMath.simulate!(model, stopTime=3.0, log=false)

println("... success of Move_FourBar_Collision2.jl!")
end


text:
using StaticArrays
using Base.Threads: @threads
using StatsBase: sample
using BenchmarkTools
using EvoTrees
using EvoTrees: find_split_static!, pred_leaf, update_hist!, find_split_narrow!

x1 = BitSet([1,2,3,4,5,8,9])
x2 = BitSet([2,5,8])
x3 = BitSet([1,3])
x4 = BitSet([4,9])
bags = [[x2, x3, x4]]
intersect(x1, x2)

x1 = [1,2,3,4,5,8,9]
x2 = [2,5,8]
x3 = [1,3]
x4 = [4,9]
bags = [[x2, x3, x4]]
intersect(x1, x2)

nrows = 100_000
ncols = 100
nbins = 32

id1_int = sample(1:nrows, Int(nrows/2), replace=false, ordered=true)
id2_int = sample(1:nrows, Int(nrows/4), replace=false, ordered=false)

id1_bit = BitSet(sample(1:nrows, Int(nrows/2), replace=false, ordered=true));
id2_bit = BitSet(sample(1:nrows, Int(nrows/4), replace=false, ordered=false));

hist = zeros(32)
x_bin = sample(UInt8.(1:nbins), nrows, replace=true, ordered=false)
value = rand(nrows)

function hist_sum1(x::AbstractVector{S}, hist::AbstractVector{T}, set::I, value::AbstractVector{T}) where {S,T,I}
    hist .*= 0
    for i in set
        hist[x[i]] += value[i]
    end
    return
end

function hist_sum2(x, hist, set, value)
    hist .= 0
    for i in set
        hist[x[i]] += value[i]
    end
    return
end
@btime hist_sum1($x_bin, $hist, $id1_int, $value)
@btime hist_sum1($x_bin, $hist, $id1_bit, $value)

@btime hist_sum2($x_bin, $hist, $id1_int, $value)
@btime hist_sum2($x_bin, $hist, $id1_bit, $value)

function zeroise1(hist)
    hist .*= 0.0
end
function zeroise2(hist)
    hist .= 0.0
end
function zeroise3(hist)
    hist .= zero(eltype(hist))
end

hist_v = zeros(32)
hist_s = zeros(SVector{1, Float64}, 32)
@btime zeroise1($hist_v)
@btime zeroise2($hist_v)
@btime zeroise3($hist_v)
@btime zeroise1($hist_s)
@btime zeroise2($hist_s)

function zeroise3(hist)
    hist .= [zero(eltype(hist))]
end

hist_s = zeros(SVector{1, Float64}, 32)
@btime zeroise1($hist_s)
@btime zeroise3($hist_s)

function inter(x1, x2)
    res = intersect(x1, x2)
    return res
end
@btime inter($id1_int, $id2_int)
@btime inter($id1_bit, $id2_bit)

@time inter(id1_int, id2_int)
@time inter(id1_bit, id2_bit)
@time inter(id1_int, id2_bit)

function pushtest(src, cond, child)
    for i in src
        if in(i, cond)
            push!(child, i)
        end
    end
end

function conv(x::Vector{Int})
    BitSet(x)
end

function conv(x::BitSet)
    Int.(x)
end

src = sample(1:nrows, nrows, replace=false, ordered=false)
cond = sample(1:nrows, 50_000, replace=false, ordered=false)
cond = BitSet(cond);
child = Vector{Int}()
@time pushtest(src, cond, child)

@time src_bit = conv(src);
@time src_int = conv(src_bit)

x1 = rand(100_000, 100)

function sum(x::Matrix{T}, hist::Matrix{T})

    res = intersect(x1, x2)
    return res
end

intersect.(Ref(x1), bags[1])

@inline function update_bags_intersect(bags, set)
    # new_bags = Vector{Vector{BitSet}}(undef, 10)
    new_bags = [Vector{BitSet}() for i in 1:length(bags)]
    @threads for feat in 1:length(bags)
        new_bags[feat] = intersect.(bags[feat], Ref(set))
        # new_bags[feat] = intersect.(Ref(set), bags[feat])
    end
    return new_bags
end

nbins = 64
bin_size = 5_000
𝑖 = 1:nbins*bin_size |> collect
set = BitSet(𝑖);
bags = [Vector{BitSet}() for i in 1:nbins]
for i in 1:length(bags)
    x = sample(𝑖, length(𝑖), replace=false)
    bags[i] = [BitSet(x[bin_size*(i-1)+1:bin_size*i]) for i in 1:nbins]
end

length(bags[1][1])

set = BitSet(sample(𝑖, 10_000, replace=false));
@time new_bags = update_bags_intersect(bags, set);
length(new_bags[2][2])


text:
export Classifier
export LogReg

import EmpiricalRisks
using Regression

# ==========================================================================
# Namespacing for convenient enumerations
module Classifier

export EmpiricalRiskLearner, LogisticRegression

using ..SupervisedLearning.SupervisedModel

using EmpiricalRisks
using Regression

type EmpiricalRiskLearner{L<:Loss}
  lossfunction::L
  regularizer::Regularizer
  _coefficients::Vector{Float64}
  _cost::Float64
  _native_model::Any
end

function EmpiricalRiskLearner{L<:Loss}(lossfunction::L, regularizer::Regularizer; box=true)
  model = EmpiricalRiskLearner{L}(lossfunction, regularizer, [.0], typemax(Float64), nothing)
  box ? SupervisedModel(model) : model
end

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

type LogisticRegression
  l1_coef::Float64
  l2_coef::Float64
  _risk_learner::Any
end

function LogisticRegression(; l1_coef = .0, l2_coef = .0)
  SupervisedModel(LogisticRegression(l1_coef, l2_coef, nothing))
end

end

typealias LogReg Classifier.LogisticRegression

# ==========================================================================
# Implement the interface

using SupervisedLearning.Classifier

cost(model::EmpiricalRiskLearner) = model._cost

function cost(model::EmpiricalRiskLearner, data::EncodedInMemoryLabeledDataSource)
  X = features(data)
  t = targets(data)
  rmodel = model._native_model.rmodel
  f = if typeof(model.regularizer) == EmpiricalRisks.ZeroReg
    Regression.RiskFun(rmodel, X, t)
  else
    Regression.RegRiskFun(rmodel, model.regularizer, X, t)
  end
  mean(value(f, model._coefficients))
end

function predict(model::EmpiricalRiskLearner{LogisticLoss}, data::EncodedInMemoryLabeledDataSource)
  pm = model._native_model.rmodel.predmodel
  pred = sign(EmpiricalRisks.predict(pm, model._coefficients, features(data)))
  labeldecode(data.encoding, pred)
end

function predictProb(model::EmpiricalRiskLearner{LogisticLoss}, data::EncodedInMemoryLabeledDataSource)
  pm = model._native_model.rmodel.predmodel
  pred = tanh(EmpiricalRisks.predict(pm, model._coefficients, features(data)))
  pred / 2 + .5
end

function accuracy(model::EmpiricalRiskLearner, data::DataSource)
  gt = groundtruth(data)
  correct = sum(predict(model, data) .== gt)
  correct / length(gt)
end

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

cost(model::LogReg) = cost(model._risk_learner)
cost(model::LogReg, data::DataSource) = cost(model._risk_learner, data)
predict(model::LogReg, data::DataSource) = predict(model._risk_learner, data)
predictProb(model::LogReg, data::DataSource) = predictProb(model._risk_learner, data)
accuracy(model::LogReg, data::DataSource) = accuracy(model._risk_learner, data)

# ==========================================================================
# Implement the train! functions

function train!{E<:SignedClassEncoding, L<:UnivariateLoss, G<:Any}(
    callback::Function,
    model::EmpiricalRiskLearner{L},
    data::EncodedInMemoryLabeledDataSource{E,G,1},
    solver::Solver.GradientDescent,
    args...;
    max_iter::Int = 1000,
    nargs...)
  @assert max_iter > 0
  X = features(data)
  t = targets(data)
  b = bias(data)
  n = nvar(data)
  rmodel = if b == 0.0
    riskmodel(LinearPred(n), model.lossfunction)
  else
    riskmodel(AffinePred(n, b), model.lossfunction)
  end
  f = if typeof(model.regularizer) == EmpiricalRisks.ZeroReg
    Regression.RiskFun(rmodel, X, t)
  else
    Regression.RegRiskFun(rmodel, model.regularizer, X, t)
  end

  # inital theta
  theta = typeof(rmodel.predmodel) == LinearPred ? zeros(n) : zeros(n+1)

  # Reset training _history
  model._native_model = f
  model._coefficients = theta
  model._cost = cost(model, data)

  # define callback
  function cb(iter, theta, cost, grad)
    model._coefficients = theta
    model._cost = cost
    callback()
  end

  # Fit model
  r = Regression.solve!(Regression.GD(), f, theta,
                        Regression.Options(maxiter = max_iter, nargs...),
                        cb)

  model._coefficients = r.sol
  model._cost = r.fval
  r.niters
end

function train!{E<:ClassEncoding, G<:Any}(
    callback::Function,
    model::LogisticRegression,
    data::EncodedInMemoryLabeledDataSource{E,G,1},
    solver::Solver.GradientDescent,
    args...;
    nargs...)
  if model.l1_coef == model.l2_coef == 0.
    # No reg
    newModel = EmpiricalRiskLearner(LogisticLoss(), EmpiricalRisks.ZeroReg(), box=false)
    model._risk_learner = newModel
    train!(callback, newModel, data, solver, args...; nargs...)
  elseif model.l1_coef == 0. && model.l2_coef != 0.
    # Only l2 reg
    newModel = EmpiricalRiskLearner(LogisticLoss(), EmpiricalRisks.SqrL2Reg(model.l2_coef), box=false)
    model._risk_learner = newModel
    train!(callback, newModel, data, solver, args...; nargs...)
  elseif model.l1_coef != 0. && model.l2_coef == 0.
    # Only l1 reg
    newModel = EmpiricalRiskLearner(LogisticLoss(), EmpiricalRisks.L1Reg(model.l1_coef), box=false)
    model._risk_learner = newModel
    train!(callback, newModel, data, solver, args...; nargs...)
  else
    throw(ArgumentError("Hyperparameter of LogisticRegression have an illegal combination"))
  end
end


text:
function findpackages(pattern = "*", dir = DEVDIR)
    collect(filter(istestablepackage, glob(pattern, dir)))
end

function istestablepackage(path)
    isdir(path) || return false
    contents = readdir(path)
    "Project.toml" ∈ contents || return false
    "test" ∈ contents || return false
    "runtests.jl" ∈ readdir(joinpath(path, "test"))
end

function sourcepaths(packagepath; include_tests = true)
    srcfolder = joinpath(packagepath, "src")
    @assert isdir(srcfolder)
    srcpaths = getfoldersandjlfiles(srcfolder)
    if include_tests
        testfolder = joinpath(packagepath, "test")
        srcpaths = vcat(srcpaths, getfoldersandjlfiles(testfolder))
    end
    return srcpaths
end

function getfoldersandjlfiles(folder)
    files1 = glob("*.jl", folder)
    files2 = glob("**/*.jl", folder)
    folders = glob("**/", folder)
    return vcat(files1, files2, folders)
end


text:
using Test
using DataStructures

using preprocessing
using JSON
using util
using Method
using aggregation_strategy

function aggregateMatrix(stgTypeStr, matrix)
    stgType = parseAggStrategy(stgTypeStr)
    nrow, ncol = size(matrix)
    agg = Aggregation(stgType, nrow, ncol)

    for cand_idx in 1:ncol
        for query_idx in 1:nrow
            update_score_matrix!(agg, query_idx, cand_idx, matrix[query_idx, cand_idx])
        end
    end

    return aggregate(agg)
end

matrix = [
   0.6 0.5 0.9 -1.0 -1.0;
   0.0 -2.0 1.0 -1.0 -1.0;
   0.0 -1.0 0.0 -1.0 -1.0;
]

@test aggregateMatrix("max", matrix)==1.0
@test aggregateMatrix("avg_query", matrix)== 0.6333333333333333
@test abs(aggregateMatrix("avg_cand", matrix) - 0.02000000000) < 0.000000001
@test abs(aggregateMatrix("avg_short", matrix) - 0.63333333333333) < 0.000000001
@test abs(aggregateMatrix("avg_long", matrix) - 0.02000000000)  < 0.000000001
@test abs(aggregateMatrix("avg_query_cand", matrix) -  0.3266666666)  < 0.000000001

@test abs(aggregateMatrix("avg_query", transpose(matrix)) - 0.02000000000) < 0.000000001
@test abs(aggregateMatrix("avg_cand", transpose(matrix)) - 0.63333333333333)  < 0.000000001
@test abs(aggregateMatrix("avg_short", transpose(matrix)) - 0.63333333333333) < 0.000000001
@test abs(aggregateMatrix("avg_long", transpose(matrix)) - 0.02000000000)  < 0.000000001
@test abs(aggregateMatrix("avg_query_cand", transpose(matrix)) -   0.3266666666)  < 0.000000001

idfs = [(0x00010f8f, 0.0007444291882413274), (0x00010f90, 0.0007444291882413274), (0x00020d3a, 0.0006782577048420984), (0x00020d3a, 0.0006782577048420984), (0x00000446, 0.12638753329252758), (0x00010f91, 0.001505401247332462), (0x00010f8f, 0.0007444291882413274), (0x00010f90, 0.0007444291882413274), (0x00020d3a, 0.0006782577048420984), (0x00020d3a, 0.0006782577048420984), (0x00000446, 0.12638753329252758), (0x00010f91, 0.001505401247332462), (0x0000c9a1, 0.0010918294760872802), (0x00012afe, 0.0009264007675892074), (0x0002b716, 0.0001323429667984582), (0x0002b717, 0.0001323429667984582), (0x0002b718, 0.0001323429667984582)]
df = zeros(500000)

for (idx, df2) in idfs
    df[idx] = df2 * 100
end

query = [[0x00010f8f, 0x00010f90, 0x00020d3a, 0x00020d3a, 0x00000446, 0x00010f91, 0x00022c31, 0x00022c32, 0x0002e06c, 0x0002e06d, 0x0002e06e, 0x0002e06f, 0x0002277c, 0x000009c1, 0x00000780, 0x0000a336, 0x0000703f, 0x00006da4, 0x00011e66, 0x00006da4, 0x00006da5, 0x0000b3a0, 0x00006da6, 0x0002675a], [0x00010f8f, 0x00010f90, 0x00020d3a, 0x00020d3a, 0x00000446, 0x00010f91, 0x0000c9a1, 0x00012afe, 0x0002b716, 0x0002b717, 0x0002b718], [0x00010f8f, 0x00010f90, 0x00020d3a, 0x00020d3a, 0x00000446, 0x00010f91, 0x00005f39, 0x00007072, 0x00014fe2, 0x00014fe3, 0x00014fe4, 0x00001818, 0x000003eb, 0x0000b296, 0x00007d53, 0x0000adc2, 0x00000a60, 0x000088ad, 0x000088ae, 0x0002436b, 0x0001012b, 0x0000794a, 0x0001e4af, 0x0000bb81]]

cand = [[0x00010f8f, 0x00010f90, 0x00020d3a, 0x00020d3a, 0x00000446, 0x00010f91],
[0x00010f8f, 0x00010f90, 0x00020d3a, 0x00020d3a, 0x00000446, 0x00010f91, 0x0000c9a1, 0x00012afe, 0x0002b716, 0x0002b717, 0x0002b718]]

tracesim = TraceSim(1.5 ,0.8, 1.0, false, 0,false, false, false)

matrix = zeros((length(query), length(cand)))

for cand_idx in 1:length(cand)
    for query_idx in 1:length(query)
        matrix[query_idx, cand_idx] = similarity(tracesim, query[query_idx], cand[cand_idx], df)
    end
end

# 237531, 0.8572043698222129
print(matrix)
print(aggregateMatrix("avg_query_cand", matrix))

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
# This script provides implementations of the
# Adams–Bashforth's methods of different orders
# for solving a system of differential equations
# and compares them to the exact solutions and
# other methods

println('\n', " "^4, "> Loading the packages...")

using LinearAlgebra # Norm
using Plots # Plotting

# Use the GR backend for plots
gr()

# Change the default font for plots
default(fontfamily="Computer Modern", dpi=300, legend=:topright)

# Define the value of ϰ
const ϰ = -1

# Integrate the three-body problem using the
# Euler's method, return the values of position
# and velocity on each step
function motion3_euler(
    r₀₁::Vector{F},
    r₀₂::Vector{F},
    r₀₃::Vector{F},
    v₀₁::Vector{F},
    v₀₂::Vector{F},
    v₀₃::Vector{F},
    h::F,
    n::I,
)::Tuple{Matrix{F},Matrix{F},Matrix{F},Matrix{F},Matrix{F},Matrix{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Prepare output matrixes
    r₁ = Matrix{F}(undef, n + 1, N)
    r₂ = Matrix{F}(undef, n + 1, N)
    r₃ = Matrix{F}(undef, n + 1, N)
    v₁ = Matrix{F}(undef, n + 1, N)
    v₂ = Matrix{F}(undef, n + 1, N)
    v₃ = Matrix{F}(undef, n + 1, N)
    r₁[1, :] = copy(r₀₁)
    r₂[1, :] = copy(r₀₂)
    r₃[1, :] = copy(r₀₃)
    v₁[1, :] = copy(v₀₁)
    v₂[1, :] = copy(v₀₂)
    v₃[1, :] = copy(v₀₃)
    # Compute the solutions
    for i in 2:(n + 1)
        ρ₁₂ = norm(r₁[i - 1, :] - r₂[i - 1, :])^3
        ρ₁₃ = norm(r₁[i - 1, :] - r₃[i - 1, :])^3
        ρ₂₃ = norm(r₂[i - 1, :] - r₃[i - 1, :])^3
        for k in 1:N
            r₁[i, k] = r₁[i - 1, k] + h * v₁[i - 1, k]
            v₁[i, k] = v₁[i - 1, k] + h * ϰ * (r₁[i - 1, k] - r₂[i - 1, k]) / ρ₁₂ +
                                      h * ϰ * (r₁[i - 1, k] - r₃[i - 1, k]) / ρ₁₃
            r₂[i, k] = r₂[i - 1, k] + h * v₂[i - 1, k]
            v₂[i, k] = v₂[i - 1, k] + h * ϰ * (r₂[i - 1, k] - r₁[i - 1, k]) / ρ₁₂ +
                                      h * ϰ * (r₂[i - 1, k] - r₃[i - 1, k]) / ρ₂₃
            r₃[i, k] = r₃[i - 1, k] + h * v₃[i - 1, k]
            v₃[i, k] = v₃[i - 1, k] + h * ϰ * (r₃[i - 1, k] - r₁[i - 1, k]) / ρ₁₃ +
                                      h * ϰ * (r₃[i - 1, k] - r₂[i - 1, k]) / ρ₂₃
        end
    end
    return r₁, r₂, r₃, v₁, v₂, v₃
end

# Integrate the three-body problem using the
# two-step Adams–Bashforth's method, return the
# values of position and velocity on each step
function motion3_ab2(
    r₀₁::Vector{F},
    r₀₂::Vector{F},
    r₀₃::Vector{F},
    v₀₁::Vector{F},
    v₀₂::Vector{F},
    v₀₃::Vector{F},
    h::F,
    n::I,
)::Tuple{Matrix{F},Matrix{F},Matrix{F},Matrix{F},Matrix{F},Matrix{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Prepare output matrixes
    r₁ = Matrix{F}(undef, n + 1, N)
    r₂ = Matrix{F}(undef, n + 1, N)
    r₃ = Matrix{F}(undef, n + 1, N)
    v₁ = Matrix{F}(undef, n + 1, N)
    v₂ = Matrix{F}(undef, n + 1, N)
    v₃ = Matrix{F}(undef, n + 1, N)
    r₁[1, :] = copy(r₀₁)
    r₂[1, :] = copy(r₀₂)
    r₃[1, :] = copy(r₀₃)
    v₁[1, :] = copy(v₀₁)
    v₂[1, :] = copy(v₀₂)
    v₃[1, :] = copy(v₀₃)
    # Compute the second value of the solution
    # by using the one-step Euler's method
    ρ₁₂ = norm(r₁[1, :] - r₂[1, :])^3
    ρ₁₃ = norm(r₁[1, :] - r₃[1, :])^3
    ρ₂₃ = norm(r₂[1, :] - r₃[1, :])^3
    for k in 1:N
        r₁[2, k] = r₁[1, k] + h * v₁[1, k]
        v₁[2, k] = v₁[1, k] + h * ϰ * (r₁[1, k] - r₂[1,k]) / ρ₁₂ +
                              h * ϰ * (r₁[1, k] - r₃[1,k]) / ρ₁₃
        r₂[2, k] = r₂[1, k] + h * v₂[1, k]
        v₂[2, k] = v₂[1, k] + h * ϰ * (r₂[1, k] - r₁[1,k]) / ρ₁₂ +
                              h * ϰ * (r₂[1, k] - r₃[1,k]) / ρ₂₃
        r₃[2, k] = r₃[1, k] + h * v₃[1, k]
        v₃[2, k] = v₃[1, k] + h * ϰ * (r₃[1, k] - r₁[1,k]) / ρ₁₃ +
                              h * ϰ * (r₃[1, k] - r₂[1,k]) / ρ₂₃
    end
    # Define a couple of independent coefficients
    k₁ = 3 / 2 * h
    k₂ = 1 / 2 * h
    # Compute the rest in two steps
    for i in 3:(n + 1)
        ρ₁₂_ₖ = norm(r₁[i - 1, :] - r₂[i - 1, :])^3
        ρ₁₃_ₖ = norm(r₁[i - 1, :] - r₃[i - 1, :])^3
        ρ₂₃_ₖ = norm(r₂[i - 1, :] - r₃[i - 1, :])^3
        ρ₁₂_ₖ₋₁ = norm(r₁[i - 2, :] - r₂[i - 2, :])^3
        ρ₁₃_ₖ₋₁ = norm(r₁[i - 2, :] - r₃[i - 2, :])^3
        ρ₂₃_ₖ₋₁ = norm(r₂[i - 2, :] - r₃[i - 2, :])^3
        # Define a couple of dependent coefficients
        for k in 1:N
            r₁[i, k] = r₁[i - 1, k] + k₁ * v₁[i - 1, k] - k₂ * v₁[i - 2, k]
            v₁[i, k] = v₁[i - 1, k] + k₁ * ϰ * (r₁[i - 1, k] - r₂[i - 1, k]) / ρ₁₂_ₖ +
                                      k₁ * ϰ * (r₁[i - 1, k] - r₃[i - 1, k]) / ρ₁₃_ₖ -
                                      k₂ * ϰ * (r₁[i - 2, k] - r₂[i - 2, k]) / ρ₁₂_ₖ₋₁ -
                                      k₂ * ϰ * (r₁[i - 2, k] - r₃[i - 2, k]) / ρ₁₃_ₖ₋₁
            r₂[i, k] = r₂[i - 1, k] + k₁ * v₂[i - 1, k] - k₂ * v₂[i - 2, k]
            v₂[i, k] = v₂[i - 1, k] + k₁ * ϰ * (r₂[i - 1, k] - r₁[i - 1, k]) / ρ₁₂_ₖ +
                                      k₁ * ϰ * (r₂[i - 1, k] - r₃[i - 1, k]) / ρ₂₃_ₖ -
                                      k₂ * ϰ * (r₂[i - 2, k] - r₁[i - 2, k]) / ρ₁₂_ₖ₋₁ -
                                      k₂ * ϰ * (r₂[i - 2, k] - r₃[i - 2, k]) / ρ₂₃_ₖ₋₁
            r₃[i, k] = r₃[i - 1, k] + k₁ * v₃[i - 1, k] - k₂ * v₃[i - 2, k]
            v₃[i, k] = v₃[i - 1, k] + k₁ * ϰ * (r₃[i - 1, k] - r₁[i - 1, k]) / ρ₁₃_ₖ +
                                      k₁ * ϰ * (r₃[i - 1, k] - r₂[i - 1, k]) / ρ₂₃_ₖ -
                                      k₂ * ϰ * (r₃[i - 2, k] - r₁[i - 2, k]) / ρ₁₃_ₖ₋₁ -
                                      k₂ * ϰ * (r₃[i - 2, k] - r₂[i - 2, k]) / ρ₂₃_ₖ₋₁
        end
    end
    return r₁, r₂, r₃, v₁, v₂, v₃
end

println('\n', " "^4, "> Integrating the three-body problem...")

# Define the initial values of the position and velocity
r₀₁ = [0.97, -0.2431]
r₀₂ = [-0.97, 0.2431]
r₀₃ = [0.0, 0.0]
v₀₁ = [0.4662, 0.4324]
v₀₂ = [0.4662, 0.4324]
v₀₃ = [-0.9324, -0.8647]

@userplot CirclePlot
@recipe function f(cp::CirclePlot)
    x, y, i = cp.args
    n = length(x)
    inds = circshift(1:n, 1 - i)
    linewidth --> range(0, 10, length=n)
    seriesalpha --> range(0, 1, length=n)
    aspect_ratio --> 1
    label --> false
    x[inds], y[inds]
end

h = 0.0001
n = UInt(62887)

r₁ₑ, r₂ₑ, r₃ₑ, _, _, _ = motion3_ab2(r₀₁, r₀₂, r₀₃, v₀₁, v₀₂, v₀₃, h, n)

# Animate the figure-8 orbit
anim = @animate for i ∈ 1:446:n
    circleplot(r₁ₑ[:, 1], r₁ₑ[:, 2], i; xlabel="x", ylabel="y")
    circleplot!(r₂ₑ[:, 1], r₂ₑ[:, 2], i)
    circleplot!(r₃ₑ[:, 1], r₃ₑ[:, 2], i)
end

# Save the animation
gif(anim, "$(@__DIR__)/../plots/anim_3body.gif"; fps=30, show_msg=false)

println(" "^6, "* The animation `anim_3body.gif` is saved. *")

println()


text:
function run_test(M,net_param,time_param,dt)

N = length(M[1,:]);
extinput = zeros(N);
u = zeros(N);

x, r, b, stim = net_param[1:N], net_param[N+1:2*N], net_param[2*N+1], net_param[2*N+2:3*N+1];
stim_on, stim_off, test_time = time_param[1], time_param[2], time_param[3];

#===== Run simulation =====#
dtheta(x_var,I_var) = 1. - cos(x_var) + I_var.*(1. + cos(x_var));
dr(r_var) = -b*r_var;

t = 0.;

utest = zeros(Int(test_time/dt),N);

for i=1:Int(test_time/dt)

      if i > Int(stim_on/dt) && i <= Int(stim_off/dt)
            extinput[:] = stim[:];
      else
            extinput[:] = zeros(N)[:];
      end

      # Define synaptic activity
      u = M*r;

      # Update neuron phase
      k1 = dt*dtheta(x,      u+extinput);
      k2 = dt*dtheta(x+k1/2, u+extinput);
      k3 = dt*dtheta(x+k2/2, u+extinput);
      k4 = dt*dtheta(x+k3,   u+extinput);
      xnext  = x + (k1+2*k2+2*k3+k4)/6;

      # Update filtered spikes
      l1 = dt*dr(r);
      l2 = dt*dr(r+l1/2);
      l3 = dt*dr(r+l2/2);
      l4 = dt*dr(r+l3);

      # Spike detection.
      idx1 = xnext - x .> 0.;
      idx2 = xnext - x .> mod(pi - mod(x,2*pi),2*pi); # distance to next pi.
      idx  = idx1.*idx2;
      ind  = collect(1:N)[idx[:]]';

      # Update each neuron's spikes
      if !isempty(ind)
          r[idx[:]] = r[idx[:]] + b;
      end

      # Update x, r
      x = x + (k1+2*k2+2*k3+k4)/6;
      r = r + (l1+2*l2+2*l3+l4)/6;
      u = M*r;
      t = t + dt;

      # save test data
      utest[i,:] = u[:];
end

return utest

end


text:
function train!(net::Network, train_data::Vector{TestDatum},
        test_data::Vector{TestDatum}, epoch::Int,
        mini_batch_size::Int, eta::Float64; verbose = false)
    @assert all((length(train_data), length(test_data)) .> 1);
    @assert all(inlength(net) .== inlength.(train_data))
    @assert all(outlength(net) .== outlength.(train_data))
    @assert all(inlength(net) .== inlength.(test_data))
    @assert all(outlength(net) .== outlength.(test_data))

    #Gradient Buffers
    #For the first layer the gradients are not usable.
    dCx_dw_buffer::NTuple{net.depth, Matrix{Float64}} =
        tuple(zeros(Float64,1,1),[Matrix{Float64}(undef, net.lsizes[l], net.lsizes[l-1]) for l in 2:net.depth]...);
    dCx_db_buffer::NTuple{net.depth, Vector{Float64}} =
        tuple(zeros(Float64,1),[Vector{Float64}(undef,n) for n in net.lsizes[2:end]]...);

    #Extra Buffers
    #Store the gradients of dC/da for all layers...
    #For the first layer the value are not usable.
    dCx_da_buffer::NTuple{net.depth, Vector{Float64}} = tuple(zeros(Float64,1),[Vector{Float64}(undef,n) for n in net.lsizes[2:end]]...)
    #Will Store a partial result that is redundant.
    #For the first layer the value are not usable.
    Bj_buffer::NTuple{net.depth, Vector{Float64}} = tuple(zeros(Float64,1),[Vector{Float64}(undef,n) for n in net.lsizes[2:end]]...);

    #Epoch
    if verbose
        acost::Float64 = 0.0;
        bcost::Float64 = 0.0;
        startime::Float64 = 0.0;
    end

    for e in 1:epoch

        #Info
        if verbose
            println("###### Epoch $e ###################");
            startime = time();
            if e == 1
                bcost = cost!(net, test_data);
            else
                bcost = acost;
            end
            println("Cost Before $bcost");
        end

        #Unsorting data
        shuffle!(train_data);

        #Mini batching
        train_subset_length = min(mini_batch_size, length(train_data));
        for tdi in train_subset_length:train_subset_length:length(train_data)

            #Mini_batch
            train_subset = train_data[tdi - train_subset_length + 1:tdi];

            #Updating gradients for this mini batch
            update_gradients!(net, train_subset, dCx_dw_buffer, dCx_db_buffer, dCx_da_buffer, Bj_buffer);

            #appliying gradient desendence
            update_ws!(net, dCx_dw_buffer, eta);
            update_bs!(net, dCx_db_buffer, eta);

        end # mini batch train

        #Evaluate Info
        if verbose
            acost = cost!(net, test_data);
            # println(net);
            print("Cost After  $acost");
            println(" ($(acost - bcost))");
            println("Accuracy $(evaluate!(net,test_data))");
            println("Time $(round(Int,time() - startime)) seconds")
            println();
        end

    end #Epoch

end

function update_gradients!(net, train_data, dCx_dw_buffer, dCx_db_buffer, dCx_da_buffer, Bj_buffer)

    #Reset gradients to 0.0
    for l in 1:net.depth
        dCx_dw_buffer[l] .= 0.0;
        dCx_da_buffer[l] .= 0.0;
    end

    #Accumulate gradients for each testdata
    for train_datum in train_data
        accumulative_backpropagation!(net, train_datum, dCx_dw_buffer, dCx_db_buffer, dCx_da_buffer, Bj_buffer);
    end

    #Averaging
    average_grads!(net, dCx_dw_buffer,dCx_db_buffer, length(train_data));

end

function average_grads!(net, dCx_dw_buffer, dCx_db_buffer, train_set_length)
    for l in 2:net.depth
        prevl = l - 1;
        for j in 1:net.lsizes[l]
            for k in 1:net.lsizes[prevl]
                dCx_dw_buffer[l][j,k] = dCx_dw_buffer[l][j,k] / train_set_length;
            end
            dCx_db_buffer[l][j] = dCx_db_buffer[l][j] / train_set_length;
        end
    end
end

function update_ws!(net, dCx_dw_buffer, eta)
    for l in 2:net.depth
        prevl = l -1;
        for j in 1:net.lsizes[l]
            for k in 1:net.lsizes[prevl]
                net.w[l][j,k] -= eta * dCx_dw_buffer[l][j,k];
            end
        end
    end
end

function update_bs!(net, dCx_db_buffer, eta)
    for l in 2:net.depth
        for j in 1:net.lsizes[l]
            net.b[l][j] -= eta * dCx_db_buffer[l][j];
        end
    end
end

function evaluate!(net, test_data)
    good_answers = 0.0;
    for test_datum in test_data
        expout = output(test_datum);
        netout = feedforward!(net, input(test_datum));

        if findmax(expout)[2] == findmax(netout)[2]
           good_answers = good_answers + 1;
        end
    end

    return good_answers/length(test_data);
end


text:
# https://github.com/JuliaPlots/UnicodePlots.jl

# don't warn on unsupported... there's just too many warnings!!
warn_on_unsupported_args(::UnicodePlotsBackend, plotattributes::KW) = nothing

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

_canvas_map() = (
    ascii = UnicodePlots.AsciiCanvas,
    block = UnicodePlots.BlockCanvas,
    braille = UnicodePlots.BrailleCanvas,
    density = UnicodePlots.DensityCanvas,
    dot = UnicodePlots.DotCanvas,
    heatmap = UnicodePlots.HeatmapCanvas,
    lookup = UnicodePlots.LookupCanvas,
)

# do all the magic here... build it all at once, since we need to know about all the series at the very beginning
function unicodeplots_rebuild(plt::Plot{UnicodePlotsBackend})
    plt.o = UnicodePlots.Plot[]
    canvas_map = _canvas_map()
    for sp in plt.subplots
        xaxis = sp[:xaxis]
        yaxis = sp[:yaxis]
        xlim = collect(axis_limits(sp, :x))
        ylim = collect(axis_limits(sp, :y))

        # we set x/y to have a single point, since we need to create the plot with some data.
        # since this point is at the bottom left corner of the plot, it shouldn't actually be shown
        x = Float64[xlim[1]]
        y = Float64[ylim[1]]

        # create a plot window with xlim/ylim set, but the X/Y vectors are outside the bounds
        canvas_type = if (ct = _canvas_type[]) == :auto
            isijulia() ? :ascii : :braille
        else
            ct
        end

        kw = (
            title = sp[:title],
            xlim = xlim,
            ylim = ylim,
            border = isijulia() ? :ascii : :solid,
            xlabel = xaxis[:guide],
            ylabel = yaxis[:guide],
            xscale = xaxis[:scale],
            yscale = yaxis[:scale],
        )

        o = UnicodePlots.Plot(x, y, canvas_map[canvas_type]; kw...)
        for series in series_list(sp)
            o = addUnicodeSeries!(sp, o, kw, series, sp[:legend] != :none)
        end

        for ann in sp[:annotations]
            x, y, val = locate_annotation(sp, ann...)
            o = UnicodePlots.annotate!(
                o, x, y, val.str;
                color = up_color(val.font.color), halign = val.font.halign, valign = val.font.valign
            )
        end

        push!(plt.o, o)  # save the object
    end
end

function up_color(col)
    if typeof(col) <: UnicodePlots.UserColorType
        color = col
    elseif typeof(col) <: RGBA
        col = convert(ARGB32, col)
        color = map(Int, (red(col).i, green(col).i, blue(col).i))
    else
        color = :auto
    end
    color
end

# add a single series
function addUnicodeSeries!(
    sp::Subplot{UnicodePlotsBackend},
    up::UnicodePlots.Plot,
    kw, series, addlegend::Bool,
)
    st = series[:seriestype]

    # get the series data and label
    x, y = if st == :straightline
        straightline_data(series)
    elseif st == :shape
        shape_data(series)
    else
        float(series[:x]), float(series[:y])
    end

    # special handling (src/interface)
    if st == :histogram2d
        kw[:xlim][:] .= kw[:ylim][:] .= 0
        return UnicodePlots.densityplot(x, y; kw...)
    elseif st == :heatmap
        rng = range(0, 1, length = length(UnicodePlots.COLOR_MAP_DATA[:viridis]))
        cmap = [(red(c), green(c), blue(c)) for c in get(get_colorgradient(series), rng)]
        return UnicodePlots.heatmap(
            series[:z].surf;
            zlabel = sp[:colorbar_title],
            colormap = cmap,
            kw...
        )
    elseif st == :spy
        return UnicodePlots.spy(series[:z].surf; kw...)
    end

    # now use the ! functions to add to the plot
    if st in (:path, :straightline, :shape)
        func = UnicodePlots.lineplot!
    elseif st == :scatter || series[:markershape] != :none
        func = UnicodePlots.scatterplot!
    else
        error("Series type $st not supported by UnicodePlots")
    end

    label = addlegend ? series[:label] : ""

    for (n, segment) in enumerate(series_segments(series, st; check = true))
        i, rng = segment.attr_index, segment.range
        lc = get_linecolor(series, i)
        up = func(up, x[rng], y[rng]; color = up_color(lc), name = n == 1 ? label : "")
    end

    for (xi, yi, str, fnt) in EachAnn(series[:series_annotations], x, y)
        up = UnicodePlots.annotate!(
            up, xi, yi, str;
            color = up_color(fnt.color), halign = fnt.halign, valign = fnt.valign
        )
    end

    return up
end

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

# since this is such a hack, it's only callable using `png`... should error during normal `show`
function png(plt::Plot{UnicodePlotsBackend}, fn::AbstractString)
    fn = addExtension(fn, "png")

    @static if Sys.isapple()
        # make some whitespace and show the plot
        println("\n\n\n\n\n\n")
        gui(plt)

        # BEGIN HACK

        # wait while the plot gets drawn
        sleep(0.5)

        # use osx screen capture when my terminal is maximized and cursor starts at the bottom (I know, right?)
        # TODO: compute size of plot to adjust these numbers (or maybe implement something good??)
        run(`screencapture -R50,600,700,420 $fn`)

        # END HACK (phew)
        return
    elseif Sys.islinux()
        run(`clear`)
        gui(plt)
        run(`import -window $(ENV["WINDOWID"]) $fn`)
        return
    end

    error(
        "Can only savepng on MacOS or Linux with UnicodePlots (though even then I wouldn't do it)",
    )
end

# -------------------------------
Base.show(plt::Plot{UnicodePlotsBackend}) = _show(stdout, MIME("text/plain"), plt)

function _show(io::IO, ::MIME"text/plain", plt::Plot{UnicodePlotsBackend})
    unicodeplots_rebuild(plt)
    nr, nc = size(plt.layout)
    lines_colored = Array{Union{Nothing,Vector{String}}}(undef, nr, nc)
    lines_uncolored = copy(lines_colored)
    l_max = zeros(Int, nr)
    buf = IOBuffer()
    cbuf = IOContext(buf, :color => true)
    sps = wmax = 0
    for r in 1:nr
        lmax = 0
        for c in 1:nc
            l = plt.layout[r, c]
            if l isa GridLayout && size(l) != (1, 1)
                @error "UnicodePlots: complex nested layout is currently unsupported !"
            else
                if get(l.attr, :blank, false)
                    lines_colored[r, c] = lines_uncolored[r, c] = nothing
                else
                    sp = plt.o[sps += 1]
                    show(cbuf, sp)
                    colored = String(take!(buf))
                    uncolored = replace(colored, r"\x1B\[[0-9;]*[a-zA-Z]" => "")
                    lines_colored[r, c] = lc = split(colored, "\n")
                    lines_uncolored[r, c] = lu = split(uncolored, "\n")
                    lmax = max(length(lc), lmax)
                    wmax = max(maximum(length.(lu)), wmax)
                end
            end
        end
        l_max[r] = lmax
    end
    empty = ' '^wmax
    for r in 1:nr
        for n in 1:l_max[r]
            for c in 1:nc
                pre = c == 1 ? '\0' : ' '
                lc = lines_colored[r, c]
                if lc === nothing || length(lc) < n
                    print(io, pre, empty)
                else
                    lu = lines_uncolored[r, c]
                    print(io, pre, lc[n], ' '^(wmax - length(lu[n])))
                end
            end
            println(io)
        end
        r < nr && println(io)
    end
    nothing
end

function _display(plt::Plot{UnicodePlotsBackend})
    unicodeplots_rebuild(plt)
    map(display, plt.o)
    nothing
end


text:
# Interactions of the Julia and Python garbage collection:
#
#    * When we wrap a Python object in a Julia object jo (e.g. a PyDict),
#      we keep an explicit PyObject reference inside jo, whose finalizer
#      decrefs the Python object when it is called.
#
#    * When we wrap a Julia object jo inside a Python object po
#      (e.g a numpy array), we add jo to the pycall_gc dictionary,
#      keyed by a weak reference to po.  The Python weak reference
#      allows us to register a callback function that is called
#      when po is deallocated, and this callback function removes
#      jo from pycall_gc so that Julia can garbage-collect it.

const pycall_gc = Dict{PyPtr,Any}()

function weakref_callback(callback::PyPtr, wo::PyPtr)
    delete!(pycall_gc, wo)
    ccall((@pysym :Py_DecRef), Cvoid, (PyPtr,), wo)
    return pyincref_(pynothing[])
end

const weakref_callback_obj = PyNULL() # weakref_callback Python method

# Python expects the PyMethodDef structure to be a constant, so
# we put it in a global to prevent gc.
const weakref_callback_meth = Ref{PyMethodDef}()

# "embed" a reference to jo in po, using the weak-reference mechanism
function pyembed(po::PyObject, jo::Any)
    # If there's a need to support immutable embedding,
    # the API needs to be changed to return the pointer.
    isimmutable(jo) && ArgumentError("pyembed: immutable argument not allowed")
    if ispynull(weakref_callback_obj)
        cf = @cfunction(weakref_callback, PyPtr, (PyPtr,PyPtr))
        weakref_callback_meth[] = PyMethodDef("weakref_callback", cf, METH_O)
        copy!(weakref_callback_obj,
              PyObject(@pycheckn ccall((@pysym :PyCFunction_NewEx), PyPtr,
                                       (Ref{PyMethodDef}, Ptr{Cvoid}, Ptr{Cvoid}),
                                       weakref_callback_meth, C_NULL, C_NULL)))
    end
    wo = @pycheckn ccall((@pysym :PyWeakref_NewRef), PyPtr, (PyPtr,PyPtr),
                         po, weakref_callback_obj)
    pycall_gc[wo] = jo
    return po
end


text:
#=
This file is part of Pseudospectra.jl.

Julia translation
copyright (c) 2017 Ralph Smith

Portions derived from EigTool
Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
of the University of Oxford, and the EigTool Developers. All rights reserved.

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
    recalc_levels(σ, ax)

construct a reasonable set of contour levels for displaying `log10(σ)`
where `σ` is a meshed field on axes `ax`.
"""
function recalc_levels(sigmin,ax)
    err = 0
    smin,smax = extrema(sigmin)
    if smax <= smallσ
        levels = Vector{typeof(smin)}()
        err = -2
        return levels,err
    end
    if smin <= smallσ
        smin=minimum(sigmin[sigmin .>= smallσ])
    end
    max_val = log10(smax)
    min_val = log10(smin)
    num_lines = 8
    scale = min(ax[4]-ax[3],ax[2]-ax[1])
    last_lev = log10(0.03*scale)
    if max_val >= last_lev
        num_lines = ceil(Int,num_lines*(last_lev-min_val)/(max_val-min_val))
        max_val = last_lev
    end
    if max_val < min_val
        max_val = log10(smin) + 0.1*log10(smax/smin)
        num_lines = 3
    end
    max_lines = num_lines
    if num_lines > 0
        stepsize = max(round((max_val-min_val)/num_lines*4)/4,0.25)
        (stepsize >= 0.75) && (stepsize = 1.0)
        if (max_val - min_val)/stepsize < max_lines
            stepsize = max(round((max_val-min_val)/num_lines*10)/10,0.1)
            (stepsize == 0.3) && (stepsize = 0.2)
            (stepsize == 0.4) && (stepsize = 0.5)
            ((stepsize >= 0.6) && (stepsize <= 0.8)) && (stepsize = 0.5)
            (stepsize >= 0.9) && (stepsize = 1.0)
        end
        if (max_val - min_val)/stepsize < ceil(max_lines/2)
            stepsize = (max_val - min_val) / max_lines
        else
            min_val = ceil(min_val/stepsize) * stepsize
            max_val = floor(max_val/stepsize) * stepsize
        end
        if stepsize > 0
            levels = collect(min_val:stepsize:max_val)
        else
            levels = [min_val,max_val]
        end
    end
    ll = length(levels)
    levels = reverse(levels[end:-max(floor(Int,ll/9),1):1],dims=1)
    (length(levels) == 1) && (levels = levels .* ones(2))
    # upstream has commented-out normality message logic here
    return levels,err
end

"""
determine a grid size which will result in a quick computation
"""
function setgridsize(n,nmin,nmax,iscplx)
    npts = round(Int,min(max(300/n^(5/6),nmin),nmax))
    iscplx && (npts = max(nmin,floor(Int,3*npts/4)))
    isodd(npts) && (npts += 1)
    npts
end

"""
format a time interval (in seconds) into a short string
"""
function prettytime(ttime)
    if ttime < 180
        total_time = ceil(Int,ttime)
        str = @sprintf("%d seconds",total_time)
    elseif ttime < 10800
        total_time = ceil(Int,ttime/60)
        str = @sprintf("%d minutes",total_time)
    else
        total_time = round(Int,ttime/360)/10
        str = @sprintf("%.1f hours",total_time)
    end
    return total_time,str
end

isvalidax(ax) = false
function isvalidax(ax::Vector{T}) where T<:Real
    (length(ax) == 4) && (ax[1] < ax[2]) && (ax[3] < ax[4])
end

"""
construct a mesh for the imaginary axis to exploit symmetry

shift_axes(ax,npts) -> y,n_mirror
"""
function shift_axes(ax,npts)
    # note: linspace arg len::Float is ok
    if ax[4] > -ax[3]
        tpts = floor(Int,(ax[4]/(ax[4]-ax[3]))*npts)
        y = collect(range(ax[4]/(2*(tpts-1)+1), stop=ax[4], length=tpts))
        n = count(x -> (x < -ax[3]),y)
        if abs(y[min(length(y),n+1)] + ax[3]) > 1e-15
            y = vcat(y[1:n],-ax[3],y[n+1:end])
        end
        num_mirror = n+1
    elseif ax[4] < -ax[3]
        bpts = floor(Int,(ax[3] / (ax[3] - ax[4]))*npts)
        y = collect(range(ax[3], stop=ax[3]/(2*(bpts-1)+1), length=bpts))
        n = count(x -> (x < -ax[4]),y)
        if abs(y[max(1,n)] + ax[4]) > 1e-15
            y = vcat(y[1:n],-ax[4],y[n+1:end])
        end
        num_mirror = -(length(y)-(n+1)+1)
    else
        if isodd(npts)
            y = collect(range(0, stop=ax[4], length=round(Int,(npts+1)/2)))
            num_mirror = ((npts+1) >> 1) -1
        else
            step = (ax[4] - ax[3]) / npts
            y = collect(range(step/2, stop=ax[4], length=round(Int,npts/2)))
            num_mirror = npts >> 1
        end
    end
    return y,num_mirror
end

function get_step_size(n,ly,routine)
    if n < 100
        step = max(1,floor(Int,ly/8))
    else
        step = min(ly,max(1,floor(Int,4*ly/n)))
    end
    # upstream decreases by factor of 4 if fast implementation is missing
    return step
end

"""
Prompt user for key character, with optional explanatory details.
Returns index of selected option.
"""
function replqdlg(query,details=""; options=["Yes","No"])
    abbrevs = [s[1] for s in options]
    nopts = length(options)
    qstr = (isempty(details) ? "): " : ",?): ")
    while true
        print(query," (",join(options,","),qstr)
        str = readline()
        c = uppercase(str[1])
        if !isempty(details) && (c == '?')
            println(details)
            continue
        end
        for i in 1:nopts
            if c == abbrevs[i]
                return i
            end
        end
        println("Please respond with one of $(abbrevs).")
    end
    false
end

dummyqdlg(query,details=""; options=["Yes","No"]) = 1

"""
compute reasonable axis limits for displaying a vector of points & environs
"""
function vec2ax(dispvec)
    ca = [minimum(real(dispvec)), maximum(real(dispvec)),
          minimum(imag(dispvec)), maximum(imag(dispvec))]
    if ca[1]==ca[2]
        ca[1] -= 0.5
        ca[2] += 0.5
    end
    if ca[3]==ca[4]
        ca[3] -= 0.5
        ca[4] += 0.5
    end
    ext = max(ca[2]-ca[1], ca[4]-ca[3]) / 3
    ax = [ca[1]-ext,ca[2]+ext,ca[3]-ext,ca[4]+ext]
    ax[1],ax[2] = tidyaxes(ax[1],ax[2],0.02)
    ax[3],ax[4] = tidyaxes(ax[3],ax[4],0.02)
    return ax
end

"""
Compute triangular factor(s) for a non-square matrix.

As used to prepare for pseudospectra computation.
"""
function rect_fact(A)
    m,n = size(A)
    if m >= 2*n
        # QR form
        F = qrfact(A[n+1:end,:])
        S = vcat(A[1:n,:],F[:R])
        T = I
    else
        # QZ form
        tmp = eye(m,n)
        F = schurfact(A[end-n+1:end,:],tmp[end-n+1:end,:])
        S = vcat(A[1:m-n,:]*F[:Z],F[:S])
        T = vcat(F[:Z][1:m-n,:],F[:T])
    end
    return S,T
end

function expandlevels(ld::LevelDesc)
    if ld.isunif
        if ld.step != 0
            levels = ld.first:ld.step:ld.last
        else
            levels = ones(2)*ld.first
        end
    else
        levels = ld.full_levels
        (length(levels) == 1) && (levels = levels[1]*ones(2))
    end
    isa(levels,AbstractRange) && (levels = collect(levels))
    levels
end

"""
tidy up axis limits

    tidyaxes(vmin,vmax,tol) -> nmin,nmax

rounds [vmin,vmax] outward to have cleaner decimal expressions
"""
function tidyaxes(vmin::Number,vmax::Number,tol)
    span = vmax - vmin
    round_to = span * tol
    lround_to = log10(round_to)
    expon = floor(lround_to)
    mant = round(10^(lround_to-expon))
    round_to = mant*10^expon
    nmin = round_to*floor(vmin/round_to)
    nmax = round_to*ceil(vmax/round_to)

    return nmin,nmax
end

macro mywarn(args...)
    msg = args[2]
    quote
        @warn($msg)
    end
end


text:
export RidgeRegressor, PCA

import MLJBase
import MLJBase: @mlj_model, metadata_model, metadata_pkg
# using Distances
using LinearAlgebra
using Tables
using ScientificTypes

import MultivariateStats

const MS = MultivariateStats

struct LinearFitresult{F} <: MLJBase.MLJType
    coefficients::Vector{F}
    intercept::F
end

const RIDGE_DESCR = "Ridge regressor with regularization parameter lambda. Learns a linear regression with a penalty on the l2 norm of the coefficients."
const PCA_DESCR = "Principal component analysis. Learns a linear transformation to project the data  on a lower dimensional space while preserving most of the initial variance."

####
#### RIDGE
####

"""
RidgeRegressor(; lambda=1.0)

$RIDGE_DESCR

## Parameters

* `lambda=1.0`: non-negative parameter for the regularization strength.
"""
@mlj_model mutable struct RidgeRegressor <: MLJBase.Deterministic
    lambda::Real = 1.0::(_ ≥ 0)
end

function MLJBase.fit(model::RidgeRegressor, verbosity::Int, X, y)
    Xmatrix   = MLJBase.matrix(X)
    features  = Tables.schema(X).names
    θ         = MS.ridge(Xmatrix, y, model.lambda)
    coefs     = θ[1:end-1]
    intercept = θ[end]

    fitresult = LinearFitresult(coefs, intercept)
    report    = NamedTuple()
    cache     = nothing

    return fitresult, cache, report
end

MLJBase.fitted_params(::RidgeRegressor, fr) =
    (coefficients=fr.coefficients, intercept=fr.intercept)

function MLJBase.predict(::RidgeRegressor, fr, Xnew)
    Xmatrix = MLJBase.matrix(Xnew)
    return Xmatrix * fr.coefficients .+ fr.intercept
end

####
#### PCA
####

const PCAFitResultType = MS.PCA

"""
PCA(; maxoutdim=nothing, method=:auto, pratio=0.99, mean=nothing)

$PCA_DESCR

## Parameters

* `maxoutdim=nothing`: maximum number of output dimensions, unconstrained if nothing.
* `method=:auto`: method to use to solve the problem, one of `:auto`, `:cov` or `:svd`
* `pratio=0.99`: ratio of variance preserved
* `mean=nothing`: if set to nothing centering will be computed and applied, if set to `0` no centering (assumed pre-centered), if a vector is passed, the centering is done with that vector.
"""
@mlj_model mutable struct PCA <: MLJBase.Unsupervised
    maxoutdim::Union{Nothing,Int} = nothing::(_ === nothing || _ ≥ 1)
    method::Symbol  = :auto::(_ in (:auto, :cov, :svd))
    pratio::Float64 = 0.99::(0.0 < _ ≤ 1.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MLJBase.fit(model::PCA, verbosity::Int, X)
    Xarray = MLJBase.matrix(X)
    mindim = minimum(size(Xarray))

    maxoutdim = model.maxoutdim === nothing ? mindim : model.maxoutdim

    # NOTE: copy/transpose
    fitresult = MS.fit(MS.PCA, permutedims(Xarray);
                       method=model.method,
                       pratio=model.pratio,
                       maxoutdim=maxoutdim,
                       mean=model.mean)

    cache = nothing
    report = (indim=size(fitresult)[1],
              outdim=size(fitresult)[2],
              mean=MS.mean(fitresult),
              principalvars=MS.principalvars(fitresult),
              tprincipalvar=MS.tprincipalvar(fitresult),
              tresidualvar=MS.tresidualvar(fitresult),
              tvar=MS.var(fitresult))

    return fitresult, cache, report
end

MLJBase.fitted_params(::PCA, fr) = (projection=fr,)

function MLJBase.transform(::PCA, fr::PCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MLJBase.matrix(X)
    Xnew   = permutedims(MS.predict(fr, permutedims(Xarray)))
    return MLJBase.table(Xnew, prototype=X)
end

####
#### METADATA
####

metadata_pkg.((RidgeRegressor, PCA),
              name="MultivariateStats",
              uuid="6f286f6a-111f-5878-ab1e-185364afe411",
              url="https://github.com/JuliaStats/MultivariateStats.jl",
              license="MIT",
              julia=true,
              is_wrapper=false)

metadata_model(RidgeRegressor,
               input=MLJBase.Table(MLJBase.Continuous),
               target=AbstractVector{MLJBase.Continuous},
               weights=false,
               descr=RIDGE_DESCR)

metadata_model(PCA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=MLJBase.Table(MLJBase.Continuous),
               weights=false,
               descr=PCA_DESCR)


text:
function hNewton(AA::HMat,A0::HMat,numIter=30)
    tol = 1e-2;
    X = hcopy(A0);
    Xold = hscale(-1,X);
    A = hcopy(AA);
    hscale!(-1,A);
    for iter = 1:numIter
        Z = hmul(X,A);
        hadddiag!(Z,2);
        X = hmul(Z,X);
        EMat = hadd(Xold,X);
        err = hnorm(EMat);
        println(iter,": ",err);
        if err < tol
            break;
        end
        Xold = hscale(-1,X);
    end
    return X;
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
# Perform a seesaw optimization searching for a good distillation scheme
# in the vicinity of a known one

using Convex
using SCS

export seesaw, seesawAlice, seesawBob;

"""
`(CA,CB, F, psucc) = seesaw(rhoAB, nA, nB, k, startCA, startCB)` or
`(CA,CB, F, psucc) = seesaw(rhoAB, nA, nB, k, startCA, startCB, P)` or
`(CA,CB, F, psucc) = seesaw(rhoAB, nA, nB, k, startCA, startCB, P, max_iter)` or
`(CA,CB, F, psucc) = seesaw(rhoAB, nA, nB, k, startCA, startCB, P, max_iter, delta)`

Performs a seesaw optimization searching for a good scheme to distill *rhoAB*
into a maximally entangled state of local dimension *k*. Good here means that we want to maximize the fidelity of the output state to the maximally entangled state, where the success probability should be the same as the existing one - or, if supplied, it should achieve a specific success probability *delta*.

Alice's map is assumed to be from
A -> hatA, F_A, and Bob's from B -> hatB, F_B, where F_A and F_B are the flag registered that yield success or failure. After the map *P* is measured on the flag registeres to determine success where *P* is the projector determining success. Example P=|11><11| on the flag registers F_A and F_B.

Input:
- *rhoAB* state to be distilled
- *nA* dimension of A
- *nB* dimension of B
- *k*  local dimension of the desired maximally entangled output state
- *startCA* choi state of Alice's map to start from (ordering of systems: hatA (output state), F_A (flag register), A (input))
- *startCB* choi state of Bob's map to start from (ordering of systems: hatB (output state), F_B (flag register), B (input))
- *P* (optional) projector onto success: Default is |11><11|
- *max_iter* (optional) maximum number of iterations (default 10)
- *delta* (optional) desired success probability if different from that of the starting map.

Output
- *CA* new map for Alice as choi state
- *CB* new map for Bob as choi state
- *F* fidelity obtained
- *psucc* success probability attained
"""

function seesaw(
	rho::AbstractMatrix,
	nA::Int, nB::Int, k::Int,
	startCA::AbstractMatrix, startCB::AbstractMatrix,
	P::AbstractMatrix = [0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1],
	max_iter::Number=10,
	delta::Number=0)

	# Validate the essential elements of the input
	@assert isQuantumState(rho) "Input is not a valid quantum state."
	nTot = nA * nB;
	(d1,d2) = size(rho);
	@assert d1 == nTot "Input state doesn't match indicated dimensions."
	@assert delta <= 1 "Success probability cannot exceed 1."
	@assert delta >= 0 "Success probability must be positive."

	# The output state we want to produce
	want = maxEnt(k);

	# Dimensions of the systems in ordering hatA, FA, A, hatB, FB, B
	dimSys = [k, k, 2, 2, nA, nB];

	# The current probability of success
				M = permutesystems(kron(kron(eye(k^2),P),rho'), [1,3,5,2,4,6], dimSys);
	psuccOrig = nA * nB * trace(M * kron(startCA, startCB));
	bestP = psuccOrig;

	# The current fidelity achieved
				M = permutesystems(kron(kron(want,P),rho'), [1,3,5,2,4,6], dimSys);
	FOrig = nA * nB * trace(M * kron(startCA, startCB))/psuccOrig;
	bestF = FOrig;

	print("\n\nSEESAW Optimization\n\n");
	print("Starting with fidelity F=",round(FOrig,3), " and p_succ=", round(psuccOrig,3),"\n\n");

	# Set the desired probability of success
	if delta > 0
		# Desired delta is supplied: make that the desired acceptable succes probability
		fixedSucc = delta;
	else
		# Delte is not supplied: we want to keep the success probability of the starting scheme
		fixedSucc = psuccOrig;
	end

	# Iteration to search for a better scheme
	currCA = startCA;
	currCB = startCB;

	for j = 1:max_iter

		# A flag whether we achieved an improvement on either side
		betterA = false;
		betterB = false;

		# Fix Bob and find a better map on Alice's side
		(newCA, FA, psuccA) = seesawAlice(rho, nA, nB, k, currCB, P, fixedSucc);

		if (FA > bestF)
			# We improved the fidelity while achieving the desired success prob.
			# Update Alice's map
			currCA = newCA;
			betterA = true;
			bestF = FA;
			bestP = psuccA;
		end

		# Fix Alice and find a better map on Bob's side
		(newCB, FB, psuccB) = seesawBob(rho, nA, nB, k, currCA, P, fixedSucc);

		if (FB > bestF)
			# We improved the fidelity while maintaining the desired success prob.
			# Update Bob's map
			currCB = newCB;
			betterB = true;
			bestF = FB;
			bestP = psuccB;
		end

		# Check whether any improvement was achieved. Otherwise just return
		if !(betterA || betterB)
			print("\nStopping early after ",j," iterations.\n");
			print("Final fidelity F=",round(bestF,3), " and psucc=", round(bestP,3),"\n");
			return(currCA, currCB, bestF, bestP);
		end
	end
	print("\nFinal fidelity F=",round(bestF,3), " and psucc=", round(bestP,3),"\n");

	return(currCA, currCB, bestF, bestP);
end

""" `(newCA, FA, psuccA) = seesawAlice(rhoAB, nA, nB, k, CB, P, fixedSucc)`

Searches for a better distillation map for Alice, optimizing the product of the fidelity and success probability, for a fixed desired success probability. Called from seesaw.

Input (assumed to be valid):
- *rhoAB* state to be distilled
- *nA* dimension of A
- *nB* dimension of B
- *k* local dimension of desired max. entangled state
- *CB* choi state of Bob's map
- *P* measurement indicating success
- *fixedSucc * desired probability of success

Output:
- *newCA* new map for Alice as choi state
- *F* fidelity of new scheme
- *psucc* success probability of new scheme
"""

function seesawAlice(
	rho::AbstractMatrix,
	nA::Int, nB::Int, k::Int,
	CB::AbstractMatrix,
	P::AbstractMatrix,
	fixedSucc::Number)

	# Dimension of the maps
	d = nA * 2 * k;

	# Desired output state
	want = maxEnt(k);

	# Check whether we want a real or complex optimization problem
	if (real(rho) == rho) && (real(CB) == CB)
		CA = Semidefinite(d);
	else
		CA = ComplexVariable(d,d);
	end

	# Dimensions in ordering hatA, Fa, A,...
	dimSys = [k, k, 2, 2, nA, nB];

	# Set up the optimization problem
	# The current probability of success
				Mwant = permutesystems(kron(kron(want,P),rho'), [1,3,5,2,4,6], dimSys);
	tv = Variable(1);
	problem = maximize(tv);
	problem.constraints += tv == nA * nB * trace(Mwant * simplekron(CA, CB));
	problem.constraints += CA in :SDP;

	# Constraints
	# CA should be a choi state
	problem.constraints += ([trace(CA) == 1]);

	CAin = partialtrace(CA,[1],[2 * k, nA]);
	problem.constraints += ([CAin == eye(nA)/nA]);

	# We want a fixed success probability
				Mid = permutesystems(kron(kron(eye(k^2),P),rho'), [1,3,5,2,4,6], dimSys);
	problem.constraints += ([nA * nB * trace(Mid * simplekron(CA, CB)) == fixedSucc]);

	# Solve the SDP
	solve!(problem, SCSSolver(verbose=false));

	# Compute the new success probability
	psuccNew = nA * nB * trace(Mid * kron(CA.value, CB));

	# The current fidelity achieved
	Fnew = nA * nB * trace(Mwant * kron(CA.value, CB))/psuccNew;

	return(CA.value, Fnew, psuccNew);
end

""" `(newCB, FB, psuccB) = seesawBob(rhoAB, nA, nB, k, CA, P, fixedSucc)`

Searches for a better distillation map for Bob, optimizing the product of the fidelity and success probability, for a fixed desired success probability. Called from seesaw.

Input (assumed to be valid):
- *rhoAB* state to be distilled
- *nA* dimension of A
- *nB* dimension of B
- *k* local dimension of desired max. entangled state
- *CA* choi state of Alice's map
- *P* measurement indicating success
- *fixedSucc* desired probability of success

Output:
- *newCB* new map for Bob as choi state
- *F* fidelity of new scheme
- *psucc* success probability of new scheme
"""

function seesawBob(
	rho::AbstractMatrix,
	nA::Int, nB::Int, k::Int,
	CA::AbstractMatrix,
	P::AbstractMatrix,
	fixedSucc::Number)

	# Dimension of the maps
	d = nB * 2 * k;

	# Desired output state
	want = maxEnt(k);

	# Check whether we want a real or complex optimization problem
	if (real(rho) == rho) && (real(CA) == CA)
		CB = Semidefinite(d);
	else
		CB = ComplexVariable(d,d);
	end

	# Dimensions in ordering hatA, Fa, A,...
	dimSys = [k, k, 2, 2, nA, nB];

	# Set up the optimization problem
	# The current probability of success
	Mwant = permutesystems(kron(kron(want,P),rho'), [1,3,5,2,4,6], dimSys);
	tv = Variable(1);
	problem = maximize(tv);
	problem.constraints += tv == nA * nB * trace(Mwant * kron(CA, CB));
	problem.constraints += CB in :SDP;

	# Constraints
	# CA should be a choi state
	problem.constraints += ([trace(CB) == 1]);

	CBin = partialtrace(CB,[1],[2 * k, nB]);
	problem.constraints += ([CBin == eye(nA)/nA]);

	# We want a fixed success probability
	Mid = permutesystems(kron(kron(eye(k^2),P),rho'), [1,3,5,2,4,6], dimSys);
	problem.constraints += ([nA * nB * trace(Mid * kron(CA, CB)) == fixedSucc]);

	# Solve the SDP
	solve!(problem, SCSSolver(verbose=false));

	# Compute the new success probability
	psuccNew = nA * nB * trace(Mid * kron(CA, CB.value));

	# The current fidelity achieved
	Fnew = nA * nB * trace(Mwant * kron(CA, CB.value))/psuccNew;

	return(CB.value, Fnew, psuccNew);
end


text:
#####################
# Handle properties #
#####################
abstract type MetaPropertiesArg end
struct AppendProperties <: MetaPropertiesArg end
struct ReplaceProperties <: MetaPropertiesArg end
struct IgnoreProperties <: MetaPropertiesArg end

function updateproperties!(::Type{AppendProperties}, s::AbstractMetaStream, sink::ImageMeta, propkeys::SVector{0,String})
    append!(sink.properties, properties(s))
    return nothing
end
function updateproperties!(::Type{AppendProperties}, s::AbstractMetaStream, sink::ImageMeta, propkeys::SVector{N,String}) where {N}
    for i in propkeys
        append!(properties(sink)[i], properties(s)[i])
    end
end

function updateproperties!(::Type{ReplaceProperties}, s::AbstractMetaStream, sink::ImageMeta, propkeys::SVector{0,String})
    properties(sink) = properties(s)
    return nothing
end
function updateproperties!(::Type{ReplaceProperties}, s::AbstractMetaStream, sink::ImageMeta, propkeys::SVector{N,String}) where {N}
    for i in propkeys
        properties(sink)[i] = properties(s)[i]
    end
    return nothing
end

updateproperties!(::Type{IgnoreProperties}, s::AbstractMetaStream, sink::ImageMeta, propkeys::SVector) = nothing

########
# Read #
########
function Base.read!(s::S, sink::ImageMeta; propchange::MetaPropertiesArg, propkeys::SVector{N,String}) where {S<:AbstractMetaStream,N}
   read!(getfield(s,:data), getfield(s,:data))
   updateproperties!(s, sink, propchange)
   return sink
end

function Base.read(s::S) where {S<:AbstractMetaStream}
    read!(s, ImageMeta(AxisArray(Array{eltype(s),ndims(s)}(undef, size(s)), AxisArrays.axisparams(s)),
                       properties(s)); propchange=IgnoreProperties, propkeys=SVector{0,String}())
end

#########
# Write #
#########
Base.write(s::S, sink::ImageMeta) where {S<:AbstractMetaStream} = write(getfield(s, :data), getfield(sink, :data))