text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

#   hadamard_sum!(du, A, volume_flux, orientation, u, equations,
#                              sparsity_pattern, skip_index=(i,j)->false)
#
# Computes the flux difference ∑_j A[i, j] * f(u_i, u_j) and accumulates the result into `du`.
# - `du`, `u` are vectors
# - `A` is the skew-symmetric flux differencing matrix.
# - `sparsity_pattern`: either `nothing` or an AbstractSparseMatrix which specifies the sparsity
#   pattern of `A`
@inline function hadamard_sum!(du, A, volume_flux, orientation, u, equations,
                               sparsity_pattern::Nothing, skip_index=(i,j)->false)

  rows, cols = axes(A)
  for i in rows
    u_i = u[i]
    du_i = du[i]
    for j in cols
      # This routine computes only the upper-triangular part of the hadamard sum (A .* F).
      # We avoid computing the lower-triangular part, and instead accumulate those contributions
      # while computing the upper-triangular part (using the fact that A is skew-symmetric and F
      # is symmetric).
      if j > i && !skip_index(i, j)
          AF_ij = A[i,j] * volume_flux(u_i, u[j], orientation, equations)
          du_i = du_i + AF_ij
          du[j] = du[j] - AF_ij
      end
    end
    du[i] = du_i
  end
end

# If `skip_index` isn't specified, set it to `nothing` and dispatch based on `sparsity_pattern`
@inline function hadamard_sum!(du, A, volume_flux, orientation, u, equations, sparsity_pattern)
  return hadamard_sum!(du, A, volume_flux, orientation, u, equations, sparsity_pattern, nothing)
end

@inline function hadamard_sum!(du, A, volume_flux, orientation, u,
                               equations, sparsity_pattern::AbstractSparseMatrix{Bool},
                               skip_index::Union{Function, Nothing})
  n = size(sparsity_pattern, 2)
  rows = rowvals(sparsity_pattern)
  for i in 1:n
    u_i = u[i]
    du_i = du[i]
    for id in nzrange(sparsity_pattern, i)
      j = rows[id]
      # This routine computes only the upper-triangular part of the hadamard sum (A .* F).
      # We avoid computing the lower-triangular part, and instead accumulate those contributions
      # while computing the upper-triangular part (using the fact that A is skew-symmetric and F
      # is symmetric).
      if j > i
        AF_ij = A[i,j] * volume_flux(u_i, u[j], orientation, equations)
        du_i = du_i + AF_ij
        du[j] = du[j] - AF_ij
      end
    end
    du[i] = du_i
  end
end

# For DGMulti implementations, we construct "physical" differentiation operators by taking linear
# combinations of reference differentiation operators scaled by geometric change of variables terms.
# We use a lazy evaluation of physical differentiation operators, so that we can compute linear
# combinations of differentiation operators on-the-fly in an allocation-free manner.
function build_lazy_physical_derivative(element, orientation,
                                        mesh::VertexMappedMesh{2}, dg, cache)
  @unpack Qrst_skew = cache
  @unpack rxJ, sxJ, ryJ, syJ = mesh.md
  if orientation == 1
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (rxJ[1,element], sxJ[1,element]))
  else # if orientation == 2
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (ryJ[1,element], syJ[1,element]))
  end
end

function build_lazy_physical_derivative(element, orientation,
                                        mesh::VertexMappedMesh{3}, dg, cache)
  @unpack Qrst_skew = cache
  @unpack rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ = mesh.md
  if orientation == 1
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (rxJ[1,element], sxJ[1,element], txJ[1,element]))
  elseif orientation == 2
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (ryJ[1,element], syJ[1,element], tyJ[1,element]))
  else # if orientation == 3
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (rzJ[1,element], szJ[1,element], tzJ[1,element]))
  end
end

# Return the contravariant basis vector corresponding to the Cartesian
# coordinate diretion `orientation` in a given `element` of the `mesh`.
# The contravariant basis vectors have entries `dx_i / dxhat_j` where
# j ∈ {1, ..., NDIMS}. Here, `x_i` and `xhat_j` are the ith physical coordinate
# and jth reference coordinate, respectively. These are geometric terms which
# appear when using the chain rule to compute physical derivatives as a linear
# combination of reference derivatives.
@inline function get_contravariant_vector(element, orientation, mesh::VertexMappedMesh{2})
  @unpack rxJ, sxJ, ryJ, syJ = mesh.md
  if orientation == 1
    return 2 * SVector(rxJ[1, element], ryJ[1, element])
  else # if orientation == 2
    return 2 * SVector(sxJ[1, element], syJ[1, element])
  end
end

@inline function get_contravariant_vector(element, orientation, mesh::VertexMappedMesh{3})
  @unpack rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ = mesh.md
  if orientation == 1
    return 2 * SVector(rxJ[1, element], ryJ[1, element], rzJ[1, element])
  elseif orientation == 2
    return 2 * SVector(sxJ[1, element], syJ[1, element], szJ[1, element])
  else # if orientation == 3
    return 2 * SVector(txJ[1, element], tyJ[1, element], tzJ[1, element])
  end
end

function compute_flux_differencing_SBP_matrices(dg::DGMulti)
  rd = dg.basis
  Qrst_hybridized, VhP, Ph = StartUpDG.hybridized_SBP_operators(rd)
  Qrst_skew = map(A -> 0.5 * (A - A'), Qrst_hybridized)
  return Qrst_skew, VhP, Ph
end

function compute_flux_differencing_SBP_matrices(dg::DGMultiFluxDiff{<:SBP})
  rd = dg.basis
  @unpack M, Drst, Pq = rd
  Qrst = map(D -> M * D, Drst)
  Qrst_skew = map(A -> 0.5 * (A - A'), Qrst)
  return Qrst_skew
end

# precompute sparsity pattern for optimized flux differencing routines for tensor product elements
function compute_sparsity_pattern(flux_diff_matrices, dg::DG,
                                  tol = 100 * eps(real(dg))) where {DG <: DGMultiFluxDiff{ApproxType, <:Union{Quad, Hex}}} where {ApproxType}
  sparsity_pattern = sum(map(A->abs.(A), droptol!.(sparse.(flux_diff_matrices), tol))) .!= 0
  return sparsity_pattern
end

compute_sparsity_pattern(flux_diff_matrices, dg::DGMulti) = nothing

# For flux differencing SBP-type approximations, store solutions in Matrix{SVector{nvars}}.
# This results in a slight speedup for `calc_volume_integral!`.
function allocate_nested_array(uEltype, nvars, array_dimensions, dg::DGMultiFluxDiff{<:SBP})
  return zeros(SVector{nvars, uEltype}, array_dimensions...)
end

function create_cache(mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:SBP}, RealT, uEltype)

  rd = dg.basis
  md = mesh.md

  # for use with flux differencing schemes
  Qrst_skew = compute_flux_differencing_SBP_matrices(dg)
  # `sparsity_pattern` is the global sparsity pattern summing
  # the individual `sparsity_patterns` of each spatial derivative
  sparsity_pattern = compute_sparsity_pattern(Qrst_skew, dg)
  sparsity_patterns = map(Qi -> compute_sparsity_pattern((Qi,), dg), Qrst_skew)

  nvars = nvariables(equations)

  # Todo: DGMulti. Factor common storage into a struct (MeshDataCache?) for reuse across solvers?
  # storage for volume quadrature values, face quadrature values, flux values
  u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
  u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
  flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
  lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators

  local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) for _ in 1:Threads.nthreads()]

  # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
  fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, rd.Nq) for _ in 1:Threads.nthreads()]

  return (; md, Qrst_skew, sparsity_pattern, sparsity_patterns,
            invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
            u_values, u_face_values, flux_face_values,
            local_values_threaded, fluxdiff_local_threaded)
end

function create_cache(mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:Polynomial}, RealT, uEltype)

  rd = dg.basis
  @unpack md = mesh

  Qrst_skew, VhP, Ph = compute_flux_differencing_SBP_matrices(dg)
  # `sparsity_pattern` is the global sparsity pattern summing
  # the individual `sparsity_patterns` of each spatial derivative
  sparsity_pattern = compute_sparsity_pattern(Qrst_skew, dg)
  sparsity_patterns = map(Qi -> compute_sparsity_pattern((Qi,), dg), Qrst_skew)

  nvars = nvariables(equations)

  # temp storage for entropy variables at volume quad points
  entropy_var_values = allocate_nested_array(uEltype, nvars, (rd.Nq, md.num_elements), dg)

  # storage for all quadrature points (concatenated volume / face quadrature points)
  num_quad_points_total = rd.Nq + rd.Nfq
  entropy_projected_u_values = allocate_nested_array(uEltype, nvars, (num_quad_points_total, md.num_elements), dg)
  projected_entropy_var_values = allocate_nested_array(uEltype, nvars, (num_quad_points_total, md.num_elements), dg)

  # For this specific solver, `prolong2interfaces` will not be used anymore.
  # Instead, this step is also performed in `entropy_projection!`. Thus, we set
  # `u_face_values` as a `view` into `entropy_projected_u_values`. We do not do
  # the same for `u_values` since we will use that with LoopVectorization, which
  # cannot handle such views as of v0.12.66, the latest version at the time of writing.
  u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
  u_face_values = view(entropy_projected_u_values, rd.Nq+1:num_quad_points_total, :)
  flux_face_values = similar(u_face_values)

  # local storage for interface fluxes, rhs, and source
  local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) for _ in 1:Threads.nthreads()]

  # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
  # The result is then transferred to rhs_local_threaded::StructArray{<:SVector} before
  # projecting it and storing it into `du`.
  fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, num_quad_points_total) for _ in 1:Threads.nthreads()]
  rhs_local_threaded = [allocate_nested_array(uEltype, nvars, (num_quad_points_total,), dg)  for _ in 1:Threads.nthreads()]

  return (; md, Qrst_skew, sparsity_pattern, sparsity_patterns,
            VhP, Ph, invJ = inv.(md.J),
            entropy_var_values, projected_entropy_var_values, entropy_projected_u_values,
            u_values, u_face_values,  flux_face_values,
            local_values_threaded, fluxdiff_local_threaded, rhs_local_threaded)
end

# TODO: DGMulti. Address hard-coding of `entropy2cons!` and `cons2entropy!` for this function.
function entropy_projection!(cache, u, mesh::VertexMappedMesh, equations, dg::DGMulti)

  rd = dg.basis
  @unpack Vq = rd
  @unpack VhP, entropy_var_values, u_values = cache
  @unpack projected_entropy_var_values, entropy_projected_u_values = cache

  apply_to_each_field(mul_by!(Vq), u_values, u)

  # TODO: DGMulti. `@threaded` crashes when using `eachindex(u_values)`.
  # See https://github.com/JuliaSIMD/Polyester.jl/issues/37 for more details.
  @threaded for i in Base.OneTo(length(u_values))
    entropy_var_values[i] = cons2entropy(u_values[i], equations)
  end

  # "VhP" fuses the projection "P" with interpolation to volume and face quadrature "Vh"
  apply_to_each_field(mul_by!(VhP), projected_entropy_var_values, entropy_var_values)

  # TODO: DGMulti. `@threaded` crashes when using `eachindex(projected_entropy_var_values)`.
  # See https://github.com/JuliaSIMD/Polyester.jl/issues/37 for more details.
  @threaded for i in Base.OneTo(length(projected_entropy_var_values))
    entropy_projected_u_values[i] = entropy2cons(projected_entropy_var_values[i], equations)
  end
end

