text:
# So we can access "binsize" before a `CompactHeap` is instantiated.
struct PowerOfTwo
    val::Int
end

poweroftwo(x::PowerOfTwo) = x
poweroftwo(x::Integer) = PowerOfTwo(x)

getbin(x::PowerOfTwo, sz) = ((sz - 1) >> x.val) + 1
binsize(x::PowerOfTwo, bin) = bin << x.val

#####
##### Iterate over a vector of `Blocks` as a vector of `FreelistPtr`.
#####

struct FreelistIterator{T}
    lists::Vector{T}
end

Base.length(x::FreelistIterator) = length(x.lists)
Base.pointer(x::FreelistIterator, i) = pointer(x.lists, i)
function Base.iterate(x::FreelistIterator, i = 1)
    i > length(x) && return nothing
    return FreelistPtr(pointer(x, i)), i + 1
end

#####
##### CompactHeap
#####

mutable struct CompactHeap <: AbstractHeap
    minallocation::PowerOfTwo

    # Where is and what kind of memory do we have?
    base::Ptr{Nothing}
    len::Int
    allocated::Int
    pool::Pool

    # Maintaining freelist status
    status::FindNextTree
    freelists::Vector{Block}
    canalloc::Bool

    # Allocator token.
    # Can be used to preserve `base` as long as the heap is still alive.
    token::Any
end

topointer(ptr::Ptr) = ptr
topointer(A::AbstractArray) = convert(Ptr{Nothing}, pointer(A))
getstate(heap::CompactHeap) = (heap.allocated, heap.len)

getfreelist(heap::CompactHeap, i) = FreelistPtr(pointer(heap.freelists, i))
freelists(heap::CompactHeap) = FreelistIterator(heap.freelists)

function CompactHeap(
    allocator::AbstractAllocator,
    sz;
    pool = Local,
    # Power of two
    minallocation = PowerOfTwo(21),
)
    minallocation = poweroftwo(minallocation)

    # Round down size to the nearest multiple of `minallocation`
    sz = (sz >> minallocation.val) << minallocation.val

    # Allocate the memory managed by this heap
    base, token = allocate(allocator, sz)
    allocated = 0

    numbins = getbin(minallocation, sz)
    status = FindNextTree(numbins)
    freelists = [Block() for _ in Base.OneTo(numbins)]

    heap = CompactHeap(
        minallocation,
        base,
        sz,
        allocated,
        pool,
        status,
        freelists,
        true,
        token,
    )

    finalizer(heap) do x
        free(x.base, x.token)
    end

    # Add an entry for the biggest freelist.
    block = Block(base)
    block.size = sz
    block.backsize = 0
    block.free = true
    block.pool = pool
    block.next = Block()
    block.previous = Block()

    push!(heap, block)
    return heap
end

getbin(heap::CompactHeap, sz) = getbin(heap.minallocation, sz)
getbin(heap::CompactHeap, block::Block) = getbin(heap, block.size)

binsize(heap::CompactHeap, bin) = binsize(heap.minallocation, bin)

basepointer(heap::CompactHeap) = heap.base
Base.sizeof(heap::CompactHeap) = heap.len
baseaddress(heap::CompactHeap) = convert(UInt, heap.base)
numbins(heap::CompactHeap) = length(heap.freelists)

# The CompactHeap maintains backsizes, so we can also walk backwards.
function walkprevious(::CompactHeap, block)
    # Invariant, the first block will always have a backsize of zero.
    backsize = block.backsize
    if iszero(backsize)
        return nothing
    else
        return Block(pointer(block) - backsize)
    end
end

"""
    sortfreelists!(heap)

Sort all of the freelists in `heap` in order of increasing address.
"""
function sortfreelists!(heap)
    for freelist in freelists(heap)
        sort!(freelist)
    end
end

function Base.push!(heap::CompactHeap, block::Block)
    bin = getbin(heap, block)
    list = @inbounds(getfreelist(heap, bin))

    isempty(list) && @inbounds(setentry!(heap.status, bin))
    push!(list, block)
    return nothing
end

popfreelist!(heap::CompactHeap, bin::Integer) = popfreelist!(heap, getfreelist(heap, bin), bin)
Base.@propagate_inbounds function popfreelist!(
    heap::CompactHeap,
    list::AbstractFreelist,
    bin::Integer,
)
    block = pop!(list)
    isempty(list) && clearentry!(heap.status, bin)
    return block
end

# Pop a block for this bin.
# If a freelist doesn't exist for this bin, split the next largest block.
#
# Assume `sz` is a multiple of the heap granularity
function Base.pop!(heap::CompactHeap, bin)
    if isset(heap.status, bin)
        return @inbounds(popfreelist!(heap, bin))
    end
    # list = @inbounds getfreelist(heap, bin)
    # if !isempty(list)
    #     return @inbounds(popfreelist!(heap, list, bin))
    # end

    # Okay, we don't have a block.
    # We need to get the next highest available block and split it.
    nextbin = @inbounds(findnext(heap.status, bin))

    # No bigger block available.
    # This means we can't fulfill this request.
    # Return `nothing` to indicate this.
    nextbin === nothing && return nothing

    # Otherwise, take this block and split it.
    # TODO: Maybe set a limit where we don't split if the over-allocation amount is
    # less than this limit.
    block = @inbounds(popfreelist!(heap, nextbin))

    # Now that we have the block, split it down to the size we want and put the remainder
    # back in the heap.
    sz = binsize(heap, bin)
    remainder = Block(pointer(block) + sz)

    remainder.free = true
    remainder.size = block.size - sz
    remainder.next = Block()
    remainder.previous = Block()

    # Since we split this block, we need to maintain the backsize pointers for efficient
    # merging.
    block.size = sz
    remainder.backsize = block.size
    next = walknext(heap, remainder)
    if next !== nothing
        next.backsize = remainder.size
    end

    push!(heap, remainder)
    return block
end

function canpop(heap::CompactHeap, bin)
    isset(heap.status, bin) && return true

    # Okay, we don't have a block.
    # We need to get the next highest available block and split it.
    nextbin = @inbounds(findnext(heap.status, bin))

    # No bigger block available - cannot allocate.
    return nextbin === nothing ? false : true
end

function remove!(heap::CompactHeap, block::Block)
    bin = getbin(heap, block)
    list = getfreelist(heap, bin)
    remove!(list, block)
    isempty(list) && @inbounds(clearentry!(heap.status, bin))
    return nothing
end

function putback!(heap::CompactHeap, block::Block)
    @check isfree(block)

    # Check the block before and after are free.
    #
    # If so, merge the blocks
    previous = walkprevious(heap, block)
    next = walknext(heap, block)

    if previous !== nothing && isfree(previous)
        remove!(heap, previous)
        previous.size += block.size

        # Purge metadata for this block to try and catch errors earlier.
        zerometa!(block)
        block = previous
    end

    if next !== nothing
        # Similar to the case with "block.previous", if the next block is being evicted,
        # we don't want to try to merge since that will break eviction logic.
        if isfree(next)
            remove!(heap, next)
            block.size += next.size

            # Set the backedge for the nextnext block
            nextnext = walknext(heap, next)
            if nextnext !== nothing
                nextnext.backsize = block.size
            end

            # Purge metadata for the next block to try and catch errors earlier.
            zerometa!(next)
        else
            next.backsize = block.size
        end
    end

    push!(heap, block)
    return nothing
end

#####
##### High level API
#####

const CAN_ALLOC_MESSAGE = """
Something went wrong. Trying to allocate from a heap that has been temporarily marked
as "no-allocate". It's likely something went wrong during eviction.
"""

function alloc(heap::CompactHeap, bytes::Integer, id::UInt)
    heap.canalloc || error(CAN_ALLOC_MESSAGE)
    iszero(bytes) && return nothing

    needed_size = max(bytes + headersize(), one(bytes) << heap.minallocation.val)
    needed_size > sizeof(heap) && return nothing
    bin = getbin(heap, needed_size)
    block = pop!(heap, bin)
    block === nothing && return nothing

    # Setup the default state for the block
    heap.allocated += sizeof(block)
    clearbits!(block)
    block.id = id
    block.pool = heap.pool

    # Preemptively mark as dirty.
    # Up to callee's to determine if this is actually clean or not.
    block.dirty = true
    ptr = pointer(block) + headersize()
    return ptr
end

@inline function canalloc(heap::CompactHeap, bytes::Integer)
    iszero(bytes) && return false

    needed_size = max(bytes + headersize(), one(bytes) << heap.minallocation.val)
    needed_size > sizeof(heap) && return false
    bin = getbin(heap, needed_size)
    return canpop(heap, bin)
end

function free(heap::CompactHeap, ptr::Ptr{Nothing})
    block = unsafe_block(ptr)
    if block.evicting == true
        block.orphaned = true
        return nothing
    end
    block.queued = false
    block.free = true
    heap.allocated -= sizeof(block)
    putback!(heap, block)
    return nothing
end

#####
##### Defrag API
#####

function nextfree(heap, block = baseblock(heap))
    while block !== nothing && !isfree(block)
        block = walknext(heap, block)
    end
    return block
end

"""
    defrag!(f::F, heap::CompactHeap; [nthreads], [queued_callback])

Defragment `heap`. Callback function `f` is given an ID and a new block.
If a block is marked as `queued`, then, by default, we can't touch it.
However, the caller may provide a zero argument `queued_callback` which must free this
block, in which case defragmentation may continue.
"""
function defrag!(
    f::F,
    heap::CompactHeap;
    nthreads = Threads.nthreads(),
    queued_callback = nothing,
) where {F}
    @checknothing freeblock = nextfree(heap)
    @checknothing block = walknext(heap, freeblock)

    while true
        # By the invariants we keep on the heap, the next block should be not-free.
        # However, it MIGHT be queued for deletion, in which case, we can't move it.
        if block.queued
            if queued_callback !== nothing
                queued_callback()
                # Recurse to begin the walk again since it's possible that blocks
                # we've already visited have now been cleaned up.
                # println("REWALKING!")
                defrag!(f, heap; nthreads, queued_callback)
                return nothing
            else
                @checknothing freeblock = nextfree(heap, block)
                @checknothing block = walknext(heap, freeblock)
            end
            continue
        end

        # The following operations need to be atomic as far as the GC is concerned.
        # This takes care of if the GC runs while we're moving data and the block we
        # are moving suddenly becomes free ...

        # Now, we need to move "block" to "freeblock" and then fixup the heap.
        id = block.id
        # safeprint("Defrag: Moving $id")
        pool = block.pool
        block.backsize = freeblock.backsize
        freeblock_size = freeblock.size
        remove!(heap, freeblock)

        # After the "memcpy", the new location for "block" is "freeblock".
        # Thus, we need to perform the callback so that references can be updated.
        aliases = pointer(freeblock) + sizeof(block) > pointer(block)
        f(id, freeblock, block)
        if aliases
            _memcpy!(pointer(freeblock), pointer(block), sizeof(block); nthreads = nothing)
        else
            _memcpy!(pointer(freeblock), pointer(block), sizeof(block); nthreads)
        end

        @check freeblock.id == id

        # Update siblings
        sibling = getsibling(freeblock)
        if sibling !== nothing
            sibling.sibling = freeblock
        end

        # Since we haven't changed the size of the block, "walknext" should still work,
        # though the actual block returned will be garbage.
        block = freeblock
        freeblock = walknext(heap, block)

        freeblock.backsize = block.size
        freeblock.size = freeblock_size
        freeblock.free = true

        # TODO: In all likelyhood, we're just going to be removing the freeblock again.
        # However, dealing with queued blocks makes this tricky so do the conservative
        # thing for now.
        putback!(heap, freeblock)

        @checknothing block = walknext(heap, freeblock)
        block.backsize = freeblock.size
    end
    return nothing
end

#####
##### Eviction API
#####

function unsafe_free!(f::F, block::Block) where {F}
    @check !isfree(block)

    # Mark this block as being evicted.
    # This will short-circuit any calls to `free` the result from the callback while
    # still keeping the block as `not-free` to be compatible with the other heap logic.
    block.evicting = true

    # Do the callback on the block and then mark the blocked as freed.
    return f(block)
end

function evict!(f::F, heap::CompactHeap, block::Block) where {F}
    if isfree(block)
        remove!(heap, block)
        block.free = false
        heap.allocated += sizeof(block)
        return nothing
    else
        return unsafe_free!(f, block)
    end
end

# Helper method to get a block from a pointer.
function evictfrom!(heap::CompactHeap, ptr::Ptr{Nothing}, sz; kw...)
    evictfrom!(heap, unsafe_block(ptr), sz; kw...)
end

function evictfrom!(heap::CompactHeap, block::Block, sz; cb = donothing)
    heap.canalloc = false
    sizefreed = 0

    last = block
    current = block
    aborted = false

    # NOTE ON ABORTING.
    # Control flow must always pass through the bottom of this function to ensure
    # that even if an eviction is aborted early that we maintain the heap in a
    # consistent state.

    # Walk forward until either we have freed enough space or reach the end of the heap.
    while true
        abort = evict!(cb, heap, current)

        # Check for early abort
        if (abort === true)
            # Cleanup eviction status
            current.evicting = false

            # Aborted on the first block.
            # If this block was NOT orphaned as a side-effect of the eviction callback,
            # then we need to rollback the current block to the last successfully
            # processed block.
            #
            # However, if this block was recaimed, then it is free but hasn't been
            # put back on the heap.
            #
            # As a result, we need to include it in the freed frontier.
            if current.orphaned
                sizefreed += current.size
            else
                if iszero(sizefreed)
                    heap.canalloc = true
                    return nothing
                end
                current = last
            end

            @goto cleanup
        end
        sizefreed += current.size

        # Update state
        last = current
        current = walknext(heap, current)
        # If we've either freed enough memory or walked off the end of the array, exit.
        (sizefreed >= sz || current === nothing) && break
    end

    # Walk backwards if we have to in order to free enough space.
    while sizefreed < sz
        candidate = walkprevious(heap, block)
        candidate === nothing && break
        # Intenionally overwrite the argument "block".
        block = candidate
        abort = evict!(cb, heap, block)
        if (abort === true)
            block.evicting = false
            if !block.orphaned
                # Roll back block to last successful eviction.
                block = walknext(heap, block)
            else
                sizefreed += block.size
            end
            @goto cleanup
        end
        sizefreed += block.size
    end

    # Resize our block and put it back in the heap.
    @label cleanup
    block.free = true
    block.orphaned = false
    block.size = sizefreed
    heap.allocated -= sizefreed

    putback!(heap, block)
    heap.canalloc = true
    return nothing
end

#####
##### Materializer
#####

# Touch all the pages in an allocation.
function materialize_os_pages!(heap::CompactHeap)
    # First, we must make sure that NOTHING is allocated in the heap - otherwise, we might
    # over write some existing data.
    if length(heap) != 1 || !isfree(first(heap))
        throw(error("Can only materialize OS pages for an empty heap!"))
    end

    # Now, we threaded touch all pages in the heap
    #
    # Take steps of 4096, which is the smallest possible page size.
    ptr = convert(Ptr{UInt8}, datapointer(baseblock(heap)))
    Polyester.@batch per = core for i = 1:4096:sizeof(heap)
        unsafe_store!(ptr, one(UInt8), i)
    end
    return nothing
end

function touchheap(heap::CompactHeap)
    v = zeros(UInt8, Threads.nthreads())
    ptr = convert(Ptr{UInt8}, datapointer(baseblock(heap)))
    Polyester.@batch per = core for i = 1:64:sizeof(heap)
        v[Threads.threadid()] += unsafe_load(ptr, i)
    end
    return sum(v)
end

#####
##### Checker
#####

# Check the various invariants of the CompactHeap
#
# 1. First Block should always have zero backsize
# 2. Any freelist entry in the freelist dict must have an entry in the `status` tree.
# 3. All non-entries in the `status` tree must not have entry in the freelist dict.
# 4. Forward and backsizes must be consistent.
# 5. Free blocks found by walking the heap must be consistent with free-blocks found walking
#    the freelists

function first_block_invariant(heap::CompactHeap)
    block = baseblock(heap)
    if !iszero(block.backsize)
        println("Base block has nonzero backsize: $(block.backsize)")
        return false
    end
    return true
end

function freelist_status_invariant(heap::CompactHeap)
    maxbin = getbin(heap, sizeof(heap))
    passed = true

    # Walk the `status` tree and correlate the entries with the freelist dict.
    for (i, mask) in enumerate(last(heap.status.runs))
        for j = 1:64
            bin = 64 * (i - 1) + j
            bin > maxbin && return passed

            if hasentryat(mask, j)
                # There must be a key in the `freelists` dictionary and it must be non-empty
                if isempty(getfreelist(heap, bin))
                    println(
                        "Heap has a `status` entry for bin $bin but no corresponding freelist bound",
                    )
                    passed = false
                end
            else
                # No entry here - there must not be an entry in the dict.
                if !isempty(getfreelist(heap, bin))
                    println("Heap as no `status` entry for $bin but found a freelist")
                    passed = false
                end
            end
        end
    end
    return passed
end

function size_invariant(heap::CompactHeap)
    passed = true

    # This certainly better not be nothing!
    # Heaps should always have at least one entry
    y = iterate(heap)
    if y === nothing
        println("Empty Heap!!")
        passed = false
    end
    block, state = y

    count = 1
    while true
        y = iterate(heap, state)
        y === nothing && break
        next, state = y

        if next.backsize != block.size
            println(
                "Block $(count + 1) has backsize $(next.backsize). Expected $(block.size)",
            )
            passed = false
        end

        if block.free && next.free
            println("Unmerged free blocks in position $(count) and $(count + 1)")
            passed = false
        end
        count += 1
        block = next
    end
    return passed
end

function free_invariant(heap::CompactHeap)
    passed = true
    counts = Dict{Int,Int}()

    for block in heap
        if isfree(block)
            bin = getbin(heap, block)
            get!(counts, bin, 0)
            counts[bin] += 1
        end
    end

    for (bin, freelist) in enumerate(freelists(heap))
        if !isempty(freelist)
            if !haskey(counts, bin)
                println("Found a freelist for bin $bin but no block of this size found")
                passed = false
            else
                if length(freelist) != counts[bin]
                    println(
                        "Length of freelist is $(length(freelist)) but counted $(counts[bin]) blocks",
                    )
                    passed = false
                end
                delete!(counts, bin)
            end
        end
    end

    # Did we miss any
    if !isempty(counts)
        println("Unaccounted blocks")
        passed = false
        for (bin, count) in counts
            println("    Bin: $bin - Count: $count")
        end
    end

    return passed
end

function check(heap::CompactHeap)
    passed = true
    if !first_block_invariant(heap)
        printstyled(stdout, "Heap failed First Block Invariant\n"; color = :red)
        passed = false
    end

    if !freelist_status_invariant(heap)
        printstyled(stdout, "Heap failed First List Invariant\n"; color = :red)
        passed = false
    end

    if !size_invariant(heap)
        printstyled(stdout, "Heap failed size invariant\n"; color = :red)
        passed = false
    end

    if !free_invariant(heap)
        printstyled(stdout, "Heap failed free invariant\n"; color = :red)
        passed = false
    end
    return passed
end


text:
#!/usr/bin/env julia

struct Literal
    name
    polarity
end
function Literal(v)
    if v < 0
        Literal(-v, false)
    else
        Literal(v, true)
    end
end

struct Cnf
    clauses
end

function firstclause(dimacs)
    line = readline(dimacs)
    while startswith(line, 'c') || startswith(line, 'p')
        line = readline(dimacs)
    end

    return line
end

function Cnf(f::String)
    clauses = []
    dimacs = open(f, "r")
    line = firstclause(dimacs)
    push!(clauses, Set())
    while line != ""
        literals = [parse(Int, l) for l in split(line)]
        i = 1
        while (i = findnext(l -> l == 0, literals, i)) != nothing
            union!(clauses[1], map(Literal, literals[1:i-1]))
            pushfirst!(clauses, Set())
            literals = literals[i+1:length(literals)]
        end
        union!(clauses[1], map(Literal, literals))

        line = readline(dimacs)
    end

    @assert isempty(clauses[1]) "dimacs not terminated with 0"
    popfirst!(clauses)

    Cnf(reverse(clauses))
end

function check(assignments, cnf)
    assigned = Set()
    # Note that unsat is a set instead of an array here. This will
    # work since if two clauses are identical, they will both be
    # satisfied by the same conditions.
    unsat = Set(cnf.clauses)

    for v in assignments
        l = Literal(v)
        @assert !(l.name in assigned) "literal is already assigned"
        union!(assigned, [l.name])

        for clause in unsat
            if l in clause
                setdiff!(unsat, [clause])
            end
        end
    end

    @assert isempty(unsat) "not all clauses are satisfied"
end

assignments = [parse(Int, line) for line in readlines(stdin)]
cnf = Cnf(ARGS[1])
check(assignments, cnf)

ext: jl

max_stars_repo_name: velentr/sat

max_issues_repo_name: velentr/sat

max_forks_repo_name: velentr/sat


text:
# This file provides reexported functions.

using ArnoldiMethod

"""
    adjacency_matrix(g[, T=Int; dir=:out])

Return a sparse adjacency matrix for a graph, indexed by `[u, v]`
vertices. Non-zero values indicate an edge between `u` and `v`. Users may
override the default data type (`Int`) and specify an optional direction.

### Optional Arguments
`dir=:out`: `:in`, `:out`, or `:both` are currently supported.

### Implementation Notes
This function is optimized for speed and directly manipulates CSC sparse matrix fields.
"""
function adjacency_matrix(g::AbstractGraph, T::DataType=Int; dir::Symbol=:out)
    nzmult = 1
    # see below - we iterate over columns. That's why we take the
    # "opposite" neighbor function. It's faster than taking the transpose
    # at the end.
    if (dir == :out)
        _adjacency_matrix(g, T, inneighbors, 1)
    elseif (dir == :in)
        _adjacency_matrix(g, T, outneighbors, 1)
    elseif (dir == :both)
        if is_directed(g)
            _adjacency_matrix(g, T, all_neighbors, 2)
        else
            _adjacency_matrix(g, T, outneighbors, 1)
        end
    else
        error("Not implemented")
    end
end

function _adjacency_matrix(g::AbstractGraph{U}, T::DataType, neighborfn::Function, nzmult::Int=1) where U
    n_v = nv(g)
    nz = ne(g) * (is_directed(g) ? 1 : 2) * nzmult
    colpt = ones(U, n_v + 1)

    rowval = sizehint!(Vector{U}(), nz)
    selfloops = Vector{U}()
    for j in 1:n_v  # this is by column, not by row.
        if has_edge(g, j, j)
            push!(selfloops, j)
        end
        dsts = neighborfn(g, j)
        colpt[j + 1] = colpt[j] + length(dsts)
        append!(rowval, sort!(dsts))
    end
    spmx = SparseMatrixCSC(n_v, n_v, colpt, rowval, ones(T, nz))

    # this is inefficient. There should be a better way of doing this.
    # the issue is that adjacency matrix entries for self-loops are 2,
    # not one(T).
    if !is_directed(g)
        for i in selfloops
            if !(T <: Bool)
                spmx[i, i] += one(T)
            end
        end
    end
    return spmx
end

"""
    laplacian_matrix(g[, T=Int; dir=:unspec])

Return a sparse [Laplacian matrix](https://en.wikipedia.org/wiki/Laplacian_matrix)
for a graph `g`, indexed by `[u, v]` vertices. `T` defaults to `Int` for both graph types.

### Optional Arguments
`dir=:unspec`: `:unspec`, `:both`, :in`, and `:out` are currently supported.
For undirected graphs, `dir` defaults to `:out`; for directed graphs,
`dir` defaults to `:both`.
"""
function laplacian_matrix(g::AbstractGraph{U}, T::DataType=Int; dir::Symbol=:unspec) where U
    if dir == :unspec
        dir = is_directed(g) ? :both : :out
    end
    A = adjacency_matrix(g, T; dir=dir)
    D = convert(SparseMatrixCSC{T,U}, Diagonal(sparse(sum(A, dims=2)[:])))
    return D - A
end

"""
    laplacian_spectrum(g[, T=Int; dir=:unspec])

Return the eigenvalues of the Laplacian matrix for a graph `g`, indexed
by vertex. Default values for `T` are the same as those in
[`laplacian_matrix`](@ref).

### Optional Arguments
`dir=:unspec`: Options for `dir` are the same as those in [`laplacian_matrix`](@ref).

### Performance
Converts the matrix to dense with ``nv^2`` memory usage.

### Implementation Notes
Use `eigs(laplacian_matrix(g);  kwargs...)` to compute some of the
eigenvalues/eigenvectors.
"""
laplacian_spectrum(g::AbstractGraph, T::DataType=Int; dir::Symbol=:unspec) = eigvals(Matrix(laplacian_matrix(g, T; dir=dir)))

"""
    adjacency_spectrum(g[, T=Int; dir=:unspec])

Return the eigenvalues of the adjacency matrix for a graph `g`, indexed
by vertex. Default values for `T` are the same as those in
[`adjacency_matrix`](@ref).

### Optional Arguments
`dir=:unspec`: Options for `dir` are the same as those in [`laplacian_matrix`](@ref).

### Performance
Converts the matrix to dense with ``nv^2`` memory usage.

### Implementation Notes
Use `eigs(adjacency_matrix(g);  kwargs...)` to compute some of the
eigenvalues/eigenvectors.
"""
function adjacency_spectrum(g::AbstractGraph, T::DataType=Int; dir::Symbol=:unspec)
    if dir == :unspec
        dir = is_directed(g) ?  :both : :out
    end
    return eigvals(Matrix(adjacency_matrix(g, T; dir=dir)))
end

"""
    incidence_matrix(g[, T=Int; oriented=false])

Return a sparse node-arc incidence matrix for a graph, indexed by
`[v, i]`, where `i` is in `1:ne(g)`, indexing an edge `e`. For
directed graphs, a value of `-1` indicates that `src(e) == v`, while a
value of `1` indicates that `dst(e) == v`. Otherwise, the value is
`0`. For undirected graphs, both entries are `1` by default (this behavior
can be overridden by the `oriented` optional argument).

If `oriented` (default false) is true, for an undirected graph `g`, the
matrix will contain arbitrary non-zero values representing connectivity
between `v` and `i`.
"""
function incidence_matrix(g::AbstractGraph, T::DataType=Int; oriented=false)
    isdir = is_directed(g)
    n_v = nv(g)
    n_e = ne(g)
    nz = 2 * n_e

    # every col has the same 2 entries
    colpt = collect(1:2:(nz + 1))
    nzval = repeat([(isdir || oriented) ? -one(T) : one(T), one(T)], n_e)

    # iterate over edges for row indices
    rowval = zeros(Int, nz)
    i = 1
    for u in vertices(g)
        for v in outneighbors(g, u)
            if isdir || u < v # add every edge only once
                rowval[2 * i - 1] = u
                rowval[2 * i] = v
                i += 1
            end
        end
    end

    spmx = SparseMatrixCSC(n_v, n_e, colpt, rowval, nzval)
    return spmx
end

"""
    spectral_distance(G₁, G₂ [, k])

Compute the spectral distance between undirected n-vertex
graphs `G₁` and `G₂` using the top `k` greatest eigenvalues.
If `k` is ommitted, uses full spectrum.

### References
- JOVANOVIC, I.; STANIC, Z., 2014. Spectral Distances of Graphs Based on their Different Matrix Representations
"""
function spectral_distance end

# can't use Traitor syntax here (https://github.com/mauro3/SimpleTraits.jl/issues/36)
@traitfn function spectral_distance(G₁::G, G₂::G, k::Integer) where {G <: AbstractGraph; !IsDirected{G}}
    A₁ = adjacency_matrix(G₁)
    A₂ = adjacency_matrix(G₂)

    λ₁ = k < nv(G₁) - 1 ? eigs(A₁, nev=k, which=LR())[1] : eigvals(Matrix(A₁))[end:-1:(end - (k - 1))]
    λ₂ = k < nv(G₂) - 1 ? eigs(A₂, nev=k, which=LR())[1] : eigvals(Matrix(A₂))[end:-1:(end - (k - 1))]

    return sum(abs, (λ₁ - λ₂))
end

# can't use Traitor syntax here (https://github.com/mauro3/SimpleTraits.jl/issues/36)
@traitfn function spectral_distance(G₁::G, G₂::G) where {G <: AbstractGraph; !IsDirected{G}}
    nv(G₁) == nv(G₂) || throw(ArgumentError("Spectral distance not defined for |G₁| != |G₂|"))
    return spectral_distance(G₁, G₂, nv(G₁))
end


text:
using Oceananigans.Solvers: solve_for_pressure!, solve_poisson_equation!

function pressure_solver_instantiates(FT, Nx, Ny, Nz, planner_flag)
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), extent=(100, 100, 100))
    solver = PressureSolver(CPU(), grid, PressureBoundaryConditions(grid), planner_flag)
    return true  # Just making sure the PressureSolver does not error/crash.
end

function divergence_free_poisson_solution(arch, FT, topology, Nx, Ny, Nz, planner_flag=FFTW.MEASURE)
    ArrayType = array_type(arch)
    grid = RegularCartesianGrid(FT, topology=topology, size=(Nx, Ny, Nz), extent=(1.0, 2.5, π))
    fbcs = TracerBoundaryConditions(grid)
    pbcs = PressureBoundaryConditions(grid)
    solver = PressureSolver(arch, grid, fbcs, planner_flag)

    # Generate right hand side from a random (divergent) velocity field.
    Ru = CellField(FT, arch, grid, UVelocityBoundaryConditions(grid))
    Rv = CellField(FT, arch, grid, VVelocityBoundaryConditions(grid))
    Rw = CellField(FT, arch, grid, WVelocityBoundaryConditions(grid))
    U = (u=Ru, v=Rv, w=Rw)

    set!(Ru, rand(Nx, Ny, Nz))
    set!(Rv, rand(Nx, Ny, Nz))
    set!(Rw, rand(Nx, Ny, Nz))

    # Adding (nothing, nothing) in case we need to dispatch on ::NFBC
    fill_halo_regions!(Ru, arch, nothing, nothing)
    fill_halo_regions!(Rv, arch, nothing, nothing)
    fill_halo_regions!(Rw, arch, nothing, nothing)

    # Compute the right hand side R = ∇⋅U
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, U.u.data, U.v.data, U.w.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    ϕ   = CellField(FT, arch, grid, pbcs)  # "pressure"
    ∇²ϕ = CellField(FT, arch, grid, pbcs)

    # Using Δt = 1 but it doesn't matter since velocities = 0.
    solve_for_pressure!(ϕ.data, solver, arch, grid, 1, datatuple(U))

    fill_halo_regions!(ϕ, arch)
    event = launch!(arch, grid, :xyz, ∇²!, grid, ϕ.data, ∇²ϕ.data, dependencies=Event(device(arch)))
    wait(device(arch), event)
    fill_halo_regions!(∇²ϕ, arch)

    return CUDA.@allowscalar interior(∇²ϕ) ≈ R
end

#####
##### Test that Poisson solver error converges as error ~ N⁻²
#####

ψ(::Type{Bounded}, n, x) = cos(n*x/2)
ψ(::Type{Periodic}, n, x) = cos(n*x)

k²(::Type{Bounded}, n) = (n/2)^2
k²(::Type{Periodic}, n) = n^2

using Oceananigans.Solvers: permute_index, unpermute_index

@kernel function permute_indices!(dst, src, solver_type, arch, grid)
    i, j, k = @index(Global, NTuple)

    i′, j′, k′ = permute_index(solver_type, arch, i, j, k, grid.Nx, grid.Ny, grid.Nz)

    @inbounds dst[i′, j′, k′] = src[i, j, k]
end

@kernel function unpermute_indices!(dst, src, solver_type, arch, grid)
    i, j, k = @index(Global, NTuple)

    i′, j′, k′ = unpermute_index(solver_type, arch, i, j, k, grid.Nx, grid.Ny, grid.Nz)

    @inbounds dst[i′, j′, k′] = src[i, j, k]
end

function analytical_poisson_solver_test(arch, N, topo; FT=Float64, mode=1)
    grid = RegularCartesianGrid(FT, topology=topo, size=(N, N, N), x=(0, 2π), y=(0, 2π), z=(0, 2π))
    solver = PressureSolver(arch, grid, TracerBoundaryConditions(grid))

    xC, yC, zC = nodes((Cell, Cell, Cell), grid, reshape=true)

    Tx, Ty, Tz = topology(grid)
    Ψ(x, y, z) = ψ(Tx, mode, x) * ψ(Ty, mode, y) * ψ(Tz, mode, z)
    f(x, y, z) = -(k²(Tx, mode) + k²(Ty, mode) + k²(Tz, mode)) * Ψ(x, y, z)

    if arch isa GPU && topo == PBB_topo
        storage = solver.storage.storage1
    else
        storage = solver.storage
    end

    buffer = similar(storage)
    buffer .= convert(array_type(arch), f.(xC, yC, zC))

    event = launch!(arch, grid, :xyz,
                    permute_indices!, storage, buffer, solver.type, arch, grid,
                    dependencies = Event(device(arch)))
    wait(device(arch), event)

    solve_poisson_equation!(solver, grid)

    event = launch!(arch, grid, :xyz,
                    unpermute_indices!, buffer, storage, solver.type, arch, grid,
                    dependencies = Event(device(arch)))
    wait(device(arch), event)

    ϕ = real(Array(buffer))

    L¹_error = mean(abs, ϕ - Ψ.(xC, yC, zC))

    return L¹_error
end

function poisson_solver_convergence(arch, topo, N¹, N²; FT=Float64)
    error¹ = analytical_poisson_solver_test(arch, N¹, topo; FT=FT)
    error² = analytical_poisson_solver_test(arch, N², topo; FT=FT)

    rate = log(error¹ / error²) / log(N² / N¹)

    Tx, Ty, Tz = topo
    @info "Convergence of L¹-normed error, $(typeof(arch)), $FT, ($(N¹)³ -> $(N²)³), topology=($Tx, $Ty, $Tz): $rate"

    return isapprox(rate, 2, rtol=5e-3)
end

#####
##### Run pressure solver tests
#####

const PPP_topo = (Periodic, Periodic, Periodic)
const PPB_topo = (Periodic, Periodic, Bounded)
const PBB_topo = (Periodic, Bounded,  Bounded)
const BBB_topo = (Bounded,  Bounded,  Bounded)

topos = (PPP_topo, PPB_topo, PBB_topo, BBB_topo)

@testset "Pressure solvers" begin
    @info "Testing pressure solvers..."

    @testset "Pressure solver instantiation" begin
        @info "  Testing pressure solver instantiation..."

        for FT in float_types
            @test pressure_solver_instantiates(FT, 32, 32, 32, FFTW.ESTIMATE)
            @test pressure_solver_instantiates(FT, 1,  32, 32, FFTW.MEASURE)
            @test pressure_solver_instantiates(FT, 32,  1, 32, FFTW.ESTIMATE)
            @test pressure_solver_instantiates(FT,  1,  1, 32, FFTW.MEASURE)
        end
    end

    @testset "Divergence-free solution [CPU]" begin
        @info "  Testing divergence-free solution [CPU]..."

        for topo in topos
            @info "    Testing $topo topology on square grids..."
            for N in [7, 16], FT in float_types
                @test divergence_free_poisson_solution(CPU(), FT, topo, N, N, N, FFTW.ESTIMATE)
                @test divergence_free_poisson_solution(CPU(), FT, topo, 1, N, N, FFTW.MEASURE)
                @test divergence_free_poisson_solution(CPU(), FT, topo, N, 1, N, FFTW.ESTIMATE)
                @test divergence_free_poisson_solution(CPU(), FT, topo, 1, 1, N, FFTW.MEASURE)
            end
        end

        Ns = [11, 16]
        for topo in topos
            @info "    Testing $topo topology on rectangular grids..."
            for Nx in Ns, Ny in Ns, Nz in Ns, FT in float_types
                @test divergence_free_poisson_solution(CPU(), FT, topo, Nx, Ny, Nz, FFTW.ESTIMATE)
            end
        end
    end

    @hascuda @testset "Divergence-free solution [GPU]" begin
        @info "  Testing divergence-free solution [GPU]..."
        for topo in (PPP_topo, PPB_topo, PBB_topo)
            @info "    Testing $topo topology on GPUs..."
            @test divergence_free_poisson_solution(GPU(), Float64, topo, 16, 16, 16)
            @test divergence_free_poisson_solution(GPU(), Float64, topo, 32, 32, 16)
            @test divergence_free_poisson_solution(GPU(), Float64, topo, 16, 32, 24)
	        @test divergence_free_poisson_solution(GPU(), Float64, topo, 5,  7,  11)
        end
    end

    @testset "Convergence to analytical solution [CPU]" begin
        @info "  Testing convergence to analytical solution [CPU]..."
        for topo in topos
            @test poisson_solver_convergence(CPU(), topo, 2^6, 2^7)
            @test poisson_solver_convergence(CPU(), topo, 67, 131)
        end
    end

    @hascuda @testset "Convergence to analytical solution [GPU]" begin
        @info "  Testing convergence to analytical solution [GPU]..."
        for topo in (PPP_topo, PPB_topo, PBB_topo)
            @test poisson_solver_convergence(GPU(), topo, 2^6, 2^7)
            @test poisson_solver_convergence(GPU(), topo, 67, 131)
        end
    end
end


text:
# Functions for contracting two arrays, both for general StridedArrays as well as
# for SymArrays. we only have a few specialized functions that we have needed up
# to now, but carefully optimize each of them

using TensorOperations
using LinearAlgebra

# Array[i]*SymArray[(i,j,k)]
# indices 1, 2, and 3 are exchangeable here
function contract(A::StridedVector{T},S::SymArray{(3,),U},n::Union{Val{1},Val{2},Val{3}}) where {T,U}
    TU = promote_type(T,U)
    @assert size(S,1) == length(A)
    res = SymArray{(2,),TU}(size(S,1),size(S,2))
    contract!(res,A,S,n)
end

# We know that $j\leq k$ (because $R$ is itself exchange symmetric)
# \begin{align}
# R_{jk} &= \sum_{i=1}^N g_i S_{ijk}
# \end{align}

# Matrix elements represented by $S_{ijk}$:
# \begin{equation}
# \begin{cases}
# S_{ijk}, S_{ikj}, S_{jik}, S_{jki}, S_{kij}, S_{kji} & i<j<k\\
# S_{ijk}, S_{jik}, S_{jki} & i<j=k\\
# S_{ijk}, S_{ikj}, S_{kij} & i=j<k\\
# S_{ijk} & i=j=k
# \end{cases}
# \end{equation}
# We only need to take the contributions that show up for each $R_{jk}$

# Array[i]*SymArray[(i,j,k)]
# indices 1, 2, and 3 are exchangeable here
function contract!(res::SymArray{(2,),TU}, A::StridedVector{T}, S::SymArray{(3,),U}, n::Union{Val{1},Val{2},Val{3}}) where {T,U,TU}
    # only loop over S once, and put all the values where they should go
    # R[j,k] = sum_i A[i] B[i,j,k]
    # S[i,j,k] with i<=j<=k represents the 6 (not always distinct) terms: Bijk, Bikj, Bjik, Bjki, Bkij, Bkji
    # since R[j,k] is also exchange symmetric, we only need to calculate j<=k
    # this means we only have to check each adjacent pair of indices and include
    # their permutation if not equal, but keeping the result indices ordered
    @assert size(S,1) == length(A)
    @assert size(S,1) == size(res,1)
    res.data .= 0
    @inbounds for (v,inds) in zip(S.data,CartesianIndices(S))
        i,j,k = Tuple(inds)
        res[j,k] += v*A[i]
        # have to include i<->j
        if i<j res[i,k] += v*A[j] end
        # have to include i<->k, but keep indices of res sorted
        if j<k res[i,j] += v*A[k] end
    end
    res
end

# Array[k]*SymArray[(i,j),k]
function contract(A::StridedVector{T},S::SymArray{(2,1),U},n::Val{3}) where {T,U}
    TU = promote_type(T,U)
    sumsize = length(A)
    @assert sumsize == size(S,3)
    res = SymArray{(2,),TU}(size(S,1),size(S,2))
    contract!(res,A,S,n)
end

# Array[k]*SymArray[(i,j),k]
function contract!(res::SymArray{(2,),TU},A::StridedVector{T},S::SymArray{(2,1),U},::Val{3}) where {T,U,TU}
    # use that S[(i,j),k] == S[I,k] (i.e., the two symmetric indices act like a "big" index)
    mul!(res.data,reshape(S.data,:,length(A)),A)
    res
end

# Array[i]*SymArray[(i,j),k)]
# since indices 1 and 2 are exchangeable here, use this
function contract(A::StridedVector{T},S::SymArray{(2,1),U},n::Union{Val{1},Val{2}}) where {T,U}
    TU = promote_type(T,U)
    @assert size(S,1) == length(A)
    # the result is a normal 2D array
    res = Array{TU,2}(undef,size(S,2),size(S,3))
    contract!(res,A,S,n)
end

# Array[i]*SymArray[(i,j),k]
# since indices 1 and 2 are exchangeable here, use this
function contract!(res::StridedArray{TU,2},A::StridedVector{T},S::SymArray{(2,1),U},::Union{Val{1},Val{2}}) where {T,U,TU}
    # only loop over S once, and put all the values where they should go
    @assert size(A,1) == size(S,1)
    @assert size(res,1) == size(S,1)
    @assert size(res,2) == size(S,3)
    res .= zero(TU)
    @inbounds for (v,inds) in zip(S.data,CartesianIndices(S))
        i1,i2,i3 = Tuple(inds)
        res[i2,i3] += v*A[i1]
        # if i1 != i2, we have to add the equal contribution from S[i2,i1,i3]
        if i1 != i2
            res[i1,i3] += v*A[i2]
        end
    end
    res
end

# Array[i]*SymArray[(i,j)]
# this is symmetric in i1 and i2
function contract(A::StridedVector{T},S::SymArray{(2,),U},n::Union{Val{1},Val{2}}) where {T,U}
    TU = promote_type(T,U)
    @assert size(S,1) == length(A)
    # the result is a normal 1D vector
    res = Vector{TU}(undef,size(S,2))
    contract!(res,A,S,n)
end

# Array[i]*SymArray[(i,j)]
# this is symmetric in i1 and i2
function contract!(res::StridedVector{TU},A::StridedVector{T},S::SymArray{(2,),U},::Union{Val{1},Val{2}}) where {T,U,TU}
    @assert size(A,1) == size(S,1)
    @assert size(res,1) == size(S,1)
    res .= zero(TU)
    # only loop over S once, and put all the values where they should go
    for (v,inds) in zip(S.data,CartesianIndices(S))
        i1,i2 = Tuple(inds)
        res[i2] += v*A[i1]
        # if i1 != i2, we have to add the contribution from S[i2,i1]
        if i1 != i2
            res[i1] += v*A[i2]
        end
    end
    res
end

# Array[i_n]*Array[i1,i2,i3,...,iN]
function contract(A::StridedVector{T},B::StridedArray{U,N},::Val{n}) where {T,U,N,n}
    TU = promote_type(T,U)
    @assert 1 <= n <= N

    resdims = size(B)[1:N .!= n]
    res = similar(B,TU,resdims)

    A = convert(AbstractArray{TU},A)
    B = convert(AbstractArray{TU},B)
    contract!(res,A,B,Val{n}())
end

mygemv!(args...) = BLAS.gemv!(args...)

function _contract_middle!(res,A,B)
    @inbounds for k=1:size(B,3)
        mul!(@view(res[:,k]), @view(B[:,:,k]), A)
    end
end

# Array[i_n]*Array[i1,i2,i3,...,iN]
function contract!(res::StridedArray{TU},A::StridedVector{TU},B::StridedArray{TU,N},::Val{n}) where {TU,N,n}
    nsum = length(A)
    @assert size(B,n) == nsum
    @assert ndims(res)+1 == ndims(B)
    ii = 0
    for jj = 1:ndims(B)
        jj==n && continue
        ii += 1
        @assert size(B,jj) == size(res,ii)
    end

    if n==1      # A[i]*B[i,...]
        mygemv!('T',one(TU),reshape(B,nsum,:),A,zero(TU),vec(res))
    elseif n==N  # B[...,i]*A[i]
        mygemv!('N',one(TU),reshape(B,:,nsum),A,zero(TU),vec(res))
    else
        rightsize = prod(size(B,i) for i=n+1:N)
        Br = reshape(B,:,nsum,rightsize)
        resr = reshape(res,:,rightsize)
        _contract_middle!(resr,A,Br)
    end
    res
end

# Array[i_n]*Array[i1,i2,i3,...,iN]
function contract!(res::StridedArray{TU,Nres},A::StridedArray{TU,NA},B::StridedArray{TU,NB},::Val{nA},::Val{nB}) where {TU,Nres,NA,NB,nA,nB}
    # use specialized code for 1D-arrays if possible
    if NA==1
        @assert nA==1
        return contract!(res,A,B,Val(nB))
    elseif NB==1
        @assert nB==1
        return contract!(res,B,A,Val(nA))
    end

    nsum = size(A,nA)
    @assert size(B,nB) == nsum
    @assert Nres == NA + NB - 2
    ii = 0
    for jj = 1:NA
        jj==nA && continue
        ii += 1
        @assert size(A,jj) == size(res,ii)
    end
    for jj = 1:NB
        jj==nB && continue
        ii += 1
        @assert size(B,jj) == size(res,ii)
    end

    if nA==NA && nB==1      # A[...,i]*B[i,...]
        resAsize = prod(ntuple(i->size(A,i),Val(NA-1)))
        mul!(reshape(res,resAsize,:),reshape(A,resAsize,nsum),reshape(B,nsum,:))
    else
        packedsize(M,::Val{n}) where n = begin
            nM1 = prod(ntuple(i->size(M,i),Val(n-1)))
            nM2 = size(M,n)
            nM3 = prod(ntuple(i->size(M,n+i),Val(ndims(M)-n)))
            (nM1, nM2, nM3)
        end
        # just use @tensor after reshaping to pack together
        NAps = packedsize(A,Val(nA))
        NBps = packedsize(B,Val(nB))
        Nresps = (NAps[1],NAps[3],NBps[1],NBps[3])
        Ap = reshape(A,NAps...)
        Bp = reshape(B,NBps...)
        resp = reshape(res,Nresps...)
        @tensor resp[iA1,iA3,iB1,iB3] = Ap[iA1,ii,iA3] * Bp[iB1,ii,iB3]
    end
    res
end

"""return the symmetry group index and the number of symmetric indices in the group"""
@inline which_symgrp(S::T,nS) where T<:SymArray = which_symgrp(T,nS)
@inline @generated function which_symgrp(::Type{<:SymArray{Nsyms}},nS) where Nsyms
    grps = ()
    for (ii,Nsym) in enumerate(Nsyms)
        grps = (grps...,ntuple(_->ii,Nsym)...)
    end

    quote
        ng = $grps[nS]
        ng, $Nsyms[ng]
    end
end

"""Check if the arguments correspond to a valid contraction. Do all "static" checks at compile time."""
@generated function check_contraction_compatibility(res::SymArray{Nsymsres,TU}, A::StridedArray{T,NA}, S::SymArray{NsymsS,U}, ::Val{nA}, ::Val{nS}) where {T,U,TU,NsymsS,Nsymsres,NA,nA,nS}
    promote_type(T,U) <: TU || error("element types not compatible: T = $T, U = $U, TU = $TU")

    contracted_group, Nsym_ctrgrp = which_symgrp(S,nS)
    NsymsA_contracted = ntuple(_->1,NA-1)
    if Nsym_ctrgrp == 1
        NsymsS_contracted = TupleTools.deleteat(NsymsS,contracted_group)
    else
        NsymsS_contracted = Base.setindex(NsymsS,Nsym_ctrgrp-1,contracted_group)
    end
    Nsymsres_check = (NsymsA_contracted...,NsymsS_contracted...)
    # assure that symmetry structure is compatible
    errmsg = "symmetry structure not compatible:"
    errmsg *= "\nNA = $NA, NsymsS = $NsymsS, nA = $nA, nS = $nS"
    errmsg *= "\nNsymsres = $Nsymsres, expected Nsymssres = $Nsymsres_check"
    Nsymsres == Nsymsres_check || error(errmsg)

    Sresinds = ((1:nS-1)...,(nS+1:sum(NsymsS))...)
    Aresinds = ((1:nA-1)...,(nA+1:NA)...)
    sizerescheck = ((i->:(sizeA[$i])).(Aresinds)...,
                    (i->:(sizeS[$i])).(Sresinds)...)
    code = quote
        sizeA = size(A)
        sizeS = size(S)
        @assert sizeA[$nA] == sizeS[$nS]
        @assert size(res) == ($(sizerescheck...),)
    end
    #display(code)
    code
end

"""
A[iAprev,icntrct,iApost]
S[iSprev,Icntrct,ISpost]
res[iAprev,iApost,iSprev,Icntrct-1,ISpost]
"""
@generated function contract_symindex!(res::Array{TU,5}, A::Array{T,3}, ::Val{sizeA13unit}, S::Array{U,3}, ::Val{sizeS13unit}, ::Val{Nsym}) where {T,U,TU,sizeA13unit,sizeS13unit,Nsym}
    # Nsym-dimensional index tuple
    iS2s = Tuple(Symbol.(:iS2_,1:Nsym))
    # combined (Nsym-1)-dimensional index for the Nsym possible permutations
    iSm2s = Tuple(Symbol.(:iSm2_,1:Nsym))

    iAsetters = Expr(:block)
    iAusers = Expr(:block)
    for n = 1:Nsym
        chk = n==1 ? true : :( $(iS2s[n-1])<$(iS2s[n]) )
        syminds = TupleTools.deleteat(iS2s,n)
        push!(iAsetters.args, :( $chk && ($(iSm2s[n]) = symgrp_sortedsub2ind($(syminds...))) ))
        push!(iAusers.args, :( $chk && (res[iA1,iA3,iS1,$(iSm2s[n]),iS3] += v * A[iA1,$(iS2s[n]),iA3]) ))
    end

    iA1max = sizeA13unit[1] ? 1 : :(size(A,1))
    iA3max = sizeA13unit[2] ? 1 : :(size(A,3))
    iS1max = sizeS13unit[1] ? 1 : :(size(S,1))
    iS3max = sizeS13unit[2] ? 1 : :(size(S,3))

    code = quote
        # size of iterated index is middle index of A
        iterdimS = SymIndexIter($Nsym,size(A,2))

        res .= zero(TU)
        @inbounds for iS3 = 1:$iS3max
            for (iS2,IS) = enumerate(iterdimS)
                ($(iS2s...),) = Tuple(IS)
                $iAsetters
                for iS1 = 1:$iS1max
                    v = S[iS1,iS2,iS3]
                    for iA3 = 1:$iA3max
                        for iA1 = 1:$iA1max
                            $iAusers
                        end
                    end
                end
            end
        end
    end
    code
end

function contract!(res::StridedArray{TU,Nres}, A::StridedArray{T,NA}, S::SymArray{NsymsS,U}, ::Val{nA}, ::Val{nS}) where {T,U,TU,NsymsS,Nres,NA,nA,nS}
    # StridedArray is equivalent to SymArray with all Nsyms equal to 1, provide this as overlay of that data
    ressym = SymArray{ntuple(i->1,Val(Nres))}(res,size(res)...)
    contract!(ressym,A,S,Val(nA),Val(nS))
    res
end

@generated function contract!(res::SymArray{Nsymsres,TU}, A::StridedArray{T,NA}, S::SymArray{NsymsS,U}, ::Val{nA}, ::Val{nS}) where {T,U,TU,NsymsS,Nsymsres,NA,nA,nS}
    contracted_group, Nsym_ctrgrp = which_symgrp(S,nS)
    sizeA13unit = (nA==1,nA==NA)
    sizeS13unit = (contracted_group==1,contracted_group==length(NsymsS))
    newsize_centered_expr(sizeA::Symbol,nA::Int,NA::Int) = begin
        sizeAs = [:( $sizeA[$ii] ) for ii=1:NA]
        t1 = nA > 1 ? :( *($(sizeAs[1:nA-1]...)) ) : 1
        t3 = nA < NA ? :( *($(sizeAs[nA+1:NA]...)) ) : 1
        :( ($t1, $sizeA[$nA], $t3) )
    end

    code = quote
        # first check that all the sizes are compatible etc
        check_contraction_compatibility(res,A,S,Val($nA),Val($nS))

        sizeA = size(A)
        sizeAp = $(newsize_centered_expr(:sizeA,nA,NA))
        Apacked = reshape(A,sizeAp)

        grpsizeS = symgrp_size.(S.Nts,Val.(NsymsS))
        sizeSp = $(newsize_centered_expr(:grpsizeS, contracted_group, length(NsymsS)))
        Spacked = reshape(S.data,sizeSp)

        if $Nsym_ctrgrp > 1
            size_respS1 = symgrp_size(S.Nts[$contracted_group],Val($(Nsym_ctrgrp-1)))
            respacked = reshape(res.data,sizeAp[1],sizeAp[3],sizeSp[1],size_respS1,sizeSp[3])
            contract_symindex!(respacked,Apacked,Val($sizeA13unit),Spacked,Val($sizeS13unit),Val($Nsym_ctrgrp))
        else
            # iS2 is a single (not symmetric) dimension and it disappears after summing
            # we have res[iS1,iS3,iA1,iA3] = S_iS1,k,iS3 * A_iA1,k,iA3
            respacked = reshape(res.data,sizeAp[1],sizeAp[3],sizeSp[1],sizeSp[3])
            @tensor respacked[iA1,iA3,iS1,iS3] = Apacked[iA1,ii,iA3] * Spacked[iS1,ii,iS3]
        end
    end
    #display(code)
    code
end


text:
using AdversarialDriving
using AutomotiveSimulator
using POMDPs, POMDPPolicies, POMDPSimulators
using AutomotiveVisualization

sut_agent = BlinkerVehicleAgent(up_left(id=1, s=25., v=15.), TIDM(Tint_TIDM_template, noisy_observations = true))
adv1 = BlinkerVehicleAgent(left_straight(id=2, s=0., v=15.0), TIDM(Tint_TIDM_template))
adv2 = BlinkerVehicleAgent(left_turnright(id=3, s=10., v=15.0), TIDM(Tint_TIDM_template))
adv3 = BlinkerVehicleAgent(right_straight(id=4, s=40., v=20.0), TIDM(Tint_TIDM_template))
adv4 = BlinkerVehicleAgent(right_turnleft(id=5, s=30., v=20.0), TIDM(Tint_TIDM_template))

mdp = AdversarialDrivingMDP(sut_agent, [adv1, adv2, adv3, adv4], Tint_roadway, 0.1)

s = rand(initialstate(mdp))
cam = StaticCamera(position=(-10,-10), zoom=10.)
c = render([mdp.roadway, s,
        FancyCar(car = s[4], color = colorant"blue")
                ], camera = cam, surface=AutomotiveVisualization.CairoPDFSurface(IOBuffer(), AutomotiveVisualization.canvas_width(cam), AutomotiveVisualization.canvas_height(cam)))
write("T_intersection.pdf", c)

null_action = actions(mdp)[1]
blinker_action = actions(mdp)[7]

# Nominal Behavior
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> null_action))
scenes_to_gif(state_hist(hist), mdp.roadway, "Tint_nominal.gif")

# Behavior with noise
blinker_pol = FunctionPolicy((s) -> !blinker(get_by_id(s, 2)) ? blinker_action : null_action)
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, blinker_pol)
scenes_to_gif(state_hist(hist), mdp.roadway, "Tint_failure.gif")

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
using Test
using DataFrames
using CSV
using ChemistryFeaturization.Utils.ElementFeatureUtils

@testset "ElementFeatureUtils" begin
    df = CSV.read(abspath(@__DIR__, "..", "test_data", "lookup_table.csv"), DataFrame)

    # fea_minmax
    @test_throws AssertionError fea_minmax("heffalump")
    @test fea_minmax("Group") == [1, 18]
    @test fea_minmax("MeaningOfLife", df) == [-1, 42]

    # default_log
    @test default_log("Block") == false # not numbers
    @test default_log("MeaningOfLife", df) == false # values span 0
    @test default_log("Valence") == false # extremal value 0
    @test default_log("Atomic mass") == true
    @test default_log("Atomic mass", threshold = 3) == false

    # default_categorical
    @test default_categorical("Block") == true
    @test get_param_vec(false, 2) == [false, false]
    @test get_param_vec([false, false], 2) == [false, false]
    @test get_param_vec([false, false], 3) == [false, false, false]
    @test get_param_vec([false, false, false], 2) == [false, false]
    @test get_param_vec([3, 2], 3, pad_val = 0) == [3, 2, 0]
    @test get_param_vec(['a'], 3, pad_val = 'b') == ['a', 'b', 'b']
    @test get_param_vec('a', 3, pad_val = 'b') == ['a', 'a', 'a']

    @test get_param_vec([false, false], 3) == [false, false, false]
    @test get_param_vec([false, false, false], 2) == [false, false]
    @test get_param_vec([3, 2], 3, pad_val = 0) == [3, 2, 0]
    @test get_param_vec(['a'], 3, pad_val = 'b') == ['a', 'b', 'b']
    @test get_param_vec('a', 3, pad_val = 'b') == ['a', 'a', 'a']

    # get_bins
    @test get_bins("Block") == ["s", "p", "d", "f"]
    @test get_bins("first_letter", df) == ["A", "C", "T"]
    @test get_bins("Group", categorical = false, nbins = 2) == 1.0:8.5:18.0
    @test get_bins("MeaningOfLife", df, categorical = true) == [-1, 0, 42]
    logbins = get_bins("Atomic mass")
    @test isapprox(logbins[2] / logbins[1], logbins[end] / logbins[end-1])
    neglogbins = get_bins("neg_nums", df, logspaced = true)
    @test neglogbins[1] == -1000
    # errors and things that should be ignored...
    @test get_bins("Block", nbins = 3) == ["s", "p", "d", "f"]
    @test get_bins("Valence", logspaced = true) == 0:1:13
    @test_throws AssertionError get_bins("Valence", categorical = false, logspaced = true)
    @test get_bins("noarsenic", df) == 1.0:0.1:2.0
    @test get_bins("noarsenic", df, categorical = true) == [1, 2]

end


text:
using OrdinaryDiffEq, Test

f(u,p,t) = u
prob = ODEProblem(f,1.0,(0.0,100eps(1.0)))
dt = 4.7eps(1.0)
cond(u,p,t) = u - exp(70eps(1.0))

c = Ref(0)
function affect!(integrator)
	c[] += 1
end
cb = ContinuousCallback(cond, affect!)
@info "Event Repeat Test 1"
sol = solve(prob, Tsit5(); adaptive=false,callback=cb,dt=dt,save_everystep=false)
@test c[] == 1

function condition_v(out, u, t, integrator)
	out[1] = u - exp(50eps(1.0))
	out[2] = u - exp(60eps(1.0))
	out[3] = u - exp(70eps(1.0))
	out[4] = u - exp(80eps(1.0))
	out[5] = u - exp(90eps(1.0))
end

c1 = Ref(0)
c2 = Ref(0)
c3 = Ref(0)
c4 = Ref(0)
c5 = Ref(0)

function affect_v!(integrator, idx)
	if idx == 1
		c1[] += 1
	elseif idx == 2
		c2[] += 1
	elseif idx == 3
		c3[] += 1
	elseif idx == 4
		c4[] += 1
	elseif idx == 5
		c5[] += 1
	end
end

cb = VectorContinuousCallback(condition_v, affect_v!, 5)
@info "Event Repeat Test 2"
sol = solve(prob, Tsit5(); adaptive=false,callback=cb,dt=dt,save_everystep=false)

@test c1[] == 1
@test c2[] == 1
@test c3[] == 1
@test c4[] == 1
@test c5[] == 1


text:
using GraphletCounting
using Test

include("DegreeOrientedDAG.jl")
include("CountTriangles.jl")

println("add some utils to run orca and retrieve certain counts for certain nodes")


text:
@testset "chainrules" begin
    Random.seed!(1)

    @testset "general: single input" begin
        for x in (1.0, -1.0, 0.0, 0.5, 10.0, -17.1, 1.5 + 0.7im)
            test_scalar(erf, x)
            test_scalar(erfc, x)
            test_scalar(erfi, x)

            test_scalar(airyai, x)
            test_scalar(airyaiprime, x)
            test_scalar(airybi, x)
            test_scalar(airybiprime, x)

            test_scalar(erfcx, x)
            test_scalar(dawson, x)

            if x isa Real
                test_scalar(invdigamma, x)
            end

            if x isa Real && 0 < x < 1
                test_scalar(erfinv, x)
                test_scalar(erfcinv, x)
            end

            if x isa Real && x > 0 || x isa Complex
                test_scalar(gamma, x)
                test_scalar(digamma, x)
                test_scalar(trigamma, x)
            end
        end
    end

    @testset "Bessel functions" begin
        for x in (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
            test_scalar(besselj0, x)
            test_scalar(besselj1, x)

            isreal(x) && x < 0 && continue

            test_scalar(bessely0, x)
            test_scalar(bessely1, x)

            for nu in (-1.5, 2.2, 4.0)
                test_frule(besseli, nu, x)
                test_rrule(besseli, nu, x)

                test_frule(besselj, nu, x)
                test_rrule(besselj, nu, x)

                test_frule(besselk, nu, x)
                test_rrule(besselk, nu, x)

                test_frule(bessely, nu, x)
                test_rrule(bessely, nu, x)

                # use complex numbers in `rrule` for FiniteDifferences
                test_frule(hankelh1, nu, x)
                test_rrule(hankelh1, nu, complex(x))

                # use complex numbers in `rrule` for FiniteDifferences
                test_frule(hankelh2, nu, x)
                test_rrule(hankelh2, nu, complex(x))
            end
        end
    end

    @testset "beta and logbeta" begin
        test_points = (1.5, 2.5, 10.5, 1.6 + 1.6im, 1.6 - 1.6im, 4.6 + 1.6im)
        for _x in test_points, _y in test_points
            # ensure all complex if any complex for FiniteDifferences
            x, y = promote(_x, _y)
            test_frule(beta, x, y)
            test_rrule(beta, x, y)

            test_frule(logbeta, x, y)
            test_rrule(logbeta, x, y)
        end
    end

    @testset "log gamma and co" begin
        # It is important that we have negative numbers with both odd and even integer parts
        test_points = (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
        for x in test_points
            for m in (0, 1, 2, 3)
                test_frule(polygamma, m, x; rtol=1e-8)
                test_rrule(polygamma, m, x; rtol=1e-8)
            end

            isreal(x) && x < 0 && continue
            test_scalar(loggamma, x)
            for a in test_points
                # ensure all complex if any complex for FiniteDifferences
                _a, _x = promote(a, x)
                test_frule(gamma, _a, _x; rtol=1e-8)
                test_rrule(gamma, _a, _x; rtol=1e-8)

                test_frule(loggamma, _a, _x)
                test_rrule(loggamma, _a, _x)
            end

            isreal(x) || continue
            test_frule(logabsgamma, x)
            test_rrule(logabsgamma, x; output_tangent=(randn(), randn()))
            for a in test_points
                isreal(a) && a > 0 || continue
                test_frule(gamma_inc, a, x, 0)
                test_rrule(gamma_inc, a, x, 0; output_tangent=(randn(), randn()))
            end
        end
    end
end


text:
"""
SelectTransform
```
    dims = [1,3,5,6,7]
    tr = SelectTransform(dims)
    X = rand(100,10)
    transform(tr,X,obsdim=2) == X[dims,:]
```
Select the dimensions `dims` that the kernel is applied to.
"""
struct SelectTransform{T<:AbstractVector{<:Int}} <: Transform
    select::T
    function SelectTransform{V}(dims::V) where {V}
        new{V}(dims)
    end
end

function SelectTransform(dims::V) where {V<:AbstractVector{T} where  {T<:Int}}
    @assert all(dims.>0) "Selective dimensions should all be positive integers"
    SelectTransform{V}(dims)
end

set!(t::SelectTransform{<:AbstractVector{T}},dims::AbstractVector{T}) where {T<:Int} = t.select .= dims

params(t::SelectTransform) = t.select

duplicate(t::SelectTransform,θ) = t

Base.maximum(t::SelectTransform) = maximum(t.select)

function transform(t::SelectTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs)
    @boundscheck maximum(t) >= size(X,feature_dim(obsdim)) ?
        throw(DimensionMismatch("The highest index $(maximum(t)) is higher then the feature dimension of X : $(size(X,feature_dim(obsdim)))")) : nothing
    @inbounds _transform(t,X,obsdim)
end

function transform(t::SelectTransform,x::AbstractVector{<:Real},obsdim::Int=defaultobs) #TODO Add test
    @assert maximum(t) <= length(x) "The highest index $(maximum(t)) is higher then the vector length : $(length(x))"
    return @inbounds view(x,t.select)
end

_transform(t::SelectTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs) = obsdim == 2 ? view(X,t.select,:) : view(X,:,t.select)


text:
# =========================================================================================
"""
    mutable struct XieModel{T}
        radius ::T                                 # radius of the origin-centered sphere
        charges::Vector{Charge{T}}                 # point charges in the sphere
        params ::Option{T}         = defaultopt(T) # system constants
    end

System model of a dielectric sphere containing multiple point charges. On construction, the
given point charge model will be translated and rescaled to fit inside an origin-centered
sphere with the specified radius.

# Special contructors
```julia
XieModel{T}(
    radius ::T,
    charges::Vector{Charge{T}},
    params ::Option{T}          = defaultopt(T);
    # kwargs
    compat ::Bool               = false
)
```
`compat` enables the compatibility mode and scales the model exactly like the reference
implementation ([[Xie16]](@ref Bibliography)). Use this flag if you intend to compare
the results to the reference.
"""
mutable struct XieModel{T}
    "radius of the origin-centered sphere"
    radius::T
    "point charges in the sphere"
    charges::Vector{Charge{T}}
    "system constants"
    params::Option{T}

    XieModel{T}(
        radius ::T,
        charges::Vector{Charge{T}},
        params ::Option{T} = defaultopt(T);
        compat ::Bool = false
    ) where T = new(radius, scalemodel(charges, radius, compat=compat), params)
end

XieModel(
    radius ::T,
    charges::Vector{Charge{T}},
    params ::Option{T} = defaultopt(T);
    compat ::Bool = false
) where T = XieModel{T}(radius, charges, params, compat=compat)

# =========================================================================================
"""
    function scalemodel{T}(
        charges::Vector{Charge{T}},
        radius ::T;
        # kwargs
        compat ::Bool              = false
    )

Translates and rescales the given point charge model to fit inside an origin-centered
sphere with the specified radius. More specifically, the function will center the given
model at the origin and scaled in a way such that the outermost point charge will be
located at 80% `radius` distance from the origin.

# Arguments
 * `compat` Enables compatibility mode and scales the model exactly like the reference
   implementation ([[Xie16]](@ref Bibliography)). Use this flag if you intend to compare
   the results to the reference.

# Return type
`Vector{Charge{T}}`
"""
function scalemodel(
        charges::Vector{Charge{T}},
        radius ::T;
        compat ::Bool               = false
    ) where T
    # center model
    cpos   = sum.([extrema(q.pos[i] for q in charges) for i in 1:3])/2
    newpos = [q.pos - cpos for q in charges]

    # compute and apply scaling factor
    h = x -> compat ? ceil(x) : x  # in compat mode, scaling factor is rounded up
    sf = .8radius / h(maximum(norm(pos) for pos in newpos))
    rmul!(newpos, sf)

    [Charge{T}(pos, q.val) for (pos, q) in zip(newpos, charges)]
end


text:
module FireBall

using ..Ahorn, Maple

function fireballFinalizer(entity::Maple.FireBall)
    x, y = Ahorn.position(entity)

    entity.data["nodes"] = [(x + 16, y)]
end

const placements = Ahorn.PlacementDict(
    "Fireball" => Ahorn.EntityPlacement(
        Maple.FireBall,
        "point",
        Dict{String, Any}(
            "amount" => 3
        ),
        fireballFinalizer
    )
)

Ahorn.nodeLimits(entity::Maple.FireBall) = 1, -1

sprite = "objects/fireball/fireball01.png"

function Ahorn.selection(entity::Maple.FireBall)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.FireBall)
    px, py = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)

        px, py = nx, ny
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.FireBall, room::Maple.Room)
    x, y = Ahorn.position(entity)
    Ahorn.drawSprite(ctx, sprite, x, y)
end

end


text:
using SnoopCompile

println("tests infer benchmark")

@snoopiBench "Juno"


text:
function readModelAndGenerateMeshMref(m,pad::Int64,newSize::Vector,domain::Vector)
########################## m,mref are in Velocity here. ###################################

mref = copy(m);
# mref[:,:,1:end-17] = getSimilarLinearModel(m[:,:,1:end-17],1.75,2.8);
mref = getSimilarLinearModel(m,1.75,2.9);

sea = abs.(m[:] .- minimum(m)) .< 5e-2;
mref[sea] = m[sea];
if newSize!=[]
	m    = expandModelNearest(m,   collect(size(m)),newSize);
	mref = expandModelNearest(mref,collect(size(mref)),newSize);
end

Minv = getRegularMesh(domain,collect(size(m))-1);

(mPadded,MinvPadded) = addAbsorbingLayer(m,Minv,pad);
(mrefPadded,MinvPadded) = addAbsorbingLayer(mref,Minv,pad);

N = prod(MinvPadded.n+1);
boundsLow  = minimum(mPadded);
boundsHigh = maximum(mPadded);

boundsLow  = ones(N)*boundsLow;
boundsLow = convert(Array{Float32},boundsLow);
boundsHigh = ones(N)*boundsHigh;
boundsHigh = convert(Array{Float32},boundsHigh);

return (mPadded,MinvPadded,mrefPadded,boundsHigh,boundsLow);
end


text:
function fg_one!(w::Vector, g::Vector, i::Int64, X1::Matrix, y::Vector, C::Real, n::Int64, n_c::Int64, idmi::Vector)
  m = length(g)
  halfnorm = 0.5 * dot(w, w)

  psis = psi_list(w, X1, y, i, n_c, idmi)
  psis_id, xi = best_psis(psis)      # most violated constraints

  dxi = calc_dconst((i,psis_id), X1, y, n_c, idmi)

  for i=1:m
    g[i] = w[i] +  C * dxi[i]
  end

  return halfnorm + C * xi
end

# train adversarial method using SGD
function train_adv_sgd(X::Matrix, y::Vector, C::Real=1.0;
  step::Float64=0.1, use_adagrad::Bool=false,
  ftol::Real=1e-6, grtol::Real=1e-6,
  show_trace::Bool=true, max_iter::Int=1000)

  n = length(y)
  # add one
  X1 = [ones(n) X]'   # transpose
  m = size(X1, 1)

  # number of class
  n_c = maximum(y)
  n_f = n_c * m   # number of features

  # parameters. init with zero
  w = rand(n_f) - 0.5
  grad = zeros(n_f)

  # prepare saved vars
  idmi = map(i -> idi(m, i), collect(1:n_c))

  # adagrad
  # rate is something you need to set beforehand
  rate = step
  square_g = zeros(n_f)   # for storing historical square of grads

  f_prev = Inf
  iter = 0

  while true
    iter = iter + 1

    ids = randperm(n)
    sum_f = 0.0

    for i in ids

      f = fg_one!(w, grad, i, X1, y, C, n, n_c, idmi)
      sum_f += f
      if show_trace
        println("iter : ", iter, ", sample : ", i, ", C : ", C, ", f : ", f, ", abs grad : ", mean(abs(grad)))
      end

      if use_adagrad
        # adagrad
        square_g += grad .^ 2
        w = w - (rate ./ sqrt(square_g)) .* grad
      else
        w = w - step * grad
      end

    end

    # discount step
    step = step * 0.95

    if iter >= max_iter
      if show_trace println("maximum iteration reached!!") end
      break
    end

    if mean(abs(grad)) < grtol
      if show_trace println("gradient breaks!!") end
      break
    end

    f = sum_f / n
    if abs(f_prev - f) < ftol
      if show_trace println("function breaks!!") end
      break
    end
    f_prev = f

  end

  # finalizing losses
  gv_aug = zeros(n)
  gv_01 = zeros(n)
  l_adv = zeros(n)
  l_01 = zeros(n)
  for i=1:n
    psis = psi_list(w, X1, y, i, n_c, idmi)
    psis_id, val = best_psis(psis)      # most violated constraints
    n_ps = length(psis_id)

    gv_aug[i] = val

    # compute probs
    p_hat = zeros(n_c)
    p_check = zeros(n_c)
    for j=1:n_c
      if j in psis_id
        p_hat[j] = ( (n_ps-1.0)*psis[j] - sum(psis[psis_id[psis_id .!= j]]) + 1.0 ) / n_ps
        p_check[j] = 1.0 / n_ps
      else
        p_hat[j] = 0.0
        p_check[j] = 0.0
      end
    end

    C01 = 1 - eye(n_c)    # 01 loss matrix
    v = p_hat' * C01 * p_check   # the result is vector size 1 not a number
    gv_01[i] = v[1]     #

    ## training loss
    l_adv[i] = 1.0 - p_hat[y[i]]
    l_01[i] = 1.0 - round(Int, indmax(p_hat) == y[i])
  end

  game_value_01 = mean(gv_01)
  game_value_augmented = mean(gv_aug)

  # create model
  adv_model = MultiAdversarialModel(w, zeros(0), Tuple{Integer, Vector}[], n_c, game_value_01, game_value_augmented, mean(l_adv), mean(l_01))

  return adv_model::MultiAdversarialModel
end

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
# slt_trim.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 0.7
# 09.12.2016
# Last Edit: 18.04.2019

"""# slt_trim

Remove columns from SLT files, source files must not have all equal numbers of columns

`slt_trim(source,destination,mindate,maxdate,maxanalogcols)`\n
* **source**::String = Source directory, `\"K:\\Data\"`
* **destination**::String = Destination directory, `\"K:\\Data\"`
* **mindate**::DateTime = Start of period to process
* **maxdate**::DateTime = End of period to process
* **maxanalogcols**::Int = Maximum number of columns that should remain"""
function slt_trim(f1::String,f2::String,mindate::DateTime,maxdate::DateTime,maxanalogcols::Int)
	# f1::String # Source Directory
	# f2::String # Destination Directory
	# mindate::DateTime
	# maxdate::DateTime
	# maxanalogcols::Int # Maximum number of analog columns

	# Check Inputs
	if !isdir(f1)
		error("First input must be a directory")
	end
	if !isdir(f2)
		error("Second input must be a directory")
	end

	# List and sort files
	(files,folders) = dirlist(f1,regex=r"\w\d{8}\.")

	# Parse File Dates
	fdates = Array{DateTime}(undef,length(files))
	i = 0
	for i=1:1:length(files)
		try
			yearstr = Meta.parse(files[i][end-14:end-11])
			daystr = Meta.parse(files[i][end-10:end-8])
			hourstr = Meta.parse(files[i][end-7:end-6])
			minutestr = Meta.parse(files[i][end-5:end-4])
			fdates[i] = DateTime(yearstr) + Dates.Day(daystr) - Dates.Day(1) + Dates.Hour(hourstr) + Dates.Minute(minutestr)
		catch
			println(files[i])
			error("Failed to parse the header of " * files[i])
		end
	end

	# Sort the Files By Date
	f = sortperm(fdates)
	fdates = fdates[f]
	files = files[f]

	# Save Config files before dates are filtered
	f = [occursin(r"\.cfg$",i) for i in files]
	cfgfiles = files[f]
	cfgdates = fdates[f]

	# Remove files beyond the given time bounds
	f = findall(mindate .<= fdates)
	fdates = fdates[f]
	files = files[f]
	f = findall(maxdate .> fdates)
	fdates = fdates[f]
	files = files[f]

	# Display Info
	println("Adjusting SLT Analog Columns")
	println("   Source:      " * f1)
	println("   Destination: " * f2)
	println("   Number of files: " * string(length(files)))

	# Sort Files by Type
	f = [occursin(r"\.slt$",i) for i in files]
	sltfiles = files[f]
	sltdates = fdates[f]
	f = [occursin(r"\.csr$",i) for i in files]
	csrfiles = files[f]
	csrdates = fdates[f]
	f = [occursin(r"\.csu$",i) for i in files]
	csufiles = files[f]
	csudates = fdates[f]
	f = [occursin(r"\.flx$",i) for i in files]
	flxfiles = files[f]
	flxdates = fdates[f]
	f = [occursin(r"\.log$",i) for i in files]
	logfiles = files[f]
	logdates = fdates[f]
	f = [occursin(r"\.csv$",i) for i in files]
	csvfiles = files[f]
	csvdates = fdates[f]

	# Load the CFG info
	configs = slt_config(cfgfiles)
	configs = [configs;deepcopy(configs[1])] # Setup a fake final row with a time well beyond any real file name
	configs[end]["Time"] = DateTime(9999)

	# Load the SLT header data
	println("\nLoading SLT Header Data")
	sltinfo = DataFrame()
	offset = 1 # File processing list offset
	while offset < length(sltfiles)
		f = findall(get.(configs,"Time",Array{DateTime}(undef,1)) .<= sltdates[offset])[end] # Find latest config
		config1 = configs[f,:] # Current config and lower bound on file list
		config2 = configs[f+1,:] # Upper bound on file list

		f = findall(get.(config1,"Time",Array{DateTime}(undef,1)) .<= sltdates .< get.(config2,"Time",Array{DateTime}(undef,1)))
		tempfiles = sltfiles[f]
		analog_count = get(config1[1],"Analog_Count",0) # Length?
		sample_frequency = Meta.parse(get(config1[1],"Frequency",0))
		for i=1:1:length(tempfiles)
			temp = slt_header(tempfiles[i],analog_count,sample_frequency)
			temp[!,:Analog_Count] .= analog_count
			temp[!,:Sample_Frequency] .= sample_frequency
			if isempty(sltinfo)
				sltinfo = temp
			else
				sltinfo = [sltinfo;temp]
			end
			offset += 1
		end
	end

	# Process the SLT Files
	println("\nProcessing SLT Files")
	for i = 1:1:size(sltinfo,1)
		temp = splitdir(sltinfo.FileName[i])[2]
		new_filename = joinpath(f2,temp)

		# Column Difference
		colldiff = length(sltinfo.Channels[i]) - maxanalogcols
		if colldiff <= 0
			println("   ## " * temp * " ##   - Too few columns, skipped")
			continue
		end
		println("   " * temp * ", " * string(colldiff) * " columns removed")

		fid = open(new_filename,"w+") # Open New File
		# Write the header
		write(fid,Int8(2(maxanalogcols + 4))) # Bytes per record = Number of columns x 2
		write(fid,Int8(sltinfo.EddyMeas_Version[i]))
		write(fid,Int8(Dates.value(Dates.Day(sltinfo.T0[i]))))
		write(fid,Int8(Dates.value(Dates.Month(sltinfo.T0[i]))))
		yr1 = Int8(floor(Float64(Dates.value(Dates.Year(sltinfo.T0[i])))/100))
		yr2 = Int8(Int(Dates.value(Dates.Year(sltinfo.T0[i]))) - 100*floor(Float64(Dates.value(Dates.Year(sltinfo.T0[i])))/100))
		write(fid,yr1)
		write(fid,yr2)
		write(fid,Int8(Dates.value(Dates.Hour(sltinfo.T0[i]))))
		write(fid,Int8(Dates.value(Dates.Minute(sltinfo.T0[i]))))

		# Write Channels and Bit Masks
		for j=1:1:sltinfo.Analog_Count[i] - colldiff
			write(fid,Int8(sltinfo.Bit_Mask[i][j])) # Bit Mask
			write(fid,Int8(sltinfo.Channels[i][j])) # Channel
		end

		# Read and Write Data
		fid0 = open(sltinfo.FileName[i]) # Open Original File
		seek(fid0,sltinfo.Start_Pos[i]) # Move file position to data
		early_end = false
		for j=1:1:sltinfo.Line_Count[i]
			for k=1:1:4 + length(sltinfo.Channels[i]) - colldiff
				if eof(fid0)
					early_end = true
				else
					temp_data = read!(fid0,Array{Int16}(undef,1))

					# Only save columns within column limit
					write(fid,temp_data)
				end
			end
			for k = 1:1:colldiff
				if eof(fid0)
					early_end = true
				else
					temp_data = read!(fid0,Array{Int16}(undef,1))
				end
			end
		end
		if early_end
			println("    Early End of File: " * string(sltinfo.Line_Count[i]) * " lines") # Temp
		end
		close(fid0) # Close Original File
		close(fid) # Close New File
	end

	# Process CFG Files
	println("\nModifying CFG Files")
	for i=cfgfiles
		println("   " * splitdir(i)[2])
		fid1 = open(i,"r")
		fid2 = open(joinpath(f2,splitdir(i)[2]),"w+")

		for j=1:1:6 + maxanalogcols
			write(fid2,readline(fid1,keep=true))
		end

		for j=1:1:6 - maxanalogcols
			readline(fid1) # Skip line
			write(fid2,"Analog " * string(j+maxanalogcols) * ", 0 , 5000 , mV , 0 , D\n")
		end

		for j=1:1:9
			write(fid2,readline(fid1,keep=true))
		end

		close(fid1)
		close(fid2)
	end

	# Process CSR Files
	println("\nCopying CSR Files")
	for i=csrfiles
		println("   " * splitdir(i)[2])
		cp(i,joinpath(f2,splitdir(i)[2]))
	end

	# Process CSU Files
	println("\nCopying CSU Files")
	for i=csufiles
		println("   " * splitdir(i)[2])
		cp(i,joinpath(f2,splitdir(i)[2]))
	end

	# Process CSV Files
	println("\nCopying CSV Files")
	for i=csvfiles
		println("   " * splitdir(i)[2])
		cp(i,joinpath(f2,splitdir(i)[2]))
	end

	# Process LOG Files
	println("\nCopying LOG Files")
	for i=logfiles
		println("   " * splitdir(i)[2])
		cp(i,joinpath(f2,splitdir(i)[2]))
	end

	# Process FLX Files
	println("\nCopying FLX Files")
	for i=flxfiles
		println("   " * splitdir(i)[2])
		cp(i,joinpath(f2,splitdir(i)[2]))
	end
end # slt_trim(f1::String,f2::String,mindate::DateTime,maxdate::DateTime,maxanalogcols::Int)

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
module ReturnTypes

using BinaryTraits: has_proper_return_type
using Test
import Base: +

struct MyInt <: Integer
    value::Int
end
+(x::MyInt, y::Integer) = x.value + y

function test()
    # concrete return type, zero ambiguity
    f1(x::Int) = x + 1
    @test has_proper_return_type(f1, Tuple{Int}, Int)
    @test has_proper_return_type(f1, Tuple{Int}, Float64) == false

    # abstract type argument, multiple possible return types
    f2(x::Integer) = x + 1
    @test has_proper_return_type(f2, Tuple{Int}, Int)
    @test has_proper_return_type(f2, Tuple{Int8}, Int)     # due to promotion rule
    @test has_proper_return_type(f2, Tuple{UInt8}, Int)    # due to promotion rule
    @test has_proper_return_type(f2, Tuple{UInt64}, UInt64)
    @test has_proper_return_type(f2, Tuple{Integer}, Int) == false

    # One more level up, int's promoting to Float64.  Complex remains complex.
    f3(x::Number) = x + 1.0
    @test has_proper_return_type(f3, Tuple{Int}, Float64)       # promoted to Float64
    @test has_proper_return_type(f3, Tuple{Complex}, Complex)   # remains complex as input
    @test has_proper_return_type(f3, Tuple{Int}, Int) == false  # negative: no longer int

    # Using custom type

    # We know this has to be an Int
    @test has_proper_return_type(+, Tuple{MyInt,Int}, Int)

    # In this case, the returned type may not be Int because we don't have the concrete
    # type for the 2nd argument.  The result could be anything.
    @test has_proper_return_type(+, Tuple{MyInt,Integer}, Int) == false

    # Same as before, we cannot say the result would be an Integer either.
    @test has_proper_return_type(+, Tuple{MyInt,Integer}, Integer) == false
end

end # module

using Test
@testset "Return Types" begin
    import .ReturnTypes
    ReturnTypes.test()
end


text:
##
@inline function getproperty(p::T,f::Symbol) where T<:AbstractNode
    if ~hasfield(T,f)
        data = getfield(p,:data)
        haskey(data,f) ? data[f][getfield(p,:id)] : 0.0
    else
        getfield(p,f)
    end
end

@inline function setproperty!(p::T,f::Symbol,x::Float64) where T<:AbstractNode
    if ~haskey(getfield(p,:data),f)
        n = length(getfield(p,:data)[:w])
        getfield(p,:data)[f] = zeros(n)
    end
    getfield(p,:data)[f][getfield(p,:id)] = x
end

@inline -(n::T,x::NTuple{3,Float64}) where T<:AbstractNode = (n.x-x[1],n.y-x[2],n.z-x[3])
@inline -(x::NTuple{3,Float64},n::T) where T<:AbstractNode = (x[1]-n.x,x[2]-n.y,x[3]-n.z)
@inline +(n::T,x::NTuple{3,Float64}) where T<:AbstractNode = (n.x+x[1],n.y+x[2],n.z+x[3])
@inline +(x::NTuple{3,Float64},n::T) where T<:AbstractNode = (x[1]+n.x,x[2]+n.y,x[3]+n.z)
@inline *(x::NTuple{N,Float64},y::NTuple{N,Float64}) where N = sum(x[i]*y[i] for i in 1:N)

## ----------------- Node ------------------
struct Node<:AbstractNode
    id::Int
    data::Dict{Symbol,Vector{Float64}}
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
import LineSearches

let
   methods = [ LBFGS(), ConjugateGradient(),
               LBFGS(extrapolate=true, linesearch = LineSearches.interpbacktrack!) ]
   msgs = ["LBFGS Default Options: ",  "CG Default Options: ",
          "LBFGS + Backtracking + Extrapolation: " ]

   println("--------------------")
   println("Rosenbrock Example: ")
   println("--------------------")
   rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
   for (method, msg) in zip(methods, msgs)
      results = Optim.optimize(rosenbrock, zeros(2), method)
      println(msg, "g_calls = ", results.g_calls, ", f_calls = ", results.f_calls)
   end

   println("--------------------------------------")
   println("p-Laplacian Example (preconditioned): ")
   println("--------------------------------------")
   plap(U; n=length(U)) = (n-1) * sum( (0.1 + diff(U).^2).^2 ) - sum(U) / (n-1)
   plap1(U; n=length(U), dU = diff(U), dW = 4 * (0.1 + dU.^2) .* dU) =
                           (n-1) * ([0.0; dW] - [dW; 0.0]) - ones(U) / (n-1)
   precond(x::Vector) = precond(length(x))
   precond(n::Number) = spdiagm( ( -ones(n-1), 2*ones(n), -ones(n-1) ),
                                 (-1,0,1), n, n) * (n+1)
   df = DifferentiableFunction( X->plap([0;X;0]),
                                (X, G)->copy!(G, (plap1([0;X;0]))[2:end-1]) )
   GRTOL = 1e-6
   N = 100
   initial_x = zeros(N)
   P = precond(initial_x)
   methods = [ LBFGS(P=P), ConjugateGradient(P=P),
         LBFGS(extrapolate=true, linesearch = LineSearches.interpbacktrack!, P=P) ]

   for (method, msg) in zip(methods, msgs)
      results = Optim.optimize(df, copy(initial_x), method)
      println(msg, "g_calls = ", results.g_calls, ", f_calls = ", results.f_calls)
   end
end


text:
module RPCJuliaSer
using Base
using ZMQ

export
    #zmq_serialize_julia
    zmq_serialize,zmq_deserialize,
    #zmq_client_julia
    launch_client,zmqcall,zmqparse,zmqsetvar,zmqgetvar,
    #zmq_server_julia
    run_server,zmqquit,parse_eval

include("zmq_serialize_julia.jl")
include("zmq_client_julia.jl")
include("zmq_server_julia.jl")

end


text:
#!/bin/env julia

using ArgParse

function parse_commandline()
    s = ArgParseSettings(description = "usage: preprocess_fasta.jl fasta_file --output_dir dir")

    @add_arg_table s begin
        "--fasta_path"
          help = "fasta file path"
          required = true
        "--output_dir"
          help = "output directory"
          required = true
    end

    return parse_args(s)
end

function chomp_fasta(fa_file)
    d = Dict{String, String}()
    open(fa_file) do file
        i = 0
        contig = ""
        list = String[]
        while !eof(file)
            i += 1
            i % 1000_000 == 0 && info("preprocessed $i lines of $fa_file ")
            line = readline(file) |> chomp
            if startswith(line, ">chr")
                if length(list) > 0
                    d[contig] = join(list, "")
                end
                contig = line[2:end]
                list = String[]
            else
                push!(list, line)
            end
        end
    end
    d
end

function save_txt(data_dir, d)
    for contig in keys(d)
        if length(contig) > 5
            continue
        end
        s = d[contig]
        l = length(s)
        s_tr = s[1:div(l,2)]
        s_te = s[div(l,2)+1:end]
        write(joinpath(data_dir, string(contig, "_train.txt")), s_tr)
        write(joinpath(data_dir, string(contig, "_test.txt")),  s_te)
    end
end

function preprocess_fasta(data_dir, fa_file)
    d = chomp_fasta(fa_file)
    save_txt(data_dir, d)
end

function main()
    parsed_args = parse_commandline()
    for (arg,val) in parsed_args
        println("$arg => $val")
    end
    data_dir = parsed_args["output_dir"]
    fa_path  = parsed_args["fasta_path"]
    preprocess_fasta(data_dir, fa_path)
end

main()


text:
"""
    UnravelledLayout

Parameters are

- `x` (a function to get the value on x)
- `y` (a function to get the value on y)

Both of these functions *must* accept a unipartite network as input, and return
a dictionary with species and a single numerical value as output. Note that `x`
and/or `y` can be `λ`s.
"""
struct UnravelledLayout{TX<:Function,TY<:Function}
    x::TX
    y::TY
end

"""
    UnravelledLayout(;x=omnivory, y=fractional_trophic_level)

By default, the two functions used are omnivory and fractional trophic level,
but this can be changed.
"""
UnravelledLayout(;x=omnivory, y=fractional_trophic_level) = UnravelledLayout(x, y)

"""
    position!(LA::LT, L::Dict{K,NodePosition}, N::T) where {LT <: UnravelledLayout, T <: AbstractEcologicalNetwork} where {K}

Position species according to the function defined in the `UnravelledLayout`.
"""
function position!(LA::LT, L::Dict{K,NodePosition}, N::T) where {LT <: UnravelledLayout, T <: AbstractEcologicalNetwork} where {K}
  X = LA.x(N)
  Y = LA.y(N)
  for (i, s) in enumerate(species(N))
    L[s].x = X[s]
    L[s].y = Y[s]
  end
end

"""
    NestedBipartiteLayout

Parameters are

- `align` (whether the two levels should be centered together)
- `relative` (whether the two levels should occupy a length equal to their relative richness)
- `spread` (the distance between the two)

Note that to see the effect of `spread`, you may have to use `aspectratio=1`; if
not, the spacing between levels will be determined by the dimensions of the
plot.
"""
struct NestedBipartiteLayout
    align::Bool
    relative::Bool
    spread::Float64
end

"""
    NestedBipartiteLayout()

By default, a `NestedBipartiteLayout` is aligned, centered, and with a spread of
`1.0`.
"""
NestedBipartiteLayout() = NestedBipartiteLayout(true, true, 1.0)
NestedBipartiteLayout(spread::Float64) = NestedBipartiteLayout(true, true, spread)

"""
    position!(LA::NestedBipartiteLayout, L::Dict{K,NodePosition}, N::T) where {T <: AbstractBipartiteNetwork} where {K}

Rank species according to their degree.
"""
function position!(LA::NestedBipartiteLayout, L::Dict{K,NodePosition}, N::T) where {T <: AbstractBipartiteNetwork} where {K}
    r_top = ordinalrank(collect(values(degree(N; dims=1))); rev=true).-1
    r_bot = ordinalrank(collect(values(degree(N; dims=2))); rev=true).-1

    r_top = r_top./maximum(r_top)
    r_bot = r_bot./maximum(r_bot)
    if LA.relative
        if richness(N; dims=2)>richness(N; dims=1)
            r_bot = r_bot .* (richness(N; dims=2)./richness(N; dims=1))
        end
        if richness(N; dims=1)>richness(N; dims=2)
            r_top = r_top .* (richness(N; dims=1)./richness(N; dims=2))
        end
    end

    if LA.align
        r_bot = r_bot .+ (0.5 - (maximum(r_bot)-minimum(r_bot))/2.0)
        r_top = r_top .+ (0.5 - (maximum(r_top)-minimum(r_top))/2.0)
    end

    d_bot = Dict(zip(keys(degree(N; dims=2)), r_bot))
    d_top = Dict(zip(keys(degree(N; dims=1)), r_top))

    d_all = merge(d_bot, d_top)

    for s in species(N)
        L[s].x = d_all[s]
        L[s].y = LA.spread * L[s].y
    end
end


text:
using DynamicPolynomials
using MultivariateMoments
using SemialgebraicSets
using CSDP

# -----------------------------------------------------------------------------
# Jiawang Nie, Zi Yang, & Guangming Zhou. (2018).
# The Saddle Point Problem of Polynomials.
# -----------------------------------------------------------------------------

# Players are inverted!

@testset "Saddle Example 6.1 (II) " begin
        @polyvar x[1:3] y[1:3]
        p =
                x[1]^3 + x[2]^3 - x[3]^3 - y[1]^3 - y[2]^3 +
                y[3]^3 +
                x[3]y[1]y[2] * (y[1] + y[2]) +
                x[2]y[1]y[3] * (y[1] + y[3]) +
                x[1]y[2]y[3] * (y[2] + y[3])
        Sx = @set sum(x) == 1 && x[1] >= 0 && x[2] >= 0 && x[3] >= 0
        Sy = @set sum(y) == 1 && y[1] >= 0 && y[2] >= 0 && y[3] >= 0
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [WeightedDiracMeasure([0, 0, 1], 1)],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [WeightedDiracMeasure([0, 0, 1], 1)],
        )
        expected_value = 0

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x0)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x0)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Saddle Example 6.1 (III)" begin
        @polyvar x[1:4] y[1:4]
        p =
                sum([(x[i]^2 * y[i]^2) for i in 1:4]) - sum([
                        (i != j) ? (x[i]x[j] + y[i]y[j]) : 0
                        for i in 1:4, j in 1:4
                ])
        Sx = @set sum(x) == 1 &&
             x[1] >= 0 &&
             x[2] >= 0 &&
             x[3] >= 0 &&
             x[4] >= 0
        Sy = @set sum(y) == 1 &&
             y[1] >= 0 &&
             y[2] >= 0 &&
             y[3] >= 0 &&
             y[4] >= 0
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [WeightedDiracMeasure([0.25, 0.25, 0.25, 0.25], 1)],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [
                        WeightedDiracMeasure([1, 0, 0, 0], 0.25),
                        WeightedDiracMeasure([0, 1, 0, 0], 0.25),
                        WeightedDiracMeasure([0, 0, 1, 0], 0.25),
                        WeightedDiracMeasure([0, 0, 0, 1], 0.25),
                ],
        )
        expected_value = -0.6875

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Saddle Example 6.1 (IV) " begin
        @polyvar x[1:3] y[1:3]
        p =
                x[1]x[2]y[1]y[2] + x[2]x[3]y[2]y[3] + x[3]x[1]y[3]y[1] -
                x[1]^2 * y[3]^2 - x[2]^2 * y[1]^2 - x[3]^2 * y[2]^2
        Sx = @set sum(x) == 1 && x[1] >= 0 && x[2] >= 0 && x[3] >= 0
        Sy = @set sum(y) == 1 && y[1] >= 0 && y[2] >= 0 && y[3] >= 0
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [
                        WeightedDiracMeasure([0, 1, 0], 1 / 3),
                        WeightedDiracMeasure([1, 0, 0], 1 / 3),
                        WeightedDiracMeasure([0, 0, 1], 1 / 3),
                ],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [WeightedDiracMeasure([1 / 3, 1 / 3, 1 / 3], 1)],
        )
        expected_value = -1 / 9

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Saddle Example 6.3 (I)  " begin
        @polyvar x[1:3] y[1:3]
        p =
                sum([(x[i] + y[i]) for i in 1:3]) - prod([(x[i] - y[i]) for i in 1:3])
        Sx = @set x[1]^2 <= 1 && x[2]^2 <= 1 && x[3]^2 <= 1
        Sy = @set y[1]^2 <= 1 && y[2]^2 <= 1 && y[3]^2 <= 1
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [
                        WeightedDiracMeasure([-1, -1, 1], 1 / 3),
                        WeightedDiracMeasure([1, -1, -1], 1 / 3),
                        WeightedDiracMeasure([-1, 1, -1], 1 / 3),
                ],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [WeightedDiracMeasure([1, 1, 1], 1)],
        )
        expected_value = 2

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Saddle Example 6.3 (II)  " begin
        @polyvar x[1:3] y[1:3]
        p =
                y'y - x'x +
                sum([(x[i]y[j] - x[j]y[i]) for i = 1:3, j = 1:3 if i < j])
        Sx = @set x[1]^2 <= 1 && x[2]^2 <= 1 && x[3]^2 <= 1
        Sy = @set y[1]^2 <= 1 && y[2]^2 <= 1 && y[3]^2 <= 1
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [
                        WeightedDiracMeasure([-1.0, -1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, 1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, -1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, 1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, -1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, 1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, 1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, -1.0, 1.0], 1 / 8),
                ],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [
                        WeightedDiracMeasure([1.0, 1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, 1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, 1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, 1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, -1.0, -1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, -1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([-1.0, -1.0, 1.0], 1 / 8),
                        WeightedDiracMeasure([1.0, -1.0, -1.0], 1 / 8),
                ],
        )
        expected_value = 0

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Saddle Example 6.4 (I)  " begin
        @polyvar x[1:3] y[1:3]
        p =
                x[1]^3 +
                x[2]^3 +
                x[3]^3 +
                y[1]^3 +
                y[2]^3 +
                y[3]^3 +
                2 * (x[1]x[2]y[1]y[2] + x[1]x[3]y[1]y[3] + x[2]x[3]y[2]y[3])
        Sx = @set x'x == 1
        Sy = @set y'y == 1
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [
                        WeightedDiracMeasure([0, -1, 0], 1 / 3),
                        WeightedDiracMeasure([-1, 0, 0], 1 / 3),
                        WeightedDiracMeasure([0, 0, -1], 1 / 3),
                ],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [
                        WeightedDiracMeasure([0, 1, 0], 1 / 3),
                        WeightedDiracMeasure([1, 0, 0], 1 / 3),
                        WeightedDiracMeasure([0, 0, 1], 1 / 3),
                ],
        )
        expected_value = 0

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Saddle Example 6.5 (I)  " begin
        @polyvar x[1:3] y[1:3]
        p =
                x[1]^2 * y[1] + 2 * x[2]^2 * y[2] + 3 * x[3]^2 * y[3] - x[1] -
                x[2] - x[3]
        Sx = @set x'x <= 1
        Sy = @set y'y <= 1
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [WeightedDiracMeasure([0.726, 0.458, 0.349], 1)],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [WeightedDiracMeasure([0.688, 0.546, 0.477], 1)],
        )
        expected_value = -0.7666

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x0)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x0)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end