text:
using Statistics
using CUDA.CUDNN:
    cudnnDropoutForward,
    cudnnDropoutForward!,
    cudnnDropoutBackward,
    cudnnDropoutSeed,
    cudnnDropoutDescriptor,
        cudnnDropoutDescriptor_t,
        cudnnCreateDropoutDescriptor,
        cudnnSetDropoutDescriptor,
        cudnnGetDropoutDescriptor,
        cudnnRestoreDropoutDescriptor,
        cudnnDestroyDropoutDescriptor,
    cudnnDropoutGetStatesSize,
    cudnnDropoutGetReserveSpaceSize

@testset "cudnn/dropout" begin
    @test cudnnDropoutDescriptor(C_NULL) isa cudnnDropoutDescriptor
    @test Base.unsafe_convert(Ptr, cudnnDropoutDescriptor(C_NULL)) isa Ptr
    @test cudnnDropoutDescriptor(0.5) isa cudnnDropoutDescriptor

    N,P = 1000, 0.7
    x = CUDA.rand(N)
    d = cudnnDropoutDescriptor(P)
    cudnnDropoutSeed[] = 1
    y = cudnnDropoutForward(x; dropout = P)
    @test isapprox(mean(Array(y).==0), P; atol = 3/sqrt(N))
    @test y == cudnnDropoutForward(x, d)
    @test y == cudnnDropoutForward!(similar(x), x; dropout = P)
    @test y == cudnnDropoutForward!(similar(x), x, d)
    cudnnDropoutSeed[] = -1
end


text:
using Test, SparseArrays, Missings

# Must be defined outside testset on v1.0
struct CubeRooter end
(::CubeRooter)(x) = cbrt(x)

@testset "Missings" begin
    x = Missings.replace([1, 2, missing, 4], 3)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test axes(x) == (1:4,)
    @test keys(x) == eachindex(x) == 1:4
    @test collect(x) == [x[i] for i in keys(x)] == 1:4
    @test collect(x) isa Vector{Int}
    x = Missings.replace([1, 2, missing, 4], 3.0)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test axes(x) == (1:4,)
    @test keys(x) == eachindex(x) == 1:4
    @test collect(x) == [x[i] for i in keys(x)] == 1:4
    @test collect(x) isa Vector{Int}
    x = Missings.replace([1 2; missing 4], 3)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (2, 2)
    @test axes(x) == (1:2, 1:2)
    @test keys(x) == CartesianIndices((1:2, 1:2))
    @test eachindex(x) == 1:4
    @test collect(x) == [x[i] for i in keys(x)] == [1 2; 3 4]
    @test collect(x) isa Matrix{Int}
    x = Missings.replace((v for v in [missing, 1, missing, 2, 4]), 0)
    @test length(x) == 5
    @test size(x) == (5,)
    if VERSION >= v"1.7.0-DEV"
        @test eachindex(x) == keys(x) == 1:5
    else
        @test_throws MethodError keys(x)
        @test_throws MethodError eachindex(x)
    end
    @test_throws MethodError x[1]
    @test eltype(x) === Any
    @test collect(x) == [0, 1, 0, 2, 4]
    @test collect(x) isa Vector{Int}

    x = Missings.fail([1, 2, 3, 4])
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test axes(x) == (1:4,)
    @test keys(x) == eachindex(x) == 1:4
    @test collect(x) == [x[i] for i in keys(x)] == [1, 2, 3, 4]
    @test collect(x) isa Vector{Int}
    x = Missings.fail([1 2; 3 4])
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (2, 2)
    @test axes(x) == (1:2, 1:2)
    @test keys(x) == CartesianIndices((1:2, 1:2))
    @test eachindex(x) == 1:4
    @test collect(x) == [x[i] for i in keys(x)] == [1 2; 3 4]
    @test collect(x) isa Matrix{Int}
    @test_throws MissingException collect(Missings.fail([1, 2, missing, 4]))
    x = Missings.fail(v for v in [1, 2, 4])
    @test eltype(x) === Any
    @test length(x) == 3
    @test size(x) == (3,)
    @test axes(x) == (1:3,)
    if VERSION >= v"1.7.0-DEV"
        @test eachindex(x) == keys(x) == 1:3
    else
        @test_throws MethodError keys(x)
        @test_throws MethodError eachindex(x)
    end
    @test_throws MethodError x[1]
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}

    x = skipmissing([1, 2, missing, 4])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = skipmissing([1  2; missing 4])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = collect(skipmissing([missing]))
    @test eltype(x) === Union{}
    @test isempty(collect(x))
    @test collect(x) isa Vector{Union{}}
    x = collect(skipmissing(Union{Int, Missing}[]))
    @test eltype(x) === Int
    @test isempty(collect(x))
    @test collect(x) isa Vector{Int}
    x = skipmissing([missing, missing, 1, 2, missing, 4, missing, missing])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = skipmissing(v for v in [missing, 1, missing, 2, 4])
    @test eltype(x) === Any
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}

    @test levels(1:1) == levels([1]) == levels([1, missing]) == levels([missing, 1]) == [1]
    @test levels(2:-1:1) == levels([2, 1]) == levels([2, missing, 1]) == [1, 2]
    @test levels([missing, "a", "c", missing, "b"]) == ["a", "b", "c"]
    @test levels([Complex(0, 1), Complex(1, 0), missing]) == [Complex(0, 1), Complex(1, 0)]
    @test levels(sparse([0 3 2])) == [0, 2, 3]
    @test typeof(levels([1])) === typeof(levels([1, missing])) === Vector{Int}
    @test typeof(levels(["a"])) === typeof(levels(["a", missing])) === Vector{String}
    @test typeof(levels(sparse([1]))) === Vector{Int}
    @test isempty(levels([missing]))
    @test isempty(levels([]))

    @test nonmissingtype(Union{Int, Missing}) == Int
    @test nonmissingtype(Any) == Any
    @test nonmissingtype(Missing) == Union{}
    @test nonmissingtype(Union{Array{Int}, Missing}) == Array{Int}

    @test isequal(missings(1), [missing])
    @test isequal(missings(Any, 1), [missing])
    @test isequal(missings(Int, 1), [missing])
    @test missings(Int, 1) isa Vector{Union{Int, Missing}}
    @test missings(Any, 1) isa Vector{Union{Any, Missing}}
    @test isequal(missings(Union{Int, Missing}, 1, 2), [missing missing])
    @test missings(Union{Int, Missing}, 1, 2) isa Matrix{Union{Int, Missing}}
    @test Union{Int, Missing}[1,2,3] == (Union{Int, Missing})[1,2,3]
    x = missings(Int, (1, 2))
    @test isa(x, Matrix{Union{Int, Missing}})
    @test isequal(x, [missing missing])
    x = missings((1, 2))
    @test isa(x, Matrix{Missing})
    @test isequal(x, [missing missing])

    @test allowmissing([1]) == [1]
    @test allowmissing([1]) isa AbstractVector{Union{Int, Missing}}
    @test allowmissing(Any[:a]) == [:a]
    @test allowmissing(Any[:a]) isa AbstractVector{Any}
    @test isequal(allowmissing([1, missing]), [1, missing])
    @test allowmissing([1, missing]) isa AbstractVector{Union{Int, Missing}}
    @test isequal(allowmissing([missing]), [missing])
    @test allowmissing([missing]) isa AbstractVector{Missing}

    @test allowmissing([1 1]) == [1 1]
    @test allowmissing([1 1]) isa AbstractArray{Union{Int, Missing}, 2}
    @test allowmissing([:a 1]) == [:a 1]
    @test allowmissing([:a 1]) isa AbstractArray{Any, 2}
    @test isequal(allowmissing([1 missing]), [1 missing])
    @test allowmissing([1 missing]) isa AbstractArray{Union{Int, Missing}, 2}
    @test isequal(allowmissing([missing missing]), [missing missing])
    @test allowmissing([missing missing]) isa AbstractArray{Missing, 2}

    @test disallowmissing(Union{Int, Missing}[1]) == [1]
    @test disallowmissing(Union{Int, Missing}[1]) isa AbstractVector{Int}
    @test disallowmissing([1]) == [1]
    @test disallowmissing([1]) isa AbstractVector{Int}
    @test disallowmissing(Any[:a]) == [:a]
    @test disallowmissing(Any[:a]) isa AbstractVector{Any}
    @test_throws MethodError disallowmissing([1, missing])
    @test_throws MethodError disallowmissing([missing])

    @test disallowmissing(Union{Int, Missing}[1 1]) == [1 1]
    @test disallowmissing(Union{Int, Missing}[1 1]) isa AbstractArray{Int, 2}
    @test disallowmissing([1 1]) == [1 1]
    @test disallowmissing([1 1]) isa AbstractArray{Int, 2}
    @test disallowmissing([:a 1]) == [:a 1]
    @test disallowmissing([:a 1]) isa AbstractArray{Any, 2}
    @test_throws MethodError disallowmissing([1 missing])
    @test_throws MethodError disallowmissing([missing missing])

    # Lifting
    ## functor
    cuberoot = CubeRooter()  # defined at top of file
    @test passmissing(cuberoot)(27) == 3.0
    @test isequal(passmissing(cuberoot)(missing), missing)
    ## type
    @test passmissing(Int)(1.0) == 1
    @test isequal(passmissing(Int)(missing), missing)
    ## function
    @test passmissing(sqrt)(4) == 2.0
    @test isequal(passmissing(sqrt)(missing), missing)
    @test isequal(passmissing(sqrt).([missing, 4]), [missing, 2.0])
    @test passmissing((x,y)->"$x $y")(1, 2) == "1 2"
    @test isequal(passmissing((x,y)->"$x $y")(missing), missing)
    if VERSION >= v"1.4.0-DEV"
        @test_throws MethodError passmissing(string)(missing, base=2)
    else
        @test_throws ErrorException passmissing(string)(missing, base=2)
    end

    @test passmissing(sin) === Missings.PassMissing{typeof(sin)}(sin)
    @test passmissing(Int) === Missings.PassMissing{Type{Int}}(Int)
    @test passmissing(cuberoot) === Missings.PassMissing{CubeRooter}(cuberoot)

    @testset "deprecated" begin
        x = [1, 2, missing, 4]
        y = ["a", "b", "c", missing]
        z = [missing, missing, 3.1, 4.5]
        l = [1, 2, 3, 4, 5]
        @test_throws ArgumentError skipmissings(x, l)
        mx, my = skipmissings(x, y)
        iobuf = IOBuffer()
        show(iobuf, MIME("text/plain"), mx)
        s = String(take!(iobuf))
        @test s == "Missings.SkipMissings{$(Vector{Union{Missing, Int}})}(" *
            "Union{Missing, $Int}[1, 2, missing, 4]) comprised of 2 iterators"
        @test collect(mx) == [1, 2]
        @test collect(mx) isa Vector{Int}
        @test reduce(+, mx) === reduce(+, collect(mx)) === sum(mx) ===
            mapreduce(identity, +, mx) === 3
        @test mapreduce(x -> x^2, +, mx) === mapreduce(x -> x^2, +, collect(mx)) === 5
        mx, my, mz = skipmissings(x, y, z)
        @test eltype(mx) == Int
        @test eltype(my) == String
        @test eltype(mz) == Float64
        @test isempty(collect(mx))
        @test sum(mx) === 0
        x = [missing 4; 2 5; 3 6]
        y = [1 4; missing 5; 3 6]
        mx, my = skipmissings(x, y)
        @test collect(mx) == [3, 4, 5, 6]
        @test mx[3] == 3
        @test_throws MissingException mx[1]
        @test reduce(+, mx) === 18
        @test isapprox(mapreduce(cos, *, collect(mx)),  mapreduce(cos, *, mx))
        @static if VERSION >= v"1.4.0-DEV"
            @inferred Union{Float64, Missing} mapreduce(cos, *, mx)
            @inferred Union{Float64, Missing} sum(mx)
            @inferred Union{Float64, Missing} reduce(+, mx)
        end

        x = [missing missing missing]
        y = [1, 2, 3]
        mx, my = skipmissings(x, y)
        @test_throws ArgumentError reduce(x -> x/2, mx)
        @test_throws ArgumentError mapreduce(x -> x/2, +, mx)
        @test_throws MethodError length(mx)
        @test IndexStyle(typeof(mx)) == IndexStyle(typeof(x))
        x = [isodd(i) ? missing : i for i in 1:64]
        y = [isodd(i) ? missing : i for i in 65:128]
        mx, my = skipmissings(x, y)
        @test sum(mx) === 1056
        @static if VERSION >= v"1.4.0-DEV"
            @inferred Union{Missing, Int} sum(mx)
            @inferred Union{Missing, Int} reduce(+, mx)
        end
    end
end

ext: jl


text:
using Combinatorics

using dataset_pb2

using BondLengthDistributions
using SmuMolecules
using SmuUtilities

const THRESHOLD = 2.0

"""Generate a BondTopology that joins each Hydrogen atom to its nearest
    heavy atom.
Args:
  bond_topology:
  distances:
Returns:
"""
function hydrogen_to_nearest_atom(bond_topology::dataset_pb2.BondTopology,
                             distances::Array{Float32,2})::Union{dataset_pb2.BondTopology, Nothing}
  result = dataset_pb2.BondTopology()
  # Not sure if this is safe or not.
  setproperty!(result, :atoms, bond_topology.atoms)

  single_bond = 1

  natoms = length(bond_topology.atoms)
  for a1 in 1:natoms
    bond_topology.atoms[a1] == dataset_pb2.BondTopology_AtomType.ATOM_H || continue

    shortest_distance = 1.0e+30
    closest_heavy_atom = -1
    for a2 in 1:natoms
      bond_topology.atoms[a2] == dataset_pb2.BondTopology_AtomType.ATOM_H && continue

      distances[a1, a2] >= THRESHOLD && continue

      if distances[a1, a2] < shortest_distance
        shortest_distance = distances[a1, a2]
        closest_heavy_atom = a2
      end
    end

    closest_heavy_atom < 0 && return nothing

    add_bond!(a1 - 1, closest_heavy_atom - 1, single_bond, result)
  end

  return result
end

# Enable sorting of BondTopology's
import Base: isless
Base.isless(bt1::BondTopology, bt2::BondTopology) = isless(bt1.score, bt2.score)

# Convenenience function for filtering scores.
nonzero = x->x>0.0

# Convenenience function for selecting non hydrogen atoms
not_hydrogen = x->x != BondTopology_AtomType.ATOM_H

"""Return all BondTopology's that are plausible.

  Given a molecule described by `bond_topology` and `geometry`, return all possible
  BondTopology that are consistent with that.
  Note that `bond_topology` will be put in a canonical form.
  Args:
    bond_length_distribution:
    bond_topology:
    geometry:
  Returns:
    TopologyMatches
"""
function bond_topologies_from_geom(
    bond_lengths::AllBondLengthDistributions,
    bond_topology::dataset_pb2.BondTopology,
    geometry::dataset_pb2.Geometry,
    matching_parameters::MatchingParameters)::dataset_pb2.TopologyMatches
  result = dataset_pb2.TopologyMatches()    # To be returned.
  length(bond_topology.atoms) == 1 && return result  # return an empty result

  canonical_bond_topology!(bond_topology)
  distances = SmuUtilities.distances(geometry)

  # First join each Hydrogen to its nearest heavy atom, thereby
  # creating a starting BondTopology from which all others can grow
  starting_bond_topology = hydrogen_to_nearest_atom(bond_topology, distances)
  starting_bond_topology === nothing && return result

  heavy_atom_indices = findall(not_hydrogen, bond_topology.atoms)
  length(heavy_atom_indices) < 2 && return result

  # For each atom pair, a list of possible bond types.
  # Key is a tuple of the two atom numbers, value is a Vector
  # with the score for each bond type.

  bonds_to_scores = Dict{Tuple{Int32,Int32}, Vector{Float32}}()

  for (i, j) in Combinatorics.combinations(heavy_atom_indices, 2)  # All pairs
    (dist = distances[i, j]) > THRESHOLD && continue

    scores = [pdf(bond_lengths, bond_topology.atoms[i], bond_topology.atoms[j], btype, dist)
              for btype in 0:4]

    any(nonzero, scores) && (bonds_to_scores[(i, j)] = scores)
  end

  @debug("bonds_to_scores $(bonds_to_scores)")
  isempty(bonds_to_scores) && return result

  found_topologies = Vector{BondTopology}()  # Will be set into `result`.

  mol = SmuMolecule(starting_bond_topology, bonds_to_scores, matching_parameters)

  search_space = generate_search_state(mol)
  for s in Iterators.product(search_space...)
    @debug("Placing state $s")
    (bt = place_bonds!(s, mol)) == nothing && continue
    is_single_fragment(bt) || continue

    canonical_bond_topology!(bt)
    same_bond_topology(bond_topology, bt) && (bt.is_starting_topology = true)
    # smiles not set
    push!(found_topologies, bt)
  end

  isempty(found_topologies) && return result

  length(found_topologies) > 1 && sort!(found_topologies, rev=true)

  setproperty!(result, :bond_topology, found_topologies)

  return result
end

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
# types needed by flux functions
# these need to be declared in a separate file from the flux functions
# themselves because they need to be defind *before* ParamType

"""
  Data for [`RoeSolver`](@ref) and related methods
"""
struct RoeFluxData{Tsol, Tres, Tmsh}
  # primal flux
  dq::Vector{Tsol}
  sat::Vector{Tres}
  roe_vars::Vector{Tres}
  euler_flux::Vector{Tres}
  v_vals::Vector{Tsol}
  nrm2::Vector{Tmsh}

  # revm
  euler_flux_bar::Vector{Tres}
  sat_bar::Vector{Tres}

  # revq
  v_vals_bar::Vector{Tres}
  dq_bar::Vector{Tres}
  roe_vars_bar::Vector{Tres}

  # diff
  roe_vars_dot::Vector{Tres}

  function RoeFluxData{Tsol, Tres, Tmsh}(numDofPerNode::Integer, dim::Integer) where {Tsol, Tres, Tmsh}

    # primal flux
    dq = zeros(Tsol, numDofPerNode)
    sat = zeros(Tres, numDofPerNode)
    roe_vars = zeros(Tres, numDofPerNode)
    euler_flux = zeros(Tres, numDofPerNode)
    v_vals = zeros(Tsol, numDofPerNode)
    nrm2 = zeros(Tmsh, dim)

    # revm
    euler_flux_bar = zeros(Tres, numDofPerNode)
    sat_bar = zeros(Tres, numDofPerNode)

    # revq
    v_vals_bar = zeros(Tres, numDofPerNode)
    dq_bar = zeros(Tres, numDofPerNode)
    roe_vars_bar = zeros(Tres, 4)

    # diff
    roe_vars_dot = zeros(Tres, 22)  # 22 needed in 3D

    obj = new(dq, sat, roe_vars, euler_flux, v_vals, nrm2,
              euler_flux_bar, sat_bar,
              v_vals_bar, dq_bar, roe_vars_bar,
              roe_vars_dot)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end

end

"""
  Data for [`calcSat`](@ref) and related methods
"""
struct CalcSatData{Tres}
  E1dq::Vector{Tres}
  E2dq::Vector{Tres}
  E3dq::Vector{Tres}
  E4dq::Vector{Tres}

  # no arrays are used for the diff version

  # revq
  E1dq_bar::Vector{Tres}
  E2dq_bar::Vector{Tres}
  E3dq_bar::Vector{Tres}
  E4dq_bar::Vector{Tres}

  function CalcSatData{Tres}(numDofPerNode::Integer) where {Tres}

    E1dq = zeros(Tres, numDofPerNode)
    E2dq = zeros(Tres, numDofPerNode)
    E3dq = zeros(Tres, numDofPerNode)
    E4dq = zeros(Tres, numDofPerNode)

    # revq
    E1dq_bar = zeros(Tres, numDofPerNode)
    E2dq_bar = zeros(Tres, numDofPerNode)
    E3dq_bar = zeros(Tres, numDofPerNode)
    E4dq_bar = zeros(Tres, numDofPerNode)

    obj = new(E1dq, E2dq, E3dq, E4dq,
              E1dq_bar, E2dq_bar, E3dq_bar, E4dq_bar)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Data needed by [`calcLFFlux`](@ref) and related method

  **Static Parameters**

   * Tres
"""
struct LFFluxData{Tres}
  fluxL::Vector{Tres}
  fluxR::Vector{Tres}
  F_dotL::Array{Tres, 2}
  F_dotR::Array{Tres, 2}
  lambda_dotL::Vector{Tres}
  lambda_dotR::Vector{Tres}

  function LFFluxData{Tres}(numDofPerNode::Integer, nd::Integer) where {Tres}

    @assert numDofPerNode == nd  # not sure what happens if this is not true
    fluxL = zeros(Tres, numDofPerNode)
    fluxR = zeros(Tres, numDofPerNode)
    F_dotL = zeros(Tres, numDofPerNode, nd)
    F_dotR = zeros(Tres, numDofPerNode, nd)
    lambda_dotL = zeros(Tres, numDofPerNode)
    lambda_dotR = zeros(Tres, numDofPerNode)

    obj = new(fluxL, fluxR, F_dotL, F_dotR, lambda_dotL, lambda_dotR)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Data needed by [`logavg_diff`](@ref)

  **Static Parameters**

   * Tl: datatype of left state
   * Tr: datatype of right state
"""
struct LogAvgData{Tl, Tr}
  xi_dotL::Vector{Tl}
  xi_dotR::Vector{Tr}
  f_dotL::Vector{Tl}
  f_dotR::Vector{Tr}
  u_dotL::Vector{Tl}
  u_dotR::Vector{Tr}
  F_dotL::Vector{Tl}
  F_dotR::Vector{Tr}

  function LogAvgData{Tl, Tr}(nd::Integer) where {Tl, Tr}
    xi_dotL = zeros(Tl, nd)
    xi_dotR = zeros(Tr, nd)
    f_dotL = zeros(Tl, nd)
    f_dotR = zeros(Tr, nd)
    u_dotL = zeros(Tl, nd)
    u_dotR = zeros(Tr, nd)
    F_dotL = zeros(Tl, nd)
    F_dotR = zeros(Tr, nd)

    return new(xi_dotL, xi_dotR, f_dotL, f_dotR, u_dotL, u_dotR, F_dotL, F_dotR)
  end
end

"""
  This struct holds all the temporary arrays needed to calculate the IR flux
"""
struct IRFluxData{Tsol}
  pL_dot::Vector{Tsol}
  pR_dot::Vector{Tsol}
  z1L_dot::Vector{Tsol}
  z2L_dot::Vector{Tsol}
  z3L_dot::Vector{Tsol}
  z4L_dot::Vector{Tsol}
  z5L_dot::Vector{Tsol}

  z1R_dot::Vector{Tsol}
  z2R_dot::Vector{Tsol}
  z3R_dot::Vector{Tsol}
  z4R_dot::Vector{Tsol}
  z5R_dot::Vector{Tsol}

  z5avg_dotL::Vector{Tsol}
  z5avg_dotR::Vector{Tsol}
  z4avg_dotL::Vector{Tsol}
  z4avg_dotR::Vector{Tsol}
  z1avg_dotL::Vector{Tsol}
  z1avg_dotR::Vector{Tsol}

  rho_hat_dotL::Vector{Tsol}
  rho_hat_dotR::Vector{Tsol}
  u_hat_dotL::Vector{Tsol}
  u_hat_dotR::Vector{Tsol}
  v_hat_dotL::Vector{Tsol}
  v_hat_dotR::Vector{Tsol}
  w_hat_dotL::Vector{Tsol}
  w_hat_dotR::Vector{Tsol}
  p1_hat_dotL::Vector{Tsol}
  p1_hat_dotR::Vector{Tsol}
  h_hat_dotL::Vector{Tsol}
  h_hat_dotR::Vector{Tsol}
  avgdata::LogAvgData{Tsol, Tsol}

  """
    Constructor for IRFluxData.

    **Inputs**

     * nd: number of derivative components
  """
  function IRFluxData{Tsol}(nd::Integer) where {Tsol}

    #TODO: consider making these views of an array to get spatial locality
    pL_dot = zeros(Tsol, nd)
    pR_dot = zeros(Tsol, nd)

    z1L_dot = zeros(Tsol, nd)
    z2L_dot = zeros(Tsol, nd)
    z3L_dot = zeros(Tsol, nd)
    z4L_dot = zeros(Tsol, nd)
    z5L_dot = zeros(Tsol, nd)

    z1R_dot = zeros(Tsol, nd)
    z2R_dot = zeros(Tsol, nd)
    z3R_dot = zeros(Tsol, nd)
    z4R_dot = zeros(Tsol, nd)
    z5R_dot = zeros(Tsol, nd)

    z5avg_dotL = zeros(Tsol, nd)
    z5avg_dotR = zeros(Tsol, nd)
    z4avg_dotL = zeros(Tsol, nd)
    z4avg_dotR = zeros(Tsol, nd)
    z1avg_dotL = zeros(Tsol, nd)
    z1avg_dotR = zeros(Tsol, nd)

    rho_hat_dotL = zeros(Tsol, nd)
    rho_hat_dotR = zeros(Tsol, nd)
    u_hat_dotL = zeros(Tsol, nd)
    u_hat_dotR = zeros(Tsol, nd)
    v_hat_dotL = zeros(Tsol, nd)
    v_hat_dotR = zeros(Tsol, nd)
    w_hat_dotL = zeros(Tsol, nd)
    w_hat_dotR = zeros(Tsol, nd)
    p1_hat_dotL = zeros(Tsol, nd)
    p1_hat_dotR = zeros(Tsol, nd)
    h_hat_dotL = zeros(Tsol, nd)
    h_hat_dotR = zeros(Tsol, nd)

    logdata = LogAvgData{Tsol, Tsol}(nd)

    obj =  new(pL_dot, pR_dot, z1L_dot, z2L_dot, z3L_dot, z4L_dot, z5L_dot,
               z1R_dot, z2R_dot, z3R_dot, z4R_dot, z5R_dot,
               z5avg_dotL, z5avg_dotR, z4avg_dotL, z4avg_dotR, z1avg_dotL,
               z1avg_dotR,
               rho_hat_dotL, rho_hat_dotR, u_hat_dotL, u_hat_dotR,
               v_hat_dotL, v_hat_dotR, w_hat_dotL, w_hat_dotR,
               p1_hat_dotL, p1_hat_dotR, h_hat_dotL,h_hat_dotR,
               logdata)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Temporary storage for [`calcHLLFlux`](@ref)
"""
struct HLLFluxData{Tres}
  fL::Vector{Tres}
  fR::Vector{Tres}

  sL_dot::Matrix{Tres}
  sR_dot::Matrix{Tres}
  fL_jac::Matrix{Tres}
  fR_jac::Matrix{Tres}

  fL_bar::Vector{Tres}
  fR_bar::Vector{Tres}

  function HLLFluxData{Tres}(numDofPerNode::Integer) where {Tres}

    fL = zeros(Tres, numDofPerNode)
    fR = zeros(Tres, numDofPerNode)

    sL_dot = zeros(Tres, numDofPerNode, 2)
    sR_dot = zeros(Tres, numDofPerNode, 2)
    fL_jac = zeros(Tres, numDofPerNode,numDofPerNode)
    fR_jac = zeros(Tres, numDofPerNode, numDofPerNode)

    fL_bar = zeros(Tres, numDofPerNode)
    fR_bar = zeros(Tres, numDofPerNode)

    return new(fL, fR,
              sL_dot, sR_dot, fL_jac, fR_jac,
              fL_bar, fR_bar)
  end
end

"""
  Temporary arrays for [`applyEntropyKernel_diagE`](@ref)
"""
struct ApplyEntropyKernel_diagEData{Tsol, Tres}
  q_avg::Vector{Tsol}

  # diff method
  q_avg_dot::Matrix{Tres}
  F::Vector{Tres}
  F_dot::Matrix{Tres}

  # for applyEntropyKernel_inner
  vL::Vector{Tsol}
  vR::Vector{Tsol}
  F_tmp::Vector{Tres}

  # diff method
  delta_w_dot::Matrix{Tsol}

  # rev method
  q_avg_bar::Vector{Tres}
  delta_w::Vector{Tsol}
  delta_w_bar::Vector{Tsol}
  A0inv::Matrix{Tsol}

  function ApplyEntropyKernel_diagEData{Tsol, Tres}(numDofPerNode::Integer, nd::Integer) where {Tsol, Tres}

    q_avg = zeros(Tsol, numDofPerNode)
    vL = zeros(Tsol, numDofPerNode)
    vR = zeros(Tsol, numDofPerNode)
    F_tmp = zeros(Tres, numDofPerNode)

    # diff  method
    q_avg_dot = zeros(Tsol, numDofPerNode, nd)
    F = zeros(Tres, numDofPerNode)
    F_dot = zeros(Tres, numDofPerNode, nd)

    delta_w_dot = zeros(Tsol, numDofPerNode, nd)

    # rev
    q_avg_bar = zeros(Tres, numDofPerNode)
    delta_w = zeros(Tsol, numDofPerNode)
    delta_w_bar = zeros(Tsol, numDofPerNode)
    A0inv = zeros(Tsol, numDofPerNode, numDofPerNode)

    obj = new(q_avg, q_avg_dot, F, F_dot, vL, vR, F_tmp, delta_w_dot,
              q_avg_bar, delta_w, delta_w_bar, A0inv)

    assertArraysUnique(obj); assertFieldsConcrete(obj)
    return obj
  end
end

"""
  Temporary storage for [`getLambdaMaxSimple`](@ref) and differentiated versions
"""
struct GetLambdaMaxSimpleData{Tsol}
  q_avg::Vector{Tsol}
  q_avg_bar::Vector{Tsol}

  function GetLambdaMaxSimpleData{Tsol}(numDofPerNode::Integer) where {Tsol}

    q_avg = zeros(Tsol, numDofPerNode)
    q_avg_bar = zeros(Tsol, numDofPerNode)

    obj = new(q_avg, q_avg_bar)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Temporary storage for [`getLambdaMax`](@ref) and differentiated versions.
"""
struct GetLambdaMaxData{Tsol}
  p_dot::Vector{Tsol}

  function GetLambdaMaxData{Tsol}(numDofPerNode::Integer) where {Tsol}

    p_dot = zeros(Tsol, numDofPerNode)

    return new(p_dot)
  end
end

"""
  Temporary storage for [`calc2RWaveSpeed`](@ref)
"""
struct TwoRWaveSpeedData{Tsol, Tres}

  pL_dot::Vector{Tsol}
  pR_dot::Vector{Tsol}
  aL_dot::Vector{Tsol}
  aR_dot::Vector{Tsol}
  u_nrmL_dot::Vector{Tres}
  u_nrmR_dot::Vector{Tres}
  p_tr_dot::Matrix{Tres}
  qfL_dot::Matrix{Tres}
  qfR_dot::Matrix{Tres}

  function TwoRWaveSpeedData{Tsol, Tres}(numDofPerNode::Integer) where {Tsol, Tres}
    pL_dot = zeros(Tsol, numDofPerNode)
    pR_dot = zeros(Tsol, numDofPerNode)
    aL_dot = zeros(Tsol, numDofPerNode)
    aR_dot = zeros(Tsol, numDofPerNode)
    u_nrmL_dot = zeros(Tres, numDofPerNode)
    u_nrmR_dot = zeros(Tres, numDofPerNode)
    p_tr_dot = zeros(Tres, 2, numDofPerNode)
    qfL_dot = zeros(Tres, 2, numDofPerNode)
    qfR_dot = zeros(Tres, 2, numDofPerNode)

    return new(pL_dot, pR_dot, aL_dot, aR_dot, u_nrmL_dot, u_nrmR_dot, p_tr_dot,
               qfL_dot, qfR_dot)
  end
end

struct CalcVorticityData{Tsol, Tres, Tmsh}
  dxidx_unscaled::Array{Tmsh, 3}
  velocities::Matrix{Tsol}
  velocity_deriv::Array{Tmsh, 3}
  velocity_deriv_xy::Array{Tmsh, 3}

  function CalcVorticityData{Tsol, Tres, Tmsh}(numNodesPerElement::Integer, dim::Integer) where {Tsol, Tres, Tmsh}

    dxidx_element = zeros(Tmsh, dim, dim, numNodesPerElement)
    velocities = zeros(Tsol, dim, numNodesPerElement)
    velocity_deriv = zeros(Tsol, dim, numNodesPerElement, dim)
    velocity_deriv_xy = zeros(Tres, dim, dim, numNodesPerElement)

    obj = new(dxidx_element, velocities, velocity_deriv, velocity_deriv_xy)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

struct ContractResEntropyVarsData{Tsol}
  w_vals::Vector{Tsol}

  function ContractResEntropyVarsData{Tsol}(numDofPerNode::Integer) where {Tsol}

    w_vals = zeros(Tsol, numDofPerNode)

    return new(w_vals)
  end
end

"""
  Temporary storage for[`getTau`](@ref)
"""
struct GetTauData{Tsol, Tres}
  AjAk::Matrix{Tsol}
  flux_term::Matrix{Tsol}
  A0inv::Matrix{Tsol}
  tmp_mat::Matrix{Tres}

  # second method
  B_d::Matrix{Tres}
  B_p::Matrix{Tres}
  A_mat_hat::Array{Tsol, 3}
  A0::Matrix{Tsol}
  tmp_mat2::Matrix{Tsol}

  function GetTauData{Tsol, Tres}(numDofPerNode::Integer, dim::Integer) where {Tsol, Tres}

    AjAk = zeros(Tsol, numDofPerNode, numDofPerNode)
    flux_term = zeros(Tsol, numDofPerNode, numDofPerNode)
    A0inv = zeros(Tsol, numDofPerNode, numDofPerNode)
    tmp_mat = zeros(Tres, numDofPerNode, numDofPerNode)

    B_d = zeros(Tres, numDofPerNode, numDofPerNode)
    B_p = zeros(Tres, numDofPerNode, numDofPerNode)
    A_mat_hat = zeros(Tsol, numDofPerNode, numDofPerNode, dim)
    A0 = zeros(Tsol, numDofPerNode, numDofPerNode)
    tmp_mat2 = zeros(Tsol, numDofPerNode, numDofPerNode)

    obj = new(AjAk, flux_term, A0inv, tmp_mat,
              B_d, B_p, A_mat_hat, A0, tmp_mat2)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

#------------------------------------------------------------------------------
# structs for face element integrals

"""
  Temporary storage for [`calcECFaceIntegral`](@ref)
"""
struct CalcECFaceIntegralData{Tres, Tmsh}
  fluxD::Matrix{Tres}
  nrmD::Matrix{Tmsh}

  # for SparseFace method
  flux_tmp::Vector{Tres}

  # diff method
  fL_dotD::Array{Tres, 3}
  fR_dotD::Array{Tres, 3}
  fL_dot::Array{Tres, 2}
  fR_dot::Array{Tres, 2}

  # rev methods
  fluxD_bar::Matrix{Tres}

  function CalcECFaceIntegralData{Tres, Tmsh}(numDofPerNode::Integer, dim::Integer) where {Tres, Tmsh}
    fluxD = zeros(Tres, numDofPerNode, dim)
    nrmD = zeros(Tmsh, dim, dim)
    flux_tmp = zeros(Tres, numDofPerNode)

    fL_dotD = zeros(Tres, numDofPerNode, numDofPerNode, dim)
    fR_dotD = zeros(Tres, numDofPerNode, numDofPerNode, dim)

    fL_dot = zeros(Tres, numDofPerNode, numDofPerNode)
    fR_dot = zeros(Tres, numDofPerNode, numDofPerNode)

    fluxD_bar = zeros(Tres, numDofPerNode, dim)

    obj = new(fluxD, nrmD, flux_tmp, fL_dotD, fR_dotD, fL_dot, fR_dot,
              fluxD_bar)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

struct CalcEntropyPenaltyIntegralData{Tsol, Tres}
  wL::Matrix{Tsol}
  wR::Matrix{Tsol}
  wL_i::Vector{Tsol}
  wR_i::Vector{Tsol}
  qL_i::Vector{Tsol}
  qR_i::Vector{Tsol}
  flux::Vector{Tres}
  A0::Matrix{Tsol}

  # sparseface method
  delta_w::Vector{Tsol}
  q_avg::Vector{Tsol}
  res_vals::Vector{Tres}

  #------------------
  # diff method
  q_avg_dot::Matrix{Tsol}
  delta_w_dot::Matrix{Tsol}
  flux_dot_i::Matrix{Tres}  # used to compute the flux

  # arrays in the format required by SBP
  flux_dotL::Array{Tres, 3}
  flux_dotR::Array{Tres, 3}

  jacLL_tmp::Array{Tres, 4}
  jacLR_tmp::Array{Tres, 4}
  jacRL_tmp::Array{Tres, 4}
  jacRR_tmp::Array{Tres, 4}

  A0invL::Matrix{Tsol}
  A0invR::Matrix{Tsol}

  #---------------
  # rev methods
  flux_bar::Vector{Tres}
  qL_bar_i::Vector{Tres}
  qR_bar_i::Vector{Tres}
  wL_bar_i::Vector{Tres}
  wR_bar_i::Vector{Tres}
  delta_w_bar::Vector{Tres}
  q_avg_bar::Vector{Tres}
  wL_bar::Matrix{Tres}
  wR_bar::Matrix{Tres}

  function CalcEntropyPenaltyIntegralData{Tsol, Tres}(numDofPerNode::Integer,
                              numNodesPerFace::Integer,
                              stencilsize::Integer, numNodesPerElement::Integer,
                              nd::Integer) where {Tsol, Tres}

    wL = zeros(Tsol, numDofPerNode, stencilsize)
    wR = zeros(Tsol, numDofPerNode, stencilsize)
    wL_i = zeros(Tsol, numDofPerNode)
    wR_i = zeros(Tsol, numDofPerNode)
    qL_i = zeros(Tsol, numDofPerNode)
    qR_i = zeros(Tsol, numDofPerNode)
    flux = zeros(Tres, numDofPerNode)
    A0 = zeros(Tsol, numDofPerNode, numDofPerNode)

    delta_w = zeros(Tsol, numDofPerNode)
    q_avg = zeros(Tsol, numDofPerNode)
    res_vals = zeros(Tres, numDofPerNode)

    # diff method
    q_avg_dot = zeros(Tsol, numDofPerNode, nd)
    delta_w_dot = zeros(Tsol, numDofPerNode, nd)
    flux_dot_i = zeros(Tres, numDofPerNode, nd)

    # arrays needed by SBP function
    flux_dotL = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerFace)
    flux_dotR = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerFace)

    jacLL_tmp = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement, numNodesPerElement)
    jacLR_tmp = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement, numNodesPerElement)
    jacRL_tmp = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement, numNodesPerElement)
    jacRR_tmp = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement, numNodesPerElement)

    A0invL = zeros(Tsol, numDofPerNode, numDofPerNode)
    A0invR = zeros(Tsol, numDofPerNode, numDofPerNode)

    # rev methods
    flux_bar = zeros(Tres, numDofPerNode)

    qL_bar_i = zeros(Tres, numDofPerNode)
    qR_bar_i = zeros(Tres, numDofPerNode)
    wL_bar_i = zeros(Tres, numDofPerNode)
    wR_bar_i = zeros(Tres, numDofPerNode)
    delta_w_bar = zeros(Tres, numDofPerNode)
    q_avg_bar = zeros(Tres, numDofPerNode)
    wL_bar = zeros(Tres, numDofPerNode, stencilsize)
    wR_bar = zeros(Tres, numDofPerNode, stencilsize)

    obj = new(wL, wR, wL_i, wR_i, qL_i, qR_i, flux, A0, delta_w, q_avg,
              res_vals,
             q_avg_dot, delta_w_dot, flux_dot_i, flux_dotL, flux_dotR,
             jacLL_tmp, jacLR_tmp, jacRL_tmp, jacRR_tmp, A0invL, A0invR,
             flux_bar, qL_bar_i, qR_bar_i, wL_bar_i, wR_bar_i, delta_w_bar,
             q_avg_bar, wL_bar, wR_bar)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Temporary storage for [`interpolateElementStaggered`](@ref)
"""
struct InterpolateElementStaggeredData{Tsol}
  wvars_s::Matrix{Tsol}
  wvars_f::Matrix{Tsol}

  function InterpolateElementStaggeredData{Tsol}(numDofPerNode::Integer,
              numNodesPerElement_s::Integer, numNodesPerElement_f) where {Tsol}

    wvars_s = zeros(Tsol, numDofPerNode, numNodesPerElement_s)
    wvars_f = zeros(Tsol, numDofPerNode, numNodesPerElement_f)

    obj = new(wvars_s, wvars_f)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Temporary storage for [`calcEulerFlux`](@ref) and differentiated versions
"""
struct CalcEulerFluxData{Tsol}
  p_dot::Vector{Tsol}

  function CalcEulerFluxData{Tsol}(numDofPerNode::Integer) where {Tsol}

    p_dot = zeros(Tsol, numDofPerNode)

    obj = new(p_dot)

    return obj
  end
end

"""
  Temporary storage for all BC functors.  This one is a bit odd, because
  we can't store these arrays in the BC functors themselves, because when
  doing boundary conditions for explicitly computed jacobians, the
  Tsol of the eqn object and the type of the solution the boundary condition
  functor operates on are different.

  This design needs to be rethought.
"""
struct BCData{Tsol, Tres}

  qg::Vector{Tsol}
  v_vals::Vector{Tsol}

  sat::Vector{Tres}
  dq::Vector{Tsol}
  roe_vars::Vector{Tsol}
  euler_flux::Vector{Tres}

  # reverse mode
  q_bar::Vector{Tres}
  qg_bar::Vector{Tres}

  function BCData{Tsol, Tres}(numDofPerNode::Integer) where {Tsol, Tres}

    qg = zeros(Tsol, numDofPerNode)
    v_vals = zeros(Tsol, numDofPerNode)
    sat = zeros(Tsol, numDofPerNode)
    dq = zeros(Tsol, numDofPerNode)
    roe_vars = zeros(Tsol, numDofPerNode)
    euler_flux = zeros(Tres, numDofPerNode)

    # reverse mode
    q_bar = zeros(Tres, numDofPerNode)
    qg_bar = zeros(Tres, numDofPerNode)

    obj = new(qg, v_vals, sat, dq, roe_vars, euler_flux,
              q_bar, qg_bar)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

#------------------------------------------------------------------------------
# entropy kernel functions for entropy-stable scheme

struct LW2Kernel{Tsol, Tres, Tmsh} <: AbstractEntropyKernel
  nrm::Array{Tmsh, 1}
  P::Array{Tmsh, 2}
  Y::Array{Tsol, 2}  # eigenvectors
  Lambda::Array{Tsol, 1}  # eigenvalues
  S2::Array{Tsol, 1}  # scaling for the eigensystem
  q_tmp::Array{Tsol, 1}
  tmp1::Array{Tres, 1}
  tmp2::Array{Tres, 1}

  function LW2Kernel{Tsol, Tres, Tmsh}(numDofPerNode::Integer, dim::Integer) where {Tsol, Tres, Tmsh}
    nrm = zeros(Tmsh, dim)
    P = zeros(Tmsh, numDofPerNode, numDofPerNode)
    Y = zeros(Tsol, numDofPerNode, numDofPerNode)
    Lambda = zeros(Tsol, numDofPerNode)
    S2 = zeros(Tsol, numDofPerNode)
    q_tmp = zeros(Tsol, numDofPerNode)
    tmp1 = zeros(Tres, numDofPerNode)
    tmp2 = zeros(Tres, numDofPerNode)

    obj =  new(nrm, P, Y, Lambda, S2, q_tmp, tmp1, tmp2)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Data for Lax-Friedrich entropy dissipation
"""
mutable struct LFKernel{Tsol, Tres, Tmsh} <: AbstractEntropyKernel
  A0::Array{Tsol, 2}

  # diff method
  A0_dot::Array{Tsol, 3}
  t1::Array{Tsol, 1}
  t1_dot::Array{Tsol, 2}
  lambda_max_dot::Array{Tres, 1}

  # revm method
  t1_bar::Array{Tres, 1}
  A0_bar::Array{Tres, 2}

  function LFKernel{Tsol, Tres, Tmsh}(numDofPerNode::Integer, nd::Integer) where {Tsol, Tres, Tmsh}

    A0 = zeros(Tsol, numDofPerNode, numDofPerNode)

    # diff method
    A0_dot = zeros(Tsol, numDofPerNode, numDofPerNode, nd)
    t1 = zeros(Tsol, numDofPerNode)
    t1_dot = zeros(Tsol, numDofPerNode, nd)
    lambda_max_dot = zeros(Tres, numDofPerNode)

    # revq method
    t1_bar = zeros(Tres, numDofPerNode)
    A0_bar = zeros(Tres, numDofPerNode, numDofPerNode)

    obj = new(A0, A0_dot, t1, t1_dot, lambda_max_dot, t1_bar, A0_bar)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

mutable struct IdentityKernel{Tsol, Tres, Tmsh} <: AbstractEntropyKernel

  function IdentityKernel{Tsol, Tres, Tmsh}() where {Tsol, Tres, Tmsh}

    return new()
  end
end

struct GetIRA0Data{Tsol}
  p_dot::Array{Tsol, 1}

  function GetIRA0Data{Tsol}(numDofPerNode::Integer) where {Tsol}

    p_dot = zeros(Tsol, numDofPerNode)

    return new(p_dot)
  end
end

#------------------------------------------------------------------------------
# structs for functions that loop over faces

"""
  Temporary storage for [`getFaceElementIntegral`](@ref) and related functions
"""
struct FaceElementIntegralData{Tsol, Tres}
  # staggered grid things
  resL_s::Matrix{Tres}
  resR_s::Matrix{Tres}
  resL_f::Matrix{Tres}
  resR_f::Matrix{Tres}
  qvars_f::Matrix{Tsol}

  function FaceElementIntegralData{Tsol, Tres}(numDofPerNode::Integer,
                                         numNodesPerElement_s::Integer,
                                         numNodesPerElement_f::Integer) where {Tsol, Tres}

    resL_s = zeros(Tres, numDofPerNode, numNodesPerElement_s)
    resR_s = zeros(Tres, numDofPerNode, numNodesPerElement_s)

    resL_f = zeros(Tres, numDofPerNode, numNodesPerElement_f)
    resR_f = zeros(Tres, numDofPerNode, numNodesPerElement_f)
    qvars_f = zeros(Tsol, numDofPerNode, numNodesPerElement_f)

    obj = new(resL_s, resR_s, resL_f, resR_f, qvars_f)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Temporary storage for [`calcFaceIntegrals`](@ref) and related functions
  (including shared face integral functions
"""
struct CalcFaceIntegralsData{Tsol, Tres}
  q_faceL::Matrix{Tsol}
  q_faceR::Matrix{Tsol}

  # diff functions
  flux_dotL::Array{Tres, 3}
  flux_dotR::Array{Tsol, 3}
  res_jacLL::Array{Tres, 4}
  res_jacLR::Array{Tres, 4}
  res_jacRL::Array{Tres, 4}
  res_jacRR::Array{Tres, 4}

  function CalcFaceIntegralsData{Tsol, Tres}(numDofPerNode::Integer,
                    numNodesPerFace::Integer, numNodesPerElement::Integer) where {Tsol, Tres}

    q_faceL = zeros(Tsol, numDofPerNode, numNodesPerFace)
    q_faceR = zeros(Tsol, numDofPerNode, numNodesPerFace)

    flux_dotL = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerFace)
    flux_dotR = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerFace)

    res_jacLL = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement,
                            numNodesPerElement)
    res_jacLR = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement,
                            numNodesPerElement)
    res_jacRL = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement,
                            numNodesPerElement)
    res_jacRR = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement,
                            numNodesPerElement)

    obj = new(q_faceL, q_faceR, flux_dotL, flux_dotR, res_jacLL, res_jacLR, res_jacRL, res_jacRR)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

"""
  Temporary storage for [`calcVolumeIntegrals`](@ref) and
  [`calcVolumeIntegralsSplitForm`](@ref).
"""
struct CalcVolumeIntegralsData{Tres, Tmsh}

  # calcVolumeIntegrals
  flux_jac::Array{Tres, 4}
  res_jac::Array{Tres, 4}
  nrm::Vector{Tmsh}

  # calcVolumeIntegralsSplitFormLinear
  nrmD::Matrix{Tmsh}
  F_d::Matrix{Tres}
  nrmD_bar::Matrix{Tmsh}
  F_d_bar::Matrix{Tres}

  S::Array{Float64, 3}  # SBP S matrix (should elementtype really be Float64?)

  # calcVolumeIntegralsSplitFormCurvilinear
  Sx::Array{Tmsh, 3}
  Sx_bar::Array{Tmsh, 3}

  # staggered grid methods
  res_s::Matrix{Tres}
  res_f::Matrix{Tres}

  function CalcVolumeIntegralsData{Tres, Tmsh}(numDofPerNode::Integer,
                  dim::Integer, numNodesPerElement_s::Integer,
                  numNodesPerElement_f::Integer, sbp::AbstractOperator) where {Tres, Tmsh}

    # calcVolumeIntegrals
    flux_jac = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement_s, dim)
    res_jac = zeros(Tres, numDofPerNode, numDofPerNode, numNodesPerElement_s, numNodesPerElement_s)
    nrm = zeros(Tmsh, dim)

    # split form stuff
    nrmD = zeros(Tmsh, dim, dim)
    F_d = zeros(Tres, numDofPerNode, dim)
    nrmD_bar = zeros(Tmsh, dim, dim)
    F_d_bar = zeros(Tres, numDofPerNode, dim)

    # the staggered grid calculation only uses the curvilinear method, so
    # make S numNodesPerElement_s, to match the sbp operator on the solution
    # grid (S is not used if they don't match>
    S = zeros(Float64, numNodesPerElement_s, numNodesPerElement_s, dim)
    Sx = zeros(Tmsh, numNodesPerElement_f, numNodesPerElement_f, dim)
    Sx_bar = zeros(Tmsh, numNodesPerElement_f, numNodesPerElement_f, dim)

    for i=1:dim
      S[:, :, i] = 0.5*(sbp.Q[:, :, i] - sbp.Q[:, :, i].')
    end

    res_s = zeros(Tres, numDofPerNode, numNodesPerElement_s)
    res_f = zeros(Tres, numDofPerNode, numNodesPerElement_f)

    obj = new(flux_jac, res_jac, nrm,
              nrmD, F_d, nrmD_bar, F_d_bar,  S, Sx, Sx_bar, res_s, res_f)

    assertArraysUnique(obj); assertFieldsConcrete(obj)

    return obj
  end
end

#------------------------------------------------------------------------------
# Stabilization structs

"""
  Data for Local Projection Stabiization (LPS).

  **Fields**

   * P: the projection matrix that captures the high frequency modes
        (I - L L^T H) in the paper
   * alpha: coefficient in front of the stabilization term
   * entropy_vars: specifies what variables to apply the stabiization to.
                   This is an abstractly-typed field, so it should be access
                   through a function barrier.

"""
struct LPSData{Tsol, Tres}
  P::Matrix{Float64}
  alpha::Float64
  entropy_vars::AbstractVariables

  t1::Matrix{Tsol}
  t2::Matrix{Tres}
  A0::Matrix{Tsol}
  dir::Vector{Tres}

  t1_jac::Array{Tsol, 3}
  t2_jac::Array{Tsol, 3}
  t3_jac::Array{Tres, 3}
  A0_diff::Array{Tsol, 3}
  q_dot::Matrix{Float64}
  lambdas::Vector{Tres}
  lambdas_jac::Matrix{Tres}
  lambda_dot_p::Vector{Tres}

  t1_bar::Matrix{Tres}
  t2_bar::Matrix{Tres}
  tmp::Vector{Tres}
  tmp_bar::Vector{Tres}
  A0_bar::Matrix{Tres}

  dir_bar::Vector{Tres}
  function LPSData{Tsol, Tres}(mesh, sbp, opts) where {Tsol, Tres}

    @unpack mesh numDofPerNode numNodesPerElement dim
    P = getLPSMatrix(sbp)
    alpha = 1.0
    #entropy_vars = ConservativeVariables()
    entropy_vars = IRVariables()

    t1 = zeros(Tsol, numDofPerNode, numNodesPerElement)
    t2 = zeros(Tres, numDofPerNode, numNodesPerElement)
    A0 = zeros(Tsol, numDofPerNode, numDofPerNode)
    dir = zeros(Tres, dim)

    t1_jac = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerElement)
    t2_jac = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerElement)
    t3_jac = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerElement)
    A0_diff = zeros(Tsol, numDofPerNode, numDofPerNode, numDofPerNode)
    q_dot = zeros(numDofPerNode, numDofPerNode)
    for i=1:numDofPerNode
      q_dot[i, i] = 1
    end

    lambdas = zeros(Tres, numNodesPerElement)
    lambdas_jac = zeros(Tres, numDofPerNode, numNodesPerElement)
    lambda_dot_p = zeros(Tres, numDofPerNode)

    t1_bar = zeros(Tres, numDofPerNode, numNodesPerElement)
    t2_bar = zeros(Tres, numDofPerNode, numNodesPerElement)
    tmp = zeros(Tres, numDofPerNode)
    tmp_bar = zeros(Tres, numDofPerNode)
    A0_bar = zeros(Tres, numDofPerNode, numDofPerNode)

    dir_bar = zeros(Tres, dim)

    return new(P, alpha, entropy_vars,
               t1, t2, A0, dir,
               t1_jac, t2_jac, t3_jac, A0_diff, q_dot, lambdas, lambdas_jac,
               lambda_dot_p,
               t1_bar, t2_bar, tmp, tmp_bar, A0_bar,
               dir_bar)
  end
end


text:
function plotOrbits(Qs,names="";index=0,kwargs...)
    if index != 0
        X = [Q[1,:][index] for Q in Qs]
        Y = [Q[2,:][index] for Q in Qs]
    else
        X = [Q[1,:] for Q in Qs]
        Y = [Q[2,:] for Q in Qs]
    end
    return plot(X,Y,label=hcat(names...);kwargs...)
end

function plotT(ts,Qs,titles;kwargs...)
    p=[]
    for (t,Q,title) in zip(ts,Qs,titles)
        push!(p,plot(t,Q',title=title);kwargs...)
    end
    return plot(p...)
end

plotT(ts,Qs,title::String;kwargs...)=plotT(ts,Qs,title,1:length(ts);kwargs...)
plotT(ts,Qs,title::String,extension;kwargs...)=plotT(ts,Qs,["$title $(extension[i])" for i in 1:length(ts)];kwargs...)

plotT(ts,Qs;kwargs...)=plotT(ts,Qs,"plot ";kwargs...)

println("Finished loading HamiltonPlot")

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
struct GridBox
    L::Float32
    D::Float32
    # Number of cells in x direction
    m::Int
    # Number of cells in y direction
    n::Int
    dx::Float64
    dy::Float64
end
function GridBox(L, D, m::Int, n::Int)
    dx = L / m
    dy = D / n
    GridBox(L, D, m, n, dx, dy)
end
GridBox(L, D, dx::AbstractFloat, dy::AbstractFloat) = GridBox(L, D, convert(Int, L / dx), convert(Int, D / dy))
GridBox(canal::CanalConfig, dx::AbstractFloat, dy::AbstractFloat) = GridBox(canal.L, canal.d, dx, dy)
GridBox(canal::CanalConfig, h::AbstractFloat) = GridBox(canal, h, h)
Base.size(gb::GridBox) = (gp.m, gb.n)

abstract type AbstractGrid{T} <: AbstractArray{T, 2} end

abstract type AbstractGridGhost{T} <: AbstractGrid{T} end
Base.size(g::AbstractGrid) = size(g.mesh)
Base.getindex(g::AbstractGrid, i::Int) = getindex(g.mesh, i)
Base.getindex(g::AbstractGrid, i::Int, j::Int) = getindex(g.mesh, i, j)
Base.setindex!(g::AbstractGrid, v,  i::Int) = setindex!(g.mesh, v, i)
Base.setindex!(g::AbstractGrid, v,  i::Int, j::Int) = setindex!(g.mesh, v, i, j)
Base.copy(g::AbstractGridGhost) = typeof(g)(copy(g.mesh), g.grid, (copy(g.ghost[1]), (copy(g.ghost[2]))), g.ghostdim)

getxs(g::AbstractGrid) = 0.5*g.grid.dx:g.grid.dx:g.grid.L-0.5*g.grid.dx
getys(g::AbstractGrid) = (-g.grid.D + g.grid.dx)*0.5:g.grid.dy:(g.grid.D - g.grid.dx)*0.5

function Base.getindex(g::AbstractGridGhost, i::Int, j::Int)
    try
        getindex(g.mesh, i, j)
    catch e
        if e isa BoundsError
            gd = g.ghostdim
            od = gd == 1 ? 2 : 1
            if e.i[g.ghostdim] == 0
                getindex(g.ghost[1], e.i[od])
            elseif e.i[g.ghostdim] == size(g)[gd] + 1
                getindex(g.ghost[2], e.i[od])
            else
                rethrow(e)
            end
        else
            rethrow(e)
        end
    end
end
function Base.setindex!(g::AbstractGridGhost, v, i::Int, j::Int)
    try
        setindex!(g.mesh, v, i, j)
    catch e
        if e isa BoundsError
            gd = g.ghostdim
            od = gd == 1 ? 2 : 1
            if e.i[g.ghostdim] == 0
                setindex!(g.ghost[1], v, e.i[od])
            elseif e.i[g.ghostdim] == size(g)[gd] + 1
                setindex!(g.ghost[2], v, e.i[od])
            else
                rethrow(e)
            end
        else
            rethrow(e)
        end
    end
end
# struct UGrid{T} <: AbstractGrid{T}
#     mesh::Matrix{T}
#     grid::GridBox
#     # ghost0::Vector{T}
#     # ghostn::Vector{T}
# end
struct UGrid{T} <: AbstractGridGhost{T}
    mesh::Matrix{T}
    grid::GridBox
    ghost::Tuple{Vector{T}, Vector{T}}
    ghostdim::Int
end

# function UGrid{T}(gb::GridBox) where T <: AbstractFloat
#     UGrid(zeros(T, gb.m, gb.n), gb, zeros(T, gb.n), zeros(T, gb.n))
# end
# function UGrid{T}(gb::GridBox) where T <: AbstractFloat
#     UGrid(zeros(T, gb.m + 1, gb.n + 2), gb)
# end
function UGrid{T}(gb::GridBox) where T <: AbstractFloat
    nx = gb.m + 1
    ny = gb.n
    mesh = zeros(T, nx, ny)
    ghost = (zeros(T, nx), zeros(T, nx))
    UGrid(mesh, gb, ghost, 2)
end

UGrid(gb::GridBox) = UGrid{Float64}(gb)
projectionV(g::UGrid, i, j) = 0.25 * (g[i, j] + g[i, j-1] + g[i+1, j - 1] + g[i+1, j])
getxs(g::UGrid) = 0:g.grid.dx:g.grid.L
getys(g::UGrid) = -g.grid.D*0.5 + 0.5*g.grid.dy:g.grid.dy:g.grid.D*0.5 - 0.5*g.grid.dy

inrange(p::UGrid) = 2:size(p)[1]-1, 1:size(p)[2]

# trueview(g::UGrid) = @view g.mesh[:, 2:end-1]

struct VGrid{T} <: AbstractGridGhost{T}
    mesh::Matrix{T}
    grid::GridBox
    ghost::Tuple{Vector{T}, Vector{T}}
    ghostdim::Int
end
function VGrid{T}(gb::GridBox) where T <: AbstractFloat
    nx = gb.m
    ny = gb.n + 1
    mesh = zeros(T, nx, ny)
    ghost = (zeros(T, ny), zeros(T, ny))
    VGrid(mesh, gb, ghost, 1)
end
# struct VGrid{T} <: AbstractGrid{T}
#     mesh::Matrix{T}
#     grid::GridBox
# end
# function VGrid{T}(gb::GridBox) where T <: AbstractFloat
#     VGrid(zeros(T, gb.m + 2, gb.n + 1), gb)
# end
VGrid(gb::GridBox) = VGrid{Float64}(gb)
projectionU(g::VGrid, i, j) = 0.25 * (g[i, j] + g[i, j+1] + g[i-1, j + 1] + g[i-1, j])

getxs(g::VGrid) = 0+g.grid.dx*0.5:g.grid.dx:g.grid.L-g.grid.dx*0.5
getys(g::VGrid) = -g.grid.D*0.5:g.grid.dy:g.grid.D*0.5

inrange(v::VGrid) = 1:size(v)[1], 2:size(v)[2]-1

# trueview(g::VGrid) = @view g.mesh[2:end-1, :]

struct PGrid{T} <: AbstractGrid{T}
    mesh::Matrix{T}
    grid::GridBox
end
function PGrid{T}(gb::GridBox) where T <: AbstractFloat
    PGrid(zeros(T, gb.m, gb.n), gb)
end
PGrid(gb::GridBox) = PGrid{Float64}(gb)

inrange(p::PGrid) = 2:size(p)[1]-1, 2:size(p)[2]-1
# getxs(g::PGrid) = 0+g.grid.dx*0.5:g.grid.dx:g.grid.L-g.grid.dx*0.5
# getys(g::PGrid) = -g.grid.D*0.5+0.5*g.grid.dy:g.grid.dy:g.grid.D*0.5-0.5*g.grid.dy

# trueview(g::PGrid) = @view g.mesh[:, :]

struct VortGrid{T} <: AbstractGrid{T}
    mesh::Matrix{T}
    grid::GridBox
end
function VortGrid{T}(gb::GridBox) where T <: AbstractFloat
    VortGrid(zeros(T, gb.m, gb.n), gb)
end
VortGrid(gb::GridBox) = VortGrid{Float64}(gb)
getxs(g::VortGrid) = 0+g.grid.dx*0.5:g.grid.dx:g.grid.L-g.grid.dx*0.5
getys(g::VortGrid) = -g.grid.D*0.5+0.5*g.grid.dy:g.grid.dy:g.grid.D*0.5-0.5*g.grid.dy

function VortGrid(u::UGrid, v::VGrid)
    vort = VortGrid(u.grid)
    for ei in eachindex(vort)
        i,j = Tuple(ei)
        dvdx = (projectionU(v, i+1, j) - projectionU(v, i, j)) / v.grid.dx
        dudy = (projectionV(u, i, j+1) - projectionV(u, i, j)) / v.grid.dy
        vort[ei] = dvdx - dudy
    end
    vort
end

struct DivGrid{T} <: AbstractGrid{T}
    mesh::Matrix{T}
    grid::GridBox
end
function DivGrid{T}(gb::GridBox) where T <: AbstractFloat
    DivGrid(zeros(T, gb.m, gb.n), gb)
end
DivGrid(gb::GridBox) = DivGrid{Float64}(gb)
getxs(g::DivGrid) = 0+g.grid.dx*0.5:g.grid.dx:g.grid.L-g.grid.dx*0.5
getys(g::DivGrid) = -g.grid.D*0.5+0.5*g.grid.dy:g.grid.dy:g.grid.D*0.5-0.5*g.grid.dy

function DivGrid(u::UGrid, v::VGrid)
    div = DivGrid(u.grid)
    for ei in eachindex(div)
        i,j = Tuple(ei)
        dudx = (u[i+1, j] - u[i, j]) / v.grid.dx
        dvdy = (v[i, j+1] - v[i, j]) / v.grid.dy
        div[ei] = dudx + dvdy
    end
    div
end

# trueview(g::VortGrid) = @view g.mesh[:, :]


text:
module BenchMeanVar
using Random
using BenchmarkTools
using TapirBenchmarks

Random.seed!(1234)

const SUITE = BenchmarkGroup()

for n in [2^20]
    SUITE[:n=>n] = s0 = BenchmarkGroup()
    xs = randn(n)
    s0[:impl=>:seq] = @benchmarkable meanvar_seq($xs)
    s0[:impl=>:threads] = @benchmarkable meanvar_threads($xs)
    s0[:impl=>:tapir] = @benchmarkable meanvar_tapir($xs)
end

end
BenchMeanVar.SUITE


text:
struct VectorvaluedDict{DICTS,S,N,T} <: CompositeDict{S,SVector{N,T}}
    dicts   ::  DICTS
    offsets ::  Vector{Int}

    function VectorvaluedDict{DICTS,S,N,T}(dicts) where {DICTS,S,N,T}
        offsets = compute_offsets(dicts)
        new(dicts, offsets)
    end
end

function VectorvaluedDict(dict::Dictionary1d, ::Val{N}) where {N}
    dicts = ntuple(d->dict, Val{N}())
    S = domaintype(dict)
    T = codomaintype(dict)
    VectorvaluedDict{typeof(dicts),S,N,T}(dicts)
end

unit_vector(k, ::Val{N}, T) where {N} = SVector{N,T}(ntuple( i->i==k, Val{N}()))

function unsafe_eval_element(dict::VectorvaluedDict{DICTS,S,N,T}, idx::MultilinearIndex, x) where {DICTS,S,N,T}
    idx1 = outerindex(idx)
    idx2 = innerindex(idx)
    z = unit_vector(idx1, Val{N}(), T) * unsafe_eval_element(component(dict, idx1), idx2, x)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to parse the table cells in HTML back-end.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _parse_cell_html(cell::T; kwargs...)

Parse the table `cell` of type `T`.

This function must return a string that will be printed to the IO.
"""
@inline function _parse_cell_html(
    cell;
    allow_html_in_cells::Bool = false,
    cell_first_line_only::Bool = false,
    compact_printing::Bool = true,
    limit_printing::Bool = true,
    linebreaks::Bool = false,
    renderer::Union{Val{:print}, Val{:show}} = Val(:print),
    kwargs...
)
    # Create the context that will be used when rendering the cell. Notice that
    # the `IOBuffer` will be neglected.
    context = IOContext(stdout, :compact => compact_printing, :limit => limit_printing)

    # Convert to string using the desired renderer.
    if renderer === Val(:show)
        if cell isa AbstractString
            cell_str = cell
        elseif showable(MIME("text/html"), cell)
            cell_str = sprint(show, MIME("text/html"), cell; context = context)
        else
            cell_str = sprint(show, cell; context = context)
        end
    else
        cell_str = sprint(print, cell; context = context)
    end

    # Check if we need to replace `\n` with `<br>`.
    replace_newline = !cell_first_line_only && linebreaks

    if cell_first_line_only
        cell_str = split(cell_str, '\n')[1]
    end

    # If the user wants HTML code inside cell, then we must not escape the HTML
    # characters.
    return _str_html_escaped(cell_str, replace_newline, !allow_html_in_cells)
end

@inline function _parse_cell_html(cell::Markdown.MD; kwargs...)
    return replace(sprint(show, MIME("text/html"), cell),"\n"=>"")
end

@inline _parse_cell_html(cell::Missing; kwargs...) = "missing"
@inline _parse_cell_html(cell::Nothing; kwargs...) = "nothing"
@inline _parse_cell_html(cell::UndefInitializer; kwargs...) = "#undef"


text:
using DescartesFEA
using Test

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


text:
using NNlib

import LinearAlgebra
import LinearAlgebra.BLAS

for (gemm, elty) in NNlib.gemm_datatype_mappings
    @eval begin

        @inline function unsafe_gemm!(
            transA::Char, transB::Char,
            m::Int, n::Int, k::Int,
            alpha::($elty), ptrA::Ptr{$elty}, lda::Int,
            ptrB::Ptr{$elty}, ldb::Int, beta::($elty),
            ptrC::Ptr{$elty}, ldc::Int)

            ccall((BLAS.@blasfunc($gemm), BLAS.libblas), Nothing,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BLAS.BlasInt}, Ref{BLAS.BlasInt}, Ref{BLAS.BlasInt},
                   Ref{$elty}, Ptr{$elty}, Ref{BLAS.BlasInt},
                   Ptr{$elty}, Ref{BLAS.BlasInt}, Ref{$elty},
                   Ptr{$elty}, Ref{BLAS.BlasInt}),
                  transA, transB, m, n, k,
                  alpha, ptrA, lda,
                  ptrB, ldb, beta,
                  ptrC, ldc)
            return nothing
        end

        @inline function unsafe_gemm_strided_batched!(
            transA::Char, transB::Char,
            m::Int, n::Int, k::Int,
            alpha::($elty), ptrA::Ptr{$elty}, lda::Int, strideA::Int,
            ptrB::Ptr{$elty}, ldb::Int, strideB::Int, beta::($elty),
            ptrC::Ptr{$elty}, ldc::Int, strideC::Int, batchCount::Int)

            strA = strideA * sizeof($elty)
            strB = strideB * sizeof($elty)
            strC = strideC * sizeof($elty)

            for i = 1:batchCount
                unsafe_gemm!(transA, transB, m, n, k,
                             alpha, ptrA, lda,
                             ptrB, ldb, beta,
                             ptrC, ldc)

                ptrA += strA
                ptrB += strB
                ptrC += strC

            end
            return nothing
        end

        @inline function gemm_strided_batched_impl!(
            transA::Char, transB::Char,
            m::Int, n::Int, k::Int,
            alpha::($elty), A::AbstractArray{$elty}, lda::Int, strideA::Int,
            B::AbstractArray{$elty}, ldb::Int, strideB::Int, beta::($elty),
            C::AbstractArray{$elty}, ldc::Int, strideC::Int, batchCount::Int)

            ptrA = pointer(A)
            ptrB = pointer(B)
            ptrC = pointer(C)

            GC.@preserve A B C begin
                unsafe_gemm_strided_batched!(
                    transA, transB, m, n, k,
                    alpha, ptrA, lda, strideA,
                    ptrB, ldb, strideB, beta,
                    ptrC, ldc, strideC, batchCount)
            end
            return C
        end

        @inline function gemm_strided_batched!(
            transA::Char, transB::Char,
            alpha::($elty), A::AbstractArray{$elty, 3},
            B::AbstractArray{$elty, 3}, beta::($elty),
            C::AbstractArray{$elty, 3})

            Base.require_one_based_indexing(A, B, C)
            BLAS.chkstride1(A, B, C)
            @assert size(A, 3) == size(C, 3) || size(A, 3) == 1 "batch size mismatch: A != C"
            @assert size(B, 3) == size(C, 3) || size(B, 3) == 1 "batch size mismatch: B != C"

            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)

            if m != size(C,1) || n != size(C,2) || ka != kb
                throw(DimensionMismatch("A has size ($m,$ka,$(size(A, 3))), B has size ($kb,$n,$(size(B, 3))), C has size $(size(C))"))
            end

            lda = max(1, stride(A,2))
            ldb = max(1, stride(B,2))
            ldc = max(1, stride(C,2))

            strideA = size(A, 3) == 1 ? 0 : stride(A, 3)
            strideB = size(B, 3) == 1 ? 0 : stride(B, 3)
            strideC = stride(C, 3)
            batchCount = size(C, 3)

            gemm_strided_batched_impl!(
                transA, transB, m, n, ka,
                alpha, A, lda, strideA,
                B, ldb, strideB, beta,
                C, ldc, strideC, batchCount)

            return C
        end

        function gemm_strided_batched(
            transA::Char, transB::Char,
            alpha::($elty), A::AbstractArray{$elty, 3},
            B::AbstractArray{$elty, 3})
            C = similar(B, (size(A, transA == 'N' ? 1 : 2), size(B, transB == 'N' ? 2 : 1), max(size(A, 3), size(B, 3))))
            return gemm_strided_batched!(transA, transB, alpha, A, B, zero($elty), C)
        end

        function gemm_strided_batched(
            transA::Char, transB::Char,
            A::AbstractArray{$elty, 3},
            B::AbstractArray{$elty, 3})
            return gemm_strided_batched(transA, transB, one($elty), A, B)
        end

        function gemm_strided_batched!(
            transA::Char, transB::Char,
            alpha::($elty), A::AbstractArray{$elty, N1},
            B::AbstractArray{$elty, N2}, beta::($elty),
            C::AbstractArray{$elty, N3},
            Ai, Aj, Bi, Bj, Ci, Cj) where {N1, N2, N3}

            # (A1, A2, ..., Ai-1, Ai, Ai+1, ..., Aj-1, Aj, ..., An)
            #  |______lda______|  |_______K/M_______|  |__batch__|
            # (B1, B2, ..., Bi-1, Bi, Bi+1, ..., Bj-1, Bj, ..., Bn)
            #  |______ldb______|  |_______K/N_______|  |__batch__|
            # (C1, C2, ..., Ci-1, Ci, Ci+1, ..., Cj-1, Cj, ..., Cn)
            #  |______ldc______|  |_______K/N_______|  |__batch__|

            Base.require_one_based_indexing(A, B, C)
            BLAS.chkstride1(A, B, C)

            sa1, sa2, sa3 = collapsed_size(A, Ai, Aj)
            sb1, sb2, sb3 = collapsed_size(B, Bi, Bj)
            sc1, sc2, sc3 = collapsed_size(C, Ci, Cj)

            @assert sa3 == sc3 || sa3 == 1 "batch size mismatch: A != C"
            @assert sb3 == sc3 || sb3 == 1 "batch size mismatch: B != C"

            m = transA == 'N' ? sa1 : sa2
            ka = transA == 'N' ? sa2 : sa1
            kb = transB == 'N' ? sb1 : sb2
            n = transB == 'N' ? sb2 : sb1

            if m != sc1 || n != sc2 || ka != kb
                throw(DimensionMismatch("A has size ($m,$ka,$sa3), B has size ($kb,$n,$sb3), C has size ($sc1, $sc2, $sc3)"))
            end

            lda = max(1, stride(A, Ai))
            ldb = max(1, stride(B, Bi))
            ldc = max(1, stride(C, Ci))

            strideA = sa3 == 1 ? 0 : stride(A, Aj)
            strideB = sb3 == 1 ? 0 : stride(B, Bj)
            strideC = stride(C, Cj)
            batchCount = sc3

            gemm_strided_batched_impl!(
                transA, transB, m, n, ka,
                alpha, A, lda, strideA,
                B, ldb, strideB, beta,
                C, ldc, strideC, batchCount)

            return C
        end

        function gemm_strided_batched(
            transA::Char, transB::Char,
            alpha::($elty), A::AbstractArray{$elty, N1},
            B::AbstractArray{$elty, N2},
            Ai, Aj, Bi, Bj) where {N1, N2}

            m = noncollapsed_size(A, Ai, Aj, transA == 'N' ? 1 : 2)
            n = noncollapsed_size(B, Bi, Bj, transB == 'N' ? 2 : 1)
            sc3 = collapsed_size(A, Ai, Aj, 3) > collapsed_size(B, Bi, Bj, 3) ?
                noncollapsed_size(A, Ai, Aj, 3) :
                noncollapsed_size(B, Bi, Bj, 3)

            Ci = length(m) + 1
            Cj = Ci + length(n)
            C = similar(B, (m..., n..., sc3...))
            return gemm_strided_batched!(transA, transB, alpha, A, B, zero($elty), C, Ai, Aj, Bi, Bj, Ci, Cj)
        end

        function gemm_strided_batched(
            transA::Char, transB::Char,
            A::AbstractArray{$elty, N1},
            B::AbstractArray{$elty, N2},
            Ai, Aj, Bi, Bj) where {N1, N2}
            return gemm_strided_batched(transA, transB, one($elty), A, B, Ai, Aj, Bi, Bj)
        end

    end
end


text:
#commented psuedocode
using Meshes
using StaticArrays
using LinearAlgebra
using CoordinateTransformations
using Tmatrix

function calculate_Tmatrix_for_any(
    face_list::Vector{Any}, #make this better typed
    point_list::Vector{Vector{Float64}},
    n_max::Int,
    k1::Complex{R},
    k2::Complex{R};
    rotationally_symmetric = false,
    symmetric_about_plane_perpendicular_z = false,
    BigFloat_precision = nothing,
) where {R <: Real}

    #spherical_list = convert_to_spherical(point_list) #format is [r, θ, ϕ], use coordinatetransformations.jl to convert

    #TODO: bundle these into a new struct and make that not break
    r_array, θ_array, ϕ_array = (
        [element[1] for element in spherical_list],
        [element[2] for element in spherical_list],
        [element[3] for element in spherical_list],
    )

    # calculate r and n̂ for the geometry

    n̂_array = point_list #copying so that size is same

    for element in n̂_array
        element = [0, 0, 0]
    end

    for face in face_list
        for (vertex, i) in enumerate(face)
            n_array[vertex] += vertex × face[(i + 1) % size(face)]
        end
    end

    for element in n̂_array
        element = normalize(element) #TODO make sure this works, use normalize!(element) later if at all
    end

    # calculate T-matrix
    k1r_array = k1 .* r_array
    k2r_array = k2 .* r_array
    T = T_matrix(
        n_max,
        k1,
        k2,
        k1r_array,
        k2r_array,
        r_array,
        θ_array,
        ϕ_array,
        n̂_array;
        HDF5_filename = HDF5_filename,
        rotationally_symmetric = rotationally_symmetric,
        symmetric_about_plane_perpendicular_z = symmetric_about_plane_perpendicular_z,
        BigFloat_precision = BigFloat_precision,
    )
    return T
end

point_list = Point(3, 8, 23), Point(1, 2, 3), Point(5, 12, 12), Point(15, 3, 7)

face_list = [connect((1, 2, 3)), connect((1, 2, 4)), connect((1, 3, 4)), connect((2, 3, 4))]

#convert_to_spherical

#TODO make generic over shape
function Spherical(point_coords::Point)
    Spherical(point_coords.coords[1], point_coords.coords[2], point_coords.coords[3])
end

polar_point_list = [Spherical(p) for p in point_list]

n_max = 1
k1_r, k1_i, k2_r, k2_i = 1e5, 1e3, 2e5, 3e3

calculate_Tmatrix_for_spheroid_SeparateRealImag(
    face_list,
    point_list,
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i;
    rotationally_symmetric = false,
    symmetric_about_plane_perpendicular_z = false,
);


text:
# Base16 Decoder
# ==============

struct Base16Decoder <: Codec
    table::CodeTable16
    state::State
    buffer::Buffer
end

Base16Decoder(table::CodeTable16) = Base16Decoder(table, State(), Buffer(1))

"""
    Base16Decoder(;ignore::String=$(escape_string(whitespace)))

Create a base16 decoding codec.

Arguments
---------
- `ignore`: ASCII characters that will be ignored while decoding
"""
function Base16Decoder(;ignore::String=whitespace)
    table = copy(BASE16_UPPER)
    ignorechars!(table, ignore)
    return Base16Decoder(table)
end

const Base16DecoderStream{S} = TranscodingStream{Base16Decoder,S} where S<:IO

"""
    Base16DecoderStream(stream::IO; kwargs...)

Create a base16 decoding stream (see `Base16Decoder` for `kwargs`).
"""
function Base16DecoderStream(stream::IO; kwargs...)
    return TranscodingStream(Base16Decoder(;kwargs...), stream)
end

function TranscodingStreams.startproc(
        codec :: Base16Decoder,
        state :: Symbol,
        error :: Error)
    start!(codec.state)
    empty!(codec.buffer)
    return :ok
end

function TranscodingStreams.process(
        codec  :: Base16Decoder,
        input  :: Memory,
        output :: Memory,
        error  :: Error)
    table = codec.table
    state = codec.state
    buffer = codec.buffer

    # Check if we can encode data.
    if !is_running(state)
        error[] = ArgumentError("decoding is already finished")
        return 0, 0, :error
    elseif output.size < 1
        # Need more output space.
        return 0, 0, :ok
    end

    # Load the first bytes.
    i = j = 0
    if buffer.size < 1 && i < input.size
        buffer[buffer.size+=1] = input[i+=1]
    end
    c1 = c2 = BASE16_CODEIGN
    if buffer.size ≥ 1
        c1 = decode(table, buffer[1])
    end
    empty!(buffer)

    # Start decoding loop.
    status = :ok
    @inbounds while true
        if c1 > 0x0f || c2 > 0x0f
            i, j, status = decode16_irregular(table, c1, c2, input, i, output, j, error)
        else
            output[j+1] = c1 << 4 | c2
            j += 1
        end
        if i + 2 ≤ input.size && j + 1 ≤ output.size
            c1 = decode(table, input[i+1])
            c2 = decode(table, input[i+2])
            i += 2
        else
            break
        end
    end

    # Epilogue.
    if status == :end || status == :error
        finish!(state)
    end
    return i, j, status
end

# Decode irregular code (e.g. non-alphabet, etc.).
function decode16_irregular(table, c1, c2, input, i, output, j, error)
    # Skip ignored chars.
    while true
        if c1 == BASE16_CODEIGN
            c1 = c2
        elseif c2 == BASE16_CODEIGN
            # pass
        else
            break
        end
        if i + 1 ≤ input.size
            c2 = decode(table, input[i+=1])
        else
            c2 = BASE16_CODEEND
            break
        end
    end

    # Write output.
    if c1 ≤ 0x0f && c2 ≤ 0x0f
        output[j+=1] = c1 << 4 | c2
        status = :ok
    elseif c1 == BASE16_CODEIGN && c2 == BASE16_CODEEND
        status = :end
    else
        error[] = DecodeError("invalid base16 data")
        status = :error
    end
    return i, j, status
end


text:
"""
    parse_node(node)

take a node and parse into Symbolics form
"""
function parse_node(node)
    tagmap[node.name](node)
end

function parse_str(str)
    doc = parsexml(str)
    parse_doc(doc)
end

function parse_doc(doc)
    node = doc.root
    parse_node(node)
end

function parse_file(fn)
    node = readxml(fn).root
    parse_node(node)
end

"""
    parse_cn(node)

parse a <cn> node
"""
function parse_cn(node)
    # elements(node) != EzXML.Node[] && error("node cant have elements rn, check if </sep> is in the node")
    # this isnt great might want to check that `sep` actually shows up
    if haskey(node, "type") && elements(node) != EzXML.Node[]
        parse_cn_w_sep(node)
    else
        Float64(Meta.parse(node.content))   # convert to Float64 for CellML compatibility
    end
end

"""
    parse_cn_w_sep(node)

parse a <cn type=".."> node

where type ∈ ["e-notation", "rational", "complex-cartesian", "complex-polar"]
"""
function parse_cn_w_sep(node)
    # node = clean_attributes(node)
    # txts = findall("text()", node)
    txts = [n for n in eachnode(node) if istext(n)]
    length(txts) != 2 && error("stop, collaborate, and listen!, problem with <cn>")
    x1, x2 = map(x -> Meta.parse(x.content), txts)
    t = node["type"]
    if t == "e-notation"
        x1 * exp10(x2)
    elseif t == "rational"
        Rational(x1, x2)
    elseif t == "complex-cartesian"
        Complex(x1, x2)
    elseif t == "complex-polar"
        x1 * exp(x2 * im)
    else
        error("$t in parse_cn_w_sep, somethings wrong")
    end
end

"""
    parse_ci(node)

parse a <ci> node
"""
function parse_ci(node)
    c = Meta.parse(strip(node.content))
    Num(Variable(c))
end

########## Parse piecewise ###################################################

function parse_piecewise(node)
    ns = elements(node)
    pieces = filter(x -> nodename(x) == "piece", ns)
    others = filter(x -> nodename(x) == "otherwise", ns)

    if length(others) > 0
        otherwise = parse_node(firstelement(others[1]))
    else
        otherwise = 0.0
    end

    return process_pieces(pieces, otherwise)
end

function process_pieces(pieces, otherwise)
    node = pieces[1]
    c = parse_node.(elements(node))
    return IfElse.ifelse(c[2] > 0.5, c[1],
                length(pieces)==1 ? otherwise : process_pieces(pieces[2:end], otherwise))
end

"""
    parse_apply(node)

parse an <apply> node into Symbolics form

how to deal w apply within apply, need to ensure we've hit bottom
"""
function parse_apply(node)
    node.name != "apply" && error("calling parse_apply requires the name of the element to be `apply`")
    elms = elements(node)
    cs = parse_node.(elms[2:end])
    applymap[elms[1].name](cs)
end

"""
    parse_bvar(node)

parse a <bvar> node
"""
function parse_bvar(node)
    es = elements(node)
    length(es) == 1 ? (parse_node(es[1]), 1) : Tuple(parse_node.(es))
end

"""
    parse_diff(x)

parse a <diff>

"""
function parse_diff(a)
    (iv, deg), x = a
    deg = trunc(Int, deg)
    # num = Num(Symbolics.Sym{Symbolics.FnType{Tuple{Real},Real}}(Symbol(x))(iv))
    D = Differential(iv)^deg
    D(x)
end

"""
    parse_lambda(node)

parse a <lambda> node

```xml
<lambda>
  <bvar> x1 </bvar><bvar> xn </bvar>
   expression-in-x1-xn
</lambda>
```
"""
function parse_lambda(node)
    es = elements(node)
    vars = findall("//x:bvar | //bvar", node, ["x" => MathML.mathml_ns])
    # vars2 = findall("//bvar", node) # works in tests
    # vars = union(vars, vars2) # FIX

    args = first.(parse_bvar.(vars))
    num = parse_apply(es[end])
    eval(build_function([num], args...)[1])
end


text:
######################################################################################
#
# Videoprocessor.jl
#
# OpenCV.jl videoplayer with trackbar control for interactive image processing
#
# Image processing functions
# BRIGHTNESS
# CONTRAST
# THRESHOLD
#
# Easy to add new filters/algorithms for live image processing
#
######################################################################################

# procesess
const BRIGHTNESS = cint(0)
const CONTRAST = cint(1)
const THRESHOLD = cint(2)
const GAUSSIAN = cint(3)

cxx"""
class Videoprocessor
{
    public:

        int framesToProcess = -1; // infinity
        bool stop = false;
        unsigned int delay = 0;  // microseconds

        class Threshold
        {
        public:
        int threshold = 120;
        int max_val = 255;
        int flag = cv::THRESH_BINARY;
              // THRESH_OTSU, THRESH_BINARY_INV,
              // THRESH_TRUNC, THRESH_TOZERO,
              // THRESH_TOZERO_INV
        std::string tn = "threshold";
        };

        class Brightness
        {
          public:
              int brightness = 50;
              int max_val = 100;
              std::string tn = "brightness";
        };

        class Contrast
        {
          public:
              int contrast = 50;
              int max_val = 100;
              std::string tn = "contrast";
        };

        // Instantiation of processing classses
        Threshold Thresh;
        int T = Thresh.threshold;
        Contrast Contr;
        int C = Contr.contrast;
        Brightness Bright;
        int B = Bright.brightness;

        // Videoplayer class methods
        void setDuration(int frames);
        void setDelay(int tsleep);
        void setParameters(int process, int param, int  max_, int flags);
        int getParameters(int process);
        void setTrackbar(int process, const std::string& w);
        cv::Mat ProcessMulti(int process, cv::Mat src, cv::Mat dst);
        cv::Mat ProcessSingle(int process, cv::Mat src, cv::Mat dst);
        int Display(std::vector<int> processes, const std::string& w, const std::string& filename, bool video, int device);
};

void Videoprocessor::setDuration(int frames)
{
     framesToProcess = frames;
}

void Videoprocessor::setDelay(int tsleep)
{
     delay = tsleep; // microseconds
}

void Videoprocessor::setParameters(int process, int  param, int  max_, int flags)
{

    switch(process)
    {

        case(0) :
            B = param;
            Bright.max_val = max_;
            break;
        case(1)  :
            C = param;
            Contr.max_val = max_;
            break;
        case(2) :
            T = param;
            Thresh.max_val = max_;
            Thresh.flag = flags;
            break;
        default :
            std::cout <<"ERROR: wrong input" << std::endl;
            break;
    }
}

int Videoprocessor::getParameters(int process)
{
      switch(process)
      {

        case(0) :
            return(B);
            break;
        case(1) :
            return(C);
            break;
        case(2) :
            return(T);
            break;
        default :
            std::cout <<"ERROR: wrong input" << std::endl;
            return(-1);
            break;
     }
}

void Videoprocessor::setTrackbar(int process, const std::string& w)
{
    //std::string w = wn;
    cv::namedWindow(w, cv::WINDOW_AUTOSIZE);

    switch(process)
    {

        case(0) :
            cv::createTrackbar(Bright.tn, w, &B, Bright.max_val);
            break;
        case(1) :
            cv::createTrackbar(Contr.tn, w, &C, Contr.max_val);
            break;
        case(2) :
            cv::createTrackbar(Thresh.tn, w, &T, Thresh.max_val);
            break;
        default :
            std::cout <<"ERROR: wrong input" << std::endl;
            break;
     }
}

cv::Mat Videoprocessor::ProcessMulti(int process, cv::Mat src, cv::Mat dst)
{

   switch(process)
   {
        case(0) :
            src.convertTo(dst, -1, 1, (B -50));
            return(dst);
            break;
        case(1) :
            src.convertTo(dst, -1, (C/50), 0);
            return(dst);
            break;
        case(2) :
            // Convert the image to Gray
            cv::cvtColor(src, dst, cv::COLOR_RGB2GRAY);
            // Gaussian filtering
            cv::GaussianBlur(dst, dst, cv::Size(5,5), 0);
            cv::threshold(dst, dst, T, Thresh.max_val, Thresh.flag);
            return(dst);
            break;
        default :
            std::cout <<"Unprocessed image!" << std::endl;
            return(dst);
            break;
    }
}

cv::Mat Videoprocessor::ProcessSingle(int process, cv::Mat src, cv::Mat dst)
{
   switch(process)
   {
        case(0) :
            src.convertTo(dst, -1, 1, (B -50));
            return(dst);
            break;
        case(1) :
            src.convertTo(dst, -1, (C/50), 0);
            return(dst);
            break;
        case(2) :
            cv::threshold(src, dst, T, Thresh.max_val, Thresh.flag);
            return(dst);
            break;
        case(3) :
             // Convert the image to Gray
            cv::cvtColor(src, dst, cv::COLOR_RGB2GRAY);
            // Gaussian filtering
            cv::GaussianBlur(dst, dst, cv::Size(5,5), 0);
            return(dst);
            break;
        default :
            std::cout <<"Unprocessed image!" << std::endl;
            return(dst);
            break;
   }
}

int Videoprocessor::Display(std::vector<int> processes, const std::string& w, const std::string& filename, bool video = true, int device = cv::CAP_ANY)
{

     cv::VideoCapture cam = cv::VideoCapture();
     cv::Mat src;

     // Open video file only if filename is provided
     if (video)
     {
         if (filename.size()==0)
         {
              // open Video device
              cam.open(device);

              //Check that capture device is open
              if (!cam.isOpened())
              {
                  std::cout<< "Can not open camera!";
                  return(-1);  // error
              }
         }
         else
         {
              // open Video file;
              cam.open(filename);
              //Check that video file is open
              if (!cam.isOpened())
              {
                  std::cout<< "Can not open video file!";
                  return(-2);  // error
              }
         }

         // Set the frame size to the camera capture WIDTH and HEIGHT
         double dWidth = cam.get(cv::CAP_PROP_FRAME_WIDTH);
         double dHeight = cam.get(cv::CAP_PROP_FRAME_HEIGHT);
         int Height = static_cast<int>(dHeight);
         int Width = static_cast<int>(dWidth);
         // Initialize the cv::Mat structures
         src = cv::Mat(Height, Width, CV_8UC3);
    }

    else   // Open and display single image
    {
         src = cv::imread(filename);

         if (src.empty())
         {
             std::cout<< "Cannot read the file!";
             return(-3);  // error
         }

         // Check that input image is RGB
         if (src.channels() != 3)
         {
            std::cout<< "Color image is required!";
            return(-4);  // error
         }
    }

    // Main processing

    // create a destination image same as src
    cv::Mat dst = cv::Mat(src.rows,src.cols, src.type());

    if (video) // video file/webstream/camera
    {
       int count = 0;  // keep a frame counter

       while(true)
       {
          count +=1;

          bool Success = cam.read(src);
          if (!Success)
          {
              break;
          }

          // select image processing operation
          for (std::size_t i = 0; i < processes.size(); i++)
          {
              src = Videoprocessor::ProcessMulti(processes[i], src, dst);
          }

          cv::imshow(w, src);

          if (count == framesToProcess)
          {
              cv::destroyAllWindows();
              break;
          }

	        // delay in microseconds
          if (delay > 0)
          { usleep(delay); }

          if (cv::waitKey(30)==27)
          {
              std::cout << "Processed " << count << " frames." << std::endl;
              cv::destroyAllWindows();
              break;
          }
       }
      cam.release();
    }

    else
    {
         cv::Mat gray = cv::Mat(src.rows,src.cols, CV_8UC1);
         cv::Mat thresh = gray.clone();

         bool converted = false;
         while(true)
         {
              // select image processing operation
              for (std::size_t i = 0; i < processes.size(); i++)
              {
                   if (processes[i] < 2)
                   {
                       dst  = Videoprocessor::ProcessSingle(processes[i], src, dst);
                   }
                   if (processes[i] ==2)
                   {
                        if (!converted)
                        {
                             gray = Videoprocessor::ProcessSingle(3, src, dst);  // make grayscale and apply Gaussian blur
                             converted = true;
                        }
                        thresh = Videoprocessor::ProcessSingle(processes[i], gray, thresh);
                   }
               }

               if (!converted)
               {
                    cv::imshow(w, dst);
               }
               else
               {
                    cv::imshow(w, thresh);
               }

               if (cv::waitKey(30) == 27)
               {
                  cv::destroyAllWindows();
                  break;
               }
         }
     }

    return (0);
}

std::vector<int> RunVideoprocessor(std::vector<int> processes, const std::string& winname, const std::string &filename,
         int frames = -1, int tsleep=0, int brightness = 30, int contrast = 30, int threshold = 120, int max_val = 255,
         int flag = cv::THRESH_BINARY, bool video = true, int device = cv::CAP_ANY)
{
   // Initialize the videoprocessor
   Videoprocessor videoprocessor;

   // Set how many frames to record, default = -1 (continuous loop);
   if (frames > 1)
   {
        videoprocessor.setDuration(frames);
   }

   // Set interframe delay in microseconds for video streams
   if (tsleep > 0)
   {
       videoprocessor.setDelay(tsleep);
   }

   // Set initial filter parameters (threshold, contrast, brightness)
   videoprocessor.setParameters(0, videoprocessor.B, videoprocessor.Bright.max_val, 0);
   videoprocessor.setParameters(1, videoprocessor.C, videoprocessor.Contr.max_val, 0);
   videoprocessor.setParameters(2, videoprocessor.T, videoprocessor.Thresh.max_val, flag);

   // Set trackbar (BRIGHTNESS, CONTRAST, THRESHOLD)
   for (std::size_t i = 0; i < processes.size(); i++)
   {
        videoprocessor.setTrackbar(processes[i], winname);
   }

   // Set display
   videoprocessor.Display(processes, winname, filename, video, device);

   // Get parameters upon termination
   std::vector<int> final_params;
   for (std::size_t i = 0; i < processes.size(); i++)
   {
        final_params.push_back(videoprocessor.getParameters(processes[i]));
   }
   return(final_params);
}
"""

# Instructions:

# 1. Create an std::vector with processing operations [BRIGHTNESS, CONTRAST, THRESHOLD]
# processes = stdvector(cint(0),cint(0))
# # Always add THRESHOLD last!
# stdpush_back(processes, BRIGHTNESS)
# stdpush_back(processes, CONTRAST)
# stdpush_back(processes, THRESHOLD)

# 2. Arguments:
# processes                     std::vector<int> (see above)
# winname                       window name (string)
# filename= Ptr{UInt8}[0][1]    file to open
# frames=-1                     number of frames to record
# delay=0                       delay in microseconds (uses system usleep), e.g., 1000 = 1s
# brightness=30                 initial brightness value
# contrast=30                   initial contrast value
# threshold=120                 initial threshold value
# max_val=255                   maximum value (paired with process, e.g., BRIGHTNESS)
# flag=THRESH_BINARY            flag for thresholding algorithm
# video=true                    if video source is used
# device=CAP_ANY                flag for video device

videoprocessor(processes, winname::String, filename="", frames=-1, delay=0, brightness=30,
    contrast=30, threshold=120, max_val=255, flag=THRESH_BINARY, video=true, device=CAP_ANY) =
    @cxx RunVideoprocessor(processes, pointer(winname), pointer(filename), frames, delay, brightness, contrast,
        threshold, max_val, flag, video, device)


text:
#=
Force calculations according to Pammer OConner 2013
where force magnitude is estimated based on the degree of bending
in the same section of whisker in every frame selected to be high
SNR (Point P)

=#

#=
x, y - whisker coordinates
theta_f - whisker angle (radians)
curv - whisker curvature
ii - index of contact
i_p - index of high SNR point - we should use the second DLC point for this (as long as it isn't past pole)
E - Elastic Modulus of Whisker
I_p - Moment of inertia of whisker at point p I_p = 1.0
=#
function calc_force(x::Array{Float64,1},y::Array{Float64,1},theta_f::Float64,curv::Float64,ii::Int,i_p::Int,curv_0=0.0; I_p = 1.0, E=1.0)

    #curv_0 intrinsic curvature of the whisker
    delta_kappa = curv - curv_0 #Change in curvature

    #x_f, y_f follicle x and y
    x_f = x[end]
    y_f = y[end]

    #x_c, y_c contact x and y
    x_c = x[ii]
    y_c = y[ii]

    #x_p, y_p point_p x and y
    x_p = x[i_p]
    y_p = y[i_p]

    theta_contact = get_theta_contact(x,y,ii)

    (theta_0,r_0) = calc_r_theta(x_f,x_c,y_f,y_c)
    (theta_p,r_p) = calc_r_theta(x_p,x_c,y_p,y_c)

    (F, M_0) = force_moment(delta_kappa,r_p,r_0,theta_p,theta_0,theta_contact,I_p = I_p, E = E)

    (F_ax, F_lat) = decompose_force(F,theta_f,theta_contact)

    (M_0,F_ax,F_lat,F,theta_contact,r_0)
end

function get_theta_contact(x,y,ii,forward_dir=false)
    theta_contact = 0.0
    try
        #We will first look outward from the point of contact. This will fail if there is not >30 pixels of length
        #beyond that point
        c_i = get_contact_ind(x,y,ii,30.0,forward_dir)
        theta_contact = contact_angle(x[ii],x[c_i],y[ii],y[c_i])
    catch
        #if this fails, we will try making a vector from
        c_i = get_contact_ind(x,y,ii,30.0,!forward_dir)
        theta_contact = contact_angle(x[c_i],x[ii],y[c_i],y[ii])
    end
    theta_contact
end

function calc_force(delta_kappa::Real,p::Tuple,c::Tuple,f::Tuple,theta_contact::Real,theta_f::Real,I_p::Real; E = 2600.0)

    (theta_0,r_0) = calc_r_theta(f[1],c[1],f[2],c[2])
    (theta_p,r_p) = calc_r_theta(p[1],c[1],p[2],c[2])

    (F, M_0) = force_moment(delta_kappa,r_p,r_0,theta_p,theta_0,theta_contact,I_p = I_p, E = E)

    (F_ax, F_lat) = decompose_force(F,theta_f,theta_contact)

    (M_0,F_ax,F_lat,F)
end

function calc_force(x::Array{Float64,1},y::Array{Float64,1},pole::Tuple,w_length::Real,i_p_loc=150.0,
    high_snr_p1=100.0,high_snr_p2=200.0)

    theta_f = get_angle(x,y,30.0) #Follicle is ind=1

    #ii - index of contact
    ii=calc_p_dist(x,y,pole[1],pole[2])[2]
    theta_contact = get_theta_contact(x,y,ii,true)

    #i_p - index of high SNR point
    i_p=culm_dist(x,y,i_p_loc,false)

    curv_ip = get_curv_ip(x,y,high_snr_p1,high_snr_p2)

    I_p = moment_of_inertia(i_p_loc,w_length)

    calc_force(curv_ip,(x[i_p],y[i_p]),(x[ii],y[ii]),(x[i_p],y[i_p]),theta_contact,theta_f,I_p)
end

function get_curv_ip(x,y,p1,p2)
    (x_i,y_i) = make_whisker_segment(x,y,p1,p2)
    curv_ip = least_squares_quad_rot(x_i,y_i)
end

#=
For taped beam (cone)
From https://journals.physiology.org/doi/full/10.1152/jn.00707.2006
x - point to calculate moment of interia
r_base - radius of whisker at base. In paper x/L = 232/2 um / 66.2 mm = 0.0017
L - length of whisker (66 mm)
=#
function moment_of_inertia(x, L=66.2, r_base=.116)
    #alpha = pi/4 * (r_base / L)^4
    alpha = pi/4 * (.116 / 66.2) ^ 4
    I = alpha * (L - x) ^ 4 # mm^4
end

#=
Calculate Force and Moment
delta_kappa - change in curvature from intrinsic
r_p - distance from point of contact to point p (pixels)
r_0 - distance from follicle to point of contact (pixels)
theta_p - angle from point_p to point of contact (radians)
theta_0 - angle from follicle to point of contact (radians)
theta_contact - angle that whisker is pointing at point of contact (radians)

E = 2.6 GPa From https://journals.physiology.org/doi/full/10.1152/jn.00707.2006 (2600 N/mm^2)

I_p varies with distance along the whisker for measurement point
if I_p is 1/4 of the way up the whisker, I_p = 4.5e-5
=#
function force_moment(delta_kappa::Real,r_p::Real,r_0::Real,theta_p::Real,
    theta_0::Real,theta_contact::Real; I_p = 4.5e-5, E = 2600.0)

    #=
    Note that F approaches infinity if 1) theta_p and theta_contact are orthogonal (very unusual)
    or 2) if r_p decreases to zero. This is possible if your high SNR point is right at the point of
    contact. Ideally, we should not be using different point Ps through because the moment of inertia and
    modulus are position dependent.
    =#
    F = abs(delta_kappa * E * I_p / (r_p * cos(theta_p - theta_contact)))

    M_0 = r_0 * F * cos(theta_0 - theta_contact)

    (F, M_0)
end

#=
F - Total Contact Force
theta_f - whisker angle (radians)
theta_contact - angle of contact (radians)
=#
function decompose_force(F::Real,theta_f::Real,theta_contact::Real)
    F_ax = F * sin(theta_f - theta_contact)
    F_lat = F * cos(theta_f - theta_contact)
    (F_ax,F_lat)
end

function calc_r_theta(x1::Real,x2::Real,y1::Real,y2::Real)

    theta_p = atan((y2 - y1),(x2 - x1))
    r_p = sqrt((x2 - x1)^2 + (y2 - y1)^2)
    (theta_p,r_p)
end

function contact_angle(x1::Real,x2::Real,y1::Real,y2::Real)
    atan(y2-y1,x2-x1)
end

#=
Finds
=#
function get_contact_ind(x::Array{T,1},y::Array{T,1},ii,s_thres = 10.0,forward_dir=true) where T
    s = 0.0

    x_0 = x[ii]
    y_0 = y[ii]

    i = ii
    while(s < s_thres)

        s = sqrt((x[i] - x_0)^2 + (y[i] - y_0)^2)

        if forward_dir
            i = i + 1
        else
            i = i - 1
        end
    end
    i
end

#=
Positive sign = retraction
negative sign = protraction
=#
function get_theta_n(c,v::Array{Float64,1},f_t::Array{Float64,1},t_c::Array{Float64,1},t=trues(length(f_t)))

    c_sign = zeros(Float64,length(v))
    t_n = zeros(Float64,length(v))

    i=1
    while (i < length(c_sign))

        if c[i]
            j = findnext(.!c,i)

            if j == nothing
                j = length(c_sign)
            end

            rr = i:(j-1) #These are the indexes corresponding to this contact epoch
            rr = rr[t[i:(j-1)]] #we remove any frames that are not tracked

            c_sign[rr] .= sign(mean(v[rr]))

            for ii in rr
                if c_sign[ii] > 0.0 #retraction
                    t_n[ii] = t_c[ii] + pi/2
                else  #protraction
                    t_n[ii] = t_c[ii] - pi/2
                end
            end

            i = j
        else
           i = i + 1
        end
    end

    (c_sign, t_n)
end

#=

=#
function calculate_all_forces(xx,yy,p,c,aa::Array{Float64,1},
    curv::Array{Float64,1},tracked=trues(length(c)); i_p_loc=50.0,moment_of_inertia_p=1.0, whisker_length=350.0)

    F_ax=zeros(Float64,length(c))
    F_lat=zeros(Float64,length(c))
    M=zeros(Float64,length(c))
    F_t=zeros(Float64,length(c))
    theta_c = zeros(Float64,length(c))
    F_calc=falses(length(c))
    r_0 = zeros(Float64,length(c))

    for i=1:length(c)
        if ((c[i])&(length(xx[i])>1))&(tracked[i])

            (new_x, new_y) = interpolate_whisker(xx[i],yy[i])

            #ii - index of contact
            ii=calc_p_dist(new_x,new_y,p[i,1],p[i,2])[2]

            #i_p - index of high SNR point
            i_p=culm_dist(new_x,new_y,i_p_loc)

            try
                (M[i],F_ax[i],F_lat[i],F_t[i],theta_c[i],r_0[i])=calc_force(new_x,new_y,aa[i],curv[i],ii,i_p,I_p=moment_of_inertia_p)
                F_calc[i] = true
                !isfinite(F_t[i]) && (F_t[i] = 0.0; F_calc[i] = false)
                !isfinite(F_lat[i]) && (F_lat[i] = 0.0; F_calc[i] = false)
                !isfinite(F_ax[i]) && (F_ax[i] = 0.0; F_calc[i] = false)
                !isfinite(M[i]) && (M[i] = 0.0; F_calc[i] = false)
            catch
            end
        end
    end

    (F_ax,F_lat,M,F_t,theta_c,F_calc,r_0)
end


text:
include("graph.jl")
include("zadatak71BFS.jl")

M = [0 1 1 0 0 0;
     1 0 0 1 0 1;
     1 0 0 0 1 0;
     0 1 0 0 0 0;
     0 0 1 0 0 1;
     0 1 0 0 1 0]
G = nonEmptyGraph(M)
path = BFS!(G, 1, 6)

println("Path = $path");


text:
using CSyntax

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))

# init window
width, height = fb_width, fb_height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# vertex data
points = GLfloat[ 0.0,  0.5, 0.0,
                  0.5, -0.5, 0.0,
                 -0.5, -0.5, 0.0]

colors = GLfloat[ 1.0, 0.0, 0.0,
                  0.0, 1.0, 0.0,
                  0.0, 0.0, 1.0]

# create VBO
points_vbo = GLuint(0)
@c glGenBuffers(1, &points_vbo)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

colors_vbo = GLuint(0)
@c glGenBuffers(1, &colors_vbo)
glBindBuffer(GL_ARRAY_BUFFER, colors_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(colors), colors, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, colors_vbo)
glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)

# create shader program
vert_shader = createshader(joinpath(@__DIR__, "geom.vert"), GL_VERTEX_SHADER)
geom_shader = createshader(joinpath(@__DIR__, "geom.geom"), GL_GEOMETRY_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "geom.frag"), GL_FRAGMENT_SHADER)
shader_prog = createprogram(vert_shader, geom_shader, frag_shader)

# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

let
updatefps = FPSCounter()
# rendering loop
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glDrawArrays(GL_POINTS, 0, 3)
    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
######################################################
#
# MMAR Plot
#
# Date:   2014-11-13
# Author: Hessel Tuinhof
#
# Simulate and plot a Brownian Motion in multifractal
# trading time (also called MMAR).
#
# Required packages: 	PyPlot
# Required functions:	MMAR.jl
#
######################################################

# load packages & functions
using PyPlot
require("../src/mmar/MMAR.jl")

# define colors
rouge = (0.79, 0.0, 0.1)

# configure parameters
m0 = 0.6
num_cascades = 10

# set seed
srand(314159265358979323846)

# plot MMAR (returns & log price)
figure()
mmar = MMAR(m0, num_cascades)
plot(mmar, linewidth = 0.7, color = rouge)
xlim(xmax = 2^num_cascades)
figure()
plot(cumsum(mmar), linewidth = 1, color = rouge)
xlim(xmax = 2^num_cascades)

max_stars_repo_name: hesseltuinhof/masterthesis

max_issues_repo_name: hesseltuinhof/masterthesis

max_forks_repo_name: hesseltuinhof/masterthesis


text:
using UnitCellGenerator
using UnitCellModelling
using FEPreProcessing
using Materials
const UCG = UnitCellGenerator
const UCM = UnitCellModelling

working_folder = mkpath(joinpath(homedir(), "WorkOuts", "scratch"))
println("\t\t Finding Effective properties of composite material using RUC of circular fibres \n")

# ============================================
#           UNIT CELL GENERATION
# ============================================
const ruc_bounds = ((-1.0, -1.0, 1.0, 1.0) .* 10.0)  #TODO getting RUC size
const ruc_info = RUC_data(
    bbox = UCG.BBox2D(ruc_bounds...,),
    ssd_ratio = 0.07,
    inclusion_distr = RANDOM,
    periodicity = true,
)
const circ_inc = Inclusion_data(
    volume_fraction = 0.50,
    shape = Circle,
    size_params = Dict(:RADIUS => Normal(2.0, 0.0),),
)
const caps_inc = Inclusion_data(
    volume_fraction = 0.2,
    shape = Capsule,
    size_params = Dict(:SMJRX => Normal(1.42, 0.0), :SMNRX => Normal(0.925, 0.0),),
)
const elliptical_inclusions_data = Inclusion_data(
    volume_fraction = 0.30,
    shape = Ellipse,
    size_params = Dict(:SMJRX => Normal(2.0, 0.0), :SMNRX => Normal(1.0, 0.0),),
)
const rectangular_inclusions_data = Inclusion_data(
    volume_fraction = 0.30,
    shape = Rectangle,
    size_params = Dict(:SMJRX => Normal(2.0, 0.0), :SMNRX => Normal(1.0, 0.0), :CRAD => Normal(0.2, 0.0),),
)
inclusions_data = UCG.generate(ruc_info, (elliptical_inclusions_data,),)

# ============================================
#       MODELLING and MESHING in GMESH
# ============================================
udc_3d = UCM.UDC3D(
    UCG.BBox3D(
        ruc_bounds[1], ruc_bounds[2], -0.5,
        ruc_bounds[3], ruc_bounds[4], 0.5,
        ),
    inclusions_data,)
ruc_model_data = UCM.create(
    udc_3d,
    mesh_periodicity=true,
    element_types = (:C3D6, :C3D8),
    geom_export_paths=(),
    extr_dir_num_ele = Int64[3,],
    extr_dir_cum_heights = Float64[1.0,],
    extr_dir_recombine_ele = true,
    min_ele_size_factor = 0.10,  #FIXME
    max_ele_size_factor = 0.15,
    mesh_opt_algorithm = "Netgen",
    show_mesh_stats = true,
    show_rve = false,
)

# all_ntags = ruc_model_data["mesh_data"]["all_node_tags"]
# all_ncoor = ruc_model_data["mesh_data"]["all_node_coordinates"]
# matrix_ele_conn = ruc_model_data["mesh_data"]["matrix_element_connectivity"]
# inclusions_ele_conn = ruc_model_data["mesh_data"]["inclusions_element_connectivity"]
# total_num_elements = ruc_model_data["mesh_data"]["num_3D_elements"]
#
# ============================================
#             FE Pre-Processing
# ============================================

matrix_material = Materials.aluminium_matrix
fibre_material = Materials.boron_fibre

uc_femodel = FEPreProcessing.unit_cell_FEModel(
    ruc_model_data["mesh_data"]["all_node_tags"],
    ruc_model_data["mesh_data"]["all_node_coordinates"],
    Dict(
        "Matrix Element Set " => (
            ruc_model_data["mesh_data"]["matrix_element_connectivity"],
            matrix_material,
        ),
        "Inclusions Element Set " => (
            ruc_model_data["mesh_data"]["inclusions_element_connectivity"],
            fibre_material,
        ),
    ),
)


text:
immutable Hypergeometric <: DiscreteUnivariateDistribution
    ns::Float64 # number of successes in population
    nf::Float64 # number of failures in population
    n::Float64  # sample size
    function Hypergeometric(s::Real, f::Real, n::Real)
        isinteger(s) && zero(s) <= s || error("ns must be a non-negative integer")
        isinteger(f) && zero(f) <= f || error("nf must be a non-negative integer")
        isinteger(n) && zero(n) < n < s + f ||
            error("n must be a positive integer <= (ns + nf)")
        new(float64(s), float64(f), float64(n))
    end
end

@_jl_dist_3p Hypergeometric hyper

function insupport(d::Hypergeometric, x::Real)
    isinteger(x) && zero(x) <= x <= d.n && (d.n - d.nf) <= x <= d.ns
end

mean(d::Hypergeometric) = d.n * d.ns / (d.ns + d.nf)

function var(d::Hypergeometric)
    N = d.ns + d.nf
    p = d.ns / N
    d.n * p * (1.0 - p) * (N - d.n) / (N - 1.0)
end


text:
module LabDataSources

using Mmap, Printf, Tables, DataFrames, CSV, CodecZlib, MAT, C3D, EzXML, OnlineStats,
    ProgressMeter, Crayons, Reexport

@reexport using DatasetManager

include("c3d.jl")
include("opensim.jl")
include("visual3d.jl")
include("dflow.jl")

end