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,
#                 flux_is_symmetric, volume_flux,
#                 orientation_or_normal_direction, u, equations)
#
# Computes the flux difference ∑_j A[i, j] * f(u_i, u_j) and accumulates the result into `du`.
# Called by `local_flux_differencing` to compute local contributions to flux differencing
# volume integrals.
#
# - `du`, `u` are vectors
# - `A` is the skew-symmetric flux differencing matrix
# - `flux_is_symmetric` is a `Val{<:Bool}` indicating if f(u_i, u_j) = f(u_j, u_i)
#
# The matrix `A` can be either dense or sparse. In the latter case, you should
# use the `adjoint` of a `SparseMatrixCSC` to mimic a `SparseMatrixCSR`, which
# is more efficient for matrix vector products.

# Version for dense operators and symmetric fluxes
@inline function hadamard_sum!(du, A,
                               flux_is_symmetric::Val{true}, volume_flux,
                               orientation_or_normal_direction, u, equations)
  row_ids, col_ids = axes(A)

  for i in row_ids
    u_i = u[i]
    du_i = du[i]
    for j in col_ids
      # 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
        u_j = u[j]
        AF_ij = A[i,j] * volume_flux(u_i, u_j, orientation_or_normal_direction, equations)
        du_i = du_i + AF_ij
        du[j] = du[j] - AF_ij
      end
    end
    du[i] = du_i
  end
end

# Version for dense operators and non-symmetric fluxes
@inline function hadamard_sum!(du, A,
                               flux_is_symmetric::Val{false}, volume_flux,
                               orientation::Integer, u, equations)
  row_ids, col_ids = axes(A)

  for i in row_ids
    u_i = u[i]
    du_i = du[i]
    for j in col_ids
      u_j = u[j]
      f_ij = volume_flux(u_i, u_j, orientation, equations)
      du_i = du_i + A[i,j] * f_ij
    end
    du[i] = du_i
  end
end

# TODO: DGMulti. Fix for curved meshes.
@inline function hadamard_sum!(du, A,
                               flux_is_symmetric::Val{false}, volume_flux,
                               normal_direction::AbstractVector, u, equations)
  row_ids, col_ids = axes(A)

  for i in row_ids
    u_i = u[i]
    du_i = du[i]
    for j in col_ids
      u_j = u[j]
      # The `normal_direction::AbstractVector` has to be passed in twice.
      # This is because on curved meshes, nonconservative fluxes are
      # evaluated using both the normal and its average at interfaces.
      f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations)
      du_i = du_i + A[i,j] * f_ij
    end
    du[i] = du_i
  end
end

# Version for sparse operators and symmetric fluxes
@inline function hadamard_sum!(du, A::LinearAlgebra.Adjoint{<:Any, <:AbstractSparseMatrixCSC},
                               flux_is_symmetric::Val{true}, volume_flux,
                               orientation_or_normal_direction, u, equations)
  A_base = parent(A) # the adjoint of a SparseMatrixCSC is basically a SparseMatrixCSR
  row_ids = axes(A, 2)
  rows = rowvals(A_base)
  vals = nonzeros(A_base)

  for i in row_ids
    u_i = u[i]
    du_i = du[i]
    for id in nzrange(A_base, 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
        u_j = u[j]
        A_ij = vals[id]
        AF_ij = A_ij * volume_flux(u_i, u_j, orientation_or_normal_direction, equations)
        du_i = du_i + AF_ij
        du[j] = du[j] - AF_ij
      end
    end
    du[i] = du_i
  end
end

# TODO: DGMulti. Fix for curved meshes.
# Version for sparse operators and non-symmetric fluxes
@inline function hadamard_sum!(du, A::LinearAlgebra.Adjoint{<:Any, <:AbstractSparseMatrixCSC},
                               flux_is_symmetric::Val{false}, volume_flux,
                               normal_direction::AbstractVector, u, equations)
  A_base = parent(A) # the adjoint of a SparseMatrixCSC is basically a SparseMatrixCSR
  row_ids = axes(A, 2)
  rows = rowvals(A_base)
  vals = nonzeros(A_base)

  for i in row_ids
    u_i = u[i]
    du_i = du[i]
    for id in nzrange(A_base, i)
      j = rows[id]
      # The `normal_direction::AbstractVector` has to be passed in twice.
      # This is because on curved meshes, nonconservative fluxes are
      # evaluated using both the normal and its average at interfaces.
      u_j = u[j]
      f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations)
      du_i = du_i + A[i,j] * f_ij
    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.
@inline function build_lazy_physical_derivative(element, orientation,
                                                mesh::VertexMappedMesh{2}, dg, cache,
                                                operator_scaling = 1.0)
  @unpack Qrst_skew = cache
  @unpack rxJ, sxJ, ryJ, syJ = mesh.md
  scaling = 2 * operator_scaling
  if orientation == 1
    return LazyMatrixLinearCombo(Qrst_skew, scaling .* (rxJ[1,element], sxJ[1,element]))
  else # if orientation == 2
    return LazyMatrixLinearCombo(Qrst_skew, scaling .* (ryJ[1,element], syJ[1,element]))
  end
end

@inline function build_lazy_physical_derivative(element, orientation,
                                                mesh::VertexMappedMesh{3}, dg, cache,
                                                operator_scaling = 1.0)
  @unpack Qrst_skew = cache
  @unpack rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ = mesh.md
  scaling = 2 * operator_scaling
  if orientation == 1
    return LazyMatrixLinearCombo(Qrst_skew, scaling .* (rxJ[1,element], sxJ[1,element], txJ[1,element]))
  elseif orientation == 2
    return LazyMatrixLinearCombo(Qrst_skew, scaling .* (ryJ[1,element], syJ[1,element], tyJ[1,element]))
  else # if orientation == 3
    return LazyMatrixLinearCombo(Qrst_skew, scaling .* (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{1})
  @unpack rxJ = mesh.md
  return 2 * SVector(rxJ[1, element]) # the 1D contravariant vector reduces to a scaling.
end

@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

# use hybridized SBP operators for general flux differencing schemes.
function compute_flux_differencing_SBP_matrices(dg::DGMulti)
  compute_flux_differencing_SBP_matrices(dg, has_sparse_operators(dg))
end

function compute_flux_differencing_SBP_matrices(dg::DGMulti, sparse_operators)
  rd = dg.basis
  Qrst_hybridized, VhP, Ph = StartUpDG.hybridized_SBP_operators(rd)
  Qrst_skew = map(A -> 0.5 * (A - A'), Qrst_hybridized)
  if sparse_operators isa Val{true}
    Qrst_skew = map(Qi -> droptol!(sparse(Qi'), 100 * eps(eltype(Qi)))', Qrst_skew)
  end
  return Qrst_skew, VhP, Ph
end

# use traditional multidimensional SBP operators for SBP approximation types.
function compute_flux_differencing_SBP_matrices(dg::DGMultiFluxDiff{<:SBP}, sparse_operators)
  rd = dg.basis
  @unpack M, Drst, Pq = rd
  Qrst = map(D -> M * D, Drst)
  Qrst_skew = map(A -> 0.5 * (A - A'), Qrst)
  if sparse_operators isa Val{true}
    Qrst_skew = map(Qi -> droptol!(sparse(Qi'), 100 * eps(eltype(Qi)))', Qrst_skew)
  end
  return Qrst_skew
end

# 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)

  # Todo: DGMulti. Factor common storage into a struct (MeshDataCache?) for reuse across solvers?
  # storage for volume quadrature values, face quadrature values, flux values
  nvars = nvariables(equations)
  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,
            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

# most general create_cache: works for `DGMultiFluxDiff{<:Polynomial}`
function create_cache(mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff, RealT, uEltype)

  rd = dg.basis
  @unpack md = mesh

  Qrst_skew, VhP, Ph = compute_flux_differencing_SBP_matrices(dg)

  # temp storage for entropy variables at volume quad points
  nvars = nvariables(equations)
  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,
            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)

  cons2entropy!(entropy_var_values, u_values, equations)

  # "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)

  entropy2cons!(entropy_projected_u_values, projected_entropy_var_values, equations)
  return nothing
end

@inline function cons2entropy!(entropy_var_values::StructArray,
                               u_values          ::StructArray,
                               equations)
  @threaded for i in eachindex(u_values)
    entropy_var_values[i] = cons2entropy(u_values[i], equations)
  end
end

@inline function entropy2cons!(entropy_projected_u_values  ::StructArray,
                               projected_entropy_var_values::StructArray,
                               equations)
  @threaded for i in eachindex(projected_entropy_var_values)
    entropy_projected_u_values[i] = entropy2cons(projected_entropy_var_values[i], equations)
  end
end

# Trait-like system to dispatch based on whether or not the SBP operators are sparse.
# Designed to be extendable to include specialized `approximation_types` too.
@inline function has_sparse_operators(dg::DGMultiFluxDiff)
  rd = dg.basis
  return has_sparse_operators(rd.elementType, rd.approximationType)
end
# For traditional SBP operators on triangles, the operators are fully dense. We avoid using
# sum factorization here, which is slower for fully dense matrices.
@inline has_sparse_operators(::Union{Tri, Tet}, approx_type::AT) where {AT <: SBP} = Val{false}()

# Polynomial-based solvers use hybridized SBP operators, which have blocks scaled by outward
# normal components. This implies that operators for different coordinate directions have
# different sparsity patterns. We default to using sum factorization (which is faster when
# operators are sparse) for all `<:Polynomial` approximation types.
@inline has_sparse_operators(element_type, approx_type::Polynomial) = Val{true}()

# SBP/GaussSBP operators on quads/hexes use tensor-product operators. Thus, sum factorization is
# more efficient and we use the sparsity structure.
@inline has_sparse_operators(::Union{Quad, Hex}, approx_type::AT) where {AT <: SBP} = Val{true}()
@inline has_sparse_operators(::Union{Quad, Hex}, approx_type::GaussSBP) = Val{true}()

# Todo: DGMulti. Dispatch on curved/non-curved mesh types, this code only works for affine meshes (accessing rxJ[1,e],...)
# Computes flux differencing contribution from each Cartesian direction over a single element.
# For dense operators, we do not use sum factorization.
@inline function local_flux_differencing!(fluxdiff_local, u_local, element_index,
                                          has_nonconservative_terms::Val{false}, volume_integral,
                                          has_sparse_operators::Val{false}, mesh,
                                          equations, dg, cache) where {Flux}
  @unpack volume_flux = volume_integral
  for dim in eachdim(mesh)
    Qi_skew = build_lazy_physical_derivative(element_index, dim, mesh, dg, cache)
    # Val{true}() indicates the volume flux is symmetric
    hadamard_sum!(fluxdiff_local, Qi_skew,
                  Val{true}(), volume_flux,
                  dim, u_local, equations)
  end
end

@inline function local_flux_differencing!(fluxdiff_local, u_local, element_index,
                                          has_nonconservative_terms::Val{true}, volume_integral,
                                          has_sparse_operators::Val{false}, mesh,
                                          equations, dg, cache) where {Flux}
  flux_conservative, flux_nonconservative = volume_integral.volume_flux
  for dim in eachdim(mesh)
    Qi_skew = build_lazy_physical_derivative(element_index, dim, mesh, dg, cache)
    # Val{true}() indicates the flux is symmetric.
    hadamard_sum!(fluxdiff_local, Qi_skew,
                  Val{true}(), flux_conservative,
                  dim, u_local, equations)

    # The final argument .5 scales the operator by 1/2 for the nonconservative terms.
    half_Qi_skew = build_lazy_physical_derivative(element_index, dim, mesh, dg, cache, 0.5)
    # Val{false}() indicates the flux is non-symmetric.
    hadamard_sum!(fluxdiff_local, half_Qi_skew,
                  Val{false}(), flux_nonconservative,
                  dim, u_local, equations)
  end
end

# Todo: DGMulti. Dispatch on curved/non-curved mesh types, this code only works for affine meshes (accessing rxJ[1,e],...)
# When the operators are sparse, we use the sum-factorization approach to
# computing flux differencing.
@inline function local_flux_differencing!(fluxdiff_local, u_local, element_index,
                                          has_nonconservative_terms::Val{false}, volume_integral,
                                          has_sparse_operators::Val{true}, mesh,
                                          equations, dg, cache)
  @unpack Qrst_skew = cache
  @unpack volume_flux = volume_integral
  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 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(element_index, dim, mesh)
    Q_skew = Qrst_skew[dim]

    # Val{true}() indicates the flux is symmetric
    hadamard_sum!(fluxdiff_local, Q_skew,
                  Val{true}(), volume_flux,
                  normal_direction, u_local, equations)
  end
end

@inline function local_flux_differencing!(fluxdiff_local, u_local, element_index,
                                          has_nonconservative_terms::Val{true}, volume_integral,
                                          has_sparse_operators::Val{true}, mesh,
                                          equations, dg, cache)
  @unpack Qrst_skew = cache
  flux_conservative, flux_nonconservative = volume_integral.volume_flux
  for dim in eachdim(mesh)
    normal_direction = get_contravariant_vector(element_index, dim, mesh)
    Q_skew = Qrst_skew[dim]

    # Val{true}() indicates the flux is symmetric
    hadamard_sum!(fluxdiff_local, Q_skew,
                  Val{true}(), flux_conservative,
                  normal_direction, u_local, equations)

    # We scale the operator by 1/2 for the nonconservative terms.
    half_Q_skew = LazyMatrixLinearCombo((Q_skew, ), (0.5, ))
    # Val{false}() indicates the flux is non-symmetric
    hadamard_sum!(fluxdiff_local, half_Q_skew,
                  Val{false}(), flux_nonconservative,
                  normal_direction, u_local, equations)
  end
end

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

  @unpack entropy_projected_u_values, Ph = cache
  @unpack fluxdiff_local_threaded, rhs_local_threaded = cache

  @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)

    local_flux_differencing!(fluxdiff_local, u_local, e,
                             have_nonconservative_terms, volume_integral,
                             has_sparse_operators(dg),
                             mesh, equations, dg, cache)

    # 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

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

  @unpack fluxdiff_local_threaded, inv_wq = cache

  @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)

    local_flux_differencing!(fluxdiff_local, u_local, e,
                             have_nonconservative_terms, volume_integral,
                             has_sparse_operators(dg),
                             mesh, equations, dg, cache)

    for i in each_quad_node(mesh, dg, cache)
      du[i, e] = du[i, e] + fluxdiff_local[i] * inv_wq[i]
    end
  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)
  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], SVector(getindex.(md.xyzq, i, e)), t, equations)
    end
  end
end

# Specializes on Polynomial (e.g., modal) DG methods with a flux differencing volume integral, 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{<:Union{Polynomial, GaussSBP}}, cache) where {Source, BC}

  @trixi_timeit timer() "reset ∂u/∂t" 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() "volume integral" calc_volume_integral!(
    du, u, mesh, have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

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

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

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

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

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

  return nothing
end

end # @muladd


text:
edges_match(x::AbstractEdge, y::AbstractEdge) = id(x) == id(y)
@communative edges_match(x::DirectedEdge, y::UndirectedEdge) = false

struct EdgeMismatch <: Exception
    e1::AbstractEdge
    e2::AbstractEdge
end

Base.showerror(io::IO, e::EdgeMismatch) = print(io, e.e1, " ≠ ", e.e2)

function check_edges_match(x, y)
    edges_match(x, y) ? nothing : throw(EdgeMismatch(x, y))
end

@communative function check_edges_match(x::DirectedEdge, y::UndirectedEdge)
    throw(EdgeMismatch(x, y))
end

nodes_match(x::AbstractNode, y::AbstractNode) = id(x) == id(y)
@communative nodes_match(x::ModifiedNode, y::ModifyingNode) = false

struct NodeMismatch <: Exception
    n1::AbstractNode
    n2::AbstractNode
end

Base.showerror(io::IO, e::NodeMismatch) = print(io, e.n1, " ≠ ", e.n2)

function check_nodes_match(x, y)
    nodes_match(x, y) ? nothing : throw(NodeMismatch(x, y))
end

function merge_(x, y)
    (merge(keep(x, Src), keep(y, Src)), merge(keep(x, Dst), keep(y, Dst)))
end

import Base.merge
"""
    merge(x, y)

Merges two (or more) edge/node and fails if not mergable.
In contrast to [`meld`](@ref), it returns always one edge/node or fails.
Works for:

    * ModifiedEdge
    * ModifyingNode
    * ModifiedNode
    * DirectedEdge
    * UndirectedEdge
    * some combinations thereof (those that are mergable)

Requires carefull implementation by subtypes.
Will unlikely work as expected out of the box for other subtypes because it is unclear what is mergable.

# Example

```jldoctest
# `StenoGraphs` does not implement any `EdgeModifier`s
julia> struct Weight{N <: Number} <: EdgeModifier w::N end

julia> merge(Edge(Node(:a), Node(:b)), Edge(Node(:a), Node(:b) * Weight(2)))
a → b * Weight{Int64}(2)

```
"""
function merge(x::DirectedEdge, y::DirectedEdge)
    check_edges_match(x, y)
    DirectedEdge(merge(keep(x, Src), keep(y, Src)), merge(keep(x, Dst), keep(y, Dst)))
end

function merge(x::UndirectedEdge, y::UndirectedEdge)
    check_edges_match(x, y)
    # sorting required
    xns = sort([keep(x, Src), keep(x, Dst)])
    yns = sort([keep(y, Src), keep(y, Dst)])
    UndirectedEdge(merge.(xns, yns)...)
end

function merge(x::Union{UndirectedEdge, DirectedEdge}, y::Union{UndirectedEdge, DirectedEdge})
    # should not merge, check will fail
    check_edges_match(x, y)
end

function merge(x::ModifiedEdge, y::ModifiedEdge)
    check_edges_match(x, y)
    ModifiedEdge(merge(keep(x, Edge), keep(y, Edge)), merge(modifiers(x), modifiers(y)))
end

@communative function merge(x::ModifiedEdge, y::Edge)
    check_edges_match(x, y)
    x
end

function merge(x::T, y::T) where {T <: Node}
    check_nodes_match(x, y)
    x
end

function merge(x::T, y::T) where {T <: ModifiedNode}
    check_nodes_match(x, y)
    ModifiedNode(keep(x, Node), merge(modifiers(x), modifiers(y)))
end

@communative function merge(x::ModifiedNode, y::Node)
    check_nodes_match(x, y)
    x
end

function merge(x::ModifyingNode, y::ModifyingNode)
    check_nodes_match(x, y)
    ModifyingNode(keep(x, Node), merge(modifiers(x), modifiers(y)))
end

@communative function merge(x::ModifyingNode, y::Node)
    check_nodes_match(x, y)
    x
end

@communative function merge(x::ModifyingNode, y::ModifiedNode)
    # should not merge, check will fail
    check_nodes_match(x, y)
end

merge(x::Union{AbstractNode, AbstractEdge}) = x

function merge(x::Union{AbstractEdge, AbstractNode}...)
    foldr(merge, x)
end

"""
    meld(x)

Meldes a vector of edge/node. That means it tries to [`merge`](@ref) elements of the vector that relate to the same nodes.
# Example

```jldoctest
# `StenoGraphs` does not implement any `EdgeModifier`s
struct Weight{N <: Number} <: EdgeModifier w::N end;
e1 = Edge(Node(:a), Node(:b));
e2 = e1 * Weight(1.0);
e3 = UndirectedEdge(Node(:a), Node(:c));
e4 = UndirectedEdge(Node(:c), Node(:a));
es = [e1, e1, e2, e2, e3, e3, e4, e4]
meld(es) == [e2, e3]

# output

true
```

"""
function meld(x)
    ids = id.(x)
    uids = unique(ids)
    out = similar(x, length(uids))
    for (i, id) in enumerate(uids)
        out[i] = merge(x[Ref(id) .== ids]...)
    end
    out
end

ext: jl


text:
function main(input)

  numcols = size(input)[2]
  numrows = size(input)[1]

  output = Int32[]

  for i=1:numrows
    row = input[i,:]
    newrow = mknewmatrixrow(row)

    # in the first loop, output will be empty
    if size(output,1) == 0
      output = newrow
    else
      output = vcat(output,newrow)
    end
  end

  output
end

function mknewmatrixrow(inputrowmatrix)

  numrows = size(inputrowmatrix)[1]
  numcols = size(inputrowmatrix)[2]

  @assert numrows == 1

  outputrow = Int32[]

  largestindex = 1

  for i=1:numcols
    if inputrowmatrix[i] > inputrowmatrix[largestindex]
      largestindex = i
    end
  end

  for i=1:numcols
    if i == largestindex
      push!(outputrow,1)
    else
      push!(outputrow,0)
    end
  end

  outputrowmatrix = reshape(outputrow,1,numcols)
  @assert size(outputrowmatrix)[2] == numcols
  outputrowmatrix
end


text:
module NonlinearOperatorsTests

using Test
using Gridap.Algebra
using Gridap.Algebra: NonlinearOperatorMock

op = NonlinearOperatorMock()

x = [1.0, 3.0]
b = [0.0,0.0]
A = [-1.0 0.0; 0.0 1.0]

test_nonlinear_operator(op,x,b,jac=A)

@test residual(op,x) ≈ b
@test jacobian(op,x) ≈ A
b1, A1 = residual_and_jacobian(op,x)
@test b1 ≈ b
@test A1 ≈ A

end # module


text:
# Testing a pretrained normalizing flow for seismic imaging
# Authors: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
"""
Example:
julia test_hint_imaging.jl
"""

using DrWatson
@quickactivate :ReliabilityAwareImaging

using InvertibleNetworks
using HDF5
using JLD2
using Random
using Statistics
using ArgParse
using ProgressMeter
using PyPlot
using Seaborn
using LinearAlgebra
using Flux: gpu, cpu
using Flux.Data: DataLoader

set_style("whitegrid")
rc("font", family="serif", size=12)
font_prop = matplotlib.font_manager.FontProperties(
    family="serif",
    style="normal",
    size=13
)
sfmt=matplotlib.ticker.ScalarFormatter(useMathText=true)
sfmt.set_powerlimits((0, 0))
matplotlib.use("Agg")

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 100
    "--epoch"
        help = "Epoch to load net params"
        arg_type = Int
        default = -1
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-4
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 2
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 16
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 128
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 12
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "seismic-imaging-posterior-learning"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
sim_name = parsed_args["sim_name"]

if parsed_args["epoch"] == -1
    parsed_args["epoch"] = parsed_args["max_epoch"]
end
epoch = parsed_args["epoch"]

# Define raw data directory
mkpath(datadir("training-data"))
data_path = datadir("training-data", "training-pairs.h5")

# Download the dataset into the data directory if it does not exist
if isfile(data_path) == false
    run(`wget https://www.dropbox.com/s/53u8ckb9aje8xv4/'
        'training-pairs.h5 -q -O $data_path`)
end

# Load seismic images and create training and testing data
file = h5open(data_path, "r")
X_train = file["dm"][:, :, :, :]
Y_train = file["rtm"][:, :, :, :]

nx, ny, nc, nsamples = size(X_train)

AN_x = ActNorm(nsamples)
AN_y = ActNorm(nsamples)
X_train = AN_x.forward(X_train)
Y_train = AN_y.forward(Y_train)

# Dimensions after wavelet squeeze to increase no. of channels
nx = Int(nx/2)
ny = Int(ny/2)
n_in = Int(nc*4)

# Create network
CH = NetworkConditionalHINT(
    n_in, n_hidden, depth, k1=3, k2=3, p1=1, p2=1
)

# Loading the experiment—only network weights and training loss
Params, fval, exp_path = load_experiment(parsed_args; return_path=true)
train_idx = wload(exp_path)["train_idx"]
put_params!(CH, convert(Array{Any,1}, Params))

# test data pairs
idx = shuffle(setdiff(1:nsamples, train_idx))[1]
X_fixed = wavelet_squeeze(X_train[:, :, :, idx:idx])
Y_fixed = wavelet_squeeze(Y_train[:, :, :, idx:idx])

# Now select single fixed sample from all Ys
Zy_fixed = CH.forward_Y(Y_fixed)

# Draw new Zx, while keeping Zy fixed
n_samples = 1000
X_post = zeros(Float32, nx, ny, n_in, n_samples)
CH = CH |> gpu

test_batchsize = 16
test_loader = DataLoader(
    (randn(Float32, nx, ny, n_in, n_samples), repeat(Zy_fixed, 1, 1, 1, n_samples)),
    batchsize=test_batchsize, shuffle=false
)

p = Progress(length(test_loader))

for (itr, (X, Y)) in enumerate(test_loader)
    Base.flush(Base.stdout)
    counter = (itr - 1)*test_batchsize + 1

    X = X |> gpu
    Y = Y |> gpu

    X_post[:, :, :, counter:(counter+size(X)[4]-1)] = (CH.inverse(X, Y)[1] |> cpu)
    ProgressMeter.next!(p)
end

X_post = wavelet_unsqueeze(X_post)
X_post[1:10, :, :, :] .= 0f0
X_post = AN_x.inverse(X_post)

# Some stats
X_post_mean = mean(X_post; dims=4)
X_post_std = std(X_post; dims=4)

X_fixed = wavelet_unsqueeze(X_fixed)
Y_fixed = wavelet_unsqueeze(Y_fixed)
Zy_fixed = wavelet_unsqueeze(Zy_fixed)

X_fixed = AN_x.inverse(X_fixed)
Y_fixed = AN_y.inverse(Y_fixed)
Y_fixed[1:10, :, :, :] .= 0f0

save_dict = @strdict max_epoch epoch lr lr_step batchsize n_hidden depth sim_name
save_path = plotsdir(sim_name, savename(save_dict; digits=6))

spacing = [20.0, 12.5]
extent = [0., size(X_fixed, 1)*spacing[1], size(X_fixed, 2)*spacing[2], 0.]/1e3

# Training loss
fig = figure("training logs", figsize=(7, 2.5))
if epoch == parsed_args["max_epoch"]
    plot(range(0, epoch, length=length(fval)), fval, color="#d48955")
else
    plot(
        range(0, epoch, length=length(fval[1:findfirst(fval .== 0f0)-1])),
        fval[1:findfirst(fval .== 0f0)-1], color="#d48955"
    )
end
title("Negative log-likelihood")
ylabel("Training objective")
xlabel("Epochs")
wsave(joinpath(save_path, "log.png"), fig)
close(fig)

# Plot the true model
fig = figure("x", figsize=(7.68, 4.8))
imshow(
    X_fixed[:, :, 1, 1], vmin=-1.5e3, vmax=1.5e3, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1, extent=extent
)
title(L"High-fidelity image, $\mathbf{x}$")
colorbar(fraction=0.03, pad=0.01, format=sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx),  "true_model.png"), fig)
close(fig)

# Plot the observed data
fig = figure("y", figsize=(7.68, 4.8))
imshow(
    Y_fixed[:, :, 1, 1], vmin=-1.5e6, vmax=1.5e6, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1, extent=extent
)
title(L"Low-fidelity reverse-time migrated image, $\mathbf{y}$")
colorbar(fraction=0.03, pad=0.01, format=sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "observed_data.png"), fig)
close(fig)

# Plot the conditional mean estimate
fig = figure("x_cm", figsize=(7.68, 4.8))
imshow(
    X_post_mean[:, :, 1, 1], vmin=-1.5e3, vmax=1.5e3, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1, extent=extent
)
title(L"Conditional mean, $\mu$")
colorbar(fraction=0.03, pad=0.01, format=sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "conditional_mean.png"), fig)
close(fig)

# Plot the pointwise standard deviation
fig = figure("x_std", figsize=(7.68, 4.8))
imshow(
    X_post_std[:, :, 1, 1], vmin=60,
    vmax=320.0, aspect=1, cmap="OrRd", resample=true,
    interpolation="lanczos", filterrad=1, extent=extent,
    norm=matplotlib.colors.LogNorm()
)
title(L"Pointwise standard deviation, $\sigma$")
cp = colorbar(fraction=0.03, pad=0.01)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "pointwise_std.png"), fig)
close(fig)

for ns= 1:10
    fig = figure("x_cm", figsize=(7.68, 4.8))
    imshow(
        X_post[:, :, 1, ns], vmin=-1.5e3, vmax=1.5e3, aspect=1, cmap="Greys",
        resample=true, interpolation="lanczos", filterrad=1, extent=extent
    )
    title(L"Posterior sample, $\mathbf{x} \sim p(\mathbf{x} \mid \mathbf{y})$")
    colorbar(fraction=0.03, pad=0.01, format=sfmt)
    grid(false)
    xlabel("Horizontal distance (km)")
    ylabel("Depth (km)")
    safesave(joinpath(save_path, string(idx), "sample.png"), fig)
    close(fig)
end


text:
module ArrowTypes

export ArrowType, NullType, PrimitiveType, ListType, FixedSizeListType, MapType, StructType, UnionType, DictEncodedType

abstract type ArrowType end

ArrowType(x::T) where {T} = ArrowType(T)
ArrowType(::Type{T}) where {T} = isprimitivetype(T) ? PrimitiveType() : StructType()

function arrowconvert end

arrowconvert(T, x) = convert(T, x)
arrowconvert(::Type{Union{T, Missing}}, x) where {T} = arrowconvert(T, x)
arrowconvert(::Type{Union{T, Missing}}, ::Missing) where {T} = missing

struct NullType <: ArrowType end

ArrowType(::Type{Missing}) = NullType()

struct PrimitiveType <: ArrowType end

ArrowType(::Type{<:Integer}) = PrimitiveType()
ArrowType(::Type{<:AbstractFloat}) = PrimitiveType()
ArrowType(::Type{Bool}) = PrimitiveType()

struct ListType <: ArrowType end

isstringtype(T) = false
isstringtype(::Type{Union{T, Missing}}) where {T} = isstringtype(T)

ArrowType(::Type{<:AbstractString}) = ListType()
isstringtype(::Type{<:AbstractString}) = true

ArrowType(::Type{Symbol}) = ListType()
isstringtype(::Type{Symbol}) = true
arrowconvert(::Type{Symbol}, x::String) = Symbol(x)
arrowconvert(::Type{String}, x::Symbol) = String(x)

ArrowType(::Type{<:AbstractArray}) = ListType()

struct FixedSizeListType <: ArrowType end

ArrowType(::Type{NTuple{N, T}}) where {N, T} = FixedSizeListType()
getsize(::Type{NTuple{N, T}}) where {N, T} = N

struct MapType <: ArrowType end

ArrowType(::Type{<:Dict}) = MapType()

struct StructType <: ArrowType end

ArrowType(::Type{<:NamedTuple}) = StructType()

@enum STRUCT_TYPES NAMEDTUPLE STRUCT # KEYWORDARGS

structtype(::Type{NamedTuple{N, T}}) where {N, T} = NAMEDTUPLE
structtype(::Type{T}) where {T} = STRUCT

struct UnionType <: ArrowType end

ArrowType(::Union) = UnionType()

struct DictEncodedType <: ArrowType end

"""
There are a couple places when writing arrow buffers where
we need to write a "dummy" value; it doesn't really matter
what we write, but we need to write something of a specific
type. So each supported writing type needs to define `default`.
"""
function default end

default(T) = zero(T)
default(::Type{Symbol}) = Symbol()
default(::Type{Char}) = '\0'
default(::Type{String}) = ""

function default(::Type{A}) where {A <: AbstractVector{T}} where {T}
    a = similar(A, 1)
    a[1] = default(T)
    return a
end

default(::Type{NTuple{N, T}}) where {N, T} = ntuple(i -> default(T), N)
default(::Type{T}) where {T <: Tuple} = Tuple(default(fieldtype(T, i)) for i = 1:fieldcount(T))
default(::Type{Dict{K, V}}) where {K, V} = Dict{K, V}()
default(::Type{NamedTuple{names, types}}) where {names, types} = NamedTuple{names}(Tuple(default(fieldtype(types, i)) for i = 1:length(names)))

const JULIA_TO_ARROW_TYPE_MAPPING = Dict{Type, Tuple{String, Type}}(
    Char => ("JuliaLang.Char", UInt32),
    Symbol => ("JuliaLang.Symbol", String),
)

istyperegistered(::Type{T}) where {T} = haskey(JULIA_TO_ARROW_TYPE_MAPPING, T)

function getarrowtype!(meta, ::Type{T}) where {T}
    arrowname, arrowtype = JULIA_TO_ARROW_TYPE_MAPPING[T]
    meta["ARROW:extension:name"] = arrowname
    meta["ARROW:extension:metadata"] = ""
    return arrowtype
end

const ARROW_TO_JULIA_TYPE_MAPPING = Dict{String, Tuple{Type, Type}}(
    "JuliaLang.Char" => (Char, UInt32),
    "JuliaLang.Symbol" => (Symbol, String),
)

function extensiontype(meta)
    if haskey(meta, "ARROW:extension:name")
        typename = meta["ARROW:extension:name"]
        if haskey(ARROW_TO_JULIA_TYPE_MAPPING, typename)
            return ARROW_TO_JULIA_TYPE_MAPPING[typename][1]
        else
            @warn "unsupported ARROW:extension:name type: \"$typename\""
        end
    end
    return nothing
end

function registertype!(juliatype::Type, arrowtype::Type, arrowname::String=string("JuliaLang.", string(juliatype)))
    # TODO: validate that juliatype isn't already default arrow type
    JULIA_TO_ARROW_TYPE_MAPPING[juliatype] = (arrowname, arrowtype)
    ARROW_TO_JULIA_TYPE_MAPPING[arrowname] = (juliatype, arrowtype)
    return
end

end # module ArrowTypes

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
using AdventOfCode2021
using Test

@testset "Day 1" begin
    @test AdventOfCode2021.Day1.num_increases([199,200,208,210,200,207,240,269,260,263]) == 7
    @test AdventOfCode2021.Day1.part1_answer() == 1564
    @test AdventOfCode2021.Day1.num_3_window_increases([199,200,208,210,200,207,240,269,260,263]) == 5
    @test AdventOfCode2021.Day1.part2_answer() == 1611
end

@testset "Day 2" begin
    course = [
        "forward" 5;
        "down" 5;
        "forward" 8;
        "up" 3;
        "down" 8;
        "forward" 2
    ]
    @test AdventOfCode2021.Day2.calculate_position_and_depth(course) == (15,10)
    @test AdventOfCode2021.Day2.part1_answer() == 1962940
    @test AdventOfCode2021.Day2.calculate_position_and_depth_with_aim(course) == (15, 60)
    @test AdventOfCode2021.Day2.part2_answer() == 1813664422
end

@testset "Day 3" begin

    diagnostics = [
        0 0 1 0 0
        1 1 1 1 0
        1 0 1 1 0
        1 0 1 1 1
        1 0 1 0 1
        0 1 1 1 1
        0 0 1 1 1
        1 1 1 0 0
        1 0 0 0 0
        1 1 0 0 1
        0 0 0 1 0
        0 1 0 1 0
    ]
    @test AdventOfCode2021.Day3.calculate_gamma_and_epsilon_rates(diagnostics) == (22, 9)
    @test AdventOfCode2021.Day3.part1_answer() == 2583164
    @test AdventOfCode2021.Day3.calculate_oxygen_and_co2_ratings(diagnostics) == (23, 10)
    @test AdventOfCode2021.Day3.part2_answer() == 2784375
end

@testset "Day 4" begin
    input = [
        "7,4,9,5,11,17,23,2,0,14,21,24,10,16,13,6,15,25,12,22,18,20,8,19,3,26,1"
        ""
        "22 13 17 11  0"
        " 8  2 23  4 24"
        "21  9 14 16  7"
        " 6 10  3 18  5"
        " 1 12 20 15 19"
        ""
        " 3 15  0  2 22"
        " 9 18 13 17  5"
        "19  8  7 25 23"
        "20 11 10 24  4"
        "14 21 16 12  6"
        ""
        "14 21 17 24  4"
        "10 16 15  9 19"
        "18  8 23 26 20"
        "22 11 13  6  5"
        " 2  0 12  3  7"
    ]
end


text:
# This file contains tests of internal mechanisms that are probably
# not of user-level interest.  One consequence is that it's OK to
# replace/break these tests if you're redesigning the internal API
# (as long as you preserve the same overall behavior).

