text:
import Base: -, +, *, filter, real, complex
import LinearAlgebra: norm, promote_leaf_eltypes

_basetype(φ::AbstractProperty) = Base.typename(typeof(φ)).wrapper

# TODO: rewrite all this with meta-programming

@inline +(φ1::AbstractProperty, φ2::AbstractProperty) = (_basetype(φ1))( φ1.val + φ2.val )
@inline -(φ1::AbstractProperty, φ2::AbstractProperty) = (_basetype(φ1))( φ1.val - φ2.val )
@inline -(φ1::AbstractProperty) = (_basetype(φ1))( -φ1.val )
@inline *(φ1::AbstractProperty, φ2::AbstractProperty) = (_basetype(φ1))( φ1.val * φ2.val )

@inline *(a::Union{Number, AbstractMatrix}, φ::AbstractProperty) = (_basetype(φ))(a * φ.val)
@inline *(φ::AbstractProperty, a::Union{Number, AbstractMatrix})  = (_basetype(φ))(φ.val * a)

*(a::AbstractVector{<: Number}, φ::AbstractProperty) = a .* Ref(φ)

Base.isapprox(φ1::AbstractProperty, φ2::AbstractProperty) = isapprox(φ1.val, φ2.val)

@inline norm(φ::AbstractProperty) = norm(φ.val)
@inline Base.length(φ::AbstractProperty) = length(φ.val)
@inline Base.size(φ::AbstractProperty) = size(φ.val)
@inline Base.zero(φ::AbstractProperty) = (_basetype(φ))(zero(φ.val))
@inline Base.zero(::Type{T}) where {T <: AbstractProperty} = zero(T())

promote_leaf_eltypes(φ::T) where {T <: AbstractProperty} = promote_leaf_eltypes(φ.val)

Base.convert(T::Type{TP}, φ::TP) where {TP <: AbstractProperty} = φ
Base.convert(T::Type, φ::AbstractProperty) = convert(T, φ.val)
Base.convert(T::Type{Any}, φ::AbstractProperty) = φ

# TODO: this was a hack; not clear this is a good idea ...
#       remove during next cleanup and see what happens?
Base.iterate(φ::AbstractProperty) = φ, nothing
Base.iterate(φ::AbstractProperty, ::Nothing) = nothing

"""
`coco_o_daa` : implements a tensor product between a coupling coefficient
(usually an `AbstractProperty`) and a gradient (usually a `DState`).
"""
function coco_o_daa(φ::AbstractProperty, b::TX) where {TX <: XState}
   SYMS = _syms(TX)
   vals = ntuple( i -> coco_o_daa(φ.val, _x(b)[SYMS[i]]), length(SYMS) )
   return TX( NamedTuple{SYMS}(vals) )
end

coco_o_daa(cc::SVector{N, <: AbstractProperty}, b::TX) where {N, TX <: XState} =
      SVector( ntuple(i -> coco_o_daa(cc[i], b), N) )

coco_o_daa(cc::Number, b::Number) = cc * b
coco_o_daa(cc::Number, b::SVector) = cc * b
coco_o_daa(cc::SVector, b::SVector) = cc * transpose(b)
coco_o_daa(cc::SMatrix{N1,N2}, b::SVector{N3}) where {N1,N2,N3} =
		reshape(cc[:] * transpose(b), Size(N1, N2, N3))
coco_o_daa(cc::SArray{Tuple{N1,N2,N3}}, b::SVector{N4}) where {N1,N2,N3,N4} =
		reshape(cc[:] * transpose(b), Size(N1, N2, N3, N4))

# TODO: is this needed or can it be removed?
#       maybe it should also be allowed for a DState?
#       it is also more of an ⊗ ... revisit during next cleanup
*(φ::AbstractProperty, b::AbstractState) = coco_o_daa(φ, b)

# default behaviour - overwritten e.g. for EuclideanVector where
# cocos are always complex independently of whether phi is real.
coco_type(T::Type{<: AbstractProperty}) = T

"""
`struct Invariant{D}` : specifies that the output of an ACE is
an invariant scalar.
"""
struct Invariant{T} <: AbstractProperty
   val::T
end

Base.show(io::IO, φ::Invariant) = print(io, "i($(φ.val))")

isrealB(::Invariant{<: Real}) = true
isrealB(::Invariant{<: Complex}) = false
isrealAA(::Invariant{<: Real}) = true
isrealAA(::Invariant{<: Complex}) = false

Invariant{T}() where {T <: Number} = Invariant{T}(zero(T))

Invariant(T::DataType = Float64) = Invariant{T}()

real(φ::Invariant) = Invariant(real(φ.val))
complex(φ::Invariant) = Invariant(complex(φ.val))
complex(::Type{Invariant{T}}) where {T} = Invariant{complex(T)}
complex(φ::AbstractVector{<: Invariant}) = complex.(φ)
+(φ::Invariant, x::Number) = Invariant(φ.val + x)
+(x::Number, φ::Invariant) = Invariant(φ.val + x)

# TODO: could generalize this if we have a valtype(::Type{<: AbstractProperty})
Base.convert(::Type{Invariant{T}}, x::Number) where {T} = Invariant(convert(T, x))

*(φ1::Invariant, φ2::Invariant) = Invariant(φ1.val * φ2.val)

write_dict(φ::Invariant{T})  where {T} =
   Dict("__id__" => "ACE_Invariant",
        "val" => φ.val,
        "T" => write_dict(T) )

read_dict(::Val{:ACE_Invariant}, D::Dict) =
      Invariant{read_dict(D["T"])}(D["val"])

function filter(φ::Invariant, grp::O3, b::Array)
   if length(b) <= 1
      return true
   end
   suml = sum( getl(grp, bi) for bi in b )
   if haskey(b[1], msym(grp))  # depends on context whether m come along?
      summ = sum( getm(grp, bi) for bi in b )
      return iseven(suml) && iszero(summ)
   end
   return iseven(suml)
end

filter(φ::Invariant, grp::O3O3, b::Array) =
      filter(φ, grp.G1, b) && filter(φ, grp.G2, b)

rot3Dcoeffs(::Invariant, T=Float64) = Rot3DCoeffs(T)

# TODO: this is a naive implementation of differentiation.
#       cf https://github.com/ACEsuit/ACE.jl/issues/27
#       for further discussion
*(φ::Invariant, dAA::SVector) = φ.val * dAA

coco_init(::Invariant{T}) where {T} = [ Invariant(complex(T(1)));; ]

coco_init(::Invariant, l, m, μ, T, A) = (
      l == m == μ == 0 ? Invariant(T(1)) : Invariant(T(0))  )

coco_zeros(::Invariant, ll, mm, kk, T, A) = Invariant(zero(T))

coco_filter(::Invariant, ll, mm) =
            iseven(sum(ll)) && (sum(mm) == 0)

coco_filter(::Invariant, ll, mm, kk) =
            iseven(sum(ll)) && (sum(mm) == sum(kk) == 0)

coco_dot(u1::Invariant, u2::Invariant) = u1.val * u2.val

# --------------------- EuclideanVector

@doc raw"""
`struct EuclideanVector{D, T}` : specifies that the output $\varphi$ of an
ACE is an equivariant $\varphi \in \mathbb{R}^{3}$, i.e., it transforms under
$O(3)$ as
```math
      \varphi \circ Q = Q \cdot \varphi,
```
where $\cdot$ denotes the standard matrix-vector product.
"""
struct EuclideanVector{T} <: AbstractProperty
   val::SVector{3, T}
end

Base.show(io::IO, φ::EuclideanVector) =
   print(io, "e[$(φ.val[1]), $(φ.val[2]), $(φ.val[3])]")

real(φ::EuclideanVector) = EuclideanVector(real.(φ.val))
complex(φ::EuclideanVector) = EuclideanVector(complex(φ.val))
complex(::Type{EuclideanVector{T}}) where {T} = EuclideanVector{complex(T)}

+(x::SVector{3}, y::EuclideanVector) = EuclideanVector(x + y.val)
Base.convert(::Type{SVector{3, T}}, φ::EuclideanVector) where {T} = convert(SVector{3, T}, φ.val)

isrealB(::EuclideanVector{T}) where {T} = (T == real(T))
isrealAA(::EuclideanVector) = false

# CO: removed this, is it needed?
# Base.getindex(φ::EuclideanVector, i::Integer) = φ.val[i]

EuclideanVector{T}() where {T <: Number} = EuclideanVector{T}(zero(SVector{3, T}))
EuclideanVector(T::DataType=Float64) = EuclideanVector{T}()

function filter(φ::EuclideanVector, grp::O3, bb::Array)
   if length(bb) == 0  # no zero-correlations allowed
      return false
   end
   if length(bb) == 1 #MS: Not sure if this should be here
                      # CO: good question - need to investigate
      return true
   end
   suml = sum( getl(grp, bi) for bi in bb )
   if haskey(bb[1], msym(grp))  # depends on context whether m come along?
      summ = sum( getm(grp, bi) for bi in bb )
      return isodd(suml) && abs(summ) <= 1
   end
   return isodd(suml)
end

rot3Dcoeffs(::EuclideanVector,T=Float64) = Rot3DCoeffsEquiv{T,1}(Dict[], ClebschGordan(T))

write_dict(φ::EuclideanVector{T}) where {T} =
      Dict("__id__" => "ACE_EuclideanVector",
             "val" => write_dict(Vector(φ.val)) )

function read_dict(::Val{:ACE_EuclideanVector}, D::Dict)
   return EuclideanVector(SVector{3}(read_dict(D["val"])))
end

# differentiation - cf #27
# *(φ::EuclideanVector, dAA::SVector) = φ.val * dAA'

coco_init(phi::EuclideanVector{CT}, l, m, μ, T, A) where {CT<:Real} = (
      (l == 1 && abs(m) <= 1 && abs(μ) <= 1)
         ? [EuclideanVector(conj(crmatrices[(l=l,m=-m,mu=-μ,i=i)])) for i=1:3]
         : coco_zeros(phi, l, m, μ, T, A)  )

#coco_init(phi::EuclideanVector{CT}, l, m, μ, T, A) where {CT<:Real} = (
#      (l == 1 && abs(m) <= 1 && abs(μ) <= 1)
#         ? [EuclideanVector(rmatrices[m,μ][:,k]) for k=1:3]
#         : coco_zeros(phi, l, m, μ, T, A)  )

# coco_init(phi::EuclideanVector{CT}
#   this is not needed, since the EuclideanVector should never give us
#   a constant basis anyhow. Still ... this could become a problem if we
#   ever want to artificially increase the AA basis in order to get some
#   savings elsewhere. Maybe need to revisit this...

coco_type(φ::EuclideanVector) = typeof(complex(φ))
coco_type(::Type{EuclideanVector{T}}) where {T} = EuclideanVector{complex(T)}

coco_zeros(φ::EuclideanVector, ll, mm, kk, T, A) = zeros(typeof(complex(φ)), 3)

coco_filter(::EuclideanVector, ll, mm) =
            isodd(sum(ll)) && (abs(sum(mm)) <= 1)

coco_filter(::EuclideanVector, ll, mm, kk) =
      abs(sum(mm)) <= 1 &&
      abs(sum(kk)) <= 1 &&
      isodd(sum(ll))

coco_dot(u1::EuclideanVector, u2::EuclideanVector) = dot(u1.val, u2.val)

include("eucl/cov_coeffs_dict.jl")

#---------------------- Equivariant matrices

struct EuclideanMatrix{T} <: AbstractProperty
   val::SMatrix{3, 3, T, 9}
end

function Base.show(io::IO, φ::EuclideanMatrix)
   # println(io, "3x3 $(typeof(φ)):")
   println(io, "e[ $(φ.val[1,1]), $(φ.val[1,2]), $(φ.val[1,3]);")
   println(io, "   $(φ.val[2,1]), $(φ.val[2,2]), $(φ.val[2,3]);")
   print(io,   "   $(φ.val[3,1]), $(φ.val[3,2]), $(φ.val[3,3]) ]")
end

real(φ::EuclideanMatrix) = EuclideanMatrix(real.(φ.val))
complex(φ::EuclideanMatrix) = EuclideanMatrix(complex(φ.val))
complex(::Type{EuclideanMatrix{T}}) where {T} = EuclideanMatrix{complex(T)}

+(x::SMatrix{3}, y::EuclideanMatrix) = EuclideanMatrix(x + y.val)
Base.convert(::Type{SMatrix{3, 3, T, 9}}, φ::EuclideanMatrix) where {T} =  convert(SMatrix{3, 3, T, 9}, φ.val)

isrealB(::EuclideanMatrix{T}) where {T} = (T == real(T))
isrealAA(::EuclideanMatrix) = false

#fltype(::EuclideanMatrix{T}) where {T} = T

EuclideanMatrix{T}() where {T <: Number} = EuclideanMatrix{T}(zero(SMatrix{3, 3, T, 9}))
EuclideanMatrix(T::DataType=Float64) = EuclideanMatrix{T}()

function filter(φ::EuclideanMatrix, grp::O3, bb::Array)
   if length(bb) == 0  # no zero-correlations allowed
      return false
   end
   if length(bb) == 1 #MS: Not sure if this should be here
      return true
   end
   suml = sum( getl(grp, bi) for bi in bb )
   if haskey(bb[1], msym(grp))  # depends on context whether m come along?
      summ = sum( getm(grp, bi) for bi in bb )
      return iseven(suml) && abs(summ) <= 2
   end
   return iseven(suml)
end

rot3Dcoeffs(::EuclideanMatrix,T=Float64) = Rot3DCoeffsEquiv{T,1}(Dict[], ClebschGordan(T))

write_dict(φ::EuclideanMatrix{T}) where {T} =
      Dict("__id__" => "ACE_EuclideanMatrix",
              "valr" => write_dict(real.(Matrix(φ.val))),
              "vali" => write_dict(imag.(Matrix(φ.val))),
                "T" => write_dict(T) )

function read_dict(::Val{:ACE_EuclideanMatrix}, D::Dict)
   T = read_dict(D["T"])
   valr = SMatrix{3, 3, T, 9}(read_dict(D["valr"]))
   vali = SMatrix{3, 3, T, 9}(read_dict(D["vali"]))
   return EuclideanMatrix{T}(valr + im * vali)
end

# differentiation - cf #27
# *(φ::EuclideanMatrix, dAA::SVector) = φ.val * dAA'

#coco_init(phi::EuclideanMatrix{CT}, l, m, μ, T, A) where {CT<:Real} = (
#      (l <= 2 && abs(m) <= l && abs(μ) <= l)
#         ? vec([EuclideanMatrix(conj.(transpose(mrmatrices[(m,μ,i,j)]))) for i=1:3 for j=1:3])
#         : coco_zeros(phi, l, m, μ, T, A)  )
coco_init(phi::EuclideanMatrix{CT}, l, m, μ, T, A) where {CT<:Real} = (
   (l <= 2 && abs(m) <= l && abs(μ) <= l)
      ? vec([EuclideanMatrix(conj.(mrmatrices[(l=l,m=-m,mu=-μ,i=i,j=j)])) for i=1:3 for j=1:3])
      : coco_zeros(phi, l, m, μ, T, A)  )

#coco_init(::EuclideanMatrix{CT}) where {CT<:Real} = [EuclideanMatrix(SMatrix{3,3,Complex{CT},9}([1.0,0,0,0,1.0,0,0,0,1.0]))]
coco_type(φ::EuclideanMatrix) = typeof(complex(φ))
coco_type(::Type{EuclideanMatrix{T}}) where {T} = EuclideanMatrix{complex(T)}

# This is slightly different from implementation in EuclideanVector!
coco_zeros(φ::EuclideanMatrix, ll, mm, kk, T, A) =  zeros(typeof(complex(φ)), 9)
#EuclideanMatrix.(zeros(SMatrix{3, 3, Complex{T}, 9},9))

coco_filter(::EuclideanMatrix, ll, mm) =
            iseven(sum(ll)) && (abs(sum(mm)) <= 2)

coco_filter(::EuclideanMatrix, ll, mm, kk) =
      abs(sum(mm)) <= 2 &&
      abs(sum(kk)) <= 2 &&
      iseven(sum(ll))

coco_dot(u1::EuclideanMatrix, u2::EuclideanMatrix) = sum(transpose(conj.( u1.val)) * u2.val)
#dot(u1.val, u2.val)

include("eucl/equi_coeffs_dict.jl")

# ---------------------
# some codes to help convert from svector, smatrix to spherical ...

@generated function __getL(::Val{N}) where {N}
   @assert isodd(N)
   quote
      $( div(N-1, 2) )
   end
end

__getL(::SVector{N}) where {N} = __getL(Val{N}())

__getL1L2(::SMatrix{N1, N2}) where {N1, N2} = __getL(Val{N1}()), __getL(Val{N2}())

# --------------------- SphericalVector

# TODO: remove the Val{L} - not needed (next cleanup)
struct SphericalVector{L, LEN, T} <: AbstractProperty
   val::SVector{LEN, T}
   _valL::Val{L}
end

SphericalVector(val::SVector) = SphericalVector(val, Val(__getL(val)))

# # differentiation - cf #27
# *(φ::SphericalVector, dAA::SVector) = φ.val * dAA'

isrealB(::SphericalVector) = false
isrealAA(::SphericalVector) = false

real(φ::SphericalVector) = SphericalVector(real(φ.val), φ._valL)

getL(φ::SphericalVector{L}) where {L} = L

# L = 0 -> (0,0)
# L = 1 -> (0,0), (1,-1), (1,0), (1,1)  -> 4
# L = 3 ->  ... + 5 -> 9
# 1 + 3 + 5 + ... + 2*L+1
# = L + 2 * (1 + ... + L) = L+1 + 2 * L * (L+1) / 2 = (L+1)^2
function SphericalVector(L::Integer; T = Float64)
   LEN = 2L+1   # length of SH basis up to L
   return SphericalVector( zero(SVector{LEN, T}), Val{L}() )
end

function SphericalVector{L, LEN, T}(x::AbstractArray) where {L, LEN, T}
   @assert length(x) == LEN
   SphericalVector{L, LEN, T}( SVector{LEN, T}(x), Val(L) )
end

SphericalVector{L, LEN, T}()  where {L, LEN, T} =
      SphericalVector( zero(SVector{LEN, T}), Val{L}() )

function filter(φ::SphericalVector, grp::O3, b::Array)
	if length(b) <= 1
		return true
	end
	suml = sum( getl(grp, bi) for bi in b )
   if haskey(b[1], msym(grp))
      summ = sum( getm(grp, bi) for bi in b )
      return iseven(suml) == iseven(getL(φ)) && abs(summ) <= getL(φ)
   end
   return iseven(suml) == iseven(getL(φ))
end

rot3Dcoeffs(::SphericalVector, T::DataType=Float64) = Rot3DCoeffs(T)

const __rotcoeff_inv = Rotations3D.Rot3DCoeffs(Invariant())

using ACE.Wigner: wigner_D_indices

# Equation (1.2) - vector value coupling coefficients
# ∫_{SO3} D^{ll}_{μμmm} D^*(Q) e^t dQ -> 2L+1 column vector
function vec_cou_coe(rotc::Rot3DCoeffs{T},
					      l::Integer, m::Integer, μ::Integer,
					      L::Integer, t::Integer) where {T,N}
	@assert 0 < t <= 2L+1
	D = wigner_D_indices(L)'   # Dt = D[:,t]  -->  # D^* ⋅ e^t
	LL = SA[l, L]
	Z = ntuple(i -> begin
			cc = (rotc(LL, SA[μ, D[i, t].m], SA[m, D[i, t].μ]).val)::T
			D[i, t].sign * cc
		end, 2*L+1)
	return SphericalVector{L, 2L+1, Complex{T}}(SVector(Z))
end

function _select_t(φ::SphericalVector{L}, l, M, K) where {L}
	D = wigner_D_indices(L)'
	tret = -1; numt = 0
	for t = 1:2L+1
		prodμt = prod( (D[i, t].μ + M) for i in 1:2L+1)  # avoid more allocations
		prodmt = prod( (D[i, t].m + K) for i in 1:2L+1)
		if prodμt == prodmt == 0
			tret = t; numt += 1
		end
	end
	# We assumed that there is only one coefficient; this will warn us if it fails
	@assert numt == 1
	return tret
end

coco_zeros(φ::TP, ll, mm, kk, T, A)  where {TP <: SphericalVector} = zero(TP)

coco_dot(u1::SphericalVector, u2::SphericalVector) =
		dot(u1.val, u2.val)

coco_filter(φ::SphericalVector{L}, ll, mm) where {L} =
		iseven(sum(ll) + L) && (abs(sum(mm)) <=  L)

coco_filter(φ::SphericalVector{L}, ll, mm, kk) where {L} =
      iseven(sum(ll) + L) && (abs(sum(mm)) <=  L) && (abs(sum(kk)) <= L)

coco_init(φ::SphericalVector{L}, l, m, μ, T, A) where {L} =
			vec_cou_coe(__rotcoeff_inv, l, m, μ, L, _select_t(φ, l, m, μ))

coco_init(φ::SphericalVector{L}) where {L} =
            L==0 ? reshape([SphericalVector(SVector(1.0+0.0im))],1,1) : []

# --------------- SphericalMatrix

struct SphericalMatrix{L1, L2, LEN1, LEN2, T, LL} <: AbstractProperty
   val::SMatrix{LEN1, LEN2, T, LL}
   _valL1::Val{L1}
   _valL2::Val{L2}
end

function Base.show(io::IO, φ::Union{SphericalVector, SphericalMatrix})
   buffer = IOBuffer()
   print(buffer, round.(φ.val, digits=4))
   str = String(take!(buffer))
   lenn = findfirst('[', str)
   print(io, "y" * str[lenn:end])
end

SphericalMatrix(val::SMatrix) = SphericalMatrix(val, Val.(__getL1L2(val))...)

# differentiation - cf #27
# actually this here appears to be the generic form how to do the
# differentiation for arbtirary order tensors.
# *(φ::SphericalMatrix{L1, L2, LEN1, LEN2}, dAA::SVector{N}
#       ) where {L1, L2, LEN1, LEN2, N} =
#       reshape(φ.val[:] * dAA', Size(LEN1, LEN2, N))

getL(φ::SphericalMatrix{L1,L2}) where {L1,L2} = L1, L2

isrealB(::SphericalMatrix) = false
isrealAA(::SphericalMatrix) = false

# L = 0 -> (0,0)
# L = 1 -> (0,0), (1,-1), (1,0), (1,1)  -> 4
# L = 3 ->  ... + 5 -> 9
# 1 + 3 + 5 + ... + 2*L+1
# = L + 2 * (1 + ... + L) = L+1 + 2 * L * (L+1) / 2 = (L+1)^2
function SphericalMatrix(L1::Integer, L2::Integer; T = Float64)
   LEN1 = 2L1+1   # length of SH basis up to L
   LEN2 = 2L2+1
   return SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )
end

# this is a mess - need to fix it, we should need just one constructor???

function SphericalMatrix{L1, L2, LEN1, LEN2, T, LL}(x::AbstractMatrix) where {L1, L2, LEN1, LEN2, T, LL}
   @assert size(x) == (LEN1, LEN2)
   SphericalMatrix{L1, L2, LEN1, LEN2, T, LL}( SMatrix{LEN1, LEN2, T}(x), Val(L1), Val(L2) )
end

SphericalMatrix{L1, L2, LEN1, LEN2, T, LL}()  where {L1, L2, LEN1, LEN2, T, LL} =
      SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )

SphericalMatrix{L1, L2, LEN1, LEN2, T}()  where {L1, L2, LEN1, LEN2, T, LL} =
		SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )

function filter(φ::SphericalMatrix, grp::O3, b::Array)
	if length(b) < 1
		return true
	end
	suml = sum( getl(grp, bi) for bi in b )
   if haskey(b[1], msym(grp))
      summ = sum( getm(grp, bi) for bi in b )
      return iseven(suml) == iseven( sum(getL(φ)) ) && abs(summ) <= sum(getL(φ))
   end
   return iseven(suml) == iseven( sum(getL(φ)) )
end

rot3Dcoeffs(::SphericalMatrix, T::DataType=Float64) = Rot3DCoeffs(T)

function mat_cou_coe(rotc::Rot3DCoeffs{T},
				   		l::Integer, m::Integer, μ::Integer,
					   	a::Integer, b::Integer,
							::Val{L1}, ::Val{L2}) where {T, L1, L2}
	@assert (0 < a <= 2L1 + 1) && (0 < b <= 2L2 + 1)
	Z = zero(MMatrix{2L1+1, 2L2+1, Complex{T}})  # zeros(2 * L1 + 1, 2 * L2 + 1)
	Dp = wigner_D_indices(L1)'
	Dq = wigner_D_indices(L2)
	LL = SA[l, L1, L2]
	for i = 1:(2 * L1 + 1)
		for j = 1:(2 * L2 + 1)
			MM = SA[μ, Dp[i,a].m, Dq[b,j].m]
			KK = SA[m, Dp[i,a].μ, Dq[b,j].μ]
			cc = (rotc(LL, MM, KK).val)::T
			Z[i,j] = Dp[i,a].sign * Dq[b,j].sign * cc
		end
	end
	return SphericalMatrix(SMatrix(Z), Val{L1}(), Val{L2}())
	# return SphericalMatrix(SMatrix{2L1+1,2L2+1,Complex{T}}(Z), Val{L1}(), Val{L2}())
end

function _select_ab(φ::SphericalMatrix{L1,L2}, M, K) where {L1,L2}
   Dp = wigner_D_indices(L1)'
   Dq = wigner_D_indices(L2)
   list_ab = Tuple{Int, Int}[]
   for a = 1:2L1+1
      for b = 1:2L2+1
			# pm = prod( ma[i] + mb[j] + K for i = 1:2L1+1, j = 1:2L2+1)
			pm = prod( Dp[i,a].m + Dq[b,j].m + K for i = 1:2L1+1, j = 1:2L2+1)
			# pμ = prod(μa[i] + μb[j] + M for i = 1:2L1+1, j = 1:2L2+1)
			pμ = prod( Dp[i,a].μ + Dq[b,j].μ + M for i = 1:2L1+1, j = 1:2L2+1)
         if pμ == pm ==0
				push!(list_ab, (a,b))
         end
      end
   end
	return list_ab
end

function coco_init(φ::SphericalMatrix{L1,L2}, l, m, μ, T, A) where{L1,L2}
   list = _select_ab(φ, m, μ)
	@assert length(list) > 0
	# MAIN CASE
   if abs(m) <= L1+L2 && abs(μ) <= L1+L2
		return [ mat_cou_coe(__rotcoeff_inv, l, m, μ, a, b, Val(L1), Val(L2))
				   for (a,b) in list ]
	end

	# @warn("SHOULDN'T BE HERE!!")
   return fill( zero(typeof(φ)), length(list) )
end

coco_init(φ::SphericalMatrix{L1,L2}) where{L1,L2} =
				  L1==L2 ? reshape([ACE.SphericalMatrix(SMatrix{2L1+1,2L2+1,ComplexF64}(I(2L1+1)),Val(L1),Val(L2))],1,1) : []

coco_zeros(φ::TP, ll, mm, kk, T, A) where{TP <: SphericalMatrix} =
            zeros(TP, length(_select_ab(φ, sum(mm), sum(kk))))

coco_filter(φ::SphericalMatrix{L1,L2}, ll, mm) where {L1,L2} =
            iseven(sum(ll)) == iseven(L1+L2) && (abs(sum(mm)) <=  L1+L2)

coco_filter(φ::SphericalMatrix{L1,L2}, ll, mm, kk) where {L1,L2} =
            iseven(sum(ll) + L1 + L2) &&
            (abs(sum(mm)) <=  L1+L2) &&
            (abs(sum(kk)) <=  L1+L2)

coco_dot(u1::SphericalMatrix, u2::SphericalMatrix) =
		dot(u1.val, u2.val)

#  Some promotion rules

Base.promote_rule(::Type{T1}, ::Type{Invariant{T2}}
                  ) where {T1 <: Number, T2 <: Number} =
      Invariant{promote_type(T1, T2)}

Base.promote_rule(::Type{T1}, ::Type{EuclideanVector{T2}}
                  ) where {T1 <: Number, T2 <: Number} =
      EuclideanVector{promote_type(T1, T2)}

Base.promote_rule(t1::Type{T1}, t2::Type{SVector{N, T2}}
                  ) where {N, T1 <: Number, T2 <: AbstractProperty} =
      SVector{N, promote_rule(t1, T2)}

# --------------------------- AD related codes

# an x -> x.val implementation with custom adjoints to sort out the
# mess created by the AbstractProperties
# maybe this feels a bit wrong, definitely a hack. What might be nicer
# is to introduce a "Dual Property" similar to the "DState"; Then we
# could have something along the lines of  DProp * Prop = scalar or
# _contract(DProp, Prop) = scalar; That would be the "systematic" and
# "disciplined" way of implementing this.

"""
`val(x) = x.val`, normally to be used if x is a property. This should be used
instead of x.val when the operation is part of a bigger expression that is
to be ADed. I.e. `val` has rrules implemented that should allow taking up
to two derivatives.

TODO: at the moment this is a bit hacky, and needs to be adjusted over time
as we learn more about how to best implement AD. The real question is what
is the correct adjoint of this operation? E.g., if
````
     val : { Invariant{T} } -> { T }
```
then should
```
   val* :  { T } -> { Invariant{T} }
```
? If this is the case, then we need a parameterised `val` function. I.e.
we need to remember what the Property is that we started from.
"""
val(x) = x.val

function _rrule_val(dp, x)     # ∂/∂x (dp * x) = dp
   @assert dp isa Number
   return NoTangent(), dp
end

rrule(::typeof(val), x) =
         val(x),
         dp -> _rrule_val(dp, x)

function rrule(::typeof(_rrule_val), dp, x)   # ∂/∂... (0 + dp * dq[2])
      @assert dp isa Number
      function second_adj(dq)
         @assert dq[1] == ZeroTangent()
         # @assert dq[2] isa Number  # TODO -> revisit this?!
         return NoTangent(), dq[2], ZeroTangent()
      end
      return _rrule_val(dp, x), second_adj
end

import ChainRulesCore: ProjectTo

(::ProjectTo{T})(φ::Invariant{T}) where {T} = val(φ)


text:
# Boosting
function boosted_tree(dtmtable::DTMTable, sett::ModelSettings)
	trnidx = dtmtable.trnidx
	validx = dtmtable.validx
	features = dtmtable.features
	weight = dtmtable.weight
	actualNumerator = dtmtable.numerator
	denominator = dtmtable.denominator
	mappings = dtmtable.mappings
	candMatWOMaxValues = dtmtable.candMatWOMaxValues
	obstrn, obsval, trn_meanobservedvalue, val_meanobservedvalue, trn_numtot, val_numtot, trn_denomtot, val_denomtot, empty_rows_after_iteration_stats, showTimeUsedByEachIteration, chosen_apply_tree_fn, startAtMean, adaptiveLearningRate, moderationvector, iterations, nameOflistofpredictorsSheet, nameOfpredictorsSheet, nameOfModelStatisticsSheet, nameOfScoresSheet, nameOfOverviewSheet, nameOfSettingsSheet, nameOfValidationSheet, xlData, showProgressBar_time, prroduceEstAndLeafMatrices = initSettingsWhichAreTheSameForBoostingAndBagging(trnidx, validx, actualNumerator, denominator, sett)
	# trn_meanobservedvalue is the mean observed RATIO, i.e. sum(numerator)/sum(denominator) for the training data
	T_Uint8_or_UInt16 = find_max_type(features)
    # obs=obstrn+obsval
    obs = size(features, 1)
    # obs=obstrn+obsval  assumes that trn and val idx form a partition!
    if size(features, 1) != obstrn + obsval
#        @warn("DTM: obstrn+obsval != size(features,1). This is experimental.")
    end

	current_error = 0.0
    moderationfactor = moderationvector[1]
	if size(moderationvector, 1) == 1
		actual_moderationvector = Array{Float64}(undef, iterations)
		fill!(actual_moderationvector, copy(moderationvector[1]))
	else
		actual_moderationvector = copy(moderationvector)
	end
	if startAtMean
		estimatedRatio = ones(length(weight)) .* trn_meanobservedvalue
	end

	estimatedNumerator = estimatedRatio .* denominator # these are for instance the estimated losses for a LR model

	estimatedNumeratorForStats = zeros(obs)
    indicatedRelativityForApplyTree_reused = zeros(obs)
	reused_fitted_leafnr_vector = zeros(Int, obs)
    sortvec_reused_trn_only = zeros(Int, length(trnidx))
    intVarsUsed = Array{Array{Array{Int,1},1}}(undef, 0);sizehint!(intVarsUsed, iterations)
	inds_considered = Array{Array{Int,1}}(undef, 0);sizehint!(inds_considered, iterations)

	res = Array{Union{Leaf{T_Uint8_or_UInt16},Node{UInt8},Node{UInt16}}}(undef, iterations)
	if prroduceEstAndLeafMatrices
		est_matrix = Array{Float64}(undef, obs, iterations + 1)
		est_matrixFromScores = copy(est_matrix)
		MatrixOfLeafNumbers = Array{Int}(undef, obs, iterations + 1)
		MatrixOfLeafNumbers[:,1] .= 0
		est_matrix[:,1] = copy(transpose(estimatedRatio))
		est_matrixFromScores[:,1] = copy(transpose(estimatedRatio))
	else
		est_matrix = Array{Float64}(undef, 0, 0)
		est_matrixFromScores = copy(est_matrix)
		MatrixOfLeafNumbers = Array{Int}(undef, 0, 0)
	end

	scores = zeros(Int, obs)
	estFromScores = zeros(obs)
	vectorOfLeafArrays = Array{Array{Leaf{T_Uint8_or_UInt16},1}}(undef, iterations + 1)
	vectorOfLeafArrays[1] = Array{Leaf{T_Uint8_or_UInt16}}(undef, 0) # the first entry is not defined

    p = ProgressMeter.Progress(iterations, 2, "Progress of Boosting Model:") # minimum update interval: x seconds (2)
    ((adaptiveLearningRate >= 1.0) || (adaptiveLearningRate <= 0.0)) ? (BoolAdaptiveLearningRate = false) : (BoolAdaptiveLearningRate = true)
	# currently not supported
		# BoolAdaptiveLearningRate&&error("DTM: this is not yet updated to reflect numerator/denominator objectives!!")
	# this header needs to be in line with the output of the function createTrnValStatsForThisIteration
	statsPerIteration = global_statsperiter_header
	currentRelativity = Array{Float64}(undef, obs)
	scoreBandLabels = createScorebandsUTF8List(sett.scorebandsstartingpoints, sett.nScores, addtotal=true)
	# this row is only here such that the variables are initialized outside the for loop
	local fristRowOfThisTable, obsPerScore, vectorWeightPerScore, numPerScore, denomPerScore, scores, scoresVAL, cumulativeStatsPerScoreBand, maxRawRelativityPerScoreSorted, rawObservedRatioPerScore, MAPPINGSmoothedEstimatePerScore

	# the following two lines are for the initialization in case any subsampling (boostrapping) of data is done
	# the idea is to preallocate many of the arrays which are used for sampling (and to save allocation/memory space)
		sampleSizeCanBeNEGATIVE = convert(Int, round(sett.subsampling_prop * obstrn))
		if abs(sampleSizeCanBeNEGATIVE) < 1
			sampleSizeCanBeNEGATIVE = convert(Int, sign(sett.subsampling_prop))
		end
		abssampleSize, sampleVector = initBootstrapSample(sampleSizeCanBeNEGATIVE) # ,trnidx,validx)

        # local stats,estimatedRatioUnsmoothed,estimatedRatioSmoothed,estimateRatioFromRelativities,resultingBT

    for iter = 1:iterations
        # currently disable timeConditional as it leads to a segfault (see https://github.com/JuliaLang/julia/issues/28536) \August 15, 2018
        # notably this segfault did not appear on 0.7beta2 (and earlier) where all tests ran perfectly fine (see travis)
        # @timeConditional(showTimeUsedByEachIteration, begin
		# Build Iteration iter
			current_mdf = moderationvector[min(iter, size(moderationvector, 1))] # mdf remains constant if a vector of size 1 was provided (instead of a vector of size iterations)

			# @show iter,"here"
			# @show size(trnidx)
			if showTimeUsedByEachIteration
			# todo tbd: this can be done more efficiently: We can avoid the tmpTree object (also we do not need to call "sometreesettings" every time!, improve this
				@time 	tmpTree = sample_data_and_build_tree!(trnidx, validx, indicatedRelativityForApplyTree_reused, candMatWOMaxValues, mappings, deepcopy(sett), actualNumerator, estimatedNumerator, weight, features, sampleSizeCanBeNEGATIVE, abssampleSize, sampleVector, T_Uint8_or_UInt16)
			else
						tmpTree = sample_data_and_build_tree!(trnidx, validx, indicatedRelativityForApplyTree_reused, candMatWOMaxValues, mappings, deepcopy(sett), actualNumerator, estimatedNumerator, weight, features, sampleSizeCanBeNEGATIVE, abssampleSize, sampleVector, T_Uint8_or_UInt16)
			end
			# indicatedRelativity
			if minimum(indicatedRelativityForApplyTree_reused) < 0.0
				# todo/tbd check this: if there are negative indicated relativites the boosting approach is not meaningful.
				@warn("The boosting model encountered negative indicated realtivities for some leaves.")
				@warn("This should not happen if you have only positive values for numerator and denominator.")
				error("Abort.")
			end
			# todo/tbd check this: what if a fitted value in a leaf is 0? what will the indicatedRelativityForApplyTree_reused be? will the boosting work as expected?
			push!(intVarsUsed, tmpTree.intVarsUsed)
			push!(inds_considered, tmpTree.inds_considered)
			res[iter] = tmpTree.rootnode # we do not need the full tree (with all its settings and so on) here, todo/tbd (someone else should) check if this is ok and creates no issues
			vectorOfLeafArrays[iter + 1] = create_leaves_array(res[iter])
			# vectorOfRulePathsToLeavesArrays[iter+1]=[x.rule_path for x in vectorOfLeafArrays[iter+1]]
			# Apply this iteration to the validation data set
			apply_tree_by_leaf!(indicatedRelativityForApplyTree_reused, reused_fitted_leafnr_vector, validx, res[iter], features)
			# update_part_of_this_vector!(validx,indicatedRelativityForApplyTree_reused,indicatedRelativity)
        if prroduceEstAndLeafMatrices
            # todo/tbd this could possibly be optimized. We are calling apply_tree_by_leaf! twice (if prroduceEstAndLeafMatrices==true), once for val and once for trn
            apply_tree_by_leaf!(indicatedRelativityForApplyTree_reused, reused_fitted_leafnr_vector, trnidx, res[iter], features)
            MatrixOfLeafNumbers[:,iter + 1] = reused_fitted_leafnr_vector # leaf_numbers(vectorOfLeafArrays[1+iter],obs) #TODO / TBD adjust this for subsampling. this may not work properly, as each tree will use a different amount of data AND DIFFERENT OBSERVATIONS! thus we cannot rely on the *.idx field of the leaves!
        end
			# Moderate the estimate, this is done for BOTH trn and val!
			_moderate!(estimatedRatio, indicatedRelativityForApplyTree_reused, current_mdf)
			if prroduceEstAndLeafMatrices
				write_column!(est_matrix, iter + 1, estimatedRatio)
			end
			update_current_rels!(currentRelativity, estimatedRatio, trn_meanobservedvalue)
			# update estimatedNumerator
			estimatedNumerator = estimatedRatio .* denominator
		# Derive Scores for this iteration #NOTE (tbd/todo: keep in mind) For large datasets (>5m rows) the sorting in construct scores may become a dominant part of the algorithm
			maxRawRelativityPerScoreSorted, MAPPINGSmoothedEstimatePerScore, vectorWeightPerScore, obsPerScore, rawObservedRatioPerScore, numPerScore, denomPerScore, nscoresPotentiallyReduced = constructANDderiveScores!(trnidx, validx, sortvec_reused_trn_only, estimatedRatio, currentRelativity, actualNumerator, denominator, weight, trn_meanobservedvalue, iter, sett)
			update_and_derive_scores!(scores, maxRawRelativityPerScoreSorted, currentRelativity)
			# todo, this could probably be improved
			# if we avoid creating small negative values in the function constructANDderiveScores altogether, then the following two function calls are obsolete
			correctSmallNegativeValues!(MAPPINGSmoothedEstimatePerScore)
			correctSmallNegativeValues!(rawObservedRatioPerScore)
		# Define estimated Numerator
			sett.smoothEstimates == "0" ? referenceForNumEstimates = rawObservedRatioPerScore : referenceForNumEstimates = MAPPINGSmoothedEstimatePerScore
			update_mapped_estFromScores!(estFromScores, referenceForNumEstimates, scores)

			if prroduceEstAndLeafMatrices
				write_column!(est_matrixFromScores, iter + 1, estFromScores)
			end
		# Derive Cumulative Statistics
			 # todo tbd there is potential for optimization here
			 # ["rawRelativities","unsmoothedPerScore","smoothedPerScore"]
			if sett.fitForStatsAndCharts == "rawRelativities"
				selectedEstimatedRatioForStats = currentRelativity .* trn_meanobservedvalue
			else
				selectedEstimatedRatioForStats = estFromScores::Vector{Float64} # smooth or unsmooth has been defined/"selected" above already
			end
        for ij = 1:length(estimatedNumeratorForStats) # a.=b.*c seemed to allocate memory ..... (see memhist of June 24 2018)
            @inbounds  estimatedNumeratorForStats[ij] = selectedEstimatedRatioForStats[ij] * denominator[ij]
        end

			statsThisIteration, singleRowWithKeyMetrics, columnOfRelativityTrn = createTrnValStatsForThisIteration(scoreBandLabels, iter, sett.scorebandsstartingpoints, view(actualNumerator, trnidx), view(denominator, trnidx), view(weight, trnidx), view(selectedEstimatedRatioForStats, trnidx), view(scores, trnidx), view(actualNumerator, validx), view(denominator, validx), view(weight, validx), view(selectedEstimatedRatioForStats, validx), view(scores, validx), sett)
			statsPerIteration = vcat(statsPerIteration, singleRowWithKeyMetrics)
			if iter == 1
				fristRowOfThisTable = sett.iterations + 2 + empty_rows_after_iteration_stats
				cumulativeStatsPerScoreBand = copy(statsThisIteration)
			else
				fristRowOfThisTable += size(statsThisIteration, 1)
				cumulativeStatsPerScoreBand = vcat(cumulativeStatsPerScoreBand, statsThisIteration)
			end
			# Add the respective chart
				loc = string(excelLetter(size(cumulativeStatsPerScoreBand, 2) + 2), fristRowOfThisTable)
				thischart = defineRelativityChart(nameOfModelStatisticsSheet, nameOfModelStatisticsSheet, loc, size(statsThisIteration, 1) - 2, columnOfRelativityTrn, fristRowOfThisTable, xscale=2.0)
				push!(xlData.charts, deepcopy(thischart))
		# end	) #end of @timeConditional
		if showProgressBar_time # TimeUsedByEachIteration
			ProgressMeter.next!(p)
		end
    end
	# GC.enable(true)	#to check: is this needed? is there any benefit to it?
	# todo/tbd: note that we have at least three different estimates should we add all of them to the output?

	# one of these vectors is actually identical to estFromScores and would not need to be recalculated
	estimatedRatioUnsmoothed = map_these_values(rawObservedRatioPerScore, scores)
	estimatedRatioSmoothed = map_these_values(MAPPINGSmoothedEstimatePerScore, scores)
	estimateRatioFromRelativities = currentRelativity .* trn_meanobservedvalue

	# fitForStatsAndCharts is an element of ["rawRelativities","unsmoothedPerScore","smoothedPerScore"]
	if sett.fitForStatsAndCharts == "rawRelativities"
		selectedEstimatedRatioForStats = estimateRatioFromRelativities
	end
	if sett.fitForStatsAndCharts == "unsmoothedPerScore"
		selectedEstimatedRatioForStats = estimatedRatioUnsmoothed
	end
	if sett.fitForStatsAndCharts == "smoothedPerScore"
		selectedEstimatedRatioForStats = estimatedRatioSmoothed
	end
	estimatedNumeratorForStats = selectedEstimatedRatioForStats .* denominator

	# add empty lines
	statsPerIteration = vcat(statsPerIteration, repeat([""], empty_rows_after_iteration_stats, size(statsPerIteration, 2)))
	# Add Charts which show stats for all iterations
		add_iteration_charts!(xlData, sett, 1)
		size(statsPerIteration, 2) > size(cumulativeStatsPerScoreBand, 2) ? cumulativeStatsPerScoreBand = hcat(cumulativeStatsPerScoreBand, repeat([""], size(cumulativeStatsPerScoreBand, 1), size(statsPerIteration, 2) - size(cumulativeStatsPerScoreBand, 2))) : statsPerIteration = hcat(statsPerIteration, repeat([""], size(statsPerIteration, 1), -size(statsPerIteration, 2) + size(cumulativeStatsPerScoreBand, 2)))
		stats = vcat(statsPerIteration, cumulativeStatsPerScoreBand)
# Create Score List
	scoreheader = ["Score" "Row Count" "Weight" "Numerator" "Denominator" "Maximal Raw Relativity" "Observed Ratio" "Smoothed Fit"]
	scoreMatrix = vcat(scoreheader, hcat(collect(1:size(maxRawRelativityPerScoreSorted, 1)), obsPerScore, vectorWeightPerScore, numPerScore, denomPerScore, maxRawRelativityPerScoreSorted, rawObservedRatioPerScore, MAPPINGSmoothedEstimatePerScore))
	scoreMatrix = hcat(scoreMatrix, repeat([""], size(scoreMatrix, 1), size(stats, 2) - size(scoreMatrix, 2)))
# Add Score Chart
	thischart = defineScoreChart(nameOfScoresSheet, nameOfScoresSheet, "B6", size(scoreMatrix, 1) - 1, 1, 3, 7, 8)
	push!(xlData.charts, deepcopy(thischart))
# Create univariate graphs
		predictorsData, predictorCharts = createPredictorData(trnidx, nameOfpredictorsSheet, mappings, candMatWOMaxValues, sett, scores, actualNumerator, denominator, weight, selectedEstimatedRatioForStats, features)
	# push predictorCharts into Excelobject
		push!(xlData.charts, predictorCharts...)
		predictorsSheet = ExcelSheet(nameOfpredictorsSheet,  DataFrame(predictorsData,:auto))
	# Populate ExcelSheet object with data
		statsSheet = ExcelSheet(nameOfModelStatisticsSheet, DataFrame(stats,:auto))
		overviewSheet = ExcelSheet(nameOfOverviewSheet, DataFrame())
		scoreMatrixDF = DataFrame(scoreMatrix,:auto)
		scoresSheet = ExcelSheet(nameOfScoresSheet, scoreMatrixDF)
		modelsettingsSheet = ExcelSheet(nameOfSettingsSheet, DataFrame( writeAllFieldsToArray(sett),:auto))
# Create two way validation charts
	twoWayValidation, validationCharts = createTwoWayValidationCharts(trnidx, validx, nameOfValidationSheet, scoreBandLabels, mappings, candMatWOMaxValues, sett, scores, actualNumerator, denominator, weight, selectedEstimatedRatioForStats, features)
	validationSheet = ExcelSheet(nameOfValidationSheet, DataFrame(twoWayValidation,:auto))
	push!(xlData.charts, validationCharts...)
# Define Excel
	xlData.sheets = [modelsettingsSheet,overviewSheet,statsSheet,scoresSheet,predictorsSheet,validationSheet]
# add boolNumeratorStats if requested
	if sett.boolNumeratorStats
		@warn("tariff est stats are not updated for DTM2: this needs review")
			errors_num_estimates = addTariffEstimationStatsAndGraphs!(xlData, trnidx, validx, actualNumerator, estimatedRatioUnsmoothed, estimatedRatioSmoothed, estimateRatioFromRelativities, est_matrix, est_matrixFromScores)
		# attach errors per iteration to modelStatistics== statsSheet (Excel file)
			statsSheetINDEX = findall(in([nameOfModelStatisticsSheet]), [x.name for x in xlData.sheets])
			@assert length(statsSheetINDEX) == 1 "Error no statsheet was found in Excel data (when trying to attach Numerator Errors to Modelstatistics)"
			statsSheetINDEX2 = statsSheetINDEX[1]
			errors_num_estimates = vcat(errors_num_estimates, repeat([""], size(stats, 1) - size(errors_num_estimates, 1), size(errors_num_estimates, 2)))
			stats = hcat(stats, errors_num_estimates)
			xlData.sheets[statsSheetINDEX2] = ExcelSheet(nameOfModelStatisticsSheet, DataFrame( stats,:auto))
	end
	z = deepcopy(stats[1:sett.iterations + 1,:])
	modelstats = DataFrame(convert(Array{Float64,2}, z[2:size(z, 1),:]),:auto)
	DataFrames.rename!(modelstats, Symbol[Symbol(x) for x in view(z, 1, :)])
# resulting BT
	# create trnidx
	@assert issorted(trnidx)
	trnidx_one_zero_full_length = map(x->UInt8(length(searchsorted(trnidx, x))), 1:length(scores))
	fp = get_feature_pools(features)
	resultingBT = BoostedTree(res, sett, intVarsUsed, candMatWOMaxValues, mappings, inds_considered, actual_moderationvector, scores, currentRelativity, maxRawRelativityPerScoreSorted, trn_meanobservedvalue, startAtMean, MAPPINGSmoothedEstimatePerScore, rawObservedRatioPerScore, est_matrix, modelstats, xlData, trnidx_one_zero_full_length, fp)
	return xlData, estimatedRatio, MatrixOfLeafNumbers, vectorOfLeafArrays, rawObservedRatioPerScore, est_matrixFromScores, stats, estimatedRatioUnsmoothed, estimatedRatioSmoothed, estimateRatioFromRelativities, resultingBT
end

ext: jl


text:
export
    NORMAL,
    TONE,
    TONE2,
    INITIALS,
    FIRST_LETTER,
    FINALS,
    FINALS_TONE,
    FINALS_TONE2,
    TONE3,
    FINALS_TONE3,
    BOPOMOFO,
    BOPOMOFO_FIRST,
    CYRILLIC,
    CYRILLIC_FIRST

export
    TONE12345,
    TONEONTOP,
    TONELESS

const NORMAL         = 0    # 中国 -> zhong guo     ⤎ 不带声调
const TONE           = 1    # 中国 -> zhōng guó     ⤎ 标准声调
const TONE2          = 2    # 中国 -> zho1ng guo2   ⤎ 声调在韵母后
const INITIALS       = 3    # 中国 -> zh g          ⤎ 只返回声母
const FIRST_LETTER   = 4    # 中国 -> z g           ⤎ 只返回拼音首字母
const FINALS         = 5    # 中国 -> ong uo        ⤎ 只返回韵母部分,不带声调
const FINALS_TONE    = 6    # 中国 -> ōng uó        ⤎ 只返回带标准声调的韵母
const FINALS_TONE2   = 7    # 中国 -> o1ng uo2      ⤎ 只返回带数字声调的韵母
const TONE3          = 8    # 中国 -> zhong1 guo2   ⤎ 数字声调在拼音末尾
const FINALS_TONE3   = 9    # 中国 -> ong1 uo2      ⤎ 只返回带数字声调的韵母,声调在末尾
const BOPOMOFO       = 10   # 中国 -> ㄓㄨㄥ ㄍㄨㄛˊ ⤎ 注音风格,带声调,第一声不标
const BOPOMOFO_FIRST = 11   # 中国 -> ㄓ ㄍ         ⤎ 仅首字母注音
const CYRILLIC       = 12   # 中国 -> чжун1 го2    ⤎ 俄语字母风格,数字声调在末尾
const CYRILLIC_FIRST = 13   # 中国 -> ч г          ⤎ 俄语字母风格,仅首字母

# most popular ones
TONE12345 = TONE3  # ni1 hao3 style, having 1 to 5 as five tones
TONEONTOP = TONE   # ní hǎo   style, having tone on top of vowel
TONELESS  = NORMAL # ni hao   style, having no tones


text:
@testset "230.kth-smallest-element-in-a-bst.jl" begin
    @test kth_smallest_in_BST(TreeNode{Int}([3, 1, 4, nothing, 2]), 1) == 1
    @test kth_smallest_in_BST(TreeNode{Int}([5, 3, 6, 2, 4, nothing, nothing, 1]), 3) == 3
end


text:
using Oceananigans.Utils: work_layout
using Oceananigans.Architectures: device
using Oceananigans.TimeSteppers: store_tracer_tendency!

import Oceananigans.TimeSteppers: store_tendencies!

""" Store source terms for `uh`, `vh`, and `h`. """
@kernel function store_solution_tendencies!(G⁻, grid, G⁰)
    i, j, k = @index(Global, NTuple)
    @unroll for t in 1:3
        @inbounds G⁻[t][i, j, k] = G⁰[t][i, j, k]
    end
end

""" Store previous source terms before updating them. """
function store_tendencies!(model::ShallowWaterModel)

    barrier = Event(device(model.architecture))

    workgroup, worksize = work_layout(model.grid, :xyz)

    store_solution_tendencies_kernel! = store_solution_tendencies!(device(model.architecture), workgroup, worksize)
    store_tracer_tendency_kernel! = store_tracer_tendency!(device(model.architecture), workgroup, worksize)

    solution_event = store_solution_tendencies_kernel!(model.timestepper.G⁻,
                                                       model.grid,
                                                       model.timestepper.Gⁿ,
                                                       dependencies=barrier)

    events = [solution_event]

    # Tracer fields
    for i in 4:length(model.timestepper.G⁻)
        @inbounds Gc⁻ = model.timestepper.G⁻[i]
        @inbounds Gc⁰ = model.timestepper.Gⁿ[i]
        tracer_event = store_tracer_tendency_kernel!(Gc⁻, model.grid, Gc⁰, dependencies=barrier)
        push!(events, tracer_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end


text:
"""
S2S(net::Function) implements the sequence to sequence model from:
Sutskever, I., Vinyals, O., & Le, Q. V. (2014). Sequence to sequence
learning with neural networks. Advances in neural information
processing systems, 3104-3112.  The @knet function `net` is duplicated
to create an encoder and a decoder.  See `S2SData` for data format.
"""
immutable S2S <: Model; encoder; decoder; params;
    function S2S(kfun::Function; o...)
        enc = Net(encoder; o..., f=kfun)
        dec = Net(decoder; o..., f=kfun)
        new(enc, dec, vcat(params(enc), params(dec)))
    end
end

params(m::S2S)=m.params
reset!(m::S2S; o...)=(reset!(m.encoder; o...);reset!(m.decoder; o...))

@knet function encoder(word; f=:lstm, hidden=0, o...)
    wvec = wdot(word; o..., out=hidden)
    hvec = f(wvec; o..., out=hidden)
end

@knet function decoder(word; f=:lstm, hidden=0, vocab=0, o...)
    hvec = encoder(word; o..., f=f, hidden=hidden)
    tvec = wdot(hvec; out=vocab)
    pvec = soft(tvec)
end

train(m::S2S, data, loss; o...)=s2s_loop(m, data, loss; trn=true, ystack=Any[], o...)

test(m::S2S, data, loss; o...)=(l=zeros(2); s2s_loop(m, data, loss; losscnt=l, o...); l[1]/l[2])

# todo: implement load/save: revive cpucopy - done.
# todo: implement predict
# data should be a sequence generator for the source side
# i.e. next(data) should be a vector of source side ints
# todo: introduce an abstract type for data
# done: we need to have our own batch matrices
# done: the batchsize for the model needs to change
# note: no need for mask if we use batchsize=1 for predict
# note: s2sdata should probably be part of train, since we'll have to write different train functions for different models anyway, how the model consumes data is not part of the data, it is part of train, will make it easier to write data generators.
# done: copyforw needs to do repmat
# done: eos(sgen) different for x2
# todo: nbest should be another parameter, for now we assume nbest=1

function predict(m::S2S, sgen; beamsize=1, dense=false, nbest=1, o...)
    clear!(m)
    p = params(m)
    ftype = eltype(p[1].out)
    eos1 = vocab1 = size(p[1].out, 2)
    eos2 = vocab2 = size(p[end].out, 1)
    # @assert eos1 == eos(sgen)
    if dense
        x1 = zeros(ftype, vocab1, 1)
        x2 = zeros(ftype, vocab2, beamsize)
    else
        x1 = sponehot(ftype, vocab1, 1)
        x2 = sponehot(ftype, vocab2, beamsize)
    end
    pred = Any[]
    for seq in sgen
        reset!(m)
        forw(m.encoder, setrow!(x1,eos1,1); trn=false, seq=true, o...)
        for x in reverse(seq)
            forw(m.encoder, setrow!(x1, x, 1); trn=false, seq=true, o...)
        end
        for i in 1:beamsize
            setrow!(x2, eos2, i)
        end
        #@show display(to_host(m.encoder.out[end]))
        initforw(m.decoder, x2)
        s2s_copyforw!(m.encoder, m.decoder, [1 for i=1:beamsize])
        push!(pred, decode(m.decoder, x2; nbest=nbest, o...))
    end
    return pred
end

# done: construct the new h, need a shuffling version of copyforw, need to be careful not to overwrite
# done: construct the new x, need to find the top scoring words without repetition, these are cumulative scores, not word scores! and cumulative words, e.g. identical words may win continuing different sequences.  we can use bad initial scores for all but the first of the initial columns.
# todo: generate the max score sequence (already have prob in y) no need to log all: we do need to log all to combine with previous logp
# done: n-best is a different parameter than beamsize, implement n-best?
# done: pop the finished sequences from the beam, stop when everything on the beam is worse than top-n in the popped sequences.
# todo: give back the answer in string
# todo: add dropout

function decode(r::Net, x; nbest=1, o...)
    (nrows,ncols) = size(x)
    eos = vocab = nrows
    # TODO: some of these should not need to be allocated every time
    global ylogp = CudaArray(eltype(x), size(x))
    global yscore = Array(eltype(x), size(x))
    # beam is going to hold the previous tokens leading up to the current columns
    global beam = [ Int[] for i=1:ncols ]
    # score is the cumulative logp of sequences on the beam, we initialize all but the first column with -Inf because initially all columns are going to be the same.
    global score = Array(eltype(x), (1, ncols))
    fill!(score, -Inf); score[1]=0
    # nbest list for completed sequences, keys: Int[], values: cumulative logp.  When scores on the beam are all worse than the ones in nbest stop.
    global nbestqueue = PriorityQueue()

    while true
        #@show x
        #@show display(to_host(r.out[end-3]))
        yprob = forw(r, x; trn=false, seq=true, o...)
        log!(yprob, ylogp)
        copysync!(yscore, 1, ylogp, 1, length(yscore))
        broadcast!(+, yscore, score, yscore)
        global ytop = topn(yscore, 2*ncols)
        global newbeam,newscore,oldcol,newn
        newbeam,newscore,oldcol,newn = similar(beam),similar(score),zeros(Int,ncols),0
        global tmp = similar(x)        # todo: avoid alloc
        for i in ytop
            (iword, ibeam) = ind2sub(yscore, i)
            # println("yscore$((beam[ibeam]..., iword))=$(yscore[i])")
            if iword == eos
                if length(nbestqueue) < nbest
                    nbestqueue[beam[ibeam]] = yscore[i]
                elseif yscore[i] > peek(nbestqueue).second
                    dequeue!(nbestqueue)
                    nbestqueue[beam[ibeam]] = yscore[i]
                end
            else
                newn += 1
                oldcol[newn] = ibeam
                newbeam[newn] = vcat(beam[ibeam],iword)
                newscore[newn] = yscore[i]
                setrow!(tmp, iword, newn)
                newn == ncols && break
            end
        end
        # todo: what if there is not enough to fill ncols?
        @assert newn == ncols
        beam,score = newbeam,newscore
        length(nbestqueue) >= nbest && score[1] < peek(nbestqueue).second && break
        s2s_copyforw!(r, r, oldcol)
        x = tmp
    end
    @assert length(nbestqueue) == nbest
    global nbestlist = Array(Vector{Int}, nbest)
    for i=nbest:-1:1
        nbestlist[i] = dequeue!(nbestqueue)
    end
    return nbestlist
end

using Base.Collections

function topn{T}(a::Array{T},n::Int)
    p = PriorityQueue{Int,T,Base.Order.ForwardOrdering}()
    plen::Int = 0
    pmin::T = typemax(T)
    @inbounds for i=1:length(a)
        if plen < n
            p[i] = a[i]
            plen += 1
            pmin = peek(p).second
        elseif a[i] > pmin
            dequeue!(p)
            p[i] = a[i]
            pmin = peek(p).second
        end
    end
    k = Array(Int,n)
    @inbounds for i=n:-1:1
        k[i] = dequeue!(p)
    end
    return k
end

function clear!(m::S2S)
    for net in (m.encoder, m.decoder)
        for i in 1:nops(net)
            if !isa(net.op[i],Par) && !isa(net.op[i],Arr)
                net.out[i] = net.out0[i] = net.dif[i] = net.dif0[i] = net.tmp[i] = nothing
            end
        end
    end
    return m
end

function s2s_loop(m::S2S, data, loss; gcheck=false, o...)
    s2s_lossreport()
    decoding = false
    reset!(m; o...)
    for (x,ygold,mask) in data
        if decoding && ygold == nothing # the next sentence started
            gcheck && break
            s2s_eos(m, data, loss; gcheck=gcheck, o...)
            reset!(m; o...)
            decoding = false
        end
        if !decoding && ygold != nothing # source ended, target sequence started
            s2s_copyforw!(m)
            decoding = true
        end
        if decoding && ygold != nothing # keep decoding target
            s2s_decode(m, x, ygold, mask, loss; o...)
        end
        if !decoding && ygold == nothing # keep encoding source
            s2s_encode(m, x; o...)
        end
    end
    s2s_eos(m, data, loss; gcheck=gcheck, o...)
end

function s2s_encode(m::S2S, x; trn=false, o...)
    forw(m.encoder, x; trn=trn, seq=true, o...)
end

function s2s_decode(m::S2S, x, ygold, mask, loss; trn=false, ystack=nothing, losscnt=nothing, o...)
    ypred = forw(m.decoder, x; trn=trn, seq=true, o...)
    ystack != nothing  && push!(ystack, (copy(ygold),copy(mask))) # TODO: get rid of alloc
    losscnt != nothing && s2s_loss(m, ypred, ygold, mask, loss; losscnt=losscnt, o...)
end

function s2s_loss(m::S2S, ypred, ygold, mask, loss; losscnt=nothing, lossreport=0, o...)
    (yrows, ycols) = size2(ygold)
    nwords = (mask == nothing ? ycols : sum(mask))
    losscnt[1] += loss(ypred,ygold;mask=mask) # loss divides total loss by minibatch size ycols.  at the end the total loss will be equal to
    losscnt[2] += nwords/ycols                # losscnt[1]*ycols.  losscnt[1]/losscnt[2] will equal totalloss/totalwords.
    # we could scale losscnt with ycols so losscnt[1] is total loss and losscnt[2] is total words, but I think that breaks gradcheck since the scaled versions are what gets used for parameter updates in order to prevent batch size from effecting step size.
    lossreport > 0 && s2s_lossreport(losscnt,ycols,lossreport)
end

s2s_time0 = s2s_time1 = s2s_inst = 0
s2s_print(a...)=(for x in a; @printf("%.2f ",x); end; println(); flush(STDOUT))

function s2s_lossreport()
    global s2s_time0, s2s_time1, s2s_inst
    s2s_inst = 0
    s2s_time0 = s2s_time1 = time_ns()
    println("time inst speed perp")
end

function s2s_lossreport(losscnt,batchsize,lossreport)
    global s2s_time0, s2s_time1, s2s_inst
    s2s_time0 == 0 && s2s_lossreport()
    losscnt == nothing && return
    losscnt[2]*batchsize < lossreport && return
    curr_time = time_ns()
    batch_time = Int(curr_time - s2s_time1)/10^9
    total_time = Int(curr_time - s2s_time0)/10^9
    s2s_time1 = curr_time
    batch_inst = losscnt[2]*batchsize
    batch_loss = losscnt[1]*batchsize
    s2s_inst += batch_inst
    s2s_print(total_time, s2s_inst, batch_inst/batch_time, exp(batch_loss/batch_inst))
    losscnt[1] = losscnt[2] = 0
end

function s2s_eos(m::S2S, data, loss; trn=false, gcheck=false, ystack=nothing, maxnorm=nothing, gclip=0, o...)
    if trn
        s2s_bptt(m, ystack, loss; o...)
        g = (gclip > 0 || maxnorm!=nothing ? gnorm(m) : 0)
        if !gcheck
            gclip=(g > gclip > 0 ? gclip/g : 0)
            update!(m; gclip=gclip, o...)
        end
    end
    if maxnorm != nothing
        w=wnorm(m)
        w > maxnorm[1] && (maxnorm[1]=w)
        g > maxnorm[2] && (maxnorm[2]=g)
    end
end

function s2s_bptt(m::S2S, ystack, loss; o...)
    while !isempty(ystack)
        (ygold,mask) = pop!(ystack)
        back(m.decoder, ygold, loss; seq=true, mask=mask, o...)
    end
    @assert m.decoder.sp == 0
    s2s_copyback!(m)
    while m.encoder.sp > 0
        back(m.encoder; seq=true, o...)
    end
end

function s2s_copyforw!(m::S2S)
    for n=1:nops(m.encoder)
        if forwref(m.encoder, n)
            # m.decoder.out0[n] == nothing && (m.decoder.out0[n] = similar(m.encoder.out[n]))
            # m.decoder.out[n] = copysync!(m.decoder.out0[n], m.encoder.out[n])
            m.decoder.out[n] = m.decoder.out0[n] = m.encoder.out[n]
        end
    end
end

function s2s_copyforw!(n1::Net, n2::Net, cols)
    ncols = length(cols)
    for n=1:nops(n1)
        if forwref(n1, n)
            o1 = n1.out[n]
            @assert o1 != nothing
            nrows = size(o1,1)
            o2 = n2.tmp[n]
            if o2 == nothing || size(o2) != (nrows, ncols)
                o2 = similar(o1, (nrows, ncols))
            end
            for i=1:length(cols)
                copysync!(o2, 1+(i-1)*nrows, o1, 1+(cols[i]-1)*nrows, nrows)
            end
            n2.tmp[n] = n2.out0[n]
            n2.out[n] = n2.out0[n] = o2
        end
    end
end

function s2s_copyback!(m::S2S)
    for n=1:nops(m.encoder)
        if forwref(m.encoder, n)
            # m.encoder.dif0[n] == nothing && (m.encoder.dif0[n] = similar(m.decoder.dif[n]))
            # m.encoder.dif[n] = copysync!(m.encoder.dif0[n], m.decoder.dif[n])
            m.encoder.dif[n] = m.encoder.dif0[n] = m.decoder.dif[n]
        end
    end
end

function Base.isequal(a::S2S,b::S2S)
    typeof(a)==typeof(b) || return false
    for n in fieldnames(a)
        if isdefined(a,n) && isdefined(b,n)
            isequal(a.(n), b.(n)) || return false
        elseif isdefined(a,n) || isdefined(b,n)
            return false
        end
    end
    return true
end

# FAQ:
#
# Q: How do we handle the transition?  Is eos fed to the encoder or the decoder?
#   It looks like the decoder from the picture.  We handle switch using the state variable.
# Q: Do we feed one best, distribution, or gold for word[t-1] to the decoder?
#   Gold during training, one-best output during testing?  (could also try the actual softmax output)
# Q: How do we handle the sequence ending?  How is the training signal passed back?
#   gold outputs, bptt first through decoder then encoder, triggered by output eos (or data=nothing when state=decode)
# Q: input reversal?
#   handled by the data generator
# Q: batching and padding of inputs of different length?
#   handled by the data generator
# Q: data format? <eos> ran dog The <eos> el perror corrio <eos> sentence next The <eos> its spanish version <eos> ...
#   so data does not need to have x/y pairs, just a y sequence for training (and testing?).
#   we could use nothings to signal eos, but we do need an actual encoding for <eos>
#   so don't replace <eos>, keep it, just insert nothings at the end of each sentence.
#   except in the very beginning, a state variable keeps track of encoding vs decoding
# Q: how about the data generator gives us x/y pairs:
#   during encoding we get x/nothing as desired input/output.
#   during decoding we get x[t-1]/x[t] as desired input/output.
#   this is more in line with what we did with rnnlm,
#   also in line with models not caring about how data is formatted, only desired input/output.
#   in this design we can distinguish encoding/decoding by looking at output
# Q: test time:
#   encode proceeds normally
#   when we reach nothing, we switch to decoding
#   when we reach nothing again, we perform bptt (need stack for gold output, similar to rnnlm)
#   the first token to decoder gets fed from the gold (presumably <eos>)
#   the subsequent input tokens come from decoder's own output
#   second nothing resets and switches to encoding, no bptt
# Q: beam search?
#   don't need this for training or test perplexity
#   don't need this for one best greedy output
#   it will make minibatching more difficult
#   this is a problem for predict (TODO)
# Q: gradient scaling?
#   this becomes a problem once we start doing minibatching and padding short sequences.  I started
#   dividing loss and its gradients with minibatch size so the choice of minibatch does not effect
#   the step sizes and require a different learning rate.  In s2s if we have a sequence minibatch of
#   size N, which consists of T token minibatches each of size N.  Some sequences in the minibatch
#   are shorter so some of the token minibatches contain padding.  Say a particular token minibatch
#   contains n<N real words, the rest is padding.  Do we normalize using n or N?  It turns out in
#   this case the correct answer is N because we do not want the word in this minibatch to have
#   higher weight compared to other words in the sequence!
# Q: do we zero out padded columns?
#   We should set padding columns to zero in the input.  Otherwise when
#   the next time step has a legitimate word in that column, the hidden
#   state column will not be zero.

### DEAD CODE:
    # if losscnt != nothing
    #     l = loss(ypred, ygold)  # returns total loss / num columns ignoring padding
    #     ycols = size(ygold,2)
    #     nzcols = 0
    #     for j=1:ycols; zerocolumn(ygold,j) || (nzcols += 1); end
    #     losscnt[1] += l # this has been scaled by (1/ycols) by the loss function
    #     losscnt[2] += (nzcols/ycols)    # this should work as long as loss for padded columns is 0 as it should be
    # end


text:
## Packages
include("../src/StencilBuilder.jl")

using BenchmarkTools
using LoopVectorization
using CUDA
using .StencilBuilder

## CPU dynamics kernel.
function dynamics_w_kernel_cpu!(
    wt, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    alpha = 9.81/300;

    @fast3d begin
        @fd (wt, u, v, w, s) wt += (
            - gradx(interpz(u) * interpx(w)) + visc * (gradx(gradx(w)))
            - grady(interpz(v) * interpy(w)) + visc * (grady(grady(w)))
            - gradz(interpz(w) * interpz(w)) + visc * (gradz(gradz(w)))
            + alpha*interpz(s) )
    end
end

## GPU dynamics kernel.
function dynamics_w_kernel_gpu!(
    wt, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    i = (blockIdx.().x - 1) * blockDim().x + threadIdx().x + is - 1
    j = (blockIdx.().y - 1) * blockDim().y + threadIdx().y + js - 1
    k = blockIdx.().z + ks - 1

    alpha = 9.81/300;

    if i <= ie && j <= je && k <= ke
        @inbounds @fd (wt, u, v, w, s) wt += (
            - gradx(interpz(u) * interpx(w)) + visc * (gradx(gradx(w)))
            - grady(interpz(v) * interpy(w)) + visc * (grady(grady(w)))
            - gradz(interpz(w) * interpz(w)) + visc * (gradz(gradz(w)))
            + alpha*interpz(s) )
    end

    return
end

## User input
itot = 384; jtot = 384; ktot = 384;
igc = 1; jgc = 1; kgc = 1;

icells = itot + 2igc; jcells = jtot + 2jgc; kcells = ktot + 2kgc
is = igc+1; ie = igc+itot; js = jgc+1; je = jgc+jtot; ks = kgc+1; ke = kgc+ktot

wt = rand(icells, jcells, kcells)
u = rand(icells, jcells, kcells)
v = rand(icells, jcells, kcells)
w = rand(icells, jcells, kcells)
s = rand(icells, jcells, kcells)
dxi = rand()
dyi = rand()
dzi = rand(kcells)
dzhi = rand(kcells)
visc = 1.

## Copy data to GPU.
wt_gpu = CuArray(wt)
u_gpu = CuArray(u)
v_gpu = CuArray(v)
w_gpu = CuArray(w)
s_gpu = CuArray(s)
dzi_gpu = CuArray(dzi)
dzhi_gpu = CuArray(dzhi)

## Benchmark functions
function benchmark_cpu(
    wt, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    dynamics_w_kernel_cpu!(
        wt, u, v, w, s,
        visc,
        dxi, dyi, dzi, dzhi,
        is, ie, js, je, ks, ke)
end

function benchmark_gpu(
    wt_gpu, u_gpu, v_gpu, w_gpu, s_gpu,
    visc,
    dxi, dyi, dzi_gpu, dzhi_gpu,
    is, ie, js, je, ks, ke)

    blocks_size = (32, 32, 1)
    blocks_num = (
        ceil(Int, (ie-is+1) / blocks_size[1]),
        ceil(Int, (je-js+1) / blocks_size[2]),
        (ke-ks+1) )

    CUDA.@sync begin
        @cuda threads=blocks_size blocks=blocks_num dynamics_w_kernel_gpu!(
            wt_gpu, u_gpu, v_gpu, w_gpu, s_gpu,
            visc,
            dxi, dyi, dzi_gpu, dzhi_gpu,
            is, ie, js, je, ks, ke)
    end
end

## Correctness test
benchmark_cpu(
    wt, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

benchmark_gpu(
    wt_gpu, u_gpu, v_gpu, w_gpu, s_gpu,
    visc,
    dxi, dyi, dzi_gpu, dzhi_gpu,
    is, ie, js, je, ks, ke)

wt_gpu_cpu = Array(wt_gpu)
println("Are wt and wt_gpu equal? ", isequal(wt, wt_gpu_cpu), " ", isapprox(wt, wt_gpu_cpu))

## Benchmarks
@btime benchmark_cpu(
    wt, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

@btime benchmark_gpu(
    wt_gpu, u_gpu, v_gpu, w_gpu, s_gpu,
    visc,
    dxi, dyi, dzi_gpu, dzhi_gpu,
    is, ie, js, je, ks, ke)


text:
using ShuHua
using Test


text:
const ALL_FILES = [FILE_A, FILE_B, FILE_C, FILE_D, FILE_E, FILE_F, FILE_G, FILE_H]
const ALL_RANKS = [RANK_1, RANK_2, RANK_3, RANK_4, RANK_5, RANK_6, RANK_7, RANK_8]

for f in ALL_FILES
    @test f == filefromchar(tochar(f))
end

for r in ALL_RANKS
    @test r == rankfromchar(tochar(r))
end

const ALL_SQUARES = [Square(f, r) for f in ALL_FILES, r in ALL_RANKS]

for s in ALL_SQUARES
    @test s == Square(file(s), rank(s))
    @test s == squarefromstring(tostring(s))
end

@test SQ_D5 - SQ_D6 == DELTA_S
@test SQ_D5 - SQ_D4 == DELTA_N
@test SQ_D5 - SQ_E5 == DELTA_W
@test SQ_D5 - SQ_C5 == DELTA_E
@test SQ_D5 - SQ_E4 == DELTA_NW
@test SQ_D5 - SQ_C4 == DELTA_NE
@test SQ_D5 - SQ_E6 == DELTA_SW
@test SQ_D5 - SQ_C6 == DELTA_SE


text:
function theme_ggplot2()
    Theme(
        Axis = (
            backgroundcolor = :gray92,
            xgridcolor = :white,
            ygridcolor = :white,
            xminorgridcolor = (:white, 0.5),
            xminorgridvisible = true,
            yminorgridcolor = (:white, 0.5),
            yminorgridvisible = true,
            leftspinevisible = false,
            rightspinevisible = false,
            bottomspinevisible = false,
            topspinevisible = false,
            xminorticksvisible = false,
            yminorticksvisible = false,
            xtickcolor = :gray21,
            ytickcolor = :gray21,
            spinewidth = 0,
            xticklabelcolor = :gray31,
            yticklabelcolor = :gray31,
            xticksize = 4,
            yticksize = 4,
            xtickwidth = 1.5,
            ytickwidth = 1.5,
            xgridwidth = 1.5,
            ygridwidth = 1.5,
            xlabelpadding = 3,
            ylabelpadding = 3
        ),
        Legend = (
            framevisible = false,
            padding = (0, 0, 0, 0),
            patchcolor = :gray92,
        ),
        Axis3 = (
            xgridcolor = :white,
            ygridcolor = :white,
            zgridcolor = :white,
            xspinesvisible = false,
            yspinesvisible = false,
            zspinesvisible = false,
            yzpanelcolor = :gray92,
            xzpanelcolor = :gray92,
            xypanelcolor = :gray92,
            xticklabelcolor = :gray31,
            yticklabelcolor = :gray31,
            zticklabelcolor = :gray31,
        ),
        Colorbar = (
            tickcolor = :white,
            tickalign = 1,
            ticklabelcolor = :gray31,
            spinewidth = 0,
            ticklabelpad = 5,
        )
    )
end


text:
abstract type AbstractDictionary{T, K} end

### Interface

# Insert a (key, record) pair.
insert!(d::AbstractDictionary, key, item) = nothing

# Remove and return one out of possible many records given a key.
remove!(d::AbstractDictionary, key) = nothing

# Remove and return an arbitrary record.
removeany!(d::AbstractDictionary) = nothing

# Return one out of possible many records given a key.
find(d::AbstractDictionary, key) = nothing

# Return the number of records in the dictionary.
length(d::AbstractDictionary) = nothing

# Check if dictionary is empty.
isempty(d::AbstractDictionary) = nothing

### Generic

#eltype(d::AbstractDictionary{T, K}) where {T, K} = T