text:
## CrystalNets main API

"""
    topological_genome(net::CrystalNet{D,T})::String where {D,T}

Compute the topological genome of a net. The topological genome is an invariant
if the net, meaning that it does not depend on its representation. It is also
the string representation of a D-periodic graph such that
`PeriodicGraph{D}(topological_genome(net))` is isomorphic to `net.pge.g` (except
possibly if the `ignore_types` option is unset).

!!! info
    Options must be passed directly within `net`.
"""
function topological_genome(net::CrystalNet{D,T})::TopologicalGenome where {D,T}
    isempty(net.types) && return TopologicalGenome(net.options.error)
    if net.options.ignore_types
        net = CrystalNet{D,T}(net.pge, fill(Symbol(""), length(net.types)), net.options)
    end
    shrunk_net, _collisions = collision_nodes(net)
    _collisions isa Nothing && return TopologicalGenome(net.pge.g, nothing, true)
    collisions::Vector{CollisionNode} = _collisions

    if !net.options.skip_minimize
        flag = true
        try
            shrunk_net, _collisions = minimize(shrunk_net, collisions)
            flag = false
        catch e
            isinterrupt(e) && rethrow()
            if T == Rational{BigInt} || !isoverfloworinexact(e)
                net.options.throw_error && rethrow()
                return TopologicalGenome(string(e)::String)
            end
        end
        if flag # not in the catch to avoid a StackOverflow of errors in case something goes wrong
            newnet = CrystalNet{D,widen(T)}(net; ignore_types=false)
            return topological_genome(newnet)
        end
        _collisions isa Nothing && return TopologicalGenome(shrunk_net.pge.g, nothing, true)
        collisions = _collisions
    end

    export_net = isempty(net.options.export_net) ? isempty(net.options._pos) ?
                 net.options.export_subnets : "" : net.options.export_net
    export_default(net, "net", net.options.name, export_net)

    return topological_genome(shrunk_net, collisions)
end

topological_genome(net::CrystalNet{0,T}) where {T} = TopologicalGenome(net.options.error)

function topological_genome(net::CrystalNet{D,T}, collisions::Vector{CollisionNode})::TopologicalGenome where {D,T}
    try
        g::PeriodicGraph{D} = topological_key(net, collisions)
        unstable = g.width[] == -2
        unstable && (g.width[] = -1)
        ne(g) == 0 && return TopologicalGenome(net.options.error)
        return TopologicalGenome(g, recognize_topology(g), unstable)
    catch e
        isinterrupt(e) && rethrow()
        if T == Rational{BigInt} || !isoverfloworinexact(e)
            net.options.throw_error && rethrow()
            return TopologicalGenome(string(e)::String)
        end
    end
    return topological_genome(CrystalNet{D,widen(T)}(net), collisions)
end

"""
    topological_genome(g::Union{String,PeriodicGraph}, options::Options=Options())
    topological_genome(g::Union{String,PeriodicGraph}; kwargs...)

Compute the topological genome of a periodic graph.
If given a topological key (as a string), it is converted to a `PeriodicGraph` first.
"""
function topological_genome(g::PeriodicGraph, options::Options)
    nets = UnderlyingNets(g, options)
    return topological_genome(nets)
end
topological_genome(s::String, options::Options) = topological_genome(PeriodicGraph(s), options)
topological_genome(g::Union{String,PeriodicGraph}; kwargs...) = topological_genome(g, Options(; throw_error=true, kwargs...))

function _loop_group!(ex::Expr, id::Symbol, net::Symbol, group::Expr)
    for j in 1:length(ex.args)
        arg = ex.args[j]
        if arg isa Symbol
            if arg === :id
                ex.args[j] = id
            elseif arg === :net
                ex.args[j] = net
            elseif arg === :group
                ex.args[j] = group
            end
        elseif arg isa Expr
            _loop_group!(arg, id, net, group)
        end
    end
    nothing
end
macro loop_group(ex)
    @toggleassert ex.head === :for
    ret = Expr(:block)
    for i in 1:3
        id = Symbol("id", i)
        net = Symbol("net", i)
        D = Symbol("D", i)
        group = :(group.$D)
        newex = quote
            if !isempty(group)
                currallowed = first(group)[2][1].options.dimensions
                if isempty(currallowed) || $i in currallowed
                    $(deepcopy(ex))
                end
            end
        end
        _loop_group!(newex, id, net, group)
        push!(ret.args, esc(newex))
    end
    return ret
end

"""
    topological_genome(group::UnderlyingNets)::Vector{Tuple{Vector{Int},String}}

Compute the topological genome of each subnet stored in `group`.

!!! info
    Options must be passed directly within the subnets.
"""
function topological_genome(group::UnderlyingNets)
    ret = Tuple{Vector{Int},TopologyResult}[]
    @loop_group for (id, net) in group
        encountered = Dict{PeriodicGraph,_Clustering}()
        subret = Vector{Tuple{_Clustering,Union{_Clustering,TopologicalGenome}}}(undef, length(net))
        for (j, subnet) in enumerate(net)
            clust = only(subnet.options.clusterings)
            refclust = get!(encountered, subnet.pge.g, clust)
            subret[j] = (clust, refclust == clust ? topological_genome(subnet) : refclust)
        end
        push!(ret, (id, TopologyResult(subret)))
    end
    return ret
end

"""
    recognize_topology(g::PeriodicGraph, arc=CRYSTAL_NETS_ARCHIVE)
    recognize_topology(genome::AbstractString, arc=CRYSTAL_NETS_ARCHIVE)

Attempt to recognize a topological genome from an archive of known genomes.

!!! warning
    This function does a simple lookup, not any kind of topology computation.
    To identify the topology of a `PeriodicGraph` or a `CrystalNet` `x`, query
    `topological_genome(x)` instead.
"""
function recognize_topology(genome::AbstractString, arc=CRYSTAL_NETS_ARCHIVE)
    get(arc, genome, nothing)
end
function recognize_topology(genome::PeriodicGraph, arc=CRYSTAL_NETS_ARCHIVE)
    recognize_topology(string(genome), arc)
end

"""
    determine_topology(path, options::Options)
    determine_topology(path; kwargs...)

Compute the topology of the structure described in the file located at `path`.
This is essentially equivalent to calling
`topological_genome(UnderlyingNets(parse_chemfile(path, options)))`.

In the case where the structure is not made of interpenetrating nets, return the topology
of the only net.
"""
function determine_topology(path, options::Options)
    genomes::Vector{Tuple{Vector{Int},TopologyResult}} =
        topological_genome(UnderlyingNets(parse_chemfile(path, options)))
    if length(genomes) == 1
        return genomes[1][2]
    end
    length(genomes) == 0 && return TopologyResult("")
    return genomes
end
determine_topology(path; kwargs...) = determine_topology(path, Options(; kwargs...))

macro ifvalidgenomereturn(opts, msg, skipcrystal=false)
    crystaldef = skipcrystal ? nothing : :(crystal = parse_chemfile(path, $opts))
    msg = lazy"(found by $msg)"
    ifprintinfo = isempty(msg) ? nothing : :(CrystalNets.@ifwarn @info $msg)

    return esc(quote
    $crystaldef
    group = UnderlyingNets(crystal)
    dim, subnets = isempty(group.D3) ? isempty(group.D2) ? (1, group.D1) : (2, group.D2) : (3, group.D3)
    for (_, nets) in subnets
        for net in nets
            net isa CrystalNet || continue
            sig = get(encountered_graphs, net.pge.g, net.pge.g)
            if haskey(encountered_genomes, net.pge.g)
                unstable, counter = encountered_genomes[sig]
                encountered_genomes[sig] = (unstable, counter + 1)
            else
                genome = topological_genome(net)::TopologicalGenome
                if dim == maxdim && genome.name !== nothing
                    $ifprintinfo
                    return genome
                end
                encountered_graphs[sig] = genome.genome
                encountered_genomes[genome.genome] = (genome.unstable, 1)
            end
        end
    end
    end)
end

"""
    guess_topology(path, options::Options)
    guess_topology(path; kwargs...)

Tries to determine the topology of the file at `path` by passing various options
(starting from the provided options if any) until finding a known topology.
If none is found, return the topological genome encountered most often through
the variation of options.
"""
function guess_topology(path, defopts::Options)
    maxdim = maximum(defopts.dimensions; init=0)
    encountered_graphs = Dict{PeriodicGraph,PeriodicGraph}()
    encountered_genomes = Dict{PeriodicGraph,Tuple{Bool,Int}}()
    if defopts.name == "unnamed"
        defopts = Options(defopts; name=splitext(splitdir(path)[2])[1])
    end

    @ifvalidgenomereturn defopts "using structure=StructureType.Guess" false

    # dryrun = Dict{Symbol,Union{Nothing,Set{Symbol}}}()
    # crystal = parse_chemfile(path, Options(defopts; dryrun))
    atoms = Set{Symbol}(crystal.types)

    if any(x -> (at = atomic_numbers[x]; (ismetal[at] | ismetalloid[at])), atoms)
        @ifvalidgenomereturn Options(defopts; structure=StructureType.MOF) "using MOF clusters"
    end

    defopts = Options(; structure=StructureType.Auto)
    # if defopts.cutoff_coeff == 0.75 # if default cutoff was used, try larger
    #     @ifvalidgenomereturn Options(defopts; cutoff_coeff=0.85) "using longer cutoff"
    # end
    # invalidatoms = union(get(dryrun, :invalidatoms, Set{Symbol}()),
    #                      get(dryrun, :suspect_smallbonds, Set{Symbol}()))
    # hashydrogens = :H ∈ atoms || :D ∈ atoms
    # if hashydrogens
    #     @ifvalidgenomereturn Options(defopts; ignore_atoms=(:H,:D)) "ignoring H"
    #     # delete!(invalidatoms, :H)
    #     # delete!(invalidatoms, :D)
    # end
    # for a in invalidatoms
    #     @ifvalidgenomereturn Options(defopts; ignore_atoms=tuple(a)) "ignoring $a"
    #     if hashydrogens
    #         @ifvalidgenomereturn Options(defopts; ignore_atoms=(a, :H, :D)) "ignoring H and $a"
    #     end
    # end
    # if any(==(:C), atoms)# && :C ∉ invalidatoms # organic solvent
    #     @ifvalidgenomereturn Options(defopts; ignore_atoms=(:C,)) "ignoring C"
    # end
    # if haskey(dryrun, :try_Input_bonds)
    #     @ifvalidgenomereturn Options(defopts; bonding=Bonding.Input) "using input bonds"
    # end
    # @ifvalidgenomereturn Options(defopts; ignore_low_occupancy=true) "removing atoms with occupancy < 0.5"
    if :Al ∈ atoms && (:P ∈ atoms || :Pc ∈ atoms) # ALPO datastructure
        # flag = false
        # for i in vertices(crystal.pge.g)
        #     crystal.types[i] === :O || continue
        #     currt = :O
        #     for x in neighbors(crystal.pge.g, i)
        #         t = crystal.types[x.v]
        #         if (t === :P) | (t === :Al)
        #             if currt === :O
        #                 currt = t
        #             elseif currt === t
        #                 flag = true
        #                 break
        #             end
        #         end
        #     end
        # end
        # if flag
            @ifvalidgenomereturn Options(defopts; ignore_homoatomic_bonds=(:Al,:P)) "removing Al-O-Al and P-O-P bonds"
        # end
        # if :Sn ∈ atoms # observed in some cases for ALPO databases
        #     @ifvalidgenomereturn Options(defopts; ignore_atoms=(:Sn,)) "ignoring Sn"
        # end
    end
    # if haskey(dryrun, :try_no_Auto_bonds)
    #     @ifvalidgenomereturn Options(defopts; bonding=Bonding.Guess) "discarding input bonds and guessing them"
    # end
    # if haskey(dryrun, :collisions)
    #     @ifvalidgenomereturn Options(defopts; authorize_pruning=false) "retaining all colliding atoms"
    # end

    # Finally, if everything fails, return the one encountered most
    maxencounter = 0
    most_plausible_genome::PeriodicGraph = PeriodicGraph{0}()
    final_unstable = false
    for (genome, (unstable, i)) in encountered_genomes
        if i > maxencounter
            maxencounter = i
            most_plausible_genome = genome
            final_unstable = unstable
        end
    end
    return most_plausible_genome == PeriodicGraph{0}() ? begin
            issetequal(defopts.dimensions, [1,2,3]) ? TopologicalGenome() :
                TopologicalGenome("could not guess (perhaps because no net with suitable dimension in $(defopts.dimensions) was found?)")
            end : TopologicalGenome(most_plausible_genome, nothing, final_unstable)
end
guess_topology(path; kwargs...) = guess_topology(path, Options(structure=StructureType.Guess; kwargs...))

"""
    determine_topology_dataset(path, save, autoclean, options::Options)
    determine_topology_dataset(path, save=true, autoclean=true; kwargs...)

Given a path to a directory containing structure input files, compute the
topology of each structure within the directory.
Return a dictionary linking each file name to the result.
The result is a [`TopologyResult`](@ref), containing the topological genome, the name if
known and the stability of the net. In case of error, the exception is reported.

It is strongly recommended to toggle warnings off (through [`toggle_warning`](@ref)) and
not to export any file since those actions may critically reduce performance,
especially for numerous files.

If `save` is set, the result is also stored in a julia serialized file located at
"\\\$path/../results_\\\$i" where `i` is the lowest integer such that this path does
not already exist at the start of the computation.
While processing, this path will be used to create a directory storing the
current state of the computation: to continue an interrupted computation, simply
pass this temporary directory as the path. If `autoclean` is set, this directory
is removed at the end if the computation was successful.

If `save` is set and `autoclean` is unset, the directory of temporary files will
be renamed into "\\\$path/../results_\\\$i.OLD\\\$j".
"""
function determine_topology_dataset(path, save, autoclean, options::Options)
    if isdirpath(path)
        path = dirname(path)
    end
    if startswith(basename(path), "results_") && isfile(joinpath(path), "data")
        resultdir = path
        path = only(readlines(joinpath(path, "data")))
        alreadydone = String[]
        for _f in readdir(resultdir; join=true)
            basename(_f) == "data" && continue
            io = open(_f, "r+")
            pos = max(filesize(io) - 2, 0)
            seek(io, pos)
            while pos > 0 && read(io, Char) != '\n'
                pos -= 1
                seek(io, pos)
            end
            truncate(io, pos+(pos>0)) # remove the last line if incomplete
            close(io)
            for l in eachline(_f)
                _splits = split(l, '/')
                isempty(last(_splits)) && continue
                push!(alreadydone, join(_splits[1:end-2], '/'))
            end
        end
        files = recursive_readdir(path)
        setdiff!(files, alreadydone)
    else
        resultdir = tmpexportname(dirname(path), "", "results", "")
        mkdir(resultdir)
        open(joinpath(resultdir, "data"), "w") do f
            println(f, path)
        end
        files = recursive_readdir(path)
    end
    resultdir::String
    files::Vector{String}

    @threads for file in files
        f = joinpath(path, file)
        # threadid() == 1 && @show f # to find infinite loops: the last one printed is probably running

        genomes::Vector{Tuple{Vector{Int},TopologyResult}} = try
            topological_genome(UnderlyingNets(parse_chemfile(f, options)))
        catch e
            (options.throw_error || isinterrupt(e)) && rethrow()
            [(Int[], TopologyResult(string(e)))]
        end
        if isempty(genomes)
            push!(genomes, (Int[], TopologyResult("")))
        end
        for (j, (_, genome)) in enumerate(genomes)
            newname = length(genomes) == 1 ? file * '/' : file * '/' * string(j)
            open(joinpath(resultdir, string(threadid())), "a") do results
                io = IOContext(results, :compact => true)
                println(io, newname, '/', genome)
            end
        end
    end

    ret = Pair{String,TopologyResult}[]
    for _f in readdir(resultdir; join=true)
        basename(_f) == "data" && continue
        for l in eachline(_f)
            isempty(l) && continue
            splits = split(l, '/')
            _genome = pop!(splits)
            isempty(splits[end]) && pop!(splits)
            push!(ret, Pair(join(splits, '/'), parse(TopologyResult, _genome)))
        end
    end
    result::Dict{String,TopologyResult} = Dict(ret)
    if save
        i = 0
        tmpresultdir = resultdir*".OLD"*string(i)
        while ispath(tmpresultdir)
            i += 1
            tmpresultdir = resultdir*".OLD"*string(i)
        end
        mv(resultdir, tmpresultdir)
        Serialization.serialize(resultdir, result)
        println("Topologies of ", path, " saved at ", resultdir)
        if autoclean
            rm(tmpresultdir; recursive=true)
        else
            println("Temporary files kept at ", tmpresultdir)
        end
    elseif autoclean
        rm(resultdir; recursive=true)
    else
        println("Temporary files kept at ", resultdir)
    end
    return result
end
function determine_topology_dataset(path, save=true, autoclean=true; kwargs...)
    opts, restore_warns = db_options(; kwargs...)
    ret = determine_topology_dataset(path, save, autoclean, opts)
    restore_warns && (DOWARN[] = true)
    ret
end

"""
    guess_topology_dataset(path, save, autoclean, options::Options)
    guess_topology_dataset(path, save=true, autoclean=true; kwargs...)

Given a path to a directory containing structure input files, guess the topology of each
structure within the directory using [`guess_topology`](@ref).
Return a dictionary linking each file name to the result.
The result is the corresponding topology name, if known, or the topological
genome preceded by an "UNKNOWN" mention otherwise. In case of error, the result
is the exception preceded by a "FAILED with" mention. Finally, if the input does
not represent a periodic structure, the result is "non-periodic".

It is strongly recommended to toggle warnings off (through [`toggle_warning`](@ref)) and
not to export any file since those actions may critically reduce performance,
especially for numerous files.

The `save` and `autoclean` arguments work identically to their counterpart for
[`determine_topology_dataset`](@ref).
"""
function guess_topology_dataset(path, save, autoclean, options::Options)
    if isdirpath(path)
        path = dirname(path)
    end
    if startswith(basename(path), "guessed_") && isfile(joinpath(path), "data")
        resultdir = path
        path = only(readlines(joinpath(path, "data")))
        alreadydone = String[]
        for _f in readdir(resultdir; join=true)
            basename(_f) == "data" && continue
            io = open(_f, "r+")
            pos = max(filesize(io) - 2, 0)
            seek(io, pos)
            while pos > 0 && read(io, Char) != '\n'
                pos -= 1
                seek(io, pos)
            end
            truncate(io, pos+(pos>0)) # remove the last line if incomplete
            close(io)
            for l in eachline(_f)
                isempty(l) && continue
                _splits = split(l, '/')
                push!(alreadydone, join(_splits[1:end-1], '/'))
            end
        end
        files = recursive_readdir(path)
        setdiff!(files, alreadydone)
    else
        resultdir = tmpexportname(dirname(path), "", "guessed", "")
        mkdir(resultdir)
        open(joinpath(resultdir, "data"), "w") do f
            println(f, path)
        end
        files = recursive_readdir(path)
    end
    resultdir::String
    files::Vector{String}

    @threads for file in files
        f = joinpath(path, file)

        genome::TopologicalGenome = try
            guess_topology(f, options)
        catch e
            (options.throw_error || isinterrupt(e)) && rethrow()
            TopologicalGenome(escape_string(string(e)::String))
        end
        open(joinpath(resultdir, string(threadid())), "a") do results
            println(results, file, '/', genome)
        end
    end

    ret = Pair{String,TopologicalGenome}[]
    for _f in readdir(resultdir; join=true)
        basename(_f) == "data" && continue
        for l in eachline(_f)
            splits = split(l, '/')
            isempty(last(splits)) && @show _f, l
            _genome = pop!(splits)
            push!(ret, Pair(join(splits, '/'), parse(TopologicalGenome, _genome)))
        end
    end
    result::Dict{String,TopologicalGenome} = Dict(ret)
    if save
        i = 0
        tmpresultdir = resultdir*".OLD"*string(i)
        while ispath(tmpresultdir)
            i += 1
            tmpresultdir = resultdir*".OLD"*string(i)
        end
        mv(resultdir, tmpresultdir)
        Serialization.serialize(resultdir, result)
        println("Guessed topologies of ", path, " saved at ", resultdir)
        if autoclean
            rm(tmpresultdir; recursive=true)
        else
            println("Temporary files kept at ", tmpresultdir)
        end
    elseif autoclean
        rm(resultdir; recursive=true)
    else
        println("Temporary files kept at ", resultdir)
    end
    return result
end
function guess_topology_dataset(path, save=true, autoclean=true; kwargs...)
    opts, restore_warns = db_options(; structure=StructureType.Guess, kwargs...)
    ret = guess_topology_dataset(path, save, autoclean, opts)
    restore_warns && (DOWARN[] = true)
    ret
end


text:
export AffineScalingLayer, getAffineScalingLayer

"""
Scales and shifts the 3D feature tensor along each dimension. This is
useful, e.g., in batch normalization.

kron(s3,kron(s2,s1)) * vec(Y) + kron(b3,kron(e2,e1)) +
kron(e3,kron(b2,e1)) + kron(e3,kron(e2,b1));
"""
mutable struct AffineScalingLayer{T} <: AbstractMeganetElement{T}
    nData  :: Array{Int,1}     # describe size of data, at least first two dim must be correct.
end
function getAffineScalingLayer(TYPE::Type, nData)
    return AffineScalingLayer{TYPE}(nData)
end

function splitWeights(this::AffineScalingLayer{T},theta_in::Array{T}) where {T <: Number}
    theta = reshape(theta_in,:,2)
    s2    = theta[:,1]
    b2    = theta[:,2]
    return s2, b2
end

function scaleChannels!(Y::Array{T},s::Array{T},b::Array{T}) where {T <: Number}
    for i=1:length(s)
        Y[:,i,:] .= Y[:,i,:].*s[i] .+ b[i]
    end
end

function apply(this::AffineScalingLayer{T},theta::Array{T},Y::Array{T},dA,doDerivative=false) where {T <: Number}

    Y   = reshape(copy(Y),this.nData[1], this.nData[2],:)
    dA  = Array{T,2}(0,0)
    nex = size(Y,3)

    s2,b2 = splitWeights(this,theta);

    scaleChannels!(Y,s2,b2);

    Yout = reshape(Y,:,nex)
    Ydata = Yout
    return Ydata, Yout, dA
end

function nTheta(this::AffineScalingLayer)
    return 2*this.nData[2]
end

function nFeatIn(this::AffineScalingLayer)
    return prod(this.nData[1:2])
end

function nFeatOut(this::AffineScalingLayer)
   return prod(this.nData[1:2])
end

function nDataOut(this::AffineScalingLayer)
    return nFeatOut(this);
end

function initTheta(this::AffineScalingLayer{T}) where {T <: Number}
    s2,b2 = splitWeights(this,ones(T,nTheta(this)))
    return  [vec(s2); convert(T,0.)*vec(b2)]
end

function Jthetamv(this::AffineScalingLayer{T},dtheta::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    Y   = reshape(copy(Y),this.nData[1], this.nData[2],:)
    nex = size(Y,3)

    ds2,db2 = splitWeights(this,dtheta)

    scaleChannels!(Y,ds2,db2)

    dY = reshape(Y,:,nex)
    dYdata = dY
    return dYdata, dY
end

function JthetaTmv(this::AffineScalingLayer{T},Z::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    Y   = reshape(Y,this.nData[1], this.nData[2],:)
    Z   = reshape(Z,this.nData[1], this.nData[2],:)

    W = Y.*Z

    dtheta = vec(sum(sum(W,1),3))
    return  [dtheta; vec(sum(sum(Z,1),3))]
end

function JYmv(this::AffineScalingLayer{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    dY   = reshape(copy(dY),this.nData[1], this.nData[2],:);
    nex = size(dY,3)

    s2,b2 = splitWeights(this,theta);
    scaleChannels!(dY,s2,b2*0)

    dY     = reshape(dY,:,nex)
    dYdata = dY
    return dYdata, dY
end

function JYTmv(this::AffineScalingLayer{T},Z::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    Z   = reshape(copy(Z),this.nData[1], this.nData[2],:)
    nex = size(Z,3)

    s2,b2 = splitWeights(this,theta)
    scaleChannels!(Z,s2,b2*0)
    return reshape(Z,:,nex)
end

ext: jl


text:
export Pinhole

struct Pinhole <: CameraModel

    Pinhole(pixelwidth::Int, pixelheight::Int, prinicipalpoint::PixelCoordinate, focallength::Vector2) = new(prinicipalpoint, focallength, pixelwidth, pixelheight)

    prinicipalpoint::PixelCoordinate # in pixels
    focallength::Vector2 # in pixels
    columns::Int
    rows::Int
end

"""
    point2pixel(model::Pinhole, pointincamera::$(Point3))

Return a transformation that converts real-world coordinates
to camera coordinates. This currently ignores any tangential
distortion between the lens and the image plane.
"""
function point2pixel(model::Pinhole, pointincamera::Point3)
    column = model.prinicipalpoint[1] + model.focallength[1] * pointincamera[1] / pointincamera[2]
    row = model.prinicipalpoint[2] - model.focallength[2] * pointincamera[3] / pointincamera[2]
    return PixelCoordinate(column, row)
end

"""
    pixel2ray(model::Pinhole, pixelcoordinate::$(PixelCoordinate))

Return a transformation that converts real-world coordinates
to camera coordinates. This currently ignores any tangential
distortion between the lens and the image plane.
"""
function pixel2ray(model::Pinhole, pixelcoordinate::PixelCoordinate)
    x = (pixelcoordinate[1] - model.prinicipalpoint[1]) / model.focallength[1]
    z = -(pixelcoordinate[2] - model.prinicipalpoint[2]) / model.focallength[2]
    return Vector3(x, 1, z)
end

lookdirection(cameramodel::Pinhole) = SVector{3}(0,1,0)
updirection(cameramodel::Pinhole) = SVector{3}(0,0,1)

columns(cameramodel::Pinhole) = cameramodel.columns
rows(cameramodel::Pinhole) = cameramodel.rows


text:
reflo = [-6.348709539607487e-14, -4.232473026404991e-14, -2.116236513202496e-14, -1.5871773849018717e-14,
-1.6929892105619966e-14, -1.058118256601248e-14, -9.069585056582124e-15, -7.935886924509359e-15,
-7.054121710674986e-15, -8.464946052809983e-15, -5.7715541269158975e-15, -5.29059128300624e-15,
-4.88362272277499e-15, -3.0231950188607086e-15, -2.821648684269994e-15, -3.967943462254679e-15,
-2.489690015532348e-15, -2.3513739035583286e-15, -2.2276173823184166e-15, -3.1743547698037436e-15,
-2.0154633459071388e-15, -2.4048142195482905e-15, 0.0, -1.7635304276687464e-15,
-1.6929892105619966e-15, 0.0, 0.0, -1.1336981320727656e-15,
-7.297367286905158e-16, 7.054121710674985e-16, 0.0, 0.0,
-6.412837918795442e-16, 0.0, -6.046390037721416e-16, 1.1756869517791643e-15,
0.0, 1.1138086911592083e-15, 0.0, 1.0581182566012479e-15,
2.06462098849024e-15, 0.0005794956745071824, 9.842960526523235e-16, 0.0,
0.0015100291800084077, 0.010170231329438926, 0.025728741385843894, 0.04005101497165152,
0.04390452127911886, 0.04525686807986621, 0.03530410264550449, 0.021670517548941896,
0.008048864519123124, 0.0010780640392045314, 0.004470323125925969, 0.016688348116340457,
0.03493460910804327, 0.05959215527770473, 0.08244270095710961, 0.10357659617076827,
0.11391336511012921, 0.11627478138119568, 0.12956428457670546, 0.13663776226609609,
0.13497016376638563, 0.13249685947909343, 0.12925392578276643, 0.12817322381349877,
0.13203388059095647, 0.13784600161681598, 0.14240723449087092, 0.13937425900528344,
0.14236086920574512, 0.15049881048630848, 0.16302536651810398, 0.17166029310408853,
0.15680191863869608, 0.14989926543222462, 0.14999356535627262, 0.15783399886553232,
0.17387732650281687, 0.18996005145684927, 0.20460513261203542, 0.2017828459833418,
0.203590444172436, 0.20965401643308929, 0.21978961885490292, 0.23334414164750214,
0.24482639625816902, 0.24450009814539006, 0.241361704385313, 0.2373974635454126,
0.23582082671730284, 0.24335346466686872, 0.25701324342128723, 0.2612141132028289,
0.2660305729207295, 0.2692404779242703, 0.263964354375628, 0.25885651869294757,
0.25754979365637176, 0.254288196603691, 0.2503362377317122, 0.2568193720912553,
0.2574604980989164, 0.2485319153927413, 0.24205644600344794, 0.23863075415041338,
0.23819906143846187, 0.23908594757055368, 0.24173305614768417, 0.24770712558590272,
0.24877543972694863, 0.24892880078238666, 0.2590256898954135, 0.27166808320495917,
0.28251928763456524, 0.29162484188457044, 0.2990287539426958, 0.3008865750588349,
0.28577157939536907, 0.2773551436197429, 0.268485802151751, 0.26277533660806096,
0.2633127557034817, 0.2663000601279467, 0.2702453624196831, 0.2753066947598586,
0.2822525773726178, 0.27968003222135107, 0.27684017791426246, 0.2737971896326203,
0.26963063318403047, 0.2643656630796935, 0.25965299728293023, 0.2585011820584089,
0.24854591958584712, 0.2430030832516569, 0.23904201480796408, 0.23985015545585645,
0.24728657601452025, 0.2569315745253772, 0.266520726664246, 0.2746487420772013,
0.28134585117926253, 0.2866414561548989, 0.2920067929281561, 0.28967430030982805,
0.2880077439312323, 0.288510829916702, 0.28919779478701874, 0.29278856854550295,
0.2997378763748437, 0.3127219617756912, 0.3078585682435972, 0.30763813418697433,
0.3123729588390938, 0.32012974068543115, 0.3302985368071382, 0.33922797405322,
0.3400150482124273, 0.33663341679143066, 0.3339007121337234, 0.3298189861614539,
0.3244127665348869, 0.32058884855731573, 0.32118808369256424, 0.32630111972700526,
0.33284953446230253, 0.32926421804558403, 0.3327651361236813, 0.3426218018642772,
0.3423855654841696, 0.3306579511745207, 0.31923898910700166, 0.3118223610448929,
0.3065518776552234, 0.3071122760309198, 0.31401481781361384, 0.32169601368595346,
0.32988432780701954, 0.3372745855961669, 0.34548506656580164, 0.35246577682985175,
0.35396308945843974, 0.3577273271576752, 0.363231053913255, 0.3638967738452223,
0.3588223472785071, 0.35554231180656304, 0.35943155749311084, 0.3666088366234249,
0.36259249062210475, 0.3487658723827629, 0.3359619736344609, 0.3223877605974673,
0.3106225302007395, 0.30296375784327684, 0.3014745016770177, 0.3046726127337049,
0.30988668050762663, 0.317728188762448, 0.3289511297424724, 0.34373311785458177,
0.36278682227654396, 0.3827612368184497, 0.4030778282669414, 0.41120559536292073,
0.41055220793579156, 0.4122571066869349, 0.40992614632401275, 0.4048223000408057,
0.400321623164694, 0.39588340931483806, 0.39313314088683077, 0.3969234373130931,
0.4060619592316241, 0.41377185912011777, 0.41714203908292624, 0.4150723889033669,
0.4100328514069579, 0.4052900191597187, 0.4002587071424406, 0.3989137429103623,
0.4031942481568127, 0.4077209702428449, 0.4099745452392779, 0.41540722380025175,
0.42153563465409927, 0.42505338976043805, 0.4259831673784225, 0.427143466428484,
0.4283866422609185, 0.43082670612005863, 0.4255250622127837, 0.4208461818605876,
0.41818437287732685, 0.4170953583386116, 0.4144990630503587, 0.41186678541389793,
0.4115996854298485, 0.4056961954038659, 0.40098900289946593, 0.39496848014113567,
0.38668681168778507, 0.3834726184828498, 0.3764527211018305, 0.36938496789459807,
0.3630576189188345, 0.3568556984778469]
refhi = [31.917204450284782, 16.932239511551998, 11.994517392552043, 9.584217728451947,
8.10628782192911, 7.0975967411252965, 6.39645666944237, 5.901264162321405,
5.495463556286324, 5.152237281011368, 4.879784949500827, 4.673054226704583,
4.482512483622468, 4.303983206813759, 4.135063903023782, 3.9739527023115193,
3.819271704133845, 3.669949208068331, 3.525139134834113, 3.4104134778162924,
3.3709627094430035, 3.322888194788645, 3.2673147703246475, 3.2051797997772273,
3.1372706686763916, 3.064253626316756, 2.986696218393756, 2.905084912646316,
2.8198721542188956, 2.7410314167421053, 2.672760581622756, 2.599502482784571,
2.542903555664134, 2.5201546408800395, 2.4909291751605025, 2.46347544713647,
2.440203948947183, 2.4095602883242973, 2.372111554685783, 2.3283681285078117,
2.35110533491961, 2.3666665009529186, 2.3725651057534756, 2.369459960314425,
2.35795131465174, 2.3439198252346114, 2.3442313576669997, 2.3521391002164744,
2.3559707670338503, 2.351344833824406, 2.3387588064720712, 2.3186719211774,
2.3021433105646447, 2.2929166351323835, 2.2837575529773844, 2.2720715966501337,
2.253451179372895, 2.228254980160326, 2.206621825538713, 2.2265381069035604,
2.237031476008932, 2.238557880222222, 2.2315443178716206, 2.223419992092579,
2.2212293720727287, 2.236721058631896, 2.244141469485805, 2.243846690327335,
2.2361721642005024, 2.2214341659753796, 2.2003537100086477, 2.1771781207007326,
2.153273120385483, 2.125323965857102, 2.098700759718762, 2.0992573478955996,
2.0977981919901754, 2.097206235543601, 2.096807455783442, 2.0954367698413843,
2.092328696672403, 2.0822702818778316, 2.0655127426248576, 2.0422953333584237,
2.0128460494907716, 1.996796743078305, 2.00646971273678, 2.01707230245634,
2.026865130062295, 2.0295172353930786, 2.0283109492751628, 2.0237732130151014,
2.016683298268531, 2.0057288593312004, 1.9936220145826855, 1.9827997217856215,
1.9705706636183393, 1.957415091422549, 1.9454072743348574, 1.9368459388294168,
1.9305387241548866, 1.9216070890633343, 1.9073061677044048, 1.8922214422747685,
1.8795088464835308, 1.8716751757925083, 1.8686801059269904, 1.8646126604756152,
1.8558735781925564, 1.8402669581694788, 1.8201274172345854, 1.8176355685089265,
1.8283580183197359, 1.8327988423536985, 1.8311219091093416, 1.8252298709241315,
1.818882200265471, 1.81559754774534, 1.8133471950068, 1.8219414394127593,
1.8277065349956902, 1.8333163356276947, 1.8388164509301608, 1.84633128829015,
1.855487774207935, 1.8620165938925903, 1.8650352280505633, 1.8683265471733517,
1.8674292766828489, 1.8601234587815316, 1.8509401695484171, 1.844702839971218,
1.8379210455386354, 1.8247614914050734, 1.8053659055638835, 1.8184014787186613,
1.8323424752884971, 1.8401158276717893, 1.8493748344785492, 1.8605959624349366,
1.8685650199728956, 1.8753642360706488, 1.8801785487437959, 1.885042242340755,
1.887752219133277, 1.8847612067297628, 1.880587087952027, 1.8760236768533078,
1.8655828807096093, 1.8489724782986583, 1.826315044579408, 1.7977299288539939,
1.7816463104590459, 1.8132756344201248, 1.8398738291418253, 1.8622333441064587,
1.8794098928940448, 1.892111677383889, 1.8990372292634377, 1.8991958175519492,
1.8972573023274302, 1.8989509741151542, 1.8977512335978843, 1.8915587466396786,
1.8875871110251927, 1.8822899594724387, 1.874052974783995, 1.8762388542760413,
1.8878371577582833, 1.897124829774482, 1.9036777374888811, 1.9046862594297183,
1.9008415274302985, 1.8986693010300657, 1.914332161615297, 1.9254872667173353,
1.933168874921557, 1.9357657226281264, 1.934925270313838, 1.9300725926121831,
1.9210576909678778, 1.907779822534061, 1.892625679521541, 1.8788674050783838,
1.866455958953833, 1.8529374204933733, 1.8593698049491787, 1.8637283539383553,
1.8616566994484653, 1.8545377325823087, 1.8442520218562986, 1.8348318206786332,
1.8249100871165258, 1.8246047978590818, 1.827126394852936, 1.8230377030675533,
1.8141675230797218, 1.8038493756323082, 1.7928085112374317, 1.7826935353028226,
1.7777073281145646, 1.786728804294096, 1.789577722515789, 1.7863448556907853,
1.781290916162904, 1.7773773528569843, 1.7715449378918218, 1.7647467432771569,
1.7553864886932242, 1.7445895103890794, 1.733718121237644, 1.7198809834856756,
1.7030749600161919, 1.6965150954625188, 1.6858632017765836, 1.693183827590814,
1.6986959961756778, 1.7005568919312795, 1.7093668548476197, 1.7292512835035259,
1.7485139381505965, 1.7695467080386964, 1.791645469388392, 1.811622337813538,
1.8275942143810395, 1.836771128131683, 1.8392204050291943, 1.8352940265818136,
1.8299966809921944, 1.8295471375152572, 1.8334033178873554, 1.8369247326917209,
1.8363162409416394, 1.8376515267951465, 1.8395185332753015, 1.8374792223593692,
1.8330884737642625, 1.8265021865684263, 1.8222991897515293, 1.8239147996140013,
1.824795151074386, 1.8216708845003273, 1.8275609800115848, 1.83711672477262,
1.8422019068491204, 1.84177423141554, 1.8357445090850546, 1.8254441682932698,
1.8118902523940805, 1.7995974466201832]

using PointPatternStatistics
using Random
window = (x=(0,1), y=(0,1))
r0 = 0:0.001:0.25
n = 2500
Random.seed!(1)
xys = [[(rand(), rand()) for _ in 1:30] for _ in 1:n]
a = [pcf(xy, window, r0)[2:end] for xy in xys]

lo, hi = globalenvelope(a)
@test refhi == hi
@test refhi == hi


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_gauge

"""
    daq_gauge(;kwargs...)

A Gauge component.
A Gauge component that points to
a value between some range.
Keyword arguments:
- `id` (String; optional): The ID used to identify this compnent in Dash callbacks
- `base` (Real; optional): Base to be used in logarithmic scale.
- `className` (String; optional): Class to apply to the root component element.
- `color` (optional): Color configuration for the gauge's track.. color has the following type: String | lists containing elements 'default', 'gradient', 'ranges'.
Those elements have the following types:
  - `default` (String; optional): Color used for current value text and other minor accents
  - `gradient` (Bool; optional): Display ranges as a gradient between given colors.
  - `ranges` (optional): Define multiple color ranges on the gauge's track.
The key determines the color of the range and
the value is the start,end of the range itself.
Ranges must be contiguous along the entirety
of the gauge's range of values.. ranges has the following type: lists containing elements 'color'.
Those elements have the following types:
  - `color` (Array of Reals; optional)
- `label` (optional): Description to be displayed alongside the control. To control styling, pass an object with label and style properties.. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the component label is positioned.
- `logarithmic` (Bool; optional): If set to true, a logarithmic scale will be
used.
- `max` (Real; optional): The maximum value of the gauge. If logarithmic,
represents the maximum exponent.
- `min` (Real; optional): The minimum value of the gauge. If logarithmic,
represents the minimum exponent.
- `scale` (optional): Configuration for the component scale.. scale has the following type: lists containing elements 'start', 'interval', 'labelInterval', 'custom'.
Those elements have the following types:
  - `start` (Real; optional): Value to start the scale from. Defaults
to min.
  - `interval` (Real; optional): Interval by which the scale goes up. Attempts
to dynamically divide min-max range by
default.
  - `labelInterval` (Real; optional): Interval by which labels are added to
scale marks. Defaults to 2 (every other
mark has a label).
  - `custom` (optional): Custom scale marks. The key determines the position
and the value determines what will show. If you want
to set the style of a specific mark point, the value
should be an object which contains style and label
properties. custom has the following type: Real | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (String; optional)
  - `label` (String; optional)
- `showCurrentValue` (Bool; optional): If true, the current value of the gauge
will be displayed
- `size` (Real; optional): The size (diameter) of the gauge in pixels
- `style` (Dict; optional): Style to apply to the root component element.
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
- `units` (String; optional): Label for the current value
- `value` (Real; optional): The value of gauge. If logarithmic, the displayed
value will be the logarithm of the inputted value.
"""
function daq_gauge(; kwargs...)
        available_props = Symbol[:id, :base, :className, :color, :label, :labelPosition, :logarithmic, :max, :min, :scale, :showCurrentValue, :size, :style, :theme, :units, :value]
        wild_props = Symbol[]
        return Component("daq_gauge", "Gauge", "dash_daq", available_props, wild_props; kwargs...)
end


text:
using Test
using SubModule

@test SubModule.f() == 1

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
#
#  The following functions allow us to get distances between points and such
#
#   getEuclideanDistance(pos1::Tuple, pos2::Tuple)
#   isclose(pos1, pos2, threshold)
#

function getEuclideanDistance(pos1::Tuple, pos2::Tuple)
    x1 = pos1[1];
    x2 = pos2[1];
    y1 = pos1[2];
    y2 = pos2[2];
    distance = sqrt((x2-x1)^2 + (y2-y1)^2)
    return distance
end;
function isclose(pos1, pos2, threshold)
    d = getEuclideanDistance(pos1, pos2)
    if d < threshold
        return true
    else
        return false
    end
end;


text:
using JuMP
using DelimitedFiles
import GLPK
import Cbc
# import Xpress

"""
For research purposes, a utility that finds the optimal stable matching using
integer programming, with commercial Xpress solver as default.
"""
function IP_stable_opt(students::Array{Int64, 2}, schools::Array{Int64, 2},
                    capacities::Array{Int64, 1})
    n, m = size(schools)
    @assert (m,) == size(capacities)
    @assert (m, n) == size(students) "Shape mismatch between schools and students"

    # Open-source, slow for this problem.
    # model = Model(GLPK.Optimizer)
    # set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)

    # Open source, fast for this problem for no reason.
    model = Model(Cbc.Optimizer)

    # Commercial option; need to have license in your path, but this problem is
    # small enough for the demo version. Fastest.
    # model = Model(Xpress.Optimizer)

    @variable(model, x[1:m, 1:n], Bin)
    @objective(model, Min, sum(x .* students))

    if sum(capacities) ≤ n
        @constraints(model, begin
            student_capacity[s in 1:n], sum(x[:, s]) ≤ 1
            school_capacity[c in 1:m], sum(x[c, :]) == capacities[c]
            stability[c in 1:m, s in 1:n], (capacities[c] * x[c, s] +
                                            capacities[c] * x[:, s]' * (students[:, s] .≤ students[c, s]) +
                                            x[c, :]' * (schools[:, c] .≤ schools[s, c])) ≥ capacities[c]
        end)
    else
        @constraints(model, begin
            student_capacity[s in 1:n], sum(x[:, s]) == 1
            school_capacity[c in 1:m], sum(x[c, :]) ≤ capacities[c]
            stability[c in 1:m, s in 1:n], (capacities[c] * x[c, s] +
                                            capacities[c] * x[:, s]' * (students[:, s] .≤ students[c, s]) +
                                            x[c, :]' * (schools[:, c] .≤ schools[s, c])) ≥ capacities[c]
        end)
    end

    optimize!(model)

    display("Total student disutility: $(objective_value(model))")
    out = value.(x)
    assn = vec(mapslices(col -> findfirst(y -> y > 0.9, col), out, dims = 1))
    return assn
end

"""
Same but without stability constraint.
"""
function LP_system_opt(students::Array{Int64, 2}, schools::Array{Int64, 2},
                    capacities::Array{Int64, 1})
    n, m = size(schools)
    @assert (m,) == size(capacities)
    @assert (m, n) == size(students) "Shape mismatch between schools and students"

    # Open source. Works great for LP case.
    model = Model(GLPK.Optimizer)
    set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)

    # Open source.
    # model = Model(Cbc.Optimizer)

    # Commercial option.
    # model = Model(Xpress.Optimizer)

    # Constraint matrix is TUM so no integrality constraint.
    @variable(model, 0 ≤ x[1:m, 1:n] ≤ 1)
    @objective(model, Min, sum(x .* students))
    if sum(capacities) ≤ n
        @constraints(model, begin
            student_capacity[s in 1:n], sum(x[:, s]) ≤ 1
            school_capacity[c in 1:m], sum(x[c, :]) == capacities[c]
        end)
    else
        @constraints(model, begin
            student_capacity[s in 1:n], sum(x[:, s]) == 1
            school_capacity[c in 1:m], sum(x[c, :]) ≤ capacities[c]
        end)
    end

    optimize!(model)
    display("Total student disutility: $(objective_value(model))")
    out = value.(x)
    @assert isempty(filter(x -> !(round(x, digits=5) in [0, 1]), out)) "Solution not integral"
    assn = vec(mapslices(col -> findfirst(y -> y > 0.9, col), out, dims = 1))
    return assn
end

students = readdlm("examples/sysopt/students.dat", Int)
schools = readdlm("examples/sysopt/schools.dat", Int)
n, m = size(schools)
capacities = ones(Int64, m)

assn_LP = LP_system_opt(students, schools, capacities)
# GLPK Simplex Optimizer, v4.64
# 200 rows, 10000 columns, 20000 non-zeros
#       0: obj =  0.000000000e+000 inf =  1.000e+002 (100)
#     199: obj =  4.704000000e+003 inf =  0.000e+000 (0)
# Perturbing LP to avoid stalling [400]...
# Removing LP perturbation [1301]...
# *  1301: obj =  2.407000000e+003 inf =  0.000e+000 (0) 8
# OPTIMAL LP SOLUTION FOUND

assn = IP_stable_opt(students, schools, capacities)
# Result - Optimal solution found
#
# Objective value:                2426.00000000
# Enumerated nodes:               0
# Total iterations:               0
# Time (CPU seconds):             15.34
# Time (Wallclock seconds):       15.34
#
# Total time (CPU seconds):       15.38   (Wallclock seconds):       15.38
#
# "Total student disutility: 2426.0"

# writedlm("examples/sysopt/LPSystemOpt.dat", assn_LP)
# writedlm("examples/sysopt/IPStableOpt.dat", assn)

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
module CachedViews

using OffsetArrays

abstract type CachedView{T, N} <: AbstractArray{T, N} end
abstract type AbstractCacheStrategy end

export NullCache, LocalWindowCache

include("NullCache.jl")
include("LocalWindowCache.jl")

end # module


text:
module GroupsCore

import Random
import Markdown
import AbstractAlgebra
import AbstractAlgebra: elem_type, gens, ngens, order, parent_type
import AbstractAlgebra: inv!, mul!
const Group = AbstractAlgebra.Group
const GroupElement = AbstractAlgebra.GroupElem

# abstract type Group end
# abstract type GroupElement end

export Group, GroupElement
export comm, istrulyequal, gens, hasorder, hasgens, ngens, order
# export one!, inv!, mul!, conj!, comm!, div_left!, div_right!

struct InterfaceNotImplemented <: Exception
    family::Symbol
    method::String
end

Base.showerror(io::IO, err::InterfaceNotImplemented) =
    print(io, "Missing method from $(err.family) interface: `$(err.method)`")

include("groups.jl")
include("group_elements.jl")

end


text:
mutable struct Philosopher
    name::String
    hungry::Bool
    righthanded::Bool
    rightforkheld::Channel
    leftforkheld::Channel
    function Philosopher(name, leftfork, rightfork)
        this = new()
        this.name = name
        this.hungry = rand([false, true]) # not specified so start as either
        this.righthanded   = (name == "Aristotle") ? false : true
        this.leftforkheld  = leftfork
        this.rightforkheld = rightfork
        this
    end
end

mutable struct FiveForkTable
    fork51::Channel
    fork12::Channel
    fork23::Channel
    fork34::Channel
    fork45::Channel
    function FiveForkTable()
        this = new()
        this.fork51 = Channel(1); put!(this.fork51, "fork") # start with one fork per channel
        this.fork12 = Channel(1); put!(this.fork12, "fork")
        this.fork23 = Channel(1); put!(this.fork23, "fork")
        this.fork34 = Channel(1); put!(this.fork34, "fork")
        this.fork45 = Channel(1); put!(this.fork45, "fork")
        this
    end
end

table = FiveForkTable();
tasks = [Philosopher("Aristotle", table.fork12, table.fork51),
         Philosopher("Kant", table.fork23, table.fork12),
         Philosopher("Spinoza", table.fork34, table.fork23),
         Philosopher("Marx", table.fork45, table.fork34),
         Philosopher("Russell", table.fork51, table.fork45)]

function dine(t,p)
    if p.righthanded
       take!(p.rightforkheld); println("$(p.name) takes right fork")
       take!(p.leftforkheld); println("$(p.name) takes left fork")
    else
       take!(p.leftforkheld); println("$(p.name) takes left fork")
       take!(p.rightforkheld); println("$(p.name) takes right fork")
    end
end

function leavetothink(t, p)
    put!(p.rightforkheld, "fork"); println("$(p.name) puts down right fork")
    put!(p.leftforkheld, "fork");  println("$(p.name) puts down left fork")
end

contemplate(t) = sleep(t)

function dophil(p, t, fullaftersecs=2.0, hungryaftersecs=10.0)
    while true
        if p.hungry
            println("$(p.name) is hungry")
            dine(table, p)
            sleep(fullaftersecs)
            p.hungry = false
            leavetothink(t, p)
        else
            println("$(p.name) is out of the dining room for now.")
            contemplate(hungryaftersecs)
            p.hungry = true
        end
    end
end

function runall(tasklist)
    for p in tasklist
        @async dophil(p, table)
    end
    while true begin sleep(5) end end
end

runall(tasks)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function ViennaModelCoefficientsLoader(arg0::JString, arg1::jdouble, arg2::jdouble, arg3::ViennaModelType)
    return ViennaModelCoefficientsLoader((JString, jdouble, jdouble, ViennaModelType), arg0, arg1, arg2, arg3)
end

function ViennaModelCoefficientsLoader(arg0::JString, arg1::jdouble, arg2::jdouble, arg3::ViennaModelType, arg4::DataProvidersManager)
    return ViennaModelCoefficientsLoader((JString, jdouble, jdouble, ViennaModelType, DataProvidersManager), arg0, arg1, arg2, arg3, arg4)
end

function ViennaModelCoefficientsLoader(arg0::jdouble, arg1::jdouble, arg2::ViennaModelType)
    return ViennaModelCoefficientsLoader((jdouble, jdouble, ViennaModelType), arg0, arg1, arg2)
end

function get_a(obj::ViennaModelCoefficientsLoader)
    return jcall(obj, "getA", Vector{jdouble}, ())
end

function get_supported_names(obj::ViennaModelCoefficientsLoader)
    return jcall(obj, "getSupportedNames", JString, ())
end

function get_zenith_delay(obj::ViennaModelCoefficientsLoader)
    return jcall(obj, "getZenithDelay", Vector{jdouble}, ())
end

function load_data(obj::ViennaModelCoefficientsLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function load_vienna_coefficients(obj::ViennaModelCoefficientsLoader)
    return jcall(obj, "loadViennaCoefficients", void, ())
end

function load_vienna_coefficients(obj::ViennaModelCoefficientsLoader, arg0::DateTimeComponents)
    return jcall(obj, "loadViennaCoefficients", void, (DateTimeComponents,), arg0)
end

function still_accepts_data(obj::ViennaModelCoefficientsLoader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
using Test
using SumOfSquares
using DynamicPolynomials

function motzkin_test(optimizer, config::MOIT.Config)
    @polyvar x y
    # Motzkin polynomial
    p = x^4*y^2 + x^2*y^4 + 1 - 3*x^2*y^2

    model = _model(optimizer)
    # We want to write `p ≥ 0` instead of `p in SOSCone()` so we need to set
    # the polymodule to `SumOfSquares` to tell `PolyJuMP` to interpret
    # `p ≥ 0` as a Sum-of-Squares constraint.
    PolyJuMP.setpolymodule!(model, SumOfSquares)
    con_ref = @constraint(model, p ≥ 0)

    JuMP.optimize!(model)
    @test JuMP.termination_status(model) == MOI.INFEASIBLE
    @test JuMP.dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE

    q = (x^2 + y^2) * p
    delete(model, con_ref)
    @constraint(model, q >= 0)

    JuMP.optimize!(model)
    @test JuMP.termination_status(model) == config.optimal_status
    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
end
sd_tests["motzkin"] = motzkin_test

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

using FEMQuad, Test
using FEMQuad: integrate_3d, get_rule

@testset "Gauss-Legendre quadratures in pyramids" begin
    # FIXME: fix tests
    f(i) = x -> sum(sum(x)^j for j=0:i)
    @test isapprox(integrate_3d(f(0), :GLPYR5B), 4.0/6.0)
    #@test isapprox(integrate_3d(f(1), :GLPYR5B), 7.0/6.0)
end


text:
##
# Demo of using DualNumbers to differentiate with respect to
# changes in parameters
##

using DualNumbers, ApproxFun
# What is the derivative of the function (without differentiating the Chebyshev expansion)?
f=Fun(x->exp(dual(x,1)),-1..1)
dualpart(f)
# check versus differentiate
norm(realpart(f)'-dualpart(f))

# What is the derivative of the first coefficient with respect to the left endpoint?
f=Fun(exp,Segment(dual(1.0,1),dual(2.0)))
dualpart(f.coefficients[1])
# check versus finite difference calculation:
h=0.00001;(Fun(exp,Segment(1.0+h,2.0)).coefficients[1]-Fun(exp,Segment(1.0,2.0)).coefficients[1])/h

# Or an ApproxFun calculation:
Fun(h->Fun(exp,Segment(1.0+h,2.0)).coefficients[1],0..1)'(0.)

# What is the derivative of the first coefficient with respect to the exponential's constant?
f=Fun(x->exp(dual(x,x)),-1..1)
dualpart(f.coefficients[1])
# check versus finite difference calculation:
h=0.00001;(Fun(x->exp((1+h)x),-1..1).coefficients[1]-Fun(exp,-1..1).coefficients[1])/h

max_stars_repo_path: Extras/Dual Numbers.jl

max_issues_repo_path: Extras/Dual Numbers.jl

max_forks_repo_path: Extras/Dual Numbers.jl


text:
"""
Topology enumeration for Ground Structure with Discrete Diameters.

Solver object to store parameter values.
"""
struct IterTopo <: GroundStructureSolver
    maxiter::Int
    timelimit::Float64 # seconds
    mastersolver
    subsolver
    writemodels::Bool

    function IterTopo(mastersolver, subsolver;
                      maxiter::Int=100, timelimit=Inf, writemodels=false)
        new(maxiter, timelimit, mastersolver, subsolver, writemodels)
    end
end

"""
Build master model using only arc selection y and flows q.

This is used to enumerate (embedded) topologies, each of which is evaluated with
a "semi-subproblem". Returns model and variables y, q.
"""
function make_semimaster(inst::Instance, topo::Topology, optimizer)
    nodes, nnodes = topo.nodes, length(topo.nodes)
    arcs, narcs = topo.arcs, length(topo.arcs)
    terms, nterms = inst.nodes, length(inst.nodes)
    termidx = termindex(nodes, terms)

    # demand for all nodes, including junctions
    dem = fill(0.0, nnodes)
    for t in 1:nterms
        dem[termidx[t]] = inst.demand[t]
    end

    # adjacency lists
    inarcs, outarcs = [Int[] for v in 1:nnodes], [Int[] for v in 1:nnodes]
    for a in 1:narcs
        tail, head = arcs[a]
        push!(inarcs[head], a)
        push!(outarcs[tail], a)
    end

    # "big-M" bound for flow on arcs
    maxflow = 0.5 * sum(abs.(inst.demand))

    # always use direct mode with SCIP
    model = JuMP.direct_model(MOI.instantiate(optimizer))

    # select arcs from topology with y
    @variable(model, y[1:narcs], Bin)

    # flow through arcs
    @variable(model, 0 <= q[1:narcs] <= maxflow)

    # mass flow balance at nodes
    @constraint(model, flowbalance[v=1:nnodes],
                sum(q[a] for a=inarcs[v]) - sum(q[a] for a=outarcs[v]) == dem[v])

    # allow flow only for active arcs
    @constraint(model, active[a=1:narcs], q[a] <= maxflow*y[a])

    # exclude antiparallel arcs.
    antiidx = antiparallelindex(topo)
    for a in 1:narcs
        if a < antiidx[a] # only at most once per pair
            @constraint(model, y[a] + y[antiidx[a]] ≤ 1)
        end
    end

    # use cost of smallest diameter as topology-based relaxation
    L = pipelengths(topo)
    cmin = inst.diameters[1].cost
    @objective(model, Min, sum(cmin * L[a] * y[a] for a=1:narcs))

    model, y, q
end

"""
Build model for a problem where the topology y and flow q are fixed, but the
diameters z are still free. This is between the master- and subproblem and can
be used
 - as a primal heuristic
 - to generate stronger no-good cuts on the y variables

Returns model, list of candidate arcs and (sparse) variables z
"""
function make_semisub(inst::Instance, topo::Topology, cand::CandSol, optimizer)
    nnodes = length(topo.nodes)
    arcs = topo.arcs
    termidx = termindex(topo.nodes, inst.nodes)
    ndiams = length(inst.diameters)
    candarcs = filter(a -> any(cand.zsol[a,:]), 1:length(arcs))
    ncandarcs = length(candarcs)
    tail = [arcs[a].tail for a in candarcs]
    head = [arcs[a].head for a in candarcs]

    πlb_min = minimum([b.lb for b in inst.pressure])^2
    πub_max = maximum([b.ub for b in inst.pressure])^2
    πlb = fill(πlb_min, nnodes)
    πub = fill(πub_max, nnodes)
    πlb[termidx] = [b.lb^2 for b in inst.pressure]
    πub[termidx] = [b.ub^2 for b in inst.pressure]

    L = pipelengths(topo)
    c = [diam.cost for diam in inst.diameters]
    Dm5 = [diam.value^(-5) for diam in inst.diameters]
    C = inst.ploss_coeff * L[candarcs] .* cand.qsol[candarcs].^2

    model = JuMP.Model(optimizer)

    @variable(model, z[1:ncandarcs, 1:ndiams], Bin)
    @variable(model, πlb[v] ≤ π[v=1:nnodes] ≤ πub[v])

    @constraint(model, ploss[a=1:ncandarcs],
                π[tail[a]] - π[head[a]] ≥ C[a] * sum(Dm5[i]*z[a,i] for i=1:ndiams))
    @constraint(model, choice[a=1:ncandarcs], sum(z[a,i] for i=1:ndiams) == 1)

    @objective(model, Min, sum(c[i] * L[candarcs[a]] * z[a,i] for a=1:ncandarcs for i=1:ndiams))

    model, candarcs, z
end

"Iteration based decomposition with semimaster and ~subproblem."
function PipeLayout.optimize(inst::Instance, topo::Topology, solver::IterTopo)
    run_semi(inst, topo, solver.mastersolver, solver.subsolver,
             maxiter=solver.maxiter, timelimit=solver.timelimit,
             writemodels=solver.writemodels)
end

function run_semi(inst::Instance, topo::Topology, mastersolver, subsolver;
                  maxiter::Int=100, timelimit=timelimit, writemodels=false)
    finaltime = time() + timelimit
    narcs = length(topo.arcs)
    ndiams = length(inst.diameters)

    # initialize
    primal, dual, status, bestsol = Inf, 0.0, MOI.OPTIMIZE_NOT_CALLED, nothing
    mastermodel, y, q = make_semimaster(inst, topo, mastersolver)

    for iter=1:maxiter
        if !stilltime(finaltime)
            status = MOI.TIME_LIMIT
            @debug "Timelimit reached."
            return Result(status, bestsol, primal, dual, iter)
        end
        @debug "Iter $(iter)"

        # resolve (relaxed) semimaster problem, build candidate solution
        writemodels && writeLP(mastermodel, "master_iter$(iter).lp", genericnames=false)
        settimelimit!(mastermodel, mastersolver, finaltime - time())
        JuMP.optimize!(mastermodel)
        status = JuMP.termination_status(mastermodel)
        if status == MOI.INFEASIBLE
            @debug "  master problem is infeasible."
            if primal == Inf
                # no solution was found
                return Result(MOI.INFEASIBLE, nothing, Inf, Inf, iter)
            else
                @assert bestsol ≠ nothing
                return Result(MOI.OPTIMAL, bestsol, primal, dual, iter)
            end
        elseif status in (MOI.NODE_LIMIT, MOI.TIME_LIMIT)
            return Result(status, bestsol, primal, dual, iter)
        elseif status == MOI.OPTIMAL
            # good, we continue below
        else
            error("Unexpected status: $(status)")
        end

        ysol, qsol = JuMP.value.(y), JuMP.value.(q)
        dual = JuMP.objective_value(mastermodel)
        @debug begin
            "  dual bound: $(dual)"
            "  cand. sol:$(findall(!iszero, qsol))"
        end

        # stopping criterion
        if dual > primal - ɛ
            @assert bestsol ≠ nothing
            @debug "  proved optimality of best solution."
            return Result(MOI.OPTIMAL, bestsol, primal, dual, iter)
        end

        # check whether candidate has tree topology
        candtopo = topology_from_candsol(topo, ysol)
        if !is_tree(candtopo)
            fullcandtopo = topology_from_candsol(topo, ysol, true)
            cycle = find_cycle(fullcandtopo)
            if length(cycle) == 0
                # TODO: see above (candidate disconnected?)
                nogood(mastermodel, y, ysol)
                @debug "  skip disconnected topology with nogood."
            else
                avoid_topo_cut(mastermodel, y, topo, cycle)
                @debug "  skip non-tree topology, cycle: $(cycle)"
            end
            continue
        end

        # solve subproblem (from scratch, no warmstart)
        zcand = fill(false, narcs, ndiams)
        zcand[ysol .> 0.5, 1] .= true
        cand = CandSol(zcand, qsol, qsol.^2)

        submodel, candarcs, z = make_semisub(inst, topo, cand, subsolver)
        writemodels && writeLP(submodel, "sub_iter$(iter).lp", genericnames=false)
        settimelimit!(submodel, subsolver, finaltime - time())
        JuMP.optimize!(submodel)
        substatus = JuMP.termination_status(submodel)
        if substatus == MOI.OPTIMAL
            # have found improving solution?
            newobj = JuMP.objective_value(submodel)
            if newobj < primal
                primal = newobj
                znew = fill(false, narcs, ndiams)
                znew[candarcs,:] = (JuMP.value.(z) .> 0.5)
                bestsol = CandSol(znew, qsol, qsol.^2)
                @debug "  found improving solution: $(primal)"
            end
        elseif substatus != MOI.INFEASIBLE
            error("Unexpected status: $(:substatus)")
        end

        # stopping criterion
        if dual > primal - ɛ
            @assert bestsol ≠ nothing
            @debug "  proved optimality of best solution."
            return Result(MOI.OPTIMAL, bestsol, primal, dual, iter)
        end

        # generate nogood cut and add to master
        nogood(mastermodel, y, ysol)
    end

    Result(MOI.ITERATION_LIMIT, bestsol, primal, dual, maxiter)
end


text:
"""
    sqrt_kalman_filter(model::StateSpaceModel{Typ}; tol::Typ = Typ(1e-5)) where Typ

Square-root Kalman filter with big Kappa initialization.
"""
function sqrt_kalman_filter(model::StateSpaceModel{Typ}; tol::Typ = Typ(1e-5)) where Typ

    time_invariant = model.mode == "time-invariant"

    # Load dimensions
    n, p, m, r = size(model)

    # Load system
    y = model.y
    Z, T, R = ztr(model)

    ensure_pos_sym!(model.H)
    ensure_pos_sym!(model.Q)

    sqrtH = cholesky(model.H).L
    sqrtQ = cholesky(model.Q).L

    # Predictive state and its sqrt-covariance
    a     = Matrix{Typ}(undef, n+1, m)
    sqrtP = Array{Typ, 3}(undef, m, m, n+1)

    # Innovation and its sqrt-covariance
    v     = Matrix{Typ}(undef, n, p)
    sqrtF = Array{Typ, 3}(undef, p, p, n)

    # Kalman gain
    K     = Array{Typ, 3}(undef, m, p, n)

    # Auxiliary matrices
    U2star = Array{Typ, 3}(undef, m, p, n)

    # Steady state initialization
    steadystate  = false
    tsteady      = n+1

    # Initial state: big Kappa initialization
    a[1, :]        = zeros(m, 1)
    sqrtP[:, :, 1]    = 1e3 .* Matrix(I, m, m)

    # Pre-allocating for performance
    zeros_pr = zeros(p, r)
    zeros_mp = zeros(m, p)
    range1   = (p + 1):(p + m)
    range2   = 1:p
    sqrtH_zeros_pr  = [sqrtH zeros_pr]
    zeros_mp_RsqrtQ = [zeros_mp R*sqrtQ]

    !isempty(model.missing_observations) && error("Treatment of missing values not implemented for SquareRootFilter, please use KalmanFilter.")

    # Square-root Kalman filter
    for t = 1:n
        v[t, :] = y[t, :] - Z[:, :, t]*a[t, :] - model.d[t, :]
        if steadystate
            sqrtF[:, :, t] = sqrtF[:, :, t-1]
            K[:, :, t] = K[:, :, t-1]
            a[t+1, :] = T * a[t, :] + K[:, :, t] * v[t, :]
            sqrtP[:, :, t+1] = sqrtP[:, :, t]
        else
            # Manipulation of auxiliary matrices
            U         = SMatrix{p + m, m + p + r}([Z[:, :, t]*sqrtP[:, :, t] sqrtH_zeros_pr;
                                                   T*sqrtP[:, :, t]         zeros_mp_RsqrtQ])
            G         = qr(U').Q
            Ustar     = U*G
            U2star[:, :, t] = Ustar[range1, range2]
            sqrtF[:, :, t]  = Ustar[range2, range2]

            # Kalman gain and predictive state update
            K[:, :, t]       = U2star[:, :, t]*inv(sqrtF[:, :, t])
            a[t+1, :]        = T*a[t, :] + K[:, :, t]*v[t, :] + model.c[t, :]
            sqrtP[:, :, t+1] = Ustar[range1, range1]

            # Checking if steady state was attained
            if time_invariant && check_steady_state(sqrtP, t, tol)
                steadystate = true
                tsteady     = t
            end
        end
    end

    # Return the auxiliary filter structre
    return SquareRootFilter(a, v, sqrtP, sqrtF, steadystate, tsteady, K)
end

"""
    filtered_state(model::StateSpaceModel{Typ}, sqrt_filter::SquareRootFilter{Typ}) where Typ

Obtain the filtered state estimates and their covariance matrices.
"""
function filtered_state(model::StateSpaceModel{Typ}, sqrt_filter::SquareRootFilter{Typ}) where Typ

    # Load dimensions data
    n, p, m, r = size(model)

    # Load system data
    Z, T, R = ztr(model)

    # Load filter data
    a = sqrt_filter.a
    v = sqrt_filter.v
    F = gram_in_time(sqrt_filter.sqrtF)
    P = gram_in_time(sqrt_filter.sqrtP)

    # Filtered state and its covariance
    att = Matrix{Typ}(undef, n, m)
    Ptt = Array{Typ, 3}(undef, m, m, n)

    for t = 1:n
        PZF = P[:, :, t] * Z[:, :, t]' * inv(F[:, :, t])
        att[t, :]    = a[t, :] + PZF * v[t, :] + model.c[t, :]
        Ptt[:, :, t] = P[:, :, t] - PZF * Z[:, :, t] * P[:, :, t]
        ensure_pos_sym!(Ptt, t)
    end

    return att, Ptt

end

"""
    sqrt_smoother(model::StateSpaceModel, sqrt_filter::SquareRootFilter) where Typ

Square-root smoother for state space model.
"""
function sqrt_smoother(model::StateSpaceModel{Typ}, sqrt_filter::SquareRootFilter{Typ}) where Typ

    # Load dimensions data
    n, p, m, r = size(model)

    # Load system data
    Z, T, R = ztr(model)

    # Load filter data
    a           = sqrt_filter.a
    v           = sqrt_filter.v
    tsteady     = sqrt_filter.tsteady
    sqrtF       = sqrt_filter.sqrtF
    sqrtP       = sqrt_filter.sqrtP
    K           = sqrt_filter.K

    # Smoothed state and its covariance
    alpha = Matrix{Typ}(undef, n, m)
    V     = Array{Typ, 3}(undef, m, m, n)
    L     = Array{Typ, 3}(undef, m, m, n)
    r     = Matrix{Typ}(undef, n, m)
    sqrtN = Array{Typ, 3}(undef, m, m, n)

    # Initialization
    sqrtN[:, :, end]  = zeros(m, m)
    r[end, :]      = zeros(m, 1)
    sqrtPsteady = sqrtP[:, :, end]
    sqrtFsteady = sqrtF[:, :, end]

    # Iterating backwards
    for t = n:-1:tsteady
        Psteady = gram(sqrtP[:, :, end])
        Fsteady = gram(sqrtF[:, :, end])
        N_t_1   = gram(sqrtN[:, :, t-1])
        L[:, :, t]   = T - K[:, :, end] * Z[:, :, t]
        r[t-1, :] = Z[:, :, t]' * inv(Fsteady) * v[t, :] + L[:, :, t]' * r[t, :]

        # QR decomposition of auxiliary matrix Nstar
        Nstar        = [Z[:, :, t]' * inv(sqrtF[:, :, end]) L[:, :, t]' * sqrtN[:, :, t]]
        G            = qr(Nstar').Q
        NstarG       = Nstar * G
        sqrtN[:, :, t-1]   = NstarG[1:m, 1:m]

        # Smoothed state and its covariance
        alpha[t, :] = a[t, :] + Psteady * r[t-1, :]
        V[:, :, t]  = Psteady - Psteady * N_t_1 * Psteady
    end

    for t = tsteady-1:-1:2
        L[:, :, t]   = T - K[:, :, t] * Z[:, :, t]
        r[t-1, :] = Z[:, :, t]' * inv(sqrtF[:, :, t] * sqrtF[:, :, t]') * v[t, :] + L[:, :, t]'*r[t, :]
        Nstar  = [Z[:, :, t]' * inv(sqrtF[:, :, t]) L[:, :, t]' * sqrtN[:, :, t]]

        # QR decomposition of auxiliary matrix Nstar
        NstarG     = Nstar*qr(Nstar').Q
        sqrtN[:, :, t-1] = NstarG[1:m, 1:m]

        # Smoothed state and its covariance
        P = gram(sqrtP[:, :, t])
        N = gram(sqrtN[:, :, t-1])
        alpha[t, :] = a[t, :] + P * r[t-1, :]
        V[:, :, t]  = P - (P * N * P)
    end

    L[:, :, 1]   = T - K[:, :, 1] * Z[:, :, 1]
    r_0    = Z[:, :, 1]' * inv(gram(sqrtF[:, :, 1])) * v[1, :] + L[:, :, 1]' * r[1, :]
    Nstar  = [Z[:, :, 1]' * inv(sqrtF[:, :, 1]) L[:, :, 1]' * sqrtN[:, :, 1]]
    G      = qr(Nstar').Q
    NstarG = Nstar*G

    sqrtN_0  = NstarG[1:m, 1:m]
    P_1 = gram(sqrtP[:, :, 1])
    alpha[1, :] = a[1, :] + P_1 * r_0
    V[:, :, 1]  = P_1 - (P_1 * gram(sqrtN_0) * P_1)

    # Return the Square Root kalman filter smoothed state
    return Smoother(alpha, V)
end

function get_log_likelihood_params(model::StateSpaceModel, filter_type::Type{SquareRootFilter{T}}) where T

    # Obtain innovation v and its variance F
    sqrt_filter = sqrt_kalman_filter(model)
    # Return v and F
    return sqrt_filter.v, gram_in_time(sqrt_filter.sqrtF)
end

function kfas(model::StateSpaceModel{T}, filter_type::Type{SquareRootFilter{T}}) where T
    # Do the SquareRootFilter
    filter_output  = sqrt_kalman_filter(model)
    smoothed_state = sqrt_smoother(model, filter_output)
    att, Ptt       = filtered_state(model, filter_output)

    return FilterOutput(filter_output.a, att, filter_output.v,
                         gram_in_time(filter_output.sqrtP), Ptt, gram_in_time(filter_output.sqrtF),
                         filter_output.steadystate, filter_output.tsteady) ,
           SmoothedState(smoothed_state.alpha, smoothed_state.V)
end


text:
module EDModels

using SparseArrays
#using ...Positions
using ...SimLib
using SpinSymmetry
using Printf: @sprintf
using XXZNumerics

export parameter_dims, initialize_model!, split_workload, do_parameters, model_name, model_fileprefix
export RandomPositionsXXZWithXField, RandomPositionsXXZWithDegeneracyLifted

abstract type Model end

"""
    parameter_dims(model) -> (dim1, dim2, ...)
"""
function parameter_dims end

"""
    initialize_model!(model, array_initializer)
"""
function initialize_model! end

"""
    split_workload(model, num_parts) -> [parameter_chunk...]

`parameter_chunk` is put through `do_parameters` to iterate all Hamiltonians in this section of the parameter space.
"""
function split_workload end#(model, num_parts(runmode))

"""
    do_parameters(f, model, parameter_chunk)

f = f(parameter_index, H)

Construct all Hamiltonians associated with the parameter_chunk and call `f` with the index and the Matrix.
"""
function do_parameters end

"""
    model_name(model)

Return a nice name that could be used as foldername :)
"""
function model_name end

"""
    model_fileprefix(model)

Return a nice name that could prefix a file with data derived from this model.
"""
function model_fileprefix end

include("Jxxz_x.jl")
include("Jxxz_Jz.jl")
#include("xxz_hx.jl")

end #module


text:
function vctest(y, X, V; bInit::Array{Float64, 1} = Float64[],
                devices::String = "CPU", nMMmax::Int = 0,
                nBlockAscent::Int = 1000, nNullSimPts::Int = 1000000,
                nNullSimNewtonIter::Int = 15, tests::String = "eLRT",
                tolX::Float64 = 1e-4,
                vcInit::Array{Float64, 1} = Float64[],
                Vform::String = "whole",
                pvalueComputings::String = "chi2",
                WPreSim::Array{Float64, 2} = [Float64[] Float64[]],
                PrePartialSumW::Array{Float64, 2} = [Float64[] Float64[]],
                PreTotalSumW::Array{Float64, 2} = [Float64[] Float64[]],
                partialSumWConst::Array{Float64, 1} = Float64[],
                totalSumWConst::Array{Float64, 1} = Float64[],
                windowSize::Int = 50,
                partialSumW::Array{Float64, 1} = Float64[],
                totalSumW::Array{Float64, 1} = Float64[],
                lambda::Array{Float64, 2} = [Float64[] Float64[]],
                W::Array{Float64, 2} = [Float64[] Float64[]],
                nPreRank::Int = 20,
                tmpmat0::Array{Float64, 2} = [Float64[] Float64[]],
                tmpmat1::Array{Float64, 2} = [Float64[] Float64[]],
                tmpmat2::Array{Float64, 2} = [Float64[] Float64[]],
                tmpmat3::Array{Float64, 2} = [Float64[] Float64[]],
                tmpmat4::Array{Float64, 2} = [Float64[] Float64[]],
                tmpmat5::Array{Float64, 2} = [Float64[] Float64[]],
                denomvec::Array{Float64, 1} = Float64[],
                d1f::Array{Float64, 1} = Float64[],
                d2f::Array{Float64, 1} = Float64[], offset::Int = 0,
                nPtsChi2::Int = 300,
                simnull::Vector{Float64} = Float64[])
  # VCTEST Fit and test for the nontrivial variance component
  #
  # [SIMNULL] = VCTEST(y,X,V) fits and then tests for $H_0:sigma_1^2=0$ in
  # the variance components model $Y \sim N(X \beta, vc0*I + vc1*V)$.
  #
  #   INPUT:
  #       y - response vector
  #       X - design matrix for fixed effects
  #       V - the variance component being tested
  #
  #   Optional input name-value pairs:
  #       'nBlockAscent' - max block ascent iterations, default is 1000
  #       'nMMmax' - max MM iterations, default is 10 (eLRT) or 1000 (eRLRT)
  #       'nSimPts'- # simulation samples, default is 10000
  #       'test'- 'eLRT'|'eRLRT'|'eScore'|'none', the requested test, it also
  #           dictates the estimation method
  #       'tolX' - tolerance in change of parameters, default is 1e-5
  #       'Vform' - format of intpu V. 'whole': V, 'half': V*V', 'eigen':
  #           V.U*diag(V.E)*V.U'. Default is 'whole'. For computational
  #           efficiency,  a low rank 'half' should be used whenever possible
  #
  #   Output:
  #       b - estimated regression coeffiicents in the mean component
  #       vc0 - estiamted variance component for I
  #       vc1 - estiamted variance component for V
  #       stats - other statistics

  # check dimensionalities
  n = length(y);
  if size(V, 1) != n
    error("vctest:wrongdimV\n", "dimension of V does not match that of X");
  end

  # set default maximum MM iteration
  if nMMmax == 0 && tests == "eLRT"
    nMMmax = 10;
  elseif nMMmax == 0 && tests == "eRLRT"
    nMMmax = 1000;
  end

  # SVD of X
  if isempty(X)
    rankX = 0;
    X = reshape(X, n, 0);
    X = convert(Array{Float64, 2}, X);
    # LRT is same as RLRT if X is null
    if tests == "eLRT"
      tests = "eRLRT";
    end
  else
    if tests == "eRLRT"
      (UX, svalX) = svd(X, thin = false);
    else
      (UX, svalX) = svd(X);
    end
    rankX = countnz(svalX .> n * eps(svalX[1]));
  end

  # eigendecomposition of V
  if Vform == "whole"
    (evalV, UV) = eig(V);
    rankV = countnz(evalV .> n * eps(sqrt.(maximum(evalV))));
    sortIdx = sortperm(evalV, rev = true);
    evalV = evalV[sortIdx[1:rankV]];
    UV = UV[:, sortIdx[1:rankV]];
    if tests == "eLRT" || tests == "eRLRT"
      wholeV = V;
    end
  elseif Vform == "half"
    (UVfull, tmpevalVfull) = svd(V, thin=false);
    evalVfull = zeros(n);
    pevalVfull = pointer(evalVfull);
    ptmpevalV = pointer(tmpevalVfull);
    BLAS.blascopy!(length(tmpevalVfull), ptmpevalV, 1, pevalVfull, 1);
    rankV = countnz(evalVfull .> n * eps(maximum(evalVfull)));
    evalVfull = evalVfull .^ 2;
    evalV = evalVfull[1:rankV];
    UV = UVfull[:, 1:rankV];
    #if tests == "eLRT" || tests == "eRLRT"
    #  wholeV = *(V, V');
    #end
  elseif Vform == "eigen"
    UV = V.U;
    evalV = V.eval;
    rankV = countnz(V.eval .> n * eps(sqrt.(maximum(V.eval))));
    sortIdx = sortperm(V.eval, rev = true);
    V.eval = V.eval[sortIdx[1:rankV]];
    V.U = V.U[:, sortIdx[1:rankV]];
    if tests == "eLRT" || tests == "eRLRT"
      wholeV = *(V.U .* reshape(V.eval, length(V.eval), 1), V.U');
    end
  end

  # obtain eigenvalues of (I-PX)V(I-PX)
  if !isempty(X) || tests == "eScore"
    #sqrtV = UV .* sqrt(evalV)';
    sqrtV = similar(UV);
    psqrtV = pointer(sqrtV);
    pUV = pointer(UV);
    BLAS.blascopy!(n*rankV, pUV, 1, psqrtV, 1);
    scale!(sqrtV, sqrt.(evalV));
  end
  if isempty(X)
    evalAdjV = evalV;
  else
    subUX = Array(Float64, n, rankX);
    psubUX = pointer(subUX);
    pUX = pointer(UX);
    BLAS.blascopy!(n*rankX, pUX, 1, psubUX, 1);
    mat1 = BLAS.gemm('T', 'N', 1.0, subUX, sqrtV);
    mat2 = BLAS.gemm('N', 'N', 1.0, subUX, mat1);
    (UAdjV, evalAdjV) = svd(sqrtV - mat2, thin = false);
    if isempty(evalAdjV)
      evalAdjV = Float64[];
    else
      evalAdjV = evalAdjV[evalAdjV .> n * eps(maximum(evalAdjV))] .^ 2;
    end
  end
  rankAdjV = length(evalAdjV);

  # fit the variance component model
  if tests == "eLRT"

    # estimates under null model
    bNull = X \ y;
    rNull = y - X * bNull;
    vc0Null = norm(rNull) ^ 2 / n;
    Xrot = UVfull' * X;
    yrot = UVfull' * y;
    loglConst = - 0.5 * n * log.(2.0 * pi);

    # set starting point
    if isempty(bInit)
      b = copy(bNull);
      r = copy(rNull);
    else
      b = copy(bInit);
      r = y - X * b;
    end
    if isempty(vcInit)
      vc0 = norm(r) ^ 2 / n;
      vc1 = 1;
      wt = 1.0 ./ sqrt.(vc1 * evalVfull + vc0);
    else
      vc0 = vcInit[1];
      # avoid sticky boundary
      if vc0 == 0
        vc0 = 1e-4;
      end
      vc1 = vcInit[2];
      if vc1 == 0
        vc1 = 1e-4;
      end
      wt = 1.0 ./ sqrt.(vc1 * evalVfull + vc0);
      Xnew = scale(wt, Xrot);
      ynew = wt .* yrot;
      b = Xnew \ ynew;
    end

    # update residuals according supplied var. components
    r = y - BLAS.gemv('N', X, b);
    rnew = BLAS.gemv('T', UVfull, r);
    #deltaRes = norm(r) ^ 2 - norm(rnew) ^ 2;
    logl = loglConst + sum(log, wt) - 0.5 * sumabs2(rnew .* wt);

    nIters = 0;
    #bOld = similar(b);
    #pbOld = pointer(bOld);
    #pb = pointer(b);
    #BLAS.blascopy!(length(b), pb, 1, pbOld, 1);
    #denvec = similar(rnew);
    denvec = Float64[];
    for iBlockAscent = 1:nBlockAscent

      nIters = iBlockAscent;

      # update variance components
      for iMM = 1:nMMmax
        denvec = 1.0 ./ (vc1 * evalVfull + vc0);
        numvec = rnew .* denvec;
        #vc0 = sqrt( (vc0 ^ 2 * sumabs2(numvec) + deltaRes) /
        #            (sumabs(denvec) + (n - rankV) / vc0) );
        vc0 = vc0 * sqrt.(sumabs2(numvec) / sumabs(denvec));
        vc1 = vc1 * sqrt.( dot(numvec, numvec .* evalVfull) /
                           sumabs(evalVfull .* denvec) );
        #wt = 1.0 ./ sqrt(vc1 * evalVfull + vc0);
      end

      # update fixed effects and residuals
      #pb = pointer(b);
      #BLAS.blascopy!(length(b), pb, 1, pbOld, 1);
      wt = sqrt.(denvec);
      Xnew = scale(wt, Xrot);
      ynew = wt .* yrot;
      b = Xnew \ ynew;
      r = y - BLAS.gemv('N', X, b);
      rnew = BLAS.gemv('T', UVfull, r);
      #deltaRes = norm(r) ^ 2 - norm(rnew) ^ 2;

      # stopping criterion
      loglOld = logl;
      logl = loglConst + sum(log, wt) -
        0.5 * BLAS.dot(length(evalVfull), rnew .^ 2, 1, denvec, 1);
      if abs(logl - loglOld) < tolX * (abs(logl) + 1.0)
        break
      end

    end

    # log-likelihood at alternative
    logLikeAlt = logl;
    # log-likelihood at null
    logLikeNull = loglConst - 0.5 * n * log.(vc0Null) -
      0.5 * sum(rNull .^ 2) / vc0Null;
    if logLikeNull >= logLikeAlt
      vc0 = vc0Null;
      vc1 = 0;
      logLikeAlt = logLikeNull;
      b = bNull;
      r = rNull;
    end

    # LRT test statistic
    statLRT = 2 * (logLikeAlt - logLikeNull);

    # obtain p-value for testing vc1=0
    vc1_pvalue = vctestnullsim(statLRT, evalV, evalAdjV, n, rankX,
                               WPreSim, device = devices,
                               nSimPts = nNullSimPts,
                               nNewtonIter = nNullSimNewtonIter,
                               test = "eLRT",
                               pvalueComputing = pvalueComputings,
                               PrePartialSumW = PrePartialSumW,
                               PreTotalSumW = PreTotalSumW,
                               partialSumWConst = partialSumWConst,
                               totalSumWConst = totalSumWConst,
                               windowSize = windowSize,
                               partialSumW = partialSumW,
                               totalSumW = totalSumW,
                               lambda = lambda, W = W,
                               nPreRank = nPreRank,
                               tmpmat0 = tmpmat0, tmpmat1 = tmpmat1,
                               tmpmat2 = tmpmat2, tmpmat3 = tmpmat3,
                               tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                               denomvec = denomvec,
                               d1f = d1f, d2f = d2f, offset = offset,
                               nPtsChi2 = nPtsChi2, simnull = simnull);

    # return values
    return b, vc0, vc1, vc1_pvalue;

  elseif tests == "eRLRT"

    if isempty(X)
      ytilde = y;
      rankBVB = rankV;
      evalBVB = evalV;
      UBVB = UV;
    else
      # obtain a basis of N(X')
      B = UX[:, rankX+1:end];
      ytilde = B' * y;

      # eigen-decomposition of B'VB and transform data
      (UBVB, evalBVB) = svd(B' * sqrtV);
      rankBVB = countnz(evalBVB .> n * eps(maximum(evalBVB)));
      evalBVB = evalBVB[1:rankBVB] .^ 2;
      UBVB = UBVB[:, 1:rankBVB];
    end
    resnew = UBVB' * ytilde;
    normYtilde2 = norm(ytilde) ^ 2;
    deltaRes = normYtilde2 - norm(resnew) ^ 2;

    # set initial point
    # TODO: is there better initial point?
    vc0Null = normYtilde2 / length(ytilde);
    if isempty(vcInit)
      vc0 = vc0Null;
      vc1 = 1.0;
    else
      vc0 = vcInit[1];
      vc1 = vcInit[2];
    end

    # MM loop for estimating variance components
    nIters = 0;
    #tmpvec = Array(Float64, rankBVB);
    #denvec = Array(Float64, rankBVB);
    #numvec = Array(Float64, rankBVB);
    tmpvec = Float64[];
    for iMM = 1:nMMmax
      nIters = iMM;
      #tmpvec = vc0 + vc1 * evalBVB;
      #tmpvec = evalBVB;
      #numvecSum = 0.0;
      #denvecSum = 0.0;
      #numvecProdSum = 0.0;
      #denvecProdSum = 0.0;
      tmpvec = BLAS.scal(rankBVB, vc1, evalBVB, 1);
      #for i = 1 : rankBVB
      #  tmpvec[i] += vc0;
      #  denvec[i] = 1 / tmpvec[i];
      #  numvec[i] = (resnew[i] * denvec[i]) ^ 2;
      #  numvecSum += numvec[i];
      #  denvecSum += denvec[i];
      #  numvecProdSum += evalBVB[i] * numvec[i];
      #  denvecProdSum += evalBVB[i] * denvec[i];
      #end
      tmpvec += vc0;
      denvec = 1 ./ tmpvec;
      numvec = (resnew .* denvec) .^ 2;
      vcOld = [vc0 vc1];
      vc0 = sqrt.( (vc0 ^ 2 * sum(numvec) + deltaRes) /
                   (sum(denvec) + (n - rankX - rankBVB) / vc0) );
      #vc0 = sqrt( (vc0 ^ 2 * numvecSum + deltaRes) /
      #             (denvecSum + (n - rankX - rankBVB) / vc0) );
      vc1 = vc1 * sqrt.( sum(evalBVB .* numvec) / sum(evalBVB .* denvec) );
      #vc1 = vc1 * sqrt( numvecProdSum / denvecProdSum );
      # stopping criterion
      if norm([vc0 vc1] - vcOld) <= tolX * (norm(vcOld) + 1)
        break;
      end
    end

    # restrictive log-likelihood at alternative

    loglConst = - 0.5 * (n - rankX) * log.(2 * pi);
    logLikeAlt =  loglConst - 0.5 * sum(log.(tmpvec)) -
      0.5 * (n - rankX - rankBVB) * log.(vc0) - 0.5 * normYtilde2 / vc0 +
      0.5 * sum(resnew .^ 2 .* (1 / vc0 - 1 ./ (tmpvec)));
    # restrictive log-likelihood at null
    logLikeNull = - 0.5 * (n - rankX) * (log.(2 * pi) + log.(vc0Null)) -
      0.5 / vc0Null * normYtilde2;
    if logLikeNull >= logLikeAlt
      vc0 = vc0Null;
      vc1 = 0;
      logLikeAlt = logLikeNull;
    end

    # RLRT test statitic
    statRLRT = 2 * (logLikeAlt - logLikeNull);

    # obtain p-value for testing vc1=0
    vc1_pvalue = vctestnullsim(statRLRT, evalV, evalAdjV, n, rankX,
                               WPreSim, device = devices,
                               nSimPts = nNullSimPts,
                               nNewtonIter = nNullSimNewtonIter,
                               test = "eRLRT",
                               pvalueComputing = pvalueComputings,
                               PrePartialSumW = PrePartialSumW,
                               PreTotalSumW = PreTotalSumW,
                               partialSumWConst = partialSumWConst,
                               totalSumWConst = totalSumWConst,
                               windowSize = windowSize,
                               partialSumW = partialSumW,
                               totalSumW = totalSumW,
                               lambda = lambda, W = W,
                               nPreRank = nPreRank,
                               tmpmat0 = tmpmat0, tmpmat1 = tmpmat1,
                               tmpmat2 = tmpmat2, tmpmat3 = tmpmat3,
                               tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                               denomvec = denomvec,
                               d1f = d1f, d2f = d2f, offset = offset,
                               nPtsChi2 = nPtsChi2, simnull = simnull);

    # estimate fixed effects
    if isempty(X)
      b = zeros(0);
    else
      Xrot = UV' * X;
      yrot = UV' * y;
      wt = 1.0 ./ sqrt.(vc1 * evalV + vc0);
      Xnew = scale(wt, Xrot);
      ynew = wt .* yrot;
      b = Xnew \ ynew;
    end

    # return values
    return b, vc0, vc1, vc1_pvalue;

  elseif tests == "eScore"

    # fit the null model
    b = X \ y;
    r = y - X * b;
    vc0 = norm(r) ^ 2 / n;
    vc1 = 0;

    # score test statistic
    #statScore = norm(r' * sqrtV) ^ 2 / norm(r) ^ 2;
    statScore = norm(BLAS.gemv('T', sqrtV, r)) ^ 2 / norm(r) ^ 2;
    #statScore = max(statScore, sum(evalV) / n);

    #=
    # obtain p-value for testing vc1=0
    vc1_pvalue = vctestnullsim(statScore, evalV, evalAdjV, n, rankX,
                               WPreSim, test = "eScore",
                               nSimPts = nNullSimPts,
                               pvalueComputing = pvalueComputings,
                               nNewtonIter = nNullSimNewtonIter,
                               device = devices,
                               PrePartialSumW = PrePartialSumW,
                               PreTotalSumW = PreTotalSumW,
                               partialSumWConst = partialSumWConst,
                               totalSumWConst = totalSumWConst,
                               windowSize = windowSize,
                               partialSumW = partialSumW,
                               totalSumW = totalSumW,
                               lambda = lambda, W = W,
                               nPreRank = nPreRank,
                               tmpmat0 = tmpmat0, tmpmat1 = tmpmat1,
                               tmpmat2 = tmpmat2, tmpmat3 = tmpmat3,
                               tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                               denomvec = denomvec,
                               d1f = d1f, d2f = d2f, offset = offset,
                               nPtsChi2 = nPtsChi2, simnull = simnull);
    =#

    if statScore <= sum(evalV) / n
      vc1_pvalue = 1.0;
    else
      w = [evalAdjV - statScore; - statScore];
      dfvec = [ones(Int, rankAdjV); n - rankX - rankAdjV];
      fun(x) = imag(prod((1 - 2 .* w * x * im) .^ (- dfvec / 2)) / x);
      (vc1_pvalue, err) = quadgk(fun, 0, Inf);
      vc1_pvalue = 0.5 + vc1_pvalue / pi;
    end

    # return values
    return b, vc0, vc1, vc1_pvalue;

  end

end


text:
using st16SseqJuliaTools
using Artifacts
using Test
@testset "taxonkit" begin
    db = artifact"taxdump"
    ids = ["535026","595","34","1350","81852","703612"]
    @test ["Bacillus", "Salmonella", "Myxococcus", "Enterococcus", "NA","Bacillus"] == taxid_to_genus(ids,db)
    @test ["Bacillus subtilis", "Salmonella enterica", "Myxococcus xanthus", "NA", "NA", "Bacillus subtilis"] == taxid_to_species(ids,db)
end


text:
using QSM: prolong!, restrict!, restrict_size
using QSM: Poisson2

@testset "Restriction" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(16, 17, 18), (19, 20, 21)]
        szc = restrict_size(sz)

        R = CartesianIndices(map(n -> 2:n-1, sz))
        Rc = CartesianIndices(map(n -> 2:n-1, szc))

        x = zeros(T, sz)
        xc = zeros(T, szc)
        yc = zeros(T, szc)
        mxc = zeros(Bool, szc)
        myc = zeros(Bool, szc)

        x[R] .= 1
        xc[Rc] .= randn(size(Rc))
        yc[Rc] .= randn(size(Rc))

        mxc[Rc] .= true
        myc[Rc] .= rand(Bool, size(Rc))

        dx = Tuple(rand(T, 3))
        Axc = Poisson2(T, mxc, dx)
        Ayc = Poisson2(T, myc, dx)

        restrict!(xc, Axc, x)
        restrict!(yc, Ayc, x)

        @test xc[Rc] ≈ ones(size(Rc))
        @test sum(xc[Rc]) ≈ sum(xc)
        @test myc.*yc ≈ myc.*xc
    end
end

@testset "Prolongation" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(16, 17, 18), (19, 20, 21)]
        szc = restrict_size(sz)

        R = CartesianIndices(map(n -> 3:n-2, sz))
        Rc = CartesianIndices(map(n -> 2:n-1, szc))

        xc = zeros(T, szc)
        x = zeros(T, sz)
        y = zeros(T, sz)
        mx = zeros(Bool, sz)
        my = zeros(Bool, sz)

        xc .= 1
        x[R] .= randn(size(R))
        y[R] .= randn(size(R))

        mx[R] .= true
        my[R] .= rand(Bool, size(R))

        dx = Tuple(rand(T, 3))
        Ax = Poisson2(T, mx, dx)
        Ay = Poisson2(T, my, dx)

        prolong!(x, Ax, xc)
        prolong!(y, Ay, xc)

        @test x[R] ≈ ones(size(R))
        @test my.*y ≈ my.*x

        R = CartesianIndices(map(n -> 2:n-1, sz))
        @test sum(x[R]) ≈ sum(x)
    end
end

@testset "R = αP^T" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(16, 17, 18), (19, 20, 21)]
        szc = restrict_size(sz)

        u = zeros(T, sz)
        v = zeros(T, sz)
        m = zeros(Bool, sz)

        uc = zeros(T, szc)
        vc = zeros(T, szc)
        mc = zeros(Bool, szc)

        R = CartesianIndices(map(n -> 2:n-1, sz))
        Rc = CartesianIndices(map(n -> 2:n-1, szc))

        m[R] .= true
        v[R] .= randn(size(R))
        u[R] .= randn(size(R))

        mc[Rc] .= true
        uc[Rc] .= randn(size(Rc))
        vc[Rc] .= randn(size(Rc))

        dx = Tuple(rand(T, 3))
        A = Poisson2(T, m, dx)

        dxc = ((sz.-1)./(szc.-1)) .* dx
        Ac = Poisson2(T, mc, dxc)

        prolong!(u, A, uc)
        restrict!(vc, Ac, v)

        R = CartesianIndices(map(n -> 1:(iseven(n) ? n-2 : n), sz))
        @test dot(uc, vc) ≈ dot(v[R], inv(8).*u[R])
    end
end


text:
function Base.hash(m::AbstractMonomial, u::UInt)
    nnz = count(!iszero, exponents(m))
    if iszero(nnz)
        hash(1, u)
    elseif isone(nnz) && isone(degree(m))
        # We want the has of `m` to match the hash of `variable(m)`
        # so we ignore the exponent one.
        hash(variable(m), u)
    else
        reduce((u, power) -> iszero(power[2]) ? u : hash(power, u), powers(m), init=u)
    end
end

function Base.hash(t::AbstractTerm, u::UInt)
    if iszero(t)
        hash(0, u)
    elseif isone(coefficient(t))
        # We want the has of `t` to match the hash of `monomial(t)`
        # so we ignore the coefficient one.
        hash(monomial(t), u)
    else
        hash(monomial(t), hash(coefficient(t), u))
    end
end

function Base.hash(p::AbstractPolynomial, u::UInt)
    if iszero(p)
        hash(0, u)
    else
        reduce((u, t) -> hash(t, u), terms(p); init=u)
    end
end


text:
@testset "$TEST $G" begin

#tests for the concrete types Graph and DiGraph
g = DiGraph(10, 20)
@test  nv(g) == 10
@test  ne(g) == 20
@test typeof(g) == DiGraph{Int}
@test vertextype(g) == Int
@test edgetype(g) == Edge{Int}

g = Graph(10, 20)
@test  nv(g) == 10
@test  ne(g) == 20
@test typeof(g) == Graph{Int}
@test edgetype(g) == Edge{Int}

w = sprand(10,10, 0.1)
@test Graph(w) == Graph{Int}(w)
@test DiGraph(w) == DiGraph{Int}(w)

@test Graph(0) == Graph() ==Graph{Int}() == Graph{Int}(0)
@test DiGraph(0) == DiGraph() == DiGraph{Int}() == DiGraph{Int}(0)

end # testset


text:
using Test
using Awaits
using Awaits: SyncSeq

function tickingfor(ctx, sec)
    dt = 0.01
    @in ctx begin
        for i in 1:Int(cld(sec, dt))
            sleep(dt)
            @check
        end
    end
end


text:
#! format: off

using SIIPExamples #hide
using PowerSimulationsDynamics
PSID = PowerSimulationsDynamics
using PowerSystems
using Sundials
using Plots
gr()

file_dir = joinpath(
    dirname(dirname(pathof(SIIPExamples))),
    "script",
    "4_PowerSimulationsDynamics_examples",
    "Data",
)
omib_sys = System(joinpath(file_dir, "omib_sys.json"))

time_span = (0.0, 30.0)
perturbation_trip = BranchTrip(1.0, "BUS 1-BUS 2-i_1")
sim = PSID.Simulation(PSID.ImplicitModel, omib_sys, pwd(), time_span, perturbation_trip)

print_device_states(sim)

x0_init = PSID.get_initial_conditions(sim)

PSID.execute!(
    sim, #simulation structure
    IDA(), #Sundials DAE Solver
    dtmax = 0.02,
); #Arguments: Maximum timestep allowed

angle = get_state_series(sim, ("generator-102-1", :δ));
Plots.plot(angle, xlabel = "time", ylabel = "rotor angle [rad]", label = "rotor angle")

volt = get_voltage_magnitude_series(sim, 102);
Plots.plot(volt, xlabel = "time", ylabel = "Voltage [pu]", label = "V_2")

sim2 = PSID.Simulation(PSID.ImplicitModel, omib_sys, pwd(), time_span, perturbation_trip)

small_sig = small_signal_analysis(sim2)

small_sig.reduced_jacobian

small_sig.eigenvalues

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl