text:
#= # for local testing, need this:
using Test
using PhyloNetworks
using Random
=#

@testset "unconstrained NNI moves" begin

str_level1 = "(((S8,S9),(((((S1,S2,S3),S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
net_level1 = readTopology(str_level1); # this network has a polytomy at node -9
# same topology as: rootatnode!(net_level1, -3). edges 1:22
str_nontreechild = "((((Ag,E))#H3,(#H1:7.159::0.056,((M:0.0)#H2:::0.996,(Ak,(#H3:0.08,#H2:0.0::0.004):0.023):0.078):2.49):2.214):0.026,((Az:2.13,As:2.027):1.697)#H1:0.0::0.944,Ap);"
net_nontreechild = readTopology(str_nontreechild);
# problem: the plot has an extra vertical segment, for a clade that's not in the major tree
# --> fix that in PhyloPlots (fixit)
str_hybridladder = "(#H2:::0.2,((C,((B)#H1)#H2:::0.8),(#H1,(A1,A2))),O);"
net_hybridladder = readTopology(str_hybridladder);

#=
using PhyloPlots
plot(net_level1, :R, showNodeNumber=true, showEdgeNumber=true)
plot(net_nontreechild, :R, showNodeNumber=true, showEdgeNumber=true)
plot(net_hybridladder, :R, showNodeNumber=true, showEdgeNumber=true)
=#

@test isnothing(nni!(net_level1, net_level1.edge[1], 0x01, true, true)) # external edge

@testset "level1 edge 3: BB undirected move $move" for move in 0x01:0x08
    undoinfo = nni!(net_level1, net_level1.edge[3], move, true, true);
    #location of v node (number -4)
    nodes = [n.number for n in net_level1.edge[20].node] #α's connections
    if move in [1, 2, 6, 8] #check that edge α connected to v
        @test -4 in nodes
    else
        @test !(-4 in nodes)
    end
    #location of u node (number -3)
    nodes = [n.number for n in net_level1.edge[1].node]
    if move in [2, 4, 5, 6] #check that edge γ connected to u
        @test -3 in nodes
    else
        @test !(-3 in nodes)
    end
    # check directionality
    if move in [3, 4, 5, 7]
        # keeping α, β or flipping uv keeps node -4 as child of edge 3
        @test PhyloNetworks.getChild(net_level1.edge[3]).number == -4
    else
        # switching α, β AND flipping uv or doing neither makes node -3 child of edge 3
        @test PhyloNetworks.getChild(net_level1.edge[3]).number == -3
    end
    # undo move
    nni!(undoinfo...);
    # confirm we're back to original topology
    @test writeTopology(net_level1) == str_level1
end # of level1 edge 3: BB undirected

@testset "level1 test edge 13: BR directed move $move" for move in 0x01:0x03
    # e.hybrid and tree parent:  BR case, 3 moves because e cannot contain the root
    # 3cycle test: α connected to γ so move 1 will create 3 cycles.
    if move == 0x01
        @test isnothing(nni!(net_level1, net_level1.edge[13], move, true, true)) # would create a 3cycle
    # move 1 would create a 3 cycle, but should work if we don't forbid 3cycles
        undoinfo = nni!(net_level1, net_level1.edge[13], move, true, false); #no3cycle=false
        nodes = [n.number for n in net_level1.edge[17].node]
        @test 8 in nodes # check that edge α connected to v
        nodes = [n.number for n in net_level1.edge[11].node]
        @test !(-11 in nodes)
        @test PhyloNetworks.getChild(net_level1.edge[13]).number == -11
        nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == "(((S8,S9),(((((S1,S2,S3),S4),(S5:0.0)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
        # confirm we're back to original topology (but now with 0.0 branch lengths)
    else
        undoinfo = nni!(net_level1, net_level1.edge[13], move, true, false);
        # test that move was made
        nodes = [n.number for n in net_level1.edge[17].node]
        if move == 0x03
            @test 8 in nodes # check that edge α connected to v
        else
            @test !(8 in nodes)
        end
        #location of u node (number -11)
        nodes = [n.number for n in net_level1.edge[11].node]
        if move == 0x03 # check that edge γ connected to u
            @test -11 in nodes
        else # check that edge δ connected to u
            @test !(-11 in nodes)
        end
        # check directionality
        @test PhyloNetworks.getChild(net_level1.edge[13]).number == 8
        nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == "(((S8,S9),(((((S1,S2,S3),S4),(S5:0.0)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"# confirm we're back to original topology
    end
end # of level1 edge 13: BR directed

@testset "level1 edge 16: BB directed move $move" for move in 0x01:0x02
    # e not hybrid, tree parent:  BB case, 2 NNIs if directed, 8 if undirected
    undoinfo = nni!(net_level1, net_level1.edge[16], move, true, true);
    nodes = [n.number for n in net_level1.edge[13].node] #β's connections
    @test -11 in nodes #check β connected to u node (number -11)
    nodes = [n.number for n in net_level1.edge[14].node] #γ's connections
    if move == 1 #check that edge γ connected to v
        @test -12 in nodes
    else #check that edge γ connected to u
        @test -11 in nodes
    end
    # check directionality node -11 child of edge 16 in both cases
    @test PhyloNetworks.getChild(net_level1.edge[16]).number == -11
    nni!(undoinfo...); # undo move
    @test writeTopology(net_level1) == "(((S8,S9),(((((S1,S2,S3),S4),(S5:0.0)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
        # confirm we're back to original topology
end #of level1 edge 16: BB directed

# BB directed has two allowed moves, so 0x03 should throw an exception
@test_throws Exception nni!(net_level1, net_level1.edge[16], 0x03);

@testset "level1 edge 17: BB directed with 4cycle move $move" for move in 0x01:0x02
    # no3cycle test
        # move 2 should fail because of 3cycle problems
        # β is connected in a 4 cycle with γ so a move that makes γ and β a pair
        # would create a 3 cycle
    if move == 0x01
        undoinfo = nni!(net_level1, net_level1.edge[17], move, true, true);
        nodes = [n.number for n in net_level1.edge[12].node] #β's connections
        @test -6 in nodes #check β connected to u node (number -6)
        nodes = [n.number for n in net_level1.edge[13].node] #γ's connections
        #check that edge γ connected to v
        @test -11 in nodes
        # check directionality node -6 child of edge 17 in both cases
        @test PhyloNetworks.getChild(net_level1.edge[17]).number == -6
        nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == "(((S8,S9),(((((S1,S2,S3),S4),(S5:0.0)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
    else
        @test isnothing(nni!(net_level1, net_level1.edge[17], move, true, true))
        undoinfo = nni!(net_level1, net_level1.edge[17], move, true, false) #should work if we dont check for 3cycles
        nodes = [n.number for n in net_level1.edge[12].node] #β's connections
        @test -6 in nodes #check β connected to u node (number -6)
        nodes = [n.number for n in net_level1.edge[13].node] #γ's connections
        #check that edge γ connected to u
        @test -6 in nodes
        # check directionality node -11 child of edge 17 in both cases
        @test PhyloNetworks.getChild(net_level1.edge[17]).number == -6
        nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == "(((S8,S9),(((((S1,S2,S3),S4),(S5:0.0)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
    end
end # of level1 edge 17: BB directed

@testset "level1 edge 18: RR (directed) move $move" for move in 0x01:0x04
    # RB case, 4 moves. uv edge cannot contain the root (always directed)
    undoinfo = nni!(net_level1, net_level1.edge[18], move, true, true);
    #test that move was made
    nodes = [n.number for n in net_level1.edge[19].node] #α's connections
    if move in [0x01, 0x03]
        @test -6 in nodes #check that edge α connected to v
    elseif move in [0x02, 0x04]
        @test !(-6 in nodes)
    end
    nodes = [n.number for n in net_level1.edge[12].node] #γ's connections
    if move in [0x03, 0x04] #check that edge γ connected to u (11)
        @test 11 in nodes
    elseif move in [0x01, 0x02]  #check that edge δ connected to u, not γ
        @test !(11 in nodes)
    end
    #check directionality (should point toward u, node 11)
    @test PhyloNetworks.getChild(net_level1.edge[18]).number == 11
    nni!(undoinfo...);
        # edge constrained at 0.0 now
    @test writeTopology(net_level1) == "(((S8,S9),(((((S1,S2,S3),S4),(S5:0.0)#H1),(#H1,(S6,S7))):0.0)#H2),(#H2,S10));"
end # of level1 edge 18: RR (directed)

@testset "non tree child net edge 3: RB (directed) move $move" for move in 0x01:0x04
    # RB case, 4 moves. uv edge cannot contain the root
    undoinfo = nni!(net_nontreechild, net_nontreechild.edge[3], move, true, true);
    #test that move was made
    #location of v node (node -5)
    nodes = [n.number for n in net_nontreechild.edge[4].node] #α's connections
    if move in [0x01, 0x03]
        #check that edge α connected to v
        @test -5 in nodes
    else
        #check that edge β connected to v, not u
        @test !(-5 in nodes)
    end
    nodes = [n.number for n in net_nontreechild.edge[2].node] #δ's connections
    if move in [0x01, 0x02] #check that edge δ connected to u
        @test 3 in nodes
    else
        @test !(3 in nodes)
    end
    #check directionality (should point toward u, node 3)
    @test PhyloNetworks.getChild(net_nontreechild.edge[3]).number == 3
    #undo move
    nni!(undoinfo...);
        # keep constrained edges at 0.0, but otherwise topology completely restored
    @test writeTopology(net_nontreechild) == "((((Ag,E):0.0)#H3,(#H1:7.159::0.056,((M:0.0)#H2:::0.996,(Ak,(#H3:0.08,#H2:0.0::0.004):0.023):0.078):2.49):2.214):0.026,((Az:2.13,As:2.027):1.697)#H1:0.0::0.944,Ap);"
end #of non tree child net edge 5: RB (directed)

@testset "hybrid ladder net edge 1: BR undirected (u at root, potential nonDAG, 3cycle) move $move" for move in 0x01:0x06
    # BR case, 6 moves. uv edge can contain the root. u at root
    # DAG check: u is root, so if α -> γ, moves 1 and 3 will create a nonDAG
        # α -> γ so moves 1 and 3 will fail.
        # moves 1 and 3 in the notes correspond to 1, 1' and 3, 3' (1, 4 and 3, 6)
    # 3cycle check: could create 3 cycle (α connected to γ) so moves 1 and 5 forbidden
    if move in [0x01, 0x03, 0x05, 0x06] # DAG check
        @test isnothing(nni!(net_hybridladder, net_hybridladder.edge[1], move, false, true))
    end
    if move in [0x01, 0x05] # 3cycle check
        @test isnothing(nni!(net_hybridladder, net_hybridladder.edge[1], move, true, true))
    end
    if move == 0x04
        undoinfo = nni!(net_hybridladder, net_hybridladder.edge[1], move, false, true);
        nodes = [n.number for n in net_hybridladder.edge[12].node] # α
        @test !(1 in nodes) # check that edge α is connected to v
        nodes = [n.number for n in net_hybridladder.edge[4].node] # δ's connections
        @test -2 in nodes # check that edge δ is connected to u
        # check directionality (edge should point toward u, node -2)
        @test PhyloNetworks.getChild(net_hybridladder.edge[1]).number == 1
        nni!(undoinfo...);
        @test writeTopology(net_hybridladder)== "(#H2:::0.2,((C,((B)#H1:0.0)#H2:::0.8),(#H1,(A1,A2))),O);" # restored but edge below hybrid node constrained at 0.0
    elseif move == 0x02
        undoinfo = nni!(net_hybridladder, net_hybridladder.edge[1], move, true, true);
        nodes = [n.number for n in net_hybridladder.edge[12].node] # α
        @test !(1 in nodes) # check that edge α not connected to v
        nodes = [n.number for n in net_hybridladder.edge[4].node] # δ's connections
        @test -2 in nodes # δ connected to u
        # check directionality (edge should point toward u, node -2)
        @test PhyloNetworks.getChild(net_hybridladder.edge[1]).number == 1
        nni!(undoinfo...);
        @test writeTopology(net_hybridladder) == "(#H2:::0.2,((C,((B)#H1:0.0)#H2:::0.8),(#H1,(A1,A2))),O);" # restored but edge below hybrid node constrained at 0.0
    end
end # of hybrid ladder net edge 1: BR undirected

@testset "hybrid ladder net edge 4: RR (directed) move $move" for move in 0x01:0x02
    # RR case, 2 moves. uv edge cannot contain the root (always directed)
    undoinfo = nni!(net_hybridladder, net_hybridladder.edge[4], move, false, true);
    #test that move was made
    nodes = [n.number for n in net_hybridladder.edge[5].node] #α's connections
    if move == 0x01 #check that edge α connected to v
        @test 4 in nodes
    else move == 0x02 #check that edge β connected to v, not α
        @test !(4 in nodes)
    end
    #check that edge δ connected to u, not γ
    nodes = [n.number for n in net_hybridladder.edge[3].node]
    @test 1 in nodes
    #check directionality (should point toward u, node 1)
    @test PhyloNetworks.getChild(net_hybridladder.edge[4]).number == 1
    #undo move
    nni!(undoinfo...);
    @test writeTopology(net_hybridladder) == "(#H2:::0.2,((C,((B)#H1:0.0)#H2:::0.8),(#H1,(A1,A2))),O);" # restored but edge below hybrid node constrained at 0.0
end #of hybrid ladder net edge 4: RR (directed)

@testset "hybrid ladder net edge 5: BR undirected move $move" for move in 0x01:0x06
    # BR case, 6 moves. uv edge can contain the root
    # 3 cycle test: α connected to γ, α -> u
    #   moves 1, 5 forbidden
    # DAG test:
    #   no path from α -> γ or β -> γ so all moves should work
    if move in [0x01, 0x05]
        @test isnothing(nni!(net_hybridladder, net_hybridladder.edge[5], move, false, true)) # 3-cycles forbidden
        undoinfo = nni!(net_hybridladder, net_hybridladder.edge[5], move, false, false) # 3-cycles allowed
    else
        undoinfo = nni!(net_hybridladder, net_hybridladder.edge[5], move, false, true);
    end
    #tests for all moves:
    nodes = [n.number for n in net_hybridladder.edge[6].node] #α
    if move in [0x01, 0x03, 0x05, 0x06]
        @test 1 in nodes #check that edge α connected to v
    else
        @test !(1 in nodes) #check that edge α not connected to v
    end
    nodes = [n.number for n in net_hybridladder.edge[4].node] #δ's connections
    if move in [0x01, 0x02, 0x04, 0x05]
        @test -4 in nodes #δ connected to u
    else
        @test !(-4 in nodes)
    end
    #check directionality
    if move in [0x01, 0x05]
        #(edge should point toward u, node -4)
        @test PhyloNetworks.getChild(net_hybridladder.edge[5]).number == -4
    else
        @test PhyloNetworks.getChild(net_hybridladder.edge[5]).number == 1
    end
    #undo move
    nni!(undoinfo...);
    @test writeTopology(net_hybridladder) == "(#H2:::0.2,((C,((B)#H1:0.0)#H2:::0.8),(#H1,(A1,A2))),O);" # restored but edge below hybrid node constrained at 0.0
end #of hybrid ladder net edge 5: BR undirected

@testset "hybrid ladder net edge 12: BB undirected (edge below root) move $move" for move in 0x01:0x08
    # BB case, 8 moves. uv edge can contain the root. no flip.
    # no3cycle: moves 1, 4, 5, 8 would create a 3 cycle because α is connected to γ
    if move in [0x01, 0x04, 0x05, 0x08]
        @test isnothing(nni!(net_hybridladder, net_hybridladder.edge[12], move, true, true))
        undoinfo = nni!(net_hybridladder, net_hybridladder.edge[12], move, true, false);
    else
        undoinfo = nni!(net_hybridladder, net_hybridladder.edge[12], move, true, true);
    end
    nodes = [n.number for n in net_hybridladder.edge[1].node] #α
    if move in [0x01, 0x02, 0x06, 0x08]
        @test -3 in nodes #check that edge α connected to v
    else
        @test !(-3 in nodes) #check that edge α not connected to v
    end
    nodes = [n.number for n in net_hybridladder.edge[11].node] #δ's connections
    if move in [0x01, 0x03, 0x07, 0x08]
        @test -2 in nodes #δ connected to u
    else
        @test !(-2 in nodes)
    end
    #check directionality
    @test PhyloNetworks.getChild(net_hybridladder.edge[12]).number == -3
    #undo move
    nni!(undoinfo...);
    @test writeTopology(net_hybridladder) == "(#H2:::0.2,((C,((B)#H1:0.0)#H2:::0.8),(#H1,(A1,A2))),O);" # restored but edge below hybrid node constrained at 0.0
end # of hybrid ladder net edge 12: BB undirected (edge below root)

@testset "test isdescendant and isconnected functions" begin
    net_level1 = readTopology(str_level1);
    @test  PhyloNetworks.isdescendant(net_level1.node[7], net_level1.node[17])  # nodes -9, -6
    @test !PhyloNetworks.isdescendant(net_level1.node[7], net_level1.node[3])   # nodes -9, -4
    @test  PhyloNetworks.isdescendant(net_level1.node[15], net_level1.node[17]) # nodes -12, -6
    @test !PhyloNetworks.isdescendant(net_level1.node[12], net_level1.node[12])
    @test  PhyloNetworks.isconnected(net_level1.node[12], net_level1.node[17])  # nodes -7, -6
    @test !PhyloNetworks.isconnected(net_level1.node[12], net_level1.node[19])  # nodes -7, -3
    # mess up the direction of some tree edges, then check descendence relationships with isdescendant_undirected
    for i in [4,5,6,7,9,10,12,17,3,20] net_level1.edge[i].isChild1 = !net_level1.edge[i].isChild1; end
    @test  PhyloNetworks.isdescendant_undirected(net_level1.node[7], net_level1.node[17], net_level1.edge[18])
    @test !PhyloNetworks.isdescendant_undirected(net_level1.node[7], net_level1.node[3], net_level1.edge[3])
    @test  PhyloNetworks.isdescendant_undirected(net_level1.node[15], net_level1.node[17], net_level1.edge[18])
    @test !PhyloNetworks.isdescendant_undirected(net_level1.node[12], net_level1.node[12], net_level1.edge[12])
end

end # of testset on unconstrained NNIs

@testset "constrained NNI moves" begin
# subsets for: species constraints; move root (species & clade constraints);
#              clade constraints (not yet: TODO)

@testset "species constraints" begin # multiple individuals from each species
str_level1_s = "(((S8,S9),((((S1,S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));" # indviduals S1A S1B S1C go on leaf 1
net_level1_s = readTopology(str_level1_s)
#=
plot(net_level1_s, :R, showNodeNumber=true, showEdgeNumber=true)
=#

# test breakedge! function
newnode, newedge = PhyloNetworks.breakedge!(net_level1_s.edge[4], net_level1_s);
@test length(net_level1_s.node) == 20
@test length(net_level1_s.edge) == 21
@test newnode.edge[1].number == 4
@test newnode.edge[2].number == 21
@test PhyloNetworks.getParent(net_level1_s.edge[4]) === newnode
@test PhyloNetworks.getChild(newedge) === newnode

# test addleaf! function
net_level1_s = readTopology(str_level1_s)
PhyloNetworks.addleaf!(net_level1_s, net_level1_s.node[4], "S1A");
@test !net_level1_s.node[findfirst([n.number == 3 for n in net_level1_s.node])].leaf
PhyloNetworks.addleaf!(net_level1_s, net_level1_s.node[4], "S1B");
PhyloNetworks.addleaf!(net_level1_s, net_level1_s.node[4], "S1C");
@test net_level1_s.edge[21].containRoot == false # check containRoot on edge 4 and exterior edges
@test net_level1_s.edge[22].containRoot == false
@test PhyloNetworks.getChild(net_level1_s.edge[21]).name == "S1A"
@test PhyloNetworks.getChild(net_level1_s.edge[22]).name == "S1B"
# test addleaf! on edge
net_level1_s = readTopology(str_level1_s)
PhyloNetworks.addleaf!(net_level1_s, net_level1_s.edge[4], "S1A");
@test length(net_level1_s.node) == 21
@test net_level1_s.node[21].leaf

# test addindividuals! function
net_level1_s = readTopology(str_level1_s)
PhyloNetworks.addindividuals!(net_level1_s, "S1", ["S1A", "S1B", "S1C"])
@test !net_level1_s.node[findfirst([n.number == 3 for n in net_level1_s.node])].leaf
@test length(net_level1_s.node[findfirst([n.number == 3 for n in net_level1_s.node])].edge) == 4
# spaces in name
net_level1_s = readTopology(str_level1_s)
@test_logs (:warn, r"^species S 1 not") PhyloNetworks.addindividuals!(net_level1_s, "S 1", ["S1A", "S1B", "S1C"])
@test writeTopology(net_level1_s) == str_level1_s # network unchanged
@test_logs (:warn, r"^Spaces in \"S1 A\" may cause errors") PhyloNetworks.addindividuals!(net_level1_s, "S1", ["S1 A", "S1B", "S1C"])
@test writeTopology(net_level1_s) == "(((S8,S9),(((((S1_A,S1B,S1C)S1,S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
# # test mapindividuals! function
net_level1_s = readTopology(str_level1_s)
# in net env
filename = joinpath(@__DIR__, "..","examples","mappingIndividuals.csv")
# filename = abspath(joinpath(dirname(Base.find_package("PhyloNetworks")), "..", "examples", "mappingIndividuals.csv"))
net_level1_i, c_species = PhyloNetworks.mapindividuals(net_level1_s, filename)
@test string(c_species[1]) == "Species constraint, on tips: S1A, S1B, S1C\n stem edge number 4\n crown node number 3"
@test c_species[1].taxonnames == ["S1A","S1B","S1C"]
@test c_species[1].taxonnums == Set([11,12,13])
@test writeTopology(net_level1_i) == "(((S8,S9),(((((S1A,S1B,S1C)S1,S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"

# test clade constraint contructor
c_clade = PhyloNetworks.TopologyConstraint(0x02, ["S1A","S1B","S1C","S4"], net_level1_i)
@test string(c_clade) == "Clade constraint, on tips: S1A, S1B, S1C, S4\n stem edge number 6\n crown node number -8"

# test errors in species constructor
@test_throws ErrorException PhyloNetworks.TopologyConstraint(0x02, ["S1A"], net_level1_i) # only 1 tip
@test_throws ErrorException PhyloNetworks.TopologyConstraint(0x02, ["S1A", "TypoTaxa"], net_level1_i) # typo
@test_throws ErrorException PhyloNetworks.TopologyConstraint(0x02, ["S1A", "S8"], net_level1_i) # not a clade

# NNIs under species constraints
Random.seed!(1234);
# no nni on stem edge for species example
@test isnothing(nni!(net_level1_i , net_level1_i.edge[4], true, true, c_species))
@testset "NNI, 1 species constraint, net level 1, edge $ei" for ei in [8,3,9]
    # 8: BR directed, 3: BB undirected, 9: BB directed, 15: RB directed
    # note: there are no cases of RR directed in net_level1_i
    undoinfo = nni!(net_level1_i , net_level1_i.edge[ei], true, true, c_species);
    @test undoinfo !== nothing
    nni!(undoinfo...);
        # restored to original network, except that edges below hybrid nodes will now have length 0.0
    @test writeTopology(net_level1_i) == "(((S8,S9),(((((S1A,S1B,S1C)S1,S4),(S5:0.0)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
end
# TODO: BR case edge 8: nni move 3 causes problems. hybrid node 6 has 0 or 2+ major hybrid parents
end # of species constraints

@testset "test move root & constraint checking under species & clade constraints" begin
# "(((S8,S9),(((((S1A,S1B,S1C)S1,S4),#H1),((S5)#H1,(S6,S7))))#H2),(#H2,S10));"
netl1_i = readTopology("(((((S1A,S1B,S1C)S1,S4),#H1),((S5)#H1,(S6,S7))));")
con = [PhyloNetworks.TopologyConstraint(0x01, ["S1A","S1B","S1C"], netl1_i),
       PhyloNetworks.TopologyConstraint(0x02, ["S5","S6","S7"], netl1_i)]
Random.seed!(765);
@test PhyloNetworks.moveroot!(netl1_i, con) # only 2 options
writeTopology(netl1_i) == "(((S1A,S1B,S1C)S1,S4),#H1,(((S5)#H1,(S6,S7))));" # now unrooted
@test PhyloNetworks.moveroot!(netl1_i, con) # only 1 option
writeTopology(netl1_i) == "((S1A,S1B,S1C)S1,S4,(#H1,(((S5)#H1,(S6,S7)))));"
netl1_i.root = 14 # back to original rooted network. This node is still of degree 2
@test !PhyloNetworks.checkspeciesnetwork!(netl1_i, con) # false: root *at* clade crown
@test netl1_i.root == 13 # now unrooted (via removedegree2nodes!), root was moved, con[2] stem edge was deleted too...
netl1_i.root = 7; directEdges!(netl1_i) # move root strictly above clade crown
con[2] = PhyloNetworks.TopologyConstraint(0x02, ["S5","S6","S7"], netl1_i)
@test PhyloNetworks.checkspeciesnetwork!(netl1_i, con) # now fine: root *above* clade crown
undoinfo = nni!(netl1_i,netl1_i.edge[8],0x01,false,false);
@test !PhyloNetworks.checkspeciesnetwork!(netl1_i, con)
nni!(undoinfo...);
@test PhyloNetworks.checkspeciesnetwork!(netl1_i, con)
undoinfo = nni!(netl1_i,netl1_i.edge[8],0x03,false,false) # creates a 2-cycle
@test netl1_i.numEdges == 13
PhyloNetworks.deletehybridedge!(netl1_i, netl1_i.edge[10])
@test netl1_i.numEdges == 10 # 2-cycle removed
netl1_i = readTopology("(((S1A,S1B,S1C),S4),#H1,((S5)#H1,(S6,S7)));")
undoinfo = nni!(netl1_i,netl1_i.edge[12],0x03,false,false) # 4-cycle now
@test nni!(netl1_i,netl1_i.edge[12],0x02,true,true) === nothing # would create a 3-cycle
end

#=
str_level1 = "(((S8,S9),(((((S1,S2,S3),S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
net_level1 = readTopology(str_level1); #polytomy at node -9 for leaves 3, 4, 5

str_nontreechild = "((((Ag,E))#H3,(#H1:7.159::0.056,((M:0.0)#H2:::0.996,(Ak,(#H3:0.08,#H2:0.0::0.004):0.023):0.078):2.49):2.214):0.026,((Az:2.13,As:2.027):1.697)#H1:0.0::0.944,Ap);"
net_nontreechild = readTopology(str_nontreechild);

str_polytomy_species = "(((S8,S9),(((((S1,S2,S3),S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));"
net_species = readTopology(str_polytomy_species);
=#

#=
@testset "clade contraints" begin
c_nontree_cladebelowhybrid = PhyloNetworks.TopologyConstraint(0x02, ["Az", "As"], net_nontreechild)
c_level1_clade = PhyloNetworks.TopologyConstraint(0x02, ["S1", "S2", "S3"], net_level1)
# this test below returns nothing but shouldn't
nni!(net_nontreechild, net_nontreechild.edge[20], true, [c_nontree_cladebelowhybrid])
# this test errors
@test isnothing(nni!(net_nontreechild, net_nontreechild.edge[18], true, [c_nontree_cladebelowhybrid]))

@test PhyloNetworks.checkspeciesnetwork!(net_nontreechild, [c_nontree_cladebelowhybrid])

PhyloNetworks.addindividuals!(net_level1_s, "S1", ["S1A", "S1B", "S1C"])
@test_throws ErrorException PhyloNetworks.checkspeciesnetwork!(net_level1_i, [c_species])
@test PhyloNetworks.cladesviolated(net_level1_i, c_species)
@test_throws ErrorException PhyloNetworks.checkspeciesnetwork!(net_level1_s, [c_level1_species])
@test PhyloNetworks.checknetwork(net_level1_i, c_species) #TODO will want to remove if we remove this function
end # of testset on checknetwork functions for species constraints

=#

end # of constrained NNI moves


text:
using Base.Test

# Setup everything.
include("CloudGraphSetup.jl")

# Build a basic type
# And check that if we encode and decode this type, it's exactly the same.
# Make a packed data test structure.
fullType = DataTest(rand(10,10), "This is a test string", rand(Int32,10,10));
# Creating a local test graph.
print("[TEST] Creating a CloudVertex from an ExVertex...")
localGraph = graph(ExVertex[], ExEdge{ExVertex}[]);
#Make an ExVertex that may be encoded
v = make_vertex(localGraph, "TestVertex");
vertex = Graphs.add_vertex!(localGraph, v);
vertex.attributes["data"] = fullType;
vertex.attributes["age"] = 64;
vertex.attributes["latestEstimate"] = [0.0,0.0,0.0];
bigData = CloudGraphs.BigData();
vertex.attributes["bigData"] = bigData;
testElement = CloudGraphs.BigDataElement("Performance test dataset.", Array(UInt8, 1000000)); #1Mb data element
append!(bigData.dataElements, [testElement]);
# Now encoding the structure to CloudGraphs vertex
cloudVertex = CloudGraphs.exVertex2CloudVertex(vertex);

using Mongoc
# function localscope(N=1000000)
  data = ("testId" => Array(UInt8, 10), "description" => "DESCRIPTION")
  @time m_oid = insert(cloudGraph.mongo.cgBindataCollection, data)
  # nothing
# end
myFavouriteKey = first(find(cloudGraph.mongo.cgBindataCollection, ("_id" => eq(m_oid))));

loopCount = 10;
vertices = [];
times = Float32[];
print("[TEST] Testing with $(loopCount) nodes...")

for i = 1:loopCount
  newNode = deepcopy(cloudVertex);
  append!(vertices, [newNode]);
end
# Add a vertex
for i = 1:loopCount
  tic();
  CloudGraphs.add_vertex!(cloudGraph, vertices[i]);
  append!(times, toc());
end
println("Time to add a single vertex = $(mean(times) * 1000.0)ms +- $(std(times)*1000.0)ms");

# Get a vertex
times = Float32[];
for i = 1:loopCount
  tic();
  CloudGraphs.get_vertex(cloudGraph, vertices[i].neo4jNode.id, false);
  append!(times, toq());
end
println("Time to get a single vertex = $(mean(times) * 1000.0)ms +- $(std(times)*1000.0)ms");

# Time to update a vertex (adding labels).
times = Float32[];
for i = 1:loopCount
  vertices[i].labels = ["PerfTesting"];
  tic();
  CloudGraphs.update_vertex!(cloudGraph, vertices[i]);
  append!(times, toq());
end
println("Time to update a single vertex = $(mean(times) * 1000.0)ms +- $(std(times)*1000.0)ms");

println("Success!")

ext: jl


text:
"""
    (p::ProjectTo{T})(dx)

Projects the tangent `dx` onto a specific tangent space.

The type `T` is meant to encode the largest acceptable space, so usually
this enforces `p(dx)::T`. But some subspaces which aren't subtypes of `T` may
be allowed, and in particular `dx::AbstractZero` always passes through.

Usually `T` is the "outermost" part of the type, and `p` stores additional
properties such as projectors for each constituent field.
Arrays have either one projector `p.element` expressing the element type for
an array of numbers, or else an array of projectors `p.elements`.
These properties can be supplied as keyword arguments on construction,
`p = ProjectTo{T}(; field=data, element=Projector(x))`. For each `T` in use,
corresponding methods should be written for `ProjectTo{T}(dx)` with nonzero `dx`.

When called on `dx::Thunk`, the projection is inserted into the thunk.
"""
struct ProjectTo{P,D<:NamedTuple}
    info::D
