text:
using GymSpaces: Discrete, TupleSpace, AbstractSpace
using Random

deck = vcat(collect(1:10), [10, 10, 10])

@inline draw_card(seed::MersenneTwister) = rand(seed, deck)
@inline draw_hand(seed::MersenneTwister) = [draw_card(seed), draw_card(seed)]
@inline usable_ace(hand) = 1 ∈ hand && sum(hand) + 10 ≤ 21
@inline sum_hand(hand) = usable_ace(hand) ? sum(hand) + 10 : sum(hand)
@inline is_bust(hand) = sum_hand(hand) > 21
@inline score(hand) = is_bust(hand) ? 0 : sum_hand(hand)
@inline is_natural(hand) = sort(hand) == [1, 10]

mutable struct BlackjackEnv <: AbstractEnv
    action_space::AbstractSpace
    observation_space::AbstractSpace
    natural::Bool
    dealer::Array
    player::Array
    seed::MersenneTwister
end

"""Simple blackjack environment

Blackjack is a card game where the goal is to obtain cards that sum to as
near as possible to 21 without going over.  They're playing against a fixed
dealer.
Face cards (Jack, Queen, King) have point value 10.
Aces can either count as 11 or 1, and it's called 'usable' at 11.
This game is placed with an infinite deck (or with replacement).
The game starts with each (player and dealer) having one face up and one
face down card.

The player can request additional cards (hit=1) until they decide to stop
(stick=0) or exceed 21 (bust).

After the player sticks, the dealer reveals their facedown card, and draws
until their sum is 17 or greater.  If the dealer goes bust the player wins.

If neither player nor dealer busts, the outcome (win, lose, draw) is
decided by whose sum is closer to 21.  The reward for winning is +1,
drawing is 0, and losing is -1.

The observation of a 3-tuple of: the players current sum,
the dealer's one showing card (1-10 where 1 is ace),
and whether or not the player holds a usable ace (0 or 1).

This environment corresponds to the version of the blackjack problem
described in Example 5.1 in Reinforcement Learning: An Introduction
by Sutton and Barto.
http://incompleteideas.net/book/the-book-2nd.html
"""
function BlackjackEnv(natural=false)
    action_space = Discrete(2)
    observation_space = TupleSpace([
        Discrete(32),
        Discrete(11),
        Discrete(2)
    ])
    seed = MersenneTwister()
    BlackjackEnv(action_space, observation_space, natural, [], [], seed)
end

seed!(env::BlackjackEnv) = (env.seed = MersenneTwister())
seed!(env::BlackjackEnv, int::Unsigned) = (env.seed = MersenneTwister(int))

@inline get_obs(env::BlackjackEnv) = (sum_hand(env.player), env.dealer[1], usable_ace(env.player))

function reset!(env::BlackjackEnv)
    env.dealer = draw_hand(env.seed)
    env.player = draw_hand(env.seed)
    get_obs(env)
end

function step!(env::BlackjackEnv, action)
    @assert action ∈ env.action_space "Invalid action"

    if action == 1
        push!(env.player, draw_card(env.seed))
        if is_bust(env.player)
            done = true
            reward = -1
        else
            done = false
            reward = 0
        end
    else
        done = true
        while sum_hand(env.dealer) < 17
            push!(env.dealer, draw_card(env.seed))
        end
        reward = Float32(cmp(score(env.player), score(env.dealer)))
        if env.natural && is_natural(env.player) && reward == 1f0
            reward = 1.5f0
        end
    end
    return get_obs(env), reward, done, Dict{Nothing, Nothing}()
end


text:
import Base.vcat, Base.hcat, Base.hvcat
struct HcatAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple
    size::Tuple{Int, Int}

    function HcatAtom(args::AbstractExpr...)
        num_rows = args[1].size[1]
        num_cols = 0
        for arg in args
            if arg.size[1] != num_rows
                error("Cannot horizontally stack expressions of varying number of rows")
            end
            num_cols += arg.size[2]
        end
        children = tuple(args...)
        return new(:hcat, hash(children), children, (num_rows, num_cols))
    end
end

function sign(x::HcatAtom)
    return sum(map(sign, x.children))
end

function monotonicity(x::HcatAtom)
    return [Nondecreasing() for c in x.children]
end

function curvature(x::HcatAtom)
    return ConstVexity()
end

function evaluate(x::HcatAtom)
    return hcat(map(evaluate, x.children)...)
end

function conic_form!(x::HcatAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        # build a list of child conic objectives and constraints
        objectives = ConicObj[]
        for child in x.children
            push!(objectives, conic_form!(child, unique_conic_forms))
        end
        # build a dict from variable ids to sizes
        variable_to_sizes = Dict{UInt64, Int}()
        for objective in objectives
            for id in keys(objective)
                if !(id in keys(variable_to_sizes))
                    if id == objectid(:constant)
                        variable_to_sizes[id] = 1
                    else
                        variable_to_sizes[id] = length(unique_conic_forms.id_to_variables[id])
                    end
                end
            end
        end

        # Suppose the child objectives for two children e1 (2 x 1) and e2 (2 x 2) look something like
        #  e1: x => 1 2 3
        #           4 5 6
        #      y => 2 4
        #           7 8
        #  e2: x => 1 1 1
        #           2 2 2
        #           3 3 3
        #           4 4 4
        # The objective of [e1 e2] will look like
        #            x => 1 2 3
        #                 4 5 6
        #                 1 1 1
        #                 2 2 2
        #                 3 3 3
        #                 4 4 4
        #            y => 2 4
        #                 7 8
        #                 0 0
        #                 0 0
        #                 0 0
        #                 0 0
        # builds the objective by aggregating a list of coefficients for each variable
        # from each child objective, and then vertically concatenating them
        objective = ConicObj()
        for (id, col_size) in variable_to_sizes
            #temp_tuple = Tuple{Value,Value}
            x1_value_list = Value[]
            x2_value_list = Value[]

            for i in 1:length(objectives)
                row_size = length(x.children[i])
                if haskey(objectives[i], id)
                    push!(x1_value_list, objectives[i][id][1])
                    push!(x2_value_list, objectives[i][id][2])
                else
                    push!(x1_value_list, spzeros(row_size, col_size))
                    push!(x2_value_list, spzeros(row_size, col_size))
                end
            end
            x1 = vcat(x1_value_list...)
            x2 = vcat(x2_value_list...)
            objective[id] = (x1,x2)

        end
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

# TODO: fix piracy!

# * `Value` is not owned by Convex.jl
# * splatting creates zero-argument functions, which again are not owned by Convex.jl

hcat(args::AbstractExpr...) = HcatAtom(args...)
hcat(args::AbstractExprOrValue...) = HcatAtom(map(arg -> convert(AbstractExpr, arg), args)...)
hcat(args::Value...) = Base.cat(args..., dims=Val(2))

# TODO: implement vertical concatenation in a more efficient way
vcat(args::AbstractExpr...) = transpose(HcatAtom(map(transpose, args)...))
vcat(args::AbstractExprOrValue...) = transpose(HcatAtom(map(arg -> transpose(convert(AbstractExpr, arg)), args)...))
vcat(args::Value...) = Base.cat(args..., dims=Val(1)) # Note: this makes general vcat slower for anyone using Convex...

function hvcat(rows::Tuple{Vararg{Int}}, args::AbstractExprOrValue...)
    nbr = length(rows)
    rs = Vector{Any}(undef, nbr)
    a = 1
    for i = 1:nbr
        rs[i] = hcat(args[a:a-1+rows[i]]...)
        a += rows[i]
    end
    return vcat(rs...)
end

ext: jl


text:
using CFTime
using Dates
using Formatting
using JLD2
using Distributed
using SharedArrays
using MPI

if !(:LogSystem in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "LogSystem.jl")))
end
using .LogSystem

if !(:ModelClockSystem in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "ModelClockSystem.jl")))
end
using .ModelClockSystem

if !(:ConfigCheck in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "ConfigCheck.jl")))
end
using .ConfigCheck

if !(:appendLine in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "AppendLine.jl")))
end

function runModel(
    OMMODULE      :: Any,
    coupler_funcs :: Any,
    config        :: Union{Dict, Nothing} = nothing,
)

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    ntask = MPI.Comm_size(comm)

    writeLog("===== [ Master Created ] =====")
    writeLog("Number of total tasks       : {:d}", ntask)
    writeLog("Number of total worker tasks: {:d}", ntask-1)

    MPI.Barrier(comm)

    is_master = rank == 0

    if is_master
        if config == nothing
            throw(ErrorException("Master thread needs to provide config parameter."))
        end

        writeLog("Validate driver config.")
        config["DRIVER"] = validateConfigEntries(config["DRIVER"], getDriverConfigDescriptor())
    end

    writeLog("Broadcast config to slaves.")
    config = MPI.bcast(config, 0, comm)

    p = config["DRIVER"]["caserun"]
    writeLog("Setting working directory to {:s}", p)
    if is_master
        if ! isdir(p)
            writeLog("Working directory does not exist. Create it.")
            mkpath(p)
        end
    end
    MPI.Barrier(comm)
    cd(p)

    local t_start = nothing
    local read_restart = nothing

    if is_master

        writeLog("Getting model start time.")
        read_restart, t_start = coupler_funcs.master_before_model_init()

        if read_restart

            writeLog("read_restart is true.")

            restart_info = JLD2.load("model_restart.jld2")
            restart_t_start = restart_info["timestamp"]

            if t_start == nothing
                writeLog("t_start == nothing. Skip checking restart time.")
            elseif t_start != restart_t_start
                throw(ErrorException(format("Model restart time inconsistent! The start time received from the coupler is {:s} while the restart file is {:s}.")))
            else
                writeLog("Model restart time is consistent.")
            end

        end

        writeLog("### Simulation time start: {:s}", Dates.format(t_start, "yyyy-mm-dd HH:MM:SS"))

    end

    writeLog("Broadcast t_start and read_restart to slaves.")
    t_start = MPI.bcast(t_start, 0, comm)
    read_restart = MPI.bcast(read_restart, 0, comm)

    # copy the start time by adding 0 seconds
    beg_datetime = t_start + Dates.Second(0)

    # Construct model clock
    clock = ModelClock("Model", beg_datetime)

    # initializing
    writeLog("===== INITIALIZING MODEL: {:s} =====", OMMODULE.name)

    OMDATA = OMMODULE.init(
        casename     = config["DRIVER"]["casename"],
        clock        = clock,
        config       = config,
        read_restart = read_restart,
    )

    if is_master
        coupler_funcs.master_after_model_init!(OMMODULE, OMDATA)
    end

    writeLog("Ready to run the model.")
    step = 0
    while true

        step += 1

        if is_master
            writeLog("Current time: {:s}", clock2str(clock))
        end

        stage = nothing
        Δt    = nothing
        write_restart = nothing

        if is_master
            stage, Δt, write_restart = coupler_funcs.master_before_model_run!(OMMODULE, OMDATA)
        end
        stage         = MPI.bcast(stage, 0, comm)
        Δt            = MPI.bcast(Δt, 0, comm)
        write_restart = MPI.bcast(write_restart, 0, comm)

        if stage == :RUN

            cost = @elapsed let

                # This advanced option is needed when deriving
                # qflux in direct method: Coupler needs to load
                # the initial data of each day and force master
                # to sync thermo variables with slaves.
                if config["DRIVER"]["compute_QFLX_direct_method"]
                    OMMODULE.syncM2S!(OMDATA)
                end

                OMMODULE.run!(
                    OMDATA;
                    Δt = Δt,
                    write_restart = write_restart,
                )
                MPI.Barrier(comm)

            end

            writeLog("Computation cost: {:.2f} secs.", cost)

            if write_restart && is_master

                driver_restart_file = "model_restart.jld2"
                writeLog("Writing restart time of driver: $(driver_restart_file)")
                JLD2.save(driver_restart_file, "timestamp", clock.time)

                archive_list_file = joinpath(
                    config["DRIVER"]["caserun"],
                    config["DRIVER"]["archive_list"],
                )

                timestamp_str = format(
                    "{:s}-{:05d}",
                    Dates.format(clock.time, "yyyy-mm-dd"),
                    floor(Int64, Dates.hour(clock.time)*3600+Dates.minute(clock.time)*60+Dates.second(clock.time)),
                )

                appendLine(archive_list_file,
                    format("cp,{},{},{}",
                        driver_restart_file,
                        config["DRIVER"]["caserun"],
                        joinpath(config["DRIVER"]["archive_root"], "rest", timestamp_str)
                    )
                )

            end

            if is_master
                coupler_funcs.master_after_model_run!(OMMODULE, OMDATA)
            end

            if config["DRIVER"]["compute_QFLX_direct_method"]
                writeLog("compute_QFLX_direct_method is true")
                OMMODULE.syncM2S!(OMDATA)
            end

            if is_master
                advanceClock!(clock, Δt)
                dropRungAlarm!(clock)
            end

            # Broadcast time to workers. Workers need time
            # because datastream needs time interpolation.
            _time = MPI.bcast(clock.time, 0, comm)
            if !is_master
                setClock!(clock, _time)
            end

        elseif stage == :END

            writeLog("Receive :END. Entering finalizing phase now.")
            break

        else

            throw(ErrorException("Unknown stage : " * string(stage)))

        end

    end

    if is_master

        OMMODULE.final(OMDATA)
        coupler_funcs.master_finalize!(OMMODULE, OMDATA)

        archive(joinpath(
            config["DRIVER"]["caserun"],
            config["DRIVER"]["archive_list"],
        ))

    end

    writeLog("Program Ends.")

end

function getDriverConfigDescriptor()

    return [
            ConfigEntry(
                "casename",
                :required,
                [String,],
            ),

            ConfigEntry(
                "caseroot",
                :required,
                [String,],
            ),

            ConfigEntry(
                "caserun",
                :required,
                [String,],
            ),

            ConfigEntry(
                "archive_root",
                :required,
                [String,],
            ),

            ConfigEntry(
                "archive_list",
                :optional,
                [String,],
                "archive_list.txt",
            ),

            ConfigEntry(
                "compute_QFLX_direct_method",
                :optional,
                [Bool,],
                false,
            ),

   ]
end

function archive(
    archive_list_file :: String,
)

    println("===== Archiving files BEGIN =====")

    for line in eachline(archive_list_file)

        args = split(line, ",")

        if length(args) == 0
            continue
        end

        action = args[1]
        args = args[2:end]

        if action in ["mv", "cp"]

            fname, src_dir, dst_dir = args

            if ! isdir(dst_dir)
                mkpath(dst_dir)
            end

            src_file = joinpath(src_dir, fname)
            dst_file = joinpath(dst_dir, fname)

            if isfile(src_file)

                if action == "mv"
                    mv(src_file, dst_file, force=true)
                    println(format("Moving file: {:s} ( {:s} => {:s} )", fname, src_dir, dst_dir))
                elseif action == "cp"
                    cp(src_file, dst_file, force=true)
                    println(format("Copying file: {:s} ( {:s} => {:s} )", fname, src_dir, dst_dir))
                end

            else
                println("File does not exist: ", src_file)
            end

        elseif action == "rm"
            fname, fdir = args
            rm(joinpath(fdir, fname), force=true)
            println(format("Removing file: {:s} in {:s}", fname, fdir))
        else
            throw(ErrorException(format("Unknown action in archive list: {:s}", action)))
        end

    end

    println("===== Archiving files END =====")

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
@testset "derivatives" begin
    @polyvar x y z

    @test @inferred(differentiate(x, x)) == 1
    @test @inferred(differentiate(x, y)) == 0
    @test @inferred(differentiate(y, x)) == 0
    @test @inferred(differentiate(x^2, x)) == 2x
    @test @inferred(differentiate(x^2, y)) == 0
    @test @inferred(differentiate(x^2 * y^3, y)) == 3 * x^2 * y^2
    @test @inferred(differentiate(x^2 * y^3, x)) == 2 * x * y^3
    @test @inferred(differentiate(x^2 * y^3, z)) == 0
    @test @inferred(differentiate(3x^2, x)) == 6x
    @test @inferred(differentiate(3x^2 * y^0, y)) == 0
    @test @inferred(differentiate(5 * x^2 + 2 * x + 1, x)) == 10x + 2

    @test @inferred(exponents(differentiate(x^0,x))==(0,))
    @test @inferred(exponents(differentiate(x^0*y*z^2,x))==(0,1,2))

    m = x^2
    @test @wrappedallocs(differentiate(m, x)) == 0
    @test @wrappedallocs(differentiate(m, y)) == 0
    m = x^2 * y^3
    @test @wrappedallocs(differentiate(m, x)) == 0
    @test @wrappedallocs(differentiate(m, y)) == 0

    @test @inferred(differentiate(1, x)) == 0
end


text:
using ITensors
using LinearAlgebra
using Random
using Test

Random.seed!(1234)

@testset "BlockSparse ITensor" begin
  @testset "Constructor" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = ITensor(QN(0), i, dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2
  end

  @testset "Construct from Array" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")

    A = [
      1.0 0.0 0.0
      0.0 2.0 3.0
      0.0 1e-10 4.0
    ]
    T = ITensor(A, i', dag(i))
    @test flux(T) == QN(0)
    @test nnzblocks(T) == 2
    @test Block(1, 1) in nzblocks(T)
    @test Block(2, 2) in nzblocks(T)
    @test T[1, 1] == 1.0
    @test T[2, 2] == 2.0
    @test T[2, 3] == 3.0
    @test T[3, 2] == 1e-10
    @test T[3, 3] == 4.0

    T = ITensor(A, i', dag(i))
    @test flux(T) == QN(0)
    @test nnzblocks(T) == 2
    @test Block(1, 1) in nzblocks(T)
    @test Block(2, 2) in nzblocks(T)
    @test T[1, 1] == 1.0
    @test T[2, 2] == 2.0
    @test T[2, 3] == 3.0
    @test T[3, 2] == 1e-10
    @test T[3, 3] == 4.0

    T = ITensor(A, i', dag(i); tol=1e-9)
    @test flux(T) == QN(0)
    @test nnzblocks(T) == 2
    @test Block(1, 1) in nzblocks(T)
    @test Block(2, 2) in nzblocks(T)
    @test T[1, 1] == 1.0
    @test T[2, 2] == 2.0
    @test T[2, 3] == 3.0
    @test T[3, 2] == 0.0
    @test T[3, 3] == 4.0

    A = [
      1e-9 0.0 0.0
      0.0 2.0 3.0
      0.0 1e-10 4.0
    ]
    T = ITensor(A, i', dag(i); tol=1e-8)
    @test flux(T) == QN(0)
    @test nnzblocks(T) == 1
    @test Block(2, 2) in nzblocks(T)
    @test T[1, 1] == 0.0
    @test T[2, 2] == 2.0
    @test T[2, 3] == 3.0
    @test T[3, 2] == 0.0
    @test T[3, 3] == 4.0

    A = [
      1e-9 2.0 3.0
      1e-9 1e-10 2e-10
      2e-9 1e-10 4e-10
    ]
    T = ITensor(A, i', dag(i); tol=1e-8)
    @test flux(T) == QN(-1)
    @test nnzblocks(T) == 1
    @test Block(1, 2) in nzblocks(T)
    @test T[1, 1] == 0.0
    @test T[1, 2] == 2.0
    @test T[1, 3] == 3.0
    @test T[2, 2] == 0.0
    @test T[2, 3] == 0.0
    @test T[3, 2] == 0.0
    @test T[3, 3] == 0.0

    A = [
      1e-9 2.0 3.0
      1e-5 1e-10 2e-10
      2e-9 1e-10 4e-10
    ]
    @test_throws ErrorException ITensor(A, i', dag(i); tol=1e-8)
  end

  @testset "Construct from Array regression test" begin
    i = Index([QN(0) => 2, QN(1) => 2])
    T = itensor([0, 0, 1, 2], i)
    @test flux(T) == QN(1)
    @test nnzblocks(T) == 1
    @test !(Block(1) in nzblocks(T))
    @test Block(2) in nzblocks(T)
    @test T[1] == 0
    @test T[2] == 0
    @test T[3] == 1
    @test T[4] == 2
    @test flux(T, 1) == QN(0)
    @test flux(T, 2) == QN(0)
    @test flux(T, 3) == QN(1)
    @test flux(T, 4) == QN(1)
    @test_throws BoundsError flux(T, 5)
    @test_throws BoundsError flux(T, 0)
    @test flux(T, Block(1)) == QN(0)
    @test flux(T, Block(2)) == QN(1)
    @test_throws BoundsError flux(T, Block(0))
    @test_throws BoundsError flux(T, Block(3))
  end

  @testset "QN ITensor Array constructor view behavior" begin
    d = 2
    i = Index([QN(0) => d ÷ 2, QN(1) => d ÷ 2])

    # no view
    A = diagm(randn(Float64, d))
    T = ITensor(A, i', dag(i); tol=1e-12)
    @test storage(T) isa NDTensors.BlockSparse{Float64}
    A[1, 1] = 2.0
    T[1, 1] ≠ 2.0

    # no view
    A = diagm(rand(Int, d))
    T = ITensor(Int, A, i', dag(i); tol=1e-12)
    @test storage(T) isa NDTensors.BlockSparse{Int}
    A[1, 1] = 2
    T[1, 1] ≠ 2

    # no view
    A = diagm(rand(Int, d))
    T = ITensor(A, i', dag(i); tol=1e-12)
    @test storage(T) isa NDTensors.BlockSparse{Float64}
    A[1, 1] = 2
    T[1, 1] ≠ 2

    # no view
    A = diagm(randn(Float64, d))
    T = ITensor(A, i', dag(i); tol=1e-12)
    @test storage(T) isa NDTensors.BlockSparse{Float64}
    A[1, 1] = 2
    T[1, 1] ≠ 2

    # no view
    A = diagm(rand(Int, d))
    T = ITensor(Int, A, i', dag(i); tol=1e-12)
    @test storage(T) isa NDTensors.BlockSparse{Int}
    A[1, 1] = 2
    T[1, 1] ≠ 2

    # no view
    A = diagm(rand(Int, d))
    T = ITensor(A, i', dag(i); tol=1e-12)
    @test storage(T) isa NDTensors.BlockSparse{Float64}
    A[1, 1] = 2
    T[1, 1] ≠ 2
  end

  @testset "Constructor Leads to No Blocks" begin
    i = Index(QN(0) => 2, QN(1) => 3; tags="i")
    j = Index(QN(1) => 2, QN(2) => 1; tags="j")
    A = ITensor(i, j)
    @test storage(A) isa NDTensors.EmptyStorage
    @test_throws ErrorException ITensor(QN(0), i, j)
  end

  @testset "ITensor iteration" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = randomITensor(i, dag(j))
    Is = eachindex(A)
    @test length(Is) == dim(A)
    sumA = 0.0
    for I in Is
      sumA += A[I]
    end
    @test sumA ≈ sum(ITensors.data(A))
    sumA = 0.0
    for a in A
      sumA += a
    end
    @test sumA ≈ sum(A)
  end

  @testset "Constructor (from Tuple)" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = ITensor(QN(0), (i, dag(j)))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2
  end

  @testset "Constructor (no flux specified)" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = ITensor(i, dag(j))

    @test flux(A) === nothing
    @test nnzblocks(A) == 0
  end

  @testset "Constructor (Tuple, no flux specified)" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = ITensor((i, dag(j)))

    @test flux(A) === nothing
    @test nnzblocks(A) == 0
  end

  @testset "No indices getindex" begin
    T = ITensor(QN())
    @test order(T) == 0
    @test flux(T) == nothing
    @test nnzblocks(T) == 1
    @test T[] == 0

    s = Index(QN(-1) => 1, QN(1) => 1)
    A = emptyITensor(s, dag(s'))
    B = emptyITensor(s', dag(s))
    A[1, 1] = 1
    B[2, 2] = 1
    C = A * B
    @test order(C) == 0
    @test flux(C) == nothing
    @test nnzblocks(C) == 0
    @test C[] == 0
  end

  @testset "Empty constructor" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")

    A = emptyITensor(i, dag(i'))

    @test nnzblocks(A) == 0
    @test nnz(A) == 0
    @test hasinds(A, i, i')
    @test isnothing(flux(A))

    A[i => 1, i' => 1] = 1.0

    @test nnzblocks(A) == 1
    @test nnz(A) == 1
    @test flux(A) == QN(0)

    A[i => 2, i' => 2] = 1.0

    @test nnzblocks(A) == 2
    @test nnz(A) == 5
    @test flux(A) == QN(0)
  end

  @testset "Check flux when setting elements" begin
    i = Index(QN(0) => 1, QN(1) => 1; tags="i")
    A = randomITensor(QN(0), i, dag(i'))
    @test_throws ErrorException A[i => 1, i' => 2] = 1.0
  end

  @testset "Random constructor" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = randomITensor(QN(1), i, dag(j))

    @test flux(A) == QN(1)
    @test nnzblocks(A) == 1

    B = randomITensor(i, dag(j))

    @test flux(B) == QN()
    @test nnzblocks(B) == 2
  end

  @testset "Complex Number Operations" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    A = randomITensor(ComplexF64, QN(0), i, dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2

    rA = real(A)
    iA = imag(A)
    @test nnzblocks(rA) == nnzblocks(A)
    @test nnzblocks(iA) == nnzblocks(A)
    @test norm(rA + 1im * iA - A) < 1E-8
    @test eltype(rA) == Float64
    @test eltype(iA) == Float64

    cA = conj(A)
    @test eltype(cA) == ComplexF64
    @test norm(cA) ≈ norm(A)

    B = randomITensor(Float64, QN(0), i, dag(j))

    cB = conj(B)
    @test eltype(cB) == Float64
    @test norm(cB) ≈ norm(B)
  end

  @testset "QN onehot" begin
    i = Index(QN(0) => 2, QN(1) => 2; tags="i")

    T = onehot(i => 1)
    @test T[i => 1] ≈ 1.0
    @test T[i => 2] ≈ 0.0
    @test T[i => 3] ≈ 0.0
    @test T[i => 4] ≈ 0.0

    T = onehot(i => 2)
    @test T[i => 1] ≈ 0.0
    @test T[i => 2] ≈ 1.0
    @test T[i => 3] ≈ 0.0
    @test T[i => 4] ≈ 0.0
  end

  @testset "setindex!" begin
    @testset "Test 1" begin
      s1 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s1")
      s2 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s2")
      A = emptyITensor(s1, s2)

      @test nnzblocks(A) == 0
      @test nnz(A) == 0
      @test hasinds(A, s1, s2)
      @test isnothing(flux(A))

      A[2, 1] = 1.0 / sqrt(2)

      @test nnzblocks(A) == 1
      @test nnz(A) == 1
      @test A[s1 => 2, s2 => 1] ≈ 1.0 / sqrt(2)
      @test flux(A) == QN("N", 1, -1)

      A[1, 2] = 1.0 / sqrt(2)

      @test nnzblocks(A) == 2
      @test nnz(A) == 2
      @test A[s1 => 2, s2 => 1] ≈ 1.0 / sqrt(2)
      @test A[s1 => 1, s2 => 2] ≈ 1.0 / sqrt(2)
      @test flux(A) == QN("N", 1, -1)
    end

    @testset "Test 2" begin
      s1 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s1")
      s2 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s2")
      A = emptyITensor(s1, s2)

      @test nnzblocks(A) == 0
      @test nnz(A) == 0
      @test hasinds(A, s1, s2)
      @test isnothing(flux(A))

      A[1, 2] = 1.0 / sqrt(2)

      @test nnzblocks(A) == 1
      @test nnz(A) == 1
      @test A[s1 => 1, s2 => 2] ≈ 1.0 / sqrt(2)
      @test flux(A) == QN("N", 1, -1)

      A[2, 1] = 1.0 / sqrt(2)

      @test nnzblocks(A) == 2
      @test nnz(A) == 2
      @test A[s1 => 2, s2 => 1] ≈ 1.0 / sqrt(2)
      @test A[s1 => 1, s2 => 2] ≈ 1.0 / sqrt(2)
      @test flux(A) == QN("N", 1, -1)
    end
  end

  @testset "Multiply by scalar" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4], "j")

    A = randomITensor(QN(0), i, dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2

    B = 2 * A

    @test flux(B) == QN(0)
    @test nnzblocks(B) == 2

    for ii in dim(i), jj in dim(j)
      @test 2 * A[i => ii, j => jj] == B[i => ii, j => jj]
    end
  end

  @testset "Check arrows when summing" begin
    s = siteinds("S=1/2", 4; conserve_qns=true)
    Tout = randomITensor(QN("Sz" => 2), s[2], s[1], s[3], s[4])
    Tin = randomITensor(QN("Sz" => 2), dag(s[1]), dag(s[2]), dag(s[3]), dag(s[4]))
    @test norm(Tout - Tout) < 1E-10 # this is ok
    @test_throws ErrorException (Tout + Tin)         # not ok
  end

  @testset "Copy" begin
    s = Index([QN(0) => 1, QN(1) => 1], "s")
    T = randomITensor(QN(0), s, s')
    cT = copy(T)
    for ss in dim(s), ssp in dim(s')
      @test T[s => ss, s' => ssp] == cT[s => ss, s' => ssp]
    end
  end

  @testset "Permute" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4], "j")

    A = randomITensor(QN(1), i, dag(j))

    @test flux(A) == QN(1)
    @test nnzblocks(A) == 1

    B = permute(A, j, i)

    @test flux(B) == QN(1)
    @test nnzblocks(B) == 1

    for ii in dim(i), jj in dim(j)
      @test A[ii, jj] == B[jj, ii]
    end
  end

  @testset "Contraction" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4], "j")

    A = randomITensor(QN(0), i, dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2

    B = randomITensor(QN(1), j, dag(i)')

    @test flux(B) == QN(1)
    @test nnzblocks(B) == 1

    C = A * B

    @test hasinds(C, i, i')
    @test flux(C) == QN(1)
    @test nnzblocks(C) == 1
  end

  @testset "Combine and uncombine" begin
    @testset "Combine no indices" begin
      i1 = Index([QN(0, 2) => 2, QN(1, 2) => 2], "i1")
      A = randomITensor(QN(), i1, dag(i1'))

      C = combiner()
      c = combinedind(C)
      @test isnothing(c)
      AC = A * C
      @test nnz(AC) == nnz(A)
      @test nnzblocks(AC) == nnzblocks(A)
      @test hassameinds(AC, A)
      @test norm(AC - A * C) ≈ 0.0
      Ap = AC * dag(C)
      @test nnz(Ap) == nnz(A)
      @test nnzblocks(Ap) == nnzblocks(A)
      @test hassameinds(Ap, A)
      @test norm(A - Ap) ≈ 0.0
    end

    @testset "Combine set direction" begin
      i1 = Index([QN(0) => 2, QN(1) => 3], "i1")
      A = randomITensor(i1', dag(i1))
      @test isnothing(ITensors.checkflux(A))
      C = combiner(dag(i1); dir=ITensors.Out)
      c = combinedind(C)
      @test dir(c) == ITensors.Out
      AC = A * C
      @test nnz(AC) == nnz(A)
      @test nnzblocks(AC) == nnzblocks(A)
      @test isnothing(ITensors.checkflux(AC))
      Ap = AC * dag(C)
      @test nnz(Ap) == nnz(A)
      @test nnzblocks(Ap) == nnzblocks(A)
      @test hassameinds(Ap, A)
      @test isnothing(ITensors.checkflux(AC))
      @test A ≈ Ap
    end

    @testset "Order 2 (IndexSet constructor)" begin
      i1 = Index([QN(0, 2) => 2, QN(1, 2) => 2], "i1")
      A = randomITensor(QN(), i1, dag(i1'))

      iss = [i1, dag(i1'), (i1, dag(i1')), (dag(i1'), i1)]

      for is in iss
        C = combiner(IndexSet(is); tags="c")
        AC = A * C
        @test nnz(AC) == nnz(A)
        Ap = AC * dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A - Ap) ≈ 0.0
      end
    end

    @testset "Order 2" begin
      i1 = Index([QN(0, 2) => 2, QN(1, 2) => 2], "i1")
      A = randomITensor(QN(), (i1, dag(i1')))

      iss = [i1, dag(i1'), (i1, dag(i1')), (dag(i1'), i1)]

      for is in iss
        C = combiner(is; tags="c")
        AC = A * C
        @test nnz(AC) == nnz(A)
        Ap = AC * dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A - Ap) ≈ 0.0
      end
    end

    @testset "Order 3, Combine 2" begin
      i = Index([QN(0) => 2, QN(1) => 2], "i")

      A = randomITensor(QN(0), i, dag(i)', dag(i)'')

      C = combiner(i, dag(i)'')
      c = combinedind(C)

      AC = A * C

      @test hasinds(AC, c, i')
      @test nnz(AC) == nnz(A)

      for b in nzblocks(AC)
        @test flux(AC, b) == QN(0)
      end

      B = ITensor(QN(0), i', c)
      @test nnz(B) == nnz(AC)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC * dag(C)

      @test norm(A - Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
      @test hassameinds(A, Ap)
    end

    @testset "Order 3" begin
      i1 = Index([QN(0, 2) => 2, QN(1, 2) => 2], "i1")
      i2 = settags(i1, "i2")
      A = randomITensor(QN(), i1, i2, dag(i1'))

      iss = [
        i1,
        i2,
        dag(i1'),
        (i1, i2),
        (i2, i1),
        (i1, dag(i1')),
        (dag(i1'), i1),
        (i2, dag(i1')),
        (dag(i1'), i2),
        (i1, i2, dag(i1')),
        (i1, dag(i1'), i2),
        (i2, i1, dag(i1')),
        (i2, dag(i1'), i1),
        (dag(i1'), i1, i2),
        (dag(i1'), i2, i1),
      ]

      for is in iss
        C = combiner(is; tags="c")
        AC = A * C
        @test nnz(AC) == nnz(A)
        Ap = AC * dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A - AC * dag(C)) ≈ 0.0
      end
    end

    @testset "Order 4" begin
      i1 = Index([QN(0, 2) => 2, QN(1, 2) => 2], "i1")
      i2 = settags(i1, "i2")
      A = randomITensor(QN(), i1, i2, dag(i1'), dag(i2'))

      iss = [
        i1,
        i2,
        dag(i1'),
        dag(i2'),
        (i1, i2),
        (i2, i1),
        (i1, dag(i1')),
        (dag(i1'), i1),
        (i1, dag(i2')),
        (dag(i2'), i1),
        (i2, dag(i1')),
        (dag(i1'), i2),
        (i2, dag(i2')),
        (dag(i2'), i2),
        (dag(i1'), dag(i2')),
        (dag(i2'), dag(i1')),
        (i1, i2, dag(i1')),
        (i1, dag(i1'), i2),
        (i2, i1, dag(i1')),
        (i2, dag(i1'), i1),
        (dag(i1'), i1, i2),
        (dag(i1'), i2, i1),
        (i1, dag(i1'), dag(i2')),
        (i1, dag(i2'), dag(i1')),
        (dag(i1'), i1, dag(i2')),
        (dag(i1'), dag(i2'), i1),
        (dag(i2'), i1, dag(i1')),
        (dag(i2'), dag(i1'), i1),
        (i1, i2, dag(i1'), dag(i2')),
        (i1, i2, dag(i2'), dag(i1')),
        (i1, dag(i1'), i2, dag(i2')),
        (i1, dag(i1'), dag(i2'), i2),
        (i1, dag(i2'), i2, dag(i1')),
        (i1, dag(i2'), dag(i1'), i2),
        (i2, i1, dag(i1'), dag(i2')),
        (i2, i1, dag(i2'), dag(i1')),
        (i2, dag(i1'), i1, dag(i2')),
        (i2, dag(i1'), dag(i2'), i1),
        (i2, dag(i2'), i1, dag(i1')),
        (i2, dag(i2'), dag(i1'), i1),
        (dag(i1'), i2, i1, dag(i2')),
        (dag(i1'), i2, dag(i2'), i1),
        (dag(i1'), i1, i2, dag(i2')),
        (dag(i1'), i1, dag(i2'), i2),
        (dag(i1'), dag(i2'), i2, i1),
        (dag(i1'), dag(i2'), i1, i2),
        (dag(i2'), i1, dag(i1'), i2),
        (dag(i2'), i1, i2, dag(i1')),
        (dag(i2'), dag(i1'), i1, i2),
        (dag(i2'), dag(i1'), i2, i1),
        (dag(i2'), i2, i1, dag(i1')),
        (dag(i2'), i2, dag(i1'), i1),
      ]

      for is in iss
        C = combiner(is; tags="c")
        AC = A * C
        @test nnz(AC) == nnz(A)
        Ap = AC * dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A - Ap) ≈ 0.0
      end
    end

    @testset "Order 4, Combine 2, Example 1" begin
      s1 = Index(
        [
          QN(("Sz", 0), ("Nf", 0)) => 1,
          QN(("Sz", +1), ("Nf", 1)) => 1,
          QN(("Sz", -1), ("Nf", 1)) => 1,
          QN(("Sz", 0), ("Nf", 2)) => 1,
        ],
        "site,n=1",
      )
      s2 = replacetags(s1, "n=1", "n=2")

      A = randomITensor(QN(), s1, s2, dag(s1)', dag(s2)')

      C = combiner(dag(s1)', dag(s2)')
      c = combinedind(C)

      AC = A * C

      @test norm(AC) ≈ norm(A)
      @test hasinds(AC, s1, s2, c)
      @test nnz(AC) == nnz(A)
      for b in nzblocks(AC)
        @test flux(AC, b) == QN()
      end

      @test nnzblocks(AC) < nnz(A)

      B = ITensor(QN(), s1, s2, c)
      @test nnz(B) == nnz(AC)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC * dag(C)

      @test hassameinds(A, Ap)
      @test norm(A - Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
    end

    @testset "Order 4, Combine 2, Example 2" begin
      s1 = Index([QN(("Nf", 0)) => 1, QN(("Nf", 1)) => 1], "site,n=1")
      s2 = replacetags(s1, "n=1", "n=2")

      A = randomITensor(QN(), dag(s2)', s2, dag(s1)', s1)

      C = combiner(dag(s2)', dag(s1)')
      c = combinedind(C)

      AC = A * C

      @test norm(AC) ≈ norm(A)
      @test hasinds(AC, s1, s2, c)
      @test nnz(AC) == nnz(A)
      for b in nzblocks(AC)
        @test flux(AC, b) == QN()
      end

      B = ITensor(QN(), s1, s2, c)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC * dag(C)

      @test hassameinds(A, Ap)
      @test norm(A - Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
    end

    @testset "Order 4, Combine 2, Example 3" begin
      s1 = Index([QN(("Nf", 0)) => 1, QN(("Nf", 1)) => 1], "site,n=1")
      s2 = replacetags(s1, "n=1", "n=2")

      A = randomITensor(QN(), dag(s1)', s2, dag(s2)', s1)

      C = combiner(dag(s2)', dag(s1)')
      c = combinedind(C)

      AC = A * C

      @test norm(AC) ≈ norm(A)
      @test hasinds(AC, s1, s2, c)
      @test nnz(AC) == nnz(A)
      for b in nzblocks(AC)
        @test flux(AC, b) == QN()
      end

      B = ITensor(QN(), s1, s2, c)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC * dag(C)

      @test hassameinds(A, Ap)
      @test norm(A - Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
    end
  end

  @testset "Check that combiner commutes" begin
    i = Index(QN(0, 2) => 2, QN(1, 2) => 2; tags="i")
    j = settags(i, "j")
    A = randomITensor(QN(0, 2), i, j, dag(i'), dag(j'))
    C = combiner(i, j)
    @test norm(A * dag(C') * C - A * C * dag(C')) ≈ 0.0
  end

  @testset "Combiner for block deficient ITensor" begin
    i = Index(QN(0, 2) => 2, QN(1, 2) => 2; tags="i")
    j = settags(i, "j")
    A = emptyITensor(i, j, dag(i'))
    A[1, 1, 1] = 1.0
    C = combiner(i, j; tags="c")
    AC = A * C
    Ap = AC * dag(C)
    @test norm(A - Ap) ≈ 0.0
    @test norm(Ap - A) ≈ 0.0
  end

  @testset "Combine Complex ITensor" begin
    s1 = Index(QN(("Sz", 1)) => 1, QN(("Sz", -1)) => 1; tags="S=1/2,Site,n=1")
    s2 = Index(QN(("Sz", 1)) => 1, QN(("Sz", -1)) => 1; tags="S=1/2,Site,n=2")

    T = randomITensor(ComplexF64, QN("Sz", 0), s1, s2)

    C = combiner(s1, s2)
    CT = C * T
    @test norm(CT) ≈ norm(T)
    TT = dag(C) * CT
    @test TT ≈ T
  end

  @testset "Combiner bug #395" begin
    i1 = Index([QN(0) => 1, QN(1) => 2], "i1")
    i2 = Index([QN(0) => 1, QN(1) => 2], "i2")
    A = randomITensor(QN(), i1, i2, dag(i1)', dag(i2)')
    CL = combiner(i1, i2)
    CR = combiner(dag(i1)', dag(i2)')
    AC = A * CR * CL
    @test AC * dag(CR) * dag(CL) ≈ A
  end

  @testset "Contract to scalar" begin
    i = Index([QN(0) => 1, QN(1) => 1], "i")
    A = randomITensor(QN(0), i, dag(i'))

    c = A * dag(A)

    @test nnz(c) == 1
    @test nnzblocks(c) == 1
    @test c[] isa Float64
    @test c[] ≈ norm(A)^2
  end

  @testset "eigen" begin
    @testset "eigen hermitian" begin
      i = Index(QN(0) => 2, QN(1) => 3, QN(2) => 4; tags="i")
      j = settags(i, "j")
      k = settags(i, "k")
      l = settags(i, "l")

      A = randomITensor(QN(), i, j, dag(k), dag(l))
      A = A * prime(dag(A), (i, j))

      F = eigen(A; ishermitian=true, tags="x")

      D, U = F
      Ut = F.Vt

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(D) isa NDTensors.DiagBlockSparse

      u = commonind(D, U)
      up = uniqueind(D, U)

      @test hastags(u, "x")
      @test plev(u) == 0
      @test hastags(up, "x")
      @test plev(up) == 1

      @test hassameinds(U, (i, j, u))
      @test hassameinds(D, (u, up))

      @test A ≈ dag(U) * D * U' atol = 1e-11
      @test A ≈ dag(U) * D * Ut atol = 1e-11
      @test A * U ≈ U' * D atol = 1e-11
      @test A * U ≈ Ut * D atol = 1e-11
    end

    @testset "eigen hermitian (truncate)" begin
      i = Index(QN(0) => 2, QN(1) => 3, QN(2) => 4; tags="i")
      j = settags(i, "j")
      k = settags(i, "k")
      l = settags(i, "l")

      A = randomITensor(QN(), i, j, dag(k), dag(l))
      A = A * prime(dag(A), (i, j))
      for i in 1:4
        A = mapprime(A * A', 2, 1)
      end
      A = A / norm(A)

      cutoff = 1e-5
      F = eigen(A; ishermitian=true, tags="x", cutoff=cutoff)

      D, U, spec = F
      Ut = F.Vt

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(D) isa NDTensors.DiagBlockSparse

      u = commonind(D, U)
      up = uniqueind(D, U)

      @test hastags(u, "x")
      @test plev(u) == 0
      @test hastags(up, "x")
      @test plev(up) == 1

      @test hassameinds(U, (i, j, u))
      @test hassameinds(D, (u, up))

      for b in nzblocks(A)
        @test flux(A, b) == QN(0)
      end
      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(D)
        @test flux(D, b) == QN(0)
      end

      Ap = dag(U) * D * U'

      @test norm(Ap - A) ≤ 1e-2
      @test norm(dag(U) * D * Ut - A) ≤ 1e-2
      @test minimum(dims(D)) == length(spec.eigs)
      @test minimum(dims(D)) < dim(i) * dim(j)

      @test spec.truncerr ≤ cutoff
      err = sqrt(1 - (Ap * dag(Ap))[] / (A * dag(A))[])
      @test err ≤ cutoff
      @test err ≈ spec.truncerr rtol = 3e-1
    end

    @testset "eigen non-hermitian" begin
      i = Index(QN(0) => 2, QN(1) => 3, QN(2) => 4; tags="i")
      j = settags(i, "j")

      A = randomITensor(QN(), i, j, dag(i'), dag(j'))

      F = eigen(A; tags="x")

      D, U = F
      Ut = F.Vt

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(D) isa NDTensors.DiagBlockSparse

      u = commonind(D, U)
      up = uniqueind(D, U)

      @test hastags(u, "x")
      @test plev(u) == 0
      @test hastags(up, "x")
      @test plev(up) == 1

      @test A ≉ U' * D * dag(U) atol = 1e-12
      @test A ≉ Ut * D * dag(U) atol = 1e-12
      @test A * U ≈ U' * D atol = 1e-12
      @test A * U ≈ Ut * D atol = 1e-12
    end

    @testset "eigen non-hermitian (general inds)" begin
      i = Index(QN(0) => 2, QN(1) => 3, QN(2) => 4; tags="i")
      j = settags(i, "j")
      ĩ, j̃ = sim(i), sim(j)

      A = randomITensor(QN(), i, j, dag(ĩ), dag(j̃))

      F = eigen(A, (i, j), (ĩ, j̃); lefttags="x", righttags="y")

      D, U = F
      Ut = F.Vt

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(D) isa NDTensors.DiagBlockSparse

      l = uniqueind(D, U)
      r = commonind(D, U)

      @test F.l == l
      @test F.r == r

      @test hastags(l, "x")
      @test plev(l) == 0
      @test hastags(r, "y")
      @test plev(r) == 0

      @test hassameinds(U, (ĩ, j̃, r))
      @test hassameinds(Ut, (i, j, l))

      @test A * U ≈ Ut * D atol = 1e-12
      @test A ≉ Ut * D * dag(U) atol = 1e-12
    end

    @testset "eigen mixed arrows" begin
      i1 = Index([QN(0) => 1, QN(1) => 2], "i1")
      i2 = Index([QN(0) => 1, QN(1) => 2], "i2")
      A = randomITensor(i1, i2, dag(i1)', dag(i2)')
      F = eigen(A, (i1, i1'), (i2', i2))
      D, U = F
      Ut = F.Vt
      @test A * U ≈ Ut * D atol = 1e-12
    end
  end

  @testset "svd" for ElT in (Float64, ComplexF64)
    @testset "svd example 1" begin
      i = Index(QN(0) => 2, QN(1) => 2; tags="i")
      j = Index(QN(0) => 2, QN(1) => 2; tags="j")
      A = randomITensor(ElT, QN(0), i, dag(j))
      for b in nzblocks(A)
        @test flux(A, b) == QN(0)
      end
      U, S, V = svd(A, i)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(0)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0)
      end
      @test U * S * V ≈ A atol = 1e-14
    end

    @testset "svd example 2" begin
      i = Index(QN(0) => 5, QN(1) => 6; tags="i")
      j = Index(QN(-1) => 2, QN(0) => 3, QN(1) => 4; tags="j")
      A = randomITensor(ElT, QN(0), i, j)
      for b in nzblocks(A)
        @test flux(A, b) == QN(0)
      end
      U, S, V = svd(A, i)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(0)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0)
      end
      @test U * S * V ≈ A atol = 1e-14
    end

    @testset "svd example 3" begin
      i = Index(QN(0) => 5, QN(1) => 6; tags="i")
      j = Index(QN(-1) => 2, QN(0) => 3, QN(1) => 4; tags="j")
      A = randomITensor(ElT, QN(0), i, dag(j))
      for b in nzblocks(A)
        @test flux(A, b) == QN(0)
      end
      U, S, V = svd(A, i)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(0)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0)
      end
      @test U * S * V ≈ A atol = 1e-14
    end

    @testset "svd example 4" begin
      i = Index(QN(0, 2) => 2, QN(1, 2) => 2; tags="i")
      j = settags(i, "j")

      A = randomITensor(ElT, QN(0, 2), i, j, dag(i'), dag(j'))

      U, S, V = svd(A, i, j)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      for b in nzblocks(A)
        @test flux(A, b) == QN(0, 2)
      end
      U, S, V = svd(A, i)
      for b in nzblocks(U)
        @test flux(U, b) == QN(0, 2)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(0, 2)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0, 2)
      end
      @test U * S * V ≈ A atol = 1e-14
    end

    @testset "svd example 5" begin
      i = Index(QN(0, 2) => 2, QN(1, 2) => 2; tags="i")
      j = settags(i, "j")

      A = randomITensor(ElT, QN(1, 2), i, j, dag(i'), dag(j'))

      U, S, V = svd(A, i, j)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      for b in nzblocks(A)
        @test flux(A, b) == QN(1, 2)
      end
      U, S, V = svd(A, i)
      for b in nzblocks(U)
        @test flux(U, b) == QN(0, 2)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(1, 2)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0, 2)
      end
      @test U * S * V ≈ A atol = 1e-14
    end

    @testset "svd example 6" begin
      i = Index(QN(0, 2) => 2, QN(1, 2) => 2; tags="i")
      j = settags(i, "j")

      A = randomITensor(ElT, QN(1, 2), i, j, dag(i'), dag(j'))

      U, S, V = svd(A, i, i')

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      for b in nzblocks(A)
        @test flux(A, b) == QN(1, 2)
      end
      U, S, V = svd(A, i)
      for b in nzblocks(U)
        @test flux(U, b) == QN(0, 2)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(1, 2)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0, 2)
      end
      @test U * S * V ≈ A atol = 1e-14
    end

    @testset "svd truncation example 1" begin
      i = Index(QN(0) => 2, QN(1) => 3; tags="i")
      j = settags(i, "j")
      A = randomITensor(ElT, QN(0), i, j, dag(i'), dag(j'))
      for i in 1:4
        A = mapprime(A * A', 2, 1)
      end
      A = A / norm(A)

      cutoff = 1e-5
      U, S, V, spec = svd(A, i, j; utags="x", vtags="y", cutoff=cutoff)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      u = commonind(S, U)
      v = commonind(S, V)

      @test hastags(u, "x")
      @test hastags(v, "y")

      @test hassameinds(U, (i, j, u))
      @test hassameinds(V, (i', j', v))

      for b in nzblocks(A)
        @test flux(A, b) == QN(0)
      end
      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(0)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0)
      end

      Ap = U * S * V

      @test norm(Ap - A) ≤ 1e-2
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i) * dim(j)

      @test spec.truncerr ≤ cutoff
      err = real(1 - (Ap * dag(Ap))[] / (A * dag(A))[])
      @test err ≤ cutoff
      @test isapprox(err, spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 2" begin
      i = Index(QN(0) => 3, QN(1) => 2; tags="i")
      j = settags(i, "j")
      A = randomITensor(ElT, QN(0), i, j, dag(i'), dag(j'))

      maxdim = 4
      U, S, V, spec = svd(A, i, j; utags="x", vtags="y", maxdim=maxdim)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      u = commonind(S, U)
      v = commonind(S, V)

      @test hastags(u, "x")
      @test hastags(v, "y")

      @test hassameinds(U, (i, j, u))
      @test hassameinds(V, (i', j', v))

      for b in nzblocks(A)
        @test flux(A, b) == QN(0)
      end
      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(0)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i) * dim(j)

      Ap = U * S * V
      err = 1 - (Ap * dag(Ap))[] / (A * dag(A))[]
      @test isapprox(err, spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 3" begin
      i = Index(QN(0) => 2, QN(1) => 3, QN(2) => 4; tags="i")
      j = settags(i, "j")
      A = randomITensor(ElT, QN(1), i, j, dag(i'), dag(j'))

      maxdim = 4
      U, S, V, spec = svd(A, i, j; utags="x", vtags="y", maxdim=maxdim)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      u = commonind(S, U)
      v = commonind(S, V)

      @test hastags(u, "x")
      @test hastags(v, "y")

      @test hassameinds(U, (i, j, u))
      @test hassameinds(V, (i', j', v))

      for b in nzblocks(A)
        @test flux(A, b) == QN(1)
      end
      for b in nzblocks(U)
        @test flux(U, b) == QN(0)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(1)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i) * dim(j)

      Ap = U * S * V
      err = 1 - (Ap * dag(Ap))[] / (A * dag(A))[]
      @test isapprox(err, spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 4" begin
      i = Index(QN(0, 2) => 3, QN(1, 2) => 4; tags="i")
      j = settags(i, "j")
      A = randomITensor(ElT, QN(1, 2), i, j, dag(i'), dag(j'))

      maxdim = 4
      U, S, V, spec = svd(A, i, j; utags="x", vtags="y", maxdim=maxdim)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      u = commonind(S, U)
      v = commonind(S, V)

      @test hastags(u, "x")
      @test hastags(v, "y")

      @test hassameinds(U, (i, j, u))
      @test hassameinds(V, (i', j', v))

      for b in nzblocks(A)
        @test flux(A, b) == QN(1, 2)
      end
      for b in nzblocks(U)
        @test flux(U, b) == QN(0, 2)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(1, 2)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0, 2)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i) * dim(j)

      Ap = U * S * V
      err = 1 - (Ap * dag(Ap))[] / (A * dag(A))[]
      @test isapprox(err, spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 5" begin
      i = Index(QN(0, 2) => 2, QN(1, 2) => 3; tags="i")
      j = settags(i, "j")
      A = randomITensor(ElT, QN(1, 2), i, j, dag(i'), dag(j'))

      maxdim = 4
      U, S, V, spec = svd(A, i, j'; utags="x", vtags="y", maxdim=maxdim)

      @test storage(U) isa NDTensors.BlockSparse
      @test storage(S) isa NDTensors.DiagBlockSparse
      @test storage(V) isa NDTensors.BlockSparse

      u = commonind(S, U)
      v = commonind(S, V)

      @test hastags(u, "x")
      @test hastags(v, "y")

      @test hassameinds(U, (i, j', u))
      @test hassameinds(V, (i', j, v))

      for b in nzblocks(A)
        @test flux(A, b) == QN(1, 2)
      end
      for b in nzblocks(U)
        @test flux(U, b) == QN(0, 2)
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN(1, 2)
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN(0, 2)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i) * dim(j)

      Ap = U * S * V
      err = 1 - (Ap * dag(Ap))[] / (A * dag(A))[]
      @test isapprox(err, spec.truncerr; rtol=1e-6)
    end

    @testset "issue #231" begin
      l = Index(
        QN("Nf", -1, -1) => 2, QN("Nf", 0, -1) => 4, QN("Nf", +1, -1) => 2; tags="CMB,Link"
      )
      s = Index(QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1; tags="Fermion,Site,n=4")
      r = Index(
        QN("Nf", 1, -1) => 2, QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 2; tags="Link,u"
      )

      A = emptyITensor(ElT, l, s, dag(r))

      insertblock!(A, Block(2, 1, 2))
      insertblock!(A, Block(1, 2, 2))
      insertblock!(A, Block(2, 2, 3))

      for b in nzblocks(A)
        @test flux(A, b) == QN()
      end

      U, S, V = svd(A, l, s)

      for b in nzblocks(U)
        @test flux(U, b) == QN()
      end
      for b in nzblocks(S)
        @test flux(S, b) == QN()
      end
      for b in nzblocks(V)
        @test flux(V, b) == QN()
      end
      @test U * S * V ≈ A atol = 1e-15
    end

    @testset "SVD no truncate bug" begin
      s = Index(
        QN("Sz", -4) => 1,
        QN("Sz", -2) => 4,
        QN("Sz", 0) => 6,
        QN("Sz", 2) => 4,
        QN("Sz", 4) => 1,
      )
      A = emptyITensor(ElT, s, s')
      insertblock!(A, Block(5, 2))
      insertblock!(A, Block(4, 3))
      insertblock!(A, Block(3, 4))
      insertblock!(A, Block(2, 5))
      randn!(A)
      U, S, V = svd(A, s)
      @test U * S * V ≈ A
    end

    @testset "SVD no truncate" begin
      s = Index(
        QN("Sz", -4) => 1,
        QN("Sz", -2) => 4,
        QN("Sz", 0) => 6,
        QN("Sz", 2) => 4,
        QN("Sz", 4) => 1,
      )
      A = emptyITensor(ElT, s, s')
      insertblock!(A, Block(5, 1))
      insertblock!(A, Block(4, 2))
      insertblock!(A, Block(3, 3))
      insertblock!(A, Block(2, 4))
      insertblock!(A, Block(1, 5))
      U, S, V = svd(A, s)
      @test dims(S) == dims(A)
      @test U * S * V ≈ A
    end

    @testset "SVD truncate zeros" begin
      s = Index(
        QN("Sz", -4) => 1,
        QN("Sz", -2) => 4,
        QN("Sz", 0) => 6,
        QN("Sz", 2) => 4,
        QN("Sz", 4) => 1,
      )
      A = emptyITensor(ElT, s, s')
      insertblock!(A, Block(5, 1))
      insertblock!(A, Block(4, 2))
      insertblock!(A, Block(3, 3))
      insertblock!(A, Block(2, 4))
      insertblock!(A, Block(1, 5))
      U, S, V = svd(A, s; cutoff=0)
      @test dims(S) == (0, 0)
      @test U * S * V ≈ A
    end
  end

  @testset "Replace Index" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(0) => 3, QN(1) => 4, QN(2) => 5], "j")

    T1 = randomITensor(QN(1), i, j)
    T2 = copy(T1)

    k = Index([QN(0) => 1, QN(1) => 2], "k")

    replaceind!(T1, i, k)
    @test hasind(T1, k)
    @test dir(inds(T1)[1]) == dir(i)

    # Check that replaceind! keeps
    # original Arrow direction
    replaceind!(T2, i, dag(k))
    @test hasind(T2, k)
    @test dir(inds(T2)[1]) == dir(i)
    @test dir(inds(T2)[1]) != dir(dag(k))
  end

  @testset "BlockSparse dag copy behavior" begin
    i = Index(QN(0) => 2, QN(1) => 2; tags="i")
    j = Index(QN(0) => 2, QN(1) => 2; tags="j")

    v1 = randomITensor(QN(1), i, j)
    orig_elt = v1[1, 3]
    cv1 = dag(v1; allow_alias=true)
    cv1[1, 3] = 123.45
    @test v1[1, 3] ≈ cv1[1, 3]

    v1 = randomITensor(QN(1), i, j)
    orig_elt = v1[1, 3]
    cv1 = dag(ITensors.AllowAlias(), v1)
    cv1[1, 3] = 123.45
    @test v1[1, 3] ≈ cv1[1, 3]

    v2 = randomITensor(QN(1), i, j)
    orig_elt = v2[1, 3]
    cv2 = dag(v2; allow_alias=false)
    cv2[1, 3] = 123.45
    @test v2[1, 3] ≈ orig_elt

    v2 = randomITensor(QN(1), i, j)
    orig_elt = v2[1, 3]
    cv2 = dag(ITensors.NeverAlias(), v2)
    cv2[1, 3] = 123.45
    @test v2[1, 3] ≈ orig_elt

    v3 = randomITensor(ComplexF64, QN(1), i, j)
    orig_elt = v3[1, 3]
    cv3 = dag(v3; allow_alias=true)
    cv3[1, 3] = 123.45
    @test v3[1, 3] ≈ orig_elt

    v3 = randomITensor(ComplexF64, QN(1), i, j)
    orig_elt = v3[1, 3]
    cv3 = dag(ITensors.AllowAlias(), v3)
    cv3[1, 3] = 123.45
    @test v3[1, 3] ≈ orig_elt

    v4 = randomITensor(ComplexF64, QN(1), i, j)
    orig_elt = v4[1, 3]
    cv4 = dag(v4; allow_alias=false)
    cv4[1, 3] = 123.45
    @test v4[1, 3] ≈ orig_elt

    v4 = randomITensor(ComplexF64, QN(1), i, j)
    orig_elt = v4[1, 3]
    cv4 = dag(ITensors.NeverAlias(), v4)
    cv4[1, 3] = 123.45
    @test v4[1, 3] ≈ orig_elt
  end

  @testset "exponentiate" begin
    @testset "Simple arrows" begin
      i1 = Index([QN(0) => 1, QN(1) => 2], "i1")
      i2 = Index([QN(0) => 1, QN(1) => 2], "i2")
      A = randomITensor(QN(), i1, i2, dag(i1)', dag(i2)')
      Aexp = exp(A)
      Amat = Array(A, i1, i2, i1', i2')
      Amatexp = reshape(exp(reshape(Amat, 9, 9)), 3, 3, 3, 3)
      @test Array(Aexp, i1, i2, i1', i2') ≈ Amatexp rtol = 1e-14
      @test flux(Aexp) == QN()
      @test length(setdiff(inds(Aexp), inds(A))) == 0

      @test exp(A, (i1, i2), (i1', i2')) ≈ Aexp rtol = 5e-14

      # test the case where indices are permuted
      A = randomITensor(QN(), i1, dag(i1)', dag(i2)', i2)
      Aexp = exp(A, (i1, i2), (i1', i2'))
      Amat = Array(A, i1, i2, i1', i2')
      Amatexp = reshape(exp(reshape(Amat, 9, 9)), 3, 3, 3, 3)
      @test Array(Aexp, i1, i2, i1', i2') ≈ Amatexp rtol = 1e-14

      # test exponentiation in the Hermitian case
      i1 = Index([QN(0) => 2, QN(1) => 2, QN(2) => 3], "i1")
      A = randomITensor(QN(), i1, dag(i1)')
      Ad = dag(swapinds(A, IndexSet(i1), IndexSet(dag(i1)')))
      Ah = A + Ad + 1e-10 * randomITensor(QN(), i1, dag(i1)')
      Amat = Array(Ah, i1', i1)
      Aexp = exp(Ah; ishermitian=true)
      Amatexp = exp(LinearAlgebra.Hermitian(Amat))
      @test Array(Aexp, i1, i1') ≈ Amatexp rtol = 5e-14
    end

    @testset "Regression test for exp of QN ITensor with missing diagonal blocks" begin
      i = Index([QN(0) => 2, QN(1) => 3])
      A = ITensor(i', dag(i))
      A[1, 1] = 1.2
      expA = exp(A; ishermitian=false)
      for n in 1:mindim(A)
        @test expA[n, n] == exp(A[n, n])
      end
      @test expA ≈ exp(dense(A))
      expA = exp(A; ishermitian=true)
      for n in 1:mindim(A)
        @test expA[n, n] == exp(A[n, n])
      end
      @test expA ≈ exp(dense(A))
    end

    @testset "Mixed arrows" begin
      i1 = Index([QN(0) => 1, QN(1) => 2], "i1")
      i2 = Index([QN(0) => 1, QN(1) => 2], "i2")
      A = randomITensor(i1, i2, dag(i1)', dag(i2)')
      expA = exp(A, (i1, i1'), (i2', i2))
      @test exp(dense(A), (i1, i1'), (i2', i2)) ≈ dense(expA)
    end

    @testset "Test contraction direction error" begin
      i = Index([QN(0) => 1, QN(1) => 1], "i")
      A = randomITensor(i', dag(i))
      A² = A' * A
      @test dense(A²) ≈ dense(A') * dense(A)
      @test_throws ErrorException A' * dag(A)
    end

    @testset "Contraction resulting in no blocks with threading bug" begin
      i = Index([QN(0) => 1, QN(1) => 1])
      A = emptyITensor(i', dag(i))
      B = emptyITensor(i', dag(i))
      A[i' => 1, i => 1] = 11.0
      B[i' => 2, i => 2] = 22.0

      using_threaded_blocksparse = ITensors.disable_threaded_blocksparse()
      C1 = A' * B
      ITensors.enable_threaded_blocksparse()
      C2 = A' * B
      if using_threaded_blocksparse
        ITensors.enable_threaded_blocksparse()
      else
        ITensors.disable_threaded_blocksparse()
      end

      @test nnzblocks(C1) == 0
      @test nnzblocks(C2) == 0
      @test nnz(C1) == 0
      @test nnz(C2) == 0
      @test C1 ≈ C2
    end

    @testset "Contraction with scalar ITensor" begin
      i = Index([QN(0) => 2, QN(1) => 2])
      A = randomITensor(i', dag(i))
      A1 = A * ITensor(1)
      A2 = ITensor(1) * A
      @test A1 ≈ A
      @test A2 ≈ A
    end
  end

  @testset "directsum" begin
    x = Index([QN(0) => 1, QN(1) => 1], "x")
    i1 = Index([QN(0) => 1, QN(1) => 2], "i1")
    j1 = Index([QN(0) => 2, QN(1) => 2], "j1")
    i2 = Index([QN(0) => 2, QN(1) => 3], "i2")
    j2 = Index([QN(0) => 3, QN(1) => 3], "j2")

    A1 = randomITensor(i1, x, j1)
    A2 = randomITensor(x, j2, i2)
    S, s = ITensors.directsum(A1 => (i1, j1), A2 => (i2, j2); tags=["sum_i", "sum_j"])

    @test hassameinds(S, (x, s...))
    @test hastags(s[1], "sum_i")
    @test hastags(s[2], "sum_j")

    for vx in 1:dim(x)
      proj = dag(onehot(x => vx))
      A1_vx = A1 * proj
      A2_vx = A2 * proj
      S_vx = S * proj
      for m in 1:dim(s[1]), n in 1:dim(s[2])
        if m ≤ dim(i1) && n ≤ dim(j1)
          @test S_vx[s[1] => m, s[2] => n] == A1_vx[i1 => m, j1 => n]
        elseif m > dim(i1) && n > dim(j1)
          @test S_vx[s[1] => m, s[2] => n] == A2_vx[i2 => m - dim(i1), j2 => n - dim(j1)]
        else
          @test S_vx[s[1] => m, s[2] => n] == 0
        end
      end
    end
  end

  @testset "Negate QN ITensor Regression Test" begin
    s = siteind("S=1/2"; conserve_qns=true)

    A = ITensor(s', dag(s))
    A[1, 1] = 1.0

    @test length(ITensors.blockoffsets(ITensors.tensor(A))) == 1
    B = -A # there was a bug where doing -A would
    # increase the number of blocks of A's storage
    @test length(ITensors.blockoffsets(ITensors.tensor(A))) == 1
  end
end

nothing


text:
module ControlExamplePlots

using ControlSystems, Plots, VisualRegressionTests, Gtk

include("genplots.jl")

end # module


text:
module Rafts

# A Raft "floats" over data within directories for interactive data exploration

Raft() = Raft(pwd())
function Raft(path::String)
    isdir(path) && return DirectoryRaft(path)
    files = readdir(path)
    # remove dot-files
    filter!(f->!beginswith(f, "."), files)

    splat_mask = [beginswith(f, "_") for f in files]

    children = Dict{Any, Raft}()
    sizehint(children, sum(splat_mask))
    for f in files[!splat_mask]
        children[f] = Raft(joinpath(path, f), false)
    end
    if all(f->isdir(f) && files)
    end
end

end # module


text:
module GeometricalProductSpecification

include("gauss_filters.jl")
include("parameters_1d.jl")

import Plots: plot

export SurfaceProfile

struct SurfaceProfile
    data::Array{Float64,1}
    resolution::Float64
    filterwavelengths::Vector{Float64}
    components::Array{Array{Float64,1},1}
    features::Array{Function,1}
    parameters::Dict{String, Float64}
    function SurfaceProfile(data, res, filterwavelengths=[2.5, 250.0, 2500.0],
                            features=[], parameters=Dict())
        components = idct.(separatefeatures(data, 1/res, filterwavelengths))
        return new(data, res, filterwavelengths, components, features, parameters)
    end
end

function plot(profile::SurfaceProfile)
    p1 = plot(profile.R, label="R")
        plot!(profile.W, label="W")
        plot!(profile.F, label="F")
    return p1
end

end #module GeometricalProductSpecifications