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

# everything related to a DG semidiscretization in 3D,
# currently limited to Lobatto-Legendre nodes

# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::TreeMesh{3}, equations,
                      dg::DG, RealT, uEltype)
  # Get cells for which an element needs to be created (i.e. all leaf cells)
  leaf_cell_ids = local_leaf_cells(mesh.tree)

  elements = init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)

  interfaces = init_interfaces(leaf_cell_ids, mesh, elements)

  boundaries = init_boundaries(leaf_cell_ids, mesh, elements)

  mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)

  cache = (; elements, interfaces, boundaries, mortars)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
  cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)

  return cache
end

# The methods below are specialized on the volume integral type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}},
                      equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype)
  NamedTuple()
end

function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, equations,
                      volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
  element_ids_dg   = Int[]
  element_ids_dgfv = Int[]

  cache = create_cache(mesh, equations,
                       VolumeIntegralFluxDifferencing(volume_integral.volume_flux_dg),
                       dg, uEltype)

  A4dp1_x = Array{uEltype, 4}
  A4dp1_y = Array{uEltype, 4}
  A4dp1_z = Array{uEltype, 4}
  fstar1_L_threaded  = A4dp1_x[A4dp1_x(undef, nvariables(equations), nnodes(dg)+1, nnodes(dg), nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar1_R_threaded  = A4dp1_x[A4dp1_x(undef, nvariables(equations), nnodes(dg)+1, nnodes(dg), nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar2_L_threaded  = A4dp1_y[A4dp1_y(undef, nvariables(equations), nnodes(dg), nnodes(dg)+1, nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar2_R_threaded  = A4dp1_y[A4dp1_y(undef, nvariables(equations), nnodes(dg), nnodes(dg)+1, nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar3_L_threaded  = A4dp1_z[A4dp1_z(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg)+1)
                             for _ in 1:Threads.nthreads()]
  fstar3_R_threaded  = A4dp1_z[A4dp1_z(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg)+1)
                             for _ in 1:Threads.nthreads()]

  return (; cache..., element_ids_dg, element_ids_dgfv, fstar1_L_threaded, fstar1_R_threaded,
            fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded)
end

function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, equations,
                      volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype)

  A4dp1_x = Array{uEltype, 4}
  A4dp1_y = Array{uEltype, 4}
  A4dp1_z = Array{uEltype, 4}
  fstar1_L_threaded  = A4dp1_x[A4dp1_x(undef, nvariables(equations), nnodes(dg)+1, nnodes(dg), nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar1_R_threaded  = A4dp1_x[A4dp1_x(undef, nvariables(equations), nnodes(dg)+1, nnodes(dg), nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar2_L_threaded  = A4dp1_y[A4dp1_y(undef, nvariables(equations), nnodes(dg), nnodes(dg)+1, nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar2_R_threaded  = A4dp1_y[A4dp1_y(undef, nvariables(equations), nnodes(dg), nnodes(dg)+1, nnodes(dg))
                             for _ in 1:Threads.nthreads()]
  fstar3_L_threaded  = A4dp1_z[A4dp1_z(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg)+1)
                             for _ in 1:Threads.nthreads()]
  fstar3_R_threaded  = A4dp1_z[A4dp1_z(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg)+1)
                             for _ in 1:Threads.nthreads()]

  return (; fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded,
            fstar3_L_threaded, fstar3_R_threaded)
end

# The methods below are specialized on the mortar type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype)
  # TODO: Taal compare performance of different types
  A3d = Array{uEltype, 3}
  fstar_upper_left_threaded  = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                                   for _ in 1:Threads.nthreads()]
  fstar_upper_right_threaded = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                                   for _ in 1:Threads.nthreads()]
  fstar_lower_left_threaded  = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                                   for _ in 1:Threads.nthreads()]
  fstar_lower_right_threaded = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                                   for _ in 1:Threads.nthreads()]
  fstar_tmp1_threaded        = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                                   for _ in 1:Threads.nthreads()]

  (; fstar_upper_left_threaded, fstar_upper_right_threaded,
     fstar_lower_left_threaded, fstar_lower_right_threaded,
     fstar_tmp1_threaded)
end

# TODO: Taal discuss/refactor timer, allowing users to pass a custom timer?

function rhs!(du, u, t,
              mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache)

  # Calculate volume integral
  @trixi_timeit timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Prolong solution to interfaces
  @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(
    cache, u, mesh, equations, dg.surface_integral, dg)

  # Calculate interface fluxes
  @trixi_timeit timer() "interface flux" calc_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg.surface_integral, dg, cache)

  # Prolong solution to boundaries
  @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(
    cache, u, mesh, equations, dg.surface_integral, dg)

  # Calculate boundary fluxes
  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(
    cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg)

  # Prolong solution to mortars
  @trixi_timeit timer() "prolong2mortars" prolong2mortars!(
    cache, u, mesh, equations, dg.mortar, dg.surface_integral, dg)

  # Calculate mortar fluxes
  @trixi_timeit timer() "mortar flux" calc_mortar_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg.mortar, dg.surface_integral, dg, cache)

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

  # Apply Jacobian from mapping to reference element
  @trixi_timeit timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @trixi_timeit timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  return nothing
end

function calc_volume_integral!(du, u,
                               mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralWeakForm,
                               dg::DGSEM, cache)

  @threaded for element in eachelement(dg, cache)
    weak_form_kernel!(du, u, element, mesh,
                      nonconservative_terms, equations,
                      dg, cache)
  end

  return nothing
end

@inline function weak_form_kernel!(du, u,
                                   element, mesh::TreeMesh{3},
                                   nonconservative_terms::Val{false}, equations,
                                   dg::DGSEM, cache, alpha=true)
  # true * [some floating point value] == [exactly the same floating point value]
  # This can (hopefully) be optimized away due to constant propagation.
  @unpack derivative_dhat = dg.basis

  for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, j, k, element)

    flux1 = flux(u_node, 1, equations)
    for ii in eachnode(dg)
      multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], flux1, equations, dg, ii, j, k, element)
    end

    flux2 = flux(u_node, 2, equations)
    for jj in eachnode(dg)
      multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], flux2, equations, dg, i, jj, k, element)
    end

    flux3 = flux(u_node, 3, equations)
    for kk in eachnode(dg)
      multiply_add_to_node_vars!(du, alpha * derivative_dhat[kk, k], flux3, equations, dg, i, j, kk, element)
    end
  end

  return nothing
end

function calc_volume_integral!(du, u,
                               mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralFluxDifferencing,
                               dg::DGSEM, cache)
  @threaded for element in eachelement(dg, cache)
    split_form_kernel!(du, u, element, mesh,
                       nonconservative_terms, equations,
                       volume_integral.volume_flux, dg, cache)
  end
end

@inline function split_form_kernel!(du, u,
                                    element, mesh::TreeMesh{3},
                                    nonconservative_terms::Val{false}, equations,
                                    volume_flux, dg::DGSEM, cache, alpha=true)
  # true * [some floating point value] == [exactly the same floating point value]
  # This can (hopefully) be optimized away due to constant propagation.
  @unpack derivative_split = dg.basis

  # Calculate volume integral in one element
  for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, j, k, element)

    # All diagonal entries of `derivative_split` are zero. Thus, we can skip
    # the computation of the diagonal terms. In addition, we use the symmetry
    # of the `volume_flux` to save half of the possible two-poitn flux
    # computations.

    # x direction
    for ii in (i+1):nnodes(dg)
      u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
      flux1 = volume_flux(u_node, u_node_ii, 1, equations)
      multiply_add_to_node_vars!(du, alpha * derivative_split[i, ii], flux1, equations, dg, i,  j, k, element)
      multiply_add_to_node_vars!(du, alpha * derivative_split[ii, i], flux1, equations, dg, ii, j, k, element)
    end

    # y direction
    for jj in (j+1):nnodes(dg)
      u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
      flux2 = volume_flux(u_node, u_node_jj, 2, equations)
      multiply_add_to_node_vars!(du, alpha * derivative_split[j, jj], flux2, equations, dg, i, j,  k, element)
      multiply_add_to_node_vars!(du, alpha * derivative_split[jj, j], flux2, equations, dg, i, jj, k, element)
    end

    # z direction
    for kk in (k+1):nnodes(dg)
      u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
      flux3 = volume_flux(u_node, u_node_kk, 3, equations)
      multiply_add_to_node_vars!(du, alpha * derivative_split[k, kk], flux3, equations, dg, i, j, k,  element)
      multiply_add_to_node_vars!(du, alpha * derivative_split[kk, k], flux3, equations, dg, i, j, kk, element)
    end
  end
end

@inline function split_form_kernel!(du, u,
                                    element, mesh::TreeMesh{3},
                                    nonconservative_terms::Val{true}, equations,
                                    volume_flux, dg::DGSEM, cache, alpha=true)
  # true * [some floating point value] == [exactly the same floating point value]
  # This can (hopefully) be optimized away due to constant propagation.
  @unpack derivative_split = dg.basis
  symmetric_flux, nonconservative_flux = volume_flux

  # Apply the symmetric flux as usual
  split_form_kernel!(du, u, element, mesh, Val(false), equations, symmetric_flux, dg, cache, alpha)

  # Calculate the remaining volume terms using the nonsymmetric generalized flux
  for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, j, k, element)

    # The diagonal terms are zero since the diagonal of `derivative_split`
    # is zero. We ignore this for now.

    # x direction
    integral_contribution = zero(u_node)
    for ii in eachnode(dg)
      u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
      noncons_flux1 = nonconservative_flux(u_node, u_node_ii, 1, equations)
      integral_contribution = integral_contribution + derivative_split[i, ii] * noncons_flux1
    end

    # y direction
    for jj in eachnode(dg)
      u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
      noncons_flux2 = nonconservative_flux(u_node, u_node_jj, 2, equations)
      integral_contribution = integral_contribution + derivative_split[j, jj] * noncons_flux2
    end

    # z direction
    for kk in eachnode(dg)
      u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
      noncons_flux3 = nonconservative_flux(u_node, u_node_kk, 3, equations)
      integral_contribution = integral_contribution + derivative_split[k, kk] * noncons_flux3
    end

    # The factor 0.5 cancels the factor 2 in the flux differencing form
    multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, dg, i, j, k, element)
  end
end

# TODO: Taal dimension agnostic
function calc_volume_integral!(du, u,
                               mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralShockCapturingHG,
                               dg::DGSEM, cache)
  @unpack element_ids_dg, element_ids_dgfv = cache
  @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral

  # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α
  alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg, cache)

  # Determine element ids for DG-only and blended DG-FV volume integral
  pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)

  # Loop over pure DG elements
  @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg)
    element = element_ids_dg[idx_element]
    split_form_kernel!(du, u, element, mesh,
                       nonconservative_terms, equations,
                       volume_flux_dg, dg, cache)
  end

  # Loop over blended DG-FV elements
  @trixi_timeit timer() "blended DG-FV" @threaded for idx_element in eachindex(element_ids_dgfv)
    element = element_ids_dgfv[idx_element]
    alpha_element = alpha[element]

    # Calculate DG volume integral contribution
    split_form_kernel!(du, u, element, mesh,
                       nonconservative_terms, equations,
                       volume_flux_dg, dg, cache, 1 - alpha_element)

    # Calculate FV volume integral contribution
    fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
               dg, cache, element, alpha_element)
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_volume_integral!(du, u,
                               mesh::TreeMesh{3},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralPureLGLFiniteVolume,
                               dg::DGSEM, cache)
  @unpack volume_flux_fv = volume_integral

  # Calculate LGL FV volume integral
  @threaded for element in eachelement(dg, cache)
    fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
               dg, cache, element, true)
  end

  return nothing
