text:
"""
Reference:
 Shengyang Sun, Guodong Zhang, Chaoqi Wang, Wenyuan Zeng, Jiaman Li, Roger Grosse (2018)
 Differentiable Compositional Kernel Learning for Gaussian Processes
"""

using Base: tail
using Flux: glorot_uniform
import Flux: functor

"""
Neural Kernel Network
modify https://github.com/FluxML/Flux.jl/tree/master/src/layers/basic.jl Chain
"""
struct NeuralKernelNetwork{T<:Tuple} <: AbstractKernel
	layers::T
	NeuralKernelNetwork{T}(ls...) where {T} = new{T}(ls)
end
NeuralKernelNetwork(ls...) = NeuralKernelNetwork{typeof(ls)}(ls...)

functor(nkn::NeuralKernelNetwork) = nkn.layers, ls->NeuralKernelNetwork(ls...)

applynkn(::Tuple{}, x) = x
applynkn(fs::Tuple, x) = applynkn(tail(fs), first(fs)(x))
applynkn(fs::Tuple, x, xo) = applynkn(tail(fs), first(fs)(x, xo))

function (nkn::NeuralKernelNetwork)(x; λ=1e-6)
	N = size(x, 2)
	x1 = applynkn(nkn.layers, x)
	K = reshape(x1, N, N)
	K + Diagonal(λ*ones(N))
end

function (nkn::NeuralKernelNetwork)(x, xo)
	M = size(x, 2)
	N = size(xo, 2)
	x1 = applynkn(nkn.layers, x, xo)
	K = reshape(x1, M, N)
	K
end

function Base.show(io::IO, kernel::NeuralKernelNetwork)
	print(io, "NeuralKernelNetwork(")
	join(io, kernel.layers, ", ")
	print(io, ")")
end

"""
Primitive layer, constituted by basic kernels
"""
struct Primitive{T<:Tuple{Vararg{<:AbstractKernel}}}
  kernels::T
  Primitive(ks...) = new{typeof(ks)}(ks)
end
functor(p::Primitive) = p.kernels, ks -> Primitive(ks...)

function (p::Primitive)(x)
	Ks = [reshape(Ker(x), 1, :) for Ker in p.kernels]
	vcat(Ks...)
end
function (p::Primitive)(x, xo)
	Ks = [reshape(Ker(x, xo), 1, :) for Ker in p.kernels]
	vcat(Ks...)
end

function Base.show(io::IO, layer::Primitive)
	print(io, "Primitive(")
	join(io, layer.kernels, ", ")
	print(io, ")")
end

"""
Linear layer: weigths & bias should be positive
"""
positive(x) = @. log(1.0 + exp(x)) + 1e-6  # plus 1e-6 to avoid this value to be 0.0
positive_back(x, dy) = @. dy*(1.0/(1.0+exp(-x)))

struct Linear{WT, BT}
	W::WT
	b::BT
end
@functor Linear
function (lin::Linear)(x)
	pos_W = positive(lin.W)
	pos_b = positive(lin.b)
	pos_W*x .+ pos_b
end
function Linear(in::Integer, out::Integer; init_W=glorot_uniform, init_b=zeros)
	Linear(init_W(out, in), init_b(out))
end

function Base.show(io::IO, layer::Linear)
	print(io, "Linear(", size(layer.W, 2), ",", size(layer.W, 1), ")")
end

"""
Product layer
"""
function Product(x; step=2)
	m, n = size(x)
	m%step == 0 || error("the first dimension of inputs must be multiple of step")
	new_x = reshape(x, step, m÷step, n)
	.*([new_x[i, :, :] for i in 1:step]...)
end

"""
common used all product and all sum
"""
# not use `prod` here, since it'll results in unknown overflow in backpropogation
allProduct(x) = .*([x[i,:] for i in 1:size(x, 1)]...)

# similar reason, not use `sum`
function allSum(x)
	m = size(x, 1)
	v = ones(eltype(x), m)
	v'*x
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
# # Conjugate Gradient
# In this tutorial we describe the basics of using the conjugate gradient iterative solvers
# At the end you should be able to
# 1. Use Conjugate Gradient to solve a linear system
# 2. Know when to not use it
# 3. Contruct a column-wise linear solver with Conjugate Gradient

# ## What is it?
# Conjugate Gradient is an iterative method for solving special kinds of linear systems:
# ```math
#  Ax = b
# ```
# via iterative methods.
# !!! warning
#     The linear operator need to be symmetric positive definite and the preconditioner must be symmetric.
# See the [wikipedia](https://en.wikipedia.org/wiki/Conjugate_gradient_method) for more details.

# ## Basic Example
# First we must load a few things
using ClimateMachine,
    ClimateMachine.LinearSolvers, ClimateMachine.ConjugateGradientSolver
using LinearAlgebra, Random

# Next we define a 3x3 symmetric positive definite linear system. (In the ClimateMachine code a symmetric positive definite system could arise from treating diffusion implicitly.)
A = [
    2.0 -1.0 0.0
    -1.0 2.0 -1.0
    0.0 -1.0 2.0
];
# We define the matrix `A` here as a global variable for convenience later.

# We can see that it is symmetric. We can check that it is positive definite by checking the spectrum
eigvals(A)

# The linear operators that are passed into the abstract iterative solvers need to be defined as functions that act on vectors. Let us do that with our matrix. We are using function closures for type stability.
function closure_linear_operator!(A)
    function linear_operator!(x, y)
        mul!(x, A, y)
    end
    return linear_operator!
end;

# We now define our linear operator using the function closure

linear_operator! = closure_linear_operator!(A)

# We now define our `b` in the linear system
b = ones(typeof(1.0), 3);

# The exact solution to the system `Ax = b` is
x_exact = [1.5, 2.0, 1.5];

# Now we can set up the ConjugateGradient struct
linearsolver = ConjugateGradient(b);
# and an initial guess for the iterative solver.
x = ones(typeof(1.0), 3);
# To solve the linear system we just need to pass to the linearsolve! function
iters = linearsolve!(linear_operator!, linearsolver, x, b)
# The variable `x` gets overwitten during the linear solve
# The norm of the error is
norm(x - x_exact) / norm(x_exact)
# The relative norm of the residual is
norm(A * x - b) / norm(b)
# The number of iterations is
iters
# Conjugate Gradient is guaranteed to converge in 3 iterations with perfect arithmetic in this case.

# ## Non-Example

# Conjugate Gradient is not guaranteed to converge with nonsymmetric matrices. Consider
A = [
    2.0 -1.0 0.0
    0.0 2.0 -1.0
    0.0 0.0 2.0
];
# We define the matrix `A` here as a global variable for convenience later.

# We can see that it is not symmetric, but it does have all positive eigenvalues
eigvals(A)

# The linear operators that are passed into the abstract iterative solvers need to be defined as functions that act on vectors. Let us do that with our matrix. We are using function closures for type stability.
function closure_linear_operator!(A)
    function linear_operator!(x, y)
        mul!(x, A, y)
    end
    return linear_operator!
end;

# We define the linear operator using our function closure

linear_operator! = closure_linear_operator!(A)

# We now define our `b` in the linear system
b = ones(typeof(1.0), 3);

# The exact solution to the system `Ax = b` is
x_exact = [0.875, 0.75, 0.5];

# Now we can set up the ConjugateGradient struct
linearsolver = ConjugateGradient(b, max_iter = 100);
# We also passed in the keyword argument "max_iter" for the maximum number of iterations of the iterative solver. By default it is assumed to be the size of the vector.
# As before we need to define an initial guess
x = ones(typeof(1.0), 3);
# To (not) solve the linear system we just need to pass to the linearsolve! function
iters = linearsolve!(linear_operator!, linearsolver, x, b)
# The variable `x` gets overwitten during the linear solve
# The norm of the error is
norm(x - x_exact) / norm(x_exact)
# The relative norm of the residual is
norm(A * x - b) / norm(b)
# The number of iterations is
iters
# Conjugate Gradient is guaranteed to converge in 3 iterations with perfect arithmetic for a symmetric positive definite matrix. Here we see that the matrix is not symmetric and it didn't converge even after 100 iterations.

# ## More Complex Example
# Here we show how to construct a column-wise iterative solver similar to what is is in the ClimateMachine code. The following is not for the faint of heart.
# We must first define a linear operator that acts like one in the ClimateMachine
function closure_linear_operator!(A, tup)
    function linear_operator!(y, x)
        alias_x = reshape(x, tup)
        alias_y = reshape(y, tup)
        for i6 in 1:tup[6]
            for i4 in 1:tup[4]
                for i2 in 1:tup[2]
                    for i1 in 1:tup[1]
                        tmp = alias_x[i1, i2, :, i4, :, i6][:]
                        tmp2 = A[i1, i2, i4, i6] * tmp
                        alias_y[i1, i2, :, i4, :, i6] .=
                            reshape(tmp2, (tup[3], tup[5]))
                    end
                end
            end
        end
    end
end;

# Now that we have this function, we can define a linear system that we will solve columnwise
# First we define the structure of our array as ```tup``` in a manner that is similar to a stacked brick topology
tup = (3, 4, 7, 2, 20, 2);
# where
# 1. tup[1] is the number of Gauss–Lobatto points in the x-direction
# 2. tup[2] is the number of Gauss–Lobatto points in the y-direction
# 3. tup[3] is the number of Gauss–Lobatto points in the z-direction
# 4. tup[4] is the number of states
# 6. tup[5] is the number of elements in the vertical direction
# 7. tup[6] is the number of elements in the other directions

# Now we define our linear operator as a random matrix.
Random.seed!(1235);
B = [
    randn(tup[3] * tup[5], tup[3] * tup[5])
    for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
];
columnwise_A = [
    B[i1, i2, i4, i6] * B[i1, i2, i4, i6]' + 10I
    for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
];
columnwise_inv_A = [
    inv(columnwise_A[i1, i2, i4, i6])
    for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
];
columnwise_linear_operator! = closure_linear_operator!(columnwise_A, tup);
columnwise_inverse_linear_operator! =
    closure_linear_operator!(columnwise_inv_A, tup);

# We define our `x` and `b` with matrix structures similar to an MPIStateArray
mpi_tup = (tup[1] * tup[2] * tup[3], tup[4], tup[5] * tup[6]);
b = randn(mpi_tup);
x = randn(mpi_tup);

# Now we solve the linear system columnwise
linearsolver = ConjugateGradient(
    x,
    max_iter = tup[3] * tup[5],
    dims = (3, 5),
    reshape_tuple = tup,
);
# The keyword arguments dims is the reduction dimension for the linear solver. In this case dims = (3,5) are the ones associated with a column. The reshape_tuple argument is to convert the shapes of the array `x` in the a form that is more easily usable for reductions in the linear solver

# Now we can solve it
iters = linearsolve!(columnwise_linear_operator!, linearsolver, x, b);
x_exact = copy(x);
columnwise_inverse_linear_operator!(x_exact, b);
# The norm of the error is
norm(x - x_exact) / norm(x_exact)
# The number of iterations is
iters
# The algorithm converges within `tup[3]*tup[5] = 140` iterations

# ## Tips
# 1. The convergence criteria should be changed, machine precision is too small and the maximum iterations is often too large
# 2. Use a preconditioner if possible
# 3. Make sure that the linear system really is symmetric and positive-definite

ext: jl


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@testset "Interpolation" begin

#########################################################################################
@bench @testset "Random dep" begin
nep = nep_gallery("dep0")
λ1,x1 = newton(nep,logger=displaylevel,maxit=40, λ=-0.75, v=ones(size(nep,1)));
@test compute_resnorm(nep,λ1,x1) < eps()*100

# Newton on interpolated dep (interpolating pep of degree 2)
intpoints = [λ1-1, λ1, λ1+1]
pep = interpolate(nep, intpoints)
λ2,x2 = newton(pep,logger=displaylevel,maxit=40, λ=-0.75, v=ones(size(nep,1)));
@test compute_resnorm(pep,λ2,x2) < eps()*100
@test abs(λ1-λ2)/abs(λ1) < eps()*1000

# Newton on interpolated dep (interpolating pep of degree 8)
intpoints = [λ1-5, λ1-1, λ1, λ1+5, λ1+1, λ1+5im, λ1+1im, λ1-5im, λ1-1im]
pep = interpolate(nep, intpoints)
λ2,x2 = newton(pep,logger=displaylevel,maxit=40, λ=-0.75, v=ones(size(nep,1)));
@test compute_resnorm(pep,λ2,x2) < eps()*100
@test abs(λ1-λ2)/abs(λ1) < eps()*1000
end

#########################################################################################
@bench @testset "Random sparse dep" begin
nep = nep_gallery("dep0_sparse", 30)

λ1,x1 = newton(nep,logger=displaylevel,maxit=40, λ=0.2, v=ones(size(nep,1)),armijo_factor=0.5);
@test compute_resnorm(nep,λ1,x1) < eps()*100

# Newton on interpolated sparse dep
intpoints = [λ1-.5, λ1, λ1+.5]
pep = interpolate(nep, intpoints)
λ2,x2 = newton(pep,logger=displaylevel,maxit=40, λ=λ1+0.1, v=x1+0.1*ones(size(nep,1)));
@test compute_resnorm(pep,λ2,x2) < eps()*1000000
@test abs(λ1-λ2)/abs(λ1) < eps()*1000

end

#########################################################################################
@bench @testset "Random pep (original pep of degree 2)" begin
nep = nep_gallery("pep0")
λ1,x1 = newton(nep,logger=displaylevel,maxit=40, λ=1, v=ones(size(nep,1)));
@test compute_resnorm(nep,λ1,x1) < eps()*100

#Newton on interpolated pep (interpolation of degree 2)
intpoints = [λ1-1, λ1, λ1+1]
pep = interpolate(nep, intpoints)
λ2,x2 = newton(pep,logger=displaylevel,maxit=40, λ=1, v=ones(size(nep,1)));
@test compute_resnorm(pep,λ2,x2) < eps()*100
@test abs(λ1-λ2)/abs(λ1) < eps()*100

@testset "Coefficient matrix differences (monomes) deg 2" begin
for i = 1:3
    @test opnorm(nep.A[i]-pep.A[i])/opnorm(nep.A[i]) < eps()*100
end
end

#Newton on interpolated pep (interpolation of degree 4) in real arithmetics
intpoints = [λ1-3, λ1-1, λ1, λ1+1, λ1+3]

pep = interpolate(Float64, nep, intpoints)
λ2,x2 = newton(pep,logger=displaylevel,maxit=40, λ=1, v=ones(size(nep,1)));
@test compute_resnorm(pep,λ2,x2) < eps()*100
@test abs(λ1-λ2)/abs(λ1) < eps()*100

@testset "Coefficient matrix differences (monomes) deg 4" begin
for i = 1:3
    @test opnorm(nep.A[i]-pep.A[i])/opnorm(nep.A[i]) < eps()*100
end
for i = 4:5
    @test opnorm(pep.A[i]) < eps()*100
end
end

end

#########################################################################################
@bench @testset "Random sparse pep (original pep of degree 2)" begin
nep=nep_gallery("pep0_sparse")
λ1,x1 =newton(nep,logger=displaylevel,maxit=40, λ=-0.75, v=ones(size(nep,1)));
@test compute_resnorm(nep,λ1,x1) < eps()*1000

#Newton on interpolated sparse pep
intpoints = [λ1-1, λ1, λ1+1.5]
pep = interpolate(nep, intpoints)
λ2,x2 =newton(pep,logger=displaylevel,maxit=40, λ=-0.75, v=ones(size(nep,1)));
@test compute_resnorm(pep,λ2,x2) < eps()*1000
@test abs(λ1-λ2)/abs(λ1) < eps()*100

@testset "Coefficient matrix differences (monomes) deg 2" begin
for i = 1:3
 @test opnorm(nep.A[i]-pep.A[i],Inf)/opnorm(nep.A[i],Inf) < eps()*100
end
end

end

end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module GeneralizedBinomial
using Nemo, PrimesIterator, NumberTheory, Products, Triangles

export ModuleGeneralizedBinomial
export Binomial, Pascal, T007318, Multinomial

"""

 P. Luschny, [Generalized Binomial](http://oeis.org/wiki/User:Peter_Luschny/ExtensionsOfTheBinomial), OEIS Wiki.

* Binomial, Pascal, T007318
"""
const ModuleGeneralizedBinomial = ""

"""

The classical binomial coefficients defined for ``n≥0`` and ``0≤k≤n`` (a.k.a. Pascal's triangle).
"""
function Pascal(n::Int, k::Int)

    (k == 0 || k == n) && return 1
    if k > div(n, 2)
        k = n - k
    end

    nk = n - k
    factors = fmpz[]
    rootN = isqrt(n)

    # Make use of Kummer's theorem.
    for prime ∈ Primes(2, n)

        if prime > nk
            push!(factors, prime)
            continue
        end

        prime > div(n, 2) && continue

        if prime > rootN
            (n % prime < k % prime) && push!(factors, prime)
            continue
        end

        r, N, K, p = 0, n, k, 1

        while N > 0
            r = N % prime < (K % prime + r) ? 1 : 0
            if r == 1
                p *= prime
            end
            N = div(N, prime)
            K = div(K, prime)
        end

        p > 1 && push!(factors, p)
    end

    ∏(factors)
end

"""

Pascal's triangle.
"""
function T007318(n::Int)
    T = zeros(QQ, div(n * (n + 1), 2))
    j = 1
    for m ∈ 0:n-1, k ∈ 0:m
        T[j] = binomial(ZZ(m), ZZ(k))
        j += 1
    end
    T
end

# See the discussion on
# [Extensions of the Binomial](http://oeis.org/wiki/User:Peter_Luschny/ExtensionsOfTheBinomial).
"""

Return the extended binomial coefficients defined for all ``n ∈ Z`` and ``k ∈ Z``. Behaves like the binomial function in Maple and Mathematica.

``\\binom{n}{k} = \\lim \\limits_{x \\rightarrow 1}(Γ(n + x) / (Γ(k + x) Γ(n - k + x))).``

"""
function Binomial(n::Int, k::Int)
    0 ≤ k ≤ n && return binomial(ZZ(n), ZZ(k))
    k ≤ n < 0 && return binomial(ZZ(-k - 1), ZZ(n - k)) * (-1)^(n - k)
    n < 0 ≤ k && return binomial(ZZ(-n + k - 1), ZZ(k)) * (-1)^k
    ZZ(0)
end

"""

Return the multinomial coefficient of a list.
"""
function Multinomial(lst::Array{Int})
    s = 0
    result = fmpz(1)
    for c ∈ lst
        s += c
        result *= binomial(ZZ(s), ZZ(c))
    end
    result
end

#START-TEST-########################################################

using Test
function test()

    @testset "Binomial" begin
        for n ∈ 0:10, k ∈ 0:n
            @test Binomial(n, k) == div(factorial(n), (factorial(n - k) * factorial(k)))
            @test Binomial(n, k) == Pascal(n, k)
        end
    end
end

function demo()
    ShowAsΔ(T007318(8))
    println()

    for n ∈ -10:10
        println([Binomial(n, k) for k ∈ -10:10])
    end
end

"""

for n ∈ 0:10000 Binomial(2*n,n) end
    0.504729 seconds (10.00 k allocations: 156.266 KiB)
for n ∈ -100:100, k ∈ -100:100 Binomial(n,k) end
    0.008669 seconds (55.55 k allocations: 867.984 KiB)
for k ∈ -10000:10000 Binomial(-5,k) end
    0.005378 seconds (40.00 k allocations: 624.969 KiB)
"""
function perf()
    @time (for n ∈ 0:10000
        Binomial(2n, n)
    end)
    @time (for n ∈ -100:100, k ∈ -100:100
        Binomial(n, k)
    end)
    @time (for k ∈ -10000:10000
        Binomial(-5, k)
    end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# This module contains the implementation of two EMOA algorithms: NSGAII and SMS-EMOA
# Implementation follows the general algorithms presented in:
# - source nsgaii  : https://doi.org/10.1007/3-540-45356-3_83
# - source sms-emoa: https://doi.org/10.1007/978-3-540-31880-4_5

# This project was part of the course Multi-objective Metaheuristics from the master 2 ORO at University of Nantes

module GeneticAlgorithms

# Libraries
using Random

# Structs, tools, debuging, indicators computation...
include("structs.jl")
include("printing.jl")
include("debug.jl")
include("utils.jl")
include("solution.jl")
include("indicator.jl")

# GA operators
include("initialisation.jl")
include("selection.jl")
include("crossover.jl")
include("mutation.jl")
include("reparation.jl")

# EMOA algorithms
include("nsgaii.jl")
include("sms_emoa.jl")

# INTRO #####################################################

debug = false       # wether of not assertions are checked (expensive, use only for debugging)

# CODE #######################################################

# returns the configuration considered best by our experiments
function get_optimal_config(prob::_MOMKP, algorithm)
    size, crossover, mutation_rate = 0,0,0

    if algorithm == NSGAII_solve || algorithm == SMS_EMOA_solve
        if prob.n < 50
            size = 75
            crossover = CROSSOVER_binary_uniform
            mutation_rate = 0.05
        elseif prob.n < 150
            size = 150
            crossover = CROSSOVER_binary_uniform
            mutation_rate = 0.05
        else
            size = 150
            crossover = CROSSOVER_two_points
            mutation_rate = 0.2
        end
    else
        println("ERROR : algorithm not recognised")
        return false
    end

    config = Config(
        size,
        crossover,
        mutation_rate,
        REPARATION_advanced!
    )

    return config
end

function solve(prob::_MOMKP; max_t = 15.0, algorithm_solve = NSGAII_solve, seeding_sols = Vector{Sol}(undef, 0))

    config = get_optimal_config(prob, algorithm_solve)
    front = algorithm_solve(prob, config, seeding_sols = seeding_sols, max_t = max_t, verbose = false)

    return front

end

end


text:
# Allow pass --progress when running this script individually to turn on progress meter
const PROGRESS = length(ARGS) > 0 && ARGS[1] == "--progress" ? true : false

using Test, AdvancedHMC, LinearAlgebra, Random
using Parameters: reconstruct
using Statistics: mean, var, cov
include("common.jl")

θ_init = rand(MersenneTwister(1), D)
ϵ = 0.2
n_steps = 10
n_samples = 12_000
n_adapts = 2_000

function test_stats(::Union{StaticTrajectory,HMCDA}, stats, n_adapts)
    for name in (:step_size, :nom_step_size, :n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :is_adapt)
        @test all(map(s -> in(name, propertynames(s)), stats))
    end
    is_adapts = getproperty.(stats, :is_adapt)
    @test is_adapts[1:n_adapts] == ones(Bool, n_adapts)
    @test is_adapts[(n_adapts+1):end] == zeros(Bool, length(stats) - n_adapts)
end

function test_stats(::NUTS, stats, n_adapts)
    for name in (:step_size, :nom_step_size, :n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :is_adapt, :max_hamiltonian_energy_error, :tree_depth, :numerical_error)
        @test all(map(s -> in(name, propertynames(s)), stats))
    end
    is_adapts = getproperty.(stats, :is_adapt)
    @test is_adapts[1:n_adapts] == ones(Bool, n_adapts)
    @test is_adapts[(n_adapts+1):end] == zeros(Bool, length(stats) - n_adapts)
end

@testset "HMC and NUTS" begin
    @testset "$metricsym" for (metricsym, metric) in Dict(
        :UnitEuclideanMetric => UnitEuclideanMetric(D),
        :DiagEuclideanMetric => DiagEuclideanMetric(D),
        :DenseEuclideanMetric => DenseEuclideanMetric(D),
    )
        @test show(metric) == nothing
        h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
        @testset "$lfsym" for (lfsym, lf) in Dict(
            :Leapfrog => Leapfrog(ϵ),
            :JitteredLeapfrog => JitteredLeapfrog(ϵ, 1.0),
            :TemperedLeapfrog => TemperedLeapfrog(ϵ, 1.05),
        )
            @testset "$τsym" for (τsym, τ) in Dict(
                :(StaticTrajectory{EndPointTS}) => StaticTrajectory{EndPointTS}(lf, n_steps),
                :(StaticTrajectory{MultinomialTS}) => StaticTrajectory{MultinomialTS}(lf, n_steps),
                :(HMCDA{EndPointTS}) => HMCDA{EndPointTS}(lf, ϵ * n_steps),
                :(HMCDA{MultinomialTS}) => HMCDA{MultinomialTS}(lf, ϵ * n_steps),
                :(NUTS{SliceTS,Original}) => NUTS{SliceTS,ClassicNoUTurn}(lf),
                :(NUTS{SliceTS,Generalised}) => NUTS{SliceTS,GeneralisedNoUTurn}(lf),
                :(NUTS{MultinomialTS,Original}) => NUTS{MultinomialTS,ClassicNoUTurn}(lf),
                :(NUTS{MultinomialTS,Generalised}) => NUTS{MultinomialTS,GeneralisedNoUTurn}(lf),
            )
                @test show(h) == nothing
                @test show(τ) == nothing
                @testset  "NoAdaptation" begin
                    Random.seed!(1)
                    samples, stats = sample(h, τ, θ_init, n_samples; verbose=false, progress=PROGRESS)
                    @test mean(samples[n_adapts+1:end]) ≈ zeros(D) atol=RNDATOL
                end

                # Skip adaptation tests with tempering
                if lf isa TemperedLeapfrog
                    @info "Adaptation tests for $τsym with $lfsym on $metricsym are skipped"
                    continue
                end

                @testset "$adaptorsym" for (adaptorsym, adaptor) in Dict(
                    :MassMatrixAdaptorOnly => MassMatrixAdaptor(metric),
                    :StepSizeAdaptorOnly => StepSizeAdaptor(0.8, τ.integrator),
                    :NaiveHMCAdaptor => NaiveHMCAdaptor(
                        MassMatrixAdaptor(metric),
                        StepSizeAdaptor(0.8, τ.integrator),
                    ),
                    :StanHMCAdaptor => StanHMCAdaptor(
                        MassMatrixAdaptor(metric),
                        StepSizeAdaptor(0.8, τ.integrator),
                    ),
                )
                    @test show(adaptor) == nothing
                    Random.seed!(1)
                    # For `MassMatrixAdaptor`, we use the pre-defined step size as the method cannot adapt the step size.
                    # For other adapatation methods that are able to adpat the step size, we use `find_good_stepsize`.
                    τ_used = adaptorsym == :MassMatrixAdaptorOnly ? τ : reconstruct(τ, integrator=reconstruct(lf, ϵ=find_good_stepsize(h, θ_init)))
                    samples, stats = sample(h, τ_used , θ_init, n_samples, adaptor, n_adapts; verbose=false, progress=PROGRESS)
                    @test mean(samples) ≈ zeros(D) atol=RNDATOL
                    test_stats(τ_used, stats, n_adapts)
                end
            end
        end
    end
end

@testset "drop_warmup" begin
    metric = DiagEuclideanMetric(D)
    h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
    τ = NUTS(Leapfrog(ϵ))
    adaptor = StanHMCAdaptor(
        MassMatrixAdaptor(metric),
        StepSizeAdaptor(0.8, τ.integrator),
    )
    samples, stats = sample(h, τ, θ_init, n_samples, adaptor, n_adapts; verbose=false, progress=false, drop_warmup=true)
    @test length(samples) == n_samples - n_adapts
    @test length(stats) == n_samples - n_adapts
    samples, stats = sample(h, τ, θ_init, n_samples, adaptor, n_adapts; verbose=false, progress=false, drop_warmup=false)
    @test length(samples) == n_samples
    @test length(stats) == n_samples
end


text:
using ComparativeAutograder

test_cases = []
for _ in 1:20
    n = rand(20:30)
    push!(test_cases, TestCase((rand(n, n), )))
end

runtestsuite(TestSuite(
    "adjacency2modularity",
    test_cases,
))


text:
# spin-unpolarized version
function calc_Focc(evals, Nelectrons, kT::Float64; is_spinpol=false, verbose=false)

    Nstates = length(evals)
    TOL = 1e-10
    MAXITER = 100

    Focc = zeros(Nstates)
    Nocc = round(Int64,Nelectrons/2)  # normally
    if verbose
        @printf("Nelectrons = %d\n", Nelectrons)
        @printf("Nocc = %d\n", Nocc)
    end

    # use bisection to find E_fermi such that
    #  sum_i Focc(i) = Nelectrons
    if Nstates > Nocc
        ilb = Nocc - 1
        iub = Nocc + 1
        # FIXME: Need to guard against the case Nocc == 1
        lb = evals[ilb]
        ub = evals[iub]
        # make sure flb < Nelectrons and fub > Nelectrons
        # Use lb and ub as guess interval for searching Fermi energy
        Focc_lb = smear_FD(evals, lb, kT, is_spinpol=is_spinpol)
        Focc_ub = smear_FD(evals, ub, kT, is_spinpol=is_spinpol)
        flb = sum(Focc_lb)
        fub = sum(Focc_ub)
        while ( (Nelectrons-flb)*(fub-Nelectrons) < 0 )
            if verbose
                @printf("calc_Focc: initial bounds are off:\n");
                @printf("flb = %18.10f, fub = %18.10f, Nelectrons = %d\n", flb, fub, Nelectrons)
            end
            if (flb > Nelectrons)
                if (ilb > 1)
                    ilb = ilb - 1
                    lb = evals[ilb]
                    Focc = smear_FD(evals, lb, kT, is_spinpol=is_spinpol)
                    flb = sum(Focc)
                else
                    @printf("ERROR in calc_Focc: cannot find a lower bound for E_fermi\n")
                    exit()
                end
            end
            #
            if (fub < Nelectrons)
                if (iub < Nstates)
                    iub = iub + 1
                    ub  = evals[iub]
                    Focc = smear_FD(evals, ub, kT, is_spinpol=is_spinpol)
                    fub = sum(Focc)
                else
                    @printf("ERROR in calc_Focc: cannot find an upper bound for E_fermi\n")
                    @printf("Try increasing the number of states\n")
                    exit()
                end
            end
        end  # while

        if verbose
            @printf("\nInitial bounds are OK: flb = %18.10f, fub = %18.10f\n", flb, fub)
        end

        E_fermi = (lb + ub)/2
        Focc = smear_FD(evals, E_fermi, kT, is_spinpol=is_spinpol)
        occsum = sum(Focc)

        for iter = 1:MAXITER
            diffNelec = abs(occsum-Nelectrons)
            if diffNelec < TOL
                if verbose
                    @printf("Found E_fermi = %18.10f, occsum = %18.10f\n", E_fermi, occsum )
                end
                return Focc, E_fermi
            end
            if verbose
                @printf("%3d %18.10f %18.10f %18.10e\n", iter, E_fermi, occsum, diffNelec)
            end
            if (occsum < Nelectrons)
                lb = E_fermi
            else
                ub = E_fermi
            end
            E_fermi = (lb + ub)/2
            Focc = smear_FD(evals, E_fermi, kT, is_spinpol=is_spinpol)
            occsum = sum(Focc)
        end #
        @printf("WARNING: Maximum iteration achieved, E_fermi is not found within specified tolarance\n")
        return Focc, E_fermi

    #
    elseif (Nstates == Nocc)
        @printf("calc_Focc: Nstates is equal to Nocc\n")
        if is_spinpol
            Focc = 2.0*ones(Nstates)
            E_fermi = evals[Nstates]
        else
            Focc    = ones(Nstates)
            E_fermi = evals[Nstates]
        end
        return Focc, E_fermi

    else
        @printf("ERROR: The number of eigenvalues in evals should be larger than Nelectrons")
        exit()
    end

end


text:
export @metal

# Match Darwin version to MacOS version only caring about M1 release and after
# Following: https://en.wikipedia.org/wiki/Darwin_(operating_system)#History
const darwin_to_macos = Dict(
                        # Catalina
                        v"19.2.0" => v"10.15.2",
                        v"19.3.0" => v"10.15.3",
                        v"19.4.0" => v"10.15.4",
                        v"19.5.0" => v"10.15.5",
                        v"19.6.0" => v"10.15.6",
                        # Big Sur
                        v"20.0.0" => v"11.0.0",
                        v"20.1.0" => v"11.0.0",
                        v"20.2.0" => v"11.1.0",
                        v"20.3.0" => v"11.2.0",
                        v"20.4.0" => v"11.3.0",
                        v"20.5.0" => v"11.4.0",
                        v"20.6.0" => v"11.5.0",
                        # Monterey
                        v"21.0.0" => v"12.0.0",
                        v"21.0.1" => v"12.0.0",
                        v"21.1.0" => v"12.0.1",
                        v"21.2.0" => v"12.1.0",
                        v"21.3.0" => v"12.2.0",
                        v"21.4.0" => v"12.3.0")

"""
    @metal [kwargs...] func(args...)

High-level interface for executing code on a GPU. The `@metal` macro should prefix a call,
with `func` a callable function or object that should return nothing. It will be compiled to
a Metal function upon first use, and to a certain extent arguments will be converted and
managed automatically using `mtlconvert`. Finally, a call to `mtlcall` is
performed, creating a command buffer in the current global command queue then committing it.

There is one supported keyword argument that influences the behavior of `@metal`.
- `launch`: whether to launch this kernel, defaults to `true`. If `false` the returned
  kernel object should be launched by calling it and passing arguments again.
"""
macro metal(ex...)
    call = ex[end]
    kwargs = ex[1:end-1]

    # destructure the kernel call
    Meta.isexpr(call, :call) || throw(ArgumentError("second argument to @metal should be a function call"))
    f = call.args[1]
    args = call.args[2:end]

    code = quote end
    vars, var_exprs = assign_args!(code, args)

    # group keyword argument
    macro_kwargs, compiler_kwargs, call_kwargs, other_kwargs =
        split_kwargs(kwargs, [:launch], [:name], [:grid, :threads])
    if !isempty(other_kwargs)
        key,val = first(other_kwargs).args
        throw(ArgumentError("Unsupported keyword argument '$key'"))
    end

    # handle keyword arguments that influence the macro's behavior
    launch = true
    for kwarg in macro_kwargs
        key,val = kwarg.args
        if key == :launch
            isa(val, Bool) || throw(ArgumentError("`launch` keyword argument to @metal should be a Bool"))
            launch = val::Bool
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end
    if !launch && !isempty(call_kwargs)
        error("@metal with launch=false does not support launch-time keyword arguments; use them when calling the kernel")
    end

    # FIXME: macro hygiene wrt. escaping kwarg values (this broke with 1.5)
    #        we esc() the whole thing now, necessitating gensyms...
    @gensym f_var kernel_f kernel_args kernel_tt kernel

    # convert the arguments, call the compiler and launch the kernel
    # while keeping the original arguments alive
    push!(code.args,
        quote
            $f_var = $f
            GC.@preserve $(vars...) $f_var begin
                local $kernel_f = $mtlconvert($f_var)
                local $kernel_args = map($mtlconvert, ($(var_exprs...),))
                local $kernel_tt = Tuple{map(Core.Typeof, $kernel_args)...}
                local $kernel = $mtlfunction($kernel_f, $kernel_tt; $(compiler_kwargs...))
                if $launch
                    $kernel($(var_exprs...); $(call_kwargs...))
                end
                $kernel
            end
         end)
    return esc(code)
end

abstract type AbstractKernel{F,TT} end

@generated function call(kernel::AbstractKernel{F,TT}, args...; call_kwargs...) where {F,TT}
    sig = Base.signature_type(F, TT)
    args = (:F, (:( args[$i] ) for i in 1:length(args))...)

    # filter out ghost arguments that shouldn't be passed
    predicate = if VERSION >= v"1.5.0-DEV.581"
        dt -> isghosttype(dt) || Core.Compiler.isconstType(dt)
    else
        dt -> isghosttype(dt)
    end
    to_pass = map(!predicate, sig.parameters)
    call_t =                  Type[x[1] for x in zip(sig.parameters,  to_pass) if x[2]]
    call_args = Union{Expr,Symbol}[x[1] for x in zip(args, to_pass)            if x[2]]

    # replace non-isbits arguments (they should be unused, or compilation would have failed)
    for (i,dt) in enumerate(call_t)
        if !isbitstype(dt)
            call_t[i] = Ptr{Any}
            call_args[i] = :C_NULL
        end
    end

    # finalize types
    call_tt = Base.to_tuple_type(call_t)

    quote
        Base.@_inline_meta

        mtlcall(kernel.fun, $call_tt, $(call_args...); call_kwargs...)
    end
end

# Why is this all necessary? MtlFunction has the lib and function handle.
struct MtlKernel
    device::MtlDevice
    lib::MtlLibrary
    fun::MtlFunction

    function MtlKernel(dev::MtlDevice, lib::MtlLibrary, fun_name::String)
        fun = MtlFunction(lib, fun_name)
        return new(dev, lib, fun)
        # TODO: Finalizer
    end
end

function MtlKernel(dev::MtlDevice, image::Vector{UInt8}, fun_name::String)
    lib = MtlLibraryFromData(dev, image)
    return MtlKernel(dev, lib, fun_name)
end

# ## host-side kernels

struct HostKernel{F,TT} <: AbstractKernel{F,TT}
    fun::MtlKernel
end

## host-side API

# Get MacOS version from Darwin version
function get_macos_v()
    machine = Sys.MACHINE
    darwin_v = VersionNumber(machine[findfirst("darwin", machine)[end]+1:end])
    # Check for unsupported/incomplete kernel version
    if !(darwin_v in keys(darwin_to_macos))
        error("Unsupported kernel version of $darwin_v")
    end
    return darwin_to_macos[darwin_v]
end

"""
    mtlfunction(f, tt=Tuple{}; kwargs...)

Low-level interface to compile a function invocation for the currently-active GPU, returning
a callable kernel object. For a higher-level interface, use [`@metal`](@ref).

The output of this function is automatically cached, i.e. you can simply call `mtlfunction`
in a hot path without degrading performance. New code will be generated automatically when
the function changes, or when different types or keyword arguments are provided.
"""
function mtlfunction(f::Core.Function, tt::Type=Tuple{}; name=nothing, kwargs...)
    dev = MtlDevice(1)
    cache = get!(()->Dict{UInt,Any}(), mtlfunction_cache, dev)
    source = FunctionSpec(f, tt, true, name)
    target = MetalCompilerTarget(macos=get_macos_v(); kwargs...)
    params = MetalCompilerParams()
    job = CompilerJob(target, source, params)
    GPUCompiler.cached_compilation(cache, job,
                                   mtlfunction_compile, mtlfunction_link)::HostKernel{f,tt}
end

const mtlfunction_cache = Dict{Any,Any}()

function mtlfunction_compile(@nospecialize(job::CompilerJob))
    # TODO: on 1.9, this actually creates a context. cache those.
    JuliaContext() do ctx
        mtlfunction_compile(job, ctx)
    end
end
function mtlfunction_compile(@nospecialize(job::CompilerJob), ctx::Context)
    mi, mi_meta = GPUCompiler.emit_julia(job)
    ir, ir_meta = GPUCompiler.emit_llvm(job, mi; ctx)
    image, asm_meta = GPUCompiler.emit_asm(job, ir; strip=false, format=LLVM.API.LLVMObjectFile) # TODO: Undo strip eventually
    return (image, entry=LLVM.name(ir_meta.entry))
end

function mtlfunction_link(@nospecialize(job::CompilerJob), compiled)
    dev = MtlDevice(1)
    kernel = MtlKernel(dev, compiled.image, compiled.entry)
    return HostKernel{job.source.f,job.source.tt}(kernel)
end

function (kernel::HostKernel)(args...; kwargs...)
    call(kernel, map(mtlconvert, args)...; kwargs...)
end

#launchKernel

##########################################
mtlcall(f::MtlKernel, types::Tuple, args...; kwargs...) =
    mtlcall(f, Base.to_tuple_type(types), args...; kwargs...)

@inline function mtlcall(f::MtlKernel, types::Type, args...; kwargs...)
    # Handle kwargs here??
    cmdq = global_queue(f.lib.device)

    # Process kernel call arguments
    grid    = 1
    threads = 1
    for kwarg in kwargs
        key = kwarg.first
        if key == :grid
            grid = kwarg.second
        elseif key == :threads
            threads = kwarg.second
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end

    cmdbuf = MtlCommandBuffer(cmdq)
    MtlComputeCommandEncoder(cmdbuf) do cce
        enqueue_function(f, args...; grid, threads, cce)
    end
    commit!(cmdbuf)
end

##############################################
# Enqueue a function for execution
#
function enqueue_function(f::MtlKernel, args...;
                grid::MtlDim=1, threads::MtlDim=1,
                cce::MtlComputeCommandEncoder)
    grid = MtlDim3(grid)
    threads = MtlDim3(threads)
    (grid.width>0 && grid.height>0 && grid.depth>0)    || throw(ArgumentError("Grid dimensions should be non-null"))
    (threads.width>0 && threads.height>0 && threads.depth>0) || throw(ArgumentError("Threadgroup dimensions should be non-null"))

    pipeline_state = MtlComputePipelineState(f.device, f.fun)
    (threads.width * threads.height * threads.depth) > pipeline_state.maxTotalThreadsPerThreadgroup && throw(ArgumentError("Max total threadgroup size should not exceed $(pipeline_state.maxTotalThreadsPerThreadgroup)"))
    # Set the function that we are currently encoding
    MTL.set_function!(cce, pipeline_state)
    # Encode all arguments
    encode_arguments!(cce, f.fun, args...)
    # Flush everything
    MTL.append_current_function!(cce, grid, threads)
    return nothing
end

#####
function encode_arguments!(cce::MtlComputeCommandEncoder, f::MtlFunction, args...)
    for (i, a) in enumerate(args)
        encode_argument!(cce, f, i, a)
    end
end

function encode_argument!(cce::MtlComputeCommandEncoder, f::MtlFunction, idx::Integer, arg::Nothing)
    @assert idx > 0
    #@check api.clSetKernelArg(k.id, cl_uint(idx-1), sizeof(CL_mem), C_NULL)
    MTL.set_bytes!(cce, sizeof(C_NULL), C_NULL, idx)
    return cce
end

function encode_argument!(cce::MtlComputeCommandEncoder, f::MtlFunction, idx::Integer, arg::MtlBuffer)
    @assert idx > 0
    set_buffer!(cce, arg, 0, idx)
    return cce
end

function encode_argument!(cce::MtlComputeCommandEncoder, f::MtlFunction, idx::Integer, arg::Core.LLVMPtr)
    @assert idx > 0

    set_buffer!(cce, MtlBuffer{Float32}(Base.bitcast(MTL.MTLBuffer, arg)), 0, idx)
    return cce
end

function encode_argument!(enc::MTL.MtlComputeCommandEncoder, f::MtlFunction, idx::Integer, val::T) where T
    @assert idx > 0 "Kernel idx must be bigger 0"
    #if does not contain a buffer we can use setbytes
    if !contains_mtlbuffer(T)
        ref, tsize = to_mtl_ref(val)
        MTL.set_bytes!(enc, ref, tsize, idx)
    else
        #otherwise, we need an argument buffer
        throw("Not implemented: If an argument contains a mtlbuffer automatic argument encoding is not yet supported.")

        # create an encoder to write into the argument buffer
        argbuf_enc = MtlArgumentEncoder(f, idx)
        # allocate the argument buffer
        argbuf = alloc(Cchar, device(enc), sizeof(argbuf_enc), storage=Shared)
        # assign the argument buffer to the encoder
        MTL.assign_argument_buffer!(argbuf_enc, argbuf, 1)

        # TODO Implement automatic conversion of struct into a assign_argument_buffer

        #for field in val
        #    MTL.set_field!(argbuf_enc, size(val), 1)
        #    set_buffer!(argbuf_enc, pointer(val), 0, 2)
        #end

        # set argubuf_enc into cce
    end
    return enc
end

# Encode MtlDeviceArray using an argument buffer
function encode_argument!(cce::MTL.MtlComputeCommandEncoder, f::MtlFunction, idx::Integer, val::MtlDeviceArray{T}) where T
    #@assert contains_mtlbuffer(typeof(val))
    # create an encoder to write into the argument buffer
    argbuf_enc = MtlArgumentEncoder(f, idx)
    # allocate the argument buffer
    argbuf = alloc(Cchar, device(cce), sizeof(argbuf_enc), storage=Shared)
    # assign the argument buffer to the arg buff encoder
    MTL.assign_argument_buffer!(argbuf_enc, argbuf, 1)

    # Convert LLVMPtr to MtlBuffer
    mtl_buf = MtlBuffer{T}(Base.bitcast(MTL.MTLBuffer, val.ptr))
    # encode the buffer into the argument buffer
    set_buffer!(argbuf_enc, mtl_buf, 0, 1)
    # Encode the size of the MtlDeviceArray into the argument buffer
    MTL.set_field!(argbuf_enc, size(val), 2)
    # Set the device array usage for read/write TODO: Handle constant arrays
    MTL.use!(cce, mtl_buf, MTL.ReadWriteUsage)

    # Set the argument buffer at given argument index
    set_buffer!(cce, argbuf, 0, idx)
    #TODO memmgmt: Leaked temporary argument buffer (argbuf)
    return cce
end


text:
using ParallelRandomFields
using Test

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


text:
using Plots, Distributions

value_α = [0.5, 1.0, 2.0, 4.0]
value_β = [0.5, 1.0, 2.0, 4.0]

l = @layout [grid(4, 4)]
plots = []
for a in value_α, b in value_β
    if a == 0.5
        title = "β=$b"
    else
        title = ""
    end
    if b == 0.5
        ylabel = "α=$a"
        yticks = [0.0, 2.5]
    else
        ylabel = ""
        yticks = nothing
    end
    if a == 4.0
        xticks = [0.0, 0.5, 1.0]
    else
        xticks = nothing
    end
    p = plot(Beta(a, b), label=false, xlims=(0.0, 1.0), ylims=(0.0, 4.5), title=title, ylabel=ylabel, xticks=xticks, yticks=yticks, grid = false)
    push!(plots, p)
end
plot(plots..., layout=l)


text:
module SurveyDesign
include("SamplingProcedures.jl")

struct SurveyDesign
end

end


text:
#ENV["MOSEKLM_LICENSE_FILE"]="/home/ryancw/mosek/mosek.lic" #needed on cluster
using JuMP, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Compat, Test, Random
include("exact_approaches.jl")
include("alt_min.jl")
include("convexpenalties.jl")

results_template = DataFrame(n=Int[], r=Int[], p=Real[], gamma=Real[], err=Real[], rnk=Int[], solvetime=Real[], nc=Int[], gap=Real[], cuts=Int[])

results_exact=similar(results_template, 0)

for ARG in ARGS
    for iter in 1:1
        p=0.2
        n=20
        r=1
        seeds=1:20
        gammaMults=[100.0, 500.0]
        array_num = parse(Int, ARG)

        theSeed=seeds[(array_num-1)÷2+1]
        gamma=gammaMults[(array_num-1)%2+1]
        Random.seed!(theSeed)
        M_L=randn(n, r)
        M_L=round.(1000*M_L)./1000
        M_R=randn(n, r)
        M_R=round.(1000*M_R)./1000
        M=M_L*M_R'

        @show maximum(M)
        numEntriesSampled=Int(round(p*n^2))
        coords=[repeat(1:n, inner=size(1:n,1)) repeat(1:n, outer=size(1:n,1))];
        entries=sample(1:n^2, Int(numEntriesSampled), replace=false)

        theIndices=zeros(n,n)
        for t=1:Int(numEntriesSampled)
            theIndices[coords[entries[t],1], coords[entries[t],2]]=1.0
        end

        k=r

        t_GD=@elapsed X_GD,LB, UB, Y_relax=getMatrixFrobeniusNorm_greedy(coords, entries, M, n, numEntriesSampled, k, gamma, true)

        @show "SDP lower bound is: " LB
        u,sigma,v,=svd(X_GD)
        rnk_GD=sum(sigma.>1e-2);
        err_GD = norm(X_GD.-M)/norm(M);

        t_AM=@elapsed X_AM,=alternatingminimization_matrixcompletion(coords, entries, M, n, numEntriesSampled, k, gamma, 50, X_GD)

        # Test solution on overall data.
        u,sigma,v,=svd(X_AM)
        @show rnk_AM=sum(sigma.>1e-2);
        @show err_AM = norm(X_AM.-M)/norm(M);

        u_0=u[:,1:k]
        v_0=(Diagonal(sigma[1:k])*v[:,1:k]')'

        t_EX=@elapsed X_EX, nc_EX, gap_EX=getMatrixFrobeniusNorm_exact(coords, entries, M, n, numEntriesSampled, k, gamma, false, u_0, v_0)
        u,v,=svd(X_EX)
        @show rnk_EX=sum(v.>1e-2);
        @show err_EX = norm(X_EX.-M)/norm(M);

        push!(results_exact, [n, r, p, gamma, err_EX, rnk_EX, t_EX, nc_EX, gap_EX, theSeed])
        CSV.write("sensitivity_analysis_grb_p"*string(p)*"_n_"*string(n)*".csv", results_exact, append=true)

        end
    end

max_stars_repo_path: Matrix Completion/createFigure6_grb.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_path: Matrix Completion/createFigure6_grb.jl

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_path: Matrix Completion/createFigure6_grb.jl

max_forks_repo_name: ryancorywright/MixedProjectionSoftware


text:
set(obj, ::typeof(@lens last(_)), val) = @set obj[lastindex(obj)] = val
set(obj, ::typeof(@lens first(_)), val) = @set obj[firstindex(obj)] = val

set(obj::Array, ::typeof(@lens eltype(_)), T::Type) = collect(T, obj)
set(::Type{<:Array{<:Any, N}}, ::typeof(@lens eltype(_)), ::Type{T}) where {N, T} =
    Array{T, N}

set(obj::Dict, l::FunctionLens, T::Type) = set(typeof(obj), l, T)(obj)
set(::Type{<:Dict}, ::typeof(@lens eltype(_)), ::Type{Pair{K, V}}) where {K, V} =
    Dict{K, V}
set(::Type{<:Dict{<:Any,V}}, ::typeof(@lens keytype(_)), ::Type{K}) where {K, V} =
    Dict{K, V}
set(::Type{<:Dict{K}}, ::typeof(@lens valtype(_)), ::Type{V}) where {K, V} =
    Dict{K, V}


text:
# MAT_HDF5.jl
# Tools for reading MATLAB HDF5 (v7.3) files in Julia
#
# Copyright (C) 2012   Timothy E. Holy
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

###########################################
## Reading and writing MATLAB .mat files ##
###########################################

module MAT_HDF5

using HDF5, SparseArrays

import Base: read, write, close
import HDF5: names, exists, HDF5ReferenceObj, HDF5BitsKind

const HDF5Parent = Union{HDF5File, HDF5Group}
const HDF5BitsOrBool = Union{HDF5BitsKind,Bool}

mutable struct MatlabHDF5File <: HDF5.DataFile
    plain::HDF5File
    toclose::Bool
    writeheader::Bool
    refcounter::Int
    compress::Bool

    function MatlabHDF5File(plain, toclose::Bool=true, writeheader::Bool=false, refcounter::Int=0,compress::Bool=false)
        f = new(plain, toclose, writeheader, refcounter,compress)
        if toclose
            finalizer(close, f)
        end
        f
    end
end

"""
    close(matfile_handle)

Close a Matlab file.
"""
function close(f::MatlabHDF5File)
    if f.toclose
        close(f.plain)
        if f.writeheader
            magic = zeros(UInt8, 512)
            identifier = "MATLAB 7.3 MAT-file" # minimal but sufficient
            GC.@preserve magic identifier begin
                magicptr = pointer(magic)
                idptr = pointer(identifier)
                unsafe_copyto!(magicptr, idptr, length(identifier))
            end
            magic[126] = 0x02
            magic[127] = 0x49
            magic[128] = 0x4d
            rawfid = open(f.plain.filename, "r+")
            write(rawfid, magic)
            close(rawfid)
        end
        f.toclose = false
    end
    nothing
end

function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool, compress::Bool)
    local f
    if ff && !wr
        error("Cannot append to a write-only file")
    end
    if !cr && !isfile(filename)
        error("File ", filename, " cannot be found")
    end
    pa = p_create(HDF5.H5P_FILE_ACCESS)
    pa["fclose_degree"] = HDF5.H5F_CLOSE_STRONG
    if cr && (tr || !isfile(filename))
        # We're truncating, so we don't have to check the format of an existing file
        # Set the user block to 512 bytes, to save room for the header
        p = p_create(HDF5.H5P_FILE_CREATE)
        p["userblock"] = 512
        f = HDF5.h5f_create(filename, HDF5.H5F_ACC_TRUNC, p.id, pa.id)
        writeheader = true
    else
        f = HDF5.h5f_open(filename, wr ? HDF5.H5F_ACC_RDWR : HDF5.H5F_ACC_RDONLY, pa.id)
        writeheader = false
    end
    close(pa)
    fid = MatlabHDF5File(HDF5File(f, filename), true, writeheader, 0,compress)
    pathrefs = "/#refs#"
    if exists(fid.plain, pathrefs)
        g = fid.plain[pathrefs]
        fid.refcounter = length(g)-1
        close(g)
    end
    fid
end

### Matlab file format specification ###

const name_type_attr_matlab = "MATLAB_class"
const empty_attr_matlab = "MATLAB_empty"
const sparse_attr_matlab = "MATLAB_sparse"
const int_decode_attr_matlab = "MATLAB_int_decode"

### Reading
function read_complex(dtype::HDF5Datatype, dset::HDF5Dataset, ::Type{Array{T}}) where T
    if !check_datatype_complex(dtype)
        close(dtype)
        error("Unrecognized compound data type when reading ", name(dset))
    end
    memtype = build_datatype_complex(T)
    sz = size(dset)
    dbuf = Array{Complex{T}}(undef, sz...)
    HDF5.h5d_read(dset.id, memtype.id, HDF5.H5S_ALL, HDF5.H5S_ALL, HDF5.H5P_DEFAULT, vec(dbuf))
    length(dbuf) == 1 ? dbuf[1] : dbuf
end

function m_read(dset::HDF5Dataset)
    if exists(dset, empty_attr_matlab)
        # Empty arrays encode the dimensions as the dataset
        dims = convert(Vector{Int}, read(dset))
        mattype = a_read(dset, name_type_attr_matlab)
        if mattype == "char"
            return ""
        elseif mattype == "struct"
            # Not sure if this check is necessary but it is checked in
            # `m_read(g::HDF5Group)`
            if exists(dset, "MATLAB_fields")
                return Dict{String,Any}(n=>[] for n in a_read(dset, "MATLAB_fields"))
            else
                return Dict{String,Any}()
            end
        else
            T = mattype == "canonical empty" ? Union{} : str2eltype_matlab[mattype]
            return Array{T}(undef, dims...)
        end
    end

    mattype = exists(dset, name_type_attr_matlab) ? a_read(dset, name_type_attr_matlab) : "cell"

    if mattype == "cell"
        # Cell arrays, represented as an array of refs
        refs = read(dset, Array{HDF5ReferenceObj})
        out = Array{Any}(undef, size(refs))
        f = file(dset)
        for i = 1:length(refs)
            dset = f[refs[i]]
            try
                out[i] = m_read(dset)
            finally
                close(dset)
            end
        end
        return out
    end

    # Regular arrays of values
    # Convert to Julia type
    T = str2type_matlab[mattype]

    # Check for a COMPOUND data set, and if so handle complex numbers specially
    dtype = datatype(dset)
    try
        class_id = HDF5.h5t_get_class(dtype.id)
        d = class_id == HDF5.H5T_COMPOUND ? read_complex(dtype, dset, T) : read(dset, T)
        length(d) == 1 ? d[1] : d
    finally
        close(dtype)
    end
end

function add!(A, x)
    for i = 1:length(A)
        @inbounds A[i] += x
    end
    A
end

# reading a struct, struct array, or sparse matrix
function m_read(g::HDF5Group)
    mattype = a_read(g, name_type_attr_matlab)
    if mattype != "struct"
        # Check if this is a sparse matrix.
        fn = names(g)
        if exists(attrs(g), sparse_attr_matlab)
            # This is a sparse matrix.
            # ir is the row indices, jc is the column boundaries.
            # We add one to account for the zero-based (MATLAB) to one-based (Julia) transition
            jc = add!(convert(Vector{Int}, read(g, "jc")), 1)
            if "data" in fn && "ir" in fn && "jc" in fn
                # This matrix is not empty.
                ir = add!(convert(Vector{Int}, read(g, "ir")), 1)
                dset = g["data"]
                T = str2type_matlab[mattype]
                try
                    dtype = datatype(dset)
                    class_id = HDF5.h5t_get_class(dtype.id)
                    try
                        data = class_id == HDF5.H5T_COMPOUND ? read_complex(dtype, dset, T) : read(dset, T)
                    finally
                        close(dtype)
                    end
                finally
                    close(dset)
                end
            else
                # This matrix is empty.
                ir = Int[]
                data = str2eltype_matlab[mattype][]
            end
            return SparseMatrixCSC(convert(Int, HDF5.a_read(g, sparse_attr_matlab)), length(jc)-1, jc, ir, data)
        else
            error("Cannot read from a non-struct group, type was $mattype")
        end
    end
    if exists(g, "MATLAB_fields")
        fn = a_read(g, "MATLAB_fields")
    else
        fn = names(g)
    end
    s = Dict{String, Any}()
    for i = 1:length(fn)
        dset = g[fn[i]]
        try
            s[fn[i]] = m_read(dset)
        finally
            close(dset)
        end
    end
    s
end

"""
    read(matfile_handle, varname) -> value

Read a variable from an opened Matlab file and return its value.

See `matopen` and `matread`.
"""
function read(f::MatlabHDF5File, name::String)
    local val
    obj = f.plain[name]
    try
        val = m_read(obj)
    finally
        close(obj)
    end
    val
end

"""
    names(matfile_handle) -> Vector{String}

Return a list of variables in an opened Matlab file.

See `matopen`.
"""
names(f::MatlabHDF5File) = filter!(x->x != "#refs#", names(f.plain))

"""
    exists(matfile_handle, varname) -> Bool

Return true if a variable is present in an opened Matlab file.

See `matopen`.
"""
exists(p::MatlabHDF5File, path::String) = exists(p.plain, path)

### Writing

# Check whether a varname is valid for MATLAB
check_valid_varname(s::AbstractString) = if match(r"^[a-zA-Z][a-zA-Z0-9_]*$", s) == nothing
    error("Invalid variable name or key \"$s\": variable names must start with a letter and contain only alphanumeric characters and underscore")
elseif length(s) > 63
    error("Invalid variable name or key \"$s\": variable names must be less than 64 characters")
end

toarray(x::Array) = x
toarray(x::Array{Bool}) = reinterpret(UInt8, x)
toarray(x::Bool) = UInt8[x]
toarray(x) = [x]

# Write the MATLAB type string for dset
m_writetypeattr(dset, ::Type{Complex{T}}) where T = m_writetypeattr(dset, T)
function m_writetypeattr(dset, T)
    if !haskey(type2str_matlab, T)
        error("Type ", T, " is not (yet) supported")
    end
    typename = type2str_matlab[T]

    # Write the attribute
    a_write(dset, name_type_attr_matlab, typename)
    if T == Bool
        a_write(dset, int_decode_attr_matlab, Int32(1))
    end
end

# Writes an empty scalar or array
function m_writeempty(parent::HDF5Parent, name::String, data::AbstractArray)
    adata = [size(data)...]
    dset, dtype = d_create(parent, name, adata)
    try
        a_write(dset, empty_attr_matlab, 0x01)
        m_writetypeattr(dset, eltype(data))
        HDF5.writearray(dset, dtype.id, adata)
    finally
        close(dset)
        close(dtype)
    end
end

# Write an array to a dataset in a MATLAB file, returning the dataset
function m_writearray(parent::HDF5Parent, name::String, adata::AbstractArray{T},compress::Bool) where {T<:HDF5BitsOrBool}
    dset, dtype = compress ? d_create(parent, name, adata,"chunk",HDF5.heuristic_chunk(adata),"compress",3) : d_create(parent, name, adata)
    try
        HDF5.writearray(dset, dtype.id, adata)
        dset
    catch e
        close(dset)
        rethrow(e)
    finally
        close(dtype)
    end
end

function m_writearray(parent::HDF5Parent, name::String, adata::AbstractArray{Complex{T}},compress:Bool) where {T<:HDF5BitsOrBool}
    dtype = build_datatype_complex(T)
    try
        stype = dataspace(adata)
        if compress
            p = p_create(HDF5.H5P_DATASET_CREATE)
            p["compress"] = 3
            p["chunk"]=HDF5.heuristic_chunk(adata)
            obj_id = HDF5.h5d_create(parent.id, name, dtype.id, stype.id,
                                     HDF5._link_properties(name),p.id,HDF5.H5P_DEFAULT)
        else
            obj_id = HDF5.h5d_create(parent.id, name, dtype.id, stype.id)
        end
        dset = HDF5Dataset(obj_id, file(parent))
        try
            arr = reshape(reinterpret(T, adata), tuple(2, size(adata)...))
            HDF5.writearray(dset, dtype.id, arr)
        catch e
            close(dset)
            rethrow(e)
        finally
            close(stype)
        end
        dset
    finally
        close(dtype)
    end
end

# Write a scalar or array
function m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, data::Union{T, Complex{T}, Array{T}, Array{Complex{T}}}) where {T<:HDF5BitsOrBool}
    if isempty(data)
        m_writeempty(parent, name, data)
        return
    end
    dset = m_writearray(parent, name, toarray(data),mfile.compress)
    try
        m_writetypeattr(dset, T)
    finally
        close(dset)
    end
end

# Write sparse arrays
function m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, data::SparseMatrixCSC{T}) where T
    g = g_create(parent, name)
    try
        m_writetypeattr(g, T)
        a_write(g, sparse_attr_matlab, UInt64(size(data, 1)))
        if !isempty(data.nzval)
            close(m_writearray(g, "data", toarray(data.nzval),mfile.compress))
            close(m_writearray(g, "ir", add!(isa(data.rowval, Vector{UInt64}) ? copy(data.rowval) : convert(Vector{UInt64}, data.rowval), typemax(UInt64)),mfile.compress))
        end
        close(m_writearray(g, "jc", add!(isa(data.colptr, Vector{UInt64}) ? copy(data.colptr) : convert(Vector{UInt64}, data.colptr), typemax(UInt64)),mfile.compress))
    finally
        close(g)
    end
end

# Write a string
function m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, str::AbstractString)
    if isempty(str)
        data = UInt64[0, 0]

        # Create the dataset
        dset, dtype = d_create(parent, name, data)
        try
            a_write(dset, name_type_attr_matlab, "char")
            a_write(dset, empty_attr_matlab, 0x01)
            HDF5.writearray(dset, dtype.id, data)
        finally
            close(dset)
            close(dtype)
        end
    else
        # Here we assume no UTF-16
        data = zeros(UInt16, 1, length(str))
        i = 1
        for c in str
            data[i] = c
            i += 1
        end

        # Create the dataset
        dset, dtype = d_create(parent, name, data)
        try
            a_write(dset, name_type_attr_matlab, "char")
            a_write(dset, int_decode_attr_matlab, Int32(2))
            HDF5.writearray(dset, dtype.id, data)
        finally
            close(dset)
            close(dtype)
        end
    end
end

# Write cell arrays
function m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, data::Array{T}) where T
    pathrefs = "/#refs#"
    fid = file(parent)
    local g
    local refs
    if !exists(fid, pathrefs)
        g = g_create(fid, pathrefs)
    else
        g = fid[pathrefs]
    end
    try
        # If needed, create the "empty" item
        if !exists(g, "a")
            edata = zeros(UInt64, 2)
            eset, etype = d_create(g, "a", edata)
            try
                HDF5.writearray(eset, etype.id, edata)
                a_write(eset, name_type_attr_matlab, "canonical empty")
                a_write(eset, "MATLAB_empty", 0x00)
            finally
                close(etype)
                close(eset)
            end
        else
            a = g["a"]
            if !exists(attrs(a), "MATLAB_empty")
                error("Must create the empty item, with name a, first")
            end
            close(a)
        end
        # Write the items to the reference group
        refs = Array{HDF5ReferenceObj}(undef, size(data))
        for i = 1:length(data)
            mfile.refcounter += 1
            itemname = string(mfile.refcounter)
            m_write(mfile, g, itemname, data[i])
            # Extract references
            tmp = g[itemname]
            refs[i] = HDF5ReferenceObj(tmp, pathrefs*"/"*itemname)
            close(tmp)
        end
    finally
        close(g)
    end
    # Write the references as the chosen variable
    cset, ctype = d_create(parent, name, refs)
    try
        HDF5.writearray(cset, ctype.id, refs)
        a_write(cset, name_type_attr_matlab, "cell")
    finally
        close(ctype)
        close(cset)
    end
end

# Check that keys are valid for a struct, and convert them to an array of ASCIIStrings
function check_struct_keys(k::Vector)
    asckeys = Vector{String}(undef, length(k))
    for i = 1:length(k)
        key = k[i]
        if !isa(key, AbstractString)
            error("Only Dicts with string keys may be saved as MATLAB structs")
        end
        check_valid_varname(key)
        asckeys[i] = convert(String, key)
    end
    asckeys
end

# Write a struct from arrays of keys and values
function m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, k::Vector{String}, v::Vector)
    g = g_create(parent, name)
    a_write(g, name_type_attr_matlab, "struct")
    for i = 1:length(k)
        m_write(mfile, g, k[i], v[i])
    end
    a_write(g, "MATLAB_fields", HDF5Vlen(k))
end

# Write Associative as a struct
m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, s::AbstractDict) =
    m_write(mfile, parent, name, check_struct_keys(collect(keys(s))), collect(values(s)))

# Write generic CompositeKind as a struct
function m_write(mfile::MatlabHDF5File, parent::HDF5Parent, name::String, s)
    if isbits(s)
        error("This is the write function for CompositeKind, but the input doesn't fit")
    end
    T = typeof(s)
    m_write(mfile, parent, name, check_struct_keys([string(x) for x in fieldnames(T)]), [getfield(s, x) for x in fieldnames(T)])
end

# Check whether a variable name is valid, then write it
"""
    write(matfile_handle, varname, value)

Write the value into an opened Matlab file as the specified variable.

See `matopen` and `matwrite`.
"""
function write(parent::MatlabHDF5File, name::String, thing)
    check_valid_varname(name)
    m_write(parent, parent.plain, name, thing)
end

## Type conversion operations ##

struct MatlabString end

const str2type_matlab = Dict(
    "canonical empty" => nothing,
    "int8"    => Array{Int8},
    "uint8"   => Array{UInt8},
    "int16"   => Array{Int16},
    "uint16"  => Array{UInt16},
    "int32"   => Array{Int32},
    "uint32"  => Array{UInt32},
    "int64"   => Array{Int64},
    "uint64"  => Array{UInt64},
    "single"  => Array{Float32},
    "double"  => Array{Float64},
    "cell"    => Array{Any},
    "char"    => MatlabString,
    "logical" => Array{Bool}
)
# These operate on the element type rather than the whole type
const str2eltype_matlab = Dict(
    "canonical empty" => nothing,
    "int8"    => Int8,
    "uint8"   => UInt8,
    "int16"   => Int16,
    "uint16"  => UInt16,
    "int32"   => Int32,
    "uint32"  => UInt32,
    "int64"   => Int64,
    "uint64"  => UInt64,
    "single"  => Float32,
    "double"  => Float64,
    "cell"    => Any,
    "char"    => MatlabString,
    "logical" => Bool
)
const type2str_matlab = Dict(
    Int8    => "int8",
    UInt8   => "uint8",
    Int16   => "int16",
    UInt16  => "uint16",
    Int32   => "int32",
    UInt32  => "uint32",
    Int64   => "int64",
    UInt64  => "uint64",
    Float32 => "single",
    Float64 => "double",
    Bool    => "logical"
)

function read(obj::HDF5Object, ::Type{MatlabString})
    T = HDF5.hdf5_to_julia(obj)
    data = read(obj, T)
    if size(data, 1) == 1
        sz = size(data)
        data = reshape(data, sz[2:end])
    end
    if ndims(data) == 1
        return String(convert(Vector{Char}, data))
    elseif ndims(data) == 2
        return datap = String[rstrip(String(convert(Vector{Char}, vec(data[i, :])))) for i = 1:size(data, 1)]
    else
        return data
    end
end
function read(obj::HDF5Object, ::Type{Bool})
    tf = read(obj, UInt8)
    tf > 0
end
function read(obj::HDF5Dataset, ::Type{Array{Bool}})
    if HDF5.isnull(obj)
        return Bool[]
    end
    # Use the low-level HDF5 API to put the data directly into a Bool array
    tf = Array{Bool}(undef, size(obj))
    HDF5.h5d_read(obj.id, HDF5.hdf5_type_id(UInt8), tf, obj.xfer)
    return tf
end

## Utilities for handling complex numbers
function build_datatype_complex(T::Type)
    memtype_id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 2*sizeof(T))
    HDF5.h5t_insert(memtype_id, "real", 0, HDF5.hdf5_type_id(T))
    HDF5.h5t_insert(memtype_id, "imag", sizeof(T), HDF5.hdf5_type_id(T))
    HDF5Datatype(memtype_id)
end

function check_datatype_complex(dtype::HDF5Datatype)
    n = HDF5.h5t_get_nmembers(dtype.id)
    if n != 2
        return false
    end
    if HDF5.h5t_get_member_name(dtype.id, 0) != "real" ||
       HDF5.h5t_get_member_name(dtype.id, 1) != "imag"
        return false
    end
    true
end

end


text: #== In theory this is a not too complicated RNN to add punctuation into some text ==#

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
struct BoundaryDiscreteModel{Dc,Dp} <:DiscreteModel{Dc,Dp}
  model::DiscreteModel{Dc,Dp}
  bgmodel::DiscreteModel
  trian::BoundaryTriangulation
end

function BoundaryDiscreteModel(
  ::Type{<:Polytope{Df}},
  bgmodel::DiscreteModel,
  bgface_to_mask::AbstractArray) where Df

  @assert Df < num_cell_dims(bgmodel)

  bgmodel_Df = DiscreteModel(Polytope{Df},bgmodel)
  model = DiscreteModelPortion(bgmodel_Df,bgface_to_mask)
  trian = BoundaryTriangulation(bgmodel,bgface_to_mask)
  BoundaryDiscreteModel(model,bgmodel,trian)
end

get_grid(model::BoundaryDiscreteModel) = get_grid(model.model)
get_grid_topology(model::BoundaryDiscreteModel) = get_grid_topology(model.model)
get_face_labeling(model::BoundaryDiscreteModel) = get_face_labeling(model.model)
get_triangulation(a::BoundaryDiscreteModel) = a.trian
Triangulation(a::BoundaryDiscreteModel) = a.trian


text:
# This file was generated, do not modify it. # hide
ŷ = predict(tuned_ensemble, rows=test)
@show rms(ŷ, y[test])

PyPlot.close_figs() # hide

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module IsotopicCalc

    using JSON, OrderedCollections

    include("functions.jl")

    export isotopicPattern, isotopicPatternProtonated, monoisotopicMass
    export isopat, isopatprot, monoisomass

end # module


text:
using Gadfly, Distributions

plot(y=rand(Normal(), 100), Geom.boxplot)


text:
const AF_INET = 2
const AF_INET6 = @static if Sys.iswindows() # AF_INET6 in <sys/socket.h>
    23
elseif Sys.isapple()
    30
elseif Sys.KERNEL ∈ (:FreeBSD, :DragonFly)
    28
elseif Sys.KERNEL ∈ (:NetBSD, :OpenBSD)
    24
else
    10
end

const NULL = 0

# /* Socket types. */
const SOCK_STREAM = 1 #  stream (connection) socket
const SOCK_DGRAM = 2 #  datagram (conn.less) socket
const SOCK_RAW = 3 #  raw socket
const SOCK_RDM = 4 #  reliably-delivered message
const SOCK_SEQPACKET = 5 #  sequential packet socket

const SOL_IP = 0
const SOL_TCP = 6
const SOL_UDP = 17
const SOL_IPV6 = 41
const SOL_ICMPV6 = 58
const SOL_SCTP = 132
const IPPROTO_SCTP = 132
const SOL_UDPLITE = 136
const SOL_RAW = 255

const SaFamily = Cint

# struct InAddr
#   s_addr :: Cuint # address in network byte order
# end

# InAddr(address :: IPv4) = InAddr(hton(address.host))

# abstract type SockAddr end

# struct SockaddrIn <: SockAddr
#   sin_family :: Cint # AF_INET
#   sin_port :: Cushort # port in network byte order
#   sin_addr :: InAddr # internet address
# end
# function SockaddrIn(address :: IPv4, port :: Integer)
#   p = UInt16(port)
#   SockaddrIn(AF_INET,
#              hton(p),
#              InAddr(address))
# end
# const SctpAssocT = Cuint

# struct SctpSndinfo
#   # the stream number to which the application wishes to send this message
#   snd_sid :: UInt16
#   # any of the following flags and is composed of bitwise OR of these values
#   snd_flags :: UInt16
#   snd_ppid :: UInt32
#   snd_context :: UInt32
#   snd_assoc_id :: SctpAssocT
# end

# struct SctpUdpencaps
#   # sue_address :: SockaddrStorage
#   sue_address :: SockaddrIn
#   sue_assoc_id :: UInt32
#   sue_port :: UInt16
# end

#  Flags that go into the sinfo->sinfo_flags field
const SCTP_DATA_LAST_FRAG    = 0x0001 #  tail part of the message could not be sent
const SCTP_DATA_NOT_FRAG    = 0x0003 #  complete message could not be sent
const SCTP_NOTIFICATION    = 0x0010 #  next message is a notification
const SCTP_COMPLETE    = 0x0020 #  next message is complete
const SCTP_EOF    = 0x0100 #  Start shutdown procedures
const SCTP_ABORT    = 0x0200 #  Send an ABORT to peer
const SCTP_UNORDERED    = 0x0400 #  Message is un-ordered
const SCTP_ADDR_OVER    = 0x0800 #  Override the primary-address
const SCTP_SENDALL    = 0x1000 #  Send this on all associations
const SCTP_EOR    = 0x2000 #  end of message signal
const SCTP_SACK_IMMEDIATELY    = 0x4000 #  Set I-Bit

# const INVALID_SINFO_FLAG    = (x) (((x) & 0xfffffff0 \
#                                     & ~(SCTP_EOF | SCTP_ABORT | SCTP_UNORDERED |\
# 				        SCTP_ADDR_OVER | SCTP_SENDALL | SCTP_EOR |\
# 					SCTP_SACK_IMMEDIATELY)) != 0)
#  for the endpoint

const SCTP_SENDV_NOINFO     = 0 |> Cuint
const SCTP_SENDV_SNDINFO    = 1 |> Cuint
const SCTP_SENDV_PRINFO     = 2 |> Cuint
const SCTP_SENDV_AUTHINFO   = 3 |> Cuint
const SCTP_SENDV_SPA        = 4 |> Cuint

const SCTP_RECVV_NOINFO     = 0 |> Cuint
const SCTP_RECVV_RCVINFO    = 1 |> Cuint
const SCTP_RECVV_NXTINFO    = 2 |> Cuint
const SCTP_RECVV_RN         = 3 |> Cuint

# user socket options: socket API defined

# read-write options

const SCTP_RTOINFO    = 0x00000001
const SCTP_ASSOCINFO    = 0x00000002
const SCTP_INITMSG    = 0x00000003
const SCTP_NODELAY    = 0x00000004
const SCTP_AUTOCLOSE    = 0x00000005
const SCTP_PRIMARY_ADDR    = 0x00000007
const SCTP_ADAPTATION_LAYER    = 0x00000008
const SCTP_DISABLE_FRAGMENTS    = 0x00000009
const SCTP_PEER_ADDR_PARAMS    = 0x0000000a
#  ancillary data/notification interest options
#  Without this applied we will give V4 and V6 addresses on a V6 socket
const SCTP_I_WANT_MAPPED_V4_ADDR    = 0x0000000d
const SCTP_MAXSEG    = 0x0000000e
const SCTP_DELAYED_SACK    = 0x0000000f
const SCTP_FRAGMENT_INTERLEAVE    = 0x00000010
const SCTP_PARTIAL_DELIVERY_POINT    = 0x00000011
#  authentication support
const SCTP_HMAC_IDENT    = 0x00000014
const SCTP_AUTH_ACTIVE_KEY    = 0x00000015
const SCTP_AUTO_ASCONF    = 0x00000018
const SCTP_MAX_BURST    = 0x00000019
#  assoc level context
const SCTP_CONTEXT    = 0x0000001a
#  explicit EOR signalling
const SCTP_EXPLICIT_EOR    = 0x0000001b
const SCTP_REUSE_PORT    = 0x0000001c

const SCTP_EVENT    = 0x0000001e
const SCTP_RECVRCVINFO    = 0x0000001f
const SCTP_RECVNXTINFO    = 0x00000020
const SCTP_DEFAULT_SNDINFO    = 0x00000021
const SCTP_DEFAULT_PRINFO    = 0x00000022
const SCTP_REMOTE_UDP_ENCAPS_PORT    = 0x00000024
const SCTP_ECN_SUPPORTED    = 0x00000025
const SCTP_PR_SUPPORTED    = 0x00000026
const SCTP_AUTH_SUPPORTED    = 0x00000027
const SCTP_ASCONF_SUPPORTED    = 0x00000028
const SCTP_RECONFIG_SUPPORTED    = 0x00000029
const SCTP_NRSACK_SUPPORTED    = 0x00000030
const SCTP_PKTDROP_SUPPORTED    = 0x00000031
const SCTP_MAX_CWND    = 0x00000032

const SCTP_ENABLE_STREAM_RESET    = 0x00000900 #  struct sctp_assoc_value

#  Pluggable Stream Scheduling Socket option
const SCTP_PLUGGABLE_SS    = 0x00001203
const SCTP_SS_VALUE    = 0x00001204

# read-only options

const SCTP_STATUS    = 0x00000100
const SCTP_GET_PEER_ADDR_INFO    = 0x00000101
#  authentication support
const SCTP_PEER_AUTH_CHUNKS    = 0x00000102
const SCTP_LOCAL_AUTH_CHUNKS    = 0x00000103
const SCTP_GET_ASSOC_NUMBER    = 0x00000104
const SCTP_GET_ASSOC_ID_LIST    = 0x00000105
const SCTP_TIMEOUTS    = 0x00000106
const SCTP_PR_STREAM_STATUS    = 0x00000107
const SCTP_PR_ASSOC_STATUS    = 0x00000108

# write-only options

const SCTP_SET_PEER_PRIMARY_ADDR    = 0x00000006
const SCTP_AUTH_CHUNK    = 0x00000012
const SCTP_AUTH_KEY    = 0x00000013
const SCTP_AUTH_DEACTIVATE_KEY    = 0x0000001d
const SCTP_AUTH_DELETE_KEY    = 0x00000016
const SCTP_RESET_STREAMS    = 0x00000901 #  struct sctp_reset_streams
const SCTP_RESET_ASSOC    = 0x00000902 #  sctp_assoc_t
const SCTP_ADD_STREAMS    = 0x00000903 #  struct sctp_add_streams


text:
# Autogenerated wrapper script for SCS_jll for armv6l-linux-musleabihf
export libscsdir, libscsindir

using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("SCS")
JLLWrappers.@declare_library_product(libscsdir, "libscsdir.so")
JLLWrappers.@declare_library_product(libscsindir, "libscsindir.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libscsdir,
        "lib/libscsdir.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscsindir,
        "lib/libscsindir.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    function gheader(gt)
Retrieve the header infomation of a packed genotype file.
"""
function gheader(gt)
    nlc, nid, dms, gt_majored = nothing, nothing, nothing, true
    open(gt, "r") do io
        header = mmap(io, Vector{Int64}, 3)
        gt_majored = (header[1] == 1)
        nlc, nid = gt_majored ? header[2:3] : header[[3, 2]]
        dms = Tuple(header[2:3])
    end
    gt_majored, nlc, nid, dms
end

"""
    function picd(gt, twop; step = 1000, tol = 1e-4)
Given the binary file of `AGH` genotype format, and a vector of
2 × allele frequencies, this function do the `p`ivoted `i`ncomplete
`C`holesky `d`ecomposition.
Finally it returns an approximate G-inverse.

The program
"""
function picd(gt::String, twop; step = 1000, tol = 1e-4)
    goffset = 24
    gt_majored, nlc, nid, dms = gheader(gt)
    gt_majored || error("The matrix must be genotype majored for speed reasons")
    blk = zeros(dms[1], step)
    # 1/√(2pq), as product is faster than divide.
    rs2pq = 1 ./ sqrt.(twop .* (1 .- twop .* .5))
    open(gt, "r") do io
        gt = mmap(io, Matrix{Int8}, dms, goffset)
        # tp = mean(gt, dims=2) # if gt_majored, or, 1.
        copyto!(blk, view(gt, :, 1:step))
        blk .-= twop
        blk .*= rs2pq
        G = blk'blk
        G, piv, rank, info = LAPACK.pstrf!('L', G, tol)
    end
end

function tst_picd(gt, twop; step = 1000, tol = 1e-4)
    rs2pq = 1 ./ sqrt.(twop .* (1 .- twop .* .5))
    blk = zeros(size(gt)[1], step)
    copyto!(blk, view(gt, :, 1:step))
    blk .-= twop
    blk .*= rs2pq
    G = blk'blk
    G, piv, rank, info = LAPACK.pstrf!('L', G, tol)
end

"""
    function cholesky_decomposition!(A)
This function is just for fun, and should not be taken seriously.
The lower triangle of matrix `A` will be replace by `L`, such that
the original `A = LL'`.
"""
function cholesky_decomposition!(A)
    issymmetric(A) || error("Not a symmetric matrix")
    m = size(A)[1]
    for i in 1:m
        for k in 1:i-1
            A[i, i] -= A[i, k] * A[i, k]
        end
        A[i, i] = sqrt(A[i, i]) # diagonals
        for j in i+1:m
            for k in 1:i-1
                A[j, i] -= A[i, k] * A[j, k]
            end
            A[j, i] /= A[i, i]
        end                     # off diagonals
    end
end


text:
mod10(x) = mod(x, 10)
mod100(x) = mod(x, 100)

function update_s(a::NamedArray, c::Vector{Symbol}, debug=false)
  vars = names(a, 1)
  s = copy(c)
  st = Symbol[]
  while !(s == st)
    st = copy(s)
    for j in s
      for i in vars
        mod10(a[i, j]) == 1 && union!(s, [i])
      end
    end
  end
  debug && println(s)
  s
end

function update_21(ar::NamedArray, marg_vars::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a21 = copy(ar)

  for i in marg_vars
    for j in names(ar, 1)
      a31 = copy(ar)
      if mod10(ar[j, i]) == 1 || mod10(ar[i, j]) == 1
        a31[i, j] = ar[j, i]
        a31[j, i] = ar[i, j]
        debug && println("Update: a31[$i, $j] = $(ar[j, i])")
        debug && println("Update: a31[$j, $i] = $(ar[i, j])")
      end
      idx = vars[findall(x -> mod10(x) == 1, a31[:, i])]
      if length(idx) > 1 && mod10(a31[j, i]) == 1
        for k in idx
          if mod10(ar[k, j]) == 0 && k != j
            debug && println([k, j, 1])
            a21[k, j] = 1
          end
        end
      end
    end
  end
  a21
end

function update_22(ar::NamedArray, s::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a22 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in s
    idx = vars[findall(x -> StructuralCausalModels.mod10(x) == 1, ar[:, i])]
    idy = vars[findall(x -> x > 99, ar[:, i])]
    if length(idx) > 0 && length(idy) >0
      for j in idx
        for k in idy
          if StructuralCausalModels.mod10(ar[j, k]) == 0 && j != k
            debug && println([j, k, 1])
            a22[j, k] = 1
          end
        end
      end
    end
  end
  a22
end

function update_23(ar::NamedArray, m::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a23 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in m
    idx = vars[findall(x -> StructuralCausalModels.mod10(x) == 1, ar[i, :])]
    idy = vars[findall(x -> StructuralCausalModels.mod100(x) > 9, ar[i, :])]
    if length(idx) > 0 && length(idy) > 0
      for j in idx
        for k in idy
          if StructuralCausalModels.mod10(ar[k, j]) == 0 && j !== k
            debug && println([k, j, 1])
            a23[k, j] = 1
          end
        end
      end
    end
  end
  a23
end

function update_24(ar::NamedArray, m::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a24 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in m
    idx = vars[findall(x -> StructuralCausalModels.mod100(x) > 9, ar[:, i])]
    idy = vars[findall(x -> x > 99, ar[:, i])]
    if length(idx) > 0 && length(idy) > 0
      for j in idx
        for k in idy
          if StructuralCausalModels.mod10(ar[j, k]) == 0 && j !== k
            debug && println([j, k, 1])
            a24[j, k] = 1
          end
        end
      end
    end
  end
  a24
end

function update_25(ar::NamedArray, m::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a25 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in m
    idx = vars[findall(x -> StructuralCausalModels.mod10(x) == 1, ar[i, :])]
    lenidx = length(idx)
    if lenidx > 1
      for j in 1:lenidx - 1
        for k in (j + 1):lenidx
          if ar[idx[j], idx[k]] < 100
            debug && println([idx[j], idx[k], 100])
            a25[idx[j], idx[k]] = 100
          end
        end
      end
    end
  end
  a25 .+ transpose(a25)
end

function update_26(ar::NamedArray, m::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a26 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in m
    idx = vars[findall(x -> StructuralCausalModels.mod10(x) == 1, ar[i, :])]
    idy = vars[findall(x -> x > 99, ar[i, :])]
    if length(idx) > 0 && length(idy) >0
      for j in idx
        for k in idy
          if ar[k, j] < 100 && j !== k
            debug && println([k, j, 100])
            a26[k, j] = 100
          end
        end
      end
    end
  end
  a26 .+ transpose(a26)
end

function update_27(ar::NamedArray, s::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a27 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in s
    idx = vars[findall(x -> x > 99, ar[:, i])]
    lenidx = length(idx)
    if lenidx > 1
      for j in 1:lenidx-1
        for k in (j+1):lenidx
          if ar[idx[j], idx[k]] < 100
            debug && println([idx[j], idx[k], 100])
            a27[idx[j], idx[k]] = 100
          end
        end
      end
    end
  end
  a27 .+ transpose(a27)
end

function update_28(ar::NamedArray, s::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a28 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in s
    idx = vars[findall(x -> mod10(x) == 1, ar[:, i])]
    lenidx = length(idx)
    if lenidx > 1
      for j in 1:lenidx-1
        for k in (j+1):lenidx
          if StructuralCausalModels.mod100(ar[idx[j], idx[k]]) < 10
            debug && println([idx[j], idx[k], 10])
            a28[idx[j], idx[k]] = 10
          end
        end
      end
    end
  end
  a28 .+ transpose(a28)
end

function update_29(ar::NamedArray, m::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a29 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in m
    idx = vars[findall(x -> StructuralCausalModels.mod10(x) == 1, ar[:, i])]
    idy = vars[findall(x -> StructuralCausalModels.mod100(x) > 9, ar[:, i])]
    if length(idx) > 0 && length(idy) >0
      for j in idx
        for k in idy
          if mod100(ar[k, j]) < 10 && j !== k
            debug && println([k, j, 10])
            a29[k, j] = 10
          end
        end
      end
    end
  end
  a29 .+ transpose(a29)
end

function update_30(ar::NamedArray, m::Vector{Symbol}, debug=false)

  vars = names(ar, 1)
  a30 = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));

  for i in m
    idx = vars[findall(x -> StructuralCausalModels.mod100(x) > 9, ar[:, i])]
    lenidx = length(idx)
    if lenidx > 1
      for j in 1:lenidx-1
        for k in (j+1):lenidx
          if StructuralCausalModels.mod100(ar[idx[j], idx[k]]) < 10
            debug && println([idx[j], idx[k], 10])
            a30[idx[j], idx[k]] = 10
          end
        end
      end
    end
  end
  a30 .+ transpose(a30)
end

function update_a(ar::NamedArray, debug=false)
  vars = names(ar, 1)
  a = copy(ar)
  for i in vars
    for j in vars
      if mod100(a[i, j]) > 9
        a[i, j] = 10
        for k in vars
          if a[k, j] == 100
            a[j, k] = 1
            debug && println([j, k, 1])
            a[k, j] = 0
            debug && println([k, j, 0])
          end
        end
      end
    end
  end
  a
end

function update_b(ar::NamedArray, s::Vector{Symbol}, debug=false)
  vars = names(ar, 1)
  a = copy(ar)
  s2 = copy(s)
  while true
    st = copy(s2)
    for j in s2
      for k in vars
        if mod10(a[k, j]) == 1
          debug && !(k in s2) && println("$k added to $s2")
          s2 = union!(s2, [k])
        end
      end
    end
    s2 == st && break
  end
  debug && println("\n$s2\n")
  for i in s2
    for j in s2
      if mod10(a[i, j]) == 1
        a[i, j] = 10
        debug && println([i, j, 10])
        a[j, i] = 10
        debug && println([j, i, 10])
      end
    end
  end
  a
end

function update_an(ar::NamedArray, debug=false)
  vars = names(ar, 1)
  an = copy(ar)
  for i in vars
    for j in vars
      a3n = copy(ar)
      if mod10(ar[j, i]) == 1 || mod10(ar[i, j]) == 1
        a3n[i, j] = ar[j, i]
        debug && println("Update: a3n[$i, $j] = $(ar[j, i])")
        a3n[j, i] = ar[i, j]
        debug && println("Update: a3n[$j, $i] = $(ar[i, j])")
      end
      idx = vars[findall(x -> StructuralCausalModels.mod10(x) == 1, a3n[:, i])]
      if length(idx) > 1 && mod10(a3n[j, i]) == 1
        for k in idx
          if mod10(ar[k, j]) == 0 && k != j
            debug && println([k, j, 1])
            an[k, j] = 1
          end
        end
      end
    end
  end
  an
end

function update_27n(ar::NamedArray, an::NamedArray, debug=false)
  vars = names(ar, 1)
  a27n = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));
  for i in vars
    idx = vars[findall(x -> x > 99, an[:, i])]
    lenidx = length(idx)
    if length(idx) > 1
      for j in 1:lenidx-1
        for k in (j + 1):lenidx
          if ar[idx[j], idx[k]] < 100 &&
            mod10(an[i, idx[j]]) == 1 ||
            an[i, idx[k]] == 1

            debug && println([idx[j], idx[k], 100])
            a27n[idx[j], idx[k]] == 100
          end
        end
      end
    end
  end
  a27n
end

function update_22n(ar::NamedArray, an::NamedArray, debug=false)
  vars = names(ar, 1)
  a22n = NamedArray(zeros(Int, size(ar)), (vars, vars), ("Rows", "Cols"));
  for i in vars
    idx = vars[findall(x -> mod10(x) == 1, an[:, i])]
    idy = vars[findall(x -> x > 99, an[:, i])]
    if length(idx) > 0 && length(idy) > 0
      for j in idx
        for k in idy
          if mod10(ar[j, k]) == 0 &&
            j != k &&
            mod10(an[i, k]) == 1

            a22n[j, k] = 1
            debug &&println([j, k, 1])
          end
        end
      end
    end
  end
  a22n
end

function update_i(ar::NamedArray, debug=false)
  vars = names(ar, 1)
  a = copy(ar)
  for i in vars
    for j in vars
      if a[i, j] == 101
        a[i, j] = 1
        debug && println([i, j, 1])
        a[j, i] = 0
        debug && println([j, i, 0])
      end
    end
  end
  a
end


text:
@testset "Tuple" begin

s = """
    (a,b) = (1,2,3)
"""
msgs = lintstr(s)
@test msgs[1].code == :W546
@test contains(msgs[1].message, "implicitly discarding values, 2 of 3 used")

s = """
    a, = (1,2,3)
"""
msgs = lintstr(s)
@test msgs[1].code == :W546
@test contains(msgs[1].message, "implicitly discarding values, 1 of 3 used")

@test isempty(lintstr("a = (1, 2, 3)"))

s = """
function f()
    (a,b,c) = (1,2,3)
    return (b,c,a)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
    (a,b,c) = (1,2)
"""
msgs = lintstr(s)
@test msgs[1].code == :E418
@test contains(msgs[1].message, "RHS is a tuple, 3 of 2 variables used")

end


text:
@testset "Recruitment" begin

@testset "function Constructor" begin
    recruitment = Recruitment( "Recruitment" )
    @test all( [recruitment.name == "Recruitment", recruitment.freq == 1.0,
        recruitment.offset == 0.0, recruitment.targetNode == "dummy",
        recruitment.minRecruitment == 0, recruitment.maxRecruitment == -1,
        !recruitment.isAdaptive,
        recruitment.recruitmentDistType === :disc,
        recruitment.recruitmentDistNodes == Dict( 0 => 1.0 ),
        recruitment.ageDistType === :disc,
        recruitment.ageDistNodes == Dict( 0 => 1.0 )] )
end  # @testset "function Constructor"

recruitment = Recruitment( "Recruitment" )

@testset "function setRecruitmentSchedule!" begin
    setRecruitmentSchedule!( recruitment, 12 )
    @test ( recruitment.freq == 12.0 ) && ( recruitment.offset == 0.0 )
    setRecruitmentSchedule!( recruitment, 6, 8 )
    @test ( recruitment.freq == 6.0 ) && ( recruitment.offset == 2.0 )
    setRecruitmentSchedule!( recruitment, 12, -9 )
    @test ( recruitment.freq == 12.0 ) && ( recruitment.offset == 3.0 )
    @test !setRecruitmentSchedule!( recruitment, -12, 6 )
end  # @testset "function setRecruitmentSchedule!"

@testset "function setRecruitmentTarget!" begin
    setRecruitmentTarget!( recruitment, "Officer" )
    @test recruitment.targetNode == "Officer"
end  # @testset "function setRecruitmentTarget!"

@testset "function setRecruitmentAdaptiveRange!" begin
    setRecruitmentAdaptiveRange!( recruitment, 20 )
    @test all( [recruitment.isAdaptive, recruitment.minRecruitment == 20,
        recruitment.maxRecruitment == -1] )
    setRecruitmentAdaptiveRange!( recruitment, -20, -20 )
    @test ( recruitment.minRecruitment == 0 ) &&
        ( recruitment.maxRecruitment == -1 )
    setRecruitmentAdaptiveRange!( recruitment, 50, 100 )
    @test ( recruitment.minRecruitment == 50 ) &&
        ( recruitment.maxRecruitment == 100 )
    @test !setRecruitmentAdaptiveRange!( recruitment, 50, 30 )
    @test ( recruitment.minRecruitment == 50 ) &&
        ( recruitment.maxRecruitment == 100 )
end  # @testset "function setRecruitmentAdaptiveRange!"

@testset "function setRecruitmentFixed!" begin
    setRecruitmentFixed!( recruitment, 50 )
    @test all( [!recruitment.isAdaptive,
        recruitment.recruitmentDistType === :disc,
        recruitment.recruitmentDistNodes == Dict( 50 => 1.0 )] )
    setRecruitmentFixed!( recruitment, 200 )
    @test recruitment.recruitmentDistNodes == Dict( 200 => 1.0 )
    @test !setRecruitmentFixed!( recruitment, -500 )
    @test recruitment.recruitmentDistNodes == Dict( 200 => 1.0 )
end  # @testset "function setRecruitmentFixed!"

@testset "function setRecruitmentDist!" begin
    setRecruitmentDist!( recruitment, :disc,
        Dict( 50 => 1.0, 75 => 2.0, 100 => 1.0 ) )
    @test ( recruitment.recruitmentDistType === :disc ) &&
        ( recruitment.recruitmentDistNodes ==
            Dict( 50 => 1.0, 75 => 2.0, 100 => 1.0 ) )
    @test all( map( ii -> recruitment.recruitmentDist() ∈ [50, 75, 100],
        1:10000 ) )
    setRecruitmentDist!( recruitment, :pUnif,
        Dict( 50 => 1.0, 70 => 2.0, 90 => 1.0 ) )
    @test ( recruitment.recruitmentDistType === :pUnif ) &&
        ( recruitment.recruitmentDistNodes ==
            Dict( 50 => 1.0, 70 => 2.0, 90 => 1.0 ) )
    @test all( map( ii -> recruitment.recruitmentDist() ∈ 50:89,
        1:10000 ) )
    setRecruitmentDist!( recruitment, :pLin,
        Dict( 50 => 1.0, 75 => 2.0, 100 => 1.0 ) )
    @test ( recruitment.recruitmentDistType === :pLin ) &&
        ( recruitment.recruitmentDistNodes ==
        Dict( 50 => 1.0, 75 => 2.0, 100 => 1.0 ) )
    @test all( map( ii -> recruitment.recruitmentDist() ∈ 50:100,
        1:10000 ) )
    @test !setRecruitmentDist!( recruitment, :test,
        Dict( 50 => 1.0, 75 => 2.0, 100 => 1.0 ) )
    @test !setRecruitmentDist!( recruitment, :disc,
        Dict( 50 => -1.0, -20 => 5.0 ) )
    @test !setRecruitmentDist!( recruitment, :disc,
        Dict( 50 => 0.0, 100 => 0.0 ) )
    @test !setRecruitmentDist!( recruitment, :pUnif,
        Dict( 50 => 0.0, 100 => 5.0 ) )
end  # @testset "function setRecruitmentDist!"

@testset "function setRecruitmentAgeFixed!" begin
    setRecruitmentAgeFixed!( recruitment, 60 )
    @test ( recruitment.ageDistType === :disc ) &&
        ( recruitment.ageDistNodes == Dict( 60.0 => 1 ) )
    @test !setRecruitmentAgeFixed!( recruitment, -60 )
end  # @testset "function setRecruitmentAgeFixed!"

@testset "function setRecruitmentAgeDist!" begin
    setRecruitmentAgeDist!( recruitment, :disc,
        Dict( 20.0 => 1.0, 30.0 => 2.0, 45.0 => 1.0 ) )
    @test ( recruitment.ageDistType === :disc ) &&
        ( recruitment.ageDistNodes ==
            Dict( 20.0 => 1.0, 30.0 => 2.0, 45.0 => 1.0 ) )
    @test all( map( age -> age ∈ [20.0, 30.0, 45.0],
        recruitment.ageDist( 10000 ) ) )
    setRecruitmentAgeDist!( recruitment, :pUnif,
        Dict( 20.0 => 1.0, 25.0 => 2.0, 40.0 => 1.0 ) )
    @test ( recruitment.ageDistType === :pUnif ) &&
        ( recruitment.ageDistNodes ==
            Dict( 20.0 => 1.0, 25.0 => 2.0, 40.0 => 1.0 ) )
    @test all( 20.0 .<= recruitment.ageDist( 10000 ) .<= 40.0 )
    setRecruitmentAgeDist!( recruitment, :pLin,
        Dict( 20.0 => 1.0, 30.0 => 2.0, 45.0 => 1.0 ) )
    @test ( recruitment.ageDistType === :pLin ) &&
        ( recruitment.ageDistNodes ==
            Dict( 20.0 => 1.0, 30.0 => 2.0, 45.0 => 1.0 ) )
    @test all( 20.0 .<= recruitment.ageDist( 10000 ) .<= 45.0 )
    @test !setRecruitmentAgeDist!( recruitment, :test,
        Dict( 20.0 => 1.0, 30.0 => 2.0, 50.0 => 1.0 ) )
    @test !setRecruitmentAgeDist!( recruitment, :disc,
        Dict( 20.0 => -1.0, -20.0 => 5.0 ) )
    @test !setRecruitmentAgeDist!( recruitment, :disc,
        Dict( 20.0 => 0.0, 50.0 => 0.0 ) )
    @test !setRecruitmentAgeDist!( recruitment, :pUnif,
        Dict( 20.0 => 0.0, 50.0 => 5.0 ) )
end  # @testset "function setRecruitmentAgeDist!"

end  # @testset "Recruitment"


text:
"""
    BitTorrent

Julia package for working with bittorrent technology. See
https://fredrikekre.github.io/BitTorrent.jl/ for documentation.
"""
module BitTorrent

import SHA

include("bencode.jl")
include("torrent-files.jl")

end # module


text:
@inline trans(b::NNlib.BatchedAdjOrTrans) = (b isa NNlib.BatchedTranspose ? static(Int('T')) : static(Int('C'))), parent(b)
@inline trans(x) = static(Int('N')), x
@inline trans(c, x) = Int(c) == static(Int('T')) ? batched_transpose(x) :
    Int(c) == static(Int('C')) ? batched_adjoint(x) : x

matmul(a, b) = matmul(a, b, true)
matmul(a::AbstractVecOrMat, b::AbstractVecOrMat, s::Number) = (a * b) .* s
function matmul(a::AbstractArray, b::AbstractArray, s::Number)
    transA, pA = trans(a)
    transB, pB = trans(b)
    A = CollapsedDimArray(pA)
    B = CollapsedDimArray(pB)
    return matmul_wrapper(transA, transB, s, A, B)
end

@inline gemm_strided_batched_wrapper(transA::Union{AbstractChar, StaticInt}, transB::Union{AbstractChar, StaticInt}, alpha::Number, A::AbstractArray, B::AbstractArray) =
    gemm_strided_batched_wrapper(transA, transB, alpha, CollapsedDimArray(A), CollapsedDimArray(B))

@inline gemm_strided_batched_wrapper(transA::Union{AbstractChar, StaticInt}, transB::Union{AbstractChar, StaticInt}, alpha::Number, A::AbstractArray, B::CollapsedDimArray) =
    gemm_strided_batched_wrapper(transA, transB, alpha, CollapsedDimArray(A), B)

@inline gemm_strided_batched_wrapper(transA::Union{AbstractChar, StaticInt}, transB::Union{AbstractChar, StaticInt}, alpha::Number, A::CollapsedDimArray, B::AbstractArray) =
    gemm_strided_batched_wrapper(transA, transB, alpha, A, CollapsedDimArray(B))

@inline function gemm_strided_batched_wrapper(transA::Union{AbstractChar, StaticInt}, transB::Union{AbstractChar, StaticInt}, alpha::Number, A::CollapsedDimArray, B::CollapsedDimArray)
    m = noncollapsed_size(A.parent, A.si, A.sj, Int(transA) == static(Int('N')) ? static(1) : static(2))
    n = noncollapsed_size(B.parent, B.si, B.sj, Int(transB) == static(Int('N')) ? static(2) : static(1))
    # batch size differ is allow only when ones batch size is one
    sc3 = isonebatch(B) ?
        noncollapsed_size(A.parent, A.si, A.sj, static(3)) :
        noncollapsed_size(B.parent, B.si, B.sj, static(3))

    T = promote_type(eltype(A), eltype(B))
    if eltype(A) == T
        pA = parent(A)
    else
        pA = convert(AbstractArray{T}, parent(A))
    end
    if eltype(B) == T
        pB = parent(B)
    else
        pB = convert(AbstractArray{T}, parent(B))
    end

    Ci = static(length(m) + 1)
    Cj = static(Ci + length(n))
    C = similar(pB, T, (m..., n..., sc3...))
    if ndims(pA) == ndims(pB) == ndims(C) == 3
        gemm_strided_batched!(as_char(transA), as_char(transB), convert(T, alpha), pA, pB, zero(T), C)
    else
        gemm_strided_batched!(as_char(transA), as_char(transB), convert(T, alpha), pA, pB, zero(T), C, A.si, A.sj, B.si, B.sj, Ci, Cj)
    end

    return CollapsedDimArray(C, Ci, Cj, A.onebatch & B.onebatch)
end

function generic_matmul(transA::Union{AbstractChar, StaticInt}, transB::Union{AbstractChar, StaticInt}, alpha::Number, A, B)
    T = promote_type(eltype(A), eltype(B))
    scale = convert(T, alpha)
    if A isa CollapsedDimArray
        m = noncollapsed_size(A.parent, A.si, A.sj, Int(transA) == static(Int('N')) ? static(1) : static(2))
        pA = collapseddim(A)
        sa3 = noncollapsed_size(A.parent, A.si, A.sj, static(3))
    else
        m = size(A, Int(transA) == static(Int('N')) ? static(1) : static(2))
        pA = A
        sa3 = size(A, 3)
    end

    if B isa CollapsedDimArray
        n = noncollapsed_size(B.parent, B.si, B.sj, Int(transB) == static(Int('N')) ? static(2) : static(1))
        pB = collapseddim(B)
        sb3 = noncollapsed_size(B.parent, B.si, B.sj, static(3))
    else
        n = size(B, Int(transB) == static(Int('N')) ? static(2) : static(1))
        pB = B
        sb3 = size(B, 3)
    end
    sc3 = isonebatch(B) ? sa3 : sb3
    Ci = static(length(m) + 1)
    Cj = static(Ci + length(n))
    outsize = (m..., n..., sc3...)
    y = scale .* batched_mul(trans(transA, pA), trans(transB, pB))
    return CollapsedDimArray(reshape(y, outsize), Ci, Cj)
end

NNlib.is_strided(ca::CollapsedDimArray) = NNlib.is_strided(parent(ca))

@inline function matmul_wrapper(transA::Union{AbstractChar, StaticInt}, transB::Union{AbstractChar, StaticInt}, alpha::Number, A::AbstractArray{TA, 3}, B::AbstractArray{TB, 3}) where {TA, TB}
    mA = size(A, Int(transA) == static(Int('N')) ? 1 : 2)
    kA = size(A, Int(transA) == static(Int('N')) ? 2 : 1)
    bA = size(A, 3)
    kB = size(B, Int(transB) == static(Int('N')) ? 1 : 2)
    nB = size(B, Int(transB) == static(Int('N')) ? 2 : 1)
    bB = size(B, 3)

    if kA != kB || (bA != bB && bA != 1 && bB != 1)
        throw(DimensionMismatch("A has dimensions ($mA,$kA,$bA) but B has dimensions ($kB,$nB,$bB)"))
    end

    if TA <: BLAS.BlasFloat && TB <: BLAS.BlasFloat && NNlib.is_strided(A) && NNlib.is_strided(B)
        return gemm_strided_batched_wrapper(transA, transB, alpha, A, B)
    else
        return generic_matmul(transA, transB, alpha, A, B)
    end
end