function calc_volume_integral!(du, u, volume_integral,
                               mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:SBP}, cache)

  @unpack fluxdiff_local_threaded, sparsity_patterns, inv_wq, Qrst_skew = cache
  @unpack volume_flux = volume_integral

  # Todo: DGMulti. Dispatch on curved/non-curved mesh types, this code only works for affine meshes (accessing rxJ[1,e],...)
  @threaded for e in eachelement(mesh, dg, cache)
    fluxdiff_local = fluxdiff_local_threaded[Threads.threadid()]
    fill!(fluxdiff_local, zero(eltype(fluxdiff_local)))
    u_local = view(u, :, e)
    for dim in eachdim(mesh)
      # There are two ways to write this flux differencing discretization on affine meshes.
      # 1. Use numerical fluxes in Cartesian directions and sum up the discrete derivative
      #    operators per coordinate direction accordingly.
      # 2. Use discrete derivative operators per coordinate direction and corresponding
      #    numerical fluxes in arbitrary (non-Cartesian) space directions.
      #
      # The first option can be implemented using
      #
      #   Q_skew = build_lazy_physical_derivative(e, dim, mesh, dg, cache)
      #   hadamard_sum!(fluxdiff_local, Q_skew, volume_flux, dim,
      #                 u_local, equations, sparsity_pattern)
      #
      # with `sparsity_pattern === cache.sparsity_pattern`.
      # However, this option makes it necessary to sum up the individual
      # `sparsity_patterns` of each reference coordinate direction. On tensor-product
      # elements such as `Quad()` or `Hex()` elements, this increases the number of
      # potentially expensive numerical flux evaluations by a factor of `ndims(mesh)`.
      # Thus, we use the second option below (which basically corresponds to the
      # well-known sum factorization on tensor product elements).
      # Note that there is basically no difference for dense derivative operators.
      normal_direction = get_contravariant_vector(e, dim, mesh)
      sparsity_pattern = sparsity_patterns[dim]
      Q_skew = Qrst_skew[dim]

      # Use a `sparsity_pattern` to dispatch `hadamard_sum!`.
      # If using `Tri()` or `Tet()` elements, the `Q_skew` matrices are dense,
      # and `sparsity_pattern === nothing`. If using `Quad()` or `Hex()` elements
      # with an `SBP` `approximation_type`, then `sparsity_pattern::AbstractSparseMatrix{Bool}`.
      hadamard_sum!(fluxdiff_local, Q_skew, volume_flux, normal_direction,
                    u_local, equations, sparsity_pattern)
    end

    for i in each_quad_node(mesh, dg, cache)
      du[i, e] = du[i, e] + fluxdiff_local[i] * inv_wq[i]
    end
  end
end

function calc_volume_integral!(du, u, volume_integral,
                               mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:Polynomial}, cache)

  rd = dg.basis
  @unpack entropy_projected_u_values, Ph, sparsity_pattern = cache
  @unpack fluxdiff_local_threaded, rhs_local_threaded = cache
  @unpack volume_flux = volume_integral

  # skips subblock of Qi_skew which we know is zero by construction
  skip_index(i,j) = i > rd.Nq && j > rd.Nq

  # Todo: DGMulti. Dispatch on curved/non-curved mesh types, this code only works for affine meshes (accessing rxJ[1,e],...)
  @threaded for e in eachelement(mesh, dg, cache)
    fluxdiff_local = fluxdiff_local_threaded[Threads.threadid()]
    fill!(fluxdiff_local, zero(eltype(fluxdiff_local)))
    u_local = view(entropy_projected_u_values, :, e)
    for i in eachdim(mesh)
      Qi_skew = build_lazy_physical_derivative(e, i, mesh, dg, cache)
      hadamard_sum!(fluxdiff_local, Qi_skew, volume_flux, i,
                    u_local, equations, sparsity_pattern, skip_index)
    end

    # convert fluxdiff_local::Vector{<:SVector} to StructArray{<:SVector} for faster
    # apply_to_each_field performance.
    rhs_local = rhs_local_threaded[Threads.threadid()]
    for i in Base.OneTo(length(fluxdiff_local))
      rhs_local[i] = fluxdiff_local[i]
    end
    apply_to_each_field(mul_by_accum!(Ph), view(du, :, e), rhs_local)
  end
end

# Specialize since `u_values` isn't computed for DGMultiFluxDiff{<:SBP} solvers.
function calc_sources!(du, u, t, source_terms,
                       mesh, equations, dg::DGMultiFluxDiff{<:SBP}, cache)

  rd = dg.basis
  md = mesh.md

  @threaded for e in eachelement(mesh, dg, cache)
    for i in each_quad_node(mesh, dg, cache)
      du[i, e] += source_terms(u[i, e], getindex.(md.xyzq, i, e), t, equations)
    end
  end
end

# Specializes on Polynomial (e.g., modal) DG methods with a flux differencing volume kernel, e.g.,
# an entropy conservative/stable discretization. For modal DG schemes, an extra `entropy_projection`
# is required (see https://doi.org/10.1016/j.jcp.2018.02.033, Section 4.3).
function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions::BC,
              source_terms::Source, dg::DGMultiFluxDiff{Polynomial}, cache) where {Source, BC}

  @trixi_timeit timer() "Reset du/dt" fill!(du, zero(eltype(du)))

  # this function evaluates the solution at volume and face quadrature points (which was previously
  # done in `prolong2interfaces` and `calc_volume_integral`)
  @trixi_timeit timer() "entropy_projection!" entropy_projection!(cache, u, mesh, equations, dg)

  @trixi_timeit timer() "calc_volume_integral!" calc_volume_integral!(du, u, dg.volume_integral,
                                                                      mesh, equations, dg, cache)

  # the following functions are the same as in VolumeIntegralWeakForm, and can be reused from dg.jl
  @trixi_timeit timer() "calc_interface_flux!" calc_interface_flux!(cache, dg.surface_integral,
                                                                    mesh, equations, dg)

  @trixi_timeit timer() "calc_boundary_flux!" calc_boundary_flux!(cache, t, boundary_conditions,
                                                                  mesh, equations, dg)

  @trixi_timeit timer() "calc_surface_integral!" calc_surface_integral!(du, u, dg.surface_integral,
                                                                        mesh, equations, dg, cache)

  @trixi_timeit timer() "invert_jacobian" invert_jacobian!(du, mesh, equations, dg, cache)

  @trixi_timeit timer() "calc_sources!" calc_sources!(du, u, t, source_terms,
                                                      mesh, equations, dg, cache)

  return nothing
end

end # @muladd


text:
# SubgridScale parameters
const AtmosSGS = CLIMAParameters.Atmos.SubgridScale

AtmosSGS.C_smag(::AEPS) = 0.21
AtmosSGS.C_drag(::AEPS) = 0.0011
AtmosSGS.inv_Pr_turb(::AEPS) = 3
AtmosSGS.Prandtl_air(::AEPS) = 71 / 100
AtmosSGS.c_a_KASM(::AEPS) = 0.10
AtmosSGS.c_e1_KASM(::AEPS) = 0.19
AtmosSGS.c_e2_KASM(::AEPS) = 0.51
AtmosSGS.c_1_KASM(ps::AEPS) = AtmosSGS.c_a_KASM(ps) * 0.76^2
AtmosSGS.c_2_KASM(ps::AEPS) = AtmosSGS.c_e2_KASM(ps) + 2 * AtmosSGS.c_1_KASM(ps)
AtmosSGS.c_3_KASM(ps::AEPS) = AtmosSGS.c_a_KASM(ps)^(3 / 2)

# EDMF parameters
const EDMF = CLIMAParameters.Atmos.EDMF

# Entrainment - detrainment model
EDMF.c_λ(::AEPS) = 0.3
EDMF.c_ε(::AEPS) = 0.13
EDMF.c_δ(::AEPS) = 0.52
EDMF.c_γ(::AEPS) = 0.075
EDMF.β(::AEPS) = 2
EDMF.μ_0(::AEPS) = 4e-4
EDMF.χ(::AEPS) = 0.25
EDMF.w_min(::AEPS) = 0.1
EDMF.lim_ϵ(::AEPS) = 1e-4
EDMF.lim_amp(::AEPS) = 10

# Subdomain model
EDMF.a_min(::AEPS) = 0.001

# Surface model
EDMF.a_surf(::AEPS) = 0.1
EDMF.κ_star²(::AEPS) = 3.75
EDMF.ψϕ_stab(::AEPS) = 8.3

# Pressure model
EDMF.α_d(::AEPS) = 10.0
EDMF.α_a(::AEPS) = 0.1
EDMF.α_b(::AEPS) = 0.12
EDMF.H_up_min(::AEPS) = 500

# Mixing length model
EDMF.c_d(::AEPS) = 0.22
EDMF.c_m(::AEPS) = 0.14
EDMF.c_b(::AEPS) = 0.4
EDMF.a1(::AEPS) = 0.2
EDMF.a2(::AEPS) = 100
EDMF.ω_pr(::AEPS) = 53.0 / 13.0
EDMF.Pr_n(::AEPS) = 0.74
EDMF.Ri_c(::AEPS) = 0.25
EDMF.smin_ub(::AEPS) = 0.1
EDMF.smin_rm(::AEPS) = 1.5

# 0-moment microphysics parameters
const Microphysics_0M = CLIMAParameters.Atmos.Microphysics_0M

Microphysics_0M.τ_precip(::AEPS) = 1000
Microphysics_0M.qc_0(::AEPS) = 5e-3
Microphysics_0M.S_0(::AEPS) = 0.02

# 1-moment microphysics parameters
const Microphysics = CLIMAParameters.Atmos.Microphysics

# general
Microphysics.C_drag(::AEPS) = 0.55
Microphysics.K_therm(::AEPS) = 2.4e-2
Microphysics.D_vapor(::AEPS) = 2.26e-5
Microphysics.ν_air(::AEPS) = 1.6e-5
Microphysics.N_Sc(ps::AEPS) = Microphysics.ν_air(ps) / Microphysics.D_vapor(ps)

# cloud water
Microphysics.τ_cond_evap(::AEPS) = 10

# cloud ice
Microphysics.τ_sub_dep(::AEPS) = 10
Microphysics.r_ice_snow(::AEPS) = 62.5 * 1e-6
Microphysics.n0_ice(::AEPS) = 1e7 * 2
Microphysics.r0_ice(::AEPS) = 1e-5
Microphysics.me_ice(::AEPS) = 3
Microphysics.m0_ice(ps::AEPS) =
    4 / 3 *
    π *
    CLIMAParameters.Planet.ρ_cloud_ice(ps) *
    Microphysics.r0_ice(ps)^Microphysics.me_ice(ps)
Microphysics.χm_ice(::AEPS) = 1
Microphysics.Δm_ice(::AEPS) = 0

# rain
Microphysics.q_liq_threshold(::AEPS) = 5e-4
Microphysics.τ_acnv_rai(::AEPS) = 1e3
Microphysics.a_vent_rai(::AEPS) = 1.5
Microphysics.b_vent_rai(::AEPS) = 0.53
Microphysics.n0_rai(::AEPS) = 8e6 * 2
Microphysics.r0_rai(::AEPS) = 1e-3
Microphysics.me_rai(::AEPS) = 3
Microphysics.ae_rai(::AEPS) = 2
Microphysics.ve_rai(::AEPS) = 0.5
Microphysics.m0_rai(ps::AEPS) =
    4 / 3 *
    π *
    CLIMAParameters.Planet.ρ_cloud_liq(ps) *
    Microphysics.r0_rai(ps)^Microphysics.me_rai(ps)
Microphysics.a0_rai(ps::AEPS) =
    π * Microphysics.r0_rai(ps)^Microphysics.ae_rai(ps)
Microphysics.χm_rai(::AEPS) = 1
Microphysics.Δm_rai(::AEPS) = 0
Microphysics.χa_rai(::AEPS) = 1
Microphysics.Δa_rai(::AEPS) = 0
Microphysics.χv_rai(::AEPS) = 1
Microphysics.Δv_rai(::AEPS) = 0

# snow
Microphysics.q_ice_threshold(::AEPS) = 1e-6
Microphysics.τ_acnv_sno(::AEPS) = 1e2
Microphysics.a_vent_sno(::AEPS) = 0.65
Microphysics.b_vent_sno(::AEPS) = 0.44
# n0_sno = μ_sno (ρ q_sno / ρ_0)^ν_sno;  ρ_0 = 1kg/m3
Microphysics.μ_sno(::AEPS) = 4.36 * 1e9
Microphysics.ν_sno(::AEPS) = 0.63
Microphysics.r0_sno(::AEPS) = 1e-3
Microphysics.me_sno(::AEPS) = 2
Microphysics.ae_sno(::AEPS) = 2
Microphysics.ve_sno(::AEPS) = 0.25
Microphysics.m0_sno(ps::AEPS) =
    1e-1 * Microphysics.r0_sno(ps)^Microphysics.me_sno(ps)
Microphysics.a0_sno(ps::AEPS) =
    0.3 * π * Microphysics.r0_sno(ps)^Microphysics.ae_sno(ps)
Microphysics.v0_sno(ps::AEPS) =
    2^(9 / 4) * Microphysics.r0_sno(ps)^Microphysics.ve_sno(ps)
Microphysics.χm_sno(::AEPS) = 1
Microphysics.Δm_sno(::AEPS) = 0
Microphysics.χa_sno(::AEPS) = 1
Microphysics.Δa_sno(::AEPS) = 0
Microphysics.χv_sno(::AEPS) = 1
Microphysics.Δv_sno(::AEPS) = 0

Microphysics.microph_scaling(::AEPS) = 1.0
Microphysics.microph_scaling_dep_sub(::AEPS) = 1.0
Microphysics.microph_scaling_melt(::AEPS) = 1.0

# interactions
Microphysics.E_liq_rai(ps::AEPS) = 0.8
Microphysics.E_liq_sno(ps::AEPS) = 0.1
Microphysics.E_ice_rai(ps::AEPS) = 1
Microphysics.E_ice_sno(ps::AEPS) = 0.1
Microphysics.E_rai_sno(ps::AEPS) = 1

ext: jl

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
using Test

using CUDAnative, CuArrays, CUDAdrv
import Adapt, LLVM

include("util.jl")

@testset "CUDAnative" begin

@test CUDAnative.functional()

CUDAnative.version()

CUDAnative.enable_timings()

# see if we have a device to run tests on
# (do this early so that the codegen tests can target the same compute capability)
global dev, cap
dev = nothing
@test length(devices()) > 0
if length(devices()) > 0
    # the API shouldn't have been initialized
    @test CuCurrentContext() == nothing

    device_callbacked = nothing
    device_callback = (dev, ctx) -> begin
        device_callbacked = dev
    end
    push!(CUDAnative.device!_listeners, device_callback)

    # now cause initialization
    Mem.alloc(Mem.Device, 1)
    @test CuCurrentContext() != nothing
    @test device(CuCurrentContext()) == CuDevice(0)
    @test device_callbacked == CuDevice(0)

    device!(CuDevice(0))
    device!(CuDevice(0)) do
        nothing
    end

    device_reset!()

    # test the device selection functionality
    if length(devices()) > 1
        device!(1) do
            @test device(CuCurrentContext()) == CuDevice(1)
        end
        @test device(CuCurrentContext()) == CuDevice(0)

        device!(1)
        @test device(CuCurrentContext()) == CuDevice(1)
    end

    # pick a suiteable device (by available memory,
    # but also by capability if testing needs to be thorough)
    candidates = [(dev=dev,
                   cap=capability(dev),
                   mem=CuContext(ctx->CUDAdrv.available_memory(), dev))
                  for dev in devices()]
    thorough = parse(Bool, get(ENV, "CI_THOROUGH", "false"))
    if thorough
        sort!(candidates, by=x->(x.cap, x.mem))
    else
        sort!(candidates, by=x->x.mem)
    end
    pick = last(candidates)
    @info("Testing using device $(name(pick.dev)) (compute capability $(pick.cap), $(Base.format_bytes(pick.mem)) available memory) on CUDA driver $(CUDAdrv.version()) and toolkit $(CUDAnative.version())")
    device!(pick.dev)
    dev = pick.dev
end
cap = CUDAnative.current_capability()

include("base.jl")
include("pointer.jl")
include("codegen.jl")

if dev === nothing
    @warn("No CUDA-capable devices available; skipping on-device tests.")
else
    if capability(dev) < v"2.0"
        @warn("native execution not supported on SM < 2.0")
    else
        include("device/codegen.jl")
        include("device/execution.jl")
        include("device/pointer.jl")
        include("device/array.jl")
        include("device/cuda.jl")

        include("examples.jl")
    end
end

CUDAnative.timings()

end


text:
using PkgDashboards
using Test

@testset "PkgDashboards.jl" begin

    db = dashboard("mfalt", output=:markdown, autoopen=false)
    @test db isa PkgDashboards.DashBoard
    db = dashboard("mfalt", output=:html, autoopen=false)
    @testset "show" begin
        @info "Testing show"
        io = IOBuffer()
        for m in [MIME"text/plain",MIME"text/markdown",MIME"text/html"]
            @test_nowarn show(io, m(), db)
        end
    end

    @test_nowarn db = dashboard("mfalt", output=:markdown, autoopen=false, stargazers=true, activity=true, githubci=true)
    @testset "pkgdashboard" begin
        @info "Testing pkgdashboard"
        @test_nowarn db = pkgdashboard("Hyperopt", output=:markdown, autoopen=false)
    end
end


