text:
module FourierFilterFlux
using Reexport
# @reexport using CUDA
using CUDA
using Zygote, Flux, Shearlab, LinearAlgebra
using AbstractFFTs, FFTW # TODO: check the license on FFTW and such
using ContinuousWavelets
using Flux
using Adapt
using RecipesBase
using Base:tail

const use_cuda = Ref(false)

import Adapt:adapt
export pad, poolSize, originalDomain, params!, formatJLD, getBatchSize
export Periodic, Pad, ConvBoundary, Sym, analytic, outType, nFrames
# layer types and constructors
export ConvFFT, waveletLayer, shearingLayer, averagingLayer
# inits
export positive_glorot_uniform, iden_perturbed_gaussian,
    uniform_perturbed_gaussian
# Analytic types
export TransformTypes, AnalyticWavelet, RealWaveletRealSignal,
    RealWaveletComplexSignal, NonAnalyticMatching
# utils
export effectiveSize, originalSize

include("boundaries.jl")

"""
    # preconstructed filters w
    ConvFFT(w::AbstractArray{T,N}, b, originalSize, σ = identity; plan=true,
    boundary = Periodic(), dType=Float32, trainable=true, OT=Float32, An=nothing)
        where {T,N}
    # randomly constructed filters
    ConvFFT(k::NTuple{N,Integer}, nOutputChannels = 5,
            σ=identity; nConvDims=2, init = Flux.glorot_normal,
            plan=true, bias = true,
            dType=Float32, OT=Float32, boundary=Periodic(),
            trainable=true, An=nothing) where N

Similar to `Conv` from Flux.jl, but does pointwise multiplication in the
Fourier domain, with boundary conditions `boundary`, and applies the
nonlinearity `σ`. It is worth noting that the bias is added in the Fourier
domain; if you don't want a bias, set `b=nothing`. For the first method, the
weights `w` should have dimension 1 greater than the size of the convolution,
while for the second, `k` gives the size of the input (including channels and
batch size), for which appropriate weights are generated in the fourier domain
according to the distribution `init`.
# Shared Arguments
- `σ=identity`: a function to apply.
- `dType::DataType=Float32`: the data type being input. By default, it assumes
                            both the filter and the signal are real, which
                            allows us to use an rfft and half the (complex)
                            coefficients we would otherwise need.
- `OT::DataType=Float32`: the output datatype, usually determined by the (space
                          domain) type of the filters. Currently assumes that
                          if this is complex, then the filters are analytic
                          (so only defined for positive frequencies).
- `plan::Bool=true`: use precomputed fft plan(s), as this is a significant cost.
                     Set this to `false` if you have variable batch/channel sizes.
- `trainable::Bool=true`: The entries are trainable as Flux objects, and so are
                          returned in a `params` call if this is `true`.
- `boundary::ConvBoundary`: determines how the edges are treated. See e.g. `Pad`.
- `An::Union{Nothing,NTuple{<:Integer}}=nothing`: only used if the output type is
                    complex, so the filters are complex. In that case, it is a
                    list of the filters which are actually real. This is
                    included to allow for averaging wavelets which are real even
                    for analytic wavelets. For both `waveletLayer` and
                    `shearingLayer`, this is the last filter, so if there are
                    18 total wavelets, `An=(18,)`.
# First constructor only

- `w::AbstractArray{T,N}`: the weights. should be `D+1`, with the last
                           dimension being the total number of filters
` `b::Union{Nothing,AbstractArray}`: if `nothing`, no bias is added. Otherwise it should be (input channels)×(output channels)
# Second constructor only
- `k::NTuple{N,Integer}`: the dimensions of the input. In the 1D case it should
                          be three entries e.g. `(132,3,100)` which is
                          (signal)×(channels)×(examples). In the 2D case it
                          should be four entries e.g. `(132,132,3,100)`, which
                          is (x)×(y)×(channels)×(examples).
- `nOutputChannels::Integer=5`: the number of filters to use.

- `bias::Bool=true`: determines whether or not to create a bias.
- `init::function=Flux.glorot_normal`: The way to initialize both the bias (if
                    defined) and the weights.  Any function that results in a
                    matrix is allowed, though I would suggest something from
                    Flux, or one of the ones defined in this package.
- `nConvDims::Integer`: the number of dimensions that we will be doing the
                    convolution over, as `k` is somewhat ambiguous.
"""
struct ConvFFT{D,OT,F,A,V,PD,P,T,An}
    σ::F
    weight::A
    bias::V
    bc::PD
    fftPlan::P
    analytic::An
end
import Base:ndims
ndims(c::ConvFFT{D}) where D = D
function ConvFFT(w::AbstractArray{T,N}, b, originalSize, σ=identity; plan=true, boundary=Periodic(),
                 dType=Float32, trainable=true, OT=Float32,
                 An=nothing) where {T,N}
    @assert length(originalSize) >= N - 1
    if dType <: Complex
        OT = dType
    end

    if length(originalSize) == N - 1
        exSz = (originalSize..., 1) # default number of channels is 1
    else
        exSz = originalSize
    end
    netSize, boundary = effectiveSize(exSz[1:N - 1], boundary)
    convDims = (1:(N - 1)...,)

    # Check that they applied the boundary condition, and if not do it ourselves
    if dType <: Complex && netSize[1] != size(w, 1)
        wtmp = ifftshift(ifft(w, convDims), convDims)
        wtmp = applyBC(wtmp, boundary, N - 1)
        w = fft(fftshift(wtmp, convDims), convDims)
        @warn("You didn't hand me a set of filters constructed with the boundary in mind. I'm going to adjust them to fit, this may not be what you intended")
    elseif dType <: Real && netSize[1] >> 1 + 1 != size(w, 1)
        wtmp = ifftshift(irfft(w, exSz[1], convDims), convDims)
        wtmp, _ = applyBC(wtmp, boundary, N - 1)
        w = rfft(fftshift(wtmp, (convDims)), convDims)
        @warn("You didn't hand me a set of filters constructed with the boundary in mind. I'm going to adjust them to fit, this may not be what you intended")
    end

    if plan
        fftPlan = makePlan(dType, OT, w, exSz, boundary)
    else
        fftPlan = nothing
    end

    if typeof(An) <: Nothing
        An = map(x -> NonAnalyticMatching(), (1:size(w)[end]...,))
    end

    return ConvFFT{N - 1,OT,typeof(σ),typeof(w),typeof(b),typeof(boundary),typeof(fftPlan),trainable,typeof(An)}(σ, w, b, boundary, fftPlan, An)
end

function makePlan(dType, OT, w, exSz, boundary)
    N = ndims(w)
    netSize, boundary = effectiveSize(exSz[1:N - 1], boundary)
    convDims = (1:(N - 1)...,)
    nullEx = Adapt.adapt(typeof(w), zeros(dType, netSize..., exSz[N:end]...))
    if dType <: Real && OT <: Real
        fftPlan = plan_rfft(real.(nullEx), convDims)
    elseif dType <: Real # output is complex, wavelets analytic
        null2 = Adapt.adapt(typeof(w), zeros(dType, netSize..., exSz[N:end]...)) .+
            0im
        fftPlan = (plan_rfft(real.(nullEx), convDims), plan_fft!(null2, convDims))
    else
        fftPlan = plan_fft!(nullEx, convDims)
    end
end

