text:
@testset "Configurations" begin
    @testset "Construction" begin
        config = Configuration([o"1s", o"2s", o"2p", o"3s", o"3p"], [2,2,6,2,6], [:closed])
        rconfig = Configuration([ro"1s", ro"2s", ro"2p", ro"3s", ro"3p"], [2,2,6,2,6], [:closed], sorted=true)

        @test orbitals(config) ==
            [o"1s", o"2s", o"2p", o"3s", o"3p"]
        @test orbitals(rconfig) ==
            [ro"1s", ro"2s", ro"2p-", ro"2p", ro"3s", ro"3p-", ro"3p"]

        @test config.occupancy == [2, 2, 6, 2, 6]
        @test rconfig.occupancy == [2, 2, 2, 4, 2, 2, 4]

        @test config.states == [:closed, :open, :open, :open, :open]
        @test rconfig.states == [:closed, :open, :open, :open, :open, :open, :open]

        @test c"1s2c 2s2 2p6 3s2 3p6" == config
        @test c"1s2c.2s2.2p6.3s2.3p6" == config
        @test c"[He]c 2s2 2p6 3s2 3p6" == config

        # We sort the configurations since the non-relativistic
        # convenience labels (2p) &c expand to two orbitals each, one
        # of which is appended to the orbitals list.
        @test rc"1s2c 2s2 2p6 3s2 3p6"s == rconfig
        @test rc"1s2c.2s2.2p6.3s2.3p6"s == rconfig
        @test rc"[He]c 2s2 2p6 3s2 3p6"s == rconfig

        @test c"[Kr]c 5s2" == Configuration([o"1s", o"2s", o"2p", o"3s", o"3p", o"3d", o"4s", o"4p", o"5s"],
                                            [2,2,6,2,6,10,2,6,2],
                                            [:closed,:closed,:closed,:closed,:closed,:closed,:closed,:closed,:open])

        @test length(c"") == 0
        @test length(rc"") == 0

        @test rc"1s ld-2 kp6" == Configuration([ro"1s", ro"ld-", ro"kp", ro"kp-"], [1, 2, 4, 2])
        @test rc"1s ld-2 kp6"s == Configuration([ro"1s", ro"kp-", ro"kp", ro"ld-"], [1, 2, 4, 2])

        @test_throws ArgumentError parse(Configuration{Orbital}, "1sc")
        @test_throws ArgumentError parse(Configuration{Orbital}, "1s 1s")
        @test_throws ArgumentError parse(Configuration{Orbital}, "[He]c 1s")
        @test_throws ArgumentError parse(Configuration{Orbital}, "1s3")
        @test_throws ArgumentError parse(Configuration{RelativisticOrbital}, "1s3")
        @test_throws ArgumentError parse(Configuration{Orbital}, "1s2 2p-2")

        @testset "Unicode occupations/states" begin
            for c in [c"1s2 2p6", c"[He]c 2p6", c"[He]i 2p6", c"1s 3d4", c"[Xe]*"]
                @test parse(Configuration{Orbital}, string(c)) == c
            end

            for c in [rc"1s2 2p6", rc"1s 2p- 2p3", rc"[He]c 2p6", rc"[He]i 2p6", rc"1s 3d4", rc"[Xe]*"]
                @test parse(Configuration{RelativisticOrbital}, string(c)) == c
            end
        end

        @testset "Spin-configurations" begin
            a = sc""
            @test a isa SpinConfiguration{<:SpinOrbital{<:Orbital}}
            @test isempty(a)

            b = rsc""
            @test b isa SpinConfiguration{<:SpinOrbital{<:RelativisticOrbital}}
            @test isempty(b)

            @test sc"1s(0,α) 2p(1,β)" == Configuration([so"1s(0,α)", so"2p(1,β)"], ones(Int, 2))
            @test sc"[He]c 2p(1,β)" == Configuration([so"1s(0,α)", so"1s(0,β)", so"2p(1,β)"], ones(Int, 3), [:closed, :closed, :open])
            @test rsc"[He]c 2p(-1/2)" == Configuration([rso"1s(-1/2)", rso"1s(1/2)", rso"2p(-1/2)"], ones(Int, 3), [:closed, :closed, :open])

            # Test parsing of Unicode m_ℓ quantum numbers
            @test sc"1s₀α 2p₁β" == sc"1s(0,α) 2p(1,β)"

            c = rsc"[Xe]*"
            @test parse(SpinConfiguration{SpinOrbital{RelativisticOrbital}}, string(c)) == c
        end

        @test fill(c"1s 2s 2p") == c"1s2 2s2 2p6"
        @test fill(rc"1s 2s 2p- 2p") == rc"1s2 2s2 2p-2 2p4"
        @test close(c"1s2") == c"1s2c"
        let c = c"1s2c 2s 2p"
            @test fill!(c) == c"1s2c 2s2 2p6"
            @test c == c"1s2c 2s2 2p6"
            close!(c)
            @test sort(c) == c"[Ne]"s
        end
        let c = rc"1s2c 2s 2p- 2p2"
            @test fill!(c) == rc"1s2c 2s2 2p-2 2p4"
            @test c == rc"1s2c 2s2 2p-2 2p4"
            close!(c)
            @test sort(c) == rc"[Ne]"s
        end
        @test_throws ArgumentError close(c"1s")
        @test_throws ArgumentError close(rc"1s2 2s")
        @test_throws ArgumentError close!(c"1s")
        @test_throws ArgumentError close!(rc"1s2 2s")

        # Tests for #19
        @test c"10s2" == Configuration([o"10s"], [2], [:open])
        @test c"9999l32" == Configuration([o"9999l"], [32], [:open])

        @testset "Hashing" begin
            let c1a = c"1s 2s", c1b = c"1s 2s", c2 = c"1s 2p"
                @test c1a == c1b
                @test isequal(c1a, c1b)
                @test c1a != c2
                @test !isequal(c1a, c2)

                @test hash(c1a) == hash(c1a)
                @test hash(c1a) == hash(c1b)

                # If hashing is not properly implemented, unique fails to detect all equal pairs
                @test unique([c1a, c1b]) == [c1a]
                @test unique([c1a, c1a]) == [c1a]
                @test unique([c1a, c1a, c1b]) == [c1a]
                @test length(unique([c1a, c2, c1a, c1b, c2])) == 2
            end
        end

        @testset "Serialization" begin
            cs = [sc"1s(0,α) 2p(1,β)", sc"[He]c 2p(1,β)", rsc"[He]c 2p(-1/2)", sc"1s₀α 2p₁β", rsc"[Xe]*",
                  c"1s2 2p6", c"[He]c 2p6", c"[He]i 2p6", c"1s 3d4", c"[Xe]*",
                  rc"1s2 2p6", rc"1s 2p- 2p3", rc"[He]c 2p6", rc"[He]i 2p6", rc"1s 3d4", rc"[Xe]*"]

            ncs = let io = IOBuffer()
                foreach(Base.Fix1(write, io), cs)

                seekstart(io)
                [read(io, Configuration)
                 for i in eachindex(cs)]
            end

            @test cs == ncs
        end
    end

    @testset "Number of electrons" begin
        @test num_electrons(c"1s") == 1
        @test num_electrons(c"[He]") == 2
        @test num_electrons(rc"[He]") == 2
        @test num_electrons(c"[Xe]") == 54
        @test num_electrons(rc"[Xe]") == 54
        @test num_electrons(peel(c"[Kr]c 5s2 5p6")) == 8

        @test num_electrons(c"[Xe]", o"1s") == 2
        @test num_electrons(rc"[Xe]", ro"1s") == 2
        @test num_electrons(c"[Ne] 3s 3p", o"3p") == 1
        @test num_electrons(rc"[Kr] Af-2 Bf5", ro"Bf") == 5
        @test num_electrons(c"[Kr]", ro"5s") == 0
        @test num_electrons(rc"[Rn]", ro"As") == 0
    end

    @testset "Empty configurations" begin
        ec = empty(c"1s2")
        rec = empty(rc"1s2")
        @test num_electrons(ec) == 0
        @test num_electrons(rec) == 0
    end

    @testset "Access subsets" begin
        @test core(c"[Kr]c 5s2") == c"[Kr]c"
        @test peel(c"[Kr]c 5s2") == c"5s2"
        @test core(c"[Kr]c 5s2c 5p6") == c"[Kr]c 5s2c"
        @test peel(c"[Kr]c 5s2c 5p6") == c"5p6"
        @test active(c"[Kr]c 5s2") == c"5s2"
        @test inactive(c"[Kr]c 5s2i") == c"5s2i"
        @test inactive(c"[Kr]c 5s2") == c""

        @test bound(c"[Kr] 5s2 5p5 ks") == c"[Kr] 5s2 5p5"
        @test continuum(c"[Kr] 5s2 5p5 ks") == c"ks"
        @test bound(c"[Kr] 5s2 5p4 ks ld") == c"[Kr] 5s2 5p4"
        @test continuum(c"[Kr] 5s2 5p4 ks ld") == c"ks ld"

        @test bound(rc"[Kr] 5s2 5p-2 5p3 ks") == rc"[Kr] 5s2 5p-2 5p3"
        @test continuum(rc"[Kr] 5s2 5p-2 5p3 ks") == rc"ks"
        @test bound(rc"[Kr] 5s2 5p-2 5p2 ks ld") == rc"[Kr] 5s2 5p-2 5p2"
        @test continuum(rc"[Kr] 5s2 5p-2 5p2 ks ld") == rc"ks ld"

        @test c"[Ne]"[1] == (o"1s",2,:closed)
        @test c"[Ne]"[1:2] == c"1s2c 2s2c"
        @test c"[Ne]"[end-1:end] == c"2s2c 2p6c"

        @test c"1s2 2p kp"[2:3] == c"2p kp"

        @test o"1s" ∈ c"[He]"
        @test ro"1s" ∈ rc"[He]"
    end

    @testset "Sorting" begin
        @test c"1s2" < c"1s 2s"
        @test c"1s2" < c"1s 2p"
        @test c"1s2" < c"ks2"
        # @test c"kp2" < c"kp kd" # Ideally this would be true, but too
        #                         # complicated to implement
    end

    @testset "Parity" begin
        @test iseven(parity(c"1s"))
        @test iseven(parity(c"1s2"))
        @test iseven(parity(c"1s2 2s"))
        @test iseven(parity(c"1s2 2s2"))
        @test iseven(parity(c"[He]c 2s"))
        @test iseven(parity(c"[He]c 2s2"))
        @test isodd(parity(c"[He]c 2s 2p"))
        @test isodd(parity(c"[He]c 2s2 2p"))
        @test iseven(parity(c"[He]c 2s 2p2"))
        @test iseven(parity(c"[He]c 2s2 2p2"))
        @test iseven(parity(c"[He]c 2s2 2p2 3d"))
        @test isodd(parity(c"[He]c 2s kp"))
    end

    @testset "Pretty printing" begin
        Xe⁺ = rc"[Kr]c 5s2 5p-2 5p3"
        map([c"1s" => "1s",
             c"1s2" => "1s²",
             c"1s2 2s2" => "1s² 2s²",
             c"1s2 2s2 2p" => "1s² 2s² 2p",
             c"1s2c 2s2 2p" => "[He]ᶜ 2s² 2p",
             c"[He]* 2s2" => "1s² 2s²",
             c"[He]c 2s2" => "[He]ᶜ 2s²",
             c"[He]i 2s2" => "1s²ⁱ 2s²",
             c"[Kr]c 5s2" => "[Kr]ᶜ 5s²",
             Xe⁺ => "[Kr]ᶜ 5s² 5p-² 5p³",
             core(Xe⁺) => "[Kr]ᶜ",
             peel(Xe⁺) => "5s² 5p-² 5p³",
             rc"[Kr] 5s2c 5p6"s => "[Kr]ᶜ 5s²ᶜ 5p-² 5p⁴",
             c"[Ne]"[end:end] => "2p⁶ᶜ",
             sort(rc"[Ne]"[end-1:end]) => "2p-²ᶜ 2p⁴ᶜ",
             c"5s2" => "5s²",
             rc"[Kr]*"s => "1s² 2s² 2p-² 2p⁴ 3s² 3p-² 3p⁴ 3d-⁴ 3d⁶ 4s² 4p-² 4p⁴",
             c"[Kr]c" =>"[Kr]ᶜ",
             c"1s2 kp" => "1s² kp",
             c"" => "∅"]) do (c,s)
                 @test "$(c)" == s
             end
    end

    @testset "Orbital substitutions" begin
        @test replace(c"1s2", o"1s"=>o"2p") == c"1s 2p"
        @test replace(c"1s2", o"1s"=>o"kp") == c"1s kp"
        @test replace(c"1s kp", o"kp"=>o"ld") == c"1s ld"
        @test_throws ArgumentError replace(c"1s2", o"2p"=>o"3p")
        @test_throws ArgumentError replace(c"1s2 2s", o"2s"=>o"1s")

        @test replace(c"1s2 2s", o"1s"=>o"2p") == c"1s 2p 2s"
        @test replace(c"1s2 2s", o"1s"=>o"2p", append=true) == c"1s 2s 2p"
        @test replace(c"1s2 2s"s, o"1s"=>o"2p") == c"1s 2s 2p"
    end

    @testset "Orbital removal" begin
        @test c"1s2" - o"1s" == c"1s"
        @test c"1s" - o"1s" == c""
        @test c"1s 2s" - o"1s" == c"2s"
        @test c"[Ne]" - o"2s" == c"[He] 2s 2p6c"
        @test -(c"1s2 2s", o"1s", 2) == c"2s"

        @test delete!(c"1s 2s", o"1s") == c"2s"
        @test delete!(c"1s2 2s", o"1s") == c"2s"
        @test delete!(c"[Ne]", o"2p") == c"1s2c 2s2c"
        @test delete!(rc"[Ne]", ro"2p-") == rc"1s2c 2s2c 2p4c"
    end

    @testset "Configuration additions" begin
        @test c"1s" + c"1s" == c"[He]*"
        @test c"1s" + c"2p" == c"1s 2p"
        @test c"1s" + c"ld" == c"1s ld"
        @test c"[Kr]*" + c"4d10" + c"5s2" + c"5p6" == c"[Xe]*"
        @test_throws ArgumentError c"[He]*" + c"1s"
        @test_throws ArgumentError c"1si" + c"1s"
    end

    @testset "Configuration juxtapositions" begin
        @test c"" ⊗ c"" == [c""]
        @test c"1s" ⊗ c"" == [c"1s"]
        @test c"" ⊗ c"1s" == [c"1s"]
        @test c"1s" ⊗ c"2s" == [c"1s 2s"]
        @test rc"1s" ⊗ [rc"2s", rc"2p-", rc"2p"] == [rc"1s 2s", rc"1s 2p-", rc"1s 2p"]
        @test [rc"1s", rc"2s"] ⊗ rc"2p-" == [rc"1s 2p-", rc"2s 2p-"]
        @test [rc"1s", rc"2s"] ⊗ [rc"2p-", rc"2p"] == [rc"1s 2p-", rc"1s 2p", rc"2s 2p-", rc"2s 2p"]
        @test [rc"1s 2s"] ⊗ [rc"2p-2", rc"2p4"] == [rc"1s 2s 2p-2", rc"1s 2s 2p4"]
        @test [rc"1s 2s"] ⊗ [rc"kp-2", rc"lp4"] == [rc"1s 2s kp-2", rc"1s 2s lp4"]
    end

    @testset "Non-relativistic orbitals" begin
        import AtomicLevels: rconfigurations_from_orbital
        @test rconfigurations_from_orbital(1, 0, 1) == [rc"1s"]
        @test rconfigurations_from_orbital(1, 0, 2) == [rc"1s2"]
        @test rconfigurations_from_orbital(2, 0, 2) == [rc"2s2"]

        @test rconfigurations_from_orbital(2, 1, 1) == [rc"2p-", rc"2p"]
        @test rconfigurations_from_orbital(2, 1, 2) == [rc"2p-2", rc"2p- 2p", rc"2p2"]
        @test rconfigurations_from_orbital(2, 1, 3) == [rc"2p-2 2p1", rc"2p- 2p2", rc"2p3"]
        @test rconfigurations_from_orbital(2, 1, 4) == [rc"2p-2 2p2", rc"2p- 2p3", rc"2p4"]
        @test rconfigurations_from_orbital(2, 1, 5) == [rc"2p-2 2p3", rc"2p- 2p4"]
        @test rconfigurations_from_orbital(2, 1, 6) == [rc"2p-2 2p4"]

        @test rconfigurations_from_orbital(4, 2, 10) == [rc"4d-4 4d6"]
        @test rconfigurations_from_orbital(4, 2, 5) == [rc"4d-4 4d1", rc"4d-3 4d2", rc"4d-2 4d3", rc"4d-1 4d4", rc"4d5"]

        @test rconfigurations_from_orbital(:k, 2, 5) == [rc"kd-4 kd1", rc"kd-3 kd2", rc"kd-2 kd3", rc"kd-1 kd4", rc"kd5"]

        @test_throws ArgumentError rconfigurations_from_orbital(1, 2, 1)
        @test_throws ArgumentError rconfigurations_from_orbital(1, 0, 3)

        @test rconfigurations_from_orbital(o"1s", 2) == [rc"1s2"]
        @test rconfigurations_from_orbital(o"2p", 2) == [rc"2p-2", rc"2p- 2p", rc"2p2"]
        @test_throws ArgumentError rconfigurations_from_orbital(o"3d", 20)

        @test rcs"1s" == [rc"1s"]
        @test rcs"2s2" == [rc"2s2"]
        @test rcs"2p4" == [rc"2p-2 2p2", rc"2p- 2p3", rc"2p4"]
        @test rcs"4d5" == [rc"4d-4 4d1", rc"4d-3 4d2", rc"4d-2 4d3", rc"4d-1 4d4", rc"4d5"]

        @test rc"1s2" ⊗ rcs"2p" == [rc"1s2 2p-", rc"1s2 2p"]
        @test rc"1s2" ⊗ rcs"kp" == [rc"1s2 kp-", rc"1s2 kp"]
    end

    @testset "Spin-orbitals" begin
        up, down = half(1),-half(1)

        @test_throws ArgumentError Configuration(spin_orbitals(o"1s"), [2,1])
        @test_throws ArgumentError Configuration(spin_orbitals(o"1s"), [1,2])

        @test spin_configurations(c"1s2") ==
            [Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"1s",0,down)], [1, 1])]
        @test spin_configurations(c"1s2") isa Vector{Configuration{SpinOrbital{Orbital{Int},Tuple{Int,HalfInt}}}}
        @test spin_configurations(c"ks2") isa Vector{Configuration{SpinOrbital{Orbital{Symbol},Tuple{Int,HalfInt}}}}
        @test spin_configurations(c"1s ks") isa Vector{Configuration{SpinOrbital{<:Orbital,Tuple{Int,HalfInt}}}}
        @test spin_configurations(sc"1s₀α") == [sc"1s₀α"]

        @test scs"1s 2p" == spin_configurations(c"1s 2p") ==
            [Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"2p",-1,up)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"2p",-1,up)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"2p",-1,down)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"2p",-1,down)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"2p",0,up)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"2p",0,up)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"2p",0,down)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"2p",0,down)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"2p",1,up)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"2p",1,up)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"2p",1,down)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"2p",1,down)], [1, 1])]

        @test rscs"1s 2p" == spin_configurations(rc"1s 2p") ==
            [Configuration([rso"1s(-1/2)", rso"2p(-3/2)"], [1, 1])
             Configuration([rso"1s(1/2)", rso"2p(-3/2)"], [1, 1])
             Configuration([rso"1s(-1/2)", rso"2p(-1/2)"], [1, 1])
             Configuration([rso"1s(1/2)", rso"2p(-1/2)"], [1, 1])
             Configuration([rso"1s(-1/2)", rso"2p(1/2)"], [1, 1])
             Configuration([rso"1s(1/2)", rso"2p(1/2)"], [1, 1])
             Configuration([rso"1s(-1/2)", rso"2p(3/2)"], [1, 1])
             Configuration([rso"1s(1/2)", rso"2p(3/2)"], [1, 1])]

        @testset "Excited spin configurations" begin
            function excited_spin_configurations(cfg, orbitals)
                ec = excited_configurations(cfg, orbitals...,
                                            max_excitations=:singles, keep_parity=false)
                spin_configurations(ec)
            end

            a = excited_spin_configurations(c"1s2", os"2[s-p]")
            @test a isa Vector{Configuration{SpinOrbital{Orbital{Int},Tuple{Int,HalfInt}}}}
            b = excited_spin_configurations(c"1s2", os"k[s-p]")
            @test b isa Vector{<:Configuration{<:SpinOrbital{<:Orbital,Tuple{Int,HalfInt}}}}
            c = excited_spin_configurations(c"ks2", os"l[s-p]")
            @test c isa Vector{Configuration{SpinOrbital{Orbital{Symbol},Tuple{Int,HalfInt}}}}

            for (u,v) in [(a,b), (a,c), (b,c)]
                @test vcat(u,v) isa Vector{<:Configuration{<:SpinOrbital{<:Orbital,Tuple{Int,HalfInt}}}}
            end

            @test b ==
                [Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"1s",0,down)], [1, 1]),

                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"ks",0,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"ks",0,down)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"kp",-1,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"kp",-1,down)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"kp",0,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"kp",0,down)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"kp",1,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"kp",1,down)], [1, 1]),

                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"ks",0,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"ks",0,down)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"kp",-1,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"kp",-1,down)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"kp",0,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"kp",0,down)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"kp",1,up)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,down), SpinOrbital(o"kp",1,down)], [1, 1])]

            d=excited_spin_configurations(rc"1s2", ros"2[s-p]")
            @test d isa Vector{Configuration{SpinOrbital{RelativisticOrbital{Int},Tuple{HalfInt}}}}
            e=excited_spin_configurations(rc"1s2", ros"k[s-p]")
            @test e isa Vector{<:Configuration{<:SpinOrbital{<:RelativisticOrbital,Tuple{HalfInt}}}}
            f=excited_spin_configurations(rc"ks2", ros"l[s-p]")
            @test f isa Vector{Configuration{SpinOrbital{RelativisticOrbital{Symbol},Tuple{HalfInt}}}}

            for (u,v) in [(d,e), (d,f), (e,f)]
                @test vcat(u,v) isa Vector{<:Configuration{<:SpinOrbital{<:RelativisticOrbital,Tuple{HalfInt}}}}
            end

            for (u,v) in Iterators.product([a,b,c], [d,e,f])
                @test vcat(u,v) isa Vector{<:Configuration{<:SpinOrbital}}
            end
        end

        @test bound(Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"ks",0,up)], [1, 1])) ==
            Configuration([SpinOrbital(o"1s",0,up),], [1,])

        @test continuum(Configuration([SpinOrbital(o"1s",0,up), SpinOrbital(o"ks",0,up)], [1, 1])) ==
            Configuration([SpinOrbital(o"ks",0,up),], [1,])

        @test all([o ∈ spin_configurations(c"1s2")[1]
                   for o in spin_orbitals(o"1s")])

        @test map(spin_configurations(c"[Kr] 5s2 5p6")[1]) do (orb,occ,state)
            orb.orb.n < 5 && state == :closed || orb.orb.n == 5 && state == :open
        end |> all

        @test_broken string.(spin_configurations(c"2s2 2p"s)) ==
            ["2s² 2p₋₁α",
             "2s² 2p₋₁β",
             "2s² 2p₀α",
             "2s² 2p₀β",
             "2s² 2p₁α",
             "2s² 2p₁β"]

        @test_broken string.(spin_configurations(c"[Kr] 5s2 5p5 ks"s)) ==
            ["[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁α ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁β ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁α ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁β ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁² ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁² ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁² ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁² ks₀β"]

        @test_broken string.(spin_configurations(c"[Kr] 5s2 5p4"s)) ==
            ["[Kr]ᶜ 5s² 5p₋₁² 5p₀²",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀α 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀β 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀α 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀β 5p₁²",
             "[Kr]ᶜ 5s² 5p₀² 5p₁²"]

        @test substitutions(spin_configurations(c"1s2")[1],
                            spin_configurations(c"1s2")[1]) == []

        @test substitutions(spin_configurations(c"1s2")[1],
                            spin_configurations(c"1s ks")[2]) ==
                                [SpinOrbital(o"1s",0,up)=>SpinOrbital(o"ks",0,up)]

        @test spin_configurations([c"1s2", c"2p6"]) == [sc"1s₀α 1s₀β", sc"2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁α 2p₁β"]
        @test spin_configurations([rc"1s2", rc"2p6"]) == [rsc"1s(-1/2) 1s(1/2)", rsc"2p(-3/2) 2p(-1/2) 2p(1/2) 2p(3/2) 2p-(-1/2) 2p-(1/2)"]
    end

    @testset "Configuration transformations" begin
        @testset "Relativistic -> non-relativistic" begin
            @test nonrelconfiguration(rc"1s2 2p-2 2s 2p2 3s2 3p-"s) == c"1s2 2s 2p4 3s2 3p"s
            @test nonrelconfiguration(rc"1s2 ks2"s) == c"1s2 ks2"s
            @test nonrelconfiguration(rc"kp-2 kp4 lp-2 lp"s) == c"kp6 lp3"s
        end
        @testset "Non-relativistic -> relativistic" begin
            @test relconfigurations(c"1s") == [rc"1s"]
            @test relconfigurations(c"1s2") == [rc"1s2"]
            @test relconfigurations(c"[He] 2p") == [rc"[He] 2p-", rc"[He] 2p"]
            @test relconfigurations(c"[Ne]"s) == [rc"[Ne]"s]
            @test relconfigurations([c"1s 2p", c"2s 2p"]) == [rc"1s 2p-", rc"1s 2p", rc"2s 2p-", rc"2s 2p"]
        end
    end

    @testset "Internal utilities" begin
        @test AtomicLevels.get_noble_core_name(c"1s2 2s2 2p6") === nothing
        @test AtomicLevels.get_noble_core_name(c"1s2c 2s2 2p6") === "He"
        @test AtomicLevels.get_noble_core_name(c"1s2c 2s2c 2p6c") === "Ne"
        @test AtomicLevels.get_noble_core_name(c"[He] 2s2 2p6c 3s2c 3p6c") === "He"
        for gas in ["Rn", "Xe", "Kr", "Ar", "Ne", "He"]
            c = parse(Configuration{Orbital}, "[$gas]")
            @test AtomicLevels.get_noble_core_name(c) === gas
            rc = parse(Configuration{RelativisticOrbital}, "[$gas]")
            @test AtomicLevels.get_noble_core_name(rc) === gas
        end
    end

    @testset "Unsorted configurations" begin
        @testset "Comparisons" begin
            # These configurations should be similar but not equal
            @test issimilar(c"1s 2s", c"1s 2si")
            @test issimilar(c"1s 2s", c"2s 1s")
            @test issimilar(c"1s 2s", c"2si 1s")

            @test c"1s 2s" != c"1s 2si"
            @test c"1s 2s" != c"2s 1s"
            @test c"1s 2s" == c"2s 1s"s
            @test c"1s 2s" != c"2si 1s"

            @test rc"1s 2s" != rc"1s 2si"
            @test rc"1s 2s" != rc"2s 1s"
            @test rc"1s 2s" == rc"2s 1s"s
            @test rc"1s 2s" != rc"2si 1s"
        end
        @testset "Substitutions" begin
            @testset "Spatial orbitals" begin
                # Sorted case
                a = Configuration([o"1s", o"2s"], [1,1], [:open, :open], sorted=true)
                @test replace(a, o"1s" => o"2p").orbitals == [o"2s", o"2p"]

                b = Configuration([o"1s", o"2s"], [1,1], [:open, :open], sorted=false)
                c = replace(b, o"1s" => o"2p")
                @test c.orbitals == [o"2p", o"2s"]
                @test "$c" == "2p 2s"

                d = c"[He]" + Configuration([o"2p"], [1], [:open], sorted=false) + c"2s"
                @test core(d) == c"[He]"
                @test peel(d) == c
            end

            @testset "Spin orbitals" begin
                orbs = [spin_orbitals(o) for o in [o"1s", o"2s"]]

                a = spin_configurations(Configuration(o"1s", 2, :open, sorted=true))[1]
                @test replace(a, orbs[1][1]=>orbs[2][1]).orbitals == [orbs[1][2], orbs[2][1]]

                b = spin_configurations(Configuration(o"1s", 2, :open, sorted=false))[1]
                c = replace(b, orbs[1][1]=>orbs[2][1])
                @test c.orbitals == [orbs[2][1], orbs[1][2]]
                @test "$c" == "2s₀α 1s₀β"
            end
        end
    end

    @testset "String representation" begin
        @test string(c"[Ar]c 3d10c 4s2c 4p6 4d10 5s2i 5p6") == "[Ar]ᶜ 3d¹⁰ᶜ 4s²ᶜ 4p⁶ 4d¹⁰ 5s²ⁱ 5p⁶"
        @test ascii(c"[Ar]c 3d10c 4s2c 4p6 4d10 5s2i 5p6") == "[Ar]c 3d10c 4s2c 4p6 4d10 5s2i 5p6"
        @test string(rc"[Ar]*") == "1s² 2s² 2p⁴ 2p-² 3s² 3p⁴ 3p-²"
        @test ascii(rc"[Ar]*") == "1s2 2s2 2p4 2p-2 3s2 3p4 3p-2"
        @test ascii(empty(c"1s2")) == "empty"
    end

    @testset "multiplicity" begin
        @test multiplicity(c"1s") == 2
        @test multiplicity(c"2s2") == 1
        @test multiplicity(c"2p") == 6
        @test multiplicity(c"2p5") == 6
        @test multiplicity(c"2p2") == 15
        @test multiplicity(rc"1s") == 2
        @test multiplicity(rc"2p-2") == 1
        @test multiplicity(rc"2p4") == 1
        @test multiplicity(rc"3d3") == 20
        @test multiplicity(c"1s 2s") == 4
        @test multiplicity(c"1s 2s 2p") == 24
        @test multiplicity(c"1s2 2s2 2p6") == 1
        @test multiplicity(rc"1s2 2s2 2p-2 2p4") == 1

        let c = c"1s 2s2 2p3"; @test multiplicity(c) == length(spin_configurations(c)) end
        let c = c"1s2 4f1 5g1"; @test multiplicity(c) == length(spin_configurations(c)) end
        let c = rc"1s 2s2 2p3 2p-"; @test multiplicity(c) == length(spin_configurations(c)) end
        let c = rc"1s2 2s2 2p-2"; @test multiplicity(c) == length(spin_configurations(c)) end
        let c = rc"ag-6 bf4"; @test multiplicity(c) == length(spin_configurations(c)) end

        for c in spin_configurations(c"1s 2s2 2p3")
            @test multiplicity(c) == 1
        end
        for c in spin_configurations(rc"1s2 2s2 2p- 2p2")
            @test multiplicity(c) == 1
        end
    end
end

ext: jl


text:
using KrylovKit
using Rimu
using Rimu.RMPI
using DistributedDictVectors
using LinearAlgebra
using Test

# Ignore all printing on ranks other than root.
if mpi_rank() != mpi_root
    redirect_stderr(devnull)
    redirect_stdout(devnull)
end

@testset "simple operations" begin
    dv1 = DVec(zip(1:100, -50.0:49.0))
    tv1 = TVec(pairs(dv1))
    dv2 = DVec(zip(1:200, -100.0:199.0))
    tv2 = TVec(pairs(dv2))

    @test_broken dv1 == tv1 && dv2 == tv2
    @test norm(dv1) ≈ norm(tv1)
    @test sum(dv2) ≈ sum(tv2)
    @test dv1 ⋅ dv2 ≈ tv1 ⋅ tv2

    @test norm(tv1 + tv2) == norm(dv1 + dv2)
    axpy!(0.5, tv1, tv2)
    axpy!(0.5, dv1, dv2)
    normalize!(tv2)
    normalize!(dv2)
    @test dot(tv2, dv2) ≈ dot(dv2, tv2) ≈ 1
end

@testset "multiply" begin
    add = BoseFS((0,0,5,0,0))
    H = HubbardMom1D(add)
    dv = DVec(add => 1.0)
    tv = TVec(add => 1.0)

    for i in 1:10
        dv = H * dv
        tv = H * tv
        @test norm(dv) ≈ norm(tv)
        normalize!(dv)
        normalize!(tv)
        @test dot(dv, tv) ≈ 1
        @test length(dv) == length(tv)
    end
end

@testset "eigsolve" begin
    add = BoseFS((0,0,15,0,0))
    H = HubbardMom1DEP(add)
    dv = DVec(add => 1.0)
    tv = TVec(add => 1.0)

    res_dv = eigsolve(H, dv, 1, :SR; issymmetric=true)[1][1]
    res_tv = eigsolve(H, tv, 1, :SR)[1][1]

    @test res_tv ≈ res_dv
end

@testset "fciqmc" begin
    add = FermiFS2C((0,0,1,0,0), (0,1,0,1,0))
    H = Transcorrelated1D(add)

    dv = DVec(add => 1.0; style=IsDynamicSemistochastic())
    tv = TVec(add => 1.0; style=IsDynamicSemistochastic())

    println("ping")
    df_dv = lomc!(H, dv; laststep=10000).df
    println("pong")
    df_tv = lomc!(H, tv; laststep=10000).df
    println("pang")

    μ_dv = shift_estimator(df_dv; skip=2000).mean
    μ_tv = shift_estimator(df_tv; skip=2000).mean

    @test μ_dv ≈ μ_tv rtol=0.01
end


text:
#######################################################################
# Diagonalize methods
#   (All but LinearAlgebraPackage `@require` some package to be loaded)
#######################################################################
abstract type AbstractDiagonalizeMethod end

struct DiagonalizeHelper{S<:AbstractDiagonalizeMethod,C}
    method::S
    codiag::C
    minprojection::Float64
end

DiagonalizeHelper(method, codiag; minprojection = 0.5) =
    DiagonalizeHelper(method, codiag, minprojection)

## Diagonalize methods ##

defaultmethod(h::Union{Hamiltonian,ParametricHamiltonian,AbstractMatrix}) = LinearAlgebraPackage()

checkloaded(package::Symbol) = isdefined(Main, package) ||
    throw(ArgumentError("Package $package not loaded, need to be `using $package`."))

## LinearAlgebra ##
struct LinearAlgebraPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

LinearAlgebraPackage(; kw...) = LinearAlgebraPackage(values(kw))

function diagonalize(matrix, method::LinearAlgebraPackage)
    ϵ, ψ = eigen!(matrix; (method.kw)...)
    return ϵ, ψ
end

similarmatrix(h, ::LinearAlgebraPackage) = similarmatrix(h, Matrix{blockeltype(h)})

## Arpack ##
struct ArpackPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

ArpackPackage(; kw...) = (checkloaded(:Arpack); ArpackPackage(values(kw)))

function diagonalize(matrix, method::ArpackPackage)
    ϵ, ψ = Main.Arpack.eigs(matrix; (method.kw)...)
    return ϵ, ψ
end

similarmatrix(h, ::ArpackPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})

## IterativeSolvers ##

struct KrylovKitPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

KrylovKitPackage(; kw...) = (checkloaded(:KrylovKit); KrylovKitPackage(values(kw)))

function diagonalize(matrix::AbstractMatrix{M}, method::KrylovKitPackage) where {M}
    ishermitian(matrix) || throw(ArgumentError("Only Hermitian matrices supported with KrylovKitPackage for the moment"))
    origin = get(method.kw, :origin, 0.0)
    howmany = get(method.kw, :howmany, 1)
    kw´ = Base.structdiff(method.kw, NamedTuple{(:origin,:howmany)}) # Remove origin option
    lmap = shiftandinvert(matrix, origin)
    T = eltypevec(matrix)
    n = size(matrix, 2)
    x0 = rand(T, n)
    ϵ, ψ, _ = Main.KrylovKit.eigsolve(x -> lmap * x, x0, howmany, :LM; kw´...)

    ϵ´ = invertandshift(ϵ, origin)
    resize!(ϵ´, howmany)

    dimh = size(matrix, 2)
    ψ´ = Matrix{T}(undef, dimh, howmany)
    for i in 1:howmany
        copyslice!(ψ´, CartesianIndices((1:dimh, i:i)), ψ[i], CartesianIndices((1:dimh,)))
    end

    return ϵ´, ψ´
end

similarmatrix(h, ::KrylovKitPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})

#######################################################################
# shift and invert methods
#######################################################################

function shiftandinvert(matrix::AbstractMatrix{Tv}, origin) where {Tv}
    cols, rows = size(matrix)
    # Shift away from real axis to avoid pivot point error in factorize
    matrix´ = diagshift!(parent(matrix), origin + im)
    F = factorize(matrix´)
    lmap = LinearMap{Tv}((x, y) -> ldiv!(x, F, y), cols, rows,
                         ismutating = true, ishermitian = false)
    return lmap
end

function diagshift!(matrix::AbstractMatrix, origin)
    matrix´ = parent(matrix)
    vals = nonzeros(matrix´)
    rowval = rowvals(matrix´)
    for col in 1:size(matrix, 2)
        found_diagonal = false
        for ptr in nzrange(matrix´, col)
            if col == rowval[ptr]
                found_diagonal = true
                vals[ptr] -= origin * I  # To respect non-scalar eltypes
                break
            end
        end
        found_diagonal || throw(error("Sparse work matrix must include the diagonal. Possible bug in `similarmatrix`."))
    end
    return matrix
end

function invertandshift(ϵ::Vector{T}, origin) where {T}
    ϵ´ = similar(ϵ, real(T))
    ϵ´ .= real(inv.(ϵ) .+ (origin + im))  # Caution: we assume a real spectrum
    return ϵ´
end

#######################################################################
# resolve_degeneracies
#######################################################################
# Tries to make states continuous at crossings. Here, ϵ needs to be sorted
function resolve_degeneracies!(ϵ, ψ, ϕs, matrix, codiag)
    issorted(ϵ, by = real) || sorteigs!(codiag.perm, ϵ, ψ)
    hasapproxruns(ϵ, codiag.degtol) || return ϵ, ψ
    ranges, ranges´ = codiag.rangesA, codiag.rangesB
    resize!(ranges, 0)
    pushapproxruns!(ranges, ϵ, 0, codiag.degtol) # 0 is an offset
    for n in codiag.matrixindices
        v = codiag.cmatrixf(ϕs, n)
        resize!(ranges´, 0)
        for (i, r) in enumerate(ranges)
            subspace = view(ψ, :, r)
            vsubspace = subspace' * v * subspace
            veigen = eigen!(Hermitian(vsubspace))
            if hasapproxruns(veigen.values, codiag.degtol)
                roffset = minimum(r) - 1 # Range offset within the ϵ vector
                pushapproxruns!(ranges´, veigen.values, roffset, codiag.degtol)
            end
            subspace .= subspace * veigen.vectors
        end
        ranges, ranges´ = ranges´, ranges
        isempty(ranges) && break
    end
    return ψ
end

# Could perhaps be better/faster using a generalized CoSort
function sorteigs!(perm, ϵ::Vector, ψ::Matrix)
    resize!(perm, length(ϵ))
    p = sortperm!(perm, ϵ, by = real)
    # permute!(ϵ, p)
    sort!(ϵ, by = real)
    Base.permutecols!!(ψ, p)
    return ϵ, ψ
end

#######################################################################
# Codiagonalizer
#######################################################################

## Codiagonalizer
## Uses velocity operators along different directions. If not enough, use finite differences
## along mesh directions
struct Codiagonalizer{T,F<:Function}
    cmatrixf::F
    matrixindices::UnitRange{Int}
    degtol::T
    rangesA::Vector{UnitRange{Int}} # Prealloc buffer for degeneray ranges
    rangesB::Vector{UnitRange{Int}} # Prealloc buffer for degeneray ranges
    perm::Vector{Int}               # Prealloc for sortperm!
end

function codiagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, mesh, cut; kw...)
    veldirs = velocitydirections(parent(h); kw...)
    meshdirs = meshdirections(mesh; kw...)
    nv = length(veldirs)
    nm = length(meshdirs)
    matrixindices = 1:(nv + nm + 1)
    degtol = sqrt(eps(real(blockeltype(h))))
    delta = meshdelta(mesh)
    delta = iszero(delta) ? degtol : delta
    cmatrixf(meshϕs, n) =
        if n <= nv
            bloch!(matrix, h, applycut(cut, meshϕs), dn -> im * veldirs[n]' * dn)
        elseif n - nv <= nm # resort to finite differences
            bloch!(matrix, h, applycut(cut, meshϕs + delta * meshdirs[n - nv]))
        else # use a fixed random matrix
            randomfill!(matrix)
        end
    return Codiagonalizer(cmatrixf, matrixindices, degtol, UnitRange{Int}[], UnitRange{Int}[], Int[])
end

function codiagonalizer(matrixf::Function, matrix::AbstractMatrix{T}, mesh; kw...) where {T}
    meshdirs = meshdirections(mesh; kw...)
    nm = length(meshdirs)
    matrixindices = 1:(nm + 1)
    degtol = sqrt(eps(real(eltype(T))))
    delta = meshdelta(mesh)
    delta = iszero(delta) ? degtol : delta
    cmatrixf(meshϕs, n) =
        if n <= nm # finite differences
            matrixf(meshϕs + delta * meshdirs[n])
        else # use a fixed random matrix
            randomfill!(matrix)
        end
    return Codiagonalizer(cmatrixf, matrixindices, degtol, UnitRange{Int}[], UnitRange{Int}[], Int[])
end

velocitydirections(::Hamiltonian{LA,L}; kw...) where {LA,L} = _directions(Val(L); kw...)
meshdirections(::Mesh{L}; kw...) where {L} = _directions(Val(L); kw...)

function _directions(::Val{L}; direlements = 0:1, onlypositive = true) where {L}
    directions = vec(SVector{L,Int}.(Iterators.product(ntuple(_ -> direlements, Val(L))...)))
    onlypositive && filter!(ispositive, directions)
    unique!(normalize, directions)
    sort!(directions, by = norm, rev = false)
    return directions
end

meshdelta(mesh::Mesh{<:Any,T}) where {T} = T(0.1) * first(minmax_edge_length(mesh))

function randomfill!(matrix::SparseMatrixCSC, seed = 1234)
    Random.seed!(seed)
    rand!(nonzeros(matrix))  ## CAREFUL: this will be non-hermitian
    return matrix
end

function randomfill!(matrix::AbstractArray{T}, seed = 1234) where {T}
    Random.seed!(seed)
    fill!(matrix, zero(T))
    for i in 1:minimum(size(matrix))
        r = rand(T)
        @inbounds matrix[i, i] = r + r'
    end
    return matrix
end


text:
@rule NormalMixture{N}(:out, Marginalisation) (q_switch::Any, q_m::NTuple{N,  UnivariateNormalDistributionsFamily}, q_p::NTuple{N, GammaDistributionsFamily }) where { N } = begin
    πs = probvec(q_switch)
    return NormalMeanPrecision(sum(πs .* mean.(q_m)), sum(πs .* mean.(q_p)))
end

@rule NormalMixture{N}(:out, Marginalisation) (q_switch::Any, q_m::NTuple{N,  MultivariateNormalDistributionsFamily}, q_p::NTuple{N, Wishart }) where { N } = begin
    πs = probvec(q_switch)
    d  = ndims(first(q_m))
    w  = mapreduce(x -> x[1] * mean(x[2]), +, zip(πs, q_p))
    xi = mapreduce(x -> x[1] * mean(x[2]) * mean(x[3]), +, zip(πs, q_p, q_m))
    return MvNormalWeightedMeanPrecision(xi, w)
end


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

pointer_to(x::Union{CppValue, CppRef}) =
    icxx""" &$x; """

as_pointer(x::CppPtr) = x
as_pointer(x::Union{CppValue, CppRef}) = pointer_to(x)

isnullptr(ptr::CppPtr) = Ptr{Void}(ptr) == Ptr{Void}(0)

delete_if_not_nullptr(ptr::CppPtr) = begin
    if ! isnullptr(ptr)
        icxx""" delete $ptr; """
    end
end

auto_deref(v::CppPtr) = icxx"*$(v);"::CppRef
auto_deref(v::CppValue) = v
auto_deref(v::CppRef) = v


text:
function diag_full(A, X0; kwargs...)
    Neig = size(X0, 2)
    Afull = Hermitian(Array(A))
    E = eigen(Afull)
    X = E.vectors[:, 1:Neig]
    λ = E.values[1:Neig]
    (λ=λ, X=X,
     residual_norms=zeros(Neig),
     iterations=0,
     converged=true,
     n_matvec=0)
end


text:
export inplaceadd!

@inline function inplaceadd!{AT, IT, N}(arr::AbstractArray{AT, N}, v, x::Vararg{IT, N})
	arr[x...] += v
end

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# clipboard copy and paste

if Sys.isapple()
    function clipboard(x)
        pbcopy_cmd = `pbcopy`

        # OSX shells, especially when run within `tmux` or `screen`, can be
        # disconnected from the global shell namespace, which causes problems
        # with clipboards.  Luckily, the `reattach-to-user-namespace` utility
        # dodges these issues quite nicely, so we automatically utilize it if
        # it is installed.
        if Sys.which("reattach-to-user-namespace") !== nothing
            pbcopy_cmd = `reattach-to-user-namespace pbcopy`
        end

        open(pipeline(pbcopy_cmd, stderr=stderr), "w") do io
            print(io, x)
        end
        nothing
    end
    function clipboard()
        pbpaste_cmd = `pbpaste`

        # See above comment in `clipboard(x)`
        if Sys.which("reattach-to-user-namespace") !== nothing
            pbpaste_cmd = `reattach-to-user-namespace pbpaste`
        end
        return read(pbpaste_cmd, String)
    end

elseif Sys.islinux() || Sys.KERNEL === :FreeBSD
    _clipboardcmd = nothing
    const _clipboard_copy = Dict(
            :xsel  => Sys.islinux() ?
                `xsel --input --clipboard` :
                `xsel -c`,
            :xclip => `xclip -silent -in -selection clipboard`,
        )
    const _clipboard_paste = Dict(
            :xsel  => Sys.islinux() ?
                `xsel --nodetach --output --clipboard` :
                `xsel -p`,
            :xclip => `xclip -quiet -out -selection clipboard`,
        )
    function clipboardcmd()
        global _clipboardcmd
        _clipboardcmd !== nothing && return _clipboardcmd
        for cmd in (:xclip, :xsel)
            success(pipeline(`which $cmd`, devnull)) && return _clipboardcmd = cmd
        end
        pkgs = @static if Sys.KERNEL === :FreeBSD
            "x11/xsel or x11/xclip"
        else
            "xsel or xclip"
        end
        error("no clipboard command found, please install $pkgs")
    end
    function clipboard(x)
        c = clipboardcmd()
        cmd = _clipboard_copy[c]
        open(pipeline(cmd, stderr=stderr), "w") do io
            print(io, x)
        end
        nothing
    end
    function clipboard()
        c = clipboardcmd()
        cmd = _clipboard_paste[c]
        return read(pipeline(cmd, stderr=stderr), String)
    end

elseif Sys.iswindows()
    function clipboard(x::AbstractString)
        if Base.containsnul(x)
            throw(ArgumentError("Windows clipboard strings cannot contain NUL character"))
        end
        x_u16 = Base.cwstring(x)
        pdata = Ptr{UInt16}(C_NULL)
        function cleanup(cause)
            errno = cause == :success ? UInt32(0) : Libc.GetLastError()
            if cause !== :OpenClipboard
                if cause !== :success && pdata != C_NULL
                    ccall((:GlobalFree, "kernel32"), stdcall, Cint, (Ptr{UInt16},), pdata)
                end
                ccall((:CloseClipboard, "user32"), stdcall, Cint, ()) == 0 && Base.windowserror(:CloseClipboard) # this should never fail
            end
            cause == :success || Base.windowserror(cause, errno)
            nothing
        end
        ccall((:OpenClipboard, "user32"), stdcall, Cint, (Ptr{Cvoid},), C_NULL) == 0 && return Base.windowserror(:OpenClipboard)
        ccall((:EmptyClipboard, "user32"), stdcall, Cint, ()) == 0 && return cleanup(:EmptyClipboard)
        # copy data to locked, allocated space
        pdata = ccall((:GlobalAlloc, "kernel32"), stdcall, Ptr{UInt16}, (Cuint, Csize_t), 2 #=GMEM_MOVEABLE=#, sizeof(x_u16))
        pdata == C_NULL && return cleanup(:GlobalAlloc)
        plock = ccall((:GlobalLock, "kernel32"), stdcall, Ptr{UInt16}, (Ptr{UInt16},), pdata)
        plock == C_NULL && return cleanup(:GlobalLock)
        ccall(:memcpy, Ptr{UInt16}, (Ptr{UInt16}, Ptr{UInt16}, Csize_t), plock, x_u16, sizeof(x_u16))
        unlock = ccall((:GlobalUnlock, "kernel32"), stdcall, Cint, (Ptr{UInt16},), pdata)
        (unlock == 0 && Libc.GetLastError() == 0) || return cleanup(:GlobalUnlock) # this should never fail
        pset = ccall((:SetClipboardData, "user32"), stdcall, Ptr{UInt16}, (Cuint, Ptr{UInt16}), 13, pdata)
        pdata != pset && return cleanup(:SetClipboardData)
        cleanup(:success)
    end
    clipboard(x) = clipboard(sprint(print, x)::String)
    function clipboard()
        function cleanup(cause)
            errno = cause == :success ? UInt32(0) : Libc.GetLastError()
            if cause !== :OpenClipboard
                ccall((:CloseClipboard, "user32"), stdcall, Cint, ()) == 0 && Base.windowserror(:CloseClipboard) # this should never fail
            end
            if cause !== :success && (cause !== :GetClipboardData || errno != 0)
                Base.windowserror(cause, errno)
            end
            ""
        end
        ccall((:OpenClipboard, "user32"), stdcall, Cint, (Ptr{Cvoid},), C_NULL) == 0 && return Base.windowserror(:OpenClipboard)
        ccall(:SetLastError, stdcall, Cvoid, (UInt32,), 0) # allow distinguishing if the clipboard simply didn't have text
        pdata = ccall((:GetClipboardData, "user32"), stdcall, Ptr{UInt16}, (Cuint,), 13)
        pdata == C_NULL && return cleanup(:GetClipboardData)
        plock = ccall((:GlobalLock, "kernel32"), stdcall, Ptr{UInt16}, (Ptr{UInt16},), pdata)
        plock == C_NULL && return cleanup(:GlobalLock)
        s = try
            # find NUL terminator (0x0000 16-bit code unit)
            len = 0
            while unsafe_load(plock, len + 1) != 0
                len += 1
            end
            # get Vector{UInt16}, transcode data to UTF-8, make a String of it
            transcode(String, unsafe_wrap(Array, plock, len))
        finally
           unlock = ccall((:GlobalUnlock, "kernel32"), stdcall, Cint, (Ptr{UInt16},), pdata)
           (unlock != 0 || Libc.GetLastError() == 0) || return cleanup(:GlobalUnlock) # this should never fail
           cleanup(:success)
        end
        return s
    end

else
    clipboard(x="") = error("`clipboard` function not implemented for $(Sys.KERNEL)")
end

"""
    clipboard(x)

Send a printed form of `x` to the operating system clipboard ("copy").
"""
clipboard(x)

"""
    clipboard() -> AbstractString

Return a string with the contents of the operating system clipboard ("paste").
"""
clipboard()

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using FileWatching.Pidfile

using Test

using Base.Filesystem: File
using FileWatching.Pidfile: iswindows,
    write_pidfile, parse_pidfile,
    isvalidpid, stale_pidfile,
    tryopen_exclusive, open_exclusive

# helper utilities
struct MemoryFile <: Base.AbstractPipe
    io::IOBuffer
    mtime::Float64
end
Base.pipe_reader(io::MemoryFile) = io.io
Base.Filesystem.mtime(io::MemoryFile) = io.mtime

# set the process umask so we can test the behavior of
# open mask without interference from parent's state
# and create a test environment temp directory
umask(new_mask) = ccall((@static iswindows() ? :_umask : :umask), Cint, (Cint,), new_mask)

# TODO: Use targeted @test_log tests instead of suppressing all logs to hide the expected warnings
Base.CoreLogging.with_logger(Base.CoreLogging.NullLogger()) do

@testset "Pidfile.jl" begin
old_umask = umask(0o002)
try
    mktempdir() do dir
        cd(dir) do

# now start tests definitions:

@testset "validpid" begin
    mypid = getpid() % Cuint
    @test isvalidpid(gethostname(), mypid)
    @test isvalidpid("", mypid)
    @test !isvalidpid("", 0 % Cuint)
    @test isvalidpid("NOT" * gethostname(), mypid)
    @test isvalidpid("NOT" * gethostname(), 0 % Cuint)
    @test isvalidpid("NOT" * gethostname(), -1 % Cuint)
    if !iswindows()
        @test isvalidpid("", 1 % Cuint)
        @test !isvalidpid("", -1 % Cuint)
        @test !isvalidpid("", -mypid)
    end
end

@testset "write_pidfile" begin
    buf = IOBuffer()
    pid, host, age = 0, "", 123
    pid2, host2, age2 = parse_pidfile(MemoryFile(seekstart(buf), time() - age))
    @test pid == pid2
    @test host == host2
    @test age ≈ age2 atol=5

    host = " host\r\n"
    write(buf, "-1 $host")
    pid2, host2, age2 = parse_pidfile(MemoryFile(seekstart(buf), time() - age))
    @test pid == pid2
    @test host == host2
    @test age ≈ age2 atol=5
    truncate(seekstart(buf), 0)

    pid, host = getpid(), gethostname()
    write_pidfile(buf, pid)
    @test read(seekstart(buf), String) == "$pid $host"
    pid2, host2, age2 = parse_pidfile(MemoryFile(seekstart(buf), time() - age))
    @test pid == pid2
    @test host == host2
    @test age ≈ age2 atol=5
    truncate(seekstart(buf), 0)

    @testset "parse_pidfile" begin
        age = 0
        @test parse_pidfile("nonexist") === (Cuint(0), "", 0.0)
        open(io -> write_pidfile(io, pid), "pidfile", "w")
        pid2, host2, age2 = parse_pidfile("pidfile")
        @test pid == pid2
        @test host == host2
        @test age ≈ age2 atol=10
        rm("pidfile")
    end
end

@assert !ispath("pidfile")
@testset "open_exclusive" begin
    f = open_exclusive("pidfile")::File
    try
        # check that f is open and read-writable
        @test isfile("pidfile")
        @test filemode("pidfile") & 0o777 == 0o444
        @test filemode(f) & 0o777 == 0o444
        @test filesize(f) == 0
        @test write(f, "a") == 1
        @test filesize(f) == 1
        @test read(seekstart(f), String) == "a"
        chmod("pidfile", 0o600)
        @test filemode(f) & 0o777 == (iswindows() ? 0o666 : 0o600)
    finally
        close(f)
    end

    # release the pidfile after a short delay
    deleted = false
    rmtask = @async begin
        sleep(3)
        rm("pidfile")
        deleted = true
    end
    isdefined(Base, :errormonitor) && Base.errormonitor(rmtask)
    @test isfile("pidfile")
    @test !deleted

    # open the pidfile again (should wait for it to disappear first)
    t = @elapsed f2 = open_exclusive(joinpath(dir, "pidfile"))::File
    try
        @test deleted
        @test isfile("pidfile")
        @test t > 2
        if t > 6
            println("INFO: watch_file optimization appears to have NOT succeeded")
        end
        @test filemode(f2) & 0o777 == 0o444
        @test filesize(f2) == 0
        @test write(f2, "bc") == 2
        @test read(seekstart(f2), String) == "bc"
        @test filesize(f2) == 2
    finally
        close(f2)
    end
    rm("pidfile")
    wait(rmtask)

    # now test with a long delay and other non-default options
    f = open_exclusive("pidfile", mode = 0o000)::File
    try
        @test filemode(f) & 0o777 == (iswindows() ? 0o444 : 0o000)
    finally
        close(f)
    end
    deleted = false
    rmtask = @async begin
        sleep(8)
        rm("pidfile")
        deleted = true
    end
    isdefined(Base, :errormonitor) && Base.errormonitor(rmtask)
    @test isfile("pidfile")
    @test !deleted
    # open the pidfile again (should wait for it to disappear first)
    t = @elapsed f2 = open_exclusive("pidfile", mode = 0o777, poll_interval = 1.0)::File
    try
        @test deleted
        @test isfile("pidfile")
        @test filemode(f2) & 0o777 == (iswindows() ? 0o666 : 0o775)
        @test write(f2, "def") == 3
        @test read(seekstart(f2), String) == "def"
        @test t > 7
    finally
        close(f2)
    end
    rm("pidfile")
    wait(rmtask)

    @testset "test for wait == false cases" begin
        f = open_exclusive("pidfile", wait=false)
        @test isfile("pidfile")
        close(f)
        rm("pidfile")

        f = open_exclusive("pidfile")::File
        deleted = false
        rmtask = @async begin
            sleep(2)
            @test Pidfile.tryrmopenfile("pidfile")
            deleted = true
        end
        isdefined(Base, :errormonitor) && Base.errormonitor(rmtask)

        t1 = time()
        @test_throws ErrorException open_exclusive("pidfile", wait=false)
        @test time()-t1 ≈ 0 atol=1

        sleep(1)
        @test !deleted

        t1 = time()
        @test_throws ErrorException open_exclusive("pidfile", wait=false)
        @test time()-t1 ≈ 0 atol=1

        wait(rmtask)
        @test deleted
        t = @elapsed f2 = open_exclusive("pidfile", wait=false)::File
        @test isfile("pidfile")
        @test t ≈ 0 atol=1
        close(f)
        close(f2)
        rm("pidfile")
    end
end

@assert !ispath("pidfile")
@testset "open_exclusive: break lock" begin
    # test for stale_age
    t = @elapsed f = open_exclusive("pidfile", poll_interval=3, stale_age=10)::File
    try
        write_pidfile(f, getpid())
    finally
        close(f)
    end
    @test t < 2
    t = @elapsed f = open_exclusive("pidfile", poll_interval=3, stale_age=1)::File
    close(f)
    @test 20 < t < 50
    rm("pidfile")

    t = @elapsed f = open_exclusive("pidfile", poll_interval=3, stale_age=10)::File
    close(f)
    @test t < 2
    t = @elapsed f = open_exclusive("pidfile", poll_interval=3, stale_age=10)::File
    close(f)
    @test 8 < t < 20
    rm("pidfile")
end

@testset "open_exclusive: other errors" begin
    error = @test_throws(Base.IOError, open_exclusive("nonexist/folder"))
    @test error.value.code == Base.UV_ENOENT

    error = @test_throws(Base.IOError, open_exclusive(""))
    @test error.value.code == Base.UV_ENOENT
end

@assert !ispath("pidfile")
@testset "mkpidlock" begin
    lockf = mkpidlock("pidfile")
    @test lockf.update === nothing
    waittask = @async begin
        sleep(3)
        cd(homedir()) do
            return close(lockf)
        end
    end
    isdefined(Base, :errormonitor) && Base.errormonitor(waittask)

    # mkpidlock with no waiting
    t = @elapsed @test_throws ErrorException mkpidlock("pidfile", wait=false)
    @test t ≈ 0 atol=1

    t = @elapsed lockf1 = mkpidlock(joinpath(dir, "pidfile"))
    @test t > 2
    @test istaskdone(waittask) && fetch(waittask)
    @test !close(lockf)
    finalize(lockf1)
    t = @elapsed lockf2 = mkpidlock("pidfile")
    @test t < 2
    @test !close(lockf1)

    # test manual breakage of the lock
    # is correctly handled
    @test Pidfile.tryrmopenfile("pidfile")
    t = @elapsed lockf3 = mkpidlock("pidfile")
    @test t < 2
    @test isopen(lockf2.fd)
    @test !close(lockf2)
    @test !isopen(lockf2.fd)
    @test isfile("pidfile")
    @test close(lockf3)
    @test !isfile("pidfile")

    # Just for coverage's sake, run a test with do-block syntax
    lock_times = Float64[]
    t_loop = @async begin
        for idx in 1:100
            t = @elapsed mkpidlock("do_block_pidfile") do
                # nothing
            end
            sleep(0.01)
            push!(lock_times, t)
        end
    end
    isdefined(Base, :errormonitor) && Base.errormonitor(t_loop)
    mkpidlock("do_block_pidfile") do
        sleep(3)
    end
    wait(t_loop)
    @test maximum(lock_times) > 2
    @test minimum(lock_times) < 1
end

@assert !ispath("pidfile")
@testset "mkpidlock update" begin
    lockf = mkpidlock("pidfile")
    @test lockf.update === nothing
    new = mtime(lockf.fd)
    @test new ≈ time() atol=1
    sleep(1)
    @test mtime(lockf.fd) == new
    touch(lockf)
    old, new = new, mtime(lockf.fd)
    @test new != old
    @test new ≈ time() atol=1
    close(lockf)

    lockf = mkpidlock("pidfile"; refresh=0.2)
    new = mtime(lockf.fd)
    @test new ≈ time() atol=1
    for i = 1:10
        sleep(0.5)
        old, new = new, mtime(lockf.fd)
        @test new != old
        @test new ≈ time() atol=1
    end
    @test isopen(lockf.update::Timer)
    close(lockf)
    @test !isopen(lockf.update::Timer)

    lockf = mkpidlock("pidfile"; stale_age=10)
    @test lockf.update isa Timer
    close(lockf.update) # simulate a finalizer running in an undefined order
    close(lockf)
end

@assert !ispath("pidfile")
@testset "mkpidlock for child" begin
    proc = open(`cat`, "w", devnull)
    lock = mkpidlock("pidfile", proc)
    @test isopen(lock.fd)
    @test isfile("pidfile")
    close(proc)
    @test success(proc)
    sleep(1) # give some time for the other task to finish releasing the lock resources
    @test !isopen(lock.fd)
    @test !isfile("pidfile")

    error = @test_throws Base.IOError mkpidlock("pidfile", proc)
    @test error.value.code == Base.UV_ESRCH
end

@assert !ispath("pidfile-2")
@testset "mkpidlock non-blocking stale lock break" begin
    # mkpidlock with no waiting
    lockf = mkpidlock("pidfile-2", wait=false)
    @test lockf.update === nothing

    sleep(1)
    t = @elapsed @test_throws ErrorException mkpidlock("pidfile-2", wait=false, stale_age=1, poll_interval=1, refresh=0)
    @test t ≈ 0 atol=1

    sleep(5)
    t = @elapsed (lockf2 = mkpidlock("pidfile-2", wait=false, stale_age=.1, poll_interval=1, refresh=0))
    @test t ≈ 0 atol=1
    close(lockf)
    close(lockf2)
end

end; end # cd(tempdir)
finally
    umask(old_umask)
end; end # testset

end # with_logger

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
function _expansion_point(a, b, c, X0::Interval, Δ)
    c₀ = mid(X0.dat)
    f_c₀ = c₀*(a*c₀ + b) + c
    x̃ = c₀ + Δ / 2 * f_c₀
    return x̃
end

function _linearize(a, b, c, x̃)
    μ = a * x̃
    α = 2 * μ + b
    β = -μ * x̃ + c
    return α, β
end

function reach_homog_QINT(; a, b, c, # right-hand side: f(x) = ax^2 + bx + c
                          X0,      # initial set
                          T,       # time span: [0, T]
                          Δ::N,    # step size for NL reach
                          δ::N,    # step size for the linear reach
                          θ,       # remainder expansion rate
                          maxiter  # maximum number of iterations
                          ) where {N}

    # total flowpipe
    RT = ReachSet{N,Interval{N, IA.Interval{N}}}
    VRT = Vector{RT}
    FT = Flowpipe{N, RT, VRT}
    Ftot = Vector{FT}()

    # initialization
    Δti = zeroI # defines the initial time
    waiting_list = WaitingList([Δti], [StateInLocation(X0, 1)])
    R̄err = Interval(-θ * Δ, θ * Δ)

    k = 1
    ksplit = 0
    while !isempty(waiting_list)

        (Δti, elem) = pop!(waiting_list)
        X0 = state(elem)

        # compute expansion point
        x̃ = _expansion_point(a, b, c, X0, Δ)

        # linearization
        α, β = _linearize(a, b, c, x̃)

        # solve linear reachability
        prob = @ivp(x' = α * x + β, x(0) ∈ X0)
        sol = post(INT(δ=δ), prob, IA.Interval(0.0, Δ), Δt0=Δti)

        # compute admissible linearization error
        kθ = 1 / (exp(α * Δ) - 1) * α * θ * Δ
        L̄ = Interval(-kθ, kθ)

        # estimate lagrange remainder TODO: optimize
        R̂lin = ConvexHullArray([set(R) for R in sol]) ⊕ R̄err
        R̂lin_int = overapproximate(R̂lin, Interval)
        L = a * (R̂lin_int.dat - x̃)^2 |> Interval

        # validate Lagrange remainder
        if L ⊆ L̄
            # if the inclusion check succeeded, add errors to each reach-set for final storage
            F = Flowpipe([ReachSet(Interval(set(R).dat + R̄err.dat), tspan(R)) for R in sol])
            push!(Ftot, F)

            # only add the final reach-set to the waiting list if we didn't hit the time horizon
            Rf = F[end]
            tf = tend(Rf)
            T <= tstart(Rf) && continue # TODO check
            elem = StateInLocation(set(Rf), 1)
            push!(waiting_list, tspan(Rf), elem)

        else
            ksplit += 1
            # if the inclusion didn't succeed, split initial condition
            # and append left and right interval segments to the waiting list
            X0l, X0r = split(X0, 2)
            iteml = StateInLocation(X0l, 1)
            itemr = StateInLocation(X0r, 1)

            push!(waiting_list, Δti, iteml)
            push!(waiting_list, Δti, itemr)
        end

        k += 1
        if k > maxiter
            error("maximum number of iterations maxiter=$maxiter reached, try increasing `maxiter`")
        end
    end # while

    return HybridFlowpipe(Ftot)
end


text:
"""
ci_pde_aa()
Duplicate part of the data for Figure 4.2 in the book.
"""
function ci_pde_aa()
n=63;
m=10;
pdata=pdeinit(n);
Vstore=zeros(n*n,3*m+3);
VstoreS=zeros(n*n,2*m+4);
aout=PDE_aa(n, m; Vstore=Vstore, pdata=pdata);
aoutS=PDE_aa(n, m; Vstore=VstoreS, pdata=pdata);
# Same results with low storage mode?
alphaS= reldiff(aout.stats.alphanorm,aoutS.stats.alphanorm)
condS= reldiff(aout.stats.condhist,aoutS.stats.condhist)
histS = norm(aoutS.history-aout.history,Inf)
pdeerrS=  condS + alphaS + histS
aout.idid || println("pde solver failed")
(aout.errcode == 0) || println("wrong error code in pde")
(pdeerrS < 1.e-8) || println("different stats   ",condS,"  ",alphaS,"  ",histS)
(length(aout.history)==21) || println("history length wrong")
aa_ok = aout.idid && (aout.errcode==0) && (length(aout.history)==21) && (pdeerrS < 1.e-8)
aa_ok && println("pde succeeds")
return aa_ok
end


text:
# ----------------------------------------------------------------
# Alternative: multh-threading
# ----------------------------------------------------------------

#=
$ export JULIA_NUM_THREADS=16
=#

function process_multithreaded(nfiles, nstates, nattr)
    result = zeros(nstates, nattr, nfiles)
    Threads.@threads for i in 1:nfiles
        read_val_file!(i, result)
    end
    return result
end

#=
julia> @time result = process_multithreaded(nfiles, nstates, nattr);
190.520584 seconds (2.95 M allocations: 44.268 GiB, 0.44% gc time)
=#


text:
include("utils.jl")
using Distances

Random.seed!(42)

@testset "Distances" begin
    M = SymmetricPositiveDefinite(3)
    A(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)]
    ptsF = [#
        [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1],
        [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1],
        A(π / 6) * [1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] * transpose(A(π / 6)),
    ]
    dist = ManifoldML.RiemannianDistance(M)
    @test evaluate(dist, reshape(ptsF[1], 9), reshape(ptsF[3], 9)) ≈ distance(M, ptsF[1], ptsF[3])

    point_matrix = reduce(hcat, map(a -> reshape(a, 9), ptsF))
    dists = pairwise(dist, point_matrix)
    for i in 1:3, j in 1:3
        @test dists[i, j] ≈ distance(M, ptsF[i], ptsF[j]) atol = 1e-15
    end
end

max_stars_repo_name: JuliaManifolds/ManifoldML

max_issues_repo_name: JuliaManifolds/ManifoldML

max_forks_repo_name: JuliaManifolds/ManifoldML


text:
"""
Module for performing a Probilistic PCA analysis

Author: Eric Ford
"""
module PPCA

using LinearAlgebra: dot, Symmetric
using Statistics: mean, std
using Flux                           # For optimization
using PDMats, BandedMatrices         # For regularizers
using Stheno #, TemporalGPs

abstract type AbstractPPCA end

struct PPCAGridded <: AbstractPPCA
  W::Array{Float64,2}
  b::Array{Float64,1}
  x::Array{Float64,2}
  lambda::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}

  function PPCAGridded(W::AA1, b::AA2, x::AA3, lambda_range::R) where { T<:Real, AA1<:AbstractArray{T,2}, AA2<:AbstractArray{T,1}, AA3<:AbstractArray{T,2}, R<:AbstractRange{T} }
     @assert size(W,1) == size(b,1)
     @assert size(W,2) == size(x,1)
     @assert length(lambda_range) == size(b,1)
     new(W, b, x, lambda_range )
  end

  function PPCAGridded(in::PPCAGridded)
     @assert size(in.W,1) == size(in.b,1)
     @assert size(in.W,2) == size(in.x,1)
     @assert length(in.lambda) == size(in.b,1)
    new(in.W,in.b,in.x,in.lambda)
  end

end

import Base.copy
copy(in::PPCAGridded) = PPCAGridded(in)

function PPCAGridded(in::Integer, out::Integer, nobs::Integer, lambda_range::R) where {R<:AbstractRange}
    @assert length(lambda_range) == out
    PPCAGridded(randn(out, in), randn(out), randn(in,nobs), lambda_range)
end

get_num_inputs(m::PPCAGridded) = size(m.W,2)
get_num_outputs(m::PPCAGridded) = size(m.W,1)
get_num_param(m::PPCAGridded) = get_num_outputs(m)*(get_num_inputs(m)+1)

Flux.@functor PPCAGridded (W,b,x)

# Overload call, so the object can be used as a function
function (m::PPCAGridded)(x::AbstractArray{T})::Array{T,2} where { T<:Real }
    m.W * m.x .+ m.b
end

struct SpectraSet1D
    lambda_obs::Array{Float64,2}
    lambda_ssb::Array{Float64,2}
    flux::Array{Float64,2}
    flux_var::Array{Float64,2}
end

get_num_outputs(s::SpectraSet1D) = size(s.flux,1)
get_num_obs(s::SpectraSet1D) = size(s.flux,2)

function chisq(x::AbstractArray{T,2}, m::PPCAGridded, ss::SpectraSet1D)::T where { T<:Real }
    yhat = m(x) # our predictions for fluxes's
    sum(((ss.flux.-yhat).^2 ./ss.flux_var))
end

chisq(m::PPCAGridded, ss::SpectraSet1D) = chisq(m.x,m,ss)

mutable struct Regularizer
    covar::AbstractPDMat{Float64}
    scale_input::Float64
    scale_output::Float64
    num_bands::Int64
    kernel::Function
end

const sqrt5 = sqrt(5)
matern52(x) = (one(x)+x*(sqrt5+x*sqrt5/3))*exp(-sqrt5*x)
expsq(x) = exp(-0.5*x^2)

function make_covar_banded(n::Integer, b::Integer, kernel::Function; lambda::Real = 1.0)
    @assert b == 5  # TODO: Generalize for arbitrary number of bands
    lags = 1:b
    covar_at_lag = kernel.(lags./(lambda*b))
    PDMat(Symmetric(BandedMatrix((0=>ones(n),1=>covar_at_lag[1]*ones(n-1),2=>covar_at_lag[2]*ones(n-2),3=>covar_at_lag[3]*ones(n-3),4=>covar_at_lag[4]*ones(n-4),5=>covar_at_lag[5]*ones(n-5)), (n,n), (b,b) )))
end

function update_lambda!(r::Regularizer, lambda::Real)
    r.scale_input = lambda
    r.covar = make_covar_banded(size(r.covar,1), r.num_bands, r.kernel, lambda=lambda)
    return r
end

function update_sigma!(r::Regularizer, sigma::Real)
    r.scale_output = sigma
    return r
end

function (r::Regularizer)(x::AbstractArray{T,1})::T where { T<:Real }
    r.scale_output*invquad(r.covar,x./r.scale_input)
end

function (r::Regularizer)(x::AbstractArray{T,2})::T where { T<:Real }
    sum(r(x[:,i]) for i in 1:size(x,2) )
end

function Matern52Regularizer(nobs::Integer, nband::Integer; lambda::Float64 = 1.0, sigma::Float64 = 1.0)
    b = make_covar_banded(nobs,nband, matern52, lambda=lambda)
    return Regularizer(b, lambda, sigma, nband, matern52)
end

function ExpSqRegularizer(nobs::Integer, nband::Integer; lambda::Float64 = 1.0, sigma::Float64 = 1.0)
    b = make_covar_banded(nobs,nband, expsq, lambda=lambda)
    return Regularizer(b, lambda, sigma, nband, expsq)
end

function regularization_L2norm(x::AbstractArray{T}) where {T<:Real}
    sum(abs2.(x))
end

function regularization_L1norm(x::AbstractArray{T}) where {T<:Real}
    sum(abs.(x))
end

function regularization_deriv(x::AbstractArray{T,2})  where { T<:Real }
    sum(abs2.(x[2:end].-x[1:end-1]))
end

function regularization_deriv2(x::AbstractArray{T,2})  where { T<:Real }
    sum(abs2.(-x[1:end-2].+2*x[2:end-1].-x[3:end]))
end

function regularization_scores(x::AbstractArray{T,2})  where { T<:Real }
    regularization_L2norm(x)
end

function regularization_mean(x::AbstractArray{T,1}; mu=zero(T) )  where { T<:Real }
    (gp_reg.scale_output == zero(T)) ? zero(T) : gp_reg(x.-mu)
end

function regularization_basis(x::AbstractArray{T,2})  where { T<:Real }
    (gp_reg.scale_output == zero(T)) ? zero(T) : gp_reg(x)
end

function regularization(x::AbstractArray{T,2}, m::PPCAGridded) #=, ss::SpectraSet1D)=#  where { T<:Real }
    regularization_scores(x) + regularization_basis(m.W) + regularization_mean(m.b, mu=1.0)
end

function regularization(m::PPCAGridded) #=, ss::SpectraSet1D) =# where { T<:Real }
    regularization_scores(m.x) + regularization_basis(m.W) + regularization_mean(m.b, mu=1.0)
end

function loss(x::AbstractArray{T,2}, m::M, ss::SpectraSet1D)  where { T<:Real, M<:AbstractPPCA }
    chisq(x,m,ss) + regularization(x,m)
end

function loss(m::M, ss::SpectraSet1D)  where { M<:AbstractPPCA }
    chisq(m.x,m,ss) + regularization(m.x,m)
end

# random initial parameters
num_outputs = 50
num_inputs = 3
num_obs = 50
snr_per_pixel = 50
snr_per_pixel = 200

lambda_min = 5000
lambda_max = 5025
#lambda_pad = (lambda_max-lambda_min)/num_outputs
#lambda_range_extra = range(lambda_min-lambda_pad,stop=lambda_max+lambda_pad,length=num_outputs+2)
lambda_range = range(lambda_min,stop=lambda_max,length=num_outputs)
# ppca = PPCAGridded(num_inputs,num_outputs, lambda_range)

# Create a PPCA model to create "true" data
ppca_mean = 1.01 .- 0.01*4.0.*((lambda_range.-mean(lambda_range))./(maximum(lambda_range)-minimum(lambda_range))).^2
ppca_basis = zeros(num_outputs,num_inputs+1)
ppca_basis[:,1] .= sin.(2*pi.*lambda_range/20.0)
ppca_basis[:,2] .= sin.(2*pi.*lambda_range/10.0)
ppca_basis[:,3] .= sin.(2*pi.*lambda_range/4.0)
ppca_basis[:,4] .= 0.0; # sin.(2*pi.*lambda_range/6.0)

x_true = randn(num_inputs+1,num_obs)        # Generate random weights for each basis vector at each epoch
ppca_true = PPCAGridded(ppca_basis, ppca_mean, x_true, lambda_range)
ytrue = ppca_true(x_true)
#ytrue = ppca_basis*x_true.+ppca_mean
yobs = ytrue .+ 1.0/snr_per_pixel * randn(num_outputs,num_obs)  # Add measurement noise
x = copy(x_true[1:num_inputs,1:num_obs])

lambda_ranges = reshape(repeat(collect(lambda_range),num_obs),num_outputs,num_obs)
doppler_factor = 1.001
spec = SpectraSet1D(lambda_ranges.*doppler_factor, lambda_ranges, yobs, ones(size(yobs))./snr_per_pixel )

@time chisq(ppca_true,spec)

using Plots

plot(ppca_true.lambda,ppca_true.W[:,1], color=:blue, label="W1", title="Basis vectors (true)")
#scatter!(ppca_true.lambda,ppca.W[:,1]./maximum(ppca.W[:,1]), ms=2, color=:blue)
plot!(ppca_true.lambda,ppca_true.W[:,2], color=:red, label="W2")
#scatter!(ppca_true.lambda,ppca.W[:,2]./maximum(ppca.W[:,2]), ms=2,color=:red)
plot!(ppca_true.lambda,ppca_true.W[:,3], color=:green, label="W3")
plot!(ppca_true.lambda,ppca_true.W[:,4], color=:cyan, label="W4")
#scatter!(ppca_true.lambda,ppca.W[:,3]./maximum(ppca.W[:,3]), ms=2, color=:green)
plot!(ppca_true.lambda,ppca_true.b, color=:black, label="b")

plot(lambda_range,ppca_true(x_true)[:,1:5],lw=2, label=:none, title="Sample Obs")
#scatter!(lambda_range,ppca(x_true)[:,1], label="Data (init)", color=:lightblue)
xlabel!("λ") # * L"\A")

# Choose Regularization method
gp_reg = Matern52Regularizer(num_outputs, 5, lambda=0.4, sigma=0.0)
#gp_reg = ExpSqRegularizer(num_outputs,5,lambda=0.25, sigma=0.0)
#gp_reg.covar[1:7,1:7]

ppca = PPCAGridded(num_inputs,num_outputs, num_obs, lambda_range) # Create a PPCA model with random initialization
update_sigma!(gp_reg, 0.0)                               # No regularization
num_training_itterations = 11  # First test that it's working
dataset = Base.Iterators.repeated((ppca,spec),num_training_itterations)
#dataset = Base.Iterators.repeated((ppca,spec),num_training_itterations)
opt = NADAM()          # Choose Optimization algorithm
loss_tol = 0.01
itterations = 0
last_loss = Inf
consecutive_passes = 0
evalcb = function ()   # Callback to show loss and stop optimization
    global last_loss, loss_tol, consecutive_passes
    cs = chisq(ppca,spec)
    reg = regularization(ppca)
    l = cs + reg
    global itterations += 1
    if itterations%10 == 1
        println( string(itterations) * ": χ^2 = ", cs, " reg = ", reg, " loss = ", l)
    end
    num_dof = size(spec.flux,1)*size(spec.flux,2) - get_num_param(ppca) - get_num_inputs(ppca)*get_num_obs(spec)
    #if cs < min( 0.1*num_dof, num_dof - 2.0*sqrt(num_dof) )
    if (last_loss < l+loss_tol) && ( l <= last_loss)
        consecutive_passes += 1
        println( string(itterations) * ": χ^2 = ", cs, " reg = ", reg, " loss = ", l)
        if consecutive_passes >=10
            println("Hurray!")
            Flux.stop()
        end
    else
        consecutive_passes = 0
    end
    last_loss = l
end
#ppca_par = Flux.params(ppca,x)

ppca_par = Flux.params(ppca)
delete!(ppca_par, ppca.x)      # Fix x at true values
ppca.x .= x_true[1:3,:]
Flux.train!(loss, ppca_par, dataset, opt, cb = evalcb)

mean(abs.(ppca.x.-x_true[1:3,:]))

update_sigma!(gp_reg, 0.0)
num_training_itterations = 10000
dataset = Base.Iterators.repeated((ppca,spec),num_training_itterations)
Flux.train!(loss, ppca_par, dataset, opt, cb = evalcb)

plt1 = plot(ppca_true.lambda,ppca_true.b, label="true", title="Mean Vector")
plot!(plt1,ppca.lambda,ppca.b, label="est")
plt2 = plot(ppca.lambda,ppca.b.-ppca_true.b, label="resid")
plot(plt1, plt2, layout =  @layout [a ; b ] )

idx_plt = 1:(floor(Int64,size(ppca_true.W,1)/2))
plt3 = plot(ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,1]./(sqrt(2)*std(ppca_true.W[idx_plt,1])), color=:blue, title="Basis Vectors", legend=:none)
scatter!(plt3, ppca.lambda[idx_plt],ppca.W[idx_plt,1]./(sqrt(2)*std(ppca.W[idx_plt,1])), ms=2, color=:blue)
plot!(plt3,ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,2]./(sqrt(2)*std(ppca_true.W[idx_plt,2])), color=:red)
scatter!(plt3,ppca.lambda[idx_plt],ppca.W[idx_plt,2]./(sqrt(2)*std(ppca.W[idx_plt,2])), ms=2,color=:red)
plot!(plt3,ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])), color=:green)
scatter!(plt3,ppca.lambda[idx_plt],ppca.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])), ms=2, color=:green)
plt4 = plot( ppca.lambda[idx_plt],ppca.W[idx_plt,1]./(sqrt(2)*std(ppca.W[idx_plt,1])).-ppca_true.W[idx_plt,1]./(sqrt(2)*std(ppca_true.W[idx_plt,1])), label="W1", color=:blue)
plot!(plt4,ppca.lambda[idx_plt],ppca.W[idx_plt,2]./(sqrt(2)*std(ppca.W[idx_plt,2])).-ppca_true.W[idx_plt,2]./(sqrt(2)*std(ppca_true.W[idx_plt,2])), label="W2", color=:red)
plot!(plt4,ppca.lambda[idx_plt],ppca.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])).-ppca_true.W[idx_plt,3]./(sqrt(2)*std(ppca_true.W[idx_plt,3])), label="W3", color=:green)
plot(plt3, plt4, layout =  @layout [a ; b ] )

std(ppca.b.-ppca_true.b), [std(ppca.W[:,i]./(sqrt(2)*std(ppca.W[:,i])).-ppca_true.W[:,i]./(sqrt(2)*std(ppca_true.W[:,i]))) for i in 1:3]

ppca_par = Flux.params(ppca)    # Make x active again
update_sigma!(gp_reg, 0.0)
num_training_itterations = 10000
dataset = Base.Iterators.repeated((ppca,spec),num_training_itterations)
Flux.train!(loss, ppca_par, dataset, opt, cb = evalcb)

std(ppca.b.-ppca_true.b), [std(ppca.W[:,i]./(sqrt(2)*std(ppca.W[:,i])).-ppca_true.W[:,i]./(sqrt(2)*std(ppca_true.W[:,i]))) for i in 1:3]

mean(abs.(ppca.x.-x_true[1:3,:]))

plt1 = plot(ppca_true.lambda,ppca_true.b, label="true", title="Mean Vector")
plot!(plt1,ppca.lambda,ppca.b, label="est")
plt2 = plot(ppca.lambda,ppca.b.-ppca_true.b, label="resid")
plot(plt1, plt2, layout =  @layout [a ; b ] )

idx_plt = 1:(floor(Int64,size(ppca_true.W,1)/2))
plt3 = plot(ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,1]./(sqrt(2)*std(ppca_true.W[idx_plt,1])), color=:blue, title="Basis Vectors", legend=:none)
scatter!(plt3, ppca.lambda[idx_plt],ppca.W[idx_plt,1]./(sqrt(2)*std(ppca.W[idx_plt,1])), ms=2, color=:blue)
plot!(plt3,ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,2]./(sqrt(2)*std(ppca_true.W[idx_plt,2])), color=:red)
scatter!(plt3,ppca.lambda[idx_plt],ppca.W[idx_plt,2]./(sqrt(2)*std(ppca.W[idx_plt,2])), ms=2,color=:red)
plot!(plt3,ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])), color=:green)
scatter!(plt3,ppca.lambda[idx_plt],ppca.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])), ms=2, color=:green)
plt4 = plot( ppca.lambda[idx_plt],ppca.W[idx_plt,1]./(sqrt(2)*std(ppca.W[idx_plt,1])).-ppca_true.W[idx_plt,1]./(sqrt(2)*std(ppca_true.W[idx_plt,1])), label="W1", color=:blue)
plot!(plt4,ppca.lambda[idx_plt],ppca.W[idx_plt,2]./(sqrt(2)*std(ppca.W[idx_plt,2])).-ppca_true.W[idx_plt,2]./(sqrt(2)*std(ppca_true.W[idx_plt,2])), label="W2", color=:red)
plot!(plt4,ppca.lambda[idx_plt],ppca.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])).-ppca_true.W[idx_plt,3]./(sqrt(2)*std(ppca_true.W[idx_plt,3])), label="W3", color=:green)
plot(plt3, plt4, layout =  @layout [a ; b ] )

std(ppca.b.-ppca_true.b), [std(ppca.W[:,i]./(sqrt(2)*std(ppca.W[:,i])).-ppca_true.W[:,i]./(sqrt(2)*std(ppca_true.W[:,i]))) for i in 1:3]

#update_lambda!(gp_reg, 0.2)
update_sigma!(gp_reg, 3)
num_training_itterations = 10000
#ppca = PPCAGridded(num_inputs,num_outputs, num_obs, lambda_range)  # Uncomment to start fresh
dataset = Base.Iterators.repeated((ppca,spec),num_training_itterations)
ppca_par = params(ppca)
itterations = 0
#loss_tol = 0.01
last_loss = Inf
Flux.train!(loss, ppca_par, dataset, opt, cb = evalcb)

std(ppca.b.-ppca_true.b), [std(ppca.W[:,i]./(sqrt(2)*std(ppca.W[:,i])).-ppca_true.W[:,i]./(sqrt(2)*std(ppca_true.W[:,i]))) for i in 1:3]

idx_plt = 1:(floor(Int64,size(ppca_true.W,1)/1))
plt7 = plot(ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,1], color=:blue, label=:none)
scatter!(plt7, ppca_true.lambda[idx_plt],ppca.W[idx_plt,1]./(sqrt(2)*std(ppca.W[idx_plt,1])), ms=2, color=:blue)
plot!(plt7,ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,2], color=:red, label=:none)
scatter!(plt7,ppca_true.lambda[idx_plt],ppca.W[idx_plt,2]./(sqrt(2)*std(ppca.W[idx_plt,2])), ms=2,color=:red)
plot!(plt7,ppca_true.lambda[idx_plt],ppca_true.W[idx_plt,3], color=:green, label=:none)
scatter!(plt7,ppca_true.lambda[idx_plt],ppca.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])), ms=2, color=:green)
plt8 = plot( ppca_true.lambda[idx_plt],ppca.W[idx_plt,1]./(sqrt(2)*std(ppca.W[idx_plt,1])).-ppca_true.W[idx_plt,1]./(sqrt(2)*std(ppca_true.W[idx_plt,1])),  color=:blue, label=:none, title="Basis Vector Residuals (w/ reg)")
plot!(plt8,ppca_true.lambda[idx_plt],ppca.W[idx_plt,2]./(sqrt(2)*std(ppca.W[idx_plt,2])).-ppca_true.W[idx_plt,2]./(sqrt(2)*std(ppca_true.W[idx_plt,2])), color=:red, label=:none)
plot!(plt8,ppca_true.lambda[idx_plt],ppca.W[idx_plt,3]./(sqrt(2)*std(ppca.W[idx_plt,3])).-ppca_true.W[idx_plt,3]./(sqrt(2)*std(ppca_true.W[idx_plt,3])),  color=:green, label=:none)
title!(plt4,"Basis Vector Residuals (w/o reg)")
plot(plt4, plt8, layout =  @layout [a ; b ] )

plt5 = plot(ppca_true.lambda,ppca_true.b, label=:none) # "true")
plot!(plt5,ppca.lambda,ppca.b, label=:none) #"est")
plt6 = plot(ppca.lambda,ppca.b.-ppca_true.b, label=:none) # "resid")
plot(plt1, plt2, plt5, plt6, layout =  @layout [a b ; c d ] )

std(ppca.b.-ppca_true.b), [std(ppca.W[:,i]./(sqrt(2)*std(ppca.W[:,i])).-ppca_true.W[:,i]./(sqrt(2)*std(ppca_true.W[:,i]))) for i in 1:3]

mean(abs.(ppca.x.-x_true[1:3,:]))

function chisq_versus_reg_const(reg_const)
    ppca_tmp = copy(ppca)
    update_sigma!(gp_reg, reg_const)
    num_training_itterations = 10000
    #ppca = PPCAGridded(num_inputs,num_outputs, lambda_range)
    #dataset = Base.Iterators.repeated((x,ppca,spec),num_training_itterations)
    dataset = Base.Iterators.repeated((ppca_tmp,spec),num_training_itterations)
    #ppca_par = params(ppca,x)
    ppca_par = params(ppca_tmp)
    itterations = 0
    loss_tol = 0.05
    last_loss = Inf
    Flux.train!(loss, ppca_par, dataset, opt, cb = evalcb)
    #cs = chisq(x,ppca,spec)
    cs = chisq(ppca.x,ppca_tmp,spec)
end

reg_const_list = [0.3, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20. ]
#update_lambda!(gp_reg, 0.2)
chisq_list = chisq_versus_reg_const.(reg_const_list)
reg_const_list[1] = 0.1 # kludge for log
chisq_list

plt_idx = 1:length(reg_const_list) - 0
scatter(log10.(reg_const_list[plt_idx]),log10.(chisq_list[plt_idx]),legend=:none)
xlabel!("log_10 Regularization Constant")
ylabel!("log_10 chi^2")

struct PPCAInterpolated <: AbstractPPCA
  W::Array{Float64,2}
  b::Array{Float64,1}
  x::Array{Float64,2}

  lambda::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
  w_lo_array::Array{Float64,2}
  idx_lo_array::Array{Int64,2}

  function PPCAInterpolated(W::AA2, b::AA1, x::AA2, lambda_range::R; num_pad::Integer = 0) where { T<:Real, AA1<:AbstractArray{T,1}, AA2<:AbstractArray{T,2}, R<:AbstractRange }
     @assert size(W,1) == size(b,1)
     @assert size(W,2) == size(x,1)
     @assert length(lambda_range) == size(b,1)
     num_out = size(b,1) + 2 * num_pad
     num_obs = size(x,2)
     #@assert size(w_lo_array,1) ==  size(b,1)
     #@assert size(w_lo_array,2) ==  size(b,1)
     #@assert size(idx_lo_array) == size(w_lo_array)
     if num_pad == 0
        return new(W, b, x, lambda_range, zeros(Int64, num_out,num_obs), zeros(num_out,num_obs) )
     else
        @error "Need to write code to automatically pad PPCAInterpoalted's grid"
     end

  end

  function PPCAInterpolated(in::PPCAInterpolated)
     @assert size(in.W,1) == size(in.b,1)
     @assert size(in.W,2) == size(in.x,1)
     @assert length(in.lambda) == size(in.b,1)
    new(in.W,in.b,in.x,in.lambda,in.w_lo_array,in.idx_lo_array)
  end
end

function PPCAInterpolated(num_in::Integer, num_out::Integer, lambda_range::R, num_obs::Integer) where {R<:AbstractRange}
    @assert length(lambda_range) == num_out
    PPCAInterpolated(randn(num_out, num_in), randn(num_out), lambda_range, num_obs )
end

get_num_inputs(m::PPCAInterpolated) = size(m.W,2)
get_num_outputs(m::PPCAInterpolated) = size(m.W,1)
get_num_obs(m::PPCAInterpolated) = size(m.x,2)
get_num_param(m::PPCAInterpolated) = get_num_outputs(m)*(get_num_inputs(m)+1)+get_num_obs(m)*get_num_inputs(m)

Flux.@functor PPCAInterpolated (W,b)

function myfindlast(x::T; r::R) where {T<:Real, R<:StepRangeLen{T,Base.TwicePrecision{T},Base.TwicePrecision{T}} }
    idx_lo = min(max(1,ceil(Int64,(x-r.ref.hi)/r.step.hi)),length(r)-1)
end

function myfindlast(x::T; r::R) where { T<:Real, R<:AbstractArray{T} }
    idx_lo = min(max(1,ceil(Int64,convert(T,(x-r.ref)/r.step))),length(r)-1)
end

function calc_interp_w_lo(x, idx_lo; r::AR) where {T<:Real, AR<:AbstractRange }
    idx_hi = idx_lo + 1
    denom = r.step
    w = (r[idx_hi]-x)/denom
end

function calc_interp_w_lo(x, idx_lo; r::AA) where {T<:Real, AA<:AbstractArray{T} }
    idx_hi = idx_lo + 1
    denom = (r[idx_hi]-r[idx_lo])
    w = (r[idx_hi]-x)/denom
end

function update_weights!(m::PPCAInterpolated, lambda::Array{T,2}) where { T<:Real }
    m.idx_lo_array .= myfindlast.(lambda, r=m.lambda)
    m.w_lo_array   .= calc_interp_w_lo.(lambda, m.idx_lo_array, r=m.lambda)
    return m
end

function predict(m::PPCAInterpolated, x::AbstractArray{T,2}, i::Integer) where {T<:Real}
    @assert length(m.w_lo_array) >= 1
    idx = view(m.idx_lo_array,:,i)
    w = view(m.w_lo_array,:,i)
    #=
    Wtmp = w.*m.W[idx,:] + (1.0.-w).*m.W[idx.+1,:]
    btmp = w.*m.b[idx]   + (1.0.-w).*m.b[idx.+1]
    Wtmp*x[:,i].+btmp
    =#
          w .*(m.W[idx,:]   *x[:,i].+m.b[idx]   ) +
    (1.0.-w).*(m.W[idx.+1,:]*x[:,i].+m.b[idx.+1])
end

function predict(m::PPCAInterpolated, i::Integer) where {T<:Real}
    predict(m,m.x,i)
end

function (m::PPCAInterpolated)(x::AbstractArray{T,2})  where {T<:Real}
    mapreduce(i->predict(m,x,i), hcat,1:get_num_obs(m))  # WARNING: mapreduce can't be Autodiffed by Flux
end

function chisq(x::AbstractArray{T,2}, m::PPCAInterpolated, ss::SpectraSet1D)  where { T<:Real }
    @assert get_num_obs(m) == get_num_obs(ss)
    χ2 = zero(T)
    for i in 1:get_num_obs(m)
        yhat = predict(m,i)
        χ2 += sum(((ss.flux[:,i].-yhat).^2 ./ss.flux_var[:,i]))
    end
    return χ2
end

chisq(m::PPCAInterpolated, ss::SpectraSet1D) = chisq(m.x,m,ss)

function regularization(x::AbstractArray{T,2}, m::PPCAInterpolated, ss::SpectraSet1D)  where { T<:Real }
    regularization_scores(x) + regularization_basis(m.W) + regularization_mean(m.b, mu=1.0)
end

regularization(m::PPCAInterpolated, ss::SpectraSet1D) = regularization(m.x,m,ss)

ppca_interp = PPCAInterpolated(ppca.W,ppca.b,ppca.x,ppca.lambda, num_pad=0 )
update_weights!(ppca_interp, spec.lambda_ssb)

maximum(abs.(ppca_interp(ppca_interp.x).-ppca(ppca.x)))

chisq(ppca.x,ppca,spec), chisq(ppca_interp.x,ppca_interp,spec)

lmin, lmax = ppca.lambda[1], ppca.lambda[end]
lambda_pert = ppca.lambda .+ convert(Float64,5*ppca.lambda.step) .* sin.(2π*(ppca.lambda.-lmin)./(lmax-lmin))
spec.lambda_ssb .= lambda_pert
update_weights!(ppca_interp, spec.lambda_ssb)

chisq(ppca.x,ppca,spec), chisq(ppca_interp.x,ppca_interp,spec)

update_sigma!(gp_reg, 4.0)                               # No regularization
num_training_itterations = 11  # First test that it's working
dataset = Base.Iterators.repeated((ppca_interp,spec),num_training_itterations)
opt = NADAM()          # Choose Optimization algorithm
loss_tol = 0.01
itterations = 0
last_loss = Inf
consecutive_passes = 0
evalcb = function ()   # Callback to show loss and stop optimization
    global last_loss, loss_tol, consecutive_passes
    cs = chisq(ppca_interp.x,ppca_interp,spec)
    reg = regularization(ppca_interp.x,ppca_interp,spec)
    l = cs + reg
    global itterations += 1
    if itterations%10 == 1
        println( string(itterations) * ": χ^2 = ", cs, " reg = ", reg, " loss = ", l)
    end
    num_dof = size(spec.flux,1)*size(spec.flux,2) - get_num_param(ppca_interp) - get_num_inputs(ppca_interp)*get_num_obs(spec)
    #if cs < min( 0.1*num_dof, num_dof - 2.0*sqrt(num_dof) )
    if (last_loss < l+loss_tol) && ( l <= last_loss)
        consecutive_passes += 1
        println( string(itterations) * ": χ^2 = ", cs, " reg = ", reg, " loss = ", l)
        if consecutive_passes >=10
            println("Hurray!")
            Flux.stop()
        end
    else
        consecutive_passes = 0
    end
    last_loss = l
end
ppcai_par = params(ppca_interp)
Flux.train!(loss, ppcai_par, dataset, opt, cb = evalcb)

num_training_itterations = 3000
dataset = Base.Iterators.repeated((x,ppca_interp,spec),num_training_itterations)
Flux.train!(loss, ppcai_par, dataset, opt, cb = evalcb)

plt1b = plot(ppca_true.lambda,ppca_true.b, label="true", title="Mean")
plot!(plt1b,ppca_interp.lambda,ppca_interp.b, label="est")
plt2b = plot(ppca_true.lambda,ppca_interp.b.-ppca_true.b, label="resid")
plot(plt1b, plt2b, layout =  @layout [a ; b ] )

plt3b = plot(ppca_true.lambda,ppca_true.W[:,1], color=:blue, title="Basis Vectors (after λ transform)")
scatter!(plt3b,ppca.lambda,ppca.W[:,1]./(sqrt(2)*std(ppca.W[:,1])), ms=2, color=:blue)
plot!(plt3b,ppca_true.lambda,ppca_true.W[:,2], color=:red)
scatter!(plt3b,ppca.lambda,ppca.W[:,2]./(sqrt(2)*std(ppca.W[:,2])), ms=2,color=:red)
plot!(plt3b,ppca_true.lambda,ppca_true.W[:,3], color=:green)
scatter!(plt3b,ppca.lambda,ppca.W[:,3]./(sqrt(2)*std(ppca.W[:,3])), ms=2, color=:green)
plt4b = plot( ppca_true.lambda,ppca.W[:,1]./(sqrt(2)*std(ppca.W[:,1])).-ppca_true.W[:,1]./(sqrt(2)*std(ppca_true.W[:,1])), label="W1", color=:blue)
plot!(plt4b,ppca_true.lambda,ppca.W[:,2]./(sqrt(2)*std(ppca.W[:,2])).-ppca_true.W[:,2]./(sqrt(2)*std(ppca_true.W[:,2])), label="W2", color=:red)
plot!(plt4b,ppca_true.lambda,ppca.W[:,3]./(sqrt(2)*std(ppca.W[:,3])).-ppca_true.W[:,3]./(sqrt(2)*std(ppca_true.W[:,3])), label="W3", color=:green)
plot(plt3b, plt4b, layout =  @layout [a ; b ] )

plt5b = plot(ppca_true.lambda,ppca_true.b, label=:none)
plot!(plt5b,ppca.lambda,ppca_interp.b, label=:none)
plt6b = plot(ppca.lambda,ppca_interp.b.-ppca_true.b, label=:none)
title!(plt1b, "Mean before transform")
title!(plt5b, "Mean after transform")
plot(plt1b, plt2b, plt5b, plt6b, layout =  @layout [a b ; c d ] )

plt7b = plot(ppca_true.lambda,ppca_true.W[:,1], color=:blue, label=:none, title="Basis Vectors")
scatter!(plt7b, ppca_true.lambda,ppca_interp.W[:,1]./(sqrt(2)*std(ppca_interp.W[:,1])), ms=2, color=:blue)
plot!(plt7b,ppca_true.lambda,ppca_true.W[:,2], color=:red, label=:none)
scatter!(plt7b,ppca_true.lambda,ppca_interp.W[:,2]./(sqrt(2)*std(ppca_interp.W[:,2])), ms=2,color=:red)
plot!(plt7b,ppca_true.lambda,ppca_true.W[:,3], color=:green, label=:none)
scatter!(plt7b,ppca_true.lambda,ppca_interp.W[:,3]./(sqrt(2)*std(ppca_interp.W[:,3])), ms=2, color=:green)
plt8b = plot( ppca_true.lambda,ppca_interp.W[:,1]./(sqrt(2)*std(ppca_interp.W[:,1])).-ppca_true.W[:,1]./(sqrt(2)*std(ppca_true.W[:,1])),  color=:blue, label=:none)
plot!(plt8b,ppca_true.lambda,ppca_interp.W[:,2]./(sqrt(2)*std(ppca_interp.W[:,2])).-ppca_true.W[:,2]./(sqrt(2)*std(ppca_true.W[:,2])), color=:red, label=:none)
plot!(plt8b,ppca_true.lambda,ppca_interp.W[:,3]./(sqrt(2)*std(ppca_interp.W[:,3])).-ppca_true.W[:,3]./(sqrt(2)*std(ppca_true.W[:,3])),  color=:green, label=:none)
title!(plt4b,"Basis Vector Residuals before transform")
title!(plt8b,"Basis Vector Residuals after transform")
plot(plt4b, plt8b, layout =  @layout [a ; b ] )

struct PPCAGPBasis <: AbstractPPCA
  W::Array{<:Stheno.AbstractGP,1}
  b::Stheno.AbstractGP
  x::Array{Float64,2}
  W_prior::Array{<:Stheno.AbstractGP,1}
  b_prior::Stheno.AbstractGP
  lambda
#=
  function PPCAGPBasis(W::AbstractArray{WelT,1}, b::BT, x::AA2, lambda_range::R; num_pad::Integer = 0) where {
        WelT<:AbsGp,1}, BT<:AbsGp, R<:AbstractRange }
     @assert length(first(W)) == length(b)
     @assert size(W,1) == size(x,1)
     #@assert length(lambda_range) == length(b)
     num_out = length(b)
     num_obs = size(x,2)
     return new(W, b, x, lambda_range )
  end

  function PPCAGPBasis(in::PPCAGPBasis)
    new(in.W,in.b,in.x,in.lambda)
  end
    =#
end

function PPCAGPBasis(num_basis::Integer, lambda_range::R, num_obs::Integer) where {R<:AbstractRange}
    @assert 3 <= length(lambda_range) <= 512
    AbsGP = Stheno.AbstractG
    num_out = length(lambda_range)
    W = Array{AbsGP,1}(undef,num_basis)
    b = GP(k, GPC())
    W_prior = Array{AbsGP,1}(GP(k, GPC()),num_basis)
    b_prior = GP(k, GPC())
    x = randn(num_basis,num_obs)
    PPCAGPBasis(W,b,x,lambda_range, W_prior, b_prior)
end

get_num_inputs(m::PPCAGPBasis) = length(m.lambda) # or length(first(m.W))
get_num_outputs(m::PPCAGPBasis) = length(m.W)
get_num_obs(m::PPCAGPBasis) = size(m.x,2)
get_num_param(m::PPCAGPBasis) = get_num_outputs(m)*(get_num_inputs(m)+1)+get_num_obs(m)*get_num_outputs(m)

Flux.@functor PPCAGPBasis (W,b)

function predict(m::PPCAGPBasis, x::AbstractArray{T,2}, i::Integer, lambda::AbstractArray{T,1} ) where {T<:Real}
    xtmp = view(x,:,i)
    Wtmp = mapreduce(j->mean(m.W[j](lambda)), hcat,1:length(m.W) )
    btmp = mean(m.b(lambda))
    ypred = Wtmp*xtmp+btmp
end

function predict(m::PPCAGPBasis, i::Integer, lambda::AbstractArray{T,1} ) where {T<:Real}
    predict(m,m.x,i, lambda)
end

function (m::PPCAGPBasis)(x::AbstractArray{T,2})  where {T<:Real}
    mapreduce(i->predict(m,x,i), hcat,1:get_num_obs(m))  # WARNING: mapreduce can't be Autodiffed by Flux
end

function chisq(x::AbstractArray{T,2}, m::PPCAGPBasis, ss::SpectraSet1D)  where { T<:Real }
    @assert get_num_obs(m) == get_num_obs(ss)
    χ2 = zero(T)
    for i in 1:get_num_obs(m)
        yhat = predict(m,i,ss.lambda_ssb[:,i])
        χ2 += sum(((ss.flux[:,i].-yhat).^2 ./ss.flux_var[:,i]))
    end
    return χ2
end

chisq(m::PPCAGPBasis, ss::SpectraSet1D) = chisq(m.x,m,ss)

function regularization_mean(m::PPCAGPBasis; mu=zero(T), lambda = m.lambda )  where { T<:Real, GPT<:Stheno.AbstractGP }
    -logpdf(m.b_prior(lambda), mean(m.b(lambda)).-mu)
end

function regularization_basis(m::PPCAGPBasis, i::Integer, lambda = m.lambda)  where { GPT<:Stheno.AbstractGP }
    -logpdf(m.W_prior[i](lambda), mean(m.W[i](lambda)) )
end

function regularization_basis(m::PPCAGPBasis)  where { GPT<:Stheno.AbstractGP }
    sum( [regularization_basis(m,i) for i in 1:get_num_outputs(m) ] )
end

function regularization(x::AbstractArray{T,2}, m::PPCAGPBasis #=, ss::SpectraSet1D =#)  where { T<:Real }
    regularization_scores(x) + regularization_basis(m) + regularization_mean(m, mu=1.0)
end

function regularization(m::PPCAGPBasis #=, ss::SpectraSet1D=# )  where { T<:Real }
    regularization_scores(m.x) + regularization_basis(m) + regularization_mean(m, mu=1.0)
end

#= Old demo code

ktmp = 1.0 * stretch(Matern52(), 1 /1.0)
bGP = GP(1.0, ktmp, GPC())
b_prior = GP(1.0, ktmp, GPC())
W_prior = [GP(1.0, ktmp, GPC()) for j in 1:size(ppca.W,2)]
bGPx = bGP(ppca.lambda, 0.01)
bGPpost = bGP | Obs( bGPx , ppca.b)
WGPpost = [bGP | Obs( bGPx , ppca.W[:,j]) for j in 1:size(ppca.W,2)]
#ppca_gp = PPCAGPBasis(ppca.W,ppca.b,ppca.x,ppca.lambda)
#update_weights!(ppca_interp, spec.lambda_ssb)
#scatter(ppca.lambda, ppca.W[:,1])
plt = plot(ppca.lambda, mean(bGPpost(ppca.lambda)), color=1, label="b")
scatter!(plt,ppca.lambda, ppca.b, label=:none)

for j in 1:size(ppca.W,2)
    plot!(plt,ppca.lambda, mean(WGPpost[j](ppca.lambda)), color = j+1,  label = "W"*string(j))
    scatter!(plt,ppca.lambda, ppca.W[:,j], color=j+1, label=:none)
end
display(plt)

ppca_gp = PPCAGPBasis(WGPpost,bGPpost,ppca.x, W_prior,b_prior,lambda_range)

plt = plot(ppca_gp.lambda, mean(ppca_gp.b(ppca_gp.lambda)), color=1, label="b")
scatter!(plt,ppca.lambda, ppca.b, label=:none)
#=
for j in 1:size(ppca.W,2)
    plot!(plt,ppca.lambda, mean(WGPpost[j](ppca.lambda)), color = j+1,  label = "W"*string(j))
    scatter!(plt,ppca.lambda, ppca.W[:,j], color=j+1, label=:none)
end
=#
display(plt)

mean(ppca_gp.b(ppca_gp.lambda))

@time cs = chisq(ppca_gp.x,ppca_gp,spec)

#update_sigma!(gp_reg, 0.0)                               # No regularization
num_training_itterations = 1  # First test that it's working
dataset = Base.Iterators.repeated((ppca_gp,spec),num_training_itterations)
opt = NADAM()          # Choose Optimization algorithm
loss_tol = 0.01
itterations = 0
last_loss = Inf
consecutive_passes = 0
evalcb = function ()   # Callback to show loss and stop optimization
    global last_loss, loss_tol, consecutive_passes
    cs = chisq(ppca_gp.x,ppca_gp,spec)
    reg = regularization(ppca_gp.x,ppca_gp)
    l = cs + reg
    global itterations += 1
    if itterations%10 == 1
        println( string(itterations) * ": χ^2 = ", cs, " reg = ", reg, " loss = ", l)
    end
    num_dof = size(spec.flux,1)*size(spec.flux,2) - get_num_param(ppca_gp) - get_num_inputs(ppca_gp)*get_num_obs(spec)
    #if cs < min( 0.1*num_dof, num_dof - 2.0*sqrt(num_dof) )
    if (last_loss < l+loss_tol) && ( l <= last_loss)
        consecutive_passes += 1
        println( string(itterations) * ": χ^2 = ", cs, " reg = ", reg, " loss = ", l)
        if consecutive_passes >=10
            println("Hurray!")
            Flux.stop()
        end
    else
        consecutive_passes = 0
    end
    last_loss = l
end
ppca_gp_par = params(ppca_gp)
Flux.train!(loss, ppca_gp_par, dataset, opt, cb = evalcb)

num_training_itterations = 10
dataset = Base.Iterators.repeated((x,ppca_gp,spec),num_training_itterations)
@time Flux.train!(loss, ppca_gp_par, dataset, opt, cb = evalcb)

plt1b = plot(ppca_true.lambda,ppca_true.b, label="true", title="Mean")
plot!(plt1b,ppca_gp.lambda,ppca_gp.b, label="est")
plt2b = plot(ppca_true.lambda,ppca_gp.b.-ppca_true.b, label="resid")
plot(plt1b, plt2b, layout =  @layout [a ; b ] )

plt3b = plot(ppca_true.lambda,ppca_true.W[:,1], color=:blue, title="Basis Vectors (after λ transform)")
scatter!(plt3b,ppca.lambda,ppca.W[:,1]./(sqrt(2)*std(ppca.W[:,1])), ms=2, color=:blue)
plot!(plt3b,ppca_true.lambda,ppca_true.W[:,2], color=:red)
scatter!(plt3b,ppca.lambda,ppca.W[:,2]./(sqrt(2)*std(ppca.W[:,2])), ms=2,color=:red)
plot!(plt3b,ppca_true.lambda,ppca_true.W[:,3], color=:green)
scatter!(plt3b,ppca.lambda,ppca.W[:,3]./(sqrt(2)*std(ppca.W[:,3])), ms=2, color=:green)
plt4b = plot( ppca_true.lambda,ppca.W[:,1]./(sqrt(2)*std(ppca.W[:,1])).-ppca_true.W[:,1]./(sqrt(2)*std(ppca_true.W[:,1])), label="W1", color=:blue)
plot!(plt4b,ppca_true.lambda,ppca.W[:,2]./(sqrt(2)*std(ppca.W[:,2])).-ppca_true.W[:,2]./(sqrt(2)*std(ppca_true.W[:,2])), label="W2", color=:red)
plot!(plt4b,ppca_true.lambda,ppca.W[:,3]./(sqrt(2)*std(ppca.W[:,3])).-ppca_true.W[:,3]./(sqrt(2)*std(ppca_true.W[:,3])), label="W3", color=:green)
plot(plt3b, plt4b, layout =  @layout [a ; b ] )

plt5b = plot(ppca_true.lambda,ppca_true.b, label=:none)
plot!(plt5b,ppca.lambda,ppca_interp.b, label=:none)
plt6b = plot(ppca.lambda,ppca_interp.b.-ppca_true.b, label=:none)
title!(plt1b, "Mean before transform")
title!(plt5b, "Mean after transform")
plot(plt1b, plt2b, plt5b, plt6b, layout =  @layout [a b ; c d ] )

plt7b = plot(ppca_true.lambda,ppca_true.W[:,1], color=:blue, label=:none, title="Basis Vectors")
scatter!(plt7b, ppca_true.lambda,ppca_interp.W[:,1]./(sqrt(2)*std(ppca_interp.W[:,1])), ms=2, color=:blue)
plot!(plt7b,ppca_true.lambda,ppca_true.W[:,2], color=:red, label=:none)
scatter!(plt7b,ppca_true.lambda,ppca_interp.W[:,2]./(sqrt(2)*std(ppca_interp.W[:,2])), ms=2,color=:red)
plot!(plt7b,ppca_true.lambda,ppca_true.W[:,3], color=:green, label=:none)
scatter!(plt7b,ppca_true.lambda,ppca_interp.W[:,3]./(sqrt(2)*std(ppca_interp.W[:,3])), ms=2, color=:green)
plt8b = plot( ppca_true.lambda,ppca_interp.W[:,1]./(sqrt(2)*std(ppca_interp.W[:,1])).-ppca_true.W[:,1]./(sqrt(2)*std(ppca_true.W[:,1])),  color=:blue, label=:none)
plot!(plt8b,ppca_true.lambda,ppca_interp.W[:,2]./(sqrt(2)*std(ppca_interp.W[:,2])).-ppca_true.W[:,2]./(sqrt(2)*std(ppca_true.W[:,2])), color=:red, label=:none)
plot!(plt8b,ppca_true.lambda,ppca_interp.W[:,3]./(sqrt(2)*std(ppca_interp.W[:,3])).-ppca_true.W[:,3]./(sqrt(2)*std(ppca_true.W[:,3])),  color=:green, label=:none)
title!(plt4b,"Basis Vector Residuals before transform")
title!(plt8b,"Basis Vector Residuals after transform")
plot(plt4b, plt8b, layout =  @layout [a ; b ] )

=#

end # module