text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function iterate!(ph::AbstractProgressiveHedging, ::AbstractProgressiveHedgingExecution)
    # Resolve all subproblems at the current optimal solution
    Q = resolve_subproblems!(ph)
    if !(Q.status ∈ AcceptableTermination)
        # Early termination log
        log!(ph; status = Q.status)
        return Q.status
    end
    ph.data.Q = Q.value
    # Update iterate
    update_iterate!(ph)
    # Update subproblems
    update_subproblems!(ph)
    # Get dual gap
    update_dual_gap!(ph)
    # Update penalty (if applicable)
    update_penalty!(ph)
    # Check optimality
    if check_optimality(ph)
        # Final log
        log!(ph; optimal = true)
        # Optimal
        return MOI.OPTIMAL
    end
    # Calculate time spent so far and check perform time limit check
    time_spent = ph.progress.tlast - ph.progress.tinit
    if time_spent >= ph.parameters.time_limit
        log!(ph; status = MOI.TIME_LIMIT)
        return MOI.TIME_LIMIT
    end
    # Log progress
    log!(ph)
    # Dont return a status as procedure should continue
    return nothing
end

function finish_initilization!(execution::AbstractProgressiveHedgingExecution, penalty::AbstractFloat)
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                execution.subworkers[w-1],
                penalty) do sw, penalty
                    for subproblem in fetch(sw)
                        initialize!(subproblem, penalty)
                    end
                end
        end
    end
    return nothing
end

function start_workers!(::AbstractProgressiveHedging, ::AbstractProgressiveHedgingExecution)
    return nothing
end

function close_workers!(::AbstractProgressiveHedging, ::AbstractProgressiveHedgingExecution)
    return nothing
end


text:
module TestBioBridgeR

using Base.Test
using BioBridgeR.APE
using BioSymbols
using BioSequences
using RCall

include("ape.jl")

end # module TestBioBridgeR

ext: jl


text:
export getchar, set_homedir, get_homedir

charset = 'a':'z'

getchar(idx) = charset[idx]

homedir = "wrong"

function set_homedir(path)
    global homedir = path
end
get_homedir() = homedir


text:
"""
    ito_correction!(I, h=1)

Applies the Itô-correction for iterated integrals to `I`.
This amounts to subtracting ``\\frac{1}{2}h`` from every element on the diagonal.

# Example
```jldoctest; setup=:(using LevyArea)
julia> M = ones(5,5);

julia> LevyArea.ito_correction!(M, 0.5)

julia> M
5×5 Matrix{Float64}:
 0.75  1.0   1.0   1.0   1.0
 1.0   0.75  1.0   1.0   1.0
 1.0   1.0   0.75  1.0   1.0
 1.0   1.0   1.0   0.75  1.0
 1.0   1.0   1.0   1.0   0.75
```

"""
function ito_correction!(I, h=1)
    m,n = size(I)
    m == n || throw(DimensionMismatch("Matrix is not square: dimensions are $(size(I))"))
    @inbounds for i=1:m
        I[i,i] -= h/2
    end
end

"""
    iterated_integrals(W::AbstractVector, h, eps=h^(3/2);
        ito_correction=true,
        error_norm=MaxL2(),
        alg=optimal_algorithm(length(W),h,eps,error_norm)
    )

Simulates an approximation of the iterated stochastic integrals
``\\int_0^h\\int_0^sdW_i(t)dW_j(s)`` for all pairs ``1\\le i, j \\le m``
of the given ``m``-dimensional Brownian motion with step size h.

# Examples
```jldoctest; setup=:(using LinearAlgebra; using LevyArea)
julia> h = 1/2;

julia> W = [1.0, 0.5]
2-element Vector{Float64}:
 1.0
 0.5

julia> diag(iterated_integrals(W, h, h^(3/2))) ≈ 0.5*W.^2 .- 0.5h
true
```
"""
function iterated_integrals(W::AbstractVector, h::Real, eps::Real=h^(3/2);
    ito_correction=true,
    error_norm::AbstractErrorNorm=MaxL2(),
    alg::AbstractIteratedIntegralAlgorithm=optimal_algorithm(length(W),h,eps,error_norm)
)
    m = length(W)
    n = terms_needed(m, h, eps, alg, error_norm)
    I = levyarea(W/√h, n, alg)
    if ito_correction
        ito_correction!(I)
    end
    I .= 1//2 .* W.*W' .+ h.*I
    return I
end

"""
    iterated_integrals(W::AbstractVector, q_12::AbstractVector, h, eps;
        ito_correction=true,
        error_norm=FrobeniusL2(),
        alg=optimal_algorithm(length(W),q_12,h,eps,error_norm)
    )

Simulates an approximation of the iterated stochastic integrals for finite-dimensional approximations of
a Q-Wiener process with covariance matrix ``Q = Q^\\frac{1}{2}*Q^\\frac{1}{2}``.
Here `q_12` is a vector of the eigenvalues of ``Q^\\frac{1}{2}``; the square root of the covariance matrix.
Equivalently these are the square roots of the eigenvalues of ``Q``.

# Examples
```jldoctest; setup=:(using LinearAlgebra; using LevyArea)
julia> h = 0.01; dim=10; q = [1/k^2 for k=1:dim];

julia> W = √h * sqrt.(q) .* randn(dim);

julia> diag(iterated_integrals(W,sqrt.(q),h,h^(3/2))) ≈ 0.5*W.^2 .- 0.5*h*q
true
```
"""
function iterated_integrals(W::AbstractVector, q_12::AbstractVector, h::Real, eps::Real;
    ito_correction=true,
    error_norm::AbstractErrorNorm=FrobeniusL2(),
    alg::AbstractIteratedIntegralAlgorithm=optimal_algorithm(length(W),q_12,h,eps,error_norm)
)
    m = length(W)
    n = terms_needed(m, q_12, h, eps, alg, error_norm)
    I = levyarea(W./q_12./√h, n, alg)
    if ito_correction
        ito_correction!(I)
    end
    I .= 1//2 .* W.*W' .+ h.*q_12'.*I.*q_12 # scale correctly
    return I
end

"""
    iterated_integrals(W::Real, h::Real, eps::Real=0.0; ito_correction=true, kwargs...)

In the case of a scalar Brownian motion the integral can be explicitly
calculated as ``\\int_0^h\\int_0^sdW(t)dW(s) = \\frac{1}{2}W(h)^2 - \\frac{1}{2}h``.

The parameter `eps` (as well as all additional keyword arguments) has no effect but is available
to provide the same interface as the multidimensional version.
"""
iterated_integrals(W::Real, h::Real, eps::Real=0.0; ito_correction=true, kwargs...) = ito_correction ? (W^2-h)/2 : W^2/2

"""
    iterated_integrals(W::Real, q_12::Real, h::Real, eps::Real; ito_correction=true, kwargs...)

In the case of a scalar Q-Wiener process with (scalar) covariance Q the integral can be explicitly
calculated as ``\\int_0^h\\int_0^sdW(t)dW(s) = \\frac{1}{2}W(h)^2 - \\frac{1}{2}hQ``.

Note that, as in the multidimensional case, the parameter `q_12` denotes the square root of the covariance.

The parameter `eps` (as well as all additional keyword arguments) has no effect but is available
to provide the same interface as the multidimensional version.
"""
iterated_integrals(W::Real, q_12::Real, h::Real, eps::Real; ito_correction=true, kwargs...) = ito_correction ? (W^2-h*q_12^2)/2 : W^2/2

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames: DataFrame
using CUDA: CUDA

@testset "MLE tests" begin

    # Binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    estimate_parameters!(r,dfb; pseudocount=1.0)
    @test log_likelihood_avg(r,dfb) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)

    estimate_parameters!(r,dfb; pseudocount=0.0)
    @test log_likelihood_avg(r,dfb) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=0.0)

    if CUDA.functional()
        # Binary dataset
        dfb_gpu = to_gpu(dfb)

        estimate_parameters!(r, dfb_gpu; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)

        estimate_parameters!(r, dfb_gpu; pseudocount=0.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=0.0)
    end

end

@testset "Weighted MLE tests" begin
    # Binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    # Weighted binary dataset
    weights = [0.6, 0.6, 0.6]
    wdfb = weigh_samples(dfb, weights)

    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)

    estimate_parameters!(r,wdfb; pseudocount=1.0)
    @test log_likelihood_avg(r,dfb) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)
    @test log_likelihood_avg(r,dfb) ≈ log_likelihood_avg(r,wdfb)

    estimate_parameters!(r,wdfb; pseudocount=0.0)
    @test log_likelihood_avg(r,dfb) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=0.0)
    @test log_likelihood_avg(r,dfb) ≈ log_likelihood_avg(r,wdfb)

    if CUDA.functional()

        # Binary dataset
        dfb_gpu = to_gpu(dfb)

        # Weighted binary dataset
        wdfb_gpu = to_gpu(wdfb)

        estimate_parameters!(r, wdfb_gpu; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ log_likelihood_avg(r, wdfb_gpu)

        estimate_parameters!(r, dfb_gpu; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ log_likelihood_avg(r, wdfb_gpu)

        estimate_parameters!(r, wdfb_gpu; pseudocount=0.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=0.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ log_likelihood_avg(r, wdfb_gpu)

    end
end

@testset "Batched MLE tests" begin
    # Binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    # Weighted binary dataset
    weights = [0.6, 0.6, 0.6]
    wdfb = weigh_samples(dfb, weights)

    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)

    # Batched dataset
    batched_wdfb = batch(wdfb, 1)

    estimate_parameters!(r,batched_wdfb; pseudocount=1.0)
    @test log_likelihood_avg(r,dfb) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)
    @test log_likelihood_avg(r,dfb) ≈ log_likelihood_avg(r,wdfb)
    @test log_likelihood_avg(r,batched_wdfb) ≈ log_likelihood_avg(r,wdfb)

    estimate_parameters!(r,batched_wdfb; pseudocount=0.0)
    @test log_likelihood_avg(r,dfb) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=0.0)
    @test log_likelihood_avg(r,dfb) ≈ log_likelihood_avg(r,wdfb)
    @test log_likelihood_avg(r,batched_wdfb) ≈ log_likelihood_avg(r,wdfb)

    if CUDA.functional()

        # Binary dataset
        dfb_gpu = to_gpu(dfb)

        # Weighted binary dataset
        wdfb_gpu = to_gpu(wdfb)

        # Batched dataset
        batched_wdfb_gpu = to_gpu(batch(wdfb, 1))

        estimate_parameters!(r, batched_wdfb_gpu; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ log_likelihood_avg(r, wdfb_gpu)

        estimate_parameters!(r, batched_wdfb_gpu; pseudocount=0.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=0.0)
        @test log_likelihood_avg(r,dfb_gpu) ≈ log_likelihood_avg(r, wdfb_gpu)

    end
end

@testset "Soft MLE test" begin
    # Batched binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false false]), :auto)
    dfb = soften(dfb, 0.001; scale_by_marginal = false)
    batched_dfb = batch(dfb, 1)

    # Batched weighted binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false false]), :auto)
    dfb = soften(dfb, 0.001; scale_by_marginal = false)
    weights = [0.6, 0.6, 0.6]
    wdfb = weigh_samples(dfb, weights)
    batched_wdfb = batch(wdfb, 1)

    # Weighted binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false false]), :auto)
    dfb = soften(dfb, 0.001; scale_by_marginal = false)
    weights = [0.6, 0.6, 0.6]
    wdfb = weigh_samples(dfb, weights)

    # Binary dataset
    dfb = DataFrame(BitMatrix([true false; true true; false false]), :auto)
    dfb = soften(dfb, 0.001; scale_by_marginal = false)

    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    estimate_parameters!(r, dfb; pseudocount=0.0)
    @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
    @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

    estimate_parameters!(r, wdfb; pseudocount=0.0)
    @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
    @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

    estimate_parameters!(r, batched_dfb; pseudocount=0.0)
    @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
    @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

    estimate_parameters!(r, batched_wdfb; pseudocount=0.0)
    @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
    @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
    @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

    if CUDA.functional()
        estimate_parameters!(r, to_gpu(dfb); pseudocount=0.0)
        @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
        @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

        estimate_parameters!(r, to_gpu(wdfb); pseudocount=0.0)
        @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
        @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

        estimate_parameters!(r, to_gpu(batched_dfb); pseudocount=0.0)
        @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
        @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6

        estimate_parameters!(r, to_gpu(batched_wdfb); pseudocount=0.0)
        @test r.children[1].children[1].log_probs[1] ≈ -0.4059652245524093 atol = 1e-6
        @test r.children[1].children[1].log_probs[2] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[1] ≈ -1.0976127839931185 atol = 1e-6
        @test r.children[1].children[2].log_probs[2] ≈ -0.4059652245524093 atol = 1e-6
    end

end

@testset "EM tests" begin
    data = DataFrame([true missing], :auto)
    vtree2 = PlainVtree(2, :balanced)
    pc = fully_factorized_circuit(StructProbCircuit, vtree2).children[1]
    uniform_parameters!(pc)
    pc.children[1].prime.log_probs .= log.([0.3, 0.7])
    pc.children[1].sub.log_probs .= log.([0.4, 0.6])
    estimate_parameters_em!(pc, data; pseudocount=0.0)
    @test all(pc.children[1].prime.log_probs .== log.([1.0, 0.0]))
    @test pc.children[1].sub.log_probs[1] .≈ log.([0.4, 0.6])[1] atol=1e-6

    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
    uniform_parameters!(r)
    estimate_parameters!(r,dfb; pseudocount=1.0)
    paras1 = ParamBitCircuit(r, dfb).params
    uniform_parameters!(r)
    estimate_parameters_em!(r, dfb; pseudocount=1.0)
    paras2 = ParamBitCircuit(r, dfb).params
    @test all(paras1 .≈ paras2)

    if CUDA.functional()
        dfb_gpu = to_gpu(dfb)

        uniform_parameters!(r)
        estimate_parameters_em!(r, dfb_gpu; pseudocount=1.0)
        paras1 = ParamBitCircuit(r, dfb).params

        uniform_parameters!(r)
        estimate_parameters!(r, dfb_gpu; pseudocount=1.0)
        paras2 = ParamBitCircuit(r, dfb).params

        @test all(paras1 .≈ paras2)

        data = DataFrame([true missing; false missing], :auto)
        r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

        uniform_parameters!(r)
        estimate_parameters_em!(r, data; pseudocount = 0.0)
        paras1 = ParamBitCircuit(r, dfb).params

        uniform_parameters!(r)
        estimate_parameters_em!(r, to_gpu(data); pseudocount = 0.0)
        paras2 = ParamBitCircuit(r, dfb).params

        @test all(paras1 .≈ paras2)
    end
end

@testset "Weighted EM tests" begin
    data = DataFrame([true missing], :auto)
    weights = DataFrame(weight = [1.0])
    wdata = hcat(data, weights)

    vtree2 = PlainVtree(2, :balanced)
    pc = fully_factorized_circuit(StructProbCircuit, vtree2).children[1]
    uniform_parameters!(pc)
    pc.children[1].prime.log_probs .= log.([0.3, 0.7])
    pc.children[1].sub.log_probs .= log.([0.4, 0.6])
    estimate_parameters_em!(pc, wdata; pseudocount=0.0)
    @test all(pc.children[1].prime.log_probs .== log.([1.0, 0.0]))
    @test pc.children[1].sub.log_probs[1] .≈ log.([0.4, 0.6])[1] atol=1e-6

    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    weights = DataFrame(weight = [0.6, 0.6, 0.6])
    wdfb = hcat(dfb, weights)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
    uniform_parameters!(r)
    estimate_parameters!(r,wdfb; pseudocount=1.0)
    paras1 = ParamBitCircuit(r, wdfb).params
    uniform_parameters!(r)
    estimate_parameters_em!(r, wdfb; pseudocount=1.0)
    paras2 = ParamBitCircuit(r, wdfb).params
    @test all(paras1 .≈ paras2)

    if CUDA.functional()
        uniform_parameters!(r)
        estimate_parameters_em!(r, to_gpu(wdfb); pseudocount=0.0)
        paras1 = ParamBitCircuit(r, wdfb).params

        uniform_parameters!(r)
        estimate_parameters!(r, to_gpu(wdfb); pseudocount=0.0)
        paras2 = ParamBitCircuit(r, wdfb).params

        @test all(paras1 .≈ paras2)

        data = DataFrame([true missing; false missing], :auto)
        weights = DataFrame(weight = [0.6, 0.6])
        wdata = hcat(data, weights)
        r = fully_factorized_circuit(ProbCircuit, num_features(wdata))

        uniform_parameters!(r)
        estimate_parameters_em!(r, wdata; pseudocount = 0.0)
        paras1 = ParamBitCircuit(r, dfb).params

        uniform_parameters!(r)
        estimate_parameters_em!(r, to_gpu(wdata); pseudocount = 0.0)
        paras2 = ParamBitCircuit(r, dfb).params

        @test all(paras1 .≈ paras2)
    end
end

@testset "Batched EM tests" begin
    dfb = DataFrame(BitMatrix([true false; true true; false true; true true]), :auto)
    batched_dfb = batch(dfb, 1)

    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    uniform_parameters!(r)
    estimate_parameters_em!(r, dfb; pseudocount=1.0)
    paras1 = ParamBitCircuit(r, dfb).params
    uniform_parameters!(r)
    estimate_parameters_em!(r, batched_dfb; pseudocount=1.0)
    paras2 = ParamBitCircuit(r, dfb).params
    @test all(paras1 .≈ paras2)

    dfb = DataFrame(BitMatrix([true false; true true; false true; true true]), :auto)
    weights = [0.6, 0.6, 0.6, 0.6]
    wdfb = weigh_samples(dfb, weights)
    batched_wdfb = batch(wdfb, 1)

    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
    uniform_parameters!(r)

    estimate_parameters!(r, batched_wdfb; pseudocount=1.0)
    paras1 = ParamBitCircuit(r, wdfb).params
    uniform_parameters!(r)
    estimate_parameters_em!(r, batched_wdfb; pseudocount=1.0)
    paras2 = ParamBitCircuit(r, wdfb).params
    @test all(paras1 .≈ paras2)

    if CUDA.functional()
        batched_dfb_gpu = to_gpu(batched_dfb)

        uniform_parameters!(r)
        estimate_parameters_em!(r, batched_dfb_gpu; pseudocount=1.0)
        paras1 = ParamBitCircuit(r, dfb).params

        uniform_parameters!(r)
        estimate_parameters!(r, batched_dfb_gpu; pseudocount=1.0)
        paras2 = ParamBitCircuit(r, dfb).params

        @test all(paras1 .≈ paras2)

        batched_wdfb_gpu = to_gpu(batched_wdfb)

        uniform_parameters!(r)
        estimate_parameters_em!(r, batched_wdfb_gpu; pseudocount=1.0)
        paras1 = ParamBitCircuit(r, dfb).params

        uniform_parameters!(r)
        estimate_parameters!(r, batched_wdfb_gpu; pseudocount=1.0)
        paras2 = ParamBitCircuit(r, dfb).params

        @test all(paras1 .≈ paras2)

        data = DataFrame([true missing; false missing], :auto)
        r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
        batched_data = batch(data, 1)

        uniform_parameters!(r)
        estimate_parameters_em!(r, batched_data; pseudocount = 0.0)
        paras1 = ParamBitCircuit(r, dfb).params

        uniform_parameters!(r)
        estimate_parameters_em!(r, to_gpu(batched_data); pseudocount = 0.0)
        paras2 = ParamBitCircuit(r, dfb).params

        @test all(paras1 .≈ paras2)
    end
end

@testset "EM update per batch tests" begin
    dfb = DataFrame(BitMatrix([true false; true true; false true; true true]), :auto)
    batched_dfb = batch(dfb, 1)

    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    uniform_parameters!(r)
    params = estimate_parameters_em!(r, batched_dfb; pseudocount=1.0, update_per_batch = true)
    @test params[5] ≈ 0.0

    if CUDA.functional()

        uniform_parameters!(r)
        params = estimate_parameters_em!(r, to_gpu(batched_dfb); pseudocount=1.0, update_per_batch = true)
        @test to_cpu(params)[5] ≈ 0.0

    end
end

@testset "Bagging tests" begin
    # Binary dataset
    dfb = DataFrame(BitMatrix([true true; true true; true true; true true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
    # bag_dfb = bagging_dataset(dfb; num_bags = 2, frac_examples = 1.0)
    bag_dfb = Vector{DataFrame}(undef, 2)
    bag_dfb[1] = dfb[[2, 1, 3, 4], :]
    bag_dfb[2] = dfb[[4, 3, 2, 1], :]

    r = compile(SharedProbCircuit, r, 2)

    params = estimate_parameters!(r, bag_dfb; pseudocount = 1.0)
    @test all(abs.(params[1] .- params[2]) .< 1e-6)

    params = estimate_parameters!(r, bag_dfb; pseudocount = 2.0)
    @test all(abs.(params[1] .- params[2]) .< 1e-6)
end

@testset "Uniform parameters tests" begin
    r = fully_factorized_circuit(ProbCircuit, 2)

    uniform_parameters!(r)
    @test all(.≈(exp.(r.children[1].children[1].log_probs), [0.5, 0.5], atol = 1e-6))
    @test all(.≈(exp.(r.children[1].children[2].log_probs), [0.5, 0.5], atol = 1e-6))

    uniform_parameters!(r; perturbation = 0.1)
    @test all(exp.(r.children[1].children[1].log_probs) .> 0.9 / (0.9 + 1.1))
    @test all(exp.(r.children[1].children[1].log_probs) .< 1.1 / (0.9 + 1.1))
    @test all(exp.(r.children[1].children[2].log_probs) .> 0.9 / (0.9 + 1.1))
    @test all(exp.(r.children[1].children[2].log_probs) .< 1.1 / (0.9 + 1.1))
end

@testset "Entropy regularization tests" begin
    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    params = estimate_parameters!(r,dfb; pseudocount=1e-6, entropy_reg = 0.1);
    @test exp.(params) ≈ [0.6523896320397116, 0.34761036796028844, 0.6523896320397116, 0.34761036796028844, 1.0] atol = 1e-4

    params = estimate_parameters!(r,dfb; pseudocount=1e-6, entropy_reg = 1.0);
    @test exp.(params) ≈ [0.5841319507970275, 0.4158680492029725, 0.5841319507970275, 0.4158680492029725, 1.0] atol = 1e-4

    params = estimate_parameters!(r,dfb; pseudocount=1e-6, entropy_reg = 5.0);
    @test exp.(params) ≈ [0.5278256829525434, 0.4721743170474565, 0.5278256829525434, 0.4721743170474565, 1.0] atol = 1e-4

    params = estimate_parameters!(r,dfb; pseudocount=1e-6, entropy_reg = 10.0);
    @test exp.(params) ≈ [0.5151599626276689, 0.484840037372331, 0.5151599626276689, 0.484840037372331, 1.0] atol = 1e-4

    if CUDA.functional()
        params = estimate_parameters!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 0.1);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.6523896320397116, 0.34761036796028844, 0.6523896320397116, 0.34761036796028844, 1.0] atol = 1e-4

        params = estimate_parameters!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 1.0);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.5841319659187809, 0.4158680380658034, 0.5841319659187809, 0.4158680380658034, 1.0] atol = 1e-4

        params = estimate_parameters!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 5.0);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.5278256829525434, 0.4721743170474565, 0.5278256829525434, 0.4721743170474565, 1.0] atol = 1e-4

        params = estimate_parameters!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 10.0);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.5151599626276689, 0.484840037372331, 0.5151599626276689, 0.484840037372331, 1.0] atol = 1e-4
    end
end

@testset "EM Entropy regularization tests" begin
    dfb = DataFrame(BitMatrix([true false; true true; false true]), :auto)
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))

    params = estimate_parameters_em!(r,dfb; pseudocount=1e-6, entropy_reg = 0.1);
    @test exp.(params) ≈ [0.6523896320397116, 0.34761036796028844, 0.6523896320397116, 0.34761036796028844, 1.0] atol = 1e-4

    params = estimate_parameters_em!(r,dfb; pseudocount=1e-6, entropy_reg = 1.0);
    @test exp.(params) ≈ [0.5841319507970275, 0.4158680492029725, 0.5841319507970275, 0.4158680492029725, 1.0] atol = 1e-4

    params = estimate_parameters_em!(r,dfb; pseudocount=1e-6, entropy_reg = 5.0);
    @test exp.(params) ≈ [0.5278256829525434, 0.4721743170474565, 0.5278256829525434, 0.4721743170474565, 1.0] atol = 1e-4

    params = estimate_parameters_em!(r,dfb; pseudocount=1e-6, entropy_reg = 10.0);
    @test exp.(params) ≈ [0.5151599626276689, 0.484840037372331, 0.5151599626276689, 0.484840037372331, 1.0] atol = 1e-4

    if CUDA.functional()
        params = estimate_parameters_em!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 0.1);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.6523896320397116, 0.34761036796028844, 0.6523896320397116, 0.34761036796028844, 1.0] atol = 1e-4

        params = estimate_parameters_em!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 1.0);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.5841319659187809, 0.4158680380658034, 0.5841319659187809, 0.4158680380658034, 1.0] atol = 1e-4

        params = estimate_parameters_em!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 5.0);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.5278256829525434, 0.4721743170474565, 0.5278256829525434, 0.4721743170474565, 1.0] atol = 1e-4

        params = estimate_parameters_em!(r, to_gpu(dfb); pseudocount=1e-6, entropy_reg = 10.0);
        params = convert(Vector{Float64}, params)
        @test exp.(params) ≈ [0.5151599626276689, 0.484840037372331, 0.5151599626276689, 0.484840037372331, 1.0] atol = 1e-4
    end
end


text:
using GalacticOptim, Optim, Test, Random

rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p  = [1.0, 100.0]

f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
l1 = rosenbrock(x0, _p)
prob = OptimizationProblem(f, x0, _p)
sol = solve(prob, SimulatedAnnealing())
@test 10*sol.minimum < l1

Random.seed!(1234)
prob = OptimizationProblem(f, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, SAMIN())
@test 10*sol.minimum < l1

using CMAEvolutionStrategy
sol = solve(prob, CMAEvolutionStrategyOpt())
@test 10*sol.minimum < l1

rosenbrock(x, p=nothing) =  (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

l1 = rosenbrock(x0)
prob = OptimizationProblem(rosenbrock, x0)
sol = solve(prob, NelderMead())
@test 10*sol.minimum < l1

cons= (x,p) -> [x[1]^2 + x[2]^2]
optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= cons)

prob = OptimizationProblem(optprob, x0)

sol = solve(prob, ADAM(0.1), maxiters = 1000)
@test 10*sol.minimum < l1

sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

sol = solve(prob, Newton())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.KrylovTrustRegion())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lcons = [-Inf], ucons = [Inf])
sol = solve(prob, IPNewton())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lcons = [-5.0], ucons = [10.0])
sol = solve(prob, IPNewton())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lcons = [-Inf], ucons = [Inf], lb = [-500.0,-500.0], ub=[50.0,50.0])
sol = solve(prob, IPNewton())
@test sol.minimum < l1

function con2_c(x,p)
    [x[1]^2 + x[2]^2, x[2]*sin(x[1])-x[1]]
end

optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= con2_c)
prob = OptimizationProblem(optprob, x0, lcons = [-Inf,-Inf], ucons = [Inf,Inf])
sol = solve(prob, IPNewton())
@test 10*sol.minimum < l1

cons_circ = (x,p) -> [x[1]^2 + x[2]^2]
optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= cons_circ)
prob = OptimizationProblem(optprob, x0, lcons = [-Inf], ucons = [0.25^2])
sol = solve(prob, IPNewton())
@test sqrt(cons(sol.minimizer,nothing)[1]) ≈ 0.25 rtol = 1e-6

optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())
prob = OptimizationProblem(optprob, x0)
sol = solve(prob, ADAM(), maxiters = 1000, progress = false)
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, Fminbox())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
@test_broken @test_nowarn sol = solve(prob, SAMIN())
@test 10*sol.minimum < l1

using NLopt
prob = OptimizationProblem(optprob, x0)
sol = solve(prob, Opt(:LN_BOBYQA, 2))
@test 10*sol.minimum < l1

sol = solve(prob, Opt(:LD_LBFGS, 2))
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, Opt(:LD_LBFGS, 2))
@test 10*sol.minimum < l1

sol = solve(prob, Opt(:G_MLSL_LDS, 2), nstart=2, local_method = Opt(:LD_LBFGS, 2), maxiters=10000)
@test 10*sol.minimum < l1

# using MultistartOptimization
# sol = solve(prob, MultistartOptimization.TikTak(100), local_method = NLopt.LD_LBFGS)
# @test 10*sol.minimum < l1

# using QuadDIRECT
# sol = solve(prob, QuadDirect(); splits = ([-0.5, 0.0, 0.5],[-0.5, 0.0, 0.5]))
# @test 10*sol.minimum < l1

using Evolutionary
sol = solve(prob, CMAES(μ =40 , λ = 100),abstol=1e-15)
@test 10*sol.minimum < l1

#=
using BlackBoxOptim
prob = GalacticOptim.OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, BBO())
@test 10*sol.minimum < l1

prob_multi = EnsembleOptimizationProblem([prob, prob])
sol_multi = solve(prob_multi, BBO(:borg_moea), FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true), ϵ=0.05, MaxSteps=500000)
@test 10*sum(sol_multi.minimum) < l1
=#

using ModelingToolkit
rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p = [1.0, 100.0]

f = OptimizationFunction(rosenbrock,ModelingToolkit.AutoModelingToolkit(),x0,_p,grad=true,hess=true)
prob = OptimizationProblem(f,x0,_p)
sol = solve(prob,Optim.Newton())


text:
using MonteCarlo
using Debugger

betas = (1.0)#(1.0, 2.0, 5.0, 6.0, 7.0)
#lattice = SquareLattice(4)
mus = (-0.4 )#vcat(-2.0:0.5:-0.5, -0.1:0.1:1.1, 1.25, 1.5, 2.0)

dqmcs = []
    counter = 0
    N = length(mus) * length(betas)
    @time for beta in betas, mu in mus
        counter += 1
        print("\r[", lpad("$counter", 2), "/$N]")
        m = HubbardModelBi(2,1, a=0.85)#(l = lattice, t = 1.0, U = 4.0, μ = mu)
        dqmc = DQMC(
            m, beta = beta, delta_tau = 0.125, safe_mult = 8,
            thermalization = 2000, sweeps = 10000, measure_rate = 1,
            recorder = Discarder()
        )
        dqmc[:occ] = occupation(dqmc, m)
        #dqmc[:PC] = pairing_correlation(dqmc, m, kernel = MonteCarlo.pc_kernel)
      # # run!(dqmc, verbose = false)

        # for simplicity we just keep the whole simulation around
      # # push!(dqmcs, dqmc)
    end
##mean(mean(dqmcs[1][:occ]))


text:
function gridworld_main()
    gridworld = create_example_gridworld()
    planner = GridworldPlanner(gridworld, 100)
    agent = GridworldAgent(gridworld, planner)
    run(agent)
end


text:
export
    PMFGaussian

struct PMFGaussian
    r::Real   # Rating
    sigma::Real
    d::Int   # Dimension of latent space
end

function gradloglik(g::PMFGaussian, x::Vector{<:Real})
    u = x[1:g.d]
    v = x[g.d+1:end]
    e = dot(u, v) - g.r
    return -e / g.sigma^2 * [v; u]
end

function loglik(g::PMFGaussian, x::Vector{<:Real})
    u = x[1:g.d]
    v = x[g.d+1:end]
    e = dot(u, v) - g.r
    return 0.5(e / g.sigma)^2
end

# --------------------------------------------------------------------------
# HELPER FUNCTIONS FOR THE BPS sampler

function pmf_base(rexp::Real, p::Poly)
    # find solutions to [p(x) = y]
    rs  = roots(p - rexp)
    # we want the first real root > 0.0
    # this requires a bit of subtlety as `roots` is not
    # super stable so there are imaginary numbers trailing
    #
    # start by dumping roots with a negative real parts
    # and roots that are "too imaginary"
    prs = [r for r in rs if (real(r) > 0.0) && (imag(r) / abs(r) < 1e-8) ]
    prs = sort(abs.(prs))
    # now pick the one with the smallest radius
    k   = searchsortedfirst(prs, 0.0)
    return k > length(prs) ? 0.0 : prs[k]
end

function pmf_caseA(rexp::Real, p::Poly)
    # == CASE A ###### intensity:
    #  |   |
    #  |  /
    #  |/
    #  |________________
    ##############################
    # the branch starts on y axis
    # can directly look for root
    ##############################
    return pmf_base(rexp, polyint(p))
end

function pmf_caseB(rexp::T, p::Poly, r::T) where T <: Real
    @assert 0.0 < r
    # == CASE B ###### intensity:
    #   |         |
    #   |        |
    #   |      /
    #   |____r__________
    ##############################################
    # the branch starts after 0 at r
    # move the branch to zero, then add t0 to root
    ##############################################
    return r + pmf_base(rexp, polyint(polyval(p, Poly([r, 1.0]))))
end

function pmf_caseC(rexp::T, p::Poly, t0::T, tp::T) where T <: Real
    @assert 0.0 < t0 < tp
    # == CASE C ###### intensity:
    #   |             |
    #   |\           |
    #   |  \       /
    #   |___t0___tp______
    ##############################################
    # partial hump then flat then branch
    # the integral has a "step", should check on
    # what step we are
    ##############################################
    tau = 0.0
    Ip = polyint(p)
    step = Ip(t0)
    if rexp <= step # solve for 0.0->t0
        tau = pmf_base(rexp, Ip)
    else # solve for after tp
        tau = tp + pmf_base(rexp-step, polyint(polyval(p, Poly([tp, 1.0]))))
    end
    return tau
end

function pmf_caseD(rexp::T, p::Poly, tm::T, t0::T, tp::T
                    ) where T <: Real
    @assert 0.0 <= tm < t0 < tp
    # == CASE D ###### intensity:
    #   |
    #   |       _          |
    #   |     /  \       /
    #   |___tm___t0____tp_____
    ##############################################
    # flat then complete hump then flat then branch
    # the integral has a two "steps", should check
    # on what step we are
    ##############################################
    # start by shifting by tm
    Ip1  = polyint(polyval(p, Poly([tm, 1.0])))
    step = Ip1(t0-tm)
    if rexp <= step # solve between tm and t0
        tau = tm + pmf_base(rexp, Ip1)
    else # solve after tp
        tau = tp + pmf_base(rexp-step, polyint(polyval(p, Poly([tp, 1.0]))))
    end
end


text:
using Coconet
Coconet.train_model(
    # Load your model
    model = nothing,
    opt = nothing,
    losses = [(5e6, 5e6)],

    # Build model to specs
    head_count = 8,
    latent_size = 32,
    N = 6,

    # Speed it up
    sequence_length = 128,
    batch_size = 5,
    epochs = 20,
    set_length = nothing,
    max_regions = nothing,

    # Difficulty curve
    scramble_rate = 0.1,
    scramble_increase = 0.1,
    scramble_max = 0.9,
    graduation_rate = 0.5,

    # Data
    dataset_name = "Tonnetz_Bach",
    word_size = 3,
)

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet


text:
cd(Pkg.dir("GtkBuilderAid"))
Pkg.add("Coverage")
using Coverage
# Push results to coveralls
Coveralls.submit(Coveralls.process_folder())


text:
# Data structure

export ModeTree

abstract AbstractModeTree

typealias ModeDict Dict{Tree,Vector{Mode}}
immutable ModeTree <: AbstractModeTree
    tree::Tree
    free_modes::ModeDict
end
immutable SquaredModeTree <: AbstractModeTree
    mtree::ModeTree
end

# Basic functions

export tree, unsquare

tree(t::ModeTree) = t.tree
tree(t::SquaredModeTree) = tree(t.mtree)

Tensors.square(t::ModeTree) = SquaredModeTree(t)
unsquare(t::AbstractModeTree) = t
unsquare(t::SquaredModeTree) = t.mtree

for f in (:getindex, :setindex!) @eval $f(x::ModeTree, args...) = $f(x.free_modes, args...) end
getindex(x::SquaredModeTree, v) = square(x.mtree[v])

# Tree construction

export ModeTree!, link!

function ModeTree!(free_modes::Vector{Mode}, children::ModeTree...)
    r = Tree!([tree(u) for u in children]...)
    D = merge!(ModeDict(r => copy(free_modes)), [u.free_modes for u in children]...)
    return ModeTree(r,D)
end
function ModeTree!(free_modes::Vector{Mode})
    r = Tree()
    D = ModeDict(r => free_modes)
    return ModeTree(r,D)
end
ModeTree!(children::ModeTree...) = ModeTree!(Mode[], children...)
ModeTree(free_modes::Vector{Mode}, children::ModeTree...) = ModeTree!(free_modes, deepcopy(children)...)
ModeTree(free_modes::Vector{Mode}) = ModeTree!(free_modes)
ModeTree(children::ModeTree...) = ModeTree(Mode[], children...)

function link!(v::ModeTree, u::ModeTree)
    link!(tree(v), tree(u))
    merge!(v.free_modes, u.free_modes)
    return v
end

# Common tree structures

export lintree

function lintree{T <: AbstractVector{Mode}}(D::AbstractVector{T})
    mtree = ModeTree(D[1])
    for i = 2:length(D)
        mtree = ModeTree!(D[i], mtree)
    end
    return mtree
end
lintree(D::AbstractVector{Mode}) = lintree([[k] for k in D])

# Other

function Tensors.mode(t::ModeTree)
    modes = Mode[]
    for D in values(t.free_modes) append!(modes, D) end
    return modes
end
Tensors.mode(t::SquaredModeTree) = square(mode(t.mtree))


text:
include("../src/NLEVP/perturbation.jl")
println("Building WavesAndEigenvalues.jl...")

if "JULIA_WAE_PERT_ORDER" in keys(ENV)
    N=ENV["JULIA_WAE_PERT_ORDER"]
    N=parse(Int,N)
else
    N=16
end
println(" for order $N...")
disk_MN(N)
println("Done!")


text:
using Trajectories
using Tables
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

t = 0:10
x = rand(11)
X = trajectory(t, x)

@test keys(X) === t
@test values(X) === x
@test eltype(X) === eltype(x)
@test eltype(X) === eltype(x)

@test Pair(X) === Pair(t, x)
@test [pair for pair in pairs(X)] == map(t->Pair(t...), collect(zip(t, x)))

@test trajectory(zip(t, x)) == X

@test get(X, 5) == x[6]
@test_throws ErrorException get(X, 5.5)
@test_throws ErrorException get(X, 11)

#@test [y for y in X] == [t, x]

for i in 5:10
    @test Trajectories._find(5:10, i) == i - 4
end
@test_throws ErrorException Trajectories._find(5:10, 4)
@test_throws ErrorException Trajectories._find(5:10, 11)

@test interpolate(Linear(), X, 0) == interpolate(Linear(), X, -0.5) == x[1]
@test interpolate(Linear(), X, 10) == interpolate(Linear(), X, 10.5) == x[end]
@test interpolate(Linear(), X, 1/3) ≈ (2x[1] + x[2])/3

@test interpolate(Left(), X, 0) == interpolate(Left(), X, 0.5) == x[1]
@test interpolate(Left(), X, 10) == interpolate(Linear(), X, 10.5) == x[end]
@test_throws BoundsError interpolate(Left(), X, -1/3)

@test piecewise(Trajectory([1,2,3], [1,2])) == ([1, 2, 2, 3], [1, 1, 2, 2])

@test searchsorted(sort([1,missing, NaN, missing]), NaN) == 2:2
@test searchsorted(sort([1,missing, NaN, missing]), missing) == 3:4

@test Tables.schema(Tables.columns(trajectory(1:10, rand(10)))) == Tables.Schema{(:t, :x),Tuple{Int,Float64}}()

Y = map(collect, X)
push!(Y, 12=>0.0)
@test keys(Y) == [t; 12]
@test values(Y) == [x; 0.0]

Y = copy(X)
Y = vcat(Y, Y)
@test collect(pairs(Y))[end÷2+1:end] == collect(pairs(X))


text:
using Test

@testset "tests" begin
  include(joinpath("advection", "runtests.jl"))
  include(joinpath("shallow_water", "runtests.jl"))
  include(joinpath("euler", "runtests.jl"))
end

@testset "experiments" begin
  experiments_dir = joinpath(@__DIR__, "..", "experiment")
  for (root, dir, files) in walkdir(experiments_dir)
    jlfiles = filter(s -> endswith(s, ".jl"), files)
    for jlfile in jlfiles
      path = joinpath(root, jlfile)
      modname = Symbol(jlfile)
      @eval module $modname
        using Test
        const _testing = true
        @testset $jlfile begin
          include($path)
        end
      end
    end
  end
end


text:
using XRayQuadrics
using LinearAlgebra
using Colors
using Plots

axis = [0;0;1]
center = [0;0;0]

R = 1
b = 1

s = Plane(center, axis)
qs = Quadric(s)

c = Cylinder(R, center, axis)
qc = Quadric(c)

p = Paraboloid(R, center, axis)
qp = Quadric(p)

h = Hyperboloid(R,b, center, axis)
qh = Quadric(h)

p1 = Plane(c.a*0.01, axis)
p2 = Plane(c.a*2, axis)

tqc = TruncatedQuadric(c, [p1, p2], true)
tqp = TruncatedQuadric(p, [p1, p2], true)
tqh = TruncatedQuadric(h, [p1, p2], true)

gr()
splot = Plots.surface(tqp, color=:blue, plot_caps=true, size=(500,500))
display(splot)
Plots.surface!(tqh, color=:red)
Plots.surface!(tqc, color=:green)
Plots.plot!([center center+3*axis][1,:], [center center+3*axis][2,:], [center center+3*axis][3,:], color=:blue, linewidth=4)


text:
using PyPlot

export plot_trajectories, plot_deviation, plot_trajectories_and_deviation,
		plot_responses, plot_responses_timeslice, plot_hatR, plot_hatR_diag,
		plot_hatTheta, cut_res_tf, cut_resp_tf

"""
	cut_res_tf(res::Result, tf)

	Only time grid up to a specified final time tf which has to be
	a part of the vector res.time

	Returns the Result struct up to time tf

"""
function cut_res_tf(res::Result, tf)
	ts = findall(x->x==tf, res.time)
	length(ts) != 1 ? error("No matching time slice!") : t = ts[1]
	return Result(res.time[1:t], res.data[:,1:t])
end

"""
	cut_resp_tf(resp::Responses, tf)

	Only time grid up to a specified final time tf which has to be
	a part of the vector res.time

	Returns the Response struct up to time tf

"""
function cut_resp_tf(resp::Responses, tf)
	ts = findall(x->x==tf, resp.time)
	length(ts) != 1 ? error("No matching time slice!") : t = ts[1]
	return Responses(resp.time[1:t], resp.data[:,1:t,1:t])
end

"""
	plot_trajectories(res::Result)

	Plots the trajectories for all species

	Use plot_trajectories(res::Result, tf) for plotting up to final time tf

"""
function plot_trajectories(res::Result)

	num_species = size(res.data,1)
	figure(figsize=(20,5))

	for i in 1:num_species
		plot(res.time, res.data[i,:], linestyle="dashed")
	end

	legend(1:num_species)
	xlabel("Time")
	ylabel("Copy Numbers")

end

function plot_trajectories(res::Result, tf)
	plot_trajectories(cut_res_tf(res, tf))
end

"""
	plot_deviation(res::Result, ref::Result)

	Plots the deviation of res for all species to the reference ref
	(in most cases e.g. Gillespie)

	deviation = | res - ref | / ref

	Use plot_deviation(res::Result, ref::Result, tf) for plotting up to
	final time tf

"""
function plot_deviation(res::Result, ref::Result)

	num_species = size(res.data,1)
	figure(figsize=(20,5))

	for i in 1:num_species
		plot(res.time, abs.(res.data[i,:] .- ref.data[i,:]) ./ ref.data[i,:], linestyle="dashed")
	end

	legend(1:num_species)
	xlabel("Time")
	ylabel("Relative deviation from reference")

end

function plot_deviation(res::Result, ref::Result, tf)
	plot_deviation(cut_res_tf(res, tf), cut_res_tf(ref, tf))
end

"""
	plot_trajectories_and_deviation(res::Result, ref::Result)

	Plots the trajectories and deviaton for all species.

	For more details, look at `plot_trajectories(...)`
	and `plot_deviation(...)`

	Use plot_trajectories_and_deviation(res::Result, ref::Result, tf) for
	plotting up to final time tf

"""
function plot_trajectories_and_deviation(res::Result, ref::Result)

	num_species = size(res.data,1)
	figure(figsize=(20,5))

	subplot(121)
	for i in 1:num_species
		plot(res.time, res.data[i,:], linestyle="dashed")
	end
	legend(1:num_species)
	xlabel("Time")
	ylabel("Copy Numbers")

	subplot(122)
	for i in 1:num_species
		plot(res.time, abs.(res.data[i,:] .- ref.data[i,:]) ./ ref.data[i,:], linestyle="dashed")
	end
	legend(1:num_species)
	xlabel("Time")
	ylabel("Relative deviation from reference")

end

function plot_trajectories_and_deviation(res::Result, ref::Result, tf)
	plot_trajectories_and_deviation(cut_res_tf(res, tf), cut_res_tf(ref, tf))
end

"""
	plot_responses(resp::Responses)

	Plots the response functions for all species

	Use plot_responses(resp::Responses, tf) for plotting up to final time tf

"""
function plot_responses(resp::Responses)

	num_species = size(resp.data,1)

	ncols = 3
	vmin = minimum(resp.data)
	vmax = maximum(resp.data)
	lims = [resp.time[1], resp.time[end], resp.time[end], resp.time[1]]

	rr = copy(resp.data)

	for i in 1:num_species
		for j in 1:length(resp.time)
			rr[i,j,:] = circshift(resp.data[i,j,end:-1:1],j)
		end
	end

	fig, axes = subplots(figsize=(20,5), ncols=ncols, nrows=ceil(Int, num_species/ncols), sharey=true)

	for i in 1:num_species
		img = axes[i].imshow(rr[i,:,:], aspect="equal", extent=lims, vmin=vmin, vmax=vmax)
		axes[i].set_xlabel(L"$\tau - \tau^\prime$")
		axes[i].set_title(L"$ R(\tau,\tau -\tau^\prime) $"*" for species $(i)")
		i == ncols ? fig.colorbar(img, ax=vec(axes)) : nothing
	end

end

function plot_responses(resp::Responses, tf)
	plot_responses(cut_resp_tf(resp, tf))
end

"""
	plot_responses_timeslice(resp::Responses, timeslice)

	Plots the response for all species at a given timeslice

	R(t=timeslice, t')

"""
function plot_responses_timeslice(resp::Responses, timeslice)

	num_species = size(resp.data,1)
	ts = findall(x->x==timeslice, resp.time)
	length(ts) != 1 ? error("No matching time slice!") : t = ts[1]
	figure(figsize=(20,5))

	for i in 1:num_species
		plot(resp.time, resp.data[i,t,:], linestyle="dashed")
	end

	legend(1:num_species)
	xlabel("Time")
	ylabel(L"$ R(\tau="*string(timeslice)*L",\tau^\prime) $")

end

"""
	plot_hatR(tspan, fields::Fields_quad2)

	Plots the field hatR2 for the quadratic order parameter and alpha squared
	case.

	Needs improvement if number of species is not divisible by 3.

"""
function plot_hatR(tspan, fields::Fields_quad2)

	hatR = fields.hatR2
	num_species = size(hatR,1)

	ncols = 3
	vmin = minimum(hatR)
	vmax = maximum(hatR)
	vmin == 0 && vmax == 0 ? vmax = 1.0 : nothing
	lims = [tspan[1], tspan[end], tspan[end], tspan[1]]

	rr = copy(hatR)

	for i in 1:num_species
		for j in 1:length(tspan)
			rr[i,j,:] = circshift(hatR[i,j,end:-1:1],j)
		end
	end

	fig, axes = subplots(figsize=(20,5), ncols=ncols, nrows=ceil(Int, num_species/ncols), sharey=true)

	for i in 1:num_species
		img = axes[i].imshow(rr[i,:,:], aspect="equal", extent=lims, vmin=vmin, vmax=vmax)
		axes[i].set_xlabel(L"$\tau - \tau^\prime$")
		axes[i].set_title(L"$\hat{R}_2(\tau,\tau -\tau^\prime)$"*" for species $(i)")
		i == ncols ? fig.colorbar(img, ax=vec(axes)) : nothing
	end

end

"""
	plot_hatR_diag(tspan, fields)

	Plots the hatR1 diagonal elements R1(τ,τ-Δt)

"""
function plot_hatR_diag(tspan, fields)

	hatR = fields.hatR1
	num_species = size(hatR,1)
	dt = tspan[2] - tspan[1]

    figure(figsize=(20,5))
    for i in 1:num_species
        plot(tspan[1:end-1], hatR[i,1:end-1]*dt)
    end
    legend(1:num_species)
    xlabel("Time")
    ylabel(L"$\hat{R}_1(\tau,\tau_-)\,\Delta t$")

end

# Plot hatTheta1
function plot_hatTheta1(tspan, fields)

	num_species = size(fields.hatTheta1,1)
	figure(figsize=(20,5))

	for i in 1:num_species
		plot(tspan[1:end-1], fields.hatTheta1[i,1:end-1])
	end
	legend(1:num_species)
	xlabel("Time")
	ylabel(L"$\tilde{\theta}_1$")

end

"""
	plot_hatTheta(tspan, fields::Fields_lin1)

	Plots the hat Theta field for linear order parameters and linear alpha

"""
function plot_hatTheta(tspan, fields::Fields_lin1)
	plot_hatTheta1(tspan, fields)
end

"""
	plot_hatTheta(tspan, fields::Fields_quad1)

	Plots the hat Theta field for quadratic order parameters and linear alpha

"""
function plot_hatTheta(tspan, fields::Fields_quad1)
	plot_hatTheta1(tspan, fields)
end

# Plot hatTheta1 and hatTheta2
function plot_hatTheta12(tspan, fields)

	num_species = size(fields.hatTheta1,1)
	figure(figsize=(20,5))

	subplot(121)
	for i in 1:num_species
		plot(tspan[1:end-1], fields.hatTheta1[i,1:end-1])
	end
	legend(1:num_species)
	xlabel("Time")
	ylabel(L"$\tilde{\theta}_1$")

	subplot(122)
	for i in 1:num_species
		plot(tspan[1:end-1], fields.hatTheta2[i,1:end-1])
	end
	legend(1:num_species)
	xlabel("Time")
	ylabel(L"$\tilde{\theta}_2$")

end

"""
	plot_hatTheta(tspan, fields::Fields_lin2)

	Plots the hat Theta fields 1 and 2 for linear order parameters and
	alpha squared

"""
function plot_hatTheta(tspan, fields::Fields_lin2)
	plot_hatTheta12(tspan, fields)
end

"""
	plot_hatTheta(tspan, fields::Fields_quad2)

	Plots the hat Theta fields 1 and 2 for quadratic order parameters and
	alpha squared

"""
function plot_hatTheta(tspan, fields::Fields_quad2)
	plot_hatTheta12(tspan, fields)
end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
using Combinatorics
using ITensors
using Random
using LinearAlgebra: diag
using Test

Random.seed!(1234)

include("util.jl")

@testset "MPS Basics" begin
  N = 10
  sites = [Index(2, "Site") for n in 1:N]
  psi = MPS(sites)
  @test length(psi) == N
  @test length(MPS()) == 0
  @test linkdims(psi) == fill(1, N - 1)
  @test isnothing(flux(psi))

  psi = MPS(sites; linkdims=3)
  @test length(psi) == N
  @test length(MPS()) == 0
  @test linkdims(psi) == fill(3, N - 1)
  @test isnothing(flux(psi))

  str = split(sprint(show, psi), '\n')
  @test str[1] == "MPS"
  @test length(str) == length(psi) + 2

  @test siteind(psi, 2) == sites[2]
  @test findfirstsiteind(psi, sites[2]) == 2
  @test findfirstsiteind(psi, sites[4]) == 4
  @test findfirstsiteinds(psi, IndexSet(sites[5])) == 5
  @test hasind(psi[3], linkind(psi, 2))
  @test hasind(psi[3], linkind(psi, 3))

  @test isnothing(linkind(psi, N))
  @test isnothing(linkind(psi, N + 1))
  @test isnothing(linkind(psi, 0))
  @test isnothing(linkind(psi, -1))
  @test linkind(psi, 3) == commonind(psi[3], psi[4])

  psi[1] = ITensor(sites[1])
  @test hasind(psi[1], sites[1])

  @testset "N=1 MPS" begin
    sites1 = [Index(2, "Site,n=1")]
    psi = MPS(sites1)
    @test length(psi) == 1
    @test siteind(psi, 1) == sites1[1]
    @test siteinds(psi)[1] == sites1[1]
  end

  @testset "Missing links" begin
    psi = MPS([randomITensor(sites[i]) for i in 1:N])
    @test isnothing(linkind(psi, 1))
    @test isnothing(linkind(psi, 5))
    @test isnothing(linkind(psi, N))
    @test maxlinkdim(psi) == 1
    @test psi ⋅ psi ≈ *(dag(psi)..., psi...)[]
  end

  @testset "productMPS" begin
    @testset "vector of string input" begin
      sites = siteinds("S=1/2", N)
      state = fill("", N)
      for j in 1:N
        state[j] = isodd(j) ? "Up" : "Dn"
      end
      psi = MPS(sites, state)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
      psi = productMPS(sites, state)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
      @test_throws DimensionMismatch MPS(sites, fill("", N - 1))
      @test_throws DimensionMismatch productMPS(sites, fill("", N - 1))
    end

    @testset "String input" begin
      sites = siteinds("S=1/2", N)
      psi = MPS(sites, "Dn")
      for j in 1:N
        sign = -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
      psi = productMPS(sites, "Dn")
      for j in 1:N
        sign = -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end

      psi = MPS(sites, "X+")
      for j in 1:N
        @test (psi[j] * op(sites, "X", j) * dag(prime(psi[j], "Site")))[] ≈ 1.0
      end
    end

    @testset "Int input" begin
      sites = siteinds("S=1/2", N)
      psi = MPS(sites, 2)
      for j in 1:N
        sign = -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
      psi = productMPS(sites, 2)
      for j in 1:N
        sign = -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
    end

    @testset "vector of int input" begin
      sites = siteinds("S=1/2", N)
      state = fill(0, N)
      for j in 1:N
        state[j] = isodd(j) ? 1 : 2
      end
      psi = MPS(sites, state)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
      psi = productMPS(sites, state)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
    end

    @testset "vector of ivals input" begin
      sites = siteinds("S=1/2", N)
      states = fill(0, N)
      for j in 1:N
        states[j] = isodd(j) ? 1 : 2
      end
      ivals = [sites[n] => states[n] for n in 1:length(sites)]
      psi = MPS(ivals)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end
      psi = productMPS(ivals)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
      end

      @testset "ComplexF64 eltype" begin
        sites = siteinds("S=1/2", N)
        psi = MPS(ComplexF64, sites, fill(1, N))
        for j in 1:N
          @test eltype(psi[j]) == ComplexF64
        end
        psi = productMPS(ComplexF64, sites, fill(1, N))
        for j in 1:N
          @test eltype(psi[j]) == ComplexF64
        end
        @test eltype(psi) == ITensor
        @test ITensors.promote_itensor_eltype(psi) == ComplexF64
      end
    end

    @testset "N=1 case" begin
      site = Index(2, "Site,n=1")
      psi = MPS([site], [1])
      @test psi[1][1] ≈ 1.0
      @test psi[1][2] ≈ 0.0
      psi = productMPS([site], [1])
      @test psi[1][1] ≈ 1.0
      @test psi[1][2] ≈ 0.0
      psi = MPS([site], [2])
      @test psi[1][1] ≈ 0.0
      @test psi[1][2] ≈ 1.0
      psi = productMPS([site], [2])
      @test psi[1][1] ≈ 0.0
      @test psi[1][2] ≈ 1.0
    end
  end

  @testset "randomMPS with chi==1" begin
    phi = randomMPS(sites)
    phic = randomMPS(ComplexF64, sites)

    @test maxlinkdim(phi) == 1
    @test maxlinkdim(phic) == 1

    @test hasind(phi[1], sites[1])
    @test norm(phi[1]) ≈ 1.0
    @test norm(phic[1]) ≈ 1.0

    @test hasind(phi[4], sites[4])
    @test norm(phi[4]) ≈ 1.0
    @test norm(phic[4]) ≈ 1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1]) * psi[1]
    for j in 2:N
      phipsi *= dag(phi[j]) * psi[j]
    end
    @test phipsi[] ≈ inner(phi, psi)

    badsites = [Index(2) for n in 1:(N + 1)]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(phi, badpsi)
  end

  @testset "broadcasting" begin
    psi = randomMPS(sites)
    orthogonalize!(psi, 1)
    @test ortho_lims(psi) == 1:1
    @test dim.(psi) == fill(2, length(psi))
    psi′ = prime.(psi)
    @test ortho_lims(psi′) == 1:length(psi′)
    @test ortho_lims(psi) == 1:1
    for n in 1:length(psi)
      @test prime(psi[n]) == psi′[n]
    end
    psi_copy = copy(psi)
    psi_copy .= addtags(psi_copy, "x")
    @test ortho_lims(psi_copy) == 1:length(psi_copy)
    @test ortho_lims(psi) == 1:1
    for n in 1:length(psi)
      @test addtags(psi[n], "x") == psi_copy[n]
    end
  end

  @testset "copy and deepcopy" begin
    s = siteinds("S=1/2", 3)
    M1 = randomMPS(s; linkdims=3)
    @test norm(M1) ≈ 1

    M2 = deepcopy(M1)
    M2[1] .*= 2 # Modifies the tensor data
    @test norm(M1) ≈ 1
    @test norm(M2) ≈ 2

    M3 = copy(M1)
    M3[1] *= 3
    @test norm(M1) ≈ 1
    @test norm(M3) ≈ 3

    M4 = copy(M1)
    M4[1] .*= 4
    @test norm(M1) ≈ 4
    @test norm(M4) ≈ 4
  end

  @testset "inner same MPS" begin
    psi = randomMPS(sites)
    psidag = dag(psi)
    #ITensors.prime_linkinds!(psidag)
    psipsi = psidag[1] * psi[1]
    for j in 2:N
      psipsi *= psidag[j] * psi[j]
    end
    @test psipsi[] ≈ inner(psi, psi)
  end

  @testset "norm MPS" begin
    psi = randomMPS(sites; linkdims=10)
    psidag = sim(linkinds, dag(psi))
    psi² = ITensor(1)
    for j in 1:N
      psi² *= psidag[j] * psi[j]
    end
    @test psi²[] ≈ psi ⋅ psi
    @test sqrt(psi²[]) ≈ norm(psi)

    psi = randomMPS(sites; linkdims=10)
    psi .*= 1:N
    @test norm(psi) ≈ factorial(N)

    psi = randomMPS(sites; linkdims=10)
    for j in 1:N
      psi[j] .*= j
    end
    # This fails because it modifies the MPS ITensors
    # directly, which ruins the orthogonality
    @test norm(psi) ≉ factorial(N)
    reset_ortho_lims!(psi)
    @test norm(psi) ≈ factorial(N)
  end

  @testset "lognorm MPS" begin
    psi = randomMPS(sites; linkdims=10)
    for j in 1:N
      psi[j] .*= j
    end
    psidag = sim(linkinds, dag(psi))
    psi² = ITensor(1)
    for j in 1:N
      psi² *= psidag[j] * psi[j]
    end
    @test psi²[] ≈ psi ⋅ psi
    @test 0.5 * log(psi²[]) ≉ lognorm(psi)
    @test lognorm(psi) ≉ log(factorial(N))
    # Need to manually change the orthogonality
    # limits back to 1:length(psi)
    reset_ortho_lims!(psi)
    @test 0.5 * log(psi²[]) ≈ lognorm(psi)
    @test lognorm(psi) ≈ log(factorial(N))
  end

  @testset "scaling MPS" begin
    psi = randomMPS(sites; linkdims=4)
    twopsidag = 2.0 * dag(psi)
    #ITensors.prime_linkinds!(twopsidag)
    @test inner(twopsidag, psi) ≈ 2.0 * inner(psi, psi)
  end

  @testset "flip sign of MPS" begin
    psi = randomMPS(sites)
    minuspsidag = -dag(psi)
    #ITensors.primelinkinds!(minuspsidag)
    @test inner(minuspsidag, psi) ≈ -inner(psi, psi)
  end

  @testset "add MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    # sum of many MPSs
    Ks = [randomMPS(sites) for i in 1:3]
    K12 = add(Ks[1], Ks[2])
    K123 = add(K12, Ks[3])
    @test inner(sum(Ks), K123) ≈ inner(K123, K123)
  end

  @testset "+ MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = psi + phi
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    # sum of many MPSs
    Ks = [randomMPS(sites) for i in 1:3]
    K12 = Ks[1] + Ks[2]
    K123 = K12 + Ks[3]
    @test inner(sum(Ks), K123) ≈ inner(K123, K123)
  end

  @testset "+ MPS with coefficients" begin
    Random.seed!(1234)

    N = 20
    conserve_qns = true

    s = siteinds("S=1/2", N; conserve_qns=conserve_qns)
    state = n -> isodd(n) ? "↑" : "↓"

    ψ₁ = randomMPS(s, state; linkdims=4)
    ψ₂ = randomMPS(s, state; linkdims=4)
    ψ₃ = randomMPS(s, state; linkdims=4)

    ψ = ψ₁ + ψ₂

    @test inner(ψ, ψ) ≈ inner_add(ψ₁, ψ₂)
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂)

    ψ = +(ψ₁, ψ₂; cutoff=0.0)

    @test inner(ψ, ψ) ≈ inner_add(ψ₁, ψ₂)
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂)

    ψ = ψ₁ + (-ψ₂)

    @test inner(ψ, ψ) ≈ inner_add((1, ψ₁), (-1, ψ₂))
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂)

    α₁ = 2.2
    α₂ = -4.1
    ψ = +(α₁ * ψ₁, α₂ * ψ₂; cutoff=1e-8)

    @test inner(ψ, ψ) ≈ inner_add((α₁, ψ₁), (α₂, ψ₂))
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂)

    α₁ = 2 + 3im
    α₂ = -4 + 1im
    ψ = α₁ * ψ₁ + α₂ * ψ₂

    @test inner(ψ, ψ) ≈ inner_add((α₁, ψ₁), (α₂, ψ₂))
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂)

    α₁ = 2 + 3im
    α₂ = -4 + 1im
    ψ = α₁ * ψ₁ + α₂ * ψ₂ + ψ₃

    @test inner(ψ, ψ) ≈ inner_add((α₁, ψ₁), (α₂, ψ₂), ψ₃)
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂) + maxlinkdim(ψ₃)

    ψ = ψ₁ - ψ₂

    @test inner(ψ, ψ) ≈ inner_add(ψ₁, (-1, ψ₂))
    @test maxlinkdim(ψ) ≤ maxlinkdim(ψ₁) + maxlinkdim(ψ₂)
  end

  sites = siteinds(2, N)
  psi = MPS(sites)
  @test length(psi) == N # just make sure this works
  @test length(siteinds(psi)) == N

  psi = randomMPS(sites)
  l0s = linkinds(psi)
  orthogonalize!(psi, N - 1)
  ls = linkinds(psi)
  for (l0, l) in zip(l0s, ls)
    @test tags(l0) == tags(l)
  end
  @test ITensors.leftlim(psi) == N - 2
  @test ITensors.rightlim(psi) == N
  orthogonalize!(psi, 2)
  @test ITensors.leftlim(psi) == 1
  @test ITensors.rightlim(psi) == 3
  psi = randomMPS(sites)
  ITensors.setrightlim!(psi, N + 1) # do this to test qr
  # from rightmost tensor
  orthogonalize!(psi, div(N, 2))
  @test ITensors.leftlim(psi) == div(N, 2) - 1
  @test ITensors.rightlim(psi) == div(N, 2) + 1

  @test isnothing(linkind(MPS(fill(ITensor(), N), 0, N + 1), 1))

  @testset "replacebond!" begin
    # make sure factorization preserves the bond index tags
    psi = randomMPS(sites)
    phi = psi[1] * psi[2]
    bondindtags = tags(linkind(psi, 1))
    replacebond!(psi, 1, phi)
    @test tags(linkind(psi, 1)) == bondindtags

    # check that replacebond! updates llim and rlim properly
    orthogonalize!(psi, 5)
    phi = psi[5] * psi[6]
    replacebond!(psi, 5, phi; ortho="left")
    @test ITensors.leftlim(psi) == 5
    @test ITensors.rightlim(psi) == 7

    phi = psi[5] * psi[6]
    replacebond!(psi, 5, phi; ortho="right")
    @test ITensors.leftlim(psi) == 4
    @test ITensors.rightlim(psi) == 6

    ITensors.setleftlim!(psi, 3)
    ITensors.setrightlim!(psi, 7)
    phi = psi[5] * psi[6]
    replacebond!(psi, 5, phi; ortho="left")
    @test ITensors.leftlim(psi) == 3
    @test ITensors.rightlim(psi) == 7
  end
end

@testset "orthogonalize! with QNs" begin
  N = 8
  sites = siteinds("S=1/2", N; conserve_qns=true)
  init_state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  psi0 = MPS(sites, init_state)
  orthogonalize!(psi0, 4)
  @test ITensors.leftlim(psi0) == 3
  @test ITensors.rightlim(psi0) == 5
end

# Helper function for making MPS
function basicRandomMPS(N::Int; dim=4)
  sites = [Index(2, "Site") for n in 1:N]
  M = MPS(sites)
  links = [Index(dim, "n=$(n-1),Link") for n in 1:(N + 1)]
  for n in 1:N
    M[n] = randomITensor(links[n], sites[n], links[n + 1])
  end
  M[1] *= delta(links[1])
  M[N] *= delta(links[N + 1])
  M[1] /= sqrt(inner(M, M))
  return M
end

@testset "MPS gauging and truncation" begin
  N = 30

  @testset "orthogonalize! method" begin
    c = 12
    M = basicRandomMPS(N)
    orthogonalize!(M, c)

    @test ITensors.leftlim(M) == c - 1
    @test ITensors.rightlim(M) == c + 1

    # Test for left-orthogonality
    L = M[1] * prime(M[1], "Link")
    l = linkind(M, 1)
    @test norm(L - delta(l, l')) < 1E-12
    for j in 2:(c - 1)
      L = L * M[j] * prime(M[j], "Link")
      l = linkind(M, j)
      @test norm(L - delta(l, l')) < 1E-12
    end

    # Test for right-orthogonality
    R = M[N] * prime(M[N], "Link")
    r = linkind(M, N - 1)
    @test norm(R - delta(r, r')) < 1E-12
    for j in reverse((c + 1):(N - 1))
      R = R * M[j] * prime(M[j], "Link")
      r = linkind(M, j - 1)
      @test norm(R - delta(r, r')) < 1E-12
    end

    @test norm(M[c]) ≈ 1.0
  end

  @testset "truncate! method" begin
    M = basicRandomMPS(N; dim=10)
    M0 = copy(M)
    truncate!(M; maxdim=5)

    @test ITensors.rightlim(M) == 2

    # Test for right-orthogonality
    R = M[N] * prime(M[N], "Link")
    r = linkind(M, N - 1)
    @test norm(R - delta(r, r')) < 1E-12
    for j in reverse(2:(N - 1))
      R = R * M[j] * prime(M[j], "Link")
      r = linkind(M, j - 1)
      @test norm(R - delta(r, r')) < 1E-12
    end

    @test inner(M, M0) > 0.1
  end
end

@testset "Other MPS methods" begin
  @testset "sample! method" begin
    N = 10
    sites = [Index(3, "Site,n=$n") for n in 1:N]
    psi = randomMPS(sites; linkdims=3)
    nrm2 = inner(psi, psi)
    psi[1] *= (1.0 / sqrt(nrm2))

    s = sample!(psi)

    @test length(s) == N
    for n in 1:N
      @test 1 <= s[n] <= 3
    end

    # Throws becase not orthogonalized to site 1:
    orthogonalize!(psi, 3)
    @test_throws ErrorException sample(psi)

    # Throws becase not normalized
    orthogonalize!(psi, 1)
    psi[1] *= (5.0 / norm(psi[1]))
    @test_throws ErrorException sample(psi)

    # Works when ortho & normalized:
    orthogonalize!(psi, 1)
    psi[1] *= (1.0 / norm(psi[1]))
    s = sample(psi)
    @test length(s) == N
  end

  @testset "randomMPS with chi > 1" begin
    N = 20
    chi = 8
    sites = siteinds(2, N)
    M = randomMPS(sites; linkdims=chi)

    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    @test norm(M[1]) ≈ 1.0

    @test maxlinkdim(M) == chi

    # Test for right-orthogonality
    R = M[N] * prime(M[N], "Link")
    r = linkind(M, N - 1)
    @test norm(R - delta(r, r')) < 1E-10
    for j in reverse(2:(N - 1))
      R = R * M[j] * prime(M[j], "Link")
      r = linkind(M, j - 1)
      @test norm(R - delta(r, r')) < 1E-10
    end

    # Complex case
    Mc = randomMPS(sites; linkdims=chi)
    @test inner(Mc, Mc) ≈ 1.0 + 0.0im
  end

  @testset "randomMPS from initial state (QN case)" begin
    N = 20
    chi = 8
    sites = siteinds("S=1/2", N; conserve_qns=true)

    # Make flux-zero random MPS
    state = [isodd(n) ? 1 : 2 for n in 1:N]
    M = randomMPS(sites, state; linkdims=chi)
    @test flux(M) == QN("Sz", 0)

    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    @test norm(M[1]) ≈ 1.0
    @test inner(M, M) ≈ 1.0

    @test maxlinkdim(M) == chi

    # Test making random MPS with different flux
    state[1] = 2
    M = randomMPS(sites, state; linkdims=chi)
    @test flux(M) == QN("Sz", -2)
    state[3] = 2
    M = randomMPS(sites, state; linkdims=chi)
    @test flux(M) == QN("Sz", -4)
  end

  @testset "Expected value and Correlations" begin
    N = 8
    m = 4

    # Non-fermionic case - spin system
    s = siteinds("S=1/2", N; conserve_qns=true)
    psi = randomMPS(s, n -> isodd(n) ? "Up" : "Dn"; linkdims=m)
    Cpm = correlation_matrix(psi, "S+", "S-")
    # Check using OpSum:
    for i in 1:N, j in i:N
      a = OpSum()
      a += "S+", i, "S-", j
      @test inner(psi, MPO(a, s), psi) ≈ Cpm[i, j]
    end
    PM = expect(psi, "S+*S-")
    @test norm(PM - diag(Cpm)) < 1E-8

    range = 3:7
    @test norm(PM[range] - expect(psi, "S+*S-"; site_range=range)) < 1E-8

    # With start_site, end_site arguments:
    s = siteinds("S=1/2", N)
    psi = randomMPS(ComplexF64, s; linkdims=m)
    ss, es = 3, 6
    Nb = es - ss + 1
    Cpm = correlation_matrix(psi, "S+", "S-"; site_range=ss:es)
    Czz = correlation_matrix(psi, "Sz", "Sz"; site_range=ss:es)
    @test size(Cpm) == (Nb, Nb)
    # Check using OpSum:
    for i in ss:es, j in i:es
      a = OpSum()
      a += "S+", i, "S-", j
      @test inner(psi, MPO(a, s), psi) ≈ Cpm[i - ss + 1, j - ss + 1]
    end

    # Fermionic case
    s = siteinds("Electron", N)
    psi = randomMPS(s; linkdims=m)
    Cuu = correlation_matrix(psi, "Cdagup", "Cup")
    # Check using OpSum:
    for i in 1:N, j in i:N
      a = OpSum()
      a += "Cdagup", i, "Cup", j
      @test inner(psi, MPO(a, s), psi) ≈ Cuu[i, j]
    end
  end

  @testset "expect regression test for in-place modification of input MPS" begin
    s = siteinds("S=1/2", 5)
    psi = randomMPS(s; linkdims=3)
    orthogonalize!(psi, 1)
    expect_init = expect(psi, "Sz")
    norm_scale = 10
    psi[1] *= norm_scale
    @test ortho_lims(psi) == 1:1
    @test norm(psi) ≈ norm_scale
    expect_Sz = expect(psi, "Sz")
    @test all(≤(1 / 2), expect_Sz)
    @test expect_Sz ≈ expect_init
    @test ortho_lims(psi) == 1:1
    @test norm(psi) ≈ norm_scale
  end

  @testset "correlation_matrix regression test for in-place modification of input MPS" begin
    s = siteinds("S=1/2", 5)
    psi = randomMPS(s; linkdims=3)
    orthogonalize!(psi, 1)
    correlation_matrix_init = correlation_matrix(psi, "Sz", "Sz")
    norm_scale = 10
    psi[1] *= norm_scale
    @test ortho_lims(psi) == 1:1
    @test norm(psi) ≈ norm_scale
    correlation_matrix_SzSz = correlation_matrix(psi, "Sz", "Sz")
    @test all(≤(1 / 2), correlation_matrix_SzSz)
    @test correlation_matrix_SzSz ≈ correlation_matrix_init
    @test ortho_lims(psi) == 1:1
    @test norm(psi) ≈ norm_scale
  end

  @testset "swapbondsites" begin
    N = 5
    sites = siteinds("S=1/2", N)
    ψ0 = randomMPS(sites)
    ψ = replacebond(ψ0, 3, ψ0[3] * ψ0[4]; swapsites=true, cutoff=1e-15)
    @test siteind(ψ, 1) == siteind(ψ0, 1)
    @test siteind(ψ, 2) == siteind(ψ0, 2)
    @test siteind(ψ, 4) == siteind(ψ0, 3)
    @test siteind(ψ, 3) == siteind(ψ0, 4)
    @test siteind(ψ, 5) == siteind(ψ0, 5)
    @test prod(ψ) ≈ prod(ψ0)
    @test maxlinkdim(ψ) == 1

    ψ = swapbondsites(ψ0, 4; cutoff=1e-15)
    @test siteind(ψ, 1) == siteind(ψ0, 1)
    @test siteind(ψ, 2) == siteind(ψ0, 2)
    @test siteind(ψ, 3) == siteind(ψ0, 3)
    @test siteind(ψ, 5) == siteind(ψ0, 4)
    @test siteind(ψ, 4) == siteind(ψ0, 5)
    @test prod(ψ) ≈ prod(ψ0)
    @test maxlinkdim(ψ) == 1
  end

  @testset "map!" begin
    N = 5
    s = siteinds("S=½", N)
    M0 = MPS(s, "↑")

    # Test map! with limits getting set
    M = orthogonalize(M0, 1)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2
    map!(prime, M)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == N + 1

    # Test map! without limits getting set
    M = orthogonalize(M0, 1)
    map!(prime, M; set_limits=false)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    # Test prime! with limits getting set
    M = orthogonalize(M0, 1)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2
    prime!(M; set_limits=true)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == N + 1

    # Test prime! without limits getting set
    M = orthogonalize(M0, 1)
    prime!(M)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2
  end

  @testset "setindex!(::MPS, _, ::Colon)" begin
    N = 4
    s = siteinds("S=½", N)
    ψ = randomMPS(s)
    ϕ = MPS(s, "↑")
    orthogonalize!(ϕ, 1)
    ψ[:] = ϕ
    @test ITensors.orthocenter(ψ) == 1
    @test inner(ψ, ϕ) ≈ 1

    ψ = randomMPS(s)
    ϕ = MPS(s, "↑")
    orthogonalize!(ϕ, 1)
    ψ[:] = ITensors.data(ϕ)
    @test ITensors.leftlim(ψ) == 0
    @test ITensors.rightlim(ψ) == N + 1
    @test inner(ψ, ϕ) ≈ 1
  end

  @testset "findsite[s](::MPS/MPO, is)" begin
    s = siteinds("S=1/2", 5)
    ψ = randomMPS(s)
    l = linkinds(ψ)

    A = randomITensor(s[4]', s[2]', dag(s[4]), dag(s[2]))

    @test findsite(ψ, s[3]) == 3
    @test findsite(ψ, (s[3], s[5])) == 3
    @test findsite(ψ, l[2]) == 2
    @test findsite(ψ, A) == 2

    @test findsites(ψ, s[3]) == [3]
    @test findsites(ψ, (s[4], s[1])) == [1, 4]
    @test findsites(ψ, l[2]) == [2, 3]
    @test findsites(ψ, (l[2], l[3])) == [2, 3, 4]
    @test findsites(ψ, A) == [2, 4]

    M = randomMPO(s)
    lM = linkinds(M)

    @test findsite(M, s[4]) == 4
    @test findsite(M, s[4]') == 4
    @test findsite(M, (s[4]', s[4])) == 4
    @test findsite(M, (s[4]', s[3])) == 3
    @test findsite(M, lM[2]) == 2
    @test findsite(M, A) == 2

    @test findsites(M, s[4]) == [4]
    @test findsites(M, s[4]') == [4]
    @test findsites(M, (s[4]', s[4])) == [4]
    @test findsites(M, (s[4]', s[3])) == [3, 4]
    @test findsites(M, (lM[2], lM[3])) == [2, 3, 4]
    @test findsites(M, A) == [2, 4]
  end

  @testset "[first]siteind[s](::MPS/MPO, j::Int)" begin
    s = siteinds("S=1/2", 5)
    ψ = randomMPS(s)
    @test siteind(first, ψ, 3) == s[3]
    @test siteind(ψ, 4) == s[4]
    @test isnothing(siteind(ψ, 4; plev=1))
    @test siteinds(ψ, 3) == IndexSet(s[3])
    @test siteinds(ψ, 3; plev=1) == IndexSet()

    M = randomMPO(s)
    @test noprime(siteind(first, M, 4)) == s[4]
    @test siteind(first, M, 4; plev=0) == s[4]
    @test siteind(first, M, 4; plev=1) == s[4]'
    @test siteind(M, 4) == s[4]
    @test siteind(M, 4; plev=0) == s[4]
    @test siteind(M, 4; plev=1) == s[4]'
    @test isnothing(siteind(M, 4; plev=2))
    @test hassameinds(siteinds(M, 3), (s[3], s[3]'))
    @test siteinds(M, 3; plev=1) == IndexSet(s[3]')
    @test siteinds(M, 3; plev=0) == IndexSet(s[3])
    @test siteinds(M, 3; tags="n=2") == IndexSet()
  end

  @testset "movesites $N sites" for N in 1:4
    s0 = siteinds("S=1/2", N)
    for perm in permutations(1:N)
      s = s0[perm]
      ψ = MPS(s, rand(("↑", "↓"), N))
      ns′ = [findfirst(==(i), s0) for i in s]
      @test ns′ == perm
      ψ′ = movesites(ψ, 1:N .=> ns′; cutoff=1e-15)
      if N == 1
        @test maxlinkdim(ψ′) == 1
      else
        @test maxlinkdim(ψ′) == 1
      end
      for n in 1:N
        @test s0[n] == siteind(ψ′, n)
      end
      @test prod(ψ) ≈ prod(ψ′)
    end
  end

  @testset "Construct MPS from ITensor" begin
    N = 5
    s = siteinds("S=1/2", N)
    l = [Index(3, "left_$n") for n in 1:2]
    r = [Index(3, "right_$n") for n in 1:2]

    #
    # MPS
    #

    A = randomITensor(s...)
    ψ = MPS(A, s)
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == N
    @test maxlinkdim(ψ) == 4

    ψ0 = MPS(s, "↑")
    A = prod(ψ0)
    ψ = MPS(A, s; cutoff=1e-15)
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == N
    @test maxlinkdim(ψ) == 1

    ψ0 = randomMPS(s; linkdims=2)
    A = prod(ψ0)
    ψ = MPS(A, s; cutoff=1e-15, orthocenter=2)
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == 2
    @test maxlinkdim(ψ) == 2

    A = randomITensor(s..., l[1], r[1])
    ψ = MPS(A, s; leftinds=l[1], orthocenter=3)
    ls = linkinds(ψ)
    @test hassameinds(ψ[1], (l[1], s[1], ls[1]))
    @test hassameinds(ψ[N], (r[1], s[N], ls[N - 1]))
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == 3
    @test maxlinkdim(ψ) == 12

    A = randomITensor(s..., l..., r...)
    ψ = MPS(A, s; leftinds=l)
    ls = linkinds(ψ)
    @test hassameinds(ψ[1], (l..., s[1], ls[1]))
    @test hassameinds(ψ[N], (r..., s[N], ls[N - 1]))
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == N
    @test maxlinkdim(ψ) == 36

    #
    # Construct from regular Julia tensor
    #
    Ajt = randn(2, 2, 2, 2, 2) # Julia tensor
    ψ = MPS(Ajt, s)
    @test prod(ψ) ≈ ITensor(Ajt, s...)

    Ajv = randn(2^N) # Julia vector
    ψ = MPS(Ajv, s)
    @test prod(ψ) ≈ ITensor(Ajv, s...)
  end

  @testset "Set range of MPS tensors" begin
    N = 5
    s = siteinds("S=1/2", N)
    ψ0 = randomMPS(s; linkdims=3)

    ψ = orthogonalize(ψ0, 2)
    A = prod(ITensors.data(ψ)[2:(N - 1)])
    randn!(A)
    ϕ = MPS(A, s[2:(N - 1)]; orthocenter=1)
    ψ[2:(N - 1)] = ϕ
    @test prod(ψ) ≈ ψ[1] * A * ψ[N]
    @test maxlinkdim(ψ) == 4
    @test ITensors.orthocenter(ψ) == 2

    ψ = orthogonalize(ψ0, 1)
    A = prod(ITensors.data(ψ)[2:(N - 1)])
    randn!(A)
    @test_throws AssertionError ψ[2:(N - 1)] = A

    ψ = orthogonalize(ψ0, 2)
    A = prod(ITensors.data(ψ)[2:(N - 1)])
    randn!(A)
    ψ[2:(N - 1), orthocenter=3] = A
    @test prod(ψ) ≈ ψ[1] * A * ψ[N]
    @test maxlinkdim(ψ) == 4
    @test ITensors.orthocenter(ψ) == 3
  end

  @testset "movesites reverse sites" begin
    N = 6
    s = siteinds("S=1/2", N)
    ψ0 = randomMPS(s)
    ψ = movesites(ψ0, 1:N .=> reverse(1:N))
    for n in 1:N
      @test siteind(ψ, n) == s[N - n + 1]
    end
  end

  @testset "movesites subsets of sites" begin
    N = 6
    s = siteinds("S=1/2", N)
    ψ = randomMPS(s)

    for i in 1:N, j in 1:N
      ns = [i, j]
      !allunique(ns) && continue
      min_ns = minimum(ns)
      ns′ = collect(min_ns:(min_ns + length(ns) - 1))
      ψ′ = movesites(ψ, ns .=> ns′; cutoff=1e-15)
      @test siteind(ψ′, min_ns) == siteind(ψ, i)
      @test siteind(ψ′, min_ns + 1) == siteind(ψ, j)
      @test maxlinkdim(ψ′) == 1
      ψ̃ = movesites(ψ′, ns′ .=> ns; cutoff=1e-15)
      for n in 1:N
        @test siteind(ψ̃, n) == siteind(ψ, n)
      end
      @test maxlinkdim(ψ̃) == 1
    end

    for i in 1:N, j in 1:N, k in 1:N
      ns = [i, j, k]
      !allunique(ns) && continue
      min_ns = minimum(ns)
      ns′ = collect(min_ns:(min_ns + length(ns) - 1))
      ψ′ = movesites(ψ, ns .=> ns′; cutoff=1e-15)
      @test siteind(ψ′, min_ns) == siteind(ψ, i)
      @test siteind(ψ′, min_ns + 1) == siteind(ψ, j)
      @test siteind(ψ′, min_ns + 2) == siteind(ψ, k)
      @test maxlinkdim(ψ′) == 1
      ψ̃ = movesites(ψ′, ns′ .=> ns; cutoff=1e-15)
      for n in 1:N
        @test siteind(ψ̃, n) == siteind(ψ, n)
      end
      @test maxlinkdim(ψ̃) == 1
    end

    #for i in 1:N, j in 1:N, k in 1:N, l in 1:N
    #  ns = [i, j, k, l]
    #  !allunique(ns) && continue
    #  min_ns = minimum(ns)
    #  ns′ = collect(min_ns:min_ns+length(ns)-1)
    #  ψ′ = movesites(ψ, ns .=> ns′; cutoff = 1e-15)
    #  @test siteind(ψ′, min_ns) == siteind(ψ, i)
    #  @test siteind(ψ′, min_ns+1) == siteind(ψ, j)
    #  @test siteind(ψ′, min_ns+2) == siteind(ψ, k)
    #  @test siteind(ψ′, min_ns+3) == siteind(ψ, l)
    #  @test maxlinkdim(ψ′) == 1
    #  ψ̃ = movesites(ψ′, ns′ .=> ns; cutoff = 1e-15)
    #  for n in 1:N
    #    @test siteind(ψ̃, n) == siteind(ψ, n)
    #  end
    #  @test maxlinkdim(ψ̃) == 1
    #end

    #for i in 1:N, j in 1:N, k in 1:N, l in 1:N, m in 1:N
    #  ns = [i, j, k, l, m]
    #  !allunique(ns) && continue
    #  min_ns = minimum(ns)
    #  ns′ = collect(min_ns:min_ns+length(ns)-1)
    #  ψ′ = movesites(ψ, ns .=> ns′; cutoff = 1e-15)
    #  for n in 1:length(ns)
    #    @test siteind(ψ′, min_ns+n-1) == siteind(ψ, ns[n])
    #  end
    #  @test maxlinkdim(ψ′) == 1
    #  ψ̃ = movesites(ψ′, ns′ .=> ns; cutoff = 1e-15)
    #  for n in 1:N
    #    @test siteind(ψ̃, n) == siteind(ψ, n)
    #  end
    #  @test maxlinkdim(ψ̃) == 1
    #end

  end

  @testset "product(::Vector{ITensor}, ::MPS)" begin
    N = 6
    s = siteinds("Qubit", N)

    I = [op("Id", s, n) for n in 1:N]
    X = [op("X", s, n) for n in 1:N]
    Y = [op("Y", s, n) for n in 1:N]
    Z = [op("Z", s, n) for n in 1:N]
    H = [op("H", s, n) for n in 1:N]
    CX = [op("CX", s, n, m) for n in 1:N, m in 1:N]
    CY = [op("CY", s, n, m) for n in 1:N, m in 1:N]
    CZ = [op("CZ", s, n, m) for n in 1:N, m in 1:N]
    CCNOT = [op("CCNOT", s, n, m, k) for n in 1:N, m in 1:N, k in 1:N]
    CSWAP = [op("CSWAP", s, n, m, k) for n in 1:N, m in 1:N, k in 1:N]
    CCCNOT = [op("CCCNOT", s, n, m, k, l) for n in 1:N, m in 1:N, k in 1:N, l in 1:N]

    v0 = [onehot(s[n] => "0") for n in 1:N]
    v1 = [onehot(s[n] => "1") for n in 1:N]

    # Single qubit
    @test product(I[1], v0[1]) ≈ v0[1]
    @test product(I[1], v1[1]) ≈ v1[1]

    @test product(H[1], H[1]) ≈ I[1]
    @test product(H[1], v0[1]) ≈ 1 / sqrt(2) * (v0[1] + v1[1])
    @test product(H[1], v1[1]) ≈ 1 / sqrt(2) * (v0[1] - v1[1])

    @test product(X[1], v0[1]) ≈ v1[1]
    @test product(X[1], v1[1]) ≈ v0[1]

    @test product(Y[1], v0[1]) ≈ im * v1[1]
    @test product(Y[1], v1[1]) ≈ -im * v0[1]

    @test product(Z[1], v0[1]) ≈ v0[1]
    @test product(Z[1], v1[1]) ≈ -v1[1]

    @test product(X[1], X[1]) ≈ I[1]
    @test product(Y[1], Y[1]) ≈ I[1]
    @test product(Z[1], Z[1]) ≈ I[1]
    @test -im * product([Y[1], X[1]], Z[1]) ≈ I[1]

    @test dag(X[1]) ≈ -product([X[1], Y[1]], Y[1])
    @test dag(Y[1]) ≈ -product([Y[1], Y[1]], Y[1])
    @test dag(Z[1]) ≈ -product([Z[1], Y[1]], Y[1])

    @test product(X[1], Y[1]) - product(Y[1], X[1]) ≈ 2 * im * Z[1]
    @test product(Y[1], Z[1]) - product(Z[1], Y[1]) ≈ 2 * im * X[1]
    @test product(Z[1], X[1]) - product(X[1], Z[1]) ≈ 2 * im * Y[1]

    @test product([Y[1], X[1]], v0[1]) - product([X[1], Y[1]], v0[1]) ≈
          2 * im * product(Z[1], v0[1])
    @test product([Y[1], X[1]], v1[1]) - product([X[1], Y[1]], v1[1]) ≈
          2 * im * product(Z[1], v1[1])
    @test product([Z[1], Y[1]], v0[1]) - product([Y[1], Z[1]], v0[1]) ≈
          2 * im * product(X[1], v0[1])
    @test product([Z[1], Y[1]], v1[1]) - product([Y[1], Z[1]], v1[1]) ≈
          2 * im * product(X[1], v1[1])
    @test product([X[1], Z[1]], v0[1]) - product([Z[1], X[1]], v0[1]) ≈
          2 * im * product(Y[1], v0[1])
    @test product([X[1], Z[1]], v1[1]) - product([Z[1], X[1]], v1[1]) ≈
          2 * im * product(Y[1], v1[1])

    #
    # 2-qubit
    #

    @test product(I[1] * I[2], v0[1] * v0[2]) ≈ v0[1] * v0[2]

    @test product(CX[1, 2], v0[1] * v0[2]) ≈ v0[1] * v0[2]
    @test product(CX[1, 2], v0[1] * v1[2]) ≈ v0[1] * v1[2]
    @test product(CX[1, 2], v1[1] * v0[2]) ≈ v1[1] * v1[2]
    @test product(CX[1, 2], v1[1] * v1[2]) ≈ v1[1] * v0[2]

    @test product(CY[1, 2], v0[1] * v0[2]) ≈ v0[1] * v0[2]
    @test product(CY[1, 2], v0[1] * v1[2]) ≈ v0[1] * v1[2]
    @test product(CY[1, 2], v1[1] * v0[2]) ≈ im * v1[1] * v1[2]
    @test product(CY[1, 2], v1[1] * v1[2]) ≈ -im * v1[1] * v0[2]

    @test product(CZ[1, 2], v0[1] * v0[2]) ≈ v0[1] * v0[2]
    @test product(CZ[1, 2], v0[1] * v1[2]) ≈ v0[1] * v1[2]
    @test product(CZ[1, 2], v1[1] * v0[2]) ≈ v1[1] * v0[2]
    @test product(CZ[1, 2], v1[1] * v1[2]) ≈ -v1[1] * v1[2]

    #
    # 3-qubit
    #

    @test product(CCNOT[1, 2, 3], v0[1] * v0[2] * v0[3]) ≈ v0[1] * v0[2] * v0[3]
    @test product(CCNOT[1, 2, 3], v0[1] * v0[2] * v1[3]) ≈ v0[1] * v0[2] * v1[3]
    @test product(CCNOT[1, 2, 3], v0[1] * v1[2] * v0[3]) ≈ v0[1] * v1[2] * v0[3]
    @test product(CCNOT[1, 2, 3], v0[1] * v1[2] * v1[3]) ≈ v0[1] * v1[2] * v1[3]
    @test product(CCNOT[1, 2, 3], v1[1] * v0[2] * v0[3]) ≈ v1[1] * v0[2] * v0[3]
    @test product(CCNOT[1, 2, 3], v1[1] * v0[2] * v1[3]) ≈ v1[1] * v0[2] * v1[3]
    @test product(CCNOT[1, 2, 3], v1[1] * v1[2] * v0[3]) ≈ v1[1] * v1[2] * v1[3]
    @test product(CCNOT[1, 2, 3], v1[1] * v1[2] * v1[3]) ≈ v1[1] * v1[2] * v0[3]

    @test product(CSWAP[1, 2, 3], v0[1] * v0[2] * v0[3]) ≈ v0[1] * v0[2] * v0[3]
    @test product(CSWAP[1, 2, 3], v0[1] * v0[2] * v1[3]) ≈ v0[1] * v0[2] * v1[3]
    @test product(CSWAP[1, 2, 3], v0[1] * v1[2] * v0[3]) ≈ v0[1] * v1[2] * v0[3]
    @test product(CSWAP[1, 2, 3], v0[1] * v1[2] * v1[3]) ≈ v0[1] * v1[2] * v1[3]
    @test product(CSWAP[1, 2, 3], v1[1] * v0[2] * v0[3]) ≈ v1[1] * v0[2] * v0[3]
    @test product(CSWAP[1, 2, 3], v1[1] * v0[2] * v1[3]) ≈ v1[1] * v1[2] * v0[3]
    @test product(CSWAP[1, 2, 3], v1[1] * v1[2] * v0[3]) ≈ v1[1] * v0[2] * v1[3]
    @test product(CSWAP[1, 2, 3], v1[1] * v1[2] * v1[3]) ≈ v1[1] * v1[2] * v1[3]

    #
    # Apply to an MPS
    #

    ψ = MPS(s, "0")
    @test prod(product(X[1], ψ)) ≈ prod(MPS(s, n -> n == 1 ? "1" : "0"))
    @test prod(product(X[1], product(X[2], ψ))) ≈
          prod(MPS(s, n -> n == 1 || n == 2 ? "1" : "0"))
    @test prod(product(X[1] * X[2], ψ)) ≈ prod(MPS(s, n -> n == 1 || n == 2 ? "1" : "0"))
    @test prod(product([X[2], X[1]], ψ)) ≈ prod(MPS(s, n -> n == 1 || n == 2 ? "1" : "0"))
    @test prod(product(CX[1, 2], ψ)) ≈ prod(MPS(s, "0"))
    @test prod(product(CX[1, 2], product(X[1], ψ))) ≈
          prod(MPS(s, n -> n == 1 || n == 2 ? "1" : "0"))
    @test prod(product(product(CX[1, 2], X[1]), ψ)) ≈
          prod(MPS(s, n -> n == 1 || n == 2 ? "1" : "0"))
    @test prod(product([X[1], CX[1, 2]], ψ)) ≈
          prod(MPS(s, n -> n == 1 || n == 2 ? "1" : "0"))

    for i in 1:N, j in 1:N
      !allunique((i, j)) && continue
      # Don't move sites back
      CXij_ψ = product([X[i], CX[i, j]], ψ; move_sites_back=false, cutoff=1e-15)
      @test maxlinkdim(CXij_ψ) == 1
      @test prod(CXij_ψ) ≈ prod(MPS(s, n -> n == i || n == j ? "1" : "0"))

      # Move sites back
      CXij_ψ = product([X[i], CX[i, j]], ψ)
      for n in 1:N
        @test siteind(CXij_ψ, n) == siteind(ψ, n)
      end
      @test prod(CXij_ψ) ≈ prod(MPS(s, n -> n == i || n == j ? "1" : "0"))
    end

    for i in 1:N, j in 1:N, k in 1:N
      ns = (i, j, k)
      !allunique(ns) && continue
      # Don't move sites back
      CCNOTijk_ψ = product(
        [X[j], X[i], CCNOT[ns...]], ψ; move_sites_back=false, cutoff=1e-15
      )
      @test maxlinkdim(CCNOTijk_ψ) == 1
      @test prod(CCNOTijk_ψ) ≈ prod(MPS(s, n -> n ∈ ns ? "1" : "0"))

      # Move sites back
      CCNOTijk_ψ = product([X[j], X[i], CCNOT[ns...]], ψ; cutoff=1e-15)
      @test maxlinkdim(CCNOTijk_ψ) == 1
      for n in 1:N
        @test siteind(CCNOTijk_ψ, n) == siteind(ψ, n)
      end
      @test prod(CCNOTijk_ψ) ≈ prod(MPS(s, n -> n ∈ ns ? "1" : "0"))
    end

    for i in 1:N, j in i:N, k in 1:N, l in k:N
      ns = (i, j, k, l)
      !allunique(ns) && continue
      # Don't move sites back
      CCCNOTijkl_ψ = product(
        [X[i], X[j], X[k], CCCNOT[ns...]], ψ; move_sites_back=false, cutoff=1e-15
      )
      @test maxlinkdim(CCCNOTijkl_ψ) == 1
      @test prod(CCCNOTijkl_ψ) ≈ prod(MPS(s, n -> n ∈ ns ? "1" : "0"))

      # Move sites back
      CCCNOTijkl_ψ = product([X[i], X[j], X[k], CCCNOT[ns...]], ψ; cutoff=1e-15)
      @test maxlinkdim(CCCNOTijkl_ψ) == 1
      for n in 1:N
        @test siteind(CCCNOTijkl_ψ, n) == siteind(ψ, n)
      end
      @test prod(CCCNOTijkl_ψ) ≈ prod(MPS(s, n -> n ∈ ns ? "1" : "0"))
    end
  end

  @testset "product" begin
    @testset "Contraction order of operations" begin
      s = siteind("Qubit")
      Q = SiteType("Qubit")
      @test product(ops([s], [("Y", 1), ("X", 1)]), setelt(s => 1)) ≈
            itensor(op("X", Q) * op("Y", Q) * [1; 0], s)
      @test product(ops([s], [("Y", 1), ("Z", 1)]), setelt(s => 1)) ≈
            itensor(op("Z", Q) * op("Y", Q) * [1; 0], s)
      @test product(ops([s], [("X", 1), ("Y", 1)]), setelt(s => 1)) ≈
            itensor(op("Y", Q) * op("X", Q) * [1; 0], s)
    end

    @testset "Simple on-site state evolution" begin
      N = 3

      pos = [("Z", 3), ("Y", 2), ("X", 1)]

      s = siteinds("Qubit", N)
      gates = ops(s, pos)
      ψ0 = MPS(s, "0")

      # Apply the gates
      ψ = product(gates, ψ0)

      # Move site 1 to position 3
      ψ′ = movesite(ψ, 1 => 3)
      @test siteind(ψ′, 1) == s[2]
      @test siteind(ψ′, 2) == s[3]
      @test siteind(ψ′, 3) == s[1]
      @test prod(ψ) ≈ prod(ψ′)

      # Move the site back
      ψ′′ = movesite(ψ′, 3 => 1)
      @test siteind(ψ′′, 1) == s[1]
      @test siteind(ψ′′, 2) == s[2]
      @test siteind(ψ′′, 3) == s[3]
      @test prod(ψ) ≈ prod(ψ′′)
    end

    @testset "More complex evolution" begin
      N = 7

      osX = [("X", n) for n in 1:N]

      osZ = [("Z", n) for n in 1:N]

      osSw = [("SWAP", n, n + 1) for n in 1:(N - 2)]

      osCx = [("CX", n, n + 3) for n in 1:(N - 3)]

      osT = [("CCX", n, n + 1, n + 3) for n in 1:(N - 3)]

      osRx = [("Rx", n, (θ=π,)) for n in 1:N]

      osXX = [("XX", (n, n + 1), (ϕ=π / 8,)) for n in 1:(N - 1)]

      #os_noise = [("noise", n, n+2, n+4) for n in 1:N-4]

      os = vcat(osX, osXX, osSw, osRx, osZ, osCx, osT)
      s = siteinds("Qubit", N)
      gates = ops(os, s)

      @testset "Pure state evolution" begin
        ψ0 = MPS(s, "0")
        ψ = product(gates, ψ0; cutoff=1e-15)
        @test maxlinkdim(ψ) == 8
        prodψ = product(gates, prod(ψ0))
        @test prod(ψ) ≈ prodψ rtol = 1e-12
      end

      M0 = MPO(s, "Id")
      maxdim = prod(dim(siteinds(M0, j)) for j in 1:N)

      @testset "Mixed state evolution" begin
        M = product(gates, M0; cutoff=1e-15, maxdim=maxdim)
        @test maxlinkdim(M) == 24 || maxlinkdim(M) == 25
        sM0 = siteinds(M0)
        sM = siteinds(M)
        for n in 1:N
          @test hassameinds(sM[n], sM0[n])
        end
        @set_warn_order 15 begin
          prodM = product(gates, prod(M0))
          @test prod(M) ≈ prodM rtol = 1e-6
        end
      end

      #@testset "Mixed state noisy evolution" begin
      #  prepend!(os, os_noise)
      #  gates = ops(os, s)
      #  M = product(gates, M0; apply_dag = true,
      #              cutoff = 1e-15, maxdim = maxdim)
      #  @test maxlinkdim(M) == 64
      #  sM0 = siteinds(M0)
      #  sM = siteinds(M)
      #  for n in 1:N
      #    @test hassameinds(sM[n], sM0[n])
      #  end
      #  @set_warn_order 16 begin
      #    prodM = product(gates, prod(M0); apply_dag = true)
      #    @test prod(M) ≈ prodM rtol = 1e-7
      #  end
      #end

      #@testset "Mixed state noisy evolution" begin
      #  prepend!(os, os_noise)
      #  gates = ops(os, s)
      #  M = product(gates, M0;
      #              apply_dag = true, cutoff = 1e-15, maxdim = maxdim-1)
      #  @test maxlinkdim(M) == 64
      #  sM0 = siteinds(M0)
      #  sM = siteinds(M)
      #  for n in 1:N
      #    @test hassameinds(sM[n], sM0[n])
      #  end
      #  @set_warn_order 16 begin
      #    prodM = product(gates, prod(M0); apply_dag = true)
      #    @test prod(M) ≈ prodM rtol = 1e-7
      #  end
      #end

    end

    @testset "Gate evolution open system" begin
      N = 8
      osX = [("X", n) for n in 1:N]
      osZ = [("Z", n) for n in 1:N]
      osSw = [("SWAP", n, n + 2) for n in 1:(N - 2)]
      osCx = [("CX", n, n + 3) for n in 1:(N - 3)]
      osT = [("CCX", n, n + 1, n + 3) for n in 1:(N - 3)]
      osRx = [("Rx", n, (θ=π,)) for n in 1:N]
      os = vcat(osX, osSw, osRx, osZ, osCx, osT)

      s = siteinds("Qubit", N)
      gates = ops(os, s)

      M0 = MPO(s, "Id")

      # Apply the gates

      s0 = siteinds(M0)

      M = apply(gates, M0; apply_dag=true, cutoff=1e-15, maxdim=500, svd_alg="qr_iteration")

      s = siteinds(M)
      for n in 1:N
        @assert hassameinds(s[n], s0[n])
      end

      @set_warn_order 18 begin
        prodM = apply(gates, prod(M0); apply_dag=true)
        @test prod(M) ≈ prodM rtol = 1e-6
      end
    end

    @testset "Gate evolution state" begin
      N = 10

      osX = [("X", n) for n in 1:N]
      osZ = [("Z", n) for n in 1:N]
      osSw = [("SWAP", n, n + 1) for n in 1:(N - 1)]
      osCx = [("CX", n, n + 1) for n in 1:(N - 1)]
      osT = [("CCX", n, n + 2, n + 4) for n in 1:(N - 4)]
      os = vcat(osX, osSw, osZ, osCx, osT)

      s = siteinds("Qubit", N)
      gates = ops(os, s)

      ψ0 = MPS(s, "0")

      # Apply the gates
      ψ = apply(gates, ψ0; cutoff=1e-15, maxdim=100)

      prodψ = apply(gates, prod(ψ0))
      @test prod(ψ) ≈ prodψ rtol = 1e-4
    end

    @testset "With fermions" begin
      N = 3

      s = siteinds("Fermion", N; conserve_qns=true)

      # Ground state |000⟩
      ψ000 = MPS(s, "0")

      # Start state |011⟩
      ψ011 = MPS(s, n -> n == 2 || n == 3 ? "1" : "0")

      # Reference state |110⟩
      ψ110 = MPS(s, n -> n == 1 || n == 2 ? "1" : "0")

      function ITensors.op(::OpName"CdagC", ::SiteType, s1::Index, s2::Index)
        return op("Cdag", s1) * op("C", s2)
      end

      os = [("CdagC", 1, 3)]
      Os = ops(os, s)

      # Results in -|110⟩
      ψ1 = product(Os, ψ011; cutoff=1e-15)

      @test inner(ψ1, ψ110) == -1

      a = OpSum()
      a += "Cdag", 1, "C", 3
      H = MPO(a, s)

      # Results in -|110⟩
      ψ2 = noprime(contract(H, ψ011; cutoff=1e-15))

      @test inner(ψ2, ψ110) == -1
    end

    @testset "Spinless fermion (gate evolution)" begin
      N = 6

      s = siteinds("Fermion", N; conserve_qns=true)

      # Starting state
      ψ0 = MPS(s, n -> isodd(n) ? "0" : "1")

      t = 1.0
      U = 1.0
      ampo = OpSum()
      for b in 1:(N - 1)
        ampo .+= -t, "Cdag", b, "C", b + 1
        ampo .+= -t, "Cdag", b + 1, "C", b
        ampo .+= U, "N", b, "N", b + 1
      end
      H = MPO(ampo, s)

      sweeps = Sweeps(6)
      maxdim!(sweeps, 10, 20, 40)
      cutoff!(sweeps, 1E-12)
      energy, ψ0 = dmrg(H, ψ0, sweeps; outputlevel=0)

      function ITensors.op(::OpName"CdagC", ::SiteType, s1::Index, s2::Index)
        return op("Cdag", s1) * op("C", s2)
      end

      function ITensors.op(
        ::OpName"CCCC", ::SiteType, s1::Index, s2::Index, s3::Index, s4::Index
      )
        return -1 * op("Cdag", s1) * op("Cdag", s2) * op("C", s3) * op("C", s4)
      end

      for i in 1:(N - 1), j in (i + 1):N
        G1 = op("CdagC", s, i, j)

        @disable_warn_order begin
          G2 = op("Cdag", s, i)
          for n in (i + 1):(j - 1)
            G2 *= op("F", s, n)
          end
          G2 *= op("C", s, j)
        end

        ampo = OpSum()
        ampo += "Cdag", i, "C", j
        G3 = MPO(ampo, s)

        A_OP = prod(product(G1, ψ0; cutoff=1e-6))

        A_OPS = noprime(G2 * prod(ψ0))

        A_MPO = noprime(prod(contract(G3, ψ0; cutoff=1e-6)))

        @test A_OP ≈ A_OPS
        @test A_OP ≈ A_MPO
      end

      for i in 1:(N - 3), j in (i + 1):(N - 2), k in (j + 1):(N - 1), l in (k + 1):N
        G1 = op("CCCC", s, i, j, k, l)
        @disable_warn_order begin
          G2 = -1 * op("Cdag", s, i)
          for n in (i + 1):(j - 1)
            G2 *= op("F", s, n)
          end
          G2 *= op("Cdag", s, j)
          for n in (j + 1):(k - 1)
            G2 *= op("Id", s, n)
          end
          G2 *= op("C", s, k)
          for n in (k + 1):(l - 1)
            G2 *= op("F", s, n)
          end
          G2 *= op("C", s, l)

          ampo = OpSum()
          ampo += "Cdag", i, "Cdag", j, "C", k, "C", l
          G3 = MPO(ampo, s)

          A_OP = prod(product(G1, ψ0; cutoff=1e-16))

          A_OPS = noprime(G2 * prod(ψ0))

          A_MPO = noprime(prod(contract(G3, ψ0; cutoff=1e-16)))
        end
        @test A_OPS ≈ A_OP rtol = 1e-12
        @test A_MPO ≈ A_OP rtol = 1e-12
      end
    end

    @testset "Spinful Fermions (Electron) gate evolution" begin
      N = 8
      s = siteinds("Electron", N; conserve_qns=true)
      ψ0 = randomMPS(s, n -> isodd(n) ? "↑" : "↓")
      t = 1.0
      U = 1.0
      ampo = OpSum()
      for b in 1:(N - 1)
        ampo .+= -t, "Cdagup", b, "Cup", b + 1
        ampo .+= -t, "Cdagup", b + 1, "Cup", b
        ampo .+= -t, "Cdagdn", b, "Cdn", b + 1
        ampo .+= -t, "Cdagdn", b + 1, "Cdn", b
      end
      for n in 1:N
        ampo .+= U, "Nupdn", n
      end
      H = MPO(ampo, s)
      sweeps = Sweeps(6)
      maxdim!(sweeps, 10, 20, 40)
      cutoff!(sweeps, 1E-12)
      energy, ψ = dmrg(H, ψ0, sweeps; outputlevel=0)

      function ITensors.op(::OpName"CCup", ::SiteType"Electron", s1::Index, s2::Index)
        return op("Adagup*F", s1) * op("Aup", s2)
      end

      for i in 1:(N - 1), j in (i + 1):N
        ampo = OpSum()
        ampo += "Cdagup", i, "Cup", j
        G1 = MPO(ampo, s)
        G2 = op("CCup", s, i, j)
        A_MPO = prod(noprime(contract(G1, ψ; cutoff=1e-8)))
        A_OP = prod(product(G2, ψ; cutoff=1e-8))
        @test A_MPO ≈ A_OP atol = 1E-4
      end
    end
  end

  @testset "dense conversion of MPS" begin
    N = 4
    s = siteinds("S=1/2", N; conserve_qns=true)
    QM = randomMPS(s, ["Up", "Dn", "Up", "Dn"]; linkdims=4)
    qsz1 = scalar(QM[1] * op("Sz", s[1]) * dag(prime(QM[1], "Site")))

    M = dense(QM)
    @test !hasqns(M[1])
    sz1 = scalar(M[1] * op("Sz", removeqns(s[1])) * dag(prime(M[1], "Site")))
    @test sz1 ≈ qsz1
  end

  @testset "inner of MPS with more than one site Index" begin
    s = siteinds("S=½", 4)
    sout = addtags.(s, "out")
    sin = addtags.(s, "in")
    sinds = IndexSet.(sout, sin)
    Cs = combiner.(sinds)
    cinds = combinedind.(Cs)
    ψ = randomMPS(cinds)
    @test norm(ψ) ≈ 1
    @test inner(ψ, ψ) ≈ 1
    ψ .*= dag.(Cs)
    @test norm(ψ) ≈ 1
    @test inner(ψ, ψ) ≈ 1
  end

  @testset "inner(::MPS, ::MPO, ::MPS) with more than one site Index" begin
    N = 8
    s = siteinds("S=1/2", N)
    a = OpSum()
    for j in 1:(N - 1)
      a .+= 0.5, "S+", j, "S-", j + 1
      a .+= 0.5, "S-", j, "S+", j + 1
      a .+= "Sz", j, "Sz", j + 1
    end
    H = MPO(a, s)
    ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=10)
    # Create MPO/MPS with pairs of sites merged
    H2 = MPO([H[b] * H[b + 1] for b in 1:2:N])
    ψ2 = MPS([ψ[b] * ψ[b + 1] for b in 1:2:N])
    @test inner(ψ, ψ) ≈ inner(ψ2, ψ2)
    @test inner(ψ', H, ψ) ≈ inner(ψ2', H2, ψ2)
    @test_throws ErrorException inner(ψ2, ψ2')
    @test_throws ErrorException inner(ψ2, H2, ψ2)
  end

  @testset "orthogonalize! on MPS with no link indices" begin
    N = 4
    s = siteinds("S=1/2", N)
    ψ = MPS([itensor(randn(ComplexF64, 2), s[n]) for n in 1:N])
    @test all(==(IndexSet()), linkinds(all, ψ))
    ϕ = orthogonalize(ψ, 2)
    @test ITensors.hasdefaultlinktags(ϕ)
    @test ortho_lims(ϕ) == 2:2
    @test ITensors.dist(ψ, ϕ) ≈ 0 atol = 1e-6
    # TODO: use this instead?
    # @test lognorm(ψ - ϕ) < -16
    @test norm(ψ - ϕ) ≈ 0 atol = 1e-6
  end

  @testset "MPO from MPS with no link indices" begin
    N = 4
    s = siteinds("S=1/2", N)
    ψ = MPS([itensor(randn(ComplexF64, 2), s[n]) for n in 1:N])
    ρ = outer(ψ, ψ)
    @test !ITensors.hasnolinkinds(ρ)
    @test inner(ρ, ρ) ≈ inner(ψ, ψ)^2
    @test inner(ψ, ρ, ψ) ≈ inner(ψ, ψ)^2

    # Deprecated syntax
    ρ = MPO(ψ)
    @test !ITensors.hasnolinkinds(ρ)
    @test inner(ρ, ρ) ≈ inner(ψ, ψ)^2
    @test inner(ψ, ρ, ψ) ≈ inner(ψ, ψ)^2
  end

  @testset "Truncate MPO with no link indices" begin
    N = 4
    s = siteinds("S=1/2", N)
    M = MPO([itensor(randn(ComplexF64, 2, 2), s[n]', dag(s[n])) for n in 1:N])
    @test ITensors.hasnolinkinds(M)
    Mt = truncate(M; cutoff=1e-15)
    @test ITensors.hasdefaultlinktags(Mt)
    @test norm(M - Mt) ≈ 0 atol = 1e-12
  end
end

nothing