text:
@testset "Float" begin
    Random.seed!(0)
    ket = rand(4)
    bra = ket'
    ρ = rand(4,4)

    DiracNotation.reset_properties()
    @test sprint(dirac, ket) == "|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩+0.177329|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket)) == "|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩+0.177329|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ket)) == "|ψ⟩ = 0.8236475079774124|00⟩+0.9103565379264364|01⟩+0.16456579813368521|10⟩+0.17732884646626457|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2])) == "|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩+0.177329|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2], header=true)) == "4-element Array{Float64,1}\n|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩+0.177329|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2], "ϕ")) == "|ϕ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩+0.177329|11⟩\n"

    DiracNotation.set_properties(newline=true)
    @test sprint(dirac, ket) == "|ψ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩\n     +0.177329|11⟩\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket)) == "|ψ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩\n     +0.177329|11⟩\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ket)) == "|ψ⟩ = 0.8236475079774124|00⟩\n     +0.9103565379264364|01⟩\n     +0.16456579813368521|10⟩\n     +0.17732884646626457|11⟩\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2])) == "|ψ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩\n     +0.177329|11⟩\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2], "ϕ")) == "|ϕ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩\n     +0.177329|11⟩\n     "

    DiracNotation.reset_properties()
    DiracNotation.set_properties(precision=2)
    @test sprint(dirac, ket) == "|ψ⟩ = 0.82|00⟩+0.91|01⟩+0.16|10⟩+0.18|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket)) == "|ψ⟩ = 0.82|00⟩+0.91|01⟩+0.16|10⟩+0.18|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ket)) == "|ψ⟩ = 0.82|00⟩+0.91|01⟩+0.16|10⟩+0.18|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2])) == "|ψ⟩ = 0.82|00⟩+0.91|01⟩+0.16|10⟩+0.18|11⟩\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2], "ϕ")) == "|ϕ⟩ = 0.82|00⟩+0.91|01⟩+0.16|10⟩+0.18|11⟩\n"

    DiracNotation.reset_properties()
    DiracNotation.set_properties(displayall=false, numhead=3)
    sprint(dirac, ket) == "|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩ +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>true), ket)) == "|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩ +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>false), ket)) == "|ψ⟩ = 0.8236475079774124|00⟩+0.9103565379264364|01⟩+0.16456579813368521|10⟩ +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2])) == "|ψ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩ +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2], "ϕ")) == "|ϕ⟩ = 0.823648|00⟩+0.910357|01⟩+0.164566|10⟩ +...\n"

    DiracNotation.reset_properties()
    DiracNotation.set_properties(displayall=false, numhead=3, newline=true)
    @test sprint(dirac, ket) == "|ψ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩ +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket)) == "|ψ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩ +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ket)) == "|ψ⟩ = 0.8236475079774124|00⟩\n     +0.9103565379264364|01⟩\n     +0.16456579813368521|10⟩ +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2])) == "|ψ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩ +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ket, [2,2], "ϕ")) == "|ϕ⟩ = 0.823648|00⟩\n     +0.910357|01⟩\n     +0.164566|10⟩ +...\n"

    DiracNotation.reset_properties()
    @test sprint(dirac, bra) == "⟨ψ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10|+0.177329⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra)) == "⟨ψ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10|+0.177329⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), bra)) == "⟨ψ| = 0.8236475079774124⟨00|+0.9103565379264364⟨01|+0.16456579813368521⟨10|+0.17732884646626457⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2])) == "⟨ψ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10|+0.177329⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2], "ϕ")) == "⟨ϕ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10|+0.177329⟨11|\n"

    DiracNotation.set_properties(newline=true)
    @test sprint(dirac, bra) == "⟨ψ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10|\n     +0.177329⟨11|\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra)) == "⟨ψ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10|\n     +0.177329⟨11|\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>false), bra)) == "⟨ψ| = 0.8236475079774124⟨00|\n     +0.9103565379264364⟨01|\n     +0.16456579813368521⟨10|\n     +0.17732884646626457⟨11|\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2])) == "⟨ψ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10|\n     +0.177329⟨11|\n     "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2], "ϕ")) == "⟨ϕ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10|\n     +0.177329⟨11|\n     "

    DiracNotation.reset_properties()
    DiracNotation.set_properties(precision=2)
    @test sprint(dirac, bra) == "⟨ψ| = 0.82⟨00|+0.91⟨01|+0.16⟨10|+0.18⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra)) == "⟨ψ| = 0.82⟨00|+0.91⟨01|+0.16⟨10|+0.18⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), bra)) == "⟨ψ| = 0.82⟨00|+0.91⟨01|+0.16⟨10|+0.18⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2])) == "⟨ψ| = 0.82⟨00|+0.91⟨01|+0.16⟨10|+0.18⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2], "ϕ")) == "⟨ϕ| = 0.82⟨00|+0.91⟨01|+0.16⟨10|+0.18⟨11|\n"

    DiracNotation.reset_properties()
    DiracNotation.set_properties(displayall=false, numhead=3)
    sprint(dirac, bra) == "⟨ψ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10| +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>true), bra)) == "⟨ψ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10| +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>false), bra)) == "⟨ψ| = 0.8236475079774124⟨00|+0.9103565379264364⟨01|+0.16456579813368521⟨10| +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2])) == "⟨ψ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10| +...\n"
    sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2], "ϕ")) == "⟨ϕ| = 0.823648⟨00|+0.910357⟨01|+0.164566⟨10| +...\n"

    DiracNotation.reset_properties()
    DiracNotation.set_properties(displayall=false, numhead=3, newline=true)
    @test sprint(dirac, bra) == "⟨ψ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra)) == "⟨ψ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), bra)) == "⟨ψ| = 0.8236475079774124⟨00|\n     +0.9103565379264364⟨01|\n     +0.16456579813368521⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2])) == "⟨ψ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), bra, [2,2], "ϕ")) == "⟨ϕ| = 0.823648⟨00|\n     +0.910357⟨01|\n     +0.164566⟨10| +...\n"

    DiracNotation.reset_properties()
    @test sprint(dirac, ρ) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|+0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|+0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|+0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ)) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|+0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|+0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|+0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ρ)) == "ρ = 0.278880109331201|00⟩⟨00|+0.3618283907762174|00⟩⟨01|+0.26003585026904785|00⟩⟨10|+0.5758873948500367|00⟩⟨11|+0.20347655804192266|01⟩⟨00|+0.9732164043865108|01⟩⟨01|+0.910046541351011|01⟩⟨10|+0.8682787096942046|01⟩⟨11|+0.042301665932029664|10⟩⟨00|+0.5858115517433242|10⟩⟨01|+0.16703619444214968|10⟩⟨10|+0.9677995536192001|10⟩⟨11|+0.06826925550564478|11⟩⟨00|+0.5392892841426182|11⟩⟨01|+0.6554484126999125|11⟩⟨10|+0.7676903325581188|11⟩⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2])) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|+0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|+0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|+0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2], header=true)) == "4×4 Array{Float64,2}\nρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|+0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|+0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|+0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2], "ϕ")) == "ϕ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|+0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|+0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|+0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n"

    DiracNotation.reset_properties()
    DiracNotation.set_properties(newline=true)
    @test sprint(dirac, ρ) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|\n   +0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|\n   +0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|\n   +0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n   "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ)) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|\n   +0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|\n   +0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|\n   +0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n   "
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ρ)) == "ρ = 0.278880109331201|00⟩⟨00|+0.3618283907762174|00⟩⟨01|+0.26003585026904785|00⟩⟨10|+0.5758873948500367|00⟩⟨11|\n   +0.20347655804192266|01⟩⟨00|+0.9732164043865108|01⟩⟨01|+0.910046541351011|01⟩⟨10|+0.8682787096942046|01⟩⟨11|\n   +0.042301665932029664|10⟩⟨00|+0.5858115517433242|10⟩⟨01|+0.16703619444214968|10⟩⟨10|+0.9677995536192001|10⟩⟨11|\n   +0.06826925550564478|11⟩⟨00|+0.5392892841426182|11⟩⟨01|+0.6554484126999125|11⟩⟨10|+0.7676903325581188|11⟩⟨11|\n   "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2])) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|\n   +0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|\n   +0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|\n   +0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n   "
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2], "ϕ")) == "ϕ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10|+0.575887|00⟩⟨11|\n   +0.203477|01⟩⟨00|+0.973216|01⟩⟨01|+0.910047|01⟩⟨10|+0.868279|01⟩⟨11|\n   +0.0423017|10⟩⟨00|+0.585812|10⟩⟨01|+0.167036|10⟩⟨10|+0.9678|10⟩⟨11|\n   +0.0682693|11⟩⟨00|+0.539289|11⟩⟨01|+0.655448|11⟩⟨10|+0.76769|11⟩⟨11|\n   "

    DiracNotation.reset_properties()
    DiracNotation.set_properties(displayall=false, numhead=3)
    @test sprint(dirac, ρ) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ)) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>false), ρ)) == "ρ = 0.278880109331201|00⟩⟨00|+0.3618283907762174|00⟩⟨01|+0.26003585026904785|00⟩⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2])) == "ρ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10| +...\n"
    @test sprint(io -> dirac(IOContext(io, :compact=>true), ρ, [2,2], [2,2], "ϕ")) == "ϕ = 0.27888|00⟩⟨00|+0.361828|00⟩⟨01|+0.260036|00⟩⟨10| +...\n"

