using Test
using OpenQuantumSystems
using LinearAlgebra
using SparseArrays

# https://github.com/JuliaIO/Suppressor.jl/blob/master/src/Suppressor.jl
macro suppress(block)
    quote
        if ccall(:jl_generating_output, Cint, ()) == 0
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            out_reader = @async read(out_rd, String)

            original_stderr = stderr
            err_rd, err_wr = redirect_stderr()
            err_reader = @async read(err_rd, String)

            # approach adapted from https://github.com/JuliaLang/IJulia.jl/pull/667/files
            logstate = Base.CoreLogging._global_logstate
            logger = logstate.logger
            if logger.stream == original_stderr
                new_logstate =
                    Base.CoreLogging.LogState(typeof(logger)(err_wr, logger.min_level))
                Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), new_logstate))
            end
        end

        try
            $(esc(block))
        finally
            if ccall(:jl_generating_output, Cint, ()) == 0
                redirect_stdout(original_stdout)
                close(out_wr)

                redirect_stderr(original_stderr)
                close(err_wr)

                if logger.stream == stderr
                    Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), logstate))
                end
            end
        end
    end
end

@testset "aggregate" begin

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))


    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.4, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1, mode2], 3, Energy)
    mol2 = Molecule([mode1, mode2], 3, Energy)
    agg = Aggregate([mol1, mol2])

    @test getNvib(agg) == [[3, 3], [3, 3]]
    @test electronicIndices(agg) == [[1, 1], [2, 1], [1, 2]]

    vibInds = [
        [[1, 1], [1, 1]],
        [[1, 1], [1, 2]],
        [[1, 1], [1, 3]],
        [[1, 1], [2, 1]],
        [[1, 1], [2, 2]],
        [[1, 1], [2, 3]],
        [[1, 1], [3, 1]],
        [[1, 1], [3, 2]],
        [[1, 1], [3, 3]],
        [[1, 2], [1, 1]],
        [[1, 2], [1, 2]],
        [[1, 2], [1, 3]],
        [[1, 2], [2, 1]],
        [[1, 2], [2, 2]],
        [[1, 2], [2, 3]],
        [[1, 2], [3, 1]],
        [[1, 2], [3, 2]],
        [[1, 2], [3, 3]],
        [[1, 3], [1, 1]],
        [[1, 3], [1, 2]],
        [[1, 3], [1, 3]],
        [[1, 3], [2, 1]],
        [[1, 3], [2, 2]],
        [[1, 3], [2, 3]],
        [[1, 3], [3, 1]],
        [[1, 3], [3, 2]],
        [[1, 3], [3, 3]],
        [[2, 1], [1, 1]],
        [[2, 1], [1, 2]],
        [[2, 1], [1, 3]],
        [[2, 1], [2, 1]],
        [[2, 1], [2, 2]],
        [[2, 1], [2, 3]],
        [[2, 1], [3, 1]],
        [[2, 1], [3, 2]],
        [[2, 1], [3, 3]],
        [[2, 2], [1, 1]],
        [[2, 2], [1, 2]],
        [[2, 2], [1, 3]],
        [[2, 2], [2, 1]],
        [[2, 2], [2, 2]],
        [[2, 2], [2, 3]],
        [[2, 2], [3, 1]],
        [[2, 2], [3, 2]],
        [[2, 2], [3, 3]],
        [[2, 3], [1, 1]],
        [[2, 3], [1, 2]],
        [[2, 3], [1, 3]],
        [[2, 3], [2, 1]],
        [[2, 3], [2, 2]],
        [[2, 3], [2, 3]],
        [[2, 3], [3, 1]],
        [[2, 3], [3, 2]],
        [[2, 3], [3, 3]],
        [[3, 1], [1, 1]],
        [[3, 1], [1, 2]],
        [[3, 1], [1, 3]],
        [[3, 1], [2, 1]],
        [[3, 1], [2, 2]],
        [[3, 1], [2, 3]],
        [[3, 1], [3, 1]],
        [[3, 1], [3, 2]],
        [[3, 1], [3, 3]],
        [[3, 2], [1, 1]],
        [[3, 2], [1, 2]],
        [[3, 2], [1, 3]],
        [[3, 2], [2, 1]],
        [[3, 2], [2, 2]],
        [[3, 2], [2, 3]],
        [[3, 2], [3, 1]],
        [[3, 2], [3, 2]],
        [[3, 2], [3, 3]],
        [[3, 3], [1, 1]],
        [[3, 3], [1, 2]],
        [[3, 3], [1, 3]],
        [[3, 3], [2, 1]],
        [[3, 3], [2, 2]],
        [[3, 3], [2, 3]],
        [[3, 3], [3, 1]],
        [[3, 3], [3, 2]],
        [[3, 3], [3, 3]],
    ]
    @test vibrationalIndices(agg) == vibInds

    aggInds = [
        [[1, 1], [[1, 1], [1, 1]]],
        [[1, 1], [[1, 1], [1, 2]]],
        [[1, 1], [[1, 1], [1, 3]]],
        [[1, 1], [[1, 1], [2, 1]]],
        [[1, 1], [[1, 1], [2, 2]]],
        [[1, 1], [[1, 1], [2, 3]]],
        [[1, 1], [[1, 1], [3, 1]]],
        [[1, 1], [[1, 1], [3, 2]]],
        [[1, 1], [[1, 1], [3, 3]]],
        [[1, 1], [[1, 2], [1, 1]]],
        [[1, 1], [[1, 2], [1, 2]]],
        [[1, 1], [[1, 2], [1, 3]]],
        [[1, 1], [[1, 2], [2, 1]]],
        [[1, 1], [[1, 2], [2, 2]]],
        [[1, 1], [[1, 2], [2, 3]]],
        [[1, 1], [[1, 2], [3, 1]]],
        [[1, 1], [[1, 2], [3, 2]]],
        [[1, 1], [[1, 2], [3, 3]]],
        [[1, 1], [[1, 3], [1, 1]]],
        [[1, 1], [[1, 3], [1, 2]]],
        [[1, 1], [[1, 3], [1, 3]]],
        [[1, 1], [[1, 3], [2, 1]]],
        [[1, 1], [[1, 3], [2, 2]]],
        [[1, 1], [[1, 3], [2, 3]]],
        [[1, 1], [[1, 3], [3, 1]]],
        [[1, 1], [[1, 3], [3, 2]]],
        [[1, 1], [[1, 3], [3, 3]]],
        [[1, 1], [[2, 1], [1, 1]]],
        [[1, 1], [[2, 1], [1, 2]]],
        [[1, 1], [[2, 1], [1, 3]]],
        [[1, 1], [[2, 1], [2, 1]]],
        [[1, 1], [[2, 1], [2, 2]]],
        [[1, 1], [[2, 1], [2, 3]]],
        [[1, 1], [[2, 1], [3, 1]]],
        [[1, 1], [[2, 1], [3, 2]]],
        [[1, 1], [[2, 1], [3, 3]]],
        [[1, 1], [[2, 2], [1, 1]]],
        [[1, 1], [[2, 2], [1, 2]]],
        [[1, 1], [[2, 2], [1, 3]]],
        [[1, 1], [[2, 2], [2, 1]]],
        [[1, 1], [[2, 2], [2, 2]]],
        [[1, 1], [[2, 2], [2, 3]]],
        [[1, 1], [[2, 2], [3, 1]]],
        [[1, 1], [[2, 2], [3, 2]]],
        [[1, 1], [[2, 2], [3, 3]]],
        [[1, 1], [[2, 3], [1, 1]]],
        [[1, 1], [[2, 3], [1, 2]]],
        [[1, 1], [[2, 3], [1, 3]]],
        [[1, 1], [[2, 3], [2, 1]]],
        [[1, 1], [[2, 3], [2, 2]]],
        [[1, 1], [[2, 3], [2, 3]]],
        [[1, 1], [[2, 3], [3, 1]]],
        [[1, 1], [[2, 3], [3, 2]]],
        [[1, 1], [[2, 3], [3, 3]]],
        [[1, 1], [[3, 1], [1, 1]]],
        [[1, 1], [[3, 1], [1, 2]]],
        [[1, 1], [[3, 1], [1, 3]]],
        [[1, 1], [[3, 1], [2, 1]]],
        [[1, 1], [[3, 1], [2, 2]]],
        [[1, 1], [[3, 1], [2, 3]]],
        [[1, 1], [[3, 1], [3, 1]]],
        [[1, 1], [[3, 1], [3, 2]]],
        [[1, 1], [[3, 1], [3, 3]]],
        [[1, 1], [[3, 2], [1, 1]]],
        [[1, 1], [[3, 2], [1, 2]]],
        [[1, 1], [[3, 2], [1, 3]]],
        [[1, 1], [[3, 2], [2, 1]]],
        [[1, 1], [[3, 2], [2, 2]]],
        [[1, 1], [[3, 2], [2, 3]]],
        [[1, 1], [[3, 2], [3, 1]]],
        [[1, 1], [[3, 2], [3, 2]]],
        [[1, 1], [[3, 2], [3, 3]]],
        [[1, 1], [[3, 3], [1, 1]]],
        [[1, 1], [[3, 3], [1, 2]]],
        [[1, 1], [[3, 3], [1, 3]]],
        [[1, 1], [[3, 3], [2, 1]]],
        [[1, 1], [[3, 3], [2, 2]]],
        [[1, 1], [[3, 3], [2, 3]]],
        [[1, 1], [[3, 3], [3, 1]]],
        [[1, 1], [[3, 3], [3, 2]]],
        [[1, 1], [[3, 3], [3, 3]]],
        [[2, 1], [[1, 1], [1, 1]]],
        [[2, 1], [[1, 1], [1, 2]]],
        [[2, 1], [[1, 1], [1, 3]]],
        [[2, 1], [[1, 1], [2, 1]]],
        [[2, 1], [[1, 1], [2, 2]]],
        [[2, 1], [[1, 1], [2, 3]]],
        [[2, 1], [[1, 1], [3, 1]]],
        [[2, 1], [[1, 1], [3, 2]]],
        [[2, 1], [[1, 1], [3, 3]]],
        [[2, 1], [[1, 2], [1, 1]]],
        [[2, 1], [[1, 2], [1, 2]]],
        [[2, 1], [[1, 2], [1, 3]]],
        [[2, 1], [[1, 2], [2, 1]]],
        [[2, 1], [[1, 2], [2, 2]]],
        [[2, 1], [[1, 2], [2, 3]]],
        [[2, 1], [[1, 2], [3, 1]]],
        [[2, 1], [[1, 2], [3, 2]]],
        [[2, 1], [[1, 2], [3, 3]]],
        [[2, 1], [[1, 3], [1, 1]]],
        [[2, 1], [[1, 3], [1, 2]]],
        [[2, 1], [[1, 3], [1, 3]]],
        [[2, 1], [[1, 3], [2, 1]]],
        [[2, 1], [[1, 3], [2, 2]]],
        [[2, 1], [[1, 3], [2, 3]]],
        [[2, 1], [[1, 3], [3, 1]]],
        [[2, 1], [[1, 3], [3, 2]]],
        [[2, 1], [[1, 3], [3, 3]]],
        [[2, 1], [[2, 1], [1, 1]]],
        [[2, 1], [[2, 1], [1, 2]]],
        [[2, 1], [[2, 1], [1, 3]]],
        [[2, 1], [[2, 1], [2, 1]]],
        [[2, 1], [[2, 1], [2, 2]]],
        [[2, 1], [[2, 1], [2, 3]]],
        [[2, 1], [[2, 1], [3, 1]]],
        [[2, 1], [[2, 1], [3, 2]]],
        [[2, 1], [[2, 1], [3, 3]]],
        [[2, 1], [[2, 2], [1, 1]]],
        [[2, 1], [[2, 2], [1, 2]]],
        [[2, 1], [[2, 2], [1, 3]]],
        [[2, 1], [[2, 2], [2, 1]]],
        [[2, 1], [[2, 2], [2, 2]]],
        [[2, 1], [[2, 2], [2, 3]]],
        [[2, 1], [[2, 2], [3, 1]]],
        [[2, 1], [[2, 2], [3, 2]]],
        [[2, 1], [[2, 2], [3, 3]]],
        [[2, 1], [[2, 3], [1, 1]]],
        [[2, 1], [[2, 3], [1, 2]]],
        [[2, 1], [[2, 3], [1, 3]]],
        [[2, 1], [[2, 3], [2, 1]]],
        [[2, 1], [[2, 3], [2, 2]]],
        [[2, 1], [[2, 3], [2, 3]]],
        [[2, 1], [[2, 3], [3, 1]]],
        [[2, 1], [[2, 3], [3, 2]]],
        [[2, 1], [[2, 3], [3, 3]]],
        [[2, 1], [[3, 1], [1, 1]]],
        [[2, 1], [[3, 1], [1, 2]]],
        [[2, 1], [[3, 1], [1, 3]]],
        [[2, 1], [[3, 1], [2, 1]]],
        [[2, 1], [[3, 1], [2, 2]]],
        [[2, 1], [[3, 1], [2, 3]]],
        [[2, 1], [[3, 1], [3, 1]]],
        [[2, 1], [[3, 1], [3, 2]]],
        [[2, 1], [[3, 1], [3, 3]]],
        [[2, 1], [[3, 2], [1, 1]]],
        [[2, 1], [[3, 2], [1, 2]]],
        [[2, 1], [[3, 2], [1, 3]]],
        [[2, 1], [[3, 2], [2, 1]]],
        [[2, 1], [[3, 2], [2, 2]]],
        [[2, 1], [[3, 2], [2, 3]]],
        [[2, 1], [[3, 2], [3, 1]]],
        [[2, 1], [[3, 2], [3, 2]]],
        [[2, 1], [[3, 2], [3, 3]]],
        [[2, 1], [[3, 3], [1, 1]]],
        [[2, 1], [[3, 3], [1, 2]]],
        [[2, 1], [[3, 3], [1, 3]]],
        [[2, 1], [[3, 3], [2, 1]]],
        [[2, 1], [[3, 3], [2, 2]]],
        [[2, 1], [[3, 3], [2, 3]]],
        [[2, 1], [[3, 3], [3, 1]]],
        [[2, 1], [[3, 3], [3, 2]]],
        [[2, 1], [[3, 3], [3, 3]]],
        [[1, 2], [[1, 1], [1, 1]]],
        [[1, 2], [[1, 1], [1, 2]]],
        [[1, 2], [[1, 1], [1, 3]]],
        [[1, 2], [[1, 1], [2, 1]]],
        [[1, 2], [[1, 1], [2, 2]]],
        [[1, 2], [[1, 1], [2, 3]]],
        [[1, 2], [[1, 1], [3, 1]]],
        [[1, 2], [[1, 1], [3, 2]]],
        [[1, 2], [[1, 1], [3, 3]]],
        [[1, 2], [[1, 2], [1, 1]]],
        [[1, 2], [[1, 2], [1, 2]]],
        [[1, 2], [[1, 2], [1, 3]]],
        [[1, 2], [[1, 2], [2, 1]]],
        [[1, 2], [[1, 2], [2, 2]]],
        [[1, 2], [[1, 2], [2, 3]]],
        [[1, 2], [[1, 2], [3, 1]]],
        [[1, 2], [[1, 2], [3, 2]]],
        [[1, 2], [[1, 2], [3, 3]]],
        [[1, 2], [[1, 3], [1, 1]]],
        [[1, 2], [[1, 3], [1, 2]]],
        [[1, 2], [[1, 3], [1, 3]]],
        [[1, 2], [[1, 3], [2, 1]]],
        [[1, 2], [[1, 3], [2, 2]]],
        [[1, 2], [[1, 3], [2, 3]]],
        [[1, 2], [[1, 3], [3, 1]]],
        [[1, 2], [[1, 3], [3, 2]]],
        [[1, 2], [[1, 3], [3, 3]]],
        [[1, 2], [[2, 1], [1, 1]]],
        [[1, 2], [[2, 1], [1, 2]]],
        [[1, 2], [[2, 1], [1, 3]]],
        [[1, 2], [[2, 1], [2, 1]]],
        [[1, 2], [[2, 1], [2, 2]]],
        [[1, 2], [[2, 1], [2, 3]]],
        [[1, 2], [[2, 1], [3, 1]]],
        [[1, 2], [[2, 1], [3, 2]]],
        [[1, 2], [[2, 1], [3, 3]]],
        [[1, 2], [[2, 2], [1, 1]]],
        [[1, 2], [[2, 2], [1, 2]]],
        [[1, 2], [[2, 2], [1, 3]]],
        [[1, 2], [[2, 2], [2, 1]]],
        [[1, 2], [[2, 2], [2, 2]]],
        [[1, 2], [[2, 2], [2, 3]]],
        [[1, 2], [[2, 2], [3, 1]]],
        [[1, 2], [[2, 2], [3, 2]]],
        [[1, 2], [[2, 2], [3, 3]]],
        [[1, 2], [[2, 3], [1, 1]]],
        [[1, 2], [[2, 3], [1, 2]]],
        [[1, 2], [[2, 3], [1, 3]]],
        [[1, 2], [[2, 3], [2, 1]]],
        [[1, 2], [[2, 3], [2, 2]]],
        [[1, 2], [[2, 3], [2, 3]]],
        [[1, 2], [[2, 3], [3, 1]]],
        [[1, 2], [[2, 3], [3, 2]]],
        [[1, 2], [[2, 3], [3, 3]]],
        [[1, 2], [[3, 1], [1, 1]]],
        [[1, 2], [[3, 1], [1, 2]]],
        [[1, 2], [[3, 1], [1, 3]]],
        [[1, 2], [[3, 1], [2, 1]]],
        [[1, 2], [[3, 1], [2, 2]]],
        [[1, 2], [[3, 1], [2, 3]]],
        [[1, 2], [[3, 1], [3, 1]]],
        [[1, 2], [[3, 1], [3, 2]]],
        [[1, 2], [[3, 1], [3, 3]]],
        [[1, 2], [[3, 2], [1, 1]]],
        [[1, 2], [[3, 2], [1, 2]]],
        [[1, 2], [[3, 2], [1, 3]]],
        [[1, 2], [[3, 2], [2, 1]]],
        [[1, 2], [[3, 2], [2, 2]]],
        [[1, 2], [[3, 2], [2, 3]]],
        [[1, 2], [[3, 2], [3, 1]]],
        [[1, 2], [[3, 2], [3, 2]]],
        [[1, 2], [[3, 2], [3, 3]]],
        [[1, 2], [[3, 3], [1, 1]]],
        [[1, 2], [[3, 3], [1, 2]]],
        [[1, 2], [[3, 3], [1, 3]]],
        [[1, 2], [[3, 3], [2, 1]]],
        [[1, 2], [[3, 3], [2, 2]]],
        [[1, 2], [[3, 3], [2, 3]]],
        [[1, 2], [[3, 3], [3, 1]]],
        [[1, 2], [[3, 3], [3, 2]]],
        [[1, 2], [[3, 3], [3, 3]]],
    ]
    @test getIndices(agg) == aggInds

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    mode3 = Mode(0.4, 3.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1, mode2, mode3], 2, Energy)
    mol2 = Molecule([mode1, mode2], 2, Energy)
    mol3 = Molecule([mode1], 2, Energy)
    agg = Aggregate([mol1, mol2, mol3])

    @test getShifts(agg) == [[1.0, 2.0, 3.0], [1.0, 2.0], [1.0]]
    @test getFrequencies(agg) == [[0.2, 0.3, 0.4], [0.2, 0.3], [0.2]]

    mode1 = Mode(0.2, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, Energy)
    mol2 = Molecule([mode1], 2, Energy)
    agg = Aggregate([mol1, mol2])
    aggInds = getIndices(agg)
    FC1 = getFranckCondonFactors(agg)
    FC2 = getFranckCondonFactors(agg, aggInds)

    @test size(FC1) == (12, 12)
    @test FC1 == FC2
    FC = [
        1.0 0.0 0.0 0.0 0.7788007830714049 0.0 0.5506953149031838 0.0 0.7788007830714049 0.5506953149031838 0.0 0.0; 
        0.0 1.0 0.0 0.0 0.0 0.7788007830714049 0.0 0.5506953149031838 -0.5506953149031837 0.3894003915357024 -0.0 0.0; 
        0.0 0.0 1.0 0.0 -0.5506953149031837 -0.0 0.3894003915357024 0.0 0.0 0.0 0.7788007830714049 0.5506953149031838; 
        0.0 0.0 0.0 1.0 -0.0 -0.5506953149031837 0.0 0.3894003915357024 -0.0 0.0 -0.5506953149031837 0.3894003915357024; 
        0.7788007830714049 0.0 -0.5506953149031837 -0.0 1.0 0.0 0.0 0.0 0.6065306597126334 0.42888194248035344 -0.42888194248035333 -0.3032653298563167; 
        0.0 0.7788007830714049 -0.0 -0.5506953149031837 0.0 1.0 0.0 0.0 -0.42888194248035333 0.30326532985631666 0.30326532985631666 -0.21444097124017664; 
        0.5506953149031838 0.0 0.3894003915357024 0.0 0.0 0.0 1.0 0.0 0.42888194248035344 0.30326532985631677 0.30326532985631666 0.2144409712401767; 
        0.0 0.5506953149031838 0.0 0.3894003915357024 0.0 0.0 0.0 1.0 -0.3032653298563167 0.2144409712401767 -0.21444097124017664 0.15163266492815833; 
        0.7788007830714049 -0.5506953149031837 0.0 -0.0 0.6065306597126334 -0.42888194248035333 0.42888194248035344 -0.3032653298563167 1.0 0.0 0.0 0.0; 
        0.5506953149031838 0.3894003915357024 0.0 0.0 0.42888194248035344 0.30326532985631666 0.30326532985631677 0.2144409712401767 0.0 1.0 0.0 0.0; 
        0.0 -0.0 0.7788007830714049 -0.5506953149031837 -0.42888194248035333 0.30326532985631666 0.30326532985631666 -0.21444097124017664 0.0 0.0 1.0 0.0; 
        0.0 0.0 0.5506953149031838 0.3894003915357024 -0.3032653298563167 -0.21444097124017664 0.2144409712401767 0.15163266492815833 0.0 0.0 0.0 1.0
    ]
    @test 1e-12 > D(FC1, FC)

    @test getAggStateEnergy(agg, [1, 1], [[1], [2]]) == 0.4

    @test getAggStateEnergy(agg, [1, 1], [[1], [2]]) == 0.4

    @test OpenQuantumSystems.elIndOrder([1, 1, 1]) == 1
    @test OpenQuantumSystems.elIndOrder([1, 1, 2]) == 4
    @test OpenQuantumSystems.elIndOrder([2, 1, 1]) == 2

    agg.coupling[2, 3] = 200
    agg.coupling[3, 2] = 200

    FCSparse = getFranckCondonFactorsSparse(agg, aggInds)
    @test 1e-12 > D(FC, Matrix(FCSparse))

    FCSparse = getFranckCondonFactorsSparse(agg)
    @test 1e-12 > D(FC, Matrix(FCSparse))

    #=
    HamSparse1 = getAggHamiltonianSparse(agg, aggInds, FCSparse)
    HamSparse2 = getAggHamiltonianSparse(agg, aggInds)
    HamSparse3 = getAggHamiltonianSparse(agg)
    @test 1e-12 > D(Ham_ref, Matrix(HamSparse1.data))
    @test 1e-12 > D(Ham_ref, Matrix(HamSparse2.data))
    @test 1e-12 > D(Ham_ref, Matrix(HamSparse3.data))
    =#


    agg = Aggregate([mol1, mol2])
    agg.coupling[2, 3] = 200
    agg.coupling[3, 2] = 200
    agg2 = Aggregate([mol1, mol2], [0.0 0.0 0.0; 0.0 0.0 200.0; 0.0 200.0 0.0])
    @test agg2.coupling == agg.coupling

    mode1 = Mode(0.2, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, [2.0, 200.0])
    mol2 = Molecule([mode1], 2, [3.0, 300.0])
    agg = Aggregate([mol1, mol2])
    aggInds = getIndices(agg)
    vibIndices = getVibIndices(agg, aggInds)
    FCFact = getFranckCondonFactors(agg, aggInds)

    agg.coupling[2, 3] = 50
    agg.coupling[3, 2] = 50

    @test vibIndices == [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]

    Ham_S_ref = [
        5.0 0.0 0.0; 
        0.0 203.5 50.0; 
        0.0 50.0 302.5
    ]
    Ham_S = getAggHamSystemSmall(agg; groundEnergy=true)
    @test 1e-12 > D(Ham_S.data, Ham_S_ref)
    
    Ham_S_ref = [5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 5.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0; 0.0 0.0 5.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 5.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0; 0.0 0.0 -0.0 -0.0 203.5 0.0 0.0 0.0 30.326532985631673 21.444097124017674 -21.444097124017667 -15.163266492815836; 0.0 0.0 -0.0 -0.0 0.0 203.5 0.0 0.0 -21.444097124017667 15.163266492815833 15.163266492815833 -10.722048562008831; 0.0 0.0 0.0 0.0 0.0 0.0 203.5 0.0 21.444097124017674 15.163266492815838 15.163266492815833 10.722048562008835; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 203.5 -15.163266492815836 10.722048562008835 -10.722048562008831 7.581633246407916; 0.0 -0.0 0.0 -0.0 30.326532985631673 -21.444097124017667 21.444097124017674 -15.163266492815836 302.5 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 21.444097124017674 15.163266492815833 15.163266492815838 10.722048562008835 0.0 302.5 0.0 0.0; 0.0 -0.0 0.0 -0.0 -21.444097124017667 15.163266492815833 15.163266492815833 -10.722048562008831 0.0 0.0 302.5 0.0; 0.0 0.0 0.0 0.0 -15.163266492815836 -10.722048562008831 10.722048562008835 7.581633246407916 0.0 0.0 0.0 302.5]
    Ham_S = getAggHamSystemBig(agg, aggInds, FCFact; groundEnergy=true)
    @test 1e-12 > D(Ham_S.data, Ham_S_ref)
    
    Ham_B_ref = [
        0.2  0.0  0.0  0.0
        0.0  0.4  0.0  0.0
        0.0  0.0  0.4  0.0
        0.0  0.0  0.0  0.6 
    ]
    Ham_B = getAggHamBathSmall(agg; groundEnergy=true)
    @test 1e-12 > D(Ham_B.data, Ham_B_ref)
    
    Ham_B_ref = [
        0.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
        0.0  0.4  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.4  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.6  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0  0.4  0.0  0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0  0.0  0.4  0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.6  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.2  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.4  0.0  0.0
        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.4  0.0
        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.6
    ]
    Ham_B = getAggHamBathBig(agg; groundEnergy=true)
    @test 1e-12 > D(Ham_B.data, Ham_B_ref)
    
    Ham_0_ref = [5.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 5.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 5.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 5.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 203.7 0.0 0.0 0.0 30.326532985631673 21.444097124017674 -21.444097124017667 -15.163266492815836; 0.0 0.0 0.0 0.0 0.0 203.9 0.0 0.0 -21.444097124017667 15.163266492815833 15.163266492815833 -10.722048562008831; 0.0 0.0 0.0 0.0 0.0 0.0 203.9 0.0 21.444097124017674 15.163266492815838 15.163266492815833 10.722048562008835; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 204.1 -15.163266492815836 10.722048562008835 -10.722048562008831 7.581633246407916; 0.0 0.0 0.0 0.0 30.326532985631673 -21.444097124017667 21.444097124017674 -15.163266492815836 302.7 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 21.444097124017674 15.163266492815833 15.163266492815838 10.722048562008835 0.0 302.9 0.0 0.0; 0.0 0.0 0.0 0.0 -21.444097124017667 15.163266492815833 15.163266492815833 -10.722048562008831 0.0 0.0 302.9 0.0; 0.0 0.0 0.0 0.0 -15.163266492815836 -10.722048562008831 10.722048562008835 7.581633246407916 0.0 0.0 0.0 303.1]
    Ham_0 = getAggHamSystemBath(agg, aggInds, FCFact; groundEnergy=true)
    @test 1e-12 > D(Ham_0.data, Ham_0_ref)

    Ham_I_ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0]
    Ham_I = getAggHamInteraction(agg, aggInds, FCFact)
    @test 1e-12 > D(Ham_I.data, Ham_I_ref)

    Ham_ref = [5.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 5.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 5.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 5.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 203.7 0.0 -0.7071067811865476 0.0 30.326532985631673 21.444097124017674 -21.444097124017667 -15.163266492815836; 0.0 0.0 0.0 0.0 0.0 203.9 0.0 -0.7071067811865476 -21.444097124017667 15.163266492815833 15.163266492815833 -10.722048562008831; 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 203.9 0.0 21.444097124017674 15.163266492815838 15.163266492815833 10.722048562008835; 0.0 0.0 0.0 0.0 0.0 -0.7071067811865476 0.0 204.1 -15.163266492815836 10.722048562008835 -10.722048562008831 7.581633246407916; 0.0 0.0 0.0 0.0 30.326532985631673 -21.444097124017667 21.444097124017674 -15.163266492815836 302.7 -0.7071067811865476 0.0 0.0; 0.0 0.0 0.0 0.0 21.444097124017674 15.163266492815833 15.163266492815838 10.722048562008835 -0.7071067811865476 302.9 0.0 0.0; 0.0 0.0 0.0 0.0 -21.444097124017667 15.163266492815833 15.163266492815833 -10.722048562008831 0.0 0.0 302.9 -0.7071067811865476; 0.0 0.0 0.0 0.0 -15.163266492815836 -10.722048562008831 10.722048562008835 7.581633246407916 0.0 0.0 -0.7071067811865476 303.1]
    Ham = getAggHamiltonian(agg, aggInds, FCFact; groundEnergy=true)
    @test 1e-12 > D(Ham.data, Ham_ref)

    #=
    aggInds_ref = getIndices(agg)
    vibindices_ref = getVibIndices(agg, aggInds_ref)
    aggIndLen_ref = length(aggInds_ref)
    basis_ref = GenericBasis([aggIndLen_ref])
    FCFact_ref = getFranckCondonFactors(agg, aggInds_ref)
    FCProd_ref =
        getFCProd(agg, FCFact_ref, aggInds_ref, vibindices_ref)
    Ham_ref = getAggHamiltonian(
        agg,
        aggInds_ref,
        FCFact_ref;
        groundEnergy = true,
    )
    Ham_0_ref = getAggHamSysBath(agg, aggInds_ref; groundEnergy = true)
    Ham_I_ref = Ham_ref - Ham_0_ref

    aggInds, vibindices, aggIndLen, basis, FCFact, FCProd, Ham, Ham_0, Ham_I =
        setupAggregate(agg;  groundEnergy = true, verbose = false)
    @test aggInds_ref == aggInds
    @test vibindices_ref == vibindices
    @test aggIndLen_ref == aggIndLen
    @test basis_ref == basis
    @test FCFact_ref == FCFact
    @test FCProd_ref == FCProd
    @test Ham_ref == Ham
    @test Ham_0_ref == Ham_0
    @test Ham_I_ref == Ham_I

    @suppress aggInds, vibindices, aggIndLen, basis, FCFact, FCProd, Ham, Ham_0, Ham_I =
        setupAggregate(agg; groundEnergy = true, verbose = true)
    @test aggInds_ref == aggInds
    @test vibindices_ref == vibindices
    @test aggIndLen_ref == aggIndLen
    @test basis_ref == basis
    @test FCFact_ref == FCFact
    @test FCProd_ref == FCProd
    @test Ham_ref == Ham
    @test Ham_0_ref == Ham_0
    @test Ham_I_ref == Ham_I
    =#
end
