text:
__precompile__()
"""
# Public API (nothing is exported)

* lookupname(str)
* matchchar(char)
* matches(str)
* longestmatches(str)
* completions(str)
"""
module Emoji_Entities

using StrTables

VER = UInt32(1)

struct Emoji_Table{T} <: AbstractEntityTable
    ver::UInt32
    tim::String
    inf::String
    base32::UInt32
    base2c::UInt32
    nam::StrTable{T}
    ind::Vector{UInt16}
    val16::Vector{UInt16}
    ind16::Vector{UInt16}
    val32::Vector{UInt16}
    ind32::Vector{UInt16}
    val2c::StrTable{T}
    ind2c::Vector{UInt16}
    max2c::UInt32
end

function __init__()
    global default =
        Emoji_Table(StrTables.load(joinpath(@__DIR__, "../data", "emoji.dat"))...)
    nothing
end

StrTables._get_val2c(tab::Emoji_Table, val) = val

function StrTables.matches(tab::Emoji_Table, vec::String)
    (isempty(vec)
     ? StrTables._empty_str_vec
     : (length(vec) == 1
        ? matchchar(tab, vec[1])
        : StrTables._get_strings(tab, vec, tab.val2c, tab.ind2c)))
end

StrTables.matches(tab::Emoji_Table, str::AbstractString) = matches(tab, String(str))

function StrTables.longestmatches(tab::Emoji_Table, vec::Vector{T}) where {T}
    isempty(vec) && return StrTables._empty_str_vec
    ch = vec[1]
    len = length(vec)
    len == 1 && return matchchar(tab, ch)
    # Get range that matches the first character, if any
    rng = StrTables.matchfirstrng(tab.val2c, string(ch))
    if !isempty(rng)
        maxlen = min(len, tab.max2c)
        # Truncate vec
        vec = vec[1:maxlen]
        # Need to find longest matching strings
        for l = 2:maxlen
            length(rng) == 1 && break
            prevrng = rng
            rng = StrTables.matchfirstrng(tab.val2c, string(vec[1:l]))
            isempty(rng) && (rng = prevrng; break)
        end
        return tab.nam[tab.ind2c[rng]]
    end
    # Fall through and check only the first character
    matchchar(tab, ch)
end

end # module


text:
using DataStructures, StringDistances

x = map(randstring, rand(5:25,100_000))
y = map(randstring, rand(5:25,100_000))
function f(out, t, x, y)
    d = Array(out, length(x))
    @inbounds for i in 1:length(x)
        d[i] = evaluate(t, x[i], y[i])
    end
end

# similar
@time f(Int, Levenshtein(), x, y)
@time f(Float64, Jaro(), x, y)

# 2x slower compared to StringDist
@time f(Int, QGram(2), x, y)
@time f(Float64, Cosine(2), x, y)
@time f(Float64, Jaccard(2), x, y)

#
@time f(Float64, RatcliffObershelp(), x, y)

#= Rcode
library(stringdist)
x <- sapply(sample(5:25,1e5,replace=TRUE), function(n) paste(sample(letters,n,replace=TRUE),collapse=""))
y <- sapply(sample(5:25,1e5,replace=TRUE), function(n) paste(sample(letters,n,replace=TRUE),collapse=""))
system.time(stringdist(x,y,method='lv', nthread = 1))
system.time(stringdist(x,y,method='jaccard', nthread = 1))
system.time(stringdist(x,y,method='cosine', nthread = 1))
system.time(stringdist(x,y,method='qgram', nthread = 1))

=#

ext: jl


text:
# Model.
function bayesian_linear_regression(x::Vector{Float64})
    σ = rand(:σ, InverseGamma(2, 3))
    β = rand(:β, Normal(0.0, 1.0))
    for i in 1 : length(x)
        rand(:y => i, Normal(β * x[i], σ))
    end
end

# Data.
data_len = 100
x = [Float64(i) for i in 1 : data_len]
y = map(x) do k
    3.0 * k + randn()
end
obs = target(map(1 : data_len) do i
                 (:y => i, ) => y[i]
             end)

@testset "Inference library smoke test 1 - Bayesian linear regression" begin

    # Importance sampling.
    is_test = () -> begin
        n_samples = 50000
        ps, nw = importance_sampling(obs, n_samples, bayesian_linear_regression, (x, ))
        zipped = zip(ps.calls, nw)

        est_σ = sum(map(zipped) do (cl, w)
                        (cl[:σ]) * w
                    end)

        est_β = sum(map(zipped) do (cl, w)
                        (cl[:β]) * w
                    end)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # MH random walk kernel.
    mh_test = () -> begin
        n_iters = 50000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = mh(target([(:σ, ), (:β, )]), cl)
            i % 50 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=2e-1)
    end

    # MH random walk kernel with proposal.
    proposal = (cl, y) -> begin
        σ = rand(:σ, Normal(6.0, 1.0))
        β = rand(:β, Normal(mean(y) / length(y), 3.0))
    end

    mh_test_with_proposal = () -> begin
        n_iters = 10000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = mh(cl, proposal, (y, ))
            i % 30 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # ADVI.
    surrogate = () -> begin
        μ₁ = learnable(:μ₁)
        μ₂ = learnable(:μ₂)
        σ₂ = learnable(:σ₂)
        σ = rand(:σ, Normal(1.0 + μ₁, 0.3))
        β = rand(:β, Normal(μ₂, σ₂ ^ 2))
    end

    advi_test = () -> begin
        ps = learnables([(:μ₁, ) => 5.0,
                         (:μ₂, ) => 1.0,
                         (:σ₂, ) => 1.0])
        @time ps, elbows, _ = advi(obs, ps,
                                   surrogate, (),
                                   bayesian_linear_regression, (x, );
                                   opt = ADAM(0.05, (0.9, 0.8)),
                                   iters = 1000)
    end

    # ADGV.
    surrogate = () -> begin
        μ₁ = learnable(:μ₁)
        μ₂ = learnable(:μ₂)
        σ₂ = learnable(:σ₂)
        σ = rand(:σ, Normal(1.0 + μ₁, 0.3))
        β = rand(:β, Normal(μ₂, σ₂ ^ 2))
    end

    adgv_test = () -> begin
        ps = learnables([(:μ₁, ) => 5.0,
                         (:μ₂, ) => 1.0,
                         (:σ₂, ) => 1.0])
        ps, elbows, _ = adgv(obs, ps, 10,
                             surrogate, (),
                             bayesian_linear_regression, (x, );
                             opt = ADAM(0.05, (0.9, 0.8)),
                             iters = 1000)
    end

    # HMC kernel.
    hmc_test = () -> begin
        n_iters = 1000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = hmc(target([(:σ, ), (:β, )]), cl)
            i % 10 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=2e-1)
    end

    # Combination kernel.
    combo_kernel_test = () -> begin
        n_iters = 1000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = mh(cl, proposal, (y, ))
            cl, _ = hmc(target([(:σ, ), (:β, )]), cl)
            i % 10 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # Boomerang kernel.
    boomerang_test = () -> begin
        n_iters = 500
        ret, cl = generate(obs, bayesian_linear_regression, x)
        sel = get_target(cl)
        sel_values, _ = get_choice_gradients(sel, cl, 1.0)
        d = length(array(sel_values, Float64))
        flow = Boomerang(sparse(I, d, d), zeros(d), 1.0)
        θ = rand(MvNormal(d, 1.0))
        calls = []
        for i in 1 : n_iters
            cl, _ = pdmk(target([(:σ, ), (:β, )]), cl, flow, θ)
            i % 10 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end)
    end

    is_test()
    mh_test()
    mh_test_with_proposal()
    hmc_test()
    combo_kernel_test()
    #@time boomerang_test()
    #@time advi_test()
    #@time adgv_test()

end


text:
# This script solves a Fredholm equation of the first type

println('\n', " "^4, "> Loading the packages...")

using LinearAlgebra # Norm
using Optim # Optimization
using Plots # Plotting
using Printf # Formatted printing
using QuadGK # Gauss–Kronrod integration

# Use the GR backend for plots
gr()

# Change the default font for plots
default(fontfamily="Computer Modern", dpi=300, legend=:outerright)

println(" "^4, "> Computing the solution (Test)...")

# Define the kernel (Test)
K(x, s) = exp(s * x)

# Define the right part of the equation (Test)
u(x) = (exp(x + 1) - 1) / (x + 1)

# Set the integration intervals
a = 0
b = 1
c = 0
d = 1

# Set the number of nodes
n = 100

# Calculate the step
h = (b - a) / n

# Calculate the nodes for the s argument
s = [ a + 0.5 * h + (i - 1) * h for i in 1:n ]

# Set the initial value of the regularization parameter
α = 0.001

# Prepare a range of nodes for the residual calculation
xr = range(c, d; length=1000)

# Compute the residual of the solution with the specified number of nodes and the regularization parameter
function calculate(n::Int, α::Float64)::Tuple{Vector{Float64},Float64}
    # Calculate the step
    h = (b - a) / n

    # Calculate the nodes for the s argument
    s = [ a + 0.5 * h + (i - 1) * h for i in 1:n ]

    # Calculate the nodes for the t argument
    t = copy(s)

    # Compute the g vector
    g = [ quadgk(x -> K(x, s[i]) * u(x), c, d; rtol=1e-8)[1] for i in 1:n ]

    # Prepare a matrix for the computation
    Bα = Matrix{Float64}(undef, n, n)

    # Compute the Bα matrix
    for i in 1:n, j in 1:n
        Bα[i, j] = quadgk(x -> K(x, s[i]) * K(x, t[j]), c, d; rtol=1e-8)[1] * h + (i == j ? α : 0)
    end

    # Compute the solution
    z = Symmetric(Bα) \ g

    # Calculate the residual
    r = norm([ sum(K.(x, s) .* z .* h) for x in xr ] .- u.(xr))

    return z, r
end

# Compute the residual of the solution with the specified regularization parameter (optimization only)
function residual(θ::Vector{Float64})::Float64
    return calculate(100, θ[1])[2]
end

# Optimize over the regularization parameter
res = Optim.optimize(
    residual,
    [α,],
    LBFGS(),
    Optim.Options(
        show_trace=false,
        extended_trace=true,
        store_trace=true,
    );
    inplace=false,
)

# Save the trace and the results
open("test_optimization.trace", "w") do io
    println(io, res.trace)
    println(
        io,
        " * Parameters:\n",
        " "^4, "α = $(res.minimum)\n",
    )
    show(io, res)
end

# Print the residual
println(
    '\n',
    " "^6, "Regularization parameter: ", res.minimizer[1], '\n',
    " "^6, "Residual: ", res.minimum, '\n',
)

println(" "^4, "> Plotting the solution (Test)...")

# Recalculate the solution
z = calculate(n, res.minimizer[1])

# Plot the solution
p = plot(s, z; label="Приближенное решение", xlabel="s", ylabel="z(s)");

# Add the solution to the plot
plot!(p, x -> exp(x), a, b; label="Точное решение");

# Save the figure
savefig(p, "$(@__DIR__)/../plots/test.pdf")

println('\n', " "^6, "* The figure `test.pdf` is saved. *", '\n')

println(" "^4, "> Plotting the residuals against the number of nodes (Test)...")

# Define the `n`s to plot against
ns = [ i for i in 10:300 ]

# Calculate the residuals for these `n`s
rs = @. getindex(calculate(ns, res.minimizer[1]), 2)

# Make a plot of the residual against `n`
p = plot(ns, rs; label="", xlabel="n", ylabel="r");

# Save the figure
savefig(p, "$(@__DIR__)/../plots/test_residual.pdf")

println('\n', " "^6, "* The figure `test_residual.pdf` is saved. *", '\n')

println(" "^4, "> Computing the solution (Assignment)...")

# Define the kernel
K(x, s) = 1 / (1 + x + s)

# Define the right part of the equation
u(x) = 1 / sqrt(2 - x) * (log((2 + x) / (1 + x)) +
       2 * log((sqrt(3) + sqrt(2 - x)) / (2 + sqrt(2 - x))))

# Compute the g vector
g = [ quadgk(x -> K(x, s[i]) * u(x), c, d; rtol=1e-8)[1] for i in 1:n ]

# Optimize over the regularization parameter
res = Optim.optimize(
    residual,
    [α,],
    LBFGS(),
    Optim.Options(
        show_trace=false,
        extended_trace=true,
        store_trace=true,
    );
    inplace=false,
)

# Save the trace and the results
open("assignment_optimization.trace", "w") do io
    println(io, res.trace)
    println(
        io,
        " * Parameters:\n",
        " "^4, "α = $(res.minimum)\n",
    )
    show(io, res)
end

# Print the residual
println(
    '\n',
    " "^6, "Regularization parameter: ", res.minimizer[1], '\n',
    " "^6, "Residual: ", res.minimum, '\n',
)

println(" "^4, "> Plotting the solution (Assignment)...")

# Recalculate the solution
z = calculate(n, res.minimizer[1])

# Plot the solution
p = plot(s, z; label="", xlabel="s", ylabel="z(s)");

# Save the figure
savefig(p, "$(@__DIR__)/../plots/assignment.pdf")

println('\n', " "^6, "* The figure `assignment.pdf` is saved. *", '\n')

println(" "^4, "> Plotting the residuals against the number of nodes (Assignment)...")

# Define the `n`s to plot against
ns = [ i for i in 10:300 ]

# Calculate the residuals for these `n`s
rs = @. getindex(calculate(ns, res.minimizer[1]), 2)

# Make a plot of the residual against `n`
p = plot(ns, rs; label="", xlabel="n", ylabel="r");

# Save the figure
savefig(p, "$(@__DIR__)/../plots/assignment_residual.pdf")

println('\n', " "^6, "* The figure `assignment_residual.pdf` is saved. *", '\n')

println(" "^4, "> Calculating the solutions for tables (Assignment)...")

# Define the header
header = "α = 1e-4" * " "^7 * "1e-5" * " "^11 * "1e-6" * " "^11 * "1e-7" *
         " "^11 * "1e-8" * " "^11 * "1e-9" * " "^11 * "1e-10" * " "^10 * "1e-11" *
         " "^10 * "1e-12" * " "^10 * "1e-13" * " "^10 * "1e-14" * " "^10 * "1e-15"

# Calculate the solutions for different numbers of nodes and regularization parameters
for n in [100, 200, 300]
    step = Int(n * 0.05)
    data = Matrix{Float64}(undef, 20, 12)
    format = Printf.Format("%.8e "^11 * "%.8e\n")
    open("$(@__DIR__)/../tables/$(n).dat", "w") do io
        println(io, header)
        for (i, α) in enumerate([1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15])
            sol, _ = calculate(n, α)
            data[:, i] = sol[1:step:end]
        end
        for i in 1:20
            Printf.format(io, format, data[i, :]...)
        end
    end
end

println('\n', " "^6, "* Tables `100.dat`, `200.dat`, and `300.dat` are saved. *", '\n')


text:
# Interconnections of linear systems

# include("StateSpace.jl");
# include("Utilities.jl");

function Negate(S::StateSpace)::StateSpace
    return StateSpace(S.A,S.B,-S.C,-S.D);
end

function Subsystem(S::StateSpace, inputs::Vector{Int64}, outputs::Vector{Int64})::StateSpace
    # Check IO dimensions.
    if (!isempty(filter(x->(x>S.nu),inputs))) error("Input index exceeds input dimension."); end
    if (!isempty(filter(x->(x>S.ny),outputs))) error("Output index exceeds input dimension."); end

    B = S.B[:,inputs];
    C = S.C[outputs,:];
    D = S.D[outputs,inputs];

    return StateSpace(S.A,B,C,D);
end

function Parallel(S1::StateSpace,S2::StateSpace)::StateSpace
    # Check io dimensions.
    if (S1.nu!=S2.nu || S1.ny!=S2.ny)
        error("Input and output dimensions must match.")
    end

    A = [S1.A zeros(S1.ns,S2.ns);zeros(S2.ns,S1.ns) S2.A];
    B = [S1.B;S2.B];
    C = [S1.C S2.C];
    D = (S1.D+S2.D);
    return StateSpace(A,B,C,D);
end

function Series(S1::StateSpace,S2::StateSpace)::StateSpace
    # Check io dimensions.
    if (S1.ny!=S2.nu)
        error("Output of S1 must match input of S2.");
    end

    A = [S1.A zeros(S1.ns,S2.ns);S2.B*S1.C S2.A];
    B = [S1.B;S2.B*S1.D];
    C = [S2.D*S1.C S2.C];
    D = S2.D*S1.D;

    return StateSpace(A,B,C,D);
end

function Feedback(S1::StateSpace,S2::StateSpace)::StateSpace
    # Check io dimensions -- or provide io options for feedback.
    if (S1.ny != S2.nu || S1.nu != S2.ny)
        error("Incompatible I/O dimensions.")
    end

    A1 = S1.A; B1 = S1.B; C1 = S1.C; D1 = S1.D;
    A2 = S2.A; B2 = S2.B; C2 = S2.C; D2 = S2.D;

    I = eye(S2.ny);
    X = I+D2*D1;

    if (det(X) == 0)
        error("Feedback interconnection is illposed.");
    end

    invX = inv(X);
    Cb = [D2*C1 C2];

    AA1 = [A1 zeros(S1.ns,S2.ns)] - B1*inv(X)Cb;
    AA2 = [B2*C1 A2];
    AA = [AA1;AA2];

    BB1 = B1*invX;
    BB2 = B2*D1;
    BB = [BB1;BB2];

    CC = [C1 zeros(S1.ny,S2.ns)] - D1*invX*Cb;
    DD = D1*inv(X);

    return StateSpace(AA,BB,CC,DD);
end

# Need to implement partial interconnections.
# Eventually, these will replace the above code.

function Parallel(S1::StateSpace,S2::StateSpace,
    in1::Vector{Int64}, in2::Vector{Int64},
    out1::Vector{Int64}, out2::Vector{Int64})::StateSpace
    #                               +------+
    #                 v1 ---------->|      |----------> z1
    #                               |  S1  |
    #                        u1 +-->|      |---+ y1
    #                           |   +------+   |
    #                  u ------>+              O------> y
    #                           |   +------+   |
    #                        u2 +-->|      |---+ y2
    #                               |  S2  |
    #                 v2 ---------->|      |----------> z2
    #                               +------+
    # New system maps (v1;u;v2) -> (z1;y;z2)
    #   -- in1 defines channel u1
    #   -- in2 defines channel u2
    #   -- out1 defines channel y1
    #   -- out2 defines channel y2

    # Check input/output indices are within signal dimensions
    if (!isempty(filter(x->(x>S1.nu),in1))) error("S1: Input index exceeds input dimension."); end
    if (!isempty(filter(x->(x>S1.ny),out1))) error("S1: Output index exceeds output dimension."); end
    if (!isempty(filter(x->(x>S2.nu),in2))) error("S2: Input index exceeds input dimension."); end
    if (!isempty(filter(x->(x>S2.ny),out2))) error("S2: Output index exceeds output dimension."); end

    # Check compatibility of interconnection
    if (length(in1)!= length(in2)) error("Size of in1 and in2 must be same."); end
    if (length(out1)!= length(out2)) error("Size of out1 and out2 must be same."); end

    # Get indices for partitions of system matrices
    in1c = setdiff(1:S1.nu,in1); out1c = setdiff(1:S1.ny,out1); # These can be empty
    in2c = setdiff(1:S2.nu,in2); out2c = setdiff(1:S2.ny,out2); # These can be empty

    nv1 = length(in1c);
    nv2 = length(in2c);
    nz1 = length(out1c);
    nz2 = length(out2c);
    ns1 = S1.ns;
    ns2 = S2.ns;

    # Partitions of S1
    B1 = S1.B[:,in1c]; # defines v1 channel
    B2 = S1.B[:,in1]; # defines u1 channel
    C1 = S1.C[out1c,:];
    C2 = S1.C[out1,:];
    D11 = S1.D[out1c,in1c];
    D12 = S1.D[out1c,in1];
    D21 = S1.D[out1,in1c];
    D22 = S1.D[out1,in1];

    # Partitions of S2
    F1 = S2.B[:,in2]; # defines u2 channel
    F2 = S2.B[:,in2c];  # defines v2 channel
    G1 = S2.C[out2,:];
    G2 = S2.C[out2c,:];
    H11 = S2.D[out2,in2];
    H12 = S2.D[out2,in2c];
    H21 = S2.D[out2c,in2];
    H22 = S2.D[out2c,in2c];

    # Construct the new system matrices
    AA = [S1.A zeros(ns1,ns2);
    zeros(ns2,ns1) S2.A];

    # If nz1, nz2, nv1, nv2 are zero, this will result in empty matrices in the partitions.
    # Creating matrices with empty submatrices is a problem.
    # Using hcat and vcat to create the matrices, eliminates this problem.
    BB1 = hcat(B1,B2,zeros(ns1,nv2));
    BB2 = hcat(zeros(ns2,nv1),F1,F2);
    BB  = vcat(BB1,BB2);

    CC1 = hcat(C1,zeros(nz1,ns2));
    CC2 = hcat(C2,G1);
    CC3 = hcat(zeros(nz2,ns1),G2);
    CC  = vcat(CC1,CC2,CC3);

    DD1 = hcat(D11, D12, zeros(nz1,nv2));
    DD2 = hcat(D21, (D22+H11), H12);
    DD3 = hcat(zeros(nz2,nv1), H21, H22);
    DD = vcat(DD1,DD2,DD3);

    return StateSpace(AA,BB,CC,DD);
end

function Series(S1::StateSpace,S2::StateSpace, out1::Vector{Int64}, in2::Vector{Int64})::StateSpace
    #                                       +------+
    #                                v2 --->|      |
    #                       +------+        |  S2  |-----> y2
    #                       |      |------->|      |
    #              u1 ----->|      |y1   u2 +------+
    #                       |  S1  |
    #                       |      |---> z1
    #                       +------+
    # New system maps (v2;u1) -> (y2;z1)
    # -- out1 defines y1
    # -- in2 defines u2

    # Check input/output indices are within signal dimensions
    if (!isempty(filter(x->(x>S1.ny),out1))) error("S1: Output index exceeds output dimension."); end
    if (!isempty(filter(x->(x>S2.nu),in2)))  error("S2: Input index exceeds input dimension."); end

    # Check compatibility of interconnection
    if (length(out1)!= length(in2)) error("Incompatible I/O dimension for series connection."); end

    # Get indices for partitions of system matrices
    out1c = setdiff(1:S1.ny,out1); # Defines z, These can be empty
    in2c = setdiff(1:S2.nu,in2);   # These can be empty

    nv2 = length(in2c);
    nz1 = length(out1c);

    ns1 = S1.ns;
    ns2 = S2.ns;

    # Partitions of S1 -- follows the diagram
    B = S1.B;
    C1 = S1.C[out1,:]; # x1->y1
    C2 = S1.C[out1c,:]; # x1->z1
    D1 = S1.D[out1,:]; # u->y1
    D2 = S1.D[out1c,:]; # u->z1

    # Partitions of S2
    F1 = S2.B[:,in2c]; # defines v2 channel
    F2 = S2.B[:,in2];  # defines u2 channel
    G = S2.C; # y2
    H1 = S2.D[:,in2c]; # v2->y2
    H2 = S2.D[:,in2];  # u2->y2

    # Construct the new system matrices
    AA = [S1.A zeros(ns1,ns2);
          F2*C1 S2.A];

    # If nz1, nv2 are zero, this will result in empty matrices in the partitions.
    # Creating matrices with empty submatrices is a problem.
    # Using hcat and vcat to create the matrices, eliminates this problem.
    BB1 = hcat(zeros(ns1,nv2),B);
    BB2 = hcat(F1,F2*D1);
    BB  = vcat(BB1,BB2);

    CC1 = hcat(H2*C1,G);
    CC2 = hcat(C2,zeros(nz1,ns2));
    CC  = vcat(CC1,CC2);

    DD1 = hcat(H1,H2*D1);
    DD2 = hcat(zeros(nz1,nv2),D2);
    DD = vcat(DD1,DD2);

    return StateSpace(AA,BB,CC,DD);
end

function Feedback(S1::StateSpace,S2::StateSpace, feedin::Vector{Int64}, feedout::Vector{Int64},FeedbackSign::Int64=-1)::StateSpace

    #                            +------+
    #                v --------->|      |--------> z
    #                            |  S1  |
    #                u --->O---->|      |----+---> y
    #                      |  u1 +------+    |
    #                      |                 |
    #                      +-----[  S2  ]<---+
    # Maps  (v;u) -> (z;y)
    #   -- feedin defines u1
    #   -- feedout defines y

    # Check input/output indices are within signal dimensions
    if (!isempty(filter(x->(x>S1.nu),feedin))) error("S1: Feedback input index exceeds input dimension."); end
    if (!isempty(filter(x->(x>S1.ny),feedout))) error("S1: Feedback output index exceeds output dimension."); end

    # Check compatibility of interconnection
    if (length(feedin)!= S2.ny || length(feedout)!=S2.nu)
        error("Incompatible I/O dimensions for feedback.");
    end

    # Get indices for partitions of system matrices
    in1c = setdiff(1:S1.nu,feedin); # defines v, can be empty
    out1c = setdiff(1:S1.ny,feedout); # defines z, can be empty

    nv = length(in1c);
    nz = length(out1c);
    ns1 = S1.ns;
    ns2 = S2.ns;
    ny = length(feedout);
    nu = length(feedin);

    # Partitions of S1
    A = S1.A;
    B1 = S1.B[:,in1c]; # defines v channel
    B2 = S1.B[:,feedin]; # defines u1 channel
    C1 = S1.C[out1c,:]; # defines z channel
    C2 = S1.C[feedout,:]; # defines y channel
    D11 = S1.D[out1c,in1c]; # v->z
    D12 = S1.D[out1c,feedin]; # u1->z
    D21 = S1.D[feedout,in1c]; # v->y
    D22 = S1.D[feedout,feedin]; # u1->y

    # System matrices of S2
    E = S2.A;
    F = S2.B;
    G = S2.C;
    H = S2.D;

    I = eye(nu);
    K = (I-FeedbackSign*H*D22);
    if (det(K) == 0)
        error("Feedback interconnection is illposed.");
    end
    invK = inv(K);

    X = FeedbackSign*invK*hcat(H*C2,G);
    Y = invK*hcat(FeedbackSign*H*D21,I);

    # Construct the new system matrices
    AA1 = hcat(A,zeros(ns1,ns2)) + B2*X;
    AA2 = hcat(F*C2,E) + F*D22*X;
    AA = vcat(AA1,AA2);

    BB1 = hcat(B1,zeros(ns1,nu)) + B2*Y;
    BB2 = hcat(F*D21,zeros(ns2,nu)) + F*D22*Y;
    BB = vcat(BB1,BB2);

    CC1 = hcat(C1,zeros(nz,ns2)) + D12*X;
    CC2 = hcat(C2,zeros(ny,ns2)) + D22*X;
    CC = vcat(CC1,CC2);

    DD1 = hcat(D11,zeros(nz,nu)) + D12*Y;
    DD2 = hcat(D21,zeros(ny,nu)) + D22*Y;
    DD = vcat(DD1,DD2);

    return StateSpace(AA,BB,CC,DD);
end

# TO DO
# 1. Linear fractional transformation, and LFT algebra.
# LinearFractionalTransformation(...) -- uppper and lower.

max_stars_repo_name: isrlab/RobustControl

max_issues_repo_name: isrlab/RobustControl

max_forks_repo_name: isrlab/RobustControl


text:
using Stheno: FiniteGP

"""
    check_consistency(rng::AbstractRNG, θ, f, x::AV, y::AV, A, z::AV)

Some basic consistency checks for the function `f(θ)::Tuple{GP, GP}`. Mainly just checks
that Zygote works properly for `f`, and correctly derives the gradients w.r.t. `θ` for
`rand`, `logpdf`, `elbo` when considering the f.d.d.s `f(x, Σ)` and observations `y`, where
`Σ = _to_psd(A)`. The first output of `f` will be the GP sampled from and whose `logpdf`
will be computed, while the second will be used as the process for the pseudo-points, whose
inputs are `z`.
"""
function check_consistency(rng::AbstractRNG, θ, f, x::AV, y::AV, A, z::AV, B)

    # Check input consistency to prevent test failures for the wrong reasons.
    @assert length(x) == length(y)

    g = (θ, x, A)->FiniteGP(first(f(θ)), x, _to_psd(A))
    function h(θ, x, A, z, B)
        v, u = f(θ)
        return FiniteGP(v, x, _to_psd(A)), FiniteGP(u, z, _to_psd(B))
    end

    #
    # mean / cov tests
    #

    # Check that we can differentiate through evaluation of the mean vector.
    adjoint_test(
        (θ, x, A)->mean(g(θ, x, A)),
        randn(rng, length(x)), θ, x, A;
        rtol=1e-4, atol=1e-4,
    )

    # Check that we can differentiate through evaluation of the covariance matrix.
    adjoint_test(
        (θ, x, A)->cov(g(θ, x, A)),
        randn(rng, length(x), length(x)), θ, x, A;
        rtol=1e-4, atol=1e-4,
    )

    adjoint_test(
        (θ, x, A, z, B)->cov(h(θ, x, A, z, B)...),
        randn(rng, length(x), length(z)), θ, x, A, z, B;
        rtol=1e-4, atol=1e-4,
    )

    # #
    # # rand / logpdf / elbo tests
    # #

    # # Check that the gradient w.r.t. the samples is correct (single-sample).
    # adjoint_test(
    #     (θ, x, A)->rand(MersenneTwister(123456), g(θ, x, A)),
    #     randn(rng, length(x)), θ, x, A;
    #     rtol=1e-4, atol=1e-4,
    # )

    # # Check that the gradient w.r.t. the samples is correct (multi-sample).
    # adjoint_test(
    #     (θ, x, A)->rand(MersenneTwister(123456), g(θ, x, A), 11),
    #     randn(rng, length(x), 11), θ, x, A;
    #     rtol=1e-4, atol=1e-4,
    # )

    # # Check adjoints for logpdf.
    # adjoint_test(
    #     (θ, x, A, y)->logpdf(g(θ, x, A), y), randn(rng), θ, x, A, y;
    #     rtol=1e-4, atol=1e-4,
    # )

    # # Check adjoint for elbo.
    # adjoint_test(
    #     (ϴ, x, A, y, z, B)->begin
    #         fx, uz = h(θ, x, A, z, B)
    #         return elbo(fx, y, uz)
    #     end,
    #     randn(rng), θ, x, A, y, z, B;
    #     rtol=1e-4, atol=1e-4,
    # )

    # #
    # # multi-process rand / logpdf / elbo tests - this stuff won't work for anything if
    # # cross-related functionality doesn't work properly
    # #

    # adjoint_test(
    #     (θ, x, A) -> begin
    #         _g = g(θ, x, A)
    #         return rand(MersenneTwister(123456), [_g, _g])
    #     end,
    #     [randn(rng, length(x)), randn(rng, length(x))], θ, x, A;
    #     rtol=1e-4, atol=1e-4,
    # )

    # adjoint_test(
    #     (θ, x, A) -> begin
    #         _g = g(θ, x, A)
    #         return rand(MersenneTwister(123456), [_g, _g], 11)
    #     end,
    #     [randn(rng, length(x), 11), randn(rng, length(x), 11)], θ, x, A;
    #     rtol=1e-4, atol=1e-4,
    # )

    # adjoint_test(
    #     (θ, x, A, y) -> begin
    #         _g = g(θ, x, A)
    #         return logpdf([_g, _g], [y, y])
    #     end,
    #     randn(rng), θ, x, A, y;
    #     rtol=1e-4, atol=1e-4,
    # )

    # adjoint_test(
    #     (θ, x, A, y) -> begin
    #         _g = g(θ, x, A)
    #         return logpdf([_g ← y, _g ← y])
    #     end,
    #     randn(rng), θ, x, A, y;
    #     rtol=1e-4, atol=1e-4,
    # )

    # adjoint_test(
    #     (ϴ, x, A, y, z, B)->begin
    #         fx, uz = h(θ, x, A, z, B)
    #         return elbo([fx, fx], [y, y], [uz, uz])
    #     end,
    #     randn(rng), θ, x, A, y, z, B;
    #     rtol=1e-4, atol=1e-4,
    # )
end

# inputs(k::Kernel, N::Int) = collect(range(-5.0, 5.0; length=N))
# function inputs(k::BlockKernel, N::Int)
#     if size(k.ks, 1) == 1
#         return BlockData(collect(range(-5.0, 5.0; length=N)))
#     elseif size(k.ks, 1) == 2
#         x = collect(range(-5.0, 5.0; length=N))
#         return BlockData([x[1:div(N, 2)], x[div(N, 2) + 1:end]])
#     else
#         error("Incorrectly-sized block kernel")
#     end
# end

function standard_1D_dense_test(rng::AbstractRNG, θ, f, x::AV, z::AV)
    g, u = f(θ)
    N, M = length(x), length(z)
    A, B = randn(rng, N, N), randn(rng, M, M)
    y = rand(rng, first(f(θ))(x, _to_psd(A)))
    check_consistency(rng, θ, f, x, y, A, z, B)
end

function standard_1D_diag_test(rng::AbstractRNG, θ, f, x::AV, z::AV)
    g, u = f(θ)
    a, b = randn(rng, length(x)), randn(rng, length(z))
    y = rand(rng, first(f(θ))(x, _to_psd(a)))
    check_consistency(rng, θ, f, x, y, a, z, b)
end

function standard_1D_isotropic_test(rng::AbstractRNG, θ, f, x::AV, z::AV)
    g, u = f(θ)
    a, b = randn(rng), randn(rng)
    y = rand(rng, first(f(θ))(x, _to_psd(a)))
    check_consistency(rng, θ, f, x, y, a, z, b)
end

function standard_1D_tests(rng::AbstractRNG, θ, f, x::AV, z::AV)

    g, u = f(θ)

    @test cov(g, x) ≈ cov(g, x)'
    @test minimum(eigvals(cov(g, x))) > -1e-9
    @test cov_diag(g, x) ≈ diag(cov(g, x))

    @test cov(g, g, x, x) ≈ cov(g, x)
    @test cov_diag(g, g, x, x) ≈ diag(cov(g, g, x, x))

    standard_1D_dense_test(rng, θ, f, x, z)
    standard_1D_diag_test(rng, θ, f, x, z)
    standard_1D_isotropic_test(rng, θ, f, x, z)
end


text:
using LinearKernelModels
using LinearAlgebra, OffsetArrays
using Test

@testset "Single stimulus" begin
    S = rand(100)
    ktrue = OffsetArray(randn(5), -4:0)
    r = compute_r(S, ktrue)
    C, d = compute_Cd(S, r, 5)
    kcalc = reshape(C \ d, axes(ktrue))
    @test kcalc ≈ ktrue
    nfailed = 0
    for i = 1:100
        rnoisy = r .+ randn(size(r))
        C, d = compute_Cd(S, rnoisy, 5)
        knoisy = reshape(C \ d, axes(ktrue))
        kerr = compute_kerr(knoisy, C, S, rnoisy)
        nfailed += any((knoisy .- ktrue).^2 .> 3*kerr)
    end
    @test nfailed <= 5

    nfailed = 0
    for i = 1:100
        S = rand(10000)
        r = compute_r(S, ktrue)
        rnoisy = r .+ randn(size(r))
        C, d = compute_Cd(S, rnoisy, 5)
        knoisy = reshape(C \ d, axes(ktrue))
        kerr = compute_kerr(knoisy, C, S, rnoisy)
        nfailed += any((knoisy .- ktrue).^2 .> 3*kerr)
    end
    @test nfailed <= 5
end

@testset "Single acausal stimulus" begin
    S = rand(100)
    ktrue = OffsetArray(randn(5), -2:2)
    r = compute_r(S, ktrue)
    C, d = compute_Cd(S, r, axes(ktrue, 1))
    kcalc = reshape(C \ d, axes(ktrue))
    @test kcalc ≈ ktrue
end

@testset "Multiple stimuli and NaNs" begin
    S = rand(100, 3)
    ktrue = OffsetArray(randn(5, 3), -4:0, 1:3)
    r = compute_r(S, ktrue)
    C, d = compute_Cd(S, r, 5)
    kcalc = reshape(C \ d, axes(ktrue))
    @test kcalc ≈ ktrue

    r[5:10] .= NaN
    r[50:52] .= NaN
    C, d = compute_Cd(S, r, 5)
    kcalc = reshape(C \ d, axes(ktrue))
    @test kcalc ≈ ktrue
end

@testset "Multiple acausal stimuli and NaNs" begin
    S = rand(100,3)
    ktrue = OffsetArray(randn(5,3), -2:2, 1:3)
    r = compute_r(S, ktrue)
    C, d = compute_Cd(S, r, -2:2)
    kcalc = reshape(C \ d, axes(ktrue))
    @test kcalc ≈ ktrue

    r[5:10] .= NaN
    r[50:52] .= NaN
    C, d = compute_Cd(S, r, -2:2)
    kcalc = reshape(C \ d, axes(ktrue))
    @test kcalc ≈ ktrue
end

@testset "Single stimulus with constraints" begin
    S = rand(100)
    ktrue = OffsetArray([rand(3); 0; 0], -2:2)
    r = compute_r(S, ktrue)
    kcalc = solve_for_kernel(S, r, ktrue.!=0)
    @test kcalc ≈ ktrue
end

@testset "Multiple stimuli with constraints" begin
    S = rand(100,3)
    ktrue = OffsetArray(randn(5,3), -2:2, 1:3)
    ktrue[rand(1:15,5,1)] .= 0
    r = compute_r(S,ktrue)
    kcalc = solve_for_kernel(S, r, ktrue.!=0)
    @test kcalc ≈ ktrue

    # With a degeneracy that is broken by the constraint
    S = zeros(20, 2)
    S[10,:] .= 1
    r = zeros(20)
    r[11:12] .= 0.5
    kmask = OffsetArray([true false; false true], -2:-1, 1:2)
    kcalc = solve_for_kernel(S, r, kmask)
    @test kcalc ≈ OffsetArray([0.5 0; 0 0.5], -2:-1, 1:2)
end

@testset "Multiple stimuli with constraints and Nans" begin
    S = rand(100,7)
    ktrue = OffsetArray(rand(5,7), -2:2, 1:7)
    ktrue[rand(1:35, 10)] .= 0
    r = compute_r(S,ktrue)
    r[5:10] .= NaN
    r[50:52] .= NaN
    kcalc = solve_for_kernel(S, r, ktrue.!=0)
    @test kcalc ≈ ktrue
end


text:
struct IsData end
struct IsAntiData end
const PolarityTrait = Union{IsData,IsAntiData}

Base.adjoint(::IsData) = IsAntiData()
Base.adjoint(::IsAntiData) = IsData()

abstract type AbstractSSQNode{T} end

mutable struct SSQDataNode{T} <: AbstractSSQNode{T}
    @atomic next::Union{AbstractSSQNode{T},Nothing}
    value::T

    SSQDataNode{T}() where {T} = new{T}(nothing)
    SSQDataNode{T}(next::Union{AbstractSSQNode{T},Nothing}, value::T) where {T} =
        new{T}(next, value)
end

mutable struct SSQWaiterNode{T} <: AbstractSSQNode{T}
    @atomic next::Union{AbstractSSQNode{T},Nothing}
    value::Waiter{T}

    SSQWaiterNode{T}() where {T} = new{T}(nothing)
    SSQWaiterNode{T}(next::Union{AbstractSSQNode{T},Nothing}, value::Waiter{T}) where {T} =
        new{T}(next, value)
end

const SSQNode{T} = Union{SSQDataNode{T},SSQWaiterNode{T}}

polarityof(::SSQDataNode) = IsData
polarityof(::SSQWaiterNode) = IsAntiData
ssqnodetype(::IsData, ::Type{T}) where {T} = SSQDataNode{T}
ssqnodetype(::IsAntiData, ::Type{T}) where {T} = SSQWaiterNode{T}

#=
mutable struct SSQNode{Polarity<:PolarityTrait,T,V<:Union{T,Waiter{T}}}
    @atomic next::Union{SSQNode{IsData,T,T},SSQNode{IsAntiData,T,Waiter{T}},Nothing}
    value::V

    SSQNode{Polarity,T}() where {Polarity,T} = new{Polarity,T,vtype(Polarity(), T)}(nothing)
    function SSQNode{Polarity,T}(
        next::Union{SSQNode{<:Any,T},Nothing},
        value::Union{T,Waiter{T}},
    ) where {Polarity,T}
        V = vtype(Polarity(), T)
        value = value::V
        return new{Polarity,T,V}(next, value)
    end
end

@inline SSQNode{Polarity,T,V}(
    next::Union{SSQNode{<:Any,T},Nothing},
    value::Union{T,Waiter{T}},
) where {Polarity,T,V} = SSQNode{Polarity,T}(next, value)::SSQNode{Polarity,T,V}

polarityof(::SSQNode{Polarity}) where {Polarity} = Polarity
=#

mutable struct DualLinkedQueue{T}
    @atomic head::SSQNode{T}
    @atomic tail::SSQNode{T}
end

Base.eltype(::Type{DualLinkedQueue{T}}) where {T} = T

DualLinkedQueue() = DualLinkedQueue{Any}()
function DualLinkedQueue{T}() where {T}
    n = SSQDataNode{T}()
    return DualLinkedQueue{T}(n, n)
end

function denqueue!(
    queue::DualLinkedQueue{T},
    x::Union{T,Waiter{T}},
    polarity::PolarityTrait,
) where {T}
    local node::Union{Nothing,SSQNode{T}} = nothing

    head = @atomic queue.head
    tail = @atomic queue.tail
    while true
        next = (@atomic head.next)::Union{Nothing,SSQNode{T}}
        head′ = @atomic queue.head
        if head !== head′  # snapshot failed
            head = head′
            tail = @atomic queue.tail
            continue
        end
        if head === tail
            last = next
            should_enqueue = true
        elseif polarity isa polarityof(next)
            last = @atomic tail.next
            should_enqueue = true
        else
            last = nothing
            should_enqueue = false
        end

        if should_enqueue
            tail′ = @atomic queue.tail
            if tail′ !== tail  # snapshot failed
                tail = tail′
                continue
            end
            if last !== nothing
                old, ok = @atomicreplace(queue.tail, tail => last)
                tail = ok ? next : old
                continue
            end
            node = if node === nothing
                ssqnodetype(polarity, T)(nothing, x)
            else
                node
            end::ssqnodetype(polarity, T)
            last, ok = @atomicreplace(tail.next, nothing => node)
            if ok
                @atomicreplace(queue.tail, tail => node)
                return nothing
            end
            last = last::SSQNode{T}  # can be any polarity
            old, ok = @atomicreplace(queue.tail, tail => last)
            tail = ok ? last : old
        else
            next = next::ssqnodetype(polarity', T)
            value = next.value
            head, ok = @atomicreplace(queue.head, head => next)
            if ok
                return Some(value)
            end
            tail = @atomic queue.tail
        end
    end
end

function Base.push!(queue::DualLinkedQueue{T}, x) where {T}
    x = convert(T, x)
    while true
        y = denqueue!(queue, x, IsData())
        if y isa Some
            w = something(y)::Waiter{T}
            tryput!(w, x) || continue
        else
            y::Nothing
        end
        return queue
    end
end

function Base.popfirst!(queue::DualLinkedQueue{T}) where {T}
    # TODO: cache Waiter
    w = Waiter{T}()
    y = denqueue!(queue, w, IsAntiData())
    if y isa Some
        return something(y)
    else
        y::Nothing
        x = fetch(w)
        return x::T
    end
end

Base.IteratorSize(::Type{<:DualLinkedQueue}) = Base.SizeUnknown()

function Base.iterate(queue::DualLinkedQueue{T}) where {T}
    head = @atomic queue.head
    node = @atomic head.next
    if node === nothing
        return nothing
    else
        if node isa SSQWaiterNode
            return nothing
        else
            return (node.value, node)
        end
    end
end

function Base.iterate(::DualLinkedQueue{T}, prev::SSQDataNode) where {T}
    node = (@atomic prev.next)::Union{Nothing,SSQDataNode}
    if node === nothing
        return nothing
    else
        return (node.value, node)
    end
end

struct NodeIterator{T}
    queue::T
end

Base.IteratorSize(::Type{<:NodeIterator}) = Base.SizeUnknown()

Base.eltype(::Type{NodeIterator{DualLinkedQueue{T}}}) where {T} = SSQNode{T}

function Base.iterate(
    iter::NodeIterator{DualLinkedQueue{T}},
    prev::SSQNode{T} = let queue = iter.queue
        @atomic queue.head
    end,
) where {T}
    node = (@atomic prev.next)::Union{Nothing,SSQNode{T}}
    if node === nothing
        return nothing
    else
        return (node, node)
    end
end

function check_invariance(queue::DualLinkedQueue{T}) where {T}
    isdata = nothing
    for node in NodeIterator(queue)
        if isdata === nothing
            isdata = node isa SSQDataNode
        elseif isdata !=′ (node isa SSQDataNode)
            return false
        end
    end
    return true
end

function summaryinfo(queue::DualLinkedQueue{T}) where {T}
    counter = Ref(0)
    isdata = Ref(true)
    for node in NodeIterator(queue)
        isdata[] = !(node isa SSQWaiterNode)
        counter[] += 1
    end
    return (; nitems = counter[], isdata = isdata[])
end

function Base.summary(io::IO, queue::DualLinkedQueue)
    show(io, MIME"text/plain"(), typeof(queue))
    nitems, isdata = summaryinfo(queue)
    s = nitems > 1 ? "s" : ""
    print(io, " with ", nitems, isdata ? " data item$s" : " waiter$s")
end

function Base.show(io::IO, ::MIME"text/plain", queue::DualLinkedQueue)
    summary(io, queue)
    eio = IOContext(io, :typeinfo => eltype(queue), :limit => true, :compact => true)
    n = 0
    for x in queue
        if n == 0
            println(io, ":")
        else
            println(io)
        end
        n += 1
        if n > 3
            print(io, "  ⋮")
            return
        end
        print(io, "  ")
        show(eio, MIME"text/plain"(), x)
    end
end


text:
import DelimitedFiles

data = DelimitedFiles.readdlm("data/day_three.txt", ' ')

function count(data, rsteps, dsteps)
    tot = 0
    for (idx, row) in enumerate(data)
        if (idx - 1) % dsteps == 0
            slope = rsteps / dsteps
            row_idx = (Int((idx - 1) * slope) % length(row)) + 1
            if row[row_idx] == '#'
                tot += 1
            end
        end
    end
    return tot
end

function solve()
    tot = 1
    for (rstep, dstep) in ((1, 1), (3, 1), (5, 1), (7, 1), (1, 2))
        tot *= count(data, rstep, dstep)
    end
    return tot
end

println(count(data, 3, 1))
println(solve())

max_stars_repo_name: ColCarroll/adventofcode

max_issues_repo_name: ColCarroll/adventofcode

max_forks_repo_name: ColCarroll/adventofcode


text:
#=##############################################################################
# DESCRIPTION
    Auxiliary functions for rotor geometry generation and calculations.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

"""
Generates the Rotor object. `pitch` is the pitch of the blades in degrees,
`n` is the number of lattices in the VLM.

`ReD` is the diameter Reynolds number based on rotational speed calculated
as ReD = (omega*R)*(2*R)/nu, and `Matip` is the rotational+freestream Mach
number at the tip. This values are used for calculate airfoil properties
through XFOIL (chord Reynolds), so ignore them if airfoil properties are
prescribed.

Give it `altReD = [RPM, J, nu]`, and it'll calculate the chord-based Reynolds
accounting for the effective velocity at every station ignoring ReD (this is
more accurate, but not needed).

NOTE: If Matip is different than zero while running XFOIL, remember to deactive
compressibility corrections when running `vlm.solvefromCCBlade()` by giving it
`sound_spd=nothing`
"""
function generate_rotor(Rtip::Real, Rhub::Real, B::Int,
                        chorddist::Array{Float64,2},
                        pitchdist::Array{Float64,2},
                        sweepdist::Array{Float64,2},
                        heightdist::Array{Float64,2},
                        airfoil_files::Array{Tuple{Float64,String,String},1};
                        # INPUT OPTIONS
                        data_path=def_data_path,
                        read_polar=vlm.ap.read_polar,
                        # PROCESSING OPTIONS
                        pitch=0.0,
                        n=10, CW=true, blade_r=1.0,
                        ReD=5*10^5, altReD=nothing, Matip=0.0,
                        ncrit=9,
                        xfoil=false,
                        rotor_file="apc10x7.jl",
                        spline_k=5, spline_s=0.001, splines_s=nothing, spline_bc="extrapolate",
                        turbine_flag=false,
                        rfl_n_lower=15, rfl_n_upper=15,
                        rediscretize_airfoils=true,
                        # OUTPUT OPTIONS
                        verbose=false, v_lvl=1,
                        plot_disc=true, figsize_factor=2/3)
    if verbose; println("\t"^v_lvl*"Generating geometry..."); end;
    n_bem = n

    # Read airfoil contours
    # Airfoils along the blade as
    # airfoil_contours=[ (pos1, contour1, polar1), (pos2, contour2, pol2), ...]
    # with contour=(x,y) and pos the position from root to tip between 0 and 1.
    # pos1 must equal 0 (root airfoil) and the last must be 1 (tip airfoil)
    airfoil_contours = []
    airfoil_path = joinpath(data_path, "airfoils")
    for (r, rfl_file, clcurve_file) in airfoil_files
        x,y = gt.readcontour(rfl_file; delim=",", path=airfoil_path, header_len=1)
        rfl = hcat(x,y)

        push!(airfoil_contours, (r, rfl, clcurve_file))
    end

    # Splines
    _spl_chord = Dierckx.Spline1D(chorddist[:, 1]*Rtip, chorddist[:, 2]*Rtip;
                                        k= size(chorddist)[1]>2 ? spline_k : 1,
                                        s=splines_s!=nothing ? splines_s[1] : spline_s, bc=spline_bc)
    _spl_theta = Dierckx.Spline1D(pitchdist[:, 1]*Rtip, pitchdist[:, 2];
                                        k= size(pitchdist)[1]>2 ? spline_k : 1,
                                        s=splines_s!=nothing ? splines_s[2] : spline_s, bc=spline_bc)
    _spl_LE_x = Dierckx.Spline1D(sweepdist[:, 1]*Rtip, sweepdist[:, 2]*Rtip;
                                        k= size(sweepdist)[1]>2 ? spline_k : 1,
                                        s=splines_s!=nothing ? splines_s[3] : spline_s, bc=spline_bc)
    _spl_LE_z = Dierckx.Spline1D(heightdist[:, 1]*Rtip, heightdist[:, 2]*Rtip;
                                        k= size(heightdist)[1]>2 ? spline_k : 1,
                                        s=splines_s!=nothing ? splines_s[4] : spline_s, bc=spline_bc)
    spl_chord(x) = Dierckx.evaluate(_spl_chord, x)
    spl_theta(x) = pitch + Dierckx.evaluate(_spl_theta, x)
    spl_LE_x(x) = Dierckx.evaluate(_spl_LE_x, x)
    spl_LE_z(x) = Dierckx.evaluate(_spl_LE_z, x)

    # Geometry for CCBlade & FLOWVLM
    r = [Rhub + i*(Rtip-Rhub)/n_bem for i in 0:n_bem] # r is discretized in n+1 sections
    chord = [spl_chord(ri) for ri in r]
    theta = [spl_theta(ri) for ri in r]
    LE_x = [spl_LE_x(ri) for ri in r]
    LE_z = [spl_LE_z(ri) for ri in r]

    if verbose; println("\t"^v_lvl*"Generating airfoils..."); end;

    airfoils = []
    Mas = xfoil ? [] : nothing
    for (pos, contour, file_name) in airfoil_contours
        x, y = contour[:,1], contour[:,2]
        # Calls XFOIL to calculate polars of each airfoil
        if xfoil
            roR = (Rhub + pos*(Rtip-Rhub))/Rtip       # (r/R) position along blade.

            if altReD!=nothing
                RPM, J, nu = altReD
                this_ReD = calc_ReD(nu, 2*Rtip, RPM, J, roR*Rtip) # Diameter Reynolds at r
                this_Re = ceil(Int, this_ReD*spl_chord(roR*Rtip)/(2Rtip)) # Chord Reynolds

                if Matip != 0
                    speedofsnd = (2*pi*RPM/60)*Rtip / Matip
                    this_Ma = calc_Veff(2*R, RPM, J, roR*Rtip) / speedofsnd
                else
                    this_Ma = Matip*roR
                end
            else
                    # Chord-based Reynolds number at
                this_Re = Int(ceil(ReD*roR*spl_chord(roR*Rtip)/(2*Rtip))) # this position
                this_Ma = Matip*roR
            end

            push!(Mas, this_Ma)

            polar = vlm.ap.runXFOIL(x, y, this_Re;
                                    alphas=[i for i in -20:1.0:20],
                                    verbose=verbose, Mach=this_Ma,
                                    iter=100, ncrit=ncrit)

        else # Reads polars from files
            if verbose; println("\t"^(v_lvl+1)*"$file_name"); end;
            polar = read_polar(file_name; path=data_path*"airfoils/", x=x, y=y)
        end

        push!(airfoils, (pos, polar))
    end

    if verbose; println("\t"^v_lvl*"Generating FLOWVLM Rotor..."); end;
    propeller = vlm.Rotor(CW, r, chord, theta, LE_x, LE_z, B, airfoils, turbine_flag)

    vlm.initialize(propeller, n; r_lat=blade_r, verif=plot_disc,
                    genblade_args=[(:spl_k,spline_k), (:spl_s,spline_s)],
                    rfl_n_lower=rfl_n_lower, rfl_n_upper=rfl_n_upper,
                    figsize_factor=figsize_factor,
                    rediscretize=rediscretize_airfoils)

    if plot_disc
        fig = figure("discretization_verif", figsize=[7*2,5*1]*figsize_factor)
        suptitle("Discretization Verification")

        subplot(121)
        plot(chorddist[:, 1], chorddist[:, 2], "ok", label="Chord data", alpha=0.75)
        plot(r/Rtip, chord/Rtip, "--or", label="Chord Spline", alpha=0.75)
        plot(sweepdist[:, 1], sweepdist[:, 2], "^k", label="LE-x data", alpha=0.75)
        plot(r/Rtip, LE_x/Rtip, "--^g", label="LE-x Spline", alpha=0.75)
        plot(heightdist[:, 1], heightdist[:, 2], "*k", label="LE-z data", alpha=0.75)
        plot(r/Rtip, LE_z/Rtip, "--*b", label="LE-z Spline", alpha=0.75)
        xlabel(L"$r/R$")
        ylabel(L"$c/R$, $x/R$, $z/R$")
        legend(loc="best", frameon=true)
        grid(true, color="0.8", linestyle="--")

        subplot(122)
        plot(pitchdist[:, 1], pitchdist[:, 2], "ok", label="Twist data", alpha=0.75)
        plot(r/Rtip, theta, "--^r", label="Twist Spline", alpha=0.75)
        xlabel(L"$r/R$")
        ylabel(L"Twist $\theta$ ($^\circ$)")
        legend(loc="best", frameon=true)
        grid(true, color="0.8", linestyle="--")

        for (i,(pos, polar)) in enumerate(airfoils)
            vlm.ap.plot(polar; geometry=true, label="pos=$pos, Re=$(vlm.ap.get_Re(polar))"*
                                        (Mas!=nothing ? ", Ma=$(round(Mas[i], digits=2))" : ""),
                    cdpolar=false, fig_id="prelim_curves", title_str="Re sweep", rfl_figfactor=figsize_factor)
        end
    end

    return propeller
end

function generate_rotor(Rtip::Real, Rhub::Real, B::Int, blade_file::String;
                        data_path=def_data_path, optargs...)

    (chorddist, pitchdist,
     sweepdist, heightdist,
     airfoil_files, spl_k,
     spl_s) = read_blade(blade_file; data_path=data_path)

    return generate_rotor(Rtip, Rhub, B, chorddist, pitchdist, sweepdist,
                            heightdist, airfoil_files;
                            data_path=data_path,
                            spline_k=spl_k, spline_s=spl_s, optargs...)
end

function generate_rotor(rotor_file::String;
                        data_path=def_data_path, optargs...)

    Rtip, Rhub, B, blade_file = read_rotor(rotor_file; data_path=data_path)

    return generate_rotor(Rtip, Rhub, B, blade_file;
                            data_path=data_path, optargs...)
end

function read_rotor(rotor_file::String; data_path=def_data_path)

    # Path to rotor files
    rotor_path = joinpath(data_path, "rotors")

    data = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, rotor_file)))
    Rtip = Meta.parse(data[1, 2])
    Rhub = Meta.parse(data[2, 2])
    B = Meta.parse(data[3, 2])
    blade_file = data[4, 2]

    return Rtip, Rhub, B, blade_file
end

function read_blade(blade_file::String; data_path=def_data_path)

    # Path to rotor files
    rotor_path = joinpath(data_path, "rotors")

    # Read blade
    files = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, blade_file)))
    chorddist = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, files[1, 2])))
    pitchdist = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, files[2, 2])))
    sweepdist = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, files[3, 2])))
    heightdist = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, files[4, 2])))
    airfoil_files = DataFrames.DataFrame!(CSV.File(joinpath(rotor_path, files[5, 2])))
    spl_k = Meta.parse(files[6, 2])
    spl_s = Meta.parse(files[7, 2])

    # Convert DataFrames to concrete types
    chorddist = Array{Float64, 2}(chorddist)
    pitchdist = Array{Float64, 2}(pitchdist)
    sweepdist = Array{Float64, 2}(sweepdist)
    heightdist = Array{Float64, 2}(heightdist)

    af = airfoil_files
    airfoil_files = [(Float64(af[i, 1]), String(af[i, 2]), String(af[i, 3]))
                        for i in 1:size(af, 1)]

    return chorddist, pitchdist, sweepdist, heightdist, airfoil_files, spl_k, spl_s
end

"""
    Given a diameter-based Reynolds number \$\\mathrm{Re}_D(r)\$ at a certain
radial position r and an advance ratio J, return the rotational speed
"""
calc_n(ReD, J, r, D, nu) = nu/( D * sqrt( (2*pi*r)^2 + (J*D)^2 ) ) * ReD
calc_RPM(args...) = 60*calc_n(args...)

calc_Veff(D, RPM, J, r) = RPM/60 * sqrt((2*pi*r)^2 + (J*D)^2)

"Reynolds number at radial station r (m)"
calc_ReD(nu, D, args...) = D/nu * calc_Veff(D, args...)

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
include("graph.jl")
include("opinion.jl")

function main()
    (W, U, n) = read_edgelist("./data/raw/cloister.edgelist")
    m = 2
    normalize_edgelist!(W)
    normalize_edgelist!(U)
    y = generate_opinions(n, m)
    α, β, γ = generate_factors(n)
    M = Matrix(β .* W - γ .* U)
    Q = get_fundamental_matrix(β, W, γ, U)
    z∞ = get_equlibrium(α, Q, y)
    return y, Q, z∞
end

main()

max_stars_repo_name: InzamamRahaman/SignedDegroot

max_issues_repo_name: InzamamRahaman/SignedDegroot

max_forks_repo_name: InzamamRahaman/SignedDegroot


text:
using NonconvexCore: flatten
using OrderedCollections, JuMP, Zygote, SparseArrays, LinearAlgebra

struct SS
    a
    b
end

xs = [
    1.0,
    [1.0],
    [1.0, 2.0],
    [1.0, [1.0, 2.0]],
    [1.0, (1.0, 2.0)],
    [1.0, OrderedDict(1 => [1.0, 2.0])],
    [[1.0], OrderedDict(1 => [1.0, 2.0])],
    [(1.0,), [1.0,], OrderedDict(1 => [1.0, 2.0])],
    [1.0 1.0; 1.0 1.0],
    rand(2, 2, 2),
    [[1.0, 2.0], [3.0, 4.0]],
    OrderedDict(1 => 1.0),
    OrderedDict(1 => [1.0]),
    OrderedDict(1 => 1.0, 2 => [2.0]),
    OrderedDict(1 => 1.0, 2 => [2.0], 3 => [[1.0, 2.0], [3.0, 4.0]]),
    JuMP.Containers.DenseAxisArray(reshape([1.0, 1.0], (2,)), 1),
    (1.0,),
    (1.0, 2.0),
    (1.0, (1.0, 2.0)),
    (1.0, [1.0, 2.0]),
    (1.0, OrderedDict(1 => [1.0, 2.0])),
    ([1.0], OrderedDict(1 => [1.0, 2.0])),
    ((1.0,), [1.0,], OrderedDict(1 => [1.0, 2.0])),
    (a = 1.0,),
    (a = 1.0, b = 2.0),
    (a = 1.0, b = (1.0, 2.0)),
    (a = 1.0, b = [1.0, 2.0]),
    (a = 1.0, b = OrderedDict(1 => [1.0, 2.0])),
    (a = [1.0], b = OrderedDict(1 => [1.0, 2.0])),
    (a = (1.0,), b = [1.0,], c = OrderedDict(1 => [1.0, 2.0])),
    sparsevec([1.0, 2.0], [1, 3], 10),
    sparse([1, 2, 2, 3], [2, 3, 1, 4], [1.0, 2.0, 3.0, 4.0], 10, 10),
    SS(1.0, 2.0),
    [SS(1.0, 2.0), 1.0],
]

for x in xs
    @show x
    xvec, unflatten = flatten(x)
    @test x == unflatten(xvec)
    J = Zygote.jacobian(xvec) do x
        unflatten(x)
        flatten(x)[1]
    end[1]
    @test logabsdet(J) == (0.0, 1.0)
end


text:
using Luna

a = 13e-6
gas = :Ar
pres = 5

τfwhm = 30e-15
λ0 = 800e-9
energy = 1e-6

L = 15e-2

grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12)

a0 = a
aL = 3a/4

afun = let a0=a0, aL=aL, L=L
    afun(z) = a0 + (aL-a0)*z/L
end

modes = (
    Capillary.MarcatiliMode(afun, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    Capillary.MarcatiliMode(afun, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false)
)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
            #  Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

linop = LinearOps.make_linop(grid, modes, λ0);

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes, :y; full=false)

statsfun = Stats.default(grid, Eω, modes, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
Plotting.time_1D(output, [5e-2, 9.8e-2])
Plotting.time_1D(output, [5e-2, 9e-2], modes=:sum, bandpass=(150e-9, 300e-9))
Plotting.spec_1D(output, [5e-2, 9.8e-2])
Plotting.spec_1D(output, [5e-2, 9.8e-2], modes=:sum)


text:
for (func, elty) in [(:CLBlastSsymm, Float32), (:CLBlastDsymm, Float64),
                     (:CLBlastCsymm, ComplexF32), (:CLBlastZsymm, ComplexF64)]
    #TODO: (:CLBlastHsymm, Float16)

    @eval function $func(layout::CLBlastLayout,
                         side::CLBlastSide, triangle::CLBlastTriangle,
                         m::Integer, n::Integer,
                         alpha::$elty,
                         a_buffer::cl.CL_mem, a_offset::Integer, a_ld::Integer,
                         b_buffer::cl.CL_mem, b_offset::Integer, b_ld::Integer,
                         beta::$elty,
                         c_buffer::cl.CL_mem, c_offset::Integer, c_ld::Integer,
                         queue::cl.CmdQueue, event::cl.Event)
        err = ccall(
            ($(string(func)), libCLBlast),
            cl.CL_int,
            (Cint, Cint, Cint, Csize_t, Csize_t, $elty, Ptr{Cvoid}, Csize_t, Csize_t,
              Ptr{Cvoid}, Csize_t, Csize_t, $elty, Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}),
            Cint(layout), Cint(side), Cint(triangle), m, n, alpha, a_buffer, a_offset, a_ld,
              b_buffer, b_offset, b_ld, beta, c_buffer, c_offset, c_ld, Ref(queue), Ref(event)
        )
        if err != cl.CL_SUCCESS
            println(stderr, "Calling function $(string($func)) failed!")
            throw(CLBlastError(err))
        end
        return err
    end

    @eval function symm!(_side::Char, uplo::Char, α::Number, A::cl.CLArray{$elty,2},
                         B::cl.CLArray{$elty,2}, β::Number, C::cl.CLArray{$elty,2};
                         queue::cl.CmdQueue=cl.queue(C))
        # check and convert arguments
        if _side == 'L'
            side = CLBlastSideLeft
        elseif _side == 'R'
            side = CLBlastSideRight
        else
            throw(ArgumentError("Side marker `side` is $(_side) but only 'L' and 'R' are allowed."))
        end
        if uplo == 'U'
            triangle = CLBlastTriangleUpper
        elseif uplo == 'L'
            triangle = CLBlastTriangleLower
        else
            throw(ArgumentError("Upper/lower marker `uplo` is $(uplo) but only 'U' and 'L' are allowed."))
        end
        m, n = size(C)
        j = size(A,1)
        if j != size(A,2)
            throw(DimensionMismatch("`A` has dimensions $(size(A)) but must be square."))
        end
        if j != (_side == 'L' ? m : n)
            throw(DimensionMismatch("`A` has size $(size(A)), `C` has size ($m,$n)."))
        end
        if size(B,2) != n
            throw(DimensionMismatch("`B` has second dimension $(size(B,2)) but needs to match second dimension of `C`, $n."))
        end
        alpha = convert($elty, α)
        beta  = convert($elty, β)
        layout = CLBlastLayoutColMajor

        # output event
        event::cl.Event = cl.Event(C_NULL)

        $func(layout, side, triangle,
              m, n,
              alpha,
              pointer(A), 0, size(A,1),
              pointer(B), 0, size(B,1),
              beta,
              pointer(C), 0, size(C,1),
              queue, event)

        # wait for kernel
        cl.wait(event)

        C
    end

end


text:
# Unit tests for dftk_integration.jl

@testset "dftk_integration.jl" begin
    # TODO
end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
using MLDatasets
using Random

function dataloader(dataname::String)
    if dataname == "MNIST" || dataname == "mnist"
        train_x, train_y = MNIST.traindata()
        test_x, test_y = MNIST.testdata()
        train_x, train_y = collect(Float64, train_x), collect(Float64, train_y)
        test_x, test_y = collect(Float64, test_x), collect(Float64, test_y)
        return train_x, train_y, test_x, test_y
    elseif dataname == "iris"
        return Iris.features(), Iris.labels()
    elseif dataname == "BostonHousing"
        return BostonHousing.features(), BostonHousing.targets()
    end
end

struct DataSplitter
    #TODO: make it possible to DataFrame split
    splitted_data::Tuple
    function DataSplitter(xs, narray, ndata, sb::Bool, seed, test_size::Int)
        if sb == true
            seed != nothing && Random.seed!(seed)
            narray = shuffle(narray)
        end
        narray = (narray[1:test_size-1], narray[test_size:end])
        splitted_data = []
        for x in xs
            push!(splitted_data, index2element(x, narray))
        end
        new(Tuple(splitted_data))
    end
end

function index2element(x, index::Tuple)
    tr, te = index
    if ndims(x) == 1
        return ([x[i] for i in tr], [x[i] for i in te])
    else
        return ([x[i, :] for i in tr], [x[i, :] for i in te])
    end
end

#ndata is dims = 1
function split_data(x...; sb = true, train_size = nothing, test_size = nothing, seed = nothing)
    x = collect(x)
    ndata = size(x[1])[1]
    narray = collect(1 : 1 : ndata)
    test_size == train_size == nothing && throw(ArgumentError("test_size and train_size wasn't specified. you must specify either one."))
    if test_size == nothing
        #number of train data is ceiled
        test_size = (0<=train_size<=1) ? ndata-ceil(train_size*ndata) : ndata-train_size
    else
        test_size = (0<=test_size<=1) ? ndata-ceil(train_size*ndata) : ndata-train_size
    end
    return DataSplitter(x, narray, ndata, sb, seed, Int(test_size)).splitted_data
end


text:
# I want to try a pure MCTS simulation that doesn't rely on a neural network and
# see how things

using Game2048

using Game2048
using Game2048: bitboard_to_array

function find_best_move(board, n=10000)
    res = zeros(Int, 4)
    for direction in DIRS
        dir_idx = Int(direction) + 1
        new_board = move(board, direction)
        if new_board == board # not a valid move
            res[dir_idx] = -1
            continue
        end

        for _ in 1:n
            # random add one tile
            new_board = add_tile(new_board)

            # simulate from this onward
            res[dir_idx] += sum(2 .^ bitboard_to_array(simulate_bb(new_board)))
        end
    end

    DIRS[argmax(res)]
end

function play_via_monte_carlo(bitboard, n)
     while true
        best_move = find_best_move(bitboard, n)
        # println(best_move)
        new_bitboard = move(bitboard, best_move)
        # display(new_bitboard)

        if new_bitboard != bitboard
            bitboard = add_tile(new_bitboard)
            if bitboard == new_bitboard
                error("wtf")
            end
        else
            # lost
            return bitboard
        end
    end
end

board = initbboard()

@time play_via_monte_carlo(board, 100000)


text:
#######################################################################
# Field
#######################################################################
struct Field{F,N<:NamedTuple,L<:Lattice}
    f::F
    lat::L
    args::N
end

Field(::Missing, lat) = missing
Field(f::F, lat) where {F<:Function} = Field(f, lat, NamedTuple())

#######################################################################
# applyfield
#######################################################################
@inline applyfield(::Missing, h, i, j, dn) = h
function applyfield(field::Field{F}, h, i, j, dn) where {F<:Function}
    sites = field.lat.unitcell.sites
    r, dr = _rdr(sites[j], sites[i]) # rsource, rtarget
    r0 = field.lat.bravais.matrix * dn
    return field.f(r + r0, dr, h)
end


text:
module TestMachines

using MLJBase
using Test
using Statistics
using ..Models
using StableRNGs
using Serialization
using ..TestUtilities

const MLJModelInterface = MLJBase.MLJModelInterface

N=50
X = (a=rand(N), b=rand(N), c=rand(N));
y = 2*X.a - X.c + 0.05*rand(N);

train, test = partition(eachindex(y), 0.7);

tree = DecisionTreeRegressor(max_depth=5)
pca = PCA()

@testset "_contains_unknown" begin
    @test MLJBase._contains_unknown(Unknown)
    @test MLJBase._contains_unknown(Tuple{Unknown})
    @test MLJBase._contains_unknown(Tuple{Unknown, Int})
    @test MLJBase._contains_unknown(Union{Tuple{Unknown}, Tuple{Int,Char}})
    @test MLJBase._contains_unknown(Union{Tuple{Int}, Tuple{Int,Unknown}})
    @test !MLJBase._contains_unknown(Int)
    @test !MLJBase._contains_unknown(Tuple{Int})
    @test !MLJBase._contains_unknown(Tuple{Char, Int})
    @test !MLJBase._contains_unknown(Union{Tuple{Int}, Tuple{Int,Char}})
    @test !MLJBase._contains_unknown(Union{Tuple{Int}, Tuple{Int,Char}})
end

t = machine(tree, X, y)
@test_throws MLJBase.NotTrainedError(t, :fitted_params) fitted_params(t)
@test_throws MLJBase.NotTrainedError(t, :report) report(t)
@test_throws MLJBase.NotTrainedError(t, :training_losses) training_losses(t)
@test_logs (:info, r"Training") fit!(t)
@test_logs (:info, r"Training") fit!(t, rows=train)
@test_logs (:info, r"Not retraining") fit!(t, rows=train)
@test_logs (:info, r"Training") fit!(t)
t.model.max_depth = 1
@test_logs (:info, r"Updating") fit!(t)

@test training_losses(t) === nothing

predict(t, selectrows(X,test));
@test rms(predict(t, selectrows(X, test)), y[test]) < std(y)

mach = machine(ConstantRegressor(), X, y)
@test_logs (:info, r"Training") fit!(mach)
yhat = predict_mean(mach, X);

n = nrows(X)
@test rms(yhat, y) ≈ std(y)*sqrt(1 - 1/n)

# test an unsupervised univariate case:
mach = machine(UnivariateStandardizer(), float.(1:5))
@test_logs (:info, r"Training") fit!(mach)
@test isempty(params(mach))

# test a frozen Machine
stand = machine(Standardizer(), source((x1=rand(10),)))
freeze!(stand)
@test_logs (:warn, r"not trained as it is frozen\.$") fit!(stand)

@testset "warnings" begin
    @test_throws DimensionMismatch machine(tree, X, y[1:end-1])

    # supervised model with bad target:
    @test_logs((:warn,
                MLJBase.alert_generic_scitype_mismatch(
                    Tuple{scitype(X), AbstractVector{Multiclass{N}}},
                    MLJBase.fit_data_scitype(tree),
                    typeof(tree)
                    )
                ),
               machine(tree, X, categorical(1:N)))

    # ordinary transformer:
    @test_logs((:warn,
                MLJBase.alert_generic_scitype_mismatch(
                    Tuple{scitype(42),},
                    MLJBase.fit_data_scitype(pca),
                    typeof(pca)
                    )
                ),
               machine(pca, 42))
    y2 = coerce(1:N, OrderedFactor);

    # bad weight vector:
    @test_logs((:warn,
                MLJBase.alert_generic_scitype_mismatch(
                    Tuple{scitype(X), scitype(y2), scitype(42)},
                    MLJBase.fit_data_scitype(ConstantClassifier()),
                    ConstantClassifier
                    )
                ),
               machine(ConstantClassifier(), X, y2, 42))
end

struct FooBar <: Model end

MLJBase.fit_data_scitype(::Type{<:FooBar}) =
    Union{Tuple{AbstractVector{Count}},
          Tuple{AbstractVector{Count},AbstractVector{Continuous}}}

struct FooBarUnknown <: Model end

@testset "machine scitype_check_level" begin

    X = [1, 2, 3, 4]
    y = rand(4)

    # with no Unknown scitypes

    model = FooBar()

    for scitype_check_level in [1, 2]
        @test_logs machine(model, X, y; scitype_check_level)
        @test_logs machine(model, X; scitype_check_level)
        @test_logs((:warn,
                    MLJBase.alert_generic_scitype_mismatch(Tuple{scitype(y)},
                                                          fit_data_scitype(model),
                                                          FooBar)),
                   machine(model, y; scitype_check_level))
    end

    scitype_check_level = 3
    @test_logs machine(model, X, y; scitype_check_level)
    @test_logs machine(model, X; scitype_check_level)
    @test_throws(ArgumentError(
        MLJBase.alert_generic_scitype_mismatch(Tuple{scitype(y)},
                                              fit_data_scitype(model),
                                              FooBar)),
                 machine(model, y; scitype_check_level))

    @test default_scitype_check_level() == 1
    default_scitype_check_level(3)
    @test default_scitype_check_level() == 3

    @test_logs machine(model, X, y)
    @test_logs machine(model, X)
    @test_throws(ArgumentError(
        MLJBase.alert_generic_scitype_mismatch(Tuple{scitype(y)},
                                              fit_data_scitype(model),
                                              FooBar)),
                 machine(model, y))
    default_scitype_check_level(1)

    # with Unknown scitypes

    model = FooBarUnknown()

    scitype_check_level = 1
    @test_logs machine(model, X, y; scitype_check_level)
    @test_logs machine(model, X; scitype_check_level)

    warning = MLJBase.WARN_UNKNOWN_SCITYPE
    for scitype_check_level in [2, 3]
        @test_logs (:warn, warning) machine(model, X, y; scitype_check_level)
        @test_logs (:warn, warning) machine(model, X; scitype_check_level)
    end

    scitype_check_level = 4
    @test_throws ArgumentError(warning) machine(model, X, y; scitype_check_level)
    @test_throws ArgumentError(warning) machine(model, X; scitype_check_level)

end

@testset "weights" begin
    yraw = ["Perry", "Antonia", "Perry", "Skater"]
    X = (x=rand(4),)
    y = categorical(yraw)
    w = [2, 3, 2, 5]

    # without weights:
    mach = machine(ConstantClassifier(), X, y)
    fit!(mach, verbosity=0)
    d1 = predict(mach, X)[1]
    d2 = MLJBase.UnivariateFinite([y[1], y[2], y[4]], [0.5, 0.25, 0.25])
    @test all([pdf(d1, c) ≈ pdf(d2, c) for c in MLJBase.classes(d1)])

    # with weights:
    mach = machine(ConstantClassifier(), X, y, w)
    fit!(mach, verbosity=0)
    d1 = predict(mach, X)[1]
    d2 = MLJBase.UnivariateFinite([y[1], y[2], y[4]], [1/3, 1/4, 5/12])
    @test all([pdf(d1, c) ≈ pdf(d2, c) for c in MLJBase.classes(d1)])
end

mutable struct Scale <: MLJBase.Static
    scaling::Float64
end

function MLJBase.transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X * s.scaling
    MLJBase.table(s.scaling * MLJBase.matrix(X), prototype=X)
end

function MLJBase.inverse_transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X / s.scaling
    MLJBase.table(MLJBase.matrix(X) / s.scaling, prototype=X)
end

@testset "static transformer machines" begin
    s = Scale(2)
    X = ones(2, 3)
    Xt = MLJBase.table(X)

    @test_throws ArgumentError machine(Scale(2), X)
    @test_throws ArgumentError machine(Scale(2), source(X))

    mach = machine(Scale(2))
    @test_logs (:info, r"Training") fit!(mach) # no-op
    state = mach.state

    R  = transform(mach, X)
    IR = inverse_transform(mach, R)
    @test IR ≈ X

    # changing rows does not alter state (and "training" is skipped):
    @test_logs (:info, r"Not retraining") fit!(mach, rows=1:3)
    @test mach.state == state

    # changing hyper-parameters *does* change state (and "training" is
    # not skipped):
    mach.model.scaling = 3.0
    @test_logs (:info, r"Updating") fit!(mach, rows=1:3)
    @test mach.state != state

    @test_throws ArgumentError transform(mach, rows=1:2)
end

mutable struct Box
    matrix::Matrix{Int}
end

## DUMMY UNSUPERVISED MODEL

mutable struct Fozy <: Unsupervised end
MLJBase.fit(model::Fozy, verbosity, X) = minimum(X.matrix), nothing, nothing
MLJBase.transform(model::Fozy, fitresult, newbox) =
    fill(fitresult, nrows(newbox.matrix))
MLJBase.MLJModelInterface.reformat(model::Fozy, user_data) =
    (Box(MLJBase.matrix(user_data)),)
MLJBase.selectrows(model::Fozy, I, X...) = (Box(X[1].matrix[I,:]),)

## BABY SUPERVISED MODEL WITH AN UPDATE METHOD

mutable struct SomeModel <: Deterministic
    n::Int
end

function MLJModelInterface.fit(model::SomeModel,
                               verbosity,
                               A,
                               y)
    n = model.n
    cache = (A \ y)'  # all coefficients
    n_features = length(cache)

    # kill all but first n coefficients:
    fitresult = vcat(cache[1:n], fill(0.0, n_features - n))

    report = (n_features=n_features, )
    return fitresult, cache, report
end

function MLJModelInterface.update(model::SomeModel,
                                  verbosity,
                                  old_fitresult,
                                  old_cache,
                                  A,  # ignored in this case
                                  y)  # ignored in this case
    n = model.n
    cache = old_cache # coefficients already obtained in `fit`
    n_features = length(cache)

    # kill all but first n coefficients:
    fitresult = vcat(cache[1:n], fill(0.0, n_features - n))

    report = (n_features=n_features, )
    return fitresult, cache, report
end

function MLJModelInterface.predict(::SomeModel, fitresult, Xnew)
    Anew = MLJBase.matrix(Xnew)
    return Anew*fitresult
end

MLJModelInterface.reformat(model::SomeModel, X, y) = (MLJBase.matrix(X), y)
MLJModelInterface.selectrows(model::SomeModel, I, A, y) =
    (view(A, I, :), view(y, I))

@testset "overloading reformat(::Model, ...), selectrows(::Model, ...)" begin

    # dummy unsupervised model:
    model = Fozy()
    args = ((x1=[10, 30, 50], x2 = [20, 40, 60]),)
    data = MLJBase.MLJModelInterface.reformat(model, args...)
    @test data[1] isa Box && data[1].matrix == [10 20; 30 40; 50 60]
    @test selectrows(model, 2:3, data...)[1].matrix == [30 40; 50 60]
    @test fit(model, 1, data...)[1] == 10
    mach = machine(model, args...)
    @test_logs (:info, r"Training") fit!(mach, rows=2:3);
    @test transform(mach, (x1 = 1:4, x2 = 1:4)) == [30, 30, 30, 30]

    # supervised model with an update method:
    rng = StableRNGs.StableRNG(123)
    A = rand(rng, 8, 3)
    y = A*[1, 2, 3]
    X =  MLJBase.table(A)
    model = SomeModel(1)
    mach = machine(model, X, y)
    @test_mach_sequence fit!(mach, rows=1:4) [(:train, mach),]
    Xnew = selectrows(X, 1:4)
    @test predict(mach, Xnew) ≈ A[1:4,1]
    # mutate the model to trigger `update` call:
    model.n=3
    @test_mach_sequence fit!(mach, rows=1:4) [(:update, mach), ]
    @test predict(mach, Xnew) ≈ y[1:4]
    # change the rows to be sampled:
    @test_mach_sequence fit!(mach) [(:train, mach),]
    @test predict(mach, Xnew) ≈ y[1:4]
end

@testset "fit! for models with reformat front-end" begin
    X = (x1=ones(5), x2=2*ones(5))
    y = categorical(collect("abaaa"))

    clf = ConstantClassifier()
    clf = ConstantClassifier(testing=true)
    mach = machine(clf, X, y, cache=true)
    # first call to fit reformats data and resamples data:
    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(mach, rows=1:3, verbosity=0))
    @test mach.data == (MLJBase.matrix(X), y)
    @test mach.resampled_data[1] == mach.data[1][1:3,:]
    @test mach.resampled_data[2] == y[1:3]
    yhat = @test_logs (:info, r"reformatting X") predict_mode(mach, X)
    @test yhat == fill('a', 5)
    yhat = @test_logs (:info, "resampling X") predict_mode(mach, rows=1:2)
    @test yhat == fill('a', 2)
    # calling fit! with new `rows` triggers resampling but no
    # reformatting:
    @test_logs((:info, "resampling X, y"),
               fit!(mach, rows=1:2, verbosity=0))

end

@testset "Test serializable method of Supervised Machine" begin
    X, y = make_regression(100, 1)
    filename = "decisiontree.jls"
    mach = machine(DecisionTreeRegressor(), X, y)
    fit!(mach, verbosity=0)
    # Check serializable function
    smach = MLJBase.serializable(mach)
    @test smach.report == mach.report
    @test smach.fitresult == mach.fitresult
    @test_throws(ArgumentError, predict(smach))
    @test_logs (:warn, MLJBase.warn_serializable_mach(predict)) predict(smach, X)

    TestUtilities.generic_tests(mach, smach)
    # Check restore! function
    Serialization.serialize(filename, smach)
    smach = Serialization.deserialize(filename)
    MLJBase.restore!(smach)

    @test smach.state == 1
    @test MLJBase.predict(smach, X) == MLJBase.predict(mach, X)
    @test fitted_params(smach) isa NamedTuple
    @test report(smach) == report(mach)

    rm(filename)

    # End to end save and reload
    MLJBase.save(filename, mach)
    smach = machine(filename)
    @test smach.state == 1
    @test predict(smach, X) == predict(mach, X)

    rm(filename)
end

@testset "Test serializable method of Unsupervised Machine" begin
    X, _ = make_regression(100, 1)
    filename = "standardizer.jls"
    mach = machine(Standardizer(), X)
    fit!(mach, verbosity=0)

    MLJBase.save(filename, mach)
    smach = machine(filename)

    @test transform(mach, X) == transform(smach, X)
    @test_throws(ArgumentError, transform(smach))

    # warning on non-restored machine
    smach = deserialize(filename)
    @test_logs (:warn, MLJBase.warn_serializable_mach(transform)) transform(smach, X)

    rm(filename)
end

@testset "Test Misc functions used in `serializable`" begin
    X, y = make_regression(100, 1)
    mach = machine(DeterministicConstantRegressor(), X, y)
    fit!(mach, verbosity=0)
    # setreport! default
    @test mach.report isa NamedTuple
    MLJBase.setreport!(mach, "toto")
    @test mach.report == "toto"
end

end # module

true


text:
using PrettyTables

function Base.show(io::IO, br::BenchmarkResult)
    hb = Highlighter(
        (br,i,j) -> (j > 1 && maximum(@view(br.results[:, i])) == br.results[j-1,i]),
        foreground = :green
    );
    pretty_table(
        io, br.sizedresults, br.tests, crop = :none, highlighters = (hb,)
    )
end

if (Sys.ARCH === :aarch64) && Sys.isapple()
    nothing
else

using Colors, ColorSchemes, Gadfly
const COLORS = [RGB(0.0,0.0,0.0),RGB(1.0,0.0,0.0)]
# const COLORS = [RGB(0.0,0.0,0.0),RGB(0.0,1.0,0.0)]
# const COLORS = distinguishable_colors(14, pushfirst!(get(ColorSchemes.Paired_12, (0.5:11.5) ./ 12), RGB(0.0,0.0,0.0)))
for i ∈ 1:12 # 11 is number of tested libs - 2
    push!(COLORS, get(ColorSchemes.cyclic_mygbm_30_95_c78_n256_s25, i/12))
    # push!(COLORS, get(ColorSchemes.vikO, (i-0.5)/12))
end
# const COLOR_MAP = Dict{String,RGB{Float64}}()
# const COLOR_MAP = Dict{String,RGB{Colors.N0f8}}()
const COLOR_MAP64 = Dict{String,RGB{Float64}}()
function getcolor(s::String)
    get!(COLOR_MAP64, s) do
        COLORS[length(COLOR_MAP64) + 1]
    end
end
replace_and(str) = replace(str, '&' => "with")

function Gadfly.plot(br::BenchmarkResult)
    res = br.sizedresults.results
    sizes = br.sizedresults.sizes
    # sizes = Vector{eltype(brsizes)}(undef, length(res))
    tests = replace_and.(@view(br.tests[2:end]))
    colors = getcolor.(tests)
    addlabel = false

    maxxval, maxxind = findmax(sizes)
    maxxtick = 10cld(maxxval,10) + (addlabel ? 20 : 0)
    xt = 0:20:maxxtick
    maxres = maximum(res)
    maxtick = 10round(Int, 0.1maxres)
    yt = if iszero(maxtick)
        maxtick = 10round(0.1maxres)
        range(0, maxres, length = 20)
    elseif maxtick < 10
        0:1:maxtick
    elseif maxtick < 20
        0:2:maxtick
    elseif maxtick < 50
        0:5:maxtick
    else
        0:10:maxtick
    end
    p = Gadfly.plot(
        Gadfly.Guide.manual_color_key("Methods", tests, colors),
        Guide.xlabel("Size"), Guide.ylabel("GFLOPS"),
        Guide.xticks(ticks=collect(xt)), Guide.yticks(ticks=collect(yt))
    )
    for i ∈ eachindex(tests)
        push!(p, layer(x = sizes, y = res[i,:], Geom.line, Theme(default_color=colors[i])))
    end
    addlabel && push!(p, layer(x = fill(maxxtick - 10, length(tests)), y = res[:,maxxind], label=tests, Geom.label(position=:centered)))
    p
end

end
# using VegaLite, IndexedTables
# function plot(br::BenchmarkResult)
#     res = vec(br.sizedresults.results)
#     brsizes = br.sizedresults.sizes
#     sizes = Vector{eltype(brsizes)}(undef, length(res))
#     ntests = length(br.tests) - 1
#     for i ∈ 0:length(brsizes)-1
#         si = brsizes[i+1]
#         for j ∈ 1:ntests
#             sizes[j + i*ntests] = si
#         end
#     end
#     names = ["$(i > 9 ? string(i) : "0$i"). $test" for (i,test) ∈ enumerate(@view(br.tests[2:end]))]
#     tests = reduce(vcat, (names for _ ∈ eachindex(brsizes)))
#     t = table((GFLOPS = res, Size = sizes, Method = tests))
#     t |> @vlplot(
#         :line,
#         x = :Size,
#         y = :GFLOPS,
#         width = 900,
#         height = 600,
#         color={
#             :Method,
#             scale={scheme="category20"}
#         }
#     )
# end


text:
struct ODECompositeSolution{T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE} <: DiffEqBase.AbstractODESolution{T,N,uType}
  u::uType
  u_analytic::uType2
  errors::EType
  t::tType
  k::rateType
  prob::P
  alg::A
  interp::IType
  alg_choice::Vector{Int}
  dense::Bool
  tslocation::Int
  destats::DE
  retcode::Symbol
end
(sol::ODECompositeSolution)(t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(t,idxs,deriv,sol.prob.p,continuity)
(sol::ODECompositeSolution)(v,t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(v,t,idxs,deriv,sol.prob.p,continuity)

function DiffEqBase.build_solution(
        prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDDEProblem},
        alg::OrdinaryDiffEqCompositeAlgorithm,t,u;
        timeseries_errors=length(u)>2,
        dense=false,dense_errors=dense,
        calculate_error = true,
        k=[],
        du=[],
        interp = !isempty(du) ? HermiteInterpolation(t,u,du) : LinearInterpolation(t,u),
        alg_choice=[1], retcode = :Default, destats=DiffEqBase.DEStats(), kwargs...)

  T = eltype(eltype(u))
  if typeof(prob.u0) <: Tuple
    N = length((size(ArrayPartition(prob.u0))..., length(u)))
  else
    N = length((size(prob.u0)..., length(u)))
  end

  f = prob.f

  if DiffEqBase.has_analytic(f)
    if !(typeof(prob.u0) <: Tuple)
      u_analytic = Vector{typeof(prob.u0)}(undef, 0)
      errors = Dict{Symbol,real(eltype(prob.u0))}()
    else
      u_analytic = Vector{typeof(ArrayPartition(prob.u0))}(undef, 0)
      errors = Dict{Symbol,real(eltype(prob.u0[1]))}()
    end

    sol = ODECompositeSolution{T,N,typeof(u),typeof(u_analytic),typeof(errors),typeof(t),typeof(k),
                               typeof(prob),typeof(alg),typeof(interp),typeof(destats)}(u,u_analytic,errors,t,k,prob,alg,interp,alg_choice,dense,0,destats,retcode)
    if calculate_error
      DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
    end
    return sol
  else
    return ODECompositeSolution{T,N,typeof(u),typeof(nothing),typeof(nothing),typeof(t),typeof(k),
                                typeof(prob),typeof(alg),typeof(interp),typeof(destats)}(u,nothing,nothing,t,k,prob,alg,interp,alg_choice,dense,0,destats,retcode)
  end
end

function DiffEqBase.solution_new_retcode(sol::ODECompositeSolution{T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE},retcode) where {T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE}
  ODECompositeSolution{T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE}(
                     sol.u,sol.u_analytic,sol.errors,sol.t,sol.k,sol.prob,
                     sol.alg,sol.interp,sol.alg_choice,sol.dense,sol.tslocation,
                     sol.destats,retcode)
end

function DiffEqBase.solution_new_tslocation(sol::ODECompositeSolution{
   T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE},tslocation) where {T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE}
   ODECompositeSolution{T,N,uType,uType2,EType,tType,rateType,P,A,IType,DE}(
                      sol.u,sol.u_analytic,sol.errors,sol.t,sol.k,sol.prob,
                      sol.alg,sol.interp,sol.alg_choice,sol.dense,tslocation,
                      sol.destats,sol.retcode)
end

function DiffEqBase.sensitivity_solution(sol::ODECompositeSolution,u,t)
    T = eltype(eltype(u))
    N = length((size(sol.prob.u0)..., length(u)))
    interp = if !sol.dense
      sol.interp
    else
      DiffEqBase.SensitivityInterpolation(t,u)
    end

    ODECompositeSolution{T,N,typeof(u),typeof(sol.u_analytic),typeof(sol.errors),
                typeof(t),Nothing,typeof(sol.prob),typeof(sol.alg),
                typeof(interp),typeof(sol.destats)}(
                u,sol.u_analytic,sol.errors,t,nothing,sol.prob,
                sol.alg,interp,sol.alg_choice,
                sol.dense,sol.tslocation,
                sol.destats,sol.retcode)
end


text:
using Test
using CorpusLoaders
using CorpusLoaders: @NestedVector
using MultiResolutionIterators

@testset "NestedVector" begin
    @test @NestedVector(Int,1)==Vector{Int}
    @test @NestedVector(Char,2)==Vector{Vector{Char}}
    @test @NestedVector(String, 3)==Vector{Vector{Vector{String}}}
end

@testset "Document" begin
    @test collect(Document(["a", "b", "c"])) == ["a", "b", "c"]
    @test collect(Document("foos",["a", "b", "c"])) == ["a", "b", "c"]

    @testset "consolidate" begin
        cdoc = consolidate(Document("foos",["a", "b", "c"]))
        @test cdoc.title == "foos"
        @test cdoc[1] == "a"
        @test cdoc[2] == "b"
        @test cdoc[3] == "c"
        @test_throws BoundsError cdoc[4]
    end

    @testset "full_consolidate" begin
        cdocs = full_consolidate([Document("foos",["a", "b", "c"]), Document("bars",["A", "B"])])
        cdoc = cdocs[1]
        @test cdoc.title == "foos"
        @test cdoc[1] == "a"
        @test cdoc[2] == "b"
        @test cdoc[3] == "c"
        @test_throws BoundsError cdoc[4]
    end

end


text:
function lanczosIterations(A, stepsize, αs, βs, vs, reorthogonalize)

    iter = length(αs)
    v1 = vs[iter]
    v = vs[iter + 1]
    β = βs[iter]

    for i = iter .+ (1:stepsize)
        w = A*v
        α = real(dot(w,v))

        axpy!(-α, v, w)
        axpy!(-β, v1, w)

        # Orthogonlize (FixMe! consider a partial reorthogonal strategy here)
        if reorthogonalize
            for j = 1:i - 1
                axpy!(-dot(vs[j], w), vs[j], w)
            end
        end

        v1 = v
        v, β = qr!(w)

        push!(αs, α)
        push!(βs, β)
        push!(vs, v)
    end
    return αs, βs, vs
end

function _teig(A,
    nvals = 1;
    maxiter = size(A, 1),
    initvec = convert(Vector{eltype(A)}, randn(size(A, 1))),
    tolconv = sqrt(eps(real(eltype(A)))),
    reorthogonalize = true,
    stepsize = max(1, div(nvals, 10)),
    debug = false)

    # LinAlg.chksquare(A)

    ### initialize αs, βs and vs ###
    # run a single iteration of the Lanczos procedure. Note that this step handles promotion Lanczos vectors as well. They will in in general be different from `initvec`.
    v1 = A*initvec
    v0 = convert(typeof(v1), initvec)
    v0, nm0 = qr!(v0)
    rmul!(v1, inv(nm0)) # because we don't assume that initvec is normalized
    α = real(dot(v1, v0))
    axpy!(-α, v0, v1)
    v1, β = qr!(v1)
    αs = [α]
    βs = [β]
    vs = typeof(v1)[v0, v1]

    ### initial iterations ###
    lanczosIterations(A, nvals - 1, αs, βs, vs, reorthogonalize)

    ### Iteration count ###
    iter = nvals

    # values0 = sort!(abs(eigvals(SymTridiagonal(αs, βs[1:end - 1]))), rev = true)
    # values1 = values0

    z = zeros(eltype(αs), 1)

    ### run itertions ###
    # continue until convergence or the maximum number of iterations has been reached.
    converged = UInt[]
    while iter <= maxiter
        lanczosIterations(A, stepsize, αs, βs, vs, reorthogonalize)
        # values, vectors = LinearAlgebra.EigenSelfAdjoint.eigQL!(SymTridiagonal(copy(αs), βs[1:end - 1]), [zeros(length(αs) - 1); ones(eltype(αs), 1)]')
        values, vectors = eigen(SymTridiagonal(αs, βs[1:end-1]))

        ### Use criteria (9.18) of Peter Arbenz' lecture notes ###
        βj = βs[end]
        converged = Int[]
        for i = 1:length(αs)
            if βj*abs(vectors[end,i]) < abs(values[i])*tolconv
                push!(converged, i)
            end
        end

        debug && @show converged

        # if length(converged) >= nvals
        #     break
        # end

        # For now we hardcode convergence for largest positive value. This should be a choice later on
        if all([length(αs) - i + 1 in converged for i = 1:nvals])
            break
        end

        iter += stepsize
        debug && @show iter
    end
    if length(converged) < nvals
        error("no convergence")
    end

    S = SymTridiagonal(αs, βs[1:end-1])
    values, vectors = eigen(S)
    return values[converged], vectors[:, converged], S, vs, βs[end]
end

teig(A,
    nvals = 1;
    maxiter = size(A, 1),
    initvec = convert(Vector{eltype(A)}, randn(size(A, 1))),
    tolconv = 1e-12,
    reorthogonalize = true,
    stepsize = max(1, div(nvals, 10)),
    debug = false) = _teig(A, nvals, maxiter = maxiter, initvec = initvec,
        tolconv = tolconv, reorthogonalize = reorthogonalize, stepsize = stepsize, debug = debug)


text:
@testset "community" begin

    #Create community
    @test_nowarn begin
        m = @agent(
            0,

            [v,l,i]::Local,
            g::Global,
            ga::GlobalArray,

            Equation = v̇ = 34*dt ,
            UpdateLocal = l += 1,
            UpdateGlobal = begin
                g += Normal(1,2) + Normal(1,2) + Uniform(1,2)
                ga[3,2] += 3
            end,
            UpdateInteraction = i.i += 1,
            UpdateLocalInteraction = i.i += 1
        )

        m = compile(m)

        com = Community(m)
        com = Community(m,N=30)
    end

    #Get variables
    @test_nowarn begin
        m = @agent(
            2,

            [id1,id2,id3]::Identity,
            [l1,l2,l3]::Local,
            [g1,g2,g3]::Global,
            [ga1,ga2,ga3]::GlobalArray,
            [m1,m2]::Medium,

            Boundary = BoundaryFlat(2,Periodic(),Periodic())
        )

        m = compile(m)

        com = Community(m,N=10,mediumN=[10,10])

        com.declaredSymbols_
        com.t
        com.N
        com.id1
        com.id2
        com.id3
        com.l1
        com.l2
        com.l3
        com.g1
        com.g2
        com.g3
        com.ga1
        com.ga2
        com.ga3
        com.simulationBox
        com.radiusInteraction
    end

    #Get variables
    @test_nowarn begin
        m = @agent(
            2,

            [id1,id2,id3]::Identity,
            [l1,l2,l3]::Local,
            [g1,g2,g3]::Global,
            [ga1,ga2,ga3]::GlobalArray,
            [m1,m2,m3]::Medium,

            Boundary = BoundaryFlat(2,Periodic(),Periodic())
        )

        m = compile(m)

        com = Community(m,N=10,mediumN=[2,2])

        com.id1[1] = 1
        com.id2[1] = 2
        com.id3[1] = 3
        com.l1[1] = 1
        com.l2[1] = 2
        com.l3[1] = 3
        com.ga1 = ones(2,2)
        com.ga2 = 2*ones(2,2)
        com.ga3 = 3*ones(2,2)
        com.g1 = 1
        com.g2 = 2
        com.g3 = 3
        com.m1[1,1] = 1
        com.m2[1,1] = 2
        com.m3[1,1] = 3
        com.simulationBox = [0. 1.;0. 1.]
        com.radiusInteraction = [0.,1.]
        com.radiusInteraction = 1.

        if com.id1[1] != 1 error() end
        if com.id2[1] != 2 error() end
        if com.id3[1] != 3 error() end
        if com.l1[1] != 1. error() end
        if com.l2[1] != 2. error() end
        if com.l3[1] != 3. error() end
        if com.ga1[1,1] != 1. error() end
        if com.ga2[1,1] != 2. error() end
        if com.ga3[1,1] != 3. error() end
        if com.g1 != 1. error() end
        if com.g2 != 2. error() end
        if com.g3 != 3. error() end
        if com.m1[1,1] != 1. error() end
        if com.m2[1,1] != 2. error() end
        if com.m3[1,1] != 3. error() end

        com.id1 .= 1
        com.id2 .= 2
        com.id3 .= 3
        com.l1 .= 1
        com.l2 .= 2
        com.l3 .= 3
        com.ga1 .= 4
        com.ga2 .= 5
        com.ga3 .= 6
        com.m1 .= 1
        com.m2 .= 2
        com.m3 .= 3

        if com.id1 != 1*ones(com.N) error() end
        if com.id2 != 2*ones(com.N) error() end
        if com.id3 != 3*ones(com.N) error() end
        if com.l1 != 1. *ones(com.N) error() end
        if com.l2 != 2. *ones(com.N) error() end
        if com.l3 != 3. *ones(com.N) error() end
        if com.ga1 != 4. *ones(2,2) error() end
        if com.ga2 != 5. *ones(2,2) error() end
        if com.ga3 != 6. *ones(2,2) error() end
        if com.m1 != 1. *ones(2,2) error() end
        if com.m2 != 2. *ones(2,2) error() end
        if com.m3 != 3. *ones(2,2) error() end
    end

    @test_throws ErrorException begin
        m = @agent(0)
        m = compile(m)
        com = Community(m,N=10)

        com.id = 1
    end

    @test_throws MethodError begin
        m = @agent(0)
        m = compile(m)
        com = Community(m,N=10)

        com.id = [1]
    end

    @test_throws ErrorException begin
        m = @agent(0)
        m = compile(m)
        com = Community(m,N=10)

        com.id = 1
    end

    @test_throws MethodError begin
        m = @agent(0)
        m = compile(m)
        com = Community(m,N=10)

        com.id = [1]
    end

    @test_throws ErrorException begin
        m = @agent(0,ga::GlobalArray)
        m = compile(m)
        com = Community(m,N=10)

        com.ga = 1
    end

    #initialisers
    m = @agent 3
    mo = compile(m)
    fCom(x) = sqrt(sum(x.^2)) < 10

    @test_nowarn initialiseCommunityCompactHexagonal(mo,[-10 10;-10 10;-10 10],.5)
    @test_nowarn initialiseCommunityCompactHexagonal(mo,[-10 10;-10 10;-10 10],.5,fExtrude = fCom)
    @test_nowarn initialiseCommunityCompactCubic(mo,[-10 10;-10 10;-10 10],.5)
    @test_nowarn initialiseCommunityCompactCubic(mo,[-10 10;-10 10;-10 10],.5,fExtrude = fCom)

end


text:
using AbInitioSoftwareBase: parentdir, extension

@testset "`load` and `save`" begin
    dict = Dict(
        "a" => 1,
        "b" => 1.0e4,
        "c" => [1, 2, "3", "4.0"],
        "4" => true,
        "5.0" => "Hello",
    )
    @testset "Save to YAML" begin
        file = "test.yaml"
        save(file, dict)
        @test load(file) == dict
        file = "test.yml"
        save(file, dict)
        @test load(file) == dict
        save("mixed.YaMl", dict)
        @test load("mixed.YaMl") == dict
    end
    @testset "Save to JSON" begin
        file = "test.json"
        save(file, dict)
        @test load(file) == dict
        save("mixed.Json", dict)
        @test load("mixed.Json") == dict
    end
    @testset "Save to TOML" begin
        file = "test.toml"
        save(file, dict)
        @test load(file) == dict
        save("mixed.ToMl", dict)
        @test load("mixed.ToMl") == dict
    end
    @testset "Other extensions" begin
        @test_throws MethodError save("x.aml", dict)
        @test_throws MethodError save("x.jon", dict)
        @test_throws MethodError save("x.Tml", dict)
    end
    @testset "Format conversions" begin
        of_format("testconvert.json", "test.toml")
        @test load("testconvert.json") == load("test.toml")
        of_format("testconvert.toml", "test.json")
        @test load("testconvert.toml") == load("test.json")
    end
end

@testset "Test `parentdir`" begin
    @test parentdir("test.yaml") == parentdir("./test.yaml") == pwd()
end

@testset "Extensions" begin
    @test extension("a.b_c.d.yaml") == "yaml"
    @test extension("a.b_c.d.YaMl") == "yaml"
end


text:
# This file was generated, do not modify it. # hide
gdf = groupby(iris, :Species)
combine(gdf, :PetalLength => mean => :MPL, :PetalLength => std => :SPL)


text:
module FundamentosDMC

using Printf
using Parameters
using StaticArrays
using FastPow
using Statistics: mean

export System, Options
export Point2D, Point3D
export minimize!
export init_velocities
export printxyz
export md, md_isokinetic, md_berendsen, md_langevin, mc
export radial_distribution, velocity_distribution

include("./Points.jl")
include("./System.jl")
include("./Options.jl")

include("./image.jl")

include("./potential.jl")
include("./kinetic.jl")
include("./remove_drift.jl")
include("./forces.jl")

include("./minimize.jl")
include("./printxyz.jl")

include("./init_velocities.jl")
include("./md-simple.jl")
include("./md-isokinetic.jl")
include("./md-berendsen.jl")
include("./md-langevin.jl")
include("./mc.jl")

include("./radial_distribution.jl")
include("./velocity_distribution.jl")

end


text:
import ..UncertainValues: UncertainScalarPopulation

using RecipesBase

#@recipe f(::Type{UncertainScalarPopulation{T}}, x::UncertainScalarPopulation{T}) where {T} =
#    rand(x, 10000)

@recipe function f(p::UncertainScalarPopulation{T}) where T
    @series begin
        rand(p, 10000)
    end
end

@recipe function f(populations::Vector{UncertainScalarPopulation{T}}) where {T}

    for p in populations
        @series begin
            p
        end
    end
end