text:
"""
    Dimension

Abstract supertype of all dimension types.

Example concrete implementations are [`X`](@ref), [`Y`](@ref), [`Z`](@ref),
[`Ti`](@ref) (Time), and the custom [`Dim`]@ref) dimension.

`Dimension`s label the axes of an `AbstractDimArray`,
or other dimensional objects, and are used to index into the array.

They may also provide an alternate index to lookup for each array axis.
This may be any `AbstractVector` matching the array axis length, or a `Val`
holding a tuple for compile-time index lookups.

`Dimension`s also have `lookup` and `metadata` fields.

`lookup` gives more details about the dimension, such as that it is
[`Categorical`](@ref) or [`Sampled`](@ref) as [`Points`](@ref) or
[`Intervals`](@ref) along some transect. DimensionalData will
attempt to guess the lookup from the passed-in index value.

Example:

```jldoctest Dimension
using DimensionalData, Dates

x = X(2:2:10)
y = Y(['a', 'b', 'c'])
ti = Ti(DateTime(2021, 1):Month(1):DateTime(2021, 12))

A = DimArray(zeros(3, 5, 12), (y, x, ti))

# output

3×5×12 DimArray{Float64,3} with dimensions:
  Y Categorical Char[a, b, c] ForwardOrdered,
  X Sampled 2:2:10 ForwardOrdered Regular Points,
  Ti Sampled DateTime("2021-01-01T00:00:00"):Month(1):DateTime("2021-12-01T00:00:00") ForwardOrdered Regular Points
[:, :, 1]
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
[and 11 more slices...]
```

For simplicity, the same `Dimension` types are also used as wrappers
in `getindex`, like:

```jldoctest Dimension
x = A[X(2), Y(3)]

# output

12-element DimArray{Float64,1} with dimensions:
  Ti Sampled DateTime("2021-01-01T00:00:00"):Month(1):DateTime("2021-12-01T00:00:00") ForwardOrdered Regular Points
and reference dimensions:
  Y Categorical Char[c] ForwardOrdered,
  X Sampled 4:2:4 ForwardOrdered Regular Points
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
```

A `Dimension` can also wrap [`Selector`](@ref).

```jldoctest Dimension
x = A[X(Between(3, 4)), Y(At('b'))]

# output

1×12 DimArray{Float64,2} with dimensions:
  X Sampled 4:2:4 ForwardOrdered Regular Points,
  Ti Sampled DateTime("2021-01-01T00:00:00"):Month(1):DateTime("2021-12-01T00:00:00") ForwardOrdered Regular Points
and reference dimensions:
  Y Categorical Char[b] ForwardOrdered
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```

`Dimension` objects may have [`lookup`](@ref) and [`metadata`](@ref) fields
to track additional information about the data and the index, and their relationship.
"""
abstract type Dimension{T} end

"""
    IndependentDim <: Dimension

Abstract supertype for independent dimensions. Thise will plot on the X axis.
"""
abstract type IndependentDim{T} <: Dimension{T} end

"""
    DependentDim <: Dimension

Abstract supertype for Dependent dimensions. These will plot on the Y axis.
"""
abstract type DependentDim{T} <: Dimension{T} end

"""
    XDim <: IndependentDim

Abstract supertype for all X dimensions.
"""
abstract type XDim{T} <: IndependentDim{T} end

"""
    YDim <: DependentDim

Abstract supertype for all Y dimensions.
"""
abstract type YDim{T} <: DependentDim{T} end

"""
    ZDim <: DependentDim

Abstract supertype for all Z dimensions.
"""
abstract type ZDim{T} <: DependentDim{T} end

"""
    TimeDim <: IndependentDim

Abstract supertype for all time dimensions.

In a `TimeDime` with `Interval` sampling the locus will automatically
be set to `Start()`. Dates and times generally refer to the start of a
month, hour, second etc., not the central point as is more common with spatial data.
`"""
abstract type TimeDim{T} <: IndependentDim{T} end

ConstructionBase.constructorof(d::Type{<:Dimension}) = basetypeof(d)
Adapt.adapt_structure(to, dim::Dimension) = rebuild(dim; val=Adapt.adapt(to, val(dim)))

const DimType = Type{<:Dimension}
const DimTuple = Tuple{<:Dimension,Vararg{<:Dimension}}
const SymbolTuple = Tuple{Symbol,Vararg{Symbol}}
const DimTypeTuple = Tuple{<:DimType,Vararg{<:DimType}}
const VectorOfDim = Vector{<:Union{Dimension,DimType,Symbol}}
const DimOrDimType = Union{Dimension,DimType,Symbol}
const AllDims = Union{Symbol,Dimension,DimTuple,SymbolTuple,DimType,DimTypeTuple,VectorOfDim}

# DimensionalData interface methods

"""
    rebuild(dim::Dimension, val) => Dimension
    rebuild(dim::Dimension; val=val(dim)) => Dimension

Rebuild dim with fields from `dim`, and new fields passed in.
"""
function rebuild(dim::D, val) where D <: Dimension
    ConstructionBase.constructorof(D)(val)
end

dims(dim::Union{Dimension,DimType,Val{<:Dimension}}) = dim
dims(dims::DimTuple) = dims
dims(::Tuple{}) = ()
dims(x) = nothing
dims(::Nothing) = error("No dims found")

val(dim::Dimension) = dim.val
refdims(x) = ()

lookup(dim::Dimension{<:AbstractArray}) = val(dim)
lookup(dim::Union{DimType,Val{<:Dimension}}) = NoLookup()

lookuptype(dim::Dimension) = typeof(lookup(dim))
lookuptype(::Type{<:Dimension{L}}) where L = L
lookuptype(x) = NoLookup

# LookupArrays methods
metadata(dim::Dimension) = metadata(lookup(dim))

index(dim::Dimension{<:AbstractArray}) = index(val(dim))
index(dim::Dimension{<:Val}) = unwrap(index(val(dim)))

name(dim::Dimension) = name(typeof(dim))
name(dim::Val{D}) where D = name(D)

label(x) = string(string(name(x)), (units(x) === nothing ? "" : string(" ", units(x))))

bounds(dim::Dimension) = bounds(val(dim))
shiftlocus(locus::Locus, dim::Dimension) = rebuild(dim, shiftlocus(locus, lookup(dim)))
maybeshiftlocus(locus::Locus, d::Dimension) = rebuild(d, maybeshiftlocus(locus, lookup(d)))

function hasselection(x, selectors::Union{DimTuple,SelTuple,Selector,Dimension})
    hasselection(dims(x), selectors)
end
hasselection(x::Nothing, selectors::Union{DimTuple,SelTuple,Selector,Dimension}) = false
function hasselection(ds::DimTuple, seldims::DimTuple)
    sorted = dims(seldims, ds)
    hasselection(dims(ds, sorted), map(val, sorted))
end
hasselection(ds::DimTuple, selectors::SelTuple) = all(map(hasselection, ds, selectors))
hasselection(ds::DimTuple, selector::Dimension) = hasselection(dims(ds, selector), selector)
function hasselection(ds::DimTuple, selector::Selector)
    throw(ArgumentError("Cannot select from multiple Dimensions with a single Selector"))
end
hasselection(dim::Dimension, seldim::Dimension) = hasselection(dim, val(seldim))
hasselection(dim::Dimension, sel::Selector) = hasselection(lookup(dim), sel)

for func in (:order, :span, :sampling, :locus)
    @eval ($func)(dim::Dimension) = ($func)(lookup(dim))
end

# Dipatch on Tuple{<:Dimension}, and map to single dim methods
for f in (:val, :index, :lookup, :metadata, :order, :sampling, :span, :bounds, :locus,
          :name, :label, :units)
    @eval begin
        $f(ds::DimTuple) = map($f, ds)
        $f(ds::Tuple{}) = ()
        $f(ds::DimTuple, i1, I...) = $f(ds, (i1, I...))
        $f(ds::DimTuple, I) = $f(dims(ds, key2dim(I)))
    end
end

@inline function selectindices(x, selectors)
    if dims(x) isa Nothing
        # This object has no dimensions and no `selectindices` method.
        # Just return whatever it is, maybe the underlying array can use it.
        return selectors
    else
        # Otherwise select indices based on the object `Dimension`s
        return selectindices(dims(x), selectors)
    end
end
@inline selectindices(ds::DimTuple, sel...) = selectindices(ds, sel)
@inline selectindices(ds::DimTuple, sel::Tuple) = selectindices(val(ds), sel)
@inline selectindices(dim::Dimension, sel) = selectindices(val(dim), sel)

# Base methods
const ArrayOrVal = Union{AbstractArray,Val}

Base.parent(d::Dimension) = val(d)
Base.eltype(d::Type{<:Dimension{T}}) where T = T
Base.eltype(d::Type{<:Dimension{A}}) where A<:AbstractArray{T} where T = T
Base.size(d::Dimension, args...) = size(val(d), args...)
Base.axes(d::Dimension, args...) = axes(val(d), args...)
Base.eachindex(d::Dimension) = eachindex(val(d))
Base.length(d::Dimension) = length(val(d))
Base.ndims(d::Dimension) = 0
Base.ndims(d::Dimension{<:AbstractArray}) = ndims(val(d))
Base.iterate(d::Dimension{<:AbstractArray}, args...) = iterate(lookup(d), args...)
Base.first(d::Dimension) = val(d)
Base.first(d::Dimension{<:AbstractArray}) = first(lookup(d))
Base.last(d::Dimension) = val(d)
Base.last(d::Dimension{<:AbstractArray}) = last(lookup(d))
Base.firstindex(d::Dimension) = 1
Base.lastindex(d::Dimension) = 1
Base.firstindex(d::Dimension{<:AbstractArray}) = firstindex(lookup(d))
Base.lastindex(d::Dimension{<:AbstractArray}) = lastindex(lookup(d))
Base.step(d::Dimension) = step(lookup(d))
Base.Array(d::Dimension{<:AbstractArray}) = collect(lookup(d))
function Base.:(==)(d1::Dimension, d2::Dimension)
    basetypeof(d1) == basetypeof(d2) && val(d1) == val(d2)
end

Base.size(dims::DimTuple) = map(length, dims)

# Extents.jl
function Extents.extent(ds::DimTuple, args...)
    extent_dims = _astuple(dims(ds, args...))
    extent_bounds = bounds(extent_dims)
    return Extents.Extent{dim2key(extent_dims)}(extent_bounds)
end

dims(extent::Extents.Extent{K}) where K = map(rebuild, key2dim(K), values(extent))
dims(extent::Extents.Extent, ds) = dims(dims(extent), ds)

# Produce a 2 * length(dim) matrix of interval bounds from a dim
dim2boundsmatrix(dim::Dimension)  = dim2boundsmatrix(lookup(dim))
function dim2boundsmatrix(lookup::LookupArray)
    samp = sampling(lookup)
    samp isa Intervals || error("Cannot create a bounds matrix for $(nameof(typeof(samp)))")
    _dim2boundsmatrix(locus(lookup), span(lookup), lookup)
end

_dim2boundsmatrix(::Locus, span::Explicit, lookup) = val(span)
function _dim2boundsmatrix(::Locus, span::Regular, lookup)
    starts = permutedims(LookupArrays._shiftindexlocus(Start(), lookup))
    ends = permutedims(LookupArrays._shiftindexlocus(End(), lookup))
    return vcat(starts, ends)
end
@noinline _dim2boundsmatrix(::Center, span::Regular{Dates.TimeType}, lookupj) =
    error("Cannot convert a Center TimeType index to Explicit automatically: use a bounds matrix e.g. Explicit(bnds)")
@noinline _dim2boundsmatrix(::Start, span::Irregular, lookupj) =
    error("Cannot convert Irregular to Explicit automatically: use a bounds matrix e.g. Explicit(bnds)")

"""
    Dim{S}(val=:)

A generic dimension. For use when custom dims are required when loading
data from a file. Can be used as keyword arguments for indexing.

Dimension types take precedence over same named `Dim` types when indexing
with symbols, or e.g. creating Tables.jl keys.

```jldoctest
using DimensionalData

dim = Dim{:custom}(['a', 'b', 'c'])

# output

Dim{:custom} Char[a, b, c]
```
"""
struct Dim{S,T} <: Dimension{T}
    val::T
end
Dim{S}(val::T) where {S,T} = Dim{S,T}(val)
function Dim{S}(val::AbstractArray; kw...) where S
    if length(kw) > 0
        val = AutoLookup(val, values(kw))
    end
    Dim{S,typeof(val)}(val)
end
Dim{S}() where S = Dim{S}(:)

name(::Type{<:Dim{S}}) where S = S
basetypeof(::Type{<:Dim{S}}) where S = Dim{S}
key2dim(s::Val{S}) where S = Dim{S}()
dim2key(::Type{D}) where D<:Dim{S} where S = S

"""
    AnonDim <: Dimension

    AnonDim()

Anonymous dimension. Used when extra dimensions are created,
such as during transpose of a vector.
"""
struct AnonDim{T} <: Dimension{T}
    val::T
end
AnonDim() = AnonDim(Colon())
AnonDim(val, arg1, args...) = AnonDim(val)

lookup(::AnonDim) = NoLookup()
metadata(::AnonDim) = NoMetadata()
name(::AnonDim) = :Anon

"""
    @dim typ [supertype=Dimension] [name::String=string(typ)]

Macro to easily define new dimensions. The supertype will be inserted
into the type of the dim. The default is simply `YourDim <: Dimension`. Making
a Dimesion inherit from `XDim`, `YDim`, `ZDim` or `TimeDim` will affect
automatic plot layout and other methods that dispatch on these types. `<: YDim`
are plotted on the Y axis, `<: XDim` on the X axis, etc.

Example:
```julia
@dim Lat YDim "latitude"
@dim Lon XDim "Longitude"
```
"""
macro dim end
macro dim(typ::Symbol, args...)
    dimmacro(typ::Symbol, :(DimensionalData.Dimension), args...)
end
macro dim(typ::Symbol, supertyp::Symbol, args...)
    dimmacro(typ, supertyp, args...)
end

function dimmacro(typ, supertype, name::String=string(typ))
    quote
        Base.@__doc__ struct $typ{T} <: $supertype{T}
            val::T
        end
        function $typ(val::AbstractArray; kw...)
            if length(kw) > 0
                val = AutoLookup(val, values(kw))
            end
            $typ{typeof(val)}(val)
        end
        $typ() = $typ(:)
        Dimensions.name(::Type{<:$typ}) = $(QuoteNode(Symbol(name)))
        Dimensions.key2dim(::Val{$(QuoteNode(typ))}) = $typ()
    end |> esc
end

# Define some common dimensions.

"""
    X <: XDim

    X(val=:)

X [`Dimension`](@ref). `X <: XDim <: IndependentDim`

## Example:
```julia
xdim = X(2:2:10)
# Or
val = A[X(1)]
# Or
mean(A; dims=X)
```
"""
@dim X XDim

"""
    Y <: YDim

    Y(val=:)

Y [`Dimension`](@ref). `Y <: YDim <: DependentDim`

## Example:
```julia
ydim = Y(['a', 'b', 'c'])
# Or
val = A[Y(1)]
# Or
mean(A; dims=Y)
```
"""
@dim Y YDim

"""
    Z <: ZDim

    Z(val=:)

Z [`Dimension`](@ref). `Z <: ZDim <: Dimension`

## Example:
```julia
zdim = Z(10:10:100)
# Or
val = A[Z(1)]
# Or
mean(A; dims=Z)
```
"""
@dim Z ZDim

"""m
    Ti <: TimeDim

    Ti(val=:)

Time [`Dimension`](@ref). `Ti <: TimeDim <: IndependentDim`

`Time` is already used by Dates, and `T` is a common type parameter,
We use `Ti` to avoid clashes.

## Example:
```julia
timedim = Ti(DateTime(2021, 1):Month(1):DateTime(2021, 12))
# Or
val = A[Ti(1)]
# Or
mean(A; dims=Ti)
```
"""
@dim Ti TimeDim "Time"

const Time = Ti # For some backwards compat


text:
"""
    module Unitful_U_str

Temporary solution to fix issue with Unitful.
For details see: https://github.com/ModiaSim/Modia.jl/issues/49

This module exports string macro @U_str. It is then possible to use:

```julia
using Unitful
using Modia    # Unitful_U_str.jl is included in Modia.jl

const T1_type = u"N*m"    # gives warning
const T2_type = u"N*m"    # gives no warning
```
"""
module Unitful_U_str

export @U_str

using Unitful

macro U_str(unit)
    ex = Meta.parse(unit)
    esc(replace_value(ex))
end

const allowed_funcs = [:*, :/, :^, :sqrt, :√, :+, :-, ://]
function replace_value(ex::Expr)
    if ex.head == :call
        ex.args[1] in allowed_funcs ||
            error("""$(ex.args[1]) is not a valid function call when parsing a unit.
             Only the following functions are allowed: $allowed_funcs""")
        for i=2:length(ex.args)
            if typeof(ex.args[i])==Symbol || typeof(ex.args[i])==Expr
                ex.args[i]=replace_value(ex.args[i])
            end
        end
        return ex
    elseif ex.head == :tuple
        for i=1:length(ex.args)
            if typeof(ex.args[i])==Symbol
                ex.args[i]=replace_value(ex.args[i])
            else
                error("only use symbols inside the tuple.")
            end
        end
        return ex
    else
        error("Expr head $(ex.head) must equal :call or :tuple")
    end
end

replace_value(sym::Symbol) = Unitful.replace_value(sym)

replace_value(literal::Number) = literal

end

ext: jl


text:
using AndTraits
const Contract = Trait{:Contract}
const Err = Trait{:Err}


text:
module TestRedis
using Test
using Redis
using StructDatabaseMapping
using Dates

DB_NUMBER = 0

include("../includes/basic_test.jl")

function test()
    test_redis()
end
params = (host=get(ENV, "REDIS_HOST", "localhost"),
          db=DB_NUMBER,
          port=parse(Int64, get(ENV, "REDIS_PORT", "6379")))
function cleanup()
    Redis.flushall(
        Redis.RedisConnection(;params...))
end
function test_redis()
    f = ()->Redis.RedisConnection(;params...)
    _test_basic_functionalities(f)
end

end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
using Distances, StructuredDistances, InteractionNetworkModels
using Plots, BenchmarkTools, Distributions, StatsPlots

d = Cityblock()

V = 9
τ = 0.2
mode = [rand() < 0.1 ? rand(1:4) : 0 for i in 1:V,j in 1:V]
γ = 1.9

model = SNF(mode, γ, d)
model = vectorise(model)

# Gibbs scan
gibbs_move_rand = GibbsRandMove(ν=1)
gibbs_move = GibbsScanMove(ν=1)

mcmc = McmcSampler(gibbs_move_rand)
mcmc_scan = McmcSampler(gibbs_move)

@time x = mcmc(model, desired_samples=10000, lag=1)
plot(x)

sum(mode)

acceptance_prob(mcmc_scan)
acceptance_prob(mcmc)
typeof(model)

# Testing posterior sampler
@time out  = mcmc_scan(model, desired_samples=100, lag=5, burn_in=500)
plot(out)
similar(model, rand(1:10, 20), 1.0)

model.mode
data = out.sample
γ_prior = Gamma(model.γ, 1)
plot(γ_prior)
G_prior = SNF(
    model.mode,
    0.1,
    model.d,
    directed=model.directed,
    self_loops=model.self_loops
)
posterior = SnfPosterior(data, G_prior, γ_prior)
mode_move = GibbsScanMove(ν=1)
aux_mcmc = McmcSampler(GibbsRandMove(ν=1), burn_in=3000, lag=50)

@time x = aux_mcmc(model, desired_samples=100)
acceptance_prob(aux_mcmc)
plot(x)

mcmc_posterior = SnfPosteriorSampler(
    mode_move, aux_mcmc, posterior,
    ε=0.1,
    aux_init_at_prev=true
)
mcmc_posterior.aux.data

@time x = mcmc_posterior(posterior, desired_samples=500, γ_init=4.0)

plot(x, model.mode)
x_sample_mat = get_sample_matrices(x)[end]

data

# Plot graph
using CairoMakie, GraphMakie, Graphs

g = SimpleDiGraph(x_sample_mat[1])

f,ax,p = graphplot(
    g,
    node_color=:blue,
    edge_plottype=:linesegments,
    showaxis=false
)
hidedecorations!(ax)  # hides ticks, grid and lables
hidespines!(ax)  # hide the frame
f

V = 15
g = SimpleDiGraph(
    [rand() < 0.5 ? true : false for i in 1:V, j in 1:V]
)

f, ax, p = graphplot(
    g,
    node_size=10,
    node_color=:red,
    edge_plottype=:linesegments
)
hidedecorations!(ax)  # hides ticks, grid and lables
hidespines!(ax)  # hide the frame
f

max_stars_repo_path: Examples/Graph Models/SNF_test.jl

max_issues_repo_path: Examples/Graph Models/SNF_test.jl

max_forks_repo_path: Examples/Graph Models/SNF_test.jl


text:
"""
    Kernel PCA object
"""
struct kPCA{T <: Real}
    X::Matrix{Union{T, Missing}} # original data
    κ::Kernel{T} # kernel function (from MLKernels)
    μ::Vector{T} # row/col means of the kernel matrix
    μ2::T        # mean of the kernel matrix
    λ::Vector{T} # eigenvalues in feature space
    α::Matrix{T} # eigenvectors in feature space

    function kPCA(
        X  :: Matrix{Union{T, Missing}},
        κ  :: Kernel{T},
        μ  :: Vector{T},
        μ2 :: T,
        λ  :: Vector{T},
        α  :: Matrix{T}
    ) where T <: Real
        @assert length(λ) == size(α, 2)
        @assert length(μ) == size(X, 2)

        new{T}(X, κ, μ, μ2, λ, α)
    end
end

indim(M::kPCA) = size(M.X, 1)
outdim(M::kPCA) = length(M.λ)

projection(M::kPCA) = M.V ./ sqrt.(M.λ)
principalvars(M::kPCA) = M.λ

"""
fit a kernel PCA.

    κ :: Kernel{T} the kernel function, see MLKernels.jl for documentation.
    X :: Array{Union{T, Missing}} the data to embed, with observations ins columns.
"""
function fit(
      :: Type{kPCA},
    κ :: Kernel{T},
    X :: Matrix{Union{T, Missing}},
    ncomp = 2
) where T

    K = kernelmatrix(Val(:col), κ, X, true)
    K2, μ, μ2 = missing_to_mean(K)
    K2 .= K2 .- μ .- μ' .+ μ2

    e = K2 |> Hermitian |> eigen
    e_ev_perm = sortperm(e.values, rev = true)[1:ncomp]

    λ = e.values[e_ev_perm]::Vector{T}
    α = e.vectors[:, e_ev_perm]::Matrix{T}
    return kPCA(X, κ, μ, μ2, λ, α)
end

"""Calculate transformation to kernel space"""
# function transform(M::kPCA{T}, x::AbstractMatrix{mT}) where mT <: Union{T, Missing} where T <: Real
function transform(M::kPCA{T}, x::AbstractMatrix{mT}) where mT where T

    ##### There is some typeinference bug that leads to illegal instructions and
    ##### reaches the unreachable...

    K = kernelmatrix(Val(:col), M.κ, M.X, x)
#     ##### instead of `kernelmatrix`:
#     # K = allocate_basematrix(Val(:col), M.X, x)
#     K = Array{T}(undef, size(M.X, 2), size(x, 2))
#     f = basefunction(M.κ)

#     # basematrix!(Val(:col), K, f, M.X, x)
#     ##### instead of `basematrix!`:
#     n, m = checkdimensions(Val(:col), K, M.X, x)
#     @inbounds for j = 1:m
#         yj = subvector(Val(:col), x, j)
#         for i = 1:n
#             xi = subvector(Val(:col), M.X, i)

# ##### this causes slightly less allocations but far worse performance at least it doesn't crash the julia session
# # # unsafe_base_evaluate
# #     s, nn = base_initiate(f, T), 0
# #     # s, nn = base_initiate(f, promote_type(eltype(x), eltype(y))), 0
# #     for i in eachindex(xi, yj)
# #         s, nn = base_aggregate(f, s, nn, xi[i], yj[i])
# #     end
# #     K[i, j] = nn > 0 ? base_return(f, s / nn) : missing
# # # end unsafe_base_evaluate

#             K[i,j] = unsafe_base_evaluate(f, xi, yj)
#         end
#     end

#     ##### end instead of `basematrix!`

#     kappamatrix!(M.κ, K)
#     ##### end instead of `kernelmatrix`

    K .= K .- M.μ .- mean(K, dims = 1) .+ M.μ2

    return (M.α' ./ sqrt.(M.λ)) * K
end

transform(M::kPCA) = sqrt.(M.λ) .* M.α'

function Base.show(io::IO, ::MIME"text/plain", M::kPCA{T}) where T
    println(io, "kPCA{", T, "}(κ: ", M.κ, ", dims: ",  indim(M), "─→", outdim(M), ")")
end

# """Kernel PCA type"""
# struct KernelPCA{T<:Real}
#     X::AbstractMatrix{T}           # fitted data or precomputed kernel
#     ker::Union{Nothing, Function}  # kernel function
#     center::KernelCenter           # kernel center
#     λ::AbstractVector{T}           # eigenvalues  in feature space
#     α::AbstractMatrix{T}           # eigenvectors in feature space
#     inv::AbstractMatrix{T}         # inverse transform coefficients
# end
# struct KernelCenter{T<:Real}
#     means::AbstractVector{T}
#     total::T
# end

# import MultivariateStats: KernelPCA, KernelCenter, transform!
# import Base.convert

# """Center kernel matrix."""
# function transform!(C::KernelCenter{T}, K::AbstractMatrix{T}) where {T<:Real}
#     r, c = size(K)
#     tot = C.total
#     means = mean(K, dims=1)
#     K .= K .- C.means .- means .+ tot
#     return K
# end

# TODO: does not work
# function convert(::Type{KernelPCA{T}}, M::kPCA) where T
#     KernelPCA(missing_to_mean_slice(M.X, dims = 2),
#               (x, y) -> kappa(M.κ, unsafe_base_evaluate(basefunction(M.κ), x, y)),
#               KernelCenter(M.μ, M.μ2),
#               M.λ,
#               M.α,
#               zeros(T, 0, 0))
# end


text:
using Symbolics
using Symbolics: symbol_to_poly, poly_to_symbol

using Test

@variables x y z

syms = [
    [1], [1, BigInt(2)], [x], [x^2 + 2], [x + (5//8)y],
    [1 - x*y*z], [(x + y + z)^5], [0, BigInt(0)y^2, Rational(1)z^3],
    [x, sin((44//31)y)*z]
]
for sym in syms
    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(sym)
    sym2 = Symbolics.poly_to_symbol(polynoms, pvar2sym, sym2term, Real)
    @test isequal(expand.(sym2), expand.(sym))
end

@test isequal(expand.(groebner_basis([x, y])), [y, x])
@test isequal(expand.(groebner_basis([y, x])), [y, x])
@test isequal(expand.(groebner_basis([x])), [x])
@test isequal(expand.(groebner_basis([x, x^2])), [x])
@test isequal(expand.(groebner_basis([BigInt(1)x + 2//3])), [x  + 2//3])

@variables x1 x2 x3 x4
@test isequal(expand.(groebner_basis([x1, x, y])), [y, x1, x])

# input unchanged
f1 = [x2, x1, x4, x3]
f2 = [x2, x1, x4, x3]
groebner_basis(f1)
@test isequal(expand.(f1), expand.(f2))

@variables x1 x2 x3 x4 x5
system = [
   x1 + x2 + x3 + x4 + x5,
   x1*x2 + x1*x3 + x1*x4 + x1*x5 + x2*x3 + x2*x4 + x2*x5 + x3*x4 + x3*x5 + x4*x5,
   x1*x2*x3 + x1*x2*x4 + x1*x2*x5 + x1*x3*x4 + x1*x3*x5 + x1*x4*x5 + x2*x3*x4 + x2*x3*x5 + x2*x4*x5 + x3*x4*x5,
   x1*x2*x3*x4 + x1*x2*x3*x5 + x1*x2*x4*x5 + x1*x3*x4*x5 + x2*x3*x4*x5,
   x1*x2*x3*x4*x5 - 1
]
truebasis = [
        x1 + x2 + x3 + x4 + x5,
        x2^2 + x2*x3 + x2*x4 + x2*x5 + x3^2 + x3*x4 + x3*x5 + x4^2 + x4*x5 + x5^2,
        x3^3 + x3*(x4^2) + x3*(x5^2) + x4^3 + x4*(x3^2) + x5*(x4^2) + x4*(x5^2) + x5^3 + x5*(x3^2) + x3*x4*x5,
        x4^4 + x4*(x5^3) + x5^4 + x5*(x4^3) + (x4^2)*(x5^2),
        x5^5 - 1
]
basis = expand.(groebner_basis(system))
@test isequal(basis, truebasis)

N = 45671930739135174346839766056203605080877915151
system = [
    x1 + x2 + x3 + x4,
    x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x3*x4,
    x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4,
    x1*x2*x3*x4 + N
]
truebasis = [
    x1 + x2 + x3 + x4,
    x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2,
    x3^3 + x3^2*x4 + x3*x4^2 + x4^3,
    x4^4 - N
]
basis = groebner_basis(system)
@test isequal(expand.(basis), truebasis)

# Groebner does not yet work with constant ideals
@test_broken groebner_basis([1])


text:
#На прямоугольном поле произвольных размеров расставить маркеры в виде "шахматных" клеток, начиная с юго-западного угла поля,
#когда каждая отдельная "шахматная" клетка имеет размер n x n клеток поля (n - это параметр функции).
#Начальное положение Робота - произвольное, конечное - совпадает с начальным. Клетки на севере и востоке могут получаться
#"обрезанными" - зависит от соотношения размеров поля и "шахматных" клеток.
include("library.jl")
function Chess_with_parameter(r, N::Integer)
    back_path = BackPath(r, (Nord, West))
    R = NChessRobot(r, Coords(0, 0), N)
    marks!(R, (Sud, Ost))
    back!(R, back_path)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
#for debug only
# using CRC32c
# cksum(t::Array)=crc32c(collect(reinterpret(UInt8,t)))
# cksum(t::Array, x::Int32)=crc32c(collect(reinterpret(UInt8,t)), x)

import Base: isless, zero, iterate, getindex, +, -, *, /, ==, !=, mod

struct OpTuple{N,T}
    v::NTuple{N,T}
end
Base.show(io::IO, ot::OpTuple) = print(io, ot.v)

(+)(v1::OpTuple,v2::OpTuple) = OpTuple(v1.v .+ v2.v)
(-)(v1::OpTuple,v2::OpTuple) = OpTuple(v1.v .- v2.v)
(-)(v::OpTuple) = OpTuple(0 .- v.v)
(==)(v1::OpTuple,v2::OpTuple) = v1.v == v2.v
(!=)(v1::OpTuple,v2::OpTuple) = v1.v != v2.v
(*)(a::Number, v::OpTuple) = OpTuple( a .* v.v)
(*)(v::OpTuple, a::Number) = OpTuple( v.v .* a)
(/)(v::OpTuple, a::Number) = OpTuple( v.v ./ a)
Base.zero(::Type{OpTuple{N,T}}) where{N,T}=OpTuple(ntuple(x->zero(T),N))
Base.iterate(v::OpTuple)=iterate(v.v)
Base.iterate(v::OpTuple, state)=iterate(v.v, state)
Base.length(::OpTuple{N}) where N=N
Base.getindex(ot::OpTuple, ind...)=getindex(ot.v,ind...)
Base.mod(v1::OpTuple,v2::OpTuple) = OpTuple(mod.(v1.v, v2.v))

# contruct nb iterators thta is the split of 1:lgtot iterator
function splititr(nb, lgtot)
    lg, r = divrem(lgtot, nb)
    return vcat(
        map(x -> ((x-1)*(lg+1)+1):x*(lg+1), 1:r),
        map(x -> ((x-1)*lg+r+1):(x*lg+r), (r+1):nb),
    )
end

function splitvec(nb, v)
    return map(x -> v[x], splititr(nb, length(v)))
end

# construct the permutation of size n that is the transposition of a and b
function transposition(a, b, n)
    p = collect(1:n)
    p[a], p[b] = b, a
    return p
end

# convert a vector to a tuple
totuple(v) = Tuple(x for x in v)

# convert a tuple to a vector
tovector(t) = [x for x in t]

# construct a tuple of size nb, with ones except at index ind the value sz
tupleshape(ind::Int, nb::Int, sz::Int) = Tuple(((x == ind) ? sz : 1) for x = 1:nb)

# construct an array with nb dims, mesh.points on dim=ind, the other dims have a size of one
function tupleshape(ind::Int, nb::Int, v::Vector{T}) where {T}
    return reshape(v, tupleshape(ind, nb, length(v)))
end
function dotprod(t_v::NTuple{N,Vector{T}}) where {N,T}
    #    res = ones(T,totuple(ones(Int,N))) # array of N dimensions with only one one.
    res = ones(T, ntuple(x -> 1, N)) # array of N dimensions with only one one.
    for (ind, v) in enumerate(t_v)
        res = res .* tupleshape(ind, N, v)
    end
    return res
end
dotprod(t_v::Tuple{Vector{T}}) where {T} = t_v[1]
function dotprodother(t_v::NTuple{N,Vector{T}}) where {N,T}
    return prod.(Iterators.product(t_v...))
end

# modulo or div for "begin to one" array
# an other argument for "begin to zero" array ...
@inline modone(ind, n) = mod(ind - 1, n) + 1
@inline modone(ind::CartesianIndex, n) = CartesianIndex(modone.(ind.I, n))
divone(ind, n) = div(ind - 1, n) + 1
gettabmod(lg) = modone.(1:5lg, lg)

# dotprod(v_v::Vector{Vector{T}}) where{T}=dotprod(totuple(v_v))
# modone(x,n)=(x-1)%n+1
# struct CircularArray{T,N} <: AbstractArray{T,N} #inherits from AbstractArray
#     x::AbstractArray{T,N}
#     tabmod::NTuple{N,Vector{Int}}
#     function CircularArray(x::AbstractArray{T,N}) where {T,N}
# #creates the type only with a vector x
#         sz = size(x)
#         tabmod = ntuple( i -> modone.(1:10*sz[i], sz[i]), N)
#         return new{T,N}(x, tabmod)
#     end
# end
# modtabmod(tm::AbstractVector{Int}, i::Int)=tm[i]

# Base.size(A::CircularArray) = 10 .* size(A.x) #important: use of Base.function

# Base.length(A::CircularArray)=length(A.x)

# function Base.getindex(A::CircularArray, I::Vararg{Int, N}) where N # implements A[I]
#     return Base.getindex(A.x,modtabmod.(A.tabmod, I)...) #this is the magic operation
# end
# function Base.getindex(A::CircularArray, I::Tuple{N, Int}) where N # implements A[I]
#      return Base.getindex(A.x,modtabmod.(A.tabmod, I)...) #this is the magic operation
# end

# Base.getindex(A::CircularArray{T,nd}, CI::CartesianIndex{nd}) where{T,nd}=getindex(A, CI.I)

# Base.getindex(A::CircularArray, I) = (A[i] for i in I) #A[1:5], for example

# function Base.setindex!(A::CircularArray,value,I::Vararg{Int, N}) where N # A[I] = value
#     I2 = size(A)
#     return Base.setindex!(A.x,value,modtabmod.(A.tabmod,I)...)
# end
# function Base.setindex!(A::CircularArray,value,I::Tuple{N, Int}) where N # A[I] = value
#     I2 = size(A)
#     return Base.setindex!(A.x,value,modtabmod.(A.tabmod,I)...)
# end
# Base.setindex!(A::CircularArray,value, CI::CartesianIndex)= setindex!(A,value,CI.I)

# Base.IndexStyle(::Type{CircularArray}) = IndexCartesian()
# a finir

function getextarray(
    tabor::Array{T2,N},
    decbeg::NTuple{N,Int},
    decend::NTuple{N,Int},
) where {T2,N}
    tab = zeros(T2, size(tabor) .+ decbeg .+ decend)
    inddecbeg = CartesianIndex(decbeg)
    sz = size(tabor)
    for ind in CartesianIndices(tab)
        tab[ind] = tabor[modone(ind - inddecbeg, sz)]
    end
    return tab
end

# struct ExtArray{T,N} <: AbstractArray{T,N}
#     tab::Array{T,N}
#     inddecbegin::CartesianIndex{N}
#     function ExtArray(
#         tabor::Array{T,N},
#         decbeg::NTuple{N,Int},
#         decend::NTuple{N,Int},
#     ) where {T,N}
#         tab = zeros(T, size(tabor) .+ decbeg .+ decend)
#         inddecbeg = CartesianIndex(decbeg)
#         sz = size(tabor)
#         for ind in CartesianIndices(tab)
#             tab[ind] = tabor[modone(ind + inddecbeg, sz)]
#         end
#         return new{T,N}(tab, inddecbeg)
#     end
# end
# Base.getindex(A::ExtArray{T,N}, ind::CartesianIndex{N}) where {T,N} =
#     A.tab[ind+A.inddecbegin]
# Base.axes(A::ExtArray{T,N}) where {T,N} =
#     ntuple(x -> (-A.inddecbegin[x]+1):(size(A.tab, x)-A.inddecbegin[x]), N)
# Base.size(A::ExtArray{T,N}) where {T,N} = size(A.tab)
# getarray(A::ExtArray) = A.tab


text:
function test_zeros(T=Float64, dims=(13, 45))
  a1, b1 = zeros(T, dims), zeros(T, dims)
  OceanTurb.@zeros T dims a2 b2
  eltype(a1) == T && a1 == a2 && b1 == b2
end

function test_run_until(stepper=:ForwardEuler)
    model = Diffusion.Model(N=3, H=1, K=1.0, stepper=stepper)

    dt, tfinal = 0.5, 1.3
    run_until!(model, dt, tfinal)

    return time(model) == tfinal
end

function test_diffusive_flux()
    model = Diffusion.Model(N=4, H=2, K=0.1)
    c0(z) = z
    model.solution.c = c0
    flux = OceanTurb.diffusive_flux(:c, model)
    return flux[2] == -0.1
end

@testset "Utils" begin
    @test test_zeros(Float64)
    @test test_zeros(Float32)
    @test test_diffusive_flux()
    for stepper in steppers
        @test test_run_until(stepper)
    end
end


text:
Cluster = Tuple{Int64,Int64,Int64}

function lt_cl(x::Cluster, y::Cluster)
    return     (x[1] < y[1]) ||
             ( (x[1] == y[1]) && (x[2] < y[2]) ) ||
             ( (x[1] == y[1]) && (x[2] == y[2]) && (x[3] < y[3]) )
end

mutable struct NDSep # A nested dissection separator, "s"

    ## Hierarchy & dofs data
    ptr::Array{Int64,1}             # A pointer array to the end of each cluster within s at the current stage of the elimination
    hrch::IT                        # A hierarchy, i.e., a tree of integer, monotone left -> right
    start::Int64                    # Where does the separator start
    depth::Int64                    # Where the clustering stand in terms of depth

    ## Data
    Piv::Matrix{Float64}            # Ass: Diagonal
    Low::Matrix{Float64}            # Ans: Lower-part
    Upp::Matrix{Float64}            # Asn: Upper-part

    ## Block information, basically pointers to the data above
    APiv::Matrix{StridedMatrix{Float64}} # APiv[i,j] = ith partition of s, jth partition of s - (self lvl, self sep, i)
    ALow::Matrix{StridedMatrix{Float64}} # ALow[n,i] = ith partition of s, nth neighbor
    AUpp::Matrix{StridedMatrix{Float64}} # AUpp[i,n] = ith partition of s, nth neighbor
    Nbr::Vector{Cluster}                 # Nbr[n]    = nth neighbor (lvl, sep, part)
    # When fully merged, size(APiv) == (1,1), size(ALow,2) == size(AUpp,1) == 1

    ## Range information, for easier backward/forward pass
    NbrRanges::Vector{UnitRange{Int64}}

    function NDSep(start, subs, hrch)
        if length(subs) >= 1
            ptr = sorted_to_ptr(subs)
        else
            ptr = [1, 1]
        end
        this = new()
        this.start = start
        @assert length(ptr) >= 2
        this.ptr = ptr
        this.hrch = hrch
        this.depth = length(this.hrch.nit)
        this.NbrRanges = Vector{UnitRange{Int64}}()
        this
    end
end

struct NDTree
    t::Vector{Vector{NDSep}}
    p::Vector{Int64}
end

function get_dofs(self::NDSep, i::Int64)
    return (self.start-1) .+ (self.ptr[i]:self.ptr[i+1]-1)
end

function get_dofs(self::NDSep)
    return (self.start-1) .+ (self.ptr[1]:self.ptr[end]-1)
end

function get_size(self::NDSep)
    return (self.ptr[end] - self.ptr[1])
end

function get_size(self::NDSep, i::Int64)
    return (self.ptr[i+1] - self.ptr[i])
end

function nparts(self::NDSep)
    np = length(self.ptr)-1;
    @assert length(self.hrch.nit[self.depth]) == np
    return np
end

function nparts_merged(self::NDSep)
    return length(self.hrch.nit[self.depth-1])
end

function ptr_merged(self::NDSep)
    ptr = self.ptr
    ptr_new = ones(Int64, nparts_merged(self)+1);
    for n in self.hrch.nit[self.depth-1]
        ptr_new[n.i+1] = ptr[n.c[end].i+1]
    end
    return ptr_new
end

# Given i, go up l level in the tree
# Return the range low:high (i \in low:high)
# corresponding to i's ancestor at i's level
function go_up_down(self::NDSep, i::Int64, l::Int64)
    @assert l >= 0
    if l == 0
        return i:i
    end
    # Go up
    h = self.hrch.nit
    p = h[end][i]
    for _ in 1:l
        p = p.p
    end
    # Go left
    left = p
    for _ in 1:l
        left = left.c[1]
    end
    # Go right
    right = p
    for _ in 1:l
        right = right.c[end]
    end
    return left.i:right.i
end

function parent(self::NDSep, i::Int64)
    self.hrch.nit[self.depth][i].p.i
end

function siblings(self::NDSep, i::Int64)
    p = self.hrch.nit[self.depth][i].p
    return (p.c[1].i:p.c[end].i, p.i)
end

function factorize(A::SparseMatrixCSC{Float64,Int64}, maxLevel::Int64, verbose::Bool=false)

    N = size(A, 1)
    @assert size(A, 2) == N
    @assert maxLevel >= 1

    # Symmetric partition
    (subseps, hrch, dofs) = mnd(A+A', maxLevel, verbose=verbose);

    ## Permute the matrix & initialize an empty tree
    perm = Vector{Int64}(undef, N)
    Tree = Vector{Vector{NDSep}}(undef, maxLevel)
    start = 1
    for lvl = 1:maxLevel
        nsep = length(dofs[lvl])
        Tree[lvl] = Vector{NDSep}(undef, nsep)
        for sep = 1:nsep
            dofs_ = dofs[lvl][sep]
            subs_ = subseps[dofs_]
            perm_ = sortperm(subs_)
            perm[start:(start+length(dofs_)-1)] .= dofs_[perm_]
            Tree[lvl][sep] = NDSep(start, subs_[perm_], hrch[lvl][sep])
            start += length(dofs_)
        end
    end
    A = A[perm,perm]

    ## Create dofs2cl map
    dofs2cl = Vector{Cluster}(undef, N)
    for (lvl, sep) in BinTreeIt(maxLevel)
        for i in dofs[lvl][sep]
            dofs2cl[i] = (lvl, sep, subseps[i])
        end
    end
    dofs2cl = dofs2cl[perm]

    ## Symbolic computation: Compute sparsity pattern, allocate & assemble
    if(verbose) @printf("Assembling...\n") end
    for (lvl, sep) in BinTreeIt(maxLevel)
        h = Tree[lvl][sep].hrch
        s = length(dofs[lvl][sep])
        np = length(h.nit[end])
        # Get the original neighbors
        (I, t1, t2) = findnz(A[:,get_dofs(Tree[lvl][sep])])
        (t1, J, t2) = findnz(A[get_dofs(Tree[lvl][sep]),:]) # FIXME: probably slow since in CSC
        nbrs = sort(unique(vcat(I,J)))
        Nbr = Set(dofs2cl[nbrs])
        # Add the children
        if lvl > 1
            union!(Nbr, Tree[lvl-1][2*sep-1].Nbr)
            union!(Nbr, Tree[lvl-1][2*sep  ].Nbr)
        end
        # Remove too small
        filter!(x -> x[1] > lvl, Nbr)
        # Subsample
        Nbr2 = copy(Nbr)
        for (l_, s_, i_) in Nbr2
            lr = go_up_down(Tree[l_][s_], i_, lvl - 1)
            for i2_ in lr
                push!(Nbr, (l_, s_, i2_))
            end
        end
        # Clean & sort
        Nbr = collect(Nbr);
        sort!(Nbr, lt=lt_cl)
        Tree[lvl][sep].Nbr = Nbr
        # Allocate full data
        nsizes = Vector{Int64}(map(x -> Tree[x[1]][x[2]].ptr[x[3]+1]-Tree[x[1]][x[2]].ptr[x[3]], Nbr))
        n = sum(nsizes)
        Low = zeros(n, s)
        Upp = zeros(s, n)
        Piv = zeros(s, s)
        Tree[lvl][sep].Low = Low
        Tree[lvl][sep].Upp = Upp
        Tree[lvl][sep].Piv = Piv
        # Get subblocks
        ALow = Matrix{StridedMatrix{Float64}}(undef, (length(Nbr), np))
        AUpp = Matrix{StridedMatrix{Float64}}(undef, (np, length(Nbr)))
        APiv = Matrix{StridedMatrix{Float64}}(undef, (np, np))
        sptr = Tree[lvl][sep].ptr
        nptr = sizes_to_ptr(nsizes)
        # Fill diag
        for i in 1:np
            idofs = get_dofs(Tree[lvl][sep], i)
            iids  = sptr[i]:(sptr[i+1]-1)
            # Diagonal
            for j in 1:np
                jdofs = get_dofs(Tree[lvl][sep], j)
                jids  = sptr[j]:(sptr[j+1]-1)
                APiv[i,j] = view(Piv, iids, jids)
                APiv[i,j][:,:] = A[idofs,jdofs]
            end
            # Lower/Upper
            for n in 1:length(Nbr)
                (lvln,sepn,in) = Nbr[n]
                ndofs = get_dofs(Tree[lvln][sepn], in)
                nids  = nptr[n]:(nptr[n+1]-1)
                ALow[n,i] = view(Low, nids, iids);
                ALow[n,i][:,:] = A[ndofs, idofs]
                AUpp[i,n] = view(Upp, iids, nids)
                AUpp[i,n][:,:] = A[idofs, ndofs]
            end
        end
        Tree[lvl][sep].ALow = ALow
        Tree[lvl][sep].AUpp = AUpp
        Tree[lvl][sep].APiv = APiv
    end

    ## Numerical factorization
    for lvl = 1:maxLevel

        if(verbose) @printf("Lvl: %d, eliminating\n", lvl) end

        ## Eliminate
        for sep = 1:length(dofs[lvl])
            s = Tree[lvl][sep]
            # Check it's fully merged
            @assert all(size(s.APiv) .== (1,1))
            # GETRF, factor pivot (no pivoting so far)
            Ass = s.APiv[1,1]
            lu!(Ass, Val(false)); # Ass -> L U
            # TRSM, solve panel in place
            for n = 1:length(s.Nbr)
                (ln,sn,in) = s.Nbr[n]
                push!(s.NbrRanges, get_dofs(Tree[ln][sn],in))
                # Lower
                Ans = s.ALow[n,1]
                BLAS.trsm!('R','U','N','N', 1.0, Ass, Ans) # Ans -> Ans U^-1
                # Upper
                Asn = s.AUpp[1,n]
                BLAS.trsm!('L','L','N','U', 1.0, Ass, Asn) # Asn -> L^-1 Asn
            end
            # GEMM, update schur complement
            for ni = 1:length(s.Nbr)
                for nj = 1:length(s.Nbr)
                    Ais = s.ALow[ni,1]
                    Asj = s.AUpp[1,nj]
                    (li,si,ii) = s.Nbr[ni]
                    (lj,sj,jj) = s.Nbr[nj]
                    if (li,si) == (lj,sj)
                        Aij = Tree[li][si].APiv[ii,jj]
                    else
                        @assert li != lj # Can't be any edge between different separators at a given level
                        if li < lj
                            n = Tree[li][si]
                            id = searchsortedfirst(n.Nbr, s.Nbr[nj], lt=lt_cl)
                            Aij = n.AUpp[ii,id]
                        else
                            n = Tree[lj][sj]
                            id = searchsortedfirst(n.Nbr, s.Nbr[ni], lt=lt_cl)
                            Aij = n.ALow[id,jj]
                        end
                    end
                    BLAS.gemm!('N', 'N', - 1.0, Ais, Asj, 1.0, Aij) # Aij -> Aij - Ais Asj
                end
            end
        end

        ## Merge, i.e., redefine all the views
        for lvl2 = lvl+1:maxLevel
            for sep2 = 1:length(dofs[lvl2])
                s = Tree[lvl2][sep2]
                ptr2 = ptr_merged(s)
                # Pivot
                for i = 1:length(ptr2)-1
                    for j = 1:length(ptr2)-1
                        s.APiv[i,j] = view(s.Piv, ptr2[i]:(ptr2[i+1]-1), ptr2[j]:(ptr2[j+1]-1))
                    end
                end
                s.APiv = s.APiv[1:length(ptr2)-1, 1:length(ptr2)-1]
                # Lower/Upper
                n_i1, n_b1, n_merged = 1, 1, 1
                while n_b1 <= length(s.Nbr)
                    # What are n's siblings ?
                    (l_, s_, p_) = s.Nbr[n_b1]
                    (n_siblg, parent) = siblings(Tree[l_][s_], p_)
                    # All the siblings should follow
                    n_i2, n_b2 = n_i1, n_b1
                    for i = 1:length(n_siblg)
                        n_p_  = n_siblg[i]
                        @assert s.Nbr[n_b2] == (l_, s_, n_p_)
                        n_b2 += 1
                        n_i2 += get_size(Tree[l_][s_], n_p_)
                    end
                    # Update the row in Low and in Nbr
                    s.Nbr[n_merged] = (l_, s_, parent)
                    for c = 1:length(ptr2)-1
                        s.ALow[n_merged, c] = view(s.Low, n_i1:(n_i2-1), ptr2[c]:(ptr2[c+1]-1))
                        s.AUpp[c, n_merged] = view(s.Upp, ptr2[c]:(ptr2[c+1]-1), n_i1:(n_i2-1))
                    end
                    # Next
                    n_merged += 1
                    n_i1, n_b1 = n_i2, n_b2
                end
                # Resize
                s.ALow = s.ALow[1:n_merged-1, 1:length(ptr2)-1]
                s.AUpp = s.AUpp[1:length(ptr2)-1, 1:n_merged-1]
                s.Nbr  = s.Nbr[1:n_merged-1]
            end
        end

        # Update ptr & depth
        for lvl2 = lvl+1:maxLevel
            for sep2 = 1:length(dofs[lvl2])
                s = Tree[lvl2][sep2]
                s.ptr = ptr_merged(s)
                s.depth -= 1
            end
        end
    end

    return NDTree(Tree, perm)

end


text:
module DocumentationModule
export Citation, JournalCitation, BookCitation, GenericCitation, summarize
import Base.show
import Base.convert
import Base.promote_rule
using Markdown
using InteractiveUtils

"""Abstract class supertyping all objects that are entries of a bibliography (for documenting functions). Each subtype has its corresponding fields and defines its own show function."""
abstract type Citation end

"""
Documentation is a generic type to which all types which are to be used with `@doc` can be automatically converted.
At the moment this means: Any type of Citation and Markdown
"""
struct Documentation
  doc::Union{Citation,Markdown.MD}
end

"""Citation class for articles in scientific jounals."""
mutable struct JournalCitation <: Citation
  authors::Array{AbstractString,1}
  journalabbr::AbstractString
  vol::UInt
  page::UInt
  year::UInt
end

function show(io::IO,mime::MIME"text/plain",cite::JournalCitation)
  show(io,mime,Markdown.parse(""" - $(join(cite.authors,", ")), $(cite.journalabbr) **$(cite.vol)**, *$(cite.page)* ($(cite.year))"""))
  print(io,"\n")
end

"""Citation class for complete Books. Refer to BookSectionCitation for more specific parts of a book."""
mutable struct BookCitation <: Citation
  authors::Array{AbstractString,1}
  title::AbstractString
  ISBN::AbstractString
end

function show(io::IO,mime::MIME"text/plain",cite::BookCitation)
  show(io,mime,Markdown.parse(""" - $(join(cite.authors,", ")): *$(cite.title)* (ISBN $(cite.ISBN))"""))
  print(io,"\n")
end

"""Citation class for anything that can not be described by any of the other subclasses of Citation. Make sure none of the others fits the purpose before using this one."""
struct GenericCitation <: Citation
  string::AbstractString
end

function show(io::IO,mime::MIME"text/plain",cite::GenericCitation)
  show(io,mime,Markdown.parse(""" - $(cite.string)"""))
  print(io,"\n")
end

function show(io::IO,mime::MIME"text/plain",docuArray::Array{Documentation,1})
  local cites = Array{Citation,1}()
  show(io,mime,Markdown.parse("""
  # Documentation
  """))
  print(io,"\n")
  for document in docuArray
    docu = document.doc
    if typeof(docu)==Markdown.MD
      show(io,mime,docu)
    elseif isa(docu,Citation)
      push!(cites,docu)
    end
  end

  print(io,"\n\n")

  show(io,mime,Markdown.parse("""
  # Citations:
  """))
  print(io,"\n")
  for citekey in cites
    show(io,mime,citekey)
  end
end

function show(io::IO,mime::MIME"text/plain",docuArray::Array{Citation,1})
  show(io,mime,Markdown.parse("""
  # Citations:
  """))
  print(io,"\n")
  for citekey in docuArray
    show(io,mime,citekey)
  end
end

function show(io::IO,mime::MIME"text/plain",docuArray::Array{T,1}) where T<:Citation
  show(io,mime,Markdown.parse("""
  # Citations:
  """))
  print(io,"\n")
  for citekey in docuArray
    show(io,mime,citekey)
  end
end

"""
Taken from base/docs/Docs.jl, this function shows a summary of the fields, type and supertypes of a DataType.
In Docs.jl the function `typesummary` is only meant to be called when no other documentation can be found for it.
If `typesummary` doesn't print "No documentation found" anymore in the future, then this function becomes obsolete.
"""
function summarize(f::DataType)
    parts = [
    """
    **Summary:**
    ```julia
    $(f.abstract ? "abstract" : f.mutable ? "type" : "immutable") $f <: $(supertype(f))
    ```
    """
    ]
    if !f.abstract && !isempty(fieldnames(f))
        pad    = maximum([length(string(f)) for f in fieldnames(f)])
        fields = ["$(rpad(f, pad)) :: $(t)" for (f, t) in zip(fieldnames(f), f.types)]
        push!(parts,
        """
        **Fields:**
        ```julia
        $(join(fields, "\n"))
        ```
        """)
    end
    if !isempty(subtypes(f))
        push!(parts,
        """
        **Subtypes:**
        ```julia
        $(join(subtypes(f), "\n"))
        ```
        """)
    end
    Markdown.parse(join(parts, "\n"))
end

# promotion and conversion
Base.promote_rule(::Type{String},::Type{Documentation})         = Documentation
Base.promote_rule(::Type{String},::Type{Citation})         	= Documentation
Base.promote_rule(::Type{String},::Type{GenericCitation}) 	= Documentation
Base.promote_rule(::Type{String},::Type{BookCitation}) 		= Documentation
Base.promote_rule(::Type{String},::Type{JournalCitation}) 	= Documentation
Base.promote_rule(::Type{Markdown.MD},::Type{Documentation}) 	= Documentation
Base.promote_rule(::Type{Markdown.MD},::Type{Citation})         = Documentation
Base.promote_rule(::Type{Markdown.MD},::Type{GenericCitation}) 	= Documentation
Base.promote_rule(::Type{Markdown.MD},::Type{BookCitation}) 	= Documentation
Base.promote_rule(::Type{Markdown.MD},::Type{JournalCitation}) 	= Documentation
Base.promote_rule(::Type{Documentation},::Type{GenericCitation}) 	= Documentation
Base.promote_rule(::Type{Documentation},::Type{BookCitation}) 		= Documentation
Base.promote_rule(::Type{Documentation},::Type{JournalCitation}) 	= Documentation
Base.convert(::Type{Documentation},str::String) = Documentation(Markdown.parse(str))
Base.convert(::Type{Documentation},md::Markdown.MD) = Documentation(md)
Base.convert(::Type{Documentation},cite::Citation) = Documentation(cite)

macro ignore_warnings(expr)
  quote
    STDERRorig = stderr
    drop = redirect_stderr()
    result = $expr
    redirect_stderr(STDERRorig)
    return result
  end
end

end # module


text:
# I map some model ids to the enzimatic cost data reported in
# Beg2007, so, I will set the cost of this reactions using this data
# The rest of the model reactions will get the average enzimatic cost also
# reported by Beg.
const _beg_enz_cost_dict = Dict()
function _load_beg_enz_cost_dict()
    empty!(_beg_enz_cost_dict)
    beg_rxns_map = load_beg_rxns_map()
    beg_enz_data = Bd.load_enz_data()
    for (i, beg_rxn) in enumerate(beg_enz_data.Enzyme)
        for model_rnx in beg_rxns_map[beg_rxn]
            if model_rnx != ""
                _beg_enz_cost_dict[model_rnx] =
                    beg_enz_data[i, Symbol("ai/xi (h g / mmol)")]
            end
        end
    end
end

function beg_enz_cost(model_rxn)
    isempty(_beg_enz_cost_dict) && _load_beg_enz_cost_dict()
    return get(_beg_enz_cost_dict, model_rxn, Bd.ave_a)
end


text:
@reexport using Distributions
import Distributions: rand
export MultivariateMix, NonMixedMultivariateDistribution, cov, cor, simulate_nobs_independent_vectors, covariance_matrix, correlation_matrix
# ## create new type not in Distributions for multivariate mixture

"""
    MultivariateMix(vecd::AbstractVector{<:UnivariateDistribution}, Γ)
An N dimensional `MultivariateMix` constructed from a vector of N `UnivariateDistribution`s, and the user-specified covariance matrix, Γ.
Since Distributions.jl does not currently allow for multivariate vectors with mixed continuous and discrete components, this type is not of super type Distributions.MultivariateDistribution.
We will pre-allocate a Vector of `ContinuousUnivariateCopula`s and `DiscreteUnivariateCopula`s to then fill in the appropriate constants c0, c1, c2, recursively.
"""
struct MultivariateMix{
    V<:AbstractVector{UnivariateDistribution},
    T<: BlasReal
}
    vecd::V
    Γ::Matrix{T}
    trΓ::Float64
    gc_obs::Vector{Union{ContinuousUnivariateCopula, DiscreteUnivariateCopula}}
    function MultivariateMix(vecd::V, Γ::Matrix{T}) where {T <: BlasReal,
        V<:AbstractVector{UnivariateDistribution}}
        n = length(vecd)
        trΓ = tr(Γ)
        gc_obs = Vector{Union{ContinuousUnivariateCopula, DiscreteUnivariateCopula}}(undef, n)
        return new{V, T}(vecd, Γ, trΓ, gc_obs)
    end
end

"""
    NonMixedMultivariateDistribution(vecd::AbstractVector{<:UnivariateDistribution}, Γ)
An N dimensional `MultivariateDistribution` constructed from a vector of N `UnivariateDistribution`s, and the user-specified covariance matrix, Γ.
We will pre-allocate a Vector of `ContinuousUnivariateCopula`s or `DiscreteUnivariateCopula`s to then fill in the appropriate constants c0, c1, c2, recursively.
"""
struct NonMixedMultivariateDistribution{
    S<:ValueSupport,
    T<:UnivariateDistribution{S},
    V<:AbstractVector{T},
}
    vecd::V
    Γ::Matrix{Float64}
    trΓ::Float64
    gc_obs::Vector{Union{DiscreteUnivariateCopula, ContinuousUnivariateCopula}}
    function NonMixedMultivariateDistribution(vecd::V, Γ::Matrix{Float64}) where
        V<:AbstractVector{T} where
        T<:UnivariateDistribution{S} where
        S<:ValueSupport
        n = length(vecd)
        trΓ = tr(Γ)
        gc_obs = Vector{Union{ContinuousUnivariateCopula, DiscreteUnivariateCopula}}(undef, n)
        return new{S, T, V}(vecd, Γ, trΓ, gc_obs)
    end
end

####
function rand(gc_vec::NonMixedMultivariateDistribution{S, T, V},
    Y::Vector{Float64},
    res::Vector{Float64}) where {S<: ValueSupport, T <: UnivariateDistribution{S}, V<:AbstractVector{T}}
    for i in 1:length(gc_vec.vecd)
        # form constants for conditional density of i given 1, ..., i-1
        gc_vec.gc_obs[i] = pdf_constants(gc_vec.Γ, res, i, gc_vec.vecd[i])
        # generate y_i given y_1, ..., y_i-1
        Y[i] = rand(gc_vec.gc_obs[i])
        # update residuals 1,..., i-1
        res[i] = update_res!(Y[i], res[i], gc_vec.gc_obs[i])
    end
    res[end] = update_res!(Y[end], res[end], gc_vec.gc_obs[end])
    Y
end

# ### we need a function that is like rand, but cannot use the distributions framework.
function rand(gc_vec::MultivariateMix{V, T},
    Y::Vector{T},
    res::Vector{T}) where {T <: BlasReal, V<:AbstractVector{UnivariateDistribution}}
    for i in 1:length(gc_vec.vecd)
        # form constants for conditional density of i given 1, ..., i-1
        gc_vec.gc_obs[i] = pdf_constants(gc_vec.Γ, res, i, gc_vec.vecd[i])
        # generate y_i given y_1, ..., y_i-1
        Y[i] = rand(gc_vec.gc_obs[i])
        # update residuals 1,..., i-1
        res[i] = update_res!(Y[i], res[i], gc_vec.gc_obs[i])
    end
    res[end] = update_res!(Y[end], res[end], gc_vec.gc_obs[end])
    Y
end

"""
    simulate_nobs_independent_vectors(multivariate_distribution::Union{NonMixedMultivariateDistribution, MultivariateMix}, n_obs::Integer)
Simulate n_obs independent realizations from multivariate copula-like distribution.
"""
function simulate_nobs_independent_vectors(
    multivariate_distribution::Union{NonMixedMultivariateDistribution, MultivariateMix},
    n_obs::Integer)
    dimension = length(multivariate_distribution.vecd)
    Y = [Vector{Float64}(undef, dimension) for i in 1:n_obs]
    res = [Vector{Float64}(undef, dimension) for i in 1:n_obs]
    for i in 1:n_obs
        rand(multivariate_distribution, Y[i], res[i])
    end
    Y
end

"""
    logpdf(d::Union{ContinuousUnivariateCopula, DiscreteUnivariateCopula}, x::Real)
Theoretical log pdf under the copula model.
"""
function logpdf(multivar::Union{MultivariateMix, NonMixedMultivariateDistribution}, Y::Vector{T}) where T <: BlasReal
    logl = zero(T)
    for i in 1:length(multivar.vecd)
        logl += logpdf(multivar.gc_obs[i], Y[i])
    end
    logl
end

"""
    cov(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix}, k::Int64, l::Int64)
Theoretical covariance of the kth and lth element in the random vector.
"""
function cov(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix}, k::Int64, l::Int64)
    cov = sqrt(var(gc_vec.vecd[k])) * sqrt(var(gc_vec.vecd[l])) * gc_vec.Γ[k, l] * inv(1 + 0.5 * gc_vec.trΓ)
    return cov
end

"""
    cor(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix}, k::Int64, l::Int64)
Theoretical correlation of the kth and lth element in the random vector.
"""
function cor(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix}, k::Int64, l::Int64)
    cor = gc_vec.Γ[k, l] / (sqrt(1 + 0.5 * gc_vec.trΓ + gc_vec.Γ[k, k]) * sqrt(1 + 0.5 * gc_vec.trΓ + gc_vec.Γ[l, l]) )
    return cor
end

"""
    covariance_matrix(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix})
Theoretical covariance matrix of a random vector.
"""
function covariance_matrix(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix})
    n = length(gc_vec.gc_obs)
    Covariance = zeros(n, n)
    for i in 1:n
        Covariance[i, i] = QuasiCopula.var(gc_vec.gc_obs[i])
        for j = i+1:n
            Covariance[j, i] = QuasiCopula.cov(gc_vec, j, i)
            Covariance[i, j] = Covariance[j, i]
        end
    end
    Covariance
end

"""
    correlation_matrix(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix})
Theoretical covariance matrix of a random vector.
"""
function correlation_matrix(gc_vec::Union{NonMixedMultivariateDistribution, MultivariateMix})
    n = length(gc_vec.gc_obs)
    Corr = zeros(n, n)
    for i in 1:n
        Corr[i, i] = 1.0
        for j = i+1:n
            Corr[j, i] = QuasiCopula.cor(gc_vec, j, i)
            Corr[i, j] = Corr[j, i]
        end
    end
    Corr
end


text:
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#
# It is included in file Modia3D/Composition/sceneOptions.jl
# in order to be used as default renderer in SceneOptions(..)
#

mutable struct DummyRenderer <: Modia3D.AbstractRenderer
  DummyRenderer(info; host="127.0.0.1",port=11000,sync=false) = new()
end


text:
using InstrumentOperator
using NCDatasets
using OrderedCollections
using YAML

# Load L1 file (replace with your example!):
L1File = "/net/fluo/data1/group/oco2/L1bSc/oco2_L1bScND_26780a_190715_B10003r_200429212407.h5"
ocoData = Dataset(L1File);

# Load corresponding met file:
metFile = "/net/fluo/data1/group/oco2/L2Met/oco2_L2MetND_26780a_190715_B10003r_200429212406.h5"
metData = Dataset(metFile);

# Load dictionary:
dictOCO2 = YAML.load_file("json/oco2.yaml"; dicttype=OrderedDict{String,Any});

# Load L1 file (could just use filenames here as well)
oco = InstrumentOperator.load_L1(dictOCO2,ocoData, metData);

# Pick some bands as tuple (or just one)
bands = (1,2);
# Indices within that band:
indices = (1:1016,20:1000);
# Geo Index (footprint,sounding):
GeoInd = [5,3000];

# Get data for that sounding:
oco_sounding = InstrumentOperator.getMeasurement(oco, bands, indices, GeoInd)


text:
### A Pluto.jl notebook ###
# v0.17.2

# using Markdown
# using InteractiveUtils

# ╔═╡ bb26cad0-4720-11ec-20c0-6538c8fcc18f
#=╠═╡ notebook_exclusive
begin
	import Pkg
	Pkg.activate("..")
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ c0e30b73-0973-47c8-b57c-3a0b9c534ec5
using BenchmarkTools

# ╔═╡ 82ace475-cc8d-4f3d-84bd-91a49c141648
using PlutoDevMacros

# ╔═╡ 8c251cd4-3050-4809-80a3-07e9aab558ab
struct TestStruct1 end

# ╔═╡ 8145c16e-0a16-4e8b-a073-8e1774273bb4
asd(x::Int) = "INTSSS"

# ╔═╡ 354605a3-297e-4647-b4cf-f0aea4987f9e
"""
	asd(::TestStruct1)
This is a magical docstring
"""
asd(::TestStruct1) = "TESTSTRUCT1"

# ╔═╡ 9e1468e5-beda-486a-9cde-9c498f77999e
@benchmark asd(TestStruct1())

# ╔═╡ Cell order:
# ╠═bb26cad0-4720-11ec-20c0-6538c8fcc18f
# ╠═c0e30b73-0973-47c8-b57c-3a0b9c534ec5
# ╠═82ace475-cc8d-4f3d-84bd-91a49c141648
# ╠═8c251cd4-3050-4809-80a3-07e9aab558ab
# ╠═8145c16e-0a16-4e8b-a073-8e1774273bb4
# ╠═354605a3-297e-4647-b4cf-f0aea4987f9e
# ╠═9e1468e5-beda-486a-9cde-9c498f77999e


text:
module FrostHelperCassetteTempoTrigger

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/CassetteTempoTrigger" FHCassetteTempoTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16, Tempo::Number=1.0, ResetOnLeave::Bool=false)

const placements = Ahorn.PlacementDict(
	"Cassette Tempo Trigger (Frost Helper)" => Ahorn.EntityPlacement(
		FHCassetteTempoTrigger,
		"rectangle"
	)
)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
# Unit Test for `SolveQuadraticProgram()`

using LinearAlgebra;
using SparseArrays;
using Plots;
using Random;
using StableRNGs;
using MAT;
using BenchmarkTools;

seedNumber = 1234;
Random.seed!(seedNumber);
# Random.seed!(StableRNG(seedNumber), seedNumber);

include("GenerateQuadraticProgram.jl");
include("SolveQuadraticProgram.jl");
include("LinearSystemSolvers.jl");

# @enum ProblemClass randomQp = 1 inequalityConstrainedQp equalityConstrainedQp optimalControl portfolioOptimization lassoOptimization huberFitting supportVectorMachine isotonicRegression
@enum DataSource dataSourceGenerated = 1 dataSourceLoaded

## Parameters

# Simulaion
numElements     = 800;
numConstraints  = 400;

dataSource      = dataSourceGenerated;
dataFileName    = "QpModel.mat";

# problemClass    = rand(instances(ProblemClass));
problemClass    = randomQp;

# Solver
numIterations   = 7500;
ρ               = 1e6
adptΡ           = true;
linSolverMode   = modeItertaive;
numItrPolish    = 0;

if (dataSource == dataSourceGenerated)
    mP, vQ, mA, vL, vU = GenerateRandomQP(problemClass, numElements, numConstraints = numConstraints);
else
    dMatFile = matread(dataFileName);
    mP = dMatFile["mP"];
    vQ = dropdims(dMatFile["vQ"]; dims = 2);
    mA = dMatFile["mA"];
    vL = dropdims(dMatFile["vL"]; dims = 2);
    vU = dropdims(dMatFile["vU"]; dims = 2);
end

# Some problems might change the actual data
numElements     = size(mP, 1);
numConstraints  = size(mA, 1);

vXX = zeros(numElements);

hLinSolInit = QDLdlInit;
hLinSol     = QDLdl!;

sBenchMark = @benchmarkable SolveQuadraticProgram!(vX, $mP, $vQ, $mA, $vL, $vU, hLinSolInit, hLinSol; numIterations = $numIterations, ρ = $ρ, adptΡ = $adptΡ, numItrPolish = $numItrPolish) setup = (vX = copy($vXX));
currTime    = time();
tuRunResult = BenchmarkTools.run_result(sBenchMark, samples = 10, evals = 1, seconds = 150); #<! A trick to get result as well
runTime     = time() - currTime;

display(tuRunResult[1]);
println("\nTotal Run Time: $runTime [Sec]");
println("The Alogrihtm Convergence Status: $(tuRunResult[2])\n\n");

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
" Saves region data.
# Arguments
- `region_dict::Dict` - Dictionary with region extents
- `region_post_dict::Dict` - Dictionary with region image overapproximations
- `filename::String` - Location to save region data
"
function save_region_data(region_dict::Dict, region_post_dict::Dict, filename::String; lazy_flag=false)
	region_data = serialize_region_data(region_dict, region_post_dict, lazy_flag=lazy_flag)
	bson(filename, Dict(:res => region_data))
end

" Saves region data to the experiment root directory.
# Arguments
- `exp_dir::String` - Experiment root directory
- `region_data::Dict` - Dictionary containing region extent, post dictionaries
"
function save_region_data(exp_dir::String, region_data::Dict)
    @info "Saving the region data to the experiment directory..."
    region_filename = "$exp_dir/regions.bson"
    save_region_data(region_data[:extents], region_data[:posts], region_filename)
end

" Loads region data.
# Arguments
- `filename::String` - Region data file to load
"
function load_region_data(filename::String)
	region_dict, region_post_dict, _ = deserialize_region_data(filename)
    region_data = Dict(:extents => region_dict, :posts => region_post_dict)
	return region_data
end

function create_region_data(space::Dict, grid_sizes::Dict)
    region_dict = discretize_set(space, grid_sizes)
    # ? THIS IS NOT CORRECT for a 3d model ? Not sure what I mean by this comment
    region_dict[region_dict.count+1] = space # Region with index L+1 corresponds to the safe set as a whole.
    return region_dict
end

function create_region_data_new(space::Dict, grid_size::Dict)
    region_dict, extent_dict = discretize_set_lazy(space, grid_size)
    x = space["x1"]
    y = space["x2"]
    region_dict[region_dict.count+1] = VPolygon([[x[1], y[1]], [x[1], y[2]],
                                    [x[2], y[2]], [x[2], y[1]]])
    extent_dict[region_dict.count+1] = space

    return region_dict, extent_dict
end

function serialize_region_data(region_dict::Dict, region_post_dict::Dict; lazy_flag=false)

	if lazy_flag
		return serialize_region_data_lazy(region_dict, region_post_dict)
	end

	num_regions = region_dict.count
	# Assuming 2D square for now:
	M = 10
	region_data_array = zeros(num_regions, M)
	# Format: [minX, maxX, minY, maxY, μlX, μUX, μlY, μUY, σUX, σUY]
	for i=1:num_regions - 1
		region_data_array[i,1:2] = region_dict[i]["x1"]
		region_data_array[i,3:4] = region_dict[i]["x2"]
		region_data_array[i,5:6] = region_post_dict[i][1]["x1"]
		region_data_array[i,7:8] = region_post_dict[i][1]["x2"]
		region_data_array[i,9:10] = region_post_dict[i][2][:]
	end

	region_data_array[num_regions,1:2] = region_dict[num_regions]["x1"]
	region_data_array[num_regions,3:4] = region_dict[num_regions]["x2"]

	return region_data_array
end

function serialize_region_data_lazy(region_dict::Dict, region_post_dict)
	num_regions = region_dict.count
	# Assuming 2D square for now:
	M = 5
	region_data_array = zeros(num_regions, M)
	# Format: [minX, maxX, minY, maxY, postmap]
	for i=1:num_regions - 1
		region_data_array[i,1:2] = region_dict[i]["x1"]
		region_data_array[i,3:4] = region_dict[i]["x2"]
		region_data_array[i,5] = 0.
	end

	region_data_array[num_regions,1:2] = region_dict[num_regions]["x1"]
	region_data_array[num_regions,3:4] = region_dict[num_regions]["x2"]

	return region_data_array
end

function deserialize_region_data(region_data_filename::String; lazy_flag=false)

	if lazy_flag
		return deserialize_region_data_lazy(region_data_filename)
	end

	region_data_array = BSON.load(region_data_filename)[:res]
	num_regions = size(region_data_array,1)
	pair_iterator = Base.product(1:num_regions, 1:num_regions)
	region_dict = Dict()
	region_post_dict = Dict()

	for i=1:num_regions-1
		region_dict[i] = Dict()
		region_dict[i]["x1"] = region_data_array[i,1:2]
		region_dict[i]["x2"] = region_data_array[i,3:4]

		region_post_dict[i] = (Dict(), region_data_array[i,9:10])
		region_post_dict[i][1]["x1"] = region_data_array[i,5:6]
		region_post_dict[i][1]["x2"] = region_data_array[i,7:8]
	end

	region_dict[num_regions] = Dict()
	region_dict[num_regions]["x1"] = region_data_array[num_regions,1:2]
	region_dict[num_regions]["x2"] = region_data_array[num_regions,3:4]

	return region_dict, region_post_dict, pair_iterator
end

function deserialize_region_data_lazy(region_data_filename::String)
	region_data_array = BSON.load(region_data_filename)[:res]
	num_regions = size(region_data_array,1)
	pair_iterator = Base.product(1:num_regions, 1:num_regions)
	region_dict = Dict()
	region_post_dict = Dict()

	for i=1:num_regions-1
		region_dict[i] = Dict()
		region_dict[i]["x1"] = region_data_array[i,1:2]
		region_dict[i]["x2"] = region_data_array[i,3:4]

		# TODO: Ignoring post map for now
	end

	region_dict[num_regions] = Dict()
	region_dict[num_regions]["x1"] = region_data_array[num_regions,1:2]
	region_dict[num_regions]["x2"] = region_data_array[num_regions,3:4]

	return region_dict, region_post_dict, pair_iterator
end


text:
#!usr/bin/env julia
#testing version: 0.6.0

__precompile__()

module Units

import Base: show, +, -, *, /, ^, sqrt, exp, log
export Physical, show, +, -, *, /, ^, sqrt, exp, log

export meter, kilogram, second, ampere, kelvin, mole, candela
export hertz, newton, pascal, joule, watt, coulomb, volt, farad, ohm, weber, tesla, henry, lumen, lux

struct Physical{T, m, kg, s, A, K, mol, cd}
    value :: T
end

Physical(x, m, kg, s, A, K, mol, cd) = Physical{typeof(x), m, kg, s, A, K, mol, cd}(x)
Physical(x; m=0, kg=0, s=0, A=0, K=0, mol=0, cd=0) = Physical{typeof(x), m, kg, s, A, K, mol, cd}(x)

@inline character_lift(x) = Dict(
    '0' => '⁰',
    '1' => '¹',
    '2' => '²',
    '3' => '³',
    '4' => '⁴',
    '5' => '⁵',
    '6' => '⁶',
    '7' => '⁷',
    '8' => '⁸',
    '9' => '⁹',
)[x]

function show(io::IO, x::Physical{T, m, kg, s, A, K, mol, cd}) where {T, m, kg, s, A, K, mol, cd}
    print(io, "$(x.value)")
    for (u, v) in [(kg, :kg), (A, :A), (m, :m), (K, :K), (mol, :mol), (cd, :cd), (s, :s)]
        if u != 0
            print(io, '⋅', v)
            if u < 0
                print(io, '⁻')
            end
            if u != 1
                print(io, map(character_lift, collect("$(abs(u))")) |> String)
            end
        end
    end
end

include("arithmetic.jl")
include("constants.jl")

end #module


text:
module test_julia_show

using Test

buf = IOBuffer()
ioc = IOContext(buf, :SHOWN_SET => 5)
Base.show_circular(ioc, 5)
@test String(take!(buf)) == "#= circular reference @-1 =#"

end # module test_julia_show


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# This file contains objective-related functions

"""
    relative_gap(model::Model)

Return the final relative optimality gap after a call to `optimize!(model)`.
Exact value depends upon implementation of MathOptInterface.RelativeGap()
by the particular solver used for optimization.
"""
function relative_gap(model::Model)::Float64
    return MOI.get(model, MOI.RelativeGap())
end

"""
    objective_bound(model::Model)

Return the best known bound on the optimal objective value after a call to
`optimize!(model)`.
"""
function objective_bound(model::Model)::Float64
    return MOI.get(model, MOI.ObjectiveBound())
end

"""
    objective_value(model::Model; result::Int = 1)

Return the objective value associated with result index `result` of the
most-recent solution returned by the solver.

See also: [`result_count`](@ref).
"""
function objective_value(model::Model; result::Int = 1)::Float64
    return MOI.get(model, MOI.ObjectiveValue(result))
end

"""
    dual_objective_value(model::Model; result::Int = 1)

Return the value of the objective of the dual problem associated with result
index `result` of the most-recent solution returned by the solver.

Throws `MOI.UnsupportedAttribute{MOI.DualObjectiveValue}` if the solver does
not support this attribute.

See also: [`result_count`](@ref).
"""
function dual_objective_value(model::Model; result::Int = 1)::Float64
    return MOI.get(model, MOI.DualObjectiveValue(result))
end

"""
    objective_sense(model::Model)::MathOptInterface.OptimizationSense

Return the objective sense.
"""
function objective_sense(model::Model)
    return MOI.get(model, MOI.ObjectiveSense())::MOI.OptimizationSense
end

"""
    set_objective_sense(model::Model, sense::MathOptInterface.OptimizationSense)

Sets the objective sense of the model to the given sense. See
[`set_objective_function`](@ref) to set the objective function. These are
low-level functions; the recommended way to set the objective is with the
[`@objective`](@ref) macro.
"""
function set_objective_sense(model::Model, sense::MOI.OptimizationSense)
    return MOI.set(model, MOI.ObjectiveSense(), sense)
end

"""
    set_objective_function(
        model::Model,
        func::Union{AbstractJuMPScalar, MathOptInterface.AbstractScalarFunction})

Sets the objective function of the model to the given function. See
[`set_objective_sense`](@ref) to set the objective sense. These are low-level
functions; the recommended way to set the objective is with the
[`@objective`](@ref) macro.
"""
function set_objective_function end

function set_objective_function(model::Model, func::MOI.AbstractScalarFunction)
    attr = MOI.ObjectiveFunction{typeof(func)}()
    if !MOI.supports(backend(model), attr)
        error(
            "The solver does not support an objective function of type ",
            typeof(func),
            ".",
        )
    end
    MOI.set(model, attr, func)
    # Nonlinear objectives override regular objectives, so if there was a
    # nonlinear objective set, we must clear it.
    if model.nlp_data !== nothing
        model.nlp_data.nlobj = nothing
    end
    return
end

function set_objective_function(model::Model, func::AbstractJuMPScalar)
    check_belongs_to_model(func, model)
    return set_objective_function(model, moi_function(func))
end

function set_objective_function(model::Model, func::Real)
    return set_objective_function(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm{Float64}[],
            Float64(func),
        ),
    )
end

function set_objective_function(model::AbstractModel, ::MutableArithmetics.Zero)
    return set_objective_function(model, 0.0)
end

function set_objective_function(model::AbstractModel, func)
    return error("The objective function `$(func)` is not supported by JuMP.")
end

"""
    set_objective(model::AbstractModel, sense::MOI.OptimizationSense, func)

The functional equivalent of the [`@objective`](@ref) macro.

Sets the objective sense and objective function simultaneously, and is
equivalent to:
```julia
set_objective_sense(model, sense)
set_objective_function(model, func)
```

## Examples

```jldoctest; setup=:(using JuMP)
model = Model()
@variable(model, x)
set_objective(model, MOI.MIN_SENSE, x)
```
"""
function set_objective(model::AbstractModel, sense::MOI.OptimizationSense, func)
    set_objective_sense(model, sense)
    return set_objective_function(model, func)
end

"""
    objective_function_type(model::Model)::AbstractJuMPScalar

Return the type of the objective function.
"""
function objective_function_type(model::Model)
    return jump_function_type(
        model,
        MOI.get(backend(model), MOI.ObjectiveFunctionType()),
    )
end

"""
    objective_function(model::Model,
                   T::Type{<:AbstractJuMPScalar}=objective_function_type(model))

Return an object of type `T` representing the objective function.
Error if the objective is not convertible to type `T`.

## Examples

```jldoctest objective_function; setup = :(using JuMP)
julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, x)
x

julia> @objective(model, Min, 2x + 1)
2 x + 1

julia> objective_function(model, AffExpr)
2 x + 1

julia> objective_function(model, QuadExpr)
2 x + 1

julia> typeof(objective_function(model, QuadExpr))
GenericQuadExpr{Float64,VariableRef}
```
We see with the last two commands that even if the objective function is affine,
as it is convertible to a quadratic function, it can be queried as a quadratic
function and the result is quadratic.

However, it is not convertible to a variable.
```jldoctest objective_function; filter = r"MathOptInterface\\."s
julia> objective_function(model, VariableRef)
ERROR: InexactError: convert(MathOptInterface.SingleVariable, MathOptInterface.ScalarAffineFunction{Float64}(MathOptInterface.ScalarAffineTerm{Float64}[MathOptInterface.ScalarAffineTerm{Float64}(2.0, MathOptInterface.VariableIndex(1))], 1.0))
[...]
```
"""
function objective_function(
    model::Model,
    FunType::Type{<:AbstractJuMPScalar} = objective_function_type(model),
)
    MOIFunType = moi_function_type(FunType)
    func =
        MOI.get(backend(model), MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
    return jump_function(model, func)
end

"""
    set_objective_coefficient(model::Model, variable::VariableRef, coefficient::Real)

Set the linear objective coefficient associated with `Variable` to `coefficient`.

Note: this function will throw an error if a nonlinear objective is set.
"""
function set_objective_coefficient(
    model::Model,
    variable::VariableRef,
    coeff::Real,
)
    if model.nlp_data !== nothing && _nlp_objective_function(model) !== nothing
        error("A nonlinear objective is already set in the model")
    end
    coeff = convert(Float64, coeff)::Float64
    obj_fct_type = objective_function_type(model)
    if obj_fct_type == VariableRef
        # Promote the objective function to be an affine expression.
        current_obj = objective_function(model)
        if index(current_obj) == index(variable)
            set_objective_function(model, coeff * variable)
        else
            set_objective_function(
                model,
                add_to_expression!(coeff * variable, current_obj),
            )
        end
    elseif obj_fct_type == AffExpr || obj_fct_type == QuadExpr
        MOI.modify(
            backend(model),
            MOI.ObjectiveFunction{moi_function_type(obj_fct_type)}(),
            MOI.ScalarCoefficientChange(index(variable), coeff),
        )
    else
        error("Objective function type not supported: $(obj_fct_type)")
    end
    return
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Random Pairs" begin
  @test isa(rand(Pair{String,Int}), Pair{String,Int})
end


text:
# yuck. there were less fucks given making this
"""

    FVMesh(
        model::SFFM.Model,
        Nodes::Array{Float64,1};
        Fil::IndexDict=IndexDict(),
    )

Constructor for a mesh for a finite volume scheme.
    Inputs:
     - `model::Model` a Model object
     - `Nodes::Array{Float64,1}` a vector specifying the cell edges
     - `Fil::IndexDict` an optional dictionary allocating the cells to the sets Fᵢᵐ
"""
struct FVMesh <: SFFM.Mesh
    Nodes::Array{Float64,1}
    order::Int
    Fil::IndexDict
end
function FVMesh(
    model::SFFM.Model,
    Nodes::Array{Float64,1},
    order::Int;
    Fil::IndexDict=IndexDict(),
)
    ## Construct the sets Fᵐ = ⋃ᵢ Fᵢᵐ, global index for sets of type m
    if isempty(Fil)
        Fil = MakeFil(model, Nodes)
    end

    return FVMesh(Nodes, order, Fil)
end

"""

    NBases(mesh::FVMesh)

Constant 1
"""
NBases(mesh::FVMesh) = 1
Order(mesh::FVMesh) = mesh.order

"""

    CellNodes(mesh::FVMesh)

The cell centres
"""
CellNodes(mesh::FVMesh) = Array(((mesh.Nodes[1:end-1] + mesh.Nodes[2:end]) / 2 )')

"""

    Basis(mesh::FVMesh)

Constant ""
"""
Basis(mesh::FVMesh) = ""

TotalNBases(mesh::FVMesh) = NIntervals(mesh)

function MakeFVFlux(mesh::SFFM.Mesh, order::Int)
    nNodes = TotalNBases(mesh)
    F = zeros(Float64,nNodes,nNodes)
    ptsLHS = Int(ceil(order/2))
    interiorCoeffs = lagrange_poly_basis(CellNodes(mesh)[1:order],mesh.Nodes[ptsLHS+1])
    for n in 2:nNodes
        evalPt = mesh.Nodes[n]
        if n-ptsLHS-1 < 0
            nodesIdx = 1:order
            nodes = CellNodes(mesh)[nodesIdx]
            coeffs = lagrange_poly_basis(nodes,evalPt)
        elseif n-ptsLHS-1+order > nNodes
            nodesIdx = (nNodes-order+1):nNodes
            nodes = CellNodes(mesh)[nodesIdx]
            coeffs = lagrange_poly_basis(nodes,evalPt)
        else
            nodesIdx =  (n-ptsLHS-1) .+ (1:order)
            coeffs = interiorCoeffs
        end
        F[nodesIdx,n-1:n] += [-coeffs coeffs]./Δ(mesh)[n-1]
    end
    F[end-order+1:end,end] += -lagrange_poly_basis(CellNodes(mesh)[end-order+1:end],mesh.Nodes[end])./Δ(mesh)[end]
    return F
end

function MakeFullGenerator(model::SFFM.Model, mesh::SFFM.FVMesh; v::Bool=false)
    N₊ = sum(model.C .>= 0)
    N₋ = sum(model.C .<= 0)
    order = Order(mesh)
    F = SFFM.MakeFVFlux(mesh, order)

    B = SparseArrays.spzeros(
        Float64,
        NPhases(model) * NIntervals(mesh) + N₋ + N₊,
        NPhases(model) * NIntervals(mesh) + N₋ + N₊,
    )
    B[N₋+1:end-N₊,N₋+1:end-N₊] = SparseArrays.kron(
            model.T,
            SparseArrays.I(NIntervals(mesh))
        )

    ## Make QBD index
    QBDidx = MakeQBDidx(model,mesh)

    # Boundary conditions
    T₋₋ = model.T[model.C.<=0,model.C.<=0]
    T₊₋ = model.T[model.C.>=0,:].*((model.C.<0)')
    T₋₊ = model.T[model.C.<=0,:].*((model.C.>0)')
    T₊₊ = model.T[model.C.>=0,model.C.>=0]
    # yuck
    begin
        nodes = CellNodes(mesh)[1:order]
        coeffs = lagrange_poly_basis(nodes,mesh.Nodes[1])
        idxdown = ((1:order).+TotalNBases(mesh)*(findall(model.C .<= 0) .- 1)')[:] .+ N₋
        B[idxdown, 1:N₋] = LinearAlgebra.kron(
            LinearAlgebra.diagm(0 => model.C[model.C.<=0]),
            -coeffs,
        )
    end
    # inLower = [
    #     SparseArrays.diagm(abs.(model.C).*(model.C.<=0))[:,model.C.<=0];
    #     SparseArrays.zeros((NIntervals(mesh)-1)*NPhases(model),N₋)
    # ]
    outLower = [
        T₋₊./Δ(mesh)[1] SparseArrays.zeros(N₋,N₊+(NIntervals(mesh)-1)*NPhases(model))
    ]
    begin
        nodes = CellNodes(mesh)[end-order+1:end]
        coeffs = lagrange_poly_basis(nodes,mesh.Nodes[end])
        idxup =
            ((1:order).+TotalNBases(mesh)*(findall(model.C .>= 0) .- 1)')[:] .+
            (N₋ + TotalNBases(mesh) - order)
        B[idxup, (end-N₊+1):end] = LinearAlgebra.kron(
            LinearAlgebra.diagm(0 => model.C[model.C.>=0]),
            coeffs,
        )
    end
    # inUpper = [
    #     SparseArrays.zeros((NIntervals(mesh)-1)*NPhases(model),N₊);
    #     (SparseArrays.diagm(abs.(model.C).*(model.C.>=0)))[:,model.C.>=0]
    # ]
    outUpper = [
        SparseArrays.zeros(N₊,N₋+(NIntervals(mesh)-1)*NPhases(model)) T₊₋./Δ(mesh)[end]
    ]

    B[1:N₋,QBDidx] = [T₋₋ outLower]
    B[end-N₊+1:end,QBDidx] = [outUpper T₊₊]
    # B[QBDidx[N₋+1:end-N₊],1:N₋] = inLower
    # B[QBDidx[N₋+1:end-N₊],(end-N₊+1):end] = inUpper
    for i = 1:NPhases(model)
        idx = ((i-1)*NIntervals(mesh)+1:i*NIntervals(mesh)) .+ N₋
        if model.C[i] > 0
            B[idx, idx] += model.C[i] * F
        elseif model.C[i] < 0
            B[idx, idx] += abs(model.C[i]) * F[end:-1:1,end:-1:1]
        end
    end

    BDict = SFFM.MakeDict(B,model,mesh)
    out = FullGenerator(BDict, B, mesh.Fil)
    v && println("FullGenerator created with keys ", keys(out))
    return out
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
export read_trace

"""
Use:    read_trace!(s::IO,
                    fh::BinaryFileHeader,
                    datatype::Type,
                    headers::AbstractArray{BinaryTraceHeader,1},
                    data::AbstractArray{<:Union{IBMFloat32, Float32}, 2},
                    trace::Int,
                    th_byte2sample::Dict{String,Int32})

Reads 'trace' from the current position of stream 's' into 'headers' and
'data'.
"""
function read_trace!(s::IO, fh::BinaryFileHeader,
                    datatype::Type,
                    headers::AbstractArray{BinaryTraceHeader,1},
                    data::AbstractArray{<:Union{IBMFloat32, Float32}, 2},
                    trace::Int,
                    th_byte2sample::Dict{String,Int32})

    # Read trace header
    setindex!(headers, read_traceheader(s, th_byte2sample), trace)

    # Read trace
    setindex!(data, read_tracedata(s, fh, datatype), :, trace)

    nothing
end

"""
Use:    read_trace!(s::IO,
                    fh::BinaryFileHeader,
                    datatype::Type,
                    headers::AbstractArray{BinaryTraceHeader,1},
                    data::AbstractArray{<:Union{IBMFloat32, Float32}, 2},
                    trace::Int,
                    keys::Array{String,1},
                    th_byte2sample::Dict{String,Int32})

Reads 'trace' from the current position of stream 's' into 'headers' and
'data'. Only the header values in 'keys' and read.
"""
function read_trace!(s::IO, fh::BinaryFileHeader,
                    datatype::Type,
                    headers::AbstractArray{BinaryTraceHeader,1},
                    data::AbstractArray{<:Union{IBMFloat32, Float32}, 2},
                    trace::Int,
                    keys::Array{String,1},
                    th_byte2sample::Dict{String,Int32})

    # Read trace header
    setindex!(headers, read_traceheader(s, keys, th_byte2sample), trace)

    # Read trace
    setindex!(data, read_tracedata(s, fh, datatype), :, trace)

    nothing
end

"""
Read a single trace from current position in stream as Float32
"""
function read_tracedata(s::IO, fh::BinaryFileHeader, dsf::Type{DT}) where {DT<:IBMFloat32}
    dummy=Vector{DT}(undef, fh.ns)
    read!(s, dummy)
    return dummy
end

function read_tracedata(s::IO, fh::BinaryFileHeader, dsf::Type{DT}) where {DT<:Float32}
    dummy=Vector{DT}(undef, fh.ns)
    read!(s, dummy)
    return bswap.(dummy)
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    isomorphismitervf2, edgeisomorphismitervf2,
    graphmatches, graphmatch, isgraphmatch,
    subgraphmatches, subgraphmatch, issubgraphmatch,
    edgesubgraphmatches, edgesubgraphmatch, isedgesubgraphmatch

mutable struct VF2State{T1<:UndirectedGraph,T2<:UndirectedGraph}
    G::T1
    H::T2

    g_core::Dict{Int,Int}
    h_core::Dict{Int,Int}
    g_term::Dict{Int,Int}
    h_term::Dict{Int,Int}
    expire::Union{UInt64,Nothing} # UInt64, nanoseconds

    mappings::Vector{Dict{Int,Int}}
    status::Symbol

    function VF2State{T1,T2}(G::T1, H::T2; timeout=nothing
            ) where {T1<:UndirectedGraph,T2<:UndirectedGraph}
        if timeout !== nothing
            expire = (time_ns() + timeout * 1_000_000_000)::UInt64
        else
            expire = nothing
        end
        return new(G, H, Dict(), Dict(), Dict(), Dict(), expire, [], :ready)
    end
end

function candidatepairs(state::VF2State; kwargs...)
    # Mandatory pair
    if haskey(kwargs, :mandatory)
        md = setdiff(keys(kwargs[:mandatory]), keys(state.g_core))
        if !isempty(md)
            n = pop!(md)
            return [(n, kwargs[:mandatory][n])]
        end
    end

    pairs = Tuple{Int,Int}[]
    g_cand = setdiff(keys(state.g_term), keys(state.g_core))
    h_cand = setdiff(keys(state.h_term), keys(state.h_core))
    if isempty(g_cand) || isempty(h_cand)
        # New connected component
        g_cand = setdiff(nodeset(state.G), keys(state.g_core))
        h_cand = setdiff(nodeset(state.H), keys(state.h_core))
    end
    if !isempty(h_cand)
        h_min = minimum(h_cand)
        for g in g_cand
            # Forbidden pair
            if haskey(kwargs, :forbidden)
                fb = kwargs[:forbidden]
                haskey(fb, g) && fb[g] == h_min && continue
            end
            push!(pairs, (g, h_min))
        end
    end
    return pairs
end

function is_feasible(state::VF2State, g, h; mode=:Subgraph, kwargs...)
    # TODO: assume no self loop
    g_nbrs = adjacencies(state.G, g)
    h_nbrs = adjacencies(state.H, h)
    for n in intersect(g_nbrs, keys(state.g_core))
        if !(state.g_core[n] in h_nbrs)
            # println("Infeasible: neighbor connectivity")
            return false
        end
    end
    for n in intersect(h_nbrs, keys(state.h_core))
        if !(state.h_core[n] in g_nbrs)
            # println("Infeasible: neighbor connectivity")
            return false
        end
    end
    g_term_count = length(setdiff(keys(state.g_term), keys(state.g_core)))
    h_term_count = length(setdiff(keys(state.h_term), keys(state.h_core)))
    if mode == :Isomorphism && g_term_count != h_term_count
        # println("Infeasible: terminal set size")
        return false
    elseif mode == :Subgraph && g_term_count < h_term_count
        # println("Infeasible: terminal set size")
        return false
    end
    g_new_count = length(setdiff(nodeset(state.G), keys(state.g_term)))
    h_new_count = length(setdiff(nodeset(state.H), keys(state.h_term)))
    if mode == :Isomorphism && g_new_count != h_new_count
        # println("Infeasible: yet unexplored nodes")
        return false
    elseif mode == :Subgraph && g_new_count < h_new_count
        # println("Infeasible: yet unexplored nodes")
        return false
    end
    # println("Syntactically feasible")
    return true
end

function is_semantic_feasible(state::VF2State, g, h; kwargs...)
    if haskey(kwargs, :nodematcher)
        if !kwargs[:nodematcher](g, h)
            # println("Infeasible: node attribute mismatch")
            return false
        end
    end
    if haskey(kwargs, :edgematcher)
        for (inc, adj) in neighbors(state.G, g)
            haskey(state.g_core, adj) || continue
            hinc = findedgekey(state.H, h, state.g_core[adj])
            if !kwargs[:edgematcher](inc, hinc)
                # println("Infeasible: edge attribute mismatch")
                return false
            end
        end
    end
    return true
end

function updatestate!(state::VF2State, g, h)
    state.g_core[g] = h
    state.h_core[h] = g
    depth = length(state.g_core)
    if !haskey(state.g_term, g)
        state.g_term[g] = depth
    end
    if !haskey(state.h_term, h)
        state.h_term[h] = depth
    end
    g_nbrset = union([adjacencies(state.G, n) for n in keys(state.g_core)]...)
    for n in setdiff(g_nbrset, keys(state.g_term))
        state.g_term[n] = depth
    end
    h_nbrset = union([adjacencies(state.H, n) for n in keys(state.h_core)]...)
    for n in setdiff(h_nbrset, keys(state.h_term))
        state.h_term[n] = depth
    end
    return
end

function restore!(state::VF2State, g, h)
    depth = length(state.g_core)
    if g !== nothing && h !== nothing
        delete!(state.g_core, g)
        delete!(state.h_core, h)
    end
    for (k, v) in state.g_term
        v == depth && delete!(state.g_term, k)
    end
    for (k, v) in state.h_term
        v == depth && delete!(state.h_term, k)
    end
end

function yieldallnodes!(state::VF2State)
    push!(state.mappings, copy(state.g_core))
    return
end

function yieldfirstnode!(state::VF2State)
    push!(state.mappings, copy(state.g_core))
    state.status = :done
    return
end

function yieldedgefunc(G, H, yieldtype)
    return function (state)
        delta_y_mismatch(G, H, state.g_core) && return
        push!(state.mappings, copy(state.g_core))
        yieldtype == :first && (state.status = :done)
        return
    end
end

function expand!(state::VF2State; kwargs...)
    state.status == :done && return
    # Recursive
    # println("depth $(length(state.g_core))")
    if length(state.g_core) == nodecount(state.H)
        # println("done $(state.g_core)")
        kwargs[:reportfunc](state)
        return
    elseif state.expire !== nothing && time_ns() > state.expire
        state.status = :timedout
        # println("Timed out")
        return
    end
    for (g, h) in candidatepairs(state; kwargs...)
        # println("candidates $(g) $(h)")
        is_feasible(state, g, h; kwargs...) || continue
        is_semantic_feasible(state, g, h; kwargs...) || continue
        updatestate!(state, g, h)
        # println("g_core $(state.g_core)")
        mp = expand!(state; kwargs...)
        # println("restored $(state.g_core)")
        restore!(state, g, h)
    end
    return
end

function isomorphismitervf2(G::T1, H::T2; timeout=nothing, yield=:all, kwargs...
        ) where {T1<:UndirectedGraph,T2<:UndirectedGraph}
    (nodecount(G) == 0 || nodecount(H) == 0) && return Dict{Int,Int}()
    state = VF2State{T1,T2}(G, H, timeout=timeout)
    yieldfunc = yield == :all ? yieldallnodes! : yieldfirstnode!
    expand!(state, reportfunc=yieldfunc; kwargs...)
    if state.status == :timedout
        throw(ErrorException("Timeout reached"))
    end

    return state.mappings
end

function edgeisomorphismitervf2(
            G::T1, H::T2; timeout=nothing, yield=:all,
            nodematcher=(g,h)->true , edgematcher=(g,h)->true,
            kwargs...
        ) where {T1<:UndirectedGraph,T2<:UndirectedGraph}
    (edgecount(G) == 0 || edgecount(H) == 0) && return Dict{Int,Int}()
    lg = linegraph(G)
    lh = linegraph(H)
    state = VF2State{LineGraph,LineGraph}(lg, lh, timeout=timeout)
    lgematch = lgedgematcher(lg, lh, nodematcher)
    lgnmatch = lgnodematcher(lg, lh, nodematcher, edgematcher)
    yieldfunc = yieldedgefunc(G, H, yield)
    expand!(state, nodematcher=lgnmatch, edgematcher=lgematch,
            reportfunc=yieldfunc; kwargs...)
    if state.status == :timedout
        throw(ErrorException("Timeout reached"))
    end
    return state.mappings
end

# Graph isomorphism

"""
    graphmatches(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Iterator

Generate isomorphism mappings between `G` and `H`. If no match found, return
nothing.
"""
graphmatches(G, H; kwargs...
    ) = isomorphismitervf2(G, H, mode=:Isomorphism; kwargs...)

"""
    graphmatch(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Dict{Int,Int}

Return an isomorphism mapping between `G` and `H`. If no match found, return
nothing.
"""
function graphmatch(G, H; kwargs...)
    res = isomorphismitervf2(G, H, mode=:Isomorphism, yield=:first; kwargs...)
    return isempty(res) ? nothing : res[1]
end

"""
    isgraphmatch(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Bool

Return true if `G` and `H` are isomorphic.
"""
isgraphmatch(G, H; kwargs...) = graphmatch(G, H; kwargs...) !== nothing

# Node induced subgraph isomorphism

"""
    subgraphmatches(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Iterator

Generate subgraph isomorphism mappings between `G` and `H`.

# Keyword arguments

- nodematcher(Function): node matcher function that takes two node indices as
arguments.
- edgematcher(Function): edge matcher function that takes two edge indices as
arguments.
- mandatory(Dict{Int,Int}): mandatory node matches (available for only VF2)
- forbidden(Dict{Int,Int}):
    forbidden node matches (available for only VF2)
"""
subgraphmatches(G, H; kwargs...
    ) = isomorphismitervf2(G, H, mode=:Subgraph; kwargs...)

"""
    subgraphmatch(
        G::AbstractGraph, H::AbstractGraph; kwargs...) -> Dict{Int,Int}

Return a subgraph isomorphism mapping between `G` and `H`. If no match found,
return nothing.
"""
function subgraphmatch(G, H; kwargs...)
    res = isomorphismitervf2(G, H, mode=:Subgraph, yield=:first; kwargs...)
    return isempty(res) ? nothing : res[1]
end

"""
    issubgraphmatch(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Bool

Return true if a node induced subgraph of `G` and `H` are isomorphic.
"""
issubgraphmatch(G, H; kwargs...) = subgraphmatch(G, H; kwargs...) !== nothing

# Edge induced subgraph isomorphism

"""
    edgesubgraphmatch(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Iterator

Generate edge induced subgraph isomorphism mappings between `G` and `H`.
"""
edgesubgraphmatches(G, H; kwargs...
    ) = edgeisomorphismitervf2(G, H, mode=:Subgraph; kwargs...)

"""
    edgesubgraphmatch(
        G::AbstractGraph, H::AbstractGraph; kwargs...) -> Dict{Int,Int}

Return a edge induced subgraph isomorphism mapping between `G` and `H`.
If no match found, return nothing.
"""
function edgesubgraphmatch(G, H; kwargs...)
    res = edgeisomorphismitervf2(G, H, mode=:Subgraph, yield=:first; kwargs...)
    return isempty(res) ? nothing : res[1]
end

"""
    isedgesubgraphmatch(G::AbstractGraph, H::AbstractGraph; kwargs...) -> Bool

Return true if a node induced subgraph of `G` and `H` are isomorphic.
"""
isedgesubgraphmatch(G, H; kwargs...
    ) = edgesubgraphmatch(G, H; kwargs...) !== nothing


text:
using FCCQuad: Fct, Cheb, getweights, fccquad_cc, fccquad_cs, fccquad_sc, fccquad_ss, all_methods, supported_types, interval_methods
using QuadGK, Test

#=
Integral of cos(13cos(2x))*52sin(x)cos(x)dx from 0 to 17pi/4
  == integral of cos(u)*du from 0 to 13 == sin(13).
Integral of sin(13cos(2x))*52sin(x)cos(x)dx from 0 to 17pi/4
  == integral of sin(u)*du from 0 to 13 == 1-cos(13).
=#
@testset "analytic integrals" begin
    println("Integral of (cos or sin)(13cos(2x))*52sin(x)cos(x)dx from 0 to 17pi/4:")
    angle(x)=13cos(2x)
    ampcos(x)=52cos(x)
    ampsin(x)=52sin(x)
    println("\tmethod \tdata type \tintegrand evaluation counts")
    for mthd in all_methods
        for T in supported_types[mthd]
            print("\t$mthd \t$T \t")
            R=real(T)
            freqs=[one(R)]
            kwargs=(:xmin=>zero(R),:xmax=>17*R(pi)/4,:method=>mthd,:T=>R)

            exact=sin(R(13))
            result=fccquad_cc(ampsin,angle,freqs;kwargs...)
            print("$(result[2]) \t")
            @test isapprox(result[1][1,1],exact)

            result=fccquad_cs(ampcos,angle,freqs;kwargs...)
            print("$(result[2]) \t")
            @test isapprox(result[1][1,1],exact)

            exact=one(R)-cos(R(13))
            result=fccquad_sc(ampsin,angle,freqs;kwargs...)
            print("$(result[2]) \t")
            @test isapprox(result[1][1,1],exact)

            result=fccquad_ss(ampcos,angle,freqs;kwargs...)
            println("$(result[2])")
            @test isapprox(result[1][1,1],exact)
        end
    end
end

#slow type I DCT (for testing purposes):
#(dct(x_0,...,x_N))_n = sum_{0<=m<=N}x_m*cos(pi*m*n/N)
#with 1/2 weight given to first and last summands.
#N can be any positive integer
function dct(x::AbstractVector,T::Type=Complex{Float64})::Vector{T} #type I DCT
    n=length(x)-1
    y=Vector{T}(undef,n+1)
    p = real(T)(pi)
    for i in 0:n
        y[i+1]=0.5x[1]
        @simd for j in 1:n-1
            @inbounds y[i+1]+=x[j+1]*cos(p*i*j/n)
        end
        y[i+1]+=x[n+1]*ifelse(i%2==0,0.5,-0.5)
    end
    y
end

@testset "Comparing fast (I)DCT to simpler implementation" begin
    x0=rand(Complex{Float64},129)
    y1=Fct.fct(x0)
    y2=dct(x0)
    @test isapprox(y1,y2)
    x1=Fct.idct(y1)
    @test isapprox(x0,x1)
end

@testset "Chebyshev series evaluation, dilation, and multiplication" begin
    f=Dict(:cube=>x->x^3, :sin=>sin, :cos=>cos, :exp=>exp)
    c=Dict(:cube=>Cheb.OddChebSeries(x->x^3),
           :sin=>Cheb.OddChebSeries(sin),
           :cos=>Cheb.EvenChebSeries(cos),
           :exp=>Cheb.ChebSeries(exp))
    x0=-0.5
    kappa=0.9
    for k in keys(f)
        @test isapprox(Cheb.chebinterp(c[k],x0),f[k](x0))
        @test isapprox(Cheb.chebinterp(Cheb.dilate(c[k],kappa),x0),f[k](kappa*x0))
        for j in keys(f)
            @test isapprox(Cheb.chebinterp(c[j]*c[k],x0),f[j](x0)*f[k](x0))
        end
    end
end

#T_n(x)
function chebpoly(n::Integer,x)
    @assert n>=0
    n<=1 && return iszero(n) ? one(x) : x
    m = div(n,2)
    iseven(n) && return 2chebpoly(m,x)^2 - one(x) #T_{2m}(x)
    a,b = chebpair(m,x) #T_m(x), T_{m+1}(x)
    2a*b - x #T_{2m+1}(x)
end

#T_n(x),T_{n+1}(x)
function chebpair(n::Integer,x)
    @assert n>=0
    n==0 && return (one(x),x)
    m = div(n,2)
    e,o = chebpair(m,x) #T_m, T_{m+1}
    if iseven(n) #n == 2m
        2e^2 - one(x), 2e*o - x #T_{2m}, T_{2m+1}
    else # n == 2m+1
        2o*e - x, 2o^2 - one(x) #T_{2m+1}, T_{2m+2}
    end
end

@testset "Chebyshev weights" begin
    function checkweights(N,w,args...)
        a = getweights(N,w,args...)
        atol2 = N*eps()^2
        gk_atol=eps()
        gk_rtol=10^4 * eps()
        for n in 0:2:N
        gk,delta = collect(quadgk(x->chebpoly(n,x)*cos(w*x),-1,1,rtol=gk_rtol,atol=gk_atol))
            pass = abs2(a[n+1]-gk) <= max(atol2,abs2(delta))
            pass || println("N=$N,n=$n,omega=$w,weight=$(a[n+1]),gk=$gk,gk_delta=$delta")
            @test pass
        end
        for n in 1:2:N
            gk,delta = collect(quadgk(x->chebpoly(n,x)*sin(w*x),-1,1,rtol=gk_rtol,atol=gk_atol))
            pass = abs2(a[n+1]-gk) <= max(atol2,abs2(delta))
            pass || println("N=$N,n=$n,omega=$w,weight=$(a[n+1]),gk=$gk,gk_delta=$delta")
            @test pass
        end
    end
    println("Testing Chebyshev weights (0<=n<=100):")
    for freq in [0,-1e-9,1e-6,-1e-3,1e-2,-1e-1,1,2,3,4,5,-99.5,200]
        println("\ttesting angular frequency $freq...")
        checkweights(100,freq)
    end
end


text:
include("rivers.jl")
include("horizon.jl")
include("resolution.jl")
include("segmenter.jl")
include("hydrodata.jl")
include("water_value.jl")
include("pricedata.jl")
include("flowdata.jl")
include("forecaster.jl")

Scenario = Int

abstract type AbstractModelIndices end

function plants(indices::AbstractModelIndices)
    :plants ∈ fieldnames(typeof(indices)) || error("Model indices do not contain hydro plants")
    return indices.plants
end

abstract type AbstractModelData end

function hydrodata(data::AbstractModelData)
    :hydrodata ∈ fieldnames(typeof(data)) || error("Model data does not contain hydro data")
    return data.hydrodata
end

modelindices(data::AbstractModelData; kw...) = error("No definition of modelindices for ", typeof(data), " with arguments: ", args...)