end
ProjectTo{P}(info::D) where {P,D<:NamedTuple} = ProjectTo{P,D}(info)

# We'd like to write
# ProjectTo{P}(; kwargs...) where {P} = ProjectTo{P}(NamedTuple(kwargs))
#
# but the kwarg dispatcher has non-trivial complexity. See rules.jl for an
# explanation of this trick.
const EMPTY_NT = NamedTuple()
ProjectTo{P}() where {P} = ProjectTo{P}(EMPTY_NT)

const Type_kwfunc = Core.kwftype(Type).instance
function (::typeof(Type_kwfunc))(kws::Any, ::Type{ProjectTo{P}}) where {P}
    return ProjectTo{P}(NamedTuple(kws))
end

Base.getproperty(p::ProjectTo, name::Symbol) = getproperty(backing(p), name)
Base.propertynames(p::ProjectTo) = propertynames(backing(p))
backing(project::ProjectTo) = getfield(project, :info)

project_type(p::ProjectTo{T}) where {T} = T

function Base.show(io::IO, project::ProjectTo{T}) where {T}
    print(io, "ProjectTo{")
    show(io, T)
    print(io, "}")
    if isempty(backing(project))
        print(io, "()")
    else
        show(io, backing(project))
    end
end

# Structs
# Generic method is to recursively make `ProjectTo`s for all their fields. Not actually
# used on unknown structs, but useful for handling many known ones in the same manner.
function generic_projector(x::T; kw...) where {T}
    fields_nt::NamedTuple = backing(x)
    fields_proj = map(_maybe_projector, fields_nt)
    # We can't use `T` because if we have `Foo{Matrix{E}}` it should be allowed to make a
    # `Foo{Diagaonal{E}}` etc. Official API for this? https://github.com/JuliaLang/julia/issues/35543
    wrapT = T.name.wrapper
    return ProjectTo{wrapT}(; fields_proj..., kw...)
end

function generic_projection(project::ProjectTo{T}, dx::T) where {T}
    sub_projects = backing(project)
    sub_dxs = backing(dx)
    return construct(T, map(_maybe_call, sub_projects, sub_dxs))
end

# Used for encoding fields, leaves alone non-diff types:
_maybe_projector(x::Union{AbstractArray,Number,Ref}) = ProjectTo(x)
_maybe_projector(x) = x
# Used for re-constructing fields, restores non-diff types:
_maybe_call(f::ProjectTo, x) = f(x)
_maybe_call(f, x) = f

"""
    ProjectTo(x)

Returns a `ProjectTo{T}` functor which projects a tangent `dx` onto the
relevant tangent space for `x`.

Custom `ProjectTo` methods are provided for many subtypes of `Number` (to e.g. ensure precision),
and `AbstractArray` (to e.g. ensure sparsity structure is maintained by tangent).
Called on unknown types it will (as of v1.5.0) simply return `identity`, thus can be safely
applied to arbitrary `rrule` arguments.

# Examples
```jldoctest
julia> pr = ProjectTo(1.5f0)  # preserves real numbers, and floating point precision
ProjectTo{Float32}()

julia> pr(3 + 4im)
3.0f0

julia> pd = ProjectTo(Diagonal([1,2,3]))  # preserves structured matrices
ProjectTo{Diagonal}(diag = ProjectTo{AbstractArray}(element = ProjectTo{Float64}(), axes = (Base.OneTo(3),)),)

julia> th = @thunk reshape(1:9,3,3);

julia> pd(th) isa Thunk
true

julia> unthunk(pd(th))
3×3 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅    ⋅
  ⋅   5.0   ⋅
  ⋅    ⋅   9.0

julia> ProjectTo([1 2; 3 4]')  # no special structure, integers are promoted to float(x)
ProjectTo{AbstractArray}(element = ProjectTo{Float64}(), axes = (Base.OneTo(2), Base.OneTo(2)))
```
"""
ProjectTo(::Any) = identity

# Generic
(::ProjectTo{T})(dx::AbstractZero) where {T} = dx
(::ProjectTo{T})(dx::NotImplemented) where {T} = dx

# Thunks
(project::ProjectTo)(dx::Thunk) = Thunk(project ∘ dx.f)

# Zero
ProjectTo(::AbstractZero) = ProjectTo{NoTangent}()  # Any x::Zero in forward pass makes this one projector,
(::ProjectTo{NoTangent})(dx) = NoTangent()          # but this is the projection only for nonzero gradients,
(::ProjectTo{NoTangent})(dx::AbstractZero) = dx     # and this one solves an ambiguity.

# Also, any explicit construction with fields, where all fields project to zero, itself
# projects to zero. This simplifies projectors for wrapper types like Diagonal([true, false]).
const _PZ = ProjectTo{<:AbstractZero}
const _PZ_Tuple = Tuple{_PZ,Vararg{_PZ}}  # 1 or more ProjectTo{<:AbstractZeros}
function ProjectTo{P}(::NamedTuple{T,<:_PZ_Tuple}) where {P,T}
    return ProjectTo{NoTangent}()
end

# Tangent
# We haven't entirely figured out when to convert Tangents to "natural" representations such as
# dx::AbstractArray (when both are possible), or the reverse. So for now we just pass them through:
(::ProjectTo{T})(dx::Tangent{<:T}) where {T} = dx

#####
##### `Base`
#####

# Bool
ProjectTo(::Bool) = ProjectTo{NoTangent}()  # same projector as ProjectTo(::AbstractZero) above

# Other never-differentiable types
for T in (:Symbol, :Char, :AbstractString, :RoundingMode, :IndexStyle)
    @eval ProjectTo(::$T) = ProjectTo{NoTangent}()
end

# Numbers
ProjectTo(::Real) = ProjectTo{Real}()
ProjectTo(::Complex) = ProjectTo{Complex}()
ProjectTo(::Number) = ProjectTo{Number}()

ProjectTo(x::Integer) = ProjectTo(float(x))
ProjectTo(x::Complex{<:Integer}) = ProjectTo(float(x))

# Preserve low-precision floats as accidental promotion is a common performance bug
for T in (Float16, Float32, Float64, ComplexF16, ComplexF32, ComplexF64)
    @eval ProjectTo(::$T) = ProjectTo{$T}()
end

# In these cases we can just `convert` as we know we are dealing with plain and simple types
(::ProjectTo{T})(dx::AbstractFloat) where {T<:AbstractFloat} = convert(T, dx)
(::ProjectTo{T})(dx::Integer) where {T<:AbstractFloat} = convert(T, dx)  #needed to avoid ambiguity
# simple Complex{<:AbstractFloat}} cases
function (::ProjectTo{T})(dx::Complex{<:AbstractFloat}) where {T<:Complex{<:AbstractFloat}}
    return convert(T, dx)
end
(::ProjectTo{T})(dx::AbstractFloat) where {T<:Complex{<:AbstractFloat}} = convert(T, dx)
function (::ProjectTo{T})(dx::Complex{<:Integer}) where {T<:Complex{<:AbstractFloat}}
    return convert(T, dx)
end
(::ProjectTo{T})(dx::Integer) where {T<:Complex{<:AbstractFloat}} = convert(T, dx)

# Other numbers, including e.g. ForwardDiff.Dual and Symbolics.Sym, should pass through.
# We assume (lacking evidence to the contrary) that it is the right subspace of numebers.
(::ProjectTo{<:Number})(dx::Number) = dx

(project::ProjectTo{<:Real})(dx::Complex) = project(real(dx))
(project::ProjectTo{<:Complex})(dx::Real) = project(complex(dx))

# Tangents: we prefer to reconstruct numbers, but only safe to try when their constructor
# understands, including a mix of Zeros & reals. Other cases, we just let through:
(project::ProjectTo{<:Number})(dx::Tangent{<:Complex}) = project(Complex(dx.re, dx.im))
(::ProjectTo{<:Number})(dx::Tangent{<:Number}) = dx

# Arrays
# If we don't have a more specialized `ProjectTo` rule, we just assume that there is
# no structure worth re-imposing. Then any array is acceptable as a gradient.

# For arrays of numbers, just store one projector:
function ProjectTo(x::AbstractArray{T}) where {T<:Number}
    return ProjectTo{AbstractArray}(; element=_eltype_projectto(T), axes=axes(x))
end
ProjectTo(x::AbstractArray{Bool}) = ProjectTo{NoTangent}()

_eltype_projectto(::Type{T}) where {T<:Number} = ProjectTo(zero(T))
_eltype_projectto(::Type{<:Irrational}) = ProjectTo{Real}()

# In other cases, store a projector per element:
function ProjectTo(xs::AbstractArray)
    elements = map(ProjectTo, xs)
    if elements isa AbstractArray{<:ProjectTo{<:AbstractZero}}
        return ProjectTo{NoTangent}()  # short-circuit if all elements project to zero
    else
        # Arrays of arrays come here, and will apply projectors individually:
        return ProjectTo{AbstractArray}(; elements=elements, axes=axes(xs))
    end
end

function (project::ProjectTo{AbstractArray})(dx::AbstractArray{S,M}) where {S,M}
    # First deal with shape. The rule is that we reshape to add or remove trivial dimensions
    # like dx = ones(4,1), where x = ones(4), but throw an error on dx = ones(1,4) etc.
    dy = if axes(dx) === project.axes
        dx
    else
        for d in 1:max(M, length(project.axes))
            if size(dx, d) != length(get(project.axes, d, 1))
                throw(_projection_mismatch(project.axes, size(dx)))
            end
        end
        reshape(dx, project.axes)
    end
    # Then deal with the elements. One projector if AbstractArray{<:Number},
    # or one per element for arrays of anything else, including arrays of arrays:
    dz = if hasproperty(project, :element)
        T = project_type(project.element)
        S <: T ? dy : map(project.element, dy)
    else
        map((f, y) -> f(y), project.elements, dy)
    end
    return dz
end

# Trivial case, this won't collapse Any[NoTangent(), NoTangent()] but that's OK.
(project::ProjectTo{AbstractArray})(dx::AbstractArray{<:AbstractZero}) = NoTangent()

# Row vectors aren't acceptable as gradients for 1-row matrices:
function (project::ProjectTo{AbstractArray})(dx::LinearAlgebra.AdjOrTransAbsVec)
    return project(reshape(vec(dx), 1, :))
end

# Zero-dimensional arrays -- these have a habit of going missing,
# although really Ref() is probably a better structure.
function (project::ProjectTo{AbstractArray})(dx::Number) # ... so we restore from numbers
    if !(project.axes isa Tuple{})
        throw(
            DimensionMismatch(
                "array with ndims(x) == $(length(project.axes)) >  0 cannot have dx::Number"
            ),
        )
    end
    return fill(project.element(dx))
end

function _projection_mismatch(axes_x::Tuple, size_dx::Tuple)
    size_x = map(length, axes_x)
    return DimensionMismatch(
        "variable with size(x) == $size_x cannot have a gradient with size(dx) == $size_dx"
    )
end

#####
##### `Base`, part II: return of the Tangent
#####

# Ref
# Note that Ref is mutable. This causes Zygote to represent its structral tangent not as a NamedTuple,
# but as `Ref{Any}((x=val,))`. Here we use a Tangent, there is at present no mutable version, but see
# https://github.com/JuliaDiff/ChainRulesCore.jl/issues/105
function ProjectTo(x::Ref)
    sub = ProjectTo(x[])  # should we worry about isdefined(Ref{Vector{Int}}(), :x)?
    return ProjectTo{Tangent{typeof(x)}}(; x=sub)
end
(project::ProjectTo{<:Tangent{<:Ref}})(dx::Tangent) = project(Ref(first(backing(dx))))
function (project::ProjectTo{<:Tangent{<:Ref}})(dx::Ref)
    dy = project.x(dx[])
    return project_type(project)(; x=dy)
end
# Since this works like a zero-array in broadcasting, it should also accept a number:
(project::ProjectTo{<:Tangent{<:Ref}})(dx::Number) = project(Ref(dx))

# Tuple and NamedTuple
function ProjectTo(x::Tuple)
    elements = map(ProjectTo, x)
    if elements isa NTuple{<:Any,ProjectTo{<:AbstractZero}}
        return ProjectTo{NoTangent}()
    else
        return ProjectTo{Tangent{typeof(x)}}(; elements=elements)
    end
end
function ProjectTo(x::NamedTuple)
    elements = map(ProjectTo, x)
    if Tuple(elements) isa NTuple{<:Any,ProjectTo{<:AbstractZero}}
        return ProjectTo{NoTangent}()
    else
        return ProjectTo{Tangent{typeof(x)}}(; elements...)
    end
end

# This method means that projection is re-applied to the contents of a Tangent.
# We're not entirely sure whether this is every necessary; but it should be safe,
# and should often compile away:
function (project::ProjectTo{<:Tangent{<:Union{Tuple,NamedTuple}}})(dx::Tangent)
    return project(backing(dx))
end

function (project::ProjectTo{<:Tangent{<:Tuple}})(dx::Tuple)
    len = length(project.elements)
    if length(dx) != len
        str = "tuple with length(x) == $len cannot have a gradient with length(dx) == $(length(dx))"
        throw(DimensionMismatch(str))
    end
    # Here map will fail if the lengths don't match, but gives a much less helpful error:
    dy = map((f, x) -> f(x), project.elements, dx)
    return project_type(project)(dy...)
end
function (project::ProjectTo{<:Tangent{<:NamedTuple}})(dx::NamedTuple)
    dy = _project_namedtuple(backing(project), dx)
    return project_type(project)(; dy...)
end

# Diffractor returns not necessarily a named tuple with all keys and of the same order as
# the projector
# Thus we can't use `map`
function _project_namedtuple(f::NamedTuple{fn,ft}, x::NamedTuple{xn,xt}) where {fn,ft,xn,xt}
    if @generated
        vals = Any[
            if xn[i] in fn
                :(getfield(f, $(QuoteNode(xn[i])))(getfield(x, $(QuoteNode(xn[i])))))
            else
                throw(
                    ArgumentError(
                        "named tuple with keys(x) == $fn cannot have a gradient with key $(xn[i])",
                    ),
                )
            end for i in 1:length(xn)
        ]
        :(NamedTuple{$xn}(($(vals...),)))
    else
        vals = ntuple(Val(length(xn))) do i
            name = xn[i]
            if name in fn
                getfield(f, name)(getfield(x, name))
            else
                throw(
                    ArgumentError(
                        "named tuple with keys(x) == $fn cannot have a gradient with key $(xn[i])",
                    ),
                )
            end
        end
        NamedTuple{xn}(vals)
    end
end

function (project::ProjectTo{<:Tangent{<:Tuple}})(dx::AbstractArray)
    for d in 1:ndims(dx)
        if size(dx, d) != get(length(project.elements), d, 1)
            throw(_projection_mismatch(axes(project.elements), size(dx)))
        end
    end
    dy = reshape(dx, axes(project.elements))  # allows for dx::OffsetArray
    dz = ntuple(i -> project.elements[i](dy[i]), length(project.elements))
    return project_type(project)(dz...)
end

#####
##### `LinearAlgebra`
#####

using LinearAlgebra: AdjointAbsVec, TransposeAbsVec, AdjOrTransAbsVec

# Row vectors
ProjectTo(x::AdjointAbsVec) = ProjectTo{Adjoint}(; parent=ProjectTo(parent(x)))
# Note that while [1 2; 3 4]' isa Adjoint, we use ProjectTo{Adjoint} only to encode AdjointAbsVec.
# Transposed matrices are, like PermutedDimsArray, just a storage detail,
# but row vectors behave differently, for example [1,2,3]' * [1,2,3] isa Number
function (project::ProjectTo{Adjoint})(dx::LinearAlgebra.AdjOrTransAbsVec)
    return adjoint(project.parent(adjoint(dx)))
end
function (project::ProjectTo{Adjoint})(dx::AbstractArray)
    if size(dx, 1) != 1 || size(dx, 2) != length(project.parent.axes[1])
        throw(_projection_mismatch((1:1, project.parent.axes...), size(dx)))
    end
    dy = eltype(dx) <: Real ? vec(dx) : adjoint(dx)
    return adjoint(project.parent(dy))
end

function ProjectTo(x::LinearAlgebra.TransposeAbsVec)
    return ProjectTo{Transpose}(; parent=ProjectTo(parent(x)))
end
function (project::ProjectTo{Transpose})(dx::LinearAlgebra.AdjOrTransAbsVec)
    return transpose(project.parent(transpose(dx)))
end
function (project::ProjectTo{Transpose})(dx::AbstractArray)
    if size(dx, 1) != 1 || size(dx, 2) != length(project.parent.axes[1])
        throw(_projection_mismatch((1:1, project.parent.axes...), size(dx)))
    end
    dy = eltype(dx) <: Number ? vec(dx) : transpose(dx)
    return transpose(project.parent(dy))
end

# Diagonal
ProjectTo(x::Diagonal) = ProjectTo{Diagonal}(; diag=ProjectTo(x.diag))
(project::ProjectTo{Diagonal})(dx::AbstractMatrix) = Diagonal(project.diag(diag(dx)))
(project::ProjectTo{Diagonal})(dx::Diagonal) = Diagonal(project.diag(dx.diag))

# Symmetric
for (SymHerm, chk, fun) in
    ((:Symmetric, :issymmetric, :transpose), (:Hermitian, :ishermitian, :adjoint))
    @eval begin
        function ProjectTo(x::$SymHerm)
            sub = ProjectTo(parent(x))
            # Because the projector stores uplo, ProjectTo(Symmetric(rand(3,3) .> 0)) isn't automatically trivial:
            sub isa ProjectTo{<:AbstractZero} && return sub
            return ProjectTo{$SymHerm}(; uplo=LinearAlgebra.sym_uplo(x.uplo), parent=sub)
        end
        function (project::ProjectTo{$SymHerm})(dx::AbstractArray)
            dy = project.parent(dx)
            # Here $chk means this is efficient on same-type.
            # If we could mutate dx, then that could speed up action on dx::Matrix.
            dz = $chk(dy) ? dy : (dy .+ $fun(dy)) ./ 2
            return $SymHerm(project.parent(dz), project.uplo)
        end
        # This is an example of a subspace which is not a subtype,
        # not clear how broadly it's worthwhile to try to support this.
        function (project::ProjectTo{$SymHerm})(dx::Diagonal)
            sub = project.parent # this is going to be unhappy about the size
            sub_one = ProjectTo{project_type(sub)}(;
                element=sub.element, axes=(sub.axes[1],)
            )
            return Diagonal(sub_one(dx.diag))
        end
    end
end

# Triangular
for UL in (:UpperTriangular, :LowerTriangular, :UnitUpperTriangular, :UnitLowerTriangular) # UpperHessenberg
    @eval begin
        ProjectTo(x::$UL) = ProjectTo{$UL}(; parent=ProjectTo(parent(x)))
        (project::ProjectTo{$UL})(dx::AbstractArray) = $UL(project.parent(dx))
        function (project::ProjectTo{$UL})(dx::Diagonal)
            sub = project.parent
            sub_one = ProjectTo{project_type(sub)}(;
                element=sub.element, axes=(sub.axes[1],)
            )
            return Diagonal(sub_one(dx.diag))
        end
    end
end

# Weird -- not exhaustive!
# one strategy is to recurse into the struct:
ProjectTo(x::Bidiagonal{T}) where {T<:Number} = generic_projector(x)
function (project::ProjectTo{Bidiagonal})(dx::AbstractMatrix)
    uplo = LinearAlgebra.sym_uplo(project.uplo)
    dv = project.dv(diag(dx))
    ev = project.ev(uplo === :U ? diag(dx, 1) : diag(dx, -1))
    return Bidiagonal(dv, ev, uplo)
end
function (project::ProjectTo{Bidiagonal})(dx::Bidiagonal)
    if project.uplo == dx.uplo
        return generic_projection(project, dx) # fast path
    else
        uplo = LinearAlgebra.sym_uplo(project.uplo)
        dv = project.dv(diag(dx))
        ev = fill!(similar(dv, length(dv) - 1), 0)
        return Bidiagonal(dv, ev, uplo)
    end
end

ProjectTo(x::SymTridiagonal{T}) where {T<:Number} = generic_projector(x)
function (project::ProjectTo{SymTridiagonal})(dx::AbstractMatrix)
    dv = project.dv(diag(dx))
    ev = project.ev((diag(dx, 1) .+ diag(dx, -1)) ./ 2)
    return SymTridiagonal(dv, ev)
end
(project::ProjectTo{SymTridiagonal})(dx::SymTridiagonal) = generic_projection(project, dx)

# another strategy is just to use the AbstractArray method
function ProjectTo(x::Tridiagonal{T}) where {T<:Number}
    notparent = invoke(ProjectTo, Tuple{AbstractArray{T2}} where {T2<:Number}, x)
    return ProjectTo{Tridiagonal}(; notparent=notparent)
end
function (project::ProjectTo{Tridiagonal})(dx::AbstractArray)
    dy = project.notparent(dx)
    return Tridiagonal(dy)
end
# Note that backing(::Tridiagonal) doesn't work, https://github.com/JuliaDiff/ChainRulesCore.jl/issues/392

#####
##### `SparseArrays`
#####

using SparseArrays
# Word from on high is that we should regard all un-stored values of sparse arrays as
# structural zeros. Thus ProjectTo needs to store nzind, and get only those.
# This implementation very naiive, can probably be made more efficient.

function ProjectTo(x::SparseVector{T}) where {T<:Number}
    return ProjectTo{SparseVector}(;
        element=ProjectTo(zero(T)), nzind=x.nzind, axes=axes(x)
    )
end
function (project::ProjectTo{SparseVector})(dx::AbstractArray)
    dy = if axes(dx) == project.axes
        dx
    else
        if size(dx, 1) != length(project.axes[1])
            throw(_projection_mismatch(project.axes, size(dx)))
        end
        reshape(dx, project.axes)
    end
    nzval = map(i -> project.element(dy[i]), project.nzind)
    return SparseVector(length(dx), project.nzind, nzval)
end
function (project::ProjectTo{SparseVector})(dx::SparseVector)
    if size(dx) != map(length, project.axes)
        throw(_projection_mismatch(project.axes, size(dx)))
    end
    # When sparsity pattern is unchanged, all the time is in checking this,
    # perhaps some simple hash/checksum might be good enough?
    samepattern = project.nzind == dx.nzind
    # samepattern = length(project.nzind) == length(dx.nzind)
    if eltype(dx) <: project_type(project.element) && samepattern
        return dx
    elseif samepattern
        nzval = map(project.element, dx.nzval)
        SparseVector(length(dx), dx.nzind, nzval)
    else
        nzind = project.nzind
        # Or should we intersect? Can this exploit sorting?
        # nzind = intersect(project.nzind, dx.nzind)
        nzval = map(i -> project.element(dx[i]), nzind)
        return SparseVector(length(dx), nzind, nzval)
    end
end

function ProjectTo(x::SparseMatrixCSC{T}) where {T<:Number}
    return ProjectTo{SparseMatrixCSC}(;
        element=ProjectTo(zero(T)),
        axes=axes(x),
        rowval=rowvals(x),
        nzranges=nzrange.(Ref(x), axes(x, 2)),
        colptr=x.colptr,
    )
end
# You need not really store nzranges, you can get them from colptr -- TODO
# nzrange(S::AbstractSparseMatrixCSC, col::Integer) = getcolptr(S)[col]:(getcolptr(S)[col+1]-1)
function (project::ProjectTo{SparseMatrixCSC})(dx::AbstractArray)
    dy = if axes(dx) == project.axes
        dx
    else
        if size(dx) != (length(project.axes[1]), length(project.axes[2]))
            throw(_projection_mismatch(project.axes, size(dx)))
        end
        reshape(dx, project.axes)
    end
    nzval = Vector{project_type(project.element)}(undef, length(project.rowval))
    k = 0
    for col in project.axes[2]
        for i in project.nzranges[col]
            row = project.rowval[i]
            val = dy[row, col]
            nzval[k += 1] = project.element(val)
        end
    end
    m, n = map(length, project.axes)
    return SparseMatrixCSC(m, n, project.colptr, project.rowval, nzval)
end

function (project::ProjectTo{SparseMatrixCSC})(dx::SparseMatrixCSC)
    if size(dx) != map(length, project.axes)
        throw(_projection_mismatch(project.axes, size(dx)))
    end
    samepattern = dx.colptr == project.colptr && dx.rowval == project.rowval
    # samepattern = length(dx.colptr) == length(project.colptr) && dx.colptr[end] == project.colptr[end]
    if eltype(dx) <: project_type(project.element) && samepattern
        return dx
    elseif samepattern
        nzval = map(project.element, dx.nzval)
        m, n = size(dx)
        return SparseMatrixCSC(m, n, dx.colptr, dx.rowval, nzval)
    else
        invoke(project, Tuple{AbstractArray}, dx)
    end
end


text:
# mechanism for registering physics modules, initial conditions, and boundary
# conditions

"""
  This dictionary maps physics module names to the module and functions needed
  to start a simulation.  Specifically, the values in the dictionary are a
  tuple of the form:

  `(PhysicsModule, _createObjects, _checkOptions)`

  See [`register_physics`](@ref) for the meaning of each element of the tuple.

  The fact that this dictionary exists is an
  implementation details, it should never be accessed directly.  Instead,
  the accessor functions below should be used.
"""
global const PhysicsModDict = Dict{String, Tuple{Module, Function, Function}}()

"""
  This function registered a new physics module with the global list of all
  known physics modules.  Every physics module should do this as part of
  module initialization.  The name, the module, and the startup function must
  be unique (ie. they must not already exist in the list).  This function
  throws and exception if they are not.

  **Inputs**

   * modname:  an String name for this entry in the list.  It is used
                to retrieve the module and startup function in the
                retrieve_physics function. Typically the name is capitalized.
   * mod:  the Module itself

   * _createObjects: function that creates the mesh, sbp, eqn, and opts
                      objects. Must have a methods with signatures
                      `_createObjects(opt::Dict)`  and

                      If this physics modules supports solving on a submesh,
                      this function should also have a method
                      `_createObjects(mesh::AbstractMesh, sbp::AbstractOperator, opts::Dict)`
                      to create a new equation object.  See [`createObjects`](@ref).

   * _checkOptions: physics-specific function for supplying default options and
                   checking options.  Note that most of the input option
                   processing is done in [`read_input`](@ref Input.read_input),
                   `_checkOptions` need only do the physics-specific part.
                   This function must have signature `_checkOptions(opts::Dict)`.

  **Outputs**

    none
"""
function register_physics(modname::String, mod::Module, _createObjects::Function, _checkOptions::Function)

  # check if this name is already used
  if haskey(PhysicsModDict, modname)
    mod_i, func_i = PhysicsModDict[modname]
    throw(ErrorException("Physics module with name $modname is already registered to module $mod_i"))
  end

  # check if the module or function is already registered
  for (key, val) in PhysicsModDict
    mod_i = val[1]
    _createObjects_i = val[2]
    _checkOptions_i = val[3]

    if mod ==  mod_i
      throw(ErrorException("Physics module $mod is already registered to name $key"))
    end

    if _createObjects == _createObjects_i
      throw(ErrorException("_createObjects function $_createObjects is already registered to physics module $mod_i with name $key"))
    end

    if _checkOptions == _checkOptions_i
      throw(ErrorException("_checkOptions function $_checkOptions is already registered to physics module $mod_i with name $key"))
    end

  end

  if _createObjects == PDESolver.createObjects
    error("PDESolver.createObjects cannot be registered as the physics module _createObjects function, this would create an infinite loop")
  end

  # if we get here, the registration is new
  PhysicsModDict[modname] = (mod, _createObjects, _checkOptions)

  # register the options checker
  registerOptionsChecker(modname, _checkOptions)

  return nothing
end  # end function

"""
  Retrieves the physics module and function registered using [`register_physics`](@ref)

  **Input**

   * modname: an String containing the name of the module supplied to
             `register_physics`

  **Outputs**

   * mod: the physics Module
   * _createObjects: function that creates the solver objects
   * _checkOptions: the options checking function
"""
function retrieve_physics(modname::String)

  if !haskey(PhysicsModDict, modname)
    # construct the error message
    err_msg = "unknown physics module name $modname, known names are:"
    for known_modname in keys(PhysicsModDict)
      err_msg *= "\n"*known_modname
    end

    throw(ErrorException(err_msg))
  end

  # if we get here, then this is a real physics module name

  return PhysicsModDict[modname]
end

"""
  This function registers a new initial condition function with the specified
  physics module.  The function must have the signature:

  ```
   ICfunc(mesh::AbstractMesh, sbp::AbstractOperator, eqn:AbstractSolutionData{Tsol},
                opts:Dict, q_vec::AbstractVector{Tsol})
  ```

  where q_vec is a vector of length numDof that will be populated with the
  initial condition.  If the function is used to provide the exact solution
  for an unsteady problem (for calculating the error via the `exact_soln_func`
  key in the options dictionary), then it should use eqn.params.t as the
  current time value.

  This function does not attempt to verify that the functor has the correct
  signature, so the user should take care to ensure it is correct for the
  physics module.

  The physics module must have an associative container called `ICDict` that
  accepts Strings as keys and functions as values.

  It is an error to register two functions under the same name or to register
  one function under two different names.

  Inputs:

    mod: physics module to register the function with
    fname: name associated with the function, used as the value for the
           for any key in the options dictionary that specifies an initial
           condition function
    func: the function itself

  Outputs:

    none
"""
function registerIC(mod::Module, fname::String, func::Function)

  # check if name is already registered
  if haskey(mod.ICDict, fname)
    throw(ErrorException("IC name $fname is already registered with physics module $mod"))
  end

  # check if function is already registered
  for (key, val) in mod.ICDict
    if val == func
      throw(ErrorException("IC function $func is already registered with physics module $mod"))
    end
  end

  # there is currently no good way of verify the signature. func.env.def.sig
  # is a tuple of argument types, but julia's type system is covarient type
  # system makes checking them difficult

  # if we get here, this name/function pair is unique
  mod.ICDict[fname] = func

  return nothing
end

"""
  This function registers a new boundary condition functor with the specified
  physics module.  The exact signature of the functor varies by physics module.
  The purpose of the functor is to calculate the flux across a boundary at
  a particular node. It usually takes in the solution values at that node, as
  well as the coordinates and mapping jacobian.  The params type is typically
  used to dispatch to different methods for 2D or 3D, if needed.  Generally,
  this functor will calculate the boundary state and then call a numerical
  flux function of some kind to compute the flux.

  This function does not attempt to verify that the functor has the correct
  signature, so the user should take care to ensure it is correct for the
  physics module.

  It is an error to register two functors with the same name or the same
  functor under two different names.  An exception is made for "reanalysisBC",
  because it is special.

  **Inputs**

   * mod: module to register the the functor with
   * fname: the name associated with this function, used as the value for any
            key in the options dictionary that specifies a boundary condition,
            for example `BC1_name`
   * func: the functor itself, or the constructor for one

  **Outputs**

    none

  **Notes For Implementers**

  This function works by utilizing `mod.BCDict`, which must be an associative
  collection mapping BC names to functors.
"""
function registerBC(mod::Module, fname::String, func::Union{BCType, Type{T} where T <: BCType})

  # check if name is already registered
  # special case for reanalysisBC because it is special
  if haskey(mod.BCDict, fname) && fname != "reanalysisBC"
    throw(ErrorException("BC name $fname is already registered with physics module $mod"))
  end

  # check if function is already registered
  for (key, val) in mod.BCDict
    if val == func && key != "reanalysisBC"
      throw(ErrorException("BC functor $func is already registered with physics module $mod"))
    end
  end

  mod.BCDict[fname] = func

  return nothing
end

#TODO: add method to register source terms


text:
@testset "diagonal" begin
    @testset "Diagonal" begin
        rng, N = MersenneTwister(123456), 3
        rrule_test(Diagonal, randn(rng, N, N), (randn(rng, N), randn(rng, N)))
        D = Diagonal(randn(rng, N))
        rrule_test(Diagonal, D, (randn(rng, N), randn(rng, N)))
        # Concrete type instead of UnionAll
        rrule_test(typeof(D), D, (randn(rng, N), randn(rng, N)))
    end
    @testset "diag" begin
        rng, N = MersenneTwister(123456), 7
        rrule_test(diag, randn(rng, N), (randn(rng, N, N), randn(rng, N, N)))
        rrule_test(diag, randn(rng, N), (Diagonal(randn(rng, N)), randn(rng, N, N)))
        rrule_test(diag, randn(rng, N), (randn(rng, N, N), Diagonal(randn(rng, N))))
        rrule_test(diag, randn(rng, N), (Diagonal(randn(rng, N)), Diagonal(randn(rng, N))))
    end
end


text:
# ---
# title: 1207. Unique Number of Occurrences
# id: problem1207
# author: Jerry Ling
# date: 2021-11-30
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/unique-number-of-occurrences/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`, write a function that returns `true` if and
# only if the number of occurrences of each value in the array is unique.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,2,2,1,1,3]
#     Output: true
#     Explanation:  The value 1 has 3 occurrences, 2 has 2 and 3 has 1. No two values have the same number of occurrences.
#
# **Example 2:**
#
#
#
#     Input: arr = [1,2]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: arr = [-3,0,1,-3,1,1,1,-3,10,0]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 1000`
#   * `-1000 <= arr[i] <= 1000`
#
#
## @lc code=start
using LeetCode

function unique_number_of_occurrences(arr::Vector{Int})
    counts = Dict{Int, Int}()
    for a in arr
        counts[a] = get(counts, a, 0) + 1
    end
    vs = values(counts)
    return length(vs) == length(unique(vs))
end

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


text:
export FieldmapNFFTOp, createInhomogeneityData_

mutable struct FieldmapNFFTOp{T,F1,F2} <:
                      AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod :: Function
  tprod :: F1
  ctprod :: F2
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
end

mutable struct InhomogeneityData
  A_k::Matrix{ComplexF64}
  C_k::Matrix{ComplexF64}
  times::Vector{Float64}
  Cmap::Matrix{ComplexF64}
  t_hat::Float64
  z_hat::ComplexF64
  method::String
end

"""
    FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory,
                        correctionmap::Array{ComplexF64,D};
                        method::String="nfft",
                        echoImage::Bool=true,
                        alpha::Float64=1.75,
                        m::Float64=3.0,
                        K=20,
                        kargs...) where D

generates a `FieldmapNFFTOp` which evaluates the MRI Fourier signal encoding operator,
including B0-inhomogeneities using time-segmented NFFTs.

# Arguments:
* `shape::NTuple{D,Int64}`             - size of image to encode/reconstruct
* `tr::Trajectory`                     - Trajectory with the kspace nodes to sample
* `correctionmap::Array{ComplexF64,D}` - fieldmap for the correction of off-resonance effects
* (`method::String="nfft"`)            - method to use for time-segmentation when correctio field inhomogeneities
* (`echoImage::Bool=false`)            - if true sampling times will only be considered relative to the echo time
                                         this results in complex valued image even for real-valued input.
* (`alpha::Float64=1.75`)              - oversampling factor for interpolation
* (`m::Float64=3.0`)                   - truncation size of interpolation kernel
* (`K=20`)                             - number of translates for LeastSquares approaches
                                         (not NFFT-approach) to time-segmentation
* (`kargs`)                            - additional keyword arguments
"""
function FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory,
                        correctionmap::Array{ComplexF64,D};
                        method::String="nfft",
                        echoImage::Bool=true,
                        alpha::Float64=1.25,
                        m::Float64=2.0,
                        K=20,
                        kargs...) where D

  nodes,times = kspaceNodes(tr), readoutTimes(tr)
  if echoImage
    times = times .- echoTime(tr)
  end
  nrow = size(nodes,2)
  ncol = prod(shape)

 # create and truncate low-rank expansion
  cparam = createInhomogeneityData_(vec(times), correctionmap; K=K, alpha=alpha, m=m, method=method)
  K = size(cparam.A_k,2)

  @debug "K = $K"

  plan = Vector{NFFTPlan{D,0,Float64}}(undef,K)
  idx = Vector{Vector{Int64}}(undef,K)
  for κ=1:K
    idx[κ] = findall(x->x!=0.0, cparam.A_k[:,κ])
    plan[κ] = NFFTPlan(nodes[:,idx[κ]], shape, 3, 1.25, precompute = NFFT.FULL)
  end

  p = [zeros(ComplexF64, ncol) for t=1:Threads.nthreads() ]
  y = [zeros(ComplexF64, nrow) for t=1:Threads.nthreads() ]
  d = [zeros(ComplexF64, length(idx[κ])) for κ=1:K ]

  mul(x::Vector{T}) where T<:ComplexF64 =
     produ(x,nrow,ncol,shape,plan,idx,cparam,isCircular(tr),p,y,d)
  ctmul(y::Vector{T}) where T<:ComplexF64 =
     ctprodu(y,shape,plan,idx,cparam,isCircular(tr),p,y,d)
  inverse(y::Vector{T}) where T<:ComplexF64 =
     inv(y,shape,plan,idx,cparam,isCircular(tr),p,y,d)

  return FieldmapNFFTOp{ComplexF64,Nothing,Function}(nrow, ncol, false, false
            , mul
            , nothing
            , ctmul, 0, 0, 0)
end

# function produ{T<:ComplexF64}(x::Vector{T}, numOfNodes::Int, numOfPixel::Int, shape::Tuple, plan::Vector{NFFTPlan{2,0,ComplexF64}}, cparam::InhomogeneityData, density::Vector{Float64}, symmetrize::Bool)
function produ(x::Vector{T}, numOfNodes::Int, numOfPixel::Int, shape::Tuple, plan,
               idx::Vector{Vector{Int64}}, cparam::InhomogeneityData,
               shutter::Bool, p, y, d) where T<:ComplexF64
  K = size(cparam.A_k,2)
  s = zeros(ComplexF64,numOfNodes)

  if shutter
    circularShutter!(reshape(x, shape), 1.0)
  end

  # Preprocessing step when time and correctionMap are centered
  if cparam.method == "nfft"
    x_ = x .* exp.(-vec(cparam.Cmap) * cparam.t_hat )
  else
    x_ = copy(x)
  end

  sp = Threads.SpinLock()
  produ_inner(K,cparam.C_k, cparam.A_k, shape, p, d, y, s, sp, plan, idx, x_)

  # Postprocessing step when time and correctionMap are centered
  if cparam.method == "nfft"
      s .*= exp.(-cparam.z_hat*(cparam.times .- cparam.t_hat) )
  end
  return s
end

function produ_inner(K, C, A, shape, p, d, y, s, sp, plan, idx, x_)
  Threads.@threads for κ=1:K
    t = Threads.threadid()
    for l=1:length(x_)
      p[t][l] = C[κ,l] * x_[l]
    end
    NFFT.nfft!(plan[κ], reshape(p[t], shape), d[κ])
    fill!(y[t], 0.0)
    for k=1:length(idx[κ])
      y[t][idx[κ][k]] = d[κ][k]
    end
    for l=1:size(A,1)
      y[t][l] *= A[l,κ]
    end
    lock(sp)
    s .+= y[t]
    unlock(sp)
  end
  return
end

# function ctprodu{T<:ComplexF64}(x::Vector{T}, shape::Tuple, plan::Vector{NFFTPlan{2,0,ComplexF64}}, cparam::InhomogeneityData, density::Vector{Float64}, symmetrize::Bool)
function ctprodu(x::Vector{T}, shape::Tuple, plan, idx::Vector{Vector{Int64}},
                 cparam::InhomogeneityData, shutter::Bool, p, y_, d) where T<:ComplexF64

  y = zeros(ComplexF64,prod(shape))
  K = size(cparam.A_k,2)

  x_ = copy(x)

  # Preprocessing step when time and correctionMap are centered
  if cparam.method == "nfft"
      x_ .*= conj.(exp.(-cparam.z_hat*(cparam.times .- cparam.t_hat)))
  end

  sp = Threads.SpinLock()
  ctprodu_inner(K,cparam.C_k, cparam.A_k, shape, p, d, y, sp, plan, idx, x_)

  if cparam.method == "nfft"
    y .*=  conj(exp.(-vec(cparam.Cmap) * cparam.t_hat))
  end

  if shutter
    circularShutter!(reshape(y, shape), 1.0)
  end

  return y
end

function ctprodu_inner(K, C, A, shape, p, d, y, sp, plan, idx, x_)
  Threads.@threads for κ=1:K
    t = Threads.threadid()
    for k=1:length(idx[κ])
      d[κ][k] = conj.(A[idx[κ][k],κ]) * x_[idx[κ][k]]
    end
    NFFT.nfft_adjoint!(plan[κ], d[κ], reshape(p[t], shape))
    for k=1:length(p[t])
      p[t][k] = conj(C[κ,k]) * p[t][k]
    end
    lock(sp)
    y .+= p[t]
    unlock(sp)
  end
  return
end

