text:
#
# error messages:
#
SOM_ERRORS = Dict(
:ERR_MPL =>
"""
    Matplotlib is not correctly installed!

    See the documentation at https://andreasdominik.github.io/SOM.jl/stable/
    for details and potential solutions.
""",
:ERR_MATRIX =>
"""
    Input data is not a numerical 2D-matrix!

    Please provide a numerical 2D-array or DataFrame
    that can be converted into an Array{Float64,2}.
""",
:ERR_COLOUR_DEF =>
"""
    Incorrect colour definition!

    Please provide a Dict of class lables and colours or
    a colour map name (as String or as Symbol).
""",
:ERR_COL_NUM =>
"""
    Wrong number of attributes!

    Number of columns of data (i.e. attributes) does not match
    number of dimensions of codes.
"""

)


text:
using Pkg; Pkg.activate(".."); Pkg.instantiate()

using OceanTurb, Printf

@use_pyplot_utils

# Instantiate a diffusion model with 100 grid points, height 1.0,
# diffusivity 0.01, and a backward Euler timestepper.
model = Diffusion.Model(
    N = 100,
    L = 1.0,
    K = 0.01,
    stepper = :BackwardEuler)

# Set the initial condition to a Gaussian cenetered at z=-0.5.
c₀(z) = exp(-(z + 0.5)^2 / 0.005)
model.solution.c = c₀

# Iterate the model forward
iterate!(model, Δt=0.01, Nt=100)

# Plot some results
fig, axs = subplots()
xlabel(L"c")
ylabel(L"z")

plot(c₀.(model.grid.zc), label=L"t=0")
gcf()

plot(model.solution.c, label=@sprintf("\$ t = %0.2f \$", time(model)))
removespines("top", "right")

legend()
gcf()

ext: jl


text:
"""
    test trajectory optimization
"""
prob_traj = prob.prob

# objective
z0 = rand(prob_traj.num_var)
tmp_o(z) = eval_objective(prob_traj, z)
∇j = zeros(prob_traj.num_var)
eval_objective_gradient!(∇j, z0, prob_traj)
@assert norm(ForwardDiff.gradient(tmp_o, z0) - ∇j) < 1.0e-10

# constraints
c0 = zeros(prob_traj.num_con)
eval_constraint!(c0, z0, prob_traj)
tmp_c(c, z) = eval_constraint!(c, z, prob_traj)
∇c_fd = ForwardDiff.jacobian(tmp_c, c0, z0)
spar = sparsity_jacobian(prob_traj)
∇c_vec = zeros(length(spar))
∇c = zeros(prob_traj.num_con, prob_traj.num_var)
eval_constraint_jacobian!(∇c_vec, z0, prob_traj)
for (i,k) in enumerate(spar)
    ∇c[k[1],k[2]] = ∇c_vec[i]
end
@assert norm(vec(∇c) - vec(∇c_fd)) < 1.0e-10
@assert sum(∇c) - sum(∇c_fd) < 1.0e-10

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
mutable struct Arrow <: Shape
    position::Tuple{Luxor.Point,Luxor.Point}
    color::Union{Blending,Coloring}
    style::Union{Symbol, Vector{Float64}}
end

function Arrow( ; position, color, style)
    return Arrow(position, color, style)
end


text:
using GeneticVariation

@testset "convert_gt" begin
    #@code_warntype convert_gt(Float64, (false, false), true, :additive)
    @inferred convert_gt(Float64, (false, false), true, :additive)
    # REF is the minor allele
    @test convert_gt(Float64, (false, false), true, :additive) == 0.0
    @test convert_gt(Float64, (false, true), true, :additive) == 1.0
    @test convert_gt(Float64, (true, false), true, :additive) == 1.0
    @test convert_gt(Float64, (true, true), true, :additive) == 2.0
    @test convert_gt(Float64, (false, false), true, :dominant) == 0.0
    @test convert_gt(Float64, (false, true), true, :dominant) == 1.0
    @test convert_gt(Float64, (true, false), true, :dominant) == 1.0
    @test convert_gt(Float64, (true, true), true, :dominant) == 1.0
    @test convert_gt(Float64, (false, false), true, :recessive) == 0.0
    @test convert_gt(Float64, (false, true), true, :recessive) == 0.0
    @test convert_gt(Float64, (true, false), true, :recessive) == 0.0
    @test convert_gt(Float64, (true, true), true, :recessive) == 1.0
    # ALT is the minor allele
    @test convert_gt(Float64, (false, false), false, :additive) == 2.0
    @test convert_gt(Float64, (false, true), false, :additive) == 1.0
    @test convert_gt(Float64, (true, false), false, :additive) == 1.0
    @test convert_gt(Float64, (true, true), false, :additive) == 0.0
    @test convert_gt(Float64, (false, false), false, :dominant) == 1.0
    @test convert_gt(Float64, (false, true), false, :dominant) == 1.0
    @test convert_gt(Float64, (true, false), false, :dominant) == 1.0
    @test convert_gt(Float64, (true, true), false, :dominant) == 0.0
    @test convert_gt(Float64, (false, false), false, :recessive) == 1.0
    @test convert_gt(Float64, (true, false), false, :recessive) == 0.0
    @test convert_gt(Float64, (false, true), false, :recessive) == 0.0
    @test convert_gt(Float64, (true, true), false, :recessive) == 0.0
end

@testset "convert_gt(vcfile)" begin
    vcffile = "test.08Jun17.d8b.vcf.gz"
    isfile(vcffile) || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        abspath(joinpath(dirname(pathof(VCFTools)), "..", "test/$vcffile")))
    #@code_warntype convert_gt(UInt8, vcffile; impute = false, center = false, scale = false)
    @inferred convert_gt(UInt8, vcffile; impute = false, center = false, scale = false)
    # convert to a matrix of UInt8, no impute/center/scale (default)
    @time A = convert_gt(UInt8, vcffile)
    @test size(A) == (191, 1356)
    @test eltype(A) == Union{UInt8, Missing}
    # convert to a matrix of Float64, impute = center = scale = true
    @time B = convert_gt(Float64, vcffile; as_minorallele = true, impute = true, center = true, scale = true)
    @test size(B) == (191, 1356)
    @test eltype(B) == Union{Float64, Missing}
    # convert to a matrix of Float32 without checking minor alleles (reading 0/1 as-is)
    # @code_warntype convert_gt(Float32, vcffile; as_minorallele = false)
    @time C = convert_gt(Float32, vcffile; as_minorallele = false)
    @test size(C) == (191, 1356)
    @test eltype(C) == Union{Float32, Missing}
    # read the first record to a UInt8 vector, no impute/center/scale (default)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    A1 = Vector{Union{UInt8, Missing}}(undef, nsamples(vcffile))
    @time copy_gt!(A1, reader)
    @test all(A1[findall(!ismissing, A1)] .== A[findall(!ismissing, A[:, 1])])
    @test all(findall(ismissing, A1) .== findall(ismissing, A[:, 1]))
    # convert next 5 records to a Float64 matrix, impute = center = scale = true
    B5 = Matrix{Union{Float64, Missing}}(undef, nsamples(vcffile), 5)
    @time copy_gt!(B5, reader; impute = true, center = true, scale = true)
    B_copy = copy(B[:, 2:6])
    @test all(B5[findall(!ismissing, B5)] .== B_copy[findall(!ismissing, B[:, 2:6])])
    @test all(B5[findall(ismissing, B5)] .== B_copy[findall(ismissing, B_copy)])
    # convert next 5 records to a Float32 matrix without checking which is minor alleles
    C5 = Matrix{Union{Float64, Missing}}(undef, nsamples(vcffile), 5)
    @time copy_gt_as_is!(C5, reader)
    @test findall(x -> x === missing, C5) == []
    close(reader)
end

@testset "convert_ht" begin
    # REF is the minor allele
    @test convert_ht(Float64, true, true) == 1.0
    @test convert_ht(Float64, false, true) == 0.0
    # ALT is the minor allele
    @test convert_ht(Float64, true, false) == 0.0
    @test convert_ht(Float64, false, false) == 1.0
end

@testset "convert_ht(vcfile)" begin
    vcffile = "test.08Jun17.d8b.vcf.gz"
    # convert to a reference haplotype panel of Float32, no checking for minor allele
    @time H1 = convert_ht(Float32, vcffile, as_minorallele = false)
    @test size(H1) == (382, 1356)
    @test typeof(H1) == Matrix{Float32}
    # convert to a reference haplotype panel of Float64, checking for minor allele
    @time H2 = convert_ht(Float64, vcffile, as_minorallele = true)
    @test size(H2) == (382, 1356)
    @test typeof(H2) == Matrix{Float64}
    # convert first record into 2 haplotype vectors and check their sum is the genotype vector
    reader = VCF.Reader(openvcf(vcffile, "r"))
    h1h2 = Matrix{Float64}(undef, 2, nrecords(vcffile))
    copy_ht_as_is!(h1h2, reader)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    g1   = Matrix{Union{Float64, Missing}}(undef, 1, nrecords(vcffile))
    copy_gt_as_is!(g1, reader)
    @test all(sum(h1h2, dims=1) .== g1)
end

@testset "convert_ds" begin
    vcffile = "test.08Jun17.d8b.vcf.gz"
    D = convert_ds(Float64, vcffile, key = "DS", impute=false, center=false, scale=false)
    @test D[1, 5] == 1.0
    @test D[1, 1345] == 0.05
    @test D[2, 13] == 1.75
    @test eltype(D) == Union{Missing, Float64}
    @test size(D) == (191, 1356)
end


text:
using MomentTensors, Test

@testset "I/O" begin
    @testset "GlobalCMT" begin
        let (m, e) = ndk("""
            PDE  2000/01/01 05:24:35.3  36.87   69.95  54.3 5.2 4.5 HINDU KUSH REGION
            B010100A         B: 13   21  45 S:  0    0   0 M:  0    0   0 CMT: 1 BOXHD:  1.0
            CENTROID:      4.0 1.0  37.15 0.13   69.72 0.14  84.6  7.2 FREE O-00000000000000
            23  2.495 0.595 -1.952 0.880 -0.544 0.633 -0.154 0.774  0.416 0.840 -5.825 1.128
            V10   4.700 11 228   2.420 79  43  -7.120  1 138   5.910 273 82  173   4 83    8
            """)
            @test m ≈ MT(1e16.*[2.495, -1.952, -0.544, -0.154, 0.416, -5.825]) atol=1e13
            @test e ≈ MT(1e16.*[0.595,  0.880,  0.633,  0.774, 0.840,  1.128]) atol=1e13
        end
        # Wrong number of columns on fourth line
        @test_throws ArgumentError ndk("""
            PDE  2000/01/01 05:24:35.3  36.87   69.95  54.3 5.2 4.5 HINDU KUSH REGION
            B010100A         B: 13   21  45 S:  0    0   0 M:  0    0   0 CMT: 1 BOXHD:  1.0
            CENTROID:      4.0 1.0  37.15 0.13   69.72 0.14  84.6  7.2 FREE O-00000000000000
            23  2.495 0.595 -1.952 0.880 -0.544 0.633 -0.154 0.774  0.416 0.840 -5.825
            V10   4.700 11 228   2.420 79  43  -7.120  1 138   5.910 273 82  173   4 83    8
            """)
        # Incorrect number of lines
        @test_throws ArgumentError ndk("""
            B010100A         B: 13   21  45 S:  0    0   0 M:  0    0   0 CMT: 1 BOXHD:  1.0
            CENTROID:      4.0 1.0  37.15 0.13   69.72 0.14  84.6  7.2 FREE O-00000000000000
            23  2.495 0.595 -1.952 0.880 -0.544 0.633 -0.154 0.774  0.416 0.840 -5.825 1.128
            V10   4.700 11 228   2.420 79  43  -7.120  1 138   5.910 273 82  173   4 83    8
            """)
        # Exponent is not an integer
        @test_throws ArgumentError ndk("""
            PDE  2000/01/01 05:24:35.3  36.87   69.95  54.3 5.2 4.5 HINDU KUSH REGION
            B010100A         B: 13   21  45 S:  0    0   0 M:  0    0   0 CMT: 1 BOXHD:  1.0
            CENTROID:      4.0 1.0  37.15 0.13   69.72 0.14  84.6  7.2 FREE O-00000000000000
            23.0  2.495 0.595 -1.952 0.880 -0.544 0.633 -0.154 0.774  0.416 0.840 -5.825 1.128
            V10   4.700 11 228   2.420 79  43  -7.120  1 138   5.910 273 82  173   4 83    8
            """)
    end

    @testset "CMTSOLUTION" begin
        @test cmtsolution("""
            PDE 1994  6  9  0 33 16.40 -13.8300  -67.5600 637.0 6.9 6.8 NORTHERN BOLIVIA
            event name:     060994A
            time shift:     29.0000
            half duration:  20.0000
            latitude:      -13.8200
            longitude:     -67.2500
            depth:         647.1000
            Mrr:      -7.590000d+27
            Mtt:       7.750000d+27
            Mpp:      -1.600000e+26
            Mrt:      -2.503000e+28
            Mrp:       4.200000e+26
            Mtp:      -2.480000e+27
            """) ≈ MT(-7.59e20, 7.75e20, -1.6e19, -2.503e21, 4.2e19, -2.48e20) atol=1e14
        # Missing a line
        @test_throws ArgumentError cmtsolution("""
            PDE 1994  6  9  0 33 16.40 -13.8300  -67.5600 637.0 6.9 6.8 NORTHERN BOLIVIA
            event name:     060994A
            time shift:     29.0000
            half duration:  20.0000
            latitude:      -13.8200
            longitude:     -67.2500
            depth:         647.1000
            Mrr:      -7.590000e+27
            Mtt:       7.750000e+27
            Mpp:      -1.600000e+26
            Mrt:      -2.503000e+28
            Mrp:       4.200000e+26
            """)
    end
end


text:
function max_independent_set(g::Graph)

    m = Model(solver=CbcSolver())

    @variable(m, b[v in vertices(g)], Bin)

    @constraint(m , nonadjacency[e in edges(g)], b[src(e)] + b[dst(e)] <= 1)

    @objective(m, Max, sum{b[v] , v in vertices(g)})

    solve(m)
    println("The max indepedent set has $(getobjectivevalue(m)) elements")
    getvalue(b)
end

function min_dominating_set(g::Graph)

    m = Model(solver=CbcSolver())

    @variable(m, b[v in vertices(g)], Bin)

    @constraint(m , dominance[u in vertices(g)], b[u] + sum{b[v], v in out_neighbors(g,u)} >= 1)

    @objective(m, Min, sum{b[v] , v in vertices(g)})

    solve(m)
    println("The min dominating set has $(getobjectivevalue(m)) elements")
    getvalue(b)
end

function min_vertex_cover(g::Graph)

    m = Model(solver=CbcSolver())

    @variable(m, b[v in vertices(g)], Bin)

    @constraint(m , cover[e in edges(g)], b[src(e)] + b[dst(e)] >= 1)

    @objective(m, Min, sum{b[v] , v in vertices(g)})

    solve(m)
    println("The min vertex cover has $(getobjectivevalue(m)) elements")
    getvalue(b)
end

function bipartite_set(g::Graph)

    m = Model(solver=CbcSolver())

    @variable(m, b[v in vertices(g)], Bin)

    @constraint(m , cover[e in edges(g)], b[src(e)] + b[dst(e)] == 1)

    solve(m)

    print("These vertices are independent : ")
    for(v in vertices(g))
        getvalue(b[v]) == 1 && print("$v, ")
    end
    println()

    getvalue(b)
end

function distances(g::Graph, root::Int)

    m = Model(solver=CbcSolver())

    @variable(m, d[v in vertices(g)], Int)

    @constraint(m , d[root] == 0)
    @constraint(m , hop[e in edges(g)], -1 <= d[src(e)] - d[dst(e)] <= 1)

    @objective(m, Max, sum{d[v] , v in vertices(g)})

    solve(m)
    getvalue(d)
end


text:
using Literate

Literate.notebook("QCBM.jl", joinpath(@__DIR__, "../notebooks"), execute=false)


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

using Markdown
using InteractiveUtils

# ╔═╡ 004d3056-5ee4-11ec-12ac-750b7ebee606
begin
	using Plots
	x=1:10
	y=rand(10)
	plot(x,y,title="first Plots demo")
end

# ╔═╡ 665707df-cbb3-4570-a4c7-c8bb98c1fdc9
begin
	using Plots;
	plot(rand(100, 4), layout = 4)
end

# ╔═╡ Cell order:
# ╠═004d3056-5ee4-11ec-12ac-750b7ebee606
# ╠═665707df-cbb3-4570-a4c7-c8bb98c1fdc9

max_stars_repo_name: BioMedAI/Pluto-Note

max_issues_repo_name: BioMedAI/Pluto-Note

max_forks_repo_name: BioMedAI/Pluto-Note


text:
# Redistribute data for load balancing after partitioning the mesh
function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations,
                            dg::DGSEM, cache, old_mpi_ranks_per_cell)
  if cache.elements.cell_ids == local_leaf_cells(mesh.tree)
    # Cell ids of the current elements are the same as the local leaf cells of the
    # newly partitioned mesh, so the solver doesn't need to be rebalanced on this rank.
    return
  end

  # Retain current solution data
  old_n_elements = nelements(dg, cache)
  old_cell_ids = copy(cache.elements.cell_ids)
  old_u_ode = copy(u_ode)
  GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
    # Use `wrap_array_native` instead of `wrap_array` since MPI might not interact
    # nicely with non-base array types
    old_u = wrap_array_native(old_u_ode, mesh, equations, dg, cache)

    @trixi_timeit timer() "reinitialize data structures" reinitialize_containers!(mesh, equations, dg, cache)

    resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
    u = wrap_array_native(u_ode, mesh, equations, dg, cache)

    # Get new list of leaf cells
    leaf_cell_ids = local_leaf_cells(mesh.tree)

    @trixi_timeit timer() "exchange data" begin
      # Collect MPI requests for MPI_Waitall
      requests = Vector{MPI.Request}()

      # Find elements that will change their rank and send their data to the new rank
      for old_element_id in 1:old_n_elements
        cell_id = old_cell_ids[old_element_id]
        if !(cell_id in leaf_cell_ids)
          # Send element data to new rank, use cell_id as tag (non-blocking)
          dest = mesh.tree.mpi_ranks[cell_id]
          request = MPI.Isend(@view(old_u[:, .., old_element_id]), dest, cell_id, mpi_comm())
          push!(requests, request)
        end
      end

      # Loop over all elements in new container and either copy them from old container
      # or receive them with MPI
      for element in eachelement(dg, cache)
        cell_id = cache.elements.cell_ids[element]
        if cell_id in old_cell_ids
          old_element_id = searchsortedfirst(old_cell_ids, cell_id)
          # Copy old element data to new element container
          @views u[:, .., element] .= old_u[:, .., old_element_id]
        else
          # Receive old element data
          src = old_mpi_ranks_per_cell[cell_id]
          request = MPI.Irecv!(@view(u[:, .., element]), src, cell_id, mpi_comm())
          push!(requests, request)
        end
      end

      # Wait for all non-blocking MPI send/receive operations to finish
      MPI.Waitall!(requests)
    end
  end # GC.@preserve old_u_ode
end

function reinitialize_containers!(mesh::TreeMesh{2}, equations, dg::DGSEM, cache)
  # Get new list of leaf cells
  leaf_cell_ids = local_leaf_cells(mesh.tree)

  # re-initialize elements container
  @unpack elements = cache
  resize!(elements, length(leaf_cell_ids))
  init_elements!(elements, leaf_cell_ids, mesh, dg.basis)

  # re-initialize interfaces container
  @unpack interfaces = cache
  resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids))
  init_interfaces!(interfaces, elements, mesh)

  # re-initialize boundaries container
  @unpack boundaries = cache
  resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids))
  init_boundaries!(boundaries, elements, mesh)

  # re-initialize mortars container
  @unpack mortars = cache
  resize!(mortars, count_required_mortars(mesh, leaf_cell_ids))
  init_mortars!(mortars, elements, mesh)

  if mpi_isparallel()
    # re-initialize mpi_interfaces container
    @unpack mpi_interfaces = cache
    resize!(mpi_interfaces, count_required_mpi_interfaces(mesh, leaf_cell_ids))
    init_mpi_interfaces!(mpi_interfaces, elements, mesh)

    # re-initialize mpi cache
    @unpack mpi_cache = cache
    init_mpi_cache!(mpi_cache, mesh, elements, mpi_interfaces, nvariables(equations), nnodes(dg))
  end
end

# Refine elements in the DG solver based on a list of cell_ids that should be refined
function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}},
                 equations, dg::DGSEM, cache, elements_to_refine)
  # Return early if there is nothing to do
  if isempty(elements_to_refine)
    return
  end

  # Determine for each existing element whether it needs to be refined
  needs_refinement = falses(nelements(dg, cache))
  needs_refinement[elements_to_refine] .= true

  # Retain current solution data
  old_n_elements = nelements(dg, cache)
  old_u_ode = copy(u_ode)
  GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
    old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)

    reinitialize_containers!(mesh, equations, dg, cache)

    resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
    u = wrap_array(u_ode, mesh, equations, dg, cache)

    # Loop over all elements in old container and either copy them or refine them
    element_id = 1
    for old_element_id in 1:old_n_elements
      if needs_refinement[old_element_id]
        # Refine element and store solution directly in new data structure
        refine_element!(u, element_id, old_u, old_element_id,
                        adaptor, equations, dg)
        element_id += 2^ndims(mesh)
      else
        # Copy old element data to new element container
        @views u[:, .., element_id] .= old_u[:, .., old_element_id]
        element_id += 1
      end
    end
    # If everything is correct, we should have processed all elements.
    # Depending on whether the last element processed above had to be refined or not,
    # the counter `element_id` can have two different values at the end.
    @assert element_id == nelements(dg, cache) + 1 || element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
  end # GC.@preserve old_u_ode

  # Sanity check
  if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && !mpi_isparallel()
    @assert ninterfaces(cache.interfaces) == ndims(mesh) * nelements(dg, cache) ("For $(ndims(mesh))D and periodic domains and conforming elements, the number of interfaces must be $(ndims(mesh)) times the number of elements")
  end

  return nothing
end

# TODO: Taal compare performance of different implementations
# Refine solution data u for an element, using L2 projection (interpolation)
function refine_element!(u::AbstractArray{<:Any,4}, element_id,
                         old_u, old_element_id,
                         adaptor::LobattoLegendreAdaptorL2, equations, dg)
  @unpack forward_upper, forward_lower = adaptor

  # Store new element ids
  lower_left_id  = element_id
  lower_right_id = element_id + 1
  upper_left_id  = element_id + 2
  upper_right_id = element_id + 3

  @boundscheck begin
    @assert old_element_id >= 1
    @assert size(old_u, 1) == nvariables(equations)
    @assert size(old_u, 2) == nnodes(dg)
    @assert size(old_u, 3) == nnodes(dg)
    @assert size(old_u, 4) >= old_element_id
    @assert     element_id >= 1
    @assert size(    u, 1) == nvariables(equations)
    @assert size(    u, 2) == nnodes(dg)
    @assert size(    u, 3) == nnodes(dg)
    @assert size(    u, 4) >= element_id + 3
  end

  # Interpolate to lower left element
  for j in eachnode(dg), i in eachnode(dg)
    acc = zero(get_node_vars(u, equations, dg, i, j, element_id))
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, old_element_id) * forward_lower[i, k] * forward_lower[j, l]
    end
    set_node_vars!(u, acc, equations, dg, i, j, lower_left_id)
  end

  # Interpolate to lower right element
  for j in eachnode(dg), i in eachnode(dg)
    acc = zero(get_node_vars(u, equations, dg, i, j, element_id))
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, old_element_id) * forward_upper[i, k] * forward_lower[j, l]
    end
    set_node_vars!(u, acc, equations, dg, i, j, lower_right_id)
  end

  # Interpolate to upper left element
  for j in eachnode(dg), i in eachnode(dg)
    acc = zero(get_node_vars(u, equations, dg, i, j, element_id))
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, old_element_id) * forward_lower[i, k] * forward_upper[j, l]
    end
    set_node_vars!(u, acc, equations, dg, i, j, upper_left_id)
  end

  # Interpolate to upper right element
  for j in eachnode(dg), i in eachnode(dg)
    acc = zero(get_node_vars(u, equations, dg, i, j, element_id))
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, old_element_id) * forward_upper[i, k] * forward_upper[j, l]
    end
    set_node_vars!(u, acc, equations, dg, i, j, upper_right_id)
  end

  return nothing
end

# Coarsen elements in the DG solver based on a list of cell_ids that should be removed
function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}},
                  equations, dg::DGSEM, cache, elements_to_remove)
  # Return early if there is nothing to do
  if isempty(elements_to_remove)
    return
  end

  # Determine for each old element whether it needs to be removed
  to_be_removed = falses(nelements(dg, cache))
  to_be_removed[elements_to_remove] .= true

  # Retain current solution data
  old_n_elements = nelements(dg, cache)
  old_u_ode = copy(u_ode)
  GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
    old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)

    reinitialize_containers!(mesh, equations, dg, cache)

    resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
    u = wrap_array(u_ode, mesh, equations, dg, cache)

    # Loop over all elements in old container and either copy them or coarsen them
    skip = 0
    element_id = 1
    for old_element_id in 1:old_n_elements
      # If skip is non-zero, we just coarsened 2^ndims elements and need to omit the following elements
      if skip > 0
        skip -= 1
        continue
      end

      if to_be_removed[old_element_id]
        # If an element is to be removed, sanity check if the following elements
        # are also marked - otherwise there would be an error in the way the
        # cells/elements are sorted
        @assert all(to_be_removed[old_element_id:(old_element_id+2^ndims(mesh)-1)]) "bad cell/element order"

        # Coarsen elements and store solution directly in new data structure
        coarsen_elements!(u, element_id, old_u, old_element_id,
                          adaptor, equations, dg)
        element_id += 1
        skip = 2^ndims(mesh) - 1
      else
        # Copy old element data to new element container
        @views u[:, .., element_id] .= old_u[:, .., old_element_id]
        element_id += 1
      end
    end
    # If everything is correct, we should have processed all elements.
    @assert element_id == nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
  end # GC.@preserve old_u_ode

  # Sanity check
  if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && !mpi_isparallel()
    @assert ninterfaces(cache.interfaces) == ndims(mesh) * nelements(dg, cache) ("For $(ndims(mesh))D and periodic domains and conforming elements, the number of interfaces must be $(ndims(mesh)) times the number of elements")
  end

  return nothing
end

# TODO: Taal compare performance of different implementations
# Coarsen solution data u for four elements, using L2 projection
function coarsen_elements!(u::AbstractArray{<:Any,4}, element_id,
                           old_u, old_element_id,
                           adaptor::LobattoLegendreAdaptorL2, equations, dg)
  @unpack reverse_upper, reverse_lower = adaptor

  # Store old element ids
  lower_left_id  = old_element_id
  lower_right_id = old_element_id + 1
  upper_left_id  = old_element_id + 2
  upper_right_id = old_element_id + 3

  @boundscheck begin
    @assert old_element_id >= 1
    @assert size(old_u, 1) == nvariables(equations)
    @assert size(old_u, 2) == nnodes(dg)
    @assert size(old_u, 3) == nnodes(dg)
    @assert size(old_u, 4) >= old_element_id + 3
    @assert     element_id >= 1
    @assert size(    u, 1) == nvariables(equations)
    @assert size(    u, 2) == nnodes(dg)
    @assert size(    u, 3) == nnodes(dg)
    @assert size(    u, 4) >= element_id
  end

  for j in eachnode(dg), i in eachnode(dg)
    acc = zero(get_node_vars(u, equations, dg, i, j, element_id))

    # Project from lower left element
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, lower_left_id) * reverse_lower[i, k] * reverse_lower[j, l]
    end

    # Project from lower right element
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, lower_right_id) * reverse_upper[i, k] * reverse_lower[j, l]
    end

    # Project from upper left element
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, upper_left_id) * reverse_lower[i, k] * reverse_upper[j, l]
    end

    # Project from upper right element
    for l in eachnode(dg), k in eachnode(dg)
      acc += get_node_vars(old_u, equations, dg, k, l, upper_right_id) * reverse_upper[i, k] * reverse_upper[j, l]
    end

    # Update value
    set_node_vars!(u, acc, equations, dg, i, j, element_id)
  end
end

# this method is called when an `ControllerThreeLevel` is constructed
function create_cache(::Type{ControllerThreeLevel}, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DG, cache)

  controller_value = Vector{Int}(undef, nelements(dg, cache))
  return (; controller_value)
end

function create_cache(::Type{ControllerThreeLevelCombined}, mesh::TreeMesh{2}, equations, dg::DG, cache)

  controller_value = Vector{Int}(undef, nelements(dg, cache))
  return (; controller_value)
end


text:
using Plots
pyplot()
using Random

include("hexbinplots_orig.jl")

function shapecoords(hx::HexBinPlots.HexHistogram)
    h,vh = HexBinPlots.make_shapes(hx)
    vmax = maximum(vh)
    xsh = Vector{Float64}()
    ysh = Vector{Float64}()
    for k in eachindex(h)
        append!(xsh,h[k].x)
        push!(xsh,h[k].x[1])
        push!(xsh,NaN)
        append!(ysh,h[k].y)
        push!(ysh,h[k].y[1])
        push!(ysh,NaN)
    end
    return (xsh, ysh, vh)
end

using Interpolations
function control_point(xi, xj, yi, yj, dist_from_mid)
    xmid = 0.5 * (xi+xj)
    ymid = 0.5 * (yi+yj)
    # get the angle of y relative to x
    theta = atan((yj-yi) / (xj-xi)) + 0.5pi
    # dist = sqrt((xj-xi)^2 + (yj-yi)^2)
    # dist_from_mid = curvature_scalar * 0.5dist
    # now we have polar coords, we can compute the position, adding to the midpoint
    (xmid + dist_from_mid * cos(theta),
     ymid + dist_from_mid * sin(theta))
end

function _add_edge_points!(xlines,ylines,xsi,ysi,xdi,ydi,curve_amount,interplen)
    xpt, ypt = control_point(xsi,xdi,ysi,ydi,curve_amount)
    xpts = [xsi, xpt, xdi]
    ypts = [ysi, ypt, ydi]
    t = range(0, stop=1, length=3)
    A = hcat(xpts, ypts)
    itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)
    tfine = range(0, 1, length=interplen)
    for t in tfine
        push!(xlines, itp(t, 1))
    end
    for t in tfine
        push!(ylines, itp(t, 2))
    end
end
# Mutating version from HexBinPlots
function _xy2counts!(counts::Dict{(Tuple{Int, Int}), Int},
                     x::AbstractArray,y::AbstractArray,
                     xsize::Real,ysize::Real,x0,y0)
    @inbounds for i in eachindex(x)
        h = convert(HexagonOffsetOddR,
                    cube_round(x[i] - x0,y[i] - y0,xsize, ysize))
        idx = (h.q, h.r)
        counts[idx] = 1 + get(counts,idx,0)
    end
    counts
end

using Statistics
using Hexagons

# Generate counts for a graph edge plot histogram
# hhparam = (xsize,ysize,x0,y0) are the hex-histogrma params
function graphedgehist(src,dst,x,y,ptscale,hhparams,curvature = 0.05)
    xlines = similar(x,0)
    ylines = similar(y,0)
    sizehint!(xlines,max(length(src)+1000,10000000)) # we'll process them in batches this big
    sizehint!(ylines,max(length(src)+1000,10000000))
    # this is the counts of the edge vectors...
    edgecounts = Dict{(Tuple{Int, Int}), Int}()
    # now, run through each edges...
    for ei in eachindex(src)
        si = src[ei]
        di = dst[ei]
        xsi = x[si]
        ysi = y[si]
        xdi = x[di]
        ydi = y[di]
        curve_amount = curvature*sign((di-si))
        edgelen = sqrt((xdi-xsi)^2 + (ydi-ysi)^2)
        interplen = max(ceil(Int,2*edgelen/ptscale),2)
        if edgelen == 0
            continue
        end
        _add_edge_points!(xlines, ylines, xsi, ysi, xdi, ydi, curve_amount, interplen)
        if length(xlines) >= 9500000
            # compress all to edgecounts
            println("Compressing counts... $(ei) / $(length(src))")
            _xy2counts!(edgecounts, xlines, ylines, hhparams...)
            resize!(xlines,0)
            resize!(ylines,0)
        end
    end
    # handle final update
    return _xy2counts!(edgecounts, xlines, ylines, hhparams...)
end

function hexbingraphplot(src,dst,x,y;nbins=500)
    xbins = nbins
    ybins = nbins
    # compute params for hex histogram bins
    xmin,xmax = extrema(x)
    ymin,ymax = extrema(y)
    xspan, yspan = xmax - xmin, ymax - ymin
    xsize, ysize = xspan / xbins, yspan / ybins
    x0, y0 = xmin - xspan / 2,ymin - yspan / 2
    hhparams = (xsize, ysize, x0, y0) # save them as a tuple
    edgescale = 6
    hsize = min((xmax-xmin)/(xbins*edgescale),(ymax-ymin)/(ybins*edgescale))
    # need to compute hhparams
    # get coords for edges
    counts = graphedgehist(src,dst,x,y,hsize,hhparams)
    # do this for the nodes...
    #counts = HexBinPlots.xy2counts_(x,y,hhparams...)
    xh,yh,vh = HexBinPlots.counts2xy_(counts,hhparams...)
    hedges = HexBinPlots.HexHistogram(xh,yh,vh,xsize,ysize,false)
    return hedges
end

function generate_nice_plot(xy, src, dst, imagefilename;
                            labels=[],
                            mymarkeralpha=0.1,
                            mymarkersize=5,
                            invpermval=false,
                            colorscale=[],
                            skipedges=false)

    if invpermval
        x = xy[:,1]
        y = xy[:,2]
        x = invperm(sortperm(x))
        y = invperm(sortperm(y))
        xy = Float64.(hcat(x,y))
    end

    if !isempty(labels)
	labels = Int.(labels)
	ulabels = unique(labels)
	d = Dict(sort(ulabels) .=> 1:length(ulabels))
	color_palette = distinguishable_colors(length(ulabels))
        color_palette[end] = RGB(1,1,1);
	colorstouse = map(i->color_palette[d[labels[i]]],1:length(labels))
    end

    xd = xy[:,1]
    yd = xy[:,2]

    plot()
    if !skipedges
        hedges = hexbingraphplot(src, dst, xd, yd, nbins=900)
        xsh, ysh, axform = shapecoords(hedges)
        avals = axform/maximum(axform)
        avals = map(x -> x <= 1/4 ? sqrt(x) : ((sqrt(1/4)-1/4) + x), avals )
        avals = 1/1.25.*avals
        plot!(xsh, ysh,
              seriestype=:shape,
              alpha=avals,
              fillcolor=:darkblue,
              legend=false,
              linealpha=0,
              framestyle=:none,
              background=nothing)
    end

    #cscale = colorscale ./ maximum(colorscale)
    #rgbvals = get(ColorSchemes.valentine,cscale)
    #ids = findall(colorscale.!=0)
    ids = collect(length(xd):-1:1)
    scatter!(xd[ids], yd[ids],
             color=colorstouse[ids],
             alpha=mymarkeralpha,
             markersize=mymarkersize,
             legend=false,
             colorbar=false,
             markerstrokewidth=0,
             framestyle=:none,
             background=nothing)

    #plot!(dpi=300,size=(800,800))
    #plot!(dpi=200,size=(600,600))
    savefig(imagefilename)
end

function generate_nice_biplot(xy1, xy2, src, dst, imagefilename;
                              labels=[],
                              mymarkeralpha1=0.1, mymarkeralpha2=0.1,
                              mymarkersize1=5, mymarkersize2=5,
                              incledges=true,
                              inclsecond=true)

    xy = [xy1; xy2]
    xd = xy[:,1]
    yd = xy[:,2]

    plot()
    if incledges
        hedges = hexbingraphplot(src, dst, xd, yd, nbins=900)
        xsh, ysh, axform = shapecoords(hedges)
        avals = axform/maximum(axform)
        avals = map(x -> x <= 1/4 ? sqrt(x) : ((sqrt(1/4)-1/4) + x), avals )
        avals = 1/1.25.*avals
        plot!(xsh, ysh,
              seriestype=:shape,
              alpha=avals,
              fillcolor=:darkblue,
              legend=false,
              linealpha=0,
              framestyle=:none,
              background=nothing)
    end

    if inclsecond
        scatter!(xy2[:, 1], xy2[:, 2],
                 color=RGB(1,1,1),
                 alpha=mymarkeralpha2,
                 markersize=mymarkersize2,
                 legend=false,
                 colorbar=false,
                 framestyle=:none,
                 background=nothing)

    end

    function get_colors()
        if isempty(labels)
            return :red
        end
	labels = Int.(labels)
	ulabels = unique(labels)
	d = Dict(sort(ulabels) .=> 1:length(ulabels))
        # skip first in palette
	color_palette = distinguishable_colors(length(ulabels) + 1)
	return map(i->color_palette[d[labels[i]] + 1], 1:length(labels))
    end
    colorstouse = get_colors()
    scatter!(xy1[:, 1], xy1[:, 2],
             color=colorstouse,
             alpha=mymarkeralpha1,
             markersize=mymarkersize1,
             legend=false,
             colorbar=false,
             framestyle=:none,
             background=nothing)

    plot!(dpi=300,size=(800,800))
    #plot!(dpi=200,size=(400,400))
    savefig(imagefilename)
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
using FASTX
using DataFrames
using CSV
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--output", "-o"
        help = "Path to the tsv file containing a list of all the human interactions (conservative) from the BIOGRID database."
        arg_type = String
        required = true
    "--biogrid_version", "-b"
        help = "Version of BIOGRID to download (default: 4.3.196)"
        arg_type = String
        required = false
        default = "4.3.196"
end

args = parse_args(ARGS, s)

DESTINATION = args["output"]
BIOGRID_VERSION = args["biogrid_version"]
BIOGRID_URL = "https://downloads.thebiogrid.org/Download/BioGRID/Release-Archive/BIOGRID-$BIOGRID_VERSION/BIOGRID-ALL-$BIOGRID_VERSION.tab3.zip"

# Download the interactions only if not already available
run(`curl $BIOGRID_URL -o biogrid.tab3.zip`)

# Unzip
run(`unzip biogrid.tab3.zip`)
run(`rm biogrid.tab3.zip`)

#setenv(`unzip biogrid.tab3.zip`, dir="/tmp")

# Load the interaction data
interactions_df = DataFrame(CSV.File("BIOGRID-ALL-$BIOGRID_VERSION.tab3.txt"))
run(`rm BIOGRID-ALL-$BIOGRID_VERSION.tab3.txt`)

# Filter out non-human
filter!(row -> row["Organism ID Interactor A"] == 9606 && row["Organism ID Interactor B"] == 9606, interactions_df)

# Only keep the interactions discovered with stringent methods (high-quality interactions)
CONSERVATIVE_SYSTEMS = Set([
      "Two-hybrid",
      "Affinity Capture-MS",
      "Affinity Capture-Western",
      "Reconstituted Complex",
      "Affinity Capture-Luminescence",
      "Co-crystal Structure",
      "Far Western",
      "FRET",
      "Protein-peptide",
      "Co-localization",
      "Affinity Capture-RNA",
      "Co-purification"
])

# Only keep interactions detected with one of the conservative systems
filter!(row -> row["Experimental System"] ∈ CONSERVATIVE_SYSTEMS, interactions_df)

# Get rid of useless columns
interactions = select(
    interactions_df[:, [
        "SWISS-PROT Accessions Interactor A", "SWISS-PROT Accessions Interactor B", "Publication Source"
    ]],
    "SWISS-PROT Accessions Interactor A" => "prot_a",
    "SWISS-PROT Accessions Interactor B" => "prot_b",
    "Publication Source" => "source"
)

# Count the number of times the interactions have been validated in different sources
sources = Dict()
pair_counts = Dict()

for pair in eachrow(interactions)
    p1 = min(pair["prot_a"], pair["prot_b"])
    p2 = max(pair["prot_a"], pair["prot_b"])
    source = pair["source"]

    pair_tuple = (p1, p2)

    if pair_tuple ∉ keys(pair_counts)
        pair_counts[pair_tuple] = 1
        sources[pair_tuple] = [source]
    elseif pair_tuple ∈ keys(pair_counts) && source ∉ sources[pair_tuple]
        pair_counts[pair_tuple] += 1
        push!(sources[pair_tuple], source)
    else
        continue
    end
end

# Only keep pairs that are present more than once (multiple lines of evidence)
duplicated_pairs = [k for (k, v) in pair_counts if v >= 2]

open(DESTINATION, "w") do io
    write(io, join(["$(pair[1])\t$(pair[2])" for pair in duplicated_pairs if pair[1] != "-" && pair[2] != "-"], "\n"))
end

max_stars_repo_name: GreenCUBIC/PeptideScreening

max_issues_repo_name: GreenCUBIC/PeptideScreening

max_forks_repo_name: GreenCUBIC/PeptideScreening


text:
module InterpretableModels
using CSV, BlackBoxOptim, DataFrames, Statistics, StatsBase
using Lasso, MLDataUtils, ROC

export SLIMScoringSystem, bboptimize_slim_scoring_system,
    isclassificationmodel, isregressionmodel, istimeseriesmodel, horizon

export LaggedFactorLITRM, MovingAverageLITRM, predict, fit!

export LassoLIMOptimizer, LIMScoringSystem

# Core and utils
include("utils.jl")
include("loss_functions.jl")
include(joinpath("models", "types.jl"))
include(joinpath("optimizers", "types.jl"))

# Optimizers
include(joinpath("optimizers", "blackboxoptim.jl"))

# Default optimizer is a black-box, heuristic one since they can fit all
# models. You can use a more model-specific optimizer to (maybe) get
# better performance.
DefaultOptimizer(m) = BlackBoxOptimizer(m)

# Models
include(joinpath("models", "slim_scoring.jl"))
include(joinpath("models", "abstract_linear_itrm.jl"))
include(joinpath("models", "lagged_factor_linear_itrm.jl"))
include(joinpath("models", "moving_average_linear_itrm.jl"))
include(joinpath("models", "lasso_simple_scoring.jl"))

# Optimizers

end # module


text:
import MeshArrays
MeshArrays.GRID_LLC90_download()

import OceanStateEstimation
OceanStateEstimation.ECCOdiags_download()
#OceanStateEstimation.ECCOdiags_add("release4")
#OceanStateEstimation.ECCOdiags_add("interp_coeffs")

import Downloads
Downloads.download(
  "https://zenodo.org/record/5784905/files/interp_coeffs_halfdeg.jld2",
  joinpath(OceanStateEstimation.ECCOdiags_path,"interp_coeffs_halfdeg.jld2");
  timeout=60000.0)

include("../tutorials/list_notebooks.jl")
notebooks=list_notebooks()
download_notebooks("jl",notebooks)
#find . -name "*.jl" |grep notebooks

max_issues_repo_name: gaelforget/MeshArrayNotebooks

max_forks_repo_name: JuliaClimate/Notebooks


text:
orders = 1:2:21#; orders = [3;7;9;11;21]
n_err_evals = 10_001 # x-vales to evaluate the error metrics with
n_sims = 5_000_000 # need lots of these for good error plots
n_boot = 1_000 # number of bootstrap samples to perform to estimate confidence intervals
qtiles = (0.05,0.95) # 90% confidence intervals
transient_time = 2.0 # time to evolve the solutions to
approx_types = (:dg,:dg,:order1,:qbdrap)


text:
let λ = Partition([5,4,2,2]), μ = Partition([2,1])
    @test λ \ μ == SkewDiagram(λ, μ)
    @test leglength(λ, μ) == leglength(λ \ μ) == 3
end
@test leglength(Partition(Int[]), Partition(Int[])) == -1

@test isrimhook(Partition([4,3,2]), Partition([2,2,2]))
@test !isrimhook(Partition([4,3,2]), Partition([2,2,1]))
@test !isrimhook(Partition([4,3,2]), Partition([1,1]))

let λ = Partition([5,4,2,1]), μ = Partition([4,3,2,2,1])
    @test partitionsequence(λ) == [1, 0, 1, 0, 1, 1, 0, 1, 0]
    @test character(λ, μ) == 0
end

let λ = μ = Partition([1])
    @test leglength(λ, μ) == -1
    @test character(λ, μ) == 1
end


text:
#= These methods are to be overloaded by Cases.jl =#
function update_surface end
function update_forcing end
function update_radiation end

function update_cloud_frac(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables)
    # update grid-mean cloud fraction and cloud cover
    aux_bulk = center_aux_bulk(state)
    aux_gm = center_aux_grid_mean(state)
    aux_en = center_aux_environment(state)
    a_up_bulk = aux_bulk.area
    @inbounds for k in real_center_indices(grid) # update grid-mean cloud fraction and cloud cover
        aux_gm.cloud_fraction[k] = aux_en.area[k] * aux_en.cloud_fraction[k] + a_up_bulk[k] * aux_bulk.cloud_fraction[k]
    end
    gm.cloud_cover = min(edmf.EnvVar.cloud_cover + sum(edmf.UpdVar.cloud_cover), 1)
end

function compute_les_Γᵣ(z::ClimaCore.Geometry.ZPoint, τᵣ::Real = 24.0 * 3600.0, zᵢ::Real = 3000.0, zᵣ::Real = 3500.0)
    # returns height-dependent relaxation timescale from eqn. 9 in `Shen et al. 2021`
    if z < zᵢ
        return 0.0
    elseif zᵢ <= z <= zᵣ
        cos_arg = pi * ((z - zᵢ) / (zᵣ - zᵢ))
        return (0.5 / τᵣ) * (1 - cos(cos_arg))
    elseif z > zᵣ
        return (1 / τᵣ)
    end
end

function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE{N_up}, grid, state, Case, gm, TS) where {N_up}
    tendencies_gm = center_tendencies_grid_mean(state)
    kc_toa = kc_top_of_atmos(grid)
    FT = eltype(grid)
    param_set = parameter_set(gm)
    prog_gm = center_prog_grid_mean(state)
    aux_gm = center_aux_grid_mean(state)
    ∇θ_liq_ice_gm = center_aux_grid_mean(state).∇θ_liq_ice_gm
    ∇q_tot_gm = center_aux_grid_mean(state).∇q_tot_gm
    aux_en = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)
    aux_up = center_aux_updrafts(state)
    aux_bulk = center_aux_bulk(state)
    aux_tc_f = face_aux_turbconv(state)
    aux_up_f = face_aux_updrafts(state)
    ρ0_f = face_ref_state(state).ρ0
    p0_c = center_ref_state(state).p0
    α0_c = center_ref_state(state).α0
    kf_surf = kf_surface(grid)
    kc_surf = kc_surface(grid)
    up = edmf.UpdVar
    en = edmf.EnvVar
    massflux_h = aux_tc_f.massflux_h
    massflux_qt = aux_tc_f.massflux_qt
    aux_tc = center_aux_turbconv(state)

    θ_liq_ice_gm_toa = prog_gm.θ_liq_ice[kc_toa]
    q_tot_gm_toa = prog_gm.q_tot[kc_toa]
    RBθ = CCO.RightBiasedC2F(; top = CCO.SetValue(θ_liq_ice_gm_toa))
    RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa))
    wvec = CC.Geometry.WVector
    ∇ = CCO.DivergenceF2C()
    @. ∇θ_liq_ice_gm = ∇(wvec(RBθ(prog_gm.θ_liq_ice)))
    @. ∇q_tot_gm = ∇(wvec(RBq(prog_gm.q_tot)))

    @inbounds for k in real_center_indices(grid)
        # Apply large-scale horizontal advection tendencies
        ts = thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k])
        Π = TD.exner(ts)

        if Case.Fo.apply_coriolis
            tendencies_gm.u[k] -= Case.Fo.coriolis_param * (aux_gm.vg[k] - prog_gm.v[k])
            tendencies_gm.v[k] += Case.Fo.coriolis_param * (aux_gm.ug[k] - prog_gm.u[k])
        end
        if rad_type(Case.Rad) <: Union{RadiationDYCOMS_RF01, RadiationLES}
            tendencies_gm.θ_liq_ice[k] += aux_gm.dTdt_rad[k] / Π
        end

        if force_type(Case.Fo) <: ForcingDYCOMS_RF01
            tendencies_gm.q_tot[k] += aux_gm.dqtdt[k]
            # Apply large-scale subsidence tendencies
            tendencies_gm.θ_liq_ice[k] -= ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
            tendencies_gm.q_tot[k] -= ∇q_tot_gm[k] * aux_gm.subsidence[k]
        end

        if force_type(Case.Fo) <: ForcingStandard
            if Case.Fo.apply_subsidence
                tendencies_gm.θ_liq_ice[k] -= ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
                tendencies_gm.q_tot[k] -= ∇q_tot_gm[k] * aux_gm.subsidence[k]
            end
            tendencies_gm.θ_liq_ice[k] += aux_gm.dTdt[k] / Π
            tendencies_gm.q_tot[k] += aux_gm.dqtdt[k]
        end

        if force_type(Case.Fo) <: ForcingLES
            H_horz_adv = aux_gm.dTdt_hadv[k] / Π
            H_fluc = aux_gm.dTdt_fluc[k] / Π

            gm_U_nudge_k = (aux_gm.u_nudge[k] - prog_gm.u[k]) / Case.Fo.nudge_tau
            gm_V_nudge_k = (aux_gm.v_nudge[k] - prog_gm.v[k]) / Case.Fo.nudge_tau

            Γᵣ = compute_les_Γᵣ(grid.zc[k])
            if Γᵣ != 0
                tau_k = 1.0 / Γᵣ
                gm_H_nudge_k = (aux_gm.H_nudge[k] - prog_gm.θ_liq_ice[k]) / tau_k
                gm_q_tot_nudge_k = (aux_gm.qt_nudge[k] - prog_gm.q_tot[k]) / tau_k
            else
                gm_H_nudge_k = 0.0
                gm_q_tot_nudge_k = 0.0
            end

            if Case.Fo.apply_subsidence
                # Apply large-scale subsidence tendencies
                gm_H_subsidence_k = -∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
                gm_QT_subsidence_k = -∇q_tot_gm[k] * aux_gm.subsidence[k]
            else
                gm_H_subsidence_k = 0.0
                gm_QT_subsidence_k = 0.0
            end

            tendencies_gm.θ_liq_ice[k] += H_horz_adv + gm_H_nudge_k + H_fluc + gm_H_subsidence_k
            tendencies_gm.q_tot[k] +=
                aux_gm.dqtdt_hadv[k] + gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k] + gm_QT_subsidence_k

            tendencies_gm.u[k] += gm_U_nudge_k
            tendencies_gm.v[k] += gm_V_nudge_k
        end
        tendencies_gm.q_tot[k] +=
            aux_bulk.qt_tendency_precip_formation[k] +
            aux_en.qt_tendency_precip_formation[k] +
            aux_tc.qt_tendency_precip_sinks[k]
        tendencies_gm.θ_liq_ice[k] +=
            aux_bulk.θ_liq_ice_tendency_precip_formation[k] +
            aux_en.θ_liq_ice_tendency_precip_formation[k] +
            aux_tc.θ_liq_ice_tendency_precip_sinks[k]
    end

    # TODO: we shouldn't need to call parent here
    a_en = aux_en.area
    w_en = aux_en_f.w
    a_en_bcs = (; bottom = CCO.SetValue(1 - sum(i -> edmf.area_surface_bc[i], 1:N_up)), top = CCO.Extrapolate())
    Ifae = CCO.InterpolateC2F(; a_en_bcs...)
    # Compute the mass flux and associated scalar fluxes
    @inbounds for i in 1:N_up
        a_up_bcs = (; bottom = CCO.SetValue(edmf.area_surface_bc[i]), top = CCO.Extrapolate())
        Ifau = CCO.InterpolateC2F(; a_up_bcs...)
        a_up = aux_up[i].area
        w_up = aux_up_f[i].w
        @. aux_up_f[i].massflux = ρ0_f * Ifau(a_up) * Ifae(a_en) * (w_up - w_en)
        aux_up_f[i].massflux[kf_surf] = 0
    end

    parent(massflux_h) .= 0
    parent(massflux_qt) .= 0
    If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
    @inbounds for i in 1:N_up
        # We know that, since W = 0 at z = 0, m = 0 also, and
        # therefore θ_liq_ice / q_tot values do not matter
        @. massflux_h += aux_up_f[i].massflux * (If(aux_up[i].θ_liq_ice) - If(aux_en.θ_liq_ice))
        @. massflux_qt += aux_up_f[i].massflux * (If(aux_up[i].q_tot) - If(aux_en.q_tot))
    end

    wvec = CC.Geometry.WVector
    ∇c = CCO.DivergenceF2C()
    # Compute the  mass flux tendencies
    # Adjust the values of the grid mean variables
    # Prepare the output
    @. aux_tc.massflux_tendency_h = -∇c(wvec(massflux_h)) * α0_c
    @. aux_tc.massflux_tendency_qt = ∇c(wvec(massflux_qt))
    @. tendencies_gm.θ_liq_ice += -∇c(wvec(massflux_h)) * α0_c
    @. tendencies_gm.q_tot += -∇c(wvec(massflux_qt)) * α0_c

    aeKHq_tot_bc = Case.Sur.rho_qtflux / aux_en.area[kc_surf]
    aeKHθ_liq_ice_bc = Case.Sur.rho_hflux / aux_en.area[kc_surf]
    aeKMu_bc = Case.Sur.rho_uflux / aux_en.area[kc_surf]
    aeKMv_bc = Case.Sur.rho_vflux / aux_en.area[kc_surf]

    tbc = (; top = CCO.SetValue(wvec(FT(0))))
    ∇aeKH_q_tot = CCO.DivergenceF2C(; bottom = CCO.SetValue(wvec(aeKHq_tot_bc)), tbc...)
    ∇aeKH_θ_liq_ice = CCO.DivergenceF2C(; bottom = CCO.SetValue(wvec(aeKHθ_liq_ice_bc)), tbc...)
    ∇aeKM_u = CCO.DivergenceF2C(; bottom = CCO.SetValue(wvec(aeKMu_bc)), tbc...)
    ∇aeKM_v = CCO.DivergenceF2C(; bottom = CCO.SetValue(wvec(aeKMv_bc)), tbc...)

    @. tendencies_gm.q_tot += -α0_c * aux_en.area * ∇aeKH_q_tot(wvec(aux_tc_f.diffusive_flux_qt))
    @. tendencies_gm.θ_liq_ice += -α0_c * aux_en.area * ∇aeKH_θ_liq_ice(wvec(aux_tc_f.diffusive_flux_h))
    @. tendencies_gm.u += -α0_c * aux_en.area * ∇aeKM_u(wvec(aux_tc_f.diffusive_flux_u))
    @. tendencies_gm.v += -α0_c * aux_en.area * ∇aeKM_v(wvec(aux_tc_f.diffusive_flux_v))

    return nothing
end

function compute_diffusive_fluxes(
    edmf::EDMF_PrognosticTKE,
    grid,
    state,
    gm::GridMeanVariables,
    Case::CasesBase,
    TS::TimeStepping,
    param_set,
)
    FT = eltype(grid)
    ρ0_f = face_ref_state(state).ρ0
    aux_bulk = center_aux_bulk(state)
    aux_tc_f = face_aux_turbconv(state)
    aux_en = center_aux_environment(state)
    KM = center_aux_turbconv(state).KM
    KH = center_aux_turbconv(state).KH
    aeKM = aux_en.area .* KM
    aeKH = aux_en.area .* KH
    kc_surf = kc_surface(grid)
    kc_toa = kc_top_of_atmos(grid)
    kf_surf = kf_surface(grid)
    prog_gm = center_prog_grid_mean(state)
    IfKH = CCO.InterpolateC2F(; bottom = CCO.SetValue(aeKH[kc_surf]), top = CCO.SetValue(aeKH[kc_toa]))
    IfKM = CCO.InterpolateC2F(; bottom = CCO.SetValue(aeKM[kc_surf]), top = CCO.SetValue(aeKM[kc_toa]))

    @. aux_tc_f.ρ_ae_KH = IfKH(aeKH) * ρ0_f
    @. aux_tc_f.ρ_ae_KM = IfKM(aeKM) * ρ0_f

    aeKHq_tot_bc = -Case.Sur.rho_qtflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf]
    aeKHθ_liq_ice_bc = -Case.Sur.rho_hflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf]
    aeKMu_bc = -Case.Sur.rho_uflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf]
    aeKMv_bc = -Case.Sur.rho_vflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf]

    ∇q_tot_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHq_tot_bc), top = CCO.SetDivergence(FT(0)))
    ∇θ_liq_ice_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHθ_liq_ice_bc), top = CCO.SetDivergence(FT(0)))
    ∇u_gm = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKMu_bc), top = CCO.SetDivergence(FT(0)))
    ∇v_gm = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKMv_bc), top = CCO.SetDivergence(FT(0)))

    wvec = CC.Geometry.WVector
    @. aux_tc_f.diffusive_flux_qt = -aux_tc_f.ρ_ae_KH * ∇q_tot_en(wvec(aux_en.q_tot))
    @. aux_tc_f.diffusive_flux_h = -aux_tc_f.ρ_ae_KH * ∇θ_liq_ice_en(wvec(aux_en.θ_liq_ice))
    @. aux_tc_f.diffusive_flux_u = -aux_tc_f.ρ_ae_KM * ∇u_gm(wvec(prog_gm.u))
    @. aux_tc_f.diffusive_flux_v = -aux_tc_f.ρ_ae_KM * ∇v_gm(wvec(prog_gm.v))

    return
end

function affect_filter!(edmf, grid, state, gm, Case, TS)
    # TODO: figure out why this filter kills the DryBubble results if called at t = 0.
    if Case.casename == "DryBubble" && !(TS.t > 0)
        return nothing
    end
    prog_en = center_prog_environment(state)
    up = edmf.UpdVar
    ###
    ### Filters
    ###
    set_edmf_surface_bc(edmf, grid, state, up, Case.Sur)
    filter_updraft_vars(edmf, grid, state, gm)

    @inbounds for k in real_center_indices(grid)
        prog_en.ρatke[k] = max(prog_en.ρatke[k], 0.0)
        prog_en.ρaHvar[k] = max(prog_en.ρaHvar[k], 0.0)
        prog_en.ρaQTvar[k] = max(prog_en.ρaQTvar[k], 0.0)
        prog_en.ρaHQTcov[k] = max(prog_en.ρaHQTcov[k], -sqrt(prog_en.ρaHvar[k] * prog_en.ρaQTvar[k]))
        prog_en.ρaHQTcov[k] = min(prog_en.ρaHQTcov[k], sqrt(prog_en.ρaHvar[k] * prog_en.ρaQTvar[k]))
    end
end

# Compute the sum of tendencies for the scheme
function ∑tendencies!(tendencies, prog, params, t)
    UnPack.@unpack edmf, grid, gm, case, aux, TS = params

    TS.t = t
    state = State(prog, aux, tendencies)

    affect_filter!(edmf, grid, state, gm, case, TS)

    gm = gm
    up = edmf.UpdVar
    en = edmf.EnvVar
    param_set = parameter_set(gm)
    en_thermo = edmf.EnvThermo

    # Update aux / pre-tendencies filters. TODO: combine these into a function that minimizes traversals
    # Some of these methods should probably live in `compute_tendencies`, when written, but we'll
    # treat them as auxiliary variables for now, until we disentangle the tendency computations.

    compute_updraft_surface_bc(edmf, grid, state, case)
    update_aux!(edmf, gm, grid, state, case, param_set, TS)

    tends_face = tendencies.face
    tends_cent = tendencies.cent
    parent(tends_face) .= 0
    parent(tends_cent) .= 0

    # causes division error in dry bubble first time step
    compute_precipitation_formation_tendencies(grid, state, up, edmf.Precip, TS.dt, param_set)

    microphysics(en_thermo, grid, state, en, edmf.Precip, TS.dt, param_set) # saturation adjustment + rain creation
    if edmf.Precip.precipitation_model == "clima_1m"
        compute_precipitation_sink_tendencies(grid, state, gm, TS)
        compute_precipitation_advection_tendencies(edmf, grid, state, gm)
    end

    # compute tendencies
    compute_gm_tendencies!(edmf, grid, state, case, gm, TS)
    compute_updraft_tendencies(edmf, grid, state, gm)

    compute_en_tendencies!(edmf, grid, state, param_set, TS, Val(:tke), Val(:ρatke))
    compute_en_tendencies!(edmf, grid, state, param_set, TS, Val(:Hvar), Val(:ρaHvar))
    compute_en_tendencies!(edmf, grid, state, param_set, TS, Val(:QTvar), Val(:ρaQTvar))
    compute_en_tendencies!(edmf, grid, state, param_set, TS, Val(:HQTcov), Val(:ρaHQTcov))

    return
end

function set_edmf_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, up, surface)
    en = edmf.EnvVar
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    ρ0_c = center_ref_state(state).ρ0
    ρ0_f = face_ref_state(state).ρ0
    prog_up = center_prog_updrafts(state)
    prog_en = center_prog_environment(state)
    prog_up_f = face_prog_updrafts(state)
    aux_bulk = center_aux_bulk(state)
    @inbounds for i in 1:(up.n_updrafts)
        prog_up[i].ρarea[kc_surf] = ρ0_c[kc_surf] * edmf.area_surface_bc[i]
        prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * edmf.h_surface_bc[i]
        prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * edmf.qt_surface_bc[i]
        prog_up_f[i].ρaw[kf_surf] = ρ0_f[kf_surf] * edmf.w_surface_bc[i]
    end

    flux1 = surface.rho_hflux
    flux2 = surface.rho_qtflux
    zLL = grid.zc[kc_surf]
    ustar = surface.ustar
    oblength = surface.obukhov_length
    α0LL = center_ref_state(state).α0[kc_surf]
    # TODO: is bulk even defined before this is called?
    ae = 1 .- aux_bulk.area # area of environment

    ρ0_ae = ρ0_c[kc_surf] * ae[kc_surf]

    prog_en.ρatke[kc_surf] = ρ0_ae * get_surface_tke(surface.ustar, grid.zc[kc_surf], surface.obukhov_length)
    prog_en.ρaHvar[kc_surf] = ρ0_ae * get_surface_variance(flux1 * α0LL, flux1 * α0LL, ustar, zLL, oblength)
    prog_en.ρaQTvar[kc_surf] = ρ0_ae * get_surface_variance(flux2 * α0LL, flux2 * α0LL, ustar, zLL, oblength)
    prog_en.ρaHQTcov[kc_surf] = ρ0_ae * get_surface_variance(flux1 * α0LL, flux2 * α0LL, ustar, zLL, oblength)
end

function compute_updraft_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, Case::CasesBase)
    kc_surf = kc_surface(grid)

    Δzi = grid.Δzi
    zLL = grid.zc[kc_surf]
    ustar = Case.Sur.ustar
    oblength = Case.Sur.obukhov_length
    α0LL = center_ref_state(state).α0[kc_surf]
    prog_gm = center_prog_grid_mean(state)
    qt_var = get_surface_variance(Case.Sur.rho_qtflux * α0LL, Case.Sur.rho_qtflux * α0LL, ustar, zLL, oblength)
    h_var = get_surface_variance(Case.Sur.rho_hflux * α0LL, Case.Sur.rho_hflux * α0LL, ustar, zLL, oblength)

    if Case.Sur.bflux > 0.0
        a_total = edmf.surface_area
        edmf.entr_surface_bc = 2.0 * Δzi
        edmf.detr_surface_bc = 0.0
        a_ = a_total / edmf.n_updrafts
        @inbounds for i in 1:(edmf.n_updrafts)
            surface_scalar_coeff =
                percentile_bounds_mean_norm(1.0 - a_total + i * a_, 1.0 - a_total + (i + 1) * a_, 1000)
            edmf.area_surface_bc[i] = a_
            edmf.w_surface_bc[i] = 0.0
            edmf.h_surface_bc[i] = prog_gm.θ_liq_ice[kc_surf] + surface_scalar_coeff * sqrt(h_var)
            edmf.qt_surface_bc[i] = prog_gm.q_tot[kc_surf] + surface_scalar_coeff * sqrt(qt_var)
        end
    else
        edmf.entr_surface_bc = 0.0
        edmf.detr_surface_bc = 2.0 * Δzi
        @inbounds for i in 1:(edmf.n_updrafts)
            edmf.area_surface_bc[i] = 0
            edmf.w_surface_bc[i] = 0.0
            edmf.h_surface_bc[i] = prog_gm.θ_liq_ice[kc_surf]
            edmf.qt_surface_bc[i] = prog_gm.q_tot[kc_surf]
        end
    end

    return
end

# Note: this assumes all variables are defined on half levels not full levels (i.e. phi, psi are not w)
# if covar_e.name is not "tke".
function get_GMV_CoVar(
    edmf::EDMF_PrognosticTKE{N_up},
    grid,
    state,
    covar_sym::Symbol,
    ϕ_sym::Symbol,
    ψ_sym::Symbol = ϕ_sym,
) where {N_up}

    is_tke = covar_sym == :tke
    tke_factor = is_tke ? 0.5 : 1
    prog_gm_c = center_prog_grid_mean(state)
    prog_gm_f = face_prog_grid_mean(state)
    aux_up_f = face_aux_updrafts(state)
    aux_en_c = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)
    aux_en = is_tke ? aux_en_f : aux_en_c
    aux_up_c = center_aux_updrafts(state)
    aux_up = is_tke ? aux_up_f : aux_up_c
    gmv_covar = getproperty(center_aux_grid_mean(state), covar_sym)
    covar_e = getproperty(center_aux_environment(state), covar_sym)
    prog_gm = is_tke ? prog_gm_f : prog_gm_c
    ϕ_gm = getproperty(prog_gm, ϕ_sym)
    ψ_gm = getproperty(prog_gm, ψ_sym)
    ϕ_en = getproperty(aux_en, ϕ_sym)
    ψ_en = getproperty(aux_en, ψ_sym)
    area_en = aux_en_c.area

    Icd = is_tke ? CCO.InterpolateF2C() : x -> x
    @. gmv_covar = tke_factor * area_en * Icd(ϕ_en - ϕ_gm) * Icd(ψ_en - ψ_gm) + area_en * covar_e
    @inbounds for i in 1:N_up
        ϕ_up = getproperty(aux_up[i], ϕ_sym)
        ψ_up = getproperty(aux_up[i], ψ_sym)
        @. gmv_covar += tke_factor * aux_up_c[i].area * Icd(ϕ_up - ϕ_gm) * Icd(ψ_up - ψ_gm)
    end
    return
end

function compute_pressure_plume_spacing(edmf::EDMF_PrognosticTKE, param_set)

    H_up_min = CPEDMF.H_up_min(param_set)
    @inbounds for i in 1:(edmf.n_updrafts)
        edmf.pressure_plume_spacing[i] =
            max(edmf.aspect_ratio * edmf.UpdVar.updraft_top[i], H_up_min * edmf.aspect_ratio)
    end
    return
end

function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE{N_up}, grid, state, gm::GridMeanVariables) where {N_up}
    param_set = parameter_set(gm)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    FT = eltype(grid)

    aux_up = center_aux_updrafts(state)
    aux_en = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)
    aux_up_f = face_aux_updrafts(state)
    tendencies_up = center_tendencies_updrafts(state)
    tendencies_up_f = face_tendencies_updrafts(state)
    ρ0_c = center_ref_state(state).ρ0
    ρ0_f = face_ref_state(state).ρ0
    au_lim = edmf.max_area

    @inbounds for i in 1:N_up
        aux_up[i].entr_turb_dyn .= aux_up[i].entr_sc .+ aux_up[i].frac_turb_entr
        aux_up[i].detr_turb_dyn .= aux_up[i].detr_sc .+ aux_up[i].frac_turb_entr
    end

    UB = CCO.UpwindBiasedProductC2F(bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
    Ic = CCO.InterpolateF2C()

    wvec = CC.Geometry.WVector
    ∂ = CCO.DivergenceF2C()
    w_bcs = (; bottom = CCO.SetValue(wvec(FT(0))), top = CCO.SetValue(wvec(FT(0))))
    LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(FT(0)))

    # Solve for updraft area fraction
    @inbounds for i in 1:N_up
        aux_up_i = aux_up[i]
        w_up = aux_up_f[i].w
        a_up = aux_up_i.area
        q_tot_up = aux_up_i.q_tot
        q_tot_en = aux_en.q_tot
        θ_liq_ice_en = aux_en.θ_liq_ice
        θ_liq_ice_up = aux_up_i.θ_liq_ice
        entr_turb_dyn = aux_up_i.entr_turb_dyn
        detr_turb_dyn = aux_up_i.detr_turb_dyn
        θ_liq_ice_tendency_precip_formation = aux_up_i.θ_liq_ice_tendency_precip_formation
        qt_tendency_precip_formation = aux_up_i.qt_tendency_precip_formation

        tends_ρarea = tendencies_up[i].ρarea
        tends_ρaθ_liq_ice = tendencies_up[i].ρaθ_liq_ice
        tends_ρaq_tot = tendencies_up[i].ρaq_tot

        @. tends_ρarea =
            -∂(wvec(LB(Ic(w_up) * ρ0_c * a_up))) + (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn) -
            (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn)

        @. tends_ρaθ_liq_ice =
            -∂(wvec(LB(Ic(w_up) * ρ0_c * a_up * θ_liq_ice_up))) +
            (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn * θ_liq_ice_en) -
            (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn * θ_liq_ice_up) + (ρ0_c * θ_liq_ice_tendency_precip_formation)

        @. tends_ρaq_tot =
            -∂(wvec(LB(Ic(w_up) * ρ0_c * a_up * q_tot_up))) + (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn * q_tot_en) -
            (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn * q_tot_up) + (ρ0_c * qt_tendency_precip_formation)

        tends_ρarea[kc_surf] = 0
        tends_ρaθ_liq_ice[kc_surf] = 0
        tends_ρaq_tot[kc_surf] = 0
    end

    # Solve for updraft velocity

    # We know that, since W = 0 at z = 0, BCs for entr, detr,
    # and buoyancy should not matter in the end
    zero_bcs = (; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
    I0f = CCO.InterpolateC2F(; zero_bcs...)
    adv_bcs = (; bottom = CCO.SetValue(wvec(FT(0))), top = CCO.SetValue(wvec(FT(0))))
    LB = CCO.LeftBiasedF2C(; bottom = CCO.SetValue(FT(0)))
    ∂ = CCO.DivergenceC2F(; adv_bcs...)

    @inbounds for i in 1:N_up
        a_up_bcs = (; bottom = CCO.SetValue(edmf.area_surface_bc[i]), top = CCO.SetValue(FT(0)))
        Iaf = CCO.InterpolateC2F(; a_up_bcs...)
        tends_ρaw = tendencies_up_f[i].ρaw
        nh_pressure = aux_up_f[i].nh_pressure
        a_up = aux_up[i].area
        w_up = aux_up_f[i].w
        w_en = aux_en_f.w
        entr_w = aux_up[i].entr_turb_dyn
        detr_w = aux_up[i].detr_turb_dyn
        buoy = aux_up[i].buoy

        @. tends_ρaw =
            -(∂(wvec(LB(Iaf(a_up) * ρ0_f * w_up * w_up)))) +
            (ρ0_f * Iaf(a_up) * w_up * (I0f(entr_w) * w_en - I0f(detr_w) * w_up)) +
            (ρ0_f * Iaf(a_up) * I0f(buoy)) +
            nh_pressure
        tends_ρaw[kf_surf] = 0
    end

    return
end

function filter_updraft_vars(edmf::EDMF_PrognosticTKE{N_up}, grid, state, gm::GridMeanVariables) where {N_up}
    param_set = parameter_set(gm)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    FT = eltype(grid)

    up = edmf.UpdVar
    prog_up = center_prog_updrafts(state)
    prog_gm = center_prog_grid_mean(state)
    aux_up = center_aux_updrafts(state)
    aux_up_f = face_aux_updrafts(state)
    prog_up_f = face_prog_updrafts(state)
    ρ0_c = center_ref_state(state).ρ0
    ρ0_f = face_ref_state(state).ρ0

    @inbounds for i in 1:N_up
        prog_up[i].ρarea .= max.(prog_up[i].ρarea, 0)
        prog_up[i].ρaθ_liq_ice .= max.(prog_up[i].ρaθ_liq_ice, 0)
        prog_up[i].ρaq_tot .= max.(prog_up[i].ρaq_tot, 0)
        @inbounds for k in real_center_indices(grid)
            prog_up[i].ρarea[k] = min(prog_up[i].ρarea[k], ρ0_c[k] * edmf.max_area)
        end
    end

    a_min = edmf.minimum_area
    @inbounds for i in 1:N_up
        @. prog_up_f[i].ρaw = max.(prog_up_f[i].ρaw, 0)
        a_up_bcs = (; bottom = CCO.SetValue(edmf.area_surface_bc[i]), top = CCO.Extrapolate())
        If = CCO.InterpolateC2F(; a_up_bcs...)
        @. prog_up_f[i].ρaw = Int(If(prog_up[i].ρarea) >= ρ0_f * a_min) * prog_up_f[i].ρaw
    end

    @inbounds for k in real_center_indices(grid)
        @inbounds for i in 1:N_up
            is_surface_center(grid, k) && continue
            prog_up[i].ρaq_tot[k] = max(prog_up[i].ρaq_tot[k], 0)
            # this is needed to make sure Rico is unchanged.
            # TODO : look into it further to see why
            # a similar filtering of ρaθ_liq_ice breaks the simulation
            if prog_up[i].ρarea[k] / ρ0_c[k] < edmf.minimum_area
                prog_up[i].ρaq_tot[k] = 0
            end
        end
    end

    Ic = CCO.InterpolateF2C()
    @inbounds for i in 1:N_up
        @. prog_up[i].ρarea = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρarea)
        @. prog_up[i].ρaθ_liq_ice = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρaθ_liq_ice)
        @. prog_up[i].ρaq_tot = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρaq_tot)
        prog_up[i].ρarea[kc_surf] = ρ0_c[kc_surf] * edmf.area_surface_bc[i]
        prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * edmf.h_surface_bc[i]
        prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * edmf.qt_surface_bc[i]
    end
    return
end

function compute_covariance_shear(
    edmf::EDMF_PrognosticTKE,
    grid,
    state,
    gm::GridMeanVariables,
    ::Val{covar_sym},
    ϕs::Val{ϕ_en_sym},
    ψs::Val{ψ_en_sym} = ϕs,
) where {covar_sym, ϕ_en_sym, ψ_en_sym}

    aux_tc = center_aux_turbconv(state)
    ρ0_c = center_ref_state(state).ρ0
    prog_gm = center_prog_grid_mean(state)
    is_tke = covar_sym == :tke
    tke_factor = is_tke ? 0.5 : 1
    k_eddy = is_tke ? aux_tc.KM : aux_tc.KH
    aux_en_2m = center_aux_environment_2m(state)
    aux_covar = getproperty(aux_en_2m, covar_sym)

    aux_en_c = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)
    aux_en = is_tke ? aux_en_f : aux_en_c
    wvec = CC.Geometry.WVector
    ϕ_en = getproperty(aux_en, ϕ_en_sym)
    ψ_en = getproperty(aux_en, ψ_en_sym)
    FT = eltype(grid)

    bcs = (; bottom = CCO.Extrapolate(), top = CCO.SetGradient(wvec(zero(FT))))
    If = CCO.InterpolateC2F(; bcs...)
    ∇c = CCO.DivergenceF2C()
    u = prog_gm.u
    v = prog_gm.v
    area_en = aux_en_c.area
    shear = aux_covar.shear

    if is_tke
        @. shear =
            tke_factor *
            2 *
            ρ0_c *
            area_en *
            k_eddy *
            (∇c(wvec(ϕ_en)) * ∇c(wvec(ψ_en)) + (∇c(wvec(If(u))))^2 + (∇c(wvec(If(v))))^2)
    else
        @. shear = tke_factor * 2 * ρ0_c * area_en * k_eddy * ∇c(wvec(If(ϕ_en))) * ∇c(wvec(If(ψ_en)))
    end
    return
end

function compute_covariance_interdomain_src(
    edmf::EDMF_PrognosticTKE{N_up},
    grid,
    state,
    ::Val{covar_sym},
    ϕs::Val{ϕ_sym},
    ψs::Val{ψ_sym} = ϕs,
) where {N_up, covar_sym, ϕ_sym, ψ_sym}

    is_tke = covar_sym == :tke
    tke_factor = is_tke ? 0.5 : 1
    aux_up = center_aux_updrafts(state)
    aux_up_f = face_aux_updrafts(state)
    aux_en_2m = center_aux_environment_2m(state)
    interdomain = getproperty(aux_en_2m, covar_sym).interdomain
    prog_up = is_tke ? aux_up_f : aux_up
    aux_en_c = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)
    aux_en = is_tke ? aux_en_f : aux_en_c
    ϕ_en = getproperty(aux_en, ϕ_sym)
    ψ_en = getproperty(aux_en, ψ_sym)
    Ic = is_tke ? CCO.InterpolateF2C() : x -> x

    parent(interdomain) .= 0
    @inbounds for i in 1:N_up
        ϕ_up = getproperty(prog_up[i], ϕ_sym)
        ψ_up = getproperty(prog_up[i], ψ_sym)
        a_up = aux_up[i].area
        @. interdomain += tke_factor * a_up * (1 - a_up) * (Ic(ϕ_up) - Ic(ϕ_en)) * (Ic(ψ_up) - Ic(ψ_en))
    end
    return
end

function compute_covariance_entr(
    edmf::EDMF_PrognosticTKE{N_up},
    grid,
    state,
    ::Val{covar_sym},
    ϕs::Val{ϕ_sym},
    ψs::Val{ψ_sym} = ϕs,
) where {covar_sym, ϕ_sym, ψ_sym, N_up}

    ρ0_c = center_ref_state(state).ρ0
    FT = eltype(grid)
    is_tke = covar_sym == :tke
    tke_factor = is_tke ? 0.5 : 1
    aux_up = center_aux_updrafts(state)
    aux_up_f = face_aux_updrafts(state)
    prog_gm_c = center_prog_grid_mean(state)
    prog_gm_f = face_prog_grid_mean(state)
    prog_gm = is_tke ? prog_gm_f : prog_gm_c
    prog_up = is_tke ? aux_up_f : aux_up
    ϕ_gm = getproperty(prog_gm, ϕ_sym)
    ψ_gm = getproperty(prog_gm, ψ_sym)
    aux_en_2m = center_aux_environment_2m(state)
    aux_covar = getproperty(aux_en_2m, covar_sym)
    aux_en_c = center_aux_environment(state)
    covar = getproperty(aux_en_c, covar_sym)
    aux_en_f = face_aux_environment(state)
    aux_en = is_tke ? aux_en_f : aux_en_c
    ϕ_en = getproperty(aux_en, ϕ_sym)
    ψ_en = getproperty(aux_en, ψ_sym)
    entr_gain = aux_covar.entr_gain
    detr_loss = aux_covar.detr_loss
    Ic = CCO.InterpolateF2C()
    Idc = is_tke ? Ic : x -> x
    # TODO: we shouldn't need `parent` call here:
    parent(entr_gain) .= 0
    parent(detr_loss) .= 0
    min_area = edmf.minimum_area

    @inbounds for i in 1:N_up
        aux_up_i = aux_up[i]
        frac_turb_entr = aux_up_i.frac_turb_entr
        eps_turb = frac_turb_entr
        detr_sc = aux_up_i.detr_sc
        entr_sc = aux_up_i.entr_sc
        w_up = aux_up_f[i].w
        prog_up_i = prog_up[i]
        ϕ_up = getproperty(prog_up_i, ϕ_sym)
        ψ_up = getproperty(prog_up_i, ψ_sym)

        a_up = aux_up_i.area

        @. entr_gain +=
            Int(a_up > min_area) *
            (tke_factor * ρ0_c * a_up * abs(Ic(w_up)) * detr_sc * (Idc(ϕ_up) - Idc(ϕ_en)) * (Idc(ψ_up) - Idc(ψ_en))) + (
                tke_factor *
                ρ0_c *
                a_up *
                abs(Ic(w_up)) *
                eps_turb *
                ((Idc(ϕ_en) - Idc(ϕ_gm)) * (Idc(ψ_up) - Idc(ψ_en)) + (Idc(ψ_en) - Idc(ψ_gm)) * (Idc(ϕ_up) - Idc(ϕ_en)))
            )

        @. detr_loss += Int(a_up > min_area) * tke_factor * ρ0_c * a_up * abs(Ic(w_up)) * (entr_sc + eps_turb) * covar

    end

    return
end

function compute_covariance_dissipation(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, param_set)
    en = edmf.EnvVar
    c_d = CPEDMF.c_d(param_set)
    aux_tc = center_aux_turbconv(state)
    ρ0_c = center_ref_state(state).ρ0
    prog_en = center_prog_environment(state)
    aux_en = center_aux_environment(state)
    aux_en_2m = center_aux_environment_2m(state)
    aux_covar = getproperty(aux_en_2m, covar_sym)
    covar = getproperty(aux_en, covar_sym)
    dissipation = aux_covar.dissipation
    area_en = aux_en.area
    tke_en = aux_en.tke
    mixing_length = aux_tc.mixing_length

    @. dissipation = ρ0_c * area_en * covar * max(tke_en, 0)^0.5 / max(mixing_length, 1.0e-3) * c_d
    return
end

function compute_en_tendencies!(
    edmf::EDMF_PrognosticTKE{N_up},
    grid::Grid,
    state,
    param_set,
    TS,
    ::Val{covar_sym},
    ::Val{prog_sym},
) where {N_up, covar_sym, prog_sym}
    kc_surf = kc_surface(grid)
    kc_toa = kc_top_of_atmos(grid)
    ρ0_c = center_ref_state(state).ρ0
    ρ0_f = face_ref_state(state).ρ0
    prog_en = center_prog_environment(state)
    aux_en_2m = center_aux_environment_2m(state)
    tendencies_en = center_tendencies_environment(state)
    prog_covar = getproperty(tendencies_en, prog_sym)
    aux_up_f = face_aux_updrafts(state)
    aux_en = center_aux_environment(state)
    covar = getproperty(aux_en, covar_sym)
    aux_covar = getproperty(aux_en_2m, covar_sym)
    aux_up = center_aux_updrafts(state)
    w_en_f = face_aux_environment(state).w
    c_d = CPEDMF.c_d(param_set)
    is_tke = covar_sym == :tke
    FT = eltype(grid)

    ρ_ae_K = face_aux_turbconv(state).ρ_ae_K
    KM = center_aux_turbconv(state).KM
    KH = center_aux_turbconv(state).KH
    aux_tc = center_aux_turbconv(state)
    aux_bulk = center_aux_bulk(state)
    D_env = aux_tc.ϕ_temporary
    ae = 1 .- aux_bulk.area
    aeK = is_tke ? ae .* KM : ae .* KH

    press = aux_covar.press
    buoy = aux_covar.buoy
    shear = aux_covar.shear
    entr_gain = aux_covar.entr_gain
    rain_src = aux_covar.rain_src

    wvec = CC.Geometry.WVector
    aeK_bcs = (; bottom = CCO.SetValue(aeK[kc_surf]), top = CCO.SetValue(aeK[kc_toa]))
    prog_bcs = (; bottom = CCO.SetGradient(wvec(FT(0))), top = CCO.SetGradient(wvec(FT(0))))

    If = CCO.InterpolateC2F(; aeK_bcs...)
    ∇f = CCO.GradientC2F(; prog_bcs...)
    ∇c = CCO.DivergenceF2C()

    mixing_length = aux_tc.mixing_length
    min_area = edmf.minimum_area
    pressure_plume_spacing = edmf.pressure_plume_spacing

    Ic = CCO.InterpolateF2C()
    area_en = aux_en.area
    tke_en = aux_en.tke

    parent(D_env) .= 0

    @inbounds for i in 1:N_up
        turb_entr = aux_up[i].frac_turb_entr
        entr_sc = aux_up[i].entr_sc
        w_up = aux_up_f[i].w
        a_up = aux_up[i].area
        # TODO: using `Int(bool) *` means that NaNs can propagate
        # into the solution. Could we somehow call `ifelse` instead?
        @. D_env += Int(a_up > min_area) * ρ0_c * a_up * Ic(w_up) * (entr_sc + turb_entr)
    end

    RB = CCO.RightBiasedC2F(; top = CCO.SetValue(FT(0)))
    @. prog_covar =
        press + buoy + shear + entr_gain + rain_src - D_env * covar -
        (ρ0_c * area_en * c_d * sqrt(max(tke_en, 0)) / max(mixing_length, 1)) * covar -
        ∇c(wvec(RB(ρ0_c * area_en * Ic(w_en_f) * covar))) + ∇c(ρ0_f * If(aeK) * ∇f(covar))

    prog_covar[kc_surf] = covar[kc_surf]

    return nothing
end

function GMV_third_m(
    edmf::EDMF_PrognosticTKE{N_up},
    grid,
    state,
    ::Val{covar_en_sym},
    ::Val{var},
    ::Val{gm_third_m_sym},
) where {N_up, covar_en_sym, var, gm_third_m_sym}

    gm_third_m = getproperty(center_aux_grid_mean(state), gm_third_m_sym)
    kc_surf = kc_surface(grid)

    aux_bulk = center_aux_bulk(state)
    aux_up_f = face_aux_updrafts(state)
    is_tke = covar_en_sym == :tke
    aux_en_c = center_aux_environment(state)
    covar_en = getproperty(aux_en_c, covar_en_sym)
    aux_en_f = face_aux_environment(state)
    aux_en = is_tke ? aux_en_f : aux_en_c
    aux_up_c = center_aux_updrafts(state)
    aux_tc = center_aux_turbconv(state)
    ϕ_gm = aux_tc.ϕ_gm
    ϕ_gm_cov = aux_tc.ϕ_gm_cov
    ϕ_en_cov = aux_tc.ϕ_en_cov
    ϕ_up_cubed = aux_tc.ϕ_up_cubed
    aux_up = is_tke ? aux_up_f : aux_up_c
    var_en = getproperty(aux_en, var)
    area_en = aux_en_c.area
    Ic = is_tke ? CCO.InterpolateF2C() : x -> x
    wvec = CC.Geometry.WVector
    ∇c = CCO.DivergenceF2C()
    w_en = aux_en_f.w

    @. ϕ_gm = area_en * Ic(var_en)
    @inbounds for i in 1:N_up
        a_up = aux_up_c[i].area
        var_up = getproperty(aux_up[i], var)
        @. ϕ_gm += a_up * Ic(var_up)
    end

    if is_tke
        parent(ϕ_en_cov) .= 0
        @inbounds for i in 1:N_up
            horiz_K_eddy = aux_up_c[i].horiz_K_eddy
            a_up = aux_up_c[i].area
            a_bulk = aux_bulk.area
            @. ϕ_en_cov += -horiz_K_eddy * ∇c(wvec(w_en)) * a_up / a_bulk
        end
    else
        @. ϕ_en_cov = covar_en
    end

    parent(ϕ_up_cubed) .= 0
    @. ϕ_gm_cov = area_en * (ϕ_en_cov + (Ic(var_en) - ϕ_gm)^2)
    @inbounds for i in 1:N_up
        a_up = aux_up_c[i].area
        var_up = getproperty(aux_up[i], var)
        @. ϕ_gm_cov += a_up * (Ic(var_up) - ϕ_gm)^2
        @. ϕ_up_cubed += a_up * Ic(var_up)^3
    end

    @. gm_third_m = ϕ_up_cubed + area_en * (Ic(var_en)^3 + 3 * Ic(var_en) * ϕ_en_cov) - ϕ_gm^3 - 3 * ϕ_gm_cov * ϕ_gm

    gm_third_m[kc_surf] = 0
    return
end


text:
module PairPlots

using RecipesBase
using Measures
using NamedTupleTools
using Tables

using Statistics
using StatsBase

using Printf
using Requires
using StaticArrays
using PolygonOps
using Contour

function corner end
# We pull in Plots only using Requires instead of as a diret dependency, so
# that upstream packages can add us as a dependency e.g. for visualizing
# a custom datastructure without needing all of Plots.
# Give the user a helpful error message if they have not imported Plots yet
corner(args...; kwargs...) = error("You must run `using Plots` before using this package.")
export corner

function hist(a, histfunc, hist_kwargs, plotpercentiles, percentiles_kwargs, titlefmt, unit)

    x, h = histfunc(vec(a), hist_kwargs.nbins)

    if length(x) == 1
        @warn "1D histgoram has only one bin"
        # RecipesBase.plot!(; framestyle=:none, title="invalid", hist_kwargs.subplot, hist_kwargs.inset, )
        # return
    end

    minx, maxx = extrema(a)
    extent = maxx - minx
    # h_scaled = h ./ maximum(h) .* extent
    # h_scaled .+= minx
    h_scaled = h
    pcs = []
    a_sorted = sort(vec(a))
    if length(plotpercentiles) > 0
        pcs = quantile(a_sorted, plotpercentiles./100, sorted=true)
    end
    # WIP: replace the title with an annotation so that the margins don't get messed up.
    if get(hist_kwargs, :title, "") != ""
        pcs_title = quantile(a_sorted, [0.16, 0.5, 0.84], sorted=true)
        med = pcs_title[2]
        low = med - pcs_title[1]
        high =  pcs_title[3] - med
        title = @eval (@sprintf($titlefmt, $(hist_kwargs.title), $(pcs_title[2]), $high, $low, $unit))
        xc = mean(x)
        yc = maximum(h_scaled) * 1.15
        hist_kwargs = (;hist_kwargs..., title)
    end
    hist_kwargs = delete(hist_kwargs, :nbins)
    percentiles_kwargs = delete(percentiles_kwargs,:title)
    # p = RecipesBase.plot()

    kw = merge((;seriestype=:step), hist_kwargs)
    if kw.seriestype==:step
        x = x .- step(x)/2
    end
    RecipesBase.plot!(x, h_scaled; kw...)
    # RecipesBase.plot!(x, h_scaled; kw.inset, kw.subplot)
    # RecipesBase.plot!(x, h_scaled; kw.inset, kw.subplot, kw.seriestype, kw.xrotation, kw.yrotation, kw.tick_direction, kw.grid, kw.yformatter, kw.xlims, yticks)

    if length(pcs) > 0
        RecipesBase.plot!(pcs; seriestype=:vline, percentiles_kwargs...)
    end

    # return p
end
function hist(a, b, histfunc, hist2d_kwargs, contour_kwargs, scatter_kwargs, plotcontours, plotscatter, filterscatter)

    x, y, H = histfunc(vec(a), vec(b), hist2d_kwargs.nbins)

    threeD = get(hist2d_kwargs, :seriestype, nothing) == :wireframe

    # Calculate levels for contours
    levels = 1 .- exp.(-0.5 .* (0.5:0.5:2.1).^2)
    ii = sortperm(reshape(H,:))
    h2flat = H[ii]
    sm = cumsum(h2flat)
    sm /= sm[end]
    if all(isnan, sm) || length(H) <= 1
        @warn "Could not compute valid contours"
        plotcontours = false
        V = [0]
    else
        V = sort(map(v0 -> h2flat[sm .≤ v0][end], levels))
        if any(==(0), diff(V))
            @warn "Too few points to create valid contours"
        end
    end

    # # Exclude histogram spaces outside the outer contour
    # masked_weights = fill(NaN, size(H))
    # # # Old method of masking out histogram squares
    # # # if plotscatter && !threeD
    # # #     mask = H .>= V[1]
    # # # else
    # # #     mask = trues(size(H))
    # # # end
    # if plotscatter &&!threeD
    #     mask = H .> V[1]
    # else
    #     mask = trues(size(H))
    # end
    # masked_weights[mask] .= H[mask]
    masked_weights = H

    levels_final = [0; V; maximum(H) * (1 + 1e-4)]

    # p = RecipesBase.plot()
    if threeD
        scatter_kwargs = merge(scatter_kwargs, (; seriestype=:scatter3d))
        # if !hasproperty(contour_kwargs, :seriestype)
        #     contour_kwargs = merge(contour_kwargs, (;seriestype=:contour3d))
        # end
        plotcontours = false
    end
    RecipesBase.plot!(x, y, masked_weights; seriestype=:heatmap, hist2d_kwargs...)

    # We place a line of zeros on all sides of the histogram
    # for the sake of cleaner contours.
    pad_x = [first(x)-step(x); x; last(x)+step(x)]
    pad_y = [first(y)-step(y); y; last(y)+step(y)]
    pad_H = [
        zeros(size(H,2))' 0 0
        zeros(size(H,1)) H zeros(size(H,1))
        zeros(size(H,2))' 0 0
    ]
    if plotcontours
        # RecipesBase.plot!(x, y, H; seriestype=:contour, levels=levels_final, colorbar=:none, contour_kwargs...)

        # calculate the outer contour manually
        c = contours(pad_x,pad_y,pad_H', V)
        points = scatter_filtering(b, a, x,y, first(Contour.levels(c)))
    elseif filterscatter
        c = contours(pad_x,pad_y,pad_H', [V[1]])
    end

    if filterscatter
        points = scatter_filtering(b, a, x,y, first(Contour.levels(c)))
    else
        points = (b,a)
    end
    if plotscatter
        # TODO: handle threeD
        if get(scatter_kwargs, :seriestype, nothing) == :scatter3d
            z = zeros(size(y))
            RecipesBase.plot!(b, a, z; scatter_kwargs...)
        else
            RecipesBase.plot!(points; seriestype=:scatter, scatter_kwargs...)
        end
    end

    if plotcontours
        # RecipesBase.plot!(x, y, H; seriestype=:contour, levels=levels_final, colorbar=:none, contour_kwargs...)
        for level in Contour.levels(c), poly in lines(level)
            xs, ys = coordinates(poly)
            RecipesBase.plot!(xs, ys; levels=levels_final, colorbar=:none, contour_kwargs...)
        end

    end
    # p
end

emptyplot() = RecipesBase.plot(framestyle=:none, background_color_inside=:transparent)

# Default histogram calculations
"""
    prepare_hist(vector, nbins)

Use the StatsBase function to return a centered histogram from `vector` with `nbins`.
Must return a tuple of bin centres, followed by bin weights (of the same length).
"""
function prepare_hist(a, nbins)
    h = fit(Histogram, vec(a); nbins)
    # x = first(only(h.edges))+step(only(h.edges))/2 : step(only(h.edges)) : last(only(h.edges))-step(only(h.edges))/2
    x = range(first(h.edges[1])+step(h.edges[1])/2, step=step(h.edges[1]), length=size(h.weights,1))
    return x, h.weights
end
"""
    prepare_hist(vector, nbins)

Use the StatsBase function to return a centered histogram from `vector` with `nbins`x`nbins`.
Must return a tuple of bin centres along the horizontal, bin centres along the vertical,
and a matrix of bin weights (of matching dimensions).
"""
function prepare_hist(a, b, nbins)
    h = fit(Histogram, (vec(a),vec(b)); nbins)
    y = range(first(h.edges[1])+step(h.edges[1])/2, step=step(h.edges[1]), length=size(h.weights,1))
    x = range(first(h.edges[2])+step(h.edges[2])/2, step=step(h.edges[2]), length=size(h.weights,2))
    return x, y, h.weights
end

# Filter the scatter plot to remove points inside the first contour
# for performance and better display
function scatter_filtering(b, a, x,y, level)

    inds = eachindex(b,a)
    outside = trues(size(inds))

    # calculate the outer contour manually
    for poly in lines(level)
        xs, ys = coordinates(poly)
        poly = SVector.(xs,ys)
        push!(poly, SVector(xs[begin], ys[begin]))
        # poly = [xs ys; xs[begin] ys[begin]]
        for i in inds
            point = SVector(b[i],a[i])
            # This logic does not seem to be 100% right
            ins = inpolygon(point, poly, in=false, on=true, out=true)
            outside[i] &= ins
        end
    end
    return b[outside], a[outside]

end

# Special support for MCMCChains.
# They have additional table fields :iteration and :chain that should not be shown
# by default.
# In future, we might add the ability to colour the output by chain number.
function __init__()
    @require MCMCChains="c7f686f2-ff18-58e9-bc7b-31028e88f75d" begin
        function corner(chains::MCMCChains.Chains, args...; kwargs...)
            props = keys(chains)
            # values = (reshape(chains[prop].data,:) for prop in props)
            values = [vec(chains[prop].data) for prop in props]
            table = namedtuple(props, values)
            return corner(table, args...; kwargs...)
        end
    end

    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("corner.jl")

    nothing
end

end # module

max_stars_repo_name: sefffal/CornerPlots

max_issues_repo_name: sefffal/CornerPlots


text:
# usage example for package developer:
#
#     import BinarizationAPI: AbstractImageBinarizationAlgorithm,
#                             binarize, binarize!

"""
    AbstractImageBinarizationAlgorithm <: AbstractImageFilter

The root type for `ImageBinarization` package.

Any concrete binarization algorithm shall subtype it to support
[`binarize`](@ref) and [`binarize!`](@ref) APIs.

# Examples

All algorithms in ImageBinarization are called in the following pattern:

```julia
# first generate an algorithm instance
f = Otsu()

# then pass the algorithm to `binarize`
img₀₁ = binarize(img, f) # `eltype(img₀₁)` is `Gray{N0f8}`

# or use in-place version `binarize!`
img₀₁ = similar(img)
binarize!(img₀₁, img, f)
```

Some algorithms also receive additional information of image as an argument to infer the
"best" parameters, e.g., `window_size` of `AdaptiveThreshold`.

```julia
# you could explicit specify the it
f = AdaptiveThreshold(window_size = 32)

# or infer the "best" default value from given image
img = testimage("cameraman")
f = AdaptiveThreshold(img)
```

For more examples, please check [`binarize`](@ref) and [`binarize!`](@ref) and concret
algorithms.
"""
abstract type AbstractImageBinarizationAlgorithm <: AbstractImageFilter end

binarize!(out::GenericGrayImage,
          img,
          f::AbstractImageBinarizationAlgorithm,
          args...; kwargs...) =
    f(out, img, args...; kwargs...)

# in-place binarize! only supports gray image input
function binarize!(img::GenericGrayImage,
                   f::AbstractImageBinarizationAlgorithm,
                   args...; kwargs...)
    tmp = copy(img)
    f(img, tmp, args...; kwargs...)
    return img
end

function binarize(::Type{T},
                  img,
                  f::AbstractImageBinarizationAlgorithm,
                  args...; kwargs...) where T
    out = Array{T}(undef, size(img))
    binarize!(out, img, f, args...; kwargs...)
    return out
end

binarize(img::AbstractArray{T},
         f::AbstractImageBinarizationAlgorithm,
         args...; kwargs...) where T <:Union{Number, Colorant} =
    binarize(Gray{eltype(T)}, img, f, args...; kwargs...)

### Docstrings

"""
    binarize!([out,] img, f::AbstractImageBinarizationAlgorithm, args...; kwargs...)

Binarize `img` using algorithm `f`.

# Output

If `out` is specified, it will be changed in place. Otherwise `img` will be changed in place.

# Examples

Just simply pass an algorithm to `binarize!`:

```julia
img₀₁ = similar(img)
binarize!(img₀₁, img, f)
```

For cases you just want to change `img` in place, you don't necessarily need to manually
allocate `img₀₁`; just use the convenient method:

```julia
binarize!(img, f)
```

See also: [`binarize`](@ref)
"""
binarize!

"""
    binarize([T::Type,] img, f::AbstractImageBinarizationAlgorithm, args...; kwargs...)

Binarize `img` using algorithm `f`.

# Output

The return image `img₀₁` is an `Array{T}`.

If `T` is not specified, then it's inferred as
`Gray{eltype(eltype(img))}`, which is `Gray{N0f8}` for img of type `Array{N0f8}` and `Array{Gray{N0f8}}`, and `Gray{Float32}` for img of type `Array{Float32}` and `Array{Gray{Float32}}`

# Examples

Just simply pass the input image and algorithm to `binarize`

```julia
img₀₁ = binarize(img, f)
```

This reads as "`binarize` image `img` using binarization algorithm `f`".

You can also explicitly specify the return type:

```julia
img₀₁_float32 = binarize(Gray{Float32}, img, f)
```

See also [`binarize!`](@ref) for in-place binarization.
"""
binarize


text:
function resolve_import_block(x::EXPR, state::State, root, usinged, markfinal=true)
    if x.head == :as
        resolve_import_block(x.args[1], state, root, usinged, markfinal)
        if x.args[2].meta === nothing
            x.args[2].meta = Meta()
        end
        if hasbinding(last(x.args[1].args)) && CSTParser.isidentifier(x.args[2])
            lhsbinding = bindingof(last(x.args[1].args))
            x.args[2].meta.binding = Binding(x.args[2], lhsbinding.val, lhsbinding.type, lhsbinding.refs)
            setref!(x.args[2], bindingof(x.args[2]))
            last(x.args[1].args).meta.binding = nothing
        end
        return
    end
    n = length(x.args)
    for i = 1:length(x.args)
        arg = x.args[i]
        if isoperator(arg) && valof(arg) == "."
            # Leading dots. Can only be leading elements.
            if root == getsymbols(state)
                root = state.scope
            elseif root isa Scope && parentof(root) !== nothing
                root = parentof(root)
            else
                return
            end
        elseif isidentifier(arg) || (i == n && (CSTParser.ismacroname(arg) || isoperator(arg)))
            root = maybe_lookup(hasref(arg) ? refof(arg) : _get_field(root, arg, state), state)
            setref!(arg, root)
            if i == n
                markfinal && _mark_import_arg(arg, root, state, usinged)
                return refof(arg)
            end
        else
            return
        end
    end
end

function resolve_import(x::EXPR, state::State, root=getsymbols(state))
    if headof(x) === :using || headof(x) === :import
        usinged = headof(x) === :using
        if length(x.args) > 0 && isoperator(headof(x.args[1])) && valof(headof(x.args[1])) == ":"
            root = resolve_import_block(x.args[1].args[1], state, root, false, false)
            for i = 2:length(x.args[1].args)
                resolve_import_block(x.args[1].args[i], state, root, usinged)
            end
        else
            for i = 1:length(x.args)
                resolve_import_block(x.args[i], state, root, usinged)
            end
        end
    end
end

function _mark_import_arg(arg, par, state, usinged)
    if par !== nothing && CSTParser.is_id_or_macroname(arg)
        if par isa Binding # mark reference to binding
            push!(par.refs, arg)
        end
        if par isa SymbolServer.VarRef
            par = SymbolServer._lookup(par, getsymbols(state), true)
            !(par isa SymbolServer.SymStore) && return
        end
        if bindingof(arg) === nothing
            if !hasmeta(arg)
                arg.meta = Meta()
            end
            arg.meta.binding = Binding(arg, par, _typeof(par, state), [])
            setref!(arg, bindingof(arg))
        end

        if usinged
            if par isa SymbolServer.ModuleStore
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), par)
            elseif par isa Binding && par.val isa SymbolServer.ModuleStore
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), par.val)
            elseif par isa Binding && par.val isa EXPR && CSTParser.defines_module(par.val)
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), scopeof(par.val))
            elseif par isa Binding && par.val isa Binding && par.val.val isa EXPR && CSTParser.defines_module(par.val.val)
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), scopeof(par.val.val))
            end
        end
    end
end

function has_workspace_package(server, name)
    haskey(server.workspacepackages, name) &&
    hasscope(getcst(server.workspacepackages[name])) &&
    haskey(scopeof(getcst(server.workspacepackages[name])).names, name) &&
    scopeof(getcst(server.workspacepackages[name])).names[name] isa Binding &&
    scopeof(getcst(server.workspacepackages[name])).names[name].val isa EXPR &&
    CSTParser.defines_module(scopeof(getcst(server.workspacepackages[name])).names[name].val)
end

function add_to_imported_modules(scope::Scope, name::Symbol, val)
    if scope.modules isa Dict
        scope.modules[name] = val
    else
        Dict(name => val)
    end
end
no_modules_above(s::Scope) = !CSTParser.defines_module(s.expr) || s.parent === nothing || no_modules_above(s.parent)
function get_named_toplevel_module(s, name)
    return nothing
end
function get_named_toplevel_module(s::Scope, name::String)
    if CSTParser.defines_module(s.expr)
        m_name = CSTParser.get_name(s.expr)
        if ((headof(m_name) === :IDENTIFIER && valof(m_name) == name) || headof(m_name) === :NONSTDIDENTIFIER && length(m_name.args) == 2 && valof(m_name.args[2]) == name) && no_modules_above(s)
            return s.expr
        end
    end
    if s.parent isa Scope
        return get_named_toplevel_module(s.parent, name)
    end
    return nothing
end
function _get_field(par, arg, state)
    arg_str_rep = CSTParser.str_value(arg)
    if par isa SymbolServer.EnvStore
        if (arg_scope = retrieve_scope(arg)) !== nothing && (tlm = get_named_toplevel_module(arg_scope, arg_str_rep)) !== nothing && hasbinding(tlm)
            return bindingof(tlm)
        # elseif has_workspace_package(state.server, arg_str_rep)
        #     return scopeof(getcst(state.server.workspacepackages[arg_str_rep])).names[arg_str_rep]
        elseif haskey(par, Symbol(arg_str_rep))
            if isempty(state.env.project_deps) || Symbol(arg_str_rep) in state.env.project_deps
                return par[Symbol(arg_str_rep)]
            end
        end
    elseif par isa SymbolServer.ModuleStore # imported module
        if Symbol(arg_str_rep) === par.name.name
            return par
        elseif haskey(par, Symbol(arg_str_rep))
            par = par[Symbol(arg_str_rep)]
            if par isa SymbolServer.VarRef # reference to dependency
                return SymbolServer._lookup(par, getsymbols(state), true)
            end
            return par
        end
        for used_module_name in par.used_modules
            used_module = maybe_lookup(par[used_module_name], state)
            if used_module !== nothing && isexportedby(Symbol(arg_str_rep), used_module)
                return used_module[Symbol(arg_str_rep)]
            end
        end
    elseif par isa Scope
        if scopehasbinding(par, arg_str_rep)
            return par.names[arg_str_rep]
        elseif par.modules !== nothing
            for used_module in values(par.modules)
                if used_module isa SymbolServer.ModuleStore && isexportedby(Symbol(arg_str_rep), used_module)
                    return maybe_lookup(used_module[Symbol(arg_str_rep)], state)
                elseif used_module isa Scope && scope_exports(used_module, arg_str_rep, state)
                    return used_module.names[arg_str_rep]
                end
            end
        end
    elseif par isa Binding
        if par.val isa Binding
            return _get_field(par.val, arg, state)
        elseif par.val isa EXPR && CSTParser.defines_module(par.val) && scopeof(par.val) isa Scope
            return _get_field(scopeof(par.val), arg, state)
        elseif par.val isa EXPR && isassignment(par.val)
            if hasref(par.val.args[2])
                return _get_field(refof(par.val.args[2]), arg, state)
            elseif is_getfield_w_quotenode(par.val.args[2])
                return _get_field(refof_maybe_getfield(par.val.args[2]), arg, state)
            end
        elseif par.val isa SymbolServer.ModuleStore
            return _get_field(par.val, arg, state)
        end
    end
    return
end


text:
# Theories can just be vectors of rules!

macro theory(e)
    e = macroexpand(__module__, e)
    e = rmlines(e)
    if isexpr(e, :block)
        Vector{Rule}(e.args .|> x -> Rule(x; mod=__module__))
    else
        error("theory is not in form begin a => b; ... end")
    end
end

"""
A Theory is either a vector of [`Rule`](@ref) or
a compiled, callable function.
"""
const Theory = Union{Vector{Rule}, Function}