end

@inline function fv_kernel!(du, u,
                            mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}},
                            nonconservative_terms, equations,
                            volume_flux_fv, dg::DGSEM, cache, element, alpha=true)
  @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded = cache
  @unpack inverse_weights = dg.basis

  # Calculate FV two-point fluxes
  fstar1_L = fstar1_L_threaded[Threads.threadid()]
  fstar2_L = fstar2_L_threaded[Threads.threadid()]
  fstar3_L = fstar3_L_threaded[Threads.threadid()]
  fstar1_R = fstar1_R_threaded[Threads.threadid()]
  fstar2_R = fstar2_R_threaded[Threads.threadid()]
  fstar3_R = fstar3_R_threaded[Threads.threadid()]

  calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u,
               mesh, nonconservative_terms, equations, volume_flux_fv, dg, element, cache)

  # Calculate FV volume integral contribution
  for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    for v in eachvariable(equations)
      du[v, i, j, k, element] += ( alpha *
                                   (inverse_weights[i] * (fstar1_L[v, i+1, j, k] - fstar1_R[v, i, j, k]) +
                                    inverse_weights[j] * (fstar2_L[v, i, j+1, k] - fstar2_R[v, i, j, k]) +
                                    inverse_weights[k] * (fstar3_L[v, i, j, k+1] - fstar3_R[v, i, j, k])) )

    end
  end

  return nothing
end

# Calculate the finite volume fluxes inside the elements (**without non-conservative terms**).
@inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u,
                              mesh::TreeMesh{3}, nonconservative_terms::Val{false}, equations,
                              volume_flux_fv, dg::DGSEM, element, cache)

  fstar1_L[:, 1,            :, :] .= zero(eltype(fstar1_L))
  fstar1_L[:, nnodes(dg)+1, :, :] .= zero(eltype(fstar1_L))
  fstar1_R[:, 1,            :, :] .= zero(eltype(fstar1_R))
  fstar1_R[:, nnodes(dg)+1, :, :] .= zero(eltype(fstar1_R))

  for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg)
    u_ll = get_node_vars(u, equations, dg, i-1, j, k, element)
    u_rr = get_node_vars(u, equations, dg, i,   j, k, element)
    flux = volume_flux_fv(u_ll, u_rr, 1, equations) # orientation 1: x direction
    set_node_vars!(fstar1_L, flux, equations, dg, i, j, k)
    set_node_vars!(fstar1_R, flux, equations, dg, i, j, k)
  end

  fstar2_L[:, :, 1           , :] .= zero(eltype(fstar2_L))
  fstar2_L[:, :, nnodes(dg)+1, :] .= zero(eltype(fstar2_L))
  fstar2_R[:, :, 1           , :] .= zero(eltype(fstar2_R))
  fstar2_R[:, :, nnodes(dg)+1, :] .= zero(eltype(fstar2_R))

  for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg)
    u_ll = get_node_vars(u, equations, dg, i, j-1, k, element)
    u_rr = get_node_vars(u, equations, dg, i, j,   k, element)
    flux = volume_flux_fv(u_ll, u_rr, 2, equations) # orientation 2: y direction
    set_node_vars!(fstar2_L, flux, equations, dg, i, j, k)
    set_node_vars!(fstar2_R, flux, equations, dg, i, j, k)
  end

  fstar3_L[:, :, :, 1           ] .= zero(eltype(fstar3_L))
  fstar3_L[:, :, :, nnodes(dg)+1] .= zero(eltype(fstar3_L))
  fstar3_R[:, :, :, 1           ] .= zero(eltype(fstar3_R))
  fstar3_R[:, :, :, nnodes(dg)+1] .= zero(eltype(fstar3_R))

  for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg)
    u_ll = get_node_vars(u, equations, dg, i, j, k-1, element)
    u_rr = get_node_vars(u, equations, dg, i, j, k,   element)
    flux = volume_flux_fv(u_ll, u_rr, 3, equations) # orientation 3: z direction
    set_node_vars!(fstar3_L, flux, equations, dg, i, j, k)
    set_node_vars!(fstar3_R, flux, equations, dg, i, j, k)
  end

  return nothing
end

# Calculate the finite volume fluxes inside the elements (**without non-conservative terms**).
@inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u,
                              mesh::TreeMesh{3}, nonconservative_terms::Val{true}, equations,
                              volume_flux_fv, dg::DGSEM, element, cache)
  volume_flux, nonconservative_flux = volume_flux_fv

  fstar1_L[:, 1,            :, :] .= zero(eltype(fstar1_L))
  fstar1_L[:, nnodes(dg)+1, :, :] .= zero(eltype(fstar1_L))
  fstar1_R[:, 1,            :, :] .= zero(eltype(fstar1_R))
  fstar1_R[:, nnodes(dg)+1, :, :] .= zero(eltype(fstar1_R))

  for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg)
    u_ll = get_node_vars(u, equations, dg, i-1, j, k, element)
    u_rr = get_node_vars(u, equations, dg, i,   j, k, element)

    # Compute conservative part
    flux = volume_flux(u_ll, u_rr, 1, equations) # orientation 1: x direction

    # Compute nonconservative part
    # Note the factor 0.5 necessary for the nonconservative fluxes based on
    # the interpretation of global SBP operators coupled discontinuously via
    # central fluxes/SATs
    flux_L = flux + 0.5 * nonconservative_flux(u_ll, u_rr, 1, equations)
    flux_R = flux + 0.5 * nonconservative_flux(u_rr, u_ll, 1, equations)

    set_node_vars!(fstar1_L, flux_L, equations, dg, i, j, k)
    set_node_vars!(fstar1_R, flux_R, equations, dg, i, j, k)
  end

  fstar2_L[:, :, 1           , :] .= zero(eltype(fstar2_L))
  fstar2_L[:, :, nnodes(dg)+1, :] .= zero(eltype(fstar2_L))
  fstar2_R[:, :, 1           , :] .= zero(eltype(fstar2_R))
  fstar2_R[:, :, nnodes(dg)+1, :] .= zero(eltype(fstar2_R))

  for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg)
    u_ll = get_node_vars(u, equations, dg, i, j-1, k, element)
    u_rr = get_node_vars(u, equations, dg, i, j,   k, element)

    # Compute conservative part
    flux = volume_flux(u_ll, u_rr, 2, equations) # orientation 2: y direction

    # Compute nonconservative part
    # Note the factor 0.5 necessary for the nonconservative fluxes based on
    # the interpretation of global SBP operators coupled discontinuously via
    # central fluxes/SATs
    flux_L = flux + 0.5 * nonconservative_flux(u_ll, u_rr, 2, equations)
    flux_R = flux + 0.5 * nonconservative_flux(u_rr, u_ll, 2, equations)

    set_node_vars!(fstar2_L, flux_L, equations, dg, i, j, k)
    set_node_vars!(fstar2_R, flux_R, equations, dg, i, j, k)
  end

  fstar3_L[:, :, :, 1           ] .= zero(eltype(fstar3_L))
  fstar3_L[:, :, :, nnodes(dg)+1] .= zero(eltype(fstar3_L))
  fstar3_R[:, :, :, 1           ] .= zero(eltype(fstar3_R))
  fstar3_R[:, :, :, nnodes(dg)+1] .= zero(eltype(fstar3_R))

  for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg)
    u_ll = get_node_vars(u, equations, dg, i, j, k-1, element)
    u_rr = get_node_vars(u, equations, dg, i, j, k,   element)

    # Compute conservative part
    flux = volume_flux(u_ll, u_rr, 3, equations) # orientation 3: z direction

    # Compute nonconservative part
    # Note the factor 0.5 necessary for the nonconservative fluxes based on
    # the interpretation of global SBP operators coupled discontinuously via
    # central fluxes/SATs
    flux_L = flux + 0.5 * nonconservative_flux(u_ll, u_rr, 3, equations)
    flux_R = flux + 0.5 * nonconservative_flux(u_rr, u_ll, 3, equations)

    set_node_vars!(fstar3_L, flux_L, equations, dg, i, j, k)
    set_node_vars!(fstar3_R, flux_R, equations, dg, i, j, k)
  end

  return nothing
end

function prolong2interfaces!(cache, u,
                             mesh::TreeMesh{3}, equations, surface_integral, dg::DG)
  @unpack interfaces = cache
  @unpack orientations = interfaces

  @threaded for interface in eachinterface(dg, cache)
    left_element  = interfaces.neighbor_ids[1, interface]
    right_element = interfaces.neighbor_ids[2, interface]

    if orientations[interface] == 1
      # interface in x-direction
      for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, j, k, interface] = u[v, nnodes(dg), j, k, left_element]
        interfaces.u[2, v, j, k, interface] = u[v,          1, j, k, right_element]
      end
    elseif orientations[interface] == 2
      # interface in y-direction
      for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, k, interface] = u[v, i, nnodes(dg), k, left_element]
        interfaces.u[2, v, i, k, interface] = u[v, i,          1, k, right_element]
      end
    else # if orientations[interface] == 3
      # interface in z-direction
      for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, j, interface] = u[v, i, j, nnodes(dg), left_element]
        interfaces.u[2, v, i, j, interface] = u[v, i, j,          1, right_element]
      end
    end
  end

  return nothing
end

function calc_interface_flux!(surface_flux_values,
                              mesh::TreeMesh{3},
                              nonconservative_terms::Val{false}, equations,
                              surface_integral, dg::DG, cache)
  @unpack surface_flux = surface_integral
  @unpack u, neighbor_ids, orientations = cache.interfaces

  @threaded for interface in eachinterface(dg, cache)
    # Get neighboring elements
    left_id  = neighbor_ids[1, interface]
    right_id = neighbor_ids[2, interface]

    # Determine interface direction with respect to elements:
    # orientation = 1: left -> 2, right -> 1
    # orientation = 2: left -> 4, right -> 3
    # orientation = 3: left -> 6, right -> 5
    left_direction  = 2 * orientations[interface]
    right_direction = 2 * orientations[interface] - 1

    for j in eachnode(dg), i in eachnode(dg)
      # Call pointwise Riemann solver
      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, j, interface)
      flux = surface_flux(u_ll, u_rr, orientations[interface], equations)

      # Copy flux to left and right element storage
      for v in eachvariable(equations)
        surface_flux_values[v, i, j, left_direction,  left_id]  = flux[v]
        surface_flux_values[v, i, j, right_direction, right_id] = flux[v]
      end
    end
  end
end

function calc_interface_flux!(surface_flux_values,
                              mesh::TreeMesh{3},
                              nonconservative_terms::Val{true}, equations,
                              surface_integral, dg::DG, cache)
  surface_flux, nonconservative_flux = surface_integral.surface_flux
  @unpack u, neighbor_ids, orientations = cache.interfaces

  @threaded for interface in eachinterface(dg, cache)
    # Get neighboring elements
    left_id  = neighbor_ids[1, interface]
    right_id = neighbor_ids[2, interface]

    # Determine interface direction with respect to elements:
    # orientation = 1: left -> 2, right -> 1
    # orientation = 2: left -> 4, right -> 3
    # orientation = 3: left -> 6, right -> 5
    left_direction  = 2 * orientations[interface]
    right_direction = 2 * orientations[interface] - 1

    for j in eachnode(dg), i in eachnode(dg)
      # Call pointwise Riemann solver
      orientation = orientations[interface]
      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, j, interface)
      flux = surface_flux(u_ll, u_rr, orientation, equations)

      # Compute both nonconservative fluxes
      noncons_left  = nonconservative_flux(u_ll, u_rr, orientation, equations)
      noncons_right = nonconservative_flux(u_rr, u_ll, orientation, equations)

      # Copy flux to left and right element storage
      for v in eachvariable(equations)
        # Note the factor 0.5 necessary for the nonconservative fluxes based on
        # the interpretation of global SBP operators coupled discontinuously via
        # central fluxes/SATs
        surface_flux_values[v, i, j, left_direction,  left_id]  = flux[v] + 0.5 * noncons_left[v]
        surface_flux_values[v, i, j, right_direction, right_id] = flux[v] + 0.5 * noncons_right[v]
      end
    end
  end

  return nothing
end

function prolong2boundaries!(cache, u,
                             mesh::TreeMesh{3}, equations, surface_integral, dg::DG)
  @unpack boundaries = cache
  @unpack orientations, neighbor_sides = boundaries

  @threaded for boundary in eachboundary(dg, cache)
    element = boundaries.neighbor_ids[boundary]

    if orientations[boundary] == 1
      # boundary in x-direction
      if neighbor_sides[boundary] == 1
        # element in -x direction of boundary
        for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations)
          boundaries.u[1, v, j, k, boundary] = u[v, nnodes(dg), j, k, element]
        end
      else # Element in +x direction of boundary
        for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations)
          boundaries.u[2, v, j, k, boundary] = u[v, 1,          j, k, element]
        end
      end
    elseif orientations[boundary] == 2
      # boundary in y-direction
      if neighbor_sides[boundary] == 1
        # element in -y direction of boundary
        for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations)
          boundaries.u[1, v, i, k, boundary] = u[v, i, nnodes(dg), k, element]
        end
      else
        # element in +y direction of boundary
        for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations)
          boundaries.u[2, v, i, k, boundary] = u[v, i, 1,          k, element]
        end
      end
    else #if orientations[boundary] == 3
      # boundary in z-direction
      if neighbor_sides[boundary] == 1
        # element in -z direction of boundary
        for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations)
          boundaries.u[1, v, i, j, boundary] = u[v, i, j, nnodes(dg), element]
        end
      else
        # element in +z direction of boundary
        for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations)
          boundaries.u[2, v, i, j, boundary] = u[v, i, j, 1,          element]
        end
      end
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::TreeMesh{3}, equations, surface_integral, dg::DG)
  @assert isempty(eachboundary(dg, cache))
end

function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple,
                             mesh::TreeMesh{3}, equations, surface_integral, dg::DG)
  @unpack surface_flux_values = cache.elements
  @unpack n_boundaries_per_direction = cache.boundaries

  # Calculate indices
  lasts = accumulate(+, n_boundaries_per_direction)
  firsts = lasts - n_boundaries_per_direction .+ 1

  # Calc boundary fluxes in each direction
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1],
                                   equations, surface_integral, dg, cache,
                                   1, firsts[1], lasts[1])
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2],
                                   equations, surface_integral, dg, cache,
                                   2, firsts[2], lasts[2])
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3],
                                   equations, surface_integral, dg, cache,
                                   3, firsts[3], lasts[3])
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4],
                                   equations, surface_integral, dg, cache,
                                   4, firsts[4], lasts[4])
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[5],
                                   equations, surface_integral, dg, cache,
                                   5, firsts[5], lasts[5])
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[6],
                                   equations, surface_integral, dg, cache,
                                   6, firsts[6], lasts[6])
end

function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,5}, t,
                                          boundary_condition, equations,
                                          surface_integral, dg::DG, cache,
                                          direction, first_boundary, last_boundary)
  @unpack surface_flux = surface_integral
  @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries

  @threaded for boundary in first_boundary:last_boundary
    # Get neighboring element
    neighbor = neighbor_ids[boundary]

    for j in eachnode(dg), i in eachnode(dg)
      # Get boundary flux
      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, j, boundary)
      if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right
        u_inner = u_ll
      else # Element is on the right, boundary on the left
        u_inner = u_rr
      end
      x = get_node_coords(node_coordinates, equations, dg, i, j, boundary)
      flux = boundary_condition(u_inner, orientations[boundary], direction, x, t, surface_flux,
                                equations)

      # Copy flux to left and right element storage
      for v in eachvariable(equations)
        surface_flux_values[v, i, j, direction, neighbor] = flux[v]
      end
    end
  end

  return nothing
end

function prolong2mortars!(cache, u,
                          mesh::TreeMesh{3}, equations,
                          mortar_l2::LobattoLegendreMortarL2,
                          surface_integral, dg::DGSEM)
  # temporary buffer for projections
  @unpack fstar_tmp1_threaded = cache

  @threaded for mortar in eachmortar(dg, cache)
    fstar_tmp1 = fstar_tmp1_threaded[Threads.threadid()]

    lower_left_element  = cache.mortars.neighbor_ids[1, mortar]
    lower_right_element = cache.mortars.neighbor_ids[2, mortar]
    upper_left_element  = cache.mortars.neighbor_ids[3, mortar]
    upper_right_element = cache.mortars.neighbor_ids[4, mortar]
    large_element       = cache.mortars.neighbor_ids[5, mortar]

    # Copy solution small to small
    if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
      if cache.mortars.orientations[mortar] == 1
        # L2 mortars in x-direction
        for k in eachnode(dg), j in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u_upper_left[2, v, j, k, mortar]  = u[v, 1, j, k, upper_left_element]
            cache.mortars.u_upper_right[2, v, j, k, mortar] = u[v, 1, j, k, upper_right_element]
            cache.mortars.u_lower_left[2, v, j, k, mortar]  = u[v, 1, j, k, lower_left_element]
            cache.mortars.u_lower_right[2, v, j, k, mortar] = u[v, 1, j, k, lower_right_element]
          end
        end
      elseif cache.mortars.orientations[mortar] == 2
        # L2 mortars in y-direction
        for k in eachnode(dg), i in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u_upper_left[2, v, i, k, mortar]  = u[v, i, 1, k, upper_left_element]
            cache.mortars.u_upper_right[2, v, i, k, mortar] = u[v, i, 1, k, upper_right_element]
            cache.mortars.u_lower_left[2, v, i, k, mortar]  = u[v, i, 1, k, lower_left_element]
            cache.mortars.u_lower_right[2, v, i, k, mortar] = u[v, i, 1, k, lower_right_element]
          end
        end
      else # orientations[mortar] == 3
        # L2 mortars in z-direction
        for j in eachnode(dg), i in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u_upper_left[2, v, i, j, mortar]  = u[v, i, j, 1, upper_left_element]
            cache.mortars.u_upper_right[2, v, i, j, mortar] = u[v, i, j, 1, upper_right_element]
            cache.mortars.u_lower_left[2, v, i, j, mortar]  = u[v, i, j, 1, lower_left_element]
            cache.mortars.u_lower_right[2, v, i, j, mortar] = u[v, i, j, 1, lower_right_element]
          end
        end
      end
    else # large_sides[mortar] == 2 -> small elements on left side
      if cache.mortars.orientations[mortar] == 1
        # L2 mortars in x-direction
        for k in eachnode(dg), j in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u_upper_left[1, v, j, k, mortar]  = u[v, nnodes(dg), j, k, upper_left_element]
            cache.mortars.u_upper_right[1, v, j, k, mortar] = u[v, nnodes(dg), j, k, upper_right_element]
            cache.mortars.u_lower_left[1, v, j, k, mortar]  = u[v, nnodes(dg), j, k, lower_left_element]
            cache.mortars.u_lower_right[1, v, j, k, mortar] = u[v, nnodes(dg), j, k, lower_right_element]
          end
        end
      elseif cache.mortars.orientations[mortar] == 2
        # L2 mortars in y-direction
        for k in eachnode(dg), i in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u_upper_left[1, v, i, k, mortar]  = u[v, i, nnodes(dg), k, upper_left_element]
            cache.mortars.u_upper_right[1, v, i, k, mortar] = u[v, i, nnodes(dg), k, upper_right_element]
            cache.mortars.u_lower_left[1, v, i, k, mortar]  = u[v, i, nnodes(dg), k, lower_left_element]
            cache.mortars.u_lower_right[1, v, i, k, mortar] = u[v, i, nnodes(dg), k, lower_right_element]
          end
        end
      else # if cache.mortars.orientations[mortar] == 3
        # L2 mortars in z-direction
        for j in eachnode(dg), i in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u_upper_left[1, v, i, j, mortar]  = u[v, i, j, nnodes(dg), upper_left_element]
            cache.mortars.u_upper_right[1, v, i, j, mortar] = u[v, i, j, nnodes(dg), upper_right_element]
            cache.mortars.u_lower_left[1, v, i, j, mortar]  = u[v, i, j, nnodes(dg), lower_left_element]
            cache.mortars.u_lower_right[1, v, i, j, mortar] = u[v, i, j, nnodes(dg), lower_right_element]
          end
        end
      end
    end

    # Interpolate large element face data to small interface locations
    if cache.mortars.large_sides[mortar] == 1 # -> large element on left side
      leftright = 1
      if cache.mortars.orientations[mortar] == 1
        # L2 mortars in x-direction
        u_large = view(u, :, nnodes(dg), :, :, large_element)
        element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1)
      elseif cache.mortars.orientations[mortar] == 2
        # L2 mortars in y-direction
        u_large = view(u, :, :, nnodes(dg), :, large_element)
        element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1)
      else # cache.mortars.orientations[mortar] == 3
        # L2 mortars in z-direction
        u_large = view(u, :, :, :, nnodes(dg), large_element)
        element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1)
      end
    else # large_sides[mortar] == 2 -> large element on right side
      leftright = 2
      if cache.mortars.orientations[mortar] == 1
        # L2 mortars in x-direction
        u_large = view(u, :, 1, :, :, large_element)
        element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1)
      elseif cache.mortars.orientations[mortar] == 2
        # L2 mortars in y-direction
        u_large = view(u, :, :, 1, :, large_element)
        element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1)
      else # cache.mortars.orientations[mortar] == 3
        # L2 mortars in z-direction
        u_large = view(u, :, :, :, 1, large_element)
        element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1)
      end
    end
  end

  return nothing
end

@inline function element_solutions_to_mortars!(mortars, mortar_l2::LobattoLegendreMortarL2, leftright, mortar,
                                               u_large::AbstractArray{<:Any,3}, fstar_tmp1)
  multiply_dimensionwise!(view(mortars.u_upper_left,  leftright, :, :, :, mortar), mortar_l2.forward_lower, mortar_l2.forward_upper, u_large, fstar_tmp1)
  multiply_dimensionwise!(view(mortars.u_upper_right, leftright, :, :, :, mortar), mortar_l2.forward_upper, mortar_l2.forward_upper, u_large, fstar_tmp1)
  multiply_dimensionwise!(view(mortars.u_lower_left,  leftright, :, :, :, mortar), mortar_l2.forward_lower, mortar_l2.forward_lower, u_large, fstar_tmp1)
  multiply_dimensionwise!(view(mortars.u_lower_right, leftright, :, :, :, mortar), mortar_l2.forward_upper, mortar_l2.forward_lower, u_large, fstar_tmp1)
  return nothing
end

function calc_mortar_flux!(surface_flux_values,
                           mesh::TreeMesh{3},
                           nonconservative_terms::Val{false}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
                           surface_integral, dg::DG, cache)
  @unpack surface_flux = surface_integral
  @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations = cache.mortars
  @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded,
           fstar_lower_left_threaded, fstar_lower_right_threaded,
           fstar_tmp1_threaded) = cache

  @threaded for mortar in eachmortar(dg, cache)
    # Choose thread-specific pre-allocated container
    fstar_upper_left  = fstar_upper_left_threaded[Threads.threadid()]
    fstar_upper_right = fstar_upper_right_threaded[Threads.threadid()]
    fstar_lower_left  = fstar_lower_left_threaded[Threads.threadid()]
    fstar_lower_right = fstar_lower_right_threaded[Threads.threadid()]
    fstar_tmp1        = fstar_tmp1_threaded[Threads.threadid()]

    # Calculate fluxes
    orientation = orientations[mortar]
    calc_fstar!(fstar_upper_left,  equations, surface_flux, dg, u_upper_left,  mortar, orientation)
    calc_fstar!(fstar_upper_right, equations, surface_flux, dg, u_upper_right, mortar, orientation)
    calc_fstar!(fstar_lower_left,  equations, surface_flux, dg, u_lower_left,  mortar, orientation)
    calc_fstar!(fstar_lower_right, equations, surface_flux, dg, u_lower_right, mortar, orientation)

    mortar_fluxes_to_elements!(surface_flux_values,
                               mesh, equations, mortar_l2, dg, cache, mortar,
                               fstar_upper_left, fstar_upper_right,
                               fstar_lower_left, fstar_lower_right,
                               fstar_tmp1)
  end

  return nothing
end

function calc_mortar_flux!(surface_flux_values,
                           mesh::TreeMesh{3},
                           nonconservative_terms::Val{true}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
                           surface_integral, dg::DG, cache)
  surface_flux, nonconservative_flux = surface_integral.surface_flux
  @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations, large_sides = cache.mortars
  @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded,
           fstar_lower_left_threaded, fstar_lower_right_threaded,
           fstar_tmp1_threaded) = cache

  @threaded for mortar in eachmortar(dg, cache)
    # Choose thread-specific pre-allocated container
    fstar_upper_left  = fstar_upper_left_threaded[Threads.threadid()]
    fstar_upper_right = fstar_upper_right_threaded[Threads.threadid()]
    fstar_lower_left  = fstar_lower_left_threaded[Threads.threadid()]
    fstar_lower_right = fstar_lower_right_threaded[Threads.threadid()]
    fstar_tmp1        = fstar_tmp1_threaded[Threads.threadid()]

    # Calculate fluxes
    orientation = orientations[mortar]
    calc_fstar!(fstar_upper_left,  equations, surface_flux, dg, u_upper_left,  mortar, orientation)
    calc_fstar!(fstar_upper_right, equations, surface_flux, dg, u_upper_right, mortar, orientation)
    calc_fstar!(fstar_lower_left,  equations, surface_flux, dg, u_lower_left,  mortar, orientation)
    calc_fstar!(fstar_lower_right, equations, surface_flux, dg, u_lower_right, mortar, orientation)

    # Add nonconservative fluxes.
    # These need to be adapted on the geometry (left/right) since the order of
    # the arguments matters, based on the global SBP operator intepretation.
    # The same interpretation (global SBP operators coupled discontinuously via
    # central fluxes/SATs) explains why we need the factor 0.5.
    # Alternatively, you can also follow the argumentation of Bohm et al. 2018
    # ("nonconservative diamond flux")
    if large_sides[mortar] == 1 # -> small elements on right side
      for j in eachnode(dg), i in eachnode(dg)
        # Pull the left and right solutions
        u_upper_left_ll,  u_upper_left_rr  = get_surface_node_vars(u_upper_left,  equations, dg, i, j, mortar)
        u_upper_right_ll, u_upper_right_rr = get_surface_node_vars(u_upper_right, equations, dg, i, j, mortar)
        u_lower_left_ll,  u_lower_left_rr  = get_surface_node_vars(u_lower_left,  equations, dg, i, j, mortar)
        u_lower_right_ll, u_lower_right_rr = get_surface_node_vars(u_lower_right, equations, dg, i, j, mortar)
        # Call pointwise nonconservative term
        noncons_upper_left  = nonconservative_flux(u_upper_left_ll,  u_upper_left_rr,  orientation, equations)
        noncons_upper_right = nonconservative_flux(u_upper_right_ll, u_upper_right_rr, orientation, equations)
        noncons_lower_left  = nonconservative_flux(u_lower_left_ll,  u_lower_left_rr,  orientation, equations)
        noncons_lower_right = nonconservative_flux(u_lower_right_ll, u_lower_right_rr, orientation, equations)
        # Add to primary and secondary temporay storage
        multiply_add_to_node_vars!(fstar_upper_left,  0.5, noncons_upper_left,  equations, dg, i, j)
        multiply_add_to_node_vars!(fstar_upper_right, 0.5, noncons_upper_right, equations, dg, i, j)
        multiply_add_to_node_vars!(fstar_lower_left,  0.5, noncons_lower_left,  equations, dg, i, j)
        multiply_add_to_node_vars!(fstar_lower_right, 0.5, noncons_lower_right, equations, dg, i, j)
      end
    else # large_sides[mortar] == 2 -> small elements on the left
      for j in eachnode(dg), i in eachnode(dg)
        # Pull the left and right solutions
        u_upper_left_ll,  u_upper_left_rr  = get_surface_node_vars(u_upper_left,  equations, dg, i, j, mortar)
        u_upper_right_ll, u_upper_right_rr = get_surface_node_vars(u_upper_right, equations, dg, i, j, mortar)
        u_lower_left_ll,  u_lower_left_rr  = get_surface_node_vars(u_lower_left,  equations, dg, i, j, mortar)
        u_lower_right_ll, u_lower_right_rr = get_surface_node_vars(u_lower_right, equations, dg, i, j, mortar)
        # Call pointwise nonconservative term
        noncons_upper_left  = nonconservative_flux(u_upper_left_rr,  u_upper_left_ll,  orientation, equations)
        noncons_upper_right = nonconservative_flux(u_upper_right_rr, u_upper_right_ll, orientation, equations)
        noncons_lower_left  = nonconservative_flux(u_lower_left_rr,  u_lower_left_ll,  orientation, equations)
        noncons_lower_right = nonconservative_flux(u_lower_right_rr, u_lower_right_ll, orientation, equations)
        # Add to primary and secondary temporay storage
        multiply_add_to_node_vars!(fstar_upper_left,  0.5, noncons_upper_left,  equations, dg, i, j)
        multiply_add_to_node_vars!(fstar_upper_right, 0.5, noncons_upper_right, equations, dg, i, j)
        multiply_add_to_node_vars!(fstar_lower_left,  0.5, noncons_lower_left,  equations, dg, i, j)
        multiply_add_to_node_vars!(fstar_lower_right, 0.5, noncons_lower_right, equations, dg, i, j)
      end
    end

    mortar_fluxes_to_elements!(surface_flux_values,
                               mesh, equations, mortar_l2, dg, cache, mortar,
                               fstar_upper_left, fstar_upper_right,
                               fstar_lower_left, fstar_lower_right,
                               fstar_tmp1)
  end

  return nothing
end

@inline function calc_fstar!(destination::AbstractArray{<:Any,3}, equations,
                             surface_flux, dg::DGSEM,
                             u_interfaces, interface, orientation)

  for j in eachnode(dg), i in eachnode(dg)
    # Call pointwise two-point numerical flux function
    u_ll, u_rr = get_surface_node_vars(u_interfaces, equations, dg, i, j, interface)
    flux = surface_flux(u_ll, u_rr, orientation, equations)

    # Copy flux to left and right element storage
    set_node_vars!(destination, flux, equations, dg, i, j)
  end

  return nothing
end

@inline function mortar_fluxes_to_elements!(surface_flux_values,
                                            mesh::TreeMesh{3}, equations,
                                            mortar_l2::LobattoLegendreMortarL2,
                                            dg::DGSEM, cache,
                                            mortar,
                                            fstar_upper_left, fstar_upper_right,
                                            fstar_lower_left, fstar_lower_right,
                                            fstar_tmp1)
  lower_left_element  = cache.mortars.neighbor_ids[1, mortar]
  lower_right_element = cache.mortars.neighbor_ids[2, mortar]
  upper_left_element  = cache.mortars.neighbor_ids[3, mortar]
  upper_right_element = cache.mortars.neighbor_ids[4, mortar]
  large_element       = cache.mortars.neighbor_ids[5, mortar]

  # Copy flux small to small
  if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
    if cache.mortars.orientations[mortar] == 1
      # L2 mortars in x-direction
      direction = 1
    elseif cache.mortars.orientations[mortar] == 2
      # L2 mortars in y-direction
      direction = 3
    else # if cache.mortars.orientations[mortar] == 3
      # L2 mortars in z-direction
      direction = 5
    end
  else # large_sides[mortar] == 2 -> small elements on left side
    if cache.mortars.orientations[mortar] == 1
      # L2 mortars in x-direction
      direction = 2
    elseif cache.mortars.orientations[mortar] == 2
      # L2 mortars in y-direction
      direction = 4
    else # if cache.mortars.orientations[mortar] == 3
      # L2 mortars in z-direction
      direction = 6
    end
  end
  surface_flux_values[:, :, :, direction, upper_left_element]  .= fstar_upper_left
  surface_flux_values[:, :, :, direction, upper_right_element] .= fstar_upper_right
  surface_flux_values[:, :, :, direction, lower_left_element]  .= fstar_lower_left
  surface_flux_values[:, :, :, direction, lower_right_element] .= fstar_lower_right

  # Project small fluxes to large element
  if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
    if cache.mortars.orientations[mortar] == 1
      # L2 mortars in x-direction
      direction = 2
    elseif cache.mortars.orientations[mortar] == 2
      # L2 mortars in y-direction
      direction = 4
    else # if cache.mortars.orientations[mortar] == 3
      # L2 mortars in z-direction
      direction = 6
    end
  else # large_sides[mortar] == 2 -> small elements on left side
    if cache.mortars.orientations[mortar] == 1
      # L2 mortars in x-direction
      direction = 1
    elseif cache.mortars.orientations[mortar] == 2
      # L2 mortars in y-direction
      direction = 3
    else # if cache.mortars.orientations[mortar] == 3
      # L2 mortars in z-direction
      direction = 5
    end
  end

  multiply_dimensionwise!(
    view(surface_flux_values, :, :, :, direction, large_element),
    mortar_l2.reverse_lower, mortar_l2.reverse_upper, fstar_upper_left, fstar_tmp1)
  add_multiply_dimensionwise!(
    view(surface_flux_values, :, :, :, direction, large_element),
    mortar_l2.reverse_upper, mortar_l2.reverse_upper, fstar_upper_right, fstar_tmp1)
  add_multiply_dimensionwise!(
    view(surface_flux_values, :, :, :, direction, large_element),
    mortar_l2.reverse_lower, mortar_l2.reverse_lower, fstar_lower_left, fstar_tmp1)
  add_multiply_dimensionwise!(
    view(surface_flux_values, :, :, :, direction, large_element),
    mortar_l2.reverse_upper, mortar_l2.reverse_lower, fstar_lower_right, fstar_tmp1)

  return nothing
end

function calc_surface_integral!(du, u, mesh::Union{TreeMesh{3}, StructuredMesh{3}},
                                equations, surface_integral, dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  @threaded for element in eachelement(dg, cache)
    for m in eachnode(dg), l in eachnode(dg)
      for v in eachvariable(equations)
        # surface at -x
        du[v, 1,          l, m, element] -= surface_flux_values[v, l, m, 1, element] * boundary_interpolation[1,          1]
        # surface at +x
        du[v, nnodes(dg), l, m, element] += surface_flux_values[v, l, m, 2, element] * boundary_interpolation[nnodes(dg), 2]
        # surface at -y
        du[v, l, 1,          m, element] -= surface_flux_values[v, l, m, 3, element] * boundary_interpolation[1,          1]
        # surface at +y
        du[v, l, nnodes(dg), m, element] += surface_flux_values[v, l, m, 4, element] * boundary_interpolation[nnodes(dg), 2]
        # surface at -z
        du[v, l, m, 1,          element] -= surface_flux_values[v, l, m, 5, element] * boundary_interpolation[1,          1]
        # surface at +z
        du[v, l, m, nnodes(dg), element] += surface_flux_values[v, l, m, 6, element] * boundary_interpolation[nnodes(dg), 2]
      end
    end
  end

  return nothing
end

function apply_jacobian!(du, mesh::TreeMesh{3},
                         equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    factor = -cache.elements.inverse_jacobian[element]

    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      for v in eachvariable(equations)
        du[v, i, j, k, element] *= factor
      end
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_sources!(du, u, t, source_terms::Nothing,
                       equations::AbstractEquations{3}, dg::DG, cache)
  return nothing
end

function calc_sources!(du, u, t, source_terms,
                       equations::AbstractEquations{3}, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_local = get_node_vars(u, equations, dg, i, j, k, element)
      x_local = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, k, element)
      du_local = source_terms(u_local, x_local, t, equations)
      add_to_node_vars!(du, du_local, equations, dg, i, j, k, element)
    end
  end

  return nothing
end

end # @muladd


text:
include( "../src/bootstrapPosteriorEstimation.jl" )
using DataFrames
using Random
using Statistics
using Test
using Optim
using MultivariateStats

println( "running bootstrap algorithm unit tests." );
predMatTraining, predMatEval, yTraining, yEval, errMatTraining = makeupPredictions();
posterior = bootstrapPosteriorEstimation( errMatTraining, 64, 10000 );
@test all( posterior .>= 0.0 );
@test all( posterior .<= 1.0 );
@test all( sum( posterior ) ≈ 1.0 );

posterior = zeros( Float64, size( predMatTraining, 2 ) );
bootstrapPosteriorEstimation!( errMatTraining, 100, 10000, posterior );
@test all( posterior .>= 0.0 );
@test all( posterior .<= 1.0 );
@test all( sum( posterior ) ≈ 1.0 );

ext: jl


text:
# phylo/phylogeny.jl
# ==================
#
# Types and methods for phylogenetic trees.
#
# Part of the Bio.Phylo module.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# The Phylogeny Type
# -------------------

# The internal structure and representation of the Phylogeny type uses a DiGraph.
# The vertices of the DiGraph are reserved for different specific purposes.

# Vertex ranges:
# * 1:tree.ntaxa => leaf vertices
# * tree.ntaxa + 1 => root vertex
# * tree.ntaxa + 1:<#nodes> => clade vertices

"""
Phylogeny represents a phylogenetic tree.

The type is parametric with two parameters `C` and `B`.

This is because it is common to want to annotate tips, clades, and branches
in a phylogeny with data to create a richer model of evolution of do other things
like dictate aesthetic values when plotting.

Type parameter `C` dictates what datatype can be stored in the phylogeny to annotate
clades and tips. Type parameter `B` dictates what datatype can be stored in the
phylogeny to annotate branches. Think `C` for clades and `B` for branches.
"""
type Phylogeny{C, B}
    graph::DiGraph
    vertexindex::Indexer{Int}
    vertexdata::Vector{C}
    edgedata::Dict{Edge, B}
    ntaxa::Int
    rooted::Bool
    rerootable::Bool
end

# Phylogeny Constructors
# -----------------------

"""
    Phylogeny{C, B}(tips::Vector{Symbol}, ::Type{C}, ::Type{B})

Construct a completely unresolved phylogeny, with no annotated data.

The phylogeny will store annotation data of type `C` for tree nodes
(clades and tips). It will also store data of type `B` associated with branches.

# Examples
```julia
Phylogeny(Void, Void, [:Human, :Chimp, :Dog, :Cat])
Phylogeny(Float64, Float64, [:Human, :Chimp, :Dog, :Cat])
```
"""
function Phylogeny{C, B}(::Type{C}, ::Type{B}, taxa::Vector{Symbol})
    ntaxa = length(taxa)
    maxVertices = maxvertices(ntaxa)
    taxaAndRoot = push!(copy(taxa), :Root)
    idxer = Indexer(taxaAndRoot, Int)
    g = DiGraph(maxVertices)
    for i in 1:ntaxa
        add_edge!(g, ntaxa + 1, i)
    end
    v = Vector{C}()
    d = Dict{Edge, B}()
    return Phylogeny{C, B}(g, idxer, v, d, ntaxa, false, true)
end

"""
    Phylogeny(taxa::Vector{Symbol})

Construct a completely unresolved phylogeny, with no annotated data.

This constructor assumes you want a phylogeny where the only data associated
with it is branch lengths (`Float64`) for the edges, and confidence values
(`Float64`) for the clades.

# Examples
```julia
Phylogeny([:Human, :Chimp, :Dog, :Cat])
```
"""
function Phylogeny(taxa::Vector{Symbol})
    return Phylogeny(Float64, Float64, taxa)
end

"""
    Phylogeny{T <: AbstractString}(taxa::Vector{T})
Construct a completely unresolved phylogeny, with no annotated data.

This constructor assumes you want a phylogeny where the only data associated
with it is branch lengths (`Float64`) for the edges, and confidence values
(`Float64`) for the clades.
"""
function Phylogeny{T <: AbstractString}(taxa::Vector{T})
    return Phylogeny(convert(Vector{Symbol}, taxa))
end

function Base.copy{C, B}(tree::Phylogeny{C, B})
    g = copy(tree.graph)
    i = copy(tree.vertexindex)
    v = copy(tree.vertexdata)
    d = copy(tree.edgedata)
    return Phylogeny{C, B}(g, i, v, d, tree.ntaxa, tree.rooted, tree.rerootable)
end

function Base.:(==)(x::Phylogeny, y::Phylogeny)
    g = x.graph == x.graph
    vi = x.vertexindex == y.vertexindex
    vd = x.vertexdata == y.vertexdata
    ed = x.edgedata == y.edgedata
    t = x.ntaxa == y.ntaxa
    r = x.rooted == y.rooted
    rr = x.rerootable == y.rerootable
    return g && vi && vd && ed && t && r && rr
end

# Miscellaneous Functions
# ------------------------

"""
    n_possible_rooted(n::Int)

Returns the number of possible rooted tree topologies for n taxa.

**Warning:** This method currently overflows when `n` > 12.
"""
function n_possible_rooted(n::Int)
    @assert n >= 2 DomainError("n must be greater or equal to 2.")
    num = factorial(2 * n - 3)
    den = (2^(n - 2)) * factorial(n - 2)
    return Int(num / den)
end

"""
    n_possible_unrooted(n::Int)

Returns the number of possible rooted tree topologies for n taxa.

**Warning:** This method currently overflows when `n` > 13.
"""
function n_possible_unrooted(n::Int)
    @assert n >= 3 DomainError("n must be greater of equal to 3.")
    return n_possible_rooted(n - 1)
end

maxvertices(n::Int) = 2 * n - 1

typealias SymString Union{Symbol, AbstractString}

"""
    isrooted(x::Phylogeny)

Test whether a Phylogeny is rooted.

# Examples
```julia
isrooted(my_phylogeny)
```
"""
function isrooted(x::Phylogeny)
    return x.rooted
end

"""
    isrerootable(x::Phylogeny)

Test whether a Phylogeny is re-rootable.

# Examples
```julia
isrerootable(my_phylogeny)
```
"""
function isrerootable(x::Phylogeny)
    return x.rerootable
end

"Get the vertex of the tree which represents the root of the tree."
root(tree::Phylogeny) = tree.ntaxa + 1

"Get the range of vertices which represent the internal clades of the tree."
clades(tree::Phylogeny) = root(tree):nv(tree.graph)

"Get the range of vertices which represent the leaves of the tree."
leaves(tree::Phylogeny) = 1:tree.ntaxa

"""
    hasparent(tree::Phylogeny)

Check every vertex in a phylogeny and test if each is linked to a parent.
"""
function hasparent(tree::Phylogeny)
    return indegree(tree.graph) .>= 1
end

"""
    haschildren(tree::Phylogeny)

Test each vertex in a phylogeny to see if it has children.
"""
function haschildren(tree::Phylogeny)
    return nchildren(tree) .> 0
end

"""
    nchildren(tree::Phylogeny)

Count the number of clades that are direct children of each clade.
"""
function nchildren(tree::Phylogeny)
    return outdegree(tree.graph)
end

"""
    isredundant(tree::Phylogeny)

Test each vertex in a phylogeny and test if it redundant.

Vertices are considered redundant if they have a parent and 1 child.
"""
function isredundant(tree::Phylogeny)
    return hasparent(tree) & nchildren(tree) .== 1
end


text:
using Librsync
if VERSION < v"0.7-"
    using Compat: equalto
end
function findPoptArg(ARGVEC, shortArg::String, longArg::String, default = "")
    longArgStr = string("--", longArg, "=")
    nLongArgStr = length(longArgStr)
    for (i, arg) in enumerate(ARGVEC)
        if arg == string("-", shortArg)
            return get(ARGVEC, i + 1, default)
        elseif startswith(arg, longArgStr)
            return arg[nLongArgStr:end]
        end
    end
    return default
end
function findPoptPosArgTwo(ARGVEC, sigpos::Int)
    @inbounds if sigpos == 0 || length(ARGVEC) == sigpos
        return "", ""
    elseif length(ARGVEC) == sigpos + 1
        return ARGVEC[sigpos + 1], ""
    elseif length(ARGVEC) > sigpos + 1
        return ARGVEC[sigpos + 1], ARGVEC[sigpos + 2]
    end

    return "", ""
end
function findPoptPosArgThree(ARGVEC, sigpos)
    #TODO: Add a Val{} argument and generate this function (just for fun...)
    @inbounds if sigpos == 0 || length(ARGVEC) == sigpos
        return "", "", ""
    elseif length(ARGVEC) == sigpos + 1
        return ARGVEC[sigpos + 1], "", ""
    elseif length(ARGVEC) == sigpos + 2
        return ARGVEC[sigpos + 1], ARGVEC[sigpos + 2], ""
    elseif length(ARGVEC) > sigpos + 2
        return ARGVEC[sigpos + 1], ARGVEC[sigpos + 2], ARGVEC[sigpos + 3]
    end

    return "", ""
end

findHash(ARGVEC) = findPoptArg(ARGVEC, "H", "hash", "blake2")
function findBlockSize(ARGVEC)::Int
    #TODO: Figure out how to make this usage of parse inferrable?
    # Test cases to be added:
    # julia -e 'using Librsync; Librsync.main()' -- -b "Testing." -H "blake"
    # julia -e 'using Librsync; Librsync.main()' -- -b "" -H "blake"
    # julia -e 'using Librsync; Librsync.main()' -- -b -H "blake"
    # julia -e 'using Librsync; Librsync.main()' -- -b 1000 -H "blake"
    rawArg = findPoptArg(ARGVEC, "b", "block-size")
    try
        return parse(Int, rawArg)
    catch
        return Librsync.RS_DEFAULT_BLOCK_LEN
    end
end
function findSumSize(ARGVEC)::Int
    #TODO: Figure out how to make this usage of parse inferrable?
    rawArg = findPoptArg(ARGVEC, "S", "sum-size")
    try
        return parse(Int, rawArg)
    catch
        return 0
    end
end
const defaultSumSize = 0
function findSumSize(ARGVEC)::Int
    #TODO: Figure out how to make this usage of parse inferrable?
    rawArg = findPoptArg(ARGVEC, "S", "sum-size")
    try
        return parse(Int, rawArg)
    catch
        return defaultSumSize
    end
end
function findInputSize(ARGVEC)::Int
    #TODO: Figure out how to make this usage of parse inferrable?
    rawArg = findPoptArg(ARGVEC, "I", "input-size")
    try
        return parse(Int, rawArg)
    catch
        return Librsync.inbuflen
    end
end
function findOutputSize(ARGVEC)::Int
    #TODO: Figure out how to make this usage of parse inferrable?
    rawArg = findPoptArg(ARGVEC, "O", "output-size")
    try
        return parse(Int, rawArg)
    catch
        return Librsync.outbuflen
    end
end
findStats(ARGVEC) = in("-s", ARGVEC) || in("--statistics", ARGVEC)
findParanoia(ARGVEC) = in("--paranoia", ARGVEC)

function parseRdiffArgs(ARGVEC)
    hashMethod = findHash(ARGVEC)
    blockSize = findBlockSize(ARGVEC)
    sumSize = findSumSize(ARGVEC)
    inputSize = findInputSize(ARGVEC)
    outputSize = findOutputSize(ARGVEC)
    stats = findStats(ARGVEC)
    paranoia = findParanoia(ARGVEC)

    hashMethod, blockSize, sumSize, inputSize, stats, paranoia
end
function findRdiffMode(ARGVEC)
    for mode in ("signature", "delta", "patch")
        sigpos = findfirst(equalto(mode), ARGVEC)
        if sigpos > 0
            return mode, sigpos
        end
    end
    return "", 0
end
function generate_sig(basis::IO, sig::IO, hashMethod, blockSize, sumSize, inputSize, stats, paranoia)
    sig_magic = Librsync.RS_BLAKE2_SIG_MAGIC
    if hashMethod == "blake2"
        sig_magic = Librsync.RS_BLAKE2_SIG_MAGIC
    elseif hashMethod == "md4"
        if sumSize == defaultSumSize
            sumSize = 8
        end
        sig_magic = Librsync.RS_MD4_SIG_MAGIC
    else
        warn("Unknown hash algorithm $(hashMethod)")
        return -1
    end

    statStruct = Ref(Librsync.construct_rs_stats_t())

    result = Librsync.rs_sig_file(basis, sig, blockSize, sumSize, sig_magic, statStruct)

    if result != Librsync.RS_DONE
        return result
    end
    if stats
        Librsync.rs_log_stats(statStruct)
    end

    return result
end

function usage()
println("""
Usage: julia rdiff.jl [OPTIONS] signature [BASIS [SIGNATURE]]
                      [OPTIONS] delta SIGNATURE [NEWFILE [DELTA]]
                      [OPTIONS] patch BASIS [DELTA [NEWFILE]]

    Options:
        -?, --help                Show this help message
        -s, --statistics          Show performance statistics
    Signature generation options:
        -H, --hash=ALG            Hash algorithm: blake2 (default), md4
    Delta-encoding options:
        -b, --block-size=BYTES    Signature block size
        -S, --sum-size=BYTES      Set signature strength
            --paranoia            Verify all rolling checksums
    IO options:
        -I, --input-size=BYTES    Input buffer size
        -O, --output-size=BYTES   Output buffer size
""")
end
function main(ARGVEC)
    mode, sigpos = findRdiffMode(ARGVEC)
    if in("-h", ARGVEC) || in("-?", ARGVEC) || in("--help", ARGVEC)
        usage()
        exit(-1)
    end
    hashMethod, blockSize, sumSize, inputSize, stats, paranoia = parseRdiffArgs(ARGVEC)
    if mode == "signature"
        basis_file, sig_file = findPoptPosArgTwo(ARGVEC, sigpos)
        if basis_file == ""
            return -1
        end
        retval = -1
        if sig_file == ""
            open(basis_file, "r") do basis_io
                retval = generate_sig(basis_io, STDOUT, hashMethod, blockSize, sumSize, inputSize, stats, paranoia)
            end
        else
            open(sig_file, "w") do sig_io
                open(basis_file, "r") do basis_io
                    retval = generate_sig(basis_io, sig_io, hashMethod, blockSize, sumSize, inputSize, stats, paranoia)
                end
            end
        end

        return retval
    elseif mode == "delta"
        @show "delta unimplemented"
        @show findPoptPosArgThree(ARGVEC, sigpos)
        return 0
    elseif mode == "patch"
        @show "patch unimplemented"
        @show findPoptPosArgThree(ARGVEC, sigpos)
        return 0
    else
        usage()
        return -1
    end
end
exit(main(ARGS))


text:
import ChainRulesCore:
    rrule, @non_differentiable, NoTangent, Tangent, backing, AbstractTangent, ZeroTangent

function create_circuit_tangent(circuit, params)
    gc = dispatch(circuit, params)
    res = recursive_create_tangent(gc)
    return res
end

# fallback
function recursive_create_tangent(c::AbstractBlock)
    if nparameters(c) == 0
        return NoTangent()
    else
        error("`Tangent` for type $(typeof(c)) is not defined!")
    end
end
# primitive blocks
unsafe_primitive_tangent(::Any) = NoTangent()
unsafe_primitive_tangent(x::Number) = x
for GT in [:RotationGate, :ShiftGate, :PhaseGate, :(Scale{<:Number})]
    @eval function recursive_create_tangent(c::$GT)
        lst = map(fieldnames(typeof(c))) do fn
            fn => unsafe_primitive_tangent(getfield(c, fn))
        end
        nt = NamedTuple(lst)
        Tangent{typeof(c),typeof(nt)}(nt)
    end
end
function recursive_create_tangent(c::TimeEvolution)
    Tangent{typeof(c)}(; H=NoTangent(), dt=c.dt, tol=NoTangent())
end

# composite blocks
unsafe_composite_tangent(::Any) = NoTangent()
unsafe_composite_tangent(c::AbstractVector{<:AbstractBlock}) = recursive_create_tangent.(c)
unsafe_composite_tangent(c::AbstractBlock) = recursive_create_tangent(c)
for GT in [
    :ChainBlock,
    :Add,
    :KronBlock,
    :RepeatedBlock,
    :PutBlock,
    :ControlBlock,
    :Subroutine,
    :CachedBlock,
    :Daggered,
    :Scale,
]
    @eval function recursive_create_tangent(c::$GT)
        lst = map(fieldnames(typeof(c))) do fn
            fn => unsafe_composite_tangent(getfield(c, fn))
        end
        nt = NamedTuple(lst)
        Tangent{typeof(c),typeof(nt)}(nt)
    end
end

extract_circuit_gradients!(c::Number, output) = push!(output, c)
extract_circuit_gradients!(::Nothing, output) = output
extract_circuit_gradients!(::NoTangent, output) = output
extract_circuit_gradients!(::ZeroTangent, output) = output
function extract_circuit_gradients!(c::AbstractVector, output)
    for ci in c
        extract_circuit_gradients!(ci, output)
    end
    return output
end
function extract_circuit_gradients!(c::Tangent, output)
    for fn in getfield(c, :backing)
        extract_circuit_gradients!(fn, output)
    end
    return output
end
function extract_circuit_gradients!(c::NamedTuple, output)
    for fn in c
        extract_circuit_gradients!(fn, output)
    end
    return output
end

function rrule(::typeof(apply), reg::AbstractArrayReg, block::AbstractBlock)
    out = apply(reg, block)
    out, function (outδ)
        (in, inδ), paramsδ = apply_back((copy(out), outδ), block)
        return (NoTangent(), inδ, create_circuit_tangent(block, paramsδ))
    end
end
function rrule(::typeof(apply), reg::AbstractArrayReg, block::Add)
    out = apply(reg, block)
    out, function (outδ)
        (in, inδ), paramsδ = apply_back((copy(out), outδ), block; in = reg)
        return (NoTangent(), inδ, create_circuit_tangent(block, paramsδ))
    end
end

function rrule(::typeof(dispatch), block::AbstractBlock, params)
    out = dispatch(block, params)
    out, function (outδ::AbstractTangent)
        g = extract_circuit_gradients!(outδ, empty(params))
        res = (NoTangent(), NoTangent(), g)
        return res
    end
end

function rrule(::typeof(expect), op::AbstractBlock, reg::AbstractArrayReg)
    B = YaoArrayRegister._asint(nbatch(reg))
    out = expect(op, reg)
    out, function (outδ)
        greg = expect_g(op, reg)
        for b = 1:B
            viewbatch(greg, b).state .*= 2 * outδ[b]
        end
        return (NoTangent(), NoTangent(), greg)
    end
end

function rrule(
    ::typeof(expect),
    op::AbstractBlock,
    reg_and_circuit::Pair{<:AbstractArrayReg,<:AbstractBlock})
    out = expect(op, reg_and_circuit)
    out,
    function (outδ)
        reg, c = reg_and_circuit
        out = copy(reg) |> c
        goutreg = copy(out) |> op
        for b = 1:YaoArrayRegister._asint(nbatch(goutreg))
            viewbatch(goutreg, b).state .*= 2 * outδ[b]
        end
        # apply backward rule
        (in, greg), gcircuit = apply_back((out, goutreg), c)
        return (
            NoTangent(),
            NoTangent(),
            Tangent{typeof(reg_and_circuit)}(;
                first = greg,
                second = create_circuit_tangent(reg_and_circuit.second, gcircuit),
            ),
        )
    end
end

function rrule(::Type{T}, block::AbstractBlock) where {T<:Matrix}
    out = T(block)
    out, function (outδ)
        paramsδ = mat_back(block, outδ)
        return (NoTangent(), create_circuit_tangent(block, paramsδ))
    end
end

function rrule(::typeof(mat), ::Type{T}, block::AbstractBlock) where {T}
    out = mat(T, block)
    out, function (outδ)
        paramsδ = mat_back(block, outδ)
        return (NoTangent(), NoTangent(), create_circuit_tangent(block, paramsδ))
    end
end

function rrule(::Type{ArrayReg{D}}, raw::AbstractArray) where {D}
    ArrayReg{D}(raw), adjy -> (NoTangent(), reshape(adjy.state, size(raw)))
end

function rrule(::Type{ArrayReg}, raw::AbstractArray)
    ArrayReg(raw), adjy -> (NoTangent(), reshape(adjy.state, size(raw)))
end

function rrule(::Type{BatchedArrayReg{D}}, raw::AbstractArray, nbatch::Int) where {D}
    BatchedArrayReg{D}(raw, nbatch), adjy -> (NoTangent(), reshape(adjy.state, size(raw)), NoTangent())
end

function rrule(::Type{BatchedArrayReg}, raw::AbstractArray, nbatch::Int)
    BatchedArrayReg(raw, nbatch), adjy -> (NoTangent(), reshape(adjy.state, size(raw)), NoTangent())
end

function rrule(::typeof(copy), reg::AbstractArrayReg)
    copy(reg), adjy -> (NoTangent(), adjy)
end

for (BT, BLOCKS) in [(:Add, :(outδ.list)) (:ChainBlock, :(outδ.blocks))]
    for ST in [:AbstractVector, :Tuple]
        @eval function rrule(::Type{BT}, source::$ST) where {BT<:$BT}
            out = BT(source)
            out, function (outδ)
                return (NoTangent(), $ST($BLOCKS))
            end
        end
    end
    @eval function rrule(::Type{BT}, args::AbstractBlock...) where {BT<:$BT}
        out = BT(args...)
        out, function (outδ)
            return (NoTangent(), $BLOCKS...)
        end
    end
end

rrule(::typeof(state), reg::AbstractArrayReg{D,T}) where {D,T} =
    state(reg), adjy -> (NoTangent(), _match_type(reg, _totype(T, adjy)))
rrule(::typeof(statevec), reg::AbstractArrayReg{D,T}) where {D,T} =
    statevec(reg), adjy -> (NoTangent(), _match_type(reg, _totype(T, adjy)))
rrule(::typeof(state), reg::AdjointArrayReg{D,T}) where {D,T} =
    state(reg), adjy -> (NoTangent(), _match_type(parent(reg), _totype(T, adjy)')')
rrule(::typeof(statevec), reg::AdjointArrayReg{D,T}) where {D,T} =
    statevec(reg), adjy -> (NoTangent(), _match_type(parent(reg), _totype(T, adjy)')')
rrule(::typeof(parent), reg::AdjointArrayReg) = parent(reg), adjy -> (NoTangent(), adjy')
rrule(::typeof(Base.adjoint), reg::AbstractArrayReg) =
    Base.adjoint(reg), adjy -> (NoTangent(), parent(adjy))

_totype(::Type{T}, x::AbstractArray{T}) where {T} = x
_totype(::Type{T}, x::AbstractArray{T2}) where {T,T2} = convert.(T, x)
_match_type(::ArrayReg{D}, mat) where D = ArrayReg{D}(mat)
_match_type(r::BatchedArrayReg{D}, mat) where D = BatchedArrayReg{D}(mat, r.nbatch)

@non_differentiable nparameters(::Any)
@non_differentiable zero_state(args...)
@non_differentiable rand_state(args...)
@non_differentiable uniform_state(args...)
@non_differentiable product_state(args...)


text:
"""
sim!(output, ruleset; init=nothing, tstpan=(1, length(output)),
     fps=fps(output), data=nothing, nreplicates=nothing)

Runs the whole simulation, passing the destination aray to
the passed in output for each time-step.

### Arguments
- `output`: An [AbstractOutput](@ref) to store frames or display them on the screen.
- `ruleset`: A Rule() containing one ore more [`AbstractRule`](@ref). These will each be run in sequence.

### Keyword Arguments
- `init`: the initialisation array. If `nothing`, the Ruleset must contain an `init` array.
- `tspan`: the timespan simulaiton will run for.
- `fps`: the frames per second to display. Will be taken from the output if not passed in.
- `nreplicates`: the number of replicates to combine in stochastic simulations
- `data`: a SimData object. Can reduce allocations when that is important.
"""
sim!(output, ruleset; init=nothing, tspan=(1, length(output)), fps=fps(output),
     nreplicates=nothing, data=nothing) = begin
    isrunning(output) && error("A simulation is already running in this output")
    setrunning!(output, true) || error("Could not start the simulation with this output")
    starttime = first(tspan)
    fspan = tspan2fspan(tspan, timestep(ruleset))
    setstarttime!(output, starttime)
    # Copy the init array from the ruleset or keyword arg
    init = chooseinit(DynamicGrids.init(ruleset), init)
    data = initdata!(data, ruleset, init, starttime, nreplicates)
    # Delete frames output by the previous simulations
    initframes!(output, init)
    setfps!(output, fps)
    # Show the first frame
    showframe(output, data, 1)
    # Let the init frame be displayed as long as a normal frame
    delay(output, 1)
    # Run the simulation
    runsim!(output, data, fspan)
end

tspan2fspan(tspan, tstep) = 1:lastindex(first(tspan):tstep:last(tspan))

# Allows attaching an init array to the ruleset, but also passing in an
# alternate array as a keyword arg (which will take preference).
chooseinit(rulesetinit, arginit) = arginit
chooseinit(rulesetinit::Nothing, arginit) = arginit
chooseinit(rulesetinit, arginit::Nothing) = rulesetinit
chooseinit(rulesetinit::Nothing, arginit::Nothing) =
    error("Include an init array: either in the ruleset or with the `init` keyword")

"""
    resume!(output, ruleset; tadd=100, kwargs...)

Restart the simulation where you stopped last time. For arguments see [`sim!`](@ref).
The keyword arg `tadd` indicates the number of frames to add, and of course an init
array will not be accepted.
"""
resume!(output, ruleset; tstop=stoptime(output), fps=fps(output), data=nothing,
        nreplicates=nothing) = begin
    length(output) > 0 || error("There is no simulation to resume. Run `sim!` first")
    isrunning(output) && error("A simulation is already running in this output")
    setrunning!(output, true) || error("Could not start the simulation with this output")
    tstart = starttime(output)
    lastframe = lastindex(tstart:timestep(ruleset):stoptime(output))
    stopframe = lastindex(tstart:timestep(ruleset):tstop)
    fspan = lastframe:stopframe
    # Use the last frame of the existing simulation as the init frame
    init = output[lastindex(output)]
    data = initdata!(data, ruleset, init, tstart, nreplicates)
    setfps!(output, fps)
    runsim!(output, data, fspan)
end

"run the simulation either directly or asynchronously."
runsim!(output, args...) =
    if isasync(output)
        @async simloop!(output, args...)
    else
        simloop!(output, args...)
    end

"""
Loop over the selected timespan, running the ruleset and displaying the output

Operations on outputs and rulesets are allways mutable and in-place.
Operations on rules and data objects are functional as they are used in inner loops
where immutability improves performance.
"""
simloop!(output, data, fspan) = begin
    settimestamp!(output, first(fspan))
    # Initialise types etc
    data = updatetime(data, 1)# |> precalcrules
    # Loop over the simulation
    for f in fspan[2:end]
        # Get a data object with updated timestep and precalculated rules
        data = updatetime(data, f) |> precalcrules
        # Run the ruleset and setup data for the next iteration
        data = sequencerules!(data)
        # Save/do something with the the current frame
        storeframe!(output, data)
        isasync(output) && yield()
        # Stick to the FPS
        delay(output, f)
        # Exit gracefully
        if !isrunning(output) || f == last(fspan)
            showframe(output, data, f)
            setrunning!(output, false)
            setstoptime!(output, currenttime(data))
            finalize!(output)
            break
        end
    end
    output
end

"""
    precalcrules(rule, data) = rule

Rule precalculation. This is a functional approach rebuilding rules recursively.
@set from Setfield.jl helps in specific rule implementations.

The default is to return the existing rule
"""
precalcrules(rule, data) = rule

precalcrules(data::SimData) = @set data.ruleset.rules = precalcrules(rules(data), data)
precalcrules(data::MultiSimData) = @set data.data = map(precalcrules, data.data)
precalcrules(rules::Tuple, data) =
    (precalcrules(rules[1], data), precalcrules(tail(rules), data)...)
precalcrules(rules::Tuple{}, data) = ()
precalcrules(chain::Chain, data) = Chain(precalcrules(val(chain), data))


text:
using MachineElements
using Test

@testset "MachineElements.jl" begin
    # Write your tests here.
end


text:
module Trees

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export tree

register_normal_element("q__tree", context = @__MODULE__)

function tree(args...; kwargs...)
  q__tree(args...; kw(kwargs)...)
end

end


text:
using Base.Test

@testset "No constant tests" begin
    nothing
end


text:
module PyFortran90Namelists

include("init.jl")
include("convert.jl")
include("tokenizer.jl")
include("parser.jl")

end # module


text:
module HelperFunctions

export generatePosDefMatrix, transposeVectorized

  # generate a random pos def matrix with eigenvalues between 0.1 and 2
  function generatePosDefMatrix(n::Int64,rng)
      X = rand(rng,n,n)
      # any real square matrix can be QP decomposed into a orthogonal matrix and an uppertriangular matrix R
      Q, R = qr(X)
      eigs = rand(rng,n).*(2.-0.1) .+ 0.1
      X = Q*diagm(eigs)*Q'
      X = 0.5*(X+X')
      return X
  end

  # gives you the position of Aij in vec(A) where A in R^mxn
  function vecPos(m,n,i,j)
    ((i > m) || (j > n)) && error("Your indizes are outside the matrix dimension.")
    return Int((j-1)*m+i)
  end

  # permutes the cols of matrix A, s.t. A*vec(B) == A[:,p]*vec(B'), with B in R^{m_b x n_b}
  function transposeVectorized(A,mb,nb)
    # size of A
    m,n = size(A)

    # find permutation vector p
    k = 0
    p = zeros(Int64,n)
    for j=1:nb, i=1:mb
      k += 1
      p[vecPos(nb,mb,j,i)] = k
    end
    #println("Permutation vector: p=$(p)")
    # permute cols of A
    return A[:,p]

  end

end #MODULE


text:
function FMMtruncation(accuracy, boxSize, k)
    #TODO: use different truncation for different boxes???
    kd = sqrt(2)*boxSize*k
    P2 = ceil(Int, kd + 1.8*accuracy^(2/3)*(kd)^(1/3))
    Q = 2*P2 + 2
    return P2, Q
end

function divideSpace(centers::Array{Float64,2}, options; drawGroups = false)
    size(centers,1) < 2 && error("divideSpace: need at least 2 points")
    # using findfirst everywhere is slower but more robust to floating errors
    (x_min,y_min) = minimum(centers, dims=1)
    (x_max,y_max) = maximum(centers, dims=1)
    epss = eps(maximum(abs.([x_min;x_max;y_min;y_max])))
    x_max += epss
    y_max += epss
    x_min -= epss
    y_min -= epss
    a = options.dx
    nx = options.nx
    if a == 0
        (nx == 0) && error("divideSpace: no parameters given")
        #number of square boxes along x specified
        a = (x_max-x_min)/nx
    else
        #size of square boxes specified
        #recenter x to fit boxes
        nx = findfirst(n -> (a*n)>=(x_max-x_min),
                        1:Int(ceil((x_max-x_min)/a))+1)
        leftover_x = a*nx - (x_max-x_min)
        x_max += 0.5*leftover_x
        x_min -= 0.5*leftover_x
    end
    #recenter y to fit boxes
    ny = findfirst(n-> (a*n)>=(y_max-y_min),
                    1:Int(ceil((y_max-y_min)/a))+1)
    leftover_y = max(a*ny - (y_max-y_min),0.0)
    y_max += 0.5*leftover_y
    y_min -= 0.5*leftover_y
    FMMgroups = Array{FMMgroup}(undef,nx,ny)
    for ix = 1:nx, iy = 1:ny
        center_x = x_min + (ix-0.5)*a
        center_y = y_min + (iy-0.5)*a
        FMMgroups[ix,iy] = FMMgroup([], [center_x center_y])
    end
    for ic = 1:size(centers,1)
        ix = findfirst(n -> ((centers[ic,1] - x_min) <= a*n), 1:nx)
        iy = findfirst(n -> ((centers[ic,2] - y_min) <= a*n), 1:ny)
        push!(FMMgroups[ix,iy], ic)
    end
    #remove empty groups
    flags = Bool[length(FMMgroups[ii].point_ids) > 0 for ii=1:length(FMMgroups)]
    FMMgroups = FMMgroups[flags]

    if drawGroups
        figure()
        for ix = 0:nx
            plot([x_min + ix*a;x_min + ix*a],[y_min;y_max],"k",linewidth=2)
        end
        for iy = 0:ny
            plot([x_min;x_max],[y_min + iy*a;y_min + iy*a],"k",linewidth=2)
        end
        plot(centers[:,1],centers[:,2],"*")
        for ig = 1:length(FMMgroups)
            plot(FMMgroups[ig].center[1],FMMgroups[ig].center[2],"rx")
        end
        xlabel("\$x\$")
        ylabel("\$y\$")
        xlim([x_min - a;x_max + a])
        ylim([y_min - a;y_max + a])
        tight_layout()
        ax = gca()
        ax.set_aspect("equal", adjustable = "box")
    end
    return FMMgroups,a
end

# function FMMaggregationMatrix!(A, k, centers, box_center, t, P)
#     #P is length of first multipole expansion, *not* FMM
#     Q = length(t)
#     d = Array(Float64,2) #silly but faster
#     for i = 1:size(centers,1)
#         d[1] = centers[i,1] - box_center[1]
#         d[2] = centers[i,2] - box_center[2]
#         A[:,(i-1)*(2*P+1) .+ (1:2*P+1)] = [exp(-1.0im*(k*(cos(t[q])*d[1] + sin(t[q])*d[2])
#             + n*(pi/2-t[q]))) for q=1:Q, n=-P:P]
#     end
# end

function FMMaggregationMatrix(k, centers, box_center, t, P)
    #P is length of first multipole expansion, *not* FMM
    Q = length(t)
    Ns = size(centers,1)
    A = Array{Complex{Float64}}(undef, Q, Ns*(2*P+1))
    d = Array{Float64}(undef, 2) #silly but faster
    for i = 1:Ns
        d[1] = centers[i,1] - box_center[1]
        d[2] = centers[i,2] - box_center[2]
        A[:,(i-1)*(2*P+1) .+ (1:2*P+1)] = [exp(-1.0im*(k*(cos(t[q])*d[1] + sin(t[q])*d[2])
            + n*(pi/2-t[q]))) for q=1:Q, n=-P:P]
    end
    return A
end

# function FMMtranslationMatrix!(T, k, x, t, P2)
#     #x = c1-c2, where c2 is the source and c1 the destination
#     Q = length(t)
#     nx = sqrt(sum(abs2,x))
#     tx = Float64[atan2(x[2],x[1]) + pi/2 - t[q] for q=1:Q]
#     bess = besselh.(-P2:P2, 1, k*nx)
#     T[:] = 0
#     for i = -P2:P2
#         T[:] += bess[i+P2+1]*exp(1.0im*i*tx)/Q
#     end
# end

function FMMtranslationMatrix(k, x, t, P2)
    #x = c1-c2, where c2 is the source and c1 the destination
    Q = length(t)
    T = zeros(Complex{Float64},Q)
    nx = sqrt(sum(abs2,x))
    tx = Float64[-atan(-x[2],-x[1]) - pi/2 + t[q] for q=1:Q]
    bess = besselh.(-P2:P2, 1, k*nx)
    T[:] .= 0
    for i = -P2:P2, j = 1:Q
        T[j] += bess[i+P2+1]*exp(1.0im*i*tx[j])/Q
    end
    return T
end

function FMMnearMatrix(k, P, groups, centers, boxSize, num)
    #for now, does not include self interactions - but does include minus sign
    Ns = size(centers,1)
    G = length(groups)
    W = 2*P+1
    Is = Array{Int}(undef, num*W^2)
    Js = Array{Int}(undef, num*W^2)
    Zs = Array{Complex{Float64}}(undef, num*W^2)
    bess = Array{Complex{Float64}}(undef, 2*P+1)
    mindist2 = 3*boxSize^2 #anywhere between 2 and 4
    ind = 0
    for ig1 = 1:G
        for ig2 = 1:G
            x = groups[ig1].center - groups[ig2].center
            sum(abs2,x) > mindist2 && continue
            #for each enclosed scatterer, build translation matrix
            for ic1 = 1:length(groups[ig1].point_ids)
                for ic2 = 1:length(groups[ig2].point_ids)
                    ig1 == ig2 && (ic1 == ic2 && continue) #no  self-interactions
                    d = centers[groups[ig1].point_ids[ic1],:] - centers[groups[ig2].point_ids[ic2],:]
                    nd = sqrt(sum(abs2,d))
                	td = atan(d[2],d[1])
                    bess[:] = besselh.(0:2*P,1,k*nd)
                	for ix = 1:2*P #lower diagonals
                		rng = ix+1:1+W:W^2-W*ix
                		Zs[ind .+ rng] .= -exp(-1im*td*ix)*(-1)^(ix)*bess[ix+1]
                	end
                	for ix = 0:2*P #central and upper diagonals
                		rng = ix*W+1:1+W:W^2-ix
                		Zs[ind .+ rng] .= -exp(1im*td*ix)*bess[ix+1]
                	end
                    for ij = 1:W
                        Is[ind .+ (1:W)] = collect((groups[ig1].point_ids[ic1]-1)*W .+ (1:W))
                        Js[ind .+ (1:W)] .= (groups[ig2].point_ids[ic2]-1)*W + ij
                        ind += W
                    end
                end
            end
        end
    end
    return sparse(Is,Js,Zs,Ns*W,Ns*W)
end

function FMMnearMatrix_upperTri(k, P, groups, centers, boxSize, num)
    #half the bessels here are re-used, since hankel part is symmetric between a->b and b->a
    #for now, does not include self interactions - but does minus sign
    Ns = size(centers,1)
    G = length(groups)
    W = 2*P+1
    Is = Array{Int}(undef, num*W^2)
    Js = Array{Int}(undef, num*W^2)
    Zs = Array{Complex{Float64}}(undef, num*W^2)
    mindist2 = 3*boxSize^2 #anywhere between 2 and 4
    ind = 0
    d = Array{Float64}(undef, 2)
    bess = Array{Complex{Float64}}(undef, 2*P+1)
    #first, between group and itself
    for ig1 = 1:G
        for ic1 = 1:length(groups[ig1].point_ids)
            for ic2 = ic1+1:length(groups[ig1].point_ids)
                d[1] = centers[groups[ig1].point_ids[ic1],1] - centers[groups[ig1].point_ids[ic2],1]
                d[2] = centers[groups[ig1].point_ids[ic1],2] - centers[groups[ig1].point_ids[ic2],2]
                nd = sqrt(sum(abs2,d))
                #first for ic2->ic1, then ic1->ic2
                td1 = atan(d[2],d[1])
                td2 = atan(-d[2],-d[1]) #adding pi instead leads to error
                bess[:] = besselh.(0:2*P,1,k*nd)
                ind2 = ind + W^2

                for ix = 1:2*P #lower diagonals
                    rng = ix+1:1+W:W^2-W*ix
                    Zs[ind .+ rng] .= -exp(-1.0im*td1*ix)*(-1)^(ix)*bess[ix+1]
                    Zs[ind2 .+ rng] .= -exp(-1.0im*td2*ix)*(-1)^(ix)*bess[ix+1]
                end
                for ix = 0:2*P #central and upper diagonals
                    rng = ix*W+1:1+W:W^2-ix
                    Zs[ind .+ rng] .= -exp(1.0im*td1*ix)*bess[ix+1]
                    Zs[ind2 .+ rng] .= -exp(1.0im*td2*ix)*bess[ix+1]
                end
                #this way we avoid allocating vectors
                for ij = 1:W
                    rng = (1:W) .+ (ij-1)*W
                    Js[ind .+ rng] .= (groups[ig1].point_ids[ic2]-1)*W + ij
                    Js[ind2 .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij
                end
                for ij = 1:W
                    rng = (1:W:W^2-W+1) .+ (ij-1)
                    Is[ind .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij
                    Is[ind2 .+ rng] .= (groups[ig1].point_ids[ic2]-1)*W + ij
                end
                ind += 2*W^2 #skip over ind2 data
            end
        end
    end
    #now, between groups
    for ig1 = 1:G
        for ig2 = ig1+1:G
            x = groups[ig1].center - groups[ig2].center
            sum(abs2,x) > mindist2 && continue
            #for each enclosed scatterer, build translation matrix
            for ic1 = 1:length(groups[ig1].point_ids)
                for ic2 = 1:length(groups[ig2].point_ids)
                    d[1] = centers[groups[ig1].point_ids[ic1],1] - centers[groups[ig2].point_ids[ic2],1]
                    d[2] = centers[groups[ig1].point_ids[ic1],2] - centers[groups[ig2].point_ids[ic2],2]
                    nd = sqrt(sum(abs2,d))
                	#first for ig2->ig1, then ig1->ig2
                    td1 = atan(d[2],d[1])
                    td2 = atan(-d[2],-d[1]) #adding pi instead leads to error
                    bess[:] = besselh.(0:2*P,1,k*nd)
                    ind2 = ind + W^2

                    for ix = 1:2*P #lower diagonals
                		rng = ix+1:1+W:W^2-W*ix
                		Zs[ind .+ rng] .= -exp(-1.0im*td1*ix)*(-1)^(ix)*bess[ix+1]
                        Zs[ind2 .+ rng] .= -exp(-1.0im*td2*ix)*(-1)^(ix)*bess[ix+1]
                	end
                	for ix = 0:2*P #central and upper diagonals
                		rng = ix*W+1:1+W:W^2-ix
                		Zs[ind .+ rng] .= -exp(1.0im*td1*ix)*bess[ix+1]
                        Zs[ind2 .+ rng] .= -exp(1.0im*td2*ix)*bess[ix+1]
                	end
                    #this way we avoid allocating vectors
                    for ij = 1:W
                        rng = (1:W) .+ (ij-1)*W
                        Js[ind .+ rng] .= (groups[ig2].point_ids[ic2]-1)*W + ij
                        Js[ind2 .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij
                    end
                    for ij = 1:W
                        rng = (1:W:W^2-W+1) .+ (ij-1)
                        Is[ind .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij
                        Is[ind2 .+ rng] .= (groups[ig2].point_ids[ic2]-1)*W + ij
                    end
                    ind += 2*W^2 #skip over ind2 data
                end
            end
        end
    end
    return sparse(Is,Js,Zs,Ns*W,Ns*W)
end

#TODO: distinguish between symmetric and non-symmetric.
function FMMbuildMatrices(k, P, P2, Q, groups, centers, boxSize; tri = true)
    G = length(groups)
    t = Float64[2*pi*j/Q for j=0:(Q-1)]
    Agg = Array{Complex{Float64},2}[]
    for ig = 1:G #preallocate here? not sure what is faster in practice
        A = FMMaggregationMatrix(k, centers[groups[ig].point_ids,:], groups[ig].center, t, P)
        push!(Agg, A)
    end
    Disagg = [Agg[ii]' for ii=1:G]
    #Array of arrays for Trans, or dict, or something else???
    Trans = Vector{Complex{Float64}}[]
    #For now, the translation from j to i is at Trans[(i-1)*G + j]
    mindist2 = 3*boxSize^2 #anywhere between 2 and 4
    d = Array{Float64}(undef, 2)
    num = 0
    for ig1 = 1:G
        for ig2 = 1:G
            if ig1 == ig2
                push!(Trans,[])
                #number of near interactions inside a group is their product minus self-interactions
                num += length(groups[ig1].point_ids)*(length(groups[ig1].point_ids)-1)
                continue
            end
            d[1] = groups[ig1].center[1] - groups[ig2].center[1]
            d[2] = groups[ig1].center[2] - groups[ig2].center[2]
            if sum(abs2,d) <= mindist2
                push!(Trans,[])
                #number of near interactions between two close (but not identical) groups is their product
                num += length(groups[ig1].point_ids)*length(groups[ig2].point_ids)
            else
                #preallocate here???
                Trans1_2 = FMMtranslationMatrix(k, d, t, P2)
                push!(Trans,Trans1_2)
            end
        end
    end
    if tri
        Znear = FMMnearMatrix_upperTri(k, P, groups, centers, boxSize, num)
    else
        Znear = FMMnearMatrix(k, P, groups, centers, boxSize, num)
    end

    FMMmatrices(Agg, Trans, Disagg, Znear, groups, P2, Q)
end