text:
using StatsBase

typealias  Graph Dict{Int, Set{Int}}

function addnode!(g::Graph, u::Int)
    if !haskey(g, u)
        g[u] = valtype(g)()
    end
    nothing
end

function delnode!(g::Graph, u::Int)
    if haskey(g, u)
        for v in g[u]
            delete!(g[v], u)
        end
        delete!(g, u)
    end
    nothing
end

function addedge!(g::Graph, u::Int, v::Int)
    addnode!(g, u)
    addnode!(g, v)
    push!(g[u], v)
    push!(g[v], u)
    nothing
end

function deledge!(g::Graph, u::Int, v::Int)
    if haskey(g, u) && haskey(g, v)
        delete!(g[u], v)
        delete!(g[v], u)
    end
    nothing
end

function ne(g::Graph)
    m = 0
    for u in keys(g)
        m += length(g[u])
    end
    return div(m, 2)
end

neighbors(g::Graph, u::Int) = g[u]

nodes(g::Graph, u::Int) = keys(g)

"""Type to record neighbor labels and their counts."""
type NeighComm
    neigh_pos::Vector{Int}
    neigh_cnt::Vector{Int}
    neigh_last::Int
end

function range_shuffle!(r, a::AbstractVector)
    @inbounds for i=r:-1:2
        j = StatsBase.randi(i)
        a[i],a[j] = a[j],a[i]
    end
end

"""Return the most frequency label."""
function pre_vote!(g, m, c::NeighComm, u, all_active_nodes)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in neighbors(g,u)
        if in(neigh, all_active_nodes)
            neigh_comm = m[neigh]
            if c.neigh_cnt[neigh_comm] < 0
                c.neigh_cnt[neigh_comm] = 0
                c.neigh_pos[c.neigh_last] = neigh_comm
                c.neigh_last += 1
            end
            c.neigh_cnt[neigh_comm] += 1
            if c.neigh_cnt[neigh_comm] > max_cnt
                max_cnt = c.neigh_cnt[neigh_comm]
            end
        end
    end
    # ties breaking randomly
    range_shuffle!(c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
        if c.neigh_cnt[lbl] == max_cnt
            return lbl
        end
    end
end

function pre_update!(g, m, c::NeighComm, active_nodes, all_active_nodes)
    while !isempty(active_nodes)
        random_order = collect(active_nodes)
        shuffle!(random_order)
        for u in random_order
            old_comm = m[u]
            new_comm = pre_vote!(g, m, c, u, all_active_nodes)
            if new_comm != old_comm
                for v in neighbors(g,u)
                    if in(v, all_active_nodes)
                        push!(active_nodes, v)
                        m[u] = new_comm
                    end
                end
            else
                delete!(active_nodes, u)
            end
        end
    end
end

"""Return the most frequency label."""
function vote!(g, m, c::NeighComm, u)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in neighbors(g,u)
        neigh_comm = m[neigh]
        if c.neigh_cnt[neigh_comm] < 0
            c.neigh_cnt[neigh_comm] = 0
            c.neigh_pos[c.neigh_last] = neigh_comm
            c.neigh_last += 1
        end
        c.neigh_cnt[neigh_comm] += 1
        if c.neigh_cnt[neigh_comm] > max_cnt
            max_cnt = c.neigh_cnt[neigh_comm]
        end
    end
    # ties breaking randomly
    range_shuffle!(c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
        if c.neigh_cnt[lbl] == max_cnt
            return lbl
        end
    end
end

function update!(g, m, c::NeighComm, active_nodes)
    num_active_nodes = 0
    while !isempty(active_nodes)
        if length(active_nodes) > num_active_nodes
            num_active_nodes = length(active_nodes)
        end
        random_order = collect(active_nodes)
        shuffle!(random_order)
        for u in random_order
            old_comm = m[u]
            new_comm = vote!(g, m, c, u)
            if new_comm != old_comm
                for v in neighbors(g,u)
                    push!(active_nodes, v)
                    m[u] = new_comm
                end
            else
                delete!(active_nodes, u)
            end
        end
    end
    num_active_nodes
end

function lpa_addnode!(g, m, u)
    if !haskey(g, u)
        addnode!(g, u)
        m[u] = u
    end
end

function lpa_deletenode!(g, m, c, u, active_lbls, pre_active_nodes, active_nodes)
    if haskey(g, u)
        empty!(pre_active_nodes)
        empty!(active_nodes)
        empty!(active_lbls)
        for v in neighbors(g,u)
            push!(active_lbls, m[v])
        end
        deletenode!(g, u)
        delete!(m, u)
        for v in nodes(g)
            if in(m[v], active_lbls)
                push!(pre_active_nodes, v)
                push!(active_nodes, v)
                m[v] = v
            end
        end
        pre_update!(g, m, c, pre_active_nodes, active_nodes)
        return update!(g, m, c, active_nodes)
    end
    return 0
end

function lpa_addedge!(g, m, c, u, v, pre_active_nodes, active_nodes)
    lpa_addnode!(g, m, u)
    lpa_addnode!(g, m, v)
    if !in(v, g[u])
        u_comm = m[u]
        v_comm = m[v]
        if u_comm != v_comm
            ku_int = 0
            ku_ext = 0
            kv_int = 0
            kv_ext = 0
            for neigh in neighbors(g,u)
                if m[neigh] == u_comm
                    ku_int += 1
                else
                    ku_ext += 1
                end
            end
            for neigh in neighbors(g,v)
                if m[neigh] == v_comm
                    kv_int += 1
                else
                    kv_ext += 1
                end
            end
            if ku_ext+1 >= ku_int || kv_ext+1 >= kv_int
                empty!(pre_active_nodes)
                empty!(active_nodes)
                for i in nodes(g)
                    i_comm = m[i]
                    if i_comm == u_comm || i_comm == v_comm
                        push!(pre_active_nodes, i)
                        push!(active_nodes, i)
                        m[i] = i
                    end
                end
                addedge!(g, u, v)
                pre_update!(g, m, c, pre_active_nodes, active_nodes)
                return update!(g, m, c, active_nodes)
            else
                addedge!(g, u, v)
            end
        else
            addedge!(g, u, v)
        end
    end
    return 0
end

function lpa_deleteedge!(g, m, c, u, v, pre_active_nodes, active_nodes)
    if haskey(g, u) && haskey(g[u], v)
        u_comm = m[u]
        v_comm = m[v]
        if u_comm == v_comm
            empty!(pre_active_nodes)
            empty!(active_nodes)
            for i in nodes(g)
                if m[i] == u_comm
                    push!(pre_active_nodes, i)
                    push!(active_nodes, i)
                    m[i] = i
                end
            end
            deleteedge!(g, u, v)
            pre_update!(g, m, c, pre_active_nodes, active_nodes)
            return update!(g, m, c, active_nodes)
        else
            deleteedge!(g, u, v)
        end
    end
    return 0
end

using ParserCombinator: Parsers.GML

function _gml_read_one_graph(gs)
    nodes = [x[:id] for x in gs[:node]]
    g = Graph{Int,Int}()
    sds = [(Int(x[:source]), Int(x[:target])) for x in gs[:edge]]
    for (s,d) in (sds)
        addedge!(g, s, d)
    end
    return g
end

function loadgml(gname::AbstractString)
    p = GML.parse_dict(readall(gname))
    for gs in p[:graph]
        return _gml_read_one_graph(gs)
    end
    error("Graph $gname not found")
end

function savegml(gname::AbstractString, g, c=Vector{Integer}())
    io = open(gname, "w")
    println(io, "graph")
    println(io, "[")
    println(io, "directed 0")
    i = 0
    nodemap = Dict{Int,Int}()
    for u in keys(g)
        nodemap[u] = i
        println(io,"\tnode")
        println(io,"\t[")
        println(io,"\t\tid $i")
        println(io,"\t\tlabel $u")
        length(c) > 0 && println(io,"\t\tvalue $(c[u])")
        println(io,"\t]")
        i += 1
    end
    for u in keys(g)
        for v in keys(g[u])
            if u < v
                println(io,"\tedge")
                println(io,"\t[")
                println(io,"\t\tsource $(nodemap[u])")
                println(io,"\t\ttarget $(nodemap[v])")
                println(io,"\t]")
            end
        end
    end
    println(io, "]")
    close(io)
    return 1
end


text:
#= Github Actions uses Intel CPUs, so it is faster to use MKL than OpenBLAS.
It is recommended to load MKL before ANY other package.=#
#= Wait until https://github.com/JuliaLinearAlgebra/MKL.jl/issues/112 is resolved
before using MKL on CI for Mac =#
#import MKL
using SafeTestsets: @safetestset

@time @safetestset "util" begin include("util.jl") end
@time @safetestset "linalg" begin include("linalg.jl") end
@time @safetestset "onehot" begin include("onehot.jl") end
@time @safetestset "rbm" begin include("rbm.jl") end
@time @safetestset "layers" begin include("layers.jl") end
@time @safetestset "cd" begin include("cd.jl") end
@time @safetestset "pseudolikelihood" begin include("pseudolikelihood.jl") end
@time @safetestset "minibatches" begin include("minibatches.jl") end
@time @safetestset "initialization" begin include("initialization.jl") end
@time @safetestset "regularize" begin include("regularize.jl") end
@time @safetestset "truncnorm" begin include("truncnorm.jl") end
@time @safetestset "optim" begin include("optim.jl") end
@time @safetestset "partition" begin include("partition.jl") end
@time @safetestset "ais" begin include("ais.jl") end
@time @safetestset "metropolis" begin include("metropolis.jl") end

@time @safetestset "zerosum" begin include("gauge/zerosum.jl") end
@time @safetestset "rescale_hidden" begin include("gauge/rescale_hidden.jl") end
@time @safetestset "centered_gradient" begin include("centered_gradient.jl") end
@time @safetestset "pcd" begin include("pcd.jl") end

if Sys.islinux() # NPZ has issues on Windows
    # compare to https://github.com/jertubiana/PGM
    @time @safetestset "pgm" begin include("compare_to_pgm/pgm.jl") end
end

ext: jl


text:
#
# General helper functionality
#

#
# Contracting index sets and getting costs
#

# TODO: make a type:
# SBitSet{T <: Unsigned}
#   data::T
# end
#
# That is like a BitSet/BitVector but with a maximum set size.
# Aliases such as:
#
# const SBitSet64 = SBitSet{UInt64}
#
# would be helpful to specify a BitSet with a maximum number of
# 64 elements in the set.
# (See https://discourse.julialang.org/t/parse-an-array-of-bits-bitarray-to-an-integer/42361/11).

# Previously we used the definition in NDTensors:
#import NDTensors: dim
import ITensors: dim

# `is` could be Vector{Int} for BitSet
function dim(is::IndexSetT, ind_dims::Vector) where {IndexSetT<:Union{Vector{Int},BitSet}}
  isempty(is) && return one(eltype(inds_dims))
  dim = one(eltype(ind_dims))
  for i in is
    dim *= ind_dims[i]
  end
  return dim
end

function dim(is::Unsigned, ind_dims::Vector{DimT}) where {DimT}
  _isemptyset(is) && return one(eltype(ind_dims))
  dim = one(eltype(ind_dims))
  i = 1
  @inbounds while !iszero(is)
    if isodd(is)
      # TODO: use Base.Checked.mul_with_overflow
      dim, overflow = Base.Checked.mul_with_overflow(dim, ind_dims[i])
      overflow && return zero(DimT)
    end
    is = is >> 1
    i += 1
  end
  return dim
end

function contraction_cost(indsTᵃ::BitSet, indsTᵇ::BitSet, dims::Vector)
  indsTᵃTᵇ = _symdiff(indsTᵃ, indsTᵇ)
  dim_a = dim(indsTᵃ, dims)
  dim_b = dim(indsTᵇ, dims)
  dim_ab = dim(indsTᵃTᵇ, dims)
  # Perform the sqrt first to avoid overflow.
  # Alternatively, use a larger integer type.
  cost = round(Int, sqrt(dim_a) * sqrt(dim_b) * sqrt(dim_ab))
  return cost, indsTᵃTᵇ
end

function contraction_cost(
  indsTᵃ::IndexSetT, indsTᵇ::IndexSetT, dims::Vector
) where {IndexSetT<:Unsigned}
  unionTᵃTᵇ = _union(indsTᵃ, indsTᵇ)
  cost = dim(unionTᵃTᵇ, dims)
  indsTᵃTᵇ = _setdiff(unionTᵃTᵇ, _intersect(indsTᵃ, indsTᵇ))
  return cost, indsTᵃTᵇ
end

#
# Convert indices into unique integer labels
#

function contraction_labels!(labels1, labels2, is1, is2)
  nextlabel = 1
  nextlabel = common_contraction_labels!(labels1, labels2, is1, is2, nextlabel)
  nextlabel = uncommon_contraction_labels!(labels1, is1, nextlabel)
  nextlabel = uncommon_contraction_labels!(labels2, is2, nextlabel)
  return labels1, labels2
end

# Compute the common contraction labels and return the next label
function common_contraction_labels!(labels1, labels2, is1, is2, label)
  N1 = length(is1)
  N2 = length(is2)
  @inbounds for n1 in 1:N1, n2 in 1:N2
    i1 = is1[n1]
    i2 = is2[n2]
    if i1 == i2
      labels1[n1] = labels2[n2] = label
      label += 1
    end
  end
  return label
end

function uncommon_contraction_labels!(labels, is, label)
  N = length(labels)
  @inbounds for n in 1:N
    if iszero(labels[n])
      labels[n] = label
      label += 1
    end
  end
  return label
end

function contraction_labels!(labels, is)
  ntensors = length(is)
  nextlabel = 1
  # Loop through each tensor pair searching for
  # common indices
  @inbounds for n1 in 1:(ntensors - 1), n2 in (n1 + 1):ntensors
    nextlabel = common_contraction_labels!(
      labels[n1], labels[n2], is[n1], is[n2], nextlabel
    )
  end
  @inbounds for n in 1:ntensors
    nextlabel = uncommon_contraction_labels!(labels[n], is[n], nextlabel)
  end
  return nextlabel - 1
end

function empty_labels(is::NTuple{N}) where {N}
  return ntuple(n -> fill(0, length(is[n])), Val(N))
end

function empty_labels(is::Vector)
  ntensors = length(is)
  labels = Vector{Vector{Int}}(undef, ntensors)
  @inbounds for n in 1:ntensors
    labels[n] = fill(0, length(is[n]))
  end
  return labels
end

function contraction_labels(is)
  labels = empty_labels(is)
  contraction_labels!(labels, is)
  return labels
end

contraction_labels(is...) = contraction_labels(is)

#
# Use a Dict as a cache to map the indices to the integer label
# This only helps with many nodes/tensors (nnodes > 30)
# TODO: determine the crossover when this is useful and use
# it in `depth_first_constructive`/`breadth_first_constructive`
#

contraction_labels_caching(is) = contraction_labels_caching(eltype(eltype(is)), is)

function contraction_labels_caching(::Type{IndexT}, is) where {IndexT}
  labels = empty_labels(is)
  return contraction_labels_caching!(labels, IndexT, is)
end

function contraction_labels_caching!(labels, ::Type{IndexT}, is) where {IndexT}
  N = length(is)
  ind_to_label = Dict{IndexT,Int}()
  label = 0
  @inbounds for n in 1:N
    isₙ = is[n]
    labelsₙ = labels[n]
    @inbounds for j in 1:length(labelsₙ)
      i = isₙ[j]
      i_label = get!(ind_to_label, i) do
        label += 1
      end
      labelsₙ[j] = i_label
    end
  end
  return label
end

#
# Compute the labels and also return a data structure storing the dims.
#

function label_dims(::Type{DimT}, is) where {DimT<:Integer}
  labels = empty_labels(is)
  nlabels = contraction_labels!(labels, is)
  dims = fill(zero(DimT), nlabels)
  @inbounds for i in 1:length(is)
    labelsᵢ = labels[i]
    isᵢ = is[i]
    @inbounds for n in 1:length(labelsᵢ)
      lₙ = labelsᵢ[n]
      if iszero(dims[lₙ])
        dims[lₙ] = dim(isᵢ[n])
      end
    end
  end
  return labels, dims
end

label_dims(is...) = label_dims(is)

# Convert a contraction sequence in pair form to tree format.
# This is used in `depth_first_constructive` to convert the output.
function pair_sequence_to_tree(pairs::Vector{Pair{Int,Int}}, N::Int)
  trees = Any[1:N...]
  for p in pairs
    push!(trees, Any[trees[p[1]], trees[p[2]]])
  end
  return trees[end]
end

#
# BitSet utilities
#

function _cmp(A::BitSet, B::BitSet)
  for (a, b) in zip(A, B)
    if !isequal(a, b)
      return isless(a, b) ? -1 : 1
    end
  end
  return cmp(length(A), length(B))
end

# Returns true when `A` is less than `B` in lexicographic order.
_isless(A::BitSet, B::BitSet) = _cmp(A, B) < 0

bitset(::Type{BitSet}, ints) = BitSet(ints)

function bitset(::Type{T}, ints) where {T<:Unsigned}
  set = zero(T)
  u = one(T)
  for i in ints
    set |= (u << (i - 1))
  end
  return set
end

# Return a vector of the positions of the nonzero bits
# Used for debugging
function findall_nonzero_bits(i::Unsigned)
  nonzeros = Int[]
  n = 1
  @inbounds while !iszero(i)
    if isodd(i)
      push!(nonzeros, n)
    end
    i = i >> 1
    n += 1
  end
  return nonzeros
end

# Return the position of the first nonzero bit
function findfirst_nonzero_bit(i::Unsigned)
  n = 0
  @inbounds while !iszero(i)
    if isodd(i)
      return n + 1
    end
    i = i >> 1
    n += 1
  end
  return n
end

_isless(s1::T, s2::T) where {T<:Unsigned} = s1 < s2
_intersect(s1::BitSet, s2::BitSet) = intersect(s1, s2)
_intersect(s1::T, s2::T) where {T<:Unsigned} = s1 & s2
_union(s1::BitSet, s2::BitSet) = union(s1, s2)
_union(s1::T, s2::T) where {T<:Unsigned} = s1 | s2
_setdiff(s1::BitSet, s2::BitSet) = setdiff(s1, s2)
_setdiff(s1::T, s2::T) where {T<:Unsigned} = s1 & (~s2)
_symdiff(s1::BitSet, s2::BitSet) = symdiff(s1, s2)
_symdiff(s1::T, s2::T) where {T<:Unsigned} = xor(s1, s2)
_isemptyset(s::BitSet) = isempty(s)
_isemptyset(s::Unsigned) = iszero(s)

# TODO: use _first instead, optimize to avoid using _set
_only(s::BitSet) = only(s)
_only(s::Unsigned) = findfirst_nonzero_bit(s)

#
# Adjacency matrix and connected components
#

# For a network of tensors T (stored as index labels), return the adjacency matrix.
function adjacencymatrix(T::Vector, alldims::Vector)
  # First break up the network into disconnected parts
  N = length(T)
  _adjacencymatrix = falses(N, N)
  for nᵢ in 1:(N - 1), nⱼ in (nᵢ + 1):N
    if dim(_intersect(T[nᵢ], T[nⱼ]), alldims) > 1
      _adjacencymatrix[nᵢ, nⱼ] = _adjacencymatrix[nⱼ, nᵢ] = true
    end
  end
  return _adjacencymatrix
end

# For a given adjacency matrix of size n x n, connectedcomponents returns
# a list of components that contains integer vectors, where every integer
# vector groups the indices of the vertices of a connected component of the
# graph encoded by A. The number of connected components is given by
# length(components).
function connectedcomponents(A::AbstractMatrix{Bool})
  n = size(A, 1)
  @assert size(A, 2) == n
  components = Vector{Vector{Int}}(undef, 0)
  assignedlist = falses((n,))
  for i in 1:n
    if !assignedlist[i]
      assignedlist[i] = true
      checklist = [i]
      currentcomponent = [i]
      while !isempty(checklist)
        j = pop!(checklist)
        for k in findall(A[j, :])
          if !assignedlist[k]
            push!(currentcomponent, k)
            push!(checklist, k)
            assignedlist[k] = true
          end
        end
      end
      push!(components, currentcomponent)
    end
  end
  return components
end

# For a network of tensors T (stored as index labels), return the connected components
# (splits up T into the connected components).
function connectedcomponents(T::Vector, alldims::Vector)
  return connectedcomponents(adjacencymatrix(T, alldims))
end


text:
function findindex(x,X)
    # For an array X of increasing values,
    # we find i, sucht that X[i]<=x<X[i+1].
    for i=1:length(X)-1
        if X[i]<= x < X[i+1]
            return i
        end
    end
    error(x," is not between ",X[1]," and ",X[end],".")
    return NaN
end

function mmod(a,b) # Stupid arrays starting at 1....
    return mod(a-1,b)+1
end

function a2sa(a) # takes an array, returns a static array
    if length(a)==0
        return SVector{0,typeof(Float64)}(a)
    else
        return SVector{length(a),typeof(a[1])}(a)
    end
end

function getperiodic(newx,newy,i,j,ishift,jshift,edgesx,edgesy)

    isperiodic = false
    nx = length(edgesx)-1
    ny = length(edgesy)-1
    i += ishift
    j += jshift
    if i<1
        i = nx
        newx = edgesx[end]
        isperiodic = true
    elseif  i>nx
        i = 1
        newx = edgesx[1]
        isperiodic = true
    end

    if j<1
        j = ny
        newy = edgesy[end]
        isperiodic = true
    elseif  j>ny
        j = 1
        newy = edgesy[1]
        isperiodic = true
    end
    return isperiodic,newx,newy, i,j
end


text:
# Testing for gzip

using GZip
using Base.Test

##########################
# test_context("GZip tests")
##########################

#for epoch in 1:10

tmp = mktempdir()

test_infile = "$JULIA_HOME/../share/julia/helpdb.jl"
test_compressed = "$tmp/helpdb.jl.gz"
test_empty = "$tmp/empty.jl.gz"

@windows_only gunzip="gunzip.exe"
@unix_only    gunzip="gunzip"

test_gunzip = true
try
    run(`which $gunzip` |> DevNull)
catch
    test_gunzip = false
end

#########################
# test_group("Compress Test1: gzip.jl")
##########################

data = open(readall, test_infile);

first_char = data[1]

gzfile = gzopen(test_compressed, "wb")
@test write(gzfile, data) == length(data.data)
@test close(gzfile) == Z_OK
@test close(gzfile) != Z_OK

#@test throws_exception(write(gzfile, data), GZError)
@test_throws write(gzfile, data)

if test_gunzip
    data2 = readall(`$gunzip -c $test_compressed`)
    @test data == data2
end

data3 = gzopen(readall, test_compressed)
@test data == data3

# Test gzfdio
raw_file = open(test_compressed, "r")
gzfile = gzdopen(fd(raw_file), "r")
data4 = readall(gzfile)
close(gzfile)
close(raw_file)
@test data == data4

# Test peek
gzfile = gzopen(test_compressed, "r")
@test peek(gzfile) == first_char
readall(gzfile)
@test peek(gzfile) == -1
close(gzfile)

# Screw up the file
raw_file = open(test_compressed, "r+")
seek(raw_file, 3) # leave the gzip magic 2-byte header
write(raw_file, zeros(Uint8, 10))
close(raw_file)

#@test throws_exception(gzopen(readall, test_compressed), GZError)
@test_throws gzopen(readall, test_compressed)

##########################
# test_group("gzip file function tests (writing)")
##########################
gzfile = gzopen(test_compressed, "wb")
write(gzfile, data) == length(data.data)
@test flush(gzfile) == Z_OK

pos = position(gzfile)
@test_throws seek(gzfile, 100)   # can't seek backwards on write
@test position(gzfile) == pos
@test skip(gzfile, 100)
@test position(gzfile) == pos + 100

#@test throws_exception(truncate(gzfile, 100), ErrorException)
#@test throws_exception(seekend(gzfile), ErrorException)
@test_throws truncate(gzfile, 100)
@test_throws seekend(gzfile)

@test close(gzfile) == Z_OK

gzopen(test_empty, "w") do io
    a = "".data
    @test gzwrite(io, pointer(a), length(a)*sizeof(eltype(a))) == int32(0)
end

##########################
# test_group("gzip file function tests (strategy read/write)")
##########################

# rewrite the test file
modes = "fhR "
if GZip.ZLIB_VERSION >= (1,2,5,2)
    modes = "fhRFT "
end
for ch in modes
    if ch == ' '
        ch = ""
    end
    for level = 0:9
        gzfile = gzopen(test_compressed, "wb$level$ch")
        @test write(gzfile, data) == length(data.data)
        @test close(gzfile) == Z_OK

        file_size = filesize(test_compressed)

        #println("wb$level$ch: ", file_size)

        if ch == 'T'
            @test(file_size == length(data.data))
        elseif ch == 'F'
            @test(file_size >= length(data.data))
        elseif level == 0
            @test(file_size > length(data.data))
        else
            @test(file_size < length(data.data))
        end

        # readline test
        gzf = gzopen(test_compressed)
        s = IOBuffer()
        while !eof(gzf)
            write(s, readline(gzf))
        end
        data2 = takebuf_string(s);

        # readuntil test
        seek(gzf, 0)
        while !eof(gzf)
            write(s, readuntil(gzf, 'a'))
        end
        data3 = takebuf_string(s);
        close(gzf)

        @test(data == data2)
        @test(data == data3)

    end
end

##########################
# test_group("gzip array/matrix tests (write/read)")
##########################

let BUFSIZE = 65536
    for level = 0:3:6
        for T in [Int8,Uint8,Int16,Uint16,Int32,Uint32,Int64,Uint64,Int128,Uint128,
                  Float32,Float64,Complex64,Complex128]

            minval = 34567
            try
                minval = min(typemax(T), 34567)
            catch
                # do nothing
            end

            # Ordered array
            b = zeros(T, BUFSIZE)
            if !isa(T, Complex)
                for i = 1:length(b)
                    b[i] = (i-1)%minval;
                end
            else
                for i = 1:length(b)
                    b[i] = (i-1)%minval - (minval-(i-1))%minval * im
                end
            end

            # Random array
            if isa(T, FloatingPoint)
                r = (T)[rand(BUFSIZE)...];
            elseif isa(T, Complex64)
                r = Int32[rand(BUFSIZE)...] + Int32[rand(BUFSIZE)...] * im
            elseif isa(T, Complex128)
                r = Int64[rand(BUFSIZE)...] + Int64[rand(BUFSIZE)...] * im
            else
                r = b[rand(1:BUFSIZE, BUFSIZE)];
            end

            # Array file
            b_array_fn = "$tmp/b_array.raw.gz"
            r_array_fn = "$tmp/r_array.raw.gz"

            gzaf_b = gzopen(b_array_fn, "w$level")
            write(gzaf_b, b)
            close(gzaf_b)

            #println("$T ($level) ordered: $(filesize(b_array_fn))")

            gzaf_r = gzopen(r_array_fn, "w$level")
            write(gzaf_r, r)
            close(gzaf_r)

            #println("$T ($level) random: $(filesize(r_array_fn))")

            b2 = zeros(T, BUFSIZE)
            r2 = zeros(T, BUFSIZE)

            b2_infile = gzopen(b_array_fn)
            read(b2_infile, b2);
            close(b2_infile)

            r2_infile = gzopen(r_array_fn)
            read(r2_infile, r2);
            close(r2_infile)

            @test b == b2
            @test r == r2
        end
    end
end

##########################
# test_group("gzip unicode tests (write/read)")
##########################

unicode_gz_file = "$tmp/unicode_test.gz"

# Use perl to generate UTF-32BE unicode data a la unicode.jl, then convert to UTF-32LE, UTF-8 via iconv
UTF32BE_path = joinpath(tmp,"UTF32BE.unicode")
UTF32LE_path = joinpath(tmp,"UTF32LE.unicode")
UTF8_path = joinpath(tmp,"UTF8.unicode")
run(`perl -e 'print pack "N*", 0xfeff, 0..0xd7ff, 0xe000..0x10ffff' ` |> UTF32BE_path)
run(`iconv -f UTF-32BE -t UTF-32LE $UTF32BE_path` |> UTF32LE_path)
run(`iconv -f UTF-32BE -t UTF-8 $UTF32BE_path` |> UTF8_path)

str1 = CharString(reinterpret(Char, read(open(UTF32LE_path), Uint32, 1112065)[2:]));
str2 = UTF8String(read(open(UTF8_path), Uint8, 4382595)[4:]);

UTF32LE_gz = gzopen(unicode_gz_file, "w")
write(UTF32LE_gz, str1)
close(UTF32LE_gz)

str1b = readall(`gunzip -c $unicode_gz_file`);
str1c = gzopen(readall, unicode_gz_file);
@test str1 == str1b
@test str1 == str1c

UTF8_gz = gzopen(unicode_gz_file, "w");
write(UTF8_gz, str2)
close(UTF8_gz)

str2b = readall(`gunzip -c $unicode_gz_file`);
str2c = gzopen(readall, unicode_gz_file);
@test str2 == str2b
@test str2 == str2c

run(`rm -Rf $tmp`)

#end  # for epoch


text:
"""
    VelocityFields(arch, grid; [u, v, w])

Return a NamedTuple with fields `u`, `v`, `w` initialized on the architecture `arch`
and `grid`. Fields may be passed via the optional keyword arguments `u`, `v`, and `w`.
"""
function VelocityFields(arch, grid; u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)),
                                    v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)),
                                    w = ZFaceField(arch, grid, WVelocityBoundaryConditions(grid)))
    return (u=u, v=v, w=w)
end

"""
    VelocityFields(arch, grid, bcs::NamedTuple)

Return a NamedTuple with fields `u`, `v`, `w` initialized on the architecture `arch`
and `grid`. Boundary conditions `bcs` may be specified via a named tuple of
`FieldBoundaryCondition`s.
"""
function VelocityFields(arch, grid, bcs::NamedTuple)
    u_bcs = :u ∈ keys(bcs) ? bcs[:u] : UVelocityBoundaryConditions(grid)
    v_bcs = :v ∈ keys(bcs) ? bcs[:v] : VVelocityBoundaryConditions(grid)
    w_bcs = :w ∈ keys(bcs) ? bcs[:w] : WVelocityBoundaryConditions(grid)

    u = XFaceField(arch, grid, u_bcs)
    v = YFaceField(arch, grid, v_bcs)
    w = ZFaceField(arch, grid, w_bcs)

    return (u=u, v=v, w=w)
end

"""
    TracerFields(arch, grid, tracer_names; kwargs...)

Return a NamedTuple with tracer fields specified by `tracer_names` initialized as
`CellField`s on the architecture `arch` and `grid`. Fields may be passed via optional
keyword arguments `kwargs` for each field.

# Examples
```julia
arch = CPU()
topology = (Periodic, Periodic, Bounded)
grid = RegularCartesianGrid(topology=topology, size=(16, 16, 16), size=(1, 2, 3))
tracers = (:T, :S, :random)
noisy_field = CellField(arch, grid, TracerBoundaryConditions(grid), randn(16, 16))
tracer_fields = TracerFields(arch, grid, tracers, random=noisy_field)
```
"""
function TracerFields(arch, grid, tracer_names; kwargs...)
    tracer_fields =
        Tuple(c ∈ keys(kwargs) ?
              kwargs[c] :
              CellField(arch, grid, TracerBoundaryConditions(grid))
              for c in tracer_names)
    return NamedTuple{tracer_names}(tracer_fields)
end

"""
    TracerFields(arch, grid, tracer_names, bcs::NamedTuple)

Return a NamedTuple with tracer fields specified by `tracer_names` initialized as
`CellField`s on the architecture `arch` and `grid`. Boundary conditions `bcs` may
be specified via a named tuple of `FieldBoundaryCondition`s.
"""
function TracerFields(arch, grid, tracer_names, bcs::NamedTuple)
    tracer_fields =
        Tuple(c ∈ keys(bcs) ?
              CellField(arch, grid, bcs[c]) :
              CellField(arch, grid, TracerBoundaryConditions(grid))
              for c in tracer_names)
    return NamedTuple{tracer_names}(tracer_fields)
end

TracerFields(arch, grid, ::Union{Tuple{}, Nothing}; kwargs...) = NamedTuple()
TracerFields(arch, grid, tracer::Symbol; kwargs...) = TracerFields(arch, grid, tuple(tracer); kwargs...)
TracerFields(arch, grid, tracers::NamedTuple; kwargs...) = tracers

tracernames(::Nothing) = ()
tracernames(name::Symbol) = tuple(name)
tracernames(names::NTuple{N, Symbol}) where N = :u ∈ names ? names[4:end] : names
tracernames(::NamedTuple{names}) where names = tracernames(names)

"""
    PressureFields(arch, grid; [pHY′, pNHS])

Return a NamedTuple with pressure fields `pHY′` and `pNHS` initialized as
`CellField`s on the architecture `arch` and `grid`. Fields may be passed via the
optional keyword arguments `pHY′` and `pNHS`.
"""
function PressureFields(arch, grid; pHY′ = CellField(arch, grid, PressureBoundaryConditions(grid)),
                                    pNHS = CellField(arch, grid, PressureBoundaryConditions(grid)))

    return (pHY′=pHY′, pNHS=pNHS)
end

"""
    PressureFields(arch, grid, bcs::NamedTuple)

Return a NamedTuple with pressure fields `pHY′` and `pNHS` initialized as
`CellField`s on the architecture `arch` and `grid`.  Boundary conditions `bcs` may
be specified via a named tuple of `FieldBoundaryCondition`s.
"""
function PressureFields(arch, grid, bcs::NamedTuple)
    pHY′_bcs = :pHY′ ∈ keys(bcs) ? bcs[:pHY′] : PressureBoundaryConditions(grid)
    pNHS_bcs = :pNHS ∈ keys(bcs) ? bcs[:pNHS] : PressureBoundaryConditions(grid)

    pHY′ = CellField(arch, grid, pHY′_bcs)
    pNHS = CellField(arch, grid, pNHS_bcs)

    return (pHY′=pHY′, pNHS=pNHS)
end

"""
    TendencyFields(arch, grid, tracer_names; kwargs...)

Return a NamedTuple with tendencies for all solution fields (velocity fields and
tracer fields), initialized on the architecture `arch` and `grid`. Optional `kwargs`
can be specified to assign data arrays to each tendency field.
"""
function TendencyFields(arch, grid, tracer_names;
                        u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)),
                        v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)),
                        w = ZFaceField(arch, grid, WVelocityBoundaryConditions(grid)),
                        kwargs...)

    velocities = (u=u, v=v, w=w)
    tracers = TracerFields(arch, grid, tracer_names; kwargs...)
    return merge(velocities, tracers)
end


text:
"""
    NaNChecker{F} <: AbstractDiagnostic

A diagnostic that checks for `NaN` values and aborts the simulation if any are found.
"""
struct NaNChecker{F} <: AbstractDiagnostic
    frequency :: Int
       fields :: F
end

"""
    NaNChecker(model; frequency, fields)

Construct a `NaNChecker` for `model`. `fields` should be a `Dict{Symbol,Field}`. A
`frequency` should be passed to indicate how often to check for NaNs (in number of
iterations).
"""
function NaNChecker(model; frequency, fields)
    return NaNChecker(frequency, fields)
end

function run_diagnostic(model, nc::NaNChecker)
    for (name, field) in nc.fields
        if any(isnan, field.data.parent)
            t, i = model.clock.time, model.clock.iteration
            error("time = $t, iteration = $i: NaN found in $name. Aborting simulation.")
        end
    end
end


text:
using LazySets, MathematicalPredicates, SparseArrays

function mna5_specification()
    # initial set: xᵢ ∈ [0.0002, 0.00025] if i ≤ 10 and xᵢ = 0 otherwise
    X0 = Hyperrectangle(low=[fill(0.0002, 10); zeros(10903)],
                        high=[fill(0.00025, 10); zeros(10903)])

    # property: x1 ≤ 0.2 && x2 ≤ 0.15
    property = is_contained_in(HPolyhedron([
        HalfSpace(sparsevec([1], [1.0], 10913), 0.2),
        HalfSpace(sparsevec([2], [1.0], 10913), 0.15)]))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using Test

@testset "iroot" begin
    for T in (Int32, Int64, BigInt)
        @test iroot(T(100), 2)   == T(10)
        @test iroot(T(101), 2)   == T(10)
        @test iroot(T(99),  2)   == T(9)
        @test iroot(T(1000), 3)  == T(10)
        @test iroot(T(1001), 3)  == T(10)
        @test iroot(T(999),  3)  == T(9)
        @test iroot(T(10000), 4) == T(10)
        @test iroot(T(10001), 4) == T(10)
        @test iroot(T(9999),  4) == T(9)
    end
    @test iroot(big(23)^50, 50)     == big(23)
    @test iroot(big(23)^50 + 1, 50) == big(23)
    @test iroot(big(23)^50 - 1, 50) == big(22)
end

@testset "ispower" begin
    for T in (Int32, Int64, BigInt)
        @test ispower(T(100))    == true
        @test ispower(T(1))      == true
        @test ispower(T(0))      == true
        @test ispower(T(12))     == false
        @test ispower(T(2^30))   == true
        @test ispower(T(5)^10)   == true
        @test ispower(T(2^30)+1) == false
        @test ispower(T(5)^10+1) == false
    end
    @test ispower(big(5)^40)     == true
    @test ispower(big(5)^40 + 1) == false
    @test ispower(6 * big(5)^40) == false
end

@testset "find_exponent" begin
    for T in (Int32, Int64, BigInt)
        @test find_exponent(T(100))    === 2
        @test find_exponent(T(1))      === 1
        @test find_exponent(T(0))      === 1
        @test find_exponent(T(12))     === 1
        @test find_exponent(T(2^30))   === 30
        @test find_exponent(T(5)^10)   === 10
        @test find_exponent(T(2^30)+1) === 1
        @test find_exponent(T(5)^10+1) === 1
    end
    @test find_exponent(big(5)^40)     === 40
    @test find_exponent(big(5)^40 + 1) === 1
    @test find_exponent(6*big(5)^40 )  === 1
end

@testset "is_probably_prime" begin
    for T in (Int32, Int64, BigInt)
        @test is_probably_prime(T(2)^7-1)  == true
        @test is_probably_prime(T(2)^13-1) == true
        @test is_probably_prime(T(2)^19-1) == true
        @test is_probably_prime(T(2)^27-1) == false
        @test is_probably_prime(T(2)^23-1) == false
        @test is_probably_prime(T(2)^30-1) == false
    end
    @test is_probably_prime(big(2)^127-1) == true
    @test is_probably_prime(big(2)^128-1) == false
end


text:
import UncertainData:
    AbstractUncertainValue,
    AbstractUncertainValueDataset,
    AbstractUncertainIndexValueDataset,
    resample

function get_random_segment(source, target, seq::RandomSequences)
    seq_length = seq.sequence_length
    N = length(source)
    if seq_length isa Int
        idx_start = rand(1:(length(source) - seq_length))
        r = idx_start:(idx_start + seq_length)
    elseif seq_length isa AbstractVector{Int}
        # seq_length is now something we can sample directly from
        seqlen = rand(seq_length)
        idx_start = rand(1:(N + 1 - seqlen))
        r = idx_start:(idx_start + seqlen - 1)
    else
        throw(ArgumentError("`resampling.sequence_length`must be an integer or a collection of integers"))
    end

    return r
end

"""
    resample_and_subset(x, r::AbstractVector{Int})

Resample `x` assuming the points in `x` are independent,
then return the points with indices `r`.
"""
function resample_and_subset(x, r::AbstractVector{Int})
    v = zeros(Float64, length(r))

    if x isa AbstractVector
        v[:] = x[r]
    elseif typeof(x) <: Union{Vector{<:AbstractUncertainValue}, AbstractUncertainValueDataset}
        v[:] = resample(x)[r]
    elseif typeof(x) <: AbstractUncertainIndexValueDataset
        v[:] = resample(x)[2][r]
    end

    return v
end

"""
    causality(source, target, test::RandomSequencesTest)

Apply a causality test on random, consecutive chunks of `x` and `y`.

The chunk length can be fixed (an integer) or a collection of
chunk lengths.

Works on any inputs `x` and `y`, both scalar vectors and
uncertain datasets.

## Example

For scalar valued vectors:

```julia
x, y = rand(100), rand(100)

# Use a cross-mapping test with default parameters (not recommended,
# make sure you understand the parameters!)
cm_test = CrossMappingTest()

# A random sequences test that applies the cross mapping test
# to 150 different random chunks of the data of lengths 80 to 90
n_realisations, chunk_lengths = 150, 80:90
rs = RandomSequences(n_realisations, chunk_lengths)

# Combine to a RandomSequencesTest
rs_test = RandomSequencesTest(cm_test, chunk_lengths)

# Compute causality statistic
causality(x, y, rs_test)
```

This also works on uncertain data, or any combination of scalar vectors
and uncertain data:

```julia
# Some example data
N = 300
sys = ar1_unidir(c_xy = 0.8)
X, Y = example_uncertain_indexvalue_datasets(sys, N, (1, 2),
    d_xval = Uniform(0.001, 0.05), d_yval = Uniform(0.001, 0.05));

# Apply cross mapping test on 100 different randomly selected chunks
# of lengths 75 to 90
r = RandomSequences(100, 75:90)
rs_test = RandomSequencesTest(CrossMappingTest(), r)
causality(X, Y, rs_test)
```
"""
function causality(
        source::Union{AbstractVector, Vector{<:AbstractUncertainValue}, AbstractUncertainValueDataset, AbstractUncertainIndexValueDataset},
        target::Union{AbstractVector, Vector{<:AbstractUncertainValue}, AbstractUncertainValueDataset, AbstractUncertainIndexValueDataset},
        test::RandomSequencesTest)

    n_sequences = test.sequences_resampling.n

    # Get the return type of the causality test
    causality_test = test.test
    RT = get_return_type(causality_test)

    results = Vector{RT}(undef, n_sequences)

    for i = 1:n_sequences
        r = get_random_segment(source, target, test.sequences_resampling)
        source_segment = resample_and_subset(source, r)
        target_segment = resample_and_subset(target, r)

        results[i] = causality(source_segment, target_segment, causality_test)
    end

    return results
end

function causality(
        source::Union{AbstractVector, Vector{<:AbstractUncertainValue}, AbstractUncertainValueDataset, AbstractUncertainIndexValueDataset},
        target::Union{AbstractVector, Vector{<:AbstractUncertainValue}, AbstractUncertainValueDataset, AbstractUncertainIndexValueDataset},
        cond::Union{AbstractVector, Vector{<:AbstractUncertainValue}, AbstractUncertainValueDataset, AbstractUncertainIndexValueDataset},
        test::RandomSequencesTest)

    n_sequences = test.sequences_resampling.n

    # Get the return type of the causality test
    causality_test = test.test
    RT = get_return_type(causality_test)

    results = Vector{RT}(undef, n_sequences)

    for i = 1:n_sequences
        r = get_random_segment(source, target, test.sequences_resampling)
        source_segment = resample_and_subset(source, r)
        target_segment = resample_and_subset(target, r)
        cond_segment = resample_and_subset(cond, r)

        results[i] = causality(source_segment, target_segment, cond_segment, causality_test)
    end

    return results
end

export causality, resample_and_subset


text:
ocpi2_psys = PiezoSystem(retrieve_amplifier("30DV300"),
                         retrieve_positioner("NanoSX800"),
                                                  retrieve_pid("ocpi-2"))

register_piezosystem(ocpi2_psys, "ocpi-2")


text:
module PkgCompiler

using Pkg
using Libdl
using ProgressMeter

pkg_compiler_dir = "$(homedir())/.julia/pkgcompiler/v$(VERSION.major).$(VERSION.minor)"

global_snoop_file =  joinpath(pkg_compiler_dir, "snoop-file.jl")
snoop_file = nothing
snoop_file_io = nothing

function __init__()
    if haskey(ENV, "PkgCompilerEnabled") && ENV["PkgCompilerEnabled"] == "1"
        start_snooping()
    end
end

function start_snooping()
    global snoop_file_io
    if snoop_file_io === nothing
        global snoop_file = "$(tempname())-snoop.jl"
        mkpath(dirname(snoop_file))
        @info("Snooping -> $(snoop_file)")
        global snoop_file_io = open(snoop_file, "w")
        ccall(:jl_dump_compiles, Cvoid, (Ptr{Cvoid},), snoop_file_io.handle)
        atexit(finish_snooping)
    else
        @warn("Snooping is already running -> $(snoop_file)")
    end
end

# This function is for testing purposes
function set_global_snoop_file(file_name::AbstractString)
    global global_snoop_file = file_name
end

# This function is for testing purposes
function set_pkg_compiler_dir(new_diectory::AbstractString)
    global pkg_compiler_dir = new_diectory
end

function _flush()
    if snoop_file_io !== nothing
        flush(snoop_file_io)
    end
end

function _update_startup_file(content::AbstractString)
    begin_text = "Begin of PkgCompiler section"
    end_text = "End of PkgCompiler section"

    file_name = "$(homedir())/.julia/config/startup.jl"
    # Create new file if not existing
    if !isfile(file_name)
        @info("No startup file ($(file_name) will be created)")
        mkpath(dirname(file_name))
        touch(file_name)
    end

    lines = readlines(open(file_name))
    l = length(lines)

    # Test if PkgCompiler is already registered in ~/.julia/config/startup.jl
    start_idx = findfirst(x -> occursin(begin_text, x), lines)
    end_idx = findlast(x -> occursin(end_text, x), lines)

    if start_idx === nothing && end_idx === nothing
        # new text will be added at the end of the startup file
        start_idx = end_idx = l+1
    end
    if start_idx === nothing || end_idx === nothing
        @error("The startup file is malformed ($(file_name))")
        return
    end

    # Replace the old part from PkgCompiler
    open(file_name, "w") do file
        for i in 1:min(start_idx-1, l)
            write(file, lines[i], "\n")
        end
        if start_idx == l+1
            write(file, "\n")
        end
        write(file, "# $(begin_text) #####\n")
        write(file, content)
        write(file, "# $(end_text) #####\n")
        for i in (end_idx+1):l
            write(file, lines[i], "\n")
        end
    end
end

function register_snooping()
    content  = "ENV[\"PkgCompilerEnabled\"] = 1\n"
    content *= "try ; using PkgCompiler ; catch e ; end\n"
    _update_startup_file(content)
    __init__()
end

function deregister_snooping()
    _update_startup_file("\n")
    global snoop_file_io
    ccall(:jl_dump_compiles, Cvoid, (Ptr{Cvoid},), C_NULL)
    close(snoop_file_io)
    snoop_file_io = nothing
end

function finish_snooping()
    global snoop_file_io
    try
        if snoop_file_io === nothing
            @warn("No active snooping file")
            return
        end
        ccall(:jl_dump_compiles, Cvoid, (Ptr{Cvoid},), C_NULL)
        close(snoop_file_io)
        snoop_file_io = nothing

        save_statements()
        rm("$(snoop_file)")
    catch e
        @warn(e)
    end
end

function save_statements()
    try
        _flush()

        lines = readlines("$(snoop_file)")
        open("$(snoop_file).txt", "w") do file
            for line in lines
                sp = split(line, "\t")
                if length(sp) == 2
                    function_str = sp[2][2:end-1]
                    write(file, function_str * "\n")
                end
            end
        end

        @info("Copy snooped function into $(global_snoop_file)")
        if isfile(global_snoop_file)
            run(`sort $(snoop_file).txt $(global_snoop_file) -o $(global_snoop_file).tmp`)
        else
            mkpath(pkg_compiler_dir)
            run(`sort $(snoop_file).txt -o $(global_snoop_file).tmp`)
        end
        run(pipeline(`uniq $(global_snoop_file).tmp`, stdout="$(global_snoop_file)"))
        rm("$(snoop_file).txt")
        rm("$(global_snoop_file).tmp")
    catch e
        @warn(e)
    end
end

# inspired by https://nextjournal.com/sdanisch/packagecompiler-talk
function build_system_image()
    save_statements()

    lib = abspath(joinpath(Sys.BINDIR, Base.PRIVATE_LIBDIR))
    img = unsafe_string(Base.JLOptions().image_file)

    systemp = joinpath(pkg_compiler_dir, "sys.a")
    sysout = joinpath(pkg_compiler_dir, "sys.$(Libdl.dlext)")
    script = joinpath(@__DIR__, "compile.jl")

    # compile image
    @info("Compile $(@__DIR__)/compile.jl -> $(systemp)")
    @info("Using snoop file: $(global_snoop_file)")
    @time run(`julia -J$(img) --color=yes --startup-file=no -O3 -g 2 --output-o $(systemp) $(script) $(global_snoop_file)`)

    # From PackageCompiler.jl
    # `cc -shared '-DJULIAC_PROGRAM_LIBNAME="/home/petr/.julia/dev/PackageCompiler/sysimg/sys.so"' -o /home/petr/.julia/dev/PackageCompiler/sysimg/sys.so -Wl,--whole-archive /home/petr/.julia/dev/PackageCompiler/sysimg/sys.a -Wl,--no-whole-archive -std=gnu99 -I/home/petr/software/julia/julia-1.2.0/include/julia -DJULIA_ENABLE_THREADING=1 -fPIC -L/home/petr/software/julia/julia-1.2.0/lib -Wl,--export-dynamic -Wl,-rpath,/home/petr/software/julia/julia-1.2.0/lib -Wl,-rpath,/home/petr/software/julia/julia-1.2.0/lib/julia -ljulia -m64 -O3`

    # Link it
    @info("Link system image")
    # Prevent compiler from stripping all symbols from the shared lib.
    if Sys.isapple()
        temp_part = `-all_load $systemp`
    else
        temp_part = `-Wl,--whole-archive $systemp -Wl,--no-whole-archive`
    end
    part_end = ``
    if Sys.isapple()
        part_end = ` -Wl,-install_name,@rpath/$sysout`
    elseif Sys.iswindows()
        part_end = ` -Wl,--export-all-symbols`
    end
    @time run(`c++ -m64 -shared -DJULIAC_PROGRAM_LIBNAME=\"$sysout\" -fPIC -L$(lib) -L$(dirname(lib)) -o $(sysout) $(temp_part) $(part_end) `);
    @info("system image is prepared, you can run it by:\njulia -J $(sysout)")
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 0e4dcd02-71f8-11eb-2ac3-612951cdb3d4
md"""
# Examen scheduling
The purpose of the first project is to create an exam schedule. This is a group work where the groups are in competition with each other. The evaluation is made on the basis of:
* the correctness of the proposed solution
* the computational cost
* the functionality
* the scalability
* the generalizability

If the results meet expectations, your work can also be used effectively in the future for planning exams or class schedules.

The absolute minimum requirements (or constraints) that must be respected are
* the availability of professors (and teaching assistants)
* the non-overlapping of subjects

possible additions are:
* distinction between an oral or written exam
* whether or not it is compulsory to conjoin several examination days for the same course if it runs over several days
* courses that are evaluated in several parts (possibly by several professors)
* preferences (e.g. certain minimum number of days of preparation for a certain course, taking into account student's input)

You may consider this a constraint satisfaction problem. Different way of solving it can and may be used (backtracking <> local search etc.).

Start with one year, then extend the functionality to one faculty and finally to both faculties.
"""

# ╔═╡ Cell order:
# ╟─0e4dcd02-71f8-11eb-2ac3-612951cdb3d4

max_stars_repo_path: Projects/Project 1 - schudeling.jl

max_issues_repo_path: Projects/Project 1 - schudeling.jl

max_forks_repo_path: Projects/Project 1 - schudeling.jl


text:
module SATModel

# const Name = Symbol

# struct Atom
#   name::Name
# end

# struct NotLit
#   name::Name
# end

# const Literal = Union{Var, NotLit}

# struct Clause
#   lits::Set{Literal}
# end

# Clause(lits::Literal...) = Clause(Set{Literal}(lits))
# atoms(clause::Clause) = map(atom, clause.lits)

# struct CNF
#   clauses::Vector{Clause}
# end

# CNF(clause::Clause...) = CNF([clause...])
# atoms(cnf::CNF) = union(map(atoms, cnf.clauses))

# "Generate a function that evaluates whether `cnf` formula is true"
# function genfunc(cnf)
#   atoms(cnf)

# end

# Base.:!(v::Var) = NotVar(v.name)
# Base.:!(v::NotVar) = Var(v.name)

# ## Example

# function test()
#   A = Var(:A)
#   B = Var(:B)
#   C = Var(:C)
#   D = Var(:D)
#   E = Var(:E)
#   F = Var(:F)
#   c1 = Clause(A, !B, !C)
#   c2 = Clause(!D, E, F)
#   cnf = CNF(c1, c2)
# end

A ∧ B = A & B
A ∨ B = A | B

function ex1(A, B, C, D, E, F)
  (A ∨ !B ∨ !C) ∧ (!D ∨ E ∨ F)
end

end # module


text:
using BPAlign

n = 100
λ = 3
s = 0.55
PG = Pair_ER(n,λ,s)

@time mp_fin = run_bp(PG,5)

M_lr = create_matrix_lr(PG,mp_fin)
M_loglr = log.(M_lr)

ov1, ov2, v1,v2 = eval_M(PG, M_loglr)
println(ov1,ov2)


text:
@doc raw"""
`ParametersCGVI`: Parameters for right-hand side function of continuous Galerkin variational Integrator.

### Parameters

* `Θ`: function of the noncanonical one-form (∂L/∂v)
* `f`: function of the force (∂L/∂q)
* `Δt`: time step
* `b`: weights of the quadrature rule
* `c`: nodes of the quadrature rule
* `x`: nodes of the basis
* `m`: mass matrix
* `a`: derivative matrix
* `r₀`: reconstruction coefficients at the beginning of the interval
* `r₁`: reconstruction coefficients at the end of the interval
* `t`: initial time
* `q`: solution of q at time t
* `p`: solution of p at time t
"""
mutable struct ParametersCGVI{DT, TT, D, S, R, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    Δt::TT

    b::Vector{TT}
    c::Vector{TT}

    x::Vector{TT}
    m::Matrix{TT}
    a::Matrix{TT}
    r₀::Vector{TT}
    r₁::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersCGVI{DT,D}(equs::ET, Δt::TT, b, c, x, m, a, r₀, r₁) where {DT, TT, D, ET <: NamedTuple}
        new{DT,TT,D,length(x),length(c),ET}(equs, Δt, b, c, x, m, a, r₀, r₁, zero(TT), zeros(DT,D), zeros(DT,D))
    end

    function ParametersCGVI{DT,D}(equs::NamedTuple, Δt::TT, basis::Basis{TT}, quadrature::QuadratureRule{TT}) where {DT,TT,D}
        # compute coefficients
        r₀ = zeros(TT, nbasis(basis))
        r₁ = zeros(TT, nbasis(basis))
        m  = zeros(TT, nnodes(quadrature), nbasis(basis))
        a  = zeros(TT, nnodes(quadrature), nbasis(basis))

        for i in eachindex(basis)
            r₀[i] = basis[zero(TT), i]
            r₁[i] = basis[one(TT), i]
            for j in eachindex(quadrature)
                m[j,i] = basis[nodes(quadrature)[j], i]
                a[j,i] = basis'[nodes(quadrature)[j], i]
            end
        end

        if get_config(:verbosity) > 1
            println()
            println("  Continuous Galerkin Variational Integrator")
            println("  ==========================================")
            println()
            println("    c  = ", nodes(quadrature))
            println("    b  = ", weights(quadrature))
            println("    x  = ", nodes(basis))
            println("    m  = ", m)
            println("    a  = ", a)
            println("    r₀ = ", r₀)
            println("    r₁ = ", r₁)
            println()
        end

        ParametersCGVI{DT,D}(equs, Δt, weights(quadrature), nodes(quadrature), grid(basis), m, a, r₀, r₁)
    end
end

function update_params!(params::ParametersCGVI, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
end

struct IntegratorCacheCGVI{ST,D,S,R} <: IODEIntegratorCache{ST,D}
    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    X::Vector{Vector{ST}}
    Q::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    F::Vector{Vector{ST}}

    function IntegratorCacheCGVI{ST,D,S,R}() where {ST,D,S,R}
        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        X = create_internal_stage_vector(ST,D,S)
        Q = create_internal_stage_vector(ST,D,R)
        P = create_internal_stage_vector(ST,D,R)
        V = create_internal_stage_vector(ST,D,R)
        F = create_internal_stage_vector(ST,D,R)

        new(q̃, p̃, ṽ, f̃, s̃, X, Q, P, V, F)
    end
end

function IntegratorCache{ST}(params::ParametersCGVI{DT,TT,D,S,R}; kwargs...) where {ST,DT,TT,D,S,R}
    IntegratorCacheCGVI{ST,D,S,R}(; kwargs...)
end

@inline CacheType(ST, params::ParametersCGVI{DT,TT,D,S,R}) where {DT,TT,D,S,R} = IntegratorCacheCGVI{ST,D,S,R}

"Continuous Galerkin Variational Integrator."
struct IntegratorCGVI{DT, TT, D, S, R,
                      BT <: Basis,
                      PT <: ParametersCGVI{DT,TT,D,S,R},
                      ST <: NonlinearSolver{DT},
                      IT <: InitialGuessIODE{TT}} <: IODEIntegrator{DT,TT}
    basis::BT
    quadrature::QuadratureRule{TT,R}

    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorCGVI(basis::BT, quadrature::QuadratureRule{TT,R}, params::ParametersCGVI{DT,TT,D,S},
                    solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,BT,ST,IT}
        new{DT, TT, D, S, R, BT, typeof(params), ST, IT}(basis, quadrature, params, solver, iguess, caches)
    end

    function IntegratorCGVI{DT,D}(equations::NamedTuple, basis::Basis{TT}, quadrature::QuadratureRule{TT,R}, Δt::TT;
                                  interpolation=HermiteExtrapolation{DT}) where {DT,TT,D,R}

        # get number of stages
        S = nbasis(basis)

        # create params
        params = ParametersCGVI{DT,D}(equations, Δt, basis, quadrature)

        # create cache dict
        caches = CacheDict(params)

        # create nonlinear solver
        solver = create_nonlinear_solver(DT, D*(S+1), params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorCGVI(basis, quadrature, params, solver, iguess, caches)
    end

    function IntegratorCGVI(equation::Union{IODE{DT}, LODE{DT}}, basis::Basis, quadrature::QuadratureRule, Δt; kwargs...) where {DT}
        IntegratorCGVI{DT, ndims(equation)}(get_functions(equation), basis, quadrature, Δt; kwargs...)
    end
end

@inline equation(integrator::IntegratorCGVI, i::Symbol) = integrator.params.equs[i]
@inline equations(integrator::IntegratorCGVI) = integrator.params.equs
@inline timestep(integrator::IntegratorCGVI) = integrator.params.Δt
@inline Base.ndims(::IntegratorCGVI{DT,TT,D}) where {DT,TT,D} = D

function IntegratorCache(int::IntegratorCGVI{DT,TT}) where {DT,TT}
    IntegratorCacheCGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}()
end

function initialize!(int::IntegratorCGVI, sol::AtomicSolutionPODE)
    sol.t̄ = sol.t - timestep(int)

    equation(int, :v̄)(sol.t, sol.q, sol.v)
    equation(int, :f̄)(sol.t, sol.q, sol.v, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̄, sol.q̄, sol.p̄, sol.v̄, sol.f̄)
end

function initial_guess!(int::IntegratorCGVI{DT,TT,D,S,R}, sol::AtomicSolutionPODE{DT,TT},
                        cache::IntegratorCacheCGVI{DT}=int.caches[DT]) where {DT,TT,D,S,R}
    int.solver.x .= 0

    for i in eachindex(int.basis)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                                sol.q, sol.p, sol.v, sol.f,
                                cache.q̃, cache.p̃,
                                cache.ṽ, cache.f̃,
                                grid(int.basis)[i], grid(int.basis)[i])

        for k in 1:D
            int.solver.x[D*(i-1)+k] = cache.q̃[k]
        end
    end

    evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                          sol.q, sol.p, sol.v, sol.f,
                          cache.q̃, cache.p̃,
                          one(TT), one(TT))

    for k in 1:D
        int.solver.x[D*S+k] = cache.p̃[k]
    end
end

"Compute stages of variational partitioned Runge-Kutta methods."
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R},
                caches::CacheDict) where {ST,DT,TT,D,S,R}
    @assert length(x) == length(b)

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache, params)

    # compute rhs b of nonlinear solver
    compute_rhs!(b, cache.X, cache.Q, cache.P, cache.F, cache.p̃, params)
end

function compute_stages!(x, cache::IntegratorCacheCGVI, params::ParametersCGVI)
    compute_stages!(x, cache.X, cache.Q, cache.V, cache.P, cache.F, cache.q̃, cache.p̃, params)
end

function compute_stages!(x, X, Q, V, P, F, q, p, params::ParametersCGVI{DT,TT,D,S,R}) where {DT,TT,D,S,R}

    # copy x to X
    for i in eachindex(X)
        for k in eachindex(X[i])
            X[i][k] = x[D*(i-1)+k]
        end
    end

    # copy x to p
    for k in eachindex(p)
        p[k] = x[D*S+k]
    end

    # compute Q
    compute_stages_q!(X, Q, q, params)

    # compute V
    compute_stages_v!(X, V, params)

    # compute P and F
    compute_stages_p!(Q, V, P, F, params)
end

function compute_stages_q!(X::Vector{Vector{ST}}, Q::Vector{Vector{ST}}, q::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    @assert R == length(Q)
    @assert S == length(X)

    local y::ST

    # compute Q
    for i in eachindex(Q)
        @assert D == length(Q[i])
        for k in eachindex(Q[i])
            y = 0
            for j in eachindex(X)
                @assert D == length(X[j])
                y += params.m[i,j] * X[j][k]
            end
            Q[i][k] = y
        end
    end

    # compute q
    for k in eachindex(q)
        y = 0
        for i in eachindex(X)
            y += params.r₁[i] * X[i][k]
        end
        q[k] = y
    end
end

function compute_stages_v!(X::Vector{Vector{ST}}, V::Vector{Vector{ST}}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    @assert R == length(V)
    @assert S == length(X)

    local y::ST

    # compute V
    for i in eachindex(V)
        @assert D == length(V[i])
        for k in eachindex(V[i])
            y = 0
            for j in eachindex(X)
                @assert D == length(X[j])
                y += params.a[i,j] * X[j][k]
            end
            V[i][k] = y / params.Δt
        end
    end
end

function compute_stages_p!(Q::Vector{Vector{ST}}, V::Vector{Vector{ST}},
                           P::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                           params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    @assert R == length(Q) == length(V) == length(P) == length(F)

    local tᵢ::TT

    # compute P=ϑ(Q,V) and F=f(Q,V)
    for i in eachindex(Q,V,P,F)
        @assert D == length(Q[i]) == length(V[i]) == length(P[i]) == length(F[i])
        tᵢ = params.t + params.Δt * params.c[i]
        params.equs[:ϑ](tᵢ, Q[i], V[i], P[i])
        params.equs[:f](tᵢ, Q[i], V[i], F[i])
    end
end

function compute_rhs!(b::Vector{ST}, X::Vector{Vector{ST}}, Q::Vector{Vector{ST}},
                                     P::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                     p::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local y::ST
    local z::ST

    # compute b = - [(P-AF)]
    for i in 1:S
        for k in 1:D
            z = 0
            for j in eachindex(P,F)
                z += params.b[j] * params.m[j,i] * F[j][k]
                z += params.b[j] * params.a[j,i] * P[j][k] / params.Δt
            end
            b[D*(i-1)+k] = z - (params.r₁[i] * p[k] - params.r₀[i] * params.p[k]) / params.Δt
        end
    end

    # compute b = - [(q-r₀Q)]
    for k in eachindex(params.q)
        y = 0
        for j in eachindex(X)
            y += params.r₀[j] * X[j][k]
        end
        b[D*S+k] = y - params.q[k]
    end
end

function update_solution!(sol::AtomicSolutionPODE, cache::IntegratorCacheCGVI)
    sol.q .= cache.q̃
    sol.p .= cache.p̃
end

function integrate_step!(int::IntegratorCGVI{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCacheCGVI{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # compute final update
    update_solution!(sol, cache)

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
using Sundials, Test

## f routine. Compute function f(t,y).

function f(t, y_nv, ydot_nv, user_data)
    y = convert(Vector, y_nv)
    ydot = convert(Vector, ydot_nv)
    ydot[1] = -0.04 * y[1] + 1.0e4 * y[2] * y[3]
    ydot[3] = 3.0e7 * y[2] * y[2]
    ydot[2] = -ydot[1] - ydot[3]
    return Sundials.ARK_SUCCESS
end

neq = 3

t0 = 0.0
t1 = 0.4
tmult = 10.0
nout = 12
y0 = [1.0, 0.0, 0.0]
reltol = 1e-4
abstol = 1e-11
userdata = nothing
h0 = 1e-4 * reltol

mem_ptr = Sundials.ARKStepCreate(C_NULL, f, t0, y0)
arkStep_mem = Sundials.Handle(mem_ptr)
Sundials.@checkflag Sundials.ARKStepSetInitStep(arkStep_mem, h0)
Sundials.@checkflag Sundials.ARKStepSetMaxErrTestFails(arkStep_mem, 20)
Sundials.@checkflag Sundials.ARKStepSetMaxNonlinIters(arkStep_mem, 8)
Sundials.@checkflag Sundials.ARKStepSetNonlinConvCoef(arkStep_mem, 1.e-7)
Sundials.@checkflag Sundials.ARKStepSetMaxNumSteps(arkStep_mem, 100000)
Sundials.@checkflag Sundials.ARKStepSetPredictorMethod(arkStep_mem, 1)

Sundials.@checkflag Sundials.ARKStepSStolerances(arkStep_mem, reltol, abstol)
A = Sundials.SUNDenseMatrix(neq, neq)
LS = Sundials.SUNLinSol_Dense(y0, A)
Sundials.@checkflag Sundials.ARKStepSetLinearSolver(arkStep_mem, LS, A)

iout = 0
tout = t1
t = [t0]

while iout < nout
    y = similar(y0)
    flag = Sundials.ARKStepEvolve(arkStep_mem, tout, y, t, Sundials.ARK_NORMAL)
    @test flag == 0
    println("T=", tout, ", Y=", y)
    global iout += 1
    global tout *= tmult
end

tmp1 = Ref(Clong(-1))
tmp2 = Ref(Clong(-1))
Sundials.@checkflag Sundials.ARKStepGetNumSteps(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumStepAttempts(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumRhsEvals(arkStep_mem, tmp1, tmp2);
Sundials.@checkflag Sundials.ARKStepGetNumLinSolvSetups(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumErrTestFails(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumNonlinSolvIters(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumNonlinSolvConvFails(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumJacEvals(arkStep_mem, tmp1);
Sundials.@checkflag Sundials.ARKStepGetNumLinRhsEvals(arkStep_mem, tmp2);


text:
module TDDecoupling

export solve_td_decoupling_distribution

using ..FlexPlan
const _FP = FlexPlan
import ..FlexPlan: _LOGGER

import PowerModels
const _PM = PowerModels

import JuMP
import Memento

include("base.jl")
include("data.jl")
include("model.jl")

end


text:
@testset "Test Reserves from Thermal Dispatch" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true)
    for p in [true, false]
        op_problem = OperationsProblem(template, c_sys5_uc; use_parameters = p)
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 648, 0, 120, 216, 72, false)
        symbols = [
            :Reserve1__VariableReserve_ReserveUp,
            :Reserve11__VariableReserve_ReserveUp,
            :Reserve2__VariableReserve_ReserveDown,
            :ORDC1__ReserveDemandCurve_ReserveUp,
        ]
        for sym in symbols
            for v in op_problem.internal.optimization_container.variables[sym]
                @test JuMP.has_lower_bound(v)
                @test JuMP.lower_bound(v) == 0.0
            end
        end
    end
end

@testset "Test Reserves from Thermal Standard UC" begin
    template = get_thermal_standard_uc_template()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true)

    for p in [true, false]
        op_problem = OperationsProblem(template, c_sys5_uc; use_parameters = p)
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 1008, 0, 480, 216, 192, true)
    end
end

@testset "Test Upwards Reserves from Renewable Dispatch" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )

    c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re"; add_reserves = true)
    for p in [true, false]
        op_problem = OperationsProblem(template, c_sys5_re; use_parameters = p)
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 360, 0, 72, 48, 72, false)
    end
end

@testset "Test Reserves from Storage" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    set_device_model!(template, GenericBattery, BookKeeping)
    set_device_model!(template, RenewableDispatch, FixedOutput)
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )

    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat"; add_reserves = true)
    for p in [true, false]
        op_problem = OperationsProblem(template, c_sys5_bat; use_parameters = p)
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 408, 0, 192, 264, 96, false)
    end
end

@testset "Test Reserves from Hydro" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )

    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd"; add_reserves = true)
    for p in [true, false]
        op_problem = OperationsProblem(template, c_sys5_hyd; use_parameters = p)
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 240, 0, 24, 96, 72, false)
    end
end

@testset "Test Reserves from with slack variables" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true)
    for p in [true, false]
        op_problem = OperationsProblem(
            template,
            c_sys5_uc;
            use_parameters = p,
            services_slack_variables = true,
            balance_slack_variables = true,
        )
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 504, 0, 120, 192, 24, false)
    end
end

@testset "Test AGC" begin
    c_sys5_reg = PSB.build_system(PSITestSystems, "c_sys5_reg")
    @test_throws ArgumentError template_agc_reserve_deployment(; dummy_arg = 0.0)

    template_agc = template_agc_reserve_deployment()
    agc_problem = OperationsProblem(AGCReserveDeployment, template_agc, c_sys5_reg)
    @test build!(agc_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    # These values might change as the AGC model is refined
    moi_tests(agc_problem, false, 720, 0, 480, 0, 384, false)
end

@testset "Test GroupReserve from Thermal Dispatch" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    set_service_model!(
        template,
        ServiceModel(StaticReserveGroup{ReserveDown}, GroupReserve),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true)
    services = get_components(Service, c_sys5_uc)
    contributing_services = Vector{Service}()
    for service in services
        if !(typeof(service) <: PSY.ReserveDemandCurve)
            push!(contributing_services, service)
        end
    end
    groupservice = StaticReserveGroup{ReserveDown}(;
        name = "init",
        available = true,
        requirement = 0.0,
        ext = Dict{String, Any}(),
    )
    add_service!(c_sys5_uc, groupservice, contributing_services)

    for p in [true, false]
        op_problem = OperationsProblem(template, c_sys5_uc; use_parameters = p)
        @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        moi_tests(op_problem, p, 648, 0, 120, 240, 72, false)
    end
end

# TODO: Test is broken
#=
@testset "Test GroupReserve Errors" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    set_service_model!(
        template,
        ServiceModel(StaticReserveGroup{ReserveDown}, GroupReserve),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true)
    services = get_components(Service, c_sys5_uc)
    contributing_services = Vector{Service}()
    for service in services
        if !(typeof(service) <: PSY.ReserveDemandCurve)
            push!(contributing_services, service)
        end
    end
    groupservice = StaticReserveGroup{ReserveDown}(;
        name = "init",
        available = true,
        requirement = 0.0,
        ext = Dict{String, Any}(),
    )
    add_service!(c_sys5_uc, groupservice, contributing_services)

    off_service = VariableReserve{ReserveUp}("Reserveoff", true, 0.6, 10)
    push!(groupservice.contributing_services, off_service)

    op_problem = OperationsProblem(template, c_sys5_uc; use_parameters = false)
    @test_logs(
        (:error, r"is not stored"),
        match_mode = :any,
        @test_throws InfrastructureSystems.InvalidValue build!(op_problem; output_dir = mktempdir(cleanup = true))
    )
end
=#

@testset "Test StaticReserve" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(template, ServiceModel(StaticReserve{ReserveUp}, RangeReserve))

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    static_reserve = StaticReserve{ReserveUp}("Reserve3", true, 30, 100)
    add_service!(c_sys5_uc, static_reserve, get_components(ThermalGen, c_sys5_uc))
    op_problem = OperationsProblem(template, c_sys5_uc)
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    @test typeof(op_problem) <: OperationsProblem
end


text:
using project1
using Test

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


text:
module ReinforcementLearningZoo

using ReinforcementLearningBase
using ReinforcementLearningCore

include("algorithms/algorithms.jl")

end # module


text:
# Computing deviation in a variety of ways

# squared L2 distance
function sqL2dist{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T})
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Input dimension mismatch"))
    r = 0.0
    for i = 1:n
        @inbounds r += abs2(a[i] - b[i])
    end
    return r
end

# L2 distance
L2dist{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T}) = sqrt(sqL2dist(a, b))

# L1 distance
function L1dist{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T})
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Input dimension mismatch"))
    r = 0.0
    for i = 1:n
        @inbounds r += abs(a[i] - b[i])
    end
    return r
end

# Linf distance
function Linfdist{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T})
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Input dimension mismatch"))
    r = 0.0
    for i = 1:n
        @inbounds v = abs(a[i] - b[i])
        if r < v
            r = v
        end
    end
    return r
end

# Generalized KL-divergence
function gkldiv{T<:FloatingPoint}(a::AbstractArray{T}, b::AbstractArray{T})
    n = length(a)
    r = 0.0
    for i = 1:n
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        if ai > 0
            r += (ai * log(ai / bi) - ai + bi)
        else
            r += bi
        end
    end
    return r::Float64
end

# MeanAD: mean absolute deviation
meanad{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T}) = L1dist(a, b) / length(a)

# MaxAD: maximum absolute deviation
maxad{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T}) = Linfdist(a, b)

# MSD: mean squared deviation
msd{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T}) = sqL2dist(a, b) / length(a)

# RMSD: root mean squared deviation
function rmsd{T<:Number}(a::AbstractArray{T}, b::AbstractArray{T}; normalize::Bool=false)
	v = sqrt(msd(a, b))
	if normalize
		amin, amax = extrema(a)
		v /= (amax - amin)
	end
	return v
end

# PSNR: peak signal-to-noise ratio
function psnr{T<:Real}(a::AbstractArray{T}, b::AbstractArray{T}, maxv::Real)
    20. * log10(maxv) - 10. * log10(msd(a, b))
end


text:
using LogDensityTestSuite, Test, Statistics, LinearAlgebra, Distributions, StatsFuns
import ForwardDiff
using LogDensityProblems: capabilities, dimension, logdensity, logdensity_and_gradient,
    LogDensityOrder
using LogDensityTestSuite: hypercube_dimension, _find_x_norm, weight, weight_and_gradient

"Test gradient with automatic differentiation."
function test_gradient(ℓ, x; atol = √eps())
    l, g = logdensity_and_gradient(ℓ, x)
    l2 = logdensity(ℓ, x)
    g2 = ForwardDiff.gradient(x -> logdensity(ℓ, x), x)
    @test l ≈ l2 atol = atol
    @test g ≈ g2 atol = atol
end

####
#### primitives
####

@testset "standard multivariate normal" begin
    K, N = 5, 1000
    ℓ = StandardMultivariateNormal(K)
    @test dimension(ℓ) == hypercube_dimension(ℓ) == K
    @test capabilities(ℓ) == LogDensityOrder(1)
    Z = samples(ℓ, N)
    @test size(Z) == (K, N)
    d = MvNormal(zeros(K), Diagonal(ones(K)))
    for x in eachcol(Z)
        test_gradient(ℓ, x)
        @test logdensity(ℓ, x) ≈ logpdf(d, x)
    end
    @test mean(Z; dims = 2) ≈ zeros(K) atol = 0.01
    @test std(Z; dims = 2) ≈ ones(K) atol = 0.02
end

####
#### transformations
####

@testset "multivariate normal using transform" begin
    K = 4
    Q = qr(reshape(range(0.1; step = 0.05, length = K * K), K, K)).Q
    D = Diagonal(range(1; step = .1, length = K))
    μ = collect(range(0.04; step = 0.2, length = K))

    function test_mvnormal(μ, A, Σ)
        K = size(A, 1)
        ℓ = shift(μ, linear(A, StandardMultivariateNormal(K)))
        d = MvNormal(μ, Σ)
        @test dimension(ℓ) == hypercube_dimension(ℓ) == K
        @test capabilities(ℓ) == LogDensityOrder(1)
        Z = samples(ℓ, 1000)
        for x in eachcol(Z)
            f, ∇f = logdensity_and_gradient(ℓ, x)
            @test f ≈ logpdf(d, x)
            @test ∇f ≈ gradlogpdf(d, x)
        end
        @test mean(Z; dims = 2) ≈ μ atol = 0.02
        @test norm(cov(Z; dims = 2) .- Σ, Inf) ≤ 0.07
    end

    @testset "MvNormal Diagonal" begin
        A = Diagonal(2 .* ones(K))
        Σ = A * A
        test_mvnormal(μ, A, Σ)
    end

    @testset "MvNormal triangular" begin
        Σ = Symmetric(Q * D * Q')
        A = cholesky(Σ, Val(false)).L
        test_mvnormal(μ, A, Σ)
    end

    @testset "MvNormal full" begin
        Σ = Symmetric(Q * D * Q')
        A = Q * Diagonal(.√diag(D))
        test_mvnormal(μ, A, Σ)
    end
end

@testset "elongate building blocks" begin
    for _ in 1:1000
        y = abs(randn())
        k = abs(randn() * 3)
        x = _find_x_norm(y, k)
        @test y ≈ x * (1 + abs2(x))^k
    end
end

@testset "elongate" begin
    K, N = 5, 1000
    ℓ = elongate(0.5, StandardMultivariateNormal(5))
    @test dimension(ℓ) == hypercube_dimension(ℓ) == K
    @test capabilities(ℓ) == LogDensityOrder(1)
    Z = samples(ℓ, N)
    @test vec(mean(Z; dims = 2)) ≈ zeros(K) norm = x -> norm(x, Inf) atol = 0.02
    for x in eachcol(Z)
        test_gradient(ℓ, x)
    end
end

####
#### mixtures
####

@testset "scalar mixture" begin
    K, N = 5, 1000
    ℓ1 = StandardMultivariateNormal(K)
    μ2 = fill(1.7, K)
    ℓ2 = shift(μ2, linear(Diagonal(fill(0.01, K)), StandardMultivariateNormal(K)))
    for α in [0, √eps(), 0.7, 1-√eps(), 1] # extreme, near-extreme, interior α
        ℓ = mix(α, ℓ1, ℓ2)
        @test dimension(ℓ) == K
        @test hypercube_dimension(ℓ) == K + 1
        @test capabilities(ℓ) == LogDensityOrder(1)
        Z = samples(ℓ, N)
        @test size(Z) == (K, N)

        # test at sample values
        for x in eachcol(Z)
            @test logdensity(ℓ, x) ≈
                logaddexp(log(α) + logdensity(ℓ1, x), log1p(-α) + logdensity(ℓ2, x))
            test_gradient(ℓ, x)
        end

        # test at random extreme values (numerical stability)
        for _ in 1:1000
            x = normalize!(randn(K), 2) .* 100
            @test logdensity(ℓ, x) ≈
                logaddexp(log(α) + logdensity(ℓ1, x), log1p(-α) + logdensity(ℓ2, x))
            test_gradient(ℓ, x)
        end

        # compare means
        @test mean(Z; dims = 2) ≈ (1 - α) .* μ2 atol = 0.02
    end

    @test_throws ArgumentError mix(0.5, ℓ1, StandardMultivariateNormal(K + 1))
end

@testset "directional mixture" begin
    K, N = 5, 1000
    ℓ1 = StandardMultivariateNormal(K)
    μ2 = fill(1.7, K)
    ℓ2 = shift(μ2, linear(Diagonal(fill(0.01, K)), StandardMultivariateNormal(K)))
    α = directional_weight(ones(K))
    ℓ = mix(α, ℓ1, ℓ2)
    @test dimension(ℓ) == K
    @test hypercube_dimension(ℓ) == K + 1
    @test capabilities(ℓ) == LogDensityOrder(1)
    Z = samples(ℓ, N)
    @test size(Z) == (K, N)
    # test at sample values
    for x in eachcol(Z)
        αx, ∇αx = weight_and_gradient(α, x)
        @test ∇αx ≈ ForwardDiff.gradient(x -> weight(α, x), x)
        test_gradient(ℓ, x)
    end
    @test_throws ArgumentError directional_weight(zeros(5))
end

####
#### diagnostics
####

@testset "diagnostics" begin
    x = range(0, 1; length = 100001)
    q = quantile_boundaries(x, 10)

    # test p-values and printing
    bc = bin_counts(q, rand(1000))
    bc.bin_counts[1] = 1        # to test printing extremes
    ps = two_sided_pvalues(bc)
    @test all(0 .≤ ps .≤ 1)
    @test print_ascii_plot(String, bc) isa String # very rudimentary
    @info "this is what a printed ascii plot looks like"
    show(stdout, bc)

    # validate p-values
    q̄ = (1:3)./4
    ps = reduce(hcat, [two_sided_pvalues(bin_counts(q, rand(1000))) for _ in 1:1000])
    for i in axes(ps, 1)
        @test quantile(ps[i, :], q̄) ≈ q̄ atol = 0.1 norm = x -> norm(x, Inf)
    end
end


text:
function update_eta(m::MixSkewT, s::S) where S
  anew = copy(collect(m.eta.alpha))
  foreach(lam -> anew[lam] +=1, s.lambda)
  return rand(Dirichlet(anew))
end

function update_lambda(m::MixSkewT, s::S) where S
  logeta = log.(s.eta)

  lambda = [let
              loc = s.mu + s.psi * s.zeta[n]
              scale = sqrt.(s.omega) / sqrt(s.v[n])
              logmix = normlogpdf.(loc, scale, m.y[n]) + logeta
              all(s.nu .< 10000) && (logmix += gammalogpdf.(s.nu/2, 2 ./ s.nu, s.v[n]))
              MCMC.wsample_logprob(logmix)
            end for n in eachindex(m.y)]

  return lambda
end

function update_v(m::MixSkewT, s::S) where S
  nu = s.nu[s.lambda]
  omega = s.omega[s.lambda]
  psi = s.psi[s.lambda]
  mu = s.mu[s.lambda]

  shape = nu/2 .+ 1
  rate = (nu + s.zeta.^2 + ((m.y - mu - psi .* s.zeta) .^ 2) ./ omega) / 2

  return @.rand(Gamma(shape, 1 / rate))
end

function update_zeta(m::MixSkewT, s::S) where S
  mu = s.mu[s.lambda]
  omega = s.omega[s.lambda]
  psi = s.psi[s.lambda]

  vnew = @. 1 / (s.v + (psi ^ 2) * s.v / omega)
  mnew = @. vnew * s.v * psi * (m.y - mu) / omega

  return @.rand(truncated(Normal(mnew, sqrt(vnew)), 0, Inf))
end

function update_mu(m::MixSkewT, s::S) where S
  m_mu, s_mu = params(m.mu)
  vkernels = zero.(s.mu)
  mkernels = zero.(s.mu)

  for i in (:C, :T)
    for n in eachindex(s.lambda)
      k = s.lambda[n]
      vkernels[k] += s.v[n]
      mkernels[k] += (m.y[n] - s.psi[k] * s.zeta[n]) * s.v[n]
    end
  end

  vnew = 1 ./ (s_mu^-2 .+ vkernels ./ s.omega)
  mnew = vnew .* (m_mu/(s_mu^2) .+ mkernels ./ s.omega)

  return randn(m.K) .* sqrt.(vnew) + mnew
end

function update_tau(m::MixSkewT, s::S) where S
  a, b = params(m.tau)
  new_shape = a + m.K * m.a_omega
  new_rate = b + sum(1 ./ s.omega)
  return rand(Gamma(new_shape, 1 / new_rate))
end

function update_sigma(m::MixSkewT, s::S) where S
  return @. scalefromaltskewt(sqrt(s.omega), s.psi)
end

function update_phi(m::MixSkewT, s::S) where S
  return @. skewfromaltskewt(sqrt(s.omega), s.psi)
end

function update_sigma_phi(m::MixSkewT, s::S) where S
  return (sigma=update_sigma(m, s), phi=update_phi(m, s))
end

function update_omega(m::MixSkewT, s::S) where S
  a, b = m.a_omega, s.tau
  akernel = zero.(s.omega)
  bkernel = zero.(s.omega)

  for n in eachindex(m.y)
    k = s.lambda[n]
    akernel[k] += 1
    bkernel[k] += s.v[n] * (m.y[n] - s.mu[k] - s.psi[k] * s.zeta[n])^2
  end

  anew = a .+ akernel / 2
  bnew = b .+ bkernel / 2

  return @.rand(InverseGamma(anew, bnew))
end

function logprob_nu(m::MixSkewT, s::S, nu::AbstractVector{<:Real}) where S
  logprior = sum(logpdf.(m.nu, nu))
  loglike = let
    nu = nu[s.lambda]
    sum(gammalogpdf.(nu / 2, 2 ./ nu, s.v))
  end
  return logprior + loglike
end

function update_psi(m::MixSkewT, s::S) where S
  m_psi, s_psi = params(m.psi)
  vkernel = zero.(s.psi)
  mkernel = zero.(s.psi)

  # Update kernels
  for n in eachindex(s.v)
    k = s.lambda[n]
    vkernel[k] += s.zeta[n]^2 * s.v[n] / s.omega[k]
    mkernel[k] += s.zeta[n] * (m.y[n] - s.mu[k]) * s.v[n] / s.omega[k]
  end

  vnew = 1 ./ (s_psi^-2 .+ vkernel)
  mnew = vnew .* (m_psi/s_psi^2 .+ mkernel)

  return randn(m.K) .* sqrt.(vnew) + mnew
end


text:
module clustering

using LinearAlgebra

export connectivity, calculate_hierarchy_matrix

function order_unique( values )
    svals = Set(reshape( values , (length(values)) ))
    rvals = []; for v in svals append!(rvals,v) end
    rvals = sort(rvals)
    return(rvals)
end

function calculate_hierarchy_matrix( distance_matrix::Array{Float64,2} , bVerbose::Bool = false )
    """ This is the saiga/pelican/panda you are looking for:
        (linkage free hierarchical clustering)
        https://github.com/richardtjornhammar/impetuous/blob/master/src/impetuous/hierarchical.py
	around line 118
    """
    NP = size(distance_matrix[:,1])[1]
    uco_v = convert( Array{Float64}, order_unique( distance_matrix ) )
    results = Dict("level"=>[] , "distance"=>[] , "hierarchy_matrix"=>Array{Int64} )
    hmat = []
    for icut in 1:length(uco_v)
        cutoff = uco_v[icut]
	outp = connectivity( distance_matrix, cutoff, bVerbose )
	append!(results["level"]           , icut   )
	append!(results["distance"]        , cutoff )
	cout = reshape( convert( Array{Int64} , outp[:,1] ) ,(1,NP))
	if icut==1
	    hmat = cout
	else
	    hmat = vcat(hmat,cout)
	end
    end
    results["hierarchy_matrix"] = hmat
    results["level"] = convert(Array{Int64,1},results["level"])
    results["distance"] = convert(Array{Float64,1},results["distance"])
    return(results)
end

function connectivity( B::Array{Float64,2} , val::Float64 , bVerbose::Bool=false )
    """
This is a cutoff based clustering algorithm. The intended use is to supply a distance matrix and a cutoff value (then becomes symmetric positive definite).  For a small distance cutoff, you should see all the parts of the system and for a large distance cutoff, you should see the entire system. It has been employed for statistical analysis work as well as the original application where it was employed to segment molecular systems.
        #
	# JULIA ADAPTATION OF MY C++ CODE THAT CAN BE FOUND IN
	# https://github.com/richardtjornhammar/RichTools/blob/master/src/cluster.cc
	# AROUND LINE 2277
	# FOR A DESCRIPTION READ PAGE 30 (16 INTERNAL NUMBERING) of:
	# https://kth.diva-portal.org/smash/get/diva2:748464/FULLTEXT01.pdf
	#
    """
    nr_sq,mr_sq = size(B)
    if nr_sq != mr_sq
        println( "::ERROR::" )
        return(-1)
    end
    N = mr_sq
    res,nvisi,s,NN,ndx,C = [],[],[],[],[],0
    #res = append!(res,[0])
    for i in 1:N
        append!(nvisi,[i])
        append!(res,[0]); append!(res,[0])
        append!(ndx,[i])
    end
    res   = convert(Array{Int64},res  )
    ndx   = convert(Array{Int64},ndx  )
    nvisi = convert(Array{Int64},nvisi)
    if bVerbose
        println(res," ",ndx," ",nvisi)
        println(nvisi[end])
    end
    while length(ndx)>0
        i = pop!(ndx)
        NN = []
        if nvisi[i]>0
            C-=1
            for j in 1:N
                if B[j,i]<=val
                    append!(NN,[j])
                end
            end
            while length(NN)>0
                k = pop!(NN)
                nvisi[k] = C
                for j in 1:N
                    if B[j,k]<=val
                        for q in 1:N
                            if nvisi[q]==j
                                append!(NN,[q])
                            end
                        end
                    end
                end
            end
        end
    end
    if bVerbose
        println("INFO ", C,"\ncluster data:")
        Nc = []
        for i in 1:-1*C
            append!(Nc,[0])
        end
    end
    for q in 1:N
        res[q*2]=q
        res[q*2-1]=nvisi[q]-C+1
        if bVerbose
	    Nc[res[q*2-1]]+=1
            println(" ",res[q*2-1],"\t",res[2*q])
        end
    end
    if bVerbose
        for i in 1:-1*C
            println("Cluster ",i," has ", Nc[i]," elements")
        end
    end
    ret = transpose(reshape(res,(2,N)))
    if bVerbose
        println( ret, size(ret) )
    end
    return ( ret )
end

end

max_stars_repo_name: richardtjornhammar/ingenuous

max_issues_repo_name: richardtjornhammar/ingenuous

max_forks_repo_name: richardtjornhammar/ingenuous


text:
using NativeFileDialog
using Test

import NativeFileDialog: getindices, getpaths, NFDPathSet

@testset "NativeFileDialog.jl" begin
    @test length(NFDPathSet()) == 0
    @test getindices(NFDPathSet()) == UInt[]
    @test getpaths(NFDPathSet()) == String[]

    paths = "a\0b"
    GC.@preserve paths begin
        indices = [0, 2]
        count = 2
        pathset = NFDPathSet(
                    pointer(paths),
                    convert(Ptr{Csize_t}, pointer(indices)),
                    convert(Csize_t, count))

        @test length(pathset) == 2
        @test getindices(pathset) == indices
        @test getpaths(pathset) == ["a", "b"]
    end
end


text:
"""
    AbstractCallback

Callback to be performed after every time step when solving the BTPDE. The callback is
initialized before time stepping using the `initialize!` function, updated after every step
using `update!`, and finalized after time stepping using `finalize!`.
"""
abstract type AbstractCallback end

"""
    Printer(; nupdate = 1)

Print time stepping information to the console.
"""
Base.@kwdef mutable struct Printer <: AbstractCallback
    verbosity::Int = 2
    nupdate::Int = 1
    n::Int = 1
end

"""
    VTKWriter(; nupdate = 1, dir = "output", filename = "solution")

Write magnetization field to a VTK file after every `nupdate` timestep. The files are stored in a
ParaView data collection file (PVD). The magnetization time series may be visualized in
ParaView by opening the file `"\$dir/\$filename.pvd"`. The compartments are labeled.
"""
Base.@kwdef mutable struct VTKWriter <: AbstractCallback
    nupdate::Int = 1
    n::Int = 1
    ifile::Int = 1
    dir::String = "output/time_stepping"
    filename::String = "solution"
    pvd::WriteVTK.CollectionFile = paraview_collection("")
end

"""
    Plotter(; nupdate = 1)

Plot the evolution of the BTPDE during time stepping. This requires loading the `GLMakie`
plotting backend (`]add GLMakie; using GLMakie`). The plot is updated every `nupdate` time
step. The resulting figure contains a plot of the time profile, total signal attenuation,
and magnetization field (complex magnitude and phase shift).
"""
Base.@kwdef mutable struct Plotter{T} <: AbstractCallback
    nupdate::Int = 1
    n::Int = 1
    t::Node{Vector{T}} = Node(T[])
    f::Node{Vector{T}} = Node(T[])
    g⃗::Node{Vector{Vec3f}} = Node([Vec3f(0, 0, 0)])
    g⃗_hist::Node{Vector{Vec3f}} = Node([Vec3f(0, 0, 0)])
    ξ::Node{Vector{Complex{T}}} = Node(Complex{T}[])
    magnitude::Matrix{Node{Vector{T}}} = [Node(Vector{T}());;]
    phase::Matrix{Node{Vector{T}}} = [Node(Vector{T}());;]
    attenuation::Node{Vector{T}} = Node(T[])
    S₀::T = one(T)
    fig::Figure = Figure()
end


text:
#=
You have an N by N board. Write a function that, given N, returns the number of possible arrangements of the board where N queens can be placed on the board without threatening each other, i.e. no two queens share the same row, column, or diagonal.
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


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

using Documenter, FEMQuad

makedocs(modules=[FEMQuad],
         format = Documenter.HTML(),
         sitename = "FEMQuad.jl",
         pages = ["index.md", "api.md"])


text:
function crawl(::Val{:$}, x)
    return :(interpolate($(esc(x))))
end

"""
Determine what JSNode should be inserted when a value is interpolated.

This function is called with the value that is to be interpolated and not the
expression itself that is being interpolated (i.e., after the crawl phase).
"""
function interpolate(x)
    return JSTerminal(_tojsstring(tojs(x)))
end

"""
    _tojsstring(x)

Coerse a value to a JSString representation.
This is required since `tojs` only returns something that is either a `JSString`
or something that is `JSON.lower`-able, but we just want a `JSString`.
"""
_tojsstring(s::JSString) = s
_tojsstring(s) = JSString(JSON.json(s))