text:
@testset "predicates" begin
    @testset "IsEqual" begin
        pred = IsEqual{(:a, :b), Tuple{Int, Float64}}((1, 2.0))
        @test IsEqual(a=1, b=2.0) === pred

        @test pred((a=1, b=2.0)) === true
        @test pred((a=1, b=2.1)) === false
        @test pred((b=2.0, a=1)) === true
        @test pred((b=2.0, a=2)) === false
        @test pred((a=1, b=2.0, c=false)) === true
        @test pred((a=1, b=2.1, c=false)) === false
        @test_throws Exception pred((a=1))

        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

        @test map(pred, t) == [false, false, true]
        @test findall(pred, t) == [3]
        @test filter(pred, t) == Table(a = [1], b = [2.0])

        t2a = accelerate(t, UniqueIndex{(:a,)})
        @test map(pred, t2a) == [false, false, true]
        @test findall(pred, t2a) == [3]
        @test filter(pred, t2a) == Table(a = [1], b = [2.0])
        t2b = accelerate(t, UniqueIndex{(:b,)})
        @test map(pred, t2b) == [false, false, true]
        @test findall(pred, t2b) == [3]
        @test filter(pred, t2b) == Table(a = [1], b = [2.0])
        t2ab = accelerate(t, UniqueIndex{(:a,:b)})
        @test map(pred, t2ab) == [false, false, true]
        @test findall(pred, t2ab) == [3]
        @test filter(pred, t2ab) == Table(a = [1], b = [2.0])
        t2ba = accelerate(t, UniqueIndex{(:b,:a)})
        @test map(pred, t2ba) == [false, false, true]
        @test findall(pred, t2ba) == [3]
        @test filter(pred, t2ba) == Table(a = [1], b = [2.0])

        t3a = accelerate(t, SortIndex{(:a,)})
        @test map(pred, t3a) == [false, false, true]
        @test findall(pred, t3a) == [3]
        @test filter(pred, t3a) == Table(a = [1], b = [2.0])
        t3b = accelerate(t, SortIndex{(:b,)})
        @test map(pred, t3b) == [false, false, true]
        @test findall(pred, t3b) == [3]
        @test filter(pred, t3b) == Table(a = [1], b = [2.0])
        t3ab = accelerate(t, SortIndex{(:a,:b)})
        @test map(pred, t3ab) == [false, false, true]
        @test findall(pred, t3ab) == [3]
        @test filter(pred, t3ab) == Table(a = [1], b = [2.0])
        t3ba = accelerate(t, SortIndex{(:b,:a)})
        @test map(pred, t3ba) == [false, false, true]
        @test findall(pred, t3ba) == [3]
        @test filter(pred, t3ba) == Table(a = [1], b = [2.0])

        t4a = accelerate(t, UniqueSortIndex{(:a,)})
        @test map(pred, t4a) == [false, false, true]
        @test findall(pred, t4a) == [3]
        @test filter(pred, t4a) == Table(a = [1], b = [2.0])
        t4b = accelerate(t, UniqueSortIndex{(:b,)})
        @test map(pred, t4b) == [false, false, true]
        @test findall(pred, t4b) == [3]
        @test filter(pred, t4b) == Table(a = [1], b = [2.0])
        t4ab = accelerate(t, UniqueSortIndex{(:a,:b)})
        @test map(pred, t4ab) == [false, false, true]
        @test findall(pred, t4ab) == [3]
        @test filter(pred, t4ab) == Table(a = [1], b = [2.0])
        t4ba = accelerate(t, UniqueSortIndex{(:b,:a)})
        @test map(pred, t4ba) == [false, false, true]
        @test findall(pred, t4ba) == [3]
        @test filter(pred, t4ba) == Table(a = [1], b = [2.0])

        t5a = accelerate(t, HashIndex{(:a,)})
        @test map(pred, t5a) == [false, false, true]
        @test findall(pred, t5a) == [3]
        @test filter(pred, t5a) == Table(a = [1], b = [2.0])
        t5b = accelerate(t, HashIndex{(:b,)})
        @test map(pred, t5b) == [false, false, true]
        @test findall(pred, t5b) == [3]
        @test filter(pred, t5b) == Table(a = [1], b = [2.0])
        t5ab = accelerate(t, HashIndex{(:a,:b)})
        @test map(pred, t5ab) == [false, false, true]
        @test findall(pred, t5ab) == [3]
        @test filter(pred, t5ab) == Table(a = [1], b = [2.0])
        t5ba = accelerate(t, HashIndex{(:b,:a)})
        @test map(pred, t5ba) == [false, false, true]
        @test findall(pred, t5ba) == [3]
        @test filter(pred, t5ba) == Table(a = [1], b = [2.0])

        t6a = accelerate(t, UniqueHashIndex{(:a,)})
        @test map(pred, t6a) == [false, false, true]
        @test findall(pred, t6a) == [3]
        @test filter(pred, t6a) == Table(a = [1], b = [2.0])
        t6b = accelerate(t, UniqueHashIndex{(:b,)})
        @test map(pred, t6b) == [false, false, true]
        @test findall(pred, t6b) == [3]
        @test filter(pred, t6b) == Table(a = [1], b = [2.0])
        t6ab = accelerate(t, UniqueHashIndex{(:a,:b)})
        @test map(pred, t6ab) == [false, false, true]
        @test findall(pred, t6ab) == [3]
        @test filter(pred, t6ab) == Table(a = [1], b = [2.0])
        t6ba = accelerate(t, UniqueHashIndex{(:b,:a)})
        @test map(pred, t6ba) == [false, false, true]
        @test findall(pred, t6ba) == [3]
        @test filter(pred, t6ba) == Table(a = [1], b = [2.0])
    end

    @testset "IsLess" begin
        pred = IsLess{(:a, :b), Tuple{Int, Float64}}((2, 4.0))
        @test IsLess(a=2, b=4.0) === pred

        @test pred((a=1, b=2.0)) === true
        @test pred((a=2, b=4.0)) === false
        @test pred((a=2, b=3.9)) === true
        @test pred((b=2.0, a=1)) === true
        @test pred((b=4.0, a=2)) === false
        @test pred((a=1, b=2.0, c=false)) === true
        @test pred((a=3, b=2.0, c=false)) === false
        @test_throws Exception pred((a=1))

        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

        @test map(pred, t) == [false, false, true]
        @test findall(pred, t) == [3]
        @test filter(pred, t) == Table(a = [1], b = [2.0])

        t3a = accelerate(t, SortIndex{(:a,)})
        @test map(pred, t3a) == [false, false, true]
        @test findall(pred, t3a) == [3]
        @test filter(pred, t3a) == Table(a = [1], b = [2.0])
        t3b = accelerate(t, SortIndex{(:b,)})
        @test map(pred, t3b) == [false, false, true]
        @test findall(pred, t3b) == [3]
        @test filter(pred, t3b) == Table(a = [1], b = [2.0])
        t3ab = accelerate(t, SortIndex{(:a,:b)})
        @test map(pred, t3ab) == [false, false, true]
        @test findall(pred, t3ab) == [3]
        @test filter(pred, t3ab) == Table(a = [1], b = [2.0])
        t3ba = accelerate(t, SortIndex{(:b,:a)})
        @test map(pred, t3ba) == [false, false, true]
        @test findall(pred, t3ba) == [3]
        @test filter(pred, t3ba) == Table(a = [1], b = [2.0])

        t4a = accelerate(t, UniqueSortIndex{(:a,)})
        @test map(pred, t4a) == [false, false, true]
        @test findall(pred, t4a) == [3]
        @test filter(pred, t4a) == Table(a = [1], b = [2.0])
        t4b = accelerate(t, UniqueSortIndex{(:b,)})
        @test map(pred, t4b) == [false, false, true]
        @test findall(pred, t4b) == [3]
        @test filter(pred, t4b) == Table(a = [1], b = [2.0])
        t4ab = accelerate(t, UniqueSortIndex{(:a,:b)})
        @test map(pred, t4ab) == [false, false, true]
        @test findall(pred, t4ab) == [3]
        @test filter(pred, t4ab) == Table(a = [1], b = [2.0])
        t4ba = accelerate(t, UniqueSortIndex{(:b,:a)})
        @test map(pred, t4ba) == [false, false, true]
        @test findall(pred, t4ba) == [3]
        @test filter(pred, t4ba) == Table(a = [1], b = [2.0])
    end

    @testset "IsLessEqual" begin
        pred = IsLessEqual{(:a, :b), Tuple{Int, Float64}}((2, 4.0))
        @test IsLessEqual(a=2, b=4.0) === pred

        @test pred((a=1, b=2.0)) === true
        @test pred((a=2, b=4.0)) === true
        @test pred((a=2, b=4.1)) === false
        @test pred((b=2.0, a=1)) === true
        @test pred((b=4.1, a=2)) === false
        @test pred((a=2, b=4.0, c=false)) === true
        @test pred((a=3, b=2.0, c=false)) === false
        @test_throws Exception pred((a=1))

        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

        @test map(pred, t) == [true, false, true]
        @test findall(pred, t) == [1, 3]
        @test filter(pred, t) == Table(a = [2, 1], b = [4.0, 2.0])

        t3a = accelerate(t, SortIndex{(:a,)})
        @test map(pred, t3a) == [true, false, true]
        @test findall(pred, t3a) == [1, 3]
        @test filter(pred, t3a) == Table(a = [2, 1], b = [4.0, 2.0])
        t3b = accelerate(t, SortIndex{(:b,)})
        @test map(pred, t3b) == [true, false, true]
        @test findall(pred, t3b) == [1, 3]
        @test filter(pred, t3b) == Table(a = [2, 1], b = [4.0, 2.0])
        t3ab = accelerate(t, SortIndex{(:a,:b)})
        @test map(pred, t3ab) == [true, false, true]
        @test findall(pred, t3ab) == [1, 3]
        @test filter(pred, t3ab) == Table(a = [2, 1], b = [4.0, 2.0])
        t3ba = accelerate(t, SortIndex{(:b,:a)})
        @test map(pred, t3ba) == [true, false, true]
        @test findall(pred, t3ba) == [1, 3]
        @test filter(pred, t3ba) == Table(a = [2, 1], b = [4.0, 2.0])

        t4a = accelerate(t, UniqueSortIndex{(:a,)})
        @test map(pred, t4a) == [true, false, true]
        @test findall(pred, t4a) == [1, 3]
        @test filter(pred, t4a) == Table(a = [2, 1], b = [4.0, 2.0])
        t4b = accelerate(t, UniqueSortIndex{(:b,)})
        @test map(pred, t4b) == [true, false, true]
        @test findall(pred, t4b) == [1, 3]
        @test filter(pred, t4b) == Table(a = [2, 1], b = [4.0, 2.0])
        t4ab = accelerate(t, UniqueSortIndex{(:a,:b)})
        @test map(pred, t4ab) == [true, false, true]
        @test findall(pred, t4ab) == [1, 3]
        @test filter(pred, t4ab) == Table(a = [2, 1], b = [4.0, 2.0])
        t4ba = accelerate(t, UniqueSortIndex{(:b,:a)})
        @test map(pred, t4ba) == [true, false, true]
        @test findall(pred, t4ba) == [1, 3]
        @test filter(pred, t4ba) == Table(a = [2, 1], b = [4.0, 2.0])
    end

    @testset "IsGreater" begin
        pred = IsGreater{(:a, :b), Tuple{Int, Float64}}((2, 4.0))
        @test IsGreater(a=2, b=4.0) === pred

        @test pred((a=3, b=6.0)) === true
        @test pred((a=2, b=4.0)) === false
        @test pred((a=2, b=4.1)) === true
        @test pred((b=2.0, a=1)) === false
        @test pred((b=4.1, a=2)) === true
        @test pred((a=1, b=2.0, c=false)) === false
        @test pred((a=3, b=2.0, c=false)) === true
        @test_throws Exception pred((a=1))

        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

        @test map(pred, t) == [false, true, false]
        @test findall(pred, t) == [2]
        @test filter(pred, t) == Table(a = [3], b = [6.0])

        t3a = accelerate(t, SortIndex{(:a,)})
        @test map(pred, t3a) == [false, true, false]
        @test findall(pred, t3a) == [2]
        @test filter(pred, t3a) == Table(a = [3], b = [6.0])
        t3b = accelerate(t, SortIndex{(:b,)})
        @test map(pred, t3b) == [false, true, false]
        @test findall(pred, t3b) == [2]
        @test filter(pred, t3b) == Table(a = [3], b = [6.0])
        t3ab = accelerate(t, SortIndex{(:a,:b)})
        @test map(pred, t3ab) == [false, true, false]
        @test findall(pred, t3ab) == [2]
        @test filter(pred, t3ab) == Table(a = [3], b = [6.0])
        t3ba = accelerate(t, SortIndex{(:b,:a)})
        @test map(pred, t3ba) == [false, true, false]
        @test findall(pred, t3ba) == [2]
        @test filter(pred, t3ba) == Table(a = [3], b = [6.0])

        t4a = accelerate(t, UniqueSortIndex{(:a,)})
        @test map(pred, t4a) == [false, true, false]
        @test findall(pred, t4a) == [2]
        @test filter(pred, t4a) == Table(a = [3], b = [6.0])
        t4b = accelerate(t, UniqueSortIndex{(:b,)})
        @test map(pred, t4b) == [false, true, false]
        @test findall(pred, t4b) == [2]
        @test filter(pred, t4b) == Table(a = [3], b = [6.0])
        t4ab = accelerate(t, UniqueSortIndex{(:a,:b)})
        @test map(pred, t4ab) == [false, true, false]
        @test findall(pred, t4ab) == [2]
        @test filter(pred, t4ab) == Table(a = [3], b = [6.0])
        t4ba = accelerate(t, UniqueSortIndex{(:b,:a)})
        @test map(pred, t4ba) == [false, true, false]
        @test findall(pred, t4ba) == [2]
        @test filter(pred, t4ba) == Table(a = [3], b = [6.0])
    end

    @testset "IsGreaterEqual" begin
        pred = IsGreaterEqual{(:a, :b), Tuple{Int, Float64}}((2, 4.0))
        @test IsGreaterEqual(a=2, b=4.0) === pred

        @test pred((a=2, b=3.9)) === false
        @test pred((a=2, b=4.0)) === true
        @test pred((a=4, b=6.0)) === true
        @test pred((b=2.0, a=1)) === false
        @test pred((b=4.1, a=2)) === true
        @test pred((a=2, b=4.0, c=false)) === true
        @test pred((a=3, b=2.0, c=false)) === true
        @test_throws Exception pred((a=1))

        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

        @test map(pred, t) == [true, true, false]
        @test findall(pred, t) == [1, 2]
        @test filter(pred, t) == Table(a = [2, 3], b = [4.0, 6.0])

        t3a = accelerate(t, SortIndex{(:a,)})
        @test map(pred, t3a) == [true, true, false]
        @test findall(pred, t3a) == [1, 2]
        @test filter(pred, t3a) == Table(a = [2, 3], b = [4.0, 6.0])
        t3b = accelerate(t, SortIndex{(:b,)})
        @test map(pred, t3b) == [true, true, false]
        @test findall(pred, t3b) == [1, 2]
        @test filter(pred, t3b) == Table(a = [2, 3], b = [4.0, 6.0])
        t3ab = accelerate(t, SortIndex{(:a,:b)})
        @test map(pred, t3ab) == [true, true, false]
        @test findall(pred, t3ab) == [1, 2]
        @test filter(pred, t3ab) == Table(a = [2, 3], b = [4.0, 6.0])
        t3ba = accelerate(t, SortIndex{(:b,:a)})
        @test map(pred, t3ba) == [true, true, false]
        @test findall(pred, t3ba) == [1, 2]
        @test filter(pred, t3ba) == Table(a = [2, 3], b = [4.0, 6.0])

        t4a = accelerate(t, UniqueSortIndex{(:a,)})
        @test map(pred, t4a) == [true, true, false]
        @test findall(pred, t4a) == [1, 2]
        @test filter(pred, t4a) == Table(a = [2, 3], b = [4.0, 6.0])
        t4b = accelerate(t, UniqueSortIndex{(:b,)})
        @test map(pred, t4b) == [true, true, false]
        @test findall(pred, t4b) == [1, 2]
        @test filter(pred, t4b) == Table(a = [2, 3], b = [4.0, 6.0])
        t4ab = accelerate(t, UniqueSortIndex{(:a,:b)})
        @test map(pred, t4ab) == [true, true, false]
        @test findall(pred, t4ab) == [1, 2]
        @test filter(pred, t4ab) == Table(a = [2, 3], b = [4.0, 6.0])
        t4ba = accelerate(t, UniqueSortIndex{(:b,:a)})
        @test map(pred, t4ba) == [true, true, false]
        @test findall(pred, t4ba) == [1, 2]
        @test filter(pred, t4ba) == Table(a = [2, 3], b = [4.0, 6.0])
    end

    @testset "In (Interval)" begin
        pred = In{(:a,), Tuple{Interval{Int}}}((2..3,))
        @test In(a=2..3) === pred

        @test pred((a=1, b=2.0)) === false
        @test pred((a=2, b=4.0)) === true
        @test pred((a=3, b=6.0)) === true
        @test pred((b=2.0, a=1)) === false
        @test pred((b=4.0, a=2)) === true
        @test pred((a=2, b=4.0, c=false)) === true
        @test pred((a=3, b=2.0, c=false)) === true
        @test_throws Exception pred((b=1.0))

        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

        @test map(pred, t) == [true, true, false]
        @test findall(pred, t) == [1, 2]
        @test filter(pred, t) == Table(a = [2, 3], b = [4.0, 6.0])

        t3a = accelerate(t, SortIndex{(:a,)})
        @test map(pred, t3a) == [true, true, false]
        @test findall(pred, t3a) == [1, 2]
        @test filter(pred, t3a) == Table(a = [2, 3], b = [4.0, 6.0])
        t3b = accelerate(t, SortIndex{(:b,)})
        @test map(pred, t3b) == [true, true, false]
        @test findall(pred, t3b) == [1, 2]
        @test filter(pred, t3b) == Table(a = [2, 3], b = [4.0, 6.0])
        t3ab = accelerate(t, SortIndex{(:a,:b)})
        @test map(pred, t3ab) == [true, true, false]
        @test findall(pred, t3ab) == [1, 2]
        @test filter(pred, t3ab) == Table(a = [2, 3], b = [4.0, 6.0])
        t3ba = accelerate(t, SortIndex{(:b,:a)})
        @test map(pred, t3ba) == [true, true, false]
        @test findall(pred, t3ba) == [1, 2]
        @test filter(pred, t3ba) == Table(a = [2, 3], b = [4.0, 6.0])

        t4a = accelerate(t, UniqueSortIndex{(:a,)})
        @test map(pred, t4a) == [true, true, false]
        @test findall(pred, t4a) == [1, 2]
        @test filter(pred, t4a) == Table(a = [2, 3], b = [4.0, 6.0])
        t4b = accelerate(t, UniqueSortIndex{(:b,)})
        @test map(pred, t4b) == [true, true, false]
        @test findall(pred, t4b) == [1, 2]
        @test filter(pred, t4b) == Table(a = [2, 3], b = [4.0, 6.0])
        t4ab = accelerate(t, UniqueSortIndex{(:a,:b)})
        @test map(pred, t4ab) == [true, true, false]
        @test findall(pred, t4ab) == [1, 2]
        @test filter(pred, t4ab) == Table(a = [2, 3], b = [4.0, 6.0])
        t4ba = accelerate(t, UniqueSortIndex{(:b,:a)})
        @test map(pred, t4ba) == [true, true, false]
        @test findall(pred, t4ba) == [1, 2]
        @test filter(pred, t4ba) == Table(a = [2, 3], b = [4.0, 6.0])
    end

    @testset "Equals" begin
        pred = Equals(:a, :b)

        @test pred((a=1, b=1)) === true
        @test pred((a=2, b=4.0)) === false
        @test pred((a=1, b=1, c=2.0)) === true
        @test pred((a=2, b=4.0, c=2.0)) === false
        @test pred((b=1, a=1)) === true
        @test pred((b=2, a=4.0)) === false

        va = [2, 3, 1]
        vb = [2, 1, 3]
        t = Table(a = va, b = vb)

        @test map(pred, t) == [true, false, false]
        @test findall(pred, t) == [1]
        @test filter(pred, t) == Table(a = [2], b = [2])

        @test Predicate(pred, (a=1,)) === IsEqual((b=1,))
        @test Predicate(pred, (b=1,)) === IsEqual((a=1,))
    end

    @testset "LessThan" begin
        pred = LessThan(:a, :b)

        @test pred((a=1, b=0)) === false
        @test pred((a=1, b=1)) === false
        @test pred((a=1, b=2)) === true
        @test pred((a=1, b=1, c=2.0)) === false
        @test pred((a=1, b=2, c=2.0)) === true
        @test pred((b=1, a=1)) === false
        @test pred((b=2, a=1)) === true

        va = [2, 3, 1]
        vb = [2, 1, 3]
        t = Table(a = va, b = vb)

        @test map(pred, t) == [false, false, true]
        @test findall(pred, t) == [3]
        @test filter(pred, t) == Table(a = [1], b = [3])

        @test Predicate(pred, (a=1,)) === IsGreaterEqual((b=1,))
        @test Predicate(pred, (b=1,)) === IsLess((a=1,))
    end

    @testset "LessEqualThan" begin
        pred = LessEqualThan(:a, :b)

        @test pred((a=1, b=0)) === false
        @test pred((a=1, b=1)) === true
        @test pred((a=1, b=2)) === true
        @test pred((a=1, b=0, c=2.0)) === false
        @test pred((a=1, b=1, c=2.0)) === true
        @test pred((a=1, b=2, c=2.0)) === true
        @test pred((b=0, a=1)) === false
        @test pred((b=1, a=1)) === true
        @test pred((b=2, a=1)) === true

        va = [2, 3, 1]
        vb = [2, 1, 3]
        t = Table(a = va, b = vb)

        @test map(pred, t) == [true, false, true]
        @test findall(pred, t) == [1, 3]
        @test filter(pred, t) == Table(a = [2, 1], b = [2, 3])

        @test Predicate(pred, (a=1,)) === IsGreater((b=1,))
        @test Predicate(pred, (b=1,)) === IsLessEqual((a=1,))
    end

    @testset "GreaterThan" begin
        pred = GreaterThan(:a, :b)

        @test pred((a=1, b=0)) === true
        @test pred((a=1, b=1)) === false
        @test pred((a=1, b=2)) === false
        @test pred((a=1, b=0, c=2.0)) === true
        @test pred((a=1, b=1, c=2.0)) === false
        @test pred((a=1, b=2, c=2.0)) === false
        @test pred((b=0, a=1)) === true
        @test pred((b=1, a=1)) === false
        @test pred((b=2, a=1)) === false

        va = [2, 3, 1]
        vb = [2, 1, 3]
        t = Table(a = va, b = vb)

        @test map(pred, t) == [false, true, false]
        @test findall(pred, t) == [2]
        @test filter(pred, t) == Table(a = [3], b = [1])

        @test Predicate(pred, (a=1,)) === IsLessEqual((b=1,))
        @test Predicate(pred, (b=1,)) === IsGreater((a=1,))
    end

    @testset "GreaterEqualThan" begin
        pred = GreaterEqualThan(:a, :b)

        @test pred((a=1, b=0)) === true
        @test pred((a=1, b=1)) === true
        @test pred((a=1, b=2)) === false
        @test pred((a=1, b=0, c=2.0)) === true
        @test pred((a=1, b=1, c=2.0)) === true
        @test pred((a=1, b=2, c=2.0)) === false
        @test pred((b=0, a=1)) === true
        @test pred((b=1, a=1)) === true
        @test pred((b=2, a=1)) === false

        va = [2, 3, 1]
        vb = [2, 1, 3]
        t = Table(a = va, b = vb)

        @test map(pred, t) == [true, true, false]
        @test findall(pred, t) == [1, 2]
        @test filter(pred, t) == Table(a = [2, 3], b = [2, 1])

        @test Predicate(pred, (a=1,)) === IsLess((b=1,))
        @test Predicate(pred, (b=1,)) === IsGreaterEqual((a=1,))
    end

    @testset "Within" begin
        pred = Within(:a, :b, 1)

        @test pred((a=1, b=-1)) === false
        @test pred((a=1, b=0)) === true
        @test pred((a=1, b=1)) === true
        @test pred((a=1, b=2)) === true
        @test pred((a=1, b=3)) === false
        @test pred((a=1, b=-1, c=2.0)) === false
        @test pred((a=1, b=1, c=2.0)) === true
        @test pred((a=1, b=3, c=2.0)) === false
        @test pred((b=-1, a=1)) === false
        @test pred((b=1, a=1)) === true
        @test pred((b=3, a=1)) === false

        va = [2, 3, 1]
        vb = [2, 1, 3]
        t = Table(a = va, b = vb)

        @test map(pred, t) == [true, false, false]
        @test findall(pred, t) == [1]
        @test filter(pred, t) == Table(a = [2], b = [2])

        @test Predicate(pred, (a=1,)) === In((b=0..2,))
        @test Predicate(pred, (b=1,)) === In((a=0..2,))
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# To be deprecated in 2.0
rank(A::AbstractMatrix, tol::Real) = rank(A,rtol=tol)
nullspace(A::AbstractVector, tol::Real) = nullspace(reshape(A, length(A), 1), rtol= tol)
nullspace(A::AbstractMatrix, tol::Real) = nullspace(A, rtol=tol)
pinv(A::AbstractMatrix{T}, tol::Real) where T = pinv(A, rtol=tol)

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
module ModelingToolkit

