text:
mutable struct TextureInfo
    index::Int
    texCoord::Union{Nothing,Int}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function TextureInfo()
        obj = new()
        obj.texCoord = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::TextureInfo, sym::Symbol)
    x = getfield(obj, sym)
    sym === :texCoord && x === nothing && return 0
    return x
end

function Base.setproperty!(obj::TextureInfo, sym::Symbol, x)
    if sym === :index
        x ≥ 0 || throw(ArgumentError("the index of the texture should be ≥ 0"))
    elseif sym === :texCoord && x !== nothing
        x ≥ 0 || throw(ArgumentError("the set index of texture's TEXCOORD attribute used for texture coordinate mapping should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{TextureInfo}) = JSON3.Mutable()
StructTypes.omitempties(::Type{TextureInfo}) = (:texCoord, :extensions, :extras)

mutable struct NormalTextureInfo
    index::Int
    texCoord::Union{Nothing,Int}
    scale::Union{Nothing,Cfloat}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function NormalTextureInfo()
        obj = new()
        obj.texCoord = nothing
        obj.scale = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::NormalTextureInfo, sym::Symbol)
    x = getfield(obj, sym)
    sym === :texCoord && x === nothing && return 0
    sym === :scale && x === nothing && return 1
    return x
end

function Base.setproperty!(obj::NormalTextureInfo, sym::Symbol, x)
    if sym === :index
        x ≥ 0 || throw(ArgumentError("the index of the texture should be ≥ 0"))
    elseif sym === :texCoord && x !== nothing
        x ≥ 0 || throw(ArgumentError("the set index of texture's TEXCOORD attribute used for texture coordinate mapping should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{NormalTextureInfo}) = JSON3.Mutable()
StructTypes.omitempties(::Type{NormalTextureInfo}) = (:texCoord, :scale, :extensions, :extras)

mutable struct OcclusionTextureInfo
    index::Int
    texCoord::Union{Nothing,Int}
    strength::Union{Nothing,Cfloat}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function OcclusionTextureInfo()
        obj = new()
        obj.texCoord = nothing
        obj.strength = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::OcclusionTextureInfo, sym::Symbol)
    x = getfield(obj, sym)
    sym === :texCoord && x === nothing && return 0
    sym === :strength && x === nothing && return 1
    return x
end

function Base.setproperty!(obj::OcclusionTextureInfo, sym::Symbol, x)
    if sym === :index
        x ≥ 0 || throw(ArgumentError("the index of the texture should be ≥ 0"))
    elseif sym === :texCoord && x !== nothing
        x ≥ 0 || throw(ArgumentError("the set index of texture's TEXCOORD attribute used for texture coordinate mapping should be ≥ 0"))
    elseif sym === :strength && x !== nothing
        0 ≤ strength ≤ 1 || throw(ArgumentError("strength should be should be ≥ 0 and ≤ 1"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{OcclusionTextureInfo}) = JSON3.Mutable()
StructTypes.omitempties(::Type{OcclusionTextureInfo}) = (:texCoord, :strength, :extensions, :extras)

mutable struct PBRMetallicRoughness
    baseColorFactor::Union{Nothing,Vector{Cfloat}}
    baseColorTexture::Union{Nothing,TextureInfo}
    metallicFactor::Union{Nothing,Cfloat}
    roughnessFactor::Union{Nothing,Cfloat}
    metallicRoughnessTexture::Union{Nothing,TextureInfo}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function PBRMetallicRoughness()
        obj = new()
        obj.baseColorFactor = nothing
        obj.baseColorTexture = nothing
        obj.metallicFactor = nothing
        obj.roughnessFactor = nothing
        obj.metallicRoughnessTexture = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::PBRMetallicRoughness, sym::Symbol)
    x = getfield(obj, sym)
    sym === :baseColorFactor && x === nothing && return Cfloat[1,1,1,1]
    sym === :metallicFactor && x === nothing && return 1
    sym === :roughnessFactor && x === nothing && return 1
    return x
end

function Base.setproperty!(obj::PBRMetallicRoughness, sym::Symbol, x)
    if sym === :baseColorFactor && x !== nothing
        length(x) == 4 || throw(ArgumentError("the material's base color factor should exactly contain 4 values(RGBA)"))
        all(0 .≤ x .≤ 1) || throw(ArgumentError("the material's base color factor value should ≥ 0 and ≤ 1"))
    elseif sym === :metallicFactor && x !== nothing
        0 ≤ x ≤ 1 || throw(ArgumentError("the metalness of the material factor should ≥ 0 and ≤ 1"))
    elseif sym === :roughnessFactor && x !== nothing
        0 ≤ roughnessFactor ≤ 1 || throw(ArgumentError("the roughness of the material factor should ≥ 0 and ≤ 1"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{PBRMetallicRoughness}) = JSON3.Mutable()
StructTypes.omitempties(::Type{PBRMetallicRoughness}) = (:baseColorFactor, :baseColorTexture, :metallicFactor, :roughnessFactor, :metallicRoughnessTexture, :extensions, :extras)

mutable struct Material
    pbrMetallicRoughness::Union{Nothing,PBRMetallicRoughness}
    normalTexture::Union{Nothing,NormalTextureInfo}
    occlusionTexture::Union{Nothing,OcclusionTextureInfo}
    emissiveTexture::Union{Nothing,TextureInfo}
    emissiveFactor::Union{Nothing,Vector{Cfloat}}
    alphaMode::Union{Nothing,String}
    alphaCutoff::Union{Nothing,Cfloat}
    doubleSided::Union{Nothing,Bool}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Material()
        obj = new()
        obj.pbrMetallicRoughness = nothing
        obj.normalTexture = nothing
        obj.occlusionTexture = nothing
        obj.emissiveTexture = nothing
        obj.emissiveFactor = nothing
        obj.alphaMode = nothing
        obj.alphaCutoff = nothing
        obj.doubleSided = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::Material, sym::Symbol)
    x = getfield(obj, sym)
    sym === :emissiveFactor && x === nothing && return Cfloat[0,0,0]
    sym === :alphaMode && x === nothing && return "OPAQUE"
    sym === :alphaCutoff && x === nothing && return 0.5
    sym === :doubleSided && x === nothing && return false
    return x
end

function Base.setproperty!(obj::Material, sym::Symbol, x)
    if sym === :emissiveFactor && x !== nothing
        length(x) == 3 || throw(ArgumentError("the RGB components of the emissive color of the material should exactly contain 3 values"))
        all(0 .≤ x .≤ 1) || throw(ArgumentError("the emissive color should be ≥ 0 and ≤ 1"))
    elseif sym === :alphaMode && x !== nothing
        x === "OPAQUE" || x === "MASK" || x === "BLEND" ||
            throw(ArgumentError("""the material's alpha rendering mode should be one of: "OPAQUE", "MASK", "BLEND" """))
    elseif sym === :alphaCutoff && x !== nothing
        x ≥ 0 || throw(ArgumentError("alphaCutoff should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Material}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Material}) = (:pbrMetallicRoughness, :normalTexture, :occlusionTexture, :emissiveTexture, :emissiveFactor, :alphaMode, :alphaCutoff, :doubleSided, :name, :extensions, :extras)


text:
"""
    read(f::FortranFile [, spec [, spec [, ...]]])
    read(f::FortranFile, rec=N, [, spec [, spec [, ...]]])

Read data from a `FortranFile`. Like the READ statement in Fortran, this
reads a complete record, regardless of how man `spec`s are given. Each
`spec` can be:
* a `DataType` for scalar values; e.g. `Int32`, `Float64`, `FString{10}`
* a tuple of `DataType` and one or more integers, for reading arrays of
  the given size; e.g. `(Int32,4,2)` reads an `Array{Int32}` of size (4,2)
* a tuple of `DataType` and a tuple of integers, as an alternative way
  of reading arrays; e.g. `(Int32,(4,2))` does the same as the previous one
* an array, for reading into pre-allocated arrays; `DataType` and size
  of the array are implied through its Julia type.

For direct-access files, the number of the record to be read must be
specified with the `rec` keyword (N=1 for the first record).

Return value:
* if no `spec` is given: `nothing` (the record is skipped over)
* if one `spec` is given: the scalar or array requested
* if more `spec`s are given: a tuple of the scalars and arrays requested, in order
"""
function Base.read( f::FortranFile, specs...)
   record = Record(f)
   result = fread(record, specs...)
   close(record)
   return result
end

function Base.read( f::FortranFile{DirectAccess}, specs... ; rec::Integer=0 )
   if rec == 0
      fthrow("direct-access files require specifying the record to be read (use rec keyword argument)")
   end
   record = Record(f, rec)
   result = fread(record, specs...)
   close(record)
   return result
end

function fread( rec::Record )
   return nothing
end

function fread( rec::Record, spec )
   data = read_spec(rec, spec)
   return data
end

function fread( rec::Record, specs... )
   data = map( spec->read_spec(rec,spec), specs)
   return data
end

# workaround for "does not support byte I/O"
function read_spec( rec::Record, spec::Type{Int8} )
   b = read_spec(rec, (Int8,1))
   return b[1]
end

function read_spec( rec::Record, spec::Type{T} ) where {T}
   rec.convert.onread( read(rec, spec) )::T
end

function read_spec( rec::Record, spec::Array{T,N} ) where {T,N}
   arr = read!(rec, spec)::Array{T,N}
   map!(rec.convert.onread, arr, arr)
   return arr
end

function read_spec(rec::Record, spec::Tuple{DataType,I} ) where {I<:Integer}
   T,n = spec
   read_spec(rec, Array{T}(undef, n))::Array{T,1}
end

function read_spec( rec::Record, spec::Tuple{DataType, Vararg{Integer,N}} ) where {N}
   T = spec[1]
   sz = spec[2:end]
   read_spec(rec, Array{T}(undef, sz...))::Array{T,N}
end

function read_spec( rec::Record, spec::Tuple{DataType, Tuple{Vararg{Integer,N}}} ) where {N}
   T,sz = spec
   read_spec(rec, Array{T}(undef, sz...))::Array{T,N}
end

ext: jl


text:
# to use this, `using Revise; using GEAS; includet("this file")`
# then `ped, snp = two_random_generations(dir)`
# use `ped, snp` for your test.

using JLD2, Serialization, DataFrames
import GEAS:sim_QTL, create_storage, breeding, alleles2gt, snp_blup
"""
    function two_random_generations()
---
This is to generate a random mating two-generation pedigree to test about SNP-BLUP.
One-trait Results genotypes, pedigree and QTL info are written to `dir`.
"""
function two_random_generations()
    par = begin
        nsib = 10
        pp = Dict(:nSire=> 100,
                  :nDam => 200, # male:female = 1:2
                  :nSib => 2nsib,
                  :nC7e => nsib, # number of sib to be challenged
                  :nG8n => 1,    # number of generations
                  :nQTL => [100, 100],
                  :h²   => (a = .5, b = .5),
                  :p8h  => .5,  # percentage-of-death after challenge
                  :e19e => 1.,  # edit_successsful_rate
                  :b4y  => false,
                  :fix  => false,
                  :w4t  => 2    # weight on the binary trait EBV
                  )
        (; pp...)               # named tuple.  contents as above
    end
    @load "dat/run/base.jld2" base
    qtl = sim_QTL(base, par.nQTL...)
    ped, snp = create_storage(base, par, qtl)
    base, ped, snp
end

function test_trg()
    base, ped, snp = two_random_generations()
    g = alleles2gt(snp.prd)
    p = ped.prd.p7e
    snp_blup(g, p, .5, base.twop)
end


text:
## evaluate.jl - evaluation of graph

function remember_size!(g::AbstractExGraph, nd::ExNode)
    rsizes = @get_or_create(g.ctx, :rsizes, Dict{Symbol,Any}())
    buff_exprs = @get_or_create(g.ctx, :buff_exprs, Dict{Symbol, Any}())
    val = getvalue(nd)
    if isa(val, AbstractArray) || isa(val, Number)
        sz = size(val)
        rsizes[varname(nd)] = sz
        if isa(val, Array)
            T = eltype(val)
            buff_expr = :(zeros($T, $sz))
        else
            T = typeof(val)
            buff_expr = :(zero($T))
        end
        buff_exprs[varname(nd)] = buff_expr
    end
end

function remember_size!(g::AbstractExGraph, nd::ExNode{:tuple}) end

function isconv(nd::ExNode)
    depwarn_eingraph(:isconv)
    idxs = nd |> getexpr |> get_indices |> flatten
    return any(i -> isa(i, Expr), idxs)
end

function mk_eval_expr(g::ExGraph, nd::ExNode)
    dep_nodes = [g[dep] for dep in dependencies(nd) if haskey(g, dep)]
    deps_vals = Dict((varname(nd), getvalue(nd)) for nd in dep_nodes)
    evex = Expr(:block)
    for dep in unique(keys(deps_vals))
        val = deps_vals[dep]
        push!(evex.args, :(local $dep = $val))
    end
    codegen = eval_codegen(@get(g.ctx, :codegen, VectorCodeGen()))
    push!(evex.args, generate_code(codegen, g, nd))
    push!(evex.args, varname(nd))
    return evex
end

function mk_eval_expr(g::ExGraph, nd::ExNode{:ctor})
    dep_nodes = [g[dep] for dep in dependencies(nd) if haskey(g, dep)]
    deps_vals = Dict((varname(nd), getvalue(nd)) for nd in dep_nodes)
    evex = Expr(:block)
    for dep in unique(keys(deps_vals))
        val = deps_vals[dep]
        push!(evex.args, :(local $dep = $val))
    end
    push!(evex.args, to_expr_kw(nd))
    push!(evex.args, varname(nd))
    return evex
end

# function mk_eval_expr(g::ExGraph, nd::ExNode{:ctor})
#     dep_nodes = [g[dep] for dep in dependencies(nd) if haskey(g, dep)]
#     deps_vals = [(varname(nd), getvalue(nd)) for nd in dep_nodes]
#     eval_ex = Expr(:block, Expr(:let, Expr(:block)))
#     block = eval_ex.args[1].args[1]
#     for (dep, val) in unique(deps_vals)
#         push!(block.args, :(local $dep = $val))
#     end
#     push!(block.args, to_expr_kw(nd))
#     push!(block.args, varname(nd))
#     return eval_ex
# end

"""
Evaluate node, i.e. fill its `val` by evaluating node's expression using
values of its dependencies.
"""
function evaluate!(g::ExGraph, nd::ExNode{:constant}; force=false)
    ex = getexpr(nd)
    if getvalue(nd) == nothing && (isa(ex, Symbol) || isa(ex, Expr))
        # constant expression - need to evaluate it
        evex = mk_eval_expr(g, nd)
        val = Core.eval(g.ctx[:mod], evex)
        setvalue!(nd, val)
    end
    remember_size!(g, nd)
    return getvalue(nd)
end

function evaluate!(g::AbstractExGraph, nd::ExNode{:input}; force=false)
    remember_size!(g, nd)
    return getvalue(nd)
end

function evaluate!(g::AbstractExGraph, nd::ExNode{:(=)}; force=false)
    if (!force && getvalue(nd) != nothing) return getvalue(nd) end
    dep = dependencies(nd)[1]
    evaluate!(g, g[dep]; force=false)
    evex = mk_eval_expr(g, nd)
    setvalue!(nd, Core.eval(g.ctx[:mod], evex))
    remember_size!(g, nd)
    return getvalue(nd)
end

# nd::Union{ExNode{:call}, ExNode{:bcast}, ExNode{:ctor},
# ExNode{:tuple}, ExNode{:opaque}}
function evaluate!(g::AbstractExGraph, nd::ExNode; force=false)
    if (!force && getvalue(nd) != nothing) return getvalue(nd) end
    deps = dependencies(nd)
    for dep in deps
        # if dep is not in graph, consider it a global constant (like π)
        if haskey(g.idx, dep)
            evaluate!(g, g[dep]; force=false)  # force affects only current node
        end
    end
    evex = mk_eval_expr(g, nd)
    setvalue!(nd, Core.eval(g.ctx[:mod], evex))
    remember_size!(g, nd)
    return getvalue(nd)
end

function evaluate!(g::AbstractExGraph, nd::ExNode{:field}; force=false)
    if (!force && getvalue(nd) != nothing) return getvalue(nd) end
    depnd = g[dependencies(nd)[1]]
    obj = getvalue(depnd)
    fld = getexpr(nd).args[2].value
    val = getfield(obj, fld)
    setvalue!(nd, val)
    remember_size!(g, nd)
    return getvalue(nd)
end

evaluate!(g::AbstractExGraph, name::Symbol; force=false) = evaluate!(g, g[name]; force=force)

function evaluate!(g::AbstractExGraph; force=false)
    for i=1:length(g)
        try
            evaluate!(g, g[i]; force=force)
        catch e
            mod = @get(g.ctx, :mod, nothing)
            @info("Failed to evaluate node (in module $mod): $(g[i])")
            throw(e)
        end
    end
    return getvalue(g[end])
end


text:
using DataDeps

register(DataDep(
  "echocardiogram",
  "http://archive.ics.uci.edu/ml/datasets/Echocardiogram",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/echocardiogram/echocardiogram.data",
  "3d7dc4b3daef837f1b3b3c042e561697d6add718a2d6694de824faf1f1e340d4",
  post_fetch_method=(path -> begin
    # There is a misplaced leading comma on line 50, we should get rid of it
    lines = readlines(path)
    lines[50] = lines[50][2:end]
    write(path, join(lines, '\n'))

    UCIData.process_dataset(path,
        target_index=13,
        feature_indices=3:9,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# Index types

"""
    ConstraintIndex{F, S}

A type-safe wrapper for `Int64` for use in referencing `F`-in-`S` constraints in
a model.
The parameter `F` is the type of the function in the constraint, and the
parameter `S` is the type of set in the constraint. To allow for deletion,
indices need not be consecutive. Indices within a constraint type (i.e. `F`-in-`S`)
must be unique, but non-unique indices across different constraint types are allowed.
If `F` is [`SingleVariable`](@ref) then the index is equal to the index of the
variable. That is for an `index::ConstraintIndex{SingleVariable}`, we always
have
```julia
index.value == MOI.get(model, MOI.ConstraintFunction(), index).variable.value
```
"""
struct ConstraintIndex{F, S}
    value::Int64
end

"""
    VariableIndex

A type-safe wrapper for `Int64` for use in referencing variables in a model.
To allow for deletion, indices need not be consecutive.
"""
struct VariableIndex
    value::Int64
end

# The default hash is slow. It's important for the performance of dictionaries
# of VariableIndices to define our own.
# https://github.com/JuliaLang/julia/issues/10208
Base.hash(v::VariableIndex, h::UInt) = hash(v.value, h)

# No need to define isequal because the default matches our implementation of
# hash.

const Index = Union{ConstraintIndex,VariableIndex}

"""
    struct InvalidIndex{IndexType<:Index} <: Exception
        index::IndexType
    end

An error indicating that the index `index` is invalid.
"""
struct InvalidIndex{IndexType<:Index} <: Exception
    index::IndexType
end

function Base.showerror(io::IO, err::InvalidIndex)
    print("The index $(err.index) is invalid. Note that an index becomes invalid after it has been deleted.")
end

"""
    is_valid(model::ModelLike, index::Index)::Bool

Return a `Bool` indicating whether this index refers to a valid object in the model `model`.
"""
is_valid(model::ModelLike, ref::Index) = false

"""
    throw_if_not_valid(model::ModelLike, index::Index)

Throw an `InvalidIndex(index)` error if `MOI.is_valid(model, index)` returns
`false`.
"""
function throw_if_not_valid(model::ModelLike, index::Index)
    if !is_valid(model, index)
        throw(InvalidIndex(index))
    end
end

"""
    struct DeleteNotAllowed{IndexType <: Index} <: NotAllowedError
        index::IndexType
        message::String
    end

An error indicating that the index `index` cannot be deleted.
"""
struct DeleteNotAllowed{IndexType <: Index} <: NotAllowedError
    index::IndexType
    message::String
end
DeleteNotAllowed(index::Index) = DeleteNotAllowed(index, "")

function operation_name(err::DeleteNotAllowed)
    return "Deleting the index $(err.index)"
end

"""
    delete(model::ModelLike, index::Index)

Delete the referenced object from the model. Throw [`DeleteNotAllowed`](@ref) if
if `index` cannot be deleted.

The following modifications also take effect if `Index` is [`VariableIndex`](@ref):
* If `index` used in the objective function, it is removed from the function,
  i.e., it is substituted for zero.
* For each `func`-in-`set` constraint of the model:
  - If `func isa SingleVariable` and `func.variable == index` then the
    constraint is deleted.
  - If `func isa VectorOfVariables` and `index in func.variable` then
    * if `length(func.variable) == 1` is one, the constraint is deleted;
    * if `length(func.variable) > 1` and `supports_dimension_update(set)` then
      then the variable is removed from `func` and `set` is replaced by
      `update_dimension(set, MOI.dimension(set) - 1)`.
    * Otherwise, a [`DeleteNotAllowed`](@ref) error is thrown.
  - Otherwise, the variable is removed from `func`, i.e., it is substituted for
    zero.
"""
delete(model::ModelLike, index::Index) = throw(DeleteNotAllowed(index))

"""
    delete{R}(model::ModelLike, indices::Vector{R<:Index})

Delete the referenced objects in the vector `indices` from the model.
It may be assumed that `R` is a concrete type.
"""
function delete(model::ModelLike, indices::Vector{<:Index})
    for index in indices
        delete(model, index)
    end
end


text:
export NeuralNetworkApproximator

using Flux

struct NeuralNetworkApproximator{T,M,O,P} <: AbstractApproximator
    model::M
    optimizer::O
    params::P
end

"""
    NeuralNetworkApproximator(;kwargs)

Use a DNN model for value estimation.

# Keyword arguments

- `model`, a Flux based DNN model.
- `optimizer`
- `parameters=params(model)`
- `kind=Q_APPROXIMATOR`, specify the type of model.
"""
function NeuralNetworkApproximator(;
    model::M,
    optimizer::O,
    parameters::P = params(model),
    kind = Q_APPROXIMATOR,
) where {M,O,P}
    NeuralNetworkApproximator{kind,M,O,P}(model, optimizer, parameters)
end

device(app::NeuralNetworkApproximator) = device(app.model)

Flux.params(app::NeuralNetworkApproximator) = app.params

(app::NeuralNetworkApproximator)(s::AbstractArray) = app.model(s)
(app::NeuralNetworkApproximator{Q_APPROXIMATOR})(s::AbstractArray, a::Int) = app.model(s)[a]
(app::NeuralNetworkApproximator{HYBRID_APPROXIMATOR})(s::AbstractArray, ::Val{:Q}) =
    app.model(s, Val(:Q))
(app::NeuralNetworkApproximator{HYBRID_APPROXIMATOR})(s::AbstractArray, ::Val{:V}) =
    app.model(s, Val(:V))
(app::NeuralNetworkApproximator{HYBRID_APPROXIMATOR})(s::AbstractArray, a::Int) =
    app.model(s, Val(:Q))[a]

RLBase.batch_estimate(app::NeuralNetworkApproximator, states::AbstractArray) =
    app.model(states)

RLBase.batch_estimate(
    app::NeuralNetworkApproximator{HYBRID_APPROXIMATOR},
    states::AbstractArray,
    ::Val{:Q},
) = app.model(states, Val(:Q))

RLBase.batch_estimate(
    app::NeuralNetworkApproximator{HYBRID_APPROXIMATOR},
    states::AbstractArray,
    ::Val{:V},
) = app.model(states, Val(:V))

RLBase.update!(app::NeuralNetworkApproximator, gs) =
    Flux.Optimise.update!(app.optimizer, app.params, gs)

Base.copyto!(dest::NeuralNetworkApproximator, src::NeuralNetworkApproximator) =
    Flux.loadparams!(dest.model, src.params)


text:
Pkg.clone("git://github.com/adriantorrie/EodDataTestXml.jl.git")
tests = ["eod_utils_external", "eod_utils_internal"]

for t in tests
    fpath = "$t.jl"
    @printf("running %s ...\n", fpath)
    include(fpath)
end


text:
using Test
using JuMP
using SemialgebraicSets
using MultivariateMoments
using PolyJuMP
using DynamicPolynomials

function zero_polynomial_in_algebraic_set_test(optimizer,
                                               config::MOI.Test.TestConfig)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)
    @variable(model, β ≤ 1)

    @polyvar x y
    cref = @constraint(model, α * x - β * y == 0, domain = @set x == y)

    @objective(model, Max, α)
    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 1.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 1.0 atol=atol rtol=rtol
    @test value(β) ≈ 1.0 atol=atol rtol=rtol
    @test value(UpperBoundRef(β)) ≈ 1.0 atol=atol rtol=rtol

    @test dual_status(model) == MOI.FEASIBLE_POINT
    @test dual(UpperBoundRef(β)) ≈ -1.0 atol=atol rtol=rtol

    μ = dual(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 2
    @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == x
    @test moment_value(moments(μ)[2]) ≈ -1.0 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y

    μ = moments(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 1
    @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == y

    F = MOI.VectorAffineFunction{Float64}
    S = PolyJuMP.ZeroPolynomialSet{typeof(@set x == y), MonomialBasis,
                                   monomialtype(x), monovectype(x)}
    @test MOI.get(model, MOI.ListOfConstraints()) == [
        (MOI.SingleVariable, MOI.LessThan{Float64}), (F, S)]
    @testset "Delete" begin
        ST = PolyJuMP.ZeroPolynomialSet{FullSpace,MonomialBasis,Monomial{true},
                                        MonomialVector{true}}
        test_delete_bridge(model, cref, 2, ((F, MOI.Zeros, 0), (F, ST, 0)))
    end
end

linear_tests["zero_polynomial_in_algebraic_set"] = zero_polynomial_in_algebraic_set_test


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for ray.jl

@testset "constructor" begin
    ray1 = Ray(Point(1f0, 2f0, 3f0), Vec(5f0, 4f0, -1f0))
    ray2 = Ray(Point(1f0, 2f0, 3f0), Vec(5f0, 4f0, -1f0))
    ray3 = Ray(Point(5f0, 1f0, 4f0), Vec(3f0, 9f0, 4f0))

    @test ray1 ≈ ray2
    @test !(ray1 ≈ ray3)
end

@testset "methods" begin
    ray = Ray(Point(1f0, 2f0, 4f0), Vec(4f0, 2f0, 1f0))
    @test ray(0f0) ≈ ray.origin
    @test ray(1f0) ≈ Point(5f0, 4f0, 5f0)
    @test ray(2f0) ≈ Point(9f0, 6f0, 6f0)

    ray = Ray(Point(1f0, 2f0, 3f0), Vec(1f0, 2f0, 3f0))
    m = Transformation([1f0 2f0 3f0 4f0;
                        5f0 6f0 7f0 8f0;
                        9f0 9f0 8f0 7f0;
                        0f0 0f0 0f0 1f0],
                       [-3.75f0 2.75f0 -1f0 0f0;
                        5.75f0 -4.75f0 2f0 1f0;
                        -2.25f0 2.25f0 -1f0 -2f0;
                        0f0 0f0 0f0 1f0])
    @test m * ray ≈ Ray(Point(18f0, 46f0, 58f0), Vec(14f0, 38f0, 51f0))
end

@testset "transform" begin
    ray = Ray(Point(1f0, 2f0, 3f0), Vec(6f0, 5f0, 4f0))
    transformation = translation(Vec(10f0, 11f0, 12f0)) * rotationX(π/2)
    transformed = transformation * ray

    @test transformed.origin ≈ Point(11f0, 8f0, 14f0)
    @test transformed.dir ≈ Vec(6f0, -4f0, 5f0)
end


text:
const tri_table = ( (1,9,4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,9,4,10,9,2,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,2,3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,3,11,1,3,10,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,9,4,3,11,9,11,10,9,-1,-1,-1,-1,-1,-1),
                    (4,12,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,12,3,9,12,1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,10,1,3,4,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,12,3,2,10,12,10,9,12,-1,-1,-1,-1,-1,-1),
                    (4,11,2,12,11,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,11,2,1,9,11,9,12,11,-1,-1,-1,-1,-1,-1),
                    (4,10,1,4,12,10,12,11,10,-1,-1,-1,-1,-1,-1),
                    (10,9,11,11,9,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,8,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,4,1,8,4,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,2,10,9,5,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,2,10,5,8,2,8,4,2,-1,-1,-1,-1,-1,-1),
                    (2,3,11,9,5,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,5,8,4,1,5,2,3,11,-1,-1,-1,-1,-1,-1),
                    (10,3,11,10,1,3,9,5,8,-1,-1,-1,-1,-1,-1),
                    (3,11,10,3,10,8,3,8,4,8,10,5,-1,-1,-1),
                    (9,5,8,4,12,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,5,8,12,3,5,3,1,5,-1,-1,-1,-1,-1,-1),
                    (10,1,2,9,5,8,3,4,12,-1,-1,-1,-1,-1,-1),
                    (5,8,12,10,5,12,10,12,3,10,3,2,-1,-1,-1),
                    (4,11,2,4,12,11,8,9,5,-1,-1,-1,-1,-1,-1),
                    (2,12,11,2,5,12,2,1,5,8,12,5,-1,-1,-1),
                    (5,8,9,10,1,12,10,12,11,12,1,4,-1,-1,-1),
                    (5,8,12,5,12,10,10,12,11,-1,-1,-1,-1,-1,-1),
                    (10,6,5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,6,5,1,9,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,6,5,2,6,1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,6,5,9,4,6,4,2,6,-1,-1,-1,-1,-1,-1),
                    (2,3,11,10,6,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,1,9,2,3,11,5,10,6,-1,-1,-1,-1,-1,-1),
                    (6,3,11,6,5,3,5,1,3,-1,-1,-1,-1,-1,-1),
                    (3,11,6,4,3,6,4,6,5,4,5,9,-1,-1,-1),
                    (10,6,5,3,4,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,12,3,1,9,12,5,10,6,-1,-1,-1,-1,-1,-1),
                    (1,6,5,1,2,6,3,4,12,-1,-1,-1,-1,-1,-1),
                    (3,2,6,3,6,9,3,9,12,5,9,6,-1,-1,-1),
                    (11,4,12,11,2,4,10,6,5,-1,-1,-1,-1,-1,-1),
                    (5,10,6,1,9,2,9,11,2,9,12,11,-1,-1,-1),
                    (6,5,1,6,1,12,6,12,11,12,1,4,-1,-1,-1),
                    (6,5,9,6,9,11,11,9,12,-1,-1,-1,-1,-1,-1),
                    (10,8,9,6,8,10,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,4,1,10,6,4,6,8,4,-1,-1,-1,-1,-1,-1),
                    (1,8,9,1,2,8,2,6,8,-1,-1,-1,-1,-1,-1),
                    (2,6,4,4,6,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,8,9,10,6,8,11,2,3,-1,-1,-1,-1,-1,-1),
                    (11,2,3,10,6,1,6,4,1,6,8,4,-1,-1,-1),
                    (9,1,3,9,3,6,9,6,8,11,6,3,-1,-1,-1),
                    (3,11,6,3,6,4,4,6,8,-1,-1,-1,-1,-1,-1),
                    (8,10,6,8,9,10,4,12,3,-1,-1,-1,-1,-1,-1),
                    (10,6,8,10,8,3,10,3,1,3,8,12,-1,-1,-1),
                    (3,4,12,1,2,9,2,8,9,2,6,8,-1,-1,-1),
                    (12,3,2,12,2,8,8,2,6,-1,-1,-1,-1,-1,-1),
                    (10,6,9,9,6,8,11,2,4,11,4,12,-1,-1,-1),
                    (6,8,1,6,1,10,8,12,1,2,1,11,12,11,1),
                    (12,11,1,12,1,4,11,6,1,9,1,8,6,8,1),
                    (12,11,6,8,12,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (11,7,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,6,11,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,1,2,6,11,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,9,4,2,10,9,6,11,7,-1,-1,-1,-1,-1,-1),
                    (2,7,6,3,7,2,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,7,6,2,3,7,4,1,9,-1,-1,-1,-1,-1,-1),
                    (10,7,6,10,1,7,1,3,7,-1,-1,-1,-1,-1,-1),
                    (6,10,9,6,9,3,6,3,7,4,3,9,-1,-1,-1),
                    (3,4,12,11,7,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,1,9,12,3,1,11,7,6,-1,-1,-1,-1,-1,-1),
                    (1,2,10,3,4,12,6,11,7,-1,-1,-1,-1,-1,-1),
                    (6,11,7,2,10,3,10,12,3,10,9,12,-1,-1,-1),
                    (7,4,12,7,6,4,6,2,4,-1,-1,-1,-1,-1,-1),
                    (1,9,12,1,12,6,1,6,2,6,12,7,-1,-1,-1),
                    (4,12,7,1,4,7,1,7,6,1,6,10,-1,-1,-1),
                    (7,6,10,7,10,12,12,10,9,-1,-1,-1,-1,-1,-1),
                    (6,11,7,5,8,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,4,1,5,8,4,7,6,11,-1,-1,-1,-1,-1,-1),
                    (2,10,1,6,11,7,9,5,8,-1,-1,-1,-1,-1,-1),
                    (11,7,6,2,10,8,2,8,4,8,10,5,-1,-1,-1),
                    (7,2,3,7,6,2,5,8,9,-1,-1,-1,-1,-1,-1),
                    (2,3,6,6,3,7,4,1,5,4,5,8,-1,-1,-1),
                    (9,5,8,10,1,6,1,7,6,1,3,7,-1,-1,-1),
                    (8,4,10,8,10,5,4,3,10,6,10,7,3,7,10),
                    (4,12,3,8,9,5,11,7,6,-1,-1,-1,-1,-1,-1),
                    (6,11,7,5,8,3,5,3,1,3,8,12,-1,-1,-1),
                    (1,2,10,5,8,9,3,4,12,6,11,7,-1,-1,-1),
                    (10,3,2,10,12,3,10,5,12,8,12,5,6,11,7),
                    (9,5,8,4,12,6,4,6,2,6,12,7,-1,-1,-1),
                    (6,2,12,6,12,7,2,1,12,8,12,5,1,5,12),
                    (1,6,10,1,7,6,1,4,7,12,7,4,9,5,8),
                    (7,6,10,7,10,12,5,8,10,8,12,10,-1,-1,-1),
                    (11,5,10,7,5,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,11,7,5,10,11,1,9,4,-1,-1,-1,-1,-1,-1),
                    (11,1,2,11,7,1,7,5,1,-1,-1,-1,-1,-1,-1),
                    (9,4,2,9,2,7,9,7,5,7,2,11,-1,-1,-1),
                    (2,5,10,2,3,5,3,7,5,-1,-1,-1,-1,-1,-1),
                    (4,1,9,2,3,10,3,5,10,3,7,5,-1,-1,-1),
                    (1,3,5,5,3,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,4,3,9,3,5,5,3,7,-1,-1,-1,-1,-1,-1),
                    (11,5,10,11,7,5,12,3,4,-1,-1,-1,-1,-1,-1),
                    (1,9,3,3,9,12,5,10,11,5,11,7,-1,-1,-1),
                    (4,12,3,1,2,7,1,7,5,7,2,11,-1,-1,-1),
                    (7,5,2,7,2,11,5,9,2,3,2,12,9,12,2),
                    (10,7,5,10,4,7,10,2,4,12,7,4,-1,-1,-1),
                    (9,12,2,9,2,1,12,7,2,10,2,5,7,5,2),
                    (4,12,7,4,7,1,1,7,5,-1,-1,-1,-1,-1,-1),
                    (7,5,9,12,7,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,11,7,8,9,11,9,10,11,-1,-1,-1,-1,-1,-1),
                    (1,8,4,1,11,8,1,10,11,7,8,11,-1,-1,-1),
                    (11,7,8,2,11,8,2,8,9,2,9,1,-1,-1,-1),
                    (11,7,8,11,8,2,2,8,4,-1,-1,-1,-1,-1,-1),
                    (2,3,7,2,7,9,2,9,10,9,7,8,-1,-1,-1),
                    (3,7,10,3,10,2,7,8,10,1,10,4,8,4,10),
                    (8,9,1,8,1,7,7,1,3,-1,-1,-1,-1,-1,-1),
                    (8,4,3,7,8,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,4,12,11,7,9,11,9,10,9,7,8,-1,-1,-1),
                    (3,1,8,3,8,12,1,10,8,7,8,11,10,11,8),
                    (2,9,1,2,8,9,2,11,8,7,8,11,3,4,12),
                    (12,3,2,12,2,8,11,7,2,7,8,2,-1,-1,-1),
                    (9,10,7,9,7,8,10,2,7,12,7,4,2,4,7),
                    (1,10,2,12,7,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,9,1,8,1,7,4,12,1,12,7,1,-1,-1,-1),
                    (8,12,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,7,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,1,9,12,8,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,2,10,12,8,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,2,10,9,4,2,12,8,7,-1,-1,-1,-1,-1,-1),
                    (11,2,3,7,12,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,3,11,4,1,9,7,12,8,-1,-1,-1,-1,-1,-1),
                    (3,10,1,3,11,10,7,12,8,-1,-1,-1,-1,-1,-1),
                    (7,12,8,3,11,4,11,9,4,11,10,9,-1,-1,-1),
                    (8,3,4,7,3,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,1,9,8,7,1,7,3,1,-1,-1,-1,-1,-1,-1),
                    (3,8,7,3,4,8,1,2,10,-1,-1,-1,-1,-1,-1),
                    (2,7,3,2,9,7,2,10,9,9,8,7,-1,-1,-1),
                    (11,8,7,11,2,8,2,4,8,-1,-1,-1,-1,-1,-1),
                    (11,8,7,2,8,11,2,9,8,2,1,9,-1,-1,-1),
                    (1,4,8,1,8,11,1,11,10,7,11,8,-1,-1,-1),
                    (8,7,11,8,11,9,9,11,10,-1,-1,-1,-1,-1,-1),
                    (7,9,5,12,9,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,7,12,4,1,7,1,5,7,-1,-1,-1,-1,-1,-1),
                    (9,7,12,9,5,7,10,1,2,-1,-1,-1,-1,-1,-1),
                    (10,5,7,10,7,4,10,4,2,12,4,7,-1,-1,-1),
                    (7,9,5,7,12,9,3,11,2,-1,-1,-1,-1,-1,-1),
                    (2,3,11,4,1,12,1,7,12,1,5,7,-1,-1,-1),
                    (5,12,9,5,7,12,1,3,10,3,11,10,-1,-1,-1),
                    (11,10,4,11,4,3,10,5,4,12,4,7,5,7,4),
                    (9,3,4,9,5,3,5,7,3,-1,-1,-1,-1,-1,-1),
                    (1,5,3,5,7,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,10,1,3,4,5,3,5,7,5,4,9,-1,-1,-1),
                    (2,10,5,2,5,3,3,5,7,-1,-1,-1,-1,-1,-1),
                    (9,2,4,9,7,2,9,5,7,7,11,2,-1,-1,-1),
                    (11,2,1,11,1,7,7,1,5,-1,-1,-1,-1,-1,-1),
                    (5,7,4,5,4,9,7,11,4,1,4,10,11,10,4),
                    (11,10,5,7,11,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,10,6,8,7,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,5,10,6,12,8,7,-1,-1,-1,-1,-1,-1),
                    (6,1,2,6,5,1,8,7,12,-1,-1,-1,-1,-1,-1),
                    (12,8,7,9,4,5,4,6,5,4,2,6,-1,-1,-1),
                    (10,6,5,11,2,3,8,7,12,-1,-1,-1,-1,-1,-1),
                    (7,12,8,2,3,11,1,9,4,5,10,6,-1,-1,-1),
                    (8,7,12,6,5,11,5,3,11,5,1,3,-1,-1,-1),
                    (4,5,9,4,6,5,4,3,6,11,6,3,12,8,7),
                    (8,3,4,8,7,3,6,5,10,-1,-1,-1,-1,-1,-1),
                    (10,6,5,1,9,7,1,7,3,7,9,8,-1,-1,-1),
                    (4,7,3,4,8,7,2,6,1,6,5,1,-1,-1,-1),
                    (7,3,9,7,9,8,3,2,9,5,9,6,2,6,9),
                    (10,6,5,11,2,7,2,8,7,2,4,8,-1,-1,-1),
                    (2,7,11,2,8,7,2,1,8,9,8,1,10,6,5),
                    (5,1,11,5,11,6,1,4,11,7,11,8,4,8,11),
                    (8,7,11,8,11,9,6,5,11,5,9,11,-1,-1,-1),
                    (7,10,6,7,12,10,12,9,10,-1,-1,-1,-1,-1,-1),
                    (4,7,12,1,7,4,1,6,7,1,10,6,-1,-1,-1),
                    (1,12,9,1,6,12,1,2,6,6,7,12,-1,-1,-1),
                    (7,12,4,7,4,6,6,4,2,-1,-1,-1,-1,-1,-1),
                    (2,3,11,10,6,12,10,12,9,12,6,7,-1,-1,-1),
                    (1,12,4,1,7,12,1,10,7,6,7,10,2,3,11),
                    (12,9,6,12,6,7,9,1,6,11,6,3,1,3,6),
                    (7,12,4,7,4,6,3,11,4,11,6,4,-1,-1,-1),
                    (6,9,10,6,3,9,6,7,3,4,9,3,-1,-1,-1),
                    (10,6,7,10,7,1,1,7,3,-1,-1,-1,-1,-1,-1),
                    (2,6,9,2,9,1,6,7,9,4,9,3,7,3,9),
                    (2,6,7,3,2,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,4,7,2,7,11,4,9,7,6,7,10,9,10,7),
                    (11,2,1,11,1,7,10,6,1,6,7,1,-1,-1,-1),
                    (1,4,9,6,7,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (11,6,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,6,11,8,6,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,6,11,12,8,6,9,4,1,-1,-1,-1,-1,-1,-1),
                    (6,12,8,6,11,12,2,10,1,-1,-1,-1,-1,-1,-1),
                    (11,8,6,11,12,8,10,9,2,9,4,2,-1,-1,-1),
                    (12,2,3,12,8,2,8,6,2,-1,-1,-1,-1,-1,-1),
                    (1,9,4,2,3,8,2,8,6,8,3,12,-1,-1,-1),
                    (10,8,6,10,3,8,10,1,3,3,12,8,-1,-1,-1),
                    (8,6,3,8,3,12,6,10,3,4,3,9,10,9,3),
                    (3,6,11,3,4,6,4,8,6,-1,-1,-1,-1,-1,-1),
                    (9,3,1,9,6,3,9,8,6,11,3,6,-1,-1,-1),
                    (10,1,2,6,11,4,6,4,8,4,11,3,-1,-1,-1),
                    (10,9,3,10,3,2,9,8,3,11,3,6,8,6,3),
                    (2,4,6,4,8,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,8,1,8,2,2,8,6,-1,-1,-1,-1,-1,-1),
                    (10,1,4,10,4,6,6,4,8,-1,-1,-1,-1,-1,-1),
                    (10,9,8,6,10,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (6,9,5,6,11,9,11,12,9,-1,-1,-1,-1,-1,-1),
                    (6,1,5,6,12,1,6,11,12,12,4,1,-1,-1,-1),
                    (1,2,10,9,5,11,9,11,12,11,5,6,-1,-1,-1),
                    (11,12,5,11,5,6,12,4,5,10,5,2,4,2,5),
                    (3,6,2,3,9,6,3,12,9,5,6,9,-1,-1,-1),
                    (1,5,12,1,12,4,5,6,12,3,12,2,6,2,12),
                    (1,3,6,1,6,10,3,12,6,5,6,9,12,9,6),
                    (10,5,6,3,12,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,6,11,4,6,3,4,5,6,4,9,5,-1,-1,-1),
                    (6,11,3,6,3,5,5,3,1,-1,-1,-1,-1,-1,-1),
                    (4,11,3,4,6,11,4,9,6,5,6,9,1,2,10),
                    (6,11,3,6,3,5,2,10,3,10,5,3,-1,-1,-1),
                    (9,5,6,9,6,4,4,6,2,-1,-1,-1,-1,-1,-1),
                    (1,5,6,2,1,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,5,6,9,6,4,10,1,6,1,4,6,-1,-1,-1),
                    (10,5,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,12,8,5,10,12,10,11,12,-1,-1,-1,-1,-1,-1),
                    (1,9,4,5,10,8,10,12,8,10,11,12,-1,-1,-1),
                    (2,11,12,2,12,5,2,5,1,8,5,12,-1,-1,-1),
                    (4,2,5,4,5,9,2,11,5,8,5,12,11,12,5),
                    (5,12,8,10,12,5,10,3,12,10,2,3,-1,-1,-1),
                    (10,8,5,10,12,8,10,2,12,3,12,2,1,9,4),
                    (12,8,5,12,5,3,3,5,1,-1,-1,-1,-1,-1,-1),
                    (12,8,5,12,5,3,9,4,5,4,3,5,-1,-1,-1),
                    (3,10,11,3,8,10,3,4,8,8,5,10,-1,-1,-1),
                    (10,11,8,10,8,5,11,3,8,9,8,1,3,1,8),
                    (4,8,11,4,11,3,8,5,11,2,11,1,5,1,11),
                    (2,11,3,9,8,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,10,2,5,2,8,8,2,4,-1,-1,-1,-1,-1,-1),
                    (5,10,2,5,2,8,1,9,2,9,8,2,-1,-1,-1),
                    (5,1,4,8,5,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,9,8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,11,9,11,12,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,1,10,4,10,12,12,10,11,-1,-1,-1,-1,-1,-1),
                    (1,2,11,1,11,9,9,11,12,-1,-1,-1,-1,-1,-1),
                    (4,2,11,12,4,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,3,12,2,12,10,10,12,9,-1,-1,-1,-1,-1,-1),
                    (4,1,10,4,10,12,2,3,10,3,12,10,-1,-1,-1),
                    (1,3,12,9,1,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,3,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,4,9,3,9,11,11,9,10,-1,-1,-1,-1,-1,-1),
                    (10,11,3,1,10,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,4,9,3,9,11,1,2,9,2,11,9,-1,-1,-1),
                    (2,11,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,4,9,10,2,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,10,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,4,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1))

v_table = []
for trimap in tri_table
    a = tuple(filter(x-> x>0, unique!([trimap...]))...)
    b = zeros(UInt8, 12)
    for i in eachindex(a)
        b[i] = UInt8(a[i])
    end
    push!(v_table, tuple(b...))
end

for v in v_table
    println(string(v)*",")
end


text:
import Benchmark

Benchmark.run_paper_benchmark([1_000, 5_000, 10_000, 50_000])

max_stars_repo_name: jump-dev/MOIPaperBenchmarks

max_issues_repo_name: jump-dev/MOIPaperBenchmarks

max_forks_repo_name: jump-dev/MOIPaperBenchmarks


text:
using RecipesBase

@recipe function f(tree::RTree; usebrl = true, showlabels = true, shownodes = false)

    A = RootedTrees._preplot_tree(tree, usebrl)

    grid --> false
    legend --> false
    yaxis := false
    yticks := false
    xaxis --> false
    xticks --> false

    @series begin
        label := "Branches"
        linecolor --> :black
        seriestype := :path
        x = A.xbranch
        y = A.ybranch

        x, y
    end

    if shownodes
        @series begin
            label := "Nodes"
            seriestype := :scatter
            x = A.xnode
            y = A.ynode

            x, y
        end
    end

    if showlabels
        treelabels = getlabel.(getnodes(tree))

        # Find suitable label offset `dx`
        dx = extrema(A.xnode)
        dx = 0.025 * (dx[2] - dx[1])

        @series begin
            label := nothing
            hover := nothing
            markeralpha := 0
            markersize := 0
            series_annotations := [(lab, 10, :left) for lab in treelabels]
            textcolor --> :black
            seriestype := :scatter
            x = A.xnode .+ dx
            y = A.ynode

            x, y
        end
    end

    ()
end

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
Base.@kwdef struct Experiment
    gradient::NamedTuple
    btpde::Union{NamedTuple,Nothing} = nothing
    btpde_midpoint::Union{NamedTuple,Nothing} = nothing
    hadc::Union{NamedTuple,Nothing} = nothing
    mf::Union{NamedTuple,Nothing} = nothing
    analytical::Union{NamedTuple,Nothing} = nothing
    karger::Union{NamedTuple,Nothing} = nothing
end


text:
module MyIterativeSolvers

isdefined(:__precompile__) && __precompile__(true)

include("common.jl")
include("krylov.jl")

#Specialized factorizations
include("factorization.jl")

#Linear solvers
include("stationary.jl")
include("cg.jl")
include("gmres.jl")
include("chebyshev.jl")

#Eigensolvers
include("simple.jl")
include("lanczos.jl")

#SVD solvers
include("lanczos-svd.jl")
include("lanczos-svd-tr.jl")

#Least-squares
include("lsqr.jl")

#Randomized algorithms
include("rlinalg.jl")
include("rsvd.jl")
include("rsvd_fnkz.jl")

end


text:
function main()
    println(join(["Hello world!"], " "))
    println(join(["Hello", "world!"], " "))
end

main()


text:
module NoncommutativeGraphs

import Base.==

using DocStringExtensions
using Subspaces
using Convex, SCS, LinearAlgebra
using Random, RandomMatrices
using LightGraphs
using Compat
import Base.show

export AlgebraShape
export S0Graph
export create_S0_S1
export random_S0Graph, empty_S0Graph, complement, vertex_graph, forget_S0
export from_block_spaces, get_block_spaces
export block_expander
export random_S0_unitary, random_S0_density
export random_S1_unitary, random_S1_density

export Ψ
export dsw_schur, dsw_schur2
export dsw, dsw_via_complement

eye(n) = Matrix(1.0*I, (n,n))

"""
The structure of a finite dimensional C*-algebra.

- For example, `[1 2; 3 4]` corresponds to S₀ = M₁⊗I₂ ⊕ M₃⊗I₄.
- For an n-dimensional non-commutative graph use `[1 n]` for S₀ = Iₙ.
- For an n-vertex classical graph use `ones(Integer, n, 2)` for S₀ = diagonals.
"""
AlgebraShape = Array{<:Integer, 2}

"""
    create_S0_S1(sig::AlgebraShape) -> Tuple{Subspace, Subspace}

Create a C*-algebra and its commutant with the given structure.
"""
function create_S0_S1(sig::AlgebraShape)
    blocks0 = []
    blocks1 = []
    @compat for row in eachrow(sig)
        if length(row) != 2
            throw(ArgumentError("row length must be 2"))
        end
        dA = row[1]
        dY = row[2]
        #println("dA=$dA dY=$dY")
        blk0 = kron(full_subspace(ComplexF64, (dA, dA)), Matrix((1.0+0.0im)*I, dY, dY))
        blk1 = kron(Matrix((1.0+0.0im)*I, dA, dA), full_subspace(ComplexF64, (dY, dY)))
        #println(blk0)
        push!(blocks0, blk0)
        push!(blocks1, blk1)
    end

    S0 = cat(blocks0..., dims=(1,2))
    S1 = cat(blocks1..., dims=(1,2))

    @assert I in S0
    @assert I in S1
    s0 = random_element(S0)
    s1 = random_element(S1)
    @assert (norm(s0 * s1 - s1 * s0) < 1e-9) "S0 and S1 don't commute"

    return S0, S1
end

"""
    S0Graph(sig::AlgebraShape, S::Subspace{ComplexF64, 2})

    S0Graph(g::AbstractGraph)

Represents an S₀-graph as defined in arxiv:1002.2514.

$(TYPEDFIELDS)
"""
struct S0Graph
    """Dimension of Hilbert space A such that S ⊆ L(A)"""
    n::Integer
    """Structure of C*-algebra S₀"""
    sig::AlgebraShape
    """Subspace that defines the graph"""
    S::Subspace{ComplexF64, 2}
    """C*-algebra S₀"""
    S0::Subspace{ComplexF64, 2}
    """Commutant of C*-algebra S₀"""
    S1::Subspace{ComplexF64, 2}
    """Block scaling array D from definition 23 of arxiv:2101.00162"""
    D::Array{Float64, 2}

    function S0Graph(sig::AlgebraShape, S::Subspace{ComplexF64, 2})
        S0, S1 = create_S0_S1(sig)
        S == S' || throw(DomainError("S is not an S0-graph"))
        S0 in S || throw(DomainError("S is not an S0-graph"))
        (S == S0 * S * S0) || throw(DomainError("S is not an S0-graph"))

        n = size(S0)[1]
        da_sizes = sig[:,1]
        dy_sizes = sig[:,2]
        n_sizes = da_sizes .* dy_sizes
        D = cat([ v*eye(n) for (n, v) in zip(n_sizes, dy_sizes ./ da_sizes) ]..., dims=(1,2))

        return new(n, sig, S, S0, S1, D)
    end

    function S0Graph(g::AbstractGraph)
        n = nv(g)
        S0 = Subspace([ (x=zeros(ComplexF64, n,n); x[i,i]=1; x) for i in 1:n ])
        S  = Subspace([ (x=zeros(ComplexF64, n,n); x[src(e),dst(e)]=1; x) for e in edges(g) ])
        S = S + S' + S0
        return S0Graph(ones(Int64, n, 2), S)
    end
end

function show(io::IO, g::S0Graph)
    print(io, "S0Graph{S0=$(g.sig) S=$(g.S)}")
end

"""
$(TYPEDSIGNATURES)

Returns the S₀-graph with S=S₀.
"""
vertex_graph(g::S0Graph) = S0Graph(g.sig, g.S0)

"""
$(TYPEDSIGNATURES)

Returns an S₀-graph with S₀=ℂI.
"""
forget_S0(g::S0Graph) = S0Graph([1 g.n], g.S)

"""
    random_S0Graph(sig::AlgebraShape) -> S0Graph

Creates a random S₀-graph with S₀ having the given structure.
"""
function random_S0Graph(sig::AlgebraShape)
    S0, S1 = create_S0_S1(sig)

    num_blocks = size(sig, 1)
    function block(col, row)
        da_col, dy_col = sig[col,:]
        da_row, dy_row = sig[row,:]
        ds = Integer(round(dy_row * dy_col / 2.0))
        F = full_subspace(ComplexF64, (da_row, da_col))
        if row == col
            R = random_hermitian_subspace(ComplexF64, ds, dy_row)
        elseif row > col
            R = random_subspace(ComplexF64, ds, (dy_row, dy_col))
        else
            R = empty_subspace(ComplexF64, (dy_row, dy_col))
        end
        return kron(F, R)
    end
    blocks = Array{Subspace{ComplexF64, 2}, 2}([
        block(col, row)
        for col in 1:num_blocks, row in 1:num_blocks
   ])

    S = hvcat(num_blocks, blocks...)
    S |= S'
    S |= S0

    return S0Graph(sig, S)
end

"""
    empty_S0Graph(sig::AlgebraShape) -> S0Graph

Creates an empty S₀-graph (i.e. S=S₀) with S₀ having the given structure.
"""
function empty_S0Graph(sig::AlgebraShape)
    S0, S1 = create_S0_S1(sig)
    return S0Graph(sig, S0)
end

"""
$(TYPEDSIGNATURES)

Returns the complement graph perp(S) + S₀.
"""
complement(g::S0Graph) = S0Graph(g.sig, perp(g.S) | g.S0)

function ==(a::S0Graph, b::S0Graph)
    return a.sig == b.sig && a.S == b.S
end

function get_block_spaces(g::S0Graph)
    num_blocks = size(g.sig, 1)
    da_sizes = g.sig[:,1]
    dy_sizes = g.sig[:,2]
    n_sizes = da_sizes .* dy_sizes

    blkspaces = Array{Subspace{ComplexF64}, 2}(undef, num_blocks, num_blocks)
    offseti = 0
    for blki in 1:num_blocks
        offsetj = 0
        for blkj in 1:num_blocks
            #@show [blki, blkj, offseti, offsetj]
            blkbasis = Array{Array{ComplexF64, 2}, 1}()
            for m in each_basis_element(g.S)
                blk = m[1+offseti:dy_sizes[blki]+offseti, 1+offsetj:dy_sizes[blkj]+offsetj]
                push!(blkbasis, blk)
            end
            blkspaces[blki, blkj] = Subspace(blkbasis)
            #println(blkspaces[blki, blkj])
            offsetj += n_sizes[blkj]
        end
        @assert offsetj == size(g.S)[2]
        offseti += n_sizes[blki]
    end
    @assert offseti == size(g.S)[1]

    return blkspaces
end

function from_block_spaces(sig::AlgebraShape, blkspaces::Array{Subspace{ComplexF64}, 2})
    S0, S1 = NoncommutativeGraphs.create_S0_S1(sig)

    num_blocks = size(sig, 1)
    function block(col, row)
        da_col, dy_col = sig[col,:]
        da_row, dy_row = sig[row,:]
        ds = Integer(round(sqrt(dy_row * dy_col) / 2.0))
        F = full_subspace(ComplexF64, (da_row, da_col))
        kron(F, blkspaces[row, col])
    end
    blocks = [
        block(col, row)
        for col in 1:num_blocks, row in 1:num_blocks
    ]

    S = hvcat(num_blocks, blocks...)
    S |= S'
    S |= S0

    return S0Graph(sig, S)
end

function block_expander(sig::AlgebraShape)
    function basis_mat(n, i)
        M = zeros(n*n)
        M[i] = 1
        return reshape(M, (n, n))
    end

    n = sum(prod(sig, dims=2))

    @compat J = cat([
        cat([ kron(eye(dA), basis_mat(dY, i)) for i in 1:dY^2 ]..., dims=3)
        for (dA, dY) in eachrow(sig)
    ]..., dims=(1,2,3))

    return reshape(J, (n^2, size(J)[3]))
end

function random_positive(n)
    U = rand(Haar(2), n)
    return Hermitian(U' * Diagonal(rand(n)) * U)
end

"""
Returns a random unitary in S₀.
"""
function random_S0_unitary(sig::AlgebraShape)
    @compat return cat([
        kron(rand(Haar(2), dA), eye(dY))
        for (dA, dY) in eachrow(sig)
    ]..., dims=(1,2))
end

"""
Returns a random density operator in S₀.
"""
function random_S0_density(sig::AlgebraShape)
    @compat ρ = cat([
        kron(random_positive(dA), eye(dY))
        for (dA, dY) in eachrow(sig)
    ]..., dims=(1,2))
    return ρ / tr(ρ)
end

"""
Returns a random unitary in the commutant of S₀.
"""
function random_S1_unitary(sig::AlgebraShape)
    @compat return cat([
        kron(eye(dA), rand(Haar(2), dY))
        for (dA, dY) in eachrow(sig)
    ]..., dims=(1,2))
end

"""
Returns a random density operator in the commutant of S₀.
"""
function random_S1_density(sig::AlgebraShape)
    @compat ρ = cat([
        kron(eye(dA), random_positive(dY))
        for (dA, dY) in eachrow(sig)
    ]..., dims=(1,2))
    return ρ / tr(ρ)
end

###############
### DSW solvers
###############

# FIXME Convex.jl doesn't support cat(args..., dims=(1,2)).  It should be added.
function diagcat(args::Convex.AbstractExprOrValue...)
    num_blocks = size(args, 1)
    return vcat([
        hcat([
            row == col ? args[row] : zeros(size(args[row], 1), size(args[col], 2))
            for col in 1:num_blocks
        ]...)
        for row in 1:num_blocks
    ]...)
end

"""
$(TYPEDSIGNATURES)

Block scaling superoperator Ψ from definition 23 of arxiv:2101.00162
"""
function Ψ(g::S0Graph, w::Union{AbstractArray{<:Number, 2}, Variable})
    n = g.n
    da_sizes = g.sig[:,1]
    dy_sizes = g.sig[:,2]
    n_sizes = da_sizes .* dy_sizes
    num_blocks = size(g.sig)[1]

    k = 0
    blocks = []
    for (dai, dyi) in zip(da_sizes, dy_sizes)
        ni = dai * dyi
        TrAi = partialtrace(w[k+1:k+ni, k+1:k+ni], 1, [dai; dyi])
        blk = dyi^-1 * kron(Array(1.0*I, dai, dai), TrAi)
        k += ni
        push!(blocks, blk)
    end
    @assert k == n
    out = diagcat(blocks...)
    @assert size(out) == (n, n)
    return out
end

"""
$(TYPEDSIGNATURES)

Schur complement form of weighted θ from theorem 14 of arxiv:2101.00162.

Returns λ, w, and Z variables (for Convex.jl) in a named tuple.

See also: [`dsw_schur2`](@ref).
"""
function dsw_schur(g::S0Graph)::NamedTuple{(:λ, :w, :Z), Tuple{Convex.AbstractExpr, Convex.AbstractExpr, Convex.AbstractExpr}}
    n = size(g.S)[1]

    Z = sum(kron(m, ComplexVariable(n, n)) for m in hermitian_basis(g.S))
    # slow:
    #Z = ComplexVariable(n^2, n^2)
    #add_constraint!(Z, Z in kron(g.S, full_subspace(ComplexF64, (n, n))))
    # slow:
    #SB = kron(g.S, full_subspace(ComplexF64, (n, n)))
    #Z = variable_in_space(SB)

    λ = Variable()
    wt = partialtrace(Z, 1, [n; n])
    wv = reshape(wt, n*n, 1)

    add_constraint!(λ, [ λ  wv' ; wv  Z ] ⪰ 0)

    return (λ=λ, w=transpose(wt), Z=Z)
end

"""
$(TYPEDSIGNATURES)

Schur complement form of weighted θ from theorem 14 of arxiv:2101.00162, optimized for the
case S₀ ≠ ℂI, at the cost of w being constrained to S₁ (the commutant of S₀).

Returns λ, w, and Z variables (for Convex.jl) in a named tuple.

See also: [`dsw_schur2`](@ref).
"""
function dsw_schur2(g::S0Graph)::NamedTuple{(:λ, :w, :Z), Tuple{Convex.AbstractExpr, Convex.AbstractExpr, Convex.AbstractExpr}}
    da_sizes = g.sig[:,1]
    dy_sizes = g.sig[:,2]
    num_blocks = size(g.sig)[1]

    blkspaces = get_block_spaces(g)

    w_blocks = Array{Convex.AbstractExpr, 1}(undef, num_blocks)
    Z_blocks = Array{Convex.AbstractExpr, 2}(undef, num_blocks, num_blocks)
    for blki in 1:num_blocks
        for blkj in 1:num_blocks
            if blkj <= blki
                if dim(blkspaces[blki, blkj]) == 0
                    Z_blocks[blki, blkj] = zeros(dy_sizes[blki]^2, dy_sizes[blkj]^2)
                else
                    blkV = sum(kron(m, ComplexVariable(dy_sizes[blki], dy_sizes[blkj]))
                        for m in each_basis_element(blkspaces[blki, blkj]))
                    Z_blocks[blki, blkj] = blkV
                    # slow:
                    #SB = kron(blkspaces[blki, blkj], full_subspace(ComplexF64, (dy_sizes[blki], dy_sizes[blkj])))
                    #Z_blocks[blki, blkj] = variable_in_space(SB)
                end
            end
            if blkj == blki
                w_blocks[blki] = partialtrace(Z_blocks[blki, blkj], 1, [dy_sizes[blki], dy_sizes[blkj]])
            end
        end
    end

    for blki in 1:num_blocks
        for blkj in 1:num_blocks
            if blkj > blki
                Z_blocks[blki, blkj] = Z_blocks[blkj, blki]'
            end
            #@show blki, blkj
            #@show size(Z_blocks[blki, blkj])
        end
    end

    #@show [ size(wi) for wi in w_blocks ]

    λ = Variable()
    Z = vcat([ hcat([Z_blocks[i,j] for j in 1:num_blocks]...) for i in 1:num_blocks ]...)
    wv = vcat([ reshape(wi, dy_sizes[i]^2, 1) for (i,wi) in enumerate(w_blocks) ]...)
    #@show size(wv)
    #@show size(Z)

    add_constraint!(λ, [ λ  wv' ; wv  Z ] ⪰ 0)

    wt = diagcat([ kron(eye(da_sizes[i]), wi) for (i,wi) in enumerate(w_blocks) ]...)
    #@show size(wt)

    return (λ=λ, w=transpose(wt), Z=Z)
end

"""
$(TYPEDSIGNATURES)

Compute weighted θ using theorem 14 of arxiv:2101.00162.

Returns optimal λ, x, and Z values in a named tuple.
If `use_diag_optimization=true` (the default) then `x ⪰ w` and `x` is in the commutant
of S₀.  By theorem 29 of arxiv:2101.00162, θ(g, w) = θ(g, x).
"""
function dsw(g::S0Graph, w::AbstractArray{<:Number, 2}; use_diag_optimization=true, eps=1e-6, verbose=0)
    if use_diag_optimization
        λ, x, Z = dsw_schur2(g)
    else
        λ, x, Z = dsw_schur(g)
    end

    problem = minimize(λ, [x ⪰ w])
    solve!(problem, () -> SCS.Optimizer(verbose=verbose, eps=eps))
    return (λ=problem.optval, x=Hermitian(evaluate(x)), Z=Hermitian(evaluate(Z)))
end

"""
$(TYPEDSIGNATURES)

Compute weighted θ via the complement graph, using theorem 29 of arxiv:2101.00162.

θ(S, w) = max{ tr(w x) : x ⪰ 0, y = Ψ(x), θ(Sᶜ, y) ≤ 1 }

Returns optimal λ, x, y, and Z in a named tuple.

If w is in the commutant of S₀ then the weights w and y saturate the inequality in
theorem 32 of arxiv:2101.00162.
"""
function dsw_via_complement(g::S0Graph, w::AbstractArray{<:Number, 2}; use_diag_optimization=true, eps=1e-6, verbose=0)
    # max{ <w,x> : Ψ(S, x) ⪯ y, ϑ(S, y) ≤ 1, y ∈ S1 }
    # equal to:
    # max{ dsw(S0, √y * w * √y) : dsw(complement(S), y) <= 1 }
    if use_diag_optimization
        λ, y, Z = dsw_schur2(g)
    else
        λ, y, Z = dsw_schur(g)
    end
    x = HermitianSemidefinite(g.n, g.n)
    problem = maximize(real(tr(w * x')), [ λ <= 1, Ψ(g, x) == y ])
    solve!(problem, () -> SCS.Optimizer(verbose=verbose, eps=eps))
    return (λ=problem.optval, x=Hermitian(evaluate(x)), y=Hermitian(evaluate(y)), Z=Hermitian(evaluate(Z)))
end

end

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
"""
Integrate common transcoders with Formats.
"""
module FormatCodecs

using Formats
using CodecZlib
using CodecBzip2
using CodecXz

function __init__()
	Formats.addcoding("application/gzip")
	Formats.addextension("application/gzip", ".gz")
	Formats.addextension("application/gzip", ".gzip")
	Formats.setdecoder("application/gzip", GzipDecompressor)
	Formats.setencoder("application/gzip", GzipCompressor)

	Formats.addcoding("application/x-bzip2")
	Formats.addextension("application/x-bzip2", ".bz2")
	Formats.addextension("application/x-bzip2", ".bzip2")
	Formats.setdecoder("application/x-bzip2", Bzip2Decompressor)
	Formats.setencoder("application/x-bzip2", Bzip2Compressor)

	Formats.addcoding("application/x-xz")
	Formats.addextension("application/x-xz", ".xz")
	Formats.setdecoder("application/x-xz", XzDecompressor)
	Formats.setencoder("application/x-xz", XzCompressor)
end

end # module


text:
@testset "module evaluation" begin

function mock_output(str, mod = nothing)
    result_doc = mock_run(str; mod = mod)
    return result_doc.chunks[1].output
end

str = """
```julia
@__MODULE__
```
"""

# in sandbox
@test occursin(r"\#+WeaveSandBox[\#\d]+", mock_output(str))

# in Main
@test strip(mock_output(str, Main)) == "Main"

end # @testset "module evaluation"

@testset "clear_module!" begin

ary = rand(1000000)
size = Base.summarysize(ary)

# simple case
m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
Core.eval(m, :(a = $ary))
Weave.clear_module!(m)
@test Base.summarysize(m) < size

# recursive case
m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
Core.eval(m, :(
    module $(gensym(:WeaveTestSubModule))
    a = $ary
    end
))
Weave.clear_module!(m)
@test Base.summarysize(m) < size

# doesn't work with constants
m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
Core.eval(m, :(const a = $ary))
Weave.clear_module!(m)
@test_broken Base.summarysize(m) < size

julia_markdown_body = """
this is just to test the `out_path` option
"""

f_in = tempname()
f_out = tempname() * ".md"
write(f_in, julia_markdown_body)
f = weave(f_in; out_path=f_out)
@test isfile(f_out)

end # @testset "clear_module!"


text:
# This file contains some extra checks that don't fall into the other
# categories, for example because they're very platform-specific.

check_os_abi(::ObjectHandle, ::Platform, rest...; verbose::Bool = false, kwargs...) = true

function check_os_abi(oh::ELFHandle, ::FreeBSD, rest...; verbose::Bool = false, kwargs...)
    if oh.ei.osabi != 0x09
        # The dynamic loader should not have problems in this case, but the
        # linker may not appreciate.  Let the user know about this.
        if verbose
            @warn "OS/ABI is not set to FreeBSD (0x09) in the file header, this may be an issue at linking time"
        end
        return false
    end
    return true
end


text:
using MathML, EzXML, Symbolics

# const ns = namespace(node)
const sbml_lists = [
"listOfFunctionDefinitions",
"listOfUnitDefinitions",
"listOfCompartments",
"listOfSpecies",
"listOfParameters",
"listOfReactions/x:reaction/x:kineticLaw/x:listOfParameters",
"listOfReactions/x:reaction/x:kineticLaw/x:listOfLocalParameters",
"listOfInitialAssignments",
"listOfRules",
"listOfConstraints",
"listOfReactions",
"listOfEvents",
]

"""
give a fn, doc, or node of an sbml model and gather metadata

this gets used to build the MTK component
"""
function sbml_to_sysinfo end

function sbml_to_sysinfo(fn::AbstractString)
    !isfile(fn) && error("isfile(fn) must be true")
    sbml_to_sysinfo(readxml(fn))
end

function sbml_to_sysinfo(doc::EzXML.Document)
    sbml_to_sysinfo(doc.root)
end

function sbml_to_sysinfo(node::EzXML.Node)
    ns = namespace(node) # will this always work
    node.name != "sbml" && error("give an sbml node")
    arr = []
    for list_name in sbml_lists
        str = "/x:sbml/x:model/x:$(list_name)"
        l = findall(str, node, ["x" => ns])
        #use children or elements?
        ls = !isempty(l) ? mapreduce(elements, vcat, l) : EzXML.Node[]
        push!(arr, list_name => ls)
    end
    Dict(arr)
end

#=function sbml_to_sysinfo(node::EzXML.Node,promotelocalparameters::Bool)
    ns = namespace(node) # will this always work
    node.name != "sbml" && error("give an sbml node")
    arr = []
    for list_name in sbml_lists
        str = "/x:$(list_name)"
        l = findall(str, node, ["x" => ns])
        #use children or elements?
        ls = !isempty(l) ? mapreduce(elements, vcat, l) : nothing
        push!(arr, list_name => ls)
    end
    Dict(arr)
end=#

"doesn't handle the `constant=Bool` attribute or units"
function build_map(ps::Vector{EzXML.Node}, name, value)
    names = @. Symbol(getindex(ps, name))
    vals = @. Meta.parse(getindex(ps, value))
    names .=> vals
end

"given a list of ezxml nodes, return the union of all unique attributes of "
function uniquekeys(node_list)
    isnothing(node_list) ? nothing :
    unique(get_name.(reduce(vcat, attributes.(node_list))))
end

"create a `DataFrame`-able for the node list, returns pairs"
function to_table(node_list)
    isnothing(node_list) && return nothing
    ks = uniquekeys(node_list)
    arr = []
    for node in node_list
        r = []
        for k in ks
            x = haskey(node, k) ? node[k] : missing
            push!(r, x)
        end
        push!(arr, r)
    end
    ks .=> eachrow(reduce(hcat, arr))
end

get_name(n::EzXML.Node) = n.name

# function build_tables(fn)
#     ts = Dict()
#     sysinfo = sbml_to_sysinfo(fn)
#     for (k, v) in sysinfo
#         @show k v
#         t = isnothing(v) ? to_table(v) : nothing
#         ts[k] = t
#     end
#     ts
# end

function make_nice(d)
    ks, vs = keys(d), values(d)
    ts = to_table.(values(d))
    td = Dict(ks .=> ts)
    filter(x -> !isnothing(x.second), td)
end


text:
module TestPhenotypes

using Test
using TMLEEpistasis
using CSV
using DataFrames

@testset "Test prepare_phenotypes" begin
    # Test for continuous phenotypes with withdrawal list
    parsed_args = Dict(
        "phenotypes-file" => joinpath("data", "ukbb", "geneatlas", "quantitative_phenotypes.csv"),
        "bridge" => joinpath("data", "ukbb", "geneatlas", "bridge.csv"),
        "withdrawal-list" => joinpath("data", "ukbb", "withdrawal_list.csv"),
        "output" => "full_phenotypes.csv",
        "phenotypes-list" => nothing
    )
    TMLEEpistasis.prepare_phenotypes(parsed_args)

    phenotypes = CSV.File(parsed_args["output"]) |> DataFrame

    @test names(phenotypes) == ["SAMPLE_ID",
                                "1408-0.0",
                                "1777-0.0",
                                "1727-0.0",
                                "1548-0.0"]
    @test phenotypes.SAMPLE_ID == [1000023, 1000030]

    rm(parsed_args["output"])

    # Test for continuous phenotypes with withdrawal phenotypes list

    parsed_args = Dict(
        "phenotypes-file" => joinpath("data", "ukbb", "geneatlas", "binary_phenotypes.csv"),
        "bridge" => joinpath("data", "ukbb", "geneatlas", "bridge.csv"),
        "output" => "full_phenotypes.csv",
        "withdrawal-list" => nothing,
        "phenotypes-list" => joinpath("data", "ukbb", "phenotypes_list.csv")
    )
    TMLEEpistasis.prepare_phenotypes(parsed_args)

    phenotypes = CSV.File(parsed_args["output"]) |> DataFrame

    @test names(phenotypes) == ["SAMPLE_ID",
                                "clinical_c_G25",
                                "clinical_c_G20"]
    @test phenotypes.SAMPLE_ID == [1000190, 1000023, 1000014, 1000030]

    rm(parsed_args["output"])

end

@testset "Test tmle_phenotypes_batches" begin
    n_phenotypes = 10
    temp_file = "temp_phenotypes.csv"
    phenotypes_names = ["x$i" for i in 1:n_phenotypes]
    CSV.write(
        temp_file,
        DataFrame(rand(10, n_phenotypes+1), vcat("SAMPLE_ID", phenotypes_names))
        )

    # Test basic mode
    parsed_args = Dict(
        "phenotypes-file" => temp_file,
        "batch-size" => "3"
    )
    TMLEEpistasis.tmle_phenotypes_batches(parsed_args)

    expected_batches = [
        ["x1", "x2", "x3"],
        ["x4", "x5", "x6"],
        ["x7", "x8", "x9"],
        ["x10"]
    ]
    for batch in 1:4
        batch_file = TMLEEpistasis.batch_filename(batch)
        @test readlines(open(batch_file)) == expected_batches[batch]
        rm(batch_file)
    end

    # Test max mode
    parsed_args = Dict(
        "phenotypes-file" => temp_file,
        "batch-size" => "max"
    )
    TMLEEpistasis.tmle_phenotypes_batches(parsed_args)

    batch_file = TMLEEpistasis.batch_filename(1)
    @test readlines(open(batch_file)) == phenotypes_names
    rm(batch_file)
    rm(temp_file)

end

end;

true


text:
module Test_regular

using LinearAlgebra
using MatrixPencils
using Test

@testset "Regular Matrix Pencils Utilities" begin

@testset "_svdlikeAE" begin

for fast in (true, false)

A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = missing;
A = copy(A2); E = copy(E2); C = copy(C2); B = missing;
Q = nothing
Z = nothing

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test rE == 0 && rA22 == 0

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,2); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

Q = nothing
Z = nothing

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test rE == 0 && rA22 == 0

n = 9;
A2 = [
    1.     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E2 = [
    0.     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B2 = [
    -1.     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C2 = [
    0.     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
T = Float64;
Q = Matrix{T}(I,n,n)
Z = Matrix{T}(I,n,n)

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE == 7 && rA22 == 0

n = 9;

A2 = [
    1.     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E2 = [
    0.     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B2 = [
    -1.     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C2 = [
    0.     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];

T = Float64;
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

Q = Matrix{T}(I,n,n)
Z = Matrix{T}(I,n,n)

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE == 7 && rA22 == 0

# TG01ED EXAMPLE PROGRAM DATA
n = 4;
A2 = [
    -1.     0     0     3
     0     0     1     2
     1     1     0     4
     0     0     0     0];
 E2 = [
     1.     2     0     0
     0     1     0     1
     3     9     6     3
     0     0     2     0];
 B2 = [
     1.     0
     0     0
     0     1
     1     1];
 C2 = [
    -1.     0     1     0
     0     1    -1     1];

T = Float64;
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

Q = Matrix{T}(I,n,n)
Z = Matrix{T}(I,n,n)

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE == 3 && rA22 == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

Q = Matrix{T}(I,n,n)
Z = Matrix{T}(I,n,n)

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE == 3 && rA22 == 1

n = 4;
A2 = [
    -3.330669073875470e-16                         0                         0                         0
    -5.336424998587595e-17     4.999999999999997e-01                         0                         0
     7.882634225314346e-01     4.926646390821473e-02     1.017643702629568e+00                         0
    -7.804970910378840e-02    -5.073231091746246e-01    -3.602232247129060e-01     1.906644727694980e+00
     ];
 E2 = [
     9.708552809960075e-34    -4.715554847257422e-19     7.905694150420953e-01                         0
    -1.002546820843656e-17    -9.625586024907403e-20     9.682458365518573e-02     4.743416490252569e-01
    -9.878387353077869e-19     2.181405497860925e-19    -3.720759787739367e-01     4.673827146373172e-02
     1.017230340487078e-17     7.512741891975197e-20    -6.045704470706158e-02    -4.812889603890242e-01
     ];
 B2 = [
    -1.171404622872256e-16    -6.123724356957946e-01    -2.056086568456411e-17     2.775557561562891e-16    -2.844390262617513e-16
    -5.000000000000002e-01     1.250000000000000e-01    -5.000000000000003e-01     4.139921271247625e-17     4.999999999999999e-01
    -4.926646390821468e-02    -4.803480231050927e-01    -4.926646390821469e-02     7.882634225314342e-01     4.926646390821479e-02
    -4.975668955366508e-01    -7.804970910378836e-02    -4.975668955366509e-01    -7.804970910378847e-02    -5.073231091746250e-01
     ];
 C2 = [
     9.999999999999998e-01                         0                         0                         0
    -9.799563198300907e-17     9.999999999999998e-01                         0                         0
     ];

T = Float64;
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

Q = Matrix{T}(I,n,n)
Z = Matrix{T}(I,n,n)

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE == 2 && rA22 == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

Q = Matrix{T}(I,n,n)
Z = Matrix{T}(I,n,n)

@time rE, rA22  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE == 2 && rA22 == 1

for Ty in (Float64, Complex{Float64})

n = 10; m = 3; p = 2;
rE = 5;
rA22 = min(2,n-rE)
n2 = rE+rA22
n3 = n-n2
A2 = [rand(Ty,rE,n);
rand(Ty,rA22,n2) zeros(Ty,rA22,n3);
rand(Ty,n3,rE) zeros(Ty,n3,n-rE) ];
E2 = [rand(Ty,rE,rE) zeros(Ty,rE,n-rE);
zeros(Ty,n-rE,n)]
Q = qr(rand(Ty,n,n)).Q
Z = qr(rand(Ty,n,n)).Q
A2 = Q*A2*Z; E2 = Q*E2*Z;
B2 = rand(Ty,n,m); C2 = rand(Ty,p,n);

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);
Q = Matrix{Ty}(I,n,n)
Z = Matrix{Ty}(I,n,n)

@time rE1, rA221  = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast, atol1=1.e-7, atol2 = 1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C)  < sqrt(eps(1.))) &&
      rE1 == rE && rA221 == rA22

end
end

end

# test to enforce the call with lanv2
@testset "ordeigvals" begin
a = rand(2,2);
ev =  eigvals(a)
evord = ordeigvals(a)
@test sort(real(ev)) ≈ sort(real(evord)) && sort(imag(ev)) ≈ sort(imag(evord))

a = rand(2,2); e = triu(rand(2,2))
ev =  eigvals(a,e)
evord = ordeigvals(a,e)[1]
@test sort(real(ev)) ≈ sort(real(evord)) && sort(imag(ev)) ≈ sort(imag(evord))

end # ordeigvals

@testset "isregular" begin

M = zeros(0,0); N = zeros(0,0);
@test isregular(M, N)

M = zeros(3,0); N = zeros(3,0);
@test !isregular(M, N)

M = zeros(0,3); N = zeros(0,3);
@test !isregular(M, N)

M = zeros(1,1); N = ones(1,1);
@test isregular(M, N)

N = zeros(1,1); M = ones(1,1);
@test isregular(M, N)

M = [1 0;0 1]; N = [0 1; 0 0];
@test isregular(M, N)

M = [  22  34  31   31  17
        45  45  42   19  29
        39  47  49   26  34
        27  31  26   21  15
        38  44  44   24  30 ];

N = [   13  26  25  17  24
        31  46  40  26  37
        26  40  19  25  25
        16  25  27  14  23
        24  35  18  21  22  ];
@test !isregular(M, N)

Ty = Float64
for Ty in (Float64, Complex{Float64})

    abstol = sqrt(eps(one(real(Ty))))

    # given structure
    mr = 2; nr = 3; ni = 2; nf = 10; ml = 5; nl = 4;
    mM = mr+ni+nf+ml;
    nM = nr+ni+nf+nl;
    M = [rand(Ty,mr,nM);
         zeros(Ty,ni,nr) Matrix{Ty}(I,ni,ni) rand(Ty,ni,nf+nl);
         zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
         zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
    N = [rand(Ty,mr,nM);
         zeros(Ty,ni,nr) diagm(1 => ones(Ty,ni-1)) rand(Ty,ni,nf+nl);
         zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
         zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
    Q = qr(rand(Ty,mM,mM)).Q;
    Z = qr(rand(Ty,nM,nM)).Q;
    M = Q*M*Z;
    N = Q*N*Z;

    @test !isregular(M, N, atol1 = abstol, atol2 = abstol)

   # given structure
   mr = 0; nr = 0; ni = 4; nf = 10; ml = 0; nl = 0;
   mM = mr+ni+nf+ml;
   nM = nr+ni+nf+nl;
   M = [rand(Ty,mr,nM);
        zeros(Ty,ni,nr) Matrix{Ty}(I,ni,ni) rand(Ty,ni,nf+nl);
        zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
        zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
   N = [rand(Ty,mr,nM);
        zeros(Ty,ni,nr) diagm(1 => ones(Ty,ni-1)) rand(Ty,ni,nf+nl);
        zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
        zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
   Q = qr(rand(Ty,mM,mM)).Q;
   Z = qr(rand(Ty,nM,nM)).Q;
   M = Q*M*Z;
   N = Q*N*Z;

   @test isregular(M, N, atol1 = abstol, atol2 = abstol)
   γ = peigvals(M,N)[1][1]
   @test !isregular(M-γ*N, 0*N, atol1 = abstol, atol2 = abstol)
   @test !isregular(N, 0*N, atol1 = abstol, atol2 = abstol)
   @test !isregular(M, N, γ, atol = abstol)
   @test !isregular(M, N, Inf, atol = abstol)

end

end # isregular

@testset "isunimodular" begin

M = zeros(0,0); N = zeros(0,0);
@test isunimodular(M, N)

M = zeros(3,0); N = zeros(3,0);
@test !isunimodular(M, N)

M = zeros(0,3); N = zeros(0,3);
@test !isunimodular(M, N)

M = zeros(1,1); N = ones(1,1);
@test !isunimodular(M, N)

N = zeros(1,1); M = ones(1,1);
@test isunimodular(M, N)

M = [1 0;0 1]; N = [0 1; 0 0];
@test isregular(M, N)

M = [  22  34  31   31  17
        45  45  42   19  29
        39  47  49   26  34
        27  31  26   21  15
        38  44  44   24  30 ];

N = [   13  26  25  17  24
        31  46  40  26  37
        26  40  19  25  25
        16  25  27  14  23
        24  35  18  21  22  ];
@test !isunimodular(M, N)

Ty = Float64
for Ty in (Float64, Complex{Float64})

    abstol = sqrt(eps(one(real(Ty))))

    # given structure
    mr = 2; nr = 3; ni = 2; nf = 10; ml = 5; nl = 4;
    mM = mr+ni+nf+ml;
    nM = nr+ni+nf+nl;
    M = [rand(Ty,mr,nM);
         zeros(Ty,ni,nr) Matrix{Ty}(I,ni,ni) rand(Ty,ni,nf+nl);
         zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
         zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
    N = [rand(Ty,mr,nM);
         zeros(Ty,ni,nr) diagm(1 => ones(Ty,ni-1)) rand(Ty,ni,nf+nl);
         zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
         zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
    Q = qr(rand(Ty,mM,mM)).Q;
    Z = qr(rand(Ty,nM,nM)).Q;
    M = Q*M*Z;
    N = Q*N*Z;

    @test !isunimodular(M, N, atol1 = abstol, atol2 = abstol)

   # given structure
   mr = 0; nr = 0; ni = 10; nf = 0; ml = 0; nl = 0;
   mM = mr+ni+nf+ml;
   nM = nr+ni+nf+nl;
   M = [rand(Ty,mr,nM);
        zeros(Ty,ni,nr) Matrix{Ty}(I,ni,ni) rand(Ty,ni,nf+nl);
        zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
        zeros(ml,nr+ni+nf) rand(Ty,ml,nl)]
   N = [zeros(Ty,mr,nM);
        zeros(Ty,ni,nr) diagm(1 => ones(Ty,ni-1)) rand(Ty,ni,nf+nl);
        zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
        zeros(ml,nr+ni+nf) zeros(Ty,ml,nl)]
   Q = qr(rand(Ty,mM,mM)).Q;
   Z = qr(rand(Ty,nM,nM)).Q;
   M = Q*M*Z;
   N = Q*N*Z;

   @test isunimodular(M, N, atol1 = abstol, atol2 = abstol)
end

end # isunimodular

@testset "fisplit" begin
A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = missing;
A = copy(A2); E = copy(E2); C = copy(C2); B = missing;

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,0)

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,2); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,0)

A2 = [
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E2 = [
    0     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B2 = [
    -1     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C2 = [
    0     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C, finite_infinite=false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && istriu(A1) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [2, 2, 2, 1, 1, 1] && blkdims == (9,0)

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C, finite_infinite=true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && istriu(A1) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [1, 1, 1, 2, 2, 2] && blkdims == (0,9)

fast = true; finite_infinite = true; Ty = Float64

for fast in (true, false)

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); E2 = zeros(Ty,3,3); B2 = zeros(Ty,3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C,fast = fast, finite_infinite = true)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && istriu(A1) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [3] && blkdims == (0,3)

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C,fast = fast, finite_infinite = false)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && istriu(A1) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [3] && blkdims == (3,0)

A2 = rand(3,3); E2 = triu(rand(Ty,3,3),1); B2 = zeros(3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C,fast = fast, finite_infinite = true)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [1] && blkdims == (2,1)

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C,fast = fast, finite_infinite = false)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [1] && blkdims == (1,2)

A2 = zeros(3,3); E2 = rand(3,3); B2 = zeros(Ty,3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C,fast = fast, finite_infinite = true)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [] && blkdims == (3,0)

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = fisplit(A,E,B,C,fast = fast, finite_infinite = false)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,3)

A2 = rand(3,3); E2 = rand(3,3); B2 = zeros(Ty,3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, C1, B1, Q, Z, ν, blkdims  = fisplit(A',E',C',B',fast = fast, finite_infinite = true)
@test norm(Q'*A2'*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2'*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(B2'*Z-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(Q'*C2'-C1)   < sqrt(eps(1.))) &&
      ν == [] && blkdims == (3,0)

@time A1, E1, C1, B1, Q, Z, ν, blkdims  = fisplit(A',E',C',B',fast = fast, finite_infinite = false)
@test norm(Q'*A2'*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2'*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(B2'*Z-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(Q'*C2'-C1)   < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,3)

end
end
end # fisplit

fast = true; finite_infinite = true; Ty = Float64

for fast in (true, false)

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); E2 = zeros(Ty,3,3); B2 = zeros(Ty,3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = sfisplit(A,E,B,C,fast = fast, finite_infinite = true)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && istriu(A1) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == Int64[] && blkdims == (3, 0, 0)

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = sfisplit(A,E,B,C,fast = fast, finite_infinite = false)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && istriu(A1) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == Int64[] && blkdims == (0,0,3)

A2 = rand(3,3); E2 = triu(rand(Ty,3,3),1); B2 = zeros(3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = sfisplit(A,E,B,C,fast = fast, finite_infinite = true)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == Int64[] && blkdims == (1,2,0)

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = sfisplit(A,E,B,C,fast = fast, finite_infinite = false)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == Int64[] && blkdims == (0, 2, 1)

A2 = zeros(3,3); E2 = rand(3,3); B2 = zeros(Ty,3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = sfisplit(A,E,B,C,fast = fast, finite_infinite = true)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,3,0)

@time A1, E1, B1, C1, Q, Z, ν, blkdims  = sfisplit(A,E,B,C,fast = fast, finite_infinite = false)
@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(C2*Z-C1)  < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,3,0)

A2 = rand(3,3); E2 = rand(3,3); B2 = zeros(Ty,3,2); C2 = zeros(4,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time A1, E1, C1, B1, Q, Z, ν, blkdims  = sfisplit(A',E',C',B',fast = fast, finite_infinite = true)
@test norm(Q'*A2'*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2'*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(B2'*Z-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(Q'*C2'-C1)   < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,3,0)

@time A1, E1, C1, B1, Q, Z, ν, blkdims  = sfisplit(A',E',C',B',fast = fast, finite_infinite = false)
@test norm(Q'*A2'*Z-A1) < sqrt(eps(1.)) &&
      norm(Q'*E2'*Z-E1) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(B2'*Z-B1) < sqrt(eps(1.))) &&
      (ismissing(C) || norm(Q'*C2'-C1)   < sqrt(eps(1.))) &&
      ν == [] && blkdims == (0,3,0)

end

end # sfisplit

end # regular testset
end # module