text:
# Newton solver options
@with_kw mutable struct NewtonOptions{T}
    r_tol::T = 1.0e-5            # primal dual residual tolerance
    max_iter::Int = 10           # primal dual iter
    β_init::T = 1.0e-5           # initial dual regularization
    live_plotting::Bool = false  # visualize the trajectory during the solve
    verbose::Bool = false
    solver::Symbol = :lu_solver
end

struct NewtonResidual{T,vq2,vu1,vd,vI,vq0,vq1}
    r#::Vector{T}                           # residual

    q2::Vector{vq2}                    # rsd objective views
    u1::Vector{vu1}                    # rsd objective views

    rd::Vector{vd}                         # rsd dynamics lagrange multiplier views
    rI::Vector{vI}                         # rsd dynamics -I views [q2]
    q0::Vector{vq0}                        # rsd dynamics q0 views
    q1::Vector{vq1}                        # rsd dynamics q1 views
end

function NewtonResidual(model::Model, env::Environment, H::Int)

    nq = model.nq # configuration
    nu = model.nu # control
    nc = model.nc # contact
    nd = nq    # implicit dynamics constraint
    nr = nq + nu + nd # size of a one-time-step block

    off = 0
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iν = SizedVector{nd}(off .+ (1:nd)); off += nd # index of the dynamics lagrange multiplier ν1
    iz = copy(iq) # index of the IP solver solution [q2]

    r = zeros(H * nr)

    q2  = [view(r, (t - 1) * nr .+ iq) for t = 1:H]
    u1  = [view(r, (t - 1) * nr .+ iu) for t = 1:H]

    rd  = [view(r, (t - 1) * nr .+ iν) for t = 1:H]
    rI  = [view(r, (t - 1) * nr .+ iz) for t = 1:H]
    q0  = [view(r, (t - 3) * nr .+ iq) for t = 3:H]
    q1  = [view(r, (t - 2) * nr .+ iq) for t = 2:H]

    T = eltype(r)

    return NewtonResidual{T, eltype.((q2, u1))...,eltype.((rd, rI, q0, q1))...}(
        r, q2, u1, rd, rI, q0, q1)
end

struct NewtonJacobian{T,Vq,Vu,VI,VIT,Vq0,Vq0T,Vq1,Vq1T,Vu1,Vu1T,Vreg}
    R#::SparseMatrixCSC{T,Int}                 # jacobian

    obj_q2::Vector{Vq}                          # obj views
    obj_u1::Vector{Vu}                          # obj views

    obj_q1q2::Vector{Vq}
    obj_q2q1::Vector{Vq}

    IV::Vector{VI}                          # dynamics -I views [q2]
    ITV::Vector{VIT}                        # dynamics -I views [q2] transposed
    q0::Vector{Vq0}                         # dynamics q0 views
    q0T::Vector{Vq0T}                       # dynamics q0 views transposed
    q1::Vector{Vq1}                         # dynamics q1 views
    q1T::Vector{Vq1T}                       # dynamics q1 views transposed
    u1::Vector{Vu1}                         # dynamics u1 views
    u1T::Vector{Vu1T}                       # dynamics u1 views transposed
    reg::Vector{Vreg}                       # dual regularization views
end

function NewtonJacobian(model::Model, env::Environment, H::Int)

    nq = model.nq # configuration
    nu = model.nu # control
    nw = model.nw # disturbance
    nc = model.nc # contact
    nd = nq # implicit dynamics constraint
    nr = nq + nu + nd # size of a one-time-step block

    off = 0
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iν = SizedVector{nd}(off .+ (1:nd)); off += nd # index of the dynamics lagrange multiplier ν1
    iz = copy(iq)                                  # index of the IP solver solution [q2]
    iθ = vcat(iq .- 2nr, iq .- nr, iu) # index of the IP solver data [q0, q1, u1]

    R = spzeros(H * nr, H * nr)

    obj_q2  = [view(R, (t - 1) * nr .+ iq, (t - 1) * nr .+ iq) for t = 1:H]
    obj_u1  = [view(R, (t - 1) * nr .+ iu, (t - 1) * nr .+ iu) for t = 1:H]

    obj_q1q2  = [view(R, (t - 1) * nr .+ iq, t * nr .+ iq) for t = 1:H-1]
    obj_q2q1  = [view(R, t * nr .+ iq, (t - 1) * nr .+ iq) for t = 1:H-1]

    IV  = [view(R, CartesianIndex.((t - 1) * nr .+ iz, (t - 1) * nr .+ iν)) for t = 1:H]
    ITV = [view(R, CartesianIndex.((t - 1) * nr .+ iν, (t - 1) * nr .+ iz)) for t = 1:H]
    q0  = [view(R, (t - 1) * nr .+ iν, (t - 3) * nr .+ iq) for t = 3:H]
    q0T = [view(R, (t - 3) * nr .+ iq, (t - 1) * nr .+ iν) for t = 3:H]
    q1  = [view(R, (t - 1) * nr .+ iν, (t - 2) * nr .+ iq) for t = 2:H]
    q1T = [view(R, (t - 2) * nr .+ iq, (t - 1) * nr .+ iν) for t = 2:H]
    u1  = [view(R, (t - 1) * nr .+ iν, (t - 1) * nr .+ iu) for t = 1:H]
    u1T = [view(R, (t - 1) * nr .+ iu, (t - 1) * nr .+ iν) for t = 1:H]
    reg = [view(R, CartesianIndex.((t - 1) * nr .+ iν, (t - 1) * nr .+ iν)) for t = 1:H] # TODO: Cartesian indices to only grab diagonals

    return NewtonJacobian{eltype(R),
        eltype.((obj_q2, obj_u1))...,
        eltype.((IV, ITV, q0, q0T, q1, q1T, u1, u1T, reg))...}(
        R, obj_q2, obj_u1,
        obj_q1q2, obj_q2q1,
        IV, ITV, q0, q0T, q1, q1T, u1, u1T, reg)
end

mutable struct NewtonIndices{nq,nu,n1,n2,n3}
    nd::Int                                   # implicit dynamics constraint
    nr::Int                                   # size of a one-time-step block
    iq::SizedArray{Tuple{nq},Int,1,1}         # configuration indices
    iu::SizedArray{Tuple{nu},Int,1,1}         # control indices
    iν::SizedArray{Tuple{n1},Int,1,1}         # implicit dynamics lagrange multiplier
    iz::SizedArray{Tuple{n2},Int,1,1}         # IP solver solution [q2]
    iθ::SizedArray{Tuple{n3},Int,1,1}         # IP solver data [q0, q1, u1]
    Iq::Vector{SizedArray{Tuple{nq},Int,1,1}} # configuration indices
    Iu::Vector{SizedArray{Tuple{nu},Int,1,1}} # control indices
    Iν::Vector{SizedArray{Tuple{n1},Int,1,1}} # implicit dynamics lagrange multiplier
    Iz::Vector{SizedArray{Tuple{n2},Int,1,1}} # IP solver solution [q2]
    Iθ::Vector{SizedArray{Tuple{n3},Int,1,1}} # IP solver data [q0, q1, u1]
end

function NewtonIndices(model::Model, env::Environment, H::Int)
    nq = model.nq # configuration
    nu = model.nu # control
    nw = model.nw # disturbance
    nc = model.nc # contact
    nd = nq # implicit dynamics constraint
    nr = nq + nu + nd # size of a one-time-step block

    off = 0
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iν = SizedVector{nd}(off .+ (1:nd)); off += nd # index of the dynamics lagrange multiplier ν1
    iz = copy(iq)                                  # index of the IP solver solution [q2]
    iθ = vcat(iq .- 2nr, iq .- nr, iu) # index of the IP solver data [q0, q1, u1]

    Iq = [(t - 1) * nr .+ iq for t = 1:H]
    Iu = [(t - 1) * nr .+ iu for t = 1:H]
    Iν = [(t - 1) * nr .+ iν for t = 1:H]
    Iz = [(t - 1) * nr .+ iz for t = 1:H]
    Iθ = [(t - 1) * nr .+ iθ for t = 1:H]

    return NewtonIndices{nq,nu,nd,nd,2nq+nu}(
        nd, nr,
        iq, iu, iν, iz, iθ,
        Iq, Iu, Iν, Iz, Iθ)
end

mutable struct Newton{T,nq,nu,nw,nz,nθ,n1,n2,n3}
    jac::NewtonJacobian{T}                          # NewtonJacobian
    res::NewtonResidual{T}                          # residual
    res_cand::NewtonResidual{T}                     # candidate residual
    Δ::NewtonResidual{T}                            # step direction in the Newton solve, it contains: q2-qH+1, u1-uH, γ1-γH, b1-bH, λd1-λdH
    ν::Vector{SizedArray{Tuple{n1},T,1,1}}          # implicit dynamics lagrange multiplier
    ν_cand::Vector{SizedArray{Tuple{n1},T,1,1}}         # candidate implicit dynamics lagrange multiplier
    traj::ContactTraj       # optimized trajectory
    traj_cand::ContactTraj # trial trajectory used in line search
    Δq::Vector{SizedArray{Tuple{nq},T,1,1}}         # difference between the traj and ref_traj
    Δu::Vector{SizedArray{Tuple{nu},T,1,1}}         # difference between the traj and ref_traj
    ind::NewtonIndices                                 # indices of a one-time-step block
    obj::Objective                              # obj function
    solver::LinearSolver
    β::T
    opts::NewtonOptions{T}                          # Newton solver options
end

function Newton(s::Simulation, H::Int, h::T,
    traj::ContactTraj, im_traj::ImplicitTrajectory;
    obj::Objective = TrackingObjective(s.model, s.env, H),
    opts::NewtonOptions = NewtonOptions(), κ::T=im_traj.ip[1].κ[1]) where T

    model = s.model
    env = s.env

    ind = NewtonIndices(model, env, H)

    nq = model.nq
    nu = model.nu
    nw = model.nw
    nz = num_var(model, env)
    nθ = num_data(model)
    nd = ind.nd

    jac = NewtonJacobian(model, env, H)

    # precompute Jacobian for pre-factorization
    window = collect(1:(H + 2))
    implicit_dynamics!(im_traj, traj, window=window) #@@@
    jacobian!(jac, im_traj, obj, H, opts.β_init, window)

    res = NewtonResidual(model, env, H)
    res_cand = NewtonResidual(model, env, H)

    Δ = NewtonResidual(model, env, H)

    ν = [zeros(SizedVector{ind.nd,T}) for t = 1:H]
    ν_cand = deepcopy(ν)

    traj = contact_trajectory(model, env, H, h, κ=κ)
    traj_cand = contact_trajectory(model, env, H, h, κ=κ)

    Δq  = [zeros(SizedVector{nq,T}) for t = 1:H]
    Δu  = [zeros(SizedVector{nu,T}) for t = 1:H]

    # regularization
    β = copy(opts.β_init)

    # linear solver
    solver = eval(opts.solver)(jac.R)

    return Newton{T,nq,nu,nw,nz,nθ,nd,nd,2nq+nu}(
        jac, res, res_cand, Δ, ν, ν_cand, traj, traj_cand,
        Δq, Δu, ind, obj, solver, β, opts)
end

function initialize_jacobian!(jac::NewtonJacobian, obj::Objective, H::Int; update_hessian::Bool=true)

    fill!(jac.R, 0.0)

    # Objective
    update_hessian && hessian!(jac, obj)

    for t = 1:H
        # Implicit dynamics
        jac.IV[t] .-= 1.0
        jac.ITV[t] .-= 1.0
    end

    return nothing
end

function update_jacobian!(jac::NewtonJacobian, im_traj::ImplicitTrajectory, obj::Objective,
    H::Int, β::T, window::Vector{Int}) where T

    # reset
    for t = 1:H
        if t >= 3
            fill!(jac.q0[t-2], 0.0)
            fill!(jac.q0T[t-2], 0.0)
        end

        if t >= 2
            fill!(jac.q1[t-1], 0.0)
            fill!(jac.q1T[t-1], 0.0)
        end

        fill!(jac.u1[t], 0.0)
        fill!(jac.u1T[t], 0.0)

        # Dual regularization
        fill!(jac.reg[t], 0.0)
    end

    nq = size(jac.q0[1])[1]

    for (i, t) in enumerate(window[1:end-2])
        if i >= 3
            jac.q0[i-2]  .+= im_traj.δq0[t][1:nq, :]
            jac.q0T[i-2] .+= im_traj.δq0[t][1:nq, :]'
        end

        if i >= 2
            jac.q1[i-1]  .+= im_traj.δq1[t][1:nq, :]
            jac.q1T[i-1] .+= im_traj.δq1[t][1:nq, :]'
        end

        jac.u1[i]  .+= im_traj.δu1[t][1:nq, :]
        jac.u1T[i] .+= im_traj.δu1[t][1:nq, :]'

        # Dual regularization
        jac.reg[i] .-= 1.0 * β * im_traj.ip[t].κ # TODO sort the κ stuff, maybe make it a prameter of this function
    end

    return nothing
