using Test
#using Arpack: eigs
using ExactDiagonalization.FEDW: NCBasis,bkron,boplus
using QuantumLattices
using LinearAlgebra: eigen,diagm
using ExactDiagonalization.Perturbation: ProjState,GSpinGenerator,jval,hamiltonianeff,SpinHeff,projstate

lattice=Lattice("L2P",[Point(PID(1),(0.0,)),Point(PID(2),(1.0,))])
config=Config{Fock{:f}}(pid->Fock{:f}(atom=pid.site%2, norbital=1, nspin=2, nnambu=2), lattice.pids)
table=Table(config,usualfockindextotuple)
t=Hopping(:t,1.0,1)
U=Hubbard(:U,4.0)
μ=Onsite(:μ,-0.5)
h0=Generator((U,μ),Bonds(lattice),config,half=false,table=table)
h1 =Generator((t,),Bonds(lattice),config,half=false,table=table)

@testset "projstate" begin
    basislowi=NCBasis{UInt}([2],[[(1,)]])
    pjslowi=projstate(expand(h0),basislowi,h0.table,pick=1:2)
    @test isa(pjslowi,ProjState)
    @test pjslowi.vals == [-0.50;-0.50]
    #show( pjslowi.vecs )
    @test pjslowi.basis.table ==  basislowi.table
    pjslowij = bkron(pjslowi,pjslowi)
   
    @test pjslowij.vals == [-1.0;-1.0;-1.0;-1.0]
    @test pjslowij.vecs ==  kron(pjslowi.vecs,pjslowi.vecs)
    #show(pjslowij.basis)

    basis_high_i = NCBasis{UInt}([2],[[(2,)]])
    basis_high_j = NCBasis{UInt}([2],[[(0,)]])

    pjs_high_i = projstate(expand(h0),basis_high_i,h0.table)
    pjs_high_j = projstate(expand(h0),basis_high_j,h0.table)
    pjs_high_ij = bkron(pjs_high_i,pjs_high_j)
    pjs_high_ji = bkron(pjs_high_j,pjs_high_i)
    pjs_high = boplus(pjs_high_ij,pjs_high_ji)

    basis_high = boplus( bkron(basis_high_i,basis_high_j),bkron(basis_high_j,basis_high_i) )
   # show(basis_high)
    pjs_high2=projstate(expand(h0),basis_high,h0.table)
    @test pjs_high.vals ==pjs_high2.vals
    #@test pjs_high.vecs == pjs_high2.vecs
    @test pjs_high.basis.table == pjs_high2.basis.table
    #show(pjs_high.basis)
end


@testset "GSpinGenerator and SpinHeff" begin
    namestri=["σᵢ⁰","σᵢˣ","σᵢʸ","σᵢᶻ"]
    namsi = tuple([Symbol(nm) for nm in namestri]...)
    σ0=[1.0 0.0;0.0 1.0]
    σx=[0.0 1.0; 1.0 0.0]
    σy=[0.0 -im;im 0.0]
    σz=[1.0 0.0; 0.0 -1.0]
    gsgi=GSpinGenerator(namsi,(σ0,σx,σy,σz))
    @test gsgi[:σᵢˣ]==σx
    @test gsgi[3]==(:σᵢʸ,σy)
    @test keys(gsgi) == namsi
    @test eltype(gsgi) == ComplexF64
    @test length(gsgi)  == 4
    @test dimension(gsgi) ==(2,2)
    namestrj=["σⱼ⁰","σⱼˣ","σⱼʸ","σⱼᶻ"]
    namsj = tuple([Symbol(nm) for nm in namestrj]...)
    gsgj = GSpinGenerator(namsj,(σ0,σx,σy,σz))
    gsgij=kron(gsgi,gsgj)
    @test gsgij.names==tuple([Symbol(i,",",j) for i in namestri for j in namestrj]...)
    @test gsgij[Symbol("σᵢˣ,σⱼʸ")]==kron(σx,σy)

    basislowij=NCBasis{UInt}([2,2],[[(1,),(1,)]])
    pjslowij=projstate(expand(h0),basislowij,h0.table,pick=1:4)
    basishighij = NCBasis{UInt}([2,2],[[(2,),(0,)],[(0,),(2,)]])
    pjshighij = projstate(expand(h0),basishighij,h0.table)
    heffterm = hamiltonianeff(pjslowij, pjshighij,expand(h1),h1.table) 
    jex=jval(heffterm,gsgij)
    @test isa(jex,SpinHeff)
    @test jex[1] == diagm([-1.0+0im,-1.0,-1.0,-1.0])
    @test jex.jeff[[1,6,11,16]]==[-1.25,0.25,0.25,0.25]
    @test keys(jex) == (Symbol("σᵢ⁰,σⱼ⁰"), Symbol("σᵢ⁰,σⱼˣ"), Symbol("σᵢ⁰,σⱼʸ"), Symbol("σᵢ⁰,σⱼᶻ"), Symbol("σᵢˣ,σⱼ⁰"), Symbol("σᵢˣ,σⱼˣ"), Symbol("σᵢˣ,σⱼʸ"), Symbol("σᵢˣ,σⱼᶻ"), Symbol("σᵢʸ,σⱼ⁰"), Symbol("σᵢʸ,σⱼˣ"), Symbol("σᵢʸ,σⱼʸ"), Symbol("σᵢʸ,σⱼᶻ"), Symbol("σᵢᶻ,σⱼ⁰"), Symbol("σᵢᶻ,σⱼˣ"), Symbol("σᵢᶻ,σⱼʸ"), Symbol("σᵢᶻ,σⱼᶻ"))
    @test values(jex) == jex.jeff
    @test jex[Symbol("σᵢ⁰,σⱼ⁰")] == -1.25
    #show(jex)
    #===============================================================#
    #second method to get ProjState

    pjslowi2 = projstate(expand(h0),NCBasis{UInt}([2],[[(1,)]]),h0.table,pick=1:2)
    pjshighi2 = projstate(expand(h0),NCBasis{UInt}([2],[[(2,)]]),h0.table)
    pjshighi0 = projstate(expand(h0),NCBasis{UInt}([2],[[(0,)]]),h0.table)

    pjslowj2 = projstate(expand(h0),NCBasis{UInt}([2],[[(1,)]]),h0.table,pick=1:2)
    pjshighj2 = projstate(expand(h0),NCBasis{UInt}([2],[[(2,)]]),h0.table)
    pjshighj0 = projstate(expand(h0),NCBasis{UInt}([2],[[(0,)]]),h0.table)

    pjslowij2 = bkron(pjslowi2,pjslowj2)
    pjshighij2 = boplus(bkron(pjshighi2,pjshighj0),bkron(pjshighi0,pjshighj2))

    heffterm2 = hamiltonianeff(pjslowij2, pjshighij2,expand(h1),h1.table) 
    jex2=jval(heffterm2,gsgij)
    @test jex2[1] == diagm([-1.0+0im,-1.0,-1.0,-1.0])
    @test jex2.jeff[[1,6,11,16]]==[-1.25,0.25,0.25,0.25]
end