end # testset "Float64"

DiracNotation.reset_properties()
@testset "Int" begin
    @test sprint(dirac, [1,0]) == "|ψ⟩ = |0⟩\n"
    @test sprint(dirac, [0,1]) == "|ψ⟩ = |1⟩\n"
    @test sprint(dirac, [1,1]) == "|ψ⟩ = |0⟩+|1⟩\n"
    @test sprint(dirac, [1,-1]) == "|ψ⟩ = |0⟩-|1⟩\n"
    @test sprint(dirac, [-1,1]) == "|ψ⟩ = -|0⟩+|1⟩\n"
    @test sprint(dirac, [-1,-1]) == "|ψ⟩ = -|0⟩-|1⟩\n"
    @test sprint(dirac, [2,3]) == "|ψ⟩ = 2|0⟩+3|1⟩\n"
    @test sprint(dirac, [-2,3]) == "|ψ⟩ = -2|0⟩+3|1⟩\n"
    @test sprint(dirac, [2,-3]) == "|ψ⟩ = 2|0⟩-3|1⟩\n"
    @test sprint(dirac, [-2,-3]) == "|ψ⟩ = -2|0⟩-3|1⟩\n"
    @test sprint(dirac, [1 0; 0 0]) == "ρ = |0⟩⟨0|\n"
    @test sprint(dirac, [0 1; 0 0]) == "ρ = |0⟩⟨1|\n"
    @test sprint(dirac, [0 0; 1 0]) == "ρ = |1⟩⟨0|\n"
    @test sprint(dirac, [0 0; 0 1]) == "ρ = |1⟩⟨1|\n"
    @test sprint(dirac, ones(Int, 2,2)) == "ρ = |0⟩⟨0|+|0⟩⟨1|+|1⟩⟨0|+|1⟩⟨1|\n"
    Random.seed!(0)
    @show rho = rand(-4:4, 2, 2)
    @test sprint(dirac, rho) == "ρ = -4|0⟩⟨0|+|0⟩⟨1|-2|1⟩⟨0|-4|1⟩⟨1|\n"
    rho[4] = 0
    @test sprint(dirac, rho) == "ρ = -4|0⟩⟨0|+|0⟩⟨1|-2|1⟩⟨0|\n"

    # sprint(dirac, [1,0])
    # sprint(dirac, [0,1])
    # sprint(dirac, [1,1])
    # sprint(dirac, [1,-1])
    # sprint(dirac, [-1,1])
    # sprint(dirac, [-1,-1])
    # sprint(dirac, [2,3])
    # sprint(dirac, [-2,3])
    # sprint(dirac, [2,-3])
    # sprint(dirac, [-2,-3])
    #
    # sprint(dirac, [1 0; 0 0])
    # sprint(dirac, [0 1; 0 0])
    # sprint(dirac, [0 0; 1 0])
    # sprint(dirac, [0 0; 0 1])
    # sprint(dirac, ones(Int, 2,2))
    #
    # Random.seed!(0)
    # @show rho = rand(-4:4, 2, 2)
    # sprint(dirac, rho)
    # rho[4] = 0
    # sprint(dirac, rho)

end # testset "Int"

@testset "Complex" begin

    @test sprint(dirac, Complex{Int}[1,0]) == "|ψ⟩ = |0⟩\n"
    @test sprint(dirac, Complex{Int}[0, 1]) == "|ψ⟩ = |1⟩\n"
    @test sprint(dirac, Complex{Int}[1, 1]) == "|ψ⟩ = |0⟩+|1⟩\n"
    @test sprint(dirac, Complex{Int}[1,1]) == "|ψ⟩ = |0⟩+|1⟩\n"
    @test sprint(dirac, Complex{Int}[1,-1]) == "|ψ⟩ = |0⟩-|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1,1]) == "|ψ⟩ = -|0⟩+|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1,-1]) == "|ψ⟩ = -|0⟩-|1⟩\n"
    @test sprint(dirac, Complex{Int}[1im,1im]) == "|ψ⟩ = im|0⟩+im|1⟩\n"
    @test sprint(dirac, Complex{Int}[1im,-1im]) == "|ψ⟩ = im|0⟩-im|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1im,1im]) == "|ψ⟩ = -im|0⟩+im|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1im,-1im]) == "|ψ⟩ = -im|0⟩-im|1⟩\n"
    @test sprint(dirac, Complex{Int}[1+1im,1+1im]) == "|ψ⟩ = (1+im)|0⟩+(1+im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[1+1im,1-1im]) == "|ψ⟩ = (1+im)|0⟩+(1-im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[1-1im,1+1im]) == "|ψ⟩ = (1-im)|0⟩+(1+im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[1-1im,1-1im]) == "|ψ⟩ = (1-im)|0⟩+(1-im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1+1im,-1+1im]) == "|ψ⟩ = (-1+im)|0⟩+(-1+im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1+1im,-1-1im]) == "|ψ⟩ = (-1+im)|0⟩+(-1-im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1-1im,-1+1im]) == "|ψ⟩ = (-1-im)|0⟩+(-1+im)|1⟩\n"
    @test sprint(dirac, Complex{Int}[-1-1im,-1-1im]) == "|ψ⟩ = (-1-im)|0⟩+(-1-im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1,0]) == "|ψ⟩ = |0⟩\n"
    @test sprint(dirac, Complex{Float64}[0, 1]) == "|ψ⟩ = |1⟩\n"
    @test sprint(dirac, Complex{Float64}[1, 1]) == "|ψ⟩ = |0⟩+|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1,1]) == "|ψ⟩ = |0⟩+|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1,-1]) == "|ψ⟩ = |0⟩-|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1,1]) == "|ψ⟩ = -|0⟩+|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1,-1]) == "|ψ⟩ = -|0⟩-|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1im,1im]) == "|ψ⟩ = im|0⟩+im|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1im,-1im]) == "|ψ⟩ = im|0⟩-im|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1im,1im]) == "|ψ⟩ = -im|0⟩+im|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1im,-1im]) == "|ψ⟩ = -im|0⟩-im|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1+1im,1+1im]) == "|ψ⟩ = (1.0+im)|0⟩+(1.0+im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1+1im,1-1im]) == "|ψ⟩ = (1.0+im)|0⟩+(1.0-im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1-1im,1+1im]) == "|ψ⟩ = (1.0-im)|0⟩+(1.0+im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[1-1im,1-1im]) == "|ψ⟩ = (1.0-im)|0⟩+(1.0-im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1+1im,-1+1im]) == "|ψ⟩ = (-1.0+im)|0⟩+(-1.0+im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1+1im,-1-1im]) == "|ψ⟩ = (-1.0+im)|0⟩+(-1.0-im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1-1im,-1+1im]) == "|ψ⟩ = (-1.0-im)|0⟩+(-1.0+im)|1⟩\n"
    @test sprint(dirac, Complex{Float64}[-1-1im,-1-1im]) == "|ψ⟩ = (-1.0-im)|0⟩+(-1.0-im)|1⟩\n"
    Random.seed!(0)
    @show x = rand(-4:4, 2,2) + im*rand(-4:4, 2,2)
    @test sprint(dirac, x) == "ρ = (-4-4im)|0⟩⟨0|+(1+4im)|0⟩⟨1|+(-2-im)|1⟩⟨0|+(-4+3im)|1⟩⟨1|\n"
    x[4] = 0
    @test sprint(dirac, x) == "ρ = (-4-4im)|0⟩⟨0|+(1+4im)|0⟩⟨1|+(-2-im)|1⟩⟨0|\n"
    x = randn(Complex{Float64}, 2,2)
    @test sprint(dirac, x) == "ρ = (-0.539384+0.281063im)|0⟩⟨0|+(-0.132634-1.1365im)|0⟩⟨1|+(0.573909-0.24491im)|1⟩⟨0|+(-1.75419+1.60954im)|1⟩⟨1|\n"

    # sprint(dirac, Complex{Int}[1,0])
    # sprint(dirac, Complex{Int}[0, 1])
    # sprint(dirac, Complex{Int}[1, 1])
    # sprint(dirac, Complex{Int}[1,1])
    # sprint(dirac, Complex{Int}[1,-1])
    # sprint(dirac, Complex{Int}[-1,1])
    # sprint(dirac, Complex{Int}[-1,-1])
    # sprint(dirac, Complex{Int}[1im,1im])
    # sprint(dirac, Complex{Int}[1im,-1im])
    # sprint(dirac, Complex{Int}[-1im,1im])
    # sprint(dirac, Complex{Int}[-1im,-1im])
    # sprint(dirac, Complex{Int}[1+1im,1+1im])
    # sprint(dirac, Complex{Int}[1+1im,1-1im])
    # sprint(dirac, Complex{Int}[1-1im,1+1im])
    # sprint(dirac, Complex{Int}[1-1im,1-1im])
    # sprint(dirac, Complex{Int}[-1+1im,-1+1im])
    # sprint(dirac, Complex{Int}[-1+1im,-1-1im])
    # sprint(dirac, Complex{Int}[-1-1im,-1+1im])
    # sprint(dirac, Complex{Int}[-1-1im,-1-1im])
    #
    # sprint(dirac, Complex{Float64}[1,0])
    # sprint(dirac, Complex{Float64}[0, 1])
    # sprint(dirac, Complex{Float64}[1, 1])
    # sprint(dirac, Complex{Float64}[1,1])
    # sprint(dirac, Complex{Float64}[1,-1])
    # sprint(dirac, Complex{Float64}[-1,1])
    # sprint(dirac, Complex{Float64}[-1,-1])
    # sprint(dirac, Complex{Float64}[1im,1im])
    # sprint(dirac, Complex{Float64}[1im,-1im])
    # sprint(dirac, Complex{Float64}[-1im,1im])
    # sprint(dirac, Complex{Float64}[-1im,-1im])
    # sprint(dirac, Complex{Float64}[1+1im,1+1im])
    # sprint(dirac, Complex{Float64}[1+1im,1-1im])
    # sprint(dirac, Complex{Float64}[1-1im,1+1im])
    # sprint(dirac, Complex{Float64}[1-1im,1-1im])
    # sprint(dirac, Complex{Float64}[-1+1im,-1+1im])
    # sprint(dirac, Complex{Float64}[-1+1im,-1-1im])
    # sprint(dirac, Complex{Float64}[-1-1im,-1+1im])
    # sprint(dirac, Complex{Float64}[-1-1im,-1-1im])
    # Random.seed!(0)
    # @show x = rand(-4:4, 2,2) + im*rand(-4:4, 2,2)
    # sprint(dirac, x)

end # testset "Complex"


text:
using HDF5

function save_Sequential(model::Sequential, path::String)
    h5open(path, "w") do file
        write(file, "num_layer", model.num_layer)
        for i in 2:model.num_layer+1
            model.layers[i].save_layer(model.layers[i], file, string(i-1))
        end
    end
end

function load_Sequential(path::String)
    model = Sequential()

    h5open(path, "r") do file
        num_layer = read(file, "num_layer")
        for i in 1:num_layer
            s = read(file, string(i))
            layer = Symbol(s)
            mod = Symbol(s*"M")
            args = eval(mod).get_args(file, string(i))
            try
                activation_function = Symbol(read(file, string(i)*"activation_function"))
                model.add_layer(model, eval(layer); args..., activation_function=eval(activation_function))
            catch KeyError
                model.add_layer(model, eval(layer); args...)
            end
            model.layers[end].load_layer(model.layers[end], file, string(i))
        end
    end
    return model
end

function save_GAN(model::GAN, path::String)
    h5open(path, "w") do file
        write(file, "num_Glayer", model.num_Glayer)
        write(file, "num_Dlayer", model.num_Dlayer)
        for i in model.G_range
            model.layers[i].save_layer(model.layers[i], file, string(i-1))
        end
        for i in model.D_range
            model.layers[i].save_layer(model.layers[i], file, string(i-1))
        end
    end
end

function load_GAN(path::String, noise_generator::Any)
    model = GAN(noise_generator)

    h5open(path, "r") do file
        num_Glayer = read(file, "num_Glayer")
        num_Dlayer = read(file, "num_Dlayer")
        for i in 1:num_Glayer
            s = read(file, string(i))
            layer = Symbol(s)
            mod = Symbol(s*"M")
            args = eval(mod).get_args(file, string(i))
            try
                activation_function = Symbol(read(file, string(i)*"activation_function"))
                model.add_Glayer(model, eval(layer); args..., activation_function=eval(activation_function))
            catch KeyError
                model.add_Glayer(model, eval(layer); args...)
            end
            model.layers[end].load_layer(model.layers[end], file, string(i))
        end

        for i in num_Glayer+1:num_Glayer+num_Dlayer
            s = read(file, string(i))
            layer = Symbol(s)
            mod = Symbol(s*"M")
            args = eval(mod).get_args(file, string(i))
            try
                activation_function = Symbol(read(file, string(i)*"activation_function"))
                model.add_Dlayer(model, eval(layer); args..., activation_function=eval(activation_function))
            catch KeyError
                model.add_Dlayer(model, eval(layer); args...)
            end
            model.layers[end].load_layer(model.layers[end], file, string(i))
        end
    end
    return model
end

ext: jl


text:
function syminit(app)

nd = app.nd;
nc = app.nc;
ncx = app.ncx;
ncu = app.ncu;
ncq = app.ncq;
ncp = app.ncp;
nce = app.nce;
nco = app.nco;
ncw = app.ncw;
ntau = length(app.tau);
nuinf = length(app.uinf);
nparam = length(app.physicsparam);

time = SymPy.symbols("time");

xdg = [SymPy.symbols("xdg$i") for i=1:ncx];

udg = [SymPy.symbols("udg$i") for i=1:nc];
udg1 = [SymPy.symbols("udg1$i") for i=1:nc];
udg2 = [SymPy.symbols("udg2$i") for i=1:nc];

wdg = [SymPy.symbols("wdg$i") for i=1:ncw];
wdg1 = [SymPy.symbols("wdg1$i") for i=1:ncw];
wdg2 = [SymPy.symbols("wdg2$i") for i=1:ncw];

odg = [SymPy.symbols("odg$i") for i=1:nco];
odg1 = [SymPy.symbols("odg1$i") for i=1:nco];
odg2 = [SymPy.symbols("odg2$i") for i=1:nco];

uhg = [SymPy.symbols("uhg$i") for i=1:ncu];
nlg = [SymPy.symbols("nlg$i") for i=1:nd];

tau = [SymPy.symbols("tau$i") for i=1:ntau];
uinf = [SymPy.symbols("uinf$i") for i=1:nuinf];
param = [SymPy.symbols("param$i") for i=1:nparam];

return xdg, udg, udg1, udg2, wdg, wdg1, wdg2, odg, odg1, odg2, uhg, nlg, tau, uinf, param, time

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using PolynomialBasis
using Test
# using Revise
using CutCellDG
include("../useful_routines.jl")

polyorder = 1
numqp = 2

basis = LagrangeTensorProductBasis(2,1)
dim,numpts = size(interpolation_points(basis))
mesh = CutCellDG.CGMesh([0.0, 0.0], [3.0, 1.0], [3, 1], numpts)

x0 = [1.5,0.0]
normal = [1.,0.]
distancefunc(x) = plane_distance_function(x,normal,x0)
levelset = CutCellDG.LevelSet(distancefunc,mesh,basis)
CutCellDG.load_coefficients!(levelset,1)
@test allapprox(levelset.interpolater.coeffs,[-1.5,-1.5,-0.5,-0.5]')
@test levelset([0.,0.]) ≈ -1.0

CutCellDG.update_coefficients!(levelset,2,[-0.6,-0.6,0.4,0.4])
@test allapprox(CutCellDG.coefficients(levelset,2),[-0.6,-0.6,0.4,0.4])
@test allapprox(CutCellDG.coefficients(levelset,1),[-1.5,-1.5,-0.6,-0.6])
@test allapprox(CutCellDG.coefficients(levelset,3),[0.4,0.4,1.5,1.5])

basis = HermiteTensorProductBasis(2)
numqp = required_quadrature_order(3)
quad = tensor_product_quadrature(2,numqp)
dim,numpts = size(interpolation_points(basis))
mesh = CutCellDG.CGMesh([0.,0.],[1.,1.],[2,2],numpts)
x0 = [0.75,0.0]
normal = [1.0,0.0]
distancefunc(x) = plane_distance_function(x,normal,x0)

levelset = CutCellDG.LevelSet(distancefunc,mesh,basis)

coeffs1 = CutCellDG.coefficients(levelset,1)
testcoeffs = [-0.75  -0.75  -0.25  -0.25
               0.0    0.0    0.0    0.0
               0.25   0.25   0.25   0.25
               0.0    0.0    0.0    0.0]
testcoeffs = vec(testcoeffs)
@test allapprox(coeffs1,testcoeffs,1e3eps())

x0 = [1.0,0.5]
normal = [1.0,1.0]/sqrt(2)

function testcubic(v)
    x,y = v
    return 8x^3*y^3 + 4x^2*y^3 + 17*x*y + 11
end

levelset = CutCellDG.LevelSet(testcubic,mesh,basis)
testcellid = 4
CutCellDG.load_coefficients!(levelset,testcellid)
cellmap = CutCellDG.cell_map(mesh,testcellid)
xrange = range(-1,stop=1,length=5)
testp = ImplicitDomainQuadrature.tensor_product_points(xrange',xrange')
spatialtestp = cellmap(testp)

vals = vec(mapslices(levelset,testp,dims=1))
testvals = vec(mapslices(testcubic,spatialtestp,dims=1))
@test allapprox(vals,testvals,1e4eps())


text:
using Plots
include(joinpath(@__DIR__, "quadrotor_broken_propeller.jl"))
include(joinpath(pwd(), "src/models/visualize.jl"))

X̄, Ū = unpack(Z̄, prob)

@show sum([Ū[t][end] for t = 1:T-1])
t_nominal = [0.0, [sum([Ū[i][end] for i = 1:t]) for t = 1:T-1]...]
plot(t_nominal, hcat(X̄...)[1:3,:]', linetype=:steppost, width = 2.0)
plot(t_nominal[1:end-1], hcat(Ū...)', linetype=:steppost, width = 2.0)

vis = Visualizer()
open(vis)
visualize!(vis, model, X̄, Δt = Ū[1][end])

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
module SymbolFactorGraphs

using LightGraphs

import LightGraphs: nv, ne, vertices, is_directed, has_edge, edges, inneighbors, outneighbors, rem_vertex!, neighbors, connected_components

import ...DistributedFactorGraphs: DFGNode

export
    SymbolFactorGraph,
    SymbolEdge
    # is_directed,
    # has_edge,
    # rem_vertex!,
    # addVariable!,
    # addFactor!

# import DistributedFactorGraphs: DFGNode
# const AbstractNodeType = DFGNode
import DistributedFactorGraphs: AbstractDFGVariable, AbstractDFGFactor
const AbstractVariableType = AbstractDFGVariable
const AbstractFactorType = AbstractDFGFactor

import Base: zero
const UnionNothingSymbol = Union{Nothing, Symbol}
zero(UnionNothingSymbol) = nothing

include("symboledge.jl")

"""
    SymbolFactorGraph{V, F}

A type representing an undirected bipartite factor graph based on `label::Symbol`.
"""
mutable struct SymbolFactorGraph{V <: AbstractVariableType, F <: AbstractFactorType} <: AbstractGraph{UnionNothingSymbol}
    ne::Int
    fadjdict::Dict{Symbol,Vector{Symbol}} # [variable src id]: (dst, dst, dst)
    variables::Dict{Symbol,V}
    factors::Dict{Symbol,F}
end

################################################################################
# overwrites from SimpleGraphs.jl

nv(g::SymbolFactorGraph) = length(g.fadjdict)
ne(g::SymbolFactorGraph) = g.ne

vertices(g::SymbolFactorGraph) = keys(g.fadjdict)

inneighbors(g::SymbolFactorGraph, v::Symbol) = (badj(g, v))

outneighbors(g::SymbolFactorGraph, v::Symbol) = (fadj(g, v))

neighbors(g::SymbolFactorGraph, v::Symbol) = (fadj(g, v))

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SymbolFactorGraph{V,F}}) where {V,F} = false

is_directed(::Type{SymbolFactorGraph}) = false

is_directed(g::SymbolFactorGraph) = false

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

add_edge!(g::SymbolFactorGraph, x, y) = add_edge!(g, edgetype(g)(x, y))
eltype(x::SymbolFactorGraph) = Symbol

fadj(g::SymbolFactorGraph) = g.fadjdict
fadj(g::SymbolFactorGraph, v::Symbol) = g.fadjdict[v]

edgetype(::SymbolFactorGraph) = SymbolEdge

# eltype(x::SymbolFactorGraph{V,F}) where {T,V,F} = T

# Graph{UInt8}(6), Graph{Int16}(7), Graph{UInt8}()
"""
    SymbolFactorGraph{V,F}(nv=0, nf=0)

Construct an empty `SymbolFactorGraph{V,F}` with `nv,nf` varaibles `sizehint!`, factors and 0 edges.
If not specified, the element type `T` is the type of `nv,nf`.

## Examples
```jldoctest
# SymbolFactorGraph{AbstractVariableType,AbstractFactorType}()
```
"""
function SymbolFactorGraph{V, F}(nv::Int=100, nf::Int=100) where {V <: AbstractVariableType, F <: AbstractFactorType}
    fadjdict = Dict{Symbol,Vector{Symbol}}()
        sizehint!(fadjdict, nv + nf)
    variables = Dict{Symbol,V}()
        sizehint!(variables, nv)
    factors = Dict{Symbol,F}()
        sizehint!(factors, nf)
    return SymbolFactorGraph{V, F}(0, fadjdict, variables, factors)
end

# SymbolFactorGraph{V, F}() where {V<:AbstractVariableType,F<:AbstractFactorType}= SymbolFactorGraph{Int, V, F}()

# SimpleGraph(UInt8)
"""
    SymbolFactorGraph(::Type{T},::Type{V}, ::Type{F})

Construct an empty `SymbolFactorGraph{T}` with 0 vertices and 0 edges.

## Examples
```jldoctest
julia> SymbolFactorGraph(Int,AbstractVariableType, AbstractFactorType)
SymbolFactorGraph{Int64,AbstractVariableType,AbstractFactorType}(0, 0, Union{Nothing, FGNode{AbstractVariableType,Int64}}[], Union{Nothing, FGNode{AbstractFactorType,Int64}}[])
```
"""
# SymbolFactorGraph(::Type{T}, ::Type{V}, ::Type{F}) where {T <: Integer, V<:AbstractVariableType,F<:AbstractFactorType} = SymbolFactorGraph{T, V, F}(zero(T), zero(T))

# SymbolFactorGraph(Int,AbstractVariableType, AbstractFactorType)

#= #############################################################################
# SimpleGraph of a SimpleGraph
"""
    SimpleGraph{T}(g::SimpleGraph)

Construct a copy of g.
If the element type `T` is specified, the vertices of `g` are converted to this type.
Otherwise the element type is the same as for `g`.

## Examples
```jldoctest
julia> g = complete_graph(5)
julia> SimpleGraph{UInt8}(g)
{5, 10} undirected simple UInt8 graph
```
"""
SimpleGraph(g::SimpleGraph) = copy(g)
=#

"""
    badj(g::SymbolFactorGraph[, v::Symbol])

Return the backwards adjacency list of a graph. If `v` is specified,
return only the adjacency list for that vertex.

"""
badj(g::SymbolFactorGraph) = fadj(g)
badj(g::SymbolFactorGraph, v::Symbol) = fadj(g, v)

"""
    adj(g[, v])

Return the adjacency list of a graph. If `v` is specified, return only the
adjacency list for that vertex.

"""
adj(g::SymbolFactorGraph) = fadj(g)
adj(g::SymbolFactorGraph, v::Symbol) = fadj(g, v)

#TODO
#=
# copy(g::SymbolFactorGraph) =  SymbolFactorGraph(g.ne, deepcopy_adjlist(g.fadjlist))

==(g::SimpleGraph, h::SimpleGraph) =
vertices(g) == vertices(h) &&
ne(g) == ne(h) &&
fadj(g) == fadj(h)
=#

function has_edge(g::SymbolFactorGraph, s::Symbol, d::Symbol)
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list_s = g.fadjdict[s]
    @inbounds list_d = g.fadjdict[d]
    if length(list_s) > length(list_d)
        d = s
        list_s = list_d
    end
    # return insorted(d, list_s)
    return in(d, list_s)
end

function has_edge(g::SymbolFactorGraph, e::SymbolEdge)
    s, d = Symbol.(Tuple(e))
    return has_edge(g, s, d)
end

"""
    add_edge!(g, e)

Add an edge `e` to graph `g`. Return `true` if edge was added successfully,
otherwise return `false`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> add_edge!(g, 1, 2)
true

julia> add_edge!(g, 2, 3)
false
```
"""
function add_edge!(g::SymbolFactorGraph, e::SymbolEdge)
    s, d = Tuple(e)
    s == d && error("Self loops not allowed")

    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjdict[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
    insert!(list, index, d)

    g.ne += 1

    @inbounds list = g.fadjdict[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    return true  # edge successfully added
end

"""
    rem_edge!(g, e)

Remove an edge `e` from graph `g`. Return `true` if edge was removed successfully,
otherwise return `false`.

### Implementation Notes
If `rem_edge!` returns `false`, the graph may be in an indeterminate state, as
there are multiple points where the function can exit with `false`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> add_edge!(g, 1, 2);

julia> rem_edge!(g, 1, 2)
true

julia> rem_edge!(g, 1, 2)
false
```
"""
function rem_edge!(g::SymbolFactorGraph, e::SymbolEdge)
    s, d = Tuple(e)
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjdict[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false  # edge not in graph
    deleteat!(list, index)

    g.ne -= 1
    s == d && return true  # selfloop

    @inbounds list = g.fadjdict[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    return true  # edge successfully removed
end

"""
    add_vertex!(g)

Add a new vertex to the graph `g`. Return `true` if addition was successful.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(Int8(typemax(Int8) - 1))
{126, 0} undirected simple Int8 graph

julia> add_vertex!(g)
true

julia> add_vertex!(g)
false
```
"""
# function add_vertex!(g::SymbolFactorGraph, node::) where T
#     (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
#     push!(g.fadjlist, Vector{T}())
#     return true
# end

function addVariable!(g::SymbolFactorGraph{V, F}, variable::V) where {V, F}
    (nv(g) + one(Int) <= nv(g)) && return false       # test for overflow
    # !isa(node,V) && (@warn "not a variable type"; return false)
    label = variable.label

    haskey(g.fadjdict, label) && (@error "Label already in fg"; return false) #TODO debug error or exception?
    # haskey(g.labelIndex, label) && error("Label already in fg")

    push!(g.fadjdict, label=>Symbol[])

    push!(g.variables, label=>variable)

    return true
end

function addFactor!(g::SymbolFactorGraph{V, F}, variableLabels::Vector{Symbol}, factor::F)::Bool where {V, F}

    haskey(g.fadjdict, factor.label) && (@error "Label $(factor.label) already in fg"; return false)

    for label in variableLabels
        !haskey(g.fadjdict, label) && (@error "Variable '$(label)' not found in graph when creating Factor '$(factor.label)'"; return false) #TODO debug error or exception?
    end

    (nv(g) + one(Int) <= nv(g)) && return false       # test for overflow
    # !isa(node,V) && (@warn "not a variable type"; return false)

    # factor._variableOrderSymbols = copy(variableLabels)

    push!(g.fadjdict, factor.label=>Symbol[])

    push!(g.factors, factor.label=>factor)

    # add the edges
    for label in variableLabels
        add_edge!(g, label, factor.label) || return false
    end

    return true
end

function edges(g::SymbolFactorGraph)

    fadjdict = fadj(g)
    n = nv(g)

    edges = SymbolEdge[]
    sizehint!(edges, n)

    for p = pairs(fadjdict)
        for l = p.second
            push!(edges,SymbolEdge(p.first=>l))
        end
    end

    edges

end

"""
    rem_vertex!(g, v)

Remove the vertex `v` from graph `g`. Return `false` if removal fails
(e.g., if vertex is not in the graph); `true` otherwise.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> rem_vertex!(g, 2)
true

julia> rem_vertex!(g, 2)
false
```
"""
function rem_vertex!(g::SymbolFactorGraph, v::Symbol)
    v in vertices(g) || return false
    n = nv(g)

    # TODO behoort nie moontelik te wees nie...
    # self_loop_n = false  # true if n is self-looped (see #820)

    # remove the in_edges from v
    srcs = copy(inneighbors(g, v))
    @inbounds for s in srcs
        rem_edge!(g, edgetype(g)(s, v))
    end

    # # remove the in_edges from the last vertex
    # neigs = copy(inneighbors(g, n))
    # @inbounds for s in neigs
    #     rem_edge!(g, edgetype(g)(s, n))
    # end
    # if v != n
    #     # add the edges from n back to v
    #     @inbounds for s in neigs
    #         if s != n  # don't add an edge to the last vertex - see #820.
    #             add_edge!(g, edgetype(g)(s, v))
    #         else
    #             self_loop_n = true
    #         end
    #     end
    # end

    # if self_loop_n
    #     add_edge!(g, edgetype(g)(v, v))
    # end
    delete!(g.fadjdict, v)
    delete!(g.variables, v)
    delete!(g.factors, v)

    return true
end

function _connected_components!(label::Dict{Symbol, T}, g::SymbolFactorGraph) where T <: UnionNothingSymbol
    nvg = nv(g)

    for u in vertices(g)
        label[u] != nothing && continue
        label[u] = u
        Q = Vector{T}()
        push!(Q, u)
        while !isempty(Q)
            src = popfirst!(Q)
            for vertex in neighbors(g, src)
                if label[vertex] == zero(T)
                    push!(Q, vertex)
                    label[vertex] = u
                end
            end
        end
    end
    return label
end

function _components(labels::Dict{Symbol, UnionNothingSymbol})

    d = Dict{Symbol,Int}()
    c = Vector{Vector{Symbol}}()
    i = 1 #one(T)
    for (v, l) in enumerate(labels)
        index = get!(d, l.second, i)
        if length(c) >= index
            push!(c[index], l.first)
        else
            push!(c, [l.first])
            i += 1
        end
    end
    return c, d
end

function connected_components(g::SymbolFactorGraph)
    # label = zeros(T, nv(g))
    label = Dict{Symbol, UnionNothingSymbol}(l=>nothing for l in vertices(g))
    _connected_components!(label, g)
    c, d = _components(label)
    return c
end

end


text:
"""
    set_param!(ref::ComponentReference, name::Symbol, value)

Set a component parameter as `set_param!(reference, name, value)`.
"""
function set_param!(ref::ComponentReference, name::Symbol, value)
    set_param!(ref.parent, ref.comp_path, name, value)
end

"""
    set_param!(ref.parent, ref.comp_name, name, value)

Set a component parameter as `reference[symbol] = value`.
"""
function Base.setindex!(ref::ComponentReference, value, name::Symbol)
    set_param!(ref.parent, ref.comp_path, name, value)
end

"""
    connect_param!(dst::ComponentReference, dst_name::Symbol, src::ComponentReference, src_name::Symbol)

Connect two components as `connect_param!(dst, dst_name, src, src_name)`.
"""
function connect_param!(dst::ComponentReference, dst_name::Symbol, src::ComponentReference, src_name::Symbol)
    connect_param!(dst.parent, dst.comp_path, dst_name, src.comp_path, src_name)
end

"""
    connect_param!(dst::ComponentReference, src::ComponentReference, name::Symbol)

Connect two components with the same name as `connect_param!(dst, src, name)`.
"""
function connect_param!(dst::ComponentReference, src::ComponentReference, name::Symbol)
    connect_param!(dst.parent, dst.comp_path, name, src.comp_path, name)
end

"""
    Base.getindex(comp_ref::ComponentReference, var_name::Symbol)

Get a variable reference as `comp_ref[var_name]`.
"""
function Base.getindex(comp_ref::ComponentReference, var_name::Symbol)
    VariableReference(comp_ref, var_name)
end

function _same_composite(ref1::AbstractComponentReference, ref2::AbstractComponentReference)
    # @info "same_composite($(ref1.comp_path), $(ref2.comp_path))"
    return ref1.comp_path.names[1] == ref2.comp_path.names[1]
end

"""
    Base.setindex!(comp_ref::ComponentReference, var_ref::VariableReference, var_name::Symbol)

Connect two components as `comp_ref[var_name] = var_ref`.
"""
function Base.setindex!(comp_ref::ComponentReference, var_ref::VariableReference, var_name::Symbol)
    _same_composite(comp_ref, var_ref)|| error("Can't connect variables defined in different composite trees")

    connect_param!(comp_ref.parent, comp_ref.comp_path, var_name, var_ref.comp_path, var_ref.var_name)
end


text:
const ContactImplicitMPC = Main
include(joinpath(@__DIR__, "..", "src/dynamics", "pushbot", "visuals.jl"))
vis = Visualizer()
# render(vis)
open(vis)

s = get_simulation("pushbot", "flat_2D_lc", "flat")
model = s.model
env = s.env

# time
h = 0.04
H = 100

# reference trajectory
ref_traj = contact_trajectory(model, env, H, h)
ref_traj.h
qref = [0.0; 0.0]
ur = zeros(model.dim.u)
γr = zeros(model.dim.c)
br = zeros(model.dim.c * friction_dim(env))
ψr = zeros(model.dim.c)
ηr = zeros(model.dim.c * friction_dim(env))
wr = zeros(model.dim.w)

# set reference
for t = 1:H
	ref_traj.z[t] = pack_z(model, env, qref, γr, br, ψr, ηr)
	ref_traj.θ[t] = pack_θ(model, qref, qref, ur, wr, model.μ_world, ref_traj.h)
end

# test reference
for t = 1:H
	r = ContactImplicitMPC.residual(model, env, ref_traj.z[t], ref_traj.θ[t], 0.0)#[model.dim.q .+ (1:model.dim.c)]
	@test norm(r) < 1.0e-4
end

# initial conditions
q0 = @SVector [0.0 * π, 0.0]
q1 = @SVector [0.0 * π, 0.0]

# simulator
sim = ContactImplicitMPC.simulator(s, q0, q1, h, H,
	ip_opts = ContactImplicitMPC.InteriorPointOptions(
		r_tol = 1.0e-6, κ_tol = 1.0e-5),
	sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

# simulate
status = ContactImplicitMPC.simulate!(sim)
@test status

# MPC
N_sample = 2
H_mpc = 40
h_sim = h / N_sample
H_sim = 1000

# barrier parameter
κ_mpc = 1.0e-4

# SLOW RECOVERY
obj = TrackingVelocityObjective(model, env, H_mpc,
	q = [Diagonal([12*(t/H_mpc)^2; 2.0*(t/H_mpc)^4]) for t = 1:H_mpc-0],
	v = [Diagonal([1; 0.01] ./ (h^2.0)) for t = 1:H_mpc-0],
	u = [Diagonal([100; 1]) for t = 1:H_mpc-0],
	γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc-0],
	b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])
# FAST RECOVERY
obj = TrackingVelocityObjective(model, env, H_mpc,
    q = [Diagonal([12*(t/H_mpc)^2; 12*(t/H_mpc)^2]) for t = 1:H_mpc-0],
	v = [Diagonal([1; 0.01] ./ (h^2.0)) for t = 1:H_mpc-0],
    u = [Diagonal([100; 1]) for t = 1:H_mpc-0],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc-0],
    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = NewtonOptions(
		r_tol = 3e-4,
		max_iter = 10,
		max_time = ref_traj.h/2, # HARD REAL TIME
		# verbose = true,
		),
    mpc_opts = LinearizedMPCOptions(),
    )

idx_d1 = 20
idx_d2 = idx_d1 + 200
idx_d3 = idx_d2 + 80
idx_d4 = idx_d3 + 200
idx_d5 = idx_d4 + 30
idx = [idx_d1, idx_d2, idx_d3, idx_d4, idx_d5]
impulses = [[-5.5; 0.0], [+5.5; 0.0], [+5.5; 0.0], [-1.5; 0.0], [-6.5; 0.0]]
d = impulse_disturbances(impulses, idx)

q1_sim = SVector{model.dim.q}([0.0, 0.0])
q0_sim = SVector{model.dim.q}([0.0, 0.0])

sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
	d = d,
    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

status = ContactImplicitMPC.simulate!(sim, verbose = true)
# @profiler status = ContactImplicitMPC.simulate!(sim, verbose = true)

################################################################################
# Timing result
################################################################################
process!(sim)
# Time budget
ref_traj.h # 0.04
# Time used on average
sim.stats.μ_dt # 0.0138
sim.stats.σ_dt # 0.0138
# Speed ratio
H_sim * h_sim / sum(sim.stats.dt) # 2.90

l = 1
lu = 1
nc = model.dim.c
plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.([q[l:l] for q in sim.traj.q])...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][lu:lu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[3,1], hcat(Vector.(vcat([fill(ref_traj.γ[i][1:nc], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[lu:lu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)

sample = 1

anim = visualize_robot!(vis, model, sim.traj, sample = 1)#, sample = sample, α=1.0)
pθ_right = generate_pusher_traj(d, sim.traj, side=:right)
pθ_left  = generate_pusher_traj(d, sim.traj, side=:left)
visualize_disturbance!(vis, model, pθ_right, anim=anim, sample=sample, offset=0.05, name=:PusherRight)
visualize_disturbance!(vis, model, pθ_left,  anim=anim, sample=sample, offset=0.05, name=:PusherLeft)

γ_max = maximum(hcat(sim.traj.γ...))
u_max = maximum(hcat(sim.traj.u...))

plot((hcat(sim.traj.γ...) ./ γ_max)', linetype = :steppost)
plot!((hcat(sim.traj.u...) ./ u_max)[2:2, :]', linetype = :steppost)
plot((hcat(sim.traj.u...) ./ u_max)[1:1, :]', linetype = :steppost)
plot(hcat(sim.traj.q...)')

plot([pθ[1] for pθ in pθ_right])
plot!([pθ[1] for pθ in pθ_left])
plot!([q[1] for q in sim.traj.q[3:end]])

# Display highlights
t_highlights = [20,23,35,216]
α_highlights = ones(4)
α_pusher = [1.0, 0.0, 0.0, 0.0]

# # Initial push
# i = 1
# t = t_highlights[i]
# α = i/(length(t_highlights)+1)
# name = Symbol("High$i")
# pusher_name = Symbol("HighPusherRight$i")
# build_robot!(vis, model, name=name, α=α_highlights[i])
# build_disturbance!(vis, model, name=pusher_name, α=α_highlights[i])
# set_robot!(vis, model, sim.traj.q[t+2], name=name)
# set_disturbance!(vis, model, pθ_left[t], name=pusher_name, offset=0.05)

# # Contact with wall
# i = 2
# t = t_highlights[i]
# α = i/(length(t_highlights)+1)
# name = Symbol("High$i")
# pusher_name = Symbol("HighPusherRight$i")
# build_robot!(vis, model, name=name, α=α_highlights[i])
# build_disturbance!(vis, model, name=pusher_name, α=α_pusher[i])
# set_robot!(vis, model, sim.traj.q[t+2], name=name)
# set_disturbance!(vis, model, pθ_left[t], name=pusher_name, offset=0.05)

# # Recovery main body
# i = 3
# t = t_highlights[i]
# α = i/(length(t_highlights)+1)
# name = Symbol("High$i")
# pusher_name = Symbol("HighPusherRight$i")
# build_robot!(vis, model, name=name, α=α_highlights[i])
# build_disturbance!(vis, model, name=pusher_name, α=α_pusher[i])
# set_robot!(vis, model, sim.traj.q[t+2], name=name)
# set_disturbance!(vis, model, pθ_left[t], name=pusher_name, offset=0.05)

# Recovery hand
i = 4
t = t_highlights[i]
α = i/(length(t_highlights)+1)
name = Symbol("High$i")
pusher_name = Symbol("HighPusherRight$i")
build_robot!(vis, model, name=name, α=α_highlights[i])
build_disturbance!(vis, model, name=pusher_name, α=α_pusher[i])
set_robot!(vis, model, sim.traj.q[t+2], name=name)
set_disturbance!(vis, model, pθ_left[t], name=pusher_name, offset=0.05)

# filename = "pushbot_multiple_pushes2"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/video/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/video/$filename.mp4",
#     "/home/simon/Documents/video/$filename.gif", overwrite=true)


text:
module QJuliaGMRes

using QJuliaInterface
using QJuliaEnums
using QJuliaBlas
using QJuliaReduce
using QJuliaSolvers

using LinearAlgebra
using MPI
using Printf

norm2   = QJuliaReduce.gnorm2
rdot    = QJuliaReduce.reDotProduct
verbose = true

const ε = eps(Float64)

# j - current row index , k - current column index

@inline function HessenbergQR(Givens::Any, H::Any, j, k; compute_last_element = false)
    #extract component
    cs   = Givens[1]
    sn   = Givens[2]
    Rm1m = Givens[3]
    γ    = Givens[4]

    if compute_last_element == true # last entry in the column (all columns), that means j = k
      if (j != k) error("Wrong row index. Must be diagonal element."); end
      inv_denom = 1.0 / sqrt(sn[k]*sn[k]+H[j+1,j]*H[j+1,j])
      old_snj   = sn[k]
      cs[j]     = old_snj * inv_denom; sn[j] = H[j+1,j] * inv_denom
      Rm1m[j,j] = cs[j]*old_snj + sn[j]*H[j+1,j]
      return
    end

    if j == 1 # initialization stage for all columns
      sn[k] = H[j,j] # use sn[j] is a temp variable
      return
    end

    Rm1m[j-1, k] = cs[j-1]*sn[k] + sn[j-1]*H[j,k]
    sn[k] = -sn[j-1]*sn[k] + cs[j-1]*H[j,k]

end

@inline function PlainArnoldiCycle(Mat::Any, Givens::Any, H::Any, Vk1::Any, k, param; do_givens_qr = true)

    k1 = k+1
    σ  = param.shift

    if k < 1; error("Wrong vector dimension in the Arnodli step."); end

    w = view(Vk1, :, k1)
    v = view(Vk1, :, k )

    Mat(w, v)

    # Apply shifts:
    w .=@. w - σ*v

    for j in 1:k
      H[j,k] = rdot(Vk1[:,k1], Vk1[:,j])
      Vk1[:,k1] .=@. Vk1[:,k1] - H[j,k] * Vk1[:,j]

      if do_givens_qr == true; HessenbergQR(Givens, H, j, k); end
    end

    H[k1,k] = norm(Vk1[:,k1])
    Vk1[:,k1] .=@. 1.0 / H[k1,k] * Vk1[:,k1]

    if do_givens_qr == true; HessenbergQR(Givens, H, k, k; compute_last_element = true); end

    return
end

@inline function PipelinedArnoldiCycle_l1(Mat::Any, Givens::Any, H::Any, Vk1::Any, k, param; do_givens_qr = false)

    k1 = k+1
    σ  = param.shift

    if k < 1; error("Wrong vector dimension in the Arnodli step."); end

    v = view(Vk1, :, k )
    u = view(Vk1, :, k1)
    z = zeros(typeof(Vk1[1,k]), length(Vk1[:,k]))

    Mat(z, v)

    # Apply shifts:
    z .=@. z - σ*v

    # Prepare subdiag. element of Hess. matrix and the last Arnoldi vector:
    H[k1, k] = norm2(z)
    u .=@. z

    for j in 1:k
      vj        = view(Vk1, :, j)
      # WARNING: vj is a subarray so rdot z, vj won't work
      H[j,k]    = dot(z, vj)

      H[k1, k] -= H[j,k]*H[j,k]
      u .=@. u - H[j,k]*vj

      if j == k; H[j,k] += σ; end
      if do_givens_qr == true; HessenbergQR(Givens, H, j, k); end
    end

    if H[k1, k] <= 0.0; error("SQRT breakdown detected at Arnoldi cycle ", k, " ."); end

    H[k1,k] = sqrt(H[k1,k])
    u .=@. u / H[k1,k]

    return
end

@inline function PipelinedArnoldiCycle_l1_step1(Mat::Any, Givens::Any, H::Any, Vk::Any, k, param; do_givens_qr = false)

    k1      = k+1
    σ       = param.shift
    nkrylov = param.nKrylov

    if k < 1; error("Wrong vector dimension in the Arnodli step."); end

    Zk1 = view(Vk, :, nkrylov+2:2nkrylov+3)

    if k == 1; Zk1[:, k] .=@. Vk[:, k]; end

    z = view(Zk1, :, k )
    u = zeros(typeof(Vk[1,k]), length(Vk[:,k]))
    w = zeros(typeof(Zk1[1,k]), length(Zk1[:,k]))

    Mat(w, z)

    znorm2 = norm2(z)

    if k == 1 #compute intital values
      znorm = sqrt(znorm2)
      Vk[:,k]   .=@. 1.0 / znorm * Vk[:,k]
      Zk1[:,k1] .=@. 1.0 / znorm * w
      return;
    end

    # Prepare subdiag. element of Hess. matrix and the last Arnoldi vector:
    H[k, k-1] = znorm2
    u .=@. z

    for j in 1:(k-1)
      vj   = view(Vk, :, j  )
      zj1  = view(Zk1,:, j+1)
      # WARNING: vj is a subarray so rdot z, vj won't work
      H[j,k-1] = dot(z, vj)

      H[k, k-1] -= H[j,k-1]*H[j,k-1]
      u .=@. u - H[j,k-1]*vj
      w .=@. w - H[j,k-1]*zj1

      if j == k; H[j,k-1] += σ; end
      if do_givens_qr == true; HessenbergQR(Givens, H, j, k); end
    end

    if H[k, k-1] <= 0.0; error("SQRT breakdown detected at Arnoldi cycle ", k, " ."); end

    H[k,k-1] = sqrt(H[k,k-1])

    Vk[:,k]   .=@. u / H[k,k-1]
    Zk1[:,k1] .=@. w / H[k,k-1]

    return
end

@inline function PipelinedArnoldiCycle_p1(Mat::Any, Givens::Any, H::Any, Vk::Any, k, param; do_givens_qr = false)

    σ       = param.shift
    nkrylov = param.nKrylov

    if k < 1; error("Wrong vector dimension in the Arnodli step."); end
    if do_givens_qr == true; error("Givens QR is not supported."); end

    Zkp1 = view(Vk, :, nkrylov+2:2nkrylov+3)

    if k == 1; Zkp1[:, k] .=@. Vk[:, k]; end
    if k == (nkrylov + 2)
      Vk[:, k-1] .=@. Vk[:, k-1] / H[k-1, k-2]
      return;
    end

    zk   = view(Zkp1, :, k  )
    zkp1 = view(Zkp1, :, k+1)
    u = zeros(typeof(Vk[1,k]), length(Vk[:,k]))

    Mat(zkp1, zk)

    if k > 2
      vkm1  = view(Vk, :, k-1)
      vkm1 .=@. vkm1 / H[k-1,k-2]
      zk   .=@. zk   / H[k-1,k-2]
      zkp1 .=@. zkp1 / H[k-1,k-2]

      for j in 1:(k-2); H[j,k-1] /= H[k-1,k-2]; end

      H[k-1,k-1] /= (H[k-1,k-2]*H[k-1,k-2])
      H[k-1,k-1] += σ #must be sigma[k-1] !
    end

    if k > 1
      vk = view(Vk, :, k)
      vk .=@. zk
      for j in 1:(k-1)
        vj    = view(Vk, :, j)
        vk   .=@. vk - H[j,k-1]*vj

        if k > nkrylov; continue; end  # in this case we don't need the next Z-vector

        zjp1  = view(Zkp1, :, j+1)
        zkp1 .=@. zkp1 -  H[j,k-1]*zjp1
      end
      H[k, k-1] = norm(vk)
    end

    if k > nkrylov; return; end

    for j in 1:k
      vj   = view(Vk, :, j)
      # WARNING: vj is a subarray so rdot z, vj won't work
      H[j,k] = dot(zkp1,vj)
    end

    return
end

const do_givens = false
#TODO: givens version is broken
const pipeline_algo_type = 3

function solver(x::AbstractArray, b::AbstractArray, Mat::Any, MatSloppy::Any, param::QJuliaSolvers.QJuliaSolverParam_qj)

    if verbose == true; println("Running GMRES in maximum " , param.nKrylov, " iters."); end

    QJuliaReduce.set_blas_global_reduction(param.global_reduction)

    if MPI.Initialized() == false; error("MPI was not inititalized, copy source field to the solution."); end

    if param.nKrylov > param.maxiter; println("Warning: search space exceeds requested maximum iter count. Some routines may fail."); end

    if (param.maxiter == 0) || (param.nKrylov == 0)
      if param.use_init_guess == false
        x .=@. 0.0
      end
      return
    end #if

    mixed = (param.dtype_sloppy != param.dtype)

    σ = param.shift
    nkrylov = param.nKrylov
    data_type = typeof(x[1]) == BigFloat ? BigFloat : Float64

    local r    = zero(typeof(x)(undef, length(x)))
    # Arnoldi vectors
    set_size    = pipeline_algo_type < 2 ? nkrylov+1 : (nkrylov+1) + (nkrylov+2)
    local Vm1   = zero(Matrix{typeof(x[1])}(undef, length(x), set_size))

    if do_givens == true
      # Givens objects:
      # R matrix from givens QR of H (use triu() for explicit upper triangular view or tril for lower triangular view)
      local Rm1m  = zero(Matrix{BigFloat}(undef, (nkrylov+1), nkrylov))
      # Givens matrix coeffs
      local cs    = zeros(BigFloat, nkrylov)
      local sn    = zeros(BigFloat, nkrylov)
    end

    # Hessenberg matrix
    local Hm1m  = zero(Matrix{data_type}(undef, (nkrylov+1), nkrylov))

    bnorm = sqrt(norm2(b))  #Save norm of b
    rnorm = 0.0     #if zero source then we will exit immediately doing no work

    v0 = view(Vm1, :, 1)

    if param.use_init_guess == true
      #r = b - Ax0
      Mat(v0, x)
      v0 .=@. v0 - σ*x
      v0 .=@. b - v0
      rnorm = sqrt(norm2(v0))
    else
      rnorm = bnorm
      v0 .=@. b
    end

    if( rnorm == 0.0 ); return; end
    if verbose == true; println("GMRES: Initial residual = ", rnorm); end

    # Scale the first basis vector:
    if pipeline_algo_type != 2;  v0 .=@. (1.0 / rnorm) * v0; end
    # Aux coeffs
    β = rnorm

    # LS problem source vector
    γ  = zeros(BigFloat, nkrylov+1); γ[1] = β

    # Single object for all components
    Gm1m =  do_givens == true ? (cs, sn, Rm1m, γ) : nothing

    #start iterations
    ss = min(param.maxiter, param.nKrylov)

    if pipeline_algo_type == 0
      println("Running plain Arnoldi cycles.")
      ArnoldiCycle = PlainArnoldiCycle
    elseif pipeline_algo_type == 1
      println("Running pipelined Arnoldi cycles.")
      ArnoldiCycle = PipelinedArnoldiCycle_l1
    elseif pipeline_algo_type == 2
      println("Running pipelined Arnoldi cycles (s1 version).")
      ArnoldiCycle = PipelinedArnoldiCycle_l1_step1
      ss += 1   #+1 iteration for this type
    elseif pipeline_algo_type == 3
      println("Running pipelined Arnoldi cycles (p1 version).")
      ArnoldiCycle = PipelinedArnoldiCycle_p1
      ss += 2   #+2 iteration for this type
    else
      error("Arnoldi step is not selected.")
    end

    for k in 1:ss
      # Arnoldi step + Givens QR:
      ArnoldiCycle(MatSloppy, Gm1m, Hm1m, Vm1, k, param; do_givens_qr = do_givens)
      # Update LS source vector and compute iter residual
      if do_givens == true
        γ[k+1] = -sn[k]*γ[k]
        γ[k]   = γ[k]*cs[k]
        rnorm  = abs(γ[k+1])
      end

      if(verbose == true && do_givens == true); println("GMRes: ", k,"  residual norm = ", rnorm); end
    end #for itern
    # compute solution vector:
    if do_givens == false
      # native QR method (access via .Q and .R )
      qrH = qr(Hm1m)
      Rmm = qrH.R
      Qγ  = transpose(qrH.Q) * γ
      γ_  = view(Qγ,1:nkrylov)
    else
      Rmm = triu(view(Rm1m, 1:nkrylov, :))
      γ_  = view(γ,1:nkrylov)
    end
    η   = Rmm \ γ_

    for i in 1:nkrylov
      x .=@. x + Vm1[:, i]*η[i]
    end

    # Compute true residual
    Mat(r, x)

    r .=@. b - r
    rnorm = sqrt(norm2(r))

    param.true_res = rnorm / bnorm
    if verbose == true; println("GMRes: Relative residual: true = ", param.true_res); end

    println("Hessenberg matrix\n\n")
    for ii in 1:(nkrylov+1)
      for jj in 1:2
        @printf(": %1.10e ",Hm1m[ii,jj] )
      end
      println("\n------------------------------\n")
    end

    QJuliaReduce.reset_blas_global_reduction()

end #MR

end #QJuliaLanMR

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
using Printf
using GLMNet: glmnet

"""
    screen_lasso(X, y, pmax)

Screening out `pmax` variables by lasso (glmnet) on `X` and `y`.
"""
function screen_lasso(X::Matrix{Float64}, y::Vector{Float64}, pmax::Int64)
    fit = glmnet(X, y)
    p = size(X, 2)
    @assert pmax <= p
    betas = fit.betas
    n_vars = sum(betas.!=0, dims=1)[:]
    S = Set{Int64}()
    for n in unique(n_vars)
        Z = betas[:, n_vars.==n]
        non_zeros = sum(Z.!=0, dims=2)[:]
        vars = (1:p)[non_zeros .> 0]
        new_vars = setdiff(vars, S)
        Z = Z[new_vars, :]
        if length(S) + length(new_vars) > pmax
            # must break ties
            z = abs.(Z[:, end])
            idx = sort(1:length(z), by = i -> z[i], rev = true)
            new_vars =  new_vars[idx[1:(pmax - length(S))]]
        end
        union!(S, new_vars)
        length(S) == pmax && break
    end
    return sort(collect(S))
end

"""
    screen_HOLP(X, y, pmax)

Screen out `pmax` variables with HOLP projection.

See: Wang, Xiangyu, and Chenlei Leng. "High dimensional ordinary least squares projection for screening variables."
     Journal of the Royal Statistical Society: Series B (Statistical Methodology) 78.3 (2016): 589-611.

"""
function screen_HOLP(X::Matrix{Float64}, y::Vector{Float64}, pmax::Int64)
    n, p = size(X)
    @assert p >= n
    @assert pmax < p
    _X = X .- mean(X, 1)
    _X = _X ./ std(_X, 1)
    _y = y - mean(y)
    β = _X' * ((_X * _X' + 10 * eye(n)) \ _y)
    return sortperm(abs.(β), rev=true)[1:pmax]
end

max_issues_repo_name: richardkwo/InvariantCausal

max_forks_repo_name: richardkwo/InvariantCausal


text:
# Julia wrapper for the sgp4 Python library:
# https://pypi.python.org/pypi/sgp4/

VERSION >= v"0.4.0-dev+6521" && __precompile__()

module SGP4

import Compat.ASCIIString

using PyCall

function __init__()
    try
        global const sgp4io = PyCall.pywrap(pyimport("sgp4.io"))
        global const earth_gravity = PyCall.pywrap(pyimport("sgp4.earth_gravity"))
    catch e
        error("Error loading sgp4 python package - check to make sure it's installed")
    end
end

export GravityModel,
       twoline2rv,
       propagate

immutable GravityModel
    model::PyObject # can be any of {wgs72old, wgs72, wgs84}
end
GravityModel(ref::AbstractString) = earth_gravity.pymember(ref)

# sgp4.io convenience functions
twoline2rv(args...) = sgp4io.twoline2rv(args...)

"""
Propagate the satellite from its epoch to the date/time specified

Returns (position, velocity) at the specified time
"""
function propagate( sat::PyObject,
                    year::Real,
                    month::Real,
                    day::Real,
                    hour::Real,
                    min::Real,
                    sec::Real )
    (pos, vel) = sat[:propagate](year,month,day,hour,min,sec)

    # check for errors
    if sat[:error] != 0
        println(sat[:error_message])
    end

    return ([pos...],[vel...])
end

function propagate( sat::PyObject,
                    t::DateTime )
    propagate(sat,
              Dates.year(t),
              Dates.month(t),
              Dates.day(t),
              Dates.hour(t),
              Dates.minute(t),
              Dates.second(t))
end

"Propagate many satellites to a common time"
function propagate( sats::Vector{PyObject},
                    year::Real,
                    month::Real,
                    day::Real,
                    hour::Real,
                    min::Real,
                    sec::Real )
    pos = zeros(3,length(sats))
    vel = zeros(3,length(sats))
    for i = 1:length(sats)
        (pos[:,i],vel[:,i]) = propagate(sats[i],year,month,day,hour,min,sec)
    end
    return (pos,vel)
end

function propagate( sats::Vector{PyObject},
                    t::DateTime )
    propagate(sats,
              Dates.year(t),
              Dates.month(t),
              Dates.day(t),
              Dates.hour(t),
              Dates.minute(t),
              Dates.second(t))
end

end #module


text:
module Pio3d
using PyCall

greet() = print("Hello World!")

"""
Read 3D data from file. Data are stored int dictionary with keys "data3d" and "voxelsize_mm".
"""
function read3d(fn)
    pyIo3d = pyimport("io3d")
    return pyIo3d."read"(fn; dataplus_format=true);
end

"""
Read 3D data from file. Data are stored int dictionary with keys "data3d" and "voxelsize_mm".
"""
function read(fn)
    return read3d(fn)
end

function datasets_join_path(pth)
    pyIo3d_datasets = pyimport("io3d.datasets")
    opth = pyIo3d_datasets."join_path"(pth, get_root=true)
    return opth
end

"""
Download test data based on input label. Available labels can be list by get_labels()
"""
function datasets_download(label)
    pyIo3d_datasets = pyimport("io3d.datasets")
#     print(pth)
    pyIo3d_datasets."download"(label)
end

"""
Write data into file. Data have to be in dict with labels "data3d" and "voxelsize_mm".
"""
function write3d(arr, fn)
    pyIo3d = pyimport("io3d")
    return pyIo3d."write"(arr, fn);
#     print(pth)
end

"""
Get list of available datasets labels.
"""
function get_labels()
    pyIo3d_datasets = pyimport("io3d.datasets")
    labels = pyIo3d_datasets."get_labels".()
    return labels
end

end # module


text:
function sample_fractions(tree, sampling_at_tips)
    function aux(subtree, sub_sampling, x)
        if subtree.n_nodes < 2
            s = pop!(sub_sampling)
            pushfirst!(x,s)
            return s
        else
            s_left = aux(subtree.offsprings[2], sub_sampling,x)
            n_left = (subtree.offsprings[2].n_nodes+1)/2
            s_right = aux(subtree.offsprings[1], sub_sampling,x)
            n_right = (subtree.offsprings[1].n_nodes+1)/2
            s = (n_left + n_right)/(n_left/s_left + n_right/s_right)
            pushfirst!(x,s)

            return s
        end
    end

    fs = Array{Float64,1}(undef,0)
    aux(tree, deepcopy(sampling_at_tips), fs)
    return fs
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA