text:
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 "aggregateOperators" begin

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

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 3, [2.0, 200.0])
    mol2 = Molecule([mode2], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    aggTools = AggregateTools(aggCore)
    aggOperators = AggregateOperators(aggCore, aggTools; groundEnergy=true)
    # @test 1e-14 > D(aggOperators.Ham_sys.data, [5.0 0.0 0.0; 0.0 203.1 50.0; 0.0 50.0 302.6])

    ref = [0.25 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.55 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.85 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.45 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.75 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.05 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.65 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.95 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.25]
    @test 1e-14 > D(aggOperators.Ham_bath.data, ref)

    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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 203.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.6]
    # @test 1e-14 > D(aggOperators.Ham_S.data, ref)

    ref = [0.25 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.55 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.85 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.45 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.75 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 1.05 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.65 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.95 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 1.25 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.25 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.55 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.85 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.45 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.75 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 1.05 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.65 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.95 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 1.25 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.25 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.55 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.85 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.45 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.75 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 1.05 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.65 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.95 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 1.25]
    @test 1e-14 > D(aggOperators.Ham_B.data, ref)

    ref = aggOperators.Ham_S.data + aggOperators.Ham_B.data
    # @test 1e-14 > D(aggOperators.Ham_0.data, ref)

    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.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 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.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 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.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 -35.6747601569905 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 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.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 -20.258948470231473 -64.32523984300951 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 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.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -64.32523984300951 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 -42.83738007849525 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -57.162619921504756 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -57.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 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 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 -48.209345019623804 2.5323685587789493 2.5323685587789493; 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 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -51.790654980376196 2.1057475427305525e-15; 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 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -51.790654980376196; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -35.6747601569905 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 -64.32523984300951 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 1.6845980341844322e-14 -64.32523984300951 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 -2.0 0.0 0.0 -2.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 -10.12947423511573 14.325239843009507 -14.325239843009502 -42.83738007849525 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -57.162619921504756 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -57.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 -2.0 0.0 0.0 -2.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 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 -48.209345019623804 -2.532368558778949 2.532368558778948 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -51.790654980376196 -1.111293314979136e-15 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -51.790654980376196 0.0 -2.0 0.0 0.0 -2.0 0.0 0.0 -2.0 0.0]
    # @test 1e-9 > D(aggOperators.Ham_I.data, ref)

    ref = aggOperators.Ham_0.data + aggOperators.Ham_I.data
    # @test 1e-14 > D(aggOperators.Ham.data, ref)

    # groundEnergy = false
    aggOperators = AggregateOperators(aggCore, aggTools; groundEnergy=false)
    # @test 1e-14 > D(aggOperators.Ham_sys.data, [5.0 0.0 0.0; 0.0 203.1 50.0; 0.0 50.0 302.6])

    ref = [0.25 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.5500000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.8499999999999996 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.4499999999999993 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.75 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0499999999999998 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.6500000000000004 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9500000000000002 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.25]
    @test 1e-14 > D(aggOperators.Ham_bath.data, ref)

    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.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 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.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 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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.6]
    # @test 1e-14 > D(aggOperators.Ham_S.data, ref)

    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.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.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 0.0 0.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.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 0.0 0.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.0]
    @test 1e-14 > D(aggOperators.Ham_B.data, ref)

    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.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.29999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.79999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 199.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.90000000000003 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.20000000000005 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.8 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.1 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.40000000000003 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.3 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.6]
    # @test 1e-14 > D(aggOperators.Ham_0.data, ref)

    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.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 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.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 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.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 -35.6747601569905 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 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.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 -20.258948470231473 -64.32523984300951 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 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.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -64.32523984300951 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 -42.83738007849525 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -57.162619921504756 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 -0.7071067811865475 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -57.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 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 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 -48.209345019623804 2.5323685587789493 2.5323685587789493; 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 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -51.790654980376196 2.1057475427305525e-15; 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 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -51.790654980376196; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -35.6747601569905 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 -64.32523984300951 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 1.6845980341844322e-14 -64.32523984300951 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 -2.0 0.0 0.0 -2.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 -10.12947423511573 14.325239843009507 -14.325239843009502 -42.83738007849525 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -57.162619921504756 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -57.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 -2.0 0.0 0.0 -2.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 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 -48.209345019623804 -2.532368558778949 2.532368558778948 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -51.790654980376196 -1.111293314979136e-15 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -51.790654980376196 0.0 -2.0 0.0 0.0 -2.0 0.0 0.0 -2.0 0.0]
    # @test 1e-9 > D(aggOperators.Ham_I.data, ref)

    ref = aggOperators.Ham_0.data + aggOperators.Ham_I.data
    # @test 1e-14 > D(aggOperators.Ham.data, ref)

    # groundEnergy = false
    Ham_sys = getAggHamSystemSmall(aggCore, aggTools; groundEnergy = false)
    ref = [0.0 0.0 0.0; 0.0 198.1 50.0; 0.0 50.0 297.6]
    # @test 1e-14 > D(Ham_sys.data, ref)

    Ham_S = getAggHamSystemBig(aggCore, aggTools; groundEnergy = false)
    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.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 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.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 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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.6]
    # @test 1e-14 > D(Ham_S.data, ref)

    Ham_bath = getAggHamBathSmall(aggCore, aggTools; groundEnergy = false)
    ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 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.7999999999999998 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]
    @test 1e-14 > D(Ham_bath.data, ref)

    Ham_B = getAggHamBathBig(aggCore, aggTools; groundEnergy = false)
    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.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.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 0.0 0.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.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 0.0 0.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.0]
    @test 1e-14 > D(Ham_B.data, ref)

    Ham_0 = getAggHamSystemBath(aggCore, aggTools; groundEnergy = false)
    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.3000000000000007 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.5999999999999996 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.1999999999999993 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.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.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.7999999999999998 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.40000000000000036 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.7000000000000002 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 1.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 198.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.29999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 198.79999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.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 199.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.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.0 0.0 0.0 0.0 0.0 0.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.90000000000003 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.20000000000005 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 297.8 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.1 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.40000000000003 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.3 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 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 298.6]
    # @test 1e-14 > D(Ham_0.data, ref)

end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
using PowerSystems

############### OMIB Benchmark (Case 1 and 5) ########################

nodes_case1 = [Bus(1 , #number
                   "Bus 1", #Name
                   "REF" , #BusType (REF, PV, PQ)
                   0, #Angle in radians
                   1.05, #Voltage in pu
                   (min=0.94, max=1.06), #Voltage limits in pu
                   69), #Base voltage in kV
                   Bus(2 , "Bus 2"  , "PV" , 0 , 1.0 , (min=0.94, max=1.06), 69)]

branch_case1 = [Line("Line1", #name
                     true, #available
                     0.0, #active power flow initial condition (from-to)
                     0.0, #reactive power flow initial condition (from-to)
                     Arc(from=nodes_case1[1], to=nodes_case1[2]), #Connection between buses
                     0.01, #resistance in pu
                     0.05, #reactance in pu
                     (from=0.0, to=0.0), #susceptance in pu
                     18.046, #rate in MW
                     1.04)]  #angle limits (-min and max)

#Trip of a single circuit of Line 1 -> Resistance and Reactance doubled.
branch_case1_fault = [Line("Line1", #name
                           true, #available
                           0.0, #active power flow initial condition (from-to)
                           0.0, #reactive power flow initial condition (from-to)
                           Arc(from=nodes_case1[1], to=nodes_case1[2]), #Connection between buses
                           0.02, #resistance in pu
                           0.1, #reactance in pu
                           (from=0.0, to=0.0), #susceptance in pu
                           18.046, #rate in MW
                           1.04)]  #angle limits (-min and max)

loads_case1 = [PowerLoad("Bus1", true, nodes_case1[2], PowerSystems.ConstantPower, 0.3, 0.01, 0.3, 0.01)]

################ Benchmark 3 Bus case (Cases 2, 3 and 4)########################

nodes_case234   = [Bus(1 , "Bus 1"  , "REF" , 0 , 1.02  , (min=0.94, max=1.06), 138),
                    Bus(2 , "Bus 2"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138),
                    Bus(3 , "Bus 3"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138)]

branch_case234  =  [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_case234[1], to=nodes_case234[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case234[1], to=nodes_case234[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line3", true, 0.0, 0.0, Arc(from=nodes_case234[2], to=nodes_case234[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

#Trip of Line 1.
branch_case234_fault = [Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case234[1], to=nodes_case234[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                        Line("Line3", true, 0.0, 0.0, Arc(from=nodes_case234[2], to=nodes_case234[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

loads_case234 =   [PowerLoad("Bus1", true, nodes_case234[1], PowerSystems.ConstantPower, 1.5, 0.8, 1.5, 0.8),
                   PowerLoad("Bus2", true, nodes_case234[2], PowerSystems.ConstantPower, 1.5, 0.7, 1.5, 0.8),
                   PowerLoad("Bus3", true, nodes_case234[3], PowerSystems.ConstantPower, 0.5, 0.3, 0.5, 0.3)]

################ Benchmark 3 Bus case TG (Case 5)########################

nodes_case5   =    [Bus(1 , "Bus 1"  , "REF" , 0 , 1.05  , (min=0.94, max=1.06), 138),
                    Bus(2 , "Bus 2"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138),
                    Bus(3 , "Bus 3"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138)]

branch_case5 =     [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_case5[1], to=nodes_case5[2]), 0.01, 0.05, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case5[2], to=nodes_case5[3]), 0.02, 0.1, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line3", true, 0.0, 0.0, Arc(from=nodes_case5[1], to=nodes_case5[3]), 0.02, 0.1, (from=0.0, to=0.0), 100, 1.04)]

#Trip of a single circuit of Line 1 -> Resistance and Reactance doubled.
branch_case5_fault =     [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_case5[1], to=nodes_case5[2]), 0.02, 0.1, (from=0.0, to=0.0), 100, 1.04),
                          Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case5[2], to=nodes_case5[3]), 0.02, 0.1, (from=0.0, to=0.0), 100, 1.04),
                          Line("Line3", true, 0.0, 0.0, Arc(from=nodes_case5[1], to=nodes_case5[3]), 0.02, 0.1, (from=0.0, to=0.0), 100, 1.04)]

loads_case5 = [PowerLoad("Bus2", true, nodes_case5[2], PowerSystems.ConstantPower, 0.3, 0.05, 0.3, 0.05),
               PowerLoad("Bus3", true, nodes_case5[3], PowerSystems.ConstantPower, 0.3, 0.05, 0.3, 0.05)]

ext: jl


text:
## Variable construction

Variable(name::Union{Symbol,AbstractString}, indices::Int...) =
    Variable("$(name)$(join(map_subscripts.(indices), "₋"))")

const SUBSCRIPTS = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉']
const SUBSCRIPT_MAP = Dict([first(string(i)) => SUBSCRIPTS[i+1] for i = 0:9])
const SUBSCRIPT_TO_INT_MAP = Dict([SUBSCRIPTS[i+1] => i for i = 0:9])
map_subscripts(index) = join(SUBSCRIPT_MAP[c] for c in string(index))

function sort_key(v::Variable)
    name = string(ModelKit.name(v))
    sub_start = findnext(c -> c in ModelKit.SUBSCRIPTS, name, 1)
    if isnothing(sub_start)
        return name, Int[]
    end
    var_base = name[1:prevind(name, sub_start)]
    str_indices = split(name[sub_start:end], "₋")
    indices = map(str_indices) do str_index
        digits = map(s -> ModelKit.SUBSCRIPT_TO_INT_MAP[s], collect(str_index))
        n = length(digits)
        sum([digits[n-k+1] * 10^(k - 1) for k = 1:n])
    end
    var_base, reverse(indices)
end
Base.isless(a::Variable, b::Variable) = isless(sort_key(a), sort_key(b))
Base.sort!(vs::AbstractVector{Variable}; kwargs...) =
    permute!(vs, sortperm(sort_key.(vs); kwargs...))

Symbol(v::Variable) = name(v)

"""
    @var variable1 variable2 ...

Declare variables with the given names and automatically create the variable bindings.
The macro supports indexing notation to create `Array`s of variables.

## Examples

```julia-repl
julia> @var a b x[1:2] y[1:2,1:3]
(a, b, Variable[x₁, x₂], Variable[y₁₋₁ y₁₋₂ y₁₋₃; y₂₋₁ y₂₋₂ y₂₋₃])

julia> a
a

julia> b
b

julia> x
2-element Array{Variable,1}:
 x₁
 x₂

julia> y
2×3 Array{Variable,2}:
 y₁₋₁  y₁₋₂  y₁₋₃
 y₂₋₁  y₂₋₂  y₂₋₃
```
"""
macro var(args...)
    vars, exprs = buildvars(args; unique = false)
    :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
    $(Expr(:tuple, esc.(vars)...)))
end

"""
    @unique_var variable1 variable2

This is similar to [`@var`](@ref) with the only difference that the macro automatically
changes the names of the variables to ensure uniqueness. However, the binding is still
to the declared name.
This is useful to ensure that there are no name collisions.

## Examples

```julia-repl
julia> @unique_var a b
(a#591, b#592)

julia> a
a#591

julia> b
b#592
```
"""
macro unique_var(args...)
    vars, exprs = buildvars(args; unique = true)
    :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
    $(Expr(:tuple, esc.(vars)...)))
end

"""
    variables(prefix::Union{Symbol,String}, indices...)

Create an `Array` of variables with the given `prefix` and indices.
The expression  `@var x[1:3, 1:2]` is equivalent to  `x = variables(:x, 1:3, 1:2)`.
"""
function variables(prefix::Union{Symbol,String}, indices...)
    map(i -> Variable(prefix, i...), Iterators.product(indices...))
end

function buildvar(var; unique::Bool = false)
    if isa(var, Symbol)
        if unique
            # drop the first two ## from the gensym var
            varname = Symbol(String(gensym(var))[3:end])
        else
            varname = var
        end
        var, :($(esc(var)) = Variable($(QuoteNode(varname))))
    else
        isa(var, Expr) || error("Expected $var to be a variable name")
        Base.Meta.isexpr(var, :ref) ||
            error("Expected $var to be of the form varname[idxset]")
        (2 ≤ length(var.args)) || error("Expected $var to have at least one index set")
        varname = var.args[1]
        if unique
            prefix = String(gensym(varname))[3:end]
        else
            prefix = string(varname)
        end
        varname, :($(esc(varname)) = variables($prefix, $(esc.(var.args[2:end])...)))
    end
end

function buildvars(args; unique::Bool = false)
    vars = Symbol[]
    exprs = []
    for arg in args
        if arg isa Expr && arg.head == :tuple
            for _arg in arg.args
                var, expr = buildvar(_arg; unique = unique)
                push!(vars, var)
                push!(exprs, expr)
            end
        else
            var, expr = buildvar(arg; unique = unique)
            push!(vars, var)
            push!(exprs, expr)
        end
    end
    vars, exprs
end

"""
    unique_variable(var::Symbol, variables, parameters)

Create a variable with a unique name which doesn't clash with `variables` or
`parameters`.
If `var` is not possible the names `var##k` for `k=0,1,...` are tried until
one is possible,
"""
function unique_variable(var::Symbol, vars::Vector{Variable}, params::Vector{Variable})
    v = Variable(var)
    k = 0
    while (v in vars || v in params)
        v = Variable("$var##$k")
        k += 1
    end
    v
end

Base.adjoint(expr::Basic) = expr
Base.conj(expr::Basic) = expr
Base.transpose(expr::Basic) = expr
Base.broadcastable(v::Basic) = Ref(v)

"""
    variables(expr::Expression; parameters = Variable[])
    variables(exprs::AbstractVector{Expression}; parameters = Variable[])

Obtain all variables used in the given expression up to the ones declared in `parameters`.

## Example

```julia-repl
julia> @var x y a;
julia> variables(x^2 + y)
2-element Array{Variable,1}:
 x
 y

julia> variables([x^2 + a, y]; parameters = [a])
2-element Array{Variable,1}:
 x
 y
```
"""
variables(expr::Basic; kwargs...) = variables([expr]; kwargs...)
function variables(exprs::AbstractArray{<:Basic}; parameters = Variable[])
    S = Set{Variable}()
    for expr in exprs
        union!(S, _variables(expr))
    end
    setdiff!(sort!(collect(S)), parameters)
end

"""
    nvariables(expr::Expression; parameters = Variable[])
    nvariables(exprs::AbstractVector{Expression}; parameters = Variable[])

Obtain the number of variables used in the given expression not counting the the ones
declared in `parameters`.
"""
nvariables(exprs::Basic; kwargs...) = length(variables(exprs; kwargs...))
nvariables(exprs::AbstractArray{<:Basic}; kwargs...) = length(variables(exprs; kwargs...))

"""
    subs(expr::Expression, subsitutions::Pair...)
    subs(exprs::AbstractVector{<:Expression}, subsitutions::Pair...)

Apply the given substitutions to the given expressions.

## Examples
```julia-repl
@var x y

julia> subs(x^2, x => y)
y ^ 2

julia> subs(x * y, [x,y] => [x+2,y+2])
(x + 2) * (y + 2)

julia> subs([x + y, x^2], x => y + 2, y => x + 2)
2-element Array{Expression,1}:
 4 + x + y
 (2 + y)^2

# You can also use the callable syntax
julia> (x * y)([x,y] => [x+2,y+2])
 (x + 2) * (y + 2)
```
"""
subs(ex::Basic, args...) = subs(ex, ExpressionMap(args...))
subs(exs::AbstractArray{<:Basic}, args...) = subs.(exs, Ref(ExpressionMap(args...)))

"""
    evaluate(expr::Expression, subs...)
    evaluate(expr::AbstractArray{Expression}, subs...)

Evaluate the given expression.

## Example

```julia-repl
julia> @var x y;

julia> evaluate(x^2, x => 2)
4

julia> evaluate(x * y, [x,y] => [2, 3])
6

julia> evaluate([x^2, x * y], [x,y] => [2, 3])
2-element Array{Int64,1}:
 4
 6

# You can also use the callable syntax
julia> [x^2, x * y]([x,y] => [2, 3])
2-element Array{Int64,1}:
 4
 6
```
"""
function evaluate(expr::AbstractArray{<:Basic}, args...)
    out = map(to_number, subs(expr, args...))
    if isabstracttype(eltype(out))
        return to_smallest_eltype(out)
    else
        out
    end
end
evaluate(expr::Basic, args...) = to_number(subs(expr, args...))
(f::Union{Basic,AbstractArray{<:Basic}})(args...) = evaluate(f, args...)

"""
    evaluate!(u, expr::AbstractArray{Expression}, subs...)

Inplace form of [`evaluate`](@ref).
"""
function evaluate!(u::AbstractArray, expr::AbstractArray{<:Basic}, args...)
    map!(to_number, u, subs(expr, args...))
end

"""
    differentiate(expr::Expression, var::Variable, k = 1)
    differentiate(expr::AbstractVector{Expression}, var::Variable, k = 1)

Compute the `k`-th derivative of `expr` with respect to the given variable `var`.

    differentiate(expr::Expression, vars::AbstractVector{Variable})

Compute the partial derivatives of `expr` with respect to the given variable variables `vars`.
Retuns a `Vector` containing the partial derivatives.

    differentiate(exprs::AbstractVector{Expression}, vars::AbstractVector{Variable})

Compute the partial derivatives of `exprs` with respect to the given variable variables `vars`.
Returns a `Matrix` where the each row contains the partial derivatives for a given expression.
"""
function differentiate(expr::Basic, vars::AbstractVector{Variable})
    [differentiate(expr, v) for v in vars]
end
function differentiate(exprs::AbstractVector{<:Basic}, var::Variable, k::Int = 1)
    [differentiate(e, var, k) for e in exprs]
end
function differentiate(exprs::AbstractVector{<:Basic}, vars::AbstractVector{Variable})
    [differentiate(e, v) for e in exprs, v in vars]
end

"""
    monomials(variables::AbstractVector, d::Integer; affine = true)

Create all monomials of a given degree in the given `variables`.

```julia-repl
julia> @var x y
(x, y)

julia> monomials([x,y], 2; affine = false)
3-element Array{Operation,1}:
 x ^ 2
 x * y
 y ^ 2
```
"""
function monomials(
    vars::AbstractVector{<:Union{Variable,Expression}},
    d::Integer;
    affine::Bool = true,
    homogeneous::Bool = !affine,
)
    n = length(vars)
    exps = monomials_exponents(n, d; affine = !homogeneous)
    map(exps) do exp
        prod(i -> vars[i]^exp[i], 1:n)
    end
end
function monomials_exponents(n, d; affine::Bool)
    if affine
        pred = x -> sum(x) ≤ d
    else
        pred = x -> sum(x) == d
    end
    E = map(Iterators.filter(pred, Iterators.product(Iterators.repeated(0:d, n)...))) do e
        collect(e)
    end

    sort!(E, lt = td_order)
    E
end

function td_order(x, y)
    sx = sum(x)
    sy = sum(y)
    sx == sy ? x > y : sx > sy
end
function monomials(
    vars::AbstractVector{<:Union{Variable,Expression}},
    D::AbstractVector{<:Integer},
)
    D = sort(D; rev = true)
    M = monomials(vars, D[1]; affine = false)
    for i = 2:length(D)
        append!(M, monomials(vars, D[i]; affine = false))
    end
    M
end

"""
    dense_poly(vars::AbstractVector{Variable}, d::Integer;
               homogeneous = false,
               coeff_name::Symbol = gensym(:c))

Create a dense polynomial of degree `d` in the given variables `variables` where
each coefficient is a parameter. Returns a tuple with the first argument being the polynomial
and the second the parameters.

```julia-repl
julia> @var x y;

julia> f, c = dense_poly([x, y], 2, coeff_name = :q);

julia> f
 q₆ + x*q₄ + x^2*q₁ + y*q₅ + y^2*q₃ + x*y*q₂

julia> c
6-element Array{Variable,1}:
 q₁
 q₂
 q₃
 q₄
 q₅
 q₆
```
"""
function dense_poly(
    vars::AbstractVector{<:Union{Variable,Expression}},
    d::Integer;
    homogeneous::Bool = false,
    coeff_name::Symbol = gensym(:c),
)
    M = monomials(vars, d; affine = !homogeneous)
    c = Variable.(coeff_name, 1:length(M))
    sum(c .* M), c
end

"""
    coeffs_as_dense_poly(f, vars, d; homogeneous = false)

Given a polynomial `f` this returns a vector `c` of coefficients such that
`subs(dense_poly(vars, d; homogeneous = homogeneous), c) == f`.

## Example

```julia
@var x[1:3]
f, c = dense_poly(x, 3, coeff_name = :c)
g = x[1]^3+x[2]^3+x[3]^3-1
gc = coeffs_as_dense_poly(g, x, 3)
subs(f, c => gc) == g
```
```
true
```
"""
function coeffs_as_dense_poly(
    f::Expression,
    vars::AbstractVector{<:Union{Variable,Expression}},
    d::Integer;
    homogeneous::Bool = false,
)
    exps = monomials_exponents(length(vars), d; affine = !homogeneous)
    D = to_dict(expand(f), vars)
    to_smallest_eltype(map(exps) do exp
        if haskey(D, exp)
            to_number(D[exp])
        else
            false
        end
    end)
end

"""
    rand_poly(T = ComplexF64, vars::AbstractVector{Variable}, d::Integer; homogeneous::Bool = false)

Create a random dense polynomial of degree `d` in the given variables `variables`.
Each coefficient is sampled independently via `randn(T)`.

```julia-repl
julia> @var x y;

julia> rand_poly(Float64, [x, y], 2)
0.788764085756728 - 0.534507647623108*x - 0.778441366874946*y -
 0.128891763280247*x*y + 0.878962738754971*x^2 + 0.550480741774464*y^2
```
"""
function rand_poly(vars::AbstractVector, d::Integer; kwargs...)
    rand_poly(ComplexF64, vars, d; kwargs...)
end
function rand_poly(T, vars::AbstractVector, d::Integer; homogeneous::Bool = false)
    M = monomials(vars, d; affine = !homogeneous)
    sum(randn(T, length(M)) .* M)
end

"""
    expand(e::Expression)

Expand a given expression.

```julia-repl
julia> @var x y
(x, y)

julia> expand((x + y) ^ 2)
2*x*y + x^2 + y^2
```
"""
expand(e::Basic) = symengine_expand(e)

struct PolynomialError <: Exception end
function Base.showerror(io::IO, ::PolynomialError)
    print(io, "Encountered an unexpected rational expression.")
end

"""
    to_dict(expr::Expression, vars::AbstractVector{Variable})

Return the coefficients of `expr` w.r.t. the given variables `vars`. Assumes that `expr`
is expanded and representing a polynomial.
Throws a `PolynomialError` if a rational expression is encountered.
"""
function to_dict(expr::Expression, vars::AbstractVector{Variable})
    mul_args, pow_args = ExprVec(), ExprVec()
    dict = Dict{Vector{Int},Expression}()

    if class(expr) == :Add
        for op in args(expr)
            to_dict_op!(dict, op, vars, mul_args, pow_args)
        end
    else
        to_dict_op!(dict, expr, vars, mul_args, pow_args)
    end

    dict
end

function to_dict_op!(dict, op, vars, mul_args, pow_args)
    cls = class(op)
    d = zeros(Int, length(vars))
    coeff = Expression(1)
    if cls == :Mul
        op_args = args!(mul_args, op)
        for arg in op_args
            arg_cls = class(arg)
            is_coeff = true
            if arg_cls == :Symbol
                for (i, v) in enumerate(vars)
                    if v == arg
                        d[i] = 1
                        is_coeff = false
                        break
                    end
                end
            elseif arg_cls == :Pow
                vec = args!(pow_args, arg)
                x = vec[1]
                k = convert(Int, vec[2])
                if k < 0
                    throw(PolynomialError())
                end
                for (i, v) in enumerate(vars)
                    if x == v
                        d[i] = k
                        is_coeff = false
                        break
                    end
                end
            end
            if is_coeff
                mul!(coeff, coeff, arg)
            end
        end
    elseif cls == :Symbol
        is_in_vars = false
        for (i, v) in enumerate(vars)
            if v == op
                d[i] = 1
                is_in_vars = true
                break
            end
        end
        if !is_in_vars
            coeff = copy(op)
        end
    elseif cls == :Pow
        vec = args!(pow_args, op)
        # check that base is one of the variables
        x = vec[1]
        is_var_pow = false
        k = convert(Int, vec[2])
        if k < 0
            throw(PolynomialError())
        end
        for (i, v) in enumerate(vars)
            if x == v
                d[i] = k
                is_var_pow = true
                break
            end
        end
        if !is_var_pow
            coeff = copy(op)
        end
    else
        coeff = copy(op)
    end

    if haskey(dict, d)
        add!(dict[d], dict[d], coeff)
    else
        dict[d] = coeff
    end
    dict
end

"""
    exponents_coefficients(f::Expression, vars::AbstractVector{Variable}; expanded = false)

Return a matrix `M` containing the exponents for all occuring terms
(one term per column) and a vector `c` containing the corresponding coefficients.
Expands the given expression `f` unless `expanded = true`.
Throws a `PolynomialError` if a rational expression is encountered.
"""
function exponents_coefficients(
    f::Expression,
    vars::AbstractVector{Variable};
    expanded::Bool = false,
    unpack_coeffs::Bool = true,
)
    expanded || (f = expand(f))
    D = to_dict(f, vars)
    # make exponents to matrix
    m = length(D)
    E = collect(keys(D))
    perm = sortperm(E; lt = td_order)
    permute!(E, perm)
    n = length(vars)
    M = zeros(Int32, n, m)
    for (j, Eⱼ) in enumerate(E)
        for (i, dᵢ) in enumerate(Eⱼ)
            M[i, j] = dᵢ
        end
    end
    if unpack_coeffs
        coeffs = to_smallest_eltype(to_number.(values(D)))
    else
        coeffs = collect(values(D))
    end
    permute!(coeffs, perm)
    M, coeffs
end

"""
    coefficients(f::Expression, vars::AbstractVector{Variable}; expanded = false)

Return all coefficients of the given polynomial `f` for the given variables `vars`.
If `expanded = true` then this assumes that the expression `f` is already expanded,
e.g., with [`expand`](@ref).
"""
function coefficients(f::Expression, vars::AbstractVector{Variable}; expanded::Bool = false)
    expanded || (f = expand(f))
    D = to_dict(f, vars)
    m = length(D)
    E = collect(keys(D))
    perm = sortperm(E; lt = td_order)
    make_to_number = true
    for c in values(D)
        if !is_number(c)
            make_to_number = false
            break
        end
    end

    if make_to_number
        permute!(to_number.(values(D)), perm)
    else
        permute!(collect(values(D)), perm)
    end
end
coefficients(f::Expression, var::Variable; kwargs...) = coefficients(f, [var]; kwargs...)

"""
    degree(f::Expression, vars = variables(f); expanded = false)

Compute the degree of the expression `f`  in `vars`.
Unless `expanded` is `true` the expression is first expanded.
"""
function degree(
    f::Expression,
    vars::AbstractVector{Variable} = variables(f);
    expanded::Bool = false,
)
    if !expanded
        f = ModelKit.expand(f)
    end
    dicts = ModelKit.to_dict(f, vars)
    maximum(sum, keys(dicts))
end

"""
    degrees(f::AbstractVector{Expression}, vars = variables(f); expanded = false)

Compute the degrees of the expressions `f` in `vars`.
Unless `expanded` is `true` the expressions are first expanded.
"""
function degrees(
    f::AbstractVector{Expression},
    vars::AbstractVector{Variable} = variables(f);
    expanded::Bool = false,
)
    if !expanded
        f = expand.(f)
    end
    dicts = to_dict.(f, Ref(vars))
    maximum.(sum, keys.(dicts))
end

function LinearAlgebra.det(M::AbstractMatrix{<:Union{Variable,Expression}})
    m, n = size(M)
    m == n || error("Cannot take determinant of non-square matrix.")
    if m > 2
        return sum(
            (-1)^(i - 1) * M[i, 1] * LinearAlgebra.det(M[1:end.!=i, 2:end]) for i = 1:m
        )
    elseif m == 2
        return M[1, 1] * M[2, 2] - M[2, 1] * M[1, 2]
    else # m == 1
        return M[1, 1]
    end
end
LinearAlgebra.det(M::LinearAlgebra.Symmetric{<:Union{Variable,Expression}}) =
    LinearAlgebra.det(Matrix(M))

function is_homogeneous(f::Expression, vars::Vector{Variable}; expanded::Bool = false)
    if !expanded
        f = expand(f)
    end
    try
        exponents = keys(to_dict(f, vars))

        d = sum(first(exponents))
        all(e -> sum(e) == d, exponents)
    catch err
        if err isa PolynomialError
            return false
        else
            rethrow(e)
        end
    end
end
################
# Optimization #
################

"""
    horner(f::Expression, vars = variables(f))

Rewrite `f` using a multi-variate horner schema.

### Example
```julia-repl
julia> @var u v c[1:3]
(u, v, Variable[c₁, c₂, c₃])

julia> f = c[1] + c[2] * v + c[3] * u^2 * v^2 + c[3]u^3 * v
c₁ + v*c₂ + u^2*v^2*c₃ + u^3*v*c₃

julia> horner(f)
c₁ + v*(c₂ + u^3*c₃ + u^2*v*c₃)
```
"""
function horner(f::Expression, vars = variables(f))
    try
        M, coeffs = exponents_coefficients(f, vars; expanded = true, unpack_coeffs = false)
        multivariate_horner(M, coeffs, vars)
    catch err
        if err isa PolynomialError
            return f
        else
            rethrow(err)
        end
    end
end

function horner(coeffs::AbstractVector, var::Variable)
    d = length(coeffs) - 1
    h = copy(coeffs[d+1])
    @inbounds for k = d:-1:1
        mul!(h, h, var)
        add!(h, h, coeffs[k])
    end
    h
end

function multivariate_horner(
    M::AbstractMatrix,
    coeffs::AbstractVector,
    vars::AbstractVector{Variable},
)
    n, m = size(M)
    if m == 1
        fast_c = coeffs[1]
        for (i, var) in enumerate(vars)
            if M[i, 1] != 0
                if M[i, 1] == 1
                    fast_c *= var
                else
                    fast_c *= var^M[i, 1]
                end
            end
        end
        return fast_c::Expression
    elseif n == 1
        d = maximum(M)
        uni_coeff_indices = Union{Nothing,Int}[nothing for _ = 0:d]
        for j = 1:size(M, 2)
            uni_coeff_indices[M[1, j]+1] = j
        end

        var_coeffs = map(uni_coeff_indices) do ind
            if isnothing(ind)
                return zero(vars[1])
            else
                coeffs[ind]
            end
        end

        horner(var_coeffs, vars[1])
    else
        counts = zeros(Int, n)
        @inbounds for j = 1:size(M, 2), i = 1:size(M, 1)
            counts[i] += M[i, j] > 0
        end

        _, var_ind = findmax(counts)
        var = vars[var_ind]

        # compute degree in vars[var_ind]
        d = 0
        for j = 1:size(M, 2)
            d = max(d, M[var_ind, j])
        end

        # find coefficients
        coeff_indices = [Int[] for _ = 0:d]
        for j = 1:size(M, 2)
            push!(coeff_indices[M[var_ind, j]+1], j)
        end

        reduced_vars_ind = Int[]
        for i = 1:n
            if counts[i] > 0 && i != var_ind
                push!(reduced_vars_ind, i)
            end
        end
        reduced_vars = view(vars, reduced_vars_ind)

        var_coeffs = map(coeff_indices) do ind
            isempty(ind) && return zero(var)
            if length(ind) == 1
                j = first(ind)
                c = coeffs[j]
                for i in reduced_vars_ind
                    if M[i, j] != 0
                        if M[i, j] == 1
                            c *= vars[i]
                        else
                            c *= vars[i]^M[i, j]
                        end
                    end
                end
                return c::Expression
            end
            M_ind = view(M, reduced_vars_ind, ind)
            coeffs_ind = view(coeffs, ind)
            multivariate_horner(M_ind, coeffs_ind, reduced_vars)::Expression
        end
        horner(var_coeffs, var)
    end
end

#########################
## System and Homotopy ##
#########################

############
## System ##
############

function check_vars_params(f, vars, params)
    vars_params = params === nothing ? vars : [vars; params]
    Δ = setdiff(variables(f), vars_params)
    isempty(Δ) || throw(
        ArgumentError(
            "Not all variables or parameters of the system are given. Missing: " *
            join(Δ, ", "),
        ),
    )
    if params !== nothing
        both = Set(vars) ∩ Set(params)
        if !isempty(both)
            error("$(join(both, ", ")) appear as parameters and variables")
        end
    end
    nothing
end

"""
    System(exprs::AbstractVector{Expression};
                variables = variables(exprssion),
                parameters = Variable[])

Create a system from the given [`Expression`](@ref)s `exprs`.
The `variables` determine also the variable ordering.
The `parameters` argument allows to declare certain [`Variable`](@ref)s as parameters.

    System(support::AbstractVector{<:AbstractMatrix{<:Integer}},
           coefficients::AbstractVector{<:AbstractVector};
           variables,
           parameters = Variable[])

Create a system from the given support and coefficients.

## Examples
```julia-repl
julia> @var x y;
julia> F = System([x^2, y^2]; variables = [y, x])
System of length 2
 2 variables: y, x

 x^2
 y^2

# Systems are callable.
# This evaluates F at y=2 and x=3
julia> F([2, 3])
2-element Array{Int64,1}:
 9
 4
```

It is also possible to declare parameters.
```julia-repl
julia> @var x y a b;
julia> F = System([x^2 + a, y^2 + b]; variables = [y, x], parameters = [a, b])
System of length 2
 2 variables: y, x
 2 parameters: a, b

 a + x^2
 b + y^2

julia> F([2, 3], [5, -2])
 2-element Array{Int64,1}:
  14
   2
```
"""
struct System
    expressions::Vector{Expression}
    variables::Vector{Variable}
    parameters::Vector{Variable}
    variable_groups::Union{Nothing,Vector{Vector{Variable}}}
    _jacobian::Ref{Union{Nothing,Matrix{Expression}}}

    function System(
        exprs::Vector{Expression},
        vars::Vector{Variable},
        params::Vector{Variable},
        variable_groups::Union{Nothing,Vector{Vector{Variable}}} = nothing,
    )
        check_vars_params(exprs, vars, params)
        if !isnothing(variable_groups)
            vars == reduce(vcat, variable_groups) ||
                throw(ArgumentError("Variable groups and variables don't match."))
        end
        new(
            exprs,
            vars,
            params,
            variable_groups,
            Ref{Union{Nothing,Matrix{Expression}}}(nothing),
        )
    end
end

System(exprs; kwargs...) = System(convert(Vector{Expression}, exprs); kwargs...)
function System(
    exprs::AbstractVector{Expression};
    parameters::Union{Nothing,Vector{Variable}} = nothing,
    variable_groups::Union{Nothing,Vector{Vector{Variable}}} = nothing,
    variables::Union{Nothing,Vector{Variable}} = nothing,
)
    vars, params = _default_vars_params(exprs, variables, parameters, variable_groups)
    System(convert(Vector{Expression}, exprs), vars, params, variable_groups)
end

function _default_vars_params(exprs, vars, parameters, variable_groups)
    all_vars = variables(exprs)
    if isnothing(variable_groups)
        parameters = something(parameters, Variable[])
        if isnothing(vars)
            setdiff(all_vars, parameters), parameters
        else
            vars, parameters
        end
    elseif isnothing(vars)
        vars = reduce(vcat, variable_groups)
        if isnothing(parameters)
            vars, setdiff(all_vars, vars)
        else
            vars, parameters
        end
    else
        if isnothing(parameters)
            vars, Variable[]
        else
            vars, parameters
        end
    end
end

function System(
    exprs::AbstractVector,
    variables::Vector{Variable};
    parameters::Vector{Variable} = Variable[],
)
    System(convert(Vector{Expression}, exprs), variables, parameters)
end
function System(
    exprs::AbstractVector,
    variables::Vector{Variable},
    parameters::Vector{Variable},
)
    System(convert(Vector{Expression}, exprs), variables, parameters)
end

function System(
    F::AbstractVector{<:MP.AbstractPolynomial};
    parameters = similar(MP.variables(F), 0),
    variables = setdiff(MP.variables(F), parameters),
    variable_groups = nothing,
)
    vars = map(variables) do v
        name, ind = MP.name_base_indices(v)
        Variable(name, ind...)
    end
    params = Variable[]
    for v in parameters
        name, ind = MP.name_base_indices(v)
        push!(params, Variable(name, ind...))
    end
    if variable_groups === nothing
        var_groups = nothing
    else
        var_groups = map(var_groups) do group
            map(group) do v
                name, ind = MP.name_base_indices(v)
                Variable(name, ind...)
            end
        end
    end
    variables_parameters = [variables; parameters]
    vars_params = [vars; params]
    G = map(F) do f
        sum(MP.terms(f)) do t
            c = MP.coefficient(t)
            c * prod(w^MP.degree(t, v) for (v, w) in zip(variables_parameters, vars_params))
        end
    end
    System(G, variables = vars, parameters = params, variable_groups = var_groups)
end

System(F::System) = F

function System(
    support::AbstractVector{<:AbstractMatrix{<:Integer}},
    coefficients::AbstractVector{<:AbstractVector};
    variables::AbstractVector{Variable},
    parameters::AbstractVector{Variable} = Variable[],
)

    System(
        map(support, coefficients) do A, c
            fi = Expression(0)
            for (k, a) in enumerate(eachcol(A))
                ModelKit.add!(fi, fi, c[k] * prod(variables .^ convert.(Int, a)))
            end
            fi
        end,
        variables,
        parameters,
    )
end

Base.hash(S::System, u::UInt64) =
    hash(S.expressions, hash(S.variables, hash(S.parameters, u)))

function Base.show(io::IO, F::System)
    if !get(io, :compact, false)
        println(io, "System of length $(length(F.expressions))")
        if isnothing(F.variable_groups)
            print(io, " $(length(F.variables)) variables: ", join(F.variables, ", "))
        else
            print(
                io,
                length(F.variables),
                " variables (",
                length(F.variable_groups),
                " groups): ",
                join("[" .* join.(F.variable_groups, Ref(", ")) .* "]", ", "),
            )
        end
        if !isempty(F.parameters)
            print(io, "\n $(length(F.parameters)) parameters: ", join(F.parameters, ", "))
        end
        print(io, "\n\n")
        for i = 1:length(F)
            print(io, " ", F.expressions[i])
            i < length(F) && print(io, "\n")
        end
    else
        print(io, "[")
        for i = 1:length(F)
            print(io, F.expressions[i])
            i < length(F) && print(io, ", ")
        end
        print(io, "]")
    end
end

"""
    evaluate(F::System, x, p = nothing)

Evaluates the system `F` at `(x, p)`.

# Example

```julia
@var x y
F = System([x^2 + 3y, (y*x+1)^3])
evaluate(F, [2, 3])
```
```
2-element Array{Int32,1}:
  13
 343
```
"""

evaluate(F::System, x::AbstractVector) = evaluate(F.expressions, F.variables => x)
function evaluate(F::System, x::AbstractVector, p::AbstractVector)
    if isempty(F.parameters)
        evaluate(F.expressions, F.variables => x)
    else
        evaluate(F.expressions, F.variables => x, F.parameters => p)
    end
end
(F::System)(x::AbstractVector, p::Nothing = nothing) = evaluate(F, x)
(F::System)(x::AbstractVector, p::AbstractVector) = evaluate(F, x, p)

"""
    jacobian(F::System)

Computes the symbolic Jacobian of the system `F`.

# Example

```julia
@var x y
F = System([x^2 + 3y, (y*x+1)^3])
jacobian(F)
```
```
2×2 Array{Expression,2}:
             2*x                3
 3*y*(1 + x*y)^2  3*x*(1 + x*y)^2
```
"""
function jacobian(F::System)
    if isnothing(F._jacobian[])
        F._jacobian[] = differentiate(F.expressions, F.variables)
    end
    F._jacobian[]::Matrix{Expression}
end

"""
    jacobian(F::System, x, p = nothing)

Evaluates the Jacobian of the system `F` at `(x, p)`.

# Example

```julia
@var x y
F = System([x^2 + 3y, (y*x+1)^3])
jacobian(F, [2, 3])
```
```
2×2 Array{Int32,2}:
   4    3
 441  294
```
"""
function jacobian(F::System, x, p = nothing)
    if p isa Nothing
        evaluate(jacobian(F), F.variables => x)
    else
        evaluate(jacobian(F), F.variables => x, F.parameters => p)
    end
end

function Base.:(==)(F::System, G::System)
    F.expressions == G.expressions &&
        F.variables == G.variables &&
        F.parameters == G.parameters
end

Base.size(F::System) = (length(F.expressions), length(F.variables))
Base.size(F::System, i::Integer) = size(F)[i]
Base.length(F::System) = length(F.expressions)
Base.copy(F::System) = Base.deepcopy(F)

"""
    nvariables(F::System)

Returns the number of variables of the given system `F`.
"""
nvariables(F::System) = length(F.variables)

"""
    nparameters(F::System)

Returns the number of parameters of the given system `F`.
"""
nparameters(F::System) = length(F.parameters)

"""
    expressions(F::System)

Returns the expressions of the given system `F`.
"""
expressions(F::System) = F.expressions

"""
    variables(F::System)

Returns the variables of the given system `F`.
"""
variables(F::System) = F.variables

"""
    variable_groups(F::System)

Returns the variable groups of the given system `F`.
"""
variable_groups(F::System) = F.variable_groups

"""
    parameters(F::System)

Returns the parameters of the given system `F`.
"""
parameters(F::System) = F.parameters

"""
    degrees(F::System)

Return the degrees of the given system.
"""
degrees(F::System) = degrees(F.expressions, F.variables)

"""
    multi_degrees(F::System)

Return the degrees with respect to the given variable groups.
"""
function multi_degrees(F::System)
    if isnothing(F.variable_groups)
        reshape(degrees(F.expressions, F.variables), 1, nvariables(F))
    else
        reduce(vcat, map(vars -> degrees(F.expressions, vars)', F.variable_groups))
    end
end

"""
    support_coefficients(F::System)

Return the support of the system and the corresponding coefficients.
"""
function support_coefficients(F::System)
    supp_coeffs = exponents_coefficients.(F.expressions, Ref(F.variables))
    first.(supp_coeffs), last.(supp_coeffs)
end

"""
    is_homogeneous(F::System)

Checks whether a given system is homogeneous w.r.t to its variables
(resp. variable groups).
"""
function is_homogeneous(F::System)
    if isnothing(F.variable_groups)
        all(f -> is_homogeneous(f, F.variables), F.expressions)
    else
        all(F.expressions) do f
            g = expand(f)
            all(vars -> is_homogeneous(g, vars; expanded = true), F.variable_groups)
        end
    end
end

Base.iterate(F::System) = iterate(F.expressions)
Base.iterate(F::System, state) = iterate(F.expressions, state)

Base.push!(F::System, f::Expression) = (push!(F.expressions, f); F)
Base.append!(F::System, f::AbstractVector{Expression}) = (append!(F.expressions, f); F)

function Base.intersect(F::System, G::System)
    exprs = [
        F.expressions
        G.expressions
    ]
    vars = [
        F.variables
        setdiff(G.variables, F.variables)
    ]
    params = [
        F.parameters
        setdiff(G.parameters, F.parameters)
    ]
    System(exprs, vars, params, F.variable_groups)
end
function Base.intersect(F::System, G::AbstractVector{<:Expression})
    exprs = [F.expressions; G]
    vars = [
        F.variables
        setdiff(variables(G), F.variables)
    ]
    System(exprs, vars, F.parameters, F.variable_groups)
end
Base.intersect(F::AbstractVector{<:Expression}, G::System) = intersect(G, F)

"""
    optimize(F::System)

Optimize the evaluation cost of the given system `F`. This applies a multivariate horner
schema to the given expressions. See also [`horner`](@ref).

### Example

```julia
julia> f
System of length 4
 4 variables: z₁, z₂, z₃, z₄

 z₁ + z₂ + z₃ + z₄
 z₂*z₁ + z₂*z₃ + z₃*z₄ + z₄*z₁
 z₂*z₃*z₁ + z₂*z₃*z₄ + z₂*z₄*z₁ + z₃*z₄*z₁
 -1 + z₂*z₃*z₄*z₁

julia> optimize(f)
System of length 4
 4 variables: z₁, z₂, z₃, z₄

 z₁ + z₂ + z₃ + z₄
 (z₂ + z₄)*z₁ + (z₂ + z₄)*z₃
 z₁*(z₃*z₄ + (z₃ + z₄)*z₂) + z₂*z₃*z₄
 -1 + z₂*z₃*z₄*z₁
```
"""
function optimize(F::System)
    System(
        horner.(F.expressions, Ref(F.variables)),
        F.variables,
        F.parameters,
        F.variable_groups,
    )
end

## Conversion from MultivariatePolynomials
function system_with_coefficents_as_params(
    F::AbstractVector{<:MP.AbstractPolynomial};
    variables = MP.variables(F),
    variable_groups = nothing,
)
    vars = map(variables) do v
        name, ind = MP.name_base_indices(v)
        Variable(name, ind...)
    end
    if variable_groups === nothing
        var_groups = nothing
    else
        var_groups = map(var_groups) do group
            map(group) do v
                name, ind = MP.name_base_indices(v)
                Variable(name, ind...)
            end
        end
    end
    param_name = gensym(:c)
    k = 1
    params = Variable[]
    target_params = MP.coefficienttype(F)[]
    G = map(F) do f
        sum(MP.terms(f)) do t
            c = Variable(param_name, k)
            k += 1
            push!(params, c)
            push!(target_params, MP.coefficient(t))
            c * prod(w^MP.degree(t, v) for (v, w) in zip(variables, vars))
        end
    end
    sys = System(G, variables = vars, parameters = params, variable_groups = var_groups)
    sys, target_params
end

##############
## Homotopy ##
##############
"""
    Homotopy(exprs, vars, t, parameters = Variable[])

Create a homotopy `H(vars,t)` from the given [`Expression`](@ref)s `exprs` where `vars` are the given
variables and `t` is the dedicated variable parameterizing the family of systems.
The `parameters` argument allows to declare certain [`Variable`](@ref)s as parameters.

## Example
```julia-repl
julia> @var x y t;

julia> H = Homotopy([x + t, y + 2t], [y, x], t)
Homotopy in t of length 2
 2 variables: y, x

 t + x
 2*t + y

julia> H([2, 3], 0)
2-element Array{Int64,1}:
 3
 2

julia> H([2, 3], 1)
2-element Array{Int64,1}:
 4
 4
```

It is also possible to declare additional variables.
```julia-repl
julia> @var x y t a b;
julia> H = Homotopy([x^2 + t*a, y^2 + t*b], [x, y], t, [a, b])
Homotopy in t of length 2
 2 variables: x, y
 2 parameters: a, b

 a*t + x^2
 b*t + y^2
julia> H([2, 3], 1, [5, 2])
2-element Array{Int64,1}:
 9
 11
```
"""
struct Homotopy
    expressions::Vector{Expression}
    variables::Vector{Variable}
    t::Variable
    parameters::Vector{Variable}

    function Homotopy(
        exprs::Vector{Expression},
        vars::Vector{Variable},
        t::Variable,
        params::Vector{Variable},
    )
        check_vars_params(exprs, [vars; t], params)
        new(exprs, vars, t, params)
    end
end

function Homotopy(
    exprs::AbstractVector{Expression},
    variables::Vector{Variable},
    t::Variable;
    parameters::Vector{Variable} = Variable[],
)
    Homotopy(convert(Vector{Expression}, exprs), variables, t, parameters)
end
function Homotopy(exprs, variables::Vector{Variable}, t::Variable; kwargs...)
    Homotopy(convert(Vector{Expression}, exprs), variables, t; kwargs...)
end

function Base.show(io::IO, H::Homotopy)
    if !get(io, :compact, false)
        println(io, "Homotopy in ", H.t, " of length ", length(H.expressions))
        print(io, " $(length(H.variables)) variables: ", join(H.variables, ", "))
        if !isempty(H.parameters)
            print(io, "\n $(length(H.parameters)) parameters: ", join(H.parameters, ", "))
        end
        print(io, "\n\n")
        for i = 1:length(H)
            print(io, " ", H.expressions[i])
            i < length(H) && print(io, "\n")
        end
    else
        print(io, "[")
        for i = 1:length(H)
            print(io, H.expressions[i])
            i < length(H) && print(io, ", ")
        end
        print(io, "]")
    end
end

evaluate(H::Homotopy, x::AbstractVector, t) =
    evaluate(H.expressions, H.variables => x, H.t => t)
function evaluate(H::Homotopy, x::AbstractVector, t, p::AbstractVector)
    evaluate(H.expressions, H.variables => x, H.t => t, H.parameters => p)
end
(H::Homotopy)(x::AbstractVector, t) = evaluate(H, x, t)
(H::Homotopy)(x::AbstractVector, t, p::AbstractVector) = evaluate(H, x, t, p)

function Base.:(==)(H::Homotopy, G::Homotopy)
    H.expressions == G.expressions &&
        H.variables == G.variables &&
        H.parameters == G.parameters
end

Base.size(H::Homotopy) = (length(H.expressions), length(H.variables))
Base.size(H::Homotopy, i::Integer) = size(H)[i]
Base.length(H::Homotopy) = length(H.expressions)

"""
    nvariables(H::Homotopy)

Returns the number of variables of the given homotopy `H`.
"""
nvariables(H::Homotopy) = length(H.variables)

"""
    nparameters(H::Homotopy)

Returns the number of parameters of the given homotopy `H`.
"""
nparameters(H::Homotopy) = length(H.parameters)

"""
    expressions(H::Homotopy)

Returns the expressions of the given homotopy `H`.
"""
expressions(H::Homotopy) = H.expressions

"""
    variables(H::Homotopy)

Returns the variables of the given homotopy `H`.
"""
variables(H::Homotopy) = H.variables

"""
    parameters(H::Homotopy)

Returns the parameters of the given homotopy `H`.
"""
parameters(H::Homotopy) = H.parameters

"""
    to_smallest_eltype(A::AbstractArray)

Convert an array to the smallest eltype such that all elements still fit.

## Example
```julia
typeof(to_smallest_eltype(Any[2,3])) == Vector{Int}
```
"""
function to_smallest_eltype(A::AbstractArray)
    T = typeof(first(A))
    for a in A
        if !(typeof(a) <: T)
            T = promote_type(T, typeof(a))
        end
    end
    convert.(T, A)
end

function optimize(H::Homotopy)
    Homotopy(horner.(H.expressions, Ref(H.variables)), H.variables, H.t, H.parameters)
end


text:
# # Introduction to rendering using RayTracer.jl
#
# In this example we will render the famous UTAH Teapot model.
# We will go through the entire rendering API. We will load
# an obj file for the scene. This needs to be downloaded manually.
#
# Run this code in your terminal to get the file:
# `wget https://raw.githubusercontent.com/McNopper/OpenGL/master/Binaries/teapot.obj`

# If you are using REPL mode you need the `ImageView.jl` package
using RayTracer, Images #, ImageView

# ## General Attributes of the Scene
#
# Specify the dimensions of the image we want to generate.
# `screen_size` is never passed into the RayTracer directly so it
# need not be a named tuple.
screen_size = (w = 256, h = 256)

# Load the teapot object from an `obj` file. We can also specify
# the scene using primitive objects directly but that becomes a
# bit involved when there are complicated objects in the scene.
scene = load_obj("teapot.obj")

# We shall define a convenience function for rendering and saving
# the images.
# For understanding the parameters passed to the individual functions
# look into the documentations of [`get_primary_rays`](@ref), [`raytrace`](@ref)
# and [`get_image`](@ref)
function generate_render_and_save(cam, light, filename)
    #src # Get the primary rays for the camera
    origin, direction = get_primary_rays(cam)

    #src # Render the scene
    color = raytrace(origin, direction, scene, light, origin, 2)

    #src # This will reshape `color` into the proper dimensions and return
    #src # an RGB image
    img = get_image(color, screen_size...)

    #src # Display the image
    #src # For REPL mode change this to `imshow(img)`
    display(img)

    #src # Save the generated image
    save(filename, img)
end

# ## Understanding the Light and Camera API
#
# ### DistantLight
#
# In this example we will be using the [`DistantLight`](@ref). This king of lighting
# is useful when we want to render a scene in which all parts of the scene
# receive the same intensity of light.
#
# For the DistantLight we need to provide three attributes:
# * Color     - Color of the Light Rays. Must be a Vec3 Object
# * Intensity - Intensity of the Light
# * Direction - The direction of light rays. Again this needs to be a Vec3 Object
#
# ### Camera
#
# We use a perspective view [`Camera`](@ref) Model in RayTracer. Let us look into the
# arguments we need to pass into the Camera constructor.
#
# * LookFrom - The position of the Camera
# * LookAt   - The point in 3D space where the Camera is pointing
# * vup      - The UP vector of the world (typically Vec3(0.0, 1.0, 0.0), i.e. the y-axis)
# * vfov     - Field of View of the Camera
# * Focus    - The focal length of the Camera
# * Width    - Width of the output image
# * Height   - Height of the output image

# ## Rendering Different Views of the Teapot
#
# Now that we know what each argument means let us render the teapot
#
# ### TOP VIEW Render

light = DistantLight(
    Vec3(1.0f0),
    100.0f0,
    Vec3(0.0f0, 1.0f0, 0.0f0)
)

cam = Camera(
    Vec3(1.0f0, 10.0f0, -1.0f0),
    Vec3(0.0f0),
    Vec3(0.0f0, 1.0f0, 0.0f0),
    45.0f0,
    1.0f0,
    screen_size...
)

generate_render_and_save(cam, light, "teapot_top.jpg")

# ```@raw html
# <p align="center">
#     <img width=256 height=256 src="../../assets/teapot_top.jpg">
# </p>
# ```

# ### SIDE VIEW Render

light = DistantLight(
    Vec3(1.0f0),
    100.0f0,
    Vec3(1.0f0, 1.0f0, -1.0f0)
)

cam = Camera(
    Vec3(1.0f0, 2.0f0, -10.0f0),
    Vec3(0.0f0, 1.0f0, 0.0f0),
    Vec3(0.0f0, 1.0f0, 0.0f0),
    45.0f0,
    1.0f0,
    screen_size...
)

generate_render_and_save(cam, light, "teapot_side.jpg")

# ```@raw html
# <p align="center">
#     <img width=256 height=256 src="../../assets/teapot_side.jpg">
# </p>
# ```

# ### FRONT VIEW Render

light = DistantLight(
    Vec3(1.0f0),
    100.0f0,
    Vec3(1.0f0, 1.0f0, 0.0f0)
)

cam = Camera(
    Vec3(10.0f0, 2.0f0, 0.0f0),
    Vec3(0.0f0, 1.0f0, 0.0f0),
    Vec3(0.0f0, 1.0f0, 0.0f0),
    45.0f0,
    1.0f0,
    screen_size...
)

generate_render_and_save(cam, light, "teapot_front.jpg")

# ```@raw html
# <p align="center">
#     <img width=256 height=256 src="../../assets/teapot_front.jpg">
# </p>
# ```

# ## Next Steps
#
# * Try Rendering complex environments with RayTracer
# * Look into the other examples in `examples/`
# * Read about inverse rendering and see the examples on that


text:
abstract MPI_areaConfig
include("MPI_2D_Config.jl")
include("helper_boundaries.jl")


text:
import Unitful.𝐌, Unitful.𝐋, Unitful.𝐓

@testset "Unit Convention" begin
    @test InteratomicPotentials.ENERGY_UNIT == aunit(𝐌 * 𝐋^2 / 𝐓^2)
    @test InteratomicPotentials.FORCE_UNIT == aunit(𝐌 * 𝐋 / 𝐓^2)

    @test InteratomicPotentials.ENERGY_TYPE == typeof(1.0aunit(𝐌 * 𝐋^2 / 𝐓^2))
    @test InteratomicPotentials.FORCE_TYPE == typeof(1.0aunit(𝐌 * 𝐋 / 𝐓^2))
end


text:
## function to perform imaging opetation
##

b3spline1d = [1/16 , 1/4 , 3/8 , 1/4 ,  1/16]

## Wavelet A Trous Algorithm (See e.g. Leon et al. 2001)
## im: image array
## levels : number of Wavelet planes
## kernel1d : 1D kernel to create the 2D kernel. Default is B3-spline but  could be other (e.g. [0,1,0])
##
function atrous(im, levels , kernel1d = b3spline1d, verbose = true)
    let
        w      = []
        ker    = Base.copy(kernel1d)
        approx = Base.copy(im)

        for lev in 1:levels
            if verbose
                println("##WT--A Trous--Plane: $lev")
            end

            kernel2d = zeros(length(ker),length(ker))
            for i in 1:(length(ker))
                kernel2d[i,:] = ker[i] .* ker
            end
            ker2 = ImageFiltering.centered(kernel2d)

            smooth = ImageFiltering.imfilter(approx, ker2)
            wi = approx .- smooth
            approx = smooth

            ker = upscale(ker)
            push!(w, wi)
        end

        ## LSP ###
        if verbose
            println("##WT--A Trous--LSP ($(levels+1))")
        end

        push!(w, approx)

        return(w)
    end
end

## upscale by a factor 2 the 1D kernel
function upscale(kerr)
    xx   = 0:length(kerr)-1
    itp  = Interpolations.interpolate(kerr, Interpolations.BSpline(Interpolations.Linear()))
    sitp = Interpolations.scale(itp, xx)
    xp   = 0:0.5:length(kerr)-1
    ker2 = 0.5 .* sitp(xp)

    return(ker2)
end

### Add the WT coefficients to get the (partial) image
##
function addWav(wt,p1,p2)
    let

        nplane = size(wt)
        if p2 == -1 p2 = nplane[1] end
        nxy = size(wt[1])
        println(nxy)
        im = zeros(nxy[1], nxy[2])

        for i in p1:p2
            im = im .+ wt[i]
        end
    return(im)
    end
end

## Hard thresholding on the WT coefficients.
## wt: wavelet transform
## d : random distribution
## σ : threshold in std from the simulated noise, if an array of length== nw[1] then used for thresholding
##
function thresholdingWav(wt, d::UnivariateDistribution , σ = 3 , kernel1d=b3spline1d)
    wtfilt = Base.copy(wt)
    nw = size(wt)
    nxy = size(wt[1])

    println("##WT--Thresholding the WT")
    if length(σ) == nw[1]
        threshold = σ
    else
        println("##WT--Simulating the noise in the WT.")
        threshold = zeros(nw[1])
        noise = rand(d , nxy[1] , nxy[2])
        wsimul = atrous(noise , nw[1]-1 , kernel1d , false)
        threshold = σ .* noiseWav(wsimul)
    end

    for p in 1:nw[1]
        for i in 1:nxy[1]
            for j in 1:nxy[2]
                if abs(wt[p][i,j]) < threshold[p]
                   wtfilt[p][i,j] = 0
                end
            end
        end
    end
    return(wtfilt)
end

## Compute std per plane for the wt
##
function noiseWav(wt)
    noise = []
    for plane in wt
        σ = std(plane)
        push!(noise,σ)
    end
    return(noise)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
"""
Import packages and files
"""
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

@info("Loading Packages...")
using Zygote, LinearAlgebra, Random
using Flux
using Random
using Plots

include("bayes_nn.jl")
"""
class_data.jl contains the data of the example
for Bayesian Neural Network with Turing language from Turinglang
    https://turing.ml/dev/tutorials/3-bayesnn/.
"""
include("turing_data.jl")    #
plot_data()

X = hcat(xs...)
Y = hcat(ts...)

"""
Set up a bnn for classification. The architecture and standard deviation are the
same as in the example of Turing in order to ensure correctness of the algorithm.
"""
# Create a regularization term and a Gaussain prior variance term.
alpha = 0.09
sig = sqrt(1.0 / alpha)

""" NNet with 2 input units, 3 units in 1st hidden layer,
2 units in 2nd hidden layer and 1 output unit and their correpsonding
activation function.
"""
ffnet = network([2, 3, 2, 1], [identity, tanh, tanh, sigmoid])

"""
Sample from the posterior distribution using Hamiltonian Monte Carlo.
"""
function main(net::network, X, Y, prior_sigma, maxiter, burnin)
    # matrices to store mcmc results
    sampled_weights = []
    sampled_bias = []
    lp = zeros(maxiter)

    for its in 1:maxiter
        hmc_bnn!(net, X, Y, prior_sigma, 0.05, 4)
        lp[its] = log_posterior(net, X, Y, prior_sigma)
        if its > burnin
            push!(sampled_weights, net.weights)
            push!(sampled_bias, net.biases)
        end
        if mod(its, 100) == 0
            println("HMC iterations: $its")
        end
    end

    return sampled_weights, sampled_bias, lp
end

Random.seed!(12345)
@time sampled_weights, sampled_bias, lp = main(ffnet, X, Y, sig, 5000, 1000)

# define test set
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))

"""
Visualise the results
"""
n_end = 4000
anim = @gif for i=1:100:n_end
    plot_data();
    bnn_posterior!(ffnet, sampled_weights[i], sampled_bias[i])
    Z = [nn_forward(ffnet, [x, y]) for x=x_range, y=y_range]
    contour!(x_range, y_range, cell2array(Z), title="Iteration $i", clim = (0,1));
end every 10;

_, idx = findmax(lp)
bnn_posterior!(ffnet, sampled_weights[idx], sampled_bias[idx])
Z = [nn_forward(ffnet, [x, y]) for x=x_range, y=y_range]
plot_data()
contour!(x_range, y_range, cell2array(Z))


text:
using POMDPModels

problem =  SimpleGridWorld()
# e greedy
policy = EpsGreedyPolicy(problem, 0.5)
a = first(actions(problem))
@inferred action(policy, FunctionPolicy(s->a::Symbol), 1, GWPos(1,1))
policy = EpsGreedyPolicy(problem, 0.0)
@test action(policy, FunctionPolicy(s->a), 1, GWPos(1,1)) == a

# softmax
policy = SoftmaxPolicy(problem, 0.5)
@test loginfo(policy, 1).temperature == 0.5
on_policy = ValuePolicy(problem)
@inferred action(policy, on_policy, 1, GWPos(1,1))

# test linear schedule
policy = EpsGreedyPolicy(problem, LinearDecaySchedule(start=1.0, stop=0.0, steps=10))
for i=1:11
    action(policy, FunctionPolicy(s->a), i, GWPos(1,1))
    @test policy.eps(i) < 1.0
    @test loginfo(policy, i).eps == policy.eps(i)
end
@test policy.eps(11) ≈ 0.0

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies