text:
using Clustering
using Distances
using FileIO
using KernelFunctions
using MLDatasets
using PrototypesCriticisms
using Random

"""
    saveasimage(d, ids; outdir="out")
    saveasimage(d; kwargs...)

Save the data instances `ids` of the dataset `d` as image files.

If `ids` is unspecified, all data instances of the dataset are saved as image files.
Assumes that each data instance actually represents an image.

# Keyword arguments
- `outdir`: the location where the image files are saved.
- `ext`: the file extension of the image files.

"""
function saveasimage(d, ids; outdir="out", ext="png")
    mkpath(outdir)
    for i in ids
        image = convert2image(d, i)
        classname = d.targets[i]
        if haskey(d.metadata, "class_names")
            classnames = d.metadata["class_names"]
            classname = classnames[d.targets[i] + 1]
        end
        save("$outdir/$("0" ^ (ndigits(length(d)) - ndigits(i)))$(i)_$classname.$ext", image)
    end
end
saveasimage(d; kwargs...) = saveasimage(d, 1:length(d); kwargs...)

"""
    printclusters(clusters; headline="", indent=2)

Print all instances of every cluster in `clusters`.

# Keyword arguments
- `headline`: the headline that is printed before the output of the clusters.
- `indent`: the number of spaces to indent the output of every cluster.
"""
function printclusters(clusters; headline="", indent=2)
    println(headline)
    for (i, instances) in enumerate(clusters)
        println("$(" " ^ indent)Cluster $(lpad(i, ndigits(length(clusters)))): $instances")
    end
end

"""
    evaluate(c, ys)

Evaluate the clustering `c` using the labels `ys`.
"""
function evaluate(c, ys)
    ris = randindex(c, ys)
    println("Rand-related indices:")
    println("  Adjusted rand index: ", ris[1])
    println("  P(agree): ", ris[2])
    println("  P(disagree): ", ris[3])
    println("  P(agree) - P(disagree): ", ris[4])
    println("Variation of information: ", Clustering.varinfo(c, ys))
    println("V-measure: ", vmeasure(c, ys))
    println("Mutual information: ", mutualinfo(c, ys))
end

"""
    main()

Run the example.
"""
function main()
    Random.seed!(42)

    # Load dataset and the corresponding pre-trained image embeddings
    @info "Load data..."
    dataset = CIFAR10()
    embedding = load("cifar-10_train_embedding_vgg-19.jld2", "embedding")'

    # Save whole data set as images
    @info "Save images..."
    saveasimage(dataset, outdir="out/cifar-10")

    # Set main program parameters
    k = 10 # Number of clusters to compute
    p = 5 # Number of prototypes to find
    c = 5 # Number of criticisms to find

    # Cluster data using k-medoids
    @info "Cluster data..."
    clustering = kmedoids(pairwise(Euclidean(), embedding), k)

    # Evaluate the clustering
    @info "Evaluate clustering..."
    targets = dataset.targets .+ 1 # Cluster indices must start with 1 instead of 0
    evaluate(clustering, targets)

    # Naive prototypes and criticisms for the clustering
    @info "Compute naive prototypes and criticisms..."
    protoids = prototypes(clustering, p)
    critids = criticisms(clustering, c)
    printclusters(protoids, headline="Prototype ids:")
    printclusters(critids, headline="Criticism ids:")
    for i = 1:k
        saveasimage(dataset, protoids[i], outdir="out/cifar-10/_protos/naive/$i")
        saveasimage(dataset, critids[i], outdir="out/cifar-10/_crits/naive/$i")
    end

    # Prototypes and criticisms for the clustering using MMD-critic
    @info "Compute MMD-critic prototypes and criticisms..."
    ys = assignments(clustering)
    clusters = Vector{Vector{Int}}(undef, k)
    for i = 1:k
        clusters[i] = findall(ys .== i)
    end
    protoids = Vector{Vector{Int}}(undef, k)
    critids = Vector{Vector{Int}}(undef, k)
    for i = 1:k
        subembedding = view(embedding, :, clusters[i])
        kernel = with_lengthscale(RBFKernel(), sqrt(size(embedding, 1)))
        clusterpids = prototypes(subembedding, kernel, p)
        protoids[i] = clusters[i][clusterpids]
        clustercids = criticisms(subembedding, kernel, clusterpids, c)
        critids[i] = clusters[i][clustercids]
    end
    printclusters(protoids, headline="Prototype ids:")
    printclusters(critids, headline="Criticism ids:")
    for i = 1:k
        saveasimage(dataset, protoids[i], outdir="out/cifar-10/_protos/mmd-critic/$i")
        saveasimage(dataset, critids[i], outdir="out/cifar-10/_crits/mmd-critic/$i")
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
function linearSplineInterpolation(fineT,fineX,T)
    if length(T)>1
        return map(T -> linearSplineInterpolation(fineT,fineX,T),T)
    end
    i = findfirst(fineT.>=T)
    if isnothing(i)
        return fineX[:,end]
    end
    if i<=1
        return fineX[:,1]
    end
    return fineX[:,i-1]+(T-fineT[i-1])/(fineT[i]-fineT[i-1])*(fineX[:,i]-fineX[:,i-1])
end

println("Finished loading HamiltonInterpolate")

ext: jl

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("utilities.jl")

_functionize_error(b, bridge_type, func, name) = ErrorException(
    "Need to apply a `$bridge_type` to a `$func` $name because the" *
    " variable is bridged but $name bridges are not supported by" *
    " `$(typeof(b))`."
)

_lazy_functionize_error(bridge_type, func, name) = ErrorException(
    "Need to apply a `$bridge_type` to a `$func` $name because the variable is" *
    " bridged but no such $name bridge type was added. Add one" *
    " with `add_bridge`."
)

MOIU.@model(StandardLPModel,
            (), (MOI.EqualTo, MOIT.UnknownScalarSet), (MOI.Nonnegatives,), (),
            (), (MOI.ScalarAffineFunction,), (MOI.VectorOfVariables,), ())
MOI.supports_constraint(::StandardLPModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{T}}) where {T} = false
MOI.supports_constraint(::StandardLPModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{T}}) where {T} = false
MOI.supports_constraint(::StandardLPModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.EqualTo{T}}) where {T} = false
MOI.supports_constraint(::StandardLPModel, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Reals}) = false
MOI.supports(::StandardLPModel{T}, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}) where {T} = false
MOI.supports(::StandardLPModel, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = false

@testset "Bridged variable in `SingleVariable` constraint with $S" for T in [Float64, Int], S in [MOIT.UnknownScalarSet{T}]
    set = S(one(T))
    @testset "No constraint bridge" begin
        model = StandardLPModel{T}()
        bridged = MOIB.Variable.Vectorize{T}(model)
        x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(one(T)))
        fx = MOI.SingleVariable(x)
        err = _functionize_error(bridged, MOIB.Constraint.ScalarFunctionizeBridge, "SingleVariable", "constraint")
        @test_throws err MOI.add_constraint(bridged, fx, set)
    end
    @testset "LazyBridgeOptimizer" begin
        model = StandardLPModel{T}()
        bridged = MOIB.LazyBridgeOptimizer(model)
        MOIB.add_bridge(bridged, MOIB.Variable.VectorizeBridge{T})
        x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(one(T)))
        fx = MOI.SingleVariable(x)
        @testset "without `Constraint.ScalarFunctionizeBridge`" begin
            err = _lazy_functionize_error(MOIB.Constraint.ScalarFunctionizeBridge, "SingleVariable", "constraint")
            @test_throws err MOI.add_constraint(bridged, fx, set)
        end
        @testset "with `Constraint.ScalarFunctionizeBridge`" begin
            MOIB.add_bridge(bridged, MOIB.Constraint.ScalarFunctionizeBridge{T})
            cx = MOI.add_constraint(bridged, fx, set)
            @test MOI.get(bridged, MOI.ConstraintFunction(), cx) == fx
            @test MOI.get(bridged, MOI.ConstraintSet(), cx) == set
            a = MOI.get(model, MOI.ListOfVariableIndices())[1]
            fa = MOI.SingleVariable(a)
            ca = MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, S}())[1]
            @test MOI.get(model, MOI.ConstraintFunction(), ca) ≈ convert(MOI.ScalarAffineFunction{T}, fa)
            @test MOI.get(model, MOI.ConstraintSet(), ca) == MOIU.shift_constant(set, -one(T))
        end
    end
end

@testset "Bridged variable in `VectorOfVariables` constraint with $T" for T in [Float64, Int]
    set = MOI.Zeros(1)
    @testset "No constraint bridge" begin
        model = StandardLPModel{T}()
        bridged = MOIB.Variable.Vectorize{T}(model)
        x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(one(T)))
        fx = MOI.VectorOfVariables([x])
        err = _functionize_error(bridged, MOIB.Constraint.VectorFunctionizeBridge, "VectorOfVariables", "constraint")
        @test_throws err MOI.add_constraint(bridged, fx, set)
    end
    @testset "LazyBridgeOptimizer" begin
        model = StandardLPModel{T}()
        bridged = MOIB.LazyBridgeOptimizer(model)
        MOIB.add_bridge(bridged, MOIB.Variable.VectorizeBridge{T})
        MOIB.add_bridge(bridged, MOIB.Constraint.ScalarizeBridge{T})
        x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(one(T)))
        fx = MOI.VectorOfVariables([x])
        @testset "without `Constraint.ScalarFunctionizeBridge`" begin
            err = _lazy_functionize_error(MOIB.Constraint.VectorFunctionizeBridge, "VectorOfVariables", "constraint")
            @test_throws err MOI.add_constraint(bridged, fx, set)
        end
        @testset "with `Constraint.ScalarFunctionizeBridge`" begin
            MOIB.add_bridge(bridged, MOIB.Constraint.VectorFunctionizeBridge{T})
            cx = MOI.add_constraint(bridged, fx, set)
            @test MOI.get(bridged, MOI.ConstraintFunction(), cx) == fx
            @test MOI.get(bridged, MOI.ConstraintSet(), cx) == set
            a = MOI.get(model, MOI.ListOfVariableIndices())[1]
            fa = MOI.SingleVariable(a)
            ca = MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}())[1]
            @test MOI.get(model, MOI.ConstraintFunction(), ca) ≈ convert(MOI.ScalarAffineFunction{T}, fa)
            @test MOI.get(model, MOI.ConstraintSet(), ca) == MOI.EqualTo(-one(T))
        end
    end
end

@testset "Bridged variable in `SingleVariable` objective function with $T" for T in [Float64, Int]
    @testset "No objective bridge" begin
        model = StandardLPModel{T}()
        bridged = MOIB.Variable.Vectorize{T}(model)
        x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(one(T)))
        fx = MOI.SingleVariable(x)
        err = _functionize_error(bridged, MOIB.Objective.FunctionizeBridge, "SingleVariable", "objective")
        @test_throws err MOI.set(bridged, MOI.ObjectiveFunction{typeof(fx)}(), fx)
    end
    @testset "LazyBridgeOptimizer" begin
        model = StandardLPModel{T}()
        bridged = MOIB.LazyBridgeOptimizer(model)
        MOIB.add_bridge(bridged, MOIB.Variable.VectorizeBridge{T})
        x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(one(T)))
        fx = MOI.SingleVariable(x)
        @testset "without `Objective.FunctionizeBridge`" begin
            err = _lazy_functionize_error(MOIB.Objective.FunctionizeBridge, "SingleVariable", "objective")
            @test_throws err MOI.set(bridged, MOI.ObjectiveFunction{typeof(fx)}(), fx)
        end
        @testset "with `Objective.FunctionizeBridge`" begin
            MOIB.add_bridge(bridged, MOIB.Objective.FunctionizeBridge{T})
            MOI.set(bridged, MOI.ObjectiveFunction{typeof(fx)}(), fx)
            @test MOI.get(bridged, MOI.ObjectiveFunctionType()) == MOI.SingleVariable
            @test MOI.get(bridged, MOI.ObjectiveFunction{MOI.SingleVariable}()) ≈ fx
            a = MOI.get(model, MOI.ListOfVariableIndices())[1]
            fa = MOI.SingleVariable(a)
            @test MOI.get(model, MOI.ObjectiveFunctionType()) == MOI.ScalarAffineFunction{T}
            @test MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}()) ≈ fa + one(T)
        end
    end
end

MOIU.@model(
    LPModel,
    (), (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan), (), (),
    (), (MOI.ScalarAffineFunction,), (), ()
)

@testset "Name test" begin
    model = LPModel{Float64}()
    bridged = MOIB.full_bridge_optimizer(model, Float64)
    MOIT.nametest(bridged)
end

@testset "Unit" begin
    model = LPModel{Float64}()
    bridged = MOIB.full_bridge_optimizer(model, Float64)
    MOIT.unittest(bridged, MOIT.TestConfig(solve=false), [
        # SOC and quadratic constraints not supported
        "solve_qcp_edge_cases", "delete_soc_variables"
    ])
end

# Model similar to SDPA format, it gives a good example because it does not
# support a lot hence need a lot of bridges
MOIU.@model(SDPAModel,
            (), (MOI.EqualTo,), (MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle), (),
            (), (MOI.ScalarAffineFunction,), (MOI.VectorOfVariables,), ())
MOI.supports_constraint(::SDPAModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{T}}) where {T} = false
MOI.supports_constraint(::SDPAModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{T}}) where {T} = false
MOI.supports_constraint(::SDPAModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.EqualTo{T}}) where {T} = false
MOI.supports_constraint(::SDPAModel{T}, ::Type{MOI.SingleVariable}, ::Type{MOI.Interval{T}}) where {T} = false
MOI.supports_constraint(::SDPAModel, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Reals}) = false
MOI.supports(::SDPAModel{T}, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}) where {T} = false
MOI.supports(::SDPAModel, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = false

@testset "Name test with SDPAModel{Float64}" begin
    model = SDPAModel{Float64}()
    bridged = MOIB.full_bridge_optimizer(model, Float64)
    MOIT.nametest(bridged)
end

@testset "Show SPDA model" begin
    model = SDPAModel{Float64}()
    bridged = MOIB.full_bridge_optimizer(model, Float64)
    # no bridges
    @test sprint(show, bridged) === raw"""
    MOIB.LazyBridgeOptimizer{SDPAModel{Float64}}
    with 0 variable bridges
    with 0 constraint bridges
    with 0 objective bridges
    with inner model SDPAModel{Float64}"""

    MOI.add_constrained_variable(bridged, MOI.LessThan(1.0))
    # add variable bridges
    @test sprint(show, bridged) == raw"""
    MOIB.LazyBridgeOptimizer{SDPAModel{Float64}}
    with 2 variable bridges
    with 0 constraint bridges
    with 0 objective bridges
    with inner model SDPAModel{Float64}"""
end

@testset "SDPA format with $T" for T in [Float64, Int]
    model = SDPAModel{T}()
    bridged = MOIB.LazyBridgeOptimizer(model)
    @testset "Variable" begin
        @testset "Nonpositives" begin
            @test !MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.Nonpositives)
            @test !MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.Nonpositives)
            MOIB.add_bridge(bridged, MOIB.Variable.NonposToNonnegBridge{T})
            @test MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.Nonpositives)
            @test MOIB.bridge_type(bridged, MOI.Nonpositives) == MOIB.Variable.NonposToNonnegBridge{T}
        end
        @testset "Zeros" begin
            @test !MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.Zeros)
            @test !MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.Zeros)
            MOIB.add_bridge(bridged, MOIB.Variable.ZerosBridge{T})
            @test MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.Zeros)
            @test MOIB.bridge_type(bridged, MOI.Zeros) == MOIB.Variable.ZerosBridge{T}
        end
        @testset "Free" begin
            @test !MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.Reals)
            @test !MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.Reals)
            @test_throws MOI.UnsupportedConstraint{MOI.VectorOfVariables, MOI.Reals} MOI.add_variable(bridged)
            @test_throws MOI.UnsupportedConstraint{MOI.VectorOfVariables, MOI.Reals} MOI.add_variables(bridged, 2)
            MOIB.add_bridge(bridged, MOIB.Variable.FreeBridge{T})
            @test MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.Reals)
            @test MOIB.bridge_type(bridged, MOI.Reals) == MOIB.Variable.FreeBridge{T}
        end
        @testset "Vectorize" begin
            @test !MOI.supports_constraint(model, MOI.SingleVariable, MOI.GreaterThan{T})
            @test !MOI.supports_constraint(bridged, MOI.SingleVariable, MOI.GreaterThan{T})
            @test !MOI.supports_constraint(model, MOI.SingleVariable, MOI.LessThan{T})
            @test !MOI.supports_constraint(bridged, MOI.SingleVariable, MOI.LessThan{T})
            @test !MOI.supports_constraint(model, MOI.SingleVariable, MOI.EqualTo{T})
            @test !MOI.supports_constraint(bridged, MOI.SingleVariable, MOI.EqualTo{T})
            MOIB.add_bridge(bridged, MOIB.Variable.VectorizeBridge{T})
            @test MOI.supports_constraint(bridged, MOI.SingleVariable, MOI.GreaterThan{T})
            @test MOIB.bridge_type(bridged, MOI.GreaterThan{T}) == MOIB.Variable.VectorizeBridge{T, MOI.Nonnegatives}
            @test MOI.supports_constraint(bridged, MOI.SingleVariable, MOI.LessThan{T})
            @test MOIB.bridge_type(bridged, MOI.LessThan{T}) == MOIB.Variable.VectorizeBridge{T, MOI.Nonpositives}
            @test MOI.supports_constraint(bridged, MOI.SingleVariable, MOI.EqualTo{T})
            @test MOIB.bridge_type(bridged, MOI.EqualTo{T}) == MOIB.Variable.VectorizeBridge{T, MOI.Zeros}
        end
        @testset "RSOCtoPSD" begin
            @test !MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
            @test !MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
            MOIB.add_bridge(bridged, MOIB.Variable.RSOCtoPSDBridge{T})
            @test !MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
            MOIB.add_bridge(bridged, MOIB.Constraint.ScalarFunctionizeBridge{T})
            @test MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
            @test MOIB.bridge_type(bridged, MOI.RotatedSecondOrderCone) == MOIB.Variable.RSOCtoPSDBridge{T}
        end
        @testset "Combining two briges" begin
            x, cx = MOI.add_constrained_variable(bridged, MOI.LessThan(one(T)))
            test_delete_bridged_variable(bridged, x, MOI.LessThan{T}, 1, (
                (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
                (MOI.VectorOfVariables, MOI.Nonpositives, 0)),
                used_bridges = 2)
        end
    end
    @testset "Constraint" begin
        @testset "Slack" begin
            @test !MOI.supports_constraint(bridged, MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone)
            MOIB.add_bridge(bridged, MOIB.Constraint.VectorSlackBridge{T})
            @test !MOI.supports_constraint(bridged, MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone)
            MOIB.add_bridge(bridged, MOIB.Constraint.ScalarizeBridge{T})
            @test MOI.supports_constraint(bridged, MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone)
            @test MOIB.bridge_type(bridged, MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone) ==
                MOIB.Constraint.VectorSlackBridge{T, MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone}
            @test MOI.supports_constraint(bridged, MOI.VectorAffineFunction{T}, MOI.Zeros)
            @test MOIB.bridge_type(bridged, MOI.VectorAffineFunction{T}, MOI.Zeros) ==
                MOIB.Constraint.ScalarizeBridge{T, MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
        end
        @testset "Vectorize" begin
            @test !MOI.supports_constraint(bridged, MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T})
            MOIB.add_bridge(bridged, MOIB.Constraint.VectorizeBridge{T})
            @test MOI.supports_constraint(bridged, MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T})
            @test MOIB.bridge_type(bridged, MOI.ScalarAffineFunction{T},
                                   MOI.GreaterThan{T}) ==
                MOIB.Constraint.VectorizeBridge{
                    T, MOI.VectorAffineFunction{T}, MOI.Nonnegatives,
                    MOI.ScalarAffineFunction{T}
                }

        end
        @testset "Quadratic" begin
            @test !MOI.supports_constraint(bridged, MOI.ScalarQuadraticFunction{T}, MOI.GreaterThan{T})
            @test !MOI.supports_constraint(bridged, MOI.ScalarQuadraticFunction{T}, MOI.LessThan{T})
            MOIB.add_bridge(bridged, MOIB.Constraint.QuadtoSOCBridge{T})
            @test MOI.supports_constraint(bridged, MOI.ScalarQuadraticFunction{T}, MOI.GreaterThan{T})
            @test MOIB.bridge_type(bridged, MOI.ScalarQuadraticFunction{T},
                                   MOI.GreaterThan{T}) == MOIB.Constraint.QuadtoSOCBridge{T}
            @test MOI.supports_constraint(bridged, MOI.ScalarQuadraticFunction{T}, MOI.LessThan{T})
            @test MOIB.bridge_type(bridged, MOI.ScalarQuadraticFunction{T},
                                   MOI.LessThan{T}) == MOIB.Constraint.QuadtoSOCBridge{T}
        end
    end
    @testset "Objective" begin
        F = MOI.ScalarQuadraticFunction{T}
        @test !MOI.supports(bridged, MOI.ObjectiveFunction{MOI.SingleVariable}())
        @test !MOI.supports(bridged, MOI.ObjectiveFunction{F}())
        MOIB.add_bridge(bridged, MOIB.Objective.SlackBridge{T})
        @test !MOI.supports(bridged, MOI.ObjectiveFunction{MOI.SingleVariable}())
        @test !MOI.supports(bridged, MOI.ObjectiveFunction{F}())
        MOIB.add_bridge(bridged, MOIB.Objective.FunctionizeBridge{T})
        @test MOI.supports(bridged, MOI.ObjectiveFunction{MOI.SingleVariable}())
        @test MOIB.bridge_type(bridged, MOI.SingleVariable) == MOIB.Objective.FunctionizeBridge{T}
        @test MOI.supports(bridged, MOI.ObjectiveFunction{F}())
        @test MOIB.bridge_type(bridged, F) == MOIB.Objective.SlackBridge{T, F, F}
    end
end

@testset "SDPA debug with $T" for T in [Float64, Int]
    model = SDPAModel{T}()
    bridged = MOIB.LazyBridgeOptimizer(model)
    function debug_string(f, args...)
        s = IOBuffer()
        f(bridged, args...; io = s)
        return String(resize!(s.data, s.size))
    end
    @testset "LessThan variables" begin
        F = MOI.SingleVariable
        S = MOI.LessThan{T}
        @test debug_string(MOIB.debug_supports_constraint, F, S) == """
Constrained variables in `MOI.LessThan{$T}` are not supported and cannot be bridged into supported constrained variables and constraints. See details below:
[1] constraint variables in `MOI.LessThan{$T}` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (1) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported
(1) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Variable.VectorizeBridge{T})
        @test debug_string(MOIB.debug_supports_constraint, F, S) == """
Constrained variables in `MOI.LessThan{$T}` are not supported and cannot be bridged into supported constrained variables and constraints. See details below:
[1] constraint variables in `MOI.LessThan{$T}` are not supported because:
  Cannot use `MOIB.Variable.VectorizeBridge{$T,MOI.Nonpositives}` because:
  [2] constrained variables in `MOI.Nonpositives` are not supported
  Cannot add free variables and then constrain them because:
  (2) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported
[2] constraint variables in `MOI.Nonpositives` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (1) `MOI.VectorOfVariables`-in-`MOI.Nonpositives` constraints are not supported
(1) `MOI.VectorOfVariables`-in-`MOI.Nonpositives` constraints are not supported because no added bridge supports bridging it.
(2) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Variable.NonposToNonnegBridge{T})
        @test debug_string(MOIB.debug_supports_constraint, F, S) == "Constrained variables in `MOI.LessThan{$T}` are supported.\n"
    end
    bridged = MOIB.LazyBridgeOptimizer(model)
    @testset "Interval constraint" begin
        F = MOI.ScalarAffineFunction{T}
        S = MOI.Interval{T}
        @test debug_string(MOIB.debug_supports_constraint, F, S) == """
`MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported and cannot be bridged into supported constrained variables and constraints. See details below:
(1) `MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Constraint.SplitIntervalBridge{T})
        @test debug_string(MOIB.debug_supports_constraint, F, S) == """
`MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported and cannot be bridged into supported constrained variables and constraints. See details below:
(1) `MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.SplitIntervalBridge{$T,MOI.ScalarAffineFunction{$T}}` because:
  (2) `MOI.ScalarAffineFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (3) `MOI.ScalarAffineFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
(2) `MOI.ScalarAffineFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported because no added bridge supports bridging it.
(3) `MOI.ScalarAffineFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Constraint.ScalarSlackBridge{T})
        @test debug_string(MOIB.debug_supports_constraint, F, S) == """
`MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported and cannot be bridged into supported constrained variables and constraints. See details below:
[1] constraint variables in `MOI.GreaterThan{$T}` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (3) `MOI.SingleVariable`-in-`MOI.GreaterThan{$T}` constraints are not supported
[2] constraint variables in `MOI.LessThan{$T}` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (5) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported
[3] constraint variables in `MOI.Interval{$T}` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (6) `MOI.SingleVariable`-in-`MOI.Interval{$T}` constraints are not supported
(1) `MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.SplitIntervalBridge{$T,MOI.ScalarAffineFunction{$T}}` because:
  (2) `MOI.ScalarAffineFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (4) `MOI.ScalarAffineFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
  Cannot use `MOIB.Constraint.ScalarSlackBridge{$T,MOI.ScalarAffineFunction{$T},MOI.Interval{$T}}` because:
  [3] constrained variables in `MOI.Interval{$T}` are not supported
(2) `MOI.ScalarAffineFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.ScalarSlackBridge{$T,MOI.ScalarAffineFunction{$T},MOI.GreaterThan{$T}}` because:
  [1] constrained variables in `MOI.GreaterThan{$T}` are not supported
(3) `MOI.SingleVariable`-in-`MOI.GreaterThan{$T}` constraints are not supported because no added bridge supports bridging it.
(4) `MOI.ScalarAffineFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.ScalarSlackBridge{$T,MOI.ScalarAffineFunction{$T},MOI.LessThan{$T}}` because:
  [2] constrained variables in `MOI.LessThan{$T}` are not supported
(5) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported because no added bridge supports bridging it.
(6) `MOI.SingleVariable`-in-`MOI.Interval{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.SplitIntervalBridge{$T,MOI.SingleVariable}` because:
  (3) `MOI.SingleVariable`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (5) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported
"""
        MOIB.add_bridge(bridged, MOIB.Variable.VectorizeBridge{T})
        @test debug_string(MOIB.debug_supports_constraint, F, S) == """
`MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported and cannot be bridged into supported constrained variables and constraints. See details below:
[1] constraint variables in `MOI.LessThan{$T}` are not supported because:
  Cannot use `MOIB.Variable.VectorizeBridge{$T,MOI.Nonpositives}` because:
  [2] constrained variables in `MOI.Nonpositives` are not supported
  Cannot add free variables and then constrain them because:
  (4) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported
[2] constraint variables in `MOI.Nonpositives` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (3) `MOI.VectorOfVariables`-in-`MOI.Nonpositives` constraints are not supported
[3] constraint variables in `MOI.Interval{$T}` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (5) `MOI.SingleVariable`-in-`MOI.Interval{$T}` constraints are not supported
(1) `MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.SplitIntervalBridge{$T,MOI.ScalarAffineFunction{$T}}` because:
  (2) `MOI.ScalarAffineFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
  Cannot use `MOIB.Constraint.ScalarSlackBridge{$T,MOI.ScalarAffineFunction{$T},MOI.Interval{$T}}` because:
  [3] constrained variables in `MOI.Interval{$T}` are not supported
(2) `MOI.ScalarAffineFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.ScalarSlackBridge{$T,MOI.ScalarAffineFunction{$T},MOI.LessThan{$T}}` because:
  [1] constrained variables in `MOI.LessThan{$T}` are not supported
(3) `MOI.VectorOfVariables`-in-`MOI.Nonpositives` constraints are not supported because no added bridge supports bridging it.
(4) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported because no added bridge supports bridging it.
(5) `MOI.SingleVariable`-in-`MOI.Interval{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.SplitIntervalBridge{$T,MOI.SingleVariable}` because:
  (6) `MOI.SingleVariable`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (4) `MOI.SingleVariable`-in-`MOI.LessThan{$T}` constraints are not supported
(6) `MOI.SingleVariable`-in-`MOI.GreaterThan{$T}` constraints are not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Variable.NonposToNonnegBridge{T})
        @test debug_string(MOIB.debug_supports_constraint, F, S) == "`MOI.ScalarAffineFunction{$T}`-in-`MOI.Interval{$T}` constraints are supported.\n"
    end
    bridged = MOIB.LazyBridgeOptimizer(model)
    @testset "Quadratic objective" begin
        F = MOI.ScalarQuadraticFunction{T}
        attr = MOI.ObjectiveFunction{F}()
        @test debug_string(MOIB.debug_supports, attr) == """
Objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported and cannot be bridged into a supported objective function by adding only supported constrained variables and constraints. See details below:
|1| objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Objective.SlackBridge{T})
        @test debug_string(MOIB.debug_supports, attr) == """
Objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported and cannot be bridged into a supported objective function by adding only supported constrained variables and constraints. See details below:
(1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported because no added bridge supports bridging it.
(2) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because no added bridge supports bridging it.
|1| objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported because:
  Cannot use `MOIB.Objective.SlackBridge{$T,MOI.ScalarQuadraticFunction{$T},MOI.ScalarQuadraticFunction{$T}}` because:
  (1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (2) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
  |2| objective function of type `MOI.SingleVariable` is not supported
|2| objective function of type `MOI.SingleVariable` is not supported because no added bridge supports bridging it.
"""
        MOIB.add_bridge(bridged, MOIB.Objective.FunctionizeBridge{T})
        MOIB.add_bridge(bridged, MOIB.Constraint.QuadtoSOCBridge{T})
        @test debug_string(MOIB.debug_supports, attr) == """
Objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported and cannot be bridged into a supported objective function by adding only supported constrained variables and constraints. See details below:
(1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.QuadtoSOCBridge{$T}` because:
  (2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
(2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported because no added bridge supports bridging it.
(3) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.QuadtoSOCBridge{$T}` because:
  (2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
|1| objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported because:
  Cannot use `MOIB.Objective.SlackBridge{$T,MOI.ScalarQuadraticFunction{$T},MOI.ScalarQuadraticFunction{$T}}` because:
  (1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (3) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
"""
        MOIB.add_bridge(bridged, MOIB.Constraint.VectorSlackBridge{T})
        @test debug_string(MOIB.debug_supports, attr) == """
Objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported and cannot be bridged into a supported objective function by adding only supported constrained variables and constraints. See details below:
[1] constraint variables in `MOI.RotatedSecondOrderCone` are not supported because no added bridge supports bridging it.
  Cannot add free variables and then constrain them because:
  (3) `MOI.VectorOfVariables`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
(1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.QuadtoSOCBridge{$T}` because:
  (2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
(2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported because:
  Cannot use `MOIB.Constraint.VectorSlackBridge{$T,MOI.VectorAffineFunction{$T},MOI.RotatedSecondOrderCone}` because:
  [1] constrained variables in `MOI.RotatedSecondOrderCone` are not supported
  (4) `MOI.VectorAffineFunction{$T}`-in-`MOI.Zeros` constraints are not supported
(3) `MOI.VectorOfVariables`-in-`MOI.RotatedSecondOrderCone` constraints are not supported because no added bridge supports bridging it.
(4) `MOI.VectorAffineFunction{$T}`-in-`MOI.Zeros` constraints are not supported because no added bridge supports bridging it.
(5) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.QuadtoSOCBridge{$T}` because:
  (2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
|1| objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported because:
  Cannot use `MOIB.Objective.SlackBridge{$T,MOI.ScalarQuadraticFunction{$T},MOI.ScalarQuadraticFunction{$T}}` because:
  (1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (5) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
"""
        MOIB.add_bridge(bridged, MOIB.Variable.RSOCtoPSDBridge{T})
        MOIB.add_bridge(bridged, MOIB.Constraint.ScalarFunctionizeBridge{T})
        @test debug_string(MOIB.debug_supports, attr) == """
Objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported and cannot be bridged into a supported objective function by adding only supported constrained variables and constraints. See details below:
(1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.QuadtoSOCBridge{$T}` because:
  (2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
(2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported because:
  Cannot use `MOIB.Constraint.VectorSlackBridge{$T,MOI.VectorAffineFunction{$T},MOI.RotatedSecondOrderCone}` because:
  (3) `MOI.VectorAffineFunction{$T}`-in-`MOI.Zeros` constraints are not supported
(3) `MOI.VectorAffineFunction{$T}`-in-`MOI.Zeros` constraints are not supported because no added bridge supports bridging it.
(4) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported because:
  Cannot use `MOIB.Constraint.QuadtoSOCBridge{$T}` because:
  (2) `MOI.VectorAffineFunction{$T}`-in-`MOI.RotatedSecondOrderCone` constraints are not supported
|1| objective function of type `MOI.ScalarQuadraticFunction{$T}` is not supported because:
  Cannot use `MOIB.Objective.SlackBridge{$T,MOI.ScalarQuadraticFunction{$T},MOI.ScalarQuadraticFunction{$T}}` because:
  (1) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.GreaterThan{$T}` constraints are not supported
  (4) `MOI.ScalarQuadraticFunction{$T}`-in-`MOI.LessThan{$T}` constraints are not supported
"""
        MOIB.add_bridge(bridged, MOIB.Constraint.ScalarizeBridge{T})
        @test debug_string(MOIB.debug_supports, attr) == "Objective function of type `MOI.ScalarQuadraticFunction{$T}` is supported.\n"
    end
end

@testset "Continuous Linear with SDPAModel{$T}" for T in [Float64, Rational{Int}]
    model = SDPAModel{T}()
    bridged = MOIB.full_bridge_optimizer(model, T)
    # For `ScalarAffineFunction`-in-`GreaterThan`,
    # `Constraint.ScalarSlackBridge` -> `Variable.VectorizeBridge`
    # is equivalent to
    # `Constraint.VectorizeBridge` -> `Constraint.VectorSlackBridge`
    # however, `Variable.VectorizeBridge` do not support modification of the
    # set hence it makes some tests of `contlineartest` fail so we disable it.
    MOIB.remove_bridge(bridged, MOIB.Constraint.ScalarSlackBridge{T})
    exclude = ["partial_start"] # `VariablePrimalStart` not supported.
    MOIT.contlineartest(bridged, MOIT.TestConfig{T}(solve=false), exclude)
end

@testset "Continuous Conic with SDPAModel{Float64}" begin
    model = SDPAModel{Float64}()
    bridged = MOIB.full_bridge_optimizer(model, Float64)
    exclude = ["exp", "dualexp", "pow", "dualpow", "logdet", "rootdets"]
    MOIT.contconictest(bridged, MOIT.TestConfig(solve=false), exclude)
end

# Model not supporting RotatedSecondOrderCone
MOIU.@model(NoRSOCModel,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.NormInfinityCone, MOI.NormOneCone,
             MOI.ExponentialCone, MOI.PositiveSemidefiniteConeTriangle),
            (MOI.PowerCone,),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

# We only use floating point types as there is √2
@testset "Constrained variables in RSOC with $T" for T in [Float64, BigFloat]
    bridged = MOIB.full_bridge_optimizer(NoRSOCModel{T}(), T)
    @test MOI.supports_constraint(bridged, MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
    # It should be selected over `MOIB.Variable.RSOCtoPSDBridge` even if they
    # are tied in terms of number of bridges because it is added first in
    # `MOIB.full_bridge_optimizer`.
    @test MOIB.bridge_type(bridged, MOI.RotatedSecondOrderCone) == MOIB.Variable.RSOCtoSOCBridge{T}
    x, cx = MOI.add_constrained_variables(bridged, MOI.RotatedSecondOrderCone(3))
    for i in 1:3
        @test MOIB.bridge(bridged, x[i]) isa MOIB.Variable.RSOCtoSOCBridge{T}
    end
    @test MOIB.bridge(bridged, cx) isa MOIB.Variable.RSOCtoSOCBridge{T}
end

# Model not supporting VectorOfVariables and SingleVariable
MOIU.@model(NoVariableModel,
            (MOI.ZeroOne, MOI.Integer),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone),
            (),
            (),
            (MOI.ScalarAffineFunction,),
            (),
            (MOI.VectorAffineFunction,))
function MOI.supports_constraint(::NoVariableModel{T}, ::Type{MOI.SingleVariable},
                                 ::Type{<:MOIU.SUPPORTED_VARIABLE_SCALAR_SETS{T}}) where T
    return false
end

# Only supports GreaterThan and Nonnegatives
MOIU.@model(GreaterNonnegModel,
            (),
            (MOI.GreaterThan,),
            (MOI.Nonnegatives,),
            (),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))
function MOI.supports_constraint(
    ::GreaterNonnegModel{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.LessThan{T}, MOI.Interval{T}}}) where T
    return false
end

MOIU.@model(ModelNoVAFinSOC,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.NormInfinityCone, MOI.NormOneCone,
             MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone),
            (MOI.PowerCone, MOI.DualPowerCone),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

MOI.supports_constraint(::ModelNoVAFinSOC{Float64},
                        ::Type{MOI.VectorAffineFunction{Float64}},
                        ::Type{MOI.SecondOrderCone}) = false

# Model supporting nothing
MOIU.@model NothingModel () () () () () () () ()
function MOI.supports_constraint(
    ::NothingModel{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T},
                   MOI.Interval{T}, MOI.Integer, MOI.ZeroOne}}) where T
    return false
end

struct BridgeAddingNoConstraint{T} <: MOI.Bridges.Constraint.AbstractBridge end
MOIB.added_constrained_variable_types(::Type{<:BridgeAddingNoConstraint}) = Tuple{DataType}[]
MOIB.added_constraint_types(::Type{<:BridgeAddingNoConstraint}) = Tuple{DataType, DataType}[]
function MOI.supports_constraint(::Type{<:BridgeAddingNoConstraint},
                                 ::Type{MOI.SingleVariable},
                                 ::Type{MOI.Integer})
    return true
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:BridgeAddingNoConstraint{T}},
                                              ::Type{MOI.SingleVariable},
                                              ::Type{MOI.Integer}) where {T}
    return BridgeAddingNoConstraint{T}
end

const LessThanIndicatorSetOne{T} = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}
MOIU.@model(ModelNoZeroIndicator,
            (MOI.ZeroOne, MOI.Integer),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval,
             MOI.Semicontinuous, MOI.Semiinteger),
            (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
             MOI.NormInfinityCone, MOI.NormOneCone,
             MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
             MOI.GeometricMeanCone, MOI.ExponentialCone, MOI.DualExponentialCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.PositiveSemidefiniteConeSquare,
             MOI.RootDetConeTriangle, MOI.RootDetConeSquare, MOI.LogDetConeTriangle,
             MOI.LogDetConeSquare),
            (MOI.PowerCone, MOI.DualPowerCone, MOI.SOS1, MOI.SOS2, LessThanIndicatorSetOne),
            (), (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

@testset "Bridge adding no constraint" begin
    mock = MOIU.MockOptimizer(NothingModel{Int}())
    bridged = MOIB.LazyBridgeOptimizer(mock)
    MOI.Bridges.add_bridge(bridged, BridgeAddingNoConstraint{Float64})
    @test MOI.Bridges.supports_bridging_constraint(bridged,
                                                   MOI.SingleVariable,
                                                   MOI.Integer)
end

@testset "Unsupported constraint with cycles" begin
    # Test that `supports_constraint` works correctly when it is not
    # supported but the bridges forms a cycle
    mock = MOIU.MockOptimizer(NothingModel{Float64}())
    bridged = MOIB.full_bridge_optimizer(mock, Float64)
    @test !MOI.supports_constraint(
        bridged, MOI.SingleVariable, MOI.GreaterThan{Float64})
    @test !MOI.supports_constraint(
        bridged, MOI.VectorAffineFunction{Float64}, MOI.Nonpositives)
end

mock = MOIU.MockOptimizer(NoRSOCModel{Float64}())
bridged_mock = MOIB.LazyBridgeOptimizer(mock)

@testset "UnsupportedConstraint when it cannot be bridged" begin
    x = MOI.add_variables(bridged_mock, 4)
    err = MOI.UnsupportedConstraint{MOI.VectorOfVariables,
                                    MOI.RotatedSecondOrderCone}()
    @test_throws err begin
        MOI.add_constraint(bridged_mock, MOI.VectorOfVariables(x),
                           MOI.RotatedSecondOrderCone(4))
    end
end

MOIB.add_bridge(bridged_mock, MOIB.Constraint.SplitIntervalBridge{Float64})
MOIB.add_bridge(bridged_mock, MOIB.Constraint.RSOCtoPSDBridge{Float64})
MOIB.add_bridge(bridged_mock, MOIB.Constraint.SOCtoPSDBridge{Float64})
MOIB.add_bridge(bridged_mock, MOIB.Constraint.RSOCBridge{Float64})

@testset "Name test" begin
    MOIT.nametest(bridged_mock)
end

@testset "Copy test" begin
    MOIT.failcopytestc(bridged_mock)
    MOIT.failcopytestia(bridged_mock)
    MOIT.failcopytestva(bridged_mock)
    MOIT.failcopytestca(bridged_mock)
    MOIT.copytest(bridged_mock, NoRSOCModel{Float64}())
end

# Test that RSOCtoPSD is used instead of RSOC+SOCtoPSD as it is a shortest path.
@testset "Bridge selection" begin
    MOI.empty!(bridged_mock)
    @test !(MOI.supports_constraint(bridged_mock,
                                    MOI.VectorAffineFunction{Float64},
                                    MOI.LogDetConeTriangle))
    x = MOI.add_variables(bridged_mock, 3)
    err = MOI.UnsupportedConstraint{MOI.VectorAffineFunction{Float64},
                                    MOI.LogDetConeTriangle}()
    @test_throws err begin
        MOIB.bridge_type(bridged_mock, MOI.VectorAffineFunction{Float64},
                         MOI.LogDetConeTriangle)
    end
    c = MOI.add_constraint(bridged_mock, MOI.VectorOfVariables(x),
                           MOI.RotatedSecondOrderCone(3))
    @test MOIB.bridge_type(
        bridged_mock, MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone) == MOIB.Constraint.RSOCtoPSDBridge{
            Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}
    @test MOIB.bridge(bridged_mock, c) isa MOIB.Constraint.RSOCtoPSDBridge
    @test bridged_mock.constraint_dist[(MOI.VectorOfVariables,
                                        MOI.RotatedSecondOrderCone)] == 1
end

@testset "Supports" begin
    full_bridged_mock = MOIB.full_bridge_optimizer(mock, Float64)
    @testset "Mismatch vector/scalar" begin
        for S in [MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros]
            @test !MOI.supports_constraint(full_bridged_mock, MOI.SingleVariable, S)
        end
        for S in [MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.EqualTo{Float64}]
            @test !MOI.supports_constraint(full_bridged_mock, MOI.VectorOfVariables, S)
        end
    end
    greater_nonneg_mock = MOIU.MockOptimizer(GreaterNonnegModel{Float64}())
    full_bridged_greater_nonneg = MOIB.full_bridge_optimizer(
        greater_nonneg_mock, Float64)
    for F in [MOI.SingleVariable, MOI.ScalarAffineFunction{Float64},
              MOI.ScalarQuadraticFunction{Float64}]
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.Interval{Float64})
        @test !MOI.supports_constraint(
            greater_nonneg_mock, F, MOI.LessThan{Float64})
        @test MOI.supports_constraint(
            full_bridged_greater_nonneg, F, MOI.LessThan{Float64})
    end
    for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
              MOI.VectorQuadraticFunction{Float64}]
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.PositiveSemidefiniteConeSquare)
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.GeometricMeanCone)
        @test !MOI.supports_constraint(
            greater_nonneg_mock, F, MOI.Nonpositives)
        @test MOI.supports_constraint(
            full_bridged_greater_nonneg, F, MOI.Nonnegatives)
    end
    for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}]
        # The bridges in this for loop do not support yet
        # VectorQuadraticFunction. See TODO's for the reason.
        # TODO: Missing vcat for quadratic for supporting quadratic.
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.RotatedSecondOrderCone)
        # TODO: Det bridges need to use MOIU.operate to support quadratic.
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.LogDetConeTriangle)
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.RootDetConeTriangle)
    end
    mock2 = MOIU.MockOptimizer(ModelNoVAFinSOC{Float64}())
    @test !MOI.supports_constraint(mock2, MOI.VectorAffineFunction{Float64},
                                   MOI.SecondOrderCone)
    full_bridged_mock2 = MOIB.full_bridge_optimizer(mock2, Float64)
    @test MOI.supports_constraint(full_bridged_mock2, MOI.VectorAffineFunction{Float64},
                                  MOI.SecondOrderCone)
    mock_indicator = MOIU.MockOptimizer(ModelNoZeroIndicator{Float64}())
    full_bridged_mock_indicator = MOIB.full_bridge_optimizer(mock_indicator, Float64)
    @test !MOI.supports_constraint(mock_indicator, MOI.VectorAffineFunction{Float64},
                                MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO, MOI.LessThan{Float64}})
    @test MOI.supports_constraint(full_bridged_mock_indicator, MOI.VectorAffineFunction{Float64},
                                MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO, MOI.LessThan{Float64}})
    @testset "Unslack" begin
        for T in [Float64, Int]
            no_variable_mock = MOIU.MockOptimizer(NoVariableModel{T}())
            full_bridged_no_variable = MOIB.full_bridge_optimizer(
                no_variable_mock, T)
            for S in [MOI.LessThan{T}, MOI.GreaterThan{T}, MOI.EqualTo{T},
                      MOI.ZeroOne, MOI.Integer]
                @test MOI.supports_constraint(
                    full_bridged_no_variable, MOI.SingleVariable, S)
            end
            for S in [MOI.Nonpositives, MOI.Nonnegatives,
                      MOI.Zeros, MOI.SecondOrderCone]
                @test MOI.supports_constraint(
                    full_bridged_no_variable, MOI.VectorOfVariables, S)
            end
        end
    end
end


text:
function plot_concentrations(C_old, savepath, title)
# Visualizing
fig, ax = subplots()
x = [C_old[1].domain.facecenters.x[1]; C_old[1].domain.cellcenters.x; C_old[1].domain.facecenters.x[end]]
  phi = [0.5*(C_old[1].value[1]+C_old[1].value[2]); C_old[1].value[2:end-1]; 0.5*(C_old[1].value[end-1]+C_old[1].value[end])]
  plot(x, phi, label=L"$C_{A},\ \nu_A=$""$nu_A"L",\ a=""$a")

  x = [C_old[2].domain.facecenters.x[1]; C_old[2].domain.cellcenters.x; C_old[2].domain.facecenters.x[end]]
    phi = [0.5*(C_old[2].value[1]+C_old[2].value[2]); C_old[2].value[2:end-1]; 0.5*(C_old[2].value[end-1]+C_old[2].value[end])]
    plot(x, phi, label=L"$C_{B},\ \nu_B=$""$nu_B"L",\ b=""$b")

    x = [C_old[3].domain.facecenters.x[1]; C_old[3].domain.cellcenters.x; C_old[3].domain.facecenters.x[end]]
      phi = [0.5*(C_old[3].value[1]+C_old[3].value[2]); C_old[3].value[2:end-1]; 0.5*(C_old[3].value[end-1]+C_old[3].value[end])]
      plot(x, phi, label=L"$C_{C},\ \nu_C=$""$nu_C"L",\ c=""$c")
ax[:legend]()
xlabel(L"$L_x$")
ylabel(L"$C_{\alpha}$")
ax[:set_title](title)
ax[:legend]()

savefig(savepath)
end

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
export getBackgroundDictionaryComplete

function getBackgroundDictionaryComplete(fSF::MPIFile, f::MPIFile, frequencies,
                                       bgFrames=nothing, numBGAverages=1)
  idxBGFrames = measBGFrameIdx(fSF)

  D = measData(fSF, idxBGFrames)
  D_ = reshape(D, size(D,1), size(D,2)*size(D,3), size(D,4))
  bgdata = reshape(D_[:,frequencies,:],length(idxBGFrames),:)

  if bgFrames != nothing
    uEmpty = getMeasurementsFD(f, frequencies=frequencies, frames=bgFrames, numAverages=numBGAverages,
                               spectralLeakageCorrection=false, bgCorrection=false)
    bgdata2 = transpose(reshape(uEmpty, :, div(length(bgFrames),numBGAverages)))

    bgdata = cat(bgdata2,bgdata,dims=1)
  end

  return svd(transpose(bgdata))
end

function getBackgroundDictionary(fSF::MPIFile, f::MPIFile, frequencies,
                                 bgDictSize::Int=2, bgFrames=nothing, numBGAverages=1)
  U,S,V = getBackgroundDictionaryComplete(fSF, f, frequencies, bgFrames, numBGAverages)

  for l=1:bgDictSize
    U[:,l] *= (S[l] / S[1])^(1/2)
  end
  #return transpose(bgdata)[:,1:bgDictSize]
  return U[:,1:bgDictSize]
end

getBackgroundDictionary(::Any, ::Any, ::Any, bgDictSize::Nothing, ::Any) = nothing

"""
Joint reconstruction of MPI signal and background. Implemented in a low-level
fashion
"""
function reconstruction(S, u::Array, bgDict::AbstractMatrix;
                        sparseTrafo = nothing, beta = 0.1, β=beta,
                        lambd=0.0, lambda=lambd, λ=lambda, progress=nothing,
                        solver = "kaczmarz",
                        weights=nothing, enforceReal=false, enforcePositive=false,
                        relativeLambda=true, backgroundCoefficients = nothing, kargs...)

  N = size(S,2)
  M = div(length(S), N)

  L = size(u)[end]
  u = reshape(u, M, L)
  c = zeros(N,L)

  #β /= size(bgDict,2)#*1000

  if sum(abs.(λ)) > 0 && solver != "fusedlasso" && relativeLambda
    trace = calculateTraceOfNormalMatrix(S,weights)

    @debug "REL λ =  $(trace / N) "

    λ *= trace / N
    setlambda(S,λ)
  end

  G = transpose( cat(Float32(1/(sqrt(λ)))*transpose(S), Float32(1/(sqrt(β)))*transpose(bgDict), dims=1) )

  constraintMask = zeros(Bool, size(G,2))
  constraintMask[1:N] .= 1

  solv = createLinearSolver(solver, G; weights=weights, λ=1.0, constraintMask=constraintMask,
                            sparseTrafo=sparseTrafo, enforceReal=enforceReal,
			                enforcePositive=enforcePositive, kargs...)

  progress==nothing ? p = Progress(L, 1, "Reconstructing data...") : p = progress
  for l=1:L

    y = solve(solv, u[:,l])
    d = y[1:N,:] ./ sqrt(λ)

    if backgroundCoefficients != nothing
      append!(backgroundCoefficients, vec(y[(N+1):end,:] ./ sqrt(β)))
    end

    if sparseTrafo != nothing
      d[:] = sparseTrafo*d # backtrafo from dual space
    end

    c[:,l] = real( d )
    next!(p)
    sleep(0.001)
  end

  return c
end


text:
# add one unit to the voltage limit
# SU-
export psp_decrementvlimit

const DECREMENT_VLIMIT = [0x53, 0x55, 0x2D, 0x0D]

"""
    psp_decrementvlimit(io)

Add one unit to the voltage limit.
"""
function psp_decrementvlimit(io_psp::IO)
  write(io_psp, DECREMENT_VLIMIT)
  return nothing
end


text:
# BioSequence
# ===========
#
# A general purpose biological sequence representation.
#
# It is cheap to create a subsequence from a sequence because sequences can
# share underlying data: creating a subsequence is just allocating a new
# sequence object defining a part of the original sequence without copy.
# Destructive operations create a copy of the underlying data if and only if it
# is shared between (sub)sequences. This is often called as copy-on-write
# strategy in computer science and should be transparent to the user.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# About Internals
# ---------------
#
# The `data` field of a `BioSequence{A}` object contains binary representation
# of a biological character sequence. Each character is encoded with an encoder
# corresponding to the alphabet `A` and compactly packed into `data`. To extract
# a character from a sequence, you should decode this binary sequence with a
# decoder that is a pair of the encoder. The length of encoded binary bits is
# fixed, and hence a character at arbitrary position can be extracted in a
# constant time. To know the exact location of a character at a position, you
# can use the `bitsof(seq, i)` function, which returns a pair of element's index
# containing binary bits and bits' offset. As a whole, character extraction
# `seq[i]` can be written as:
#
#     j = bitindex(seq, i)
#     decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
#
#  index :        index(j) - 1       index(j)       index(j) + 1
#   data : ....|xxxxx...........|xxXxxxxxxxxxxxxx|............xxxx|....
# offset :                          |<-offset(j)-|
#  width :      |<---- 64 ---->| |<---- 64 ---->| |<---- 64 ---->|
#
#  * '.' : unused (4 bits/char)
#  * 'x' : used
#  * 'X' : used and pointed by index `i`

"""
Biological sequence data structure indexed by an alphabet type `A`.
"""
type BioSequence{A<:Alphabet} <: Sequence
    data::Vector{UInt64}  # encoded character sequence data
    part::UnitRange{Int}  # interval within `data` defining the (sub)sequence
    shared::Bool          # true if and only if `data` is shared between sequences

    function BioSequence(data::Vector{UInt64}, part::UnitRange{Int}, shared::Bool)
        return new(data, part, shared)
    end
end

typealias DNASequence       BioSequence{DNAAlphabet{4}}
typealias RNASequence       BioSequence{RNAAlphabet{4}}
typealias AminoAcidSequence BioSequence{AminoAcidAlphabet}
typealias CharSequence      BioSequence{CharAlphabet}

# String Decorators
# -----------------

remove_newlines(s) = replace(s, r"\r|\n", "")

macro dna_str(seq, flags...)
    return DNASequence(remove_newlines(seq))
end

macro rna_str(seq, flags...)
    return RNASequence(remove_newlines(seq))
end

macro aa_str(seq, flags...)
    return AminoAcidSequence(remove_newlines(seq))
end

macro char_str(seq, flags...)
    return CharSequence(remove_newlines(seq))
end

# Constructors
# ------------

function seq_data_len{A}(::Type{A}, len::Integer)
    return cld(len, div(64, bitsof(A)))
end

@compat function (::Type{BioSequence{A}}){A<:Alphabet}(len::Integer)
    return BioSequence{A}(Vector{UInt64}(seq_data_len(A, len)), 1:len, false)
end

BioSequence(::Type{DNANucleotide}) = DNASequence()
BioSequence(::Type{RNANucleotide}) = RNASequence()
BioSequence(::Type{AminoAcid}) = AminoAcidSequence()
BioSequence(::Type{Char}) = CharSequence()

function BioSequence()
    return BioSequence{VoidAlphabet}(Vector{UInt64}(0), 0:-1, false)
end

@compat function (::Type{BioSequence{A}}){A<:Alphabet}(
        src::Union{AbstractString,AbstractVector},
        startpos::Integer=1,
        stoppos::Integer=length(src))
    len = stoppos - startpos + 1
    seq = BioSequence{A}(len)
    return encode_copy!(seq, 1, src, startpos, len)
end

# create a subsequence
function BioSequence{A,T<:Integer}(other::BioSequence{A}, part::UnitRange{T})
    checkbounds(other, part)
    start = other.part.start + part.start - 1
    stop = start + length(part) - 1
    subseq = BioSequence{A}(other.data, start:stop, true)
    other.shared = true
    return subseq
end

@compat function (::Type{BioSequence{A}}){A}(other::BioSequence{A}, part::UnitRange)
    return BioSequence(other, part)
end

# concatenate chunks
@compat function (::Type{BioSequence{A}}){A}(chunks::BioSequence{A}...)
    len = 0
    for chunk in chunks
        len += length(chunk)
    end
    seq = BioSequence{A}(len)
    offset = 1
    for chunk in chunks
        copy!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

function Base.repeat{A}(chunk::BioSequence{A}, n::Integer)
    seq = BioSequence{A}(length(chunk) * n)
    offset = 1
    for _ in 1:n
        copy!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

# operators for concat and repeat
@compat begin
    Base.:*{A}(chunk::BioSequence{A}, chunks::BioSequence{A}...) =
        BioSequence{A}(chunk, chunks...)
    Base.:^(chunk::BioSequence, n::Integer) = repeat(chunk, n)
end

# conversion between different alphabet size
for A in [DNAAlphabet, RNAAlphabet]
    # 4 bits => 2 bits
    @eval function Base.convert(::Type{BioSequence{$(A{2})}}, seq::BioSequence{$(A{4})})
        # TODO: make it faster with bit-parallel algorithm
        newseq = BioSequence{$(A{2})}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end

    # 2 bits => 4 bits
    @eval function Base.convert(::Type{BioSequence{$(A{4})}}, seq::BioSequence{$(A{2})})
        newseq = BioSequence{$(A{4})}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end
end

# conversion between DNA and RNA
for (A1, A2) in [(DNAAlphabet, RNAAlphabet), (RNAAlphabet, DNAAlphabet)], n in (2, 4)
    # NOTE: assumes that binary representation is identical between DNA and RNA
    @eval function Base.convert(::Type{BioSequence{$(A1{n})}},
                                seq::BioSequence{$(A2{n})})
        newseq = BioSequence{$(A1{n})}(seq.data, seq.part, true)
        seq.shared = true
        return newseq
    end
end

# from a vector
function Base.convert{A<:DNAAlphabet}(::Type{BioSequence{A}},
                                      seq::AbstractVector{DNANucleotide})
    return BioSequence{A}(seq, 1, endof(seq))
end
function Base.convert{A<:RNAAlphabet}(::Type{BioSequence{A}},
                                      seq::AbstractVector{RNANucleotide})
    return BioSequence{A}(seq, 1, endof(seq))
end
function Base.convert(::Type{AminoAcidSequence}, seq::AbstractVector{AminoAcid})
    return AminoAcidSequence(seq, 1, endof(seq))
end

# to a vector
Base.convert(::Type{Vector}, seq::BioSequence) = collect(seq)
function Base.convert{A<:DNAAlphabet}(::Type{Vector{DNANucleotide}},
                                      seq::BioSequence{A})
    return collect(seq)
end
function Base.convert{A<:RNAAlphabet}(::Type{Vector{RNANucleotide}},
                                      seq::BioSequence{A})
    return collect(seq)
end
Base.convert(::Type{Vector{AminoAcid}}, seq::AminoAcidSequence) = collect(seq)

# from/to a string
function Base.convert{S<:AbstractString}(::Type{S}, seq::BioSequence)
    return convert(S, [Char(x) for x in seq])
end
Base.convert{S<:AbstractString,A}(::Type{BioSequence{A}}, seq::S) = BioSequence{A}(seq)

# Summaries
# ---------

Base.summary{A<:DNAAlphabet}(seq::BioSequence{A}) = string(length(seq), "nt ", "DNA Sequence")
Base.summary{A<:RNAAlphabet}(seq::BioSequence{A}) = string(length(seq), "nt ", "RNA Sequence")
Base.summary(seq::AminoAcidSequence) = string(length(seq), "aa ", "Amino Acid Sequence")
Base.summary(seq::CharSequence) = string(length(seq), "char ", "Char Sequence")

# Basic Operators
# ---------------

alphabet{A}(::Type{BioSequence{A}}) = alphabet(A)

Base.length(seq::BioSequence) = length(seq.part)
Base.eltype{A}(::Type{BioSequence{A}}) = eltype(A)

function Base.count(f::Function, seq::BioSequence)
    n = 0
    for x in seq
        if f(x)
            n += 1
        end
    end
    return n
end

function Base.map(f::Function, seq::BioSequence)
    return map!(f, copy(seq))
end

function Base.filter(f::Function, seq::BioSequence)
    return filter!(f, copy(seq))
end

function Base.checkbounds(seq::BioSequence, range::UnitRange)
    if 1 ≤ range.start && range.stop ≤ endof(seq)
        return true
    end
    throw(BoundsError(seq, range))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector{Bool})
    if length(seq) == length(locs)
        return true
    end
    throw(BoundsError(seq, locs))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector)
    for i in locs
        checkbounds(seq, i)
    end
    return true
end

function checkdimension(from::Integer, to::Integer)
    if from == to
        return true
    end
    throw(DimensionMismatch(string(
        "attempt to assign ",
        from, " elements to ",
        to,   " elements")))
end

function checkdimension(seq::BioSequence, locs::AbstractVector)
    return checkdimension(length(seq), length(locs))
end

function checkdimension(seq::BioSequence, locs::AbstractVector{Bool})
    return checkdimension(length(seq), sum(locs))
end

# creates a bit mask for given number of bits `n`
mask(n::Integer) = (UInt64(1) << n) - 1
mask{A<:Alphabet}(::Type{A}) = mask(bitsof(A))

# assumes `i` is positive and `bitsof(A)` is a power of 2
@inline function bitindex{A}(seq::BioSequence{A}, i::Integer)
    return BitIndex((i + first(seq.part) - 2) << trailing_zeros(bitsof(A)))
end

@inline function inbounds_getindex{A}(seq::BioSequence{A}, i::Integer)
    j = bitindex(seq, i)
    @inbounds return decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
end

Base.getindex(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)
Base.sub(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A},
                                      other::BioSequence{A},
                                      locs::AbstractVector{T})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    for (i, x) in zip(locs, other)
        unsafe_setindex!(seq, x, i)
    end
    return seq
end

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A},
                                      other::BioSequence{A},
                                      locs::UnitRange{T})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    return copy!(seq, locs.start, other, 1)
end

function Base.setindex!{A}(seq::BioSequence{A},
                           other::BioSequence{A},
                           locs::AbstractVector{Bool})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    i = j = 0
    while (i = findnext(locs, i + 1)) > 0
        unsafe_setindex!(seq, other[j+=1], i)
    end
    return seq
end

function Base.setindex!{A}(seq::BioSequence{A}, other::BioSequence{A}, ::Colon)
    return setindex!(seq, other, 1:endof(seq))
end

function Base.setindex!(seq::BioSequence, x, i::Integer)
    checkbounds(seq, i)
    orphan!(seq)
    return unsafe_setindex!(seq, x, i)
end

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A}, x, locs::AbstractVector{T})
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    for i in locs
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function Base.setindex!{A}(seq::BioSequence{A}, x, locs::AbstractVector{Bool})
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    i = j = 0
    while (i = findnext(locs, i + 1)) > 0
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

Base.setindex!{A}(seq::BioSequence{A}, x, ::Colon) = setindex!(seq, x, 1:endof(seq))

# this is "unsafe" because of no bounds check and no orphan! call
@inline function unsafe_setindex!{A}(seq::BioSequence{A}, x, i::Integer)
    bin = enc64(seq, x)
    return encoded_setindex!(seq, bin, i)
end

@inline function encoded_setindex!{A}(seq::BioSequence{A}, bin::UInt64, i::Integer)
    j, r = bitindex(seq, i)
    @inbounds seq.data[j] = (bin << r) | (seq.data[j] & ~(mask(A) << r))
    return seq
end

function Base.resize!{A}(seq::BioSequence{A}, size::Integer)
    if size < 0
        throw(ArgumentError("size must be non-negative"))
    end
    orphan!(seq, size)
    resize!(seq.data, seq_data_len(A, size + seq.part.start - 1))
    seq.part = seq.part.start:seq.part.start+size-1
    return seq
end

Base.empty!(seq::BioSequence) = resize!(seq, 0)

function Base.push!{A}(seq::BioSequence{A}, x)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    encoded_setindex!(seq, bin, endof(seq))
    return seq
end

function Base.unshift!{A}(seq::BioSequence{A}, x)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    copy!(seq, 2, seq, 1, length(seq) - 1)
    encoded_setindex!(seq, bin, 1)
    return seq
end

function Base.pop!(seq::BioSequence)
    if isempty(seq)
        throw(ArgumentError("sequence must be non-empty"))
    end
    x = seq[end]
    deleteat!(seq, endof(seq))
    return x
end

function Base.shift!(seq::BioSequence)
    if isempty(seq)
        throw(ArgumentError("sequence must be non-empty"))
    end
    x = seq[1]
    deleteat!(seq, 1)
    return x
end

function Base.insert!{A}(seq::BioSequence{A}, i::Integer, x)
    checkbounds(seq, i)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    copy!(seq, i + 1, seq, i, endof(seq) - i)
    encoded_setindex!(seq, bin, i)
    return seq
end

function Base.deleteat!{A,T<:Integer}(seq::BioSequence{A}, range::UnitRange{T})
    checkbounds(seq, range)
    copy!(seq, range.start, seq, range.stop + 1, length(seq) - range.stop)
    resize!(seq, length(seq) - length(range))
    return seq
end

function Base.deleteat!(seq::BioSequence, i::Integer)
    checkbounds(seq, i)
    copy!(seq, i, seq, i + 1, length(seq) - i)
    resize!(seq, length(seq) - 1)
    return seq
end

function Base.append!{A}(seq::BioSequence{A}, other::BioSequence{A})
    resize!(seq, length(seq) + length(other))
    copy!(seq, endof(seq) - length(other) + 1, other, 1)
    return seq
end

function encode_copy!{A}(dst::BioSequence{A},
                         src::Union{AbstractVector,AbstractString})
    return encode_copy!(dst, 1, src, 1)
end

function encode_copy!{A}(dst::BioSequence{A},
                         doff::Integer,
                         src::Union{AbstractVector,AbstractString},
                         soff::Integer)
    return encode_copy!(dst, doff, src, soff, length(src) - soff + 1)
end

function encode_copy!{A}(dst::BioSequence{A},
                         doff::Integer,
                         src::Union{AbstractVector,AbstractString},
                         soff::Integer,
                         len::Integer)
    if soff != 1 && !isascii(src)
        throw(ArgumentError("source offset ≠ 1 is not supported for non-ASCII string"))
    end

    checkbounds(dst, doff:doff+len-1)
    if length(src) < soff + len - 1
        throw(ArgumentError("source string does not contain $len elements from $soff"))
    end

    orphan!(dst)
    next = bitindex(dst, doff)
    stop = bitindex(dst, doff + len)
    i = soff
    while next < stop
        x = UInt64(0)
        j = index(next)
        while index(next) == j && next < stop
            char, i = Base.next(src, i)
            x |= enc64(dst, convert(Char, char)) << offset(next)
            next += bitsof(A)
        end
        dst.data[j] = x
    end
    return dst
end

function Base.copy!{A}(seq::BioSequence{A}, doff::Integer,
                       src::Vector{UInt8},  soff::Integer, len::Integer)
    datalen = seq_data_len(A, len)
    return seq
end

function Base.copy!{A}(dst::BioSequence{A}, src::BioSequence{A})
    return copy!(dst, 1, src, 1)
end

function Base.copy!{A}(dst::BioSequence{A}, doff::Integer,
                       src::BioSequence{A}, soff::Integer)
    return copy!(dst, doff, src, soff, length(src) - soff + 1)
end

function Base.copy!{A}(dst::BioSequence{A}, doff::Integer,
                       src::BioSequence{A}, soff::Integer, len::Integer)
    checkbounds(dst, doff:doff+len-1)
    checkbounds(src, soff:soff+len-1)

    if dst.shared || (dst === src && doff > soff)
        orphan!(dst, length(dst), true)
    end

    id = bitindex(dst, doff)
    is = bitindex(src, soff)
    rest = len * bitsof(A)

    while rest > 0
        # move `k` bits from `src` to `dst`
        x = dst.data[index(id)]
        y = src.data[index(is)]
        if offset(id) < offset(is)
            y >>= offset(is) - offset(id)
            k = min(64 - offset(is), rest)
        else
            y <<= offset(id) - offset(is)
            k = min(64 - offset(id), rest)
        end
        m = mask(k) << offset(id)
        dst.data[index(id)] = y & m | x & ~m

        id += k
        is += k
        rest -= k
    end

    return dst
end

function Base.map!(f::Function, seq::BioSequence)
    orphan!(seq)
    for i in 1:endof(seq)
        unsafe_setindex!(seq, f(inbounds_getindex(seq, i)), i)
    end
    return seq
end

function Base.filter!{A}(f::Function, seq::BioSequence{A})
    orphan!(seq)

    len = 0
    next = bitindex(seq, 1)
    j = index(next)
    datum::UInt64 = 0
    for i in 1:endof(seq)
        x = inbounds_getindex(seq, i)
        if f(x)
            datum |= enc64(seq, x) << offset(next)
            len += 1
            next += bitsof(A)
            if index(next) != j
                seq.data[j] = datum
                datum = 0
                j = index(next)
            end
        end
    end
    if offset(next) > 0
        seq.data[j] = datum
    end
    resize!(seq, len)

    return seq
end

enc64{A}(::BioSequence{A}, x) = UInt64(encode(A, convert(eltype(A), x)))

# Replace a BioSequence's data with a copy, copying only what's needed.
# The user should never need to call this, as it has no outward effect on the
# sequence.
function orphan!{A}(seq::BioSequence{A}, size::Integer=length(seq), force::Bool=false)
    if !seq.shared && !force
        return seq
    end

    j, r = bitindex(seq, 1)
    data = Vector{UInt64}(seq_data_len(A, size))

    if !isempty(seq)
        x = seq.data[j] >> r
        l = min(endof(data), index(bitindex(seq, endof(seq))) - j + 1)
        @inbounds @simd for i in 1:l-1
            y = seq.data[j+i]
            data[i] = x | y << (64 - r)
            x = y >> r
        end
        data[l] = x
    end

    seq.data = data
    seq.part = 1:length(seq)
    seq.shared = false
    return seq
end

function Base.similar{A}(seq::BioSequence{A}, len::Integer=length(seq))
    return BioSequence{A}(len)
end

# actually, users don't need to create a copy of a sequence.
function Base.copy{A}(seq::BioSequence{A})
    # NOTE: no need to set `seq.shared = true` here
    # since `newseq` will be `orphan!`ed soon.
    newseq = BioSequence{A}(seq.data, 1:endof(seq), true)
    orphan!(newseq, length(seq), true)  # force orphan!
    @assert newseq.data !== seq.data
    return newseq
end

# Transformations
# ---------------

function Base.reverse!(seq::BioSequence)
    orphan!(seq)
    for i in 1:div(endof(seq), 2)
        x = inbounds_getindex(seq, i)
        unsafe_setindex!(seq, inbounds_getindex(seq, endof(seq) - i + 1), i)
        unsafe_setindex!(seq, x, endof(seq) - i + 1)
    end
    return seq
end

Base.reverse(seq::BioSequence) = reverse!(copy(seq))

@generated function Base.reverse{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    n = bitsof(A)
    if n == 2
        nucrev = :nucrev2
    elseif n == 4
        nucrev = :nucrev4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        data = Vector{UInt64}(seq_data_len(A, length(seq)))
        next = bitindex(seq, endof(seq))
        stop = bitindex(seq, 0)
        i = 0
        while next - stop > 0
            r = offset(next) + $n
            x = seq.data[index(next)] << (64 - r)
            next -= r
            if next - stop > 0
                x |= seq.data[index(next)] >> r
                next -= 64 - r
            end
            data[i+=1] = $nucrev(x)
        end
        return BioSequence{A}(data, 1:length(seq), false)
    end
end

@inline function nucrev2(x::UInt64)
     x = (x & 0x3333333333333333) <<  2 | (x & 0xCCCCCCCCCCCCCCCC) >>  2
     x = (x & 0x0F0F0F0F0F0F0F0F) <<  4 | (x & 0xF0F0F0F0F0F0F0F0) >>  4
     x = (x & 0x00FF00FF00FF00FF) <<  8 | (x & 0xFF00FF00FF00FF00) >>  8
     x = (x & 0x0000FFFF0000FFFF) << 16 | (x & 0xFFFF0000FFFF0000) >> 16
     x = (x & 0x00000000FFFFFFFF) << 32 | (x & 0xFFFFFFFF00000000) >> 32
     return x
end

@inline function nucrev4(x::UInt64)
     x = (x & 0x0F0F0F0F0F0F0F0F) <<  4 | (x & 0xF0F0F0F0F0F0F0F0) >>  4
     x = (x & 0x00FF00FF00FF00FF) <<  8 | (x & 0xFF00FF00FF00FF00) >>  8
     x = (x & 0x0000FFFF0000FFFF) << 16 | (x & 0xFFFF0000FFFF0000) >> 16
     x = (x & 0x00000000FFFFFFFF) << 32 | (x & 0xFFFFFFFF00000000) >> 32
     return x
end

"""
    complement!(seq)

Make a complement sequence of `seq` in place.

Ambiguous nucleotides are left as-is.
"""
@generated function complement!{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    n = bitsof(A)
    if n == 2
        nuccomp = :nuccomp2
    elseif n == 4
        nuccomp = :nuccomp4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        orphan!(seq)
        next = bitindex(seq, 1)
        stop = bitindex(seq, endof(seq) + 1)
        @inbounds while next < stop
            seq.data[index(next)] = $nuccomp(seq.data[index(next)])
            next += 64
        end
        return seq
    end
end

nuccomp2(x::UInt64) = ~x

@inline function nuccomp4(x::UInt64)
    if x & 0xCCCCCCCCCCCCCCCC != 0
        # ignore ambiguous nucleotides
        m = (x & 0x8888888888888888) | (x & 0x4444444444444444) << 1
        m |= m >> 1 | m >> 2 | m >> 3
        m = ~m
        return (~x & 0x3333333333333333 & m) | (x & ~m)
    else
        return  ~x & 0x3333333333333333
    end
end

"""
    complement(seq)

Make a complement sequence of `seq`.

Ambiguous nucleotides are left as-is.
"""
function complement{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    return complement!(copy(seq))
end

"""
    reverse_complement!(seq)

Make a reversed complement sequence of `seq` in place.

Ambiguous nucleotides are left as-is.
"""
function reverse_complement!{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    return complement!(reverse!(seq))
end

"""
    reverse_complement(seq)

Make a reversed complement sequence of `seq`.

Ambiguous nucleotides are left as-is.
"""
function reverse_complement{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    return reverse_complement!(copy(seq))
end

# Ambiguous nucleotides iterator
# ------------------------------

immutable AmbiguousNucleotideIterator{A<:Union{DNAAlphabet,RNAAlphabet}}
    seq::BioSequence{A}
end

ambiguous_positions(seq::BioSequence) = AmbiguousNucleotideIterator(seq)

Base.start(it::AmbiguousNucleotideIterator) = find_next_ambiguous(it.seq, 1)
Base.done(it::AmbiguousNucleotideIterator, nextpos) = nextpos == 0
function Base.next(it::AmbiguousNucleotideIterator, nextpos)
    return nextpos, find_next_ambiguous(it.seq, nextpos + 1)
end

if VERSION > v"0.5-"
    Base.iteratorsize(::AmbiguousNucleotideIterator) = Base.SizeUnknown()
end

function find_next_ambiguous{A<:Union{DNAAlphabet{2},RNAAlphabet{2}}}(
        seq::BioSequence{A}, i::Integer)
    # no ambiguity
    return 0
end

function find_next_ambiguous{A<:Union{DNAAlphabet{4},RNAAlphabet{4}}}(
        seq::BioSequence{A}, i::Integer)
    if i > endof(seq)
        return 0
    end

    # extract ambiguity bits using following bit masks:
    #   * 0x44... = 0b01000100...
    #   * 0x88... = 0b10001000...
    #   * 0xCC... = 0b11001100...
    next = bitindex(seq, max(i, 1))
    stop = bitindex(seq, endof(seq) + 1)
    while next < stop
        m = mask(min(stop - next, 64)) << offset(next)
        x = seq.data[index(next)] & m
        if x & 0xCCCCCCCCCCCCCCCC != 0
            x = (x & 0x8888888888888888) | (x & 0x4444444444444444) << 1
            return (index(next) - 1) << 4 + (trailing_zeros(x) + 1) >> 2
        end
        next += 64 - offset(next)
    end

    # found no ambiguity
    return 0
end

# Mismatch counting
# -----------------

"""
    mismatches(seq1::BioSequence, seq2::BioSequence[, compatible=false])

Return the number of mismatches between `seq1` and `seq2`.

If `seq1` and `seq2` are of differing lengths, only the first `min(length(seq1),
length(seq2))` nucleotides are compared.  When `compatible` is `true`, sequence
symbols are comapred using `iscompatible`; otherwise using `==`.
"""
function mismatches{A<:Alphabet}(
        seq1::BioSequence{A},
        seq2::BioSequence{A},
        compatible::Bool=false)
    if ((bitsof(A) == 2 || bitsof(A) == 4) && !compatible) ||
        A == DNAAlphabet{2} ||
        A == RNAAlphabet{2}
        return bitparallel_mismatches(seq1, seq2)
    end

    mis = 0
    if compatible
        for (x, y) in zip(seq1, seq2)
            if !iscompatible(x, y)
                mis += 1
            end
        end
    else
        for (x, y) in zip(seq1, seq2)
            if x != y
                mis += 1
            end
        end
    end
    return mis
end

@generated function bitparallel_mismatches{A}(a::BioSequence{A}, b::BioSequence{A})
    n = bitsof(A)
    if n == 2
        bitpar_mismatches = :bitpar_mismatches2
    elseif n == 4
        bitpar_mismatches = :bitpar_mismatches4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        if length(a) > length(b)
            return bitparallel_mismatches(b, a)
        end
        @assert length(a) ≤ length(b)

        nexta = bitindex(a, 1)
        nextb = bitindex(b, 1)
        stopa = bitindex(a, endof(a) + 1)
        mismatches = 0

        # align reading position of `a.data` so that `offset(nexta) == 0`
        if nexta < stopa && offset(nexta) != 0
            x = a.data[index(nexta)] >> offset(nexta)
            y = b.data[index(nextb)] >> offset(nextb)
            if offset(nextb) > offset(nexta)
                y |= b.data[index(nextb)+1] << (64 - offset(nextb))
            end
            k = 64 - offset(nexta)
            m = mask(k)
            mismatches += $bitpar_mismatches(x & m, y & m)
            nexta += k
            nextb += k
        end
        @assert offset(nexta) == 0

        if offset(nextb) == 0  # data are aligned with each other
            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                mismatches += $bitpar_mismatches(x, y)
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                m = mask(stopa - nexta)
                mismatches += $bitpar_mismatches(x & m, y & m)
            end
        elseif nexta < stopa
            y = b.data[index(nextb)]
            nextb += 64

            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                z = b.data[index(nextb)]
                y = y >> offset(nextb) | z << (64 - offset(nextb))
                mismatches += $bitpar_mismatches(x, y)
                y = z
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                x = a.data[index(nexta)]
                y = y >> offset(nextb)
                if 64 - offset(nextb) < stopa - nexta
                    y |= a.data[index(nextb)] << (64 - offset(nextb))
                end
                m = mask(stopa - nexta)
                mismatches += $bitpar_mismatches(x & m, y & m)
            end
        end

        return mismatches
    end
end

# bit-parallel mismatch count algorithm for 2 and 4-bit encoding
@inline function bitpar_mismatches2(x::UInt64, y::UInt64)
    xyxor = x $ y
    mismatches = UInt64(0)
    mismatches |=  xyxor & 0x5555555555555555
    mismatches |= (xyxor & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(mismatches)
end

@inline function bitpar_mismatches4(x::UInt64, y::UInt64)
    xyxor = x $ y
    mismatches = UInt64(0)
    mismatches |=  xyxor & 0x1111111111111111
    mismatches |= (xyxor & 0x2222222222222222) >> 1
    mismatches |= (xyxor & 0x4444444444444444) >> 2
    mismatches |= (xyxor & 0x8888888888888888) >> 3
    return count_ones(mismatches)
end

# Shuffle
# -------

function Base.shuffle(seq::BioSequence)
    return shuffle!(copy(seq))
end

function Base.shuffle!(seq::BioSequence)
    orphan!(seq)
    # Fisher-Yates shuffle
    for i in 1:endof(seq)-1
        j = rand(i:endof(seq))
        seq[i], seq[j] = seq[j], seq[i]
    end
    return seq
end


text:
using Test: @test, @testset, @inferred
using LinearAlgebra: I, norm, Symmetric, Diagonal, LowerTriangular
using WhitenedRBMs: Affine, whitening_transform, whitening_transform!

@testset "affine" begin
    x = randn(5, 7)
    t = Affine(randn(5,5), randn(5))
    s = Affine(randn(5,5), randn(5))

    @test Affine(5).A ≈ I
    @test iszero(Affine(5).u)
    @test Affine(5) * x ≈ x

    @test t * x ≈ t.A * (x .- t.u)
    @test t \ (t * x) ≈ t * (t \ x) ≈ x

    @test (s * t) * x ≈ s * (t * x)
    @test (s \ t) * x ≈ s \ (t * x)
    @test (t \ t) * x ≈ x

    @test (s + t) * x ≈ s * x + t * x
    @test (s - t) * x ≈ s * x - t * x

    @test (2 * t) * x ≈ 2 * (t * x)
    @test (t * 2) * x ≈ t * (2 * x)
    @test (-t) * x ≈ -(t * x)

    @test inv(t) * x ≈ t \ x
    @test inv(t) * (t * x) ≈ (inv(t) * t) * x ≈ x

    @test (inv(t) * t).A ≈ I
    @test norm((inv(t) * t).u) < 1e-10
end

@testset "whitening_transform" begin
    μ = randn(5)
    C = randn(5,5)
    C = C * C'
    affine = whitening_transform(μ, C)
    @test norm(affine * μ) < 1e-10
    @test affine.A * C * affine.A' ≈ I
    @test affine.A' * affine.A ≈ inv(C)

    affine = whitening_transform(μ)
    @test affine.A ≈ I
    @test affine.u ≈ μ
    @test norm(affine * μ) < 1e-10

    affine = whitening_transform(C)
    @test affine.A * C * affine.A' ≈ I
    @test iszero(affine.u)

    affine = whitening_transform(C) * whitening_transform(μ)
    @test norm(affine * μ) < 1e-10
    @test affine.A * C * affine.A' ≈ I
end

@testset "whitening_transform!" begin
    μ = randn(5)
    C = randn(5,5)
    C = C * C'
    affine = whitening_transform!(Affine(LowerTriangular(randn(5,5)), randn(5)), μ, C)
    @test norm(affine * μ) < 1e-10
    @test affine.A * C * affine.A' ≈ I
    @test affine.A' * affine.A ≈ inv(C)

    affine = whitening_transform!(Affine(Diagonal(randn(5)), randn(5)), μ)
    @test affine.A ≈ I
    @test affine.u ≈ μ
    @test norm(affine * μ) < 1e-10

    affine = whitening_transform!(Affine(LowerTriangular(randn(5,5)), randn(5)), C)
    @test affine.A * C * affine.A' ≈ I
    @test iszero(affine.u)

    C = Diagonal(rand(5))
    affine = whitening_transform!(Affine(Diagonal(randn(5)), randn(5)), μ, C)
    @test norm(affine * μ) < 1e-10
    @test affine.A * C * affine.A' ≈ I
    @test affine.A' * affine.A ≈ inv(C)
end


text:
using Test
using NeXLSpectrum

using DataDeps

register(DataDep("map[15]Artifact",
    """
    Dataset:    HyperSpectrum Test Data Set
    Author:     Nicholas W. M. Ritchie
    License:    CC-SA 3.0
    """,
    "https://drive.google.com/uc?export=download&id=1C93kn9-EIXXMDPcqJ9E4Xt4j9qfs5eeX",
    "c38ea6be06c6ef7eba427c50ecd44bc7dc4c63b5f8ef21897a8a51771edb7473",
    post_fetch_method=DataDeps.unpack
))
ENV["DATADEPS_ALWAYS_ACCEPT"]="true"

include("window.jl")
include("detefficiency.jl")
include("spectrum.jl")
include("emsa.jl")
include("brukerspx.jl")
include("filterfit.jl")
include("hyperspectrum.jl")
include("line.jl")
include("multikev.jl")
include("standardize.jl")
include("multidet.jl")

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
"""
    struct InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor}

Represents an infinite matrix product state
The state is stored in the centergauge where
    state.AL[i]*state.CR[i] = state.AC[i] = state.CR[i-1]*state.AR[i]
"""
struct InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor}
    AL::PeriodicArray{A,1}
    AR::PeriodicArray{A,1}
    CR::PeriodicArray{B,1}
    AC::PeriodicArray{A,1}
end

Base.size(arr::InfiniteMPS,i) = size(arr.AL,i)
Base.length(arr::InfiniteMPS) = size(arr,1)
Base.eltype(arr::InfiniteMPS) = eltype(arr.AL[1])
Base.copy(m::InfiniteMPS) = InfiniteMPS(copy(m.AL),copy(m.AR),copy(m.CR),copy(m.AC));
Base.repeat(m::InfiniteMPS,i::Int) = InfiniteMPS(repeat(m.AL,i),repeat(m.AR,i),repeat(m.CR,i),repeat(m.AC,i));
Base.similar(st::InfiniteMPS) = InfiniteMPS(similar(st.AL),similar(st.AR),similar(st.CR),similar(st.AC))

virtualspace(psi::InfiniteMPS, n::Integer) = _firstspace(psi.AL[n+1])

function InfiniteMPS(pspaces::AbstractArray{S,1},Dspaces::AbstractArray{S,1};eltype=Defaults.eltype,kwargs...) where S
    InfiniteMPS([TensorMap(rand,eltype,Dspaces[mod1(i-1,length(Dspaces))]*pspaces[i],Dspaces[i]) for i in 1:length(pspaces)];kwargs...)
end

#allow users to pass in simple arrays
function InfiniteMPS(A::AbstractArray{T,1}; tol::Float64 = Defaults.tolgauge, maxiter::Int64 = Defaults.maxiter,
    cguess = PeriodicArray([isomorphism(Matrix{eltype(A[i])},space(A[i],1),space(A[i],1)) for i in 1:length(A)]),leftgauged = false) where T<:GenericMPSTensor
    pA = PeriodicArray(A[:])
    #perform the left gauge fixing, remember only Al
    if leftgauged
        ALs = pA;
        deltal = 0;
    else
        ALs, cguess, deltal = uniform_leftorth(pA; tol = tol, maxiter = maxiter, cguess = cguess)
    end

    #perform the right gauge fixing from which we obtain the center matrix
    ARs, Cs, deltar  = uniform_rightorth(ALs ; tol = tol, maxiter = maxiter, cguess = cguess)
    CRs = circshift(Cs,-1);
    ACs=similar(ARs)
    for loc = 1:size(A,1)
        ACs[loc] = ALs[loc]*CRs[loc]
    end

    return InfiniteMPS(ALs,ARs,CRs,ACs)
end

"
    l_RR(state,location)
    Left dominant eigenvector of the AR-AR transfermatrix
"
l_RR(state::InfiniteMPS,loc::Int=1) = @tensor toret[-1;-2]:=state.CR[loc-1][1,-2]*conj(state.CR[loc-1][1,-1])

"
    l_RL(state,location)
    Left dominant eigenvector of the AR-AL transfermatrix
"
l_RL(state::InfiniteMPS,loc::Int=1) = state.CR[loc-1]

"
    l_LR(state,location)
    Left dominant eigenvector of the AL-AR transfermatrix
"
l_LR(state::InfiniteMPS,loc::Int=1) = state.CR[loc-1]'

"
    l_LL(state,location)
    Left dominant eigenvector of the AL-AL transfermatrix
"
l_LL(state::InfiniteMPS{A},loc::Int=1) where A= isomorphism(Matrix{eltype(A)}, space(state.AL[loc],1),space(state.AL[loc],1))

"
    r_RR(state,location)
    Right dominant eigenvector of the AR-AR transfermatrix
"
r_RR(state::InfiniteMPS{A},loc::Int=length(state)) where A= isomorphism(Matrix{eltype(A)},domain(state.AR[loc]),domain(state.AR[loc]))

"
    r_RL(state,location)
    Right dominant eigenvector of the AR-AL transfermatrix
"
r_RL(state::InfiniteMPS,loc::Int=length(state)) = state.CR[loc]'

"
    r_LR(state,location)
    Right dominant eigenvector of the AL-AR transfermatrix
"
r_LR(state::InfiniteMPS,loc::Int=length(state)) = state.CR[loc]

"
    r_LL(state,location)
    Right dominant eigenvector of the AL-AL transfermatrix
"
r_LL(state::InfiniteMPS,loc::Int=length(state))= @tensor toret[-1;-2]:=state.CR[loc][-1,1]*conj(state.CR[loc][-2,1])

function TensorKit.dot(a::InfiniteMPS,b::InfiniteMPS;krylovdim = 30)
    init = TensorMap(rand,ComplexF64,space(a.AL[1],1),space(b.AL[1],1))
    num = lcm(length(a),length(b))
    (vals,vecs,convhist) = eigsolve(x->transfer_left(x,b.AL[1:num],a.AL[1:num]),init,1,:LM,Arnoldi(krylovdim=krylovdim))
    convhist.converged == 0 && @info "dot mps not converged"
    return vals[1]
end


text:
using dopingFVM
using Test

@testset "dopingFVM.jl" begin
    # Write your tests here.
end


text:
using NonstationaryProcesses
using Plots

lorenz = lorenzSim(
    X0 = [0.0, -0.01, 9.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = (10.0, 28.0, 8/3), # Sprott's recomendation
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.01,
    tmax = 1000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-12))

𝔩𝔬𝔯𝔢𝔫𝔷 = synchronisephase(lorenz, rampInterval, (0.0, 500.0, 250.0, 1000.0))()

a = animatespectrum(𝔩𝔬𝔯𝔢𝔫𝔷, downsample=100, trail=5000, colorgradient=cgrad([:black, :crimson]), phasogram=true, nperseg=1000, dpi=100)
gif(a, "./corruptlorenzsym.gif", fps=24)


text:
"""
    hermit(X, Y, Yp)

Computes interpolating polynomial passing through the points whoose coordinates are provided in `X[]` and `Y[]`,
as well as the values of differentiated function in `Yp[]`, using Hermite interpolation method.
Returns coefficiants of the Hermite polynomial.

    hermit(X::Array{Float64}, Y::Array::Float64, Yp::Array{Float64})

### OUTPUT:
  * `P` - coefficiants of the interpolating polynomial of the type `Array{Float64}`,
"""
function hermit(X::Array{Float64}, Y::Array{Float64}, Yp::Array{Float64})
    n = length(X)
    Z = zeros(n*2)
    Q = zeros(2*n+1, 2*n+1)
    for i in 1:n
        Z[2*i-1] = X[i]
        Z[2*i] = X[i]
        Q[2*i-1, 1] = Y[i]
        Q[2*i, 1] = Y[i]
        Q[2*i, 2] = Yp[i]

        if i > 1
            Q[2*i, 2] = (Q[2*i, 1] - Q[2*i-1, 1])/(Z[2*i] - Z[2*i-1])
        end
    end

    for i in 3:2*n+2
        for j in 3:i
            Q[i, j] = (Q[i, j-1] - Q[i-1, j-1])/(Z[i] - Z[(i+1)-j])
        end
    end

    P = []

    for i in i:2*n+1
        push!(P, Q[i, i])
    end

    return P
end


text:
@testset "ROMEO compile tests" begin

path = tempname()
compile_romeo(path; audit=true)

@test !isempty(readdir(path))

phasefile = joinpath(pwd(), "data", "small", "Phase.nii")
magfile = joinpath(pwd(), "data", "small", "Mag.nii")

function test_romeo(args)
    file = tempname()
    args = [args..., "-o", file]
    name = "romeo" * (Sys.iswindows() ? ".exe" : "")
    romeofile = joinpath(path, "bin", name)
    @test isfile(romeofile)
    cmd = `$romeofile $args`
    @test success(run(cmd))
end

args = [phasefile, "-m", magfile]
test_romeo(args)
end


text:
#!/usr/local/bin/julia

function bubbleSort(array::Array{Int64, 1})
    array_length = length(array)
    for iteration in 1:array_length
        for item in 1:array_length-iteration-1
            if array[item] > array[item+1]
                array[item], array[item+1] = array[item+1], array[item]
            end
        end
    end
    return array
end

# julia> a = [12,24,48,11,23,45,9,7,5,78]
# julia> bubbleSort(a)
# 10-element Array{Int64,1}:
#   5
#   7
#   9
#  11
#  12
#  23
#  24
#  45
#  48
#  78


text:
export kpow, kf1, kse, kseα, γe, rquad

"""`logpowk(logx, α)` `p(x)^(1/α)`

Compu

# Inputs
- `α ∈ [1, Inf]` - temperature 1 -> p(x)
- `logx` - logdensity - pₓ(x)

# Returns :
"""
logkpow(logx, α) = (1/α) * logx

"Real+ -> [0, 1]"
kf1(x, β = 0.0001) = x / (x + β)
kf1β(β) = d -> kf1(d, β)

"(Log) Squared exponential kernel `α = 1/2l^2`, higher α is lower temperature  "
kse(d, α = 1000.0) = - α * d
kseα(α) = d -> kse(d, α)

"γ exponential"
γe(x, l = 1, γ = 1) = exp(-(x/l)^γ)

"Rational Quadratic"
rquad(x, l = 1, α = 1) = (1 + (x^2)/(2*α*l^2))^(-α)

"Power law kernel"
kpow(d, α = 1.0, k = 2) = -k * log(d) + log(α)

burr(x, c = 1, k = 40) =  log(c) + log(k) +  (c - 1) * log(x) - (k + 1)*log(1 + x^c)

"""
`logtemps(n, k = 10)`
`n` logarithmically spaced temperatures from log(1) (ground state) .. "
"""
logtemps(n, k = 10) = exp.(k * range(-2.0, stop = 1.0, length = n))

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect