text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, BlockBandedMatrices,  LinearAlgebra, Test
import ApproxFunBase: Multiplication,InterlaceOperator, Block, ∞,
        testfunctional, testbandedoperator, testraggedbelowoperator, testinfoperator, testblockbandedoperator
import ApproxFunOrthogonalPolynomials: JacobiZ

@testset "Operator" begin
    @testset "Evaluation" begin
        testfunctional(Evaluation(Ultraspherical(1),0.1))
        d = -4 .. 4
        f = Fun(exp, Ultraspherical(1,d))
        @test f(-4) ≈ Number(ldirichlet()*f) ≈ Number(Evaluation(Ultraspherical(1,d),-4)*f)

        testfunctional(Evaluation(Chebyshev(),0.1,1))
        testfunctional(Evaluation(Chebyshev(),0.1,1)-Evaluation(Chebyshev(),0.1,1))

        let f = Fun(cos)
            @test (Evaluation(Chebyshev(),0.1,1)*f)(0.1)  ≈ f'(0.1)
        end
    end

    @testset "Conversion" begin
        C=Conversion(Ultraspherical(1),Ultraspherical(2))
        testbandedoperator(C)

        @test Matrix(C[1:5,1:5])  ≈   [1.0 0.0 -0.3333333333333333 0.0  0.0
                                              0.0 0.5  0.0               -0.25 0.0
                                              0.0 0.0  0.3333333333333333 0.0 -0.2
                                              0.0 0.0  0.0                0.25 0.0
                                              0.0 0.0  0.0                0.0  0.2]
        d=ChebyshevInterval()
        A=Conversion(Chebyshev(d),Ultraspherical(2,d))
        x = Fun()
        f=Fun(exp)
        @test AbstractMatrix(view(A.op, Block.(1:3), Block.(1:3))) isa BlockSkylineMatrix
        testbandedoperator(A)

        @test norm(A\Fun(x.*f,rangespace(A))-(x.*f)) < 100eps()

        @test (I : Chebyshev() → Chebyshev()) * Fun() ≈ Fun()
    end

    @testset "Derivative" begin
        testbandedoperator(Derivative(Ultraspherical(1)))

        d=Interval(-10.,5.);
        S=Chebyshev(d)

        @test norm(Fun(Fun(Fun(exp,S),Ultraspherical(2,d)),S)-Fun(exp,S)) < 100eps()

        @test copy(view(Derivative(Ultraspherical(1)),1:2,1:2))[1,2] ≈ Derivative(Ultraspherical(1))[1,2]
        @test exp(0.1) ≈ (Derivative()*Fun(exp,Ultraspherical(1)))(0.1)
        D=Derivative(d)
        f = Fun(exp,d)
        @test norm((Conversion(Chebyshev(d),Ultraspherical(2,d))\(D^2*f))-f'') < 100eps()

        f=Fun(x->x^2)
        D=Derivative(domain(f))
        @test norm(D*f-f')<100eps()
    end

    @testset "Toeplitz" begin
        C=ToeplitzOperator([1.,2.,3.],[4.,5.,6.])

        @time testbandedoperator(C)

        @test Matrix(C[1:5,1:5])  ≈  [4.0 5.0 6.0 0.0  0.0
                                             1.0 4.0 5.0 6.0 0.0
                                             2.0 1.0 4.0 5.0 6.0
                                             3.0 2.0 1.0 4.0 5.0
                                             0.0 3.0 2.0 1.0 4.0]

         testbandedoperator(HankelOperator([1.,2.,3.,4.,5.,6.,7.]))
    end

    @testset "Multiplication" begin
        testbandedoperator(Multiplication(Fun(Chebyshev(),[1.,2.,3.]),Chebyshev()))

        x=Fun(identity)
        X=Multiplication(x,space(x))
        f = Fun(exp)
        testbandedoperator(X)

        @test norm(X*f-(x.*f)) < 100eps()

        A=Conversion(Chebyshev(),Ultraspherical(2))*X
        @time testbandedoperator(A)

        @test norm((mul_coefficients(A,f.coefficients))-coefficients(x.*f,rangespace(A))) < 100eps()
    end

    @testset "Integral" begin
        f=Fun(exp)
        d=domain(f)
        Q=Integral(d)
        D=Derivative(d)

        @time testbandedoperator(Q)

        @test norm((Q+I)*f-(integrate(f)+f)) < 100eps()
        @test norm((Q)*f-(integrate(f))) < 100eps()
    end

    @testset "Special functions" begin
        x=Fun(identity)
        @test norm(cos(x)-Fun(cos))<10eps()
        @test norm(sin(x)-Fun(sin))<10eps()
        @test norm(exp(x)-Fun(exp))<10eps()
        @test norm(sin(x)./x-Fun(x->sinc(x/π)))<100eps()
    end

    @testset "Permutation" begin
        P=ApproxFunBase.PermutationOperator([2,1])
        testbandedoperator(P)
        @test P[1:4,1:4] ≈ [0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0]
    end

    @testset "Mixed" begin
        d = ChebyshevInterval()
        D = Derivative(d)
        x = Fun(identity,d)
        A = D*(x*D)
        B = D+x*D^2
        C = x*D^2+D

        testbandedoperator(A)
        testbandedoperator(B)
        testbandedoperator(C)
        @time testbandedoperator(x*D)

        f=Fun(exp)
        @test (A.ops[end]*f)(0.1) ≈ f'(0.1)
        @test ((x*D)*f)(0.1) ≈ 0.1*f'(0.1)
        @test (A*f)(0.1) ≈ f'(0.1)+0.1*f''(0.1)
        @test (B*f)(0.1) ≈ f'(0.1)+0.1*f''(0.1)
        @test (C*f)(0.1) ≈ f'(0.1)+0.1*f''(0.1)

        testbandedoperator(A-B)
        testbandedoperator(B-A)
        testbandedoperator(A-C)

        @test norm((A-B)[1:10,1:10]|>Matrix) < eps()
        @test norm((B-A)[1:10,1:10]|>Matrix) < eps()
        @test norm((A-C)[1:10,1:10]|>Matrix) < eps()
        @test norm((C-A)[1:10,1:10]|>Matrix) < eps()
        @test norm((C-B)[1:10,1:10]|>Matrix) < eps()
        @test norm((B-C)[1:10,1:10]|>Matrix) < eps()
    end

    @testset "Cached" begin
        @test cache(Derivative(Chebyshev(),2))[1,1] == 0

        S=Chebyshev()
        D=Derivative(S)
        @time for padding = [true,false]
          co=ApproxFunBase.CachedOperator(D,ApproxFunBase.RaggedMatrix(Float64[],Int[1],0),(0,0),domainspace(D),rangespace(D),bandwidths(D),padding) #initialise with empty RaggedMatrix
          @test co[1:20,1:10] == D[1:20,1:10]
          @test size(co.data) == (20,10)
          ApproxFunBase.resizedata!(co,10,30)
          @test size(co.data)[2] == 30 && size(co.data)[1] ≥ 20
        end
    end

    @testset "Reverse" begin
        testbandedoperator(ApproxFunBase.Reverse(Chebyshev()))
        testbandedoperator(ApproxFunBase.ReverseOrientation(Chebyshev()))

        @test ApproxFunBase.Reverse(Chebyshev())*Fun(exp) ≈ Fun(x->exp(-x))
        @test ApproxFunBase.ReverseOrientation(Chebyshev())*Fun(exp) ≈ Fun(exp,Segment(1,-1))
    end

    @testset "Sub-operator re-view bug" begin
        D = Derivative(Chebyshev())
        S = view(D[:, 2:end], Block.(3:4), Block.(2:4))
        @test parent(S) == D
        @test parentindices(S) == (3:4,2:4)
        @test bandwidths(S)  == (-2,2)
    end

    @testset "Sub-operators" begin
        f = Fun(exp)
        D = Derivative(Chebyshev())
        testbandedoperator(D[:, 2:end])

        u = D[:,2:end] \ f
        @test u(0.1) ≈ exp(0.1)-coefficient(f,1)

        D̃ = Derivative(space(u))
        @test bandwidths(D̃) == (0,0)
        testbandedoperator(D̃)
        @test norm(u'-f) < 10eps()

        testbandedoperator(D[1:end,2:end])
        u = D[1:end,2:end] \ f
        @test u(0.1) ≈ exp(0.1)-f.coefficients[1]

        testbandedoperator(D[1:∞,2:∞])
        u = D[1:∞,2:∞] \ f
        @test u(0.1) ≈ exp(0.1)-f.coefficients[1]
    end

    @testset "InterlaceOperator" begin
        A = InterlaceOperator(Diagonal([Matrix(I,2,2),Derivative(Chebyshev())]))

        @test A[Block(1):Block(2), Block(1):Block(2)] isa BlockBandedMatrix
        @test Matrix(view(A, Block(1), Block(1))) == A[1:3,1:3]
        @test Matrix(view(A, Block(1):Block(2), Block(1):Block(2))) == A[1:4,1:4]
        testblockbandedoperator(A)
    end

    @testset "Projection and subspaces" begin
        S=Chebyshev()
        SS = S|(2:5)
        @test ApproxFunBase.block(SS,3) == Block(4)

        for C in (Operator(I,S)[3:end,:], Operator(I,S)[3:end,1:end])
            @test ApproxFunBase.domaindimension(domainspace(C)) == 1
            @test union(S,domainspace(C)) == S

            B=Dirichlet(S)

            Ai=[B;C]

            @test ApproxFunBase.colstop(Ai,1) == 2

            x=Fun()
            f=exp(x)
            u=[B;C]\[[0.,0.],f]

            @test abs(u(-1)) ≤ 10eps()
            @test abs(u(1)) ≤ 10eps()

            f=(1-x^2)*exp(x)
            u=[B;C]\[[0.,0.],f]

            @test u ≈ f
        end
    end

    @testset "Zero operator has correct bandwidths" begin
        Z=ApproxFunBase.ZeroOperator(Chebyshev())
        @test ApproxFunBase.bandwidths(Z) == ApproxFunBase.bandwidths(Z+Z)
    end

    @testset "hcat of functionals (#407)" begin
        x = Fun()
        B = [1 ldirichlet()]
        @test (B*[1;x])[1] == Fun(ConstantSpace(ApproxFunBase.Point(-1.0)),[0.0])
    end

    @testset "views of views" begin
        A = Derivative(Chebyshev()) + I
        B = A[1:2:∞,1:2:∞]
        C = B[2:∞,3:∞]
        @test A[3:2:∞,5:2:∞] == C
    end

    @testset "Multiplication functions" begin
        x = Fun()
        M = Multiplication(x, Chebyshev())
        @test exp(M).f == Multiplication(exp(x), Chebyshev()).f
        @test (-M).f == Multiplication(-x, Chebyshev()).f
        @test exp(-M).f == Multiplication(exp(-x), Chebyshev()).f
        @test (M/3).f == (3\M).f == Multiplication(x/3, Chebyshev()).f
        @test (M*3).f == (3*M).f == Multiplication(x*3, Chebyshev()).f
    end

    @testset "lastindex" begin
        Z = Operator(I,Chebyshev())
        S = view(Z,1:10,1:10)
        @test lastindex(S,1) == lastindex(S,2) == 10
        @test lastindex(S) == 100
        @test S[end,end] ≈ 1
        @test S[end-1,end] ≈ 0
    end

    @testset "Chebyshev Neumann" begin
        E = Evaluation(Chebyshev(),rightendpoint,1)
        @test E[1:4] == [0,1,4,9]
    end

    @testset "Chebyshev Dirichlet" begin
        S=Chebyshev()
        @time for io in (
                [Dirichlet(S);Derivative(Chebyshev());lneumann(S)],
                [Dirichlet(S);Derivative(Chebyshev())+Fun(cos);lneumann(S)],
                [Dirichlet(S);Derivative(Chebyshev())],
                [Dirichlet(S);Derivative(Chebyshev())+Fun(cos)],
                [Derivative(Chebyshev());Dirichlet(S)],
                [Derivative(Chebyshev())+Fun(cos);Dirichlet(S)])
            testraggedbelowoperator(io)
        end
    end

    @testset "JacobiZ" begin
        testbandedoperator(JacobiZ(Chebyshev(), 0.1))
    end
end


text:
@node Categorical Stochastic [ out, p ]

@average_energy Categorical (q_out::Categorical, q_p::Any) = -sum(probvec(q_out) .* mean(log, q_p))

ext: jl


text:
using Test
using LinearAlgebra
using QCMaterial

import Random
import PyCall: pyimport

#<1|2> = <2|1>^{\dagger}
@testset "overlap" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas_right = [pi/4, pi/6]
    thetas_left = [pi/3, pi/2]
    add_parametric_RY_gate!(c_, 1, pi/2)
    add_parametric_RZ_gate!(c_, 1, pi/2)
    c = QulacsVariationalQuantumCircuit(c_)
    #thetas = get_thetas(c)
    update_quantum_state!(c, state)
    @test QCMaterial.overlap(c, state, thetas_right, thetas_left)≈ conj(QCMaterial.overlap(c,state,thetas_left,thetas_right))
end

#B_12 = real(dot(∂1, state)*dot(state, ∂2))
#test fubini-study test B12
@testset "overlap2" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [pi/4, pi/6]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RZ_gate!(c_, 1, thetas[2])
    c = QulacsVariationalQuantumCircuit(c_)
    #thetas = get_thetas(c)
    update_quantum_state!(c, state)
    thetas_left = copy(thetas)
    delta_theta = 1e-16
    thetas_left[1] += delta_theta
    overlap1 = (overlap(c,state,thetas_left, thetas)-overlap(c,state,thetas,thetas))/delta_theta
    thetas_right = copy(thetas)
    thetas_right[2] +=delta_theta
    overlap2 = (overlap(c,state,thetas,thetas_right)-overlap(c,state,thetas,thetas))/delta_theta
    @test overlap1*overlap2≈ 0.0
end

#B_12 = real(dot(∂1, state)*dot(state, ∂2))
@testset "overlap3" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [pi/4, pi/6]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RZ_gate!(c_, 1, thetas[2])
    c = QulacsVariationalQuantumCircuit(c_)
    update_quantum_state!(c, state)
    thetas_left = copy(thetas)
    delta_theta = 1e-16
    thetas_left[1] += delta_theta
    thetas_right = copy(thetas)
    thetas_right[2] +=delta_theta
    res = real((
                (overlap(c, state, thetas_left, thetas)
                - overlap(c, state, thetas, thetas))
                *(overlap(c, state, thetas, thetas_right)
                - overlap(c, state, thetas, thetas))
            )/delta_theta^2)
    @test res≈ 0.0
end

#dot(state,∂2)
#@testset "overlap4_2" begin
#    delta_theta = 1e-3
#    n_qubit = 1
#    thetas = [0.2pi, 0.3pi]
#    state  = QulacsQuantumState(n_qubit)
#    state_delta = QulacsQuantumState(n_qubit)
#    c_ = QulacsParametricQuantumCircuit(n_qubit)
#    c_delta = QulacsParametricQuantumCircuit(n_qubit)
#    add_parametric_RY_gate!(c_, 1, thetas[1])
#    add_parametric_RZ_gate!(c_, 1, thetas[2])
#    add_parametric_RY_gate!(c_delta, 1, thetas[1])
#    add_parametric_RZ_gate!(c_delta, 1, thetas[2]+delta_theta)
#    c = QulacsVariationalQuantumCircuit(c_)
#    c_delta = QulacsVariationalQuantumCircuit(c_delta)
#    #update_quantum_state!(c, state)
#    println("state=", get_vector(state))
#    #update_quantum_state!(c_delta, state_delta)
#    @test get_vector(state) ≈ [exp(im*thetas[2]/2)*cos(thetas[1]/2), -exp(-im*thetas[2]/2)*sin(thetas[1]/2)] #state_vector is correct!!

    #<state|∂2> ver2 => |∂2> = (|state_delta> - |state>) /delta_theta => (<state|state_delta> - <state|state>) /delta_theta
#    overlap_delta_state2 = inner_product(state, state_delta)
#    overlap_state_state2=  inner_product(state, state)
#    overlap_state_partial_ver3 = (overlap_delta_state2 - overlap_state_state2)/delta_theta
#    #@test overlap_state_partial_ver3 ≈   im/2*cos(thetas[1]) atol=0.001

#    thetas_right = copy(thetas)
#    thetas_right[2] +=delta_theta
#    overlap_delta = overlap(c, state, thetas, thetas_right)
#    println("overlap_delta=", overlap_delta) #0.9999998750000023 + 0.00022571683411611032im
#    overlap_state = overlap(c, state, thetas, thetas)
#    #println("overlap_state=",overlap_state)
#    res = (
#        (overlap(c, state, thetas, thetas_right)
#       - overlap(c, state, thetas, thetas))
#       /delta_theta)
#    @test res ≈ im*cos(thetas[1])/2 atol=0.001
#end

#dot(state,∂2)
@testset "overlap4" begin
    delta_theta = 1e-3
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    state_delta = QulacsQuantumState(n_qubit)
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    c_delta = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [0.2pi, 0.3pi]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RZ_gate!(c_, 1, thetas[2])
    add_parametric_RY_gate!(c_delta, 1, thetas[1])
    add_parametric_RZ_gate!(c_delta, 1, thetas[2]+delta_theta)
    c = QulacsVariationalQuantumCircuit(c_)
    c_delta = QulacsVariationalQuantumCircuit(c_delta)
    #update_quantum_state!(c, state)
    println("state=", get_vector(state))
    #state=ComplexF64[0.8473975608908426 + 0.43177062311338915im, -0.2753361580731583 + 0.14029077970429507im]
    #update_quantum_state!(c_delta, state_delta)
    println("state_delta=", get_vector(state_delta))
    #state_delta=ComplexF64[0.8471815696635883 + 0.4321942679048537im,
                            #-0.27526597826920984 + 0.1404284302412484im]
    #@test get_vector(state) ≈ [exp(im*thetas[2]/2)*cos(thetas[1]/2), -exp(-im*thetas[2]/2)*sin(thetas[1]/2)] #state_vector is correct!!

    #println("state_delta2=", (get_vector(state_delta) ))
    #println("dif_state2=", (get_vector(state_delta) - get_vector(state)))
    #∂2
    state_partial2 = (get_vector(state_delta) - get_vector(state))/delta_theta
    #println("partial_state2=", (get_vector(state_delta) - get_vector(state))/delta_theta)
    #@test (get_vector(state_delta) - get_vector(state))/delta_theta ≈ [im/2*exp(im*thetas[2]/2)*cos(thetas[1]/2), im/2*exp(-im*thetas[2]/2)*sin(thetas[1]/2)] atol=0.001

    #<state|∂2> 1.∂2 2.<state|∂2>
    #overlap_state_partial =  dot(get_vector(state), state_partial2)
    #println("overlap_state_partial=", overlap_state_partial)
    #@test overlap_state_partial ≈ im/2*cos(thetas[1]) atol=0.001

    #<state|∂2> ver2 => |∂2> = (|state_delta> - |state>) /delta_theta => (<state|state_delta> - <state|state>) /delta_theta
    #get_vetcor ver
    #overlap_delta_state = dot(get_vector(state), get_vector(state_delta))
    #overlap_state_state = dot(get_vector(state), get_vector(state))
    #overlap_state_partial_ver2 = (overlap_delta_state - overlap_state_state)/delta_theta
    #println("overlap_state_partial_ver2=", overlap_state_partial_ver2)
    #@test overlap_state_partial_ver2 ≈   im/2*cos(thetas[1]) atol=0.001

    #<state|∂2> ver2 => |∂2> = (|state_delta> - |state>) /delta_theta => (<state|state_delta> - <state|state>) /delta_theta
    #qulacs overlap
    #overlap_delta_state2 = inner_product(state, state_delta)
    #println("overlap_delta=", overlap_delta_state2) #0.9999998750000025 + 0.00040450848033295184im
    #overlap_state_state2=  inner_product(state, state)
    #println("overlap_state=", overlap_state_state2)
    #overlap_state_partial_ver3 = (overlap_delta_state2 - overlap_state_state2)/delta_theta
    #println("overlap_state_partial_ver3=", overlap_state_partial_ver3)
    #@test overlap_state_partial_ver3 ≈   im/2*cos(thetas[1]) atol=0.001

    ####################################
    #以下が問題のコード overlapらへん
    #おそらく、thetas[2]にdelta_thetaがうまく足されていないのが原因ではないか？
    #################################

    thetas_right = copy(thetas)
    println("before_thetas_right=",thetas_right)
    thetas_right[2] +=delta_theta
    println("after_thetas_right=",thetas_right)
    println("thetas_right=", thetas_right) #thetas_right=[0.6283185307179586, 0.9434777960769379]
    println("thetas=", thetas)                   #thetas=[0.6283185307179586, 0.9424777960769379]

    ##########################################
    #delta_thetaはきちんと足されていることは確認できた。
    ###########################################

    #######################3
    #state_left_delta=(0.436165  + im*0.731763,
                        #-0.466638 + im*0.237764))
    #tate_right_delta= (0.435799 +   im * 0.731981)
                       #(-0.466519 + im * 0.237997)

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

    #####################
    #ここでミスが発見！！
    ####################
    overlap_delta = overlap(c, state, thetas, thetas_right)
    println("overlap_delta=", overlap_delta) #0.9999998750000023 + 0.00022571683411611032im
    overlap_state = overlap(c, state, thetas, thetas)
    #println("overlap_state=",overlap_state)
    res = (
        (overlap(c, state, thetas, thetas_right)
       - overlap(c, state, thetas, thetas))
       /delta_theta)
    #res_state_state = (
    #    (overlap(c, state, thetas, thetas)
    #   - overlap(c, state, thetas, thetas))
    #   )
    #println("res_state_state=,", res_state_state)
    @test res≈ im*cos(thetas[1])/2 atol=0.001
end

#B_12 = real(dot(∂2, state)*dot(state, ∂2))
@testset "overlap3_2" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [pi/4, pi/6]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RZ_gate!(c_, 1, thetas[2])
    c = QulacsVariationalQuantumCircuit(c_)
    #update_quantum_state!(c, state)
    #@test get_vector(state) ≈ [exp(im*thetas[2]/2)*cos(thetas[1]/2), -exp(-im*thetas[2]/2)*sin(thetas[1]/2)] #state_vector is correct!!
    thetas_left = copy(thetas)
    delta_theta = 1e-5
    thetas_left[2] += delta_theta
    thetas_right = copy(thetas)
    thetas_right[2] +=delta_theta
    res = real((
                (overlap(c, state, thetas_left, thetas)
                - overlap(c, state, thetas, thetas))
                *(overlap(c, state, thetas, thetas_right)
                - overlap(c, state, thetas, thetas))
            )/delta_theta^2)
    @test res ≈ (cos(thetas[1]))^2/4 atol=0.001
end

@testset "vqs.A" begin
    n_qubit = 2
    n_electron = 1

    c = UCCQuantumCircuit(n_qubit)
    # a_1^dagger a_2 - a^2 a_1 -> 0.5i (X1 Y2 - X2 Y1)
    generator = gen_t1(1, 2)
    add_parametric_circuit_using_generator!(c, generator, 0.0)

    state0 = QulacsQuantumState(n_qubit, 0b01)

    new_thetas = copy(c.thetas)

    new_thetas .+= 1e-8
    c_new = copy(c)
    update_circuit_param!(c_new, new_thetas)
    A = compute_A(c, state0, 1e-2)
    @test A ≈ [0.25] atol=1e-5
end

@testset "vqs.C" begin
    n_qubit = 2
    n_electron = 1

    ham = OFQubitOperator("X1 X2", 1.0)

    c = UCCQuantumCircuit(n_qubit)
    # a_1^dagger a_2 - a^2 a_1 -> 0.5i (X1 Y2 - X2 Y1)
    generator = gen_t1(1, 2)
    add_parametric_circuit_using_generator!(c, generator, 0.0)

    state0 = QulacsQuantumState(n_qubit, 0b01)

    C = compute_C(ham, c, state0, 1e-2)
    @test C ≈ [0.5] atol=1e-5
end

@testset "vqs.thetadot" begin
    n_qubit = 2
    n_electron = 1
    ham = OFQubitOperator("X1 X2", 1.0)
    c = UCCQuantumCircuit(n_qubit)

    # a_1^dagger a_2 - a^2 a_1 -> 0.5i (X1 Y2 - X2 Y1)
    generator = gen_t1(1, 2)
    add_parametric_circuit_using_generator!(c, generator, 0.0)

    state0 = QulacsQuantumState(n_qubit,0b01)
    thetadot = compute_thetadot(ham, c, state0, 1e-2)
    @test thetadot ≈ [2] atol=1e-5
end

@testset "vqs.imag_time_evolve_single_qubit" begin
    """
    Hamiltonian H = Z1
    No time evolution occurs.
    exp(- tau * H) |1> = A |1>,
    where A is a normalization constant.
    """
    n_qubit = 1
    d_theta = 0.01

    ham = OFQubitOperator("Z1", 1.0)
    state0 = QulacsQuantumState(n_qubit, 0b1)

    c = QulacsParametricQuantumCircuit(n_qubit)
    add_parametric_RY_gate!(c, 1, 0.0)
    vc = QulacsVariationalQuantumCircuit(c)

    taus = [0.0, 1.0]
    thetas_tau = imag_time_evolve(ham, vc, state0, taus, d_theta)[1]

    for i in eachindex(taus)
        @test all(thetas_tau[i] .≈ [0.0])
    end
end

@testset "vqs.imag_time_evolve_single_qubit2" begin
    """
    Hamiltonian H = Z1
    No time evolution occurs.
    exp(- tau * H) (|0> + |1>)/sqrt(2) = R_Y(theta(tau)) |0>,
    where A is a normalization constant and
    theta(tau) = acos(exp(-tau)/sqrt(exp(-2tau) + exp(2*tau))).
    theta = pi/2 (tau=0) => pi (tau=+infty)
    """
    n_qubit = 1
    d_theta = 0.01

    ham = OFQubitOperator("Z1", 1.0)
    state0 = QulacsQuantumState(n_qubit, 0b0)

    c = QulacsParametricQuantumCircuit(n_qubit)
    add_parametric_RY_gate!(c, 1, 0.5*pi)
    vc = QulacsVariationalQuantumCircuit(c)

    taus = collect(range(0.0, 1, length=100))
    thetas_tau = imag_time_evolve(ham, vc, state0, taus, d_theta, verbose=true)[1]

    theta_extact(τ) = 2*acos(exp(-τ)/sqrt(exp(-2*τ) + exp(2*τ)))

    thetas_ref = theta_extact.(taus)
    thetas_res = [thetas_tau[i][1] for i in eachindex(taus)]
    @test isapprox(thetas_ref, thetas_res, rtol=0.01)
end

@testset "vqs.imag_time_evolve_two_qubits" begin
    """
    1site hubbard under magnetic field
    Hamiltonian H = Un_(1up)n_(2dn) - hn_(1up) + hn_(2dn)
    No time evolution occurs.
    exp(- tau * H) (|1up> + |1dn>)/sqrt(2) = exp(theta/2*(c_1up^(dag)c_1dn-c_1dn^(dag)c_1up)) |1up>,

    """
    n_qubit = 2
    d_theta = 0.01

    U = 1.0
    h = 1.0

    ham = FermionOperator()
    up = up_index(1)
    down = down_index(1)
    ham += FermionOperator("$(up)^ $(down)^ $(up) $(down)", -U)
    ham += FermionOperator("$(up)^  $(up) ", -h)
    ham += FermionOperator("$(down)^ $(down)", h)
    ham = jordan_wigner(ham)

    state0 = QulacsQuantumState(n_qubit, 0b01)

    c = QulacsParametricQuantumCircuit(n_qubit)
    # add_parametric_RYY_gate! -> add_parametric_multi_Pauli_rotation_gate!
    target = [1,2]
    #X->pauli_X, Y->pauli_Y, Z->pauli_Z
    pauli_ids = [pauli_Y, pauli_Y]
    add_parametric_multi_Pauli_rotation_gate!(c, target, pauli_ids, 0.5*pi)
    #QulacsVariationalQuantumCircuitに対してadd_S_gate!を定義する。
    #QulacsVariationalQuantumCircuitとQulacsQuantumCircuitは両方ともQuantumCircuitを親に持つが、兄弟間同士で継承(多重継承)できないことに注意する。
    add_S_gate!(c, 2)
    vc = QulacsVariationalQuantumCircuit(c)

    taus = collect(range(0.0, 1, length=100))
    thetas_tau = imag_time_evolve(ham, vc, state0, taus, d_theta)[1]

    theta_extact(τ) = 2*acos(exp(h*τ)/sqrt(exp(2*h*τ) + exp(-2*h*τ)))

    #generate vector whose elements are theta_exact(τ)
    thetas_ref = theta_extact.(taus)
    #i is tau and num of parameter is 1
    thetas_res = [thetas_tau[i][1] for i in eachindex(taus)]
    @test isapprox(thetas_ref, thetas_res, rtol=0.01)
end

#compute_gtaを使った１サイトハバードのグリーン関数の一部の計算（一つの基底状態のみの計算）
##ノルムを考慮した
@testset "vqs.compute_gtau_norm_1site_hubbard" begin
    """
    """
    n_qubit　= 2
    nsite = 1
    t = 0.0
    U = 2.0
    µ = 1.0
    up = up_index(1)
    down = down_index(1)
    d_theta = 0.01

    ham_op = generate_ham_1d_hubbard(t, U, nsite, μ)
    ham_op = jordan_wigner(ham_op)

    c_op = FermionOperator("$(down) ", 1.0)
    c_op = jordan_wigner(c_op)
    cdagg_op = FermionOperator("$(down)^ ", 1.0)
    cdagg_op = jordan_wigner(cdagg_op)

    c_ = QulacsParametricQuantumCircuit(n_qubit)
    #add_X_gate!(c_, 1)
    #add_X_gate!(c_, 2)
    target = [1,2]
    pauli_ids = [pauli_Y, pauli_Y]
    add_parametric_multi_Pauli_rotation_gate!(c_, target, pauli_ids, 0.01*pi)
    vc = QulacsVariationalQuantumCircuit(c_)

    state_gs = QulacsQuantumState(n_qubit,0b01)　#state_gs is the ground state after VQE
    state0_ex = QulacsQuantumState(n_qubit,0b11)

    taus = collect(range(0.0, 1, length=5))
    beta = taus[end]

    Gfunc_ij_list, norm = compute_gtau(ham_op, c_op, cdagg_op, vc,  state_gs, state0_ex, taus, d_theta)
    Gfunc_ij_list_exact(τ) = -exp(-U * τ + µ * τ)
    Gfunc_ij_list_ref = Gfunc_ij_list_exact.(taus)
    @test isapprox(Gfunc_ij_list_ref, Gfunc_ij_list, rtol=0.1)
end

#compute_gtaを使った2サイトハバード(U=0,かつ1サイトのポテンシャルあり)のグリーン関数の一部の計算（一つの基底状態のみの計算）
@testset "vqs.compute_gtau_2site_hubbard_U=0_tau_plus" begin
    """
    """
    nsite = 2
    n_qubit　= 2*nsite
    t = 1.0
    ε = 1.0
    up1 = up_index(1)
    dn1 = down_index(1)
    up2 = up_index(2)
    dn2 = down_index(2)
    d_theta = 1e-5

    #Hamiltonian
    ham_op = FermionOperator()
    ham_op += FermionOperator("$(up1)^ $(up2)", -t)
    ham_op += FermionOperator("$(up2)^ $(up1)", -t)
    ham_op += FermionOperator("$(dn1)^ $(dn2)", -t)
    ham_op += FermionOperator("$(dn2)^ $(dn1)", -t)

    ham_op += FermionOperator("$(up1)^ $(up1)", -ε)
    ham_op += FermionOperator("$(dn1)^ $(dn1)", -ε)

    n_electron = 2
    @assert mod(n_electron, 2) == 0
    sparse_mat = get_number_preserving_sparse_operator(ham_op, n_qubit, n_electron);
    enes_ed = eigvals(sparse_mat.toarray());
    ham_op = jordan_wigner(ham_op)

    #vc = QulacsVariationalQuantumCircuit(c)
    vc = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=false, Doubles=false, uccgsd=false, p_uccgsd=false)

    #Perform VQE
    function cost(theta_list)
        state0_gs = create_hf_state(n_qubit, n_electron)
        #state0_gs = QulacsQuantumState(n_qubit,0b0000)
        update_circuit_param!(vc, theta_list)
        update_quantum_state!(vc, state0_gs)
        get_expectation_value(ham_op, state0_gs)
    end

    theta_init = rand(num_theta(vc))
    cost_history = Float64[]
    init_theta_list = theta_init
    push!(cost_history, cost(init_theta_list))

    method = "BFGS"
    options = Dict("disp" => true, "maxiter" =>200, "gtol" =>1e-5)
    callback(x) = push!(cost_history, cost(x))
    Random.seed!(1)
    scipy_opt = pyimport("scipy.optimize")
    opt = scipy_opt.minimize(cost, init_theta_list, method=method, callback=callback)

    Eigval_min = minimum(enes_ed)

    #c^{dag},c
    left_op = FermionOperator("$(up1) ", 1.0)
    left_op = jordan_wigner(left_op)
    right_op = FermionOperator("$(up1)^ ", 1.0)
    right_op = jordan_wigner(right_op)

    ##Ansatz -> apply_qubit_op & imag_time_evolve
    #c_ex = QulacsParametricQuantumCircuit(n_qubit)

    #add_X_gate!(c_ex, 1)
    #add_X_gate!(c_ex, 2)
    #add_X_gate!(c_ex, 3)
    #target = [2,4]
    #pauli_ids = [pauli_Y, pauli_Y]
    #add_parametric_multi_Pauli_rotation_gate!(c_ex, target, pauli_ids, 0.2*pi)
    #add_S_gate!(c_ex, 4)

    #vc_ex = QulacsVariationalQuantumCircuit(c_ex)
    # we use uccgsd
    vc_ex = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=false, Doubles=false)
    #state_gs = QulacsQuantumState(n_qubit,0b0000)
    state_gs = create_hf_state(n_qubit, n_electron)
    update_quantum_state!(vc, state_gs)
    E_gs_debug = get_expectation_value(ham_op, state_gs)
    norm_gs = inner_product(state_gs, state_gs)
    n_electron_ex = 3
    state0_ex = create_hf_state(n_qubit,n_electron_ex)

    taus = collect(range(0.0, 0.02, length=4))
    beta = taus[end]

    k = (2 * t)/(ε + (ε^2 + 4 * t^2)^0.5)
    s = (2 * t)/(ε - (ε^2 + 4 * t^2)^0.5)
    D = (1 + k^2)^0.5
    E = (1 + s^2)^0.5
    ε_1 = (-ε - (ε^2 + 4*t^2)^0.5) / 2
    ε_2 = (-ε + (ε^2 + 4*t^2)^0.5) / 2
    coef_1 = (k - s) / (E^2 * D)
    #coef_3=0(s*k = -1)
    #coef_3 = (s - k) * (s*k + 1) / (E^3 * D^2)

    E_G = 2*ε_1

    Gfunc_ij_exact(τ) = -exp(τ * (-ε_2)) * coef_1^2
    Gfunc_ij_list_ref = Gfunc_ij_exact.(taus)

    Gfunc_ij_list, norm = compute_gtau(ham_op, left_op, right_op, vc_ex,  state_gs, state0_ex, taus, d_theta)
    @test isapprox(Gfunc_ij_list_ref, Gfunc_ij_list, rtol=0.01)
end

@testset "vqs.compute_gtau_2site_hubbard_U=0_tau_minus" begin
    """
    """
    nsite = 2
    n_qubit　= 2*nsite
    t = 1.0
    ε = 1.0
    up1 = up_index(1)
    dn1 = down_index(1)
    up2 = up_index(2)
    dn2 = down_index(2)
    d_theta = 1e-5

    #Hamiltonian
    ham_op = FermionOperator()
    ham_op += FermionOperator("$(up1)^ $(up2)", -t)
    ham_op += FermionOperator("$(up2)^ $(up1)", -t)
    ham_op += FermionOperator("$(dn1)^ $(dn2)", -t)
    ham_op += FermionOperator("$(dn2)^ $(dn1)", -t)

    ham_op += FermionOperator("$(up1)^ $(up1)", -ε)
    ham_op += FermionOperator("$(dn1)^ $(dn1)", -ε)

    n_electron_gs = 2
    @assert mod(n_electron_gs, 2) == 0
    sparse_mat = get_number_preserving_sparse_operator(ham_op, n_qubit, n_electron_gs);
    enes_ed = eigvals(sparse_mat.toarray());
    ham_op = jordan_wigner(ham_op)

    #vc = QulacsVariationalQuantumCircuit(c)
    vc = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=false, Doubles=false, uccgsd=false, p_uccgsd=false)

    #Perform VQE
    function cost(theta_list)
        state0_gs = create_hf_state(n_qubit, n_electron_gs)
        #state0_gs = QulacsQuantumState(n_qubit,0b0000)
        update_circuit_param!(vc, theta_list)
        update_quantum_state!(vc, state0_gs)
        get_expectation_value(ham_op, state0_gs)
    end

    theta_init = rand(num_theta(vc))
    cost_history = Float64[]
    init_theta_list = theta_init
    push!(cost_history, cost(init_theta_list))

    method = "BFGS"
    options = Dict("disp" => true, "maxiter" =>200, "gtol" =>1e-5)
    callback(x) = push!(cost_history, cost(x))
    Random.seed!(1)
    scipy_opt = pyimport("scipy.optimize")
    opt = scipy_opt.minimize(cost, init_theta_list, method=method, callback=callback)

    Eigval_min = minimum(enes_ed)

    #c^{dag},c
    right_op = FermionOperator("$(up1) ", 1.0)
    right_op = jordan_wigner(right_op)
    left_op = FermionOperator("$(up1)^ ", 1.0)
    left_op = jordan_wigner(left_op)

    vc_ex = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=false, Doubles=false)
    update_circuit_param!(vc_ex, rand(num_theta(vc_ex)))

    #state_gs = QulacsQuantumState(n_qubit,0b0000)
    state_gs = create_hf_state(n_qubit, n_electron_gs)
    update_quantum_state!(vc, state_gs)

    #delete?
    E_gs_vqe = get_expectation_value(ham_op, state_gs)
    #norm_gs = inner_product(state_gs, state_gs)
    n_electron_ex = 1
    state0_ex = create_hf_state(n_qubit,n_electron_ex)

    #@assert mod(n_electron_ex, 1) == 0
    taus = collect(range(0.0, 0.02, length=4))
    beta = taus[end]

    k = (2 * t)/(ε + (ε^2 + 4 * t^2)^0.5)
    s = (2 * t)/(ε - (ε^2 + 4 * t^2)^0.5)
    D = (1 + k^2)^0.5
    E = (1 + s^2)^0.5
    ε_1 = (-ε - (ε^2 + 4*t^2)^0.5) / 2
    ε_2 = (-ε + (ε^2 + 4*t^2)^0.5) / 2
    E_G = 2*ε_1
    E_1 = 2*ε_1 + ε

    Gfunc_ij_exact(τ) = exp(τ * (2*ε_1 - ε_1)) * s^2 / E^2
    #Gfunc_ij_exact(τ) = exp(2 * ε_1 *τ ) *(s^2* exp(ε*τ)+1) * s^2 / E^4
    Gfunc_ij_list_ref = Gfunc_ij_exact.(taus)
    #println("Gfunc_ij_list_ref=", Gfunc_ij_list_ref)

    Gfunc_ij_list, norm = compute_gtau(ham_op, left_op, right_op, vc_ex,  state_gs, state0_ex, taus, d_theta)
    #println("Gfunc_ij_list=", Gfunc_ij_list)
    Gfunc_ij_list *= -1
    @test isapprox(Gfunc_ij_list_ref, Gfunc_ij_list, rtol=0.01)
end

@testset "vqs.fubini_Paper_simple" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    println(get_vector(state))
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [0.2pi, 0.3pi]
    #add_parametric_RZ_gate!(c_, 1, -thetas[2])
    add_parametric_RY_gate!(c_, 1, thetas[1])
    #add_parametric_RZ_gate!(c_, 1, thetas[2])
    #add_Z_gate!(c_,1)
    c = QulacsVariationalQuantumCircuit(c_)
    update_quantum_state!(c, state)
    println(get_vector(state))
    @test get_vector(state) ≈ [cos(thetas[1]/2), -sin(thetas[1]/2)] #"state_vector is NOT correct!!"

    A = compute_A(c, state, 1e-2)
    println("A=", A)
    B = compute_B(c, state, 1e-5)
    println("B=", B)
    B2 = compute_B2(c, state, 1e-8)
    println("B2=", B2)
    B2x = compute_B2x(c, state, 1e-8)
    println("B2=", B2x)
    F = compute_F(c, state,1e-3)
    println("F=", F)
    F_exact = [1/4]
    @test isapprox(F, F_exact, rtol=0.01)
end

#@testset "vqs.fubini_Paper" begin
#    n_qubit = 1
#    state  = QulacsQuantumState(n_qubit)
#    println(get_vector(state))
#    c_ = QulacsParametricQuantumCircuit(n_qubit)
#    thetas = [0.2pi, 0.3pi]
#    add_parametric_RZ_gate!(c_, 1, -thetas[2])
#    add_parametric_RY_gate!(c_, 1, thetas[1])
#    add_parametric_RZ_gate!(c_, 1, thetas[2])
#    add_Z_gate!(c_,1)
#    #add_Z_gate!(c_,1)
#    c = QulacsVariationalQuantumCircuit(c_)
#    update_quantum_state!(c, state)
#    println(get_vector(state))
#    @test get_vector(state) ≈ [cos(thetas[1]/2), exp(-im*thetas[2])*sin(thetas[1]/2)] #"state_vector is NOT correct!!"
#    #without zgate ComplexF64[0.9510565162951536 + 0.0im, -0.18163563200134025 + 0.24999999999999994im]
    #add zgate     ComplexF64[0.9510565162951536 + 0.0im, -0.24999999999999994 - 0.18163563200134025im]

#    A = compute_A(c, state, 1e-2)
#    println("A=", A)
#    B = compute_B(c, state, 1e-5#
#    println("B=", B)
#    B2 = compute_B2(c, state, 1e-8)
#    println("B2=", B2)
#    B2x = compute_B2x(c, state, 1e-8)
#    println("B2=", B2x)
#    F = compute_F(c, state,1e-3)
#    println("F=", F)
#    F_exact = [1/4 0.0;0.0 (sin(2*thetas[1]))^2]
#    @test isapprox(F, F_exact, rtol=0.01)
#end

@testset "vqs.B" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    #println(get_vector(state))
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [0.2pi, 0.3pi]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RZ_gate!(c_, 1, thetas[2])
    c = QulacsVariationalQuantumCircuit(c_)
    #update_quantum_state!(c, state)
    #println(get_vector(state))
    #@test get_vector(state) ≈ [exp(im*thetas[2]/2)*cos(thetas[1]/2), -exp(-im*thetas[2]/2)*sin(thetas[1]/2)] #state_vector is correct!!
    A = compute_A(c, state, 1e-3)
    println("A=", A)
    B = compute_B3(c, state, 1e-3)
    println("B=", B)
    #B_right_22 = 0.404508497187474𝑖
    B_exact = [0.0 0.0;
               0.0 1/4*(cos(thetas[1]))^2]
    @test isapprox(B, B_exact, rtol=0.001)
end

@testset "vqs.fubini" begin
    n_qubit = 1
    state  = QulacsQuantumState(n_qubit)
    println(get_vector(state))
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [0.2pi, 0.3pi]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RZ_gate!(c_, 1, thetas[2])
    c = QulacsVariationalQuantumCircuit(c_)
    #update_quantum_state!(c, state)
    #println(get_vector(state))
    #@test get_vector(state) ≈ [exp(im*thetas[2]/2)*cos(thetas[1]/2), -exp(-im*thetas[2]/2)*sin(thetas[1]/2)] #state_vector is correct!!

    A = compute_A(c, state, 1e-3)
    #println("A=", A)
    B = compute_B3(c, state, 1e-3)
    #println("B=", B)
    #B2 = compute_B2(c, state, 1e-8)
    #println("B2=", B2)
    #B2x = compute_B2x(c, state, 1e-8)
    #println("B2=", B2x)
    #F = compute_F(c, state,1e-3)
    #println("F=", F)
    F_exact = [1/4 0.0;
               0.0 1/4 - 1/4*(cos(thetas[1]))^2]
    @test isapprox(A-B, F_exact, rtol=0.01)
end

#How about 2 qubit?
@testset "vqs.fubini2" begin
    n_qubit = 2
    state  = QulacsQuantumState(n_qubit)
    println(get_vector(state))
    c_ = QulacsParametricQuantumCircuit(n_qubit)
    thetas = [0.2pi, 0.3pi, 0.4pi, 0.5pi]
    add_parametric_RY_gate!(c_, 1, thetas[1])
    add_parametric_RY_gate!(c_, 2, thetas[2])
    add_CNOT_gate!(c_, 1, 2)
    add_parametric_RY_gate!(c_, 1, thetas[3])
    add_parametric_RY_gate!(c_, 2, thetas[4])
    c = QulacsVariationalQuantumCircuit(c_)
    #update_quantum_state!(c, state)
    #println(get_vector(state))
    #@test get_vector(state) ≈ [exp(im*thetas[2]/2)*cos(thetas[1]/2), -exp(-im*thetas[2]/2)*sin(thetas[1]/2)] #state_vector is correct!!
    A = compute_A(c, state, 1e-5)
    println("A=", A)
    B = compute_B3(c, state, 1e-5)
    println("B=", B)
    println("F=A-B=")
    F_exact = [1/4 0.0 -sin(thetas[2])/4 0.0;
               0.0 1/4 0.0 cos(thetas[1])/4;
               -sin(thetas[2])/4 0.0 1/4 (sin(thetas[1]-thetas[2])+sin(thetas[1]+thetas[2]))/8;
               0.0 cos(thetas[1])/4 (sin(thetas[1]-thetas[2])+sin(thetas[1]+thetas[2]))/8 1/4]
    #A=ComplexF64
    #[0.24999780023904347 + 0.0im 0.0 + 0.0im 0.2377653629537235 + 0.0im 0.0 + 0.0im;
    #0.0 + 0.0im 0.2500000206850927 + 0.0im 0.0 + 0.0im 0.0772548691685415 + 0.0im;
    #0.2377653629537235 + 0.0im 0.0 + 0.0im 0.2500000206850927 + 0.0im 0.07347344954666822 + 0.0im;
    #0.0 + 0.0im 0.0772548691685415 + 0.0im 0.07347344954666822 + 0.0im 0.24999780023904347 + 0.0im]

    #[0.25 0.0 0.2377641290737884 0.0;
    #0.0 0.25 0.0 0.07725424859373686;
    #0.2377641290737884 0.0 0.25 0.0;
    #0.0 0.07725424859373686 0.0 0.25]
    @test (A-B)[1,1] ≈ F_exact[1,1] atol=0.001
    #@test (A-B) ≈ complex(F_exact) atol=0.001
    @test isapprox((A-B), F_exact, atol=0.01)
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
using Curves
using Test

@testset "Curves.jl" begin
    @testset "Curves" begin
        include("test_curves.jl")
    end

    @testset "tenors" begin
        include("test_tenors.jl")
    end

    @testset "use-case" begin
        # Test use case in Readme
        # construct zero interest rate curve
        c_zero_base = Curve(["2D", "1w", "1M", "3M", "6M", "12M"], [0.5, 0.7, 0.75, 0.83, 1.1, 1.5])
        # define zero rate shifts (e.g. for stress testing or sensitivities)
        c_shifts = Curve([2, 185, 360], [0.1, -0.1, 0.2])
        # shift curve
        c_shifted = c_zero_base + c_shifts
        # calculate discount factors for the unshifted and shifted curves
        c_base_df=apply((x,y) -> exp(-x*y/100/365), c_zero_base, logy=true)
        c_shifted_df = apply((x,y) -> exp(-x*y/100/365), c_shifted, logy=true)
        # calculate log-returns of discount factors
        log_ret = log(c_shifted_df/c_base_df)
        # apply log returns to the base curve - this should give the shifted curve back
        curve_scenario = *(c_base_df, exp(log_ret), logy=true)
        @test curve_scenario ≈ c_shifted_df
    end

end


text:
module TestAqua

import Aqua
import ThreadsX
using Test

# Default `Aqua.test_all(ThreadsX)` does not work due to ambiguities
# in upstream packages.
Aqua.test_all(ThreadsX; ambiguities = false)

@testset "Method ambiguity" begin
    Aqua.test_ambiguities(ThreadsX)
end

end  # module


text:
import Base.Ordering
import Base.Forward
import Base.Reverse
import DataStructures.eq
import Base.lt
import Base.ForwardOrdering
import Base.ReverseOrdering
import DataStructures.IntSemiToken

struct CaseInsensitive <: Ordering
end

lt(::CaseInsensitive, a, b) = isless(lowercase(a), lowercase(b))
eq(::CaseInsensitive, a, b) = isequal(lowercase(a), lowercase(b))

@noinline my_assert(stmt) = stmt ? nothing : throw(AssertionError("assertion failed"))

function my_primes(N)
    w = Vector{Bool}(undef, N)
    fill!(w,true)
    for k = 2 : N
        if w[k]
            w[2 * k : k : N] .= false
        end
    end
    p = Int[]
    for k = 2 : N
        if w[k]
            push!(p,k)
        end
    end
    p
end

## Function checkcorrectness checks a balanced tree for correctness.

function checkcorrectness(t::DataStructures.BalancedTree23{K,D,Ord},
                          allowdups=false)  where {K,D,Ord <: Ordering}
    dsz = size(t.data, 1)
    tsz = size(t.tree, 1)
    r = t.rootloc
    bfstreenodes = Vector{Int}()
    tdpth = t.depth
    intree = BitSet()
    levstart = Vector{Int}(undef, tdpth)
    push!(bfstreenodes, r)
    levstart[1] = 1
    push!(intree, r)
    for curdepth = 2 : tdpth
        levstart[curdepth] = size(bfstreenodes, 1) + 1
        for l = levstart[curdepth - 1] : levstart[curdepth] - 1
            anc = bfstreenodes[l]
            c1 = t.tree[anc].child1
            if in(c1, intree)
                println("anc = $anc c1 = $c1")
                throw(ErrorException("Tree contains loops 1"))
            end
            push!(bfstreenodes, c1)
            push!(intree, c1)
            c2 = t.tree[anc].child2
            if in(c2, intree)
                throw(ErrorException("Tree contains loops 2"))
            end
            push!(bfstreenodes, c2)
            push!(intree, c2)
            c3 = t.tree[anc].child3
            if c3 > 0
                if in(c3, intree)
                    throw(ErrorException("Tree contains loops 3"))
                end
                push!(bfstreenodes, c3)
                push!(intree, c3)
            end
        end
    end
    bfstreesize = size(bfstreenodes, 1)
    dataused = BitSet()
    minkeys = Vector{K}(undef, bfstreesize)
    maxkeys = Vector{K}(undef, bfstreesize)
    for s = levstart[tdpth] : bfstreesize
        anc = bfstreenodes[s]
        c1 = t.tree[anc].child1
        if s == levstart[tdpth]
            if c1 != 1
                throw(ErrorException("Leftmost data descendant should be node 1"))
            end
        else
            minkeys[s] = t.data[c1].k
        end
        c2 = t.tree[anc].child2
        c3 = t.tree[anc].child3
        lastchild = c3 > 0 ? c3 : c2
        if s == bfstreesize
            if lastchild != 2
                throw(ErrorException("Rightmost data descendant should be node 2"))
            end
        else
            maxkeys[s] = t.data[lastchild].k
        end
        if s > levstart[tdpth] &&
            (lt(t.ord, minkeys[s], maxkeys[s - 1]) ||
             (!lt(t.ord, maxkeys[s-1],minkeys[s]) && !allowdups))
            println("tdpth = ", tdpth, " s = ", s,
                    " maxkeys[s-1] = ", maxkeys[s-1],
                    " minkeys[s] = ", minkeys[s])
            throw(ErrorException("Data nodes out of order"))
        end
        if s < bfstreesize || c3 > 0
            if !eq(t.ord, t.tree[anc].splitkey1, t.data[c2].k)
                throw(ErrorException("Splitkey1 of leaf should match key of 2nd child"))
            end
        end
        if s < bfstreesize && c3 > 0
            if !eq(t.ord, t.tree[anc].splitkey2, t.data[c3].k)
                throw(ErrorException("Splitkey2 of leaf should match key of 1st child"))
            end
        end
        if t.data[c1].parent != anc || t.data[c2].parent != anc ||
            (c3 > 0 && t.data[c3].parent != anc)
            println("c1 = $c1 c2 = $c2 c3 = $c3 anc = $anc")
            println("t.data[c1].parent = $(t.data[c1].parent) t.data[c2].parent = $(t.data[c2].parent)")
            if c3 > 0
                println("t.data[c3].parent = $(t.data[c3].parent)")
            end
            throw(ErrorException("Incorrect parent node for data child"))
        end
        push!(dataused, c1)
        push!(dataused, c2)
        if c3 > 0
            push!(dataused, c3)
        end
    end
    for curdepth = tdpth - 1 : -1 : 1
        cp = levstart[curdepth + 1]
        for s = levstart[curdepth] : levstart[curdepth + 1] - 1
            anc = bfstreenodes[s]
            c1 = t.tree[anc].child1
            my_assert(c1 == bfstreenodes[cp])
            if s > levstart[curdepth]
                mk1 = minkeys[cp]
            end
            cp += 1
            if t.tree[c1].parent != anc
                throw(ErrorException("Parent/child1 links do not match"))
            end
            c2 = t.tree[anc].child2
            my_assert(c2 == bfstreenodes[cp])
            mk2 = minkeys[cp]
            cp += 1
            if t.tree[c2].parent != anc
                throw(ErrorException("Parent/child2 links do not match"))
            end
            c3 = t.tree[anc].child3
            my_assert(s == levstart[curdepth] ||
                    lt(t.ord,mk1,mk2) || (!lt(t.ord,mk2,mk1) && allowdups))
            if c3 > 0
                if t.tree[c3].parent != anc
                    throw(ErrorException("Parent/child3 links do not match"))
                end
                mk3 = minkeys[cp]
                cp += 1
                my_assert(lt(t.ord,mk2, mk3) ||
                        !lt(t.ord,mk3,mk2) && allowdups)
            end
            if s > levstart[curdepth]
                minkeys[s] = mk1
            end
            if !eq(t.ord, t.tree[anc].splitkey1, mk2)
                throw(ErrorException("Minkey2 not equal to minimum key among descendants of child2"))
            end
            if c3 > 0 && !eq(t.ord, t.tree[anc].splitkey2, mk3)
                throw(ErrorException("Minkey3 not equal to minimum key among descendants of child3"))
            end
        end
    end
    freedata = BitSet()
    for i = 1 : size(t.freedatainds,1)
        fdi = t.freedatainds[i]
        if in(fdi, freedata)
            throw(ErrorException("t.freedatainds has repeated element $i"))
        end
        if fdi < 1 || fdi > dsz
            throw(ErrorException("t.freedatainds entry out of range"))
        end
        push!(freedata, fdi)
    end
    if last(t.useddatacells) > dsz
        throw(ErrorException("t.useddatacells has indices larger than t.data size"))
    end
    for i = 1 : dsz
        if (in(i, dataused) && !in(i, t.useddatacells)) ||
            (!in(i,dataused) && in(i, t.useddatacells))
            throw(ErrorException("Mismatch between actual data cells used and useddatacells array"))
        end
        if (in(i, freedata) && in(i, dataused)) ||
            (!in(i,freedata) && !in(i, dataused))
            throw(ErrorException("Mismatch between t.freedatainds and t.useddatacells"))
        end
    end
    freetree = BitSet()
    for i = 1 : size(t.freetreeinds,1)
        tfi = t.freetreeinds[i]
        if in(tfi, freetree)
            throw(ErrorException("Free tree index repeated twice"))
        end
        if tfi < 1 || tfi > tsz
            throw(ErrorException("Free tree index out of range"))
        end
        push!(freetree, tfi)
    end
    for i = 1 : tsz
        if (!in(i, intree) && !in(i, freetree)) ||
            (in(i, intree) && in(i, freetree))
            throw(ErrorException("Mismatch between t.freetreeinds and actual cells used"))
        end
    end
end

function testSortedDictBasic()
    # a few basic tests of SortedDict to start
    m1 = SortedDict((Dict{String,String}()), Forward)
    my_assert(typeof(m1) == SortedDict{String, String, ForwardOrdering})
    kdarray = ["hello", "jello", "alpha", "beta", "fortune", "random",
               "july", "wednesday"]
    checkcorrectness(m1.bt, false)
    for i = 1 : div(size(kdarray,1), 2)
        k = kdarray[i*2-1]
        d = kdarray[i*2]
        #println("- inserting: k = $k d = $d")
        m1[k] = d
        checkcorrectness(m1.bt, false)
        # fulldump(m1.bt)
    end
    i1 = startof(m1)
    count = 0
    while i1 != pastendsemitoken(m1)
        count += 1
        k,d = deref((m1,i1))
        #println("+ reading: k = $k, d = $d")
        i1 = advance((m1,i1))
        checkcorrectness(m1.bt, false)
    end
    my_assert(count == 4)
    true
end

function testSortedDictMethods()
    # test all methods of SortedDict here except loops
    m0 = SortedDict(Dict{Int, Float64}())
    my_assert(typeof(m0) == SortedDict{Int,Float64,ForwardOrdering})
    m1 = SortedDict(Dict(8=>32.0, 12=>33.1, 6=>18.2))
    my_assert(typeof(m1) == SortedDict{Int,Float64,ForwardOrdering})
    m01 = SortedDict(8=>32.0, 12=>33.1, 6=>18.2)
    my_assert(typeof(m01) == SortedDict{Int,Float64,ForwardOrdering})
    m11 = SortedDict((8=>32.0, 12=>33.1, 6=>18.2))
    my_assert(typeof(m11) == SortedDict{Int,Float64,ForwardOrdering})
    m02 = SortedDict{Int,Float64}()
    my_assert(typeof(m02) == SortedDict{Int,Float64,ForwardOrdering})
    m03 = SortedDict([(1,2.0), (3,4.0), (5,6.0)])
    my_assert(typeof(m03) == SortedDict{Int,Float64,ForwardOrdering})
    m04 = SortedDict{Int,Float64}(Pair[1=>1, 2=>2.0])
    my_assert(typeof(m04) == SortedDict{Int,Float64,ForwardOrdering})
    m05 = SortedDict{Int,Float64}(Reverse, Pair[1=>1, 2=>2.0])
    my_assert(typeof(m05) == SortedDict{Int,Float64,ReverseOrdering{ForwardOrdering}})
    m06a = SortedDict(Pair[1=>2.0, 3=>'a'])
    my_assert(typeof(m06a) == SortedDict{Any,Any,ForwardOrdering})
    m06b = SortedDict([(1,2.0), (3,'a')])
    my_assert(typeof(m06b) == SortedDict{Int,Any,ForwardOrdering})
    m07a = SortedDict(Pair[1.0=>2, 2=>3])
    my_assert(typeof(m07a) == SortedDict{Any,Any,ForwardOrdering})
    m07b = SortedDict([(1.0,2), (2,3)])
    my_assert(typeof(m07b) == SortedDict{Real,Int,ForwardOrdering})
    m08a = SortedDict(Pair[1.0=>2, 2=>'a'])
    my_assert(typeof(m08a) == SortedDict{Any,Any,ForwardOrdering})
    m08b = SortedDict([(1.0,2), (2,'a')])
    my_assert(typeof(m08b) == SortedDict{Real,Any,ForwardOrdering})
    m09a = SortedDict(Pair{Int}[1=>2, 3=>'a'])
    my_assert(typeof(m09a) == SortedDict{Int,Any,ForwardOrdering})
    m09b = SortedDict([(1,2), (3,'a')])
    my_assert(typeof(m09a) == SortedDict{Int,Any,ForwardOrdering})

    my_assert(m0 == m02)
    my_assert(isequal(m0, m02))
    my_assert(m1 == m01)
    my_assert(isequal(m1, m01))
    my_assert(m1 == m11)
    my_assert(isequal(m1, m11))

    # Test Exceptions
    @test_throws ArgumentError SortedDict([1,2,3,4])
    @test_throws ArgumentError SortedDict{Int,Int}([1,2,3,4])

    expected = ([6,8,12], [18.2, 32.0, 33.1])
    checkcorrectness(m1.bt, false)
    ii = startof(m1)
    m2 = packdeepcopy(m1)
    m3 = packcopy(m1)
    p = first(m1)
    my_assert(p[1] == 6 && p[2] == 18.2)
    my_assert(in(Pair(8,32.0),m3))
    my_assert(!in(Pair(8,32.1),m3))
    push!(m1, 12 => 33.2)
    checkcorrectness(m1.bt, false)
    my_assert(length(m1) == 3)
    my_assert(m1[12] == 33.2)
    push!(m1, 12 => 33.1)
    my_assert(length(m1) == 3)
    my_assert(m1[12] == 33.1)
    for j = 1 : 3
        my_assert(ii != pastendsemitoken(m1))
        pr = deref((m1,ii))
        my_assert(pr[1] == expected[1][j] && pr[2] == expected[2][j])
        checkcorrectness(m1.bt, false)
        oldii = ii
        ii = advance((m1,ii))
        delete!((m1,oldii))
    end
    checkcorrectness(m1.bt, false)
    checkcorrectness(m2.bt, false)
    my_assert(length(m2) == 3)
    ii = startof(m2)
    for j = 1 : 3
        pr = deref((m2,ii))
        my_assert(pr[1] == expected[1][j] && pr[2] == expected[2][j])
        ii = advance((m2,ii))
    end

    checkcorrectness(m3.bt, false)
    my_assert(length(m3) == 3)
    ii = startof(m3)
    for j = 1 : 3
        pr = deref((m3,ii))
        my_assert(pr[1] == expected[1][j] && pr[2] == expected[2][j])
        ii = advance((m3,ii))
    end

    my_assert(isempty(m1))
    my_assert(length(m1) == 0)
    N = 1000
    for i = N : -1 : 2
        m1[i] = convert(Float64,i) ^ 2
        if i % 200 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    my_assert(!isempty(m1))
    my_assert(length(m1) == N - 1)
    for i = 2 : N
        d = pop!(m1, i)
        my_assert(d == convert(Float64,i)^2)
        if i % 200 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    my_assert(isempty(m1))
    my_assert(length(m1) == 0)
    N = 1000
    for i = 2:N
        m1[i] = convert(Float64,i) ^ 2
    end
    my_assert(!isempty(m1))
    my_assert(length(m1) == N - 1)
    for i = 2 : N
        d = pop!(m1, i, -1.0)
        my_assert(d == convert(Float64,i)^2)
        d2 = pop!(m1, i, -1.0)
        my_assert(d2 == -1.0)
        d3 = pop!(m1, i, nothing)
        my_assert(d3 == nothing)

        if i % 200 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    my_assert(isempty(m1))
    my_assert(length(m1) == 0)
    for i = N : -1 : 2
        m1[i] = convert(Float64,i) ^ 2
        if i % 200 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    ii = lastindex(m1)
    for i = 1 : N - 1
        pr = deref((m1,ii))
        my_assert(pr[1] == N + 1 - i && pr[2] == convert(Float64,pr[1]) ^ 2)
        ii = regress((m1,ii))
    end
    lastprime = 1
    while true
        ii = searchsortedafter(m1, lastprime)
        if ii == pastendsemitoken(m1)
            break
        end
        j = deref_key((m1,ii))
        for k = j * 2 : j : N
            p = findkey(m1, k)
            if p != pastendsemitoken(m1)
                delete!((m1,p))
            end
        end
        lastprime = j
    end
    checkcorrectness(m1.bt, false)
    my_assert(ii == pastendsemitoken(m1))
    my_assert(status((m1,ii)) == 3)
    my_assert(status((m1,SDSemiToken(-1))) == 0)
    t = 0
    u = 0.0
    for pr in m1
        t += pr[1]
        u += pr[2]
    end
    numprimes = length(m1)
    pn = convert(Array{Int64}, my_primes(N))
    my_assert(t == sum(pn))
    my_assert(u == sum(pn.^2))
    ij = lastindex(m1)
    my_assert(deref_key((m1,ij)) == last(pn) &&
       convert(Float64, last(pn)^2) ==  deref_value((m1,ij)))
    m1[6] = 49.0
    my_assert(length(m1) == numprimes + 1)
    my_assert(m1[6] == 49.0)
    b, i6 = insert!(m1, 6, 50.0)
    my_assert(length(m1) == numprimes + 1)
    my_assert(!b)
    p = deref((m1,i6))
    my_assert(p[1] == 6 && p[2] == 50.0)
    m1[i6] = 9.0
    p = deref((m1,i6))
    my_assert(p[1] == 6 && p[2] == 9.0)
    my_assert(m1[i6] == 9.0)
    b2, i7 = insert!(m1, 8, 51.0)
    my_assert(b2)
    my_assert(length(m1) == numprimes + 2)
    p = deref((m1,i7))
    my_assert(p[1] == 8 && p[2] == 51.0)
    delete!((m1,i7))
    z = pop!(m1, 6)
    checkcorrectness(m1.bt, false)
    my_assert(z == 9.0)
    i8 = startof(m1)
    p = deref((m1,i8))
    my_assert(p[1] == 2 && p[2] == 4.0)
    my_assert(i8 != beforestartsemitoken(m1))
    my_assert(status((m1,i8)) == 1)
    i9 = regress((m1,i8))
    my_assert(i9 == beforestartsemitoken(m1))
    my_assert(status((m1,i9)) == 2)
    i10 = findkey(m1, 17)
    i11 = regress((m1,i10))
    my_assert(deref_key((m1,i11)) == 13)
    i12 = searchsortedfirst(m1, 47)
    i13 = searchsortedfirst(m1, 48)
    my_assert(deref_key((m1,i12)) == 47)
    my_assert(deref_key((m1,i13)) == 53)
    i14 = searchsortedafter(m1, 47)
    i15 = searchsortedafter(m1, 48)
    my_assert(deref_key((m1,i14)) == 53)
    my_assert(deref_key((m1,i15)) == 53)
    i16 = searchsortedlast(m1, 47)
    i17 = searchsortedlast(m1, 48)
    my_assert(deref_key((m1,i16)) == 47)
    my_assert(deref_key((m1,i17)) == 47)
    ww = my_primes(N)
    cc = last(m1)
    my_assert(cc[1] == last(ww))
    wwx = first(m1)
    my_assert(wwx[1] == 2)
    tpr = eltype(m1)
    my_assert(tpr == Pair{Int,Float64})
    tpr2 = eltype(typeof(m1))
    my_assert(tpr2 == Pair{Int,Float64})
    kt = keytype(m1)
    my_assert(kt == Int)
    kt2 = keytype(typeof(m1))
    my_assert(kt2 == Int)
    vt = valtype(m1)
    my_assert(vt == Float64)
    vt2 = valtype(typeof(m1))
    my_assert(vt2 == Float64)
    my_assert(ordtype(m1) == ForwardOrdering)
    my_assert(ordtype(typeof(m1)) == ForwardOrdering)

    co = orderobject(m1)
    my_assert(co == Forward)
    my_assert(haskey(m1, 71))
    my_assert(!haskey(m1, 77))
    my_assert(get(m1, 70, nothing) == nothing)
    my_assert(get(m1, 70, -45.2) == -45.2)
    my_assert(get(m1, 83, -45.2) == convert(Float64,83)^2)
    h = get!(m1, 5, 27.0)
    my_assert(h == 25.0)
    h = get!(m1, 6, 27.0)
    my_assert(h == 27.0)
    my_assert(m1[6] == 27.0)
    my_assert(length(m1) == length(ww) + 1)
    pop!(m1, 6)
    my_assert(getkey(m1,7, 9) == 7)
    my_assert(getkey(m1,7.0, 9) == 7)
    my_assert(getkey(m1,12, 9) == 9)
    my_assert(getkey(m1,12, nothing) == nothing)
    delete!(m1, 17)
    my_assert(length(m1) == length(ww) - 1)
    my_assert(deref_key((m1,advance((m1,findkey(m1,13))))) == 19)
    empty!(m1)
    checkcorrectness(m1.bt, false)
    my_assert(isempty(m1))
    c1 = SortedDict(Dict("Eggplants"=>3,
                        "Figs"=>9,
                        "Apples"=>7))
    my_assert(typeof(c1) == SortedDict{String, Int, ForwardOrdering})
    c2 = SortedDict(Dict("Eggplants"=>6,
                        "Honeydews"=>19,
                        "Melons"=>11))
    my_assert(!isequal(c1,c2))
    c3 = merge(c1, c2)
    checkcorrectness(c3.bt, false)
    c4 = SortedDict(Dict("Apples"=>7,
                        "Figs"=>9,
                        "Eggplants"=>6,
                        "Melons"=>11,
                        "Honeydews"=>19))
    my_assert(isequal(c3,c4))
    c5 = SortedDict(Dict("Apples"=>7))
    my_assert(!isequal(c4,c5))
    merge!(c1,c2)
    checkcorrectness(c1.bt, false)
    my_assert(isequal(c3,c1))
    merge!(c3,c3)
    my_assert(isequal(c3,c1))
    checkcorrectness(c3.bt, false)

    # issue #216
    my_assert(DataStructures.isordered(SortedDict{Int, String}))

    # check for get! and get
    dfc =  SortedDict{Int, Vector{Int}}()
    x1 = get!(dfc,1,[1])
    my_assert(x1 == [1])
    my_assert(x1 === dfc[1])
    my_assert(x1 === get!(dfc, 1, [1000]))
    my_assert(x1 === get(dfc, 1, [1000]))

    x2 = get!(()->[2], dfc, 2)
    my_assert(x2 == [2])
    my_assert(x2 === dfc[2])
    my_assert(x2 === get!(()->[1000], dfc, 2))
    my_assert(x2 === get(()->[1000], dfc, 2))

    my_assert([42] == get(()->[42], dfc, 3))
    my_assert(!haskey(dfc, 3))
    my_assert([43] == get(dfc, 4, [43]))
    my_assert(!haskey(dfc, 4))
    true
end

function bitreverse(i)
    zeroi = zero(i)
    onei = one(i)
    twoi = onei + onei
    r = zeroi
    for j = 1 : 32
        r *= twoi
        r += (i & onei)
        i = div(i,twoi)
    end
    r
end

function testSortedDictLoops()
    z = 0x00000000
    T = typeof(z)
    ## Test the loops
    zero1 = zero(z)
    one1 = one(z)
    two1 = one1 + one1
    m1 = SortedDict{T,T}()
    N = 1000
    for l = 1 : N
        lUi = convert(T, l)
        m1[bitreverse(lUi)] = lUi
    end
    count = 0
    for (stok,k,v) in semitokens(inclusive(m1, startof(m1), lastindex(m1)))
        for (stok2,k2,v2) in semitokens(exclusive(m1, startof(m1), pastendsemitoken(m1)))
            c = compare(m1,stok,stok2)
            if c < 0
                my_assert(deref_key((m1,stok)) < deref_key((m1,stok2)))
            elseif c == 0
                my_assert(deref_key((m1,stok)) == deref_key((m1,stok2)))
            else
                my_assert(deref_key((m1,stok)) > deref_key((m1,stok2)))
            end
            count += 1
        end
    end
    my_assert(eltype(semitokens(exclusive(m1, startof(m1), pastendsemitoken(m1)))) ==
           Tuple{IntSemiToken, T, T})
    my_assert(count == N^2)
    N = 1000
    sk = zero1
    sv = zero1
    skhalf = zero1
    svhalf = zero1
    for l = 1 : N
        lUi = convert(T, l)
        brl = bitreverse(lUi)
        sk += brl
        m1[brl] = lUi
        sv += lUi
        if brl < div(N, 2)
            skhalf += brl
            svhalf += lUi
        end
    end
    count = 0
    sk2 = zero1
    sv2 = zero1
    for (k,v) in m1
        sk2 += k
        sv2 += v
        count += 1
    end
    my_assert(count == N)
    my_assert(sk2 == sk)
    my_assert(sv == sv2)
    sk2 = zero1
    for k in keys(m1)
        sk2 += k
    end
    my_assert(eltype(keys(m1)) == T)
    my_assert(sk2 == sk)

    sk2a = zero1

    for k in eachindex(m1)
        sk2a += k
    end
    my_assert(eltype(eachindex(m1)) == T)
    my_assert(sk2a == sk)

    sk2b = zero1
    for st in onlysemitokens(m1)
        sk2b += deref_key((m1,st))
    end
    my_assert(sk2b == sk)

    sv2 = zero1
    for v in values(m1)
        sv2 += v
    end
    my_assert(eltype(values(m1)) == T)

    my_assert(sv == sv2)
    count = 0
    for (st,k) in semitokens(keys(m1))
        my_assert(deref_key((m1,st)) == k)
        count += 1
    end
    my_assert(eltype(semitokens(keys(m1))) == Tuple{IntSemiToken, T})

    my_assert(count == N)
    count = 0
    for (st,v) in semitokens(values(m1))
        my_assert(deref_value((m1,st)) == v)
        count += 1
    end
    my_assert(count == N)
    my_assert(eltype(semitokens(values(m1))) == Tuple{IntSemiToken, T})

    pos1 = searchsortedfirst(m1, div(N,2))
    sk2 = zero1
    for k in keys(exclusive(m1, startof(m1), pos1))
        sk2 += k
    end
    my_assert(sk2 == skhalf)
    my_assert(eltype(keys(exclusive(m1, startof(m1), pos1))) == T)

    sk2a = zero1
    for k in eachindex(exclusive(m1, startof(m1), pos1))
        sk2a += k
    end
    my_assert(sk2a == skhalf)
    my_assert(eltype(eachindex(exclusive(m1, startof(m1), pos1))) == T)

    sv2 = zero1
    for v in values(exclusive(m1, startof(m1), pos1))
        sv2 += v
    end
    my_assert(sv2 == svhalf)
    count = 0
    for (k,v) in exclusive(m1, pastendsemitoken(m1), pastendsemitoken(m1))
        count += 1
    end
    my_assert(eltype(keys(exclusive(m1, startof(m1), pos1))) == T)
    my_assert(count == 0)

    count = 0
    for (k,v) in inclusive(m1, startof(m1), beforestartsemitoken(m1))
        count += 1
    end
    my_assert(count == 0)
    my_assert(eltype(keys(inclusive(m1, startof(m1), beforestartsemitoken(m1)))) == T)

    count = 0
    sk5 = zero1
    for k in eachindex(inclusive(m1, startof(m1), startof(m1)))
        sk5 += k
        count += 1
    end
    my_assert(count == 1 && sk5 == deref_key((m1,startof(m1))))
    my_assert(eltype(eachindex(inclusive(m1, startof(m1), startof(m1)))) == T)

    factors = SortedMultiDict{Int,Int}()
    N = 1000
    len = 0
    sum1 = 0
    sum2 = 0
    for factor = 1 : N
        for multiple = factor : factor : N
            insert!(factors, multiple, factor)
            sum1 += multiple
            sum2 += factor
            len += 1
        end
    end

    sum1a = 0
    sum2a = 0

    for (k,v) in factors
        sum1a += k
        sum2a += v
    end
    my_assert(sum1a == sum1 && sum2a == sum2)

    sum1a = 0
    sum2a = 0
    for st in eachindex(factors)
        (k,v) = deref((factors,st))
        sum1a += k
        sum2a += v
    end
    my_assert(sum1a == sum1 && sum2a == sum2)

    sum1a = 0
    sum2a = 0
    for st in onlysemitokens(factors)
        (k,v) = deref((factors,st))
        sum1a += k
        sum2a += v
    end
    my_assert(sum1a == sum1 && sum2a == sum2)
    my_assert(eltype(onlysemitokens(factors)) == IntSemiToken)

    sum2 = 0
    for (k,v) in inclusive(factors,
                           searchsortedfirst(factors,70),
                           searchsortedlast(factors,70))
        sum2 += v
    end

    my_assert(sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70)
    my_assert(eltype(inclusive(factors,
                           searchsortedfirst(factors,70),
                           searchsortedlast(factors,70))) == Pair{Int,Int})

    sum2 = 0
    for st in eachindex(inclusive(factors,
                                  searchsortedfirst(factors,70),
                                  searchsortedlast(factors,70)))
        v = deref_value((factors,st))
        sum2 += v
    end

    my_assert(sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70)
    my_assert(eltype(eachindex(inclusive(factors,
                                     searchsortedfirst(factors,70),
                                     searchsortedlast(factors,70)))) == IntSemiToken)

    sum3 = 0
    for (k,v) in exclusive(factors,
                           searchsortedfirst(factors,60),
                           searchsortedfirst(factors,61))
        sum3 += v
    end
    my_assert(sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60)
    my_assert(eltype(exclusive(factors,
                           searchsortedfirst(factors,70),
                           searchsortedlast(factors,70))) == Pair{Int,Int})

    sum3 = 0
    for st in eachindex(exclusive(factors,
                                  searchsortedfirst(factors,60),
                                  searchsortedfirst(factors,61)))
        v = deref_value((factors,st))
        sum3 += v
    end
    my_assert(sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60)
    my_assert(eltype(eachindex(exclusive(factors,
                                     searchsortedfirst(factors,70),
                                     searchsortedlast(factors,70)))) == IntSemiToken)

    sum4 = 0
    for k in keys(factors)
        sum4 += k
    end
    my_assert(sum4 == sum1)
    my_assert(eltype(keys(factors)) == Int)

    sum5 = 0
    for v in values(factors)
        sum5 += v
    end

    my_assert(sum5 == sum2a)
    my_assert(eltype(values(factors)) == Int)

    sum2 = 0
    for k in keys(inclusive(factors,
                            searchsortedfirst(factors,70),
                            searchsortedlast(factors,70)))
        sum2 += k
    end
    my_assert(sum2 == 70 * 8)
    my_assert(eltype(keys(inclusive(factors,
                                searchsortedfirst(factors,70),
                                searchsortedlast(factors,70)))) ==  Int)

    sum3 = 0
    for k in keys(exclusive(factors,
                            searchsortedfirst(factors,60),
                            searchsortedfirst(factors,61)))
        sum3 += k
    end
    my_assert(sum3 == 60 * 12)
    my_assert(eltype(keys(exclusive(factors,
                                searchsortedfirst(factors,60),
                                searchsortedfirst(factors,61)))) == Int)

    sum2 = 0
    for v in values(inclusive(factors,
                              searchsortedfirst(factors,70),
                              searchsortedlast(factors,70)))
        sum2 += v
    end
    my_assert(sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70)
    my_assert(eltype(values(inclusive(factors,
                                  searchsortedfirst(factors,60),
                                  searchsortedfirst(factors,61)))) == Int)

    sum3 = 0
    for v in values(exclusive(factors,
                              searchsortedfirst(factors,60),
                              searchsortedfirst(factors,61)))
        sum3 += v
    end
    my_assert(sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60)
    my_assert(eltype(values(exclusive(factors,
                                  searchsortedfirst(factors,60),
                                  searchsortedfirst(factors,61)))) == Int)

    sum1b = 0
    sum2b = 0
    for (st,k,v) in semitokens(factors)
        my_assert(deref_value((factors,st)) == v)
        sum1b += k
        sum2b += v
    end
    my_assert(sum1b == sum1a && sum2b == sum2a)
    my_assert(eltype(semitokens(factors)) == Tuple{IntSemiToken, Int, Int})

    sum2 = 0
    for (st,k,v) in semitokens(inclusive(factors,
                                         searchsortedfirst(factors,70),
                                         searchsortedlast(factors,70)))
        my_assert(deref_value((factors,st)) == v)
        sum2 += v
    end
    my_assert(sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70)
    my_assert(eltype(semitokens(inclusive(factors,
                                         searchsortedfirst(factors,70),
                                         searchsortedlast(factors,70)))) ==
        Tuple{IntSemiToken, Int, Int})

    sum3 = 0
    for (st,k,v) in semitokens(exclusive(factors,
                                         searchsortedfirst(factors,60),
                                         searchsortedfirst(factors,61)))
        my_assert(deref_value((factors,st)) == v)
        sum3 += v
    end
    my_assert(sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60)
    my_assert(eltype(semitokens(exclusive(factors,
                                      searchsortedfirst(factors,60),
                                      searchsortedfirst(factors,61)))) ==
         Tuple{IntSemiToken, Int, Int})

    sum4 = 0
    for (st,k) in semitokens(keys(factors))
        my_assert(deref_key((factors,st)) == k && mod(k,deref_value((factors,st))) == 0)
        sum4 += k
    end
    my_assert(sum4 == sum1)
    my_assert(eltype(semitokens(keys(factors))) == Tuple{IntSemiToken,Int})

    sum5 = 0
    for (st,v) in semitokens(values(factors))
        my_assert(deref_value((factors,st)) == v)
        sum5 += v
    end
    my_assert(sum5 == sum2a)
    my_assert(eltype(semitokens(values(factors))) == Tuple{IntSemiToken,Int})

    sum2 = 0
    for (st,k) in semitokens(keys(inclusive(factors,
                                           searchsortedfirst(factors,70),
                                           searchsortedlast(factors,70))))
        my_assert(deref_key((factors,st)) == k && mod(k,deref_value((factors,st))) == 0)
        sum2 += k
    end
    my_assert(sum2 == 70 * 8)
    my_assert(eltype(semitokens(keys(inclusive(factors,
                                           searchsortedfirst(factors,70),
                                           searchsortedlast(factors,70))))) ==
        Tuple{IntSemiToken, Int})

    sum3 = 0
    for (st,k) in semitokens(keys(inclusive(factors,
                                            searchsortedfirst(factors,60),
                                            searchsortedlast(factors,60))))
        my_assert(deref_key((factors,st)) == k && mod(k,deref_value((factors,st))) == 0)
        sum3 += k
    end
    my_assert(sum3 == 60 * 12)
    my_assert(eltype(semitokens(keys(inclusive(factors,
                                            searchsortedfirst(factors,60),
                                            searchsortedlast(factors,60))))) ==
        Tuple{IntSemiToken, Int})

    sum2 = 0
    for (st,v) in semitokens(values(inclusive(factors,
                                              searchsortedfirst(factors,70),
                                              searchsortedlast(factors,70))))
        my_assert(deref_value((factors,st)) == v)
        sum2 += v
    end
    my_assert(sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70)
    my_assert(eltype(semitokens(values(inclusive(factors,
                                              searchsortedfirst(factors,70),
                                              searchsortedlast(factors,70))))) ==
      Tuple{IntSemiToken, Int})

    sum3 = 0
    for (st,v) in semitokens(values(exclusive(factors,
                                              searchsortedfirst(factors,60),
                                              searchsortedfirst(factors,61))))
        my_assert(deref_value((factors,st)) == v)
        sum3 += v
    end
    my_assert(sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60)
    my_assert(eltype(semitokens(values(exclusive(factors,
                                              searchsortedfirst(factors,60),
                                              searchsortedfirst(factors,61))))) ==
       Tuple{IntSemiToken, Int})

    s = SortedSet([39, 24, 2, 14, 45, 107, 66])
    my_assert(typeof(s) == SortedSet{Int, ForwardOrdering})
    sum1 = 0
    for k in s
        sum1 += k
    end
    my_assert(sum1 == sum([39, 24, 2, 14, 45, 107, 66]))

    sum1 = 0
    for (st,k) in semitokens(s)
        my_assert(deref((s,st)) == k)
        sum1 += k
    end
    my_assert(sum1 == sum([39, 24, 2, 14, 45, 107, 66]))
    my_assert(eltype(semitokens(s)) == Tuple{IntSemiToken, Int})

    sum1 = 0
    for st in onlysemitokens(s)
        k = deref((s,st))
        sum1 += k
    end
    my_assert(sum1 == sum([39, 24, 2, 14, 45, 107, 66]))
    my_assert(eltype(onlysemitokens(s)) == IntSemiToken)

    sum1 = 0
    for st in eachindex(s)
        k = deref((s,st))
        sum1 += k
    end
    my_assert(sum1 == sum([39, 24, 2, 14, 45, 107, 66]))
    my_assert(eltype(eachindex(s)) == IntSemiToken)

    sum2 = 0
    for k in inclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))
        sum2 += k
    end
    my_assert(sum2 == 24 + 39 + 45 + 66)
    my_assert(eltype(inclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))) == Int)

    sum2 = 0
    for st in eachindex(inclusive(s,
                                  searchsortedfirst(s, 24),
                                  searchsortedfirst(s, 66)))
        sum2 += deref((s,st))
    end
    my_assert(sum2 == 24 + 39 + 45 + 66)
    my_assert(eltype(eachindex(inclusive(s,
                                     searchsortedfirst(s, 24),
                                     searchsortedfirst(s, 66)))) == IntSemiToken)

    sum2 = 0
    for (st,k) in semitokens(inclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))
        my_assert(deref((s,st)) == k)
        sum2 += k
    end
    my_assert(sum2 == 24 + 39 + 45 + 66)
    my_assert(eltype(semitokens(inclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))) ==
      Tuple{IntSemiToken, Int})

    sum2 = 0
    for st in onlysemitokens(inclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))
        sum2 += deref((s,st))
    end
    my_assert(sum2 == 24 + 39 + 45 + 66)
    my_assert(eltype(onlysemitokens(inclusive(s,
                                          searchsortedfirst(s, 24),
                                          searchsortedfirst(s, 66)))) == IntSemiToken)

    sum3 = 0
    for k in exclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))
        sum3 += k
    end
    my_assert(sum3 == 24 + 39 + 45)
    my_assert(eltype(exclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))) == Int)

    sum3 = 0
    for (st,k) in semitokens(exclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))
        my_assert(deref((s,st)) == k)
        sum3 += k
    end
    my_assert(sum3 == 24 + 39 + 45)
    my_assert(eltype(semitokens(exclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))) ==
     Tuple{IntSemiToken, Int})

    sum3 = 0
    for st in eachindex(exclusive(s,
                                  searchsortedfirst(s, 24),
                                  searchsortedfirst(s, 66)))
        sum3 += deref((s,st))
    end
    my_assert(sum3 == 24 + 39 + 45)
    my_assert(eltype(eachindex(exclusive(s,
                                     searchsortedfirst(s, 24),
                            searchsortedfirst(s, 66)))) == IntSemiToken)
    true
end

function testSortedDictOrderings()
    ## Test use of alternative orderings in test5
    keylist = ["Apple", "aPPle", "berry", "CHerry", "Dairy", "diary"]
    vallist = [6,9,-4,2,1,8]
    m = SortedDict{String,Int}()
    for j = 1:6
        m[keylist[j]] = vallist[j]
    end
    checkcorrectness(m.bt, false)
    expectedord1 = [1,4,5,2,3,6]
    count = 0
    for p in m
        count += 1
        my_assert(p[1] == keylist[expectedord1[count]] &&
                p[2] == vallist[expectedord1[count]])
    end
    my_assert(count == 6)
    m2 = SortedDict((Dict{String, Int}()), Reverse)
    for j = 1 : 6
        m2[keylist[j]] = vallist[j]
    end
    checkcorrectness(m2.bt, false)
    expectedord2 = [6,3,2,5,4,1]
    count = 0
    for p in m2
        count += 1
        my_assert(p[1] == keylist[expectedord2[count]] &&
                p[2] == vallist[expectedord2[count]])
    end
    my_assert(count == 6)
    m3 = SortedDict((Dict{String, Int}()), CaseInsensitive())
    for j = 1 : 6
        m3[keylist[j]] = vallist[j]
    end
    my_assert("BERRY" in keys(m3))
    my_assert(!("BERRY" in collect(keys(m3))))
    checkcorrectness(m3.bt, false)
    expectedord3 = [2,3,4,5,6]
    count = 0
    for p in m3
        count += 1
        my_assert(p[1] == keylist[expectedord3[count]] &&
                p[2] == vallist[expectedord3[count]])
    end
    my_assert(count == 5)
    m3empty = empty(m3)
    my_assert(eltype(m3empty) == Pair{String, Int} &&
       orderobject(m3empty) == CaseInsensitive() &&
       length(m3empty) == 0 && ordtype(m3empty) == CaseInsensitive &&
       ordtype(typeof(m3empty)) == CaseInsensitive)
    m4 = SortedDict((Dict{String,Int}()), Lt((x,y) -> isless(lowercase(x),lowercase(y))))
    for j = 1 : 6
        m4[keylist[j]] = vallist[j]
    end
    checkcorrectness(m4.bt, false)
    count = 0
    for p in m4
        count += 1
        my_assert(p[1] == keylist[expectedord3[count]] &&
                p[2] == vallist[expectedord3[count]])
    end
    my_assert(count == 5)
    true
end

function testSortedMultiDict()
    test_dict = Dict('a'=>1, 'b'=>2, 'c'=>3)

    # Test all methods of SortedMultiDict except loops
    factors = SortedMultiDict{Int,Int}()
    my_assert(typeof(factors) == SortedMultiDict{Int,Int,ForwardOrdering})
    factors2 = SortedMultiDict(test_dict)
    my_assert(typeof(factors2) == SortedMultiDict{Char,Int,ForwardOrdering})
    factors3 = SortedMultiDict{Char,Int}(test_dict)
    my_assert(typeof(factors3) == SortedMultiDict{Char,Int,ForwardOrdering})
    factors4 = SortedMultiDict{Char,Float32}(Reverse, test_dict)
    my_assert(typeof(factors4) == SortedMultiDict{Char,Float32,ReverseOrdering{ForwardOrdering}})

    test_pair_array = Pair{Char}['a'=>1, 'b'=>2, 'c'=>3]
    factors5 = SortedMultiDict(test_pair_array)
    my_assert(typeof(factors5) == SortedMultiDict{Char,Any,ForwardOrdering})

    #@test factors2 == factors3   # Broken!  TODO: fix me...
    my_assert(isequal(factors2, factors3))

    N = 1000
    checkcorrectness(factors.bt, true)
    len = 0
    for factor = 1 : N
        for multiple = factor : factor : N
            insert!(factors, multiple, factor)
            len += 1
        end
    end
    my_assert(length(factors) == len)
    my_assert(Pair(70,2) in factors)
    my_assert(Pair(70,14) in factors)
    my_assert(!(Pair(70,15) in factors))
    my_assert(!(Pair(N+1,15) in factors))
    my_assert(eltype(factors) == Pair{Int,Int})
    my_assert(eltype(typeof(factors)) == Pair{Int,Int})

    my_assert(keytype(factors) == Int)
    my_assert(keytype(typeof(factors)) == Int)
    my_assert(valtype(factors) == Int)
    my_assert(valtype(typeof(factors)) == Int)
    my_assert(ordtype(factors) == ForwardOrdering)
    my_assert(ordtype(typeof(factors)) == ForwardOrdering)

    push!(factors, 70 => 3)
    my_assert(length(factors) == len+1)
    my_assert(Pair(70,3) in factors)
    i = searchsortedfirst(factors, 70)
    dcount = 0
    for (s,k,v) in semitokens(inclusive(factors, i, lastindex(factors)))
        my_assert(k == 70)
        if v == 3
            delete!((factors,s))
            dcount += 1
            break
        end
    end
    my_assert(dcount == 1)

    my_assert(orderobject(factors) == Forward)
    my_assert(haskey(factors, 60))
    my_assert(!haskey(factors, -1))
    my_assert(60 in keys(factors))
    my_assert(!(-1 in keys(factors)))
    checkcorrectness(factors.bt, true)
    i = startof(factors)
    i = advance((factors,i))
    my_assert(deref((factors,i)) == Pair(2,1))
    my_assert(deref_key((factors,i)) == 2)
    my_assert(deref_value((factors,i)) == 1)
    my_assert(factors[i] == 1)
    factors[i] = 7
    my_assert(deref((factors,i)) == Pair(2,7))
    factors[i] = 1
    i = regress((factors,i))
    i = regress((factors,i))
    my_assert(i == beforestartsemitoken(factors))
    pr = first(factors)
    my_assert(pr == Pair(1,1))
    pr2 = last(factors)
    my_assert(pr2 == Pair(N,N))
    i = searchsortedfirst(factors,77)
    my_assert(deref((factors,i)) == Pair(77,1))
    i = searchsortedlast(factors,77)
    my_assert(deref((factors,i)) == Pair(77,77))
    i = searchsortedafter(factors,77)
    my_assert(deref((factors,i)) == Pair(78,1))
    expected = [1,2,4,5,8,10,16,20,40,80]
    i1,i2 = searchequalrange(factors, 80)
    i = i1
    for e in expected
        my_assert(deref_value((factors,i)) == e)
        i = advance((factors,i))
    end
    my_assert(compare(factors,i,i2) != 0)
    my_assert(compare(factors,regress((factors,i)),i2) == 0)
    my_assert(compare(factors,i,i1) != 0)
    insert!(factors, 80, 6)
    my_assert(length(factors) == len + 1)
    checkcorrectness(factors.bt, true)
    expected1 = deepcopy(expected)
    push!(expected1, 6)
    i1,i2 = searchequalrange(factors, 80)
    i = i1
    for e in expected1
        my_assert(deref_value((factors,i)) == e)
        i = advance((factors,i))
    end
    my_assert(compare(factors,i2,regress((factors,i))) == 0)
    my_assert(compare(factors,i,i1) != 0)
    delete!((factors,i2))
    my_assert(length(factors) == len)
    checkcorrectness(factors.bt, true)
    i1,i2 = searchequalrange(factors, 80)
    i = i1
    for e in expected
        my_assert(deref_value((factors,i)) == e)
        i = advance((factors,i))
    end
    my_assert(compare(factors,regress((factors,i)),i2) == 0)
    my_assert(!isempty(factors))
    empty!(factors)
    checkcorrectness(factors.bt, true)
    my_assert(length(factors) == 0)
    my_assert(isempty(factors))
    i = startof(factors)
    my_assert(i == pastendsemitoken(factors))
    i = lastindex(factors)
    my_assert(i == beforestartsemitoken(factors))
    i1,i2 = searchequalrange(factors, N + 2)
    my_assert(i1 == pastendsemitoken(factors))
    my_assert(i2 == beforestartsemitoken(factors))
    m1 = SortedMultiDict("apples"=>2.0, "apples"=>1.0, "bananas"=>1.5)
    my_assert(typeof(m1) == SortedMultiDict{String, Float64, ForwardOrdering})
    checkcorrectness(m1.bt, true)
    m2 = SortedMultiDict("bananas"=>1.5, "apples"=>2.0, "apples"=>1.0)
    checkcorrectness(m2.bt, true)
    m3 = SortedMultiDict("apples"=>1.0, "apples"=>2.0, "bananas"=>1.5)
    checkcorrectness(m3.bt, true)
    my_assert(isequal(m1,m2))
    my_assert(!isequal(m1,m3))
    my_assert(!isequal(m1, SortedMultiDict("apples"=>2.0)))
    stok = insert!(m2, "cherries", 6.1)
    checkcorrectness(m2.bt, true)
    my_assert(!isequal(m1,m2))
    delete!((m2,stok))
    checkcorrectness(m2.bt, true)
    my_assert(isequal(m1,m2))
    m4 = deepcopy(m2)
    checkcorrectness(m4.bt, true)
    my_assert(isequal(m1,m4))
    m5 = packcopy(m2)
    checkcorrectness(m5.bt, true)
    my_assert(isequal(m1,m5))
    m6 = packdeepcopy(m2)
    checkcorrectness(m6.bt, true)
    my_assert(isequal(m1,m6))

    m1 = SortedMultiDict(zip(["bananas", "apples", "cherries", "cherries", "oranges"],
                             [1.0, 2.0, 3.0, 4.0, 5.0]))
    my_assert(typeof(m1) == SortedMultiDict{String, Float64, ForwardOrdering})
    m2 = SortedMultiDict(zip(["apples", "cherries", "cherries", "bananas", "plums"],
                             [6.0, 7.0, 8.0, 9.0, 10.0]))
    m3 = SortedMultiDict(zip(["apples", "apples", "bananas", "bananas",
                              "cherries", "cherries", "cherries", "cherries",
                              "oranges", "plums"],
                             [2.0, 6.0, 1.0, 9.0, 3.0, 4.0, 7.0, 8.0, 5.0, 10.0]))
    m3empty = empty(m3)
    my_assert((eltype(m3empty) == Pair{String, Float64}) &&
        orderobject(m3empty) == Forward &&
        length(m3empty) == 0)
    m4 = merge(m1, m2)
    my_assert(isequal(m3, m4))
    m5 = merge(m2, m1)
    my_assert(!isequal(m3, m5))
    merge!(m1, m2)
    my_assert(isequal(m1, m3))
    m7 = SortedMultiDict{Int,Int}()
    n1 = 10000
    for k = 1 : n1
        insert!(m7, k, k+1)
    end
    for k = 1 : n1
        insert!(m7, k, k+2)
    end
    for k = 1 : n1
        i1, i2 = searchequalrange(m7, k)
        count = 0
        for (key,v) in inclusive(m7, i1, i2)
            count += 1
            my_assert(key == k)
            my_assert(v == k + count)
        end
        my_assert(count == 2)
        count = 0
        for (key,v) in inclusive(m7, searchequalrange(m7,k))
            count += 1
            my_assert(key == k)
            my_assert(v == k + count)
        end
        my_assert(count == 2)
    end
    # issue #216
    my_assert(DataStructures.isordered(SortedMultiDict{Int, String}))
    # issue #773
    s = SortedMultiDict{Int, Int}()
    insert!(s, 4, 41)
    insert!(s, 3, 31)
    insert!(s, 2, 21)
    insert!(s, 2, 22)
    insert!(s, 2, 23)
    insert!(s, 2, 24)
    insert!(s, 2, 25)
    insert!(s, 2, 26)
    insert!(s, 1, 11)
    insert!(s, 1, 12)
    st1 = insert!(s, 1, 13)
    st2 = insert!(s, 1, 14)
    st3 = insert!(s, 1, 15)
    st4 = insert!(s, 1, 16)
    st5 = insert!(s, 1, 17)
    st6 = insert!(s, 1, 18)
    delete!((s, st6))
    delete!((s, st5))
    delete!((s, st4))
    delete!((s, st3))
    delete!((s, st2))
    delete!((s, st1))
    insert!(s, 1, 19)
    checkcorrectness(s.bt, true)
    true
end

function testSortedSet()
    # Test SortedSet
    N = 1000
    sm = 0.0

    m = SortedSet(Float64[])
    my_assert(typeof(m) == SortedSet{Float64, ForwardOrdering})
    mm = SortedSet{Float64}()
    my_assert(typeof(mm) == SortedSet{Float64, ForwardOrdering})
    #@test m == mm   # Broken!  TODO: Fix me...
    my_assert(isequal(m, mm))

    my_assert(typeof(SortedSet()) == SortedSet{Any, ForwardOrdering})
    my_assert(typeof(SortedSet{Float64}()) == SortedSet{Float64, ForwardOrdering})
    my_assert(typeof(SortedSet(Reverse)) == SortedSet{Any, ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedSet{Float64}(Reverse)) == SortedSet{Float64, ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedSet([1,2,3])) == SortedSet{Int, ForwardOrdering})
    my_assert(typeof(SortedSet{Float32}([1,2,3])) == SortedSet{Float32, ForwardOrdering})
    my_assert(typeof(SortedSet(Reverse, [1,2,3])) == SortedSet{Int, ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedSet{Float32}(Reverse, [1,2,3])) == SortedSet{Float32, ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedSet([1,2,3], Reverse)) == SortedSet{Int, ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedSet{Float32}([1,2,3], Reverse)) == SortedSet{Float32, ReverseOrdering{ForwardOrdering}})

    # @test_throws ArgumentError SortedSet(Reverse, Reverse)
    # @test_throws ArgumentError SortedSet{Int}(Reverse, Reverse)

    smallest = 10.0
    largest = -10.0
    for j = 1 : N
        u = j * exp(1)
        ui = u - floor(u)
        push!(m, ui)
        sm += ui
        smallest = min(smallest,ui)
        largest = max(largest,ui)
    end
    isnew,st = insert!(m, 72.5)
    my_assert(isnew)
    my_assert(deref((m,st)) == 72.5)
    delete!((m,st))
    isnew,st = insert!(m, 73.5)
    my_assert(isnew)
    my_assert(deref((m,st)) == 73.5)
    delete!(m, 73.5)
    checkcorrectness(m.bt, false)
    count = 0
    sm2 = 0.0
    prev = -1.0
    for k in m
        sm2 += k
        count += 1
        my_assert(k >= prev)
    end
    my_assert(abs(sm2 - sm) <= 1e-10)
    my_assert(count == N)
    my_assert(length(m) == N)
    ii2 = searchsortedfirst(m, 0.5)
    i3 = startof(m)
    v = first(m)
    my_assert(v == smallest)
    my_assert(deref((m,i3)) == v)
    i4 = lastindex(m)
    w = last(m)
    my_assert(w == largest)
    my_assert(deref((m,i4)) == w)
    i5 = beforestartsemitoken(m)
    my_assert(advance((m,i5)) == i3)
    i6 = pastendsemitoken(m)
    my_assert(regress((m,i6)) == i4)
    my_assert(advance((m,i5)) != i4)
    my_assert(regress((m,i6)) != i3)
    j1 = searchsortedfirst(m,0.5)
    j2 = searchsortedlast(m,0.5)
    j3 = searchsortedafter(m,0.5)
    my_assert(deref((m,j1)) > 0.5)
    my_assert(deref((m,j2)) < 0.5)
    my_assert(advance((m,j2)) == j1)
    my_assert(j1 == j3)
    k1 = searchsortedfirst(m,smallest)
    k2 = searchsortedlast(m,smallest)
    k3 = searchsortedafter(m,smallest)
    my_assert(deref((m,k1)) == smallest)
    my_assert(deref((m,k2)) == smallest)
    my_assert(deref((m,regress((m,k3)))) == smallest)
    secondsmallest = deref((m,k3))
    sk = searchsortedfirst(m,0.4)
    ek = searchsortedfirst(m,0.6)
    dcount = 0
    for (st,k) in semitokens(exclusive(m,sk,ek))
        delete!((m,st))
        dcount += 1
        if dcount % 20 == 0
            checkcorrectness(m.bt, false)
        end
    end
    newcount = 0
    for k in m
        newcount += 1
        my_assert(k < 0.4 || k > 0.6)
    end
    my_assert(newcount == N - dcount)
    my_assert(smallest in m)
    my_assert(haskey(m,smallest))
    my_assert(!(0.5 in m))
    my_assert(!haskey(m,0.5))
    my_assert(eltype(m) == Float64)
    my_assert(eltype(typeof(m)) == Float64)
    my_assert(keytype(m) == Float64)
    my_assert(keytype(typeof(m)) == Float64)
    my_assert(orderobject(m) == Forward)
    my_assert(ordtype(m) == ForwardOrdering)
    my_assert(ordtype(typeof(m)) == ForwardOrdering)
    pop!(m, smallest)
    checkcorrectness(m.bt, false)
    my_assert(length(m) == N - dcount - 1)
    key1 = pop!(m)
    my_assert(key1 == secondsmallest)
    my_assert(length(m) == N - dcount - 2)
    checkcorrectness(m.bt, false)
    my_assert(!isempty(m))
    empty!(m)
    my_assert(isempty(m))
    m1 = SortedSet(["blue", "orange", "red"])
    my_assert(typeof(m1) == SortedSet{String, ForwardOrdering})
    m2 = SortedSet(["orange", "blue", "red"])
    m3 = SortedSet(["orange", "yellow", "red"])
    m3empty = empty(m3)
    my_assert(typeof(m3empty) == SortedSet{String, ForwardOrdering})
    my_assert(eltype(m3empty) == String &&
       length(m3empty) == 0)
    let m3empty_char = empty(m3, Char)
        @test typeof(m3empty_char) == SortedSet{Char, ForwardOrdering}
        Base.emptymutable(m3) == m3empty
        Base.emptymutable(m3, Char) == m3empty_char
    end
    m3_reverse = SortedSet(DataStructures.FasterReverse(), ["orange", "yellow", "red"])
    let m3empty_reverse_char = empty(m3_reverse, Char)
        @test typeof(m3empty_reverse_char) == SortedSet{Char, DataStructures.FasterReverse}
    end
    my_assert(isequal(m1,m2))
    my_assert(!isequal(m1,m3))
    my_assert(!isequal(m1, SortedSet(["blue"])))
    m4 = packcopy(m3)
    my_assert(typeof(m4) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m3,m4))
    m5 = packdeepcopy(m4)
    my_assert(typeof(m5) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m3,m4))
    m6 = deepcopy(m5)
    my_assert(typeof(m6) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m3,m5))
    checkcorrectness(m1.bt, false)
    checkcorrectness(m2.bt, false)
    checkcorrectness(m3.bt, false)
    checkcorrectness(m4.bt, false)
    checkcorrectness(m5.bt, false)
    checkcorrectness(m5.bt, false)
    m7 = union(m1, ["yellow"])
    my_assert(typeof(m7) == SortedSet{String, ForwardOrdering})
    m8 = union(m3, SortedSet(["blue"]))
    my_assert(typeof(m8) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m7,m8))
    my_assert(!isequal(m1,m8))
    union!(m1, ["yellow"])
    my_assert(isequal(m1,m8))
    m8a = intersect(m8)
    my_assert(typeof(m8a) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m8a,m8))
    m9 = intersect(m8, SortedSet(["yellow", "red", "white"]))
    my_assert(typeof(m9) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m9, SortedSet(["red", "yellow"])))
    m9a = intersect(m8, SortedSet(["yellow", "red", "white"]), m8)
    my_assert(typeof(m9a) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m9a, SortedSet(["red", "yellow"])))
    m10 = symdiff(m8,  SortedSet(["yellow", "red", "white"]))
    my_assert(typeof(m10) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m10, SortedSet(["white", "blue", "orange"])))
    m11 = symdiff(m8, SortedSet(["yellow", "red", "blue", "orange",
                                 "zinc"]))
    my_assert(typeof(m11) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m11, SortedSet(["zinc"])))
    m12 = symdiff(SortedSet(["yellow", "red", "blue", "orange",
                                 "zinc"]), m8)
    my_assert(isequal(m12, SortedSet(["zinc"])))
    my_assert(typeof(m12) == SortedSet{String, ForwardOrdering})
    m13 = setdiff(m8, SortedSet(["yellow", "red", "white"]))
    my_assert(typeof(m13) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m13, SortedSet(["blue", "orange"])))
    m14 = setdiff(m8, SortedSet(["blue"]))
    my_assert(typeof(m14) == SortedSet{String, ForwardOrdering})
    my_assert(isequal(m14, SortedSet(["orange", "yellow", "red"])))
    my_assert(issubset(["yellow", "blue"], m8))
    my_assert(!issubset(["blue", "green"], m8))
    setdiff!(m8, SortedSet(["yellow", "red", "white"]))
    my_assert(isequal(m8, SortedSet(["blue", "orange"])))
    true
end

# test the constructors of SortedDict and SortedMultiDict

function testSortedDictConstructors()
    sd1 = SortedDict("w" => 64, "p" => 12)
    my_assert(typeof(sd1) == SortedDict{String, Int, ForwardOrdering})
    my_assert(length(sd1) == 2 && first(sd1) == ("p"=>12) &&
        last(sd1) == ("w"=>64))
    sd2 = SortedDict(Reverse, "w" => 64, "p" => 12)
    my_assert(typeof(sd2) == SortedDict{String, Int, ReverseOrdering{ForwardOrdering}})
    my_assert(length(sd2) == 2 && last(sd2) == ("p"=>12) &&
        first(sd2) == ("w"=>64))
    sd3 = SortedDict(("w"=>64, "p"=>12))
    my_assert(typeof(sd3) == SortedDict{String, Int, ForwardOrdering})
    my_assert(length(sd3) == 2 && first(sd3) == ("p"=>12) &&
        last(sd3) == ("w"=>64))
    sd4 = SortedDict(("w"=>64, "p"=>12), Reverse)
    my_assert(typeof(sd4) == SortedDict{String, Int, ReverseOrdering{ForwardOrdering}})
    my_assert(length(sd4) == 2 && last(sd4) == ("p"=>12) &&
        first(sd4) == ("w"=>64))

    my_assert(typeof(SortedDict()) == SortedDict{Any,Any,ForwardOrdering})
    my_assert(typeof(SortedDict(CaseInsensitive())) == SortedDict{Any,Any,CaseInsensitive})
    my_assert(typeof(SortedDict{Int,Int}(Reverse)) == SortedDict{Int,Int,ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedDict{Int,Int}(Reverse, 1=>2)) == SortedDict{Int,Int,ReverseOrdering{ForwardOrdering}})
    my_assert(typeof(SortedDict{Int,Int}(1=>2)) == SortedDict{Int,Int,ForwardOrdering})

    # @test_throws ArgumentError SortedDict(Reverse, Reverse)
    true
end

function testSortedMultiDictConstructors()
    sm1 = SortedMultiDict("w" => 64, "p" => 12, "p" => 9)
    my_assert(typeof(sm1) == SortedMultiDict{String, Int, ForwardOrdering})
    my_assert(length(sm1) == 3 && first(sm1) == ("p"=>12) &&
        last(sm1) == ("w"=>64))
    sm2 = SortedMultiDict(Reverse, "w" => 64, "p" => 12, "p" => 9)
    my_assert(typeof(sm2) == SortedMultiDict{String, Int, ReverseOrdering{ForwardOrdering}})
    my_assert(length(sm2) == 3 && last(sm2) == ("p"=>9) &&
        first(sm2) == ("w"=>64))
    sm3 = SortedMultiDict(("w"=>64, "p"=>12, "p"=> 9))
    my_assert(typeof(sm3) == SortedMultiDict{String, Int, ForwardOrdering})
    my_assert(length(sm3) == 3 && first(sm3) == ("p"=>12) &&
        last(sm3) == ("w"=>64))
    sm4 = SortedMultiDict(("w"=> 64, "p"=>12, "p"=>9), Reverse)
    my_assert(typeof(sm4) == SortedMultiDict{String, Int, ReverseOrdering{ForwardOrdering}})
    my_assert(length(sm4) == 3 && last(sm4) == ("p"=>9) &&
        first(sm4) == ("w"=>64))

    my_assert(typeof(SortedMultiDict()) == SortedMultiDict{Any,Any,ForwardOrdering})
    my_assert(typeof(SortedMultiDict(CaseInsensitive())) == SortedMultiDict{Any,Any,CaseInsensitive})
    my_assert(typeof(SortedMultiDict{Int,Int}(Reverse)) == SortedMultiDict{Int,Int,ReverseOrdering{ForwardOrdering}})
    # @test typeof(SortedMultiDict{Int,Int}(Reverse, 1=>2)) == SortedMultiDict{Int,Int,ReverseOrdering{ForwardOrdering}}
    # @test typeof(SortedMultiDict{Int,Int}(1=>2)) == SortedMultiDict{Int,Int,ForwardOrdering}
    true
end

@testset "SortedContainers" begin
    @test testSortedDictBasic()
    @test testSortedDictMethods()
    @test testSortedDictLoops()
    @test testSortedDictOrderings()
    @test testSortedMultiDict()
    @test testSortedSet()
    @test testSortedDictConstructors()
    @test testSortedMultiDictConstructors()

    # test all the errors of sorted containers
    m = SortedDict(Dict("a" => 6, "bb" => 9))
    @test_throws KeyError println(m["b"])
    m2 = SortedDict{String,Int}()
    @test_throws BoundsError println(first(m2))
    @test_throws BoundsError println(last(m2))
    i1 = findkey(m,"a")
    delete!((m,i1))
    i2 = findkey(m,"bb")
    @test_throws BoundsError iterate(inclusive(m,i1,i2))
    @test_throws BoundsError iterate(exclusive(m,i1,i2))
    @test m === delete!(m,"a") # Okay to delete! nonexistent keys
    @test_throws KeyError pop!(m,"a")
    m3 = SortedDict((Dict{String, Int}()), Reverse)
    @test_throws ArgumentError isequal(m2, m3)
    @test_throws BoundsError m[i1]
    @test_throws BoundsError regress((m,beforestartsemitoken(m)))
    @test_throws BoundsError advance((m,pastendsemitoken(m)))
    m1 = SortedMultiDict{Int,Int}()
    @test_throws ArgumentError SortedMultiDict([1,2,3])
    @test_throws ArgumentError SortedMultiDict(Forward, [1,2,3])
    @test_throws ArgumentError SortedMultiDict{Int,Int}([1,2,3])
    @test_throws ArgumentError SortedMultiDict{Int,Int}(Forward, [1,2,3])
    @test_throws ArgumentError SortedMultiDict(Forward, Reverse)
    @test_throws ArgumentError isequal(SortedMultiDict("a"=>1), SortedMultiDict("b"=>1.0))
    @test_throws ArgumentError isequal(SortedMultiDict(["a"=>1],Reverse), SortedMultiDict(["b"=>1]))
    @test_throws MethodError SortedMultiDict{Char,Int}(Forward, ["aa"=>2, "bbb"=>5])
    @test_throws MethodError SortedMultiDict(Forward, [("aa",2)=>2, "bbb"=>5])
    @test_throws BoundsError first(m1)
    @test_throws BoundsError last(m1)

    s = SortedSet([3,5])
    @test s === delete!(s,7) # Okay to delete! nonexistent keys
    @test_throws KeyError pop!(s, 7)
    pop!(s)
    pop!(s)
    @test_throws BoundsError pop!(s)
    @test_throws BoundsError first(s)
    @test_throws BoundsError last(s)
    @test_throws ArgumentError isequal(SortedSet(["a"]), SortedSet([1]))
    @test_throws ArgumentError isequal(SortedSet(["a"]), SortedSet(["b"],Reverse))
    @test_throws ArgumentError (("a",6) in m)
    @test_throws ArgumentError ((2,5) in m1)

    s = SortedSet([10,30,50])
    @test pop!(s,10) == 10
    @test pop!(s,30,-1) == 30
    @test pop!(s,30, nothing) == nothing
    @test pop!(s,50, nothing) == 50
    @test pop!(s,50, nothing) == nothing
    @test isempty(s)

    # Test AbstractSet/AbstractDict interface
    for m in [SortedSet([1,2]), SortedDict(1=>2, 2=>3), SortedMultiDict(1=>2, 1=>3)]
        # copy()
        let m1 = copy(m)
            @test isequal(m1, m)
            @test typeof(m1) === typeof(m)
        end
        let m1 = Base.copymutable(m)
            @test isequal(m1, m)
            @test typeof(m1) === typeof(m)
        end
    end
end


text:
ells = "spdfghiklmnoqrtuvwxyz"

export ells