end

function jacobian!(jac::NewtonJacobian, im_traj::ImplicitTrajectory, obj::Objective,
    H::Int, β::T, window) where T

    initialize_jacobian!(jac, obj, H)
    update_jacobian!(jac, im_traj, obj, H, β, window)

    return nothing
end

function residual!(res::NewtonResidual, core::Newton,
    ν::Vector, im_traj::ImplicitTrajectory, traj::ContactTraj, ref_traj::ContactTraj,
    window::Vector{Int})

    # unpack
    opts = core.opts
    obj = core.obj
    res.r .= 0.0

    nq = size(res.q0[1])[1]

    # Objective
    gradient!(res, obj, core, traj, ref_traj)

    for (i, t) in enumerate(window[1:end-2])
        # Implicit dynamics
        res.rd[i] .+= im_traj.d[t][1:nq]

        # Minus Identity term #∇qk1, ∇γk, ∇bk
        res.rI[i] .-= ν[i]

        i >= 3 ? res.q0[i-2] .+= im_traj.δq0[t][1:nq, :]' * ν[i] : nothing
        i >= 2 ? res.q1[i-1] .+= im_traj.δq1[t][1:nq, :]' * ν[i] : nothing
        res.u1[i] .+= im_traj.δu1[t][1:nq, :]' * ν[i]
    end

    return nothing
end

function delta!(Δx::SizedArray{Tuple{nx},T,1,1}, x, x_ref) where {nx,T}
    Δx .= x
    Δx .-= x_ref
    return nothing
end

#TODO: add minus function

function update_traj!(traj_cand::ContactTraj, traj::ContactTraj,
        ν_cand::Vector, ν::Vector, Δ::NewtonResidual{T}, α::T) where T

    # H = traj_cand.H

    # for t = 1:H
    #     traj_cand.q[t+2] .= traj.q[t+2] .- α .* Δ.q2[t]
    #     traj_cand.u[t] .= traj.u[t] .- α .* Δ.u1[t]

    #     ν_cand[t]  .= ν[t] .- α .* Δ.rd[t]
    # end

    # update_z!(traj_cand)
    # update_θ!(traj_cand)

    # return nothing
end

function copy_traj!(traj::ContactTraj, traj_cand::ContactTraj, H::Int)
    Ht = traj.H
    Hs = traj_cand.H # MAYBE BREAKING TEST

    @assert Hs >= H
    @assert Ht >= H

    traj.κ .= traj_cand.κ

    for t in eachindex(1:H + 2)
        traj.q[t] .= traj_cand.q[t]
    end

    for t in eachindex(1:H)
        traj.u[t] .= traj_cand.u[t]
        traj.w[t] .= traj_cand.w[t]
        traj.z[t] .= traj_cand.z[t]
        traj.θ[t] .= traj_cand.θ[t]
    end

    return nothing
end

function reset!(core::Newton, ref_traj::ContactTraj;
    warm_start::Bool = false, initial_offset::Bool = false,
    q0 = ref_traj.q[1], q1 = ref_traj.q[2])

    # H = ref_traj.H
    H_mpc = core.traj.H
    opts = core.opts

    # Reset β value
    core.β = opts.β_init

    if !warm_start
		# Reset duals
		for t = 1:H_mpc
			fill!(core.ν[t], 0.0)
			fill!(core.ν_cand[t], 0.0)
		end

		# Set up trajectory
        copy_traj!(core.traj, ref_traj, core.traj.H)

        if initial_offset
		    rd = -0.03 * [[1, 1]; ones(model.nq - 2)]
		    core.traj.q[1] .+= rd
		    core.traj.q[2] .+= rd
			update_θ!(core.traj, 1)
			update_θ!(core.traj, 2)
		end
	end
    core.traj.q[1] .= deepcopy(q0)
    core.traj.q[2] .= deepcopy(q1)

    update_θ!(core.traj, 1)
    update_θ!(core.traj, 2)

    # initialized residual Jacobian
    initialize_jacobian!(core.jac, core.obj, core.traj.H)

	# Set up traj cand
	core.traj_cand = deepcopy(core.traj)

	return nothing
end

function newton_solve!(core::Newton, s::Simulation,
    window::Vector{Int},
    im_traj::ImplicitTrajectory, ref_traj::ContactTraj;
    warm_start::Bool = false, initial_offset::Bool = false,
    q0 = ref_traj.q[1], q1 = ref_traj.q[2])

	reset!(core, ref_traj, warm_start = warm_start,
        initial_offset = initial_offset, q0 = q0, q1 = q1)

    # Compute implicit dynamics about traj
	implicit_dynamics!(im_traj, core.traj, window=window)

    # Compute residual
    residual!(core.res, core, core.ν, im_traj, core.traj, ref_traj, window)

    r_norm = norm(core.res.r, 1)

    for l = 1:core.opts.max_iter
        # check convergence
        # println("lbefor = ", l, "  norm = ", r_norm / length(core.res.r)) #@@@
        r_norm / length(core.res.r) < core.opts.r_tol && break
        # Compute NewtonJacobian
        update_jacobian!(core.jac, im_traj, core.obj, core.traj.H, core.β, window)

        # Compute Search Direction
        linear_solve!(core.solver, core.Δ.r, core.jac.R, core.res.r)

        # line search the step direction
        α = 1.0
        iter = 0

        # candidate step
        update_traj!(core.traj_cand, core.traj, core.ν_cand, core.ν, core.Δ, α)

        # Compute implicit dynamics for candidate
		implicit_dynamics!(im_traj, core.traj_cand, window=window)

        # Compute residual for candidate
        residual!(core.res_cand, core, core.ν_cand, im_traj, core.traj_cand, ref_traj, window)
        r_cand_norm = norm(core.res_cand.r, 1)

        while r_cand_norm^2.0 >= (1.0 - 0.001 * α) * r_norm^2.0
            α = 0.5 * α

            iter += 1
            if iter > 6
                break
            end

            update_traj!(core.traj_cand, core.traj, core.ν_cand, core.ν, core.Δ, α)

            # Compute implicit dynamics about trial_traj
			implicit_dynamics!(im_traj, core.traj_cand, window=window)

            residual!(core.res_cand, core, core.ν_cand, im_traj, core.traj_cand, ref_traj, window)
            r_cand_norm = norm(core.res_cand.r, 1)
        end

        # update
        update_traj!(core.traj, core.traj, core.ν, core.ν, core.Δ, α)
        core.res.r .= core.res_cand.r
        r_norm = r_cand_norm
        # println("lafter = ", l, "  norm = ", r_norm / length(core.res.r)) #@@@

        # regularization update
        if iter > 6
            core.β = min(core.β * 1.3, 1.0e2)
        else
            core.β = max(1.0e1, core.β / 1.3)
        end

        # print status
        core.opts.verbose && println(" l: ", l ,
                "     r̄: ", scn(norm(core.res_cand.r, 1) / length(core.res_cand.r), digits = 0),
                "     r: ", scn(norm(core.res.r, 1) / length(core.res.r), digits = 0),
                "     Δ: ", scn(norm(core.Δ.r, 1) / length(core.Δ.r), digits = 0),
                "     α: ", -Int(round(log(α))),
                "     κ: ", scn(im_traj.ip[1].κ[1], digits = 0))
    end

    return nothing
end


text:
"""
    load_sigmf(f; mmap=true, integrity_check=true)

f is a string pointing to the location of a SigMF file. In practice,
few data sets are proper ".sigmf" files, but rather pairs of ".sigmf-data" and ".sigmf-meta" files.
In this case you may point to them with the ".sigmf" extension and the loader will
attempt to correct the extensions to load both.

Returns a `SigMFFile`.
"""
function load(f; mmap=true, integrity_check=false)

    exists = isfile(f)
    if !exists
        @warn "$f not found, attempting with .sigm-meta and .sigm-data extensions"
        metafile = f*"-meta"
        datafile = f*"-data"
        dataexists = isfile(datafile)
        metaexists = isfile(metafile)
        !dataexists && !metaexists && error("Could not locate data and meta SigMF files")
    else
        error("TODO: TAR")
    end

    meta = open(metafile, "r") do f
        JSON3.read(f, SigMFMeta)
    end

    # determine data type from meta
    datatype, endianness = datatype_sigmf2julia(meta.meta_global.datatype)

    dataio = open(datafile, "r")
    data = mmap ? Mmap.mmap(dataio) : read(dataio)

    # integrity checks
    hassha = meta.meta_global.sha512 !== nothing
    integrity_check && !hassha && @warn "SigMF metadata does not contain a SHA512 hash, skipping integrity checks"
    integrity_check = integrity_check && !hassha
    integrity_check && sha512(dataio) != meta.meta_global.sha512 && error("SigMF data integrity check failed")

    SigMFFile(meta, reinterpret(datatype, data))
end

"""
    SigMF.save(f, ::SigMFFile; compress=false)
"""
function save(f, sigmf::SigMFFile; compress=false)
    # write meta
    save(f*"-meta", sigmf.meta)

    # write data
    open(f*"-data", "w") do io
        write(io, sigmf.data)
    end
end

function save(f, sigmf::SigMFMeta)
    open(f, "w") do f
        JSON3.write(f, sigmf)
    end
end

ext: jl


text:
# julia test/try7.jl /Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/images/artificial_reference.jpg /Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/images/artificial_moving-affine.jpg /Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/landmarks/artificial_moving-affine.csv  /Users/MrTrololord/Google_Drive/cvut/bakalarka/LAP_julia/test/outputs/ sparse_pflap

function move_landmarks(locations, flow)
    return map(loc -> loc .+ (real.(flow)[loc...], imag.(flow)[loc...]), locations)
end

flow[1,1]
move_landmarks([(1,1)], flow)

args_dict = Dict()

target = load_image_gray("/Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/images/artificial_reference.jpg")
source = load_image_gray("/Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/images/artificial_moving-affine.jpg")

args_dict["reg_alg"] = sparse_pflap_psnr
timer = TimerOutput(string(args_dict["reg_alg"]));

showflow(flow)

# target, source, flow = gen_init(:spaghetti, img_args=[(12,10)], flow_args=[2])

method_kwargs =Dict(:timer => timer, :display => true, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, time_in_secs, (figs,) = time_reg_alg(args_dict["reg_alg"],
                                                         target,
                                                         source,
                                                         timer=timer,
                                                         method_kwargs=method_kwargs);
println()

showflow(flow_est)
imgoverlay(target, source_reg)
imgoverlay(target, source)

figs[1,1,3]

target, source, flow = gen_init()

args_dict["reg_alg"] = pflap
timer = TimerOutput(string(args_dict["reg_alg"]))
method_kwargs = Dict(:timer => timer, :display => false, :max_repeats => 2, :match_source_histogram => false)
flow_est, source_reg, timer, time_in_secs = time_reg_alg(args_dict["reg_alg"],
                                                         target,
                                                         source,
                                                         timer=timer,
                                                         method_kwargs=Dict(:display => true, :timer => timer))
println()

imgshow(source_reg)
imgshow(target)
imgoverlay(target, source_reg)

#
# y = [1 2; 3 4]
# kern = centered([0 1 0;1 -4 1;0 1 0])
#
# imfilter(y, kern, Fill(0,kern))
#
# import LAP_julia: estimation_noise_variance
#
# LAP_julia.estimation_noise_variance(target)
#
# img, imgw, flow = gen_init()
#
# showflow(flow)
# imgshow(imgw)
#
# LAP_julia.mapped_out(flow.*4)
#
# imgw_copy = deepcopy(imgw)
# imgw_copy[LAP_julia.mapped_out(flow)] .= 0
# imgshow(imgw_copy)
# imgshow(imgw, origin_bot = true)
#

assess_psnr(target, source)

assess_psnr(img, imgw)
assess_psnr(img, source_reg)

imgshow(source_reg)

imgshow(img)
imgshow(imgw)
imgoverlay(img,imgw)
imgoverlay(source_reg, img)

timer=TimerOutput("pf lap");
method_kwargs =Dict(:timer => timer, :display => false, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, results = test_registration_alg(pflap, img, imgw, flow, method_kwargs=method_kwargs, timer=timer);

showflow(flow, figtitle="truth flow")

figs[1,1,1]
figs[6,1,1]

img, imgw, flow = gen_init(flow_args=[30]);

timer=TimerOutput("pflap reg");
method_kwargs =Dict(:timer => timer, :display => false, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, results = test_registration_alg(pflap, img, imgw, flow, method_kwargs=method_kwargs, timer=timer);

timer=TimerOutput("sparse pflap reg"); # fix
method_kwargs =Dict(:timer => timer, :display => true, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, results = test_registration_alg(sparse_pflap, img, imgw, flow, method_kwargs=method_kwargs, timer=timer);

timer=TimerOutput("sparse pflap psnr reg"); # fix
method_kwargs =Dict(:timer => timer, :display => true, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, results, (figs,) = test_registration_alg(sparse_pflap_psnr, img, imgw, flow, method_kwargs=method_kwargs, timer=timer);

timer=TimerOutput("sparse lap reg");
method_kwargs =Dict(:timer => timer, :display => false)
flow_est, source_reg, timer, results = test_registration_alg(sparse_lap, img, imgw, flow, method_args=[30], method_kwargs=method_kwargs, timer=timer);

timer=TimerOutput("lap reg");
method_kwargs =Dict(:timer => timer, :display => false)
flow_est, source_reg, timer, results = test_registration_alg(lap, img, imgw, flow, method_args=[30], method_kwargs=method_kwargs, timer=timer);

imgshow(img, origin_bot = true, figtitle="target")
imgshow(imgw, origin_bot = true, figtitle="source")
imgshow(source_reg, origin_bot = true, figtitle="source_reg")

imgshow(img .- source_reg, origin_bot = true, figtitle="diff")

showflow(flow.*(-1), figtitle="target flow")
showflow(flow_est, figtitle="estimated flow")
showflow((flow_est .+ flow)[30:end-30, 30:end-30], figtitle="flow diff")

###

#anhir loc table subset

base_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/"
loc_table = CSV.read(base_path * "location_table.csv")

train_rows = loc_table[loc_table[:status] .== "training", :]

#first k rows:
k = 15
subset = train_rows[1:k, :]

CSV.write(joinpath(base_path, "subset_table.csv"), subset)

##
using LAP_julia, TimerOutputs

target, source = gen_anhir(mutate=true, diag_pixels=500)

imgshow(target)
imgshow(source)

timer=TimerOutput("pflap reg");
method_kwargs =Dict(:timer => timer, :display => false)
flow_est, source_reg, timer, results = time_reg_alg(pflap, target, source, method_kwargs=method_kwargs, timer=timer);

timer=TimerOutput("sparse pflap reg");
method_kwargs =Dict(:timer => timer, :display => false)
flow_est, source_reg, timer, results = time_reg_alg(sparse_pflap_psnr, target, source, method_kwargs=method_kwargs, timer=timer);

imgshow(source_reg)

showflow(flow_est)

A = rand(1234,1234)

@benchmark LAP_julia.resize_to_diag_size(A, 500)

flow_est

## tmp

using ArgParse, CSV, TimerOutputs, FileIO, Colors, DataFrames
using LAP_julia

function save_string_to_file(str, file_path)
    open(file_path, "w") do io
        write(io, str)
    end
end

function load_image_gray(img_path)
    img = Float64.(Gray.(load(img_path)))
end

function save_shift_landmarks(source_landmarks_path, flow, args_dict)
    source_landmarks_df = CSV.read(source_landmarks_path) |> DataFrame
    locations = [source_landmarks_df[:, "Y"] source_landmarks_df[:, "X"]]
    shifted_landmarks = LAP_julia.move_landmarks(locations, flow)
    shifted_landmarks_df = DataFrame(Column1 = 0:(size(locations,1)-1), X = shifted_landmarks[:,2], Y = shifted_landmarks[:,1])
    save_path = joinpath(args_dict["output_path"], args_dict["land_warped_fname"])
    CSV.write(save_path, shifted_landmarks_df)
end

function save_time(time_in_secs, args_dict)
    save_string_to_file(string(time_in_secs), joinpath(args_dict["output_path"], args_dict["time_fname"]))
end

args_dict = Dict("img_path" => "/Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/images/artificial_reference.jpg",
                 "imgw_path" => "/Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/images/artificial_moving-affine.jpg",
                 "moving_landmark_path" => "/Users/MrTrololord/Google_Drive/cvut/bakalarka/BIRL/data-images/landmarks/artificial_moving-affine.csv",
                 "output_path" => "/Users/MrTrololord/Google_Drive/cvut/bakalarka/LAP_julia/test/outputs/",
                 "reg_alg" => "sparse_pflap_psnr",
                 "time_fname" => "exec_time_in_secs.txt",
                 "land_warped_fname" => "warped_landmarks.csv",
                 "diag_pix" => 500)

base_path = "mount/birl_small/"
args_dict = Dict("img_path" => joinpath(base_path, "images/COAD_01/scale-5pc/S3.jpg"),
                 "imgw_path" => joinpath(base_path, "images/COAD_01/scale-5pc/S1.jpg"),
                 "moving_landmark_path" => joinpath(base_path, "landmarks/COAD_01/scale-5pc/S3.csv"),
                 "output_path" => joinpath(base_path, "tmp/"),
                 "reg_alg" => "pflap",
                 "time_fname" => "exec_time_in_secs.txt",
                 "land_warped_fname" => "warped_landmarks.csv",
                 "diag_pix" => 500)

println("Running with the arguments:")

println("******ARGS:********")
args_dict["reg_alg"] = getfield(LAP_julia, Symbol(args_dict["reg_alg"]))

target = load_image_gray(args_dict["img_path"])
source = load_image_gray(args_dict["imgw_path"])

resized_source, source_resize_ratio = resize_to_diag_size(source, args_dict["diag_pix"])
resized_target = imresize(target, size(resized_source))

# Take the 3 run of the reg alg, allowing for julia compilation.
for _ in 1:3
    timer = TimerOutput(string(args_dict["reg_alg"]))
    flow_est, source_reg, timer, time_in_secs = time_reg_alg(args_dict["reg_alg"],
                                                                resized_target,
                                                                resized_source,
                                                                timer=timer,
                                                                method_kwargs=Dict(:display => false, :timer => timer, :match_source_histogram => true))
    global flow_est, source_reg, timer, time_in_secs
    println()

end

# save time
save_time(time_in_secs, args_dict)

# rescale flow_est
flow_est = imresize(flow_est.*(1/source_resize_ratio), size(source))
@assert size(flow_est) == size(source) size(flow_est), size(source)

# move landmarks
save_shift_landmarks(args_dict["moving_landmark_path"], flow_est, args_dict)

# save source_reg
colored_source = load_image(args_dict["imgw_path"])
source_reg = warp_img(colored_source, -real.(flow_est), -imag.(flow_est))
save_path = joinpath(args_dict["output_path"], "source_reg.jpg")
save(save_path, source_reg)

addpoints(locations)
addpoints(shifted_landmarks)

colorview(Gray, target)


text:
using JSON

## We start with definitions of instruction mnemonics, broken down by category:
const instruction_categories = JSON.parsefile(joinpath(@__DIR__, "instructions.json");
                                              dicttype=Dict{String,Vector{String}})

# Turn instructions "inside out", so e.g. we have "vzeroall" => "avx"
mnemonics_by_category = Dict(
    inst => cat for (cat, insts) in instruction_categories for inst in insts
)

"""
    instruction_mnemonics(path::AbstractString, platform::AbstractPlatform)

Dump a binary object with `objdump`, returning a list of instruction mnemonics
for further analysis with `analyze_instruction_set()`.

Note that this function only really makes sense for x86/x64 binaries.  Don't
run this on armv7l, aarch64, ppc64le etc... binaries and expect it to work.

This function returns the list of mnemonics as well as the counts of each,
binned by the mapping defined within `instruction_categories`.
"""
function instruction_mnemonics(path::AbstractString, platform::AbstractPlatform)
    # The outputs we are calculating
    counts = Dict(k => 0 for k in keys(instruction_categories))
    mnemonics = String[]

    ur = preferred_runner()(
        abspath(dirname(path));
        cwd="/workspace/",
        platform=platform,
        verbose=false,
    )
    output = IOBuffer()

    # Run objdump to disassemble the input binary
    if Sys.isbsd(platform)
        objdump_cmd = "llvm-objdump -d $(basename(path))"
    else
        objdump_cmd = "\${target}-objdump -d $(basename(path))"
    end
    run_interactive(ur, `/bin/bash -c "$(objdump_cmd)"`; stdout=output, stderr=devnull)
    seekstart(output)

    for line in eachline(output)
        # First, ensure that this line of output is 3 fields long at least
        fields = filter(x -> !isempty(strip(x)), split(line, '\t'))
        if length(fields) < 3
            continue
        end

        # Grab the mnemonic for this line as the first word of the 3rd field
        m = split(fields[3])[1]
        push!(mnemonics, m)

        # For each mnemonic, find it in mnemonics_by_category, if we can, and
        # increment the appropriate `counts` member:
        if haskey(mnemonics_by_category, m)
            counts[mnemonics_by_category[m]] += 1
        else
            counts["unknown"] += 1
        end
    end

    # Return both the list of mnemonics as well as the binned counts
    return mnemonics, counts
end

function generic_march(p::AbstractPlatform)
    return first(first(Base.BinaryPlatforms.arch_march_isa_mapping[arch(p)]))
end

"""
    minimum_march(counts::Dict, p::AbstractPlatform)

This function returns the minimum instruction set required, depending on
whether the object file being pointed to is a 32-bit or 64-bit one:

* For 32-bit object files, this returns one of ["i686", "prescott"]

* For 64-bit object files, this returns one of ["x86_64", "avx", "avx2", "avx512"]
"""
function minimum_march(counts::Dict, p::AbstractPlatform)
    if arch(p) == "x86_64"
        avx512_instruction_categories = (
            "pku", "rdseed", "adcx", "clflush", "xsavec",
            "xsaves", "clwb", "avx512evex", "avex512vex",
        )
        avx2_instruction_categories = (
            "movbe", "avx2", "rdwrfsgs", "fma", "bmi1", "bmi2", "f16c",
        )
        # note that the extensions mmx, sse, and sse2 are part of the generic x86-64 architecture
        avx_instruction_categories = (
            "sse3", "ssse3", "sse4", "avx", "aes", "pclmulqdq",
        )
        if any(get.(Ref(counts), avx512_instruction_categories, 0) .> 0)
            return "avx512"
        elseif any(get.(Ref(counts), avx2_instruction_categories, 0) .> 0)
            return "avx2"
        elseif any(get.(Ref(counts), avx_instruction_categories, 0) .> 0)
            return "avx"
        end
    elseif arch(p) == "i686"
        if counts["sse3"] > 0
            return "prescott"
        end
    elseif arch(p) == "aarch64"
        # TODO: Detect instructions for aarch64 extensions
    elseif arch(p) == "armv6l"
        # We're just always going to assume we're running the single armv6l that Julia runs on.
    elseif arch(p) == "armv7l"
        # TODO: Detect NEON and vfpv4 instructions
    elseif arch(p) == "powerpc64le"
        # TODO Detect POWER9/10 instructions
    end
    return generic_march(p)
end

"""
    analyze_instruction_set(oh::ObjectHandle, platform::AbstractPlatform; verbose::Bool = false)

Analyze the instructions within the binary located at the given path for which
minimum instruction set it requires, taking note of groups of instruction sets
used such as `avx`, `sse4.2`, `i486`, etc....

Some binary files (such as libopenblas) contain multiple versions of functions,
internally determining which version to call by using the `cpuid` instruction
to determine processor support.  In an effort to detect this, we make note of
any usage of the `cpuid` instruction, disabling our minimum instruction set
calculations if such an instruction is found, and notifying the user of this
if `verbose` is set to `true`.

Note that this function only really makes sense for x86/x64 binaries.  Don't
run this on armv7l, aarch64, ppc64le etc... binaries and expect it to work.
"""
function analyze_instruction_set(oh::ObjectHandle, platform::AbstractPlatform; verbose::Bool = false)
    # Get list of mnemonics
    mnemonics, counts = instruction_mnemonics(path(oh), platform)

    # Analyze for minimum instruction set
    min_march = minimum_march(counts, platform)

    # If the binary uses `cpuid`, we can't know what it's doing, so just
    # return the most conservative ISA and warn the user if `verbose` is set.
    if counts["cpuid"] > 0
        if verbose && generic_march(platform) != min_march
            msg = replace("""
            $(basename(path(oh))) contains a `cpuid` instruction; refusing to
            analyze for minimum instruction set, as it may dynamically select
            the proper instruction set internally.  Would have chosen
            $(min_march), instead choosing $(generic_march(p)).
            """, '\n' => ' ')
            @warn(strip(msg))
        end
        return generic_march(platform)
    end

    # Otherwise, return `min_march` and let 'em know!
    return min_march
end


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using CAGD
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
if todisplay
    using ViewerGL
	GL = ViewerGL
	include("views.jl")
end

V = 2 * [
    0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0
    0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
    0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
]
EV = [
    [1, 2], [1, 3], [2, 4], [3, 4],
    [1, 5], [2, 6], [3, 7], [4, 8],
    [5, 6], [5, 7], [6, 8], [7, 8]
]
FE = [
    [1, 2, 3, 4], [9, 10, 11, 12],
    [1, 5, 6, 9], [2, 5, 7, 10], [3, 6, 8, 11], [4, 7, 8, 12]
]
CF = [
    [1, 2, 3, 4, 5, 6]
]

m1 = CAGD.Model(V)
CAGD.addModelCells!(m1, 1, EV, signed = true)
CAGD.addModelCells!(m1, 2, FE, signed = true)
CAGD.addModelCells!(m1, 3, CF, signed = false)

m2 = deepcopy(m1)
m2.G .+= [1.0; 1.0; 1.0]

model = CAGD.uniteModels(m1, m2)

if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model)

FCcycles, bicon_comps = CAGD.tgw(congr_model, 3)

externals = [bc[CAGD.getExternalBoundary(congr_model, 3, FCcycles, bc)] for bc in bicon_comps]
internals = [setdiff(bicon_comps[i], externals[i]) for i = 1 : length(externals)]

FC = FCcycles[:, internals...]
JK = ones(Int8, length(externals))
FC = [FCcycles * SparseArrays.sparse(externals, JK, JK, size(FCcycles, 2), 1) FC]

CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model)  end

# Each face is percurred two times
sum(gift_model.T[3], dims = 1)

boolmatrix, arranged_model = CAGD.bool3(model)

arranged_model == gift_model


text:
struct GitHub <: Remote
    repo_name::String
    username::String
    token::String
end

GitHub(repo_name) =
    GitHub(repo_name, settings("username"), settings("github_token"))

base_url(g::GitHub) = "https://api.github.com"
headers(g::GitHub) = Dict(
    "Authorization" => "token $(g.token)",
    "User-Agent" => user_agent
)

function repos(g::GitHub)
    @info "getting github repos"
    talk_to(HTTP.get, g, "/user/repos?per_page=100") |> json_parse
end

function create(g::GitHub)
    @info "creating github"
    talk_to(HTTP.post, g, "/user/repos", json(Dict(
        "name" => g.repo_name)))
end

function get_keys(g::GitHub)
    @info "getting github keys"
    talk_to(HTTP.get, g, "/repos/$(g.username)/$(g.repo_name)/keys") |> json_parse
end

function add_key(g::GitHub, name, value; read_only = false)
    @info "creating github key"
    talk_to(HTTP.post, g, "/repos/$(g.username)/$(g.repo_name)/keys", json(Dict(
        "title" => name,
        "key" => value,
        "read_only" => read_only
    )))
end

function delete_key(g::GitHub, key_id)
    @info "deleting github key"
    talk_to(HTTP.delete, g, "/repos/$(g.username)/$(g.repo_name)/keys/$key_id")
end

function delete_keys(g::GitHub, name)
    foreach(
        github_key ->
            if github_key["title"] == name
                delete_key(g, github_key["id"])
            end,
        get_keys(g)
    )
end

# TODO: check more than the first 100 repositories
function exists(g::GitHub)
    any(repos(g)) do repo
        repo["name"] == g.repo_name
    end
end

function delete(g::GitHub)
    @info "deleting github"
    talk_to(HTTP.delete, g, "/repos/$(g.username)/$(g.repo_name)")
end


text:
using Pkg
Pkg.activate(".")
using Vlasiator
using BenchmarkTools

#file = "bulk.singleprecision.vlsv" # 80 KB rho
#file = "bulk.0000003.vlsv"        # 900 KB rho
file = "bulk1.0001000.vlsv"       # 32 MB rho
#file = "test/data/bulk.2d.vlsv"

@time meta = load(file);

@benchmark meta = load($file)

meta = load(file);

var = "proton/vg_rho";

@benchmark rho = meta[$var]

@benchmark rho_unsorted = readvariable($meta, $var, false)

println("")


text:
module RigidBodyTreeInspector

using LinearAlgebra
using RigidBodyDynamics
using RigidBodyDynamics.OdeIntegrators
import DrakeVisualizer: Visualizer,
                        settransform!,
                        setgeometry!,
                        addgeometry!,
                        batch,
                        GeometryData,
                        Triad,
                        HyperEllipsoid,
                        center
using StaticArrays: SVector, SMatrix
using CoordinateTransformations: AffineMap, IdentityTransformation, AngleAxis,
                                 LinearMap, RodriguesVec, Quat, compose,
                                 Translation, RotMatrix, RotXYZ, transform_deriv
using GeometryTypes: AbstractGeometry,
                     AbstractMesh,
                     HyperRectangle,
                     HyperSphere,
                     Vec, Point
using ColorTypes: RGBA
import Interpolations: interpolate, Linear, Gridded
import Base: convert, *, one
import Rotations
import LoopThrottle: @throttle
using MechanismGeometries: VisualElement, Skeleton, AbstractGeometrySource,
                           visual_elements, URDFVisuals

const rbd = RigidBodyDynamics

rbd.parse_urdf

export manipulate,
       inspect,
       inspect!,
       Visualizer,
       settransform!,
       addgeometry!,
       setgeometry!,
       animate,
       parse_urdf,
       create_geometry,
       DrakeVisualizerSink,
       visual_elements,
       Skeleton,
       URDFVisuals

include("manipulate.jl")
using .Manipulate
include("visualizer.jl")
include("animate.jl")
include("ode_callback.jl")

end


text:
"""
$(TYPEDEF)

It is used to handle the following situation. Assume you want to solve `F(x)=0` with a Newton algorithm but you want to avoid the process to return some already known solutions ``roots_i``. The deflation operator penalizes these roots ; the algorithm works very well despite its simplicity. You can use `DeflationOperator` to define a function `M(u)` used to find, with Newton iterations, the zeros of the following function
``F(u) \\cdot Π_i(dot(u - roots_i, u - roots_i)^{-p} + shift) := F(u) \\cdot M(u)``. The fields of the struct `DeflationOperator` are as follows:

$(TYPEDFIELDS)

The deflation operator is is ``M(u) = \\prod_{i=1}^{n_{roots}}(shift + norm(u-roots_i)^{-p})``
where ``norm(u) = dot(u,u)``.

Given `defOp::DeflationOperator`, one can access its roots as `defOp[n]` as a shortcut for `defOp.roots[n]`. Also, one can add (resp.remove) a new root by using `push!(defOp, newroot)` (resp. `pop!(defOp)`). Finally `length(defOp)` is a shortcut for `length(defOp.roots)`
"""
struct DeflationOperator{Tp <: Real, Tdot, T <: Real, vectype}
	"power `p`. You can use an  `Int` for example"
	power::Tp

	"function, this function has to be bilinear and symmetric for the linear solver to work well"
	dot::Tdot

	"shift"
	α::T

	"roots"
	roots::Vector{vectype}

	# internal, to reduce allocations during computation
	tmp::vectype
end

# constructors
DeflationOperator(p::Real, α::Real, roots::Vector{vectype}) where vectype = DeflationOperator(p, dot, α, roots, _copy(roots[1]))
DeflationOperator(p::Real, dt, α::Real, roots::Vector{vectype}) where vectype = DeflationOperator(p, dt, α, roots, _copy(roots[1]))

# methods to deal with DeflationOperator
Base.eltype(df::DeflationOperator{Tp, Tdot, T, vectype}) where {Tp, Tdot, T, vectype} = T
Base.push!(df::DeflationOperator{Tp, Tdot, T, vectype}, v::vectype) where {Tp, Tdot, T, vectype} = push!(df.roots, v)
Base.pop!(df::DeflationOperator) = pop!(df.roots)
Base.getindex(df::DeflationOperator, inds...) = getindex(df.roots, inds...)
Base.length(df::DeflationOperator) = length(df.roots)
Base.deleteat!(df::DeflationOperator, id) = deleteat!(df.roots, id)
Base.empty!(df::DeflationOperator) = empty!(df.roots)
Base.firstindex(df::DeflationOperator) = 1
Base.lastindex(df::DeflationOperator) = length(df)
Base.show(io::IO, df::DeflationOperator) = println(io, "Deflation operator with ", length(df.roots)," roots")

# Compute M(u)
function (df::DeflationOperator{Tp, Tdot, T, vectype})(u::vectype) where {Tp, Tdot, T, vectype}
	if length(df.roots) == 0
		return T(1)
	end
	ρ = u -> T(1) / df.dot(u, u)^df.power + df.α
	# compute u - df.roots[1]
	copyto!(df.tmp, u);	axpy!(T(-1), df.roots[1], df.tmp)
	out::T = ρ(df.tmp)
	for ii in 2:length(df.roots)
		copyto!(df.tmp, u); axpy!(T(-1), df.roots[ii], df.tmp)
		out *= ρ(df.tmp)
	end
	return out
end

# Compute dM(u)⋅du. We use a tmp for storing
function (df::DeflationOperator{Tp, Tdot, T, vectype})(::Val{:dMwithTmp}, tmp, u::vectype, du) where {Tp, Tdot, T, vectype}
	if length(df) == 0
		return T(0)
	end
	δ = T(1e-8)
	copyto!(tmp, u); axpy!(δ, du, tmp)
	return (df(tmp) - df(u)) / δ
end

(df::DeflationOperator)(u, du) = df(Val(:dMwithTmp), similar(u), u, du)

"""
	pb = DeflatedProblem(F, J, M::DeflationOperator)

This creates a deflated problem ``M(u) \\cdot F(u) = 0`` where `M` is a `DeflationOperator` which encodes the penalization term. `J` is the jacobian of `F`. Can be used to call `newton` and `continuation`.
"""
struct DeflatedProblem{Tp, Tdot, T, vectype, TF, TJ}
	F::TF
	J::TJ
	M::DeflationOperator{Tp, Tdot, T, vectype}
end
@inline length(prob::DeflatedProblem) = length(prob.M)

"""
Return the deflated function M(u) * F(u) where M(u) ∈ R
"""
function (df::DeflatedProblem{Tp, Tdot, T, vectype, TF, TJ})(u::vectype, par) where {Tp, Tdot, T, vectype, TF, TJ}
	out = df.F(u, par)
	rmul!(out, df.M(u))
	return out
end

"""
Return the jacobian of the deflated function M(u) * F(u) where M(u) ∈ R
"""
function (df::DeflatedProblem{Tp, Tdot, T, vectype, TF, TJ})(u::vectype, par, du) where {Tp, Tdot, T, vectype, TF, TJ}
	# dF(u)⋅du * M(u) + F(u) dM(u)⋅du
	# out = dF(u)⋅du * M(u)
	out = apply(df.J(u, par), du)
	M = df.M(u)
	rmul!(out, M)
	# we add the remaining part
	if length(df) > 0
		F = df.F(u, par)
		# F(u) dM(u)⋅du, out .+= df.M(u, du) .* F
		axpy!(df.M(u, du), F, out)
	end
	return out
end
###################################################################################################
# Implement the Jacobian operator of the Deflated problem

# this is used to define a custom linear solver
struct DeflatedLinearSolver <: AbstractLinearSolver end

"""
Implement the linear solver for the deflated problem.
"""
function (dfl::DeflatedLinearSolver)(J, rhs)
	# the expression of the Functional is now
	# F(u) * Π_i(dot(u - root_i, u - root_i)^{-power} + shift) := F(u) * M(u)
	# the expression of the differential is
	# dF(u)⋅du * M(u) + F(u) dM(u)⋅du

	# the point at which to compute the Jacobian
	u = J[1]
	p = J[2]

	# deflated Problem composite type
	defPb = J[3]
	linsolve = J[4]

	Fu = defPb.F(u, p)
	Mu = defPb.M(u)
	Ju = defPb.J(u, p)

	if length(defPb.M) == 0
		h1, _, it1 = linsolve(Ju, rhs)
		return h1, true, (it1, 0)
	end

	# linear solve for the deflated problem. We note that Mu ∈ R
	# hence dM(u)⋅du is a scalar. We now solve the following linear problem
	# M(u) * dF(u)⋅h + F(u) dM(u)⋅h = rhs
	h1, h2, _, (it1, it2) = linsolve(Ju, rhs, Fu)

	# We look for the expression of dM(u)⋅h
	# the solution is then h = Mu * h1 - z h2 where z has to be determined
	# z1 = dM(h)⋅h1
	tmp = similar(u)
	z1 = defPb.M(Val(:dMwithTmp), tmp, u, h1)

	# z2 = dM(h)⋅h2
	z2 = defPb.M(Val(:dMwithTmp), tmp, u, h2)

	z = z1 / (Mu + z2)

	# we extract the type of defPb
	_T = eltype(defPb.M)

	# return (h1 - z * h2) / Mu, true, (it1, it2)
	copyto!(tmp, h1)
	axpy!(-z, h2, tmp)
	rmul!(tmp, _T(1) / Mu)
	return tmp, true, (it1, it2)
end
####################################################################################################
"""
$(SIGNATURES)

This is the deflated version of the Krylov-Newton Solver for `F(x, p0) = 0` with Jacobian `J(x, p0)`. We refer to [`newton`](@ref) for more information. It penalises the roots saved in `defOp.roots`. The other arguments are as for `newton`. See [`DeflationOperator`](@ref) for more information.

# Arguments
Compared to [`newton`](@ref), the only different arguments are
- `defOp::DeflationOperator` deflation operator
- `linsolver` linear solver used to invert the Jacobian of the deflated functional. We have a custom solver `DeflatedLinearSolver()` with requires solving two linear systems `J⋅x = rhs`. For other linear solvers, a matrix free method is used for the deflated functional.

# Output:
- solution:
- history of residuals
- flag of convergence
- number of iterations

# Simplified call
When `J` is not passed. It then computed with finite differences. The call is as follows:

	newton(F, x0, p0, options, defOp; kwargs...)
"""
function newton(F, J, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, linsolver::DeflatedLinearSolver = DeflatedLinearSolver(); kwargs...) where {T, L, E, Tp, Tdot, vectype}
	# we create the new functional
	deflatedPb = DeflatedProblem(F, J, defOp)

	# and its jacobian
	Jacdf = (u, p) -> (u, p, deflatedPb, options.linsolver)

	# change the linear solver
	opt_def = @set options.linsolver = linsolver

	return newton(deflatedPb, Jacdf, x0, p0, opt_def; kwargs...)
end

function newton(F, J, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, linsolver::AbstractLinearSolver; kwargs...) where {Tp, T, Tdot, vectype, L, E}
	# we create the new functional
	deflatedPb = DeflatedProblem(F, J, defOp)

	# change the linear solver
	opt_def = @set options.linsolver = linsolver

	return newton(deflatedPb, (x,p) -> (dx -> deflatedPb(x,p,dx)), x0, p0, opt_def; kwargs...)
end

# simplified call when no Jacobian is given
function newton(F, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, linsolver = DeflatedLinearSolver(); kwargs...) where {Tp, T, Tdot, vectype, L, E}
	J = (u, p) -> finiteDifferences(z -> F(z,p), u)
	return newton(F, J, x0, p0, options, defOp, linsolver; kwargs...)
end

"""
$(TYPEDEF)

This specific Newton-Kyrlov method first tries to converge to a solution `sol0` close the guess `x0`. It then attempts to converge to the guess `x1` while avoiding the previous solution `sol0`. This is very handy for branch switching. The method is based on a deflated Newton-Krylov solver.
"""
function newton(F, J, x0::vectype, x1::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator = DeflationOperator(2, dot, 1.0, Vector{vectype}(), _copy(x0)); kwargs...) where {T, vectype, L, E}
	res0 = newton(F, J, x0, p0, options; kwargs...)
	@assert res0[3] "Newton did not converge to the trivial solution x0."
	push!(defOp, res0[1])
	res1 = newton(F, J, x1, p0, (@set options.maxIter = 10options.maxIter), defOp; kwargs...)
	@assert res1[3] "Deflated Newton did not converge to the non-trivial solution ( i.e. on the bifurcated branch)."
	return res1, res0
end


text:
using Images
using ImageView
using MAT
using LinearAlgebra
using Statistics
using Plots

vars = matread("../../DATA/allFaces.mat")

n = Int(vars["n"])
m = Int(vars["m"])
faces = vars["faces"]
nfaces = Int.(vars["nfaces"])

# We use the first 36 people for training data
trainingFaces = faces[:,1:sum(nfaces[1:36])];
avgFace = mean(trainingFaces, dims = 2) # size n*m by 1

# Compute eigenfaces on mean-subtracted training data
X = trainingFaces - avgFace * ones(1, size(trainingFaces)[2])
U, S, V = svd(X)

imshow(reshape(avgFace, n, m)) # Show avg face
imshow(reshape(U[:,1], n, m))  # Show first eigenface

# Now show eigenface reconstruction of image that was omitted from test set

testFace = faces[:,1+sum(nfaces[1:36])] # First face of person 37
imshow(reshape(testFace, n, m))

testFaceMS = testFace - avgFace
for r = [25, 50, 100, 200, 400, 800, 1600]
    reconFace = avgFace + (U[:,1:r]*(U[:,1:r]'*testFaceMS))
    imshow(reshape(reconFace, n, m))
end

# Project person 2 and 7 onto PC5 and PC6

P1num = 2; # Person number 2
P2num = 7; # Person number 7

P1 = faces[:,1+sum(nfaces[1:P1num-1]):sum(nfaces[1:P1num])]
P2 = faces[:,1+sum(nfaces[1:P2num-1]):sum(nfaces[1:P2num])]

P1 = P1 - avgFace * ones(1, size(P1)[2])
P2 = P2 - avgFace * ones(1, size(P2)[2])

PCAmodes = [5, 6] # Project onto PCA modes 5 and 6
PCACoordsP1 = U[:,PCAmodes]'*P1
PCACoordsP2 = U[:,PCAmodes]'*P2

p1 = scatter(PCACoordsP1[1,:], PCACoordsP1[2,:], mc = :blue, legend = false)
scatter!(PCACoordsP2[1,:], PCACoordsP2[2,:], mc = :red)


text:
@testset "KubeError" begin
    msg = "Error from server (NotFound): pods \"localhost\" not found"
    @test sprint(Base.showerror, KubeError(msg)) == "KubeError: $msg"
end

@testset "create_pod" begin
    @testset "non-pod" begin
        @test_throws ArgumentError create_pod(Dict("kind" => "Job"))
    end
end

@testset "wait_for_running_pod" begin
    function stateful_get_pod_patch(states)
        i = 0
        @patch function get_pod(name)
            i <= length(states) && (i += 1)
            Dict("status" => Dict("phase" => states[i]))
        end
    end

    @testset "basic" begin
        states = ("Pending", "Running", "Succeeded")
        pod = apply(stateful_get_pod_patch(states)) do
            wait_for_running_pod("foo"; timeout=10)
        end

        @test pod["status"]["phase"] == "Running"
    end

    # Note: It does not seems possible to reliably test timeout without using Mocking.
    @testset "timeout" begin
        patch = @patch get_pod(name) = Dict("status" => Dict("phase" => "Pending"))

        apply(patch) do
            @test_throws K8sClusterManagers.TimeoutException wait_for_running_pod("foo"; timeout=1)
        end
    end

    @testset "failure" begin
        states = ("Failed",)
        apply(stateful_get_pod_patch(states)) do
            @test_throws ErrorException wait_for_running_pod("foo"; timeout=1)
        end
    end
end

@testset "pod_status" begin
    function gen_resource(status)
        state, reason = status
        state_dict = reason !== nothing ? Dict("reason" => reason) : Dict()
        container_status = Dict("state" => Dict(state => state_dict))

        return Dict("kind" => "Pod",
                    "status" => Dict("containerStatuses" => [container_status]))
    end

    @testset "running" begin
        r = gen_resource("running" => nothing)
        @test pod_status(r) == ("running" => nothing)
    end

    @testset "terminated: Completed" begin
        r = gen_resource("terminated" => "Completed")
        @test pod_status(r) == ("terminated" => "Completed")
    end

    @testset "terminated: OOMKilled" begin
        r = gen_resource("terminated" => "OOMKilled")
        @test pod_status(r) == ("terminated" => "OOMKilled")
    end

    @testset "invalid kind" begin
        @test_throws ArgumentError pod_status(Dict("kind" => "Unknown"))
    end

    @testset "multiple container statuses" begin
        r = gen_resource("running" => nothing)
        r["status"]["containerStatuses"] = repeat(r["status"]["containerStatuses"], 2)
        @test_throws ArgumentError pod_status(r)
    end

    @testset "multiple states" begin
        r = gen_resource("running" => nothing)
        r["status"]["containerStatuses"][1]["state"]["extra"] = Dict()
        @test_throws ArgumentError pod_status(r)
    end
end

@testset "worker_pod_spec" begin
    kwargs = (; manager_name="test", image="julia", cmd=`julia`)
    pod = K8sClusterManagers.worker_pod_spec(; kwargs...)

    @test keys(pod) == Set(["apiVersion", "kind", "metadata", "spec"])
    @test pod["apiVersion"] == "v1"
    @test pod["kind"] == "Pod"

    @test keys(pod["metadata"]) == Set(["generateName", "labels"])
    @test pod["metadata"]["generateName"] == "test-worker-"
    @test keys(pod["metadata"]["labels"]) == Set(["manager"])
    @test pod["metadata"]["labels"]["manager"] == "test"

    @test pod["spec"]["restartPolicy"] == "Never"
    @test length(pod["spec"]["containers"]) == 1

    worker = pod["spec"]["containers"][1]
    @test keys(worker) == Set(["name", "image", "command", "resources"])
    @test worker["name"] == "worker"
    @test worker["image"] == "julia"
    @test worker["command"] == ["julia"]
    @test worker["resources"]["requests"]["cpu"] == DEFAULT_WORKER_CPU
    @test worker["resources"]["requests"]["memory"] == DEFAULT_WORKER_MEMORY
    @test worker["resources"]["limits"]["cpu"] == DEFAULT_WORKER_CPU
    @test worker["resources"]["limits"]["memory"] == DEFAULT_WORKER_MEMORY
end

@testset "isk8s" begin
    pod_id = "pode773d78d-9f0d-4003-a6e8-9bef75b89298/" *
             "b52d083fb1f5b45d998895ab758d10fa36b2a53f3bf3128d7cf1d6d36bc67bd6"
    cgroup = """
        11:devices:/kubepods/$pod_id
        10:blkio:/kubepods/$pod_id
        9:hugetlb:/kubepods/$pod_id
        8:net_cls,net_prio:/kubepods/$pod_id
        7:memory:/kubepods/$pod_id
        6:cpuset:/kubepods/$pod_id
        5:perf_event:/kubepods/$pod_id
        4:cpu,cpuacct:/kubepods/$pod_id
        3:println()ids:/kubepods/$pod_id
        2:freezer:/kubepods/$pod_id
        1:name=systemd:/kubepods/$pod_id
        """

    patches = [
        @patch isfile(p) = p == "/proc/self/cgroup"
        @patch open(f, p) = f(IOBuffer(cgroup))
    ]

    apply(patches) do
        @test isk8s()
    end
end


text:
@enum EventKind::UInt8 begin

# the underlying values are also priorites in case the times are equal
# therefore the order of definition implies priority

  QuarantinedEvent
  DetectionEvent
  #DetectedFromQuarantineEvent
  HomeTreatmentSuccessEvent
  QuarantineEndEvent

  GoHospitalEvent

  #DetectedOutsideQuarantineEvent
  #DetectedFromTracingEvent

  TracedEvent
  ReleasedEvent

  # the progression events have low priority to let the immediate actions execute
  BecomeInfectiousEvent

  OutsideInfectionEvent
  TransmissionEvent

  MildSymptomsEvent
  HomeTreatmentEvent

  SevereSymptomsEvent
  CriticalSymptomsEvent
  RecoveryEvent
  DeathEvent

  ScreeningEvent

  InvalidEvent # should not be executed

end

struct Event
  time::TimePoint                 # 4 bytes
  subject_id::PersonIdx           # 4 bytes
  source_id::PersonIdx            # 4 bytes
  event_kind::EventKind           # 1 byte
  extra::UInt8
  strain::StrainKind

  #contact_kind::ContactKind       # 1 byte
  #extension::Bool                 # 1 byte
  #detection_kind::DetectionKind   # 1 byte
                                  # alignment = 4 bytes

  Event() = new(0.0, 0, 0, InvalidEvent, 0, NullStrain)
  Event(::Val{E}, time::Real, subject::Integer) where E = new(time, subject, 0, E, 0, NullStrain)
  Event(::Val{OutsideInfectionEvent}, time::Real, subject::Integer, strain::StrainKind) = new(time, subject, 0, OutsideInfectionEvent, UInt8(OutsideContact), strain)
  Event(::Val{TransmissionEvent}, ::Real, ::Integer) = error("source and contact kind needed for transmission event")
  Event(::Val{TransmissionEvent}, time::Real, subject::Integer, source::Integer, contact_kind::ContactKind, strain::StrainKind) = new(time, subject, source, TransmissionEvent, UInt8(contact_kind), strain)
  Event(::Val{QuarantinedEvent}, time::Real, subject::Integer, extension::Bool) = new(time, subject, 0, QuarantinedEvent, UInt8(extension), NullStrain)
  Event(::Val{TracedEvent}, ::Real, ::Integer) = error("source and tracing kind must be given for TracedEvent")
  Event(::Val{TracedEvent}, time::Real, subject::Integer, source::Integer, tracing_kind::TracingKind) = new(time, subject, source, TracedEvent, UInt8(tracing_kind), NullStrain)
  Event(::Val{DetectionEvent}, ::Real, ::Integer) = error("detection kind must be given for detection event")
  Event(::Val{DetectionEvent}, time::Real, subject::Integer, detectionkind::DetectionKind) = new(time, subject, 0, DetectionEvent, UInt8(detectionkind), NullStrain)
  Event(::Val{ScreeningEvent}, time::Real) = new(time, 0, 0, ScreeningEvent, 0, NullStrain)
end

time(event::Event) = event.time
subject(event::Event) = event.subject_id
source(event::Event) = event.source_id
kind(event::Event) = event.event_kind

contactkind(event::Event) = istransmission(event) ? ContactKind(event.extra) : NoContact
detectionkind(event::Event) = isdetection(event) ? DetectionKind(event.extra) : NoDetection
extension(event::Event) = isquarantine(event) ? Bool(event.extra) : false
tracingkind(event::Event) = istracing(event) ? TracingKind(event.extra) : NotTraced

strainkind(event::Event) = istransmission(event) ? event.strain : NullStrain

function show(io::IO, event::Event)
  print(io, time(event), ":", kind(event), " ", subject(event))
  if TransmissionEvent == kind(event) || OutsideInfectionEvent == kind(event)
    print(io, " <= ", source(event), " ", contactkind(event), " ", strainkind(event))
  elseif QuarantinedEvent == kind(event)
    print(io, " extension=", extension(event))
  elseif TracedEvent == kind(event)
    print(io, " <= ", source(event))
  elseif InvalidEvent == kind(event)
    print(io, " <= ", source(event), " ", kind(event), " ", event.extra)
  end
end

isdetection(ek::EventKind) = ek == DetectionEvent
isdetection(e::Event) = isdetection(kind(e))
istransmission(ek::EventKind) = ek == TransmissionEvent || ek == OutsideInfectionEvent
istransmission(e::Event) = istransmission(kind(e))
isquarantine(ek::EventKind) = ek == QuarantinedEvent || ek == QuarantineEndEvent
isquarantine(e::Event) = isquarantine(kind(e))
istracing(ek::EventKind) = ek == TracedEvent
istracing(e::Event) = istracing(kind(e))
ishospitalization(ek::EventKind) = ek == GoHospitalEvent
ishospitalization(e::Event) = ishospitalization(kind(e))
isdeath(ek::EventKind) = ek == DeathEvent
isdeath(e::Event) = isdeath(kind(e))


text:
#
# Optimizations:
#  - replace leftmap, rightmap with sorted vectors
#

###########################
# SiteOp                  #
###########################

struct SiteOp{N}
  name::String
  site::NTuple{N,Int}
  params::NamedTuple
end
# Change NamedTuple() to (;) when we drop older Julia versions
SiteOp(name::String, site::Tuple) = SiteOp(name, site, NamedTuple())
SiteOp(name::String, site::Int...) = SiteOp(name, site)
function SiteOp(name::String, site_params::Union{Int,NamedTuple}...)
  return SiteOp(name, Base.front(site_params), last(site_params))
end
SiteOp(name::String, params::NamedTuple, site::Tuple) = SiteOp(name, site, params)
SiteOp(name::String, params::NamedTuple, site::Int...) = SiteOp(name, site, params)

convert(::Type{SiteOp}, op::Pair{String,Int}) = SiteOp(first(op), last(op))

name(s::SiteOp) = s.name
site(s::SiteOp) = only(s.site)
sites(s::SiteOp) = s.site
params(s::SiteOp) = s.params

site_or_sites(s::SiteOp{1}) = site(s)
site_or_sites(s::SiteOp) = sites(s)

string_site_or_sites(s::SiteOp{1}) = string(site(s))
string_site_or_sites(s::SiteOp) = string(sites(s))[2:(end - 1)]

show(io::IO, s::SiteOp) = print(io, "\"$(name(s))\"($(string_site_or_sites(s)))")

(s1::SiteOp == s2::SiteOp) = (s1.site == s2.site && s1.name == s2.name)

function isless(s1::SiteOp, s2::SiteOp)
  if site(s1) != site(s2)
    return site(s1) < site(s2)
  end
  return name(s1) < name(s2)
end

###########################
# OpTerm                  #
###########################

const OpTerm = Vector{SiteOp}

function (o1::OpTerm == o2::OpTerm)
  (length(o1) == length(o2)) || return false
  @inbounds for n in 1:length(o1)
    (o1[n] != o2[n]) && return false
  end
  return true
end

function isless(o1::OpTerm, o2::OpTerm)
  if length(o1) != length(o2)
    return length(o1) < length(o2)
  end
  for n in 1:length(o1)
    if o1[n] != o2[n]
      return (o1[n] < o2[n])
    end
  end
  return false
end

mult(t1::OpTerm, t2::OpTerm) = isempty(t2) ? t1 : vcat(t1, t2)

function isfermionic(t::OpTerm, sites)::Bool
  p = +1
  for op in t
    if has_fermion_string(name(op), sites[site(op)])
      p *= -1
    end
  end
  return (p == -1)
end

###########################
# MPOTerm                 #
###########################

mutable struct MPOTerm
  coef::ComplexF64
  ops::OpTerm
end
coef(op::MPOTerm) = op.coef
ops(op::MPOTerm) = op.ops

copy(t::MPOTerm) = MPOTerm(coef(t), copy(ops(t)))

function (t1::MPOTerm == t2::MPOTerm)
  return coef(t1) ≈ coef(t2) && ops(t1) == ops(t2)
end

function isless(t1::MPOTerm, t2::MPOTerm)
  if ops(t1) == ops(t2)
    if coef(t1) ≈ coef(t2)
      return false
    else
      ct1 = coef(t1)
      ct2 = coef(t2)
      #"lexicographic" ordering on  complex numbers
      return real(ct1) < real(ct2) || (real(ct1) ≈ real(ct2) && imag(ct1) < imag(ct2))
    end
  end
  return ops(t1) < ops(t2)
end

function MPOTerm(c::Number, op1::String, ops_rest...)
  ops = (op1, ops_rest...)
  starts = findall(x -> x isa String, ops)
  N = length(starts)
  vop = OpTerm(undef, N)
  for n in 1:N
    start = starts[n]
    stop = (n == N) ? lastindex(ops) : (starts[n + 1] - 1)
    vop[n] = SiteOp(ops[start:stop]...)
  end
  return MPOTerm(c, vop)
end

function MPOTerm(op1::String, ops...)
  return MPOTerm(one(Float64), op1, ops...)
end

function MPOTerm(ops::Vector{Pair{String,Int}})
  return MPOTerm(Iterators.flatten(ops)...)
end

function Base.show(io::IO, op::MPOTerm)
  c = coef(op)
  if iszero(imag(c))
    print(io, "$(real(c)) ")
  elseif iszero(real(c))
    print(io, "$(imag(c))im ")
  else
    print(io, "($c) ")
  end
  for o in ops(op)
    print(io, "\"$(name(o))\"($(string_site_or_sites(o))) ")
  end
end

############################
## OpSum                 #
############################

"""
An `OpSum` represents a sum of operator
terms.

Often it is used to create matrix
product operator (`MPO`) approximation
of the sum of the terms in the `OpSum` oject.
Each term is a product of local operators
specified by names such as `"Sz"` or `"N"`,
times an optional coefficient which
can be real or complex.

Which local operator names are available
is determined by the function `op`
associated with the `TagType` defined by
special Index tags, such as `"S=1/2"`, `"S=1"`,
`"Fermion"`, and `"Electron"`.
"""
mutable struct OpSum
  data::Vector{MPOTerm}
  OpSum(terms::Vector{MPOTerm}) = new(terms)
end

length(os::OpSum) = length(data(os))
getindex(os::OpSum, I...) = data(os)[I...]

const AutoMPO = OpSum

"""
    OpSum()

Construct an empty `OpSum`.
"""
OpSum() = OpSum(Vector{MPOTerm}())

data(ampo::OpSum) = ampo.data
setdata!(ampo::OpSum, ndata) = (ampo.data = ndata)

push!(ampo::OpSum, term) = push!(data(ampo), term)

Base.:(==)(ampo1::OpSum, ampo2::OpSum) = data(ampo1) == data(ampo2)

Base.copy(ampo::OpSum) = OpSum(copy(data(ampo)))

function Base.deepcopy(ampo::OpSum)
  return OpSum(map(copy, data(ampo)))
end

Base.size(ampo::OpSum) = size(data(ampo))

"""
    add!(ampo::OpSum,
         op1::String, i1::Int)

    add!(ampo::OpSum,
         coef::Number,
         op1::String, i1::Int)

    add!(ampo::OpSum,
         op1::String, i1::Int,
         op2::String, i2::Int,
         ops...)

    add!(ampo::OpSum,
         coef::Number,
         op1::String, i1::Int,
         op2::String, i2::Int,
         ops...)

    +(ampo:OpSum, term::Tuple)

Add a single- or multi-site operator
term to the OpSum `ampo`. Each operator
is specified by a name (String) and a
site number (Int). The second version
accepts a real or complex coefficient.

The `+` operator version of this function
accepts a tuple with entries either
(String,Int,String,Int,...) or
(Number,String,Int,String,Int,...)
where these tuple values are the same
as valid inputs to the `add!` function.
For inputting a very large number of
terms (tuples) to an OpSum, consider
using the broadcasted operator `.+=`
which avoids reallocating the OpSum
after each addition.

# Examples
```julia
ampo = OpSum()

add!(ampo,"Sz",2,"Sz",3)

ampo += ("Sz",3,"Sz",4)

ampo += (0.5,"S+",4,"S-",5)

ampo .+= (0.5,"S+",5,"S-",6)
```
"""
add!(os::OpSum, t::MPOTerm) = push!(os, t)

add!(os::OpSum, args...) = add!(os, MPOTerm(args...))

"""
    subtract!(ampo::OpSum,
              op1::String, i1::Int,
              op2::String, i2::Int,
              ops...)

    subtract!(ampo::OpSum,
              coef::Number,
              op1::String, i1::Int,
              op2::String, i2::Int,
              ops...)

Subtract a multi-site operator term
from the OpSum `ampo`. Each operator
is specified by a name (String) and a
site number (Int). The second version
accepts a real or complex coefficient.
"""
subtract!(os::OpSum, args...) = add!(os, -MPOTerm(args...))

-(t::MPOTerm) = MPOTerm(-coef(t), ops(t))

function (ampo::OpSum + term::Tuple)
  ampo_plus_term = copy(ampo)
  add!(ampo_plus_term, term...)
  return ampo_plus_term
end

function (ampo::OpSum + term::Vector{Pair{String,Int64}})
  ampo_plus_term = copy(ampo)
  add!(ampo_plus_term, term)
  return ampo_plus_term
end

function (ampo::OpSum - term::Tuple)
  ampo_plus_term = copy(ampo)
  subtract!(ampo_plus_term, term...)
  return ampo_plus_term
end

#
# ampo .+= ("Sz",1) syntax using broadcasting
#

struct OpSumStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:OpSum}) = OpSumStyle()

struct OpSumAddTermStyle <: Broadcast.BroadcastStyle end

Base.broadcastable(ampo::OpSum) = ampo

Base.BroadcastStyle(::OpSumStyle, ::Broadcast.Style{Tuple}) = OpSumAddTermStyle()

Broadcast.instantiate(bc::Broadcast.Broadcasted{OpSumAddTermStyle}) = bc

function Base.copyto!(ampo, bc::Broadcast.Broadcasted{OpSumAddTermStyle,<:Any,typeof(+)})
  add!(ampo, bc.args[2]...)
  return ampo
end

#
# ampo .-= ("Sz",1) syntax using broadcasting
#

function Base.copyto!(ampo, bc::Broadcast.Broadcasted{OpSumAddTermStyle,<:Any,typeof(-)})
  subtract!(ampo, bc.args[2]...)
  return ampo
end

function Base.show(io::IO, ampo::OpSum)
  println(io, "OpSum:")
  for term in data(ampo)
    println(io, "  $term")
  end
end

##################################
# MatElem (simple sparse matrix) #
##################################

struct MatElem{T}
  row::Int
  col::Int
  val::T
end

#function Base.show(io::IO,m::MatElem)
#  print(io,"($(m.row),$(m.col),$(m.val))")
#end

function toMatrix(els::Vector{MatElem{T}})::Matrix{T} where {T}
  nr = 0
  nc = 0
  for el in els
    nr = max(nr, el.row)
    nc = max(nc, el.col)
  end
  M = zeros(T, nr, nc)
  for el in els
    M[el.row, el.col] = el.val
  end
  return M
end

function Base.:(==)(m1::MatElem{T}, m2::MatElem{T})::Bool where {T}
  return (m1.row == m2.row && m1.col == m2.col && m1.val == m2.val)
end

function Base.isless(m1::MatElem{T}, m2::MatElem{T})::Bool where {T}
  if m1.row != m2.row
    return m1.row < m2.row
  elseif m1.col != m2.col
    return m1.col < m2.col
  end
  return m1.val < m2.val
end

struct QNMatElem{T}
  rowqn::QN
  colqn::QN
  row::Int
  col::Int
  val::T
end

function Base.:(==)(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T}
  return (
    m1.row == m2.row &&
    m1.col == m2.col &&
    m1.val == m2.val &&
    m1.rowqn == m2.rowqn &&
    m1.colqn == m2.colqn
  )
end

function Base.isless(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T}
  if m1.rowqn != m2.rowqn
    return m1.rowqn < m2.rowqn
  elseif m1.colqn != m2.colqn
    return m1.colqn < m2.colqn
  elseif m1.row != m2.row
    return m1.row < m2.row
  elseif m1.col != m2.col
    return m1.col < m2.col
  end
  return m1.val < m2.val
end

isempty(op_qn::Pair{OpTerm,QN}) = isempty(op_qn.first)

# the key type is OpTerm for the dense case
# and is Pair{OpTerm,QN} for the QN conserving case
function posInLink!(linkmap::Dict{K,Int}, k::K)::Int where {K}
  isempty(k) && return -1
  pos = get(linkmap, k, -1)
  if pos == -1
    pos = length(linkmap) + 1
    linkmap[k] = pos
  end
  return pos
end

function determineValType(terms::Vector{MPOTerm})
  for t in terms
    (!isreal(coef(t))) && return ComplexF64
  end
  return Float64
end

function computeSiteProd(sites, ops::OpTerm)::ITensor
  i = site(ops[1])
  T = op(sites[i], ops[1].name)
  for j in 2:length(ops)
    (site(ops[j]) != i) && error("Mismatch of site number in computeSiteProd")
    opj = op(sites[i], ops[j].name)
    T = product(T, opj)
  end
  return T
end

function remove_dups!(v::Vector{T}) where {T}
  N = length(v)
  (N == 0) && return nothing
  sort!(v)
  n = 1
  u = 2
  while u <= N
    while u < N && v[u] == v[n]
      u += 1
    end
    if v[u] != v[n]
      v[n + 1] = v[u]
      n += 1
    end
    u += 1
  end
  resize!(v, n)
  return nothing
end #remove_dups!

function svdMPO(ampo::OpSum, sites; kwargs...)::MPO
  mindim::Int = get(kwargs, :mindim, 1)
  maxdim::Int = get(kwargs, :maxdim, 10000)
  cutoff::Float64 = get(kwargs, :cutoff, 1E-15)

  N = length(sites)

  ValType = determineValType(data(ampo))

  Vs = [Matrix{ValType}(undef, 1, 1) for n in 1:N]
  tempMPO = [MatElem{MPOTerm}[] for n in 1:N]

  crosses_bond(t::MPOTerm, n::Int) = (site(ops(t)[1]) <= n <= site(ops(t)[end]))

  rightmap = Dict{OpTerm,Int}()
  next_rightmap = Dict{OpTerm,Int}()

  for n in 1:N
    leftbond_coefs = MatElem{ValType}[]

    leftmap = Dict{OpTerm,Int}()
    for term in data(ampo)
      crosses_bond(term, n) || continue

      left::OpTerm = filter(t -> (site(t) < n), ops(term))
      onsite::OpTerm = filter(t -> (site(t) == n), ops(term))
      right::OpTerm = filter(t -> (site(t) > n), ops(term))

      bond_row = -1
      bond_col = -1
      if !isempty(left)
        bond_row = posInLink!(leftmap, left)
        bond_col = posInLink!(rightmap, mult(onsite, right))
        bond_coef = convert(ValType, coef(term))
        push!(leftbond_coefs, MatElem(bond_row, bond_col, bond_coef))
      end

      A_row = bond_col
      A_col = posInLink!(next_rightmap, right)
      site_coef = 1.0 + 0.0im
      if A_row == -1
        site_coef = coef(term)
      end
      if isempty(onsite)
        if !using_auto_fermion() && isfermionic(right, sites)
          push!(onsite, SiteOp("F", n))
        else
          push!(onsite, SiteOp("Id", n))
        end
      end
      el = MatElem(A_row, A_col, MPOTerm(site_coef, onsite))
      push!(tempMPO[n], el)
    end
    rightmap = next_rightmap
    next_rightmap = Dict{OpTerm,Int}()

    remove_dups!(tempMPO[n])

    if n > 1 && !isempty(leftbond_coefs)
      M = toMatrix(leftbond_coefs)
      U, S, V = svd(M)
      P = S .^ 2
      truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
      tdim = length(P)
      nc = size(M, 2)
      Vs[n - 1] = Matrix{ValType}(V[1:nc, 1:tdim])
    end
  end

  llinks = Vector{Index{Int}}(undef, N + 1)
  llinks[1] = Index(2, "Link,l=0")

  H = MPO(sites)

  # Constants which define MPO start/end scheme
  rowShift = 2
  colShift = 2
  startState = 2
  endState = 1

  for n in 1:N
    VL = Matrix{ValType}(undef, 1, 1)
    if n > 1
      VL = Vs[n - 1]
    end
    VR = Vs[n]
    tdim = size(VR, 2)

    llinks[n + 1] = Index(2 + tdim, "Link,l=$n")

    finalMPO = Dict{OpTerm,Matrix{ValType}}()

    ll = llinks[n]
    rl = llinks[n + 1]

    idTerm = [SiteOp("Id", n)]
    finalMPO[idTerm] = zeros(ValType, dim(ll), dim(rl))
    idM = finalMPO[idTerm]
    idM[1, 1] = 1.0
    idM[2, 2] = 1.0

    defaultMat() = zeros(ValType, dim(ll), dim(rl))

    for el in tempMPO[n]
      A_row = el.row
      A_col = el.col
      t = el.val
      (abs(coef(t)) > eps()) || continue

      M = get!(finalMPO, ops(t), defaultMat())

      ct = convert(ValType, coef(t))
      if A_row == -1 && A_col == -1 #onsite term
        M[startState, endState] += ct
      elseif A_row == -1 #term starting on site n
        for c in 1:size(VR, 2)
          z = ct * VR[A_col, c]
          M[startState, colShift + c] += z
        end
      elseif A_col == -1 #term ending on site n
        for r in 1:size(VL, 2)
          z = ct * conj(VL[A_row, r])
          M[rowShift + r, endState] += z
        end
      else
        for r in 1:size(VL, 2), c in 1:size(VR, 2)
          z = ct * conj(VL[A_row, r]) * VR[A_col, c]
          M[rowShift + r, colShift + c] += z
        end
      end
    end

    s = sites[n]
    H[n] = ITensor()
    for (op, M) in finalMPO
      T = itensor(M, ll, rl)
      H[n] += T * computeSiteProd(sites, op)
    end
  end

  L = ITensor(llinks[1])
  L[startState] = 1.0

  R = ITensor(llinks[N + 1])
  R[endState] = 1.0

  H[1] *= L
  H[N] *= R

  return H
end #svdMPO

function qn_svdMPO(ampo::OpSum, sites; kwargs...)::MPO
  mindim::Int = get(kwargs, :mindim, 1)
  maxdim::Int = get(kwargs, :maxdim, 10000)
  cutoff::Float64 = get(kwargs, :cutoff, 1E-15)

  N = length(sites)

  ValType = determineValType(data(ampo))

  Vs = [Dict{QN,Matrix{ValType}}() for n in 1:(N + 1)]
  tempMPO = [QNMatElem{MPOTerm}[] for n in 1:N]

  crosses_bond(t::MPOTerm, n::Int) = (site(ops(t)[1]) <= n <= site(ops(t)[end]))

  rightmap = Dict{Pair{OpTerm,QN},Int}()
  next_rightmap = Dict{Pair{OpTerm,QN},Int}()

  # A cache of the ITensor operators on a certain site
  # of a certain type
  op_cache = Dict{Pair{String,Int},ITensor}()

  for n in 1:N
    leftbond_coefs = Dict{QN,Vector{MatElem{ValType}}}()

    leftmap = Dict{Pair{OpTerm,QN},Int}()
    for term in data(ampo)
      crosses_bond(term, n) || continue

      left::OpTerm = filter(t -> (site(t) < n), ops(term))
      onsite::OpTerm = filter(t -> (site(t) == n), ops(term))
      right::OpTerm = filter(t -> (site(t) > n), ops(term))

      function calcQN(term::OpTerm)
        q = QN()
        for st in term
          op_tensor = get(op_cache, name(st) => site(st), nothing)
          if op_tensor === nothing
            op_tensor = op(sites[site(st)], name(st))
            op_cache[name(st) => site(st)] = op_tensor
          end
          q -= flux(op_tensor)
        end
        return q
      end
      lqn = calcQN(left)
      sqn = calcQN(onsite)

      bond_row = -1
      bond_col = -1
      if !isempty(left)
        bond_row = posInLink!(leftmap, left => lqn)
        bond_col = posInLink!(rightmap, mult(onsite, right) => lqn)
        bond_coef = convert(ValType, coef(term))
        q_leftbond_coefs = get!(leftbond_coefs, lqn, MatElem{ValType}[])
        push!(q_leftbond_coefs, MatElem(bond_row, bond_col, bond_coef))
      end

      rqn = sqn + lqn
      A_row = bond_col
      A_col = posInLink!(next_rightmap, right => rqn)
      site_coef = 1.0 + 0.0im
      if A_row == -1
        site_coef = coef(term)
      end
      if isempty(onsite)
        if !using_auto_fermion() && isfermionic(right, sites)
          push!(onsite, SiteOp("F", n))
        else
          push!(onsite, SiteOp("Id", n))
        end
      end
      el = QNMatElem(lqn, rqn, A_row, A_col, MPOTerm(site_coef, onsite))
      push!(tempMPO[n], el)
    end
    rightmap = next_rightmap
    next_rightmap = Dict{Pair{OpTerm,QN},Int}()

    remove_dups!(tempMPO[n])

    if n > 1 && !isempty(leftbond_coefs)
      for (q, mat) in leftbond_coefs
        M = toMatrix(mat)
        U, S, V = svd(M)
        P = S .^ 2
        truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
        tdim = length(P)
        nc = size(M, 2)
        Vs[n][q] = Matrix{ValType}(V[1:nc, 1:tdim])
      end
    end
  end

  #
  # Make MPO link indices
  #
  d0 = 2
  llinks = Vector{QNIndex}(undef, N + 1)
  # Set dir=In for fermionic ordering, avoid arrow sign
  # <fermions>:
  linkdir = using_auto_fermion() ? In : Out
  llinks[1] = Index(QN() => d0; tags="Link,l=0", dir=linkdir)
  for n in 1:N
    qi = Vector{Pair{QN,Int}}()
    if !haskey(Vs[n + 1], QN())
      # Make sure QN=zero is first in list of sectors
      push!(qi, QN() => d0)
    end
    for (q, Vq) in Vs[n + 1]
      cols = size(Vq, 2)
      if q == QN()
        # Make sure QN=zero is first in list of sectors
        insert!(qi, 1, q => d0 + cols)
      else
        if using_auto_fermion() # <fermions>
          push!(qi, (-q) => cols)
        else
          push!(qi, q => cols)
        end
      end
    end
    # Set dir=In for fermionic ordering, avoid arrow sign
    # <fermions>:
    llinks[n + 1] = Index(qi...; tags="Link,l=$n", dir=linkdir)
  end

  H = MPO(N)

  # Constants which define MPO start/end scheme
  startState = 2
  endState = 1

  for n in 1:N
    finalMPO = Dict{Tuple{QN,OpTerm},Matrix{ValType}}()

    ll = llinks[n]
    rl = llinks[n + 1]

    function defaultMat(ll, rl, lqn, rqn)
      #ldim = qnblockdim(ll,lqn)
      #rdim = qnblockdim(rl,rqn)
      ldim = blockdim(ll, lqn)
      rdim = blockdim(rl, rqn)
      return zeros(ValType, ldim, rdim)
    end

    idTerm = [SiteOp("Id", n)]
    finalMPO[(QN(), idTerm)] = defaultMat(ll, rl, QN(), QN())
    idM = finalMPO[(QN(), idTerm)]
    idM[1, 1] = 1.0
    idM[2, 2] = 1.0

    for el in tempMPO[n]
      t = el.val
      (abs(coef(t)) > eps()) || continue
      A_row = el.row
      A_col = el.col

      M = get!(finalMPO, (el.rowqn, ops(t)), defaultMat(ll, rl, el.rowqn, el.colqn))

      # rowShift and colShift account for
      # special entries in the zero-QN sector
      # of the MPO
      rowShift = (el.rowqn == QN()) ? 2 : 0
      colShift = (el.colqn == QN()) ? 2 : 0

      ct = convert(ValType, coef(t))
      if A_row == -1 && A_col == -1 #onsite term
        M[startState, endState] += ct
      elseif A_row == -1 #term starting on site n
        VR = Vs[n + 1][el.colqn]
        for c in 1:size(VR, 2)
          z = ct * VR[A_col, c]
          M[startState, colShift + c] += z
        end
      elseif A_col == -1 #term ending on site n
        VL = Vs[n][el.rowqn]
        for r in 1:size(VL, 2)
          z = ct * conj(VL[A_row, r])
          M[rowShift + r, endState] += z
        end
      else
        VL = Vs[n][el.rowqn]
        VR = Vs[n + 1][el.colqn]
        for r in 1:size(VL, 2), c in 1:size(VR, 2)
          z = ct * conj(VL[A_row, r]) * VR[A_col, c]
          M[rowShift + r, colShift + c] += z
        end
      end
    end

    s = sites[n]
    H[n] = ITensor()
    for (q_op, M) in finalMPO
      op_prod = q_op[2]
      Op = computeSiteProd(sites, op_prod)

      rq = q_op[1]
      sq = flux(Op)
      cq = rq - sq

      if using_auto_fermion()
        # <fermions>:
        # MPO is defined with Index order
        # of (rl,s[n]',s[n],cl) where rl = row link, cl = col link
        # so compute sign that would result by permuting cl from
        # second position to last position:
        if fparity(sq) == 1 && fparity(cq) == 1
          Op .*= -1
        end
      end

      rn = qnblocknum(ll, rq)
      cn = qnblocknum(rl, cq)

      #TODO: wrap following 3 lines into a function
      _block = Block(rn, cn)
      T = BlockSparseTensor(ValType, [_block], (dag(ll), rl))
      #blockview(T, _block) .= M
      T[_block] .= M

      IT = itensor(T)
      H[n] += IT * Op
    end
  end

  L = ITensor(llinks[1])
  L[startState] = 1.0

  R = ITensor(dag(llinks[N + 1]))
  R[endState] = 1.0

  H[1] *= L
  H[N] *= R

  return H
end #qn_svdMPO

function sorteachterm!(ampo::OpSum, sites)
  ampo = copy(ampo)
  isless_site(o1::SiteOp, o2::SiteOp) = site(o1) < site(o2)
  N = length(sites)
  for t in data(ampo)
    Nt = length(t.ops)
    prevsite = N + 1 #keep track of whether we are switching
    #to a new site to make sure F string
    #is only placed at most once for each site

    # Sort operators in t by site order,
    # and keep the permutation used, perm, for analysis below
    perm = Vector{Int}(undef, Nt)
    sortperm!(perm, t.ops; alg=InsertionSort, lt=isless_site)

    t.ops = t.ops[perm]

    # Identify fermionic operators,
    # zeroing perm for bosonic operators,
    # and inserting string "F" operators
    parity = +1
    for n in Nt:-1:1
      currsite = site(t.ops[n])
      fermionic = has_fermion_string(name(t.ops[n]), sites[site(t.ops[n])])
      if !using_auto_fermion() && (parity == -1) && (currsite < prevsite)
        # Put local piece of Jordan-Wigner string emanating
        # from fermionic operators to the right
        # (Remaining F operators will be put in by svdMPO)
        t.ops[n] = SiteOp("$(name(t.ops[n]))*F", site(t.ops[n]))
      end
      prevsite = currsite

      if fermionic
        parity = -parity
      else
        # Ignore bosonic operators in perm
        # by zeroing corresponding entries
        perm[n] = 0
      end
    end
    if parity == -1
      error("Parity-odd fermionic terms not yet supported by AutoMPO")
    end

    # Keep only fermionic op positions (non-zero entries)
    filter!(!iszero, perm)
    # and account for anti-commuting, fermionic operators
    # during above sort; put resulting sign into coef
    t.coef *= parity_sign(perm)
  end
  return ampo
end

function sortmergeterms!(ampo::OpSum)
  sort!(data(ampo))

  # Merge (add) terms with same operators
  da = data(ampo)
  ndata = MPOTerm[]
  last_term = copy(da[1])
  for n in 2:length(da)
    if ops(da[n]) == ops(last_term)
      last_term.coef += coef(da[n])
    else
      push!(ndata, last_term)
      last_term = copy(da[n])
    end
  end
  push!(ndata, last_term)

  setdata!(ampo, ndata)
  return ampo
end

"""
    MPO(ampo::OpSum,sites::Vector{<:Index};kwargs...)

Convert an OpSum object `ampo` to an
MPO, with indices given by `sites`. The
resulting MPO will have the indices
`sites[1], sites[1]', sites[2], sites[2]'`
etc. The conversion is done by an algorithm
that compresses the MPO resulting from adding
the OpSum terms together, often achieving
the minimum possible bond dimension.

# Examples
```julia
ampo = OpSum()
ampo += ("Sz",1,"Sz",2)
ampo += ("Sz",2,"Sz",3)
ampo += ("Sz",3,"Sz",4)

sites = siteinds("S=1/2",4)
H = MPO(ampo,sites)
```
"""
function MPO(ampo::OpSum, sites::Vector{<:Index}; kwargs...)::MPO
  length(data(ampo)) == 0 && error("OpSum has no terms")

  ampo = deepcopy(ampo)
  sorteachterm!(ampo, sites)
  sortmergeterms!(ampo)

  if hasqns(sites[1])
    return qn_svdMPO(ampo, sites; kwargs...)
  end
  return svdMPO(ampo, sites; kwargs...)
end


text:
module TestMOIWrapper

using Test
using MathOptInterface
import Cbc

const MOI = MathOptInterface

const OPTIMIZER = Cbc.Optimizer()
MOI.set(OPTIMIZER, MOI.Silent(), true)

function test_SolverName()
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "COIN Branch-and-Cut (Cbc)"
end

function test_supports_default_copy_to()
    @test !MOI.Utilities.supports_allocate_load(OPTIMIZER, false)
    @test !MOI.Utilities.supports_allocate_load(OPTIMIZER, true)
    @test !MOI.Utilities.supports_default_copy_to(OPTIMIZER, false)
    @test !MOI.Utilities.supports_default_copy_to(OPTIMIZER, true)
end

const CACHED = MOI.Utilities.CachingOptimizer(
    MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    OPTIMIZER,
)
const BRIDGED = MOI.Bridges.full_bridge_optimizer(CACHED, Float64)

const CONFIG = MOI.Test.TestConfig(duals = false, infeas_certificates = false)

function test_basic_constraint_tests()
    return MOI.Test.basic_constraint_tests(CACHED, CONFIG)
end

function test_Unit()
    return MOI.Test.unittest(
        BRIDGED,
        CONFIG,
        [
            # TODO(odow): implement attributes:
            "number_threads",

            # INFEASIBLE_OR_UNBOUNDED instead of DUAL_INFEASIBLE
            "solve_unbounded_model",

            # No quadratics
            "delete_soc_variables",
            "solve_qcp_edge_cases",
            "solve_qp_edge_cases",
        ],
    )
end

function test_solvername()
    @test MOI.get(CACHED, MOI.SolverName()) == "COIN Branch-and-Cut (Cbc)"
end

function test_default_objective_test()
    return MOI.Test.default_objective_test(CACHED)
end

function test_default_status_test()
    return MOI.Test.default_status_test(CACHED)
end

function test_nametest()
    return MOI.Test.nametest(CACHED)
end

function test_validtest()
    return MOI.Test.validtest(CACHED)
end

function test_emptytest()
    return MOI.Test.emptytest(BRIDGED)
end

function test_orderedindicestest()
    return MOI.Test.orderedindicestest(CACHED)
end

function test_ContinuousLinear()
    return MOI.Test.contlineartest(BRIDGED, CONFIG)
end

function test_IntegerLinear()
    return MOI.Test.intlineartest(
        BRIDGED,
        CONFIG,
        [
            # Cbc does not support indicator constraints.
            "indicator1",
            "indicator2",
            "indicator3",
            "indicator4",
            # SOS issues
            "int2",
        ],
    )
end

function test_Testparams()
    # Note: we generate a non-trivial problem to ensure that Cbc struggles to
    # find a solution at the root node.
    knapsack_model = MOI.Utilities.Model{Float64}()
    N = 100
    x = MOI.add_variables(knapsack_model, N)
    MOI.add_constraint.(knapsack_model, MOI.SingleVariable.(x), MOI.ZeroOne())
    MOI.add_constraint(
        knapsack_model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1 + sin(i) for i in 1:N], x),
            0.0,
        ),
        MOI.LessThan(10.0),
    )
    MOI.set(
        knapsack_model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([cos(i) for i in 1:N], x),
            0.0,
        ),
    )
    model = Cbc.Optimizer(
        maxSol = 1,
        presolve = "off",
        cuts = "off",
        heur = "off",
        logLevel = 0,
    )
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    MOI.empty!(model)
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
end

function test_TestPrimalStatus()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.add_constraint(model, x, MOI.LessThan(0.0))
    cbc = Cbc.Optimizer()
    MOI.copy_to(cbc, model)
    MOI.optimize!(cbc)
    return MOI.get(cbc, MOI.PrimalStatus()) == MOI.NO_SOLUTION
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$name", "test_")
            continue
        end
        @testset "$(name)" begin
            getfield(@__MODULE__, name)()
        end
    end
    return
end

end

TestMOIWrapper.runtests()


text:
include("scaling_nthreads_target.jl")
using Setfield: @set!

function manual_symmetrize!(A, B)
    @argcheck size(A) == size(B)
    for j in axes(A, 2), i in axes(A, 1)
        @inbounds A[i, j] = (B[i, j] + B[j, i]) / 2
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    let benchmarks = BENCHMARKS

        # sort: use `alg = Base.QuickSort`
        @set! benchmarks[:sort].run = (p, r) -> sort!(r; alg = p.alg.value)
        @set! benchmarks[:sort].paramaxes.alg =
            [(label = "Base.QuickSort", value = Base.QuickSort)]
        @set! benchmarks[:sort].paramaxes.basesize = [nothing]
        @set! benchmarks[:sort].paramaxes.smallsize = [nothing]

        # foreach_symmetrize: use manual implementation
        @set! benchmarks[:foreach_symmetrize].run = (p, r) -> manual_symmetrize!(r.A, r.B)

        @set! benchmarks[:sum_sin].run = (p, r) -> sum(sin, r)
        @set! benchmarks[:findfirst].run = (p, r) -> findfirst(==(-1), r)
        @set! benchmarks[:unique].run = (p, r) -> unique(r)

        main(benchmark_definitions = benchmarks, tags = ["nthreads-baseline"])
    end
end