using DiffEqBase
import MacroTools: splitdef, combinedef
import IterTools: product

abstract type Expression <: Number end
abstract type AbstractOperation <: Expression end
abstract type AbstractOperator <: Expression end
abstract type AbstractComponent <: Expression end
abstract type AbstractConnection <: Expression end
abstract type AbstractSystem end
abstract type AbstractDomain end

include("domains.jl")
include("variables.jl")

Base.promote_rule(::Type{T},::Type{T2}) where {T<:Number,T2<:Expression} = Expression
Base.one(::Type{T}) where T<:Expression = Constant(1)
Base.zero(::Type{T}) where T<:Expression = Constant(0)

function caclulate_jacobian end

include("operations.jl")
include("operators.jl")
include("systems/diffeqs/diffeqsystem.jl")
include("systems/diffeqs/first_order_transform.jl")
include("systems/nonlinear/nonlinear_system.jl")
include("function_registration.jl")
include("connections.jl")
include("simplify.jl")
include("utils.jl")

export Operation, Expression, AbstractOperator, AbstractComponent, AbstractDomain
export @register
end # module


text:
using AstroBase
import EarthOrientation
using ERFA
using Test

EarthOrientation.update()

@testset "AstroBase" begin
    include("util.jl")
    include("earth_attitude.jl")
    include("bodies.jl")
    include("ephemerides.jl")
    include("two_body.jl")
    # include("n_body.jl")
    include("frames.jl")
    include("coords.jl")
end


text:
"""
    stimmap(residuals, angles)

Calculate the standardized trajectory intensity mean (STIM) map. The inputs are a cube of residuals and the corresponding parallactic angles.

This method seeks to improve upon the typical student-t S/N tests ([`snr`](@ref), [`significance`](@ref)) by calculating statistics in the temporal domain instead of the spatial domain. This is why the full residual cube is required rather than a reduced frame.

In particular, the STIM map is robust to detections with multiple objects or extended sources within the same annuli, which results in very high noise estimates using spatial methods. The STIM map also performs better at small angular separations, since the temporal domain has no limitations from limited resolution elements.

*Pairet et al. 2019* derives a detection threshold of `τ = 0.5` for the STIM map. The detection threshold can be calculated for a specific dataset using [`stim_threshold`](@ref).

# Examples

```julia
julia> cube, angles = # load data

julia> S = subtract(PCA(10), cube, angles);

julia> sm = stimmap(S, angles);
```

# References
1. [Pairet et al. 2019](http://adsabs.harvard.edu/abs/2019MNRAS.487.2262P) "STIM map: detection map for exoplanets imaging beyond asymptotic Gaussian residual speckle noise"

# See Also
[`stim_threshold`](@ref)
"""
function stimmap(residuals::AbstractArray{T,3}, angles) where T
    return stim(derotate(residuals, angles); dims=1)
end

"""
    stim_threshold([stimmap, ] residuals, angles)

Calculate the detection threshold for the standardized trajectory intensity mean (STIM) map. This method uses the same residual cube as [`stimmap`](@ref) but adds an additional step of estimating the residual noise by derotating the residuals with the *opposite* parallactic angles.

If the STIM map has already been calculated, it can be passed in, otherwise it will be calculated in addition to the noise map. Note this will not return the STIM map, only the threshold.

The threshold is derived in section 5.1 of *Pairet et al. 2019* as the ratio of the number of pixels above the approximated noise map. They found a value of `τ = 0.5` to be typical.

# Examples

```julia
julia> cube, angles = # load data

julia> S = subtract(PCA(10), cube, angles);

julia> sm = stimmap(S, angles);

julia> τ = stim_threshold(sm, S, angles);
```

# References
1. [Pairet et al. 2019](http://adsabs.harvard.edu/abs/2019MNRAS.487.2262P) "STIM map: detection map for exoplanets imaging beyond asymptotic Gaussian residual speckle noise"

# See Also
[`stimmap`](@ref)
"""
function stim_threshold(stimmap, residuals, angles)
    # estimate noise map by getting STIM map of cube
    # rotated with opposite angles
    d_opp = Metrics.stim(derotate(residuals, -angles); dims=1)
    # return ratio of values above the noise
    n_ϵ = count(stimmap .> d_opp)
    n = length(stimmap)
    return n_ϵ / n
end

function stim_threshold(residuals, angles)
    # calculate stimmap first
    d = stim(derotate(residuals, angles); dims=1)
    return stim_threshold(d, residuals, angles)
end

"""
    Metrics.stim(cube; dims)

Calculates the STIM map of a derotated cube along the given `dims`. `dims` should correspond to the temporal axis of the cube. The STIM statistic is the slice mean divided by the slice standard deviation. Invalid values will become 0.
"""
function stim(cube; dims)
    μ = mean(cube, dims=dims)
    σ = std(cube, dims=dims, mean=μ)
    d = zero(μ)
    mask = @. !iszero(σ)
    @. d[mask] = μ[mask] / σ[mask]
    return dropdims(d, dims=dims)
end


text:
SetID = Int64

"不确定有穷自动机"
struct FiniteAutomata
    𝐊::StateSet             # 状态集
    𝚺::SymbolTabel          # 字母表
    f::Transition{StateSet} # 状态映射
    𝐒::StateSet             # 开始状态集合
    𝐙::StateSet             # 接受状态集合
end

"可达状态集合"
function move(fa::FiniteAutomata, 𝐈::StateSet, α::Symbol)
    states = StateSet()
    for s in 𝐈
        union!(states, fa.f[s][α])
    end
    states
end

"ϵ-closure ϵ-闭包"
function ϵ_colsure(fa::FiniteAutomata, 𝐈::StateSet)
    states = StateSet(𝐈)
    queue = collect(State, 𝐈)
    while !isempty(queue)
        s = popfirst!(queue)
        for next in get(fa.f[s], :ϵ, StateSet())
            if next ∉ states
                push!(states, next)
                push!(queue, next)
            end
        end
    end
    states
end

"是否为确定的有穷自动机"
function isDeterministic(fa::FiniteAutomata)
    for (_, t) in fa.f
        for (_, s) in t
            "若可达状态数大于`1`,则该状态机为不确定的"
            length(s) > 1 && return false
        end
    end
    true
end

"确定的有穷自动机"
struct DeterministicFiniteAutomata
    𝐊::StateSet             # 状态集合
    𝚺::SymbolTabel          # 字母表
    f::Transition{State}    # 状态转移
    𝐒::State                # 开始状态
    𝐙::StateSet             # 接受状态集合
end

function createSym!(tb::Dict{Symbol,StateSet}, counter::Ref{Int64}, sts::StateSet)
    sym = Symbol("St_$(counter)")
    tb[sym] = sts
    counter += 1
    return sym
end

"子集法产生新的状态集"
function subset(fa::FiniteAutomata)
    counter = 1
    table = Dict{State,StateSet}()
    f = Transition{State}()
    𝐂₀ = ϵ_colsure(fa, fa.𝐊)
    s₀ = createSym!(table, Ref(counter), 𝐂₀)
    𝐂 = Tuple{StateSet,State}[(𝐂₀, s₀),]
    𝐓 = Set{StateSet}()
        while !isempty(𝐂)
        (T, s) = popfirst!(𝐂)
        push!(𝐓, T)
        for α in fa.𝚺
            𝐔 = ϵ_colsure(fa, move(fa, T, α))
            𝐔 ∈ 𝐓 && continue
            sym = createSym!(table, Ref(counter), 𝐔)
            !haskey(f, s) && f[s] = Dict{Symbol,State}()
            f[s][α] = sym
            push!(𝐂, (𝐔, sym))
        end
    end
    table, f
end

"Transform NFA to DFA"
function determine(fa::FiniteAutomata)
    𝐊, f = subset(fa)
    S = Symbol("St_1")
    𝐙 = filter(keys(𝐊)) do s
        any(∈(fa.𝐙), 𝐊[s])
    end
    DeterministicFiniteAutomata(keys(𝐊), fa.𝚺, f, S, 𝐙)
end

function remove_useless!(dfa::DeterministicFiniteAutomata)
    queue = State[dfa.𝐒,]
    𝐊 = StateSet()
    push!(𝐊, dfa.𝐒)
    while !isempty(queue)
        s = popfirst!(queue)
        for α in dfa.𝚺
            !haskey(dfa.f[s], α) && continue
            ns = dfa.f[s][α]
            ns ∈ 𝐊 && continue
            push!(queue, ns)
            push!(𝐊, ns)
        end
    end
    empty!(dfa.𝐊)
    union!(dfa.𝐊, 𝐊)
end

