using Test
using Quiqbox

@testset "TwoBody.jl tests" begin
    errT = 1e-10
    nucs = ["H", "H"]
    cens = [[-0.7, 0.0, 0.0], [ 0.7, 0.0, 0.0]]
    bs = genBasisFunc.(cens, fill("6-31G", 2)) |> flatten
    bsiz = basisSize(bs) |> sum
    vec = [[1.1401448653048112 0.6114018855091738 0.38206058180946545 0.4478418840074861;
            0.6114018855091738 0.626369765524018 0.34091363711214673 0.4981842896930728;
            0.38206058180946545 0.34091363711214673 0.6748720646269153 0.41372200192631736;
            0.4478418840074861 0.4981842896930728 0.41372200192631736 0.5200082440216406],
           [0.6114018855091738 0.34515151828171575 0.21634301778438714 0.253673516795947;
            0.34515151828171575 0.3803285686123761 0.20500253723707199 0.30420736880783067;
            0.21634301778438714 0.20500253723707199 0.4137220019263174 0.25314595728085076;
            0.253673516795947 0.30420736880783067 0.25314595728085076 0.3221534573785314],
           [0.38206058180946545 0.21634301778438714 0.17120945350064398 0.16981096699715015;
            0.21634301778438714 0.238582175642025 0.16981096699715015 0.20918696249164206;
            0.17120945350064398 0.16981096699715015 0.38206058180946545 0.21634301778438708;
            0.16981096699715015 0.20918696249164206 0.21634301778438708 0.2385821756420249],
           [0.4478418840074861 0.253673516795947 0.16981096699715015 0.19017743024985695;
            0.253673516795947 0.28111489545152707 0.1642082324755541 0.2316959448024114;
            0.16981096699715015 0.1642082324755541 0.34091363711214673 0.20500253723707199;
            0.19017743024985695 0.2316959448024114 0.20500253723707199 0.25160103976844717],
           [0.6114018855091738 0.34515151828171575 0.21634301778438714 0.253673516795947;
            0.34515151828171575 0.3803285686123761 0.20500253723707199 0.30420736880783067;
            0.21634301778438714 0.20500253723707199 0.4137220019263174 0.25314595728085076;
            0.253673516795947 0.30420736880783067 0.25314595728085076 0.3221534573785314],
           [0.626369765524018 0.3803285686123761 0.238582175642025 0.28111489545152707;
            0.3803285686123761 0.4829562753692258 0.25160103976844717 0.3918319845719349;
            0.238582175642025 0.25160103976844717 0.5200082440216407 0.32215345737853146;
            0.28111489545152707 0.3918319845719349 0.32215345737853146 0.43088268760281656],
           [0.34091363711214673 0.20500253723707199 0.16981096699715015 0.1642082324755541;
            0.20500253723707199 0.25160103976844717 0.19017743024985695 0.2316959448024114;
            0.16981096699715015 0.19017743024985695 0.4478418840074862 0.253673516795947;
            0.1642082324755541 0.2316959448024114 0.253673516795947 0.281114895451527],
           [0.4981842896930728 0.30420736880783067 0.20918696249164206 0.2316959448024114;
            0.30420736880783067 0.3918319845719349 0.2316959448024114 0.33726555786356255;
            0.20918696249164206 0.2316959448024114 0.4981842896930728 0.30420736880783067;
            0.2316959448024114 0.33726555786356255 0.30420736880783067 0.3918319845719349],
           [0.38206058180946545 0.21634301778438714 0.17120945350064398 0.16981096699715015;
            0.21634301778438714 0.238582175642025 0.16981096699715015 0.20918696249164206;
            0.17120945350064398 0.16981096699715015 0.38206058180946545 0.21634301778438708;
            0.16981096699715015 0.20918696249164206 0.21634301778438708 0.2385821756420249],
           [0.34091363711214673 0.20500253723707199 0.16981096699715015 0.1642082324755541;
            0.20500253723707199 0.25160103976844717 0.19017743024985695 0.2316959448024114;
            0.16981096699715015 0.19017743024985695 0.4478418840074862 0.253673516795947;
            0.1642082324755541 0.2316959448024114 0.253673516795947 0.281114895451527],
           [0.6748720646269153 0.4137220019263174 0.38206058180946545 0.34091363711214673;
            0.4137220019263174 0.5200082440216407 0.4478418840074862 0.4981842896930728;
            0.38206058180946545 0.4478418840074862 1.1401448653048112 0.6114018855091738;
            0.34091363711214673 0.4981842896930728 0.6114018855091738 0.626369765524018],
           [0.41372200192631736 0.25314595728085076 0.21634301778438708 0.20500253723707199;
            0.25314595728085076 0.32215345737853146 0.253673516795947 0.30420736880783067;
            0.21634301778438708 0.253673516795947 0.6114018855091738 0.34515151828171575;
            0.20500253723707199 0.30420736880783067 0.34515151828171575 0.3803285686123761],
           [0.4478418840074861 0.253673516795947 0.16981096699715015 0.19017743024985695;
            0.253673516795947 0.28111489545152707 0.1642082324755541 0.2316959448024114;
            0.16981096699715015 0.1642082324755541 0.34091363711214673 0.20500253723707199;
            0.19017743024985695 0.2316959448024114 0.20500253723707199 0.25160103976844717],
           [0.4981842896930728 0.30420736880783067 0.20918696249164206 0.2316959448024114;
            0.30420736880783067 0.3918319845719349 0.2316959448024114 0.33726555786356255;
            0.20918696249164206 0.2316959448024114 0.4981842896930728 0.30420736880783067;
            0.2316959448024114 0.33726555786356255 0.30420736880783067 0.3918319845719349],
           [0.41372200192631736 0.25314595728085076 0.21634301778438708 0.20500253723707199;
            0.25314595728085076 0.32215345737853146 0.253673516795947 0.30420736880783067;
            0.21634301778438708 0.253673516795947 0.6114018855091738 0.34515151828171575;
            0.20500253723707199 0.30420736880783067 0.34515151828171575 0.3803285686123761],
           [0.5200082440216406 0.3221534573785314 0.2385821756420249 0.25160103976844717;
            0.3221534573785314 0.43088268760281656 0.281114895451527 0.3918319845719349;
            0.2385821756420249 0.281114895451527 0.626369765524018 0.3803285686123761;
            0.25160103976844717 0.3918319845719349 0.3803285686123761 0.4829562753692258]]
    eeI = reshape(vec |> flatten, 4,4,4,4,1)

    randIdx = rand(1:4, 4)
    randVal = eeI[randIdx...,1] * ones(1,1,1,1,1)
    @test isapprox(eeInteraction(bs[randIdx]...), randVal, atol=errT)

    eeI1, uniqueIdx1 = Quiqbox.eeInteractionsCore(bs, outputUniqueIndices=true)
    isapprox(eeInteractions(bs), eeI1[:,:,:,:,1], atol=1e-15)
    @test isapprox(eeI1, eeI, atol=errT)
    uniqueInts1 = [eeI1[i...] for i in uniqueIdx1]
    @test unique(eeI1) ⊆ uniqueInts1
    @test isapprox(sort([eeI[i,j,k,l,1] for i = 1:size(eeI,1) 
                                        for j = 1:i 
                                        for k = 1:i 
                                        for l = 1:(k==i ? j : k)]), 
                   sort(uniqueInts1), atol=errT)
    
    eeI2, uniqueIdx2 = Quiqbox.eeInteractionsCore(bs[[1,3,2,4]], outputUniqueIndices=true)
    uniqueInts2 = [eeI2[i...] for i in uniqueIdx2]
    @test isapprox(sort(uniqueInts1), sort(uniqueInts2), atol=errT)
end