text:
include("auxiliary-iterators.jl")

abstract type AbstractCTMIterable end
abstract type AbstractCTMState end
include("algorithms/nosym.jl")
include("algorithms/rotsym.jl")
include("algorithms/transpconjsymm.jl")

initializeC(A::DTensor{T}, χ) where T = DTensor(rand(T,χ, χ) .- 1)
initializeT(A::DTensor{T}, χ, s = 1) where T = DTensor(rand(T,χ, size(A,s), χ) .-1)

function initializeC(A::DASTensor{T,4}, χ) where T
    C = checked_similar_from_indices(nothing, T, (1,2), A, :Y)
    s1, s2 = sizes(C)
    s1 .= s2 .= χ
    initwithrand!(C)
    return C
 end

function initializeT(A::DASTensor{S,4}, χ, s=1) where S
    T = checked_similar_from_indices(nothing, S, (1,2), (s,), (1,3,2), A, A, :N, :Y)
    s1, s2, s3 = sizes(T)
    s1 .= s3 .= χ
    initwithrand!(T)
    return T
 end

isconverged(state, tol, minit) = state.n_it[] >= minit && !isempty(state.diffs) &&  state.diffs[end] < tol

function ctm(A::AbstractTensor{S,4}, Asz::AbstractTensor{S,4}, χ::Int;
                                    Cinit::Union{Nothing, AbstractTensor{S,2}} = nothing,
                                    Tinit::Union{Nothing, AbstractTensor{S,3}} = nothing,
                                    kwargs...
                                    ) where S
    if isrotsym(A)
        iter = rotsymctmiterable(A, χ, Cinit, Tinit)
    elseif istcsym(A)
        iter = transconjctmiterable(A, χ, Cinit, Tinit)
    else
        iter = ctmiterable(A, χ, Cinit, Tinit)
    end
    obs(state) = real(magnetisation(state, A, Asz))

    ctm_kernel(iter; obs = "mag" => obs, kwargs...)
end

"""
    ctm_kernel(iter::AbstractCTMIterable [, state ; tol=1e-13, maxit=10^5, period=100, verbose=true, obs=nothing)
iterates `iter` according to the options given. If a `state` is provided,
iteration starts with that state, otherwise `iter` is used as the start.
Keywords:
    - `tol`: if `tol` > 0, it is used to decide when the algorithm has converged;
       if the last element of the field `diffs` of the iterator is smaller than `tol`,
       the algorithm stops
    - `maxit`: maximum number of iterations
    - `period`: If `verbose==true`, information will be printed everye `period` iterations
    - `verbose`: if true, every `period` iterations the number of steps, the time taken since the
      start of the calculation, the last value of the `diffs` field and possibly an `obs` will be printed.
    - `obs`: a `Pair` of a string and a function on the state which will print the value of that function
      if `verbose == true`
"""
function ctm_kernel(iter::AbstractCTMIterable,
                    state::Union{Nothing, AbstractCTMState} = nothing;
                    tol::Real = 1e-13,
                    maxit::Integer = 5000,
                    minit::Integer = 1,
                    period::Integer = 100,
                    verbose::Bool = true,
                    obs = nothing)
    stop(state) = isconverged(state,tol,minit)
    function disp(state)
        ns = lpad(state[2].n_it[],7)
        ts = lpad(@sprintf("%.3E", state[1]/1e9), 13)
        ds = lpad(@sprintf("%.3E", state[2].diffs[end]),13)
        print(ns, " |", ts, " |", ds, " |")

        if !isnothing(obs)
            x = @sprintf("%.3E", obs[2](state[2]))
            print(lpad(x, 13)," |")
        end
        println()
    end

    #initialize
    isnothing(state) && (state = iterate(iter)[1])
    iter = rest(iter, state)

    iter = halt(iter, stop)
    iter = take(iter, maxit)

    if verbose
        print(lpad("n",6),"  |")
        print(lpad("t(s)", 12),"  |")
        print(lpad("diffs", 12),"  |")
        !isnothing(obs) && print(lpad(obs[1], 12),"  |")
        println()

        iter = sample(iter, period)
        iter = stopwatch(iter)
        iter = tee(iter, disp)
        (_, state) = loop(iter)
    else
        state = loop(iter)
    end

    return  state
end


text:
module IndirectArrays

export IndirectArray

"""
    IndirectArray(index, values)

creates an array `A` where the values are looked up in the value table,
`values`, using the `index`.  Concretely, `A[i...] =
values[index[i...]]`.

`values` can be an `AbstractVector` (useful when the values in `index` are
contiguous) or an `AbstractDict` (useful when they are not).
When there are not very many such distinct values, "frozen" LittleDicts from
[OrderedCollections](https://github.com/JuliaCollections/OrderedCollections.jl)
are recommended for performance.
"""
struct IndirectArray{T,N,I,A<:AbstractArray{I,N},V<:Union{AbstractVector{T},AbstractDict{I,T}}} <: AbstractArray{T,N}
    index::A
    values::V

    @inline function IndirectArray{T,N,I,A,V}(index, values) where {T,N,I,A,V}
        # The typical logic for testing bounds and then using
        # @inbounds will not check whether index is inbounds for
        # values. So we had better check this on construction.
        if isa(values, AbstractVector)
            I <: Integer || error("with a vector `values`, the index must have integer eltype")
            @boundscheck checkbounds(values, index)
        else
            @boundscheck begin
                uindex = unique(index)
                all(idx -> haskey(values, idx), index) || Base.throw_boundserror(index, values)
            end
        end
        new{T,N,I,A,V}(index, values)
    end
end
const IndirectArrayVec{T,N,I,A,V<:AbstractVector} = IndirectArray{T,N,I,A,V}

Base.@propagate_inbounds IndirectArray(index::AbstractArray{<:Integer,N}, values::AbstractVector{T}) where {T,N} =
    IndirectArray{T,N,eltype(index),typeof(index),typeof(values)}(index, values)
Base.@propagate_inbounds IndirectArray(index::AbstractArray{I,N}, values::AbstractDict{I,T}) where {I,T,N} =
    IndirectArray{T,N,I,typeof(index),typeof(values)}(index, values)

function IndirectArray{T}(A::AbstractArray) where {T}
    values = unique(A)
    index = convert(Array{T}, indexin(A, values))
    return IndirectArray(index, values)
end
IndirectArray(A::AbstractArray) = IndirectArray{UInt8}(A)

Base.size(A::IndirectArray) = size(A.index)
Base.axes(A::IndirectArray) = axes(A.index)
Base.IndexStyle(::Type{IndirectArray{T,N,I,A,V}}) where {T,N,I,A,V} = IndexStyle(A)

function Base.copy(A::IndirectArray)
    @inbounds ret = IndirectArray(copy(A.index), copy(A.values))
    return ret
end

@inline function Base.getindex(A::IndirectArray, i::Int)
    @boundscheck checkbounds(A.index, i)
    @inbounds idx = A.index[i]
    @inbounds ret = A.values[idx]
    ret
end

@inline function Base.getindex(A::IndirectArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A.index, I...)
    @inbounds idx = A.index[I...]
    @inbounds ret = A.values[idx]
    ret
end

@inline function Base.setindex!(A::IndirectArrayVec, x, i::Int)
    @boundscheck checkbounds(A.index, i)
    idx = findfirst(isequal(x), A.values)
    if idx === nothing
        push!(A.values, x)
        A.index[i] = length(A.values)
    else
        A.index[i] = idx
    end
    return A
end

@inline function Base.push!(A::IndirectArrayVec{T,1} where T, x)
    idx = findfirst(isequal(x), A.values)
    if idx === nothing
        push!(A.values, x)
        push!(A.index, length(A.values))
    else
        push!(A.index, idx)
    end
    return A
end

function Base.append!(A::IndirectArrayVec{T,1}, B::IndirectArray{T,1}) where T
    if A.values == B.values
        append!(A.index, B.index)
    else # pretty inefficient but let's get something going
        for b in B
            push!(A, b)
        end
    end
    return A
end

function Base.append!(A::IndirectArrayVec{<:Any,1}, B::AbstractVector)
    for b in B
        push!(A, b)
    end
    return A
end

end # module

ext: jl


text:
mutable struct SRC0 <: SRCType
end

function (obj::SRC0)(
              src::AbstractVector,
              xyz::AbstractVector,
              params::ParamType{3},
              t)

  error("SRC0 should never be called")
end

#
# Given PRIMITIVE variable and its 1st order and
# 2nd order derivatives, calculate the source terms
#
function calcMmsSource(params::ParamType{3},
                       q::AbstractArray{Tsrc, 1},
                       q_x::AbstractArray{Tsrc, 2},
                       q_xx::AbstractArray{Tsrc, 3},
                       src::AbstractArray{Tsrc, 1}) where Tsrc
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma - 1

  p   = q[1]*q[5]/gamma
  E   = q[5]/(gamma*gamma_1) + 0.5*(q[2]*q[2] + q[3]*q[3] + q[4]*q[4])

  p_x   = 1.0/gamma*(q_x[1,1]*q[5] + q[1]*q_x[1,5])
  p_y   = 1.0/gamma*(q_x[2,1]*q[5] + q[1]*q_x[2,5])
  p_z   = 1.0/gamma*(q_x[3,1]*q[5] + q[1]*q_x[3,5])
  E_x   = q_x[1,5]/(gamma*gamma_1) + (q[2]*q_x[1,2] + q[3]*q_x[1,3] + q[4]*q_x[1,4])
  E_y   = q_x[2,5]/(gamma*gamma_1) + (q[2]*q_x[2,2] + q[3]*q_x[2,3] + q[4]*q_x[2,4])
  E_z   = q_x[3,5]/(gamma*gamma_1) + (q[2]*q_x[3,2] + q[3]*q_x[3,3] + q[4]*q_x[3,4])

  #
  # contribution from inviscid terms
  #
  src[1]  = q_x[1,1]*q[2] + q[1]*q_x[1,2]
  src[1] += q_x[2,1]*q[3] + q[1]*q_x[2,3]
  src[1] += q_x[3,1]*q[4] + q[1]*q_x[3,4]

  src[2]  = q_x[1,1]*q[2]*q[2] + q[1]*q_x[1,2]*q[2] + q[1]*q[2]*q_x[1,2] + p_x
  src[2] += q_x[2,1]*q[2]*q[3] + q[1]*q_x[2,2]*q[3] + q[1]*q[2]*q_x[2,3]
  src[2] += q_x[3,1]*q[2]*q[4] + q[1]*q_x[3,2]*q[4] + q[1]*q[2]*q_x[3,4]

  src[3]  = q_x[1,1]*q[3]*q[2] + q[1]*q_x[1,3]*q[2] + q[1]*q[3]*q_x[1,2]
  src[3] += q_x[2,1]*q[3]*q[3] + q[1]*q_x[2,3]*q[3] + q[1]*q[3]*q_x[2,3] + p_y
  src[3] += q_x[3,1]*q[3]*q[4] + q[1]*q_x[3,3]*q[4] + q[1]*q[3]*q_x[3,4]

  src[4]  = q_x[1,1]*q[4]*q[2] + q[1]*q_x[1,4]*q[2] + q[1]*q[4]*q_x[1,2]
  src[4] += q_x[2,1]*q[4]*q[3] + q[1]*q_x[2,4]*q[3] + q[1]*q[4]*q_x[2,3]
  src[4] += q_x[3,1]*q[4]*q[4] + q[1]*q_x[3,4]*q[4] + q[1]*q[4]*q_x[3,4] + p_z

  src[5]  = q_x[1,1]*E*q[2] + q[1]*E_x*q[2] + q[1]*E*q_x[1,2] + p_x*q[2] + p*q_x[1,2]
  src[5] += q_x[2,1]*E*q[3] + q[1]*E_y*q[3] + q[1]*E*q_x[2,3] + p_y*q[3] + p*q_x[2,3]
  src[5] += q_x[3,1]*E*q[4] + q[1]*E_z*q[4] + q[1]*E*q_x[3,4] + p_z*q[4] + p*q_x[3,4]

  if !params.isViscous
    return nothing
  end

  muK = Array{typeof(q[5])}(2)
  getMuK(q[5], muK)
  rmu = muK[1]
  rK = muK[2]
  two3rd = 2.0/3.0

  #
  # in case someone forgets the lower diagonal of q_xx
  #
  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[3,1,2] = q_xx[1,3,2]
  q_xx[3,2,2] = q_xx[2,3,2]
  q_xx[2,1,3] = q_xx[1,2,3]
  q_xx[3,1,3] = q_xx[1,3,3]
  q_xx[3,2,3] = q_xx[2,3,3]
  q_xx[2,1,4] = q_xx[1,2,4]
  q_xx[3,1,4] = q_xx[1,3,4]
  q_xx[3,2,4] = q_xx[2,3,4]

  # tau
  txx = rmu * two3rd * (2*q_x[1,2] - q_x[2,3] - q_x[3,4])
  tyy = rmu * two3rd * (2*q_x[2,3] - q_x[1,2] - q_x[3,4])
  tzz = rmu * two3rd * (2*q_x[3,4] - q_x[1,2] - q_x[2,3])
  txy = rmu * (q_x[2,2] + q_x[1,3])
  txz = rmu * (q_x[3,2] + q_x[1,4])
  tyz = rmu * (q_x[2,4] + q_x[3,3])
  tyx = txy
  tzx = txz
  tzy = tyz

  # gradient of tau
  txx_x = rmu*two3rd*(2*q_xx[1,1,2] - q_xx[2,1,3] - q_xx[3,1,4])
  txx_y = rmu*two3rd*(2*q_xx[1,2,2] - q_xx[2,2,3] - q_xx[3,2,4])
  txx_z = rmu*two3rd*(2*q_xx[1,3,2] - q_xx[2,3,3] - q_xx[3,3,4])

  tyy_x = rmu*two3rd*(2*q_xx[2,1,3] - q_xx[1,1,2] - q_xx[3,1,4])
  tyy_y = rmu*two3rd*(2*q_xx[2,2,3] - q_xx[1,2,2] - q_xx[3,2,4])
  tyy_z = rmu*two3rd*(2*q_xx[2,3,3] - q_xx[1,3,2] - q_xx[3,3,4])

  tzz_x = rmu*two3rd*(2*q_xx[3,1,4] - q_xx[1,1,2] - q_xx[2,1,3])
  tzz_y = rmu*two3rd*(2*q_xx[3,2,4] - q_xx[1,2,2] - q_xx[2,2,3])
  tzz_z = rmu*two3rd*(2*q_xx[3,3,4] - q_xx[1,3,2] - q_xx[2,3,3])

  txy_x = rmu*(q_xx[2,1,2] + q_xx[1,1,3])
  txy_y = rmu*(q_xx[2,2,2] + q_xx[1,2,3])
  txy_z = rmu*(q_xx[2,3,2] + q_xx[1,3,3])
  txz_x = rmu*(q_xx[3,1,2] + q_xx[1,1,4])
  txz_y = rmu*(q_xx[3,2,2] + q_xx[1,2,4])
  txz_z = rmu*(q_xx[3,3,2] + q_xx[1,3,4])
  tyz_x = rmu*(q_xx[2,1,4] + q_xx[3,1,3])
  tyz_y = rmu*(q_xx[2,2,4] + q_xx[3,2,3])
  tyz_z = rmu*(q_xx[2,3,4] + q_xx[3,3,3])
  tyx_x = txy_x
  tyx_y = txy_y
  tyx_z = txy_z
  tzx_x = txz_x
  tzx_y = txz_y
  tzx_z = txz_z
  tzy_x = tyz_x
  tzy_y = tyz_y
  tzy_z = tyz_z

  # these coefficients are from nondimensionalization.
  Pr = 0.72
  c1 = params.euler_params.Ma/params.Re
  c2 = c1/(Pr*gamma_1)
  src[2] -= c1*(txx_x + txy_y + txz_z)
  src[3] -= c1*(tyx_x + tyy_y + tyz_z)
  src[4] -= c1*(tzx_x + tzy_y + tzz_z)
  src[5] -= c1*(txx_x*q[2] + txx*q_x[1,2] + txy_x*q[3] + txy*q_x[1,3] + txz_x*q[4] + txz*q_x[1,4])
  src[5] -= c1*(tyx_y*q[2] + tyx*q_x[2,2] + tyy_y*q[3] + tyy*q_x[2,3] + tyz_y*q[4] + tyz*q_x[2,4])
  src[5] -= c1*(tzx_z*q[2] + tzx*q_x[3,2] + tzy_z*q[3] + tzy*q_x[3,3] + tzz_z*q[4] + tzz*q_x[3,4])
  src[5] -= c2*rK*(q_xx[1,1,5] + q_xx[2,2,5] + q_xx[3,3,5])

  return nothing
end

function calcMmsSource(params::ParamType{2},
                       q::AbstractArray{Tsrc, 1},
                       q_x::AbstractArray{Tsrc, 2},
                       q_xx::AbstractArray{Tsrc, 3},
                       src::AbstractArray{Tsrc, 1}) where Tsrc
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma - 1

  p   = q[1]*q[4]/gamma
  E   = q[4]/(gamma*gamma_1) + 0.5*(q[2]*q[2] + q[3]*q[3])

  p_x   = 1.0/gamma*(q_x[1,1]*q[4] + q[1]*q_x[1,4])
  p_y   = 1.0/gamma*(q_x[2,1]*q[4] + q[1]*q_x[2,4])
  E_x   = q_x[1,4]/(gamma*gamma_1) + (q[2]*q_x[1,2] + q[3]*q_x[1,3])
  E_y   = q_x[2,4]/(gamma*gamma_1) + (q[2]*q_x[2,2] + q[3]*q_x[2,3])

  #
  # contribution from inviscid terms
  #
  src[1]  = q_x[1,1]*q[2] + q[1]*q_x[1,2]
  src[1] += q_x[2,1]*q[3] + q[1]*q_x[2,3]

  src[2]  = q_x[1,1]*q[2]*q[2] + q[1]*q_x[1,2]*q[2] + q[1]*q[2]*q_x[1,2] + p_x
  src[2] += q_x[2,1]*q[2]*q[3] + q[1]*q_x[2,2]*q[3] + q[1]*q[2]*q_x[2,3]

  src[3]  = q_x[1,1]*q[3]*q[2] + q[1]*q_x[1,3]*q[2] + q[1]*q[3]*q_x[1,2]
  src[3] += q_x[2,1]*q[3]*q[3] + q[1]*q_x[2,3]*q[3] + q[1]*q[3]*q_x[2,3] + p_y

  src[4]  = q_x[1,1]*E*q[2] + q[1]*E_x*q[2] + q[1]*E*q_x[1,2] + p_x*q[2] + p*q_x[1,2]
  src[4] += q_x[2,1]*E*q[3] + q[1]*E_y*q[3] + q[1]*E*q_x[2,3] + p_y*q[3] + p*q_x[2,3]

  if !params.isViscous
    return nothing
  end

  #
  # in case someone forgets the lower diagonal of q_xx
  #
  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[2,1,3] = q_xx[1,2,3]

  muK = Array{typeof(q[4])}(2)
  getMuK(q[4], muK)
  rmu = muK[1]
  rK = muK[2]
  two3rd = 2.0/3.0

  txx = rmu * two3rd * (2*q_x[1,2] - q_x[2,3])
  tyy = rmu * two3rd * (2*q_x[2,3] - q_x[1,2])
  txy = rmu * (q_x[2,2] + q_x[1,3])
  tyx = txy

  txx_x = rmu*two3rd*(2*q_xx[1,1,2] - q_xx[2,1,3])
  txx_y = rmu*two3rd*(2*q_xx[1,2,2] - q_xx[2,2,3])
  tyy_x = rmu*two3rd*(2*q_xx[2,1,3] - q_xx[1,1,2])
  tyy_y = rmu*two3rd*(2*q_xx[2,2,3] - q_xx[1,2,2])

  txy_x = rmu*(q_xx[2,1,2] + q_xx[1,1,3])
  txy_y = rmu*(q_xx[2,2,2] + q_xx[1,2,3])
  tyx_x = txy_x
  tyx_y = txy_y

  Pr = 0.72
  c1 = params.euler_params.Ma/params.Re
  c2 = c1/(Pr*gamma_1)
  src[2] -= c1*(txx_x + txy_y)
  src[3] -= c1*(tyx_x + tyy_y)
  src[4] -= c1*(txx_x*q[2] + txx*q_x[1,2] + txy_x*q[3] + txy*q_x[1,3])
  src[4] -= c1*(tyx_y*q[2] + tyx*q_x[2,2] + tyy_y*q[3] + tyy*q_x[2,3])
  src[4] -= c2*rK*(q_xx[1,1,4] + q_xx[2,2,4])

  return nothing
end
#
#TODO: As far as we know `q`, `q_x`, `q_xx`, we can move
# the rest into a function, which is less error prone, and
# will save a lot of space
#
mutable struct SRCPolynomial <: SRCType
end
function (obj::SRCPolynomial)(
              src::AbstractVector,
              xyz::AbstractVector,
              params::ParamType{3},
              t)

  Tdim = 3
  sigma = 0.01
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma_1
  aoa = params.euler_params.aoa
  beta = params.sideslip_angle
  q = zeros(typeof(src[1]), Tdim+2)
  q_x = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef = zeros(typeof(src[1]), Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma * cos(beta) * cos(aoa)
  qRef[3] = params.euler_params.Ma * sin(beta) * -1
  qRef[4] = params.euler_params.Ma * cos(beta) * sin(aoa)
  qRef[5] = 1.0

  x = xyz[1]
  y = xyz[2]
  z = xyz[3]

  q[1] = (x-x*x)*(y-y*y)* (z - z*z)
  q[2] = (x-x*x)*(y-y*y)* (z - z*z)
  q[3] = (x-x*x)*(y-y*y)* (z - z*z)
  q[4] = (x-x*x)*(y-y*y)* (z - z*z)
  q[5] = (x-x*x)*(y-y*y)* (z - z*z)

  q_x[1,1] = (1.0 - 2.0*x) * (y - y*y) * (z - z*z)
  q_x[1,2] = (1.0 - 2.0*x) * (y - y*y) * (z - z*z)
  q_x[1,3] = (1.0 - 2.0*x) * (y - y*y) * (z - z*z)
  q_x[1,4] = (1.0 - 2.0*x) * (y - y*y) * (z - z*z)
  q_x[1,5] = (1.0 - 2.0*x) * (y - y*y) * (z - z*z)

  q_x[2,1] = (x - x*x) * (1.0 - 2.0*y) * (z - z*z)
  q_x[2,2] = (x - x*x) * (1.0 - 2.0*y) * (z - z*z)
  q_x[2,3] = (x - x*x) * (1.0 - 2.0*y) * (z - z*z)
  q_x[2,4] = (x - x*x) * (1.0 - 2.0*y) * (z - z*z)
  q_x[2,5] = (x - x*x) * (1.0 - 2.0*y) * (z - z*z)

  q_x[3,1] = (x - x*x) * (y - y*y) * (1.0 - 2.0*z)
  q_x[3,2] = (x - x*x) * (y - y*y) * (1.0 - 2.0*z)
  q_x[3,3] = (x - x*x) * (y - y*y) * (1.0 - 2.0*z)
  q_x[3,4] = (x - x*x) * (y - y*y) * (1.0 - 2.0*z)
  q_x[3,5] = (x - x*x) * (y - y*y) * (1.0 - 2.0*z)

  q[1] = (sigma*q[1] + 1.0)*qRef[1]
  q[2] = (sigma*q[2] + 1.0)*qRef[2]
  q[3] = (sigma*q[3] + 1.0)*qRef[3]
  q[4] = (sigma*q[4] + 1.0)*qRef[4]
  q[5] = (sigma*q[5] + 1.0)*qRef[5]

  q_x[1,1] = q_x[1,1]*qRef[1]*sigma
  q_x[2,1] = q_x[2,1]*qRef[1]*sigma
  q_x[3,1] = q_x[3,1]*qRef[1]*sigma

  q_x[1,2] = q_x[1,2]*qRef[2]*sigma
  q_x[2,2] = q_x[2,2]*qRef[2]*sigma
  q_x[3,2] = q_x[3,2]*qRef[2]*sigma

  q_x[1,3] = q_x[1,3]*qRef[3]*sigma
  q_x[2,3] = q_x[2,3]*qRef[3]*sigma
  q_x[3,3] = q_x[3,3]*qRef[3]*sigma

  q_x[1,4] = q_x[1,4]*qRef[4]*sigma
  q_x[2,4] = q_x[2,4]*qRef[4]*sigma
  q_x[3,4] = q_x[3,4]*qRef[4]*sigma

  q_x[1,5] = q_x[1,5]*qRef[5]*sigma
  q_x[2,5] = q_x[2,5]*qRef[5]*sigma
  q_x[3,5] = q_x[3,5]*qRef[5]*sigma

  if !params.isViscous
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end
  #
  # contribution from viscous terms
  #
  q_xx[1,1,2] = -2.0 * (y - y*y) * (z - z*z)
  q_xx[1,1,3] = -2.0 * (y - y*y) * (z - z*z)
  q_xx[1,1,4] = -2.0 * (y - y*y) * (z - z*z)
  q_xx[1,1,5] = -2.0 * (y - y*y) * (z - z*z)

  q_xx[2,2,2] = -2.0 * (x - x*x) * (z - z*z)
  q_xx[2,2,3] = -2.0 * (x - x*x) * (z - z*z)
  q_xx[2,2,4] = -2.0 * (x - x*x) * (z - z*z)
  q_xx[2,2,5] = -2.0 * (x - x*x) * (z - z*z)

  q_xx[3,3,2] = -2.0 * (x - x*x) * (y - y*y)
  q_xx[3,3,3] = -2.0 * (x - x*x) * (y - y*y)
  q_xx[3,3,4] = -2.0 * (x - x*x) * (y - y*y)
  q_xx[3,3,5] = -2.0 * (x - x*x) * (y - y*y)

  q_xx[1,2,2] = (1.0 - 2.0*x) * (1.0 - 2.0*y) * (z - z*z)
  q_xx[1,2,3] = (1.0 - 2.0*x) * (1.0 - 2.0*y) * (z - z*z)
  q_xx[1,2,4] = (1.0 - 2.0*x) * (1.0 - 2.0*y) * (z - z*z)

  q_xx[1,3,2] = (1.0 - 2.0*x) * (1.0 - 2.0*z) * (y - y*y)
  q_xx[1,3,3] = (1.0 - 2.0*x) * (1.0 - 2.0*z) * (y - y*y)
  q_xx[1,3,4] = (1.0 - 2.0*x) * (1.0 - 2.0*z) * (y - y*y)

  q_xx[2,3,2] = (1.0 - 2.0*y) * (1.0 - 2.0*z) * (x - x*x)
  q_xx[2,3,3] = (1.0 - 2.0*y) * (1.0 - 2.0*z) * (x - x*x)
  q_xx[2,3,4] = (1.0 - 2.0*y) * (1.0 - 2.0*z) * (x - x*x)

  q_xx[1,1,2] *= sigma*qRef[2]
  q_xx[2,2,2] *= sigma*qRef[2]
  q_xx[3,3,2] *= sigma*qRef[2]
  q_xx[1,2,2] *= sigma*qRef[2]
  q_xx[1,3,2] *= sigma*qRef[2]
  q_xx[2,3,2] *= sigma*qRef[2]
  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[3,1,2] = q_xx[1,3,2]
  q_xx[3,2,2] = q_xx[2,3,2]

  q_xx[1,1,3] *= sigma*qRef[3]
  q_xx[2,2,3] *= sigma*qRef[3]
  q_xx[3,3,3] *= sigma*qRef[3]
  q_xx[1,2,3] *= sigma*qRef[3]
  q_xx[1,3,3] *= sigma*qRef[3]
  q_xx[2,3,3] *= sigma*qRef[3]
  q_xx[2,1,3] = q_xx[1,2,3]
  q_xx[3,1,3] = q_xx[1,3,3]
  q_xx[3,2,3] = q_xx[2,3,3]

  q_xx[1,1,4] *= sigma*qRef[4]
  q_xx[2,2,4] *= sigma*qRef[4]
  q_xx[3,3,4] *= sigma*qRef[4]
  q_xx[1,2,4] *= sigma*qRef[4]
  q_xx[1,3,4] *= sigma*qRef[4]
  q_xx[2,3,4] *= sigma*qRef[4]
  q_xx[2,1,4] = q_xx[1,2,4]
  q_xx[3,1,4] = q_xx[1,3,4]
  q_xx[3,2,4] = q_xx[2,3,4]

  q_xx[1,1,5] *= sigma*qRef[5]
  q_xx[2,2,5] *= sigma*qRef[5]
  q_xx[3,3,5] *= sigma*qRef[5]

  calcMmsSource(params, q, q_x, q_xx, src)

  return nothing
end

function (obj::SRCPolynomial)(
              src::AbstractVector,
              coords::AbstractVector,
              params::ParamType{2},
              t)
  Tdim = 2
  sigma = 0.01
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma_1
  aoa = params.euler_params.aoa
  q = zeros(typeof(src[1]), Tdim+2)
  q_x = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef = zeros(typeof(src[1]), Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0
  x = coords[1]
  y = coords[2]

  q[1] = (x-x*x)*(y-y*y)
  q[2] = (x-x*x)*(y-y*y)
  q[3] = (x-x*x)*(y-y*y)
  q[4] = (x-x*x)*(y-y*y)
  q_x[1,1] = (1.0 - 2.0*x) * (y - y*y)
  q_x[1,2] = (1.0 - 2.0*x) * (y - y*y)
  q_x[1,3] = (1.0 - 2.0*x) * (y - y*y)
  q_x[1,4] = (1.0 - 2.0*x) * (y - y*y)
  q_x[2,1] = (x - x*x) * (1.0 - 2.0*y)
  q_x[2,2] = (x - x*x) * (1.0 - 2.0*y)
  q_x[2,3] = (x - x*x) * (1.0 - 2.0*y)
  q_x[2,4] = (x - x*x) * (1.0 - 2.0*y)

  q[1] = (sigma*q[1] + 1.0)*qRef[1]
  q[2] = (sigma*q[2] + 1.0)*qRef[2]
  q[3] = (sigma*q[3] + 1.0)*qRef[3]
  q[4] = (sigma*q[4] + 1.0)*qRef[4]
  q_x[1,1] = q_x[1,1]*qRef[1]*sigma
  q_x[2,1] = q_x[2,1]*qRef[1]*sigma
  q_x[1,2] = q_x[1,2]*qRef[2]*sigma
  q_x[2,2] = q_x[2,2]*qRef[2]*sigma
  q_x[1,3] = q_x[1,3]*qRef[3]*sigma
  q_x[2,3] = q_x[2,3]*qRef[3]*sigma
  q_x[1,4] = q_x[1,4]*qRef[4]*sigma
  q_x[2,4] = q_x[2,4]*qRef[4]*sigma

  if !params.isViscous
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end
  #
  # contribution from viscous terms
  #
  q_xx[1,1,2] = -2.0 * (y - y*y)
  q_xx[2,2,2] = -2.0 * (x - x*x)
  q_xx[1,2,2] = (1.0 - 2.0*x) * (1.0 - 2.0*y)

  q_xx[1,1,3] = -2.0 * (y - y*y)
  q_xx[2,2,3] = -2.0 * (x - x*x)
  q_xx[1,2,3] = (1.0 - 2.0*x) * (1.0 - 2.0*y)

  q_xx[1,1,4] = -2.0 * (y - y*y)
  q_xx[2,2,4] = -2.0 * (x - x*x)

  q_xx[1,1,2] *= sigma*qRef[2]
  q_xx[1,2,2] *= sigma*qRef[2]
  q_xx[2,2,2] *= sigma*qRef[2]

  q_xx[1,1,3] *= sigma*qRef[3]
  q_xx[1,2,3] *= sigma*qRef[3]
  q_xx[2,2,3] *= sigma*qRef[3]

  q_xx[1,1,4] *= sigma*qRef[4]
  q_xx[2,2,4] *= sigma*qRef[4]

  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[2,1,3] = q_xx[1,2,3]

  calcMmsSource(params, q, q_x, q_xx, src)

  return nothing
end

mutable struct SRCChannel <: SRCType
end
function (obj::SRCChannel)(
              src::AbstractVector,
              xyz::AbstractVector,
              params::ParamType{3},
              t)
  Tdim = 3
  sigma = 0.01
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma_1
  aoa = params.euler_params.aoa
  beta = params.sideslip_angle
  qRef = zeros(typeof(src[1]), Tdim+2)
  q = zeros(typeof(src[1]), Tdim+2)
  q_x = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma * cos(beta) * cos(aoa)
  qRef[3] = params.euler_params.Ma * sin(beta) * -1
  qRef[4] = params.euler_params.Ma * cos(beta) * sin(aoa)
  qRef[5] = 1.0

  x = xyz[1]
  y = xyz[2]
  z = xyz[3]

  q[1] = (1 + sigma *x*y*z) * qRef[1]
  ux = sin(pi*x) + 1
  uy = sin(pi*y) + 1
  uz = sin(pi*z) + 1
  q[2]  = (1 + sigma * ux * uy * uz) * qRef[2]
  vx = sin(pi*x) + 1
  vy = sin(pi*y) + 1
  vz = sin(pi*z) + 1
  q[3]  = (1 + sigma * ux * uy * uz) * qRef[3]
  wx = sin(pi*x) + 1
  wy = sin(pi*y) + 1
  wz = sin(pi*z) + 1
  q[4]  = wx * wy * wz * qRef[4]
  q[4]  = (1 + sigma * ux * uy * uz) * qRef[4]
  q[5]  = qRef[5]

  q_x[1,1] = qRef[1] * y*z*sigma
  q_x[2,1] = qRef[1] * x*z*sigma
  q_x[3,1] = qRef[1] * x*y*sigma

  q_x[1,2] = qRef[2] * pi * cos(pi*x) * uy * uz * sigma
  q_x[2,2] = qRef[2] * pi * cos(pi*y) * ux * uz * sigma
  q_x[3,2] = qRef[2] * pi * cos(pi*z) * ux * uy * sigma

  q_x[1,3] = qRef[3] * pi * cos(pi*x) * uy * uz * sigma
  q_x[2,3] = qRef[3] * pi * cos(pi*y) * ux * uz * sigma
  q_x[3,3] = qRef[3] * pi * cos(pi*z) * ux * uy * sigma

  q_x[1,4] = qRef[4] * pi * cos(pi*x) * uy * uz * sigma
  q_x[2,4] = qRef[4] * pi * cos(pi*y) * ux * uz * sigma
  q_x[3,4] = qRef[4] * pi * cos(pi*z) * ux * uy * sigma

  q_x[1,5] = 0
  q_x[2,5] = 0
  q_x[3,5] = 0

  if !params.isViscous
    q[2] += 0.2 * qRef[2]
    q[3] += 0.2 * qRef[3]
    q[4] += 0.2 * qRef[4]
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end
  #
  # contribution from viscous terms
  #
  q_xx[1,1,2] = -sigma*qRef[2]*pi*pi * sin(pi*x) * uy * uz
  q_xx[2,2,2] = -sigma*qRef[2]*pi*pi * sin(pi*y) * ux * uz
  q_xx[3,3,2] = -sigma*qRef[2]*pi*pi * sin(pi*z) * ux * uy
  q_xx[1,2,2] =  sigma*qRef[2]*pi*pi * cos(pi*x) * cos(pi*y) * uz
  q_xx[1,3,2] =  sigma*qRef[2]*pi*pi * cos(pi*x) * cos(pi*z) * uy
  q_xx[2,3,2] =  sigma*qRef[2]*pi*pi * cos(pi*y) * cos(pi*z) * ux
  q_xx[2,1,2] =  q_xx[1,2,2]
  q_xx[3,1,2] =  q_xx[1,3,2]
  q_xx[3,2,2] =  q_xx[2,3,2]

  q_xx[1,1,3] = -sigma*qRef[3]*pi*pi * sin(pi*x) * uy * uz
  q_xx[2,2,3] = -sigma*qRef[3]*pi*pi * sin(pi*y) * ux * uz
  q_xx[3,3,3] = -sigma*qRef[3]*pi*pi * sin(pi*z) * ux * uy
  q_xx[1,2,3] =  sigma*qRef[3]*pi*pi * cos(pi*x) * cos(pi*y) * uz
  q_xx[1,3,3] =  sigma*qRef[3]*pi*pi * cos(pi*x) * cos(pi*z) * uy
  q_xx[2,3,3] =  sigma*qRef[3]*pi*pi * cos(pi*y) * cos(pi*z) * ux
  q_xx[2,1,3] =  q_xx[1,2,3]
  q_xx[3,1,3] =  q_xx[1,3,3]
  q_xx[3,2,3] =  q_xx[2,3,3]

  q_xx[1,1,4] = -sigma*qRef[4]*pi*pi * sin(pi*x) * uy * uz
  q_xx[2,2,4] = -sigma*qRef[4]*pi*pi * sin(pi*y) * ux * uz
  q_xx[3,3,4] = -sigma*qRef[4]*pi*pi * sin(pi*z) * ux * uy
  q_xx[1,2,4] =  sigma*qRef[4]*pi*pi * cos(pi*x) * cos(pi*y) * uz
  q_xx[1,3,4] =  sigma*qRef[4]*pi*pi * cos(pi*x) * cos(pi*z) * uy
  q_xx[2,3,4] =  sigma*qRef[4]*pi*pi * cos(pi*y) * cos(pi*z) * ux
  q_xx[2,1,4] =  q_xx[1,2,4]
  q_xx[3,1,4] =  q_xx[1,3,4]
  q_xx[3,2,4] =  q_xx[2,3,4]

  q_xx[1,1,5] = 0
  q_xx[2,2,5] = 0
  q_xx[3,3,5] = 0

  calcMmsSource(params, q, q_x, q_xx, src)
  return nothing
end

function (obj::SRCChannel)(
              src::AbstractVector,
              coords::AbstractVector,
              params::ParamType{2},
              t)
  Tdim = 2
  pi = 3.14159265358979323846264338
  sigma = 0.1
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma_1
  aoa = params.euler_params.aoa
  q = zeros(typeof(src[1]), Tdim+2)
  q_x = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef = zeros(typeof(src[1]), Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0
  x = coords[1]
  y = coords[2]

  # Exact solution in form of primitive variables

  q[1] = qRef[1] * (exp(sin(0.5*pi*(x+y))) *sigma +  1.0)
  ux  = (exp(x) * sin(pi*x) * sigma + 1) * qRef[2]
  uy   = exp(y) * sin(pi*y)
  q[2] = ux * uy
  vx  = (exp(x) * sin(pi*x) * sigma + 1) * qRef[3]
  vy   = exp(y) * sin(pi*y)
  q[3] = vx * vy
  q[4] = (1 + sigma* exp(0.1*x+0.1*y)) * qRef[4]

  #
  # contribution from inviscid terms
  #
  q_x[1,1] = qRef[1] * exp(sin(0.5*pi*(x+y))) * 0.5*pi*cos(0.5*pi*(x+y)) * sigma
  q_x[2,1] = qRef[1] * exp(sin(0.5*pi*(x+y))) * 0.5*pi*cos(0.5*pi*(x+y)) * sigma
  ux_x = qRef[2] * exp(x) * (sin(pi*x) + pi* cos(pi*x)) * sigma
  uy_y = exp(y) * (sin(pi*y) + pi * cos(pi*y))
  q_x[1,2] = ux_x * uy
  q_x[2,2] = ux * uy_y
  vx_x = qRef[3] * exp(x) * (sin(pi*x) + pi* cos(pi*x)) * sigma
  vy_y = exp(y) * (sin(pi*y) + pi * cos(pi*y))
  q_x[1,3] = vx_x * vy
  q_x[2,3] = vx * vy_y
  q_x[1,4] = exp(0.1*x+0.1*y) * qRef[4] * sigma * 0.1
  q_x[2,4] = exp(0.1*x+0.1*y) * qRef[4] * sigma * 0.1
  # q_x[1,4] = 0.0
  # q_x[2,4] = 0.0

  if !params.isViscous
    q[2] += qRef[2] * 0.2
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end

  ux_xx  = exp(x) * (sin(pi*x) + pi*cos(pi*x))
  ux_xx += exp(x) * (cos(pi*x) - pi*sin(pi*x)) * pi
  ux_xx *= sigma * qRef[2]
  uy_yy  = exp(y) * (sin(pi*y) + pi*cos(pi*y))
  uy_yy += exp(y) * (cos(pi*y) - pi*sin(pi*y)) * pi
  q_xx[1,1,2] = ux_xx * uy
  q_xx[1,2,2] = ux_x * uy_y
  q_xx[2,2,2] = ux * uy_yy
  q_xx[2,1,2] = q_xx[1,2,2]

  vx_xx  = exp(x) * (sin(pi*x) + pi*cos(pi*x))
  vx_xx += exp(x) * (cos(pi*x) - pi*sin(pi*x)) * pi
  vx_xx *= sigma * qRef[3]
  vy_yy  = exp(y) * (sin(pi*y) + pi*cos(pi*y))
  vy_yy += exp(y) * (cos(pi*y) - pi*sin(pi*y)) * pi
  q_xx[1,1,3] = vx_xx * vy
  q_xx[1,2,3] = vx_x * vy_y
  q_xx[2,2,3] = vx * vy_yy
  q_xx[2,1,3] = q_xx[1,2,3]

  q_xx[1,1,4] = exp(0.1*x+0.1*y) * qRef[4] * sigma * 0.01
  q_xx[2,2,4] = exp(0.1*x+0.1*y) * qRef[4] * sigma * 0.01
  # q_xx[1,1,4] = 0.0
  # q_xx[2,2,4] = 0.0

  calcMmsSource(params, q, q_x, q_xx, src)

  return nothing
end

#######################
#  ___________
# |           |
# |    ---    |
# |   |   |   |
# |    ---    |
# |           |
#  -----------
#
#######################
# inner block = nonslip wall
# outer block = freestream

# This MMS is not working very well.
# (I don't remember clearly, but it's not converging).
mutable struct SRCDoubleSquare <: SRCType
end
function (obj::SRCDoubleSquare)(
              src::AbstractVector,
              coords::AbstractVector,
              params::ParamType{2},
              t)
  Tdim = 2
  pi = 3.14159265358979323846264338
  sigma = 0.01
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma - 1
  aoa = params.euler_params.aoa
  q = zeros(typeof(src[1]), Tdim+2)
  q_x = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef = zeros(typeof(src[1]), Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0
  x = coords[1]
  y = coords[2]
  #
  # Exact solution in form of primitive variables
  #
  si = [0.5, 1.5]
  a = [-2.375, 16.875, -45.0, 55.0, -30.0, 6.0]
  gx    = 0.0
  gy    = 0.0
  gx_x  = 0.0
  gy_y  = 0.0
  gx_xx = 0.0
  gy_yy = 0.0

  if x >= si[1] && x < si[2]
    gx = a[1] + a[2]*x + a[3]*x*x + a[4]*x^3 + a[5]*x^4 + a[6]*x^5
    gx_x = a[2] + 2*a[3]*x + 3*a[4]*x*x + 4*a[5]*x*x*x + 5*a[6]*x^4
    gx_xx = 2*a[3] + 6*a[4]*x + 12*a[5]*x*x + 20*a[6]*x*x*x
  elseif x >= si[2]
    gx = 1.0
  elseif x <= -si[1] && x > -si[2]
    gx = a[1] - a[2]*x + a[3]*x*x - a[4]*x^3 + a[5]*x^4 - a[6]*x^5
    gx_x = -a[2] + 2*a[3]*x - 3*a[4]*x*x + 4*a[5]*x*x*x - 5*a[6]*x^4
    gx_xx = 2*a[3] - 6*a[4]*x + 12*a[5]*x*x - 20*a[6]*x*x*x
  elseif x <= -si[2]
    gx = 1.0
  end
  if y >= si[1] && y < si[2]
    gy = a[1] + a[2]*y + a[3]*y*y + a[4]*y^3 + a[5]*y^4 + a[6]*y^5
    gy_y = a[2] + 2*a[3]*y + 3*a[4]*y*y + 4*a[5]*y*y*y + 5*a[6]*y^4
    gy_yy = 2*a[3] + 6*a[4]*y + 12*a[5]*y*y + 20*a[6]*y*y*y
  elseif y >= si[2]
    gy = 1.0
  elseif y <= -si[1] && y > -si[2]
    gy = a[1] - a[2]*y + a[3]*y*y - a[4]*y^3 + a[5]*y^4 - a[6]*y^5
    gy_y = -a[2] + 2*a[3]*y - 3*a[4]*y*y + 4*a[5]*y*y*y - 5*a[6]*y^4
    gy_yy = 2*a[3] - 6*a[4]*y + 12*a[5]*y*y - 20*a[6]*y*y*y
  elseif y <= -si[2]
    gy = 1.0
  end

  q[1] = qRef[1]
  q[2] = qRef[2] * (gx + gy - gx*gy)
  q[3] = qRef[3] * (gx + gy - gx*gy)
  q[4] = qRef[4]

  q_x[1,1] = 0.0
  q_x[2,1] = 0.0
  q_x[1,2] = qRef[2] * (gx_x - gx_x * gy)
  q_x[2,2] = qRef[2] * (gy_y - gx * gy_y)
  q_x[1,3] = qRef[3] * (gx_x - gx_x * gy)
  q_x[2,3] = qRef[3] * (gy_y - gx * gy_y)
  q_x[1,4] = 0.0
  q_x[2,4] = 0.0

  if !params.isViscous
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end

  q_xx[1,1,2] = qRef[2] * (gx_xx - gx_xx * gy)
  q_xx[1,2,2] = qRef[2] * (-gx_x * gy_y)
  q_xx[2,2,2] = qRef[2] * (gy_yy - gx * gy_yy)
  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[1,1,3] = qRef[3] * (gx_xx - gx_xx * gy)
  q_xx[1,2,3] = qRef[3] * (-gx_x * gy_y)
  q_xx[2,2,3] = qRef[3] * (gy_yy - gx * gy_yy)
  q_xx[2,1,3] = q_xx[1,2,3]
  q_xx[1,1,4] = 0.0
  q_xx[2,2,4] = 0.0

  calcMmsSource(params, q, q_x, q_xx, src)

  return nothing
end

mutable struct SRCTrigonometric <: SRCType
end
function (obj::SRCTrigonometric)(
              src::AbstractVector,
              xyz::AbstractVector,
              params::ParamType{3},
              t)
  Tdim = 3
  sigma = 0.0001
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma - 1
  aoa = params.euler_params.aoa
  beta = params.sideslip_angle
  qRef = zeros(typeof(src[1]), Tdim+2)
  q    = zeros(typeof(src[1]), Tdim+2)
  q_x  = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma * cos(beta) * cos(aoa)
  qRef[3] = params.euler_params.Ma * sin(beta) * -1
  qRef[4] = params.euler_params.Ma * cos(beta) * sin(aoa)
  qRef[5] = 1.0

  xyz1 = 1 * pi * xyz
  xyz2 = 2 * pi * xyz
  xyz4 = 4 * pi * xyz
  sin_val_1 = sin(xyz1)
  cos_val_1 = cos(xyz1)
  sin_val_2 = sin(xyz2)
  cos_val_2 = cos(xyz2)
  sin_val_4 = sin(xyz4)
  cos_val_4 = cos(xyz4)

  q[1] = sin_val_2[1] * sin_val_2[2] * sin_val_2[3]
  q[2] = sin_val_4[1] * sin_val_4[2] * sin_val_4[3]
  q[3] = sin_val_2[1] * sin_val_2[2] * sin_val_2[3]
  q[4] = sin_val_1[1] * sin_val_1[2] * sin_val_1[3]
  q[5] = (1.0 - cos_val_4[1]) * (1.0 - cos_val_4[2]) * (1.0 - cos_val_4[3])

  q_x[1,1] = 2*pi * cos_val_2[1] * sin_val_2[2] * sin_val_2[3]
  q_x[1,2] = 4*pi * cos_val_4[1] * sin_val_4[2] * sin_val_4[3]
  q_x[1,3] = 2*pi * cos_val_2[1] * sin_val_2[2] * sin_val_2[3]
  q_x[1,4] =   pi * cos_val_1[1] * sin_val_1[2] * sin_val_1[3]
  q_x[1,5] = 4*pi * sin_val_4[1] * (1.0 - cos_val_4[2]) * (1.0 - cos_val_4[3])

  q_x[2,1] = 2*pi * cos_val_2[2] * sin_val_2[1] * sin_val_2[3]
  q_x[2,2] = 4*pi * cos_val_4[2] * sin_val_4[1] * sin_val_4[3]
  q_x[2,3] = 2*pi * cos_val_2[2] * sin_val_2[1] * sin_val_2[3]
  q_x[2,4] =   pi * cos_val_1[2] * sin_val_1[1] * sin_val_1[3]
  q_x[2,5] = 4*pi * sin_val_4[2] * (1.0 - cos_val_4[1]) * (1.0 - cos_val_4[3])

  q_x[3,1] = 2*pi * cos_val_2[3] * sin_val_2[1] * sin_val_2[2]
  q_x[3,2] = 4*pi * cos_val_4[3] * sin_val_4[1] * sin_val_4[2]
  q_x[3,3] = 2*pi * cos_val_2[3] * sin_val_2[1] * sin_val_2[2]
  q_x[3,4] =    pi* cos_val_1[3] * sin_val_1[1] * sin_val_1[2]
  q_x[3,5] = 4*pi * sin_val_4[3] * (1.0 - cos_val_4[1]) * (1.0 - cos_val_4[2])

  q[1] = (sigma*q[1] + 1.0)*qRef[1]
  q[2] = (sigma*q[2] + 1.0)*qRef[2]
  q[3] = (sigma*q[3] + 1.0)*qRef[3]
  q[4] = (sigma*q[4] + 1.0)*qRef[4]
  q[5] = (sigma*q[5] + 1.0)*qRef[5]

  q_x[1,1] = q_x[1,1]*qRef[1]*sigma
  q_x[2,1] = q_x[2,1]*qRef[1]*sigma
  q_x[3,1] = q_x[3,1]*qRef[1]*sigma

  q_x[1,2] = q_x[1,2]*qRef[2]*sigma
  q_x[2,2] = q_x[2,2]*qRef[2]*sigma
  q_x[3,2] = q_x[3,2]*qRef[2]*sigma

  q_x[1,3] = q_x[1,3]*qRef[3]*sigma
  q_x[2,3] = q_x[2,3]*qRef[3]*sigma
  q_x[3,3] = q_x[3,3]*qRef[3]*sigma

  q_x[1,4] = q_x[1,4]*qRef[4]*sigma
  q_x[2,4] = q_x[2,4]*qRef[4]*sigma
  q_x[3,4] = q_x[3,4]*qRef[4]*sigma

  q_x[1,5] = q_x[1,5]*qRef[5]*sigma
  q_x[2,5] = q_x[2,5]*qRef[5]*sigma
  q_x[3,5] = q_x[3,5]*qRef[5]*sigma

  if !params.isViscous
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end

  q_xx[1,1,2] = -16*pi*pi * sin_val_4[1] * sin_val_4[2] * sin_val_4[3]
  q_xx[1,1,3] = - 4*pi*pi * sin_val_2[1] * sin_val_2[2] * sin_val_2[3]
  q_xx[1,1,4] =    -pi*pi * sin_val_1[1] * sin_val_1[2] * sin_val_1[3]
  q_xx[1,1,5] =  16*pi*pi * cos_val_4[1] * (1.0 - cos_val_4[2]) * (1.0 - cos_val_4[3])

  q_xx[2,2,2] = -16*pi*pi * sin_val_4[2] * sin_val_4[1] * sin_val_4[3]
  q_xx[2,2,3] =  -4*pi*pi * sin_val_2[2] * sin_val_2[1] * sin_val_2[3]
  q_xx[2,2,4] =    -pi*pi * sin_val_1[2] * sin_val_1[1] * sin_val_1[3]
  q_xx[2,2,5] =  16*pi*pi * cos_val_4[2] * (1.0 - cos_val_4[1]) * (1.0 - cos_val_4[3])

  q_xx[3,3,2] = -16*pi*pi * sin_val_4[3] * sin_val_4[1] * sin_val_4[2]
  q_xx[3,3,3] =  -4*pi*pi * sin_val_2[3] * sin_val_2[1] * sin_val_2[2]
  q_xx[3,3,4] =    -pi*pi * sin_val_1[3] * sin_val_1[1] * sin_val_1[2]
  q_xx[3,3,5] =  16*pi*pi * cos_val_4[3] * (1.0 - cos_val_4[1]) * (1.0 - cos_val_4[2])

  q_xx[1,2,2] = 16*pi*pi * cos_val_4[1] * cos_val_4[2] * sin_val_4[3]
  q_xx[1,2,3] =  4*pi*pi * cos_val_2[1] * cos_val_2[2] * sin_val_2[3]
  q_xx[1,2,4] =    pi*pi * cos_val_1[1] * cos_val_1[2] * sin_val_1[3]

  q_xx[1,3,2] = 16*pi*pi * cos_val_4[1] * cos_val_4[3] * sin_val_4[2]
  q_xx[1,3,3] =  4*pi*pi * cos_val_2[1] * cos_val_2[3] * sin_val_2[2]
  q_xx[1,3,4] =    pi*pi * cos_val_1[1] * cos_val_1[3] * sin_val_1[2]

  q_xx[2,3,2] = 16*pi*pi * cos_val_4[2] * cos_val_4[3] * sin_val_4[1]
  q_xx[2,3,3] =  4*pi*pi * cos_val_2[2] * cos_val_2[3] * sin_val_2[1]
  q_xx[2,3,4] =    pi*pi * cos_val_1[2] * cos_val_1[3] * sin_val_1[1]

  q_xx[1,1,2] *= sigma*qRef[2]
  q_xx[2,2,2] *= sigma*qRef[2]
  q_xx[3,3,2] *= sigma*qRef[2]
  q_xx[1,2,2] *= sigma*qRef[2]
  q_xx[1,3,2] *= sigma*qRef[2]
  q_xx[2,3,2] *= sigma*qRef[2]
  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[3,1,2] = q_xx[1,3,2]
  q_xx[3,2,2] = q_xx[2,3,2]

  q_xx[1,1,3] *= sigma*qRef[3]
  q_xx[2,2,3] *= sigma*qRef[3]
  q_xx[3,3,3] *= sigma*qRef[3]
  q_xx[1,2,3] *= sigma*qRef[3]
  q_xx[1,3,3] *= sigma*qRef[3]
  q_xx[2,3,3] *= sigma*qRef[3]
  q_xx[2,1,3] = q_xx[1,2,3]
  q_xx[3,1,3] = q_xx[1,3,3]
  q_xx[3,2,3] = q_xx[2,3,3]

  q_xx[1,1,4] *= sigma*qRef[4]
  q_xx[2,2,4] *= sigma*qRef[4]
  q_xx[3,3,4] *= sigma*qRef[4]
  q_xx[1,2,4] *= sigma*qRef[4]
  q_xx[1,3,4] *= sigma*qRef[4]
  q_xx[2,3,4] *= sigma*qRef[4]
  q_xx[2,1,4] = q_xx[1,2,4]
  q_xx[3,1,4] = q_xx[1,3,4]
  q_xx[3,2,4] = q_xx[2,3,4]

  q_xx[1,1,5] *= sigma*qRef[5]
  q_xx[2,2,5] *= sigma*qRef[5]
  q_xx[3,3,5] *= sigma*qRef[5]

  calcMmsSource(params, q, q_x, q_xx, src)

  return nothing
end

function (obj::SRCTrigonometric)(
              src::AbstractVector,
              coords::AbstractVector,
              params::ParamType{2},
              t)
  Tdim = 2
  sigma = 0.01
  gamma = params.euler_params.gamma
  gamma_1 = gamma - 1.0
  aoa = params.euler_params.aoa
  q    = zeros(typeof(src[1]), Tdim+2)
  q_x  = zeros(typeof(src[1]), Tdim, Tdim+2)
  q_xx = zeros(typeof(src[1]), Tdim, Tdim, Tdim+2)
  qRef = zeros(typeof(src[1]), Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0
  x = coords[1]
  y = coords[2]
  x2 = 2*x*pi
  y2 = 2*y*pi
  x3 = 3*x*pi
  y3 = 3*y*pi
  x4 = 4*x*pi
  y4 = 4*y*pi
  sx2 = sin(x2)
  sy2 = sin(y2)
  sx3 = sin(x3)
  sy3 = sin(y3)
  sx4 = sin(x4)
  sy4 = sin(y4)
  cx2 = cos(x2)
  cx3 = cos(x3)
  cx4 = cos(x4)
  cy2 = cos(y2)
  cy3 = cos(y3)
  cy4 = cos(y4)
  #
  # Exact solution in form of primitive variables
  #
  q[1] = sx2 * sy2
  q[2] = sx4 * sy4
  q[3] = sx3 * sy3
  q[4] = (1.0 - cx2) * (1.0 - cy2)
  q_x[1,1] = 2*pi * cx2 * sy2
  q_x[2,1] = 2*pi * sx2 * cy2
  q_x[1,2] = 4*pi* cx4 * sy4
  q_x[2,2] = 4*pi* sx4 * cy4
  q_x[1,3] = 3*pi* cx3 * sy3
  q_x[2,3] = 3*pi* cy3 * sx3
  q_x[1,4] = 2*pi* sx2 * (1.0 - cy2)
  q_x[2,4] = 2*pi* sy2 * (1.0 - cx2)

  q[1] = (sigma*q[1] + 1.0)*qRef[1]
  q[2] = (sigma*q[2] + 1.0)*qRef[2]
  q[3] = (sigma*q[3] + 1.0)*qRef[3]
  q[4] = (sigma*q[4] + 1.0)*qRef[4]

  q_x[1,1] = q_x[1,1]*qRef[1]*sigma
  q_x[2,1] = q_x[2,1]*qRef[1]*sigma
  q_x[1,2] = q_x[1,2]*qRef[2]*sigma
  q_x[2,2] = q_x[2,2]*qRef[2]*sigma
  q_x[1,3] = q_x[1,3]*qRef[3]*sigma
  q_x[2,3] = q_x[2,3]*qRef[3]*sigma
  q_x[1,4] = q_x[1,4]*qRef[4]*sigma
  q_x[2,4] = q_x[2,4]*qRef[4]*sigma

  if !params.isViscous
    calcMmsSource(params, q, q_x, q_xx, src)
    return nothing
  end

  q_xx[1,1,2] = -16*pi*pi * sx4 * sy4
  q_xx[2,2,2] = -16*pi*pi * sx4 * sy4
  q_xx[1,2,2] =  16*pi*pi * cx4 * cy4
  q_xx[1,1,3] = -9*pi*pi * sx3 * sy3
  q_xx[2,2,3] = -9*pi*pi * sx3 * sy3
  q_xx[1,2,3] =  9*pi*pi * cx3 * cy3
  q_xx[1,1,4] =  4*pi*pi * cx2 * (1.0 - cy2)
  q_xx[2,2,4] =  4*pi*pi * cy2 * (1.0 - cx2)
  q_xx[1,1,2] = q_xx[1,1,2] * qRef[2] * sigma
  q_xx[1,2,2] = q_xx[1,2,2] * qRef[2] * sigma
  q_xx[2,2,2] = q_xx[2,2,2] * qRef[2] * sigma
  q_xx[2,1,2] = q_xx[1,2,2]
  q_xx[1,1,3] = q_xx[1,1,3] * qRef[3] * sigma
  q_xx[1,2,3] = q_xx[1,2,3] * qRef[3] * sigma
  q_xx[2,2,3] = q_xx[2,2,3] * qRef[3] * sigma
  q_xx[2,1,3] = q_xx[1,2,3]
  q_xx[1,1,4] = q_xx[1,1,4] * qRef[4] * sigma
  q_xx[2,2,4] = q_xx[2,2,4] * qRef[4] * sigma

  calcMmsSource(params, q, q_x, q_xx, src)
  return nothing
end

global const SRCDict = Dict{String, SRCType}(
  "SRC0" => SRC0(),
  "SRCTrigonometric" => SRCTrigonometric(),
  "SRCDoubleSquare" => SRCDoubleSquare(),
  "SRCPolynomial" => SRCPolynomial(),
  "SRCChannel" => SRCChannel()
)

"""
### NavierStokesMod.getSRCFunctors

  This function gets the functor specified by opts["SRCname"] and stores
  it to the equation object.  Currently one 1 source functor is allowed.

"""
function getSRCFunctors(mesh::AbstractMesh, sbp::AbstractOperator,
                        eqn::NSData, opts)

  # currently we only allow 1 source functor
  eqn.src_func = SRCDict[opts["SRCname"]]
  return nothing
end


text:
using Alea, FactCheck

facts("tensorise") do

end

facts("smolyak") do

  (G, factors) = smolyak_grid(3, 2)
  on01 = (x::Vector{Float64},w::Vector{Float64}) -> ((x+1)/2,w)::Rule1d
  rules = [on01(gauss_rule(LegendrePolynomials(), i)...)::Rule1d for i=1:5]
  (x,w)=collect_rules(G, factors, rules)

  I=multiindex(3, 5)
  x = reshape(x, (1, size(x)...))
  II=vec(squeeze(prod(x.^I, 2), 2)*reshape(w, (size(w)...,1)))
  II2=squeeze(prod(1./(I+1),2),2)
  ex1 = BitArray{1}(map(≈,II,II2))
  ex2 = vec(sum(floor(I/2),2)).<2
  @fact ex1 --> ex2
end

#=
@show smolyak_grid(3, 4)

polys = [HermitePolynomials(), LegendrePolynomials(), LaguerrePolynomials(0.0)]
@show polys
rules = [Alea.gauss_rule(P, 3) for P in polys]
(I, factors) = smolyak_grid(3, 4)
@show collect_rules(I, factors, [Alea.gauss_rule(HermitePolynomials(), i) for i=1:4])

(I, factors) = smolyak_grid(2, 7)
@show collect_rules(I, factors, [Alea.gauss_rule(HermitePolynomials(), i) for i=1:7])
=#

#(tensorise_nodes([r[1] for r in rules]), tensorise_weights([r[2] for r in rules]))

#=
w = collect(([2.0,3,4],[5,6]))
w = collect(([2.0,3,4],[5.0,6],[0.1,0.01]))
@show tensorise_weights(w)
@show tensorise_weights([])
@show tensorise_nodes(w)
@show tensorise_nodes([])

using FactCheck
using Alea

polys = [HermitePolynomials(), LegendrePolynomials(), LaguerrePolynomials(0.0)]
@show polys

rules = [Alea.gauss_rule(P, 3) for P in polys]
@show [r[1] for r in rules]
@show [r[2] for r in rules]

#@show tensorise_rulesX(rules)
@show rules[1][1]
@show rules[2][1]
@show rules[3][1]

@show collect((rules[1][1], rules[2][1], rules[3][1]))
@show [Alea.gauss_rule(P, 3)[1] for P in polys]

@show tensorise_rules(rules)
=#


text:
using CPUTime
using Base.Test

function add_and_sleep()
    x = 0
    for i in 1:1e7
        x += i
    end
    sleep(1)
end

function time_difference()
    tic()
    CPUtic()
    add_and_sleep()
    return toq() - CPUtoq()
end

time_difference()  # compilation overhead should not affect results, but just in case
time_diff = time_difference()
eval(parse("@test abs($time_diff - 1.0) <= .01"))


text:
"""
    blockappend!(dest::BlockVector, sources...) -> dest

Append blocks from `sources` to `dest`.  The number of blocks in `dest` are
increased by `sum(blocklength, sources)`.

This function avoids copying the elements of the blocks in `sources` when
these blocks are compatible with `dest`.  Importantly, this means that
mutating `sources` afterwards alters the items in `dest` and it may even
break the invariance of `dest` if the length of `sources` are changed.

The blocks in `dest` must not alias with `sources` or components of them.
For example, the result of `blockappend!(x, x)` is undefined.

# Examples
```jldoctest
julia> using BlockArrays

julia> blockappend!(mortar([[1], [2, 3]]), mortar([[4, 5]]))
3-blocked 5-element BlockArray{Int64,1}:
 1
 ─
 2
 3
 ─
 4
 5
```
"""
blockappend!(dest::BlockVector, s1, s2, sources...) =
    foldl(blockappend!, (s1, s2, sources...); init = dest)

function blockappend!(dest::BlockVector{<:Any,T}, src::BlockVector{<:Any,T}) where {T}
    append!(dest.blocks, src.blocks)
    offset = last(dest.axes[1]) + 1 - src.axes[1].first
    append!(dest.axes[1].lasts, (n + offset for n in src.axes[1].lasts))
    return dest
end

function blockappend!(
    dest::BlockVector{<:Any,<:AbstractArray{T}},
    src::PseudoBlockVector{<:Any,T},
) where {T}
    if blocklength(src) == 1
        return _blockpush!(dest, src.blocks)
    else
        return blockappend_fallback!(dest, src)
    end
end

blockappend!(
    dest::BlockVector{<:Any,<:AbstractArray{T}},
    src::T,
) where {T<:AbstractVector} = _blockpush!(dest, src)

blockappend!(dest::BlockVector{<:Any,<:Any}, src::AbstractVector) =
    blockappend_fallback!(dest, src)

blockappend_fallback!(dest::BlockVector{<:Any,<:AbstractArray{T}}, src) where {T} =
    blockappend!(dest, mortar([convert(T, @view src[b]) for b in blockaxes(src, 1)]))

"""
    blockpush!(dest::BlockVector, blocks...) -> dest

Push `blocks` to the end of `dest`.

This function avoids copying the elements of the `blocks` when these blocks
are compatible with `dest`.  Importantly, this means that mutating `blocks`
afterwards alters the items in `dest` and it may even break the invariance
of `dest` if the length of `blocks` are changed.

# Examples
```jldoctest
julia> using BlockArrays

julia> blockpush!(mortar([[1], [2, 3]]), [4, 5], [6])
4-blocked 6-element BlockArray{Int64,1}:
 1
 ─
 2
 3
 ─
 4
 5
 ─
 6
```
"""
blockpush!(dest::BlockVector, blocks...) = foldl(blockpush!, blocks; init = dest)

blockpush!(dest::BlockVector{<:Any,<:AbstractArray{T}}, block::T) where {T} =
    _blockpush!(dest, block)

blockpush!(dest::BlockVector, block) = _blockpush!(dest, _newblockfor(dest, block))

_newblockfor(dest, block) =
    if Iterators.IteratorSize(block) isa Union{Base.HasShape,Base.HasLength}
        copyto!(eltype(dest.blocks)(undef, length(block)), block)
    else
        foldl(push!, block; init = eltype(dest.blocks)(undef, 0))
    end

function _blockpush!(dest, block)
    push!(dest.blocks, block)
    push!(dest.axes[1].lasts, last(dest.axes[1]) + length(block))
    return dest
end

"""
    blockpushfirst!(dest::BlockVector, blocks...) -> dest

Push `blocks` to the beginning of `dest`.  See also [`blockpush!`](@ref).

This function avoids copying the elements of the `blocks` when these blocks
are compatible with `dest`.  Importantly, this means that mutating `blocks`
afterwards alters the items in `dest` and it may even break the invariance
of `dest` if the length of `blocks` are changed.

# Examples
```jldoctest
julia> using BlockArrays

julia> blockpushfirst!(mortar([[1], [2, 3]]), [4, 5], [6])
4-blocked 6-element BlockArray{Int64,1}:
 4
 5
 ─
 6
 ─
 1
 ─
 2
 3
```
"""
blockpushfirst!(A::BlockVector, b1, b2, blocks...) =
    foldl(blockpushfirst!, reverse((b1, b2, blocks...)); init = A)

blockpushfirst!(dest::BlockVector{<:Any,<:AbstractArray{T}}, block::T) where {T} =
    _blockpushfirst!(dest, block)

blockpushfirst!(dest::BlockVector{<:Any,<:Any}, block) =
    _blockpushfirst!(dest, _newblockfor(dest, block))

function _blockpushfirst!(dest, block)
    pushfirst!(dest.blocks, block)
    dest.axes[1].lasts .+= length(block) - 1 + dest.axes[1].first
    pushfirst!(dest.axes[1].lasts, length(block))
    return dest
end

"""
    blockpop!(A::BlockVector) -> block

Pop a `block` from the end of `dest`.

# Examples
```jldoctest
julia> using BlockArrays

julia> A = mortar([[1], [2, 3]]);

julia> blockpop!(A)
2-element Array{Int64,1}:
 2
 3

julia> A
1-blocked 1-element BlockArray{Int64,1}:
 1
```
"""
function blockpop!(A::BlockVector)
    block = pop!(A.blocks)
    pop!(A.axes[1].lasts)
    return block
end

"""
    blockpopfirst!(dest::BlockVector) -> block

Pop a `block` from the beginning of `dest`.

# Examples
```jldoctest
julia> using BlockArrays

julia> A = mortar([[1], [2, 3]]);

julia> blockpopfirst!(A)
1-element Array{Int64,1}:
 1

julia> A
1-blocked 2-element BlockArray{Int64,1}:
 2
 3
```
"""
function blockpopfirst!(A::BlockVector)
    block = popfirst!(A.blocks)
    n = popfirst!(A.axes[1].lasts)
    A.axes[1].lasts .-= n
    return block
end

"""
    append!(dest::BlockVector, sources...)

Append items from `sources` to the last block of `dest`.

The blocks in `dest` must not alias with `sources` or components of them.
For example, the result of `append!(x, x)` is undefined.

# Examples
```jldoctest
julia> using BlockArrays

julia> append!(mortar([[1], [2, 3]]), mortar([[4], [5]]))
2-blocked 5-element BlockArray{Int64,1}:
 1
 ─
 2
 3
 4
 5
```
"""
Base.append!(dest::BlockVector, sources...) = foldl(append!, sources; init = dest)

Base.append!(dest::BlockVector, src) = append_itr!(dest, Base.IteratorSize(src), src)

function append_itr!(dest::BlockVector, ::Union{Base.HasShape,Base.HasLength}, src)
    block = dest.blocks[end]
    li = lastindex(block)
    resize!(block, length(block) + length(src))
    # Equivalent to `i = li; for x in src; ...; end` but (maybe) faster:
    foldl(src, init = li) do i, x
        Base.@_inline_meta
        i += 1
        @inbounds block[i] = x
        return i
    end
    da, = dest.axes
    da.lasts[end] += length(src)
    return dest
end

function append_itr!(dest::BlockVector, ::Base.SizeUnknown, src)
    block = dest.blocks[end]
    # Equivalent to `n = 0; for x in src; ...; end` but (maybe) faster:
    n = foldl(src, init = 0) do n, x
        push!(block, x)
        return n + 1
    end
    da, = dest.axes
    da.lasts[end] += n
    return dest
end

# remove empty blocks at the end
function _squash_lasts!(A::BlockVector)
    while !isempty(A.blocks) && isempty(A.blocks[end])
        pop!(A.blocks)
        pop!(A.axes[1].lasts)
    end
end

# remove empty blocks at the beginning
function _squash_firsts!(A::BlockVector)
    while !isempty(A.blocks) && isempty(A.blocks[1])
        popfirst!(A.blocks)
        popfirst!(A.axes[1].lasts)
    end
end

"""
    pop!(A::BlockVector)

Pop the last element from the last non-empty block.  Remove all empty
blocks at the end.
"""
function Base.pop!(A::BlockVector)
    isempty(A) && throw(Argument("array must be nonempty"))
    _squash_lasts!(A)
    x = pop!(A.blocks[end])
    lasts = A.axes[1].lasts
    if isempty(A.blocks[end])
        pop!(A.blocks)
        pop!(lasts)
    else
        lasts[end] -= 1
    end
    return x
end

"""
    popfirst!(A::BlockVector)

Pop the first element from the first non-empty block.  Remove all empty
blocks at the beginning.
"""
function Base.popfirst!(A::BlockVector)
    isempty(A) && throw(Argument("array must be nonempty"))
    _squash_firsts!(A)
    x = popfirst!(A.blocks[1])
    ax, = A.axes
    if isempty(A.blocks[1])
        popfirst!(A.blocks)
        popfirst!(ax.lasts)
    else
        ax.lasts[1] -= 1
    end
    return x
end

"""
    push!(dest::BlockVector, items...)

Push items to the end of the last block.
"""
Base.push!(dest::BlockVector, items...) = append!(dest, items)

"""
    pushfirst!(A::BlockVector, items...)

Push items to the beginning of the first block.
"""
function Base.pushfirst!(A::BlockVector, items...)
    pushfirst!(A.blocks[1], items...)
    A.axes[1].lasts .+= length(items)
    return A
end


text:
@testset "Errors" begin
  @testset "Learning" begin
    x = rand(1:2, 1000)
    y = rand(1:2, 1000)
    X = rand(2, 1000)
    d = georef(DataFrame(x=x, y=y), PointSet(X))
    p = LearningProblem(d, d, ClassificationTask(:x, :y))
    s = PointwiseLearn(DummyModel())

    # dummy classifier → 0.5 misclassification rate
    for m in [CrossValidation(10), BlockCrossValidation(0.1),
              LeaveBallOut(0.1), DensityRatioValidation(10)]
      e = error(s, p, m)
      @test isapprox(e[:y], 0.5, atol=0.06)
    end
  end
end


text:
# Mandatory
export polyhedron, loadpolyhedron!
export hrepiscomputed, hrep, vrepiscomputed, vrep
export volume, volume_simplex, unscaled_volume_simplex, surface

"""
    polyhedron(rep::Representation{T})

Creates a polyhedron from the representation `rep` using the default library included in the Polyhedra package.
"""
polyhedron(rep::Representation{T}) where {T} = polyhedron(rep, default_library(FullDim(rep), T))

"""
    hrepiscomputed(p::Polyhedron)

Returns whether the H-representation of this polyhedron has been computed.
"""
function hrepiscomputed end

"""
    hrep(p::Polyhedron)

Returns an H-representation for the polyhedron `p`.
"""
hrep(p::Polyhedron) = error("`hrep` not implemented for `$(eltype(p))`")

"""
    vrepiscomputed(p::Polyhedron)

Returns whether the V-representation of this polyhedron has been computed.
"""
function vrepiscomputed end

"""
    vrep(p::Polyhedron)

Returns a V-representation for the polyhedron `p`.
"""
vrep(p::Polyhedron) = error("`vrep` not implemented for `$(eltype(p))`")

"""
    volume(p::Polyhedron{T}) where {T}

Returns the `fulldim(p)`-dimensional hyper-volume of the polyhedron `p`.
Returns `Inf` or `-one(T)` if it is infinite depending on whether the type `T` has an infinite value.
"""
function volume(p::Polyhedron)
    return sum(unscaled_volume_simplex(Δ) for Δ in triangulation(p)) / factorial(fulldim(p))
end

function unscaled_volume_simplex(s)
    A = Matrix{coefficient_type(s)}(undef, fulldim(s), npoints(s) - 1)
    v = first(points(s))
    for (i, p) in enumerate(Iterators.drop(points(s), 1))
        A[:, i] = p - v
    end
    return abs(det(A))
end
function volume_simplex(s)
    return unscaled_volume_simplex(s) / factorial(fulldim(s))
end

"""
    surface(p::Polyhedron{T}) where {T}

Returns the `fulldim(p)-1`-dimensional hyper-volume of the surface of the polyhedron `p`.
Returns `Inf` or `-one(T)` if it is infinite depending on whether the type `T` has an infinite value.
"""
function surface end

loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::Type{Val{:ine}}) = error("loadpolyhedron! not implemented for .ine")
loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::Type{Val{:ext}}) = error("loadpolyhedron! not implemented for .ext") # FIXME ExtFileVRepresentation or just ExtFile

loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::Symbol) = loadpolyhedron!(p, filename, Val{extension})

function loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::AbstractString)
    s = something(findfirst(isequal(extension), ["ext", "ine"]), 0)
    if s == 0
        error("Invalid extension $extension, please give 'ext' for V-representation or 'ine' for H-representation")
    end
    loadpolyhedron!(p, filename, [:ext, :ine][s])
end


text:
using ..Atmos
using ..Atmos: thermo_state, turbulence_tensors
using ..SubgridScaleParameters: inv_Pr_turb
using ..PlanetParameters
using ..MoistThermodynamics

Base.@kwdef mutable struct AtmosCollectedDiagnostics
    zvals::Union{Nothing,Matrix} = nothing
end
const CollectedDiagnostics = AtmosCollectedDiagnostics()

function init(atmos::AtmosModel,
              mpicomm::MPI.Comm,
              dg::DGModel,
              Q::MPIStateArray,
              starttime::String)
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    CollectedDiagnostics.zvals = zeros(Nqk, nvertelem)

    if Array ∈ typeof(Q).parameters
        localvgeo = grid.vgeo
    else
        localvgeo = Array(grid.vgeo)
    end

    @visitQ nhorzelem nvertelem Nqk Nq begin
        z = localvgeo[ijk,grid.x3id,e]
        CollectedDiagnostics.zvals[k,ev] = z
    end
end

# thermodynamic variables of interest
function vars_thermo(FT)
  @vars begin
    q_liq::FT
    q_ice::FT
    q_vap::FT
    T::FT
    θ_liq_ice::FT
    θ_dry::FT
    θ_v::FT
    e_int::FT
    h_m::FT
    h_t::FT
  end
end
num_thermo(FT) = varsize(vars_thermo(FT))
thermo_vars(array) = Vars{vars_thermo(eltype(array))}(array)

function compute_thermo!(bl, state, aux, k, ijk, ev, e, thermoQ)
    e_tot = state.ρe / state.ρ
    ts = thermo_state(bl.moisture, bl.orientation, state, aux)
    e_int = internal_energy(ts)
    Phpart = PhasePartition(ts)

    th = thermo_vars(thermoQ[ijk,e])
    th.q_liq     = Phpart.liq
    th.q_ice     = Phpart.ice
    th.q_vap     = vapor_specific_humidity(ts)
    th.T         = air_temperature(ts)
    th.θ_liq_ice = liquid_ice_pottemp(ts)
    th.θ_dry     = dry_pottemp(ts)
    th.θ_v       = virtual_pottemp(ts)
    th.e_int     = e_int

    # Moist and total henthalpy
    R_m          = gas_constant_air(ts)
    th.h_m       = e_int + R_m * th.T
    th.h_t       = e_tot + R_m * th.T

    return nothing
end

# horizontal averages
function vars_horzavg(FT)
  @vars begin
    ρ::FT
    ρu::FT
    ρv::FT
    ρw::FT
    e_tot::FT
    ρq_tot::FT
    q_liq::FT
    q_vap::FT
    θ_liq_ice::FT
    θ_dry::FT
    θ_v::FT
    e_int::FT
    h_m::FT
    h_t::FT
    qt_sgs::FT
    ht_sgs::FT
  end
end
num_horzavg(FT) = varsize(vars_horzavg(FT))
horzavg_vars(array) = Vars{vars_horzavg(eltype(array))}(array)

function compute_horzsums!(atmos::AtmosModel, state, diffusive_flx, aux, k, ijk,
                           ev, e, Nqk, nvertelem, MH, localaux, thermoQ,
                           horzsums, repdvsr, LWP, currtime)
    th = thermo_vars(thermoQ[ijk,e])
    hs = horzavg_vars(horzsums[k,ev])
    hs.ρ         += MH * state.ρ
    hs.ρu        += MH * state.ρu[1]
    hs.ρv        += MH * state.ρu[2]
    hs.ρw        += MH * state.ρu[3]
    hs.e_tot     += MH * state.ρe
    hs.q_liq     += MH * th.q_liq
    hs.q_vap     += MH * th.q_vap
    hs.θ_liq_ice += MH * th.θ_liq_ice
    hs.θ_dry     += MH * th.θ_dry
    hs.θ_v       += MH * th.θ_v
    hs.e_int     += MH * th.e_int
    hs.h_m       += MH * th.h_m
    hs.h_t       += MH * th.h_t

    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        hs.ρq_tot += MH * state.moisture.ρq_tot
    end

    ν, _       = turbulence_tensors(atmos.turbulence, state, diffusive_flx, aux, currtime)
    D_t        = (ν isa Real ? ν : diag(ν)) * inv_Pr_turb

    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        d_q_tot    = (-D_t) .* diffusive_flx.moisture.∇q_tot
        hs.qt_sgs += MH * state.ρ * d_q_tot[end]
    end

    d_h_tot    = -D_t .* diffusive_flx.∇h_tot
    hs.ht_sgs += MH * state.ρ * d_h_tot[end]

    repdvsr[Nqk*(ev-1)+k] += MH

    # liquid water path
    # this condition is also going to be used to get the number of points that
    # exist on a horizontal plane provided all planes have the same number of
    # points
    # TODO adjust for possibility of non equivalent horizontal slabs
    if ev == floor(nvertelem/2) && k == floor(Nqk/2)
        # TODO: uncomment the line below after rewriting the LWP assignment below using aux.∫dz...?
        # aux = extract_aux(dg, localaux, ijk, e)
        LWP[1] += MH * (localaux[ijk,1,e] + localaux[ijk,2,e])
    end

    return nothing
end

function compute_diagnosticsums!(atmos, state, k, ijk, ev, e, MH,
                                 thermoQ, horzavgs, dsums)
    zvals = CollectedDiagnostics.zvals
    th = thermo_vars(thermoQ[ijk,e])
    ha = horzavg_vars(horzavgs[k,ev])
    ds = diagnostic_vars(dsums[k,ev])

    u     = state.ρu[1] / state.ρ
    v     = state.ρu[2] / state.ρ
    w     = state.ρu[3] / state.ρ
    ũ     = ha.ρu / ha.ρ
    ṽ     = ha.ρv / ha.ρ
    w̃     = ha.ρw / ha.ρ
    ẽ     = ha.e_tot / ha.ρ
    q̃_tot = ha.ρq_tot / ha.ρ
    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        q_tot = state.moisture.ρq_tot / state.ρ
    end

    # vertical coordinate
    ds.z        += MH * zvals[k,ev]

    # state and functions of state
    ds.u        += MH * ũ
    ds.v        += MH * ṽ
    ds.w        += MH * w̃
    ds.e_tot    += MH * ẽ
    ds.q_tot    += MH * ha.ρq_tot / ha.ρ
    ds.q_liq    += MH * ha.q_liq
    ds.thd      += MH * ha.θ_dry
    ds.thl      += MH * ha.θ_liq_ice
    ds.thv      += MH * ha.θ_v
    ds.e_int    += MH * ha.e_int
    ds.h_m      += MH * ha.h_m
    ds.h_t      += MH * ha.h_t
    ds.qt_sgs   += MH * ha.qt_sgs
    ds.ht_sgs   += MH * ha.ht_sgs

    # vertical fluxes
    ds.vert_eddy_mass_flx += MH * (w - w̃) * (state.ρ - ha.ρ)
    ds.vert_eddy_u_flx    += MH * (w - w̃) * (u - ha.ρu / ha.ρ)
    ds.vert_eddy_v_flx    += MH * (w - w̃) * (v - ha.ρv / ha.ρ)
    ds.vert_eddy_ql_flx   += MH * (w - w̃) * (th.q_liq - ha.q_liq)
    ds.vert_eddy_qv_flx   += MH * (w - w̃) * (th.q_vap - ha.q_vap)
    ds.vert_eddy_thd_flx  += MH * (w - w̃) * (th.θ_dry - ha.θ_dry)
    ds.vert_eddy_thv_flx  += MH * (w - w̃) * (th.θ_v - ha.θ_v)
    ds.vert_eddy_thl_flx  += MH * (w - w̃) * (th.θ_liq_ice - ha.θ_liq_ice)
    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        ds.vert_eddy_qt_flx   += MH * (w - w̃) * (q_tot - q̃_tot)
        ds.vert_qt_flx        += MH * w * q_tot
    end

    # variances
    ds.uvariance += MH * (u - ũ)^2
    ds.vvariance += MH * (v - ṽ)^2
    ds.wvariance += MH * (w - w̃)^2

    # skewness
    ds.wskew     += MH * (w - w̃)^3

    # turbulent kinetic energy
    ds.TKE       += MH * 0.5 * (ds.uvariance + ds.vvariance + ds.wvariance)

    return nothing
end

"""
    collect(bl, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function collect(atmos::AtmosModel, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q

    mpirank = MPI.Comm_rank(mpicomm)

    # extract grid information
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get the state, auxiliary and geo variables onto the host if needed
    if Array ∈ typeof(Q).parameters
        localQ    = Q.realdata
        localaux  = dg.auxstate.realdata
        localvgeo = grid.vgeo
        localdiff = dg.diffstate.realdata
    else
        localQ    = Array(Q.realdata)
        localaux  = Array(dg.auxstate.realdata)
        localvgeo = Array(grid.vgeo)
        localdiff = Array(dg.diffstate.realdata)
    end
    FT = eltype(localQ)

    nstate    = num_state(bl, FT)
    nauxstate = num_aux(bl, FT)
    ndiff     = num_diffusive(bl, FT)

    # vertical coordinates and thermo variables
    thermoQ = [zeros(FT, num_thermo(FT)) for _ in 1:npoints, _ in 1:nrealelem]

    # divisor for horizontal averages
    l_repdvsr = zeros(FT, Nqk * nvertelem)

    # horizontal sums and the liquid water path
    l_LWP = zeros(FT, 1)
    horzsums = [zeros(FT, num_horzavg(FT)) for _ in 1:Nqk, _ in 1:nvertelem]

    # compute thermo variables and horizontal sums in a single pass
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, localQ, ijk, e)
        aux   = extract_aux(dg, localaux, ijk, e)

        compute_thermo!(bl, state, aux, k, ijk, ev, e, thermoQ)

        diffusive_flx = extract_diffusion(dg, localdiff, ijk, e)
        MH = localvgeo[ijk,grid.MHid,e]
        compute_horzsums!(bl, state, diffusive_flx, aux, k, ijk, ev, e, Nqk,
                          nvertelem, MH, localaux, thermoQ, horzsums, l_repdvsr,
                          l_LWP, currtime)
    end

    # compute the full number of points on a slab
    repdvsr = MPI.Allreduce(l_repdvsr, +, mpicomm)

    # compute the horizontal and LWP averages
    horzavgs = [zeros(FT, num_horzavg(FT)) for _ in 1:Nqk, _ in 1:nvertelem]
    for ev in 1:nvertelem
        for k in 1:Nqk
            hsum = MPI.Allreduce(horzsums[k,ev], +, mpicomm)
            horzavgs[k,ev][1:end] = hsum ./ repdvsr[Nqk*(ev-1)+k]
        end
    end
    LWP = zero(FT)
    LWP = MPI.Reduce(l_LWP[1], +, 0, mpicomm)
    if mpirank == 0
        LWP /= repdvsr[1]
    end

    # compute the diagnostics using the previous computed values
    dsums = [zeros(FT, num_diagnostic(FT)) for _ in 1:Nqk, _ in 1:nvertelem]
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, localQ, ijk, e)
        MH = localvgeo[ijk,grid.MHid,e]
        compute_diagnosticsums!(bl, state, k, ijk, ev, e, MH,
                                thermoQ, horzavgs, dsums)
    end
    davgs = [zeros(FT, num_diagnostic(FT)) for _ in 1:Nqk, _ in 1:nvertelem]
    for ev in 1:nvertelem
        for k in 1:Nqk
            dsum = MPI.Reduce(dsums[k,ev], +, 0, mpicomm)
            if mpirank == 0
                davgs[k,ev][1:end] = dsum ./ repdvsr[Nqk*(ev-1)+k]
            end
        end
    end

    # write LWP
    if mpirank == 0
        current_time = string(currtime)
        starttime = Settings.starttime
        outdir = Settings.outdir
        jldopen(joinpath(outdir,
                "liquid_water_path-$(starttime).jld2"), "a+") do file
            file[current_time] = LWP
        end
    end

    return davgs
end # function collect

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using FrankenTuples
using Test

@testset "It's alive!" begin
    x = ftuple(1, 2, a=3, b=4)
    @test x == FrankenTuple((1, 2), (a=3, b=4))
    @test x == FrankenTuple(x)
    @test !isempty(x)
    @test length(x) == 4
    @test Tuple(x) == (1, 2)
    @test NamedTuple(x) == (a=3, b=4)
    @test sprint(show, x) == "FrankenTuple((1, 2), (a = 3, b = 4))"
    @test x == @ftuple (1, 2; a=3, b=4)
    @test x == @ftuple (1, a=3, b=4, 2)
    @test keys(x) == (1, 2, :a, :b)
    @test values(x) == (1, 2, 3, 4)
    @test eltype(x) == Int
    @test first(x) == 1
    @test Base.tail(x) == ftuple(2, a=3, b=4)

    y = FrankenTuple{Tuple{Float64,Float64},(:a,:b),Tuple{Float64,Float64}}(x)
    @test y == ftuple(1.0, 2.0, a=3.0, b=4.0)

    t = ftuple(1, 2)
    @test t == FrankenTuple((1, 2))
    @test length(t) == 2
    @test Tuple(t) == (1, 2)
    @test NamedTuple(t) == NamedTuple()
    @test t == convert(FrankenTuple{Tuple{Int,Int},(),Tuple{}}, (1, 2))
    @test sprint(show, t) == "FrankenTuple((1, 2), NamedTuple())"
    @test t == @ftuple (1, 2)
    @test keys(t) == keys(Tuple(t))
    @test values(t) == (1, 2)
    @test first(t) == 1
    @test Base.tail(t) == ftuple(2)

    nt = ftuple(a=3, b=4)
    @test nt == FrankenTuple((a=3, b=4))
    @test length(nt) == 2
    @test Tuple(nt) == ()
    @test NamedTuple(nt) == (a=3, b=4)
    @test nt == convert(FrankenTuple{Tuple{},(:a,:b),Tuple{Int,Int}}, (a=3, b=4))
    @test sprint(show, nt) == "FrankenTuple((), (a = 3, b = 4))"
    @test nt == @ftuple (a=3, b=4)
    @test keys(nt) == (:a, :b)
    @test values(nt) == (3, 4)
    @test first(nt) == 3
    @test Base.tail(nt) == ftuple(b=4)

    e = empty(x)
    @test e == FrankenTuple()
    @test isempty(e)
    @test length(e) == 0
    @test sprint(show, e) == "FrankenTuple()"
    @test e == @ftuple ()
    @test keys(e) == values(e) == ()
    @test_throws ArgumentError first(e)
    @test_throws ArgumentError Base.tail(e)

    @test eltype(ftuple(1, 2.0, a=3, b=4.0f0)) == Real
end

@testset "Indexing and iteration" begin
    x = ftuple(1, 2, a=3, b=4)
    @test x.a == x[:a] == 3
    @test x.b == x[:b] == 4
    for i = 1:4
        @test x[i] == i
    end
    for (i, t) in enumerate(x)
        @test t == i
    end
    @test collect(pairs(x)) == [1 => 1, 2 => 2, :a => 3, :b => 4]
    @test firstindex(x) == 1
    @test lastindex(x) == 4
end

@testset "☎" begin
    @test ftcall(sum, ftuple(abs2, [1 2; 3 4]; dims=2)) == reshape([5, 25], (2, 1))
    @test ftcall(string, ftuple()) == ""
    @test ftcall(+, ftuple(1, 2)) == 3
    @test ftcall((; k...)->sum(values(k.data)), ftuple(a=3.0, b=0x4)) == 7.0
end

f(x::Int; y=3) = x + y
g(; b, c, a) = a + b + c
h(x::String; a, kwargs...) = x * a

@testset "hasmethod" begin
    @test hasmethod(f, typeof(ftuple(1, y=2)))
    @test hasmethod(f, typeof(ftuple(1)))  # Agreement with using a plain Tuple
    @test !hasmethod(f, typeof(ftuple(1, a=3)))
    @test hasmethod(f, FrankenTuple{Tuple{Int},(:y,)})  # Omitting NamedTuple types
    @test hasmethod(g, typeof(ftuple(a=1, b=2, c=3)))
    @test hasmethod(g, typeof(ftuple(a=1, b=2)))
    @test !hasmethod(g, FrankenTuple{Tuple{},(:a,:b,:d)})
    @test hasmethod(g, FrankenTuple{Tuple{},(:a,:b,:c)})
    @test hasmethod(h, FrankenTuple{Tuple{String},(:a,:b,:c,:d)})
    @test !hasmethod(h, FrankenTuple{Tuple{Int},(:a,)})

    @test !hasmethod(f, FrankenTuple{Tuple{Int},(:y,)}, world=typemin(UInt))
end

@testset "map" begin
    @test FrankenTuple((11,22), (x=33,y=44)) == map(+, FrankenTuple((1,2), (x=3,y=4)), FrankenTuple((10, 20), (x=30,y=40)))
    @test_throws Exception map(+, FrankenTuple((1,2), (x=3,y=4)), FrankenTuple((10, 20), (x=30,)))
    @test_throws Exception map(+, FrankenTuple((1,2), (x=3,y=4)), FrankenTuple((10,), (x=30,y=40)))
end


text:
# IFAC benchmarks
@testset "IFAC Benchmarks" begin
    @safetestset "Binary Distillation Column" begin include( "ifac/prob9001-binaryDistillationColumn.jl" ) end
    @safetestset "Drum Boiler"                begin include( "ifac/prob9002-drumBoiler.jl" )               end
    @safetestset "B767 Flutter Condition"     begin include( "ifac/prob9006-airplaneFlutter.jl" )          end
end

# Miscellaneous benchmarks
@testset "Miscellaneous Benchmarks" begin
    @safetestset "Ball on plate system" begin include( "ballOnPlate.jl" ) end
    @safetestset "Jones-Morari system"  begin include( "jonesMorari.jl" ) end
end

# Mass-Spring system benchmarks
@testset "Mass-Spring Systems" begin
    @safetestset "Linear Mass-Spring-Damper Chain"        begin include( "massSpringSystems/linearMassSpringDamper.jl" )  end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9b8c8d1a-481e-11eb-1b85-91264e100b12
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
end

# ╔═╡ 7441759c-4853-11eb-3d63-2be1f95f59fe
using Plots

# ╔═╡ 704d34fc-4859-11eb-2d95-45c4d5246b26
begin
	using Flux
	A = TabularVApproximator(;n_state=10, init=0.0, opt=InvDecay(1.0))
	A
end

# ╔═╡ 82338e22-4805-11eb-0b54-49f734001f0c
md"""
## Introduction

In the following notebooks, we'll mainly use the [ReinforcementLearning.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl) package to demonstrate how to generate figures in the book [Reinforcement Learning: An Introduction(2nd)](http://incompleteideas.net/book/the-book-2nd.html). In case you haven't used this package before, you can visit [juliareinforcementlearning.org](https://juliareinforcementlearning.org/) for a detailed introduction. Besides, we'll also explain some basic concepts gradually when we meet them for the first time.

[ReinforcementLearning.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl) contains a collection of tools to describe and solve problems we want to handle in the reinforcement learning field. Though we are mostly interested in traditional tabular methods here, it also contains many state-of-the-art algorithms. To use it, we can simply run the following code here:
"""

# ╔═╡ bfe260f4-481e-11eb-28e8-eb3cca823162
md"""
!!! note
	As you might have noticed, it takes quite a long time to load this package for the first time (the good news is that it will be largely decreased after `Julia@v1.6`). Once loaded, things should be very fast.

"""

# ╔═╡ 147a9204-4841-11eb-04fd-4f760dff4bc8
md"""
## Tic-Tac-Toe

In chapter 1.5, a simple game of [Tic-Tac-Toe](https://en.wikipedia.org/wiki/Tic-tac-toe) is introduced to illustrate the general idea of reinforcement learning. Before looking into the details about how to implement the Monte Carlo based policy, let's take a look at how the Tic-Tac-Toe environment is defined in [ReinforcementLearning.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl).
"""

# ╔═╡ 46d5181a-4845-11eb-222c-3d8b78f89dc6
env = TicTacToeEnv()

# ╔═╡ d9a93b80-4827-11eb-39cf-c589ebcd092b
md"""
There is many important information provided above. First, the `TicTacToeEnv` is a [zero-sum](https://en.wikipedia.org/wiki/Zero-sum_game), two player environment. `Sequential` means each player takes an action alternatively. Since all the players can observe the same information (the board), it is an environment of [perfect information](https://en.wikipedia.org/wiki/Perfect_information). Note that each player only gets the reward at the end of the game (`-1`, `0`, or `1`). So we call the `RewardStyle` of this kind of environments `TerminalReward()`. At each step, only part of the actions are valid (the position of the board), so we say the `ActionStyle` of this env is of `FullActionSet()`.

All this information is provided by traits. In later chapters, we'll see how to define these traits for new customized environments. Now let's get familiar with some basic interfaces defined with the environment first.
"""

# ╔═╡ fad9dddc-4827-11eb-1b79-db762523f380
state(env) |> Text  # Here we use Text to show the string properly in Pluto

# ╔═╡ 483e774c-4849-11eb-1d48-c3ba0f63a784
md"""
By default, the state of the `TicTacToeEnv` is represented as a `String`. But some other types are also provided.

First, the number of possible valid boards is numerable. You might be interested in this [blog](http://www.occasionalenthusiast.com/tag/tic-tac-toe/) to figure out the exact number. So we can use an integer to represent the state.
"""

# ╔═╡ ef6c5598-4849-11eb-1b76-d9aff5a62abb
state(env, Observation{Int}())

# ╔═╡ f07794e2-484a-11eb-2cc8-e318bf15d469
md"""
Another way to express the state is to use a one-hot bool Array to represent the state. The first and the second dimension mean the row and column of the board. The third dimension means the role in each cell (`x`, `o` or empty). This kind of representation is usually more friendly to neural network based algorithms.
"""

# ╔═╡ b8521fe6-4837-11eb-0e53-ab179c6c3ea8
state(env, Observation{BitArray{3}}())

# ╔═╡ 9e42c688-4829-11eb-0e7a-f93d399dbf09
md"""
Each environment is a functional object which takes in an action and changes its internal state correspondingly.
"""

# ╔═╡ ff62343e-484b-11eb-2a41-29768287e2b4
begin
	env(1)
	env(2)
	env(3)
	env(4)
	env(5)
	env(6)
	env(7)
	state(env) |> Text
end

# ╔═╡ 52ac24ce-484c-11eb-25ed-bf7a4bd7f9fb
md"""
Now we can check some basic information about the `env`:
"""

# ╔═╡ 62a66722-484c-11eb-32d0-7b2695b67e6b
is_terminated(env), [reward(env, p) for p in players(env)]

# ╔═╡ 83ceeb54-484c-11eb-0499-33bafa9f569b
md"""
As you can see, the game is now terminated, and the reward of each player is $(join([reward(env, p) for p in players(env)], ", ")).
"""

# ╔═╡ c842c388-48b9-11eb-3474-b177c74c02c9
reset!(env)

# ╔═╡ cd5ddb32-48b9-11eb-0a3f-41c723469e9a
state(env) |> Text

# ╔═╡ 45d38868-484d-11eb-016a-eff04cc182a8
md"""
## RandomPolicy

Among all the policies provided in `ReinforcementLearning.jl`, the simplest one is `RandomPolicy`.
"""

# ╔═╡ 1697cfde-482d-11eb-1da9-5537e88a70c5
policy = RandomPolicy()

# ╔═╡ 177f3ee6-482b-11eb-36c3-c97e945cbd0b
[policy(env) for _ in 1:10]

# ╔═╡ 5b589dbc-482b-11eb-2395-317aca83d1f4
state(env) |> Text

# ╔═╡ 505e0fe6-482b-11eb-1df7-45ebb6c65c9b
md"""
As you can see, each policy is also a functional object which takes in a `env` and returns an action.

!!! note
	A policy will never change the internal state of an `env`. The `env` will only get updated after executing `env(action)`.
"""

# ╔═╡ 19b8290e-484e-11eb-33f3-9d2329d15a46
md"""
## Policy + Environment

Now that we have understood the two main concepts in reinforcement learning, we can roll out several simulations and collect some information we are interested in.
"""

# ╔═╡ c6045bfe-484e-11eb-2b9c-6f230a90cb03
run(policy, env, StopAfterEpisode(1))

# ╔═╡ 297a9614-484f-11eb-24a4-bd80771c1f1c
md"""
The above `run` function defined in `ReinforcementLearning.jl` is quite straight forward. The `policy` generates an action at each time step and feeds it into the `env`. The process continues until the end of an episode. Here `StopAfterEpisode(1)` is a built in stop condition. You can also see many other stop conditions in the [doc](https://juliareinforcementlearning.org/ReinforcementLearning.jl/latest/rl_core/#ReinforcementLearningCore.StopAfterEpisode).

You are encouraged to read the [source code](https://github.com/JuliaReinforcementLearning/ReinforcementLearningCore.jl/blob/master/src/core/run.jl) of this function. It's pretty simple (less than 30 lines) and easy to understand. I'll wait you here until you are finished.

If you have finished reading it, you'll notice that one important argument is missing in the above function call, the `hook`. Now we'll add the fourth argument to collect the reward of each player in every episode.
"""

# ╔═╡ 6f58815c-4852-11eb-18e8-e55b276ba228
multi_agent_policy = MultiAgentManager(
	(
		NamedPolicy(p=>RandomPolicy())
		for p in players(env)
	)...
)

# ╔═╡ 7ff2ce0a-4852-11eb-3f97-c1bc0b1bc056
multi_agent_hook = MultiAgentHook(
	(
		p => TotalRewardPerEpisode()
		for p in players(env)
	)...
)

# ╔═╡ e8e0ad64-484e-11eb-3eba-e7b323105de1
run(multi_agent_policy,	env, StopAfterEpisode(10), multi_agent_hook)

# ╔═╡ bb374d92-4852-11eb-1268-976e9e5532f4
md"""
Besides that we added a fourth argument in the `run`, another important difference is that we used a `MultiAgentManager` policy instead of a simple `RandomPolicy`. The reason is that we want different players to use a separate policy and then we can collect different information of them separately.

Now let's take a look at the total reward of each player in the above 10 episodes:
"""

# ╔═╡ 6f7bddcc-4853-11eb-2e38-edc41fcbd1ec
x, o = players(env)

# ╔═╡ 97e68c62-4853-11eb-2e19-43fc6a6a4618
begin
	plot(multi_agent_hook[x][]; label="x")
	plot!(multi_agent_hook[o][]; label="o")
end

# ╔═╡ b7ded1bc-4855-11eb-02a7-2f77d10767b3
md"""
## Tackling the Tic-Tac-Toe Problem with Monte Carlo Prediction

Actually the Monte Carlo method mentioned in the first chapter to solve the Tic-Tac-Toe problem is explained in Chapter 5.1, so it's ok if you don't fully understand it right now for the first time. Just keep reading and turn back to this chapter once you have finished chapter 5.

The intuition behind Monte Carlo prediction is that we use a table to record the estimated gain at each step. If such estimation is accurate, then we can simply look one step ahead and compare which state is of the largest gain. Then we just take the action which will lead to that state to maximize our reward.

### TabularVApproximator

In `ReinforcementLearning.jl`, the fundamental components used to estimate state values (or state-action values, or any other values we are interested in) are called *AbstractApproximator*. In this simple example, we'll use the `TabularVApproximator` to estimate the state value.
"""

# ╔═╡ 3889b408-49b4-11eb-02b6-cfd5fbb885b4
md"""
A `TabularVApproximator` uses a `Vector` underlying to estimate the state value. Here the **Tabular** means that the state is required to be an `Integer`. The code above created a `TabularVApproximator` which has a table of `10` states, `1:10`, and each state value is initialized to `0.0`. We can get the estimation by:
"""

# ╔═╡ fa9d93ac-49b4-11eb-105d-958e3031f2eb
A(1)

# ╔═╡ 0d999a3c-49b5-11eb-1b74-355c4603b71a
md"""
And we can update the estimation via:
"""

# ╔═╡ 88464c8a-49b5-11eb-3327-81e653872ecf
update!(A, 2 => A(2) - 5.)

# ╔═╡ ad7f07ca-49bb-11eb-3a96-4d649425d821
md"""
Let's spend some time on what's happening here. The `update!` method provided in `ReinforcementLearning.jl` adopts the similar interface defined in [`Flux.Optimise`](https://github.com/FluxML/Flux.jl/blob/master/src/optimise/Optimise.jl). Here the `2 => A(2) - 5.)` is something similar to `Grad`. And then we apply the optimiser (`A.optimizer`) to calculate the error we want to apply on the approximator. The default optimizer in `TabularApproximator` is `InvDecay(1.0)`. The effect is the same with calculating the mean value of all the examples.
"""

# ╔═╡ e4d58428-49bc-11eb-19aa-35759a992cac
begin
	examples = 1:10
	for x in examples
		update!(A, 1 => A(1) - x)
	end
	A(1) == (0. #= the init value is 0. =# + sum(examples)) / (1+length(examples))
end

# ╔═╡ 4c7c124a-49bd-11eb-0fa9-99fe3631af10
md"""
By using the same interfaces defined in `Flux.jl`, we get a bunch of different optimizers free. You can take a look at all the supported `optimizers` defined in [optimisers.jl](https://github.com/FluxML/Flux.jl/blob/master/src/optimise/optimisers.jl). In future notebooks, we'll see how to apply some other optimizers beyond `InvDecay`.
"""

# ╔═╡ f60a5790-49bd-11eb-16b4-df2b3a682d86
md"""
### MonteCarloLearner

In `ReinforcementLearning.jl`, many different variants of Monte Carlo learners are supported. We'll skip the implementation detail for now. But to convince you that implementing such an algorithm is quite simple and straightforward in this package, we can take a look at the code snippet and compare it with the pseudocode on the book.

```julia
function _update!(
    ::EveryVisit,
    ::TabularVApproximator,
    ::NoSampling,
    L::MonteCarloLearner,
    t::AbstractTrajectory,
)
    S, R = t[:states], t[:rewards]
    V, G, γ = L.approximator, 0.0, L.γ
    for i in length(R):-1:1
        s, r = S[i], R[i]
        G = γ * G + r
        update!(V, s => V(s) - G)
    end
end
```
"""

# ╔═╡ 266f60c2-49c0-11eb-1ad9-f1c363389e38
M = MonteCarloLearner(;
	approximator = TabularVApproximator(;n_state=5478, init=0., opt=InvDecay(1.0)),
	γ = 1.0,
	kind = FIRST_VISIT,
	sampling = NO_SAMPLING
)

# ╔═╡ 55b6eadc-49d7-11eb-2080-37fd57265468
md"""
The `MonteCarloLearner` is just a wrapper of the `TabularVApproximator` here. We can feed the `env` to it and get the the state estimation.
"""

# ╔═╡ c04e438c-49d8-11eb-3888-61e46bd24f25
M(1)

# ╔═╡ 14422efe-49dc-11eb-3c4b-81a6e239b512
M(env)

# ╔═╡ ee3ddc56-49db-11eb-2150-1dc2740b1f2f
md"""
Uh-oh! The `MethodError` above tells us that the `MonteCarloLearner` expects the input to be an `Int` but gets a `String`. The reason is that the default state style of `env` is `Observation{String}()`. So we need to set the default state style to `Observation{Int}()`.
"""

# ╔═╡ 7b979bd4-49dc-11eb-26b6-4d45d01a2301
E = DefaultStateStyleEnv{Observation{Int}()}(env)

# ╔═╡ a6176830-49dc-11eb-3e04-112e6594bfbe
state(E)

# ╔═╡ 8b87797e-49de-11eb-0f22-d79e3399bd1b
md"""
Now that we have the `MonteCarloLearner` to estimate state values, how can we use it to generate actions? A simple approach is that, we try each valid action and see what state it leads to. Then we simply select the one which results in the state of the largest state value.

!!! note
	Obviously this approach only applies to deterministic environments with a small action set.

But usually this approach is more suitable for **testing** instead of **training**. Because always selecting the best action may fall into local minima and hurt the performance. You'll see more discussions on **exploration** and **exploitation** in the next chapter. Here we use the `EpsilonGreedy` method to select action.
"""

# ╔═╡ d74aecfc-4a4e-11eb-14da-75f749a4f743
explorer = EpsilonGreedyExplorer(0.1)

# ╔═╡ fedae1de-4a4e-11eb-37ed-0d32687a9b9b
begin
	values = [1, 2, 3, 1]
	N = 1000
	actions = [explorer(values) for _ in 1:N]
	plot([sum(actions .== i)/N for i in 1:length(values)])
end

# ╔═╡ d5742096-4a4f-11eb-32da-f387bcc0ca25
md"""
The above figure shows that with $\epsilon$ set to `0.1`, we have the probability of $\epsilon$ to select an action randomly and $1-\epsilon$ to select the action of the highest value.

Below, we define a function to illustrate how to select an action given an environment.
"""

# ╔═╡ a89a0f98-4a4c-11eb-2226-1b332c23ed34
function select_action(env, V)
	A = legal_action_space(env)
    values = map(A) do a
        V(child(env, a))
    end
    A[explorer(values)]
end

# ╔═╡ dc31b262-4a50-11eb-3709-d7652d590bcd
md"""
Combining the `MonteCarloLearner` and the `select_action` function, we get our policy:
"""

# ╔═╡ 42995fe4-49e0-11eb-0a60-cfe8518ee90c
P = VBasedPolicy(;learner=M, mapping=select_action)

# ╔═╡ 752cf5ec-49e0-11eb-1440-d7ce675d8fc3
reset!(E)

# ╔═╡ 79865c8c-49e0-11eb-30b8-ddf82be2f473
P(E)

# ╔═╡ ae65a714-49e0-11eb-01d4-39bacc61d24a
md"""
That's it!

Now we have a policy and we can use it to roll our simulations.
"""

# ╔═╡ d11cc120-49e0-11eb-2588-e11f814a0548
run(P, E, StopAfterEpisode(10))

# ╔═╡ f03c55e8-49e0-11eb-2dc3-a3018c810c5e
md"""
## Training

One main question we haven't answered is, **how to train the policy?**

Well, the usage is similar to the above one, the only difference is now we wrap the `policy` in an `Agent`, which is also an `AbstractPolicy`. An `Agent` is `policy + trajectory`, or people usually call it *experience replay buffer*.
"""

# ╔═╡ d04a563a-49e1-11eb-1326-b1993ee989ec
policies = MultiAgentManager(
	(
		Agent(
			policy = NamedPolicy(
				p => VBasedPolicy(;
					learner=MonteCarloLearner(;
						approximator = TabularVApproximator(;
							n_state=length(state_space(E)),
							init=0.,
							opt=InvDecay(1.0)
						),
						γ = 1.0,
						kind = FIRST_VISIT,
						sampling = NO_SAMPLING
					),
					mapping = select_action
				)
			),
			trajectory =VectorSARTTrajectory(
					;state=Int,
					action=Union{Int, NoOp},
					reward=Int,
					terminal=Bool
			)
		)
		for p in players(E)
	)...
)

# ╔═╡ 8ed8eab4-4a51-11eb-2dfb-e5a1bf1cfac3
run(policies, E, StopAfterEpisode(100_000))

# ╔═╡ 187e6fb6-4a52-11eb-3d84-e36d0f64a0c4
md"""
The interface is almost the same with the above one. Let me explain what's happening here.

First, the `policies` is a `MultiAgentManager` and the environment is a `Sequential` one. So at each step, the agent manager forwards the `env` to its inner policies. Here each inner policy is an `Agent`. Then it will collect necessary information at each step (here we used a `VectorSARTTrajectory` to tell the agent to collect `state`, `action`, `reward`, and `is_terminated` info). Finally, after each episode, the agent sends the `trajectory` to the inner `VBasedPolicy`, and then forwards it to the `MonteCarloLearner` to update the `TabularVApproximator`. Thanks to **multiple dispatch** each step above is fully customizable.
"""

# ╔═╡ 52392820-4a54-11eb-03eb-cf339de8e3dc
md"""
## Testing

Now we have our policy trained. We can happily run several self-plays and see the result. The only necessary change is to drop out the `Agent` wrapper.
"""

# ╔═╡ d5662b30-4a54-11eb-3ec6-23f78e2437b1
test_policies = MultiAgentManager([p.policy for (x, p) in (policies.agents)]...)

# ╔═╡ ff86531e-4a55-11eb-2b13-f50d077ec7b0
md"""
Here we use a hook to collect the reward of each episode for every player:
"""

# ╔═╡ 15c37092-4a56-11eb-26a3-716552d8deba
hook = MultiAgentHook(
	(
		p => TotalRewardPerEpisode()
		for p in players(E)
	)...
)

# ╔═╡ e225bf54-4a55-11eb-0ba6-f3854fc816a2
run(test_policies, E, StopAfterEpisode(100), hook)

# ╔═╡ 33261234-4a56-11eb-28b0-f5067e4c948e
begin
	plot(hook[x][]; label="x")
	plot!(hook[o][]; label="o")
end

# ╔═╡ 696c441c-4a56-11eb-39c0-e1d1465ca4f7
md"""
As you can see, in most cases, we reach a **tie**. But why are there still some cases that don't reach a tie?

Is it because we have not trained for enough time?

*Maybe*. But a more possible reason is that we're still using the `EpsilonGreedyExplorer` in the testing mode.

I leave it as an exercise to change the `select_action` to a greedy version and confirm that our trained policy reaches a tie everytime.
"""

# ╔═╡ Cell order:
# ╟─82338e22-4805-11eb-0b54-49f734001f0c
# ╠═9b8c8d1a-481e-11eb-1b85-91264e100b12
# ╟─bfe260f4-481e-11eb-28e8-eb3cca823162
# ╟─147a9204-4841-11eb-04fd-4f760dff4bc8
# ╠═46d5181a-4845-11eb-222c-3d8b78f89dc6
# ╟─d9a93b80-4827-11eb-39cf-c589ebcd092b
# ╠═fad9dddc-4827-11eb-1b79-db762523f380
# ╟─483e774c-4849-11eb-1d48-c3ba0f63a784
# ╠═ef6c5598-4849-11eb-1b76-d9aff5a62abb
# ╟─f07794e2-484a-11eb-2cc8-e318bf15d469
# ╠═b8521fe6-4837-11eb-0e53-ab179c6c3ea8
# ╟─9e42c688-4829-11eb-0e7a-f93d399dbf09
# ╠═ff62343e-484b-11eb-2a41-29768287e2b4
# ╟─52ac24ce-484c-11eb-25ed-bf7a4bd7f9fb
# ╠═62a66722-484c-11eb-32d0-7b2695b67e6b
# ╟─83ceeb54-484c-11eb-0499-33bafa9f569b
# ╠═c842c388-48b9-11eb-3474-b177c74c02c9
# ╠═cd5ddb32-48b9-11eb-0a3f-41c723469e9a
# ╟─45d38868-484d-11eb-016a-eff04cc182a8
# ╠═1697cfde-482d-11eb-1da9-5537e88a70c5
# ╠═177f3ee6-482b-11eb-36c3-c97e945cbd0b
# ╠═5b589dbc-482b-11eb-2395-317aca83d1f4
# ╟─505e0fe6-482b-11eb-1df7-45ebb6c65c9b
# ╟─19b8290e-484e-11eb-33f3-9d2329d15a46
# ╠═c6045bfe-484e-11eb-2b9c-6f230a90cb03
# ╟─297a9614-484f-11eb-24a4-bd80771c1f1c
# ╠═6f58815c-4852-11eb-18e8-e55b276ba228
# ╠═7ff2ce0a-4852-11eb-3f97-c1bc0b1bc056
# ╠═e8e0ad64-484e-11eb-3eba-e7b323105de1
# ╟─bb374d92-4852-11eb-1268-976e9e5532f4
# ╠═6f7bddcc-4853-11eb-2e38-edc41fcbd1ec
# ╠═7441759c-4853-11eb-3d63-2be1f95f59fe
# ╠═97e68c62-4853-11eb-2e19-43fc6a6a4618
# ╟─b7ded1bc-4855-11eb-02a7-2f77d10767b3
# ╠═704d34fc-4859-11eb-2d95-45c4d5246b26
# ╟─3889b408-49b4-11eb-02b6-cfd5fbb885b4
# ╠═fa9d93ac-49b4-11eb-105d-958e3031f2eb
# ╟─0d999a3c-49b5-11eb-1b74-355c4603b71a
# ╠═88464c8a-49b5-11eb-3327-81e653872ecf
# ╟─ad7f07ca-49bb-11eb-3a96-4d649425d821
# ╠═e4d58428-49bc-11eb-19aa-35759a992cac
# ╟─4c7c124a-49bd-11eb-0fa9-99fe3631af10
# ╟─f60a5790-49bd-11eb-16b4-df2b3a682d86
# ╠═266f60c2-49c0-11eb-1ad9-f1c363389e38
# ╟─55b6eadc-49d7-11eb-2080-37fd57265468
# ╠═c04e438c-49d8-11eb-3888-61e46bd24f25
# ╠═14422efe-49dc-11eb-3c4b-81a6e239b512
# ╟─ee3ddc56-49db-11eb-2150-1dc2740b1f2f
# ╠═7b979bd4-49dc-11eb-26b6-4d45d01a2301
# ╠═a6176830-49dc-11eb-3e04-112e6594bfbe
# ╟─8b87797e-49de-11eb-0f22-d79e3399bd1b
# ╠═d74aecfc-4a4e-11eb-14da-75f749a4f743
# ╠═fedae1de-4a4e-11eb-37ed-0d32687a9b9b
# ╟─d5742096-4a4f-11eb-32da-f387bcc0ca25
# ╠═a89a0f98-4a4c-11eb-2226-1b332c23ed34
# ╟─dc31b262-4a50-11eb-3709-d7652d590bcd
# ╠═42995fe4-49e0-11eb-0a60-cfe8518ee90c
# ╠═752cf5ec-49e0-11eb-1440-d7ce675d8fc3
# ╠═79865c8c-49e0-11eb-30b8-ddf82be2f473
# ╟─ae65a714-49e0-11eb-01d4-39bacc61d24a
# ╠═d11cc120-49e0-11eb-2588-e11f814a0548
# ╟─f03c55e8-49e0-11eb-2dc3-a3018c810c5e
# ╠═d04a563a-49e1-11eb-1326-b1993ee989ec
# ╠═8ed8eab4-4a51-11eb-2dfb-e5a1bf1cfac3
# ╟─187e6fb6-4a52-11eb-3d84-e36d0f64a0c4
# ╟─52392820-4a54-11eb-03eb-cf339de8e3dc
# ╠═d5662b30-4a54-11eb-3ec6-23f78e2437b1
# ╟─ff86531e-4a55-11eb-2b13-f50d077ec7b0
# ╠═15c37092-4a56-11eb-26a3-716552d8deba
# ╠═e225bf54-4a55-11eb-0ba6-f3854fc816a2
# ╠═33261234-4a56-11eb-28b0-f5067e4c948e
# ╟─696c441c-4a56-11eb-39c0-e1d1465ca4f7