function ConvFFT(k::NTuple{N,Integer}, nOutputChannels=5,
                 σ=identity; nConvDims=2, init=Flux.glorot_normal,
                 plan=true, bias=true,
                 dType=Float32, OT=Float32, boundary=Periodic(),
                 trainable=true, An=nothing) where N

    effSize, boundary = effectiveSize(k[1:nConvDims], boundary)
    if dType <: Real && OT <: Real
        effSize = (effSize[1] >> 1 + 1, effSize[2:end]...)
    end
    w = init(effSize..., nOutputChannels) .+
        im .* init(effSize..., nOutputChannels)

    if bias == true
        b = init(k[(nConvDims + 1):end - 1]..., nOutputChannels)
    else
        b = nothing
    end

    ConvFFT(w, b, k, σ, plan=plan, boundary=boundary, dType=dType,
            trainable=trainable, OT=OT, An=An)
end

function Base.show(io::IO, l::ConvFFT)
    # stored as a brief description
    if typeof(l.fftPlan) <: Tuple
        sz = l.fftPlan[2]
    else
        sz = l.fftPlan.sz
    end
    es = originalSize(sz[1:ndims(l.weight) - 1], l.bc)
    print(io, "ConvFFT[input=($(es), " *
          "nfilters = $(size(l.weight)[end]), " *
          "σ=$(l.σ), " *
          "bc=$(l.bc)]")
end

function Base.show(io::IO, l::ConvFFT{D,OT,A,B,C,PD,P}) where {D,OT,A,B,C,PD,P <: Tuple}
    if typeof(l.fftPlan[1]) <: Tuple
        sz = l.fftPlan[1][2]
    else
        sz = l.fftPlan[1].sz
    end
    es = originalSize(sz[1:ndims(l.weight) - 1], l.bc)
    print(io, "ConvFFT[input=($(es), " *
          "nfilters = $(size(l.weight)[end]), " *
          "σ=$(l.σ), " *
          "bc=$(l.bc)]")
end

analytic(p::ConvFFT) = p.analytic != nothing

outType(p::ConvFFT{D,OT}) where {D,OT} = OT
nFrames(p::ConvFFT) = size(p.weight)[end]

include("transforms.jl")
include("Utils.jl")
include("convFFTConstructors.jl")

function __init__()
  use_cuda[] = CUDA.functional() # Can be overridden after load with `Flux.use_cuda[] = false`
  if CUDA.functional()
    if !CUDA.has_cudnn()
      @warn "CUDA.jl found cuda, but did not find libcudnn. Some functionality will not be available."
end
  end
end

end # module


text:
using StratiGraphics
using GeoStatsBase
using Plots; gr()
using VisualRegressionTests
using Test, Pkg, Random

using ImageMagick

# workaround for GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

include("dummysolver.jl")

@testset "StratiGraphics.jl" begin
  if visualtests
    Random.seed!(2019)
    proc   = GeoStatsProcess(Dummy())
    env    = Environment([proc, proc], [0.5 0.5; 0.5 0.5], ExponentialDuration(1.0))
    record = simulate(env, LandState(zeros(50,50)), 10)
    strata = Strata(record)

    @plottest plot(strata) joinpath(datadir,"strata.png") !istravis 0.1

    Random.seed!(2019)
    problem = SimulationProblem(RegularGrid{Float64}(50,50,20), :strata => Float64, 3)
    solver₁ = StratSim(:strata => (environment=env,))
    solver₂ = StratSim(:strata => (environment=env,fillbase=0))
    solver₃ = StratSim(:strata => (environment=env,fillbase=0,filltop=0))
    solvers = [solver₁, solver₂, solver₃]

    solutions = [solve(problem, solver) for solver in solvers]
    snames = ["voxel1","voxel2","voxel3"]

    for (solution, sname) in zip(solutions, snames)
      reals = solution[:strata]
      @plottest begin
        plts = [heatmap(rotr90(real[1,:,:])) for real in reals]
        plot(plts..., layout=(3,1))
      end joinpath(datadir,sname*".png") !istravis
    end
  end
end

ext: jl


text:
# Mutation operators
# ==================

# Evolutionary mutations
# ======================

"""
    gaussian(x, s::IsotropicStrategy)

Performs Gaussian isotropic mutation of the recombinant `x` given the strategy `s`  by adding Gaussian noise as follows:

- ``x_i^\\prime = x_i + s.\\sigma \\mathcal{N}_i(0,1)``

"""
function gaussian(recombinant::AbstractVector, s::IsotropicStrategy)
    vals = randn(length(recombinant)) * s.σ
    recombinant += vals
    return recombinant
end

"""
    gaussian(x, s::AnisotropicStrategy)

Performs Gaussian anisotropic mutation of the recombinant `x` given the strategy `s`  by adding Gaussian noise as follows:

- ``x_i^\\prime = x_i + s.\\sigma_i \\mathcal{N}_i(0,1)``

"""
function gaussian(recombinant::AbstractVector, s::AnisotropicStrategy)
    @assert length(s.σ) == length(recombinant) "Parameter `σ` must be defined for every dimension of objective parameter"
    vals = randn(length(recombinant)) .* s.σ
    recombinant += vals
    return recombinant
end

"""
    cauchy(x, s::IsotropicStrategy)

Performs isotropic mutation of the recombinant `x` given the strategy `s`  by adding a noise from the Cauchy distribution as follows:

- ``x_i^\\prime = x_i + s.\\sigma_i \\delta_i``

where ``\\delta`` is a Cauchy random variable with the scale parameter ``t = 1`` [^2].

"""
function cauchy(recombinant::AbstractVector, s::IsotropicStrategy)
    l = length(recombinant)
    vals = s.σ * randn(l)./(randn(l).+eps())
    recombinant += vals
    return recombinant
end

# Strategy mutation operators
# ===========================

"""
    gaussian(s::IsotropicStrategy)

Performs in-place mutation of the isotropic strategy `s` modifying its mutated strategy parameter ``\\sigma`` with Gaussian noise as follows:

- ``\\sigma^\\prime = \\sigma \\exp(\\tau_0 \\mathcal{N}(0,1))``

"""
function gaussian(s::IsotropicStrategy)
    s.σ *= exp(s.τ₀*randn())
    return s
end

"""
    gaussian(s::AnisotropicStrategy)

Performs in-place mutation of the anisotropic strategy `s` modifying its mutated strategy parameter ``\\sigma`` with Gaussian noise as follows:

- ``\\sigma_i^\\prime = \\sigma_i \\exp(\\tau_0 \\mathcal{N}(0,1) + \\tau_i \\mathcal{N}(0,1))``

"""
function gaussian(s::AnisotropicStrategy)
    s.σ .*= exp.(s.τ₀*randn())*exp.(s.τ*randn(length(s.σ)))
    return s
end

# Genetic mutations
# =================

# Binary mutations
# ---------------------

"""
    flip(recombinant)

Returns an in-place mutated binary `recombinant` with a bit flips at random positions.
"""
function flip(recombinant::T) where {T <: BitVector}
    s = length(recombinant)
    pos = rand(1:s)
    recombinant[pos] = !recombinant[pos]
    return recombinant
end

"""
    bitinversion(recombinant)

Returns an in-place mutated binary `recombinant` with its bits inverted.
"""
bitinversion(recombinant::T) where {T <: BitVector} = map!(!, recombinant, recombinant)

# Real-valued mutations
# ---------------------

"""
    uniform(r = 1.0)

Returns an in-place real valued mutation function that performs the uniform distributed mutation [^1].

The mutation operator randomly chooses a number ``z`` in from the uniform distribution on the interval ``[-r,r]``, the mutation range.
The mutated individual is given by

- ``x_i^\\prime = x_i + z_i``

"""
function uniform(r::Real = 1.0)
    function mutation(recombinant::T) where {T <: AbstractVector}
        d = length(recombinant)
        recombinant += 2r.*rand(d).-r
        return recombinant
    end
    return mutation
end

"""
    gaussian(σ = 1.0)

Returns an in-place real valued mutation function that performs the normal distributed mutation [^1].

The mutation operator randomly chooses a number ``z`` in from the normal distribution ``\\mathcal{N}(0,\\sigma)`` with standard deviation ``\\sigma``.
The mutated individual is given by

- ``x_i^\\prime = x_i + z_i``

"""
function gaussian(σ::Real = 1.0)
    function mutation(recombinant::T) where {T <: AbstractVector}
        d = length(recombinant)
        recombinant += σ.*randn(d)
        return recombinant
    end
    return mutation
end

"""
    domainrange(valrange, m = 20)

Returns an in-place real valued mutation function that performs the BGA mutation scheme with the mutation range `valrange` and the mutation probability `1/m` [^1].
"""
function domainrange(valrange::Vector, m::Int = 20)
    prob = 1.0 / m
    function mutation(recombinant::T) where {T <: AbstractVector}
        d = length(recombinant)
        @assert length(valrange) == d "Range matrix must have $(d) columns"
        δ = zeros(m)
        for i in 1:length(recombinant)
            for j in 1:m
                δ[j] = (rand() < prob) ? δ[j] = 2.0^(-j) : 0.0
            end
            if rand() > 0.5
                recombinant[i] += sum(δ)*valrange[i]
            else
                recombinant[i] -= sum(δ)*valrange[i]
            end
        end
        return recombinant
    end
    return mutation
end

"""
    PM(lower, upper, p = 2)

Returns an in-place real valued mutation function that performs the Power Mutation (PM) scheme within `lower` and `upper` bound, and an index of
mutation `p`[^3].

*Note:* The implementation is a degenerate case of Mixed Integer Power Mutation ([`MIPM`](@ref))
"""
function PM(lower::Vector, upper::Vector, p::Float64 = 5.0) # index of distribution p
    return mipmmutation(lower, upper, p)
end

"""
    MIPM(lower, upper, p_real = 10, p_int = 4)

Returns an in-place real valued mutation function that performs the Mixed Integer Power Mutation (MI-PM) scheme within `lower` and `upper` bound, and an index of mutation `p_real` for real value and `p_int` for integer values[^4].
"""
function MIPM(lowerBounds::Vector, upperBounds::Vector, p_real::Float64 = 10.0, p_int::Float64 = 4.0) # index of distribution p
    return mipmmutation(lowerBounds, upperBounds, p_real, p_int)
end

function mipmmutation(lowerBounds::Vector, upperBounds::Vector, p_real::Float64, p_int::Union{Nothing, Float64} = nothing)
    function pm_mutation(x, l, u, s, d)
        x̄ = d < rand() ? x - s * (x - l) : x + s * (u - x)
        if isa(x, Integer)
            if isinteger(x̄)
                Int(x̄)
            else
                floor(Int, x̄) + (rand() > 0.5)
            end
        else
            x̄
        end
    end
    function mutation(recombinant::T) where {T <: Vector}
        d = length(recombinant)
        @assert length(lowerBounds) == d "Bounds vector must have $(d) columns"
        @assert length(upperBounds) == d "Bounds vector must have $(d) columns"
        @assert !(p_int === nothing && any(isa(x, Integer) for x in recombinant)) "Need to set p_int for integer variables"
        u = rand()
        P = (isa(x, Integer) ? p_int : p_real for x in recombinant)
        S = u .^ (1 ./ P) # random var following power distribution
        D = (recombinant - lowerBounds) ./ (upperBounds - lowerBounds)
        broadcast!(pm_mutation, recombinant, recombinant, lowerBounds, upperBounds, S, D)
        return recombinant
    end
    return mutation
end

# Combinatorial mutations (applicable to binary vectors)
# ------------------------------------------------------

"""
    inversion(recombinant)

Returns an in-place mutated individual with a random arbitrary length segment of the genome in the reverse order.
"""
function inversion(recombinant::T) where {T <: AbstractVector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    l = round(Int,(to - from)/2)
    for i in 0:(l-1)
        swap!(recombinant, from+i, to-i)
    end
    return recombinant
end

function insertion(recombinant::T) where {T <: AbstractVector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    val = recombinant[from]
    deleteat!(recombinant, from)
    return insert!(recombinant, to, val)
end

"""
    swap2(recombinant)

Returns an in-place mutated individual with a two random elements of the genome are swapped.
"""
function swap2(recombinant::T) where {T <: AbstractVector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    swap!(recombinant, from, to)
    return recombinant
end

"""
    scramble(recombinant)

Returns an in-place mutated individual with elements, on a random arbitrary length segment of the genome, been scrambled.
"""
function scramble(recombinant::T) where {T <: AbstractVector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    diff = to - from + 1
    if diff > 1
        patch = recombinant[from:to]
        idx = randperm(diff)
        # println("$(from)-$(to), P: $(patch), I: $(idx)")
        for i in 1:(diff)
            recombinant[from+i-1] = patch[idx[i]]
        end
    end
    return recombinant
end

"""
    shifting(recombinant)

Returns an in-place mutated individual with a random arbitrary length segment of the genome been shifted to an arbitrary position.
"""
function shifting(recombinant::T) where {T <: AbstractVector}
    l = length(recombinant)
    from, to, where = sort(rand(1:l, 3))
    patch = recombinant[from:to]
    diff = where - to
    if diff > 0
        # move values after tail of patch to the patch head position
        #println([from, to, where, diff])
        for i in 1:diff
            recombinant[from+i-1] = recombinant[to+i]
        end
        # place patch values in order
        start = from + diff
        for i in 1:length(patch)
            recombinant[start+i-1] = patch[i]
        end
    end
    return recombinant
end

# Differential Evolution
# ======================

"""
    differentiation(recombinant, mutators; F = 1.0)

Returns an in-place differentially mutated individual ``x^\\prime`` from `recombinant` ``x``  by `mutators` ``\\{\\xi_1, \\ldots, \\xi_n \\}`` as follows

- ``x^\\prime = x + \\sum_{i=1}^{n/2} F (\\xi_{2i-1} - \\xi_{2i})``

"""
function differentiation(recombinant::T, mutators::AbstractVector{T}; F::Real = 1.0) where {T <: AbstractVector}
    m = length(mutators)
    @assert m%2 == 0 "Must be even number of target mutators"
    for i in 1:2:m
        recombinant .+= F.*(mutators[i] .- mutators[i+1])
    end
    return recombinant
end

# Utils
# =====
function swap!(v::T, from::Int, to::Int) where {T <: AbstractVector}
    val = v[from]
    v[from] = v[to]
    v[to] = val
end

function mutationwrapper(gamutation::Function)
    wrapper(recombinant::T, s::S) where {T <: AbstractVector, S <: AbstractStrategy} =  gamutation(recombinant)
    return wrapper
end


text:
struct Binomial <: Distribution{Int} end

"""
	binom(n::Integer, p::Real)
Sample an `Int` from the Binomial distribution with parameters `n` (number of
trials) and `p` (probability of success in each trial).
"""
const binom = Binomial()

function logpdf(::Binomial, x::Integer, n::Integer, p::Real)
    Distributions.logpdf(Distributions.Binomial(n, p), x)
end

function logpdf_grad(::Binomial, x::Integer, n::Integer, p::Real)
	(nothing, nothing, (x / p - (n - x) / (1 - p)))
end

function random(::Binomial, n::Integer, p::Real)
	rand(Distributions.Binomial(n, p))
end

(::Binomial)(n, p) = random(Binomial(), n, p)

has_output_grad(::Binomial) = false
has_argument_grads(::Binomial) = (false, true)
is_discrete(::Binomial) = true

export binom

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
module ProgressiveAligner
  push!(LOAD_PATH, dirname(@__FILE__()))

  #export DataReader,
  #      DataWriter,
  #      ProfileAligner,
  #      Clustering

  include("DataReader.jl")
  include("DataWriter.jl")
  include("ProfileAligner.jl")
  include("Clustering.jl")
end # module


text:
using LaTeXStrings, Test

tst1 = L"an equation: $\alpha^2$"
tst1s = String(tst1)
@test tst1 == "an equation: \$\\alpha^2\$" == tst1s
@test L"\alpha^2" == "\$\\alpha^2\$"
@test repr("text/latex", tst1) == repr("application/x-latex", tst1) == tst1s
@test sprint(show, "text/latex", tst1) == tst1s
@test sprint(show, "text/plain", tst1) == "L\"an equation: \$\\alpha^2\$\""
@test latexstring(tst1s) == tst1 == LaTeXString(tst1s)

@test ccall(:strlen, Csize_t, (Cstring,), tst1) == ccall(:strlen, Csize_t, (Ptr{UInt8},), tst1) == sizeof(tst1)

for f in (length, eachindex, pointer, collect, ncodeunits, codeunits)
    @test f(tst1) == f(tst1.s)
end
@test nextind(tst1, 1) == prevind(tst1, 3) == 2
@test isvalid(tst1, 1)
@test findnext(isequal('o'), L"foo", 1) == 3
@test tst1[1] == 'a'
@test codeunit(tst1, 1) == UInt8('a')
@test L"foo"[1] == '$'
@test tst1[1:2] == tst1[0x01:0x02] == tst1[[1,2]] == "an"

# issue #23 — will change if #17 is addressed
@test L"x" * L"y" == "\$x\$\$y\$"

# show should return nothing
@test show(IOBuffer(), "application/x-latex", tst1) === nothing
@test show(IOBuffer(), "text/latex", tst1) === nothing

@testset "interpolation" begin
    @test L"" == latexstring("")
    @test L"%" == latexstring("%")
    @test L"$" == latexstring("\$")
    @test L"%$" == latexstring("")

    for x in ["foo", 'c', 7, 3.1]
        @test L"%$x" == latexstring(x)
        @test L"%%$(x)foo" == latexstring("%", x, "foo")
        @test L"%$(x)%$x" == latexstring(x, x)
    end
    @test L"%$(1+2)" == latexstring(3)
    @test L"%$(raw\"$$\")" == latexstring(raw"$$")

    if VERSION >= v"1.6.0-DEV.22"
        @test L"%$(@__FILE__)" == latexstring(@__FILE__)
    end

    @test_throws ErrorException getproperty(LaTeXStrings, Symbol("@L_str"))(
        LineNumberNode(@__LINE__, @__FILE__), @__MODULE__,
        "%\$(",
    )
end

using Documenter
DocMeta.setdocmeta!(LaTeXStrings, :DocTestSetup, :(using LaTeXStrings); recursive=true)
doctest(LaTeXStrings; manual = false)


text:
#  IRStab.jl: functions for stabilizing the IR flux
#             These functions are called from bc_solvers.jl to produce
#             fluxes that include dissipation terms.

include("IR_stab_diff.jl")

"""
  Computes dq/dv, where q are the conservative variables and v are the
  IR entropy variables.  This is equiavlent to calcA0 scaled by gamma_1,
  but computed from the conservative variables, which is much less expensive.

  Methods are available for 2 and 3 dimensions.
  A0 is overwritten with the result
"""
@inline function getIRA0(params::ParamType{2},
                     q::AbstractArray{Tsol,1},
                     A0::AbstractArray{Tsol, 2}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1
  p = calcPressure(params, q)

  rho = q[1]
  rhou = q[2]
  rhov = q[3]
  rhoe = q[4]

  rhoinv = 1/rho

  h = (rhoe + p)*rhoinv  # this isn't really h, but including the factor of
                         # 1/rho is convenient
  a2 = gamma*p*rhoinv  # speed of sound

  A0[1,1] = rho
  A0[2,1] = rhou
  A0[3,1] = rhov
  A0[4,1] = rhoe

  A0[1,2] = rhou
  A0[2,2] = rhou*rhou*rhoinv + p
  A0[3,2] = rhou*rhov*rhoinv
  A0[4,2] = rhou*h

  A0[1,3] = rhov
  A0[2,3] = rhou*rhov/rho
  A0[3,3] = rhov*rhov*rhoinv + p
  A0[4,3] = rhov*h

  A0[1,4] = rhoe
  A0[2,4] = h*rhou
  A0[3,4] = h*rhov
  A0[4,4] = rho*h*h - a2*p/gamma_1

  return nothing
end

@inline function getIRA0(params::ParamType{3},
                     q::AbstractArray{Tsol,1},
                     A0::AbstractArray{Tsol, 2}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1
  p = calcPressure(params, q)

  rho = q[1]
  rhou = q[2]
  rhov = q[3]
  rhow = q[4]
  rhoe = q[5]

  rhoinv = 1/rho

  h = (rhoe + p)*rhoinv
  a2 = gamma*p*rhoinv  # speed of sound

  A0[1,1] = rho
  A0[2,1] = rhou
  A0[3,1] = rhov
  A0[4,1] = rhow
  A0[5,1] = rhoe

  A0[1,2] = rhou
  A0[2,2] = rhou*rhou*rhoinv + p
  A0[3,2] = rhou*rhov*rhoinv
  A0[4,2] = rhou*rhow*rhoinv
  A0[5,2] = rhou*h

  A0[1,3] = rhov
  A0[2,3] = rhou*rhov/rho
  A0[3,3] = rhov*rhov*rhoinv + p
  A0[4,3] = rhov*rhow*rhoinv
  A0[5,3] = rhov*h

  A0[1,4] = rhow
  A0[2,4] = rhow*rhou*rhoinv
  A0[3,4] = rhow*rhov*rhoinv
  A0[4,4] = rhow*rhow*rhoinv + p
  A0[5,4] = rhow*h

  A0[1,5] = rhoe
  A0[2,5] = h*rhou
  A0[3,5] = h*rhov
  A0[4,5] = h*rhow
  A0[5,5] = rho*h*h - a2*p/gamma_1

  return nothing
end

"""
  Gets inverse matrix of [`getIRA0`](@ref)
"""
function getIRA0inv(params::ParamType{2}, q::AbstractArray{Tsol, 1},
                    A0inv::AbstractMatrix{Tsol}) where {Tsol}

  gamma = params.gamma
  gamma_1 = params.gamma_1
  gamma_1i = 1/gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2)/q[1]
  k1_dot1 = -k1/q[1]
  k1_dot2 = q[2]/q[1]
  k1_dot3 = q[3]/q[1]
#  k1_dot4 = 0

  rho_int = q[4] - k1
  rho_int_dot1 = -k1_dot1
  rho_int_dot2 = -k1_dot2
  rho_int_dot3 = -k1_dot3
  rho_int_dot4 =  1

  s = log(gamma_1*rho_int/(q[1]^gamma))
  s_dot1 = rho_int_dot1/rho_int - gamma/q[1]
  s_dot2 = rho_int_dot2/rho_int
  s_dot3 = rho_int_dot3/rho_int
  s_dot4 = rho_int_dot4/rho_int

  fac = 1.0/rho_int
  t1 = -fac*fac
  fac_dot1 = t1*rho_int_dot1
  fac_dot2 = t1*rho_int_dot2
  fac_dot3 = t1*rho_int_dot3
  fac_dot4 = t1*rho_int_dot4

  t2 = (rho_int*(gamma + 1 - s) - q[4])

  A0inv[1, 1] = (rho_int_dot1*(gamma + 1) - (rho_int*s_dot1 + s*rho_int_dot1))*fac*gamma_1i + t2*fac_dot1*gamma_1i
  A0inv[2, 1] = q[2]*fac_dot1*gamma_1i
  A0inv[3, 1] = q[3]*fac_dot1*gamma_1i
  A0inv[4, 1] = -(fac + q[1]*fac_dot1)*gamma_1i

  A0inv[1, 2] = ( (rho_int_dot2*(gamma + 1) - (rho_int*s_dot2 +
                   s*rho_int_dot2))*fac + fac_dot2*t2)*gamma_1i
  A0inv[2, 2] = (fac + q[2]*fac_dot2)*gamma_1i
  A0inv[3, 2] = q[3]*fac_dot2*gamma_1i
  A0inv[4, 2] = -(q[1]*fac_dot2)*gamma_1i

  A0inv[1, 3] = ( (rho_int_dot3*(gamma + 1) - (rho_int*s_dot3 +
                   s*rho_int_dot3))*fac + fac_dot3*t2)*gamma_1i
  A0inv[2, 3] = q[2]*fac_dot3*gamma_1i
  A0inv[3, 3] = (fac + q[3]*fac_dot3)*gamma_1i
  A0inv[4, 3] = -(q[1]*fac_dot3)*gamma_1i

  A0inv[1, 4] = ( (rho_int_dot4*(gamma + 1) - (rho_int*s_dot4 +
                   s*rho_int_dot4) - 1)*fac + fac_dot4*t2)*gamma_1i
  A0inv[2, 4] = q[2]*fac_dot4*gamma_1i
  A0inv[3, 4] = q[3]*fac_dot4*gamma_1i
  A0inv[4, 4] = -(q[1]*fac_dot4)*gamma_1i

  return nothing
end

function getIRA0inv(params::ParamType{3}, q::AbstractArray{Tsol, 1},
                    A0inv::AbstractMatrix{Tsol}) where {Tsol}

  gamma = params.gamma
  gamma_1 = params.gamma_1
  gamma_1i = 1/gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2 + q[4]^2)/q[1]
  k1_dot1 = -k1/q[1]
  k1_dot2 = q[2]/q[1]
  k1_dot3 = q[3]/q[1]
  k1_dot4 = q[4]/q[1]
#  k1_dot4 = 0

  rho_int = q[5] - k1
  rho_int_dot1 = -k1_dot1
  rho_int_dot2 = -k1_dot2
  rho_int_dot3 = -k1_dot3
  rho_int_dot4 = -k1_dot4
  rho_int_dot5 =  1

  s = log(gamma_1*rho_int/(q[1]^gamma))
  s_dot1 = rho_int_dot1/rho_int - gamma/q[1]
  s_dot2 = rho_int_dot2/rho_int
  s_dot3 = rho_int_dot3/rho_int
  s_dot4 = rho_int_dot4/rho_int
  s_dot5 = rho_int_dot5/rho_int

  fac = 1.0/rho_int
  t1 = -fac*fac
  fac_dot1 = t1*rho_int_dot1
  fac_dot2 = t1*rho_int_dot2
  fac_dot3 = t1*rho_int_dot3
  fac_dot4 = t1*rho_int_dot4
  fac_dot5 = t1*rho_int_dot5

  t2 = (rho_int*(gamma + 1 - s) - q[5])

  A0inv[1, 1] = (rho_int_dot1*(gamma + 1) - (rho_int*s_dot1 + s*rho_int_dot1))*fac*gamma_1i + t2*fac_dot1*gamma_1i
  A0inv[2, 1] = q[2]*fac_dot1*gamma_1i
  A0inv[3, 1] = q[3]*fac_dot1*gamma_1i
  A0inv[4, 1] = q[4]*fac_dot1*gamma_1i
  A0inv[5, 1] = -(fac + q[1]*fac_dot1)*gamma_1i

  A0inv[1, 2] = ( (rho_int_dot2*(gamma + 1) - (rho_int*s_dot2 +
                   s*rho_int_dot2))*fac + fac_dot2*t2)*gamma_1i
  A0inv[2, 2] = (fac + q[2]*fac_dot2)*gamma_1i
  A0inv[3, 2] = q[3]*fac_dot2*gamma_1i
  A0inv[4, 2] = q[4]*fac_dot2*gamma_1i
  A0inv[5, 2] = -(q[1]*fac_dot2)*gamma_1i

  A0inv[1, 3] = ( (rho_int_dot3*(gamma + 1) - (rho_int*s_dot3 +
                   s*rho_int_dot3))*fac + fac_dot3*t2)*gamma_1i
  A0inv[2, 3] = q[2]*fac_dot3*gamma_1i
  A0inv[3, 3] = (fac + q[3]*fac_dot3)*gamma_1i
  A0inv[4, 3] = q[4]*fac_dot3*gamma_1i
  A0inv[5, 3] = -(q[1]*fac_dot3)*gamma_1i

  A0inv[1, 4] = (rho_int_dot4*(gamma + 1) - (rho_int*s_dot4 + s*rho_int_dot4))*fac*gamma_1i + t2*fac_dot4*gamma_1i
  A0inv[2, 4] = q[2]*fac_dot4*gamma_1i
  A0inv[3, 4] = q[3]*fac_dot4*gamma_1i
  A0inv[4, 4] = (fac + q[4]*fac_dot4)*gamma_1i
  A0inv[5, 4] = -(q[1]*fac_dot4)*gamma_1i

  A0inv[1, 5] = ( (rho_int_dot5*(gamma + 1) - (rho_int*s_dot5 +
                   s*rho_int_dot5) - 1)*fac + fac_dot5*t2)*gamma_1i
  A0inv[2, 5] = q[2]*fac_dot5*gamma_1i
  A0inv[3, 5] = q[3]*fac_dot5*gamma_1i
  A0inv[4, 5] = q[4]*fac_dot5*gamma_1i
  A0inv[5, 5] = -(q[1]*fac_dot5)*gamma_1i

  return nothing
end

function getLambdaMax(params::ParamType{Tdim},
    qL::AbstractVector{Tsol}, qR::AbstractVector{Tsol},
    dir::AbstractVector{Tmsh}) where {Tsol, Tmsh, Tdim}
# compute lambda_max approximation from Carpenter's Entropy Stable Collocation
# Schemes paper

  gamma = params.gamma
  Tres = promote_type(Tsol, Tmsh)
  UnL = zero(Tres)
  UnR = zero(Tres)
  rhoLinv = 1/qL[1]
  rhoRinv = 1/qR[1]
  dA = zero(Tmsh)

  pL = calcPressure(params, qL)
  pR = calcPressure(params, qR)

  aL = sqrt(gamma*pL*rhoLinv)  # speed of sound
  aR = sqrt(gamma*pR*rhoRinv)  # speed of sound
#  Un = nx*q_avg[2]*rhoinv + ny*q_avg[3]*rhoinv + nz*q_avg[4]*rhoinv
#  dA = sqrt(nx*nx + ny*ny + nz*nz)
  for i=1:Tdim
    UnL += dir[i]*qL[i+1]*rhoLinv
    UnR += dir[i]*qR[i+1]*rhoRinv
    dA += dir[i]*dir[i]
  end

  dA = sqrt(dA)
  aL *= dA
  aR *= dA

#  println("UnL = ", UnL, ", UnR = ", UnR, ", aL = ", aL, ", aR = ", aR)

  lambda_max = 0.5*(UnL^4 + aL^4 + UnR^4 + aR^4)
  lambda_max = lambda_max^(1/4)

  return lambda_max
end

"""
  Calculates the maximum magnitude eigenvalue of the Euler flux
  jacobian at the arithmatic average of two states.

  This functions works in both 2D and 3D
  Inputs:
    params:  ParamType, conservative variable
    qL: left state
    qR: right state
    dir: direction vector (does *not* have to be unit vector)

  Outputs:
    lambda_max: eigenvalue of maximum magnitude

"""
function getLambdaMaxSimple(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol}, qR::AbstractVector{Tsol},
                      dir::AbstractVector{Tmsh}) where {Tsol, Tmsh, Tdim}
# calculate maximum eigenvalue at simple average state

  gamma = params.gamma
  Tres = promote_type(Tsol, Tmsh)

  q_avg = params.get_lambda_max_simple_data.q_avg

  for i=1:length(q_avg)
    q_avg[i] = 0.5*(qL[i] + qR[i])
  end

  return getLambdaMax(params, q_avg, dir)
end

function getLambdaMaxSimple_revm(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol}, qR::AbstractVector{Tsol},
                      dir::AbstractVector{Tmsh}, dir_bar::AbstractVector,
                      lambda_max_bar::Number
                     ) where {Tsol, Tmsh, Tdim}
# calculate maximum eigenvalue at simple average state

  gamma = params.gamma
  Tres = promote_type(Tsol, Tmsh)

  q_avg = params.get_lambda_max_simple_data.q_avg

  for i=1:length(q_avg)
    q_avg[i] = 0.5*(qL[i] + qR[i])
  end

  return getLambdaMax_revm(params, q_avg, dir, dir_bar, lambda_max_bar)
end

function getLambdaMaxSimple_revq(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol}, qL_bar::AbstractVector,
                      qR::AbstractVector{Tsol}, qR_bar::AbstractVector,
                      dir::AbstractVector{Tmsh}, lambda_max_bar) where {Tsol, Tmsh, Tdim}
# calculate maximum eigenvalue at simple average state

  gamma = params.gamma
  @unpack params.get_lambda_max_simple_data q_avg q_avg_bar
  fill!(q_avg_bar, 0)

  for i=1:length(q_avg)
    q_avg[i] = 0.5*(qL[i] + qR[i])
  end

  val = getLambdaMax_revq(params, q_avg, q_avg_bar, dir, lambda_max_bar)

  for i=1:length(q_avg)
    qL_bar[i] += 0.5*q_avg_bar[i]
    qR_bar[i] += 0.5*q_avg_bar[i]
  end

  return val
end

function getLambdaMaxRoe(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol}, qR::AbstractVector{Tsol},
                      dir::AbstractVector{Tmsh}) where {Tsol, Tmsh, Tdim}
# compute lambda_max approximation from Carpenter's Entropy Stable Collocation
# Schemes paper

  gamma = params.gamma
  Tres = promote_type(Tsol, Tmsh)
  Un = zero(Tres)
  dA = zero(Tmsh)
  rhoLinv = 1/qL[1]
  rhoRinv = 1/qR[1]

  pL = calcPressure(params, qL)
  pR = calcPressure(params, qR)

  aL = sqrt(gamma*pL*rhoLinv)  # speed of sound
  aR = sqrt(gamma*pR*rhoRinv)  # speed of sound
  a = 0.5*(aL + aR)
#  Un = nx*q_avg[2]*rhoinv + ny*q_avg[3]*rhoinv + nz*q_avg[4]*rhoinv
#  dA = sqrt(nx*nx + ny*ny + nz*nz)
  for i=1:Tdim
    Un += dir[i]*0.5*(qL[i+1]*rhoLinv + qR[i+1]*rhoLinv)
    dA += dir[i]*dir[i]
  end

  rhoA = absvalue(Un) + dA*a
  lambda1 = Un + dA*a
  lambda2 = Un - dA*a
  lambda3 = Un
  sat_Vn = convert(Tsol, 0.025)
  sat_Vl = convert(Tsol, 0.025)
  tau = 1
  lambda1 = (tau*max(absvalue(lambda1),sat_Vn *rhoA) - lambda1)
  lambda2 = (tau*max(absvalue(lambda2),sat_Vn *rhoA) - lambda2)
  lambda3 = (tau*max(absvalue(lambda3),sat_Vl *rhoA) - lambda3)
  lambda_max1 = max(absvalue(lambda1), absvalue(lambda2))
  lambda_max1 = max(absvalue(lambda_max1), absvalue(lambda3))

  return lambda_max1
end


text:
"""
    ABCDParams

A structure holding parameters for ABCD graph generator. Fields:
* w::Vector{Int32}:             a sorted in descending order list of vertex degrees
* s::Vector{Int32}:             a sorted in descending order list of cluster sizes
* μ::Union{Float64, Nothing}: mixing parameter
* ξ::Union{Float64, Nothing}: background graph fraction
* isCL::Bool:                 if `true` a Chung-Lu model is used, otherwise configuration model
* islocal::Bool:              if `true` mixing parameter restriction is cluster local, otherwise
                              it is only global

Exactly one of ξ and μ must be passed as `Float64`. Also if `ξ` is passed then
`islocal` must be `false`.

The base ABCD graph is generated when ξ is passed and `isCL` is set to `false`.
"""
struct ABCDParams
    w::Vector{Int32}
    s::Vector{Int32}
    μ::Union{Float64, Nothing}
    ξ::Union{Float64, Nothing}
    isCL::Bool
    islocal::Bool

    function ABCDParams(w, s, μ, ξ, isCL, islocal)
        length(w) == sum(s) || throw(ArgumentError("inconsistent data"))
        if !isnothing(μ)
            0 ≤ μ ≤ 1 || throw(ArgumentError("inconsistent data on μ"))
        end
        if !isnothing(ξ)
            0 ≤ ξ ≤ 1 || throw(ArgumentError("inconsistent data ξ"))
            if islocal
                throw(ArgumentError("when ξ is provided local model is not allowed"))
            end
        end
        if isnothing(μ) && isnothing(ξ)
            throw(ArgumentError("inconsistent data: either μ or ξ must be provided"))
        end

        if !(isnothing(μ) || isnothing(ξ))
            throw(ArgumentError("inconsistent data: only μ or ξ may be provided"))
        end

        new(sort(w, rev=true),
            sort(s, rev=true),
            μ, ξ, isCL, islocal)
    end
end

function randround(x)
    d = floor(Int32, x)
    d + (rand() < x - d)
end

function populate_clusters(params::ABCDParams)
    w, s = params.w, params.s
    if isnothing(params.ξ)
        mul = 1.0 - params.μ
    else
        n = length(w)
        ϕ = 1.0 - sum((sl/n)^2 for sl in s)
        mul = 1.0 - params.ξ*ϕ
    end
    @assert length(w) == sum(s)
    @assert 0 ≤ mul ≤ 1
    @assert issorted(w, rev=true)
    @assert issorted(s, rev=true)

    slots = copy(s)
    clusters = Int32[]
    j = 0
    for (i, vw) in enumerate(w)
        while j + 1 ≤ length(s) && mul * vw + 1 ≤ s[j + 1]
            j += 1
        end
        j == 0 && throw(ArgumentError("could not find a large enough cluster for vertex of weight $vw"))
        wts = Weights(view(slots, 1:j))
        wts.sum == 0 && throw(ArgumentError("could not find an empty slot for vertex of weight $vw"))
        loc = sample(1:j, wts)
        push!(clusters, loc)
        slots[loc] -= 1
    end
    clusters
end

function CL_model(clusters, params)
    @assert params.isCL
    w, s, μ = params.w, params.s, params.μ
    cluster_weight = zeros(Int32, length(s))
    for i in axes(w, 1)
        cluster_weight[clusters[i]] += w[i]
    end
    total_weight = sum(cluster_weight)
    @assert total_weight == sum(w)
    if params.islocal
        ξl = @. μ / (1.0 - cluster_weight / total_weight)
        maximum(ξl) >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
    else
        if isnothing(params.ξ)
            ξg = μ / (1.0 - sum(x -> x^2, cluster_weight) / total_weight^2)
            ξg >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
        else
            ξg = params.ξ
        end
    end

    wf = float.(w)
    edges = Set{Tuple{Int32, Int32}}()
    mutex = ReentrantLock()
    @threads for tid in 1:nthreads()
        local thr_edges = Set{Tuple{Int32, Int32}}[]
        for i in tid:nthreads():length(s)
            local local_edges = Set{Tuple{Int32, Int32}}()
            local idxᵢ = findall(==(i), clusters)
            @debug "tid:$(tid) start CL_model for i:$(i) size:$(length(idxᵢ))"
            local wᵢ = wf[idxᵢ]
            local ξ = params.islocal ? ξl[i] : ξg
            local m = randround((1-ξ) * sum(wᵢ) / 2)
            local ww = Weights(wᵢ)
            while length(local_edges) < m
                local a = sample(idxᵢ, ww, m - length(local_edges))
                local b = sample(idxᵢ, ww, m - length(local_edges))
                for (p, q) in zip(a, b)
                    p != q && push!(local_edges, minmax(p, q))
                end
            end
            push!(thr_edges, local_edges)
            @debug "tid:$(tid) end CL_model for i:$(i) size:$(length(idxᵢ))"
        end
        @debug "tid:$(tid) synch CL_model"
        lock(mutex)
        union!(edges, thr_edges...)
        unlock(mutex)
        @debug "tid:$(tid) end CL_model"
    end
    wwt = if params.islocal
        Weights([ξl[clusters[i]]*x for (i,x) in enumerate(wf)])
    else
        Weights(ξg * wf)
    end
    while 2*length(edges) < total_weight
        a = sample(axes(w, 1), wwt, randround(total_weight / 2) - length(edges))
        b = sample(axes(w, 1), wwt, randround(total_weight / 2) - length(edges))
        for (p, q) in zip(a, b)
            p != q && push!(edges, minmax(p, q))
        end
    end
    edges
end

function config_model(clusters, params)
    @assert !params.isCL
    w, s, μ = params.w, params.s, params.μ

    cluster_weight = zeros(Int32, length(s))
    for i in axes(w, 1)
        cluster_weight[clusters[i]] += w[i]
    end
    total_weight = sum(cluster_weight)
    @assert total_weight == sum(w)
    if params.islocal
        ξl = @. μ / (1.0 - cluster_weight / total_weight)
        maximum(ξl) >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
        w_internal_raw = [w[i] * (1 - ξl[clusters[i]]) for i in axes(w, 1)]
    else
        if isnothing(params. ξ)
            ξg = μ / (1.0 - sum(x -> x^2, cluster_weight) / total_weight^2)
            ξg >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
        else
            ξg = params.ξ
        end
        w_internal_raw = [w[i] * (1 - ξg) for i in axes(w, 1)]
    end

    clusterlist = [Int32[] for i in axes(s, 1)]
    for i in axes(clusters, 1)
        push!(clusterlist[clusters[i]], i)
    end
    # order by cluster size
    # idx = sortperm([length(cluster) for cluster in clusterlist], rev=true)
    # clusterlist = clusterlist[idx]

    edges::Vector{Set{Tuple{Int32, Int32}}} = []
    unresolved_collisions, length_global_recycle = 0, 0
    w_internal = Vector{Int32}(undef, length(w_internal_raw))
    mutex = ReentrantLock()
    @threads for tid in 1:nthreads()
      local thr_clusters::Vector{Vector{Int32}} = []
      local thr_weights::Vector{Vector{Int32}} = []

      for c in tid:nthreads():length(s)
        local cluster = clusterlist[c]
        local w_cluster = Vector{Int32}(undef, length(cluster))
        local maxw_idx = argmax(view(w_internal_raw, cluster))
        local wsum = 0
        for i in axes(cluster, 1)
            if i != maxw_idx
                neww = randround(w_internal_raw[cluster[i]])
                w_cluster[i] = neww
                wsum += neww
            end
        end
        local maxw = floor(Int32, w_internal_raw[cluster[maxw_idx]])
        w_cluster[maxw_idx] = maxw + (isodd(wsum) ? iseven(maxw) : isodd(maxw))

        push!(thr_clusters, cluster)
        push!(thr_weights, w_cluster)
      end
      @debug "tid $(tid) getting 1st lock"
      lock(mutex)
      @debug "tid $(tid) got 1st lock"
      foreach((cluster,w_cluster)->w_internal[cluster]=w_cluster, thr_clusters, thr_weights)
      @debug "tid $(tid) releasing 1st lock"
      unlock(mutex)
    end

    global_edges = Set{Tuple{Int32, Int32}}()
    recycle = Tuple{Int32,Int32}[]
    w_global = w - w_internal
    stubs::Vector{Int32} = Vector{Int32}(undef, sum(w_global))

    ch = Channel{Int}(1+length(s))
    foreach(cid->put!(ch, cid), 0:length(s))
    close(ch)

    @threads for tid in 1:nthreads()
        local thr_edges   = Set{Tuple{Int32, Int32}}[]
        local thr_recycle = Vector{Tuple{Int32,Int32}}[]

        for cid in ch
            local local_recycle = Tuple{Int32,Int32}[]
            if cid == 0 # global/background graph
                sizehint!(global_edges, length(stubs)>>1)
                # populate stubs
                local v::Vector{Int} = similar(w_global, Int, length(w_global)+1)
                v[1] = 0; cumsum!(view(v, 2:length(v)), w_global)
                foreach(i->stubs[v[i]+1:v[i+1]].=i, axes(w_global,1));
                @assert sum(w) == length(stubs) + sum(w_internal)
                shuffle!(stubs)

                for i in 1:2:length(stubs)
                    e = minmax(stubs[i], stubs[i+1])
                    if (e[1] == e[2]) || (e in global_edges)
                        push!(local_recycle, e)
                    else
                        push!(global_edges, e)
                    end
                end
                length_global_recycle = length(local_recycle)
                if length_global_recycle > 0
                    @warn "global collisions: $(length_global_recycle) fraction: $(2 * length_global_recycle / total_weight)"
                    @debug "dups1 are $(local_recycle)"
                end
            else # community graphs
                local cluster = clusterlist[cid]
                local w_cluster = w_internal[cluster]

                @debug "tid $(tid) cluster $(length(cluster)) w_cluster $(sum(w_cluster))"
                local local_stubs::Vector{Int32} = Vector{Int32}(undef, sum(w_cluster))
                # populate stubs
                local v = similar(w_cluster, Int, length(w_cluster)+1)
                v[1] = 0; cumsum!(view(v, 2:length(v)), w_cluster)
                foreach((i,j)->local_stubs[v[i]+1:v[i+1]].=j, axes(cluster,1), cluster);
                @assert sum(w_cluster) == length(local_stubs)

                shuffle!(local_stubs)
                local local_edges = Set{Tuple{Int32, Int32}}()
                sizehint!(local_edges, length(local_stubs)>>1)
                for i in 1:2:length(local_stubs)
                    e = minmax(local_stubs[i], local_stubs[i+1])
                    if (e[1] == e[2]) || (e in local_edges)
                        push!(local_recycle, e)
                    else
                        push!(local_edges, e)
                    end
                end
                local last_recycle = length(local_recycle)
                local recycle_counter = last_recycle
                while !isempty(local_recycle)
                    recycle_counter -= 1
                    if recycle_counter < 0
                        if length(local_recycle) < last_recycle
                            last_recycle = length(local_recycle)
                            recycle_counter = last_recycle
                        else
                            break
                        end
                    end
                    local p1 = popfirst!(local_recycle)
                    local from_recycle = 2 * length(local_recycle) / length(local_stubs)
                    local success = false
                    for _ in 1:2:length(local_stubs)
                        local p2 = if rand() < from_recycle
                            used_recycle = true
                            recycle_idx = rand(axes(local_recycle, 1))
                            local_recycle[recycle_idx]
                        else
                            used_recycle = false
                            rand(local_edges)
                        end
                        if rand() < 0.5
                            local newp1 = minmax(p1[1], p2[1])
                            local newp2 = minmax(p1[2], p2[2])
                        else
                            local newp1 = minmax(p1[1], p2[2])
                            local newp2 = minmax(p1[2], p2[1])
                        end
                        if newp1 == newp2
                            good_choice = false
                        elseif (newp1[1] == newp1[2]) || (newp1 in local_edges)
                            good_choice = false
                        elseif (newp2[1] == newp2[2]) || (newp2 in local_edges)
                            good_choice = false
                        else
                            good_choice = true
                        end
                        if good_choice
                            if used_recycle
                                local_recycle[recycle_idx], local_recycle[end] = local_recycle[end], local_recycle[recycle_idx]
                                pop!(local_recycle)
                            else
                                pop!(local_edges, p2)
                            end
                            success = true
                            push!(local_edges, newp1)
                            push!(local_edges, newp2)
                            break
                        end
                    end
                    success || push!(local_recycle, p1)
                end
                push!(thr_edges, local_edges)
            end
            push!(thr_recycle, local_recycle)
        end
        @debug "tid $(tid) getting 2nd lock"
        lock(mutex)
        @debug "tid $(tid) got 2nd lock"
        append!(edges, thr_edges)
        append!(recycle, thr_recycle...)
        @debug "tid $(tid) releasing 2nd lock"
        unlock(mutex)
    end

    unresolved_collisions = length(recycle) - length_global_recycle
    if unresolved_collisions > 0
        @warn "Unresolved cluster collisions: $(unresolved_collisions) fraction: $(2 * unresolved_collisions / total_weight)"
    end

    @debug "GLOBAL resolving dups"
    @debug "intersect $(length(global_edges)) global_edges with $(typeof(edges)) $([length(e) for e in edges])"
    dups = [Set{Tuple{Int32, Int32}}() for _ in axes(edges, 1)]
    @threads for i in axes(edges, 1)
        if length(global_edges) > length(edges[i])
            dups[i] = intersect(global_edges, edges[i])
        else
            dups[i] = intersect(edges[i], global_edges)
        end
    end
    append!(recycle, dups...)
    setdiff!(global_edges, dups...)
    dups = Nothing
    @debug "dups2 are $(recycle)"
    last_recycle = length(recycle)
    recycle_counter = last_recycle
    while !isempty(recycle)
        recycle_counter -= 1
        if recycle_counter < 0
            if length(recycle) < last_recycle
                last_recycle = length(recycle)
                recycle_counter = last_recycle
            else
                break
            end
        end
        p1 = pop!(recycle)
        from_recycle = 2 * length(recycle) / length(stubs)
        p2 = if rand() < from_recycle
            i = rand(axes(recycle, 1))
            recycle[i], recycle[end] = recycle[end], recycle[i]
            pop!(recycle)
        else
            x = rand(global_edges)
            pop!(global_edges, x)
        end
        if rand() < 0.5
            newp1 = minmax(p1[1], p2[1])
            newp2 = minmax(p1[2], p2[2])
        else
            newp1 = minmax(p1[1], p2[2])
            newp2 = minmax(p1[2], p2[1])
        end
        for newp in (newp1, newp2)
            if (newp[1] == newp[2]) || (newp in global_edges) || any(cluster->newp in cluster, edges)
                push!(recycle, newp)
            else
                push!(global_edges, newp)
            end
        end
    end
    push!(edges, global_edges)
    if isempty(recycle)
        @debug "$(length(global_edges)) global_edges $(length(stubs)) stubs"
        @assert length(global_edges) == length(stubs)/2 + unresolved_collisions
    else
        @debug "dups3 are $(recycle)"
        last_recycle = length(recycle)
        recycle_counter = last_recycle
        cw = Weights([length(c) for c in edges])
        while !isempty(recycle)
            recycle_counter -= 1
            if recycle_counter < 0
                if length(recycle) < last_recycle
                    last_recycle = length(recycle)
                    recycle_counter = last_recycle
                else
                    break
                end
            end
            p1 = pop!(recycle)
            local cluster = sample(edges, cw)
            x = rand(cluster)
            p2 = pop!(cluster, x)
            if rand() < 0.5
                newp1 = minmax(p1[1], p2[1])
                newp2 = minmax(p1[2], p2[2])
            else
                newp1 = minmax(p1[1], p2[2])
                newp2 = minmax(p1[2], p2[1])
            end
            for newp in (newp1, newp2)
                if (newp[1] == newp[2]) || any(c->newp in c, edges)
                    push!(recycle, newp)
                else
                    push!(global_edges, newp)
                end
            end
        end
        if !isempty(recycle)
            unresolved_collisions = length(recycle)
            @warn "Very hard graph! Failed to generate $(unresolved_collisions) edges; fraction: $(2 * unresolved_collisions / total_weight)"
        end
    end
    return ChainedVector([edgeset.dict.keys[edgeset.dict.slots.==0x1] for edgeset in edges])
end

"""
    gen_graph(params::ABCDParams)

Generate ABCD graph following parameters specified in `params`.

Return a named tuple containing a set of edges of the graph and a list of cluster
assignments of the vertices.
The ordering of vertices and clusters is in descending order (as in `params`).
"""
function gen_graph(params::ABCDParams)
    clusters = populate_clusters(params)
    edges = params.isCL ? CL_model(clusters, params) : config_model(clusters, params)
    (edges=edges, clusters=clusters)
end