A = rand(3,4)
AI.checksame_inds(A, A)
@test_throws DimensionMismatch AI.checksame_inds(A, A')
@test isa(@inferred(AI.storageorder(A)), AI.FirstToLast)
R = reshape(A, (2,3,2))
@test isa(@inferred(AI.storageorder(R)), AI.FirstToLast)
S = view(A, 1:3, 1:4)
R = reshape(S, (2,3,2))
@test isa(@inferred(AI.storageorder(R)), AI.FirstToLast)
B = PermutedDimsArray(A, [2,1])
@test isa(@inferred(AI.storageorder(B)), AI.OtherOrder{(2,1)})
R = reshape(B, (2,2,3))
@test isa(@inferred(AI.storageorder(R)), AI.NoOrder)

@test @inferred(AI.contiguous_index((:, :))) == AI.Contiguous()
@test @inferred(AI.contiguous_index((:, 3))) == AI.Contiguous()
@test @inferred(AI.contiguous_index((3, :))) == AI.NonContiguous()
@test @inferred(AI.contiguous_index((3, 3))) == AI.Contiguous()
@test @inferred(AI.contiguous_index((:, :, 3))  ) == AI.Contiguous()
@test @inferred(AI.contiguous_index((:, :, 1:2))) == AI.Contiguous()
@test @inferred(AI.contiguous_index((:, 3, :))  ) == AI.NonContiguous()
@test @inferred(AI.contiguous_index((:, 3, 1:2))) == AI.NonContiguous()

@test isless(CartesianIndex((1,1)), CartesianIndex((2,1)))
@test isless(CartesianIndex((1,1)), CartesianIndex((1,2)))
@test isless(CartesianIndex((2,1)), CartesianIndex((1,2)))
@test !isless(CartesianIndex((1,2)), CartesianIndex((2,1)))

A = sparse([1,4,3],[1,1,2],[0.2,0.4,0.6])
@test isa(A, AI.ContiguousCSC)
@test isa(view(A, :, 1), AI.ContiguousCSC)
@test isa(view(A, 1:2, 2), AI.ContiguousCSC)
@test isa(view(A, :, 1:2), AI.ContiguousCSC)
@test !isa(view(A, 1:2, 1:2), AI.ContiguousCSC)
@test !isa(view(A, :, [1,2]), AI.ContiguousCSC)
@test !isa(view(A, [1,2], 1), AI.ContiguousCSC)


text:
module ThreeJSTest

using ThreeJS
using FactCheck
using Colors
using Compat

import Patchwork: Elem
# write your own tests here
facts("Testing General Functions") do
    context("Testing outerdiv") do
        @fact outerdiv() -->
            Elem(
                :div,
                style=Dict(:width=>"100%", :height=>"600px")
            )
        @fact outerdiv("90%","200px") -->
            Elem(
                :div,
                style=Dict(:width=>"90%", :height=>"200px")
            )
    end
    context("Testing initscene") do
        @fact initscene() --> Elem(
            :"three-js",
            attributes = Dict(:bgcolor => "#FFFFFF")
        )
        @fact initscene(bgcolor=colorant"black") --> Elem(
            :"three-js",
            attributes = Dict(:bgcolor => "#000000")
        )
    end
end

facts("Testing Render Elem Outputs") do
    context("Testing mesh") do
        @fact mesh(10.0, 10.0, 10.0) -->
            Elem(
                :"three-js-mesh",
                attributes = Dict(
                    :x => 10.0,
                    :y => 10.0,
                    :z => 10.0,
                    :rx => 0.0,
                    :ry => 0.0,
                    :rz => 0.0
                )
            )
        @fact mesh(10.0, 10.0, 10.0, rx = 20.0, ry = 15.0, rz = 240.0) -->
             Elem(
                :"three-js-mesh",
                attributes = Dict(
                    :x => 10.0,
                    :y => 10.0,
                    :z => 10.0,
                    :rx => 20.0,
                    :ry => 15.0,
                    :rz => 240.0
                )
            )
    end
    context("Testing geometry tags") do
        @fact box(10.0, 11.0, 12.0) -->
            Elem(
                :"three-js-box",
                attributes = Dict(:width => 10.0, :height => 11.0, :depth => 12.0)
            )
        @fact sphere(10.0) -->
            Elem(:"three-js-sphere", attributes = Dict(:r => 10.0))
        @fact pyramid(10.0, 12.0) -->
            Elem(
                :"three-js-pyramid",
                attributes = Dict(:base => 10.0, :height => 12.0)
            )
        @fact cylinder(10.0, 11.0, 12.0) -->
            Elem(
                :"three-js-cylinder",
                attributes = Dict(:top => 10.0, :bottom => 11.0, :height => 12.0)
            )
        @fact torus(12.0, 2.0) -->
            Elem(:"three-js-torus", attributes = Dict(:r => 12.0, :tube => 2.0))
        @fact plane(12.0, 2.0) -->
            Elem(:"three-js-plane", attributes = Dict(:w => 12.0, :h => 2.0))
        colormap = Colors.colormap("RdBu")
        ys = [
                0.0 1.0 2.0 3.0
                1.0 2.0 3.0 4.0
                2.0 3.0 4.0 5.0
             ]
        hexcolors = map(z -> "#"*hex(colormap[ceil(Int,(5-z)/5 * (100-1)+1)]), ys)
        colors = Array{Color}(3,4)
        map!(z -> colormap[ceil(Int,(5-z)/5 * (100-1)+1)], colors, ys)
        @fact parametric(2, 3, 0:2, 0:3, (x, y) -> x + y) -->
            Elem(
                :"three-js-parametric",
                attributes = Dict(
                    :slices => 2,
                    :stacks => 3,
                    :x => [
                        0.0 0.0 0.0 0.0
                        1.0 1.0 1.0 1.0
                        2.0 2.0 2.0 2.0
                    ],
                    :y => ys,
                    :z => [
                        0.0 1.0 2.0 3.0
                        0.0 1.0 2.0 3.0
                        0.0 1.0 2.0 3.0
                    ],
                    :vertexcolors => hexcolors
                    )
            )
        @fact meshlines(2, 3, 0:2, 0:3, (x, y) -> x + y) -->
            map(x->x << linematerial(Dict(:color=>"white", :colorkind=>"vertex")),
                [
                line(
                    [(0.0, 0.0, 0.0), (1.0, 1.0, 0.0), (2.0, 2.0, 0.0)],
                    vertexcolors= colors[:, 1]
                ),
                line(
                    [(0.0, 1.0, 1.0), (1.0, 2.0, 1.0), (2.0, 3.0, 1.0)],
                    vertexcolors= colors[:, 2]
                ),
                line(
                    [(0.0, 2.0, 2.0), (1.0, 3.0, 2.0), (2.0, 4.0, 2.0)],
                    vertexcolors= colors[:, 3]
                ),
                line(
                    [(0.0, 3.0, 3.0), (1.0, 4.0, 3.0), (2.0, 5.0, 3.0)],
                    vertexcolors= colors[:, 4]
                ),
                line(
                    [(0.0, 0.0, 0.0), (0.0, 1.0, 1.0), (0.0, 2.0, 2.0), (0.0, 3.0, 3.0)],
                    vertexcolors= reshape(transpose(colors[1, :]), 4)),
                line(
                    [(1.0, 1.0, 0.0), (1.0, 2.0, 1.0), (1.0, 3.0, 2.0), (1.0, 4.0, 3.0)],
                    vertexcolors= reshape(transpose(colors[2, :]), 4)),
                line(
                    [(2.0, 2.0, 0.0), (2.0, 3.0, 1.0), (2.0, 4.0, 2.0), (2.0, 5.0, 3.0)],
                    vertexcolors= reshape(transpose(colors[3, :]), 4))
                ]
            )
        @fact dodecahedron(4.0) -->
            Elem(:"three-js-dodecahedron", attributes = Dict(:r => 4.0))
        @fact icosahedron(4.0) -->
            Elem(:"three-js-icosahedron", attributes = Dict(:r => 4.0))
        @fact octahedron(4.0) -->
            Elem(:"three-js-octahedron", attributes = Dict(:r => 4.0))
        @fact tetrahedron(4.0) -->
            Elem(:"three-js-tetrahedron", attributes = Dict(:r => 4.0))
        verts = [(i, i, i) for i = 1.0:3.0]
    end
    context("Testing geometry") do
        verts = [(i, i, i) for i = 1.0:3.0]
        faces = [(1 , 2, 3)]
        @fact geometry(verts, faces) -->
            Elem(
                :"three-js-geometry",
                attributes = Dict(
                    :x => [1.0, 2.0, 3.0],
                    :y => [1.0, 2.0, 3.0],
                    :z => [1.0, 2.0, 3.0],
                    :faces => [0, 1, 2]
                )
            )
    end
    context("Testing light tags") do
        @fact pointlight(10.0, 11.0, 12.0) -->
            Elem(
                :"three-js-light",
                attributes = Dict(
                :x => 10.0, :y => 11.0 , :z => 12.0,
                :kind => "point",
                :color => "#FFFFFF",
                :intensity => 1.0,
                :distance => 0.0
                )
            )
        @fact pointlight(
            10.0, 11.0, 12.0,
            intensity = 2.0, distance = 1000.0, color = colorant"red"
            ) --> Elem(
                :"three-js-light",
                attributes = Dict(
                    :x => 10.0, :y => 11.0 , :z => 12.0,
                    :kind => "point",
                    :color => "#FF0000",
                    :intensity => 2.0,
                    :distance => 1000.0
                )
        )
        @fact spotlight(10.0, 11.0, 12.0) -->
            Elem(
                :"three-js-light",
                attributes = Dict(
                    :x => 10.0, :y => 11.0 , :z => 12.0,
                    :kind => "spot",
                    :color => "#FFFFFF",
                    :intensity => 1.0,
                    :distance => 0.0,
                    :angle => 60.0,
                    :exponent => 8.0,
                    :shadow => false
                )
            )
        @fact spotlight(
            10.0, 11.0, 12.0,
            intensity = 2.0, distance = 1000.0, color = colorant"red",
            angle = 45.0, exponent = 10.0, shadow = true
            ) --> Elem(
                :"three-js-light",
                attributes = Dict(
                    :x => 10.0, :y => 11.0 , :z => 12.0,
                    :kind => "spot",
                    :color => "#FF0000",
                    :intensity => 2.0,
                    :distance => 1000.0,
                    :angle => 45.0,
                    :exponent => 10.0,
                    :shadow => true
                )
            )
        @fact ambientlight() -->
            Elem(
                :"three-js-light",
                attributes = Dict(:kind => "ambient", :color => "#FFFFFF")
            )
        @fact ambientlight(colorant"red") -->
            Elem(
                :"three-js-light",
                attributes = Dict(:kind => "ambient", :color => "#FF0000")
            )
    end
    context("Testing camera tag") do
        @fact camera(10.0, 11.0, 12.0) -->
            Elem(
                :"three-js-camera",
                attributes = Dict(
                    :x => 10.0,
                    :y => 11.0,
                    :z => 12.0,
                    :fov => 45.0,
                    :aspect => nothing,
                    :near => 1.0,
                    :far => 1000.0
                )
            )
    end
    context("Testing line") do
        verts = [Tuple{Float64, Float64, Float64}((i, i, i)) for i = 1.0:3.0]
        colors = Color[colorant"red", colorant"blue", colorant"green"]
        vertswithcolors = [Tuple{Float64, Float64, Float64, Color}(
            (verts[i][1], verts[i][2], verts[i][3], colors[i])
        ) for i=1:3]
        @fact line(
        verts, x = 10.0, y = 10.0, z = 10.0, rx = 20.0, ry = 15.0, rz = 240.0,
        kind = "pieces", vertexcolors = colors) -->
             Elem(
                :"three-js-line",
                attributes = Dict(
                    :kind => "pieces",
                    :x => 10.0,
                    :y => 10.0,
                    :z => 10.0,
                    :rx => 20.0,
                    :ry => 15.0,
                    :rz => 240.0,
                    :xs => [1.0, 2.0, 3.0],
                    :ys => [1.0, 2.0, 3.0],
                    :zs => [1.0, 2.0, 3.0],
                    :vertexcolors => map(x->"#"*hex(x),colors)
                )
            )
        @fact line(
        vertswithcolors, x = 10.0, y = 10.0, z = 10.0, rx = 20.0, ry = 15.0,
        rz = 240.0, kind = "pieces") -->
             Elem(
                :"three-js-line",
                attributes = Dict(
                    :kind => "pieces",
                    :x => 10.0,
                    :y => 10.0,
                    :z => 10.0,
                    :rx => 20.0,
                    :ry => 15.0,
                    :rz => 240.0,
                    :xs => [1.0, 2.0, 3.0],
                    :ys => [1.0, 2.0, 3.0],
                    :zs => [1.0, 2.0, 3.0],
                    :vertexcolors => map(x->"#"*hex(x),colors)
                )
            )
        @fact line(verts) -->
             Elem(
                :"three-js-line",
                attributes = Dict(
                    :kind => "strip",
                    :x => 0.0,
                    :y => 0.0,
                    :z => 0.0,
                    :rx => 0.0,
                    :ry => 0.0,
                    :rz => 0.0,
                    :xs => [1.0, 2.0, 3.0],
                    :ys => [1.0, 2.0, 3.0],
                    :zs => [1.0, 2.0, 3.0],
                    :vertexcolors => Color[]
                )
            )
    end
    context("Testing line material") do
        @fact linematerial() --> Elem(:"three-js-line-material", attributes=Dict())
        @fact linematerial(Dict(:kind=>"dashed", :color=>"red")) -->
            Elem(
            :"three-js-line-material",
            attributes = Dict(
                :kind => "dashed",
                :color => "red"
            )
        )
    end
    context("Testing pointclouds") do
        verts = [Tuple{Float64, Float64, Float64}((i, i, i)) for i = 1.0:3.0]
        colors = Color[colorant"red", colorant"blue", colorant"green"]
        vertswithcolors = [
            Tuple{Float64, Float64, Float64, Color}(
                (verts[i][1], verts[i][2], verts[i][3], colors[i])
            ) for i=1:3
        ]
        @fact pointcloud(
        verts, x = 10.0, y = 10.0, z = 10.0, rx = 20.0, ry = 15.0, rz = 240.0,
        vertexcolors = colors) -->
             Elem(
                :"three-js-points",
                attributes = Dict(
                    :x => 10.0,
                    :y => 10.0,
                    :z => 10.0,
                    :rx => 20.0,
                    :ry => 15.0,
                    :rz => 240.0,
                    :xs => [1.0, 2.0, 3.0],
                    :ys => [1.0, 2.0, 3.0],
                    :zs => [1.0, 2.0, 3.0],
                    :vertexcolors => map(x->"#"*hex(x),colors)
                )
            )
        @fact pointcloud(
        vertswithcolors, x = 10.0, y = 10.0, z = 10.0, rx = 20.0, ry = 15.0,
        rz = 240.0) -->
             Elem(
                :"three-js-points",
                attributes = Dict(
                    :x => 10.0,
                    :y => 10.0,
                    :z => 10.0,
                    :rx => 20.0,
                    :ry => 15.0,
                    :rz => 240.0,
                    :xs => [1.0, 2.0, 3.0],
                    :ys => [1.0, 2.0, 3.0],
                    :zs => [1.0, 2.0, 3.0],
                    :vertexcolors => map(x->"#"*hex(x),colors)
                )
            )
        @fact pointcloud(verts) -->
             Elem(
                :"three-js-points",
                attributes = Dict(
                    :x => 0.0,
                    :y => 0.0,
                    :z => 0.0,
                    :rx => 0.0,
                    :ry => 0.0,
                    :rz => 0.0,
                    :xs => [1.0, 2.0, 3.0],
                    :ys => [1.0, 2.0, 3.0],
                    :zs => [1.0, 2.0, 3.0],
                    :vertexcolors => Color[]
                )
            )
    end
    context("Testing point material") do
        @fact pointmaterial() --> Elem(:"three-js-point-material", attributes=Dict())
        @fact pointmaterial(
            Dict(:color => "red", :size=> 10,
            :attenuation => false, :colorkind => "vertex")
        ) -->
            Elem(
            :"three-js-point-material",
            attributes = Dict(
                :colorkind => "vertex",
                :color => "red",
                :size => 10,
            )
        )
    end
    context("Testing material") do
        @fact material() --> Elem(:"three-js-material", attributes = Dict())
        @fact material(Dict(:kind=>"basic", :color=>"red", :wireframe=>false)) -->
            Elem(
            :"three-js-material",
            attributes = Dict(
                :kind => "basic",
                :color => "red"
            )
        )
    end
    context("Testing grid") do
        @fact grid(10.0, 1.0) --> Elem(
            :"three-js-grid",
            attributes = Dict(
                :size => 10.0,
                :step => 1.0,
                :x => 0.0,
                :y => 0.0,
                :z => 0.0,
                :rx => 0.0,
                :ry => 0.0,
                :rz => 0.0,
                :colorgrid => "#000000",
                :colorcenter => "#000000"
            )
        )
        @fact grid(
            10.0, 2.0, x = 10.0, colorgrid = colorant"red",
            colorcenter = colorant"green", ry = 10.0
        ) --> Elem(
            :"three-js-grid",
            attributes = Dict(
                :size => 10.0,
                :step => 2.0,
                :x => 10.0,
                :y => 0.0,
                :z => 0.0,
                :rx => 0.0,
                :ry => 10.0,
                :rz => 0.0,
                :colorgrid => "#FF0000",
                :colorcenter => "#008000"
            )
        )
    end
    context("Testing Text") do
        @fact text(0.0, 0.0, 0.0, "Hello!") --> Elem(
            :"three-js-text",
            attributes = Dict(:x => 0.0, :y => 0.0, :z => 0.0, :content => "Hello!")
        )
    end
end

facts("Testing property helpers") do
    context("Testing meshcolor") do
        @fact meshcolor(colorant"red") --> (:color, "#FF0000")
    end
    context("Testing kinds") do
        @fact lambert() --> (:kind, "lambert")
        @fact basic() --> (:kind, "basic")
        @fact phong() --> (:kind, "phong")
        @fact normal() --> (:kind, "normal")
    end
    context("Testing wireframe") do
        @fact wireframe() --> (:wireframe, true)
        @fact wireframe(false) --> (:wireframe, false)
    end
    context("Testing visibility") do
        @fact visible() --> (:hidden, false)
        @fact visible(false) --> (:hidden, true)
    end
end

FactCheck.exitstatus()
end


text:
# Scans methods in Base and its sub-modules (excluding Core) and counts the number of instances
# of each unique fragment

using BasicBlockRewriter
using MethodAnalysis
using PyPlot: PyPlot, plt

"""
    collect_methods(; exclude=())

Collect all methods in the session, excluding modules listed in `exclude`.

    collect_methods(SomeModule; exclude=())

Collect all methods in `SomeModule` or one of its sub-modules, excluding modules listed in `exclude`.
"""
function collect_methods(top::Union{Tuple{},Tuple{Module}}=(); exclude=())
    meths = Method[]
    visit(top...) do item
        isa(item, Module) && return item ∉ exclude
        if isa(item, Method)
            push!(meths, item)
            return false
        end
        return true
    end
    return meths
end
collect_methods(top::Module; kwargs...) = collect_methods((top,); kwargs...)

function catalog_fragments!(fragdict, m::Method)
    try
        for frag in fragments(m)
            fragdict[frag] = get(fragdict, frag, 0) + 1
        end
    catch err
        if isa(err, ErrorException) && err.msg == "Code for this Method is not available."
            @warn "No code available for $m"
            return fragdict
        end
        @error "$err when processing $m"
    end
    return fragdict
end

meths = collect_methods(; exclude=child_modules(Core))
fragdict = Dict{Vector{Any},Int}()
for m in meths
    catalog_fragments!(fragdict, m)
end

if lowercase(get(ENV, "CI", "false")) != "true"
    fig, axs = plt.subplots(1, 3, figsize=(8, 3))
    ax = axs[1]
    ax.hist(map(length, collect(keys(fragdict))))
    ax.set_yscale("log")
    ax.set_xlabel("# statements")
    ax.set_ylabel("# fragments")

    ax = axs[2]
    bins = [1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]
    ax.hist(collect(values(fragdict)), bins)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_xlabel("# callers")
    ax.set_ylabel("# fragments")

    ax = axs[3]
    kv = collect(fragdict)
    ax.scatter(length.(first.(kv)), last.(kv))
    ax.set_yscale("log")
    ax.set_xlabel("# statements")
    ax.set_ylabel("# callers")

    fig.tight_layout()
end


text:
"""
    Exo4Real{T <: Real} = Parabolic{Exo3Real{T}}

A 4-exo-real number is a parabolic Cayley-Dickson construct with 3-exo-real numbers.
"""
const Exo4Real{T <: Real} = Parabolic{Exo3Real{T}}

function Exo4Real(z::Exo3Real{T}) where T <: Real
    Exo4Real{T}(z, zero(Exo3Real{T}))
end

function Exo4Real(z::Exo2Real{T}) where T <: Real
    Exo4Real{T}(
        Exo3Real{T}(z, zero(Exo2Real{T})),
        zero(Exo3Real{T})
    )
end

function Exo4Real(z::Exo1Real{T}) where T <: Real
    Exo4Real{T}(
        Exo3Real{T}(
            Exo2Real{T}(z, zero(Exo1Real{T})),
            zero(Exo2Real{T})),
        zero(Exo3Real{T})
    )
end

function Exo4Real(v::T...) where T <: Real
    n = length(v)
    if n <= 8
        return Exo4Real{T}(Exo3Real(v...), zero(Exo3Real{T}))
    end

    if n > 16
        return Exo4Real{T}(Exo3Real(v[1:8]...), Exo3Real(v[9:16]...))
    end

    Exo4Real{T}(Exo3Real(v[1:8]...), Exo3Real(v[9:end]...))
end

function selfstar(z::Exo4Real)
    (z + star(z) + dagger(z) + star(dagger(z))) / 4
end

function anti_selfstar(z::Exo4Real)
    (z - star(z) + dagger(z) - star(dagger(z))) / 4
end

"""
    unreal(z::Exo4Real)

The unreal part of a 4-exo-real number is a 15-vector.
"""
function unreal(z::Exo4Real)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Exo4Real)
    print(io, "[1: ")
    print(io, z.l.l.l.l)
    print(io, ", a: ")
    print(io, z.l.l.l.r)
    print(io, ", b: ")
    print(io, z.l.l.r.l)
    print(io, ", ab: ")
    print(io, z.l.l.r.r)
    print(io, ", c: ")
    print(io, z.l.r.l.l)
    print(io, ", ac: ")
    print(io, z.l.r.l.r)
    print(io, ", bc: ")
    print(io, z.l.r.r.l)
    print(io, ", (ab)c: ")
    print(io, z.l.r.r.r)
    print(io, ", d: ")
    print(io, z.r.l.l.l)
    print(io, ", ad: ")
    print(io, z.r.l.l.r)
    print(io, ", bd: ")
    print(io, z.r.l.r.l)
    print(io, ", (ab)d: ")
    print(io, z.r.l.r.r)
    print(io, ", cd: ")
    print(io, z.r.r.l.l)
    print(io, ", (ac)d: ")
    print(io, z.r.r.l.r)
    print(io, ", (bc)d: ")
    print(io, z.r.r.r.l)
    print(io, ", ((ab)c)d: ")
    print(io, z.r.r.r.r)
    print(io, "]")
end

function random(::Type{Exo4Real{T}}) where T <: Real
    Exo4Real{T}(random(Exo3Real{T}), random(Exo3Real{T}))
end


text:
struct Landau

    ϵ  :: Float64
    kx :: Float64

end

function distribution( mesh1  :: UniformMesh,
                       mesh2  :: UniformMesh,
		       landau :: Landau)

    x = mesh1.points
    v = mesh2.points
    ϵ, kx = landau.ϵ, landau.kx

    (1.0.+ϵ*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.*v))

end


text:
# Generate uniform random data

using Revise
using SMLMSim
using SMLMData
using PlotlyJS
using Plots

## create a blinking fluorophore
γ=1e5
q=[0 50
   1e-2 0]
f=SMLMSim.GenericFluor(γ,q)

# pattern is a single molecule
pattern=SMLMSim.Point2D()
ρ=1.0
xsize=25.6
ysize=25.6
smd_true=SMLMSim.uniform2D(ρ,pattern,xsize,ysize)
plt=PlotlyJS.plot(scattergl(x=smd_true.x, y=smd_true.y, mode="markers"))
display(plt)

# generate the kinetic model
nframes=2000
framerate=50.0
smd_model=SMLMSim.kineticmodel(smd_true,f,nframes,framerate;ndatasets=10)
plt=Plots.histogram(smd_model.photons;xlabel="photons",ylabel="counts")
display(plt)

# make noisy coordinates
σ_psf=1.3
smd_noisy=SMLMSim.noise(smd_model,σ_psf)
plt=PlotlyJS.plot(scattergl(x=smd_noisy.x, y=smd_noisy.y, mode="markers"))
display(plt)


text:
# Check if any binary operator are overly complex
function flagBinOperatorComplexity(tree::Node, ::Val{op}, options::Options)::Bool where {op}
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        return flagBinOperatorComplexity(tree.l, Val(op), options)
    else
        if tree.op == op
            overly_complex::Bool = (
                    ((options.bin_constraints[op][1]::Int > -1) &&
                     (countNodes(tree.l) > options.bin_constraints[op][1]::Int))
                      ||
                    ((options.bin_constraints[op][2]::Int > -1) &&
                     (countNodes(tree.r) > options.bin_constraints[op][2]::Int))
                )
            if overly_complex
                return true
            end
        end
        return (flagBinOperatorComplexity(tree.l, Val(op), options) || flagBinOperatorComplexity(tree.r, Val(op), options))
    end
end

# Check if any unary operators are overly complex
function flagUnaOperatorComplexity(tree::Node, ::Val{op}, options::Options)::Bool where {op}
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        if tree.op == op
            overly_complex::Bool = (
                      (options.una_constraints[op]::Int > -1) &&
                      (countNodes(tree.l) > options.una_constraints[op]::Int)
                )
            if overly_complex
                return true
            end
        end
        return flagUnaOperatorComplexity(tree.l, Val(op), options)
    else
        return (flagUnaOperatorComplexity(tree.l, Val(op), options) || flagUnaOperatorComplexity(tree.r, Val(op), options))
    end
end


text:
using Images, SIUnits
using Base.Test

if !isdefined(:workdir)
    const workdir = joinpath(tempdir(), "Images")
end
if !isdefined(:writedir)
    const writedir = joinpath(workdir, "write")
end

if !isdir(workdir)
    mkdir(workdir)
end
if !isdir(writedir)
    mkdir(writedir)
end

# Gray, raw
img = imread(joinpath(dirname(@__FILE__), "io", "small.nrrd"))
@assert colorspace(img) == "Gray"
@assert ndims(img) == 3
@assert colordim(img) == 0
@assert eltype(img) == Float32
outname = joinpath(writedir, "small.nrrd")
imwrite(img, outname)
imgc = imread(outname)
@assert img.data == imgc.data

img = imread(joinpath(dirname(@__FILE__), "io", "units.nhdr"))
ps = pixelspacing(img)
@test_approx_eq ps[1]/(0.1*Milli*Meter) 1
@test_approx_eq ps[2]/(0.2*Milli*Meter) 1
@test_approx_eq ps[3]/(1*Milli*Meter) 1

# Gray, compressed (gzip)
img = imread(joinpath(dirname(@__FILE__), "io", "smallgz.nrrd"))
@assert colorspace(img) == "Gray"
@assert ndims(img) == 3
@assert colordim(img) == 0
@assert eltype(img) == Float32
outname = joinpath(writedir, "smallgz.nrrd")
imwrite(img, outname)
imgc = imread(outname)
@assert img.data == imgc.data


text:
##########################
# Gates
##########################

type Gate{T,N}
    name::AbstractString
    op::AbstractOp{T,N}
    Gate{T,N}(name::AbstractString,op::AbstractOp{T,N}) = new(name,op)
end

Gate{N}(name::AbstractString,op::MatrixOp{N})=Gate{AbstractMatrix,N}(name,op)

Hadamard = Gate("Hadamard gate",OP_Hadamard)
X        = Gate("Paili X",OP_sigmax)
Y        = Gate("Paili Y",OP_sigmay)
Z        = Gate("Paili Z",OP_sigmaz)

function call{T,N}(gate::Gate{T,N})
    return gate.op()
end

function call{T,N}(gate::Gate{T,N},x)
    return gate.op(x)
end

"""
An AbstractGateUnit should at least contain one element called `pos`,
for definition of `pos`, please refer to the docs of `GateUnit`
"""
abstract AbstractGateUnit{N}

"""
GateUnit
---

GateUnit is for the gate unit in a circuit

pos:
the first number in pos is the the column number
the second number in pos is the bits' IDs that is realted to this gate
eg.

'''
  --block 1--|----block 2----|
1 -----------|---------------|
2 ----[A]----|---------------|
3 -----------|---------------|
4 -----------|-----[   ]-----|
5 -----------|-----[ B ]-----|
6 -----------|-----[   ]-----|
'''

The position of gate B is (2,4,5,6)
"""

type GateUnit{T,N}<:AbstractGateUnit{N}
    gate::Gate{T,N}
    pos::Vector{Int}
    # TODO : time layer?
    # time_layer::Real

    GateUnit{T,N}(gate::Gate{T,N},pos::Vector{Int})=new(gate,pos)
end

GateUnit{T,N}(gate::Gate{T,N},pos::Vector{Int}) = GateUnit{T,N}(gate,pos)
GateUnit{T,N}(gate::Gate{T,N},pos::Tuple{Int}) = GateUnit{T,N}(gate,collect(pos))
GateUnit{T,N}(gate::Gate{T,N},pos::Int...) = GateUnit{T,N}(gate,collect(pos))

type CtrlGateUnit{T,N}<:AbstractGateUnit{N}
    gate::Gate{T,N}
    pos::Vector{Int}
    ctrl::Vector{Int}
end

# CtrlGateUnit{T,N}(gate::Gate{T,N},pos::Vector{Int},ctrl::Vector{Int}) = CtrlGateUnit{T,N}(gate,pos,ctrl)

bitnum{N}(unit::AbstractGateUnit{N}) = N


text:
# [[file:../fresa.org::testmultiobjective][testmultiobjective]]
using Fresa
nx = 2
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> ones(length(x)))
# initial point in domain
x = rand(nx)
# objective function
f = x -> ( [sin(x[1]-x[2]); cos(x[1]+x[2])], 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x,f)]
# now invoke Fresa to solve the problem
pareto, population = Fresa.solve(f, p0, domain;
                                 #fitnesstype = :hadamard,
                                 #fitnesstype = :borda,
                                 fitnesstype = :nondominated,
                                 ngen=200,
                                 npop=20,
                                 plotvectors=true,
                                 tolerance=0.01)

println("**** Pareto front:")
println("#+plot: ind:2 deps:(3) with:points")
println(population[pareto])
#using BenchmarkTools
#@benchmark
# testmultiobjective ends here


text:
# ==========================================================================================================
# Backend funcion

formula(trm::TableRegressionModel) = trm.mf.f

function deviances(trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}, <: AbstractArray{T}};
                    type::Int = 1, kwargs...) where {T <: BlasReal}
    isa(trm.model.pp, DensePredChol) || throw(
            ArgumentError("Methods for other PredChol types is not implemented; use model with DensesPredChol instead."))

    assign, f = trm.mm.assign, formula(trm).rhs
    # Determine null model by link function
    ## model with 0 factors as null
    ## if 0 factor is not allowed,
    ### eg, ~ 1 + A + B + ...
    ### ~ 1 as null
    ### eg, ~ 0 + A + A & B
    ### ~ A as null
    # start value
    start = isnullable(trm.model) ? 0 : first(assign)
    err1 = ArgumentError("Invalid set of model specification for ANOVA; not enough variables provided.")
    err2 = ArgumentError("Invalid set of model specification for ANOVA; try adding variables without zeros.")
    if type == 1
        todel = union(start, assign)
        # ~ 0 + A, ~ 1
        start > 0 && (length(todel) > 1 || throw(err1))
        devs = zeros(T, length(todel) + 1)
        @inbounds for id in eachindex(todel)
            devs[id] = deviance(trm, @view(todel[id + 1:end]); kwargs...)
        end
        devs = -diff(devs)
    elseif type == 2
        # Problems with Fullrank promotion occur
        ## ~ 0 + A (categorical) + B (categorical)
        ## A is aliased with 1, but 1 is not in formula.
        ## A is promoted to full rank.
        ## B is also aliased with 1, but it is not promoted.
        ## Doing type 2 and 3 anova cause problem
        ## when deleting promoted categorical variables.
        # Some warning or checks?
        todel = unique(assign)
        if start > 0
            # ~ 0 + A + A & B, all terms are related to A, ~ A as null
            selectcoef(f, Val(first(todel))) == Set(todel) && popfirst!(todel)
            # ~ 0 + A, ~ 1
            if isempty(todel)
                throw(err1)
            # ~ 0 + 2 categorical
            # Do not inspect InteractionTerm
            elseif length(todel) == 2
                all(t -> isa(t, CategoricalTerm), f.terms[todel]) && throw(err2)
            end
        end
        devs = zeros(T, length(todel) + 1)
        # cache fitted
        dict_devs = Dict{Set{Int}, T}()
        @inbounds for (id, del) in enumerate(todel)
            delcoef = selectcoef(f, Val(del))
            dev1 = get(dict_devs, delcoef, (push!(dict_devs, delcoef => deviance(trm, delcoef; kwargs...)); dict_devs[delcoef]))
            delete!(delcoef, del)
            dev2 = get(dict_devs, delcoef, (push!(dict_devs, delcoef => deviance(trm, delcoef; kwargs...)); dict_devs[delcoef]))
            devs[id] = dev1 - dev2
        end
        devs[end] = deviance(trm, 0; kwargs...)
    else
        todel = unique(assign)
        if start > 0
            # ~ 0 + A, ~ 1
            if length(todel) <= 1
                throw(err1)
            # ~ 0 + 2 categorical or
            # ~ 1 + 1 categorical
            # Do not inspect InteractionTerm
            elseif length(todel) == 2
                if isa(f.terms[todel][1], Union{CategoricalTerm, InterceptTerm{true}})
                    isa(f.terms[todel][1], CategoricalTerm) && throw(err2)
                end
            end
        end
        devs = zeros(T, length(todel) + 1)
        @inbounds for (id, del) in enumerate(todel)
            devs[id] = deviance(trm, del; kwargs...)
        end
        devr = deviance(trm, 0; kwargs...)
        devs .-= devr
        devs[end] = devr
    end
    # every method end with deviance(trm, 0; kwargs...), ie full model fit.
    installbeta!(trm.model.pp) # ensure model unchanged
    tuple(devs...)
end

submm(trm::TableRegressionModel, exclude) =
    trm.mm.m[:, map(!in(exclude), trm.mm.assign)]
submm(trm::TableRegressionModel, exclude::Int) =
    trm.mm.m[:, map(!=(exclude), trm.mm.assign)]

# Backend for LinearModel
# Only used by type 2 SS
function deviance(trm::TableRegressionModel{<: LinearModel}, exclude)
    p = trm.model.pp

    # reschema to calculate new X, reset beta
    X = submm(trm, exclude)
    p.beta0 = p.delbeta = repeat([0], size(X, 2))
    p.scratchbeta = similar(p.beta0)

    # cholesky
    F = Hermitian(float(X'X))
    p.chol = updatechol(p.chol, F)

    # reset scratch
    p.scratchm1 = similar(X)
    p.scratchm2 = similar(p.chol.factors)

    isempty(trm.model.rr.wts) ? delbeta!(p, X, trm.model.rr.y) : delbeta!(p, X, trm.model.rr.y, trm.model.rr.wts)
    updateμ!(trm.model.rr, linpred(p, X))
    # installbeta is ommited
end

updatechol(::CholeskyPivoted, F::AbstractMatrix{<: BlasReal}) =
    cholesky!(F, Val(true), tol = -one(eltype(F)), check = false)

updatechol(::Cholesky, F::AbstractMatrix{<: BlasReal}) = cholesky!(F)

# Backend for GeneralizedLinearModel
function deviance(trm::TableRegressionModel{<: GeneralizedLinearModel{<: GlmResp{<: GLM.FPVector, <: GLM.UnivariateDistribution, <: Link}}}, exclude;
                    verbose::Bool = false,
                    maxiter::Integer = 30,
                    minstepfac::Real = 0.001,
                    atol::Real = 1e-6,
                    rtol::Real = 1e-6,
                    kwargs...)

    # Check arguments
    maxiter >= 1       || throw(ArgumentError("maxiter must be positive"))
    0 < minstepfac < 1 || throw(ArgumentError("minstepfac must be in (0, 1)"))

    # Extract fields and set convergence flag
    cvg, p, r = false, trm.model.pp, trm.model.rr
    lp = r.mu

    # subset X, reset beta
    X = submm(trm, exclude)
    p.beta0 = p.delbeta = repeat([0], size(X, 2))
    p.scratchbeta = similar(p.beta0)
    # cholesky
    F = Hermitian(float(X'X))
    p.chol = updatechol(p.chol, F)

    # reset scratch
    p.scratchm1 = similar(X)
    p.scratchm2 = similar(p.chol.factors)

    # Initialize β, μ, and compute deviance
    delbeta!(p, X, GLM.wrkresp(r), r.wrkwt)
    linpred!(lp, p, X)
    updateμ!(r, lp)
    installbeta!(p)
    devold = deviance(trm.model)

    for i = 1:maxiter
        f = 1.0 # line search factor
        local dev

        # Compute the change to β, update μ and compute deviance
        try
            delbeta!(p, X, r.wrkresid, r.wrkwt)
            linpred!(lp, p, X)
            updateμ!(r, lp)
            dev = deviance(trm.model)
        catch e
            isa(e, DomainError) ? (dev = Inf) : rethrow(e)
        end

        # Line search
        ## If the deviance isn't declining then half the step size
        ## The rtol*dev term is to avoid failure when deviance
        ## is unchanged except for rouding errors.
        while dev > devold + rtol*dev
            println("ok")
            f /= 2
            f > minstepfac || error("step-halving failed at beta0 = $(p.beta0)")
            try
                updateμ!(r, linpred(p, X, f))
                dev = deviance(trm.model)
            catch e
                isa(e, DomainError) ? (dev = Inf) : rethrow(e)
            end
        end
        installbeta!(p, f)

        # Test for convergence
        verbose && println("Iteration: $i, deviance: $dev, diff.dev.:$(devold - dev)")
        if devold - dev < max(rtol*devold, atol)
            cvg = true
            break
        end
        @assert isfinite(dev)
        devold = dev
    end
    cvg || throw(ConvergenceException(maxiter))
    devold
end

# Linear prediction
function delbeta!(p::DensePredChol{<: BlasReal, <: Cholesky},
                    X::Matrix{<: BlasReal},
                    r::Vector{<: BlasReal})
    ldiv!(p.chol, mul!(p.delbeta, transpose(X), r))
    p
end
# β = (X'X)⁻¹X'y

function delbeta!(p::DensePredChol{<: BlasReal, <: CholeskyPivoted},
                    X::Matrix{<: BlasReal},
                    r::Vector{<: BlasReal})
    ch = p.chol
    beta = mul!(p.delbeta, adjoint(X), r)
    rnk = rank(ch)
    if rnk == length(beta)
        ldiv!(ch, beta)
    elseif size(ch.factors) == (0, 0) # Rank undtermined when the matrix is 0x0
        ldiv!(ch, beta)
    else
        permute!(beta, ch.piv)
        for k = (rnk + 1):length(beta)
            beta[k] = -zero(T)
        end
        LAPACK.potrs!(ch.uplo, view(ch.factors, 1:rnk, 1:rnk), view(beta, 1:rnk))
        invpermute!(beta, ch.piv)
    end
    p
end

# exterior weights
function delbeta!(p::DensePredChol{<: BlasReal, <: Cholesky},
                    X::Matrix{<: BlasReal},
                    r::Vector{<: BlasReal},
                    wt::Vector{<: BlasReal})
    scr = mul!(p.scratchm1, Diagonal(wt), X)
    cholesky!(Hermitian(mul!(cholfactors(p.chol), transpose(scr), X), :U))
    mul!(p.delbeta, transpose(scr), r)
    ldiv!(p.chol, p.delbeta)
    p
end

function delbeta!(p::DensePredChol{<: BlasReal, <: CholeskyPivoted},
                    X::Matrix{<: BlasReal},
                    r::Vector{<: BlasReal},
                    wt::Vector{<: BlasReal})
    cf = cholfactors(p.chol)
    piv = p.chol.piv
    cf .= mul!(p.scratchm2, adjoint(mul!(p.scratchm1, Diagonal(wt), X)), X)[piv, piv]
    cholesky!(Hermitian(cf, Symbol(p.chol.uplo)))
    ldiv!(p.chol, mul!(p.delbeta, transpose(p.scratchm1), r))
    p
end

# experimental
function delbeta!(p::SparsePredChol{<: BlasReal},
                    X::Matrix{<: BlasReal},
                    r::Vector{<: BlasReal},
                    wt::Vector{<: BlasReal})
    scr = mul!(p.scratch, Diagonal(wt), X)
    XtWX = X'*scr
    c = p.chol = cholesky(Symmetric{eltype(XtWX),typeof(XtWX)}(XtWX, 'L'))
    p.delbeta = c \ mul!(p.delbeta, adjoint(scr), r)
end

# skip updating beta0
linpred(p::LinPred, X::Matrix, f::Real = 1.0) =
    linpred!(Vector{eltype(p.X)}(undef, size(p.X, 1)), p, X, f)

function linpred!(out, p::LinPred, X::Matrix, f::Real = 1.0)
    mul!(out, X, iszero(f) ? p.beta0 : broadcast!(muladd, p.scratchbeta, f, p.delbeta, p.beta0))
end

# Create nestedmodels
function nestedmodels(trm::TableRegressionModel{<: LinearModel}; null::Bool = true, kwargs...)
    null = null && isnullable(trm.model.pp.chol)
    assign = unique(trm.mm.assign)
    pop!(assign)
    dropcollinear, range = null ? (false, union(0, assign)) : (true, assign)
    wts = trm.model.rr.wts
    trms = map(range) do id
        # create sub-formula, modify schema, create mf and mm
        mf, mm = subtablemodel(trm, id)
        y = response(mf)
        TableRegressionModel(fit(trm.mf.model, mm.m, y; wts, dropcollinear, kwargs...), mf, mm)
    end
    (trms..., trm)
end

function nestedmodels(trm::TableRegressionModel{<: GeneralizedLinearModel}; null::Bool = true, kwargs...)
    null = null && isnullable(trm.model)
    distr = trm.model.rr.d
    link = typeof(trm.model.rr).parameters[3]()
    wts = trm.model.rr.wts
    offset = trm.model.rr.offset
    assign = unique(trm.mm.assign)
    pop!(assign)
    range = null ? union(0, assign) : assign
    trms = map(range) do id
        # create sub-formula, modify schema, create mf and mm
        mf, mm = subtablemodel(trm, id)
        y = response(mf)
        TableRegressionModel(fit(trm.mf.model, mm.m, y, distr, link; wts, offset, kwargs...), mf, mm)
    end
    (trms..., trm)
end

nestedmodels(::Type{LinearModel}, formula, data; null::Bool = true, kwargs...) =
    nestedmodels(lm(formula, data; kwargs...); null, kwargs...)
nestedmodels(::Type{GeneralizedLinearModel}, formula, data,
                distr::UnivariateDistribution, link::Link = canonicallink(distr);
                null::Bool = true, kwargs...) =
    nestedmodels(glm(formula, data, distr, link; kwargs...); null, kwargs...)

# backend for implementing nestedmodels or dropterms (in future)
function subtablemodel(trm::TableRegressionModel, id; reschema::Bool = false)
    # create sub-formula, modify schema, create mf and mm
    f = formula(trm)
    subf = subformula(f.lhs, f.rhs, id; reschema)
    if reschema
        #=
        terms = setdiff(getterms(f.rhs), getterms(subf.rhs))
        filter!(!=(Symbol("1")), terms)
        pair = collect(pairs(trm.mf.data))
        filter!(x->!in(x.first, terms), pair)
        =#
        mf = ModelFrame(subf, trm.mf.data;
                        model = trm.mf.model, contrasts = extract_contrasts(f))
    else
        #schema = deepcopy(trm.mf.schema)
        #for term in terms
        #    pop!(schema.schema, Term(term))
        #end
        mf = ModelFrame(subf, schema, trm.mf.data, trm.mf.model)
    end
    mm = ModelMatrix(mf)
    mf, mm
end

# Null model for CholeskyPivoted is unstable
## For nestedmodels
isnullable(::CholeskyPivoted) = false
isnullable(::Cholesky) = true
isnullable(::InverseLink) = false
isnullable(::InverseSquareLink) = false
isnullable(::Link) = true

## For deviances
isnullable(::LinearModel) = true
isnullable(model::GeneralizedLinearModel) = isnullable(typeof(model.rr).parameters[3]())


text:
## gmm/types.jl  The type that implements a GMM.
## (c) 2013--2015 David A. van Leeuwen

## Some remarks on the dimension.  There are three main indexing variables:
## - The gaussian index
## - The data point
## - The feature dimension
## Often data is stores in 2D arrays, and computations can be done efficiently in
## matrix multiplications.  For this it is nice to have the data in standard row,column order
## however, we can't have these consistently over all three indexes.

## My approach is do have:
## - The data index (i) always be a row-index
## - The feature dimenssion index (k) always to be a column index
## - The Gaussian index (j) to be mixed, depending on how it is used
## For full covar Σ, we need 2 feature dims.  For individuals covars to be
## consecutive in memory, the gaussian index should be _last_.

## force Emacs utf-8: αβγδεζηθικλμνξοπρστυφχψω

using Printf

"""
`History`, a type to record the history of how a GMM is built.
"""
struct History
    """timestamp"""
    t::Float64
    """description"""
    s::AbstractString
end
History(s::AbstractString) = History(time(), s)

"""
`GaussianMixture`, an abstract type for a mixture of full-covariance or diagonal-covariance Gaussian
distributions
"""
abstract type GaussianMixture{T,CT}; end

## support for two kinds of covariance matrix
## Full covariance is represented by inverse cholesky of the covariance matrix,
## i.e., Σ^-1 = ci * ci'
DiagCov{T} = AbstractArray{T,2}
FullCov{T} = Vector{UpperTriangular{T,Matrix{T}}}
CovType{T} = Union{DiagCov{T}, FullCov{T}}

VecOrMat{T} = Union{Vector{T},AbstractArray{T,2}}
MatOrVecMat{T} = Union{AbstractArray{T,2}, Vector{AbstractArray{T,2}}}

## GMMs can be of type FLoat32 or Float64, and diagonal or full
"""
`GMM` is the type that stores information of a Guassian Mixture Model.  Currently two main covariance
types are supported: full covarariance and diagonal covariance.
"""
mutable struct GMM{T<:AbstractFloat, CT<:CovType{T}} <: GaussianMixture{T,CT}
    "number of Gaussians"
    n::Int
    "dimension of Gaussian"
    d::Int
    "weights (size n)"
    w::Vector{T}
    "means (size n x d)"
    μ::Matrix{T}
    "covariances (size n x d for diagonal, or n x (d^2) for full)"
    Σ::CT
    "history"
    hist::Vector{History}
    "number of points used to train the GMM"
    nx::Int
    function GMM{T,CT}(w::Vector{T}, μ::AbstractArray{T,2}, Σ::CT,
                               hist::Vector, nx::Int) where{T, CT}
        n = length(w)
        isapprox(1, sum(w)) || error("weights do not sum to one")
        d = size(μ, 2)
        n == size(μ, 1) || error("Inconsistent number of means")
        if isa(Σ, Matrix)
            (n,d) == size(Σ) || error("Inconsistent covar dimension")
        else
            n == length(Σ) || error(@sprintf("Inconsistent number of covars %d != %d", n, length(Σ)))
            for (i,S) in enumerate(Σ)
                (d,d) == size(S) || error(@sprintf("Inconsistent dimension for %d", i))
##                isposdef(S) || error(@sprintf("Covariance %d not positive definite", i))
            end
        end
        new(n, d, w, μ, Σ, hist, nx)
    end
end
GMM(w::Vector{T}, μ::AbstractArray{T,2}, Σ::Union{DiagCov{T},FullCov{T}},
                      hist::Vector, nx::Int) where {T<:AbstractFloat} = GMM{T, typeof(Σ)}(w, μ, Σ, hist, nx)

## Variational Bayes GMM types.

## Please note our pedantic use of the Greek letter ν (nu), don't confuse this with Latin v!
## The index-0 "₀" is part of the identifier.
"""
`GMMprior` is a type that holds the prior for training GMMs using Variational Bayes.
"""
struct GMMprior{T<:AbstractFloat}
    "effective prior number of observations"
    α₀::T
    β₀::T
    "prior on the mean μ"
    m₀::Vector{T}
    "scale of precision Λ"
    ν₀::T
    "prior of the precision Λ"
    W₀::Matrix{T}
end

## In Variational Bayes, the GMM is not specified by point estimates of the paramters,
## but distributions over these parameters.
## These are Dirichlet for the weights and Gaussian-Wishart for the mean and precision.
## These distributions have parameters themselves, and these are stored in this type...
"""
`VGMM` is the type that is used to store a GMM in the Variational Bayes training.
"""
mutable struct VGMM{T} <: GaussianMixture{T,Any}
    "number of Gaussians"
    n::Int
    "dimension of Gaussian"
    d::Int
    "The prior used in this VGMM"
    π::GMMprior{T}
    "Dirichlet, size n"
    α::Vector{T}
    "scale of precision, size n"
    β::Vector{T}
    "means of means, size n * d"
    m::Matrix{T}
    "no. degrees of freedom, size n"
    ν::Vector{T}
    "scale matrix for precision? size n * (d * d)"
    W::FullCov{T}
    "history"
    hist::Vector{History}
end

## UBM-centered and scaled stats.
## This structure currently is useful for dotscoring, so we've limited the
## order to 1.  Maybe we can make this more general allowing for uninitialized second order
## stats?

## We store the stats in a (ng * d) structure, i.e., not as a super vector yet.
## Perhaps in ivector processing a supervector is easier.
"""
`CSstats` a type holding centered and scaled zeroth and first order GMM statistics
"""
struct CSstats{T<:AbstractFloat}
    "zeroth order stats"
    n::Vector{T}          # zero-order stats, ng
    "first order stats"
    f::Matrix{T}          # first-order stats, ng * d
    function CSstats{T}(n::Vector, f::Matrix) where{T}
        @assert size(n,1)==size(f, 1)
        new(n,f)
    end
end
CSstats(n::Vector{T}, f::Matrix{T}) where {T<:AbstractFloat} = CSstats{T}(n, f)
## special case for tuple (why would I need this?)
CSstats(t::Tuple) = CSstats(t[1], t[2])

## Cstats is a type of centered but un-scaled stats, necessary for i-vector extraction
"""
`Cstats`, a type holding centered zeroth, first and second order GMM statistics
"""
struct Cstats{T<:AbstractFloat, CT<:VecOrMat}
    "zeroth order stats"
    N::Vector{T}
    "first order stats"
    F::Matrix{T}
    "second order stats"
    S::CT
    function Cstats{T,CT}(n::Vector{T}, f::Matrix{T}, s::MatOrVecMat{T}) where{T, CT}
        size(n,1) == size(f,1) || error("Inconsistent size 0th and 1st order stats")
        if size(n) == size(s)   # full covariance stats
            all([size(f,2) == size(ss,1) == size(ss,2) for ss in s]) || error("inconsistent size 1st and 2nd order stats")
       else
            size(f) == size(s) || error("inconsistent size 1st and 2nd order stats")
        end
        new(n, f, s)
    end
end
Cstats(n::Vector{T}, f::Matrix{T}, s::MatOrVecMat{T}) where {T<:AbstractFloat} = Cstats{T,typeof(s)}(n, f, s)
Cstats(t::Tuple) = Cstats(t...)

## A data handle, either in memory or on disk, perhaps even mmapped but I haven't seen any
## advantage of that.  It contains a list of either files (where the data is stored)
## or data units.  The point is, that in processing, these units can naturally be processed
## independently.

## The API is a dictionary of functions that help loading the data into memory
## Compulsory is: :load, useful is: :size
"""
`Data` is a type for holding an array of feature vectors (i.e., matrices), or references to
files on disk.  The data is automatically loaded when needed, e.g., by indexing.
"""
struct Data{T,VT<:Union{Matrix,AbstractString}}
    list::Vector{VT}
    API::Dict{Symbol,Function}
    function Data{T,VT}(list::Union{Vector{VT},Vector{Matrix{T}}}, API::Dict{Symbol,Function}) where{T,VT}
        return new(list,API)
    end
end
Data(list::Vector{Matrix{T}}) where {T} = Data{T, eltype(list)}(list, Dict{Symbol,Function}())
Data(list::Vector{S}, t::DataType, API::Dict{Symbol,Function}) where {S<:AbstractString} = Data{t, S}(list, API)

DataOrMatrix{T} = Union{Data{T}, Matrix{T}}


text:
using IonBase
using IonCLI
using Comonicon: Comonicon, PATH, Options, BuildTools

configs = Options.read_configs(IonCLI)
@info "building application"
BuildTools.build_application(IonCLI, configs)
IonBase.copy_assets(joinpath(pkgdir(IonCLI), "build", "ion"))
# pack tarball
tarball = BuildTools.tarball_name(IonCLI, configs.name; application = true)
@info "creating application tarball $tarball"
cd(PATH.project(IonCLI, configs.application.path)) do
    run(`tar -czvf $tarball $(configs.name)`)
end