function state2stid(𝐏::Vector{StateSet})
    set_tabel = Dict{State,SetID}()
    for (i, st) in 𝐏
        for s ∈ st
            set_tabel[s] = i
        end
    end
    set_tabel
end

function segmentation!(𝐏::Vector{StateSet}, f::Transition{State}, 𝚺::SymbolTabel)
    set_tabel = state2stid(𝐏)
    i = 1
    while i <= length(𝐏)
        stateset = 𝐏[i]
        length(stateset) == 1 && i += 1 && continue
        cnt = 0
        for α ∈ 𝚺
            new_set = Dict{SetID,StateSet}()
            for state in stateset
                next_state = get(f[state], α, state)
                state_id = set_tabel[next_state]
                !haskey(new_set, state_id) && new_set[state_id] = StateSet()
                push!(new_set[state_id], state)
            end
            if length(new_set) > 1
                for set in values(new_set)
                    splice!(𝐏, i, set)
                    foreach(x -> set_tabel[x] = length(𝐏), set)
                end
                break
            end
            cnt += 1
        end
        cnt == length(𝚺) && i += 1
    end
    𝐏
end

function reset!(f::Transition{State}, src::State, dst::State)
    for k in keys(f)
        if k == src
            delete!(f, k)
        else
            for e in keys(f[k])
                v = f[k][e]
                v == src && f[k][e] = dst
            end
        end
    end
    f
end

"DFA化简（最小化）"
function simplify!(dfa::DeterministicFiniteAutomata)
    remove_useless!(dfa)
    # 分割法
    𝐏 = StateSet[setdiff(dfa.𝐊, dfa.𝐙),dfa.𝐙]
    segmentation!(𝐏, dfa.f, dfa.𝚺)
    for set in 𝐏
        state = pop!(set)
        while !isempty(set)
            removed = pop!(set)
            reset!(dfa.f, removed, state)
        end
    end
    dfa
end

max_stars_repo_name: yeruoforever/SimpleGrammar

max_issues_repo_name: yeruoforever/SimpleGrammar

max_forks_repo_name: yeruoforever/SimpleGrammar


text:
#=
disk-phantom.jl
generate random disk phantoms
2019-07-03, Jeff Fessler, University of Michigan
=#

export disk_phantom_params

using Random: seed!

"""
    params = disk_phantom_params( ; ...)

Generate `ndisk × 6` ellipse phantom parameters
for a head-sized disk plus many disks within it,
designed so that the disks have some minimum separation `minsep`
to avoid overlap and to simplify patch-based model fitting.

# Options
- `fov::Real = 240` image field of view in mm
- `rhead::Function = () -> 100` background radius for "head" [mm]
- `muhead::Function = () -> 1000` "μ" (intensity) value for background head disk
- `rdisk::Function = () -> 10 + 10 * rand()` random disk radii [10,20]
- `mudisk::Function = () -> 100 + 200 * rand()` "μ" values for disks [100,300]
- `ndisk::Function = () -> 10` # of random disks
- `minsep::Real = 8` minimum disk separation in mm
- `maxtry::Int = 500` give up on adding more disks if this is reached
- `warn::Bool = false` warn if maxtry reached?
- `seed::Int = 0` if nonzero then use this seed

The function options can be replaced with rand() for other Distributions.
"""
function disk_phantom_params( ;
    fov::RealU = 240,
#   rhead::RealU = 100, # background radius for "head" [mm]
    rhead::Function = () -> 100, # background radius for "head" [mm]
    muhead::Function = () -> 1000, # "μ" value for background head
#   rmin::RealU = 10, # min radius for disks
#   rmax::RealU = 20, # max radius for disks
    rdisk::Function = () -> 10 + 10 * rand(), # disk radius [10,20]
#   mumin::RealU = 100, # range of "μ" values for disks
#   mumax::RealU = 300,
    mudisk::Function = () -> 100 + 200 * rand(), # "μ" values for disks [100,300]
    ndisk::Function = () -> rand(10:10), # of disks
    minsep::RealU = 8, # minimum separation in mm
    # radial position of disk centers, uses rhead and rmin:
    cdisk::Function = () -> sqrt(rand()) * (100 - 10 - minsep),
    maxtry::Int = 500, # give up on adding more disks if this is reached
    warn::Bool = false, # warn if maxtry reached?
    seed::Int = 0,
)

    ndisk = ndisk()
    rhead = rhead()
    μhead = muhead()

    params = zeros(Float32, ndisk+1, 6)
    params[end,:] = [0, 0, rhead, rhead, 0, μhead]

#   randu = (a, b ; f = x->x) -> a + f(rand()) * (b - a)
    seper = (a,b,r1,x,y,r2) -> sqrt(abs2(x - a) + abs2(y - b)) - r1 - r2

    (seed != 0) && seed!(seed)

    for id=1:ndisk
        trial = 0
        for ii=1:maxtry
#           rc = randu(0, rhead - rmin - minsep, f = x->x^0.5)
            rc = cdisk()
            phi = rand() * 2π
            rad = rdisk()
            val = mudisk()
            trial = [rc*cos(phi), rc*sin(phi), rad, rad, 0, val]

            id == 1 && break # first one is always fine

            # see if trial is too close to another one
            c = 1:(id-1) # check these
            sep = seper.(params[c,1], params[c,2], params[c,3],
                trial[1], trial[2], trial[3])

            (minimum(sep) > minsep) && (rhead - (rc + rad) > minsep) && break
            ii == maxtry && warn && @warn("need more tries")
            ii == maxtry && return params[[1:id; end],:]
        end

#=
        if id > 1
            c = 1:(id-1) # check these
            sep = seper(params[c,1], params[c,2], params[c,3],
                trial[1], trial[2], trial[3])
            @show minimum(sep)
        end
        @show id, trial
=#

        params[id,:] = trial
    end

    return params
end


text:
# This file contains the implementation of different methods for the
# `get_fallback` function. These methods can be used by solver wrappers as
# fallbacks for implementing the `get` method when the solver API does not
# provide the required result. For instance, if the solver does not provide the
# value of the constraints, the solver wrapper can write
# ```julia
# function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
#                  ci::MOI.ConstraintIndex)
#     return MOIU.get_fallback(model, attr, ci)
# end
# ```

"""
    get_fallback(model::MOI.ModelLike, ::MOI.ObjectiveValue)

Compute the objective function value using the `VariablePrimal` results and
the `ObjectiveFunction` value.
"""
function get_fallback(model::MOI.ModelLike, ::MOI.ObjectiveValue)
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    f = MOI.get(model, MOI.ObjectiveFunction{F}())
    # TODO do not include constant if primal solution is a ray
    return evalvariables(vi -> MOI.get(model, MOI.VariablePrimal(), vi), f)
end

"""
    get_fallback(model::MOI.ModelLike, ::MOI.ConstraintPrimal,
                 constraint_index::MOI.ConstraintIndex)

Compute the value of the function of the constraint of index `constraint_index`
using the `VariablePrimal` results and the `ConstraintFunction` values.
"""
function get_fallback(model::MOI.ModelLike, ::MOI.ConstraintPrimal,
                      idx::MOI.ConstraintIndex)
    f = MOI.get(model, MOI.ConstraintFunction(), idx)
    # TODO do not include constant if primal solution is a ray
    return evalvariables(vi -> MOI.get(model, MOI.VariablePrimal(), vi), f)
end

################ Constraint Dual for Variable-wise constraints #################
#
# In the primal we have
#   min a_0' x + b_0
#       A_i  x + b_i in C_i for all i
# In the dual we have
#   max b_0 - sum b_i' y
#       a_0 - sum A_i* y_i = 0
#                      y_i in C_i* for all i
# where A_i* is the adjoint operator of the linear operator A_i. That is, A*
# is the linear operator such that
# ⟨A x, y⟩_{C_i} = ⟨x, A* y⟩_Rn
# where
# * ⟨., .⟩_Rn is the standard scalar product over Rn: ⟨., .⟩_Rn and
# * ⟨., .⟩_{C_i} is the scalar product `set_dot` defined for the set C_i
#
# Suppose we want to get the constraint variable of a variable-wise constraint:
#   A_j x in C_j
# where A_j is zero except on a submatrix which is the identity. We have
# A_j* y_j = a_0 - sum_(i != j) A_i* y_i
# Thus to get the dual y_j, we simply have to compute the right-hand side and
# then invert A_j*. To get the kth element of A_i* y_i we need to compute
# ⟨e_k, A_i* y_i⟩_Rn = ⟨A_i e_k, y_i⟩_{C_i}. A_i e_k is computed using
# `variable_coefficient` and then it is combined with the dual y_i with
# `set_dot`.
# Once A_j* y_j is obtained, we invert A_j* with `dot_coefficients`.

function variable_coefficient(func::MOI.ScalarAffineFunction{T},
                              vi::MOI.VariableIndex) where T
    coef = zero(T)
    for term in func.terms
        if term.variable_index == vi
            coef += term.coefficient
        end
    end
    return coef
end
function variable_coefficient(func::MOI.VectorAffineFunction{T},
                              vi::MOI.VariableIndex) where T
    coef = zeros(T, MOI.output_dimension(func))
    for vector_term in func.terms
        term = vector_term.scalar_term
        if term.variable_index == vi
            coef[vector_term.output_index] += term.coefficient
        end
    end
    return coef
end

"""
    variable_dual(model::MOI.ModelLike,
                  attr::MOI.ConstraintDual,
                  vi::MOI.VariableIndex,
                  ci::MOI.ConstraintIndex{<:Union{MOI.ScalarAffineFunction,
                                                  MOI.VectorAffineFunction})

Return dual of the constraint of index `ci` multiplied by the coefficient of
`vi` in the `MOI.ConstraintFunction`.
"""
function variable_dual(model::MOI.ModelLike,
                       attr::MOI.ConstraintDual,
                       vi::MOI.VariableIndex,
                       ci::MOI.ConstraintIndex{<:MOI.VectorAffineFunction})
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    set = MOI.get(model, MOI.ConstraintSet(), ci)
    coef = variable_coefficient(func, vi)
    dual = MOI.get(model, attr, ci)
    return set_dot(coef, dual, set)
end
function variable_dual(model::MOI.ModelLike,
                       attr::MOI.ConstraintDual,
                       vi::MOI.VariableIndex,
                       ci::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction})
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    coef = variable_coefficient(func, vi)
    dual = MOI.get(model, attr, ci)
    return coef * dual
end

"""
    variable_dual(model::MOI.ModelLike,
                  attr::MOI.ConstraintDual,
                  ci::MOI.ConstraintIndex,
                  vi::MOI.VariableIndex,
                  F::Type{<:MOI.AbstractFunction},
                  S::Type{<:MOI.AbstractSet})

Return sum of the the dual of the `F`-in-`S` constraints except `ci` multiplied
by the coefficient of `vi` in the `MOI.ConstraintFunction`. It errors if another
variable-wise constraint different than `ci` uses `vi`.
"""
function variable_dual(model::MOI.ModelLike,
                       attr::MOI.ConstraintDual,
                       ci::MOI.ConstraintIndex,
                       vi::MOI.VariableIndex,
                       F::Type{<:MOI.AbstractFunction},
                       S::Type{<:MOI.AbstractSet})
    dual = 0.0
    for constraint_index in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        dual += variable_dual(model, attr, vi, constraint_index)
    end
    return dual
end
function variable_dual(model::MOI.ModelLike,
                       ::MOI.ConstraintDual,
                       ci::MOI.ConstraintIndex,
                       vi::MOI.VariableIndex,
                       F::Type{<:Union{MOI.SingleVariable,
                                       MOI.VectorOfVariables}},
                       S::Type{<:MOI.AbstractSet})
    for constraint_index in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        if constraint_index != ci
            func = MOI.get(model, MOI.ConstraintFunction(), constraint_index)
            if (F == MOI.SingleVariable && func.variable == vi) ||
               (F == MOI.VectorOfVariables && vi in func.variables)
               error("Fallback getter for variable constraint dual does not",
                     "support other variable-wise constraints on the variable.",
                     "Please report this issue to the solver wrapper package.")
            end
        end
    end
    return 0.0
end
function variable_dual(::MOI.ModelLike,
                       ::MOI.ConstraintDual,
                       ::MOI.ConstraintIndex,
                       ::MOI.VariableIndex,
                       ::Type{<:Union{MOI.ScalarQuadraticFunction,
                                      MOI.VectorQuadraticFunction}},
                       ::Type{<:MOI.AbstractSet})
    error("Fallback getter for variable constraint dual only supports affine",
          "constraint functions.",
          "Please report this issue to the solver wrapper package.")
end

"""
    variable_dual(model::MOI.ModelLike,
                  attr::MOI.ConstraintDual,
                  ci::MOI.ConstraintIndex,
                  vi::MOI.VariableIndex)

Return the dual of the variable `vi` by using the duals of constraints
of index different than `ci`. It errors if another variable-wise constraint
different than `ci` uses `vi`.
"""
function variable_dual(model::MOI.ModelLike,
                       attr::MOI.ConstraintDual,
                       ci::MOI.ConstraintIndex,
                       vi::MOI.VariableIndex)
    status = MOI.get(model, MOI.DualStatus())
    ray = status == MOI.InfeasibilityCertificate ||
          status == MOI.NearlyInfeasibilityCertificate
    dual = 0.0
    if !ray
        sense = MOI.get(model, MOI.ObjectiveSense())
        # Dual definition for maximization problem corresponds to dual
        # definition for minimization problem with flipped objectived in MOI
        sign = sense == MOI.MaxSense ? -1.0 : 1.0
        F = MOI.get(model, MOI.ObjectiveFunctionType())
        obj_attr = MOI.ObjectiveFunction{F}()
        if F == MOI.SingleVariable
            if MOI.get(model, obj_attr).variable == vi
                dual += sign
            end
        elseif F <: MOI.ScalarAffineFunction
            f = MOI.get(model, obj_attr)
            dual += sign * variable_coefficient(f, vi)
        else
            error("Fallback getter for variable constraint dual only supports",
                  "affine objective function.",
                  "Please report this issue to the solver wrapper package.")
        end
    end
    for FS in MOI.get(model, MOI.ListOfConstraints())
        dual -= variable_dual(model, attr, ci, vi, FS[1], FS[2])
    end
    return dual
end

"""
    variable_dual(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                  ci::MOI.ConstraintIndex{F},
                  func::F) where F <: Union{MOI.SingleVariable,
                                            MOI.VectorOfVariables}

Return the dual of the constraint of index `ci` for which the value of the
`MOI.ConstraintFunction` attribute is `func`.
"""
function variable_dual(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                       ci::MOI.ConstraintIndex{MOI.SingleVariable},
                       func::MOI.SingleVariable)
    return variable_dual(model, attr, ci, func.variable)
end
function variable_dual(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                       ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
                       func::MOI.VectorOfVariables)
    dual = map(vi -> variable_dual(model, attr, ci, vi), func.variables)
    set = MOI.get(model, MOI.ConstraintSet(), ci)
    return dot_coefficients(dual, set)
end

"""
    get_fallback(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{Union{MOI.SingleVariable,
                                               MOI.VectorOfVariables}})

Compute the dual of the constraint of index `ci` using the `ConstraintDual` of
other constraints and the `ConstraintFunction` values. Throws an error if some
constraints are quadratic or if there is one another `MOI.SingleVariable`-in-`S`
or `MOI.VectorOfVariables`-in-`S` constraint with one of the variables in the
function of the constraint `ci`.
"""
function get_fallback(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                      ci::MOI.ConstraintIndex{<:Union{MOI.SingleVariable,
                                                      MOI.VectorOfVariables}})
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return variable_dual(model, attr, ci, func)
end

# Scalar product. Any vector set defined that does not use the standard scalar
# product between vectors of ``R^n`` should redefine `set_dot` and
# `dot_coefficients`.

"""
    set_dot(x::Vector, y::Vector, set::AbstractVectorSet)

Return the scalar product between a vector `x` of the set `set` and a vector
`y` of the dual of the set `s`.
"""
function set_dot(x::Vector, y::Vector, set::MOI.AbstractVectorSet)
    return dot(x, y)
end

function triangle_dot(x::Vector{T}, y::Vector{T}, dim::Int, offset::Int) where T
    result = zero(T)
    k = offset
    for i in 1:dim
        for j in 1:i
            k += 1
            if i == j
                result += x[k] * y[k]
            else
                result += 2 * x[k] * y[k]
            end
        end
    end
    return result
end

function set_dot(x::Vector, y::Vector,
                 set::MOI.PositiveSemidefiniteConeTriangle)
    return triangle_dot(x, y, set.side_dimension, 0)
end

function set_dot(x::Vector, y::Vector, set::Union{MOI.LogDetConeTriangle,
                                                  MOI.RootDetConeTriangle})
    return x[1] * y[1] + triangle_dot(x, y, set.side_dimension, 1)
end

"""
    dot_coefficients(a::Vector, set::AbstractVectorSet)

Return the vector `b` such that for all vector `x` of the set `set`,
`set_dot(b, x, set)` is equal to `dot(a, x)`.
"""
function dot_coefficients(a::Vector, set::MOI.AbstractVectorSet)
    return a
end

function triangle_coefficients!(b::Vector{T}, dim::Int, offset::Int) where T
    k = offset
    for i in 1:dim
        for j in 1:i
            k += 1
            if i != j
                b[k] /= 2
            end
        end
    end
end

function dot_coefficients(a::Vector, set::MOI.PositiveSemidefiniteConeTriangle)
    b = copy(a)
    triangle_coefficients!(b, set.side_dimension, 0)
    return b
end

function dot_coefficients(a::Vector, set::Union{MOI.LogDetConeTriangle,
                                                MOI.RootDetConeTriangle})
    b = copy(a)
    triangle_coefficients!(b, set.side_dimension, 1)
    return b
end


text:
using .NumericalFluxes: EntropyConservative

include("ESDGModel_kernels.jl")

"""
    ESDGModel

Contain type and functor that is used to evaluated the tendency for a entropy
stable / conservative DGSEM discretization. Major fundamental difference between
this and the more vanilla DGSEM is that the first order flux derivatives in the
balance laws are evaluated using "flux-differencing". Namely, the following identities are used:
```math
    ∂x f(q(x)) = 2∂x F(q(x), q(y))|_{x = y},
    A(q(x)) ∂x q(x) = 2∂x D(q(x), q(y))|_{x = y},
```
where the numerical conservative flux `F` and numerical fluctuation flux `D`
satisfy the following consistency and symmetry properties
```math
    F(q, p) = F(p, q),
    F(q, q) = f(q),
    D(q, p) = B(q, p)(q - p),
    2B(q, q) = A(p).
```
For the scheme to be entropy stable (and not just consistent) other properties
of the numerical flux are also required. In particular, consider a balance laws
of the form
```math
    ∂t q + ∑_{j=1:d} (∂xj fj(q) + Aj(q) ∂xj q) = g(q, x, t),
```
where `q` is the state vector, `fj` is the conservative flux, and `Aj`
nonconservative variable coefficient matrix, and `g` is the production field.
Let there exists a scalar companion balance law of the form
```math
    ∂t η(q) + ∑_{j=1:d} ∂xj ζj(q) = Π(q, x, t),
    Π(q, x, t) = β(q)^T g(q, x, t),
    β(q) = ∂q η(q).
```
Then for the scheme to be entropy stable it is requires that the numerical flux
`H(q, p) = F(q, p) + D(q, p)` satisfy the following Tadmor-shuffle:
```math
    β(q)^T Hj(q, p) - β(p)^T Hj(p, q) <= ψj(q) - ψ(p),
    ψj(q) = β(q)^T fj(q) - ζj(q);
```
when the equality is satisfied the scheme is called entropy conservative. For
balance laws without a nonconservative term, `ψj` is the entropy potential.
"""
struct ESDGModel{BL, SA, VNFFO, SNFFO} <: SpaceDiscretization
    "definition of the physics being considered, primary dispatch type"
    balance_law::BL
    "all the grid related information (connectivity, metric terms, etc.)"
    grid::DiscontinuousSpectralElementGrid
    "auxiliary state are quantities needed to evaluate the physics that are not
    explicitly time stepped by the ode solvers"
    state_auxiliary::SA
    "first order, two-point flux to be used for volume derivatives"
    volume_numerical_flux_first_order::VNFFO
    "first order, two-point flux to be used for surface integrals"
    surface_numerical_flux_first_order::SNFFO
end

"""
    ESDGModel(
        balance_law,
        grid;
        state_auxiliary = create_state(balance_law, grid, Auxiliary()),
        volume_numerical_flux_first_order = EntropyConservative(),
        surface_numerical_flux_first_order = EntropyConservative(),
    )

Construct a `ESDGModel` type from a given `grid` and `balance_law` using the
`volume_numerical_flux_first_order` and `surface_numerical_flux_first_order`
two-point fluxes. If the two-point fluxes satisfy the appropriate Tadmor shuffle
then semi-discrete scheme will be entropy stable (or conservative).
"""
function ESDGModel(
    balance_law,
    grid;
    # FIXME: this probably should be done differently
    state_auxiliary = (
        aux = create_state(balance_law, grid, Auxiliary());
        init_state(aux, balance_law, grid, EveryDirection(), Auxiliary())
    ),
    volume_numerical_flux_first_order = EntropyConservative(),
    surface_numerical_flux_first_order = EntropyConservative(),
)

    ESDGModel(
        balance_law,
        grid,
        state_auxiliary,
        volume_numerical_flux_first_order,
        surface_numerical_flux_first_order,
    )
end

"""
    (esdg::ESDGModel)(
        tendency::MPIStateArray,
        state_prognostic::MPIStateArray,
        param::Nothing,
        t,
        α = true,
        β = false,
    )

Compute the entropy stable tendency from the model `esdg`.

    tendency .= α .* dQdt(state_prognostic, param, t) .+ β .* tendency
"""
function (esdg::ESDGModel)(
    tendency::MPIStateArray,
    state_prognostic::MPIStateArray,
    ::Nothing,
    t,
    α = true,
    β = false,
)

    balance_law = esdg.balance_law
    @assert number_states(balance_law, GradientFlux(), Int) == 0

    grid = esdg.grid
    topology = grid.topology

    # Currently only support two polynomial orders
    info = basic_launch_info(esdg)

    device = info.device # array_device(state_prognostic)
    workgroup = (info.Nq[1], info.Nq[2], info.Nqk)
    ndrange = (info.Nq[1] * info.nrealelem, info.Nq[2], info.Nqk)
    nrealelem = length(topology.realelems)

    state_auxiliary = esdg.state_auxiliary

    # XXX: When we do stacked meshes and IMEX this will change
    communicate = true

    exchange_state_prognostic = NoneEvent()

    comp_stream = Event(device)

    ########################
    # tendency Computation #
    ########################
    if communicate
        exchange_state_prognostic = MPIStateArrays.begin_ghost_exchange!(
            state_prognostic;
            dependencies = comp_stream,
        )
    end

    # volume tendency
    comp_stream = esdg_volume_tendency!(device, workgroup)(
        balance_law,
        Val(1),
        Val(info),
        esdg.volume_numerical_flux_first_order,
        tendency.data,
        state_prognostic.data,
        state_auxiliary.data,
        grid.vgeo,
        grid.D[1],
        α,
        β,
        true, # add_source
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )
    comp_stream = esdg_volume_tendency!(device, workgroup)(
        balance_law,
        Val(2),
        Val(info),
        esdg.volume_numerical_flux_first_order,
        tendency.data,
        state_prognostic.data,
        state_auxiliary.data,
        grid.vgeo,
        grid.D[2],
        α,
        true,
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )
    if info.dim == 3
        comp_stream = esdg_volume_tendency!(device, workgroup)(
            balance_law,
            Val(3),
            Val(info),
            esdg.volume_numerical_flux_first_order,
            tendency.data,
            state_prognostic.data,
            state_auxiliary.data,
            grid.vgeo,
            grid.D[3],
            α,
            true,
            ndrange = ndrange,
            dependencies = (comp_stream,),
        )
    end

    # interfaces: Horizontal => Nfp_v and Vertical => Nfp_h
    # mirror surface tendency: interior
    Nfp = info.Nfp_v
    ndrange = Nfp * info.ninteriorelem
    comp_stream = dgsem_interface_tendency!(device, (Nfp,))(
        balance_law,
        Val(info),
        HorizontalDirection(),
        esdg.surface_numerical_flux_first_order,
        nothing,
        tendency.data,
        state_prognostic.data,
        nothing,
        nothing,
        state_auxiliary.data,
        grid.vgeo,
        grid.sgeo,
        t,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        grid.interiorelems,
        α;
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )

    Nfp = info.Nfp_h
    ndrange = Nfp * info.ninteriorelem
    comp_stream = dgsem_interface_tendency!(device, (Nfp,))(
        balance_law,
        Val(info),
        VerticalDirection(),
        esdg.surface_numerical_flux_first_order,
        nothing,
        tendency.data,
        state_prognostic.data,
        nothing,
        nothing,
        state_auxiliary.data,
        grid.vgeo,
        grid.sgeo,
        t,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        grid.interiorelems,
        α;
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )

    if communicate
        exchange_state_prognostic = MPIStateArrays.end_ghost_exchange!(
            state_prognostic;
            dependencies = exchange_state_prognostic,
        )
    end

    # mirror surface tendency: exterior
    Nfp = info.Nfp_v
    ndrange = Nfp * info.nexteriorelem
    comp_stream = dgsem_interface_tendency!(device, (Nfp,))(
        balance_law,
        Val(info),
        HorizontalDirection(),
        esdg.surface_numerical_flux_first_order,
        nothing,
        tendency.data,
        state_prognostic.data,
        nothing,
        nothing,
        state_auxiliary.data,
        grid.vgeo,
        grid.sgeo,
        t,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        grid.exteriorelems,
        α;
        ndrange = ndrange,
        dependencies = (comp_stream, exchange_state_prognostic),
    )

    Nfp = info.Nfp_h
    ndrange = Nfp * info.nexteriorelem
    comp_stream = dgsem_interface_tendency!(device, (Nfp,))(
        balance_law,
        Val(info),
        VerticalDirection(),
        esdg.surface_numerical_flux_first_order,
        nothing,
        tendency.data,
        state_prognostic.data,
        nothing,
        nothing,
        state_auxiliary.data,
        grid.vgeo,
        grid.sgeo,
        t,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        grid.exteriorelems,
        α;
        ndrange = ndrange,
        dependencies = (comp_stream, exchange_state_prognostic),
    )

    # The synchronization here through a device event prevents CuArray based and
    # other default stream kernels from launching before the work scheduled in
    # this function is finished.
    wait(device, comp_stream)
end


text:
@test Tsat(0.10) ≈ 0.372_755_919E3
@test Tsat(0.10u"MPa") ≈ 0.372_755_919E3u"K"

@test Tsat(1.00) ≈ 0.453_035_632E3
@test Tsat(1.00u"MPa") ≈ 0.453_035_632E3u"K"

@test Tsat(10.0) ≈ 0.584_149_488E3
@test Tsat(10.0u"MPa") ≈ 0.584_149_488E3u"K"


text:
# This file includes testset to define new sink types

@testset "NewSinkDefinitionTestset" begin
    # New sink types must be of subtypes of `AbstractSink`.
    @test_throws Exception @eval @def_sink struct Mysink{T,S}
        field1::T
        field2::S
    end

    @test_throws Exception @eval @def_source struct Mysink{T,S} <: SomeDummyType
        readout::RO = t -> sin(t)
        output::OP = Outport()
    end
end


text:
# all functions must receive client and event arguments!
using Discord

function handleping(c::Client, e::MessageCreate)
    if e.message.content == "ping"
        reply(c, e.message, "pong")
    end
end

function handlepong(c::Client, e::MessageCreate)
    if e.message.content == "pong"
        reply(c, e.message, "ping")
    end
end

max_stars_repo_name: dzonimn/HoJDiscordBot

max_issues_repo_name: dzonimn/HoJDiscordBot

max_forks_repo_name: dzonimn/HoJDiscordBot


text:
#bellcurve.jl
#Algorithm 1.1 of Random Eigenvalues by Alan Edelman

#Experiment:  Generate random samples from the normal distribution
#Plot:        Histogram of random samples
#Theory:      The normal distribution curve

## Experiment
t = 1000000
dx = .2
v = randn(t)
x = -4:dx:4
grid,count = hist(v,x)

## Theory
y = exp(-x.^2/2)/sqrt(2*pi)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count/(t*dx), step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(x, y, "color", "blue", "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "bellcurve.png")
end


text:
# example modified from  mthelm85
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie, ForwardDiff
GLMakie.activate!()
let
    f(x,y) = -5*x*y*exp(-x^2-y^2)
    x = -1:0.05:1.0
    y = -1:0.05:1.0
    z = [f(i,j) for i in x, j in y]
    # This is the same function as above, just modified so that it will
    # work with ForwardDiff
    g(x,y) = [-5*x*y*exp(-x^2-y^2)]
    J(xx,yy) = ForwardDiff.jacobian(x -> g(x[1], x[2]), [xx,yy])
    field(i,j) = Point(J(i,j)[1], J(i,j)[2])

    zmin, zmax = minimum(z), maximum(z)
    cmap = :viridis
    function plot()
        fig = Figure(resolution = (1200,600))
        ax1 = Axis3(fig, aspect = (1,1,1), perspectiveness = 0.5,
            elevation = π/3.5, azimuth = 0.1π, )
        ax2 = Axis(fig, aspect = DataAspect(), xlabel = "x", ylabel = "y")
        surface!(ax1, x, y, z, colormap = cmap, colorrange = (zmin, zmax))
        contour3d!(ax1, x, y, z .+ 0.005, levels = 15, linewidth = 2, color=(:white, 0.5))
        wireframe!(ax1, x, y, z, overdraw = false, transparency = true, color = (:black, 0.1))
        streamplot!(ax1, field, -1..1, -1..1, colormap = cmap, gridsize= (40,40),
            arrow_size = 0.05,linewidth=1, transformation = (:xy, -zmax))
        streamplot!(ax2, field, -1..1, -1..1, colormap = cmap, gridsize= (40,40),
                arrow_size = 0.03,linewidth=1)
        fig[1,1] = ax1
        fig[1,2] = ax2
        fig
    end
    fig = with_theme(plot, theme_dark())
    save(joinpath(@__DIR__, "output", "theme_dark_Contour3d.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie","ForwardDiff"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
function mass_matrix_outer_entries!(
    mass_matrix,
    dynamic_device::PSY.DynamicInverter{C, O, IC, DC, P, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}
    @debug "Using default mass matrix entries $O $(PSY.get_name(dynamic_device))"
end

function mdl_outer_ode!(
    device_states,
    output_ode,
    f0,
    ω_sys,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
    )
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]
    Ir_filter = device_states[external_ix[3]]
    Ii_filter = device_states[external_ix[4]]

    #Obtain inner variables for component
    ω_pll = get_inner_vars(dynamic_device)[ω_freq_estimator_var]

    #Get Active Power Controller parameters
    outer_control = PSY.get_outer_control(dynamic_device)
    active_power_control = PSY.get_active_power(outer_control)
    Ta = PSY.get_Ta(active_power_control) #VSM Inertia constant
    kd = PSY.get_kd(active_power_control) #VSM damping constant
    kω = PSY.get_kω(active_power_control) #Frequency droop gain
    ωb = 2 * pi * f0 #Rated angular frequency

    #Get Reactive Power Controller parameters
    reactive_power_control = PSY.get_reactive_power(outer_control)
    kq = PSY.get_kq(reactive_power_control) #Reactive power droop gain
    ωf = PSY.get_ωf(reactive_power_control) #Reactive power filter cutoff frequency

    #Obtain external parameters
    p_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index]
    ω_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][ω_ref_index]
    V_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]
    q_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][Q_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
    )

    #Define internal states for Outer Control
    internal_states = @view device_states[local_ix]
    ω_oc = internal_states[1]
    θ_oc = internal_states[2]
    qm = internal_states[3]

    #Obtain additional expressions
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Compute 3 states ODEs
    output_ode[local_ix[1]] =
        (p_ref / Ta - p_elec_out / Ta - kd * (ω_oc - ω_pll) / Ta - kω * (ω_oc - ω_ref) / Ta)
    output_ode[local_ix[2]] = ωb * (ω_oc - ω_sys)
    output_ode[local_ix[3]] = (ωf * (q_elec_out - qm))

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ_oc
    get_inner_vars(dynamic_device)[ω_oc_var] = ω_oc
    get_inner_vars(dynamic_device)[V_oc_var] = V_ref + kq * (q_ref - qm)
end

function mdl_outer_ode!(
    device_states,
    output_ode,
    f0,
    ω_sys,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.ActivePowerDroop, PSY.ReactivePowerDroop},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerDroop, PSY.ReactivePowerDroop},
    )
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]
    Ir_filter = device_states[external_ix[3]]
    Ii_filter = device_states[external_ix[4]]

    #Obtain inner variables for component
    V_tR = get_inner_vars(dynamic_device)[VR_inv_var]
    V_tI = get_inner_vars(dynamic_device)[VI_inv_var]

    #Get Active Power Controller parameters
    outer_control = PSY.get_outer_control(dynamic_device)
    active_power_control = PSY.get_active_power(outer_control)
    Rp = PSY.get_Rp(active_power_control) #Droop Gain
    ωz = PSY.get_ωz(active_power_control) #Frequency cutoff frequency
    ωb = 2 * pi * f0 #Rated angular frequency

    #Get Reactive Power Controller parameters
    reactive_power_control = PSY.get_reactive_power(outer_control)
    kq = PSY.get_kq(reactive_power_control) #Reactive power droop gain
    ωf = PSY.get_ωf(reactive_power_control) #Reactive power filter cutoff frequency

    #Obtain external parameters
    p_ref = PSY.get_ext(dynamic_device)[PSID.CONTROL_REFS][PSID.P_ref_index]
    ω_ref = PSY.get_ext(dynamic_device)[PSID.CONTROL_REFS][PSID.ω_ref_index]
    V_ref = PSY.get_ext(dynamic_device)[PSID.CONTROL_REFS][PSID.V_ref_index]
    q_ref = PSY.get_ext(dynamic_device)[PSID.CONTROL_REFS][PSID.Q_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerDroop, PSY.ReactivePowerDroop},
    )

    #Define internal states for frequency estimator
    internal_states = @view device_states[local_ix]
    θ_oc = internal_states[1]
    pm = internal_states[2]
    qm = internal_states[3]

    #Obtain additional expressions
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Compute Frequency from Droop
    ω_oc = ω_ref + Rp * (p_ref - pm)

    #Compute 3 states ODEs
    output_ode[local_ix[1]] = ωb * (ω_oc - ω_sys)
    output_ode[local_ix[2]] = (ωz * (p_elec_out - pm))
    output_ode[local_ix[3]] = (ωf * (q_elec_out - qm))

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ_oc
    get_inner_vars(dynamic_device)[ω_oc_var] = ω_oc
    get_inner_vars(dynamic_device)[V_oc_var] = V_ref + kq * (q_ref - qm)
end

function mdl_outer_ode!(
    device_states,
    output_ode,
    f0,
    ω_sys,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.ActivePowerPI, PSY.ReactivePowerPI},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerPI, PSY.ReactivePowerPI},
    )
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]
    Ir_filter = device_states[external_ix[3]]
    Ii_filter = device_states[external_ix[4]]

    #Obtain inner variables for component
    θ_pll = get_inner_vars(dynamic_device)[θ_freq_estimator_var]
    ω_pll = get_inner_vars(dynamic_device)[ω_freq_estimator_var]

    #Get Active Power Controller parameters
    outer_control = PSY.get_outer_control(dynamic_device)
    active_power_control = PSY.get_active_power(outer_control)
    Kp_p = PSY.get_Kp_p(active_power_control) #Proportional Gain
    Ki_p = PSY.get_Ki_p(active_power_control) #Integral Gain
    ωz = PSY.get_ωz(active_power_control) #Frequency cutoff frequency

    #Get Reactive Power Controller parameters
    reactive_power_control = PSY.get_reactive_power(outer_control)
    Kp_q = PSY.get_Kp_q(reactive_power_control) #Proportional Gain
    Ki_q = PSY.get_Ki_q(reactive_power_control) #Integral Gain
    ωf = PSY.get_ωf(reactive_power_control) #Reactive power filter cutoff frequency

    #Obtain external parameters
    p_ref = PSY.get_ext(dynamic_device)[PSID.CONTROL_REFS][PSID.P_ref_index]
    q_ref = PSY.get_ext(dynamic_device)[PSID.CONTROL_REFS][PSID.Q_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerPI, PSY.ReactivePowerPI},
    )

    #Define internal states for frequency estimator
    internal_states = @view device_states[local_ix]
    σp_oc = internal_states[1]
    p_oc = internal_states[2]
    σq_oc = internal_states[3]
    q_oc = internal_states[4]

    #Obtain additional expressions
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Compute 4 states ODEs
    output_ode[local_ix[1]] = p_ref - p_oc
    output_ode[local_ix[2]] = ωz * (p_elec_out - p_oc)
    output_ode[local_ix[3]] = q_ref - q_oc
    output_ode[local_ix[4]] = ωf * (q_elec_out - q_oc)

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ_pll
    get_inner_vars(dynamic_device)[ω_oc_var] = ω_pll
    get_inner_vars(dynamic_device)[Iq_oc_var] = Kp_p * (p_ref - p_oc) + Ki_p * σp_oc
    get_inner_vars(dynamic_device)[Id_oc_var] = Kp_q * (q_ref - q_oc) + Ki_q * σq_oc
end


text: @warn "`neural/stochastic_layers.jl` is not tested."