function createInhomogeneityData_(times::Vector,
                                  correctionmap::Array{ComplexF64,D};
                                  K::Int64=20,
                                  alpha::Float64=1.75,
                                  m::Float64 = 4.0,
                                  method="nfft") where D

    if method == "const"
      A = get_AC_coefficients_one_term(K,vec(times),vec(correctionmap))
      C = get_C_coefficients_hist_lsqr(K,vec(times),vec(correctionmap))
    elseif method== "linear"
      A = get_AC_coefficients_two_terms(K,vec(times),vec(correctionmap))
      C = get_C_coefficients_hist_lsqr(K,vec(times),vec(correctionmap))
    elseif method == "nfft"
      A,C = get_AC_coefficients_nfft(K,vec(times),vec(correctionmap), alpha, m)
    elseif method == "leastsquare"
      A,C = get_AC_coefficients_lsqr(K,vec(times),vec(correctionmap))
    elseif method == "hist"
      A = get_A_coefficients_hist_lsqr(K,vec(times),vec(correctionmap))
      C = get_C_coefficients_hist_lsqr(K,vec(times),vec(correctionmap))
    else
      error("approximation scheme $(interp) is not yet implemented")
    end
    if size(A,2) != size(C,1)
      error("Consistency check failed! A and C are not compatible")
    end

    t_hat = (times[1] + times[end])/2
    z_hat = minimum(real(correctionmap)) + maximum(real(correctionmap))
    z_hat += 1im*(minimum(imag(correctionmap)) + maximum(imag(correctionmap)))
    z_hat *= 0.5

    return InhomogeneityData(A,C,vec(times),correctionmap,t_hat,z_hat,method)
end

function adjoint(op::FieldmapNFFTOp{T}) where T
  return LinearOperator{T}(op.ncol, op.nrow, op.symmetric, op.hermitian,
                        op.ctprod, nothing, op.prod)
end


text:
module Firefly

export findpeak,
       findpeaks

using KernelDensity: kde
using MCMCChains: Chains
using StatsBase: rle

include("distributions.jl")

###############################################################################

"""
    findpeak(::AbstractVector)
    findpeak(::Chains, [::Symbol])

Return a naïve maximum a posteriori (MAP) estimate from a chain.

This uses kernel density estimation (KDE) to estimate the continuous posterior distribution from the input chain and simply returns its max argument. This is naïve because we should not consider KDE peaks to be true "modes" or MAP estimates. Estimating accurate modes of a posterior sample is still an active area of statistical research.

# Examples
```jldoctest
julia> samples = randn(2000);

julia> findpeak(samples) # should be close to 0
-0.24696046398697033
```
"""
function findpeak(samples::AbstractVector)
    k = kde(samples)
    return k.x[argmax(k.density)]
end

findpeak(chain::Chains) = findpeak(Array(chain))
findpeak(chain::Chains, s::Symbol) = findpeak(chain[s])

"""
    findpeaks(::AbstractVector, [n])
    findpeaks(::Chains, [::Symbol], [n])

Return the first `n` naïve maximum a posteriori (MAP) estimates from a chain.

This uses kernel density estimation (KDE) to estimate the continuous posterior distribution from the input sample and finds maxima via the finite-difference estimation of the derivative. This is naïve because we should not consider KDE peaks to be true "modes" or MAP estimates. Estimating accurate modes of a posterior sample is still an active area of statistical research.

# Examples
```jldoctest
julia> samples = append!(randn(1000), randn(800) .+ 10);

julia> findpeaks(samples, 2) # should be close to [0, 10]
2-element Array{Float64,1}:
 -0.059655173778625246
  9.987428704498656
```
"""
function findpeaks(samples::AbstractVector)
    k = kde(samples)
    # the sign of the difference will tell use whether we are increasing or decreasing
    # using rle gives us the points at which the sign switches (local extreema)
    runs = rle(sign.(diff(k.density)))
    # if we start going up, first extreme will be maximum, else minimum
    start = runs[1][1] == 1 ? 1 : 2
    # find the peak indices at the local minima
    peak_idx = cumsum(runs[2])[start:2:end]
    sorted_idx = sortperm(k.density[peak_idx], rev = true)
    return k.x[peak_idx[sorted_idx]]
end

findpeaks(samples::AbstractVector, n::Integer) = findpeaks(samples)[1:n]
findpeaks(chain::Chains) = findpeaks(Array(chain))
findpeaks(chain::Chains, n::Integer) = findpeaks(Array(chain), n)
findpeaks(chain::Chains, s::Symbol) = findpeaks(chain[s])
findpeaks(chain::Chains, s::Symbol, n::Integer) = findpeaks(chain[s], n)

end


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test

const OPTIMIZER = SCIP.Optimizer(display_verblevel=0)
const CONFIG_DIRECT = MOIT.Config(
    atol=1e-5, rtol=1e-5,
    exclude=Any[
        MOI.ConstraintDual, MOI.ConstraintName, MOI.DualObjectiveValue, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus,
    ],
)

@testset "MOI unit tests" begin
    excluded = copy(MOI_BASE_EXCLUDED)
    append!(excluded, [
        "test_linear_integration", # Can not delete variable while model contains constraints
        "test_basic_VectorOfVariables_SecondOrderCone",
        "test_conic_SecondOrderCone_nonnegative_post_bound",
        "test_variable_delete_SecondOrderCone",
        "test_modification_func_scalaraffine_",
        "test_modification_func_vectoraffine_",
    ])
    MOIT.runtests(
        OPTIMIZER,
        CONFIG_DIRECT,
        warn_unsupported=false,
        exclude = excluded,
    )

end

function indicator4_test(model::MOI.ModelLike, config)
    atol = config.atol
    rtol = config.rtol
    # equivalent to indicator1_test with left-hand-side partially in LHS constant
    # linear problem with indicator constraint and
    # max  2x1 + 3x2
    # s.t. x1 + x2 <= 10
    #      z1 ==> x2 - 1 <= 7
    #      z2 ==> x2 + x1/5 + 1 <= 10
    #      z1 + z2 >= 1

    MOI.empty!(model)
    @test MOI.is_empty(model)

    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{Float64}})
    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1 = MOI.add_variable(model)
    z2 = MOI.add_variable(model)
    MOI.add_constraint(model, z1, MOI.ZeroOne())
    MOI.add_constraint(model, z2, MOI.ZeroOne())
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, -1.0]
    )
    iset1 = MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(7.0))
    MOI.add_constraint(model, f1, iset1)

    f2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(0.2, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 1.0],
    )
    iset2 = MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(10.0))

    MOI.add_constraint(model, f2, iset2)

    # Additional regular constraint.
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, x2)], 0.0),
        MOI.LessThan(10.0),
    )

    # Disjunction z1 ⋁ z2
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, z1), MOI.ScalarAffineTerm(1.0, z2)], 0.0),
        MOI.GreaterThan(1.0),
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0], [x1, x2]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 28.75 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), x1) ≈ 1.25 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), x2) ≈ 8.75 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), z1) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), z2) ≈ 1.0 atol=atol rtol=rtol
end

@testset "MOI Integer Linear" begin
    # replace with MOIT when MathOptInterface.jl#929 merged
    indicator4_test(OPTIMIZER, CONFIG_DIRECT)
end


text:
function solve_cvrp(cvrp::CVRP, parameters=AlgorithmParameters(); maximum_number_of_vehicles=typemax(Cint), verbose=true, use_dist_mtx=false)
    n = cvrp.dimension
    x = cvrp.coordinates[:, 1]
    y = cvrp.coordinates[:, 2]
    serv_time = zeros(size(x))
    dem = cvrp.demand
    vehicleCapacity = cvrp.capacity

    if use_dist_mtx
        dist_mtx = zeros(n, n)
        for i in 1:n, j in 1:n
            dist_mtx[i, j] = cvrp.weights[i, j]
        end
        # need to input dist_mtx' instead of dist_mtx
        # Julia: column-first indexing
        # C: row-first indexing
        return c_api_solve_cvrp_dist_mtx(n, x, y, Matrix(dist_mtx'), serv_time, dem, vehicleCapacity, maximum_number_of_vehicles, parameters, verbose)
    else
        return c_api_solve_cvrp(n, x, y, serv_time, dem, vehicleCapacity, maximum_number_of_vehicles, parameters, verbose)
    end
end

function solve_cvrp(cvrp_file_path::AbstractString, parameters=AlgorithmParameters(); maximum_number_of_vehicles=typemax(Cint), verbose=true)
    cvrp = CVRPLIB.readCVRP(cvrp_file_path)
    return solve_cvrp(cvrp, parameters; maximum_number_of_vehicles=maximum_number_of_vehicles, verbose=verbose, use_dist_mtx=true)
end

function solve_cvrp(x::Vector, y::Vector, service_time::Vector, demand::Vector, vehicle_capacity::Integer, n_vehicles::Integer, parameters=AlgorithmParameters(); verbose=true)
    @assert length(x) == length(y) == length(service_time) == length(demand)
    return c_api_solve_cvrp(length(x), x, y, service_time, demand, vehicle_capacity, n_vehicles, parameters, verbose)
end

function solve_cvrp(
    dist_mtx::Matrix, service_time::Vector, demand::Vector, vehicle_capacity::Integer, n_vehicles::Integer, parameters=AlgorithmParameters();
    verbose=true, x_coords=zeros(length(demand)), y_coords=zeros(length(demand))
)
    @assert length(x_coords) == length(y_coords) == length(service_time) == length(demand)
    return c_api_solve_cvrp_dist_mtx(length(x_coords), x_coords, y_coords, Matrix(dist_mtx'), service_time, demand, vehicle_capacity, n_vehicles, parameters, verbose)
end


text:
function test_iders()
    model = MetNets.toy_model()

    too_big = 100
    negative = -100
    @test_throws ErrorException MetNets.metindex(model, too_big)
    @test_throws ErrorException MetNets.metindex(model, negative)

    @test_throws ErrorException MetNets.rxnindex(model, too_big)
    @test_throws ErrorException MetNets.rxnindex(model, negative)

    for (i, rxn) in enumerate(MetNets.rxns(model))
        @test MetNets.rxnindex(model, i) == i
        @test MetNets.rxnindex(model, rxn) == i
    end

    for (i, met) in enumerate(MetNets.mets(model))
        @test MetNets.metindex(model, i) == i
        @test MetNets.metindex(model, met) == i
    end

    @test_throws ErrorException MetNets.metindex(model, "met_not_in_model")
    @test_throws ErrorException MetNets.rxnindex(model, "rxn_not_in_model")

end
test_iders()


text:
"""
$(TYPEDEF)

An abstract type the subtypes of which are associated with functions method
overloaded for for new extensions. An instance of the `DefaultExt <:ExtensionType`
structure to the `Optimizer` in the `ext_type` field.
"""
abstract type ExtensionType end
struct DefaultExt <: ExtensionType end
MOIU.map_indices(::Function, x::ExtensionType) = x
MOIU.map_indices(::Function, x::DefaultExt) = x

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions