text:
using Test, MLToolkit
using MLToolkit.DistributionsX: u2gumbel
using Distributions: Dirichlet
using Statistics: mean
using Tracker: gradient, data
using FiniteDifferences: central_fdm
using StatsFuns: logit, logistic
using Flux: gpu

@testset "Gumbel" begin
    n_randtests = 5
    n_samples = 5_000
    d, n = 2, 10
    atol = 0.015

    @testset "GumbelSoftmax" begin
        for _ = 1:n_randtests
            # Vector `p`
            p = rand(Dirichlet(ones(d))) |> gpu
            gs = GumbelSoftmax(p)
            @test mean(gs) == p
            x = rand(gs, n_samples)
            @test vec(mean(x; dims=2)) ≈ mean(gs) atol=atol * d

            # Matrix `p`
            p = rand(Dirichlet(ones(d)), n) |> gpu
            gs = GumbelSoftmax(p)
            xs = [rand(gs) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(gs) atol=atol * d * n

            logpdf(gs, xs[1])
            @warn "Correctness of `logpdf(::GumbelSoftmax, x)` is not tested."
        end
    end

    @testset "GumbelSoftmax2D" begin
        for _ = 1:n_randtests
            # Vector `p`
            p1 = rand()
            gs = GumbelSoftmax2D(p1)
            @test mean(gs) == [p1, 1 - p1]
            x = rand(gs, n_samples)
            @test vec(mean(x; dims=2)) ≈ mean(gs) atol=atol * 2

            # Matrix `p`
            p1 = rand(n) |> gpu
            gs = GumbelSoftmax2D(p1)
            @test mean(gs) == [p1'; 1 .- p1']
            xs = [rand(gs) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(gs) atol=atol * 2 * n
        end
    end

    @testset "track(u2gumbel, u)" begin
        for _ = 1:n_randtests
            u = rand() / 2 + 1 / 4
            f(x) = sum(u2gumbel(x))
            g_ad, = data.(gradient(f, u))
            g_fd = central_fdm(5, 1)(f, u)
            @test g_ad ≈ g_fd
        end
    end

    @testset "GumbelBernoulli" begin
        for _ = 1:n_randtests
            p = rand(d, n) |> gpu
            gb = GumbelBernoulli(p)
            @test mean(gb) == p
            xs = [rand(gb) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(gb) atol=atol * d * n

            # Consistency between two versions of log density
            x = rand(gb)
            @test logpdf(gb, x) ≈ logpdfCoV(gb, x)
        end
    end

    @testset "GumbelBernoulliLogit" begin
        for _ = 1:n_randtests
            p = rand(d, n) |> gpu
            gbl = GumbelBernoulliLogit(logit.(p))
            @test mean(gbl) ≈ p
            xs = [rand(gbl) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(gbl) atol=atol * d * n
        end
    end
end


text:
using Bokeh
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end

start = now()
x = map(d -> start + Dates.Hour(d), 1:(30 * 24))
println(typeof(x))
# x: Array{DateTime,1}
day = 1000 * 3600 * 24
y = 0.5sin(int(x) * pi / day) + cos(int(x) * pi / (day * 7))
# note x needs to be an array of type DateTime or Date
plot(x, y)
showplot("dates_xy.html")

start = Date(1980, 1, 1)
x = map(d -> start + Dates.Day(d), 1:(365 * 30))
println(typeof(x))
# x: Array{Date,1}
y = 0.5sin(int(x) * pi / 365) + cos(int(x) * pi / (365 * 7))

plot(y, x)
showplot("dates_yx.html")

ext: jl

max_stars_repo_name: Laeeth/dbokeh

max_issues_repo_name: Laeeth/dbokeh

max_forks_repo_name: Laeeth/dbokeh


text:
mutable struct Record
    header::Union{Missing, String}
    frequencies::Union{Missing, AbstractMatrix}
    sequence::Union{Missing, BioSequence}
end

function Record()
    return Record(missing, missing, missing)
end

function Record(data::Vector{UInt8})
    return convert(Record, data)
end

function Base.convert(::Type{Record}, data::Vector{UInt8})
    record = Record()
    stream = TranscodingStreams.NoopStream(IOBuffer(data))
    cs, linenum, found = readrecord!(stream, record, (1, 1))
    if !found || !allspace(stream)
        throw(ArgumentError("invalid record"))
    end
    return record
end

function Record(str::AbstractString)
    return Record(Vector{UInt8}(str))
end

function Base.convert(::Type{Record}, str::AbstractString)
    return Record(Vector{UInt8}(str))
end

function Base.copy(record::Record)
    return Record(
        record.header,
        record.frequencies,
        record.sequence
        )
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "      header: ", header(record))
        println(io, " frequencies: ", frequencies(record))
        println(io, "    sequence: ", sequence(record))
    else
        println(io, " <not filled>")
    end
end

function initialize!(record::Record)
    record.header = missing
    record.frequencies = missing
    record.sequence = missing
    return record
end

function header(record::Record)
    return record.header
end

function frequencies(record::Record)
    return record.frequencies
end

function sequence(record::Record)
    return record.sequence
end

function BioGenerics.isfilled(record::Record)
    return !ismissing(record.header) || !ismissing(record.frequencies) || !ismissing(record.sequence)
end

function hasmissing(record::Record)
    return ismissing(record.header) || ismissing(record.frequencies) || ismissing(record.sequence)
end


text:
"""
A Julia interface module for SciPy

- [GitHub repository: AtsushiSakai/SciPy.jl: Julia interface for SciPy](https://github.com/AtsushiSakai/SciPy.jl)
- [Stable Doc](https://atsushisakai.github.io/SciPy.jl/stable)
- [Dev Doc](https://atsushisakai.github.io/SciPy.jl/dev)

"""
module SciPy

using PyCall
using InteractiveUtils

export print_configulations
export cluster, constants, fft, integrate, interpolate, io, linalg, ndimage, odr
export optimize, signal, sparse, spatial, stats, special
# raw PyObject modules
export pystats, pyoptimize, pyinterpolate, pyspatial

const scipy = PyNULL()

include("util.jl")

"""
scipy.cluster module

- [Clustering package (scipy.cluster) Reference Guide](https://docs.scipy.org/doc/scipy/reference/cluster.html)

# Examples

You can do k-means clustering using this module:

```julia-repl
julia> features  = [[ 1.9  2.3];
                    [ 1.5 2.5];
                    [ 0.8 0.6];
                    [ 0.4 1.8];
                    [ 0.1 0.1];
                    [ 0.2 1.8];
                    [ 2.0 0.5];
                    [ 0.3 1.5];
                    [ 1.0 1.0]]
9×2 Array{Float64,2}:
 1.9  2.3
 1.5  2.5
 0.8  0.6
 0.4  1.8
 0.1  0.1
 0.2  1.8
 2.0  0.5
 0.3  1.5
 1.0  1.0

julia> whitened = SciPy.cluster.vq.whiten(features)
9×2 Array{Float64,2}:
 2.7396    2.91001
 2.16284   3.16306
 1.15351   0.759134
 0.576757  2.2774
 0.144189  0.126522
 0.288379  2.2774
 2.88379   0.632612
 0.432568  1.89784
 1.44189   1.26522

julia> SciPy.cluster.vq.kmeans(whitened, [whitened[1,:] whitened[3,:]] )
([1.1174670798453024 1.8345740800894272; 2.8837860125040065 0.6326117517549749], 1.073399
3090584457)

```

"""
const cluster = PyNULL()

"""
scipy.constants module

- [Constants (scipy.constants) Reference Guide](https://docs.scipy.org/doc/scipy/reference/constants.html)

# Examples

You can access each constants:

```julia-repl
julia> SciPy.constants.golden
1.618033988749895

julia> SciPy.constants.physical_constants["electron mass"]
(9.10938356e-31, "kg", 1.1e-38)

julia> SciPy.constants.convert_temperature([-40, 40.0], "Celsius", "Kelvin")

2-element Array{Float64,1}:
 233.14999999999998
 313.15

```
"""
const constants = PyNULL()

"""
scipy.fft module

- [Discrete Fourier transforms (scipy.fft) Reference Guide](https://docs.scipy.org/doc/scipy/reference/fft.html)

# Examples

You can use FFT (Fast Fourier Transform):

```julia-repl

julia> SciPy.fft.fft(exp.(π/8 * collect(1:8)))
8-element Array{Complex{Float64},1}:
   68.17385416403044 - 0.0im
   1.408601300061675 + 31.248171041435185im
 -10.268363617931092 + 15.207165888808841im
 -12.695027025520982 + 6.493878653648949im
 -13.216494113330363 - 0.0im
 -12.695027025520982 - 6.493878653648949im
 -10.268363617931092 - 15.207165888808841im
   1.408601300061675 - 31.248171041435185im

```
"""
const fft = PyNULL()

"""
scipy.integrate module

- [Integration and ODEs (scipy.integrate) Reference Guide](https://docs.scipy.org/doc/scipy/reference/integrate.html)

# Examples

You can compute a definite integral:

```julia-repl
julia> f(x) = x^2
f (generic function with 1 method)

julia> SciPy.integrate.quad(f, 0, 4)
(21.333333333333336, 2.368475785867001e-13)

```
"""
const integrate = PyNULL()

const pyinterpolate = PyNULL()
include("interpolate.jl")

"""
scipy.io module

- [Input and output (scipy.io) Reference Guide](https://docs.scipy.org/doc/scipy/reference/io.html)

# Examples

You can save a MATLAB-style .mat file:

```julia-repl
julia> mdic = Dict([("a", 100), ("label", "experiment")])
Dict{String,Any} with 2 entries:
  "label" => "experiment"
  "a"     => 100

julia> SciPy.io.savemat("sample_data.mat", mdic)
```
"""
const io = PyNULL()

"""
scipy.linalg module

- [Linear algebra (scipy.linalg) Reference Guide](https://docs.scipy.org/doc/scipy/reference/linalg.html)

# Examples

You can compute the matrix exponential using Pade approximation:

```julia-repl
julia> SciPy.linalg.expm(zeros((2,2)))
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
```
"""
const linalg = PyNULL()

"""
scipy.ndimage module

- [Multi-dimensional image processing (scipy.ndimage) Reference Guide](https://docs.scipy.org/doc/scipy/reference/ndimage.html)

# Examples

You can compute a multidimensional convolution:

```julia-repl
julia> k = [[1 1 1];[1 1 0];[1 0 0]]
3×3 Array{Int64,2}:
 1  1  1
 1  1  0
 1  0  0

julia> a = [[1 2 0 0];
            [5 3 0 4];
            [0 0 0 7];
            [9 3 0 0]]
4×4 Array{Int64,2}:
 1  2  0  0
 5  3  0  4
 0  0  0  7
 9  3  0  0

julia> SciPy.ndimage.convolve(a, k, mode="constant", cval=0.0)
4×4 Array{Int64,2}:
 11  10   7   4
 10   3  11  11
 15  12  14   7
 12   3   7   0
```
"""
const ndimage = PyNULL()

"""
scipy.odr module

- [Orthogonal distance regression (scipy.odr) Reference Guide](https://docs.scipy.org/doc/scipy/reference/odr.html)

# Examples

You can calculate orthogonal distance regression with an exponential model:

```julia-repl
julia> x = collect(0.0:5.0);

julia> y = -10.0 .+ exp.(0.5*x);

julia> data = SciPy.odr.Data(x, y)
PyObject <scipy.odr.odrpack.Data object at 0x7fe5fda4ccc0>

julia> data = SciPy.odr.Data(x, y);

julia> odr_obj = SciPy.odr.ODR(data, SciPy.odr.exponential);

julia> output = odr_obj.run();

julia> println(output.beta)
[-10.0, 0.5]
```
"""
const odr = PyNULL()

const pyoptimize = PyNULL()
include("optimize.jl")

"""
scipy.signal module

- [Signal processing (scipy.signal) Reference Guide](https://docs.scipy.org/doc/scipy/reference/signal.html)

# Examples

You can compute the Kaiser parameter beta, given the attenuation a:

```julia-repl
julia> SciPy.signal.kaiser_beta(65)
6.20426

```
"""
const signal = PyNULL()

"""
scipy.sparse module

- [Sparse matrices (scipy.sparse) Reference Guide](https://docs.scipy.org/doc/scipy/reference/sparse.html)

# Examples

You can do sparse matrix calculation:

```julia-repl
julia> A = SciPy.sparse.csc_matrix([[1.0 0.0];
                                   [1.0 2.0]]);

julia> Ainv = SciPy.sparse.linalg.inv(A);

julia> A.dot(Ainv).todense()
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

```
"""
const sparse = PyNULL()

const pyspatial = PyNULL()
include("spatial.jl")

"""
scipy.special module

- [Special functions (scipy.special) Reference Guide](https://docs.scipy.org/doc/scipy/reference/special.html)

"""
const special = PyNULL()

const pystats = PyNULL()
include("stats.jl")

"""
Module initialization function
"""
function __init__()
    copy!(scipy, pyimport_conda("scipy", "scipy"))
    copy!(cluster, pyimport_conda("scipy.cluster", "scipy"))
    copy!(constants, pyimport_conda("scipy.constants", "scipy"))
    copy!(fft, pyimport_conda("scipy.fft", "scipy"))
    copy!(integrate, pyimport_conda("scipy.integrate", "scipy"))
    copy!(io, pyimport_conda("scipy.io", "scipy"))
    copy!(linalg, pyimport_conda("scipy.linalg", "scipy"))
    copy!(ndimage, pyimport_conda("scipy.ndimage", "scipy"))
    copy!(odr, pyimport_conda("scipy.odr", "scipy"))
    copy!(signal, pyimport_conda("scipy.signal", "scipy"))
    copy!(sparse, pyimport_conda("scipy.sparse", "scipy"))
    copy!(special, pyimport_conda("scipy.special", "scipy"))

    copy!(pyspatial, pyimport_conda("scipy.spatial", "scipy"))
    copy!(pystats, pyimport_conda("scipy.stats", "scipy"))
    copy!(pyoptimize, pyimport_conda("scipy.optimize", "scipy"))
    copy!(pyinterpolate, pyimport_conda("scipy.interpolate", "scipy"))
end

"""
Print configulations:
- Julia version
- Python version
- Python path
- scipy version
"""
function print_configulations()
    println("==== SciPy.jl configulations ====")
    println("-- Julia settings --")
    versioninfo()
    println("-- PyCall.jl settings --")
    @show PyCall.pyversion
    @show PyCall.libpython
    @show PyCall.pyprogramname
    @show PyCall.conda
    println("-- SciPy settings --")
    @show scipy.__version__
    @show scipy.version.full_version
    @show scipy.version.git_revision
    return nothing
end

end


text:
using Spectra
using Test
using Statistics
using Random

# based on the example of the smooth function...

# the x axis
x = collect(0:0.01:10)

# a fake signal: perfect y
y_perfect = x.^2

# fake noisy signal
y = y_perfect + randn(size(x,1))*2

y_sv = smooth(x,y,method="savgol",window_length=9,polyorder=3)
y_whit = smooth(x,y,method="whittaker", Lambda=10.0^1)

ese_noise = sum((y .- y_perfect).^2)
ese_sg = sum((y_sv.-y_perfect).^2)
ese_whit = sum((y_whit.-y_perfect).^2)

@test ese_sg < ese_noise
@test ese_whit < ese_noise


text:
@testset "Crossover" begin

    v4n(graph::CompGraph, want) = v4n(vertices(graph), want)
    v4n(vs, want) = vs[findfirst(v -> want == name(v), vs)]

    teststrat() = NaiveGAflux.default_crossoverswap_strategy(v -> 1)

    @testset "CrossoverSwap" begin
        using NaiveGAflux: crossoverswap!, separablefrom
        using Functors: fmap

        iv(np) = denseinputvertex("$np.in", 3)
        dv(in, outsize, name) = fluxvertex(name, Dense(nout(in), outsize), in)

        @testset "Simple Dense swap" begin

            function g(sizes, np)
                vi = iv(np)
                vo = foldl(enumerate(sizes); init=vi) do vin, (i,s)::Tuple
                    dv(vin, s, "$np.dv$i")
                end
                return CompGraph(vi, vo)
            end

            @testset "Single vertex swap" begin
                ga = g(4:6, "a")
                gb = g(1:3, "b")

                @test crossoverswap!(vertices(ga)[3], vertices(gb)[3]; strategy=teststrat) == (true,true)

                @test name.(vertices(ga)) == ["a.in", "a.dv1", "b.dv2", "a.dv3"]
                @test name.(vertices(gb)) == ["b.in", "b.dv1", "a.dv2", "b.dv3"]

                @test nout.(vertices(ga)) == [3,4,2,6]
                @test nout.(vertices(gb)) == [3,1,5,3]

                @test size(ga(ones(3, 2))) == (6, 2)
                @test size(gb(ones(3, 2))) == (3, 2)
            end

            @testset "Path swap" begin
                ga = g(4:8, "a")
                gb = g(1:5, "b")

                @test crossoverswap!(vertices(ga)[2], vertices(ga)[5], vertices(gb)[3], vertices(gb)[4]; strategy=teststrat) == (true,true)

                @test nout.(vertices(ga)) == [3,2,3,8]
                @test nout.(vertices(gb)) == [3,1,4,5,6,7,4,5]

                @test size(ga(ones(3, 2))) == (8, 2)
                @test size(gb(ones(3, 2))) == (5, 2)
            end

            @testset "Swap same graph is noop" begin
                ga = g(4:7, "a")

                gb = fmap(x -> x isa String ? "b" * x[2:end] : x, ga)

                indata = reshape(collect(Float32, 1:3*2),3,:)
                orgout = ga(indata)
                @test orgout == gb(indata)

                @test name.(vertices(ga)) == ["a.in", "a.dv1", "a.dv2", "a.dv3", "a.dv4"]
                @test name.(vertices(gb)) == ["b.in", "b.dv1", "b.dv2", "b.dv3", "b.dv4"]

                vsa = vertices(ga)
                vsb = vertices(gb)
                @test crossoverswap!(vsa[2], vsa[4], vsb[2], vsb[4]) == (true,true)

                @test orgout == ga(indata) == gb(indata)
            end
        end

        @testset "Swap add and conc with conc larger" begin
            function g(sizes, np, cfun)
                vi = iv(np)
                vs = map((i, s) -> dv(vi, s, "$np.dv$i"), eachindex(sizes[1:end-1]), sizes[1:end-1])
                vm = cfun("$np.merge", vs...)
                return CompGraph(vi, dv(vm, sizes[end], "$np.dv$(length(sizes))"))
            end

            ga = g(3:7, "a", concat)
            gb = g(fill(3, 4), "b", (name, vs...) -> +(name >> vs[1], vs[2:end]...))

            @test crossoverswap!(vertices(ga)[end-1], vertices(gb)[end-1]; strategy=teststrat) == (true, true)

            @test nout.(vertices(ga)) == [3, 6, 6, 6, 6, 6, 7]
            @test nout.(vertices(gb)) == [3, 3, 3, 3, 9, 3]

            @test size(ga(ones(3, 2))) == (7, 2)
            @test size(gb(ones(3, 2))) == (3, 2)
        end

        @testset "Swap add and conc with add larger" begin
            function g(sizes, np, cfun)
                vi = iv(np)
                vs = map((i, s) -> dv(vi, s, "$np.dv$i"), eachindex(sizes[1:end-1]), sizes[1:end-1])
                vm = cfun("$np.merge", vs...)
                return CompGraph(vi, dv(vm, sizes[end], "$np.dv$(length(sizes))"))
            end

            ga = g(3:6, "a", concat)
            gb = g(fill(3, 6), "b", (name, vs...) -> +(name >> vs[1], vs[2:end]...))

            @test crossoverswap!(vertices(ga)[end-1], vertices(gb)[end-1]; strategy=teststrat) == (true, true)

            @test nout.(vertices(ga)) == [3, 5, 5, 5, 5, 6]
            @test nout.(vertices(gb)) == [3, 3, 3, 3, 3, 3, 15, 3]

            @test size(ga(ones(3, 2))) == (6, 2)
            @test size(gb(ones(3, 2))) == (3, 2)
        end

        @testset "Swap conc and dense" begin
            function g(sizes, np, cfun)
                vi = iv(np)
                vs = map((i, s) -> dv(vi, s, "$np.dv$i"), eachindex(sizes[1:end]), sizes[1:end])
                vm = cfun("$np.merge", vs...)
                dvn = dv(vm, 2, "$np.dvn")
                dvo = dv(dvn, 4, "$np.dvo")
                return CompGraph(vi, dvo)
            end

            ga = g(2:4, "a", concat)
            gb = g(2:4, "b", concat)

            crossoverswap!(vertices(ga)[end-1], vertices(gb)[end-2]; mergefun = v -> concat("extramerge", v), strategy=teststrat)

            @test name.(vertices(ga)) == ["a.in", "a.dv1", "a.dv2", "a.dv3", "a.merge", "b.merge", "a.dvo"]
            @test name.(vertices(gb)) == ["b.in", "b.dv1", "b.dv2", "b.dv3", "extramerge", "a.dvn", "b.dvn", "b.dvo"]

            @test nout.(vertices(ga)) == [3, 2, 3, 4, 9, 9, 4]
            @test nout.(vertices(gb)) == [3, 2, 3, 4, 9, 2, 2, 4]

            @test size(ga(ones(3, 2))) == (4, 2)
            @test size(gb(ones(3, 2))) == (4, 2)
        end

        @testset "Swap residual" begin
            function g(sizes, np, withinv)
                vi = iv(np)
                dv0 = dv(vi, 3, "$np.dv0")
                dv1 = dv(dv0, sizes[1], "$np.dv1")
                vo = foldl(enumerate(sizes[2:end]); init=dv1) do vin, (i,s)::Tuple
                    dv(vin, s, "$np.dv$i")
                end
                dvn = dv(vo, sizes[1], "$np.dvn")
                vn = withinv ? invariantvertex(identity, dvn; traitdecoration= named("$np.inv")) : dvn
                add = "$np.add" >> dv1 + vn
                dvo = dv(add, 4, "$np.dvo")
                return CompGraph(vi, dvo)
            end

            @testset "Invariant before add" begin
                ga = g(4:6, "a", true)
                gb = g(5:6, "b", true)

                vsa = vertices(ga)
                vsb = vertices(gb)

                @test crossoverswap!(vsa[4], vsa[end-2], vsb[4], vsb[end-2]; strategy=teststrat) == (true, true)
                @test size(ga(ones(Float32, 3, 2))) == (4,2)
                @test size(gb(ones(Float32, 3, 2))) == (4,2)
            end
        end

        @testset "Dummy shall pass size validation" begin
            import NaiveNASflux: named, validated, traitconf
            # Evilness from this case: We want to swap out tv1 to dv2, but both tv1 and dv2 is input
            # to a elementwise addition. NaiveNASlib handles this, but if the dummy vertex does not
            # adjust its output size NaiveNASlib will think that the graph is size inconsistent as
            # it does not know that we intend to remove the dummy vertex.

            function g(np, sz)
                vi = iv(np)
                dv1 = dv(vi, sz, "$np.dv1")
                tv1 = invariantvertex(identity, dv1; traitdecoration=named("$np.tv1") ∘ validated())
                tv2 = invariantvertex(identity, tv1; traitdecoration=named("$np.tv2") ∘ validated())
                dv2 = dv(tv2, nout(tv2), "$np.dv2")
                add = traitconf(named("$np.add") ∘ validated()) >> dv2 + tv1
                dv3 = dv(add, 2, "$np.dv3")
                return CompGraph(vi, dv3), tv1, dv2
            end

            ga, vin1, vout1 = g("a", 2)
            gb, vin2, vout2 = g("b", 4)

            @test crossoverswap!(vin1, vout1, vin2, vout2; strategy=teststrat) == (true, true)
        end

        @testset "Find swappable path" begin

            @testset "Residual on input" begin
                function g(np)
                    vi = iv(np)
                    dv1 = dv(vi, 4, "$np.dv1")
                    dv2 = dv(dv1, nout(vi), "$np.dv2")
                    add1 = "$np.add1" >> vi + dv2
                    dv3 = dv(add1, 4, "$np.dv3")
                    dv4 = dv(dv3, 3, "$np.dv4")
                    return CompGraph(vi, dv4)
                end

                ga = g("a")

                indata = randn(MersenneTwister(1), 3, 2)
                outa = ga(indata)

                vsa = vertices(ga)
                aswap = separablefrom(v4n(ga, "a.dv3"), ga.inputs)
                @test name.(aswap) == ["a.dv3", "a.add1"]
            end

            @testset "Dual branched graph with optional extra connection" begin
                function g(np, bconnect = false)
                    vi = iv(np)
                    dv1 = dv(vi, 4, "$np.dv1")

                    dva1 = dv(dv1, 3, "$np.dva1")
                    dvaa1 = dv(dva1, 5, "$np.dvaa1")
                    dvaa2 = dv(dvaa1, 4, "$np.dvaa2")
                    dvab1 = dv(dva1, 2, "$np.dvab1")
                    dvab2 = dv(dvab1, 4, "$np.dvab2")
                    add_aa_bb = "$np.add_aa_bb" >> dvaa2 + dvab2

                    dvb1 = dv(dv1, 5, "$np.dvb1")
                    vba1 = bconnect ? concat("$np.conc_dvb1_dvab2", dvb1, dvab2) : dv(dvb1, nout(dvb1) + nout(dvab2), "$np.dvba1")
                    dvbb1 = dv(dv1, 4, "$np.dvbb1")
                    conc_ba_bb = concat("$np.conc_ba_bb", vba1, dvbb1)

                    conc_a_b = concat("$np.conc_a_b", add_aa_bb, conc_ba_bb)
                    return CompGraph(vi, dv(conc_a_b, 4, "$np.out"))
                end

                ga = g("a")

                indata = randn(MersenneTwister(1), 3, 2)
                outa = ga(indata)

                vsa = vertices(ga)
                aswap = separablefrom(v4n(ga, "a.add_aa_bb"), ga.inputs)
                @test name.(vsa) == name.(vertices(ga))
                @test name.(aswap) == ["a.add_aa_bb", "a.dva1"]

                @test ga(indata) == outa

                gb = g("b", true)
                outb = gb(indata)
                vsb = vertices(gb)
                @test name.(separablefrom(v4n(gb, "b.add_aa_bb"), gb.inputs)) == ["b.add_aa_bb"]
                @test name.(vsb) == name.(vertices(gb))

                @test gb(indata) == outb

                bswap = separablefrom(v4n(gb, "b.dvbb1"), gb.inputs)
                @test name.(vsb) == name.(vertices(gb))
                @test name.(bswap) == ["b.dvbb1"]

                @test gb(indata) == outb

                @test crossoverswap!(aswap[end], aswap[1], bswap[end], bswap[1]; strategy=teststrat) == (true, true)

                @test name.(vertices(ga)) == ["a.in", "a.dv1", "b.dvbb1", "a.dvb1", "a.dvba1", "a.dvbb1", "a.conc_ba_bb", "a.conc_a_b", "a.out"]
                @test size(ga(ones(3,2))) == (4, 2)

                @test name.(vertices(gb)) == ["b.in", "b.dv1", "b.dva1", "b.dvaa1", "b.dvaa2", "b.dvab1", "b.dvab2", "b.add_aa_bb", "b.dvb1", "b.conc_dvb1_dvab2", "a.dva1", "a.dvaa1", "a.dvaa2", "a.dvab1", "a.dvab2", "a.add_aa_bb", "b.conc_ba_bb", "b.conc_a_b", "b.out"]
                @test size(gb(ones(3,2))) == (4, 2)
            end
        end

        @testset "Swapping preserves edge order" begin
            import NaiveGAflux: stripoutedges!, stripinedges!, addoutedges!, addinedges!
            function g(np)
                vi = iv(np)
                dv1 = dv(vi, 2, "$np.dv1")
                dv2 = dv(vi, 3, "$np.dv2")
                dv3 = dv(vi, 4, "$np.dv3")
                ca1 = concat("$np.ca1", dv1, dv2, dv1, dv3, dv1)
                dva1 = dv(ca1, 4, "$np.dva1")
                dva2 = dv(dva1, 3, "$np.dva2")
                dvb1 = dv(dv2, 3, "$np.dvb1")
                c2 = concat("$np.c2", dvb1, dva1, dvb1, dva1)
                c3 = concat("$np.c3", dva1, dva1, dvb1)
                c4 = concat("$np.c4", c2,c3,dva2)
                out = dv(c4, 4, "$np.out")
                return CompGraph(vi, out)
            end
            @testset "Strip and add out edges" begin
                gg = g("a")
                dva1 = v4n(gg, "a.dva1")

                indata = randn(MersenneTwister(0), 3,2)
                graphout_before = gg(indata)

                expected1 = name.(outputs(dva1))
                expected2 = mapreduce(vo -> name.(inputs(vo)), vcat, unique(outputs(dva1)))

                dummy = stripoutedges!(dva1)
                @test outputs(dva1) == []
                addoutedges!(dva1, dummy, teststrat)

                actual1 = name.(outputs(dva1))
                actual2 = mapreduce(vo -> name.(inputs(vo)), vcat, unique(outputs(dva1)))

                @test actual1 == expected1
                @test actual2 == expected2

                @test gg(indata) == graphout_before
            end

            @testset "Strip and add in edges" begin
                gg = g("a")
                ca1 = v4n(gg, "a.ca1")

                indata = randn(MersenneTwister(0), 3,2)
                graphout_before = gg(indata)

                expected = name.(inputs(ca1))

                ins = stripinedges!(ca1)
                @test mapreduce(inputs, vcat, inputs(ca1);init=[]) == []
                addinedges!(ca1, ins, teststrat)

                actual = name.(inputs(ca1))

                @test actual == expected

                @test gg(indata) == graphout_before
            end

            g_org = g("a")
            g_new = deepcopy(g_org)

            indata = randn(MersenneTwister(0), 3, 2)
            out_org = g_org(indata)

            swappable_new = separablefrom(v4n(g_new, "a.dva1"))
            @test name.(swappable_new) == ["a.dva1", "a.ca1"]
            @test name.(vertices(g_org)) == name.(vertices(g_new))

            @test g_org(indata) == g_new(indata) == out_org

            vs_org = vertices(g_org)
            nouts_org = nout.(vs_org)

            swappable_org = separablefrom(v4n(g_org, "a.dva1"))
            @test crossoverswap!(swappable_org[end], swappable_org[1], swappable_new[end], swappable_new[1]; strategy=teststrat) == (true, true)

            @test name.(vertices(g_org)) == name.(vertices(g_new)) == name.(vs_org)
            @test nout.(vertices(g_org)) == nout.(vertices(g_new)) == nouts_org
            @test g_org(indata) == g_new(indata) == out_org
        end

        @testset "Partial success" begin
            idv(in, outsize, name) = fluxvertex(name, Dense(nout(in), outsize), in; traitfun=t -> NamedTrait(name, Immutable()))

            function g(np, mergesize, mergeop)
                vi = iv(np)
                dv1 = dv(vi, mergesize, "$np.dv1")
                dv2 = idv(dv1, mergesize, "$np.dv2")
                m1 = mergeop("$np.m1", dv1, vi, dv2)
                dv3 = dv(m1, 5, "$np.dv3")
                return CompGraph(vi, dv3)
            end

            @testset "Fail first" begin
                ga = g("a", 5, concat)
                gb = g("b", 3, (vname, vs...) -> +(vname >> vs[1], vs[2:end]...))

                @test @test_logs (:warn, "Failed to align sizes for vertices a.dv2 and b.dv2.dummy for crossover. Attempt aborted!") crossoverswap!(vertices(ga)[end-1], vertices(gb)[end-1]; strategy=teststrat) == (false, true)

                @test name.(vertices(gb)) == ["b.in", "b.dv1", "b.dv2", "a.m1", "b.dv3"]
                @test size(gb(ones(3,2))) == (5,2)
            end

            @testset "Fail second" begin
                ga = g("a", 3, (vname, vs...) -> +(vname >> vs[1], vs[2:end]...))
                gb = g("b", 5, concat)

                @test @test_logs (:warn, "Failed to align sizes for vertices b.dv2 and a.dv2.dummy for crossover. Attempt aborted!") crossoverswap!(vertices(ga)[end-1], vertices(gb)[end-1]; strategy=teststrat) == (true, false)

                @test name.(vertices(ga)) == ["a.in", "a.dv1", "a.dv2", "b.m1", "a.dv3"]
                @test size(ga(ones(3,2))) == (5,2)
            end
        end

        @testset "Swap invariant layers" begin
            bv(in, name) = fluxvertex(name, BatchNorm(nout(in)), in)
            pv(in, name) = fluxvertex(name, MaxPool((3,3); pad=(1,1)), in)
            cv(in, outsize, name) = fluxvertex(name, Conv((3,3), nout(in) => outsize; pad=(1,1)), in)

            @testset "Linear graph" begin
                function g(np)
                    vi = iv(np)
                    cv1 = cv(vi, 4, "$np.cv1")
                    bv1 = bv(cv1, "$np.bv1")
                    cv2 = cv(bv1, 5, "$np.cv2")
                    pv1 = pv(cv2, "$np.pv1")
                    cv3 = cv(pv1, 2, "$np.cv3")
                    return CompGraph(vi, cv3)
                end

                @testset "Swap $v1 to $v2" for (v1,v2) in (
                    ("bv1", "bv1"),
                    ("pv1", "pv1"),
                    ("bv1", "pv1")
                    )

                    vs(g, np) = v4n(g, "$np.$v1"), v4n(g, "$np.$v2")

                    ga = g("a")
                    gb = g("b")

                    @test crossoverswap!(vs(ga, "a")..., vs(gb, "b")...; strategy=teststrat) == (true, true)

                    @test "a.$v1" ∉ name.(vertices(ga))
                    @test "b.$v1" ∈ name.(vertices(ga))

                    @test "a.$v2" ∉ name.(vertices(ga))
                    @test "b.$v2" ∈ name.(vertices(ga))

                    @test "b.$v1" ∉ name.(vertices(gb))
                    @test "a.$v1" ∈ name.(vertices(gb))

                    @test "b.$v2" ∉ name.(vertices(gb))
                    @test "a.$v2" ∈ name.(vertices(gb))

                    @test size(ga(ones(Float32, 4,4,3,2))) == (2,2,2,2)
                    @test size(gb(ones(Float32, 4,4,3,2))) == (2,2,2,2)
                end

                @testset "Self swap is noop" begin
                    g_org = g("a")

                    bn = layer(v4n(g_org, "a.bv1"))
                    bn.γ .= randn(Float32, nout(bn))

                    indata = randn(Float32, 4,4,3,2)
                    out_org = g_org(indata)

                    g_new = deepcopy(g_org)

                    @test crossoverswap!(v4n(g_org, "a.bv1"), v4n(g_new, "a.bv1"); strategy=teststrat) == (true, true)

                    @test g_org(indata) == out_org
                    @test g_new(indata) == out_org
                end
            end

            @testset "Swap double transparent parlayer" begin

                function g(sizes, np)
                    vi = iv(np)
                    cv1 = cv(vi, sizes[1], "$np.cv1")
                    bv1 = bv(cv1, "$np.bv1")
                    bv2 = bv(bv1, "$np.bv2")
                    cv2 = cv(bv2, sizes[2], "$np.cv2")
                    cv3 = cv(cv2, 4, "$np.cv3")
                    return CompGraph(vi, cv3)
                end

                ga = g((3,2), "a")
                gb = g((6,2), "b")

                strat = () -> NaiveGAflux.default_crossoverswap_strategy(v -> 1)
                @test crossoverswap!(vertices(ga)[3], vertices(ga)[4], vertices(gb)[3], vertices(gb)[4]; strategy = strat) == (true,true)

                @test name.(vertices(ga)) == ["a.in", "a.cv1", "b.bv1", "b.bv2", "a.cv2", "a.cv3"]
                @test name.(vertices(gb)) == ["b.in", "b.cv1", "a.bv1", "a.bv2", "b.cv2", "b.cv3"]

                @test size(ga(ones(Float32, 4,4,3,1))) == (4,4,4,1)
                @test size(gb(ones(Float32, 4,4,3,1))) == (4,4,4,1)
            end
        end
    end

    @testset "VertexCrossover" begin
        import NaiveGAflux: GlobalPool, sameactdims, default_pairgen

        iv(np) = denseinputvertex("$np.in", 3)
        dv(in, outsize, name) = fluxvertex(name, Dense(nout(in), outsize), in)
        gv(in, name) = invariantvertex(GlobalPool(MaxPool), in; traitdecoration = named(name))
        cv(in, outsize, name) = fluxvertex(name, Conv((3,3), nout(in) => outsize; pad=(1,1)), in)

        @testset "Linear graph" begin
            function g(np, cvs, dvs)
                vi = iv(np)
                cvn = foldl(enumerate(cvs); init=vi) do vin, (i,s)::Tuple
                    cv(vin, s, "$np.cv$i")
                end
                pv1 = gv(cvn, "$np.pv1")
                dvn = foldl(enumerate(dvs); init=pv1) do vin, (i,s)::Tuple
                    dv(vin, s, "$np.dv$i")
                end
                return CompGraph(vi, dvn)
            end

            @testset "regraph" begin
                import NaiveGAflux: regraph

                g_org = g("test", 3:4, 5:6)

                v = vertices(g_org)[4]

                g_new = regraph(v, 1, 1)
                @test name.(vertices(g_org)) == name.(vertices(g_new))

                @test_throws AssertionError regraph(v, 2 ,1)
                @test_throws AssertionError regraph(v, 1, 2)
            end

            @testset "Generate pairs" begin
                ga = g("a", (4,5), (3,2))

                vsa = vertices(ga)[2:end]
                @test sameactdims.(vsa[end], vsa) == sameactdims.(vsa[end-1], vsa) == Bool[0, 0, 0, 1, 1]
                @test sameactdims.(vsa[1], vsa) == sameactdims.(vsa[2], vsa) == Bool[1, 1, 0, 0, 0]
                @test sameactdims.(vsa[3], vsa) == Bool[0,0,1,0,0]

                @testset "Match same size" begin
                    gb = g("b", (4,5), (3,2))
                    vsb = vertices(gb)[2:end]

                    @test default_pairgen(vsa, vsb; ind1=1) == (1, 1)
                    @test default_pairgen(vsa, vsb; ind1=3) == (3, 3)
                    @test default_pairgen(vsa, vsb; ind1=length(vsa)+1) === nothing
                end

                @testset "Match larger" begin
                    gb = g("b", 2:7, (8,))
                    vsb = vertices(gb)[2:end]

                    @test default_pairgen(vsa, vsb; ind1=1) == (1, 2)
                    @test default_pairgen(vsa, vsb; ind1=3) == (3, 7) # Only global pool matches with global pool
                    @test default_pairgen(vsa, vsb; ind1=5) == (5, 8)
                end

                @testset "Match smaller" begin
                    gb = g("b", (2,), (3,))
                    vsb = vertices(gb)[2:end]

                    @test default_pairgen(vsa, vsb; ind1=1) == (1, 1)
                    @test default_pairgen(vsa, vsb; ind1=2) == (2, 1)
                    @test default_pairgen(vsa, vsb; ind1=3) == (3, 2) # Only global pool matches with global pool
                    @test default_pairgen(vsa, vsb; ind1=4) == (4, 3)
                    @test default_pairgen(vsa, vsb; ind1=5) == (5, 3)
                end
            end

            @testset "Test swap single" begin
                import NaiveGAflux: crossover, default_pairgen, crossoverswap
                ga = g("a", 3:4, 5:6)
                gb = g("b", 2:4, 5:7)

                pairgen_outer(vs1,vs2; ind1=1) = ind1 > 3 ? nothing : default_pairgen(vs1, vs2; ind1=min(2*ind1, length(vs1)))
                pairgen_inner(vs1,vs2) = default_pairgen(vs1, vs2; ind1 = length(vs1)-1)

                crossfun(args...) = crossoverswap(args...;pairgen = pairgen_inner, strategy=teststrat)

                anames = name.(vertices(ga))
                bnames = name.(vertices(gb))

                ga_new, gb_new = crossover(ga, gb; pairgen=pairgen_outer, crossoverfun=crossfun)

                # Originals not impacted
                @test anames == name.(vertices(ga))
                @test bnames == name.(vertices(gb))

                @test name.(vertices(ga_new)) == ["a.in", "b.cv2", "b.cv3", "a.pv1", "b.dv2", "b.dv3"]
                @test name.(vertices(gb_new)) == ["b.in", "b.cv1", "a.cv1", "a.cv2", "b.pv1", "b.dv1", "a.dv1", "a.dv2"]

                @test size(ga_new(ones(Float32, 4,4,3,2))) == (7, 2)
                @test size(gb_new(ones(Float32, 4,4,3,2))) == (7, 2)
            end

            @testset "VertexCrossover" begin
                import NaiveGAflux: default_pairgen

                pairgen_outer(vs1,vs2; ind1=1) = default_pairgen(vs1, vs2; ind1=2ind1)
                pairgen_inner(vs1, vs2) = default_pairgen(vs1, vs2; ind1=length(vs1))

                c = VertexCrossover(CrossoverSwap(;pairgen=pairgen_inner, strategy=teststrat); pairgen=pairgen_outer)

                ga = g("a", 3:6, 7:8)
                gb = g("b", 2:3, 4:10)

                anames = name.(vertices(ga))
                bnames = name.(vertices(gb))

                ga_new, gb_new = c((ga,gb))

                # Originals not impacted
                @test anames == name.(vertices(ga))
                @test bnames == name.(vertices(gb))

                @test name.(vertices(ga_new)) == ["a.in", "a.cv1", "b.cv2", "a.cv3", "a.cv4", "a.pv1", "b.dv6", "a.dv2"]
                @test name.(vertices(gb_new)) == ["b.in", "b.cv1", "a.cv2", "b.pv1", "b.dv1", "b.dv2", "b.dv3", "b.dv4", "b.dv5", "a.dv1", "b.dv7"]

                NaiveNASflux.forcemutation(ga_new)
                NaiveNASflux.forcemutation(gb_new)

                lnout(v) = nout(layer(v))
                lnin(v) = nin(layer(v))

                @test nout(v4n(ga_new, "b.cv2")) == lnout(v4n(ga_new, "b.cv2"))
                @test nin(v4n(ga_new, "a.cv3")) == lnin(v4n(ga_new, "a.cv3"))

                @test nout(v4n(ga_new, "b.dv6")) == lnout(v4n(ga_new, "b.dv6"))
                @test nin(v4n(ga_new, "a.dv2")) == lnin(v4n(ga_new, "a.dv2"))

                @test nout(v4n(gb_new, "a.cv2")) == lnout(v4n(gb_new, "a.cv2"))
                @test nin(v4n(gb_new, "b.dv1")) == lnin(v4n(gb_new, "b.dv1"))

                @test nout(v4n(gb_new, "a.dv1")) == lnout(v4n(gb_new, "a.dv1"))
                @test nin(v4n(gb_new, "b.dv7")) == lnin(v4n(gb_new, "b.dv7"))

                @test size(ga_new(ones(Float32, 4,4,3,2))) == (8,2)
                @test size(gb_new(ones(Float32, 4,4,3,2))) == (10,2)
            end
        end

        @testset "Branchy graph" begin
            function forks(vbase, vf, vs, np)
                fin = map(enumerate(vs)) do (i, s)
                    vf(vbase, s, "$np$i")
                end

                return map(enumerate(fin)) do (i, fstart)
                    foldl(enumerate(vs[1:end-i]); init=fstart) do vin, (i,s)::Tuple
                        vf(vin, s, "$(name(fstart)).$i")
                    end
                end
            end

            function g(np, cvs, dvs)
                vi = iv(np)
                cvbase = cv(vi, 3, "$np.cvbbase")

                cvfa = forks(cvbase, cv, cvs, "$np.cvfa")
                m1 = concat("$np.m1", cvfa...)

                cvfb = map(enumerate(forks(m1, cv, cvs, "$np.cvfb"))) do (i, vin)
                    cv(vin, minimum(cvs), "$(name(vin)).out")
                end
                m2 = +("$np.m2" >> cvfb[1], cvfb[2:end]...)

                pv1 = gv(m2, "$np.pv1")
                dvn = foldl(enumerate(dvs); init=pv1) do vin, (i,s)::Tuple
                    dv(vin, s, "$np.dv$i")
                end
                return CompGraph(vi, dvn)
            end

            @testset "Crossover with MutationProbability" begin
                import NaiveGAflux: default_pairgen

                pairgen_inner(vs1, vs2) = default_pairgen(vs1, vs2; ind1=length(vs1))

                c = VertexCrossover(MutationProbability(CrossoverSwap(;pairgen=pairgen_inner, strategy=teststrat), Probability(0.2, MockRng([0.3, 0.1, 0.3]))))

                ga = g("a", 3:6, 7:8)
                gb = g("b", 2:3, 4:10)

                anames = name.(vertices(ga))
                bnames = name.(vertices(gb))

                ga_new, gb_new = c((ga,gb))

                # Originals not impacted
                @test anames == name.(vertices(ga))
                @test bnames == name.(vertices(gb))

                @test unique(first.(split.(name.(vertices(ga_new)), '.'))) == ["a" ,"b"]
                @test unique(first.(split.(name.(vertices(gb_new)), '.'))) == ["b" ,"a"]

                @test size(ga_new(ones(Float32, 4,4,3,2))) == (8,2)
                @test size(gb_new(ones(Float32, 4,4,3,2))) == (10,2)
            end
        end

        @testset "Crossover conv size mismatch" begin
            pv(in, name) = fluxvertex(name, MaxPool((2,2); stride=(2,2)), in)
            function g(np)
                vi = iv(np)
                cv1 = cv(vi, 4, "$np.cv1")
                pv1 = pv(cv1, "$np.pv1")
                cv2 = cv(pv1, 2, "$np.cv2")

                cva1 = cv(cv2, 3, "$np.cva1")
                cva2 = cv(cva1, 2, "$np.cva2")

                cvb1 = cv(cv2, 4, "$np.cvb1")
                cvb2 = cv(cvb1, 2, "$np.cvb2")

                mv = "$np.mv" >> cva2 + cvb2
                cv3 = cv(mv, 3, "$np.cv3")
                pv2 = pv(cv3, "$np.pv2")
                cv4 = cv(pv2, 4, "$np.cv4")
                cv5 = cv(cv4, 3, "$np.cv5")

                return CompGraph(vi, cv5)
            end

            @testset "Swap branch for graph" begin

                ga = g("a")
                gb = g("b")

                pairgen_outer(vs1, vs2; ind1=1) = ind1 == 1 ? (5, 9) : nothing

                # Anonymous type to avoid name collisions with other tests (e.g. structs named DummyRng, DummyRng1, etc..)
                selectfirst = nr -> 1:nr
                Random.randn(rng::typeof(selectfirst), nr) = selectfirst(nr)
                function pairgen_inner(vs1, vs2)
                    # Test that the testcase is working :)
                    @test name.(vs1) == ["a.cva1", "a.cva2"]
                    @test name(first(vs2)) == "b.cv1"
                    @test name(last(vs2)) == "b.cv3"
                    return NaiveGAflux.default_inputs_pairgen(vs1, vs2, 10; ind1=2, rng=selectfirst)
                end

                c = VertexCrossover(CrossoverSwap(;pairgen=pairgen_inner, strategy=teststrat); pairgen=pairgen_outer)

                ga_new, gb_new = c((ga,gb))

                @test name.(vertices(ga_new)) == ["a.in", "a.cv1", "a.pv1", "a.cv2", "a.cva1", "b.cv2", "b.cva1", "b.cva2", "b.cvb1", "b.cvb2", "b.mv", "b.cv3", "a.cvb1", "a.cvb2", "a.mv", "a.cv3", "a.pv2", "a.cv4", "a.cv5"]
                @test name.(vertices(gb_new)) ==  ["b.in", "b.cv1", "b.pv1", "a.cva2", "b.pv2", "b.cv4", "b.cv5"]

                @test size(ga(ones(Float32, 10,10,3,1))) == size(ga_new(ones(Float32, 10,10,3,1)))
                @test size(gb(ones(Float32, 10,10,3,1))) == size(gb_new(ones(Float32, 10,10,3,1)))
            end

            @testset "Swap branch with single vertex" begin

                ga = g("a")
                gb = g("b")

                pairgen_outer(vs1, vs2; ind1=1) = ind1 == 1 ? (5, 11) : nothing

                # Anonymous type to avoid name collisions with other tests (e.g. structs named DummyRng, DummyRng1, etc..)
                selectfirst = nr -> 1:nr
                Random.randn(rng::typeof(selectfirst), nr) = selectfirst(nr)
                function pairgen_inner(vs1, vs2)
                    # Test that the testcase is working :)
                    @test name.(vs1) == ["a.cva1", "a.cva2"]
                    @test name(first(vs2)) == "b.cv1"
                    @test name(last(vs2)) == "b.cv4"
                    return NaiveGAflux.default_inputs_pairgen(vs1, vs2, 10; ind1=2, rng=selectfirst)
                end
                c = VertexCrossover(CrossoverSwap(;pairgen=pairgen_inner, strategy=teststrat); pairgen=pairgen_outer)

                ga_new, gb_new = c((ga,gb))

                @test name.(vertices(ga_new)) == ["a.in", "a.cv1", "a.pv1", "a.cv2", "a.cva1", "b.cv4", "a.cvb1", "a.cvb2", "a.mv", "a.cv3", "a.pv2", "a.cv4", "a.cv5"]
                @test name.(vertices(gb_new)) == ["b.in", "b.cv1", "b.pv1", "b.cv2", "b.cva1", "b.cva2", "b.cvb1", "b.cvb2", "b.mv", "b.cv3", "b.pv2", "a.cva2", "b.cv5"]

                @test size(ga(ones(Float32, 10,10,3,1))) == size(ga_new(ones(Float32, 10,10,3,1)))
                @test size(gb(ones(Float32, 10,10,3,1))) == size(gb_new(ones(Float32, 10,10,3,1)))
            end
        end
    end

    @testset "OptimizerCrossover" begin
        using NaiveGAflux.Flux.Optimise

        prts(o) = typeof(o)
        prts(o::Optimiser) = "$(typeof(o))$(prts.(Tuple(o.os)))"

        @testset "Swap optimizers $(prts(o1)) and $(prts(o2))" for (o1, o2) in (
            (ADAM(), Momentum()),
            (Optimiser(Descent(), WeightDecay()), Optimiser(Momentum(), Nesterov())),
            )
            oc = OptimizerCrossover()
            ooc = OptimizerCrossover(oc)
            @test prts.(oc((o1,o2))) == prts.(ooc((o1,o2))) == prts.((o2, o1))
            @test prts.(oc((o2,o1))) == prts.(ooc((o2,o1))) == prts.((o1, o2))
        end

        @testset "Don't swap shielded" begin
            o1 = ShieldedOpt(Descent())
            o2 = ShieldedOpt(Momentum())
            @test OptimizerCrossover()((o1,o2)) == (o1,o2)
        end

        @testset "Cardinality difference" begin

            @testset "Single opt vs Optimiser" begin
                oc = OptimizerCrossover()
                @test prts.(oc((Descent(), Optimiser(Momentum(), WeightDecay())))) == prts.((Momentum(), Optimiser(Descent(), WeightDecay())))
            end

            @testset "Different size Optimisers" begin
                oc = OptimizerCrossover()
                o1 = Optimiser(Descent(), WeightDecay(), Momentum())
                o2 = Optimiser(ADAM(), ADAMW(), NADAM(), RADAM())

                o1n,o2n = oc((o1,o2))

                @test prts(o1n) == prts(Optimiser(ADAM(), ADAMW(), NADAM()))
                @test prts(o2n) == prts(Optimiser(Descent(), WeightDecay(), Momentum(), RADAM()))
            end
        end

        @testset "LogMutation and MutationProbability" begin
            mplm(c) = MutationProbability(LogMutation(((o1,o2)::Tuple) -> "Crossover between $(prts(o1)) and $(prts(o2))", c), Probability(0.2, MockRng([0.3, 0.1, 0.3])))
            oc = OptimizerCrossover() |> mplm |> OptimizerCrossover

            o1 = Optimiser(Descent(), WeightDecay(), Momentum())
            o2 = Optimiser(ADAM(), ADAGrad(), AdaMax())

            o1n,o2n = @test_logs (:info, "Crossover between WeightDecay and ADAGrad") oc((o1,o2))

            @test typeof.(o1n.os) == [Descent, ADAGrad, Momentum]
            @test typeof.(o2n.os) == [ADAM, WeightDecay, AdaMax]
        end

        @testset "Learningrate crossover" begin
            import NaiveGAflux: learningrate
            @testset "Single opt" begin
                oc = LearningRateCrossover()
                o1,o2 = oc((Descent(0.1), Momentum(0.2)))

                @test typeof(o1) == Descent
                @test o1.eta == 0.2

                @test typeof(o2) == Momentum
                @test o2.eta == 0.1
            end

            @testset "Shielded opt" begin
                oc = LearningRateCrossover()
                o1,o2 = oc((ShieldedOpt(Descent(0.1)), Momentum(0.2)))

                @test typeof(o1) == ShieldedOpt{Descent}
                @test o1.opt.eta == 0.1

                @test typeof(o2) == Momentum
                @test o2.eta == 0.2
            end

            @testset "Optimiser" begin
                oc = LearningRateCrossover()
                o1 = Optimiser(Descent(0.1), Momentum(0.2), WeightDecay(0.1))
                o2 = Optimiser(ADAM(0.3), RADAM(0.4), NADAM(0.5), Nesterov(0.6))

                o1n,o2n = oc((o1,o2))

                @test prts(o1n) == prts(o1)
                @test prts(o2n) == prts(o2)

                @test learningrate.(o1n.os[1:end-1]) == [0.3, 0.4]
                @test learningrate.(o2n.os) == [0.1, 0.2, 0.5, 0.6]

            end
        end
    end
end


text:
module TestExamples1DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d")

@testset "MHD Multicomponent" begin

  @trixi_testset "elixir_mhdmulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [1.69266396e-01, 3.04958381e-01, 0.00000000e+00, 3.45735220e-01, 6.60825891e-17, 3.17662584e-01, 0.00000000e+00, 2.18373660e-02, 4.36747319e-02, 8.73494638e-02],
      linf = [5.08655909e-01, 1.25176699e+00, 0.00000000e+00, 1.99579851e+00, 1.11022302e-16, 1.68777938e+00, 0.00000000e+00, 1.00173010e-01, 2.00346020e-01, 4.00692040e-01])
  end

  @trixi_testset "elixir_mhdmulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
      l2   = [1.66364836e-01, 2.87812904e-01, 0.00000000e+00, 3.17156609e-01, 6.60825891e-17, 3.16727982e-01, 0.00000000e+00, 1.87387435e-02, 3.74774869e-02, 7.49549738e-02],
      linf = [5.25859519e-01, 1.11130058e+00, 0.00000000e+00, 1.46798145e+00, 1.11022302e-16, 1.86689645e+00, 0.00000000e+00, 7.63433974e-02, 1.52686795e-01, 3.05373590e-01])
  end

  @trixi_testset "elixir_mhdmulti_eoc.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_eoc.jl"),
      l2   = [1.70884468e-05, 3.33778117e-04, 3.33778117e-04, 3.60816702e-05, 4.13046273e-17, 3.33548882e-04, 3.33548882e-04, 1.30233539e-05, 2.60467077e-05, 5.20934154e-05],
      linf = [2.34173129e-05, 1.23807753e-03, 1.23807753e-03, 1.07512068e-04, 1.11022302e-16, 1.22960686e-03, 1.22960686e-03, 2.01157231e-05, 4.02314461e-05, 8.04628923e-05])
    end

  @trixi_testset "elixir_mhdmulti_briowu_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_briowu_shock_tube.jl"),
      l2   = [1.93221943e-01, 3.56599349e-01, 0.00000000e+00, 3.63446684e-01, 6.33435897e-16, 3.66051482e-01, 0.00000000e+00, 5.65287156e-02, 1.13057431e-01],
      linf = [4.33927958e-01, 1.09948372e+00, 0.00000000e+00, 1.03943772e+00, 2.99760217e-15, 1.50960600e+00, 0.00000000e+00, 1.96228180e-01, 3.92456360e-01])
    end

end

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Private functions and macros.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                  Constants
################################################################################

# Crayon used to reset all the styling.
const _reset_crayon = Crayon(reset = true)

################################################################################
#                             Auxiliary functions
################################################################################

"""
    _crop_str(str, crop_size)

Return a cropped string of `str` with size `crop_size`. Notice that if the last
characted before the crop does not fit due to its width, then blank spaces are
added.

"""
function _crop_str(str, crop_size)
    # If the crop_size is large than the screen size, then just return the
    # current string.
    crop_size ≥ textwidth(str) && return str

    # Process every character of the string.
    cstr  = ""
    lcstr = 0

    for c in str
        sc = textwidth(c)

        # Check if we can add the character without passing the crop size.
        if lcstr + sc ≤ crop_size
            cstr  *= c
            lcstr += sc

        # Fill the remaining spaces with blank characters.
        else
            Δ      = crop_size - lcstr
            cstr  *= " "^Δ
            break
        end
    end

    return cstr
end

"""
    _str_aligned(data::AbstractString, alignment::Symbol, field_size::Integer)

This function returns the string `data` with alignment `alignment` in a field
with size `field_size`. `alignment` can be `:l` or `:L` for left alignment, `:c`
or `:C` for center alignment, or `:r` or `:R` for right alignment. It defaults
to `:r` if `alignment` is any other symbol.

If the string is larger than `field_size`, then it will be cropped and `⋯` will
be added as the last character.

"""
function _str_aligned(data::AbstractString, alignment::Symbol,
                      field_size::Integer)

    lstr = textwidth(data)
    Δ    = field_size - lstr

    # If the length is larger than the field, then we will crop the string.
    if Δ < 0
        data  = _crop_str(data, lstr + Δ - 1)
        data *= "…"

        # In this case, the string has the same size of the field.
        Δ = 0
    end

    # TODO: If Δ is 0, can we just return the string?

    if alignment == :l || alignment == :L
        return data * " "^Δ
    elseif alignment == :c || alignment == :C
        left  = div(Δ,2)
        right = Δ-left
        return " "^left * data * " "^right
    else
        return " "^Δ * data
    end
end

"""
    _str_line_breaks(str::AbstractString, autowrap::Bool = false, width::Int = 0)

Split the string `str` into substring, each one meaning one new line. If
`autowrap` is `true`, then the text will be wrapped so that it fits the column
with the width `width`.

"""
function _str_line_breaks(str::AbstractString, autowrap::Bool = false, width::Int = 0)
    # Check for errors.
    autowrap && (width <= 0) &&
    error("If `autowrap` is true, then the width must not be positive.")

    # Get the tokens for each line.
    tokens_raw = _str_escaped.(split(str, '\n'))

    # If the user wants to auto wrap the text, then we must check if
    # the tokens must be modified.
    if autowrap
        tokens = String[]

        for token in tokens_raw
            sub_tokens = String[]
            length_tok = length(token)

            # Get the list of valid indices to handle UTF-8 strings. In this
            # case, the n-th character of the string can be accessed by
            # `token[tok_ids[n]]`.
            tok_ids = collect(eachindex(token))

            if length_tok > width
                # First, let's analyze from the beginning of the token up to the
                # field width.
                #
                # k₀ is the character that will start the sub-token.
                # k₁ is the character that will end the sub-token.
                k₀ = 1
                k₁ = k₀ + width - 1

                while k₀ <= length_tok

                    # Check if the remaining string fit in the available space.
                    if k₁ == length_tok
                        push!(sub_tokens, token[tok_ids[k₀:k₁]])

                    else
                        # If the remaining string does not fit into the
                        # available space, then we search for spaces to crop.
                        Δ = 0
                        for k = k₁:-1:k₀
                            if token[tok_ids[k]] == ' '
                                # If a space is found, then select `k₁` as this
                                # character and use `Δ` to remove it when
                                # printing, so that we hide the space.
                                k₁ = k
                                Δ  = 1

                                break
                            end
                        end

                        push!(sub_tokens, token[tok_ids[k₀:k₁-Δ]])
                    end

                    # Move to the next analysis window.
                    k₀ = k₁+1
                    k₁ = clamp(k₀ + width - 1, 0, length_tok)
                end
                push!(tokens, sub_tokens...)
            else
                push!(tokens, token)
            end
        end

        return tokens
    else
        return tokens_raw
    end
end

################################################################################
#                              Printing Functions
################################################################################

"""
    _draw_continuation_row(screen, io, tf, text_crayon, border_crayon, cols_width, vlines)

Draw the continuation row when the table has filled the vertical space
available. This function prints in each column the character `⋮` centered.

"""
function _draw_continuation_row(screen, io, tf, text_crayon, border_crayon,
                                cols_width, vlines)

    num_cols = length(cols_width)

    0 ∈ vlines && _p!(screen, io, border_crayon, tf.column)

    @inbounds for j = 1:num_cols
        data_ij_str = _str_aligned("⋮", :c, cols_width[j] + 2)
        _p!(screen, io, text_crayon, data_ij_str)

        flp = j == num_cols

        _pc!(j ∈ vlines, screen, io, border_crayon, tf.column, " ", flp)
        _eol(screen) && break
    end

    _nl!(screen, io)

    return nothing
end

"""
    _draw_line!(screen, io, left, intersection, right, row, border_crayon, cols_width)

Draw a vertical line in `io` using the information in `screen`.

"""
function _draw_line!(screen, io, left, intersection, right, row, border_crayon,
                     cols_width, vlines)

    num_cols = length(cols_width)

    0 ∈ vlines && _p!(screen, io, border_crayon, left)

    @inbounds for i = 1:num_cols
        # Check the alignment and print.
        _p!(screen, io, border_crayon, row^(cols_width[i]+2)) && break

        i != num_cols &&
            _pc!(i ∈ vlines, screen, io, border_crayon, intersection, row)
    end

    _pc!(num_cols ∈ vlines, screen, io, border_crayon, right, row, true)
    _nl!(screen, io)
end

"""
    _eol(screen)

Return `true` if the cursor is at the end of line or `false` otherwise.

"""
_eol(screen) = (screen.size[2] > 0) && (screen.col >= screen.size[2])

"""
    _nl!(screen, io)

Add a new line into `io` using the screen information in `screen`.

"""
function _nl!(screen, io)
    # Store the largest column that was printed.
    screen.max_col < screen.col && (screen.max_col = screen.col)
    screen.row += 1
    screen.col  = 0
    println(io)
end

"""
    _p!(screen, io, crayon, str, final_line_print = false)

Print `str` into `io` using the Crayon `crayon` with the screen information in
`screen`. The parameter `final_line_print` must be set to `true` if this is the
last string that will be printed in the line. This is necessary for the
algorithm to select whether or not to include the continuation character.

"""
function _p!(screen, io, crayon, str, final_line_print = false)
    # Get the size of the string.
    #
    # TODO: We might reduce the number of allocations by avoiding calling
    # `textwidth`, since we have the size of all fields.
    lstr = textwidth(str)

    # `sapp` is a string to be appended to `str`. This is used to add `⋯` if the
    # text must be wrapped. Notice that `lapp` is the length of `sapp`.
    sapp = ""
    lapp = 0

    # When printing a line, we must verify if the screen has bounds. This is
    # done by looking if the horizontal size is positive.
    @inbounds if screen.size[2] > 0

        # If we are at the end of the line, then just return.
        _eol(screen) && return true

        Δ = screen.size[2] - (lstr + screen.col)

        # Check if we can print the entire string.
        if Δ <= 0
            # If we cannot, then create a wrapped string considering how many
            # columns are left.
            if lstr + Δ - 2 > 0
                # Here we crop the string considering the available space
                # reserving 2 characters for the line continuation indicator.
                str  = _crop_str(str, lstr + Δ - 2)
                lstr = lstr + Δ - 2
                sapp = " ⋯"
                lapp = 2
            elseif screen.size[2] - screen.col == 2
                # If there are only 2 characters left, then we must only print
                # " ⋯".
                str  = ""
                lstr = 0
                sapp = " ⋯"
                lapp = 2
            elseif screen.size[2] - screen.col == 1
                # If there are only 1 character left, then we must only print
                # "⋯".
                str  = ""
                lstr = 0
                sapp = "⋯"
                lapp = 1
            else
                # This should never be reached.
                @error("Internal error!")
                return true
            end
        elseif !final_line_print
            # Here we must verify if this is the final printing on this line. If
            # it is, then we should just check if the entire string fits on the
            # available size. Otherwise, we must see if, after printing the
            # current string, we will have more than 1 space left. If not, then
            # we just add the continuation character sequence.

            if Δ == 1
                str   = lstr > 1 ? _crop_str(str, lstr - 1) : ""
                lstr -= 1
                sapp  = " ⋯"
                lapp  = 2
            end
        end
    end

    # Print the with correct formating.
    if screen.has_color
        print(io, crayon, str, _reset_crayon, sapp)
    else
        print(io, str, sapp)
    end

    # Update the current columns.
    screen.col += lstr + lapp

    # Return if we reached the end of line.
    return _eol(screen)
end

"""
    _pc!(cond, screen, io, crayon, str_true, str_false, final_line_print = false)

If `cond == true` then print `str_true`. Otherwise, print `str_false`. Those
strings will be printed into `io` using the Crayon `crayon` with the screen
information in `screen`. The parameter `final_line_print` must be set to `true`
if this is the last string that will be printed in the line. This is necessary
for the algorithm to select whether or not to include the continuation
character.

"""
function _pc!(cond, screen, io, crayon, str_true, str_false,
              final_line_print = false)
    if cond
        return _p!(screen, io, crayon, str_true, final_line_print)
    else
        return _p!(screen, io, crayon, str_false, final_line_print)
    end
end


text:
# Welded beam design problem
# from Kalyanmoy Deb, "An effcient constraint handling method for genetic algorithms", 2000

using Evolutionary, Test, Random

beam(x) = 1.10471*x[2]*x[1]^2+0.04811*x[3]*x[4]*(14+x[2])
lx = [0.125, 0.1, 0.1, 0.125]
ux = [5.0, 10.0, 10.0, 5.0]

P = 6000; L = 14; E = 30e+6; G = 12e+6;
τₘₐₓ = 13600; σₘₐₓ = 30000; δₘₐₓ = 0.25;
M(x) = P*(L+x[2]/2)
R(x) = sqrt(0.25*(x[2]^2+(x[1]+x[3])^2))
J(x) = 2/sqrt(2)*x[1]*x[2]*(x[2]^2/12+0.25*(x[1]+x[3])^2)
P_c(x) = (4.013*E/(6*L^2))*x[3]*x[4]^3*(1-0.25*x[3]*sqrt(E/G)/L)
τ₁(x) = P/(sqrt(2)*x[1]*x[2])
τ₂(x) = M(x)*R(x)/J(x)
τ(x) = sqrt(τ₁(x)^2+τ₁(x)*τ₂(x)*x[2]/R(x)+τ₂(x)^2)
σ(x) = 6*P*L/(x[4]*x[3]^2)
δ(x) = 4*P*L^3/(E*x[4]*x[3]^3)

con_c1(x) = [τ(x), σ(x), x[4]-x[1], 0.10471*x[1]^2+0.04811*x[3]*x[4]*(14.0+x[2]), P_c(x), δ(x)]
lc = [-Inf, -Inf, 0.0, -Inf,   P, -Inf]
uc = [τₘₐₓ, σₘₐₓ, Inf,  5.0, Inf, δₘₐₓ]

con_c2(x) = [τ(x)-τₘₐₓ, σ(x)-σₘₐₓ, x[1]-x[4], 0.10471*x[1]^2+0.04811*x[3]*x[4]*(14.0+x[2])-5.0, P-P_c(x), δ(x)-δₘₐₓ]
lc = fill(-Inf, 6)
uc = fill(0.0, 6)

opts = Evolutionary.Options(iterations=2000, successive_f_tol=40)
mthd = CMAES(mu = 90, c_mu=0.3)
mthd = GA(populationSize = 120, ɛ = 0.03, crossoverRate=0.8, mutationRate=0.01, selection=rouletteinv, crossover=LX(0.0,4.0), mutation = PM(lx,ux,1.0))
mthd = DE(n=2)

c = PenaltyConstraints([1e3, 1e3, 1e3, 1e3, 1e2, 1e1, 1e3, 1e1, 1e1, 1e0], lx, ux, lc, uc, con_c1)
c = PenaltyConstraints([1e3, 1e3, 1e3, 1e3, 1e1, 1e1, 1e3, 1e1, 1e1, 1e0], Evolutionary.ConstraintBounds(lx,ux,lc,[]), con_c2)
c = WorstFitnessConstraints(Evolutionary.ConstraintBounds(lx,ux,lc,uc), con_c1)
c = WorstFitnessConstraints(Evolutionary.ConstraintBounds(lx,ux,lc,uc), con_c2)

Random.seed!(0);
result = Evolutionary.optimize(beam, c, mthd, opts)
Evolutionary.minimizer(result) |> con_c1
Evolutionary.minimizer(result) |> con_c2
[Evolutionary.optimize(beam, c, mthd, opts) |> minimum for i in 1:100] |> minimum


text:
# type
mutable struct DQMCStack{GreensEltype<:Number, HoppingEltype<:Number} <: AbstractDQMCStack
  eye_flv::Matrix{Float64}
  eye_full::Matrix{Float64}
  ones_vec::Vector{Float64}

  u_stack::Array{GreensEltype, 3}
  d_stack::Matrix{Float64}
  t_stack::Array{GreensEltype, 3}

  Ul::Matrix{GreensEltype}
  Ur::Matrix{GreensEltype}
  Dl::Vector{Float64}
  Dr::Vector{Float64}
  Tl::Matrix{GreensEltype}
  Tr::Matrix{GreensEltype}

  greens::Matrix{GreensEltype}
  greens_temp::Matrix{GreensEltype}
  log_det::Float64 # contains logdet of greens_{mc.p.slices+1} === greens_1
                            # after we calculated a fresh greens in propagate()

  U::Matrix{GreensEltype}
  D::Vector{Float64}
  T::Matrix{GreensEltype}
  tmp::Matrix{GreensEltype}

  ranges::Array{UnitRange, 1}
  n_elements::Int
  current_slice::Int # running internally over 0:mc.p.slices+1, where 0 and mc.p.slices+1 are artifcial to prepare next sweep direction.
  direction::Int

  # # -------- Global update backup
  # gb_u_stack::Array{GreensEltype, 3}
  # gb_d_stack::Matrix{Float64}
  # gb_t_stack::Array{GreensEltype, 3}

  # gb_greens::Matrix{GreensEltype}
  # gb_log_det::Float64

  # gb_conf::Array{Float64, 3}
  # # --------

  # preallocated, reused arrays
  curr_U::Matrix{GreensEltype}
  eV::Matrix{GreensEltype}

  # hopping matrices
  hopping_matrix_exp::Matrix{HoppingEltype} # mu included
  hopping_matrix_exp_inv::Matrix{HoppingEltype} # mu included

  # checkerboard hopping matrices
  checkerboard::Matrix{Int} # src, trg, bondid
  groups::Vector{UnitRange}
  n_groups::Int
  chkr_hop_half::Vector{SparseMatrixCSC{HoppingEltype, Int64}}
  chkr_hop_half_inv::Vector{SparseMatrixCSC{HoppingEltype, Int64}}
  chkr_hop_half_dagger::Vector{SparseMatrixCSC{HoppingEltype, Int64}}
  chkr_hop::Vector{SparseMatrixCSC{HoppingEltype, Int64}} # without prefactor 0.5 in matrix exponentials
  chkr_hop_inv::Vector{SparseMatrixCSC{HoppingEltype, Int64}}
  chkr_hop_dagger::Vector{SparseMatrixCSC{HoppingEltype, Int64}}
  chkr_mu_half::SparseMatrixCSC{HoppingEltype, Int64}
  chkr_mu_half_inv::SparseMatrixCSC{HoppingEltype, Int64}
  chkr_mu::SparseMatrixCSC{HoppingEltype, Int64}
  chkr_mu_inv::SparseMatrixCSC{HoppingEltype, Int64}

  DQMCStack{GreensEltype, HoppingEltype}() where {GreensEltype<:Number, HoppingEltype<:Number} = begin
    # @assert isleaftype(GreensEltype);
    # @assert isleaftype(HoppingEltype);
    @assert isconcretetype(GreensEltype);
    @assert isconcretetype(HoppingEltype);
    new()
  end
end

# type helpers
geltype(::Type{DQMCStack{G,H}}) where {G,H} = G
heltype(::Type{DQMCStack{G,H}}) where {G,H} = H
geltype(mc::DQMC{M, CB, CT, S}) where {M, CB, CT, S} = geltype(S)
heltype(mc::DQMC{M, CB, CT, S}) where {M, CB, CT, S} = heltype(S)

# type initialization
function initialize_stack(mc::DQMC)
  GreensEltype = geltype(mc)
  HoppingEltype = heltype(mc)
  N = mc.model.l.sites
  flv = mc.model.flv

  mc.s.eye_flv = Matrix{Float64}(I, flv,flv)
  mc.s.eye_full = Matrix{Float64}(I, flv*N,flv*N)
  mc.s.ones_vec = ones(flv*N)

  mc.s.n_elements = convert(Int, mc.p.slices / mc.p.safe_mult) + 1

  mc.s.u_stack = zeros(GreensEltype, flv*N, flv*N, mc.s.n_elements)
  mc.s.d_stack = zeros(Float64, flv*N, mc.s.n_elements)
  mc.s.t_stack = zeros(GreensEltype, flv*N, flv*N, mc.s.n_elements)

  mc.s.greens = zeros(GreensEltype, flv*N, flv*N)
  mc.s.greens_temp = zeros(GreensEltype, flv*N, flv*N)

  mc.s.Ul = Matrix{GreensEltype}(I, flv*N, flv*N)
  mc.s.Ur = Matrix{GreensEltype}(I, flv*N, flv*N)
  mc.s.Tl = Matrix{GreensEltype}(I, flv*N, flv*N)
  mc.s.Tr = Matrix{GreensEltype}(I, flv*N, flv*N)
  mc.s.Dl = ones(Float64, flv*N)
  mc.s.Dr = ones(Float64, flv*N)

  mc.s.U = zeros(GreensEltype, flv*N, flv*N)
  mc.s.D = zeros(Float64, flv*N)
  mc.s.T = zeros(GreensEltype, flv*N, flv*N)
  mc.s.tmp = zeros(GreensEltype, flv*N, flv*N)

  # # Global update backup
  # mc.s.gb_u_stack = zero(mc.s.u_stack)
  # mc.s.gb_d_stack = zero(mc.s.d_stack)
  # mc.s.gb_t_stack = zero(mc.s.t_stack)
  # mc.s.gb_greens = zero(mc.s.greens)
  # mc.s.gb_log_det = 0.
  # mc.s.gb_conf = zero(mc.conf)

  mc.s.ranges = UnitRange[]

  for i in 1:mc.s.n_elements - 1
    push!(mc.s.ranges, 1 + (i - 1) * mc.p.safe_mult:i * mc.p.safe_mult)
  end

  mc.s.curr_U = zero(mc.s.U)
  mc.s.eV = zeros(GreensEltype, flv*N, flv*N)

  # mc.s.hopping_matrix_exp = zeros(HoppingEltype, flv*N, flv*N)
  # mc.s.hopping_matrix_exp_inv = zeros(HoppingEltype, flv*N, flv*N)
  nothing
end

# hopping
function init_hopping_matrices(mc::DQMC{M,CB}, m::Model) where {M, CB<:Checkerboard}
  init_hopping_matrix_exp(mc, m)
  CB <: CheckerboardTrue && init_checkerboard_matrices(mc, m)
  nothing
end
function init_hopping_matrix_exp(mc::DQMC, m::Model)
  N = m.l.sites
  flv = m.flv
  dtau = mc.p.delta_tau

  T = hopping_matrix(mc, m)
  size(T) == (flv*N, flv*N) || error("Hopping matrix should have size "*
                                "$((flv*N, flv*N)) but has size $(size(T)) .")
  mc.s.hopping_matrix_exp = exp(-0.5 * dtau * T)
  mc.s.hopping_matrix_exp_inv = exp(0.5 * dtau * T)
  nothing
end

# checkerboard
rem_eff_zeros!(X::AbstractArray) = map!(e -> abs.(e)<1e-15 ? zero(e) : e,X,X)
function init_checkerboard_matrices(mc::DQMC, m::Model)
  s = mc.s
  l = m.l
  flv = m.flv
  H = heltype(mc)
  N = m.l.sites
  dtau = mc.p.delta_tau
  mu = m.mu

  s.checkerboard, s.groups, s.n_groups = build_checkerboard(l)
  n_grps = s.n_groups
  cb = s.checkerboard

  T = reshape(hopping_matrix(mc, m), (N, flv, N, flv))

  s.chkr_hop_half = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)
  s.chkr_hop_half_inv = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)
  s.chkr_hop = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)
  s.chkr_hop_inv = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)

  for (g, gr) in enumerate(s.groups)
    Tg = zeros(H, N, flv, N, flv)
    for i in gr
      src, trg = cb[1:2,i]
      for f1 in 1:flv, f2 in 1:flv
        Tg[trg, f1, src, f2] = T[trg, f1, src, f2]
      end
    end

    Tgg = reshape(Tg, (N*flv, N*flv))
    s.chkr_hop_half[g] = sparse(rem_eff_zeros!(exp(-0.5 * dtau * Tgg)))
    s.chkr_hop_half_inv[g] = sparse(rem_eff_zeros!(exp(0.5 * dtau * Tgg)))
    s.chkr_hop[g] = sparse(rem_eff_zeros!(exp(- dtau * Tgg)))
    s.chkr_hop_inv[g] = sparse(rem_eff_zeros!(exp(dtau * Tgg)))
  end

  s.chkr_hop_half_dagger = adjoint.(s.chkr_hop_half)
  s.chkr_hop_dagger = adjoint.(s.chkr_hop)

  mus = diag(reshape(T, (N*flv, N*flv)))
  s.chkr_mu_half = spdiagm(0 => exp.(-0.5 * dtau * mus))
  s.chkr_mu_half_inv = spdiagm(0 => exp.(0.5 * dtau * mus))
  s.chkr_mu = spdiagm(0 => exp.(-dtau * mus))
  s.chkr_mu_inv = spdiagm(0 => exp.(dtau * mus))

  # hop_mat_exp_chkr = foldl(*,s.chkr_hop_half) * sqrt.(s.chkr_mu)
  # r = effreldiff(s.hopping_matrix_exp,hop_mat_exp_chkr)
  # r[find(x->x==zero(x),hop_mat_exp_chkr)] = 0.
  # println("Checkerboard - Exact ≈ ", round(maximum(absdiff(s.hopping_matrix_exp,hop_mat_exp_chkr)), 4))
  nothing
end

# stack construction
"""
Build stack from scratch.
"""
function build_stack(mc::DQMC)
  mc.s.u_stack[:, :, 1] = mc.s.eye_full
  mc.s.d_stack[:, 1] = mc.s.ones_vec
  mc.s.t_stack[:, :, 1] = mc.s.eye_full

  @inbounds for i in 1:length(mc.s.ranges)
    add_slice_sequence_left(mc, i)
  end

  mc.s.current_slice = mc.p.slices + 1
  mc.s.direction = -1

  nothing
end
"""
Updates stack[idx+1] based on stack[idx]
"""
function add_slice_sequence_left(mc::DQMC, idx::Int)
  copyto!(mc.s.curr_U, mc.s.u_stack[:, :, idx])

  # println("Adding slice seq left $idx = ", mc.s.ranges[idx])
  for slice in mc.s.ranges[idx]
    multiply_slice_matrix_left!(mc, mc.model, slice, mc.s.curr_U)
  end

  mc.s.curr_U *= spdiagm(0 => mc.s.d_stack[:, idx])
  mc.s.u_stack[:, :, idx + 1], mc.s.d_stack[:, idx + 1], T = udt(mc.s.curr_U)
  mc.s.t_stack[:, :, idx + 1] =  T * mc.s.t_stack[:, :, idx]
end
"""
Updates stack[idx] based on stack[idx+1]
"""
function add_slice_sequence_right(mc::DQMC, idx::Int)
  copyto!(mc.s.curr_U, mc.s.u_stack[:, :, idx + 1])

  for slice in reverse(mc.s.ranges[idx])
    multiply_daggered_slice_matrix_left!(mc, mc.model, slice, mc.s.curr_U)
  end

  mc.s.curr_U *=  spdiagm(0 => mc.s.d_stack[:, idx + 1])
  mc.s.u_stack[:, :, idx], mc.s.d_stack[:, idx], T = udt(mc.s.curr_U)
  mc.s.t_stack[:, :, idx] = T * mc.s.t_stack[:, :, idx + 1]
end

# Green's function calculation
"""
Calculates G(slice) using mc.s.Ur,mc.s.Dr,mc.s.Tr=B(slice)' ... B(M)' and
mc.s.Ul,mc.s.Dl,mc.s.Tl=B(slice-1) ... B(1)
"""
function calculate_greens(mc::DQMC)
    # U, D, T = udt(mc.s.greens) after this
    mc.s.U, mc.s.D, mc.s.T = udt_inv_one_plus(
        UDT(mc.s.Ul, mc.s.Dl, mc.s.Tl),
        UDT(mc.s.Ur, mc.s.Dr, mc.s.Tr),
        tmp = mc.s.U, tmp2 = mc.s.T, tmp3 = mc.s.tmp,
        internaluse = true
    )
    mul!(mc.s.tmp, mc.s.U, Diagonal(mc.s.D))
    mul!(mc.s.greens, mc.s.tmp, mc.s.T)
    mc.s.greens
end

"""
Only reasonable immediately after calculate_greens()!
"""
function calculate_logdet(mc::DQMC)
    mc.s.log_det = real(
        log(complex(det(mc.s.U))) +
        sum(log.(mc.s.D)) +
        log(complex(det(mc.s.T)))
    )
  # mc.s.log_det = real(logdet(mc.s.U) + sum(log.(mc.s.D)) + logdet(mc.s.T))
end

# Green's function propagation
@inline function wrap_greens!(mc::DQMC, gf::Matrix, curr_slice::Int, direction::Int)
  if direction == -1
    multiply_slice_matrix_inv_left!(mc, mc.model, curr_slice - 1, gf)
    multiply_slice_matrix_right!(mc, mc.model, curr_slice - 1, gf)
  else
    multiply_slice_matrix_left!(mc, mc.model, curr_slice, gf)
    multiply_slice_matrix_inv_right!(mc, mc.model, curr_slice, gf)
  end
  nothing
end
# @inline function wrap_greens(mc::DQMC, gf::Matrix,slice::Int,direction::Int)
#   temp = copy(gf)
#   wrap_greens!(mc, temp, slice, direction)
#   return temp
# end
function propagate(mc::DQMC)
  if mc.s.direction == 1
    if mod(mc.s.current_slice, mc.p.safe_mult) == 0
      mc.s.current_slice +=1 # slice we are going to
      if mc.s.current_slice == 1
        mc.s.Ur[:, :], mc.s.Dr[:], mc.s.Tr[:, :] = mc.s.u_stack[:, :, 1], mc.s.d_stack[:, 1], mc.s.t_stack[:, :, 1]
        mc.s.u_stack[:, :, 1] = mc.s.eye_full
        mc.s.d_stack[:, 1] = mc.s.ones_vec
        mc.s.t_stack[:, :, 1] = mc.s.eye_full
        mc.s.Ul[:,:], mc.s.Dl[:], mc.s.Tl[:,:] = mc.s.u_stack[:, :, 1], mc.s.d_stack[:, 1], mc.s.t_stack[:, :, 1]

        calculate_greens(mc) # greens_1 ( === greens_{m+1} )
        calculate_logdet(mc)

      elseif 1 < mc.s.current_slice <= mc.p.slices
        idx = Int((mc.s.current_slice - 1)/mc.p.safe_mult)

        mc.s.Ur[:, :], mc.s.Dr[:], mc.s.Tr[:, :] = mc.s.u_stack[:, :, idx+1], mc.s.d_stack[:, idx+1], mc.s.t_stack[:, :, idx+1]
        add_slice_sequence_left(mc, idx)
        mc.s.Ul[:,:], mc.s.Dl[:], mc.s.Tl[:,:] = mc.s.u_stack[:, :, idx+1], mc.s.d_stack[:, idx+1], mc.s.t_stack[:, :, idx+1]

        if mc.p.all_checks
          mc.s.greens_temp = copy(mc.s.greens)
        end

        wrap_greens!(mc, mc.s.greens_temp, mc.s.current_slice - 1, 1)

        calculate_greens(mc) # greens_{slice we are propagating to}

        if mc.p.all_checks
          greensdiff = maximum(abs.(mc.s.greens_temp - mc.s.greens)) # OPT: could probably be optimized through explicit loop
          if greensdiff > 1e-7
            @printf("->%d \t+1 Propagation instability\t %.1e\n", mc.s.current_slice, greensdiff)
          end
        end

      else # we are going to mc.p.slices+1
        idx = mc.s.n_elements - 1
        add_slice_sequence_left(mc, idx)
        mc.s.direction = -1
        mc.s.current_slice = mc.p.slices+1 # redundant
        propagate(mc)
      end

    else
      # Wrapping
      wrap_greens!(mc, mc.s.greens, mc.s.current_slice, 1)
      mc.s.current_slice += 1
    end

  else # mc.s.direction == -1
    if mod(mc.s.current_slice-1, mc.p.safe_mult) == 0
      mc.s.current_slice -= 1 # slice we are going to
      if mc.s.current_slice == mc.p.slices
        mc.s.Ul[:, :], mc.s.Dl[:], mc.s.Tl[:, :] = mc.s.u_stack[:, :, end], mc.s.d_stack[:, end], mc.s.t_stack[:, :, end]
        mc.s.u_stack[:, :, end] = mc.s.eye_full
        mc.s.d_stack[:, end] = mc.s.ones_vec
        mc.s.t_stack[:, :, end] = mc.s.eye_full
        mc.s.Ur[:,:], mc.s.Dr[:], mc.s.Tr[:,:] = mc.s.u_stack[:, :, end], mc.s.d_stack[:, end], mc.s.t_stack[:, :, end]

        calculate_greens(mc) # greens_{mc.p.slices+1} === greens_1
        calculate_logdet(mc) # calculate logdet for potential global update

        # wrap to greens_{mc.p.slices}
        wrap_greens!(mc, mc.s.greens, mc.s.current_slice + 1, -1)

      elseif 0 < mc.s.current_slice < mc.p.slices
        idx = Int(mc.s.current_slice / mc.p.safe_mult) + 1
        mc.s.Ul[:, :], mc.s.Dl[:], mc.s.Tl[:, :] = mc.s.u_stack[:, :, idx], mc.s.d_stack[:, idx], mc.s.t_stack[:, :, idx]
        add_slice_sequence_right(mc, idx)
        mc.s.Ur[:,:], mc.s.Dr[:], mc.s.Tr[:,:] = mc.s.u_stack[:, :, idx], mc.s.d_stack[:, idx], mc.s.t_stack[:, :, idx]

        if mc.p.all_checks
          mc.s.greens_temp = copy(mc.s.greens)
        end

        calculate_greens(mc)

        if mc.p.all_checks
          greensdiff = maximum(abs.(mc.s.greens_temp - mc.s.greens)) # OPT: could probably be optimized through explicit loop
          if greensdiff > 1e-7
            @printf("->%d \t-1 Propagation instability\t %.1e\n", mc.s.current_slice, greensdiff)
          end
        end

        wrap_greens!(mc, mc.s.greens, mc.s.current_slice + 1, -1)

      else # we are going to 0
        idx = 1
        add_slice_sequence_right(mc, idx)
        mc.s.direction = 1
        mc.s.current_slice = 0 # redundant
        propagate(mc)
      end

    else
      # Wrapping
      wrap_greens!(mc, mc.s.greens, mc.s.current_slice, -1)
      mc.s.current_slice -= 1
    end
  end
  nothing
end

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
fig = Figure()
ax = Axis(fig[1, 1])

fid = h5open(filename, "r")

stringslider = Slider(fig[2, 1], range=1:120, startvalue=0)
slidervalue = stringslider.value
field = @lift(read(fid["moisture"][string($slidervalue)]))
colorrange = @lift(extrema($field))
heatmap!(ax, field, colormap = :balance, colorrange = colorrange)

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
module CombinatorialSpaces

using Reexport
using Requires

include("ArrayUtils.jl")
include("CombinatorialMaps.jl")
include("ExteriorCalculus.jl")
include("SimplicialSets.jl")
include("DiscreteExteriorCalculus.jl")
include("MeshInterop.jl")

@reexport using .SimplicialSets
@reexport using .DiscreteExteriorCalculus
@reexport using .MeshInterop

function __init__()
  @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("MeshGraphics.jl")
end

end


text:
"""
Determines if of Image Type
"""
function isa_image(x::Type{T}) where T<:Matrix
    eltype(T) <: Union{Colorant, Colors.Fractional}
end
isa_image(x::Matrix) = isa_image(typeof(x))

isa_image(x) = false

# Splits a dictionary in two dicts, via a condition
function Base.split(condition::Function, associative::AbstractDict)
    A = similar(associative)
    B = similar(associative)
    for (key, value) in associative
        if condition(key, value)
            A[key] = value
        else
            B[key] = value
        end
    end
    A, B
end

function assemble_robj(data, program, bb, primitive, pre_fun, post_fun)
    transp = get(data, :transparency, Node(false))
    overdraw = get(data, :overdraw, Node(false))
    pre = if pre_fun != nothing
        _pre_fun = GLAbstraction.StandardPrerender(transp, overdraw)
        function ()
            _pre_fun()
            pre_fun()
        end
    else
        GLAbstraction.StandardPrerender(transp, overdraw)
    end
    robj = RenderObject(data, program, pre, nothing, bb, nothing)
    post = if haskey(data, :instances)
        GLAbstraction.StandardPostrenderInstanced(data[:instances], robj.vertexarray, primitive)
    else
        GLAbstraction.StandardPostrender(robj.vertexarray, primitive)
    end
    robj.postrenderfunction = if post_fun != nothing
        () -> begin
            post()
            post_fun()
        end
    else
        post
    end
    robj
end

function assemble_shader(data)
    shader = data[:shader]
    delete!(data, :shader)
    default_bb = Node(GeometryTypes.centered(AABB))
    bb  = get(data, :boundingbox, default_bb)
    if bb == nothing || isa(bb, Node{Nothing})
        bb = default_bb
    end
    glp = get(data, :gl_primitive, GL_TRIANGLES)
    robj = assemble_robj(
        data, shader, bb, glp,
        get(data, :prerender, nothing),
        get(data, :postrender, nothing)
    )
    Context(robj)
end

function y_partition_abs(area, amount)
    a = round(Int, amount)
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, r.w, a),
            SimpleRectangle{Int}(0, a, r.w, r.h - a)
        )
    end
    return lift(first, p), lift(last, p)
end
function x_partition_abs(area, amount)
    a = round(Int, amount)
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, a, r.h),
            SimpleRectangle{Int}(a, 0, r.w - a, r.h)
        )
    end
    return lift(first, p), lift(last, p)
end

function y_partition(area, percent)
    amount = percent / 100.0
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, r.w, round(Int, r.h*amount)),
            SimpleRectangle{Int}(0, round(Int, r.h*amount), r.w, round(Int, r.h*(1-amount)))
        )
    end
    return lift(first, p), lift(last, p)
end
function x_partition(area, percent)
    amount = percent / 100.0
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, round(Int, r.w*amount), r.h ),
            SimpleRectangle{Int}(round(Int, r.w*amount), 0, round(Int, r.w*(1-amount)), r.h)
        )
    end
    return lift(first, p), lift(last, p)
end

glboundingbox(mini, maxi) = AABB{Float32}(Vec3f0(mini), Vec3f0(maxi)-Vec3f0(mini))
function default_boundingbox(main, model)
    main == nothing && return Node(AABB{Float32}(Vec3f0(0), Vec3f0(1)))
    const_lift(*, model, AABB{Float32}(main))
end
AABB(a::GPUArray) = AABB{Float32}(gpu_data(a))
AABB{T}(a::GPUArray) where {T} = AABB{T}(gpu_data(a))

points2f0(positions::Vector{T}, range::AbstractRange) where {T} = Point2f0[Point2f0(range[i], positions[i]) for i=1:length(range)]

extrema2f0(x::Array{T,N}) where {T<:Intensity,N} = Vec2f0(extrema(reinterpret(Float32,x)))
extrema2f0(x::Array{T,N}) where {T,N} = Vec2f0(extrema(x))
extrema2f0(x::GPUArray) = extrema2f0(gpu_data(x))
function extrema2f0(x::Array{T,N}) where {T<:Vec,N}
    _norm = map(norm, x)
    Vec2f0(minimum(_norm), maximum(_norm))
end

function mix_linearly(a::C, b::C, s) where C<:Colorant
    RGBA{Float32}((1-s)*comp1(a)+s*comp1(b), (1-s)*comp2(a)+s*comp2(b), (1-s)*comp3(a)+s*comp3(b), (1-s)*alpha(a)+s*alpha(b))
end

color_lookup(cmap, value, mi, ma) = color_lookup(cmap, value, (mi, ma))
function color_lookup(cmap, value, color_norm)
    mi,ma = color_norm
    scaled = clamp((value-mi)/(ma-mi), 0, 1)
    index = scaled * (length(cmap)-1)
    i_a, i_b = floor(Int, index)+1, ceil(Int, index)+1
    mix_linearly(cmap[i_a], cmap[i_b], scaled)
end

"""
Converts index arrays to the OpenGL equivalent.
"""
to_index_buffer(x::GLBuffer) = x
to_index_buffer(x::TOrSignal{Int}) = x
to_index_buffer(x::VecOrSignal{UnitRange{Int}}) = x
to_index_buffer(x::TOrSignal{UnitRange{Int}}) = x
"""
For integers, we transform it to 0 based indices
"""
to_index_buffer(x::AbstractVector{I}) where {I <: Integer} = indexbuffer(Cuint.(x .- 1))
function to_index_buffer(x::Node{<: AbstractVector{I}}) where I <: Integer
    indexbuffer(lift(x-> Cuint.(x .- 1), x))
end

"""
If already GLuint, we assume its 0 based (bad heuristic, should better be solved with some Index type)
"""
function to_index_buffer(x::VectorTypes{I}) where I <: Union{GLuint, Face{2, GLIndex}}
    indexbuffer(x)
end

to_index_buffer(x) = error(
    "Not a valid index type: $(typeof(x)).
    Please choose from Int, Vector{UnitRange{Int}}, Vector{Int} or a signal of either of them"
)

"""
Creates a moving average and discards values to close together.
If discarded return (false, p), if smoothed, (true, smoothed_p).
"""
function moving_average(p, cutoff,  history, n = 5)
    if length(history) > 0
        if norm(p - history[end]) < cutoff
            return false, p # don't keep point
        end
    end
    if length(history) == 5
        # maybe better to just keep a current index
        history[1:5] = circshift(view(history, 1:5), -1)
        history[end] = p
    else
        push!(history, p)
    end
    true, sum(history) ./ length(history)# smooth
end

function layoutlinspace(n::Integer)
    if n == 1
        1:1
    else
        range(1/n, stop=1, length=n)
    end
end
xlayout(x::Int) = zip(layoutlinspace(x), Iterators.repeated(""))
function xlayout(x::AbstractVector{T}) where T <: AbstractFloat
    zip(x, Iterators.repeated(""))
end

function xlayout(x::AbstractVector)
    zip(layoutlinspace(length(x)), x)
end
function ylayout(x::AbstractVector)
    zip(layoutlinspace(length(x)), x)
end
function ylayout(x::AbstractVector{T}) where T <: Tuple
    sizes = map(first, x)
    values = map(last, x)
    zip(sizes, values)
end

function layout_rect(area, lastw, lasth, w, h)
    wp = widths(area)
    xmin = wp[1] * lastw
    ymin = wp[2] * lasth
    xmax = wp[1] * w
    ymax = wp[2] * h
    xmax = max(xmin, xmax)
    xmin = min(xmin, xmax)
    ymax = max(ymin, ymax)
    ymin = min(ymin, ymax)
    AbstractPlotting.IRect(xmin, ymin, xmax - xmin, ymax - ymin)
end


text:
#! format: off

using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder

using Cbc #solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

sys = build_system(PSITestSystems, "test_RTS_GMLC_sys")

for line in get_components(Line, sys)
    if (get_base_voltage(get_from(get_arc(line))) >= 230.0) &&
       (get_base_voltage(get_to(get_arc(line))) >= 230.0)
        #if get_area(get_from(get_arc(line))) != get_area(get_to(get_arc(line)))
        @info "Changing $(get_name(line)) to MonitoredLine"
        convert_component!(MonitoredLine, line, sys)
    end
end

template = template_unit_commitment(transmission = PTDFPowerModel)

set_device_model!(template, MonitoredLine, StaticBranch)

set_device_model!(template, Line, StaticBranchUnbounded)

uc_prob = OperationsProblem(template, sys, horizon = 24, optimizer = solver)
build!(uc_prob, output_dir = mktempdir())

solve!(uc_prob)

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


text:
# Autogenerated wrapper script for snappy_jll for aarch64-linux-gnu-cxx11
export libsnappy

using Zlib_jll
using LZO_jll
JLLWrappers.@generate_wrapper_header("snappy")
JLLWrappers.@declare_library_product(libsnappy, "libsnappy.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, LZO_jll)
    JLLWrappers.@init_library_product(
        libsnappy,
        "lib/libsnappy.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Distributions, Statistics, LinearAlgebra

import Base: size, length, show

import Statistics: mean, var, std

export ETKF

""""
Define system ETKF which performs the data assimilation
using the ensemble transform Kalman filter (Bishop, et al. 2001)

Fields:
 - 'f' : propagation function

 - 'A' : inflation

 - 'G' : filtering function acting on the state

 - 'm' : measurement function based on state

 - 'z' : real measurement function

 - 'ϵ' : measurement noise distribution

 - 'bounds' : bounds on certain states

 - 'isinflated' : Bool = true if state is inflated,
     = false otherwise

 - 'isfiltered' : Bool = true if state has to be filtered,
     = false otherwise

 - 'isaugmented' : Bool = true if measurement function is nonlinear,
     = false otherwise
"""
mutable struct ETKF{N, NZ}

    # "Ensemble of states"
    # ENS::EnsembleState{N, NS, TS}

    "Propagation function"
    f::PropagationFunction

    "Covariance Inflation"
    A::Union{InflationType, RecipeInflation}

    "Filter function"
    G::FilteringFunction

    "Measurement function based on state"
    m::MeasurementFunction

    "Real measurement function"
    z::RealMeasurementFunction

    "Measurement noise distribution"
    ϵ::AdditiveInflation{NZ}

    "Boolean: is state vector inflated"
    isinflated::Bool

    "Boolean: is state vector filtered"
    isfiltered::Bool

    "Boolean: is state vector augmented"
    isaugmented::Bool

    # "Bounds on certain state"
    # bounds

end

"""
    Define action of ETKF on EnsembleState
"""
function (etkf::ETKF{N, NZ})(t::Float64,
         Δt::Float64,
         ens::EnsembleState{N, TS}) where {N, NZ, TS}

    "Propagate each ensemble member"
    etkf.f(t, ens)

    "Covariance inflation if 'isinflated==true' "
    if etkf.isinflated ==true
        etkf.A(ens)
    end

    x̄ᶠ = mean(deepcopy(ens))

    Hx̄ᶠ = etkf.m(t, x̄ᶠ)

    d = etkf.z(t) - Hx̄ᶠ
    "Compute measurement"
    mens = EnsembleState(N, zeros(NZ))

        for (i, s) in enumerate(ens.S)
            mens.S[i] = etkf.m(t, deepcopy(s))
        end
    HXᶠ = hcat(deepcopy(mens))

    "Compute deviation from measurement of the mean"
    Ŝ = mean(deepcopy(ens))

    HX′ᶠ  = HXᶠ.- Hx̄ᶠ

    R = cov(etkf.ϵ)
    C = inv(R)*HX′ᶠ

    A₂ = (N-1)*I + (HX′ᶠ)'*C

    EVDA₂ = Eigen(A₂)
    U = EVDA₂.vectors
    DΣ  = EVDA₂.values
    X′ᶠ = hcat(ens) .- x̄ᶠ
    D = C'*d
    w₁ = U'*D
    w₂  = w₁ ./ DΣ
    w̄ = U*w₂
    W₁′ = zeros(N, N)
    for j=1:N
        W₁′[:, j] .= √(DΣ[j])* U[:,j]
    end
    W′ = W₁′*U'
    W = W′ .+ w̄
    ens += cut(X′ᶠ*W)

    " Compute a posteriori covariance"
    ensfluc = EnsembleState(N, ens.S[1])
    deviation(ensfluc, ens)

    A′ = hcat(ensfluc)

    return t+Δt, ens, A′*A′'

end

# Create constructor for ETKF

function ETKF(N, NZ, f, A, G, m, z, ϵ, isinflated, isfiltered, isaugmented)
    return ETKF{N, NZ}(f, A, G, m, z, ϵ, isinflated, isfiltered, isaugmented)
end

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF


text:
using NNlib, ImageFiltering, Images
using DrWatson
using Statistics
using BenchmarkTools
@quickactivate "MEngProject"
imgLoc = datadir("img","stairs_100gs.png")

img = convert(Array{Float32,2}, load(imgLoc))

function reshape2d_4d(img::AbstractArray)
    reshape(img, size(img)[1], size(img)[2], 1, 1)
end

img = reshape2d_4d(img)

@benchmark mean(img[:,:,1,1])

Statistics.mean(img[:,:,1,1])

a = similar(img)
a = [1.,1.,1.,1.]

@benchmark mean!(a, img[:,:,1,1])

a

a = reshape(
    Array{eltype(img)}(undef, 100, 100 * 2),
    100,
    100,
    2,
    1,
)

a[:,:,1:1,:] .= img
a[:,:,2:2,:] .= img

b = Array{eltype(img)}(undef,  2)

for k in 1:2
b[k] = mean(a[:,:,k,:])
end

b


text:
const commands = Dict{String, String}(
	"0000" => "Module",
	"0001" => "End",
	"0010" => "Sub",
	"0011" => "Main",
 	"0100" => "Console.WriteLine",
 	"0101" =>  "Console.ReadLine",
 	"0110" => "Dim",
 	"0111" => "Function",
 	"1000" => "Return",
 	"1001" => "As",
 	"1010" => "String",
 	"1011" => "args",
 	"1100" => "Using",
 	"1101" => "System",
 	"1110" => "For",
 	"1111" => "While"
)

@doc raw"""
Take the code to execute and replace the keys of commands in
the execute_code with the corresponding value and return
the value
"""
function execute_binary_c(execute_code::String)::String
	for (key, value) in commands
		execute_code = replace(execute_code, key => value)
	end
	return execute_code
end

@doc raw"""
Write the output_code into the dist folder
"""
function write_to_file(execute_code::String):Nothing
	filepath::String = joinpath(pwd(), "dist")
	if !isdir(filepath)
		mkdir(filepath)
	end

	open(joinpath(filepath, "output.vb"), "w") do output_writer
		write(output_writer, execute_code)
	end

	println("Succesfully compiled into $filepath")
end

max_stars_repo_name: Technetium-Inc/BinaryC

max_issues_repo_name: Technetium-Inc/BinaryC

max_forks_repo_name: Technetium-Inc/BinaryC


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

using JuliaFEM, LinearAlgebra, Test

# 2d heat problem (one element)

X = Dict(
    1 => [0.0,0.0],
    2 => [1.0,0.0],
    3 => [1.0,1.0],
    4 => [0.0,1.0])

# define volume element
element1 = Element(Quad4, (1, 2, 3, 4))

update!(element1, "geometry", X)
update!(element1, "thermal conductivity", 6.0)
update!(element1, "heat source", 12.0)

# define boundary element for flux
element2 = Element(Seg2, (1, 2))
update!(element2, "geometry", X)
# linear ramp from 0 -> 6 in time 0 -> 1
update!(element2, "heat flux", 0.0 => 0.0)
update!(element2, "heat flux", 1.0 => 6.0)

# define heat problem and add elements to problem
problem = Problem(PlaneHeat, "one element heat problem", 1)
add_elements!(problem, element1, element2)

# Set constant source f=12 with k=6. Accurate solution is
# T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
# when boundary flux not active (at t=0)
time = 0.0
assemble!(problem, time)
A = Matrix(problem.assembly.K)
b = Vector(problem.assembly.f)
A_expected = [
     4.0 -1.0 -2.0 -1.0
    -1.0  4.0 -1.0 -2.0
    -2.0 -1.0  4.0 -1.0
    -1.0 -2.0 -1.0  4.0]
free_dofs = [1, 2]
@test isapprox(A, A_expected)
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])

# Set constant flux g=6 on boundary. Accurate solution is
# u(x,y) = x which equals T=1 on boundary.
# at time t=1.0 all loads should be on.
empty!(problem.assembly)
time = 1.0
assemble!(problem, time)
A = Matrix(problem.assembly.K)
b = Vector(problem.assembly.f)
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0])


text:
module GeodesicTracer

using SciMLBase
using OrdinaryDiffEq
using DiffEqCallbacks

using StaticArrays
using DocStringExtensions
using Parameters

import GeodesicBase: AbstractMetricParams, geodesic_eq, constrain, on_chart, inner_radius

import ForwardDiff

include("callbacks.jl")
include("problem.jl")
include("tracer.jl")
include("constraints.jl")
include("utility.jl")
include("auto-diff.jl")

"""
    tracegeodesics(
        m::AbstractMetricParams{T},
        init_positions, init_velocities,
        time_domain::Tuple{T,T}
        ;
        μ = 0.0f0,
        callbacks=Nothing,
        solver=Tsit5(),
        solver_opts...
    )

Integrate a geodesic for metric parameterised by `m`, for some initial positions and velocities.
The positions and velocities may be

  - a single position and velocity in the form of a vector of numbers
  - a collection of positions and velocities, as either a vector of vectors, or as a matrix

The matrix specification reads each corresponding column as the initial position and velocity. When a collection of
positions and velocities is supplied, this method dispatched `EnsembleProblem`, offering `ensemble` as a `solver_opts`,
specifying the ensemble method to use.

`solver_opts` are the common solver options in DifferentialEquations.
"""
function tracegeodesics(
    m::AbstractMetricParams{T},
    init_positions,
    init_velocities,
    time_domain::Tuple{T,T};
    solver = Tsit5(),
    μ = 0.0,
    solver_opts...
) where {T}
    __tracegeodesics(
        m,
        init_positions,
        # ensure everything already normalised
        constrain_all(m, init_positions, init_velocities, T(μ)),
        time_domain,
        solver;
        callback = nothing,
        abstol = 1e-8,
        reltol = 1e-8,
        solver_opts...
    )
end

export tracegeodesics

end # module


text:
using Base.depwarn

Base.@deprecate by(d::AbstractDataFrame, cols, s::Vector{Symbol}) aggregate(d, cols, map(eval, s))
Base.@deprecate by(d::AbstractDataFrame, cols, s::Symbol) aggregate(d, cols, eval(s))
Base.@deprecate nullable!(colnames::Array{Symbol,1}, df::AbstractDataFrame) nullable!(df, colnames)
Base.@deprecate nullable!(colnums::Array{Int,1}, df::AbstractDataFrame) nullable!(df, colnums)
import Base: keys, values, insert!
Base.@deprecate keys(df::AbstractDataFrame) names(df)
Base.@deprecate values(df::AbstractDataFrame) DataFrames.columns(df)
Base.@deprecate insert!(df::DataFrame, df2::AbstractDataFrame) merge!(df, df2)

import DataArrays: array, DataArray
Base.@deprecate array(df::AbstractDataFrame) convert(Array, df)
Base.@deprecate array(r::DataFrameRow) convert(Array, r)
if VERSION < v"0.4.0-"
    Base.@deprecate DataArray(df::AbstractDataFrame) convert(DataArray, df)
end
Base.@deprecate DataArray(df::AbstractDataFrame, T::DataType) convert(DataArray{T}, df)