text:
#
# Dense storage
#
using LinearAlgebra: BlasFloat

struct Dense{ElT,VecT<:AbstractVector} <: TensorStorage{ElT}
  data::VecT
  function Dense{ElT,VecT}(data::AbstractVector) where {ElT,VecT<:AbstractVector{ElT}}
    return new{ElT,VecT}(data)
  end
end

function Dense(data::VecT) where {VecT<:AbstractVector{ElT}} where {ElT}
  return Dense{ElT,VecT}(data)
end

function Dense{ElR}(data::AbstractVector{ElT}) where {ElR,ElT}
  return ElT == ElR ? Dense(data) : Dense(ElR.(data))
end

# Construct from a set of indices
function Dense{ElT,VecT}(inds) where {ElT,VecT<:AbstractVector{ElT}}
  return Dense(VecT(dim(inds)))
end

Dense{ElT}(dim::Integer) where {ElT} = Dense(zeros(ElT, dim))

Dense{ElT}(::UndefInitializer, dim::Integer) where {ElT} = Dense(Vector{ElT}(undef, dim))

Dense(::Type{ElT}, dim::Integer) where {ElT} = Dense{ElT}(dim)

Dense(x::ElT, dim::Integer) where {ElT<:Number} = Dense(fill(x, dim))

Dense(dim::Integer) = Dense(Float64, dim)

Dense(::Type{ElT}, ::UndefInitializer, dim::Integer) where {ElT} = Dense{ElT}(undef, dim)

Dense(::UndefInitializer, dim::Integer) = Dense(Float64, undef, dim)

Dense{ElT}() where {ElT} = Dense(ElT[])
Dense(::Type{ElT}) where {ElT} = Dense{ElT}()

setdata(D::Dense, ndata) = Dense(ndata)

#
# Random constructors
#

function randn(::Type{StoreT}, dim::Integer) where {StoreT<:Dense}
  return Dense(randn(eltype(StoreT), dim))
end

copy(D::Dense) = Dense(copy(data(D)))

Base.real(::Type{Dense{ElT,Vector{ElT}}}) where {ElT} = Dense{real(ElT),Vector{real(ElT)}}

function complex(::Type{Dense{ElT,Vector{ElT}}}) where {ElT}
  return Dense{complex(ElT),Vector{complex(ElT)}}
end

similar(D::Dense) = Dense(similar(data(D)))

similar(D::Dense, length::Int) = Dense(similar(data(D), length))

function similar(
  ::Type{<:Dense{ElT,VecT}}, length::Int
) where {ElT,VecT<:AbstractVector{ElT}}
  return Dense(similar(Vector{ElT}, length))
end

function similartype(::Type{StoreT}, ::Type{ElT}) where {StoreT<:Dense,ElT}
  return Dense{ElT,similartype(datatype(StoreT), ElT)}
end

# TODO: make these more general, move to tensorstorage.jl
datatype(::Type{<:Dense{<:Any,DataT}}) where {DataT} = DataT

function similar(::Type{StorageT}, ::Type{ElT}, length::Int) where {StorageT<:Dense,ElT}
  return Dense(similar(datatype(StorageT), ElT, length))
end

similar(D::Dense, ::Type{T}) where {T<:Number} = Dense(similar(data(D), T))

zeros(DenseT::Type{<:Dense{ElT}}, inds) where {ElT} = zeros(DenseT, dim(inds))

# TODO: should this do something different for SubArray?
zeros(::Type{<:Dense{ElT}}, dim::Int) where {ElT} = Dense{ElT}(dim)

function promote_rule(
  ::Type{<:Dense{ElT1,VecT1}}, ::Type{<:Dense{ElT2,VecT2}}
) where {ElT1,VecT1,ElT2,VecT2}
  ElR = promote_type(ElT1, ElT2)
  VecR = promote_type(VecT1, VecT2)
  return Dense{ElR,VecR}
end

# This is to get around the issue in Julia that:
# promote_type(Vector{ComplexF32},Vector{Float64}) == Vector{T} where T
function promote_rule(
  ::Type{<:Dense{ElT1,Vector{ElT1}}}, ::Type{<:Dense{ElT2,Vector{ElT2}}}
) where {ElT1,ElT2}
  ElR = promote_type(ElT1, ElT2)
  VecR = Vector{ElR}
  return Dense{ElR,VecR}
end

# This is for type promotion for Scalar*Dense
function promote_rule(
  ::Type{<:Dense{ElT1,Vector{ElT1}}}, ::Type{ElT2}
) where {ElT1,ElT2<:Number}
  ElR = promote_type(ElT1, ElT2)
  VecR = Vector{ElR}
  return Dense{ElR,VecR}
end

function convert(::Type{<:Dense{ElR,VecR}}, D::Dense) where {ElR,VecR}
  return Dense(convert(VecR, data(D)))
end

#
# DenseTensor (Tensor using Dense storage)
#

const DenseTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Dense}

DenseTensor(::Type{ElT}, inds) where {ElT} = tensor(Dense(ElT, dim(inds)), inds)

# Special convenience function for Int
# dimensions
DenseTensor(::Type{ElT}, inds::Int...) where {ElT} = DenseTensor(ElT, inds)

DenseTensor(inds) = tensor(Dense(dim(inds)), inds)

DenseTensor(inds::Int...) = DenseTensor(inds)

function DenseTensor(::Type{ElT}, ::UndefInitializer, inds) where {ElT}
  return tensor(Dense(ElT, undef, dim(inds)), inds)
end

function DenseTensor(::Type{ElT}, ::UndefInitializer, inds::Int...) where {ElT}
  return DenseTensor(ElT, undef, inds)
end

DenseTensor(::UndefInitializer, inds) = tensor(Dense(undef, dim(inds)), inds)

DenseTensor(::UndefInitializer, inds::Int...) = DenseTensor(undef, inds)

# For convenience, direct Tensor constructors default to Dense
Tensor(::Type{ElT}, inds...) where {ElT} = DenseTensor(ElT, inds...)

Tensor(inds...) = Tensor(Float64, inds...)

function Tensor(::Type{ElT}, ::UndefInitializer, inds...) where {ElT}
  return DenseTensor(ElT, undef, inds...)
end

Tensor(::UndefInitializer, inds...) = DenseTensor(undef, inds...)

Tensor(A::Array{<:Number,N}, inds::Dims{N}) where {N} = tensor(Dense(vec(A)), inds)

#
# Random constructors
#

function randomDenseTensor(::Type{ElT}, inds) where {ElT}
  return tensor(randn(Dense{ElT}, dim(inds)), inds)
end

function randomDenseTensor(::Type{ElT}, inds::Int...) where {ElT}
  return randomDenseTensor(ElT, inds)
end

randomDenseTensor(inds) = randomDenseTensor(Float64, inds)

randomDenseTensor(inds::Int...) = randomDenseTensor(Float64, inds)

function randomTensor(::Type{ElT}, inds) where {ElT}
  return randomDenseTensor(ElT, inds)
end

function randomTensor(::Type{ElT}, inds::Int...) where {ElT}
  return randomDenseTensor(ElT, inds...)
end

randomTensor(inds) = randomDenseTensor(Float64, inds)

randomTensor(inds::Int...) = randomDenseTensor(Float64, inds)

# Basic functionality for AbstractArray interface
IndexStyle(::Type{<:DenseTensor}) = IndexLinear()

# Override CartesianIndices iteration to iterate
# linearly through the Dense storage (faster)
iterate(T::DenseTensor, args...) = iterate(storage(T), args...)

function _zeros(TensorT::Type{<:DenseTensor}, inds)
  return tensor(zeros(storagetype(TensorT), dim(inds)), inds)
end

function zeros(TensorT::Type{<:DenseTensor}, inds)
  return _zeros(TensorT, inds)
end

# To fix method ambiguity with zeros(::Type, ::Tuple)
function zeros(TensorT::Type{<:DenseTensor}, inds::Dims)
  return _zeros(TensorT, inds)
end

function zeros(TensorT::Type{<:DenseTensor}, inds::Tuple{})
  return _zeros(TensorT, inds)
end

# To fix method ambiguity with similar(::AbstractArray,::Type)
function similar(T::DenseTensor, ::Type{ElT}) where {ElT}
  return tensor(similar(storage(T), ElT), inds(T))
end

_similar(T::DenseTensor, inds) = similar(typeof(T), inds)

similar(T::DenseTensor, inds) = _similar(T, inds)

# To fix method ambiguity with similar(::AbstractArray,::Tuple)
similar(T::DenseTensor, inds::Dims) = _similar(T, inds)

#
# Single index
#

@propagate_inbounds function getindex(T::DenseTensor{<:Number}, I::Integer...)
  Base.@_inline_meta
  return getindex(data(T), Base._sub2ind(T, I...))
end

@propagate_inbounds function getindex(T::DenseTensor{<:Number}, I::CartesianIndex)
  Base.@_inline_meta
  return getindex(T, I.I...)
end

@propagate_inbounds function setindex!(
  T::DenseTensor{<:Number}, x::Number, I::Vararg{Integer}
)
  Base.@_inline_meta
  setindex!(data(T), x, Base._sub2ind(T, I...))
  return T
end

@propagate_inbounds function setindex!(
  T::DenseTensor{<:Number}, x::Number, I::CartesianIndex
)
  Base.@_inline_meta
  setindex!(T, x, I.I...)
  return T
end

#
# Linear indexing
#

@propagate_inbounds @inline getindex(T::DenseTensor, i::Integer) = storage(T)[i]

@propagate_inbounds @inline function setindex!(T::DenseTensor, v, i::Integer)
  return (storage(T)[i] = v; T)
end

#
# Slicing
# TODO: this doesn't allow colon right now
# Create a DenseView that stores a Dense and an offset?
#

## @propagate_inbounds function _getindex(
##   T::DenseTensor{ElT,N}, I::CartesianIndices{N}
## ) where {ElT,N}
##   storeR = Dense(vec(@view array(T)[I]))
##   indsR = Tuple(I[end] - I[1] + CartesianIndex(ntuple(_ -> 1, Val(N))))
##   return tensor(storeR, indsR)
## end
##
## @propagate_inbounds function getindex(T::DenseTensor{ElT,N}, I...) where {ElT,N}
##   return _getindex(T, CartesianIndices(I))
## end

@propagate_inbounds function getindex(T::DenseTensor, I...)
  AI = @view array(T)[I...]
  storeR = Dense(vec(AI))
  indsR = size(AI)
  return tensor(storeR, indsR)
end

# Reshape a DenseTensor using the specified dimensions
# This returns a view into the same Tensor data
function reshape(T::DenseTensor, dims)
  dim(T) == dim(dims) || error("Total new dimension must be the same as the old dimension")
  return tensor(storage(T), dims)
end

# This version fixes method ambiguity with AbstractArray reshape
function reshape(T::DenseTensor, dims::Dims)
  dim(T) == dim(dims) || error("Total new dimension must be the same as the old dimension")
  return tensor(storage(T), dims)
end

function reshape(T::DenseTensor, dims::Int...)
  return tensor(storage(T), tuple(dims...))
end

convert(::Type{Array}, T::DenseTensor) = reshape(data(storage(T)), dims(inds(T)))

# Create an Array that is a view of the Dense Tensor
# Useful for using Base Array functions
array(T::DenseTensor) = convert(Array, T)

function Array{ElT,N}(T::DenseTensor{ElT,N}) where {ElT,N}
  return copy(array(T))
end

function Array(T::DenseTensor{ElT,N}) where {ElT,N}
  return Array{ElT,N}(T)
end

# If the storage data are regular Vectors, use Base.copyto!
function copyto!(
  R::Tensor{<:Number,N,<:Dense{<:Number,<:Vector}},
  T::Tensor{<:Number,N,<:Dense{<:Number,<:Vector}},
) where {N}
  RA = array(R)
  TA = array(T)
  RA .= TA
  return R
end

# If they are something more complicated like views, use Strided copyto!
function copyto!(R::DenseTensor{<:Number,N}, T::DenseTensor{<:Number,N}) where {N}
  RA = array(R)
  TA = array(T)
  @strided RA .= TA
  return R
end

# TODO: call permutedims!(R,T,perm,(r,t)->t)?
function permutedims!(
  R::DenseTensor{<:Number,N}, T::DenseTensor{<:Number,N}, perm::NTuple{N,Int}
) where {N}
  RA = array(R)
  TA = array(T)
  @strided RA .= permutedims(TA, perm)
  return R
end

function apply!(R::DenseTensor, T::DenseTensor, f::Function=(r, t) -> t)
  RA = array(R)
  TA = array(T)
  @strided RA .= f.(RA, TA)
  return R
end

# Version that may overwrite the result or promote
# and return the result
function permutedims!!(R::DenseTensor, T::DenseTensor, perm::Tuple, f::Function=(r, t) -> t)
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  #RA = array(R)
  #TA = array(T)
  RA = ReshapedArray(data(RR), dims(RR), ())
  TA = ReshapedArray(data(T), dims(T), ())
  if !is_trivial_permutation(perm)
    @strided RA .= f.(RA, permutedims(TA, perm))
  else
    # TODO: specialize for specific functions
    RA .= f.(RA, TA)
  end
  return RR
end

# TODO: move to tensor.jl?
function permutedims(T::Tensor{<:Number,N}, perm::NTuple{N,Int}) where {N}
  Tp = similar(T, permute(inds(T), perm))
  Tp = permutedims!!(Tp, T, perm)
  return Tp
end

function permutedims(::Tensor, ::Tuple{Vararg{Int}})
  return error("Permutation size must match tensor order")
end

# TODO: move to tensor.jl?
function (x::Number * T::Tensor)
  return tensor(x * storage(T), inds(T))
end
(T::Tensor * x::Number) = x * T

function (T::Tensor / x::Number)
  return tensor(storage(T) / x, inds(T))
end

function permutedims!(
  R::DenseTensor{<:Number,N}, T::DenseTensor{<:Number,N}, perm, f::Function
) where {N}
  if nnz(R) == 1 && nnz(T) == 1
    R[1] = f(R[1], T[1])
    return R
  end
  RA = array(R)
  TA = array(T)
  @strided RA .= f.(RA, permutedims(TA, perm))
  return R
end

function outer!(
  R::DenseTensor{ElR}, T1::DenseTensor{ElT1}, T2::DenseTensor{ElT2}
) where {ElR,ElT1,ElT2}
  if ElT1 != ElT2
    # TODO: use promote instead
    # T1,T2 = promote(T1,T2)

    ElT1T2 = promote_type(ElT1, ElT2)
    if ElT1 != ElT1T2
      # TODO: get this working
      # T1 = ElR.(T1)
      T1 = one(ElT1T2) * T1
    end
    if ElT2 != ElT1T2
      # TODO: get this working
      # T2 = ElR.(T2)
      T2 = one(ElT1T2) * T2
    end
  end

  v1 = data(T1)
  v2 = data(T2)
  RM = reshape(R, length(v1), length(v2))
  #RM .= v1 .* transpose(v2)
  #mul!(RM, v1, transpose(v2))
  _gemm!('N', 'T', one(ElR), v1, v2, zero(ElR), RM)
  return R
end

export backend_auto, backend_blas, backend_generic

@eval struct GemmBackend{T}
  (f::Type{<:GemmBackend})() = $(Expr(:new, :f))
end
GemmBackend(s) = GemmBackend{Symbol(s)}()
macro GemmBackend_str(s)
  return :(GemmBackend{$(Expr(:quote, Symbol(s)))})
end

const gemm_backend = Ref(:Auto)
function backend_auto()
  return gemm_backend[] = :Auto
end
function backend_blas()
  return gemm_backend[] = :BLAS
end
function backend_generic()
  return gemm_backend[] = :Generic
end

@inline function auto_select_backend(
  ::Type{<:StridedVecOrMat{<:BlasFloat}},
  ::Type{<:StridedVecOrMat{<:BlasFloat}},
  ::Type{<:StridedVecOrMat{<:BlasFloat}},
)
  return GemmBackend(:BLAS)
end

@inline function auto_select_backend(
  ::Type{<:AbstractVecOrMat}, ::Type{<:AbstractVecOrMat}, ::Type{<:AbstractVecOrMat}
)
  return GemmBackend(:Generic)
end

function _gemm!(
  tA, tB, alpha, A::TA, B::TB, beta, C::TC
) where {TA<:AbstractVecOrMat,TB<:AbstractVecOrMat,TC<:AbstractVecOrMat}
  if gemm_backend[] == :Auto
    _gemm!(auto_select_backend(TA, TB, TC), tA, tB, alpha, A, B, beta, C)
  else
    _gemm!(GemmBackend(gemm_backend[]), tA, tB, alpha, A, B, beta, C)
  end
end

# BLAS matmul
function _gemm!(
  ::GemmBackend{:BLAS},
  tA,
  tB,
  alpha,
  A::AbstractVecOrMat,
  B::AbstractVecOrMat,
  beta,
  C::AbstractVecOrMat,
)
  #@timeit_debug timer "BLAS.gemm!" begin
  return BLAS.gemm!(tA, tB, alpha, A, B, beta, C)
  #end # @timeit
end

# generic matmul
function _gemm!(
  ::GemmBackend{:Generic},
  tA,
  tB,
  alpha::AT,
  A::AbstractVecOrMat,
  B::AbstractVecOrMat,
  beta::BT,
  C::AbstractVecOrMat,
) where {AT,BT}
  mul!(C, tA == 'T' ? transpose(A) : A, tB == 'T' ? transpose(B) : B, alpha, beta)
  return C
end

function outer!!(R::Tensor, T1::Tensor, T2::Tensor)
  outer!(R, T1, T2)
  return R
end

# TODO: call outer!!, make this generic
function outer(T1::DenseTensor{ElT1}, T2::DenseTensor{ElT2}) where {ElT1,ElT2}
  array_outer = vec(array(T1)) * transpose(vec(array(T2)))
  inds_outer = unioninds(inds(T1), inds(T2))
  return tensor(Dense{promote_type(ElT1, ElT2)}(vec(array_outer)), inds_outer)
end
const ⊗ = outer

# TODO: move to tensor.jl?
function contraction_output_type(
  ::Type{TensorT1}, ::Type{TensorT2}, ::Type{IndsR}
) where {TensorT1<:Tensor,TensorT2<:Tensor,IndsR}
  return similartype(promote_type(TensorT1, TensorT2), IndsR)
end

function contraction_output(
  ::TensorT1, ::TensorT2, indsR::IndsR
) where {TensorT1<:DenseTensor,TensorT2<:DenseTensor,IndsR}
  TensorR = contraction_output_type(TensorT1, TensorT2, IndsR)
  return similar(TensorR, indsR)
end

# TODO: move to tensor.jl?
function contraction_output(T1::Tensor, labelsT1, T2::Tensor, labelsT2, labelsR)
  indsR = contract_inds(inds(T1), labelsT1, inds(T2), labelsT2, labelsR)
  R = contraction_output(T1, T2, indsR)
  return R
end

# Version where output labels aren't supplied
function contract(T1::Tensor, labelsT1, T2::Tensor, labelsT2)
  labelsR = contract_labels(labelsT1, labelsT2)
  return contract(T1, labelsT1, T2, labelsT2, labelsR)
end

# TODO: move to tensor.jl?
function contract(T1::Tensor, labelsT1, T2::Tensor, labelsT2, labelsR)
  #@timeit_debug timer "dense contract" begin
  # TODO: put the contract_inds logic into contraction_output,
  # call like R = contraction_ouput(T1,labelsT1,T2,labelsT2)
  #indsR = contract_inds(inds(T1),labelsT1,inds(T2),labelsT2,labelsR)
  R = contraction_output(T1, labelsT1, T2, labelsT2, labelsR)
  # contract!! version here since the output R may not
  # be mutable (like UniformDiag)
  R = contract!!(R, labelsR, T1, labelsT1, T2, labelsT2)
  return R
  #end @timeit
end

# Move to tensor.jl? Is this generic for all storage types?
function contract!!(
  R::Tensor, labelsR, T1::Tensor, labelsT1, T2::Tensor, labelsT2, α::Number=1, β::Number=0
)
  NR = ndims(R)
  N1 = ndims(T1)
  N2 = ndims(T2)
  if (N1 ≠ 0) && (N2 ≠ 0) && (N1 + N2 == NR)
    # Outer product
    (α ≠ 1 || β ≠ 0) && error(
      "contract!! not yet implemented for outer product tensor contraction with non-trivial α and β",
    )
    # TODO: permute T1 and T2 appropriately first (can be more efficient
    # then permuting the result of T1⊗T2)
    # TODO: implement the in-place version directly
    R = outer!!(R, T1, T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR, labelsRp)
    if !is_trivial_permutation(perm)
      Rp = reshape(R, (inds(T1)..., inds(T2)...))
      R = permutedims!!(R, copy(Rp), perm)
    end
  else
    if α ≠ 1 || β ≠ 0
      R = _contract!!(R, labelsR, T1, labelsT1, T2, labelsT2, α, β)
    else
      R = _contract!!(R, labelsR, T1, labelsT1, T2, labelsT2)
    end
  end
  return R
end

# Move to tensor.jl? Overload this function
# for immutable storage types
function _contract!!(
  R::Tensor, labelsR, T1::Tensor, labelsT1, T2::Tensor, labelsT2, α::Number=1, β::Number=0
)
  if α ≠ 1 || β ≠ 0
    contract!(R, labelsR, T1, labelsT1, T2, labelsT2, α, β)
  else
    contract!(R, labelsR, T1, labelsT1, T2, labelsT2)
  end
  return R
end

Strided.StridedView(T::DenseTensor) = StridedView(convert(Array, T))

# Both are scalar-like tensors
function _contract_scalar!(
  R::DenseTensor{ElR},
  labelsR,
  T1::Number,
  labelsT1,
  T2::Number,
  labelsT2,
  α=one(ElR),
  β=zero(ElR),
) where {ElR}
  if iszero(β)
    R[1] = α * T1 * T2
  elseif iszero(α)
    R[1] = β * R[1]
  else
    R[1] = α * T1 * T2 + β * R[1]
  end
  return R
end

# Trivial permutation
# Version where R and T have different element types, so we can't call BLAS
# Instead use Julia's broadcasting (maybe consider Strided in the future)
function _contract_scalar_noperm!(
  R::DenseTensor{ElR}, T::DenseTensor, α, β=zero(ElR)
) where {ElR}
  Rᵈ = data(R)
  Tᵈ = data(T)
  if iszero(β)
    if iszero(α)
      fill!(Rᵈ, 0)
    else
      Rᵈ .= α .* Tᵈ
    end
  elseif isone(β)
    if iszero(α)
      # No-op
      # Rᵈ .= Rᵈ
    else
      Rᵈ .= α .* Tᵈ .+ Rᵈ
    end
  else
    if iszero(α)
      # Rᵈ .= β .* Rᵈ
      BLAS.scal!(length(Rᵈ), β, Rᵈ, 1)
    else
      Rᵈ .= α .* Tᵈ .+ β .* Rᵈ
    end
  end
  return R
end

# Trivial permutation
# Version where R and T are the same element type, so we can
# call BLAS
function _contract_scalar_noperm!(
  R::DenseTensor{ElR}, T::DenseTensor{ElR}, α, β=zero(ElR)
) where {ElR}
  Rᵈ = data(R)
  Tᵈ = data(T)
  if iszero(β)
    if iszero(α)
      fill!(Rᵈ, 0)
    else
      # Rᵈ .= α .* T₂ᵈ
      BLAS.axpby!(α, Tᵈ, β, Rᵈ)
    end
  elseif isone(β)
    if iszero(α)
      # No-op
      # Rᵈ .= Rᵈ
    else
      # Rᵈ .= α .* Tᵈ .+ Rᵈ
      BLAS.axpy!(α, Tᵈ, Rᵈ)
    end
  else
    if iszero(α)
      # Rᵈ .= β .* Rᵈ
      BLAS.scal!(length(Rᵈ), β, Rᵈ, 1)
    else
      # Rᵈ .= α .* Tᵈ .+ β .* Rᵈ
      BLAS.axpby!(α, Tᵈ, β, Rᵈ)
    end
  end
  return R
end

# Non-trivial permutation
function _contract_scalar_perm!(
  Rᵃ::AbstractArray{ElR}, Tᵃ::AbstractArray, perm, α, β=zero(ElR)
) where {ElR}
  if iszero(β)
    if iszero(α)
      fill!(Rᵃ, 0)
    else
      @strided Rᵃ .= α .* permutedims(Tᵃ, perm)
    end
  elseif isone(β)
    if iszero(α)
      # Rᵃ .= Rᵃ
      # No-op
    else
      @strided Rᵃ .= α .* permutedims(Tᵃ, perm) .+ Rᵃ
    end
  else
    if iszero(α)
      # Rᵃ .= β .* Rᵃ
      BLAS.scal!(length(Rᵃ), β, Rᵃ, 1)
    else
      Rᵃ .= α .* permutedims(Tᵃ, perm) .+ β .* Rᵃ
    end
  end
  return Rᵃ
end

function drop_singletons(::Order{N}, labels, dims) where {N}
  labelsᵣ = ntuple(zero, Val(N))
  dimsᵣ = labelsᵣ
  nkeep = 1
  for n in 1:length(dims)
    if dims[n] > 1
      labelsᵣ = @inbounds setindex(labelsᵣ, labels[n], nkeep)
      dimsᵣ = @inbounds setindex(dimsᵣ, dims[n], nkeep)
      nkeep += 1
    end
  end
  return labelsᵣ, dimsᵣ
end

function _contract_scalar_maybe_perm!(
  ::Order{N}, R::DenseTensor{ElR,NR}, labelsR, T::DenseTensor, labelsT, α, β=zero(ElR)
) where {ElR,NR,N}
  labelsRᵣ, dimsRᵣ = drop_singletons(Order(N), labelsR, dims(R))
  labelsTᵣ, dimsTᵣ = drop_singletons(Order(N), labelsT, dims(T))
  perm = getperm(labelsRᵣ, labelsTᵣ)
  if is_trivial_permutation(perm)
    # trivial permutation
    _contract_scalar_noperm!(R, T, α, β)
  else
    # non-trivial permutation
    Rᵣ = ReshapedArray(data(R), dimsRᵣ, ())
    Tᵣ = ReshapedArray(data(T), dimsTᵣ, ())
    _contract_scalar_perm!(Rᵣ, Tᵣ, perm, α, β)
  end
  return R
end

function _contract_scalar_maybe_perm!(
  R::DenseTensor{ElR,NR}, labelsR, T::DenseTensor, labelsT, α, β=zero(ElR)
) where {ElR,NR}
  N = count(≠(1), dims(R))
  _contract_scalar_maybe_perm!(Order(N), R, labelsR, T, labelsT, α, β)
  return R
end

# XXX: handle case of non-trivial permutation
function _contract_scalar_maybe_perm!(
  R::DenseTensor{ElR,NR},
  labelsR,
  T₁::DenseTensor,
  labelsT₁,
  T₂::DenseTensor,
  labelsT₂,
  α=one(ElR),
  β=zero(ElR),
) where {ElR,NR}
  if nnz(T₁) == 1
    _contract_scalar_maybe_perm!(R, labelsR, T₂, labelsT₂, α * T₁[1], β)
  elseif nnz(T₂) == 1
    _contract_scalar_maybe_perm!(R, labelsR, T₁, labelsT₁, α * T₂[1], β)
  else
    error("In _contract_scalar_perm!, one tensor must be a scalar")
  end
  return R
end

# At least one of the tensors is size 1
function _contract_scalar!(
  R::DenseTensor{ElR},
  labelsR,
  T1::DenseTensor,
  labelsT1,
  T2::DenseTensor,
  labelsT2,
  α=one(ElR),
  β=zero(ElR),
) where {ElR}
  if nnz(T1) == nnz(T2) == 1
    _contract_scalar!(R, labelsR, T1[1], labelsT1, T2[1], labelsT2, α, β)
  else
    _contract_scalar_maybe_perm!(R, labelsR, T1, labelsT1, T2, labelsT2, α, β)
  end
  return R
end

function contract!(
  R::DenseTensor{ElR,NR},
  labelsR,
  T1::DenseTensor{ElT1,N1},
  labelsT1,
  T2::DenseTensor{ElT2,N2},
  labelsT2,
  α::Elα=one(ElR),
  β::Elβ=zero(ElR),
) where {Elα,Elβ,ElR,ElT1,ElT2,NR,N1,N2}
  # Special case for scalar tensors
  if nnz(T1) == 1 || nnz(T2) == 1
    _contract_scalar!(R, labelsR, T1, labelsT1, T2, labelsT2, α, β)
    return R
  end

  if using_tblis() && ElR <: LinearAlgebra.BlasReal && (ElR == ElT1 == ElT2 == Elα == Elβ)
    #@timeit_debug timer "TBLIS contract!" begin
    contract!(Val(:TBLIS), R, labelsR, T1, labelsT1, T2, labelsT2, α, β)
    #end
    return R
  end

  if N1 + N2 == NR
    outer!(R, T1, T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR, labelsRp)
    if !is_trivial_permutation(perm)
      permutedims!(R, copy(R), perm)
    end
    return R
  end

  props = ContractionProperties(labelsT1, labelsT2, labelsR)
  compute_contraction_properties!(props, T1, T2, R)

  if ElT1 != ElT2
    # TODO: use promote instead
    # T1, T2 = promote(T1, T2)

    ElT1T2 = promote_type(ElT1, ElT2)
    if ElT1 != ElR
      # TODO: get this working
      # T1 = ElR.(T1)
      T1 = one(ElT1T2) * T1
    end
    if ElT2 != ElR
      # TODO: get this working
      # T2 = ElR.(T2)
      T2 = one(ElT1T2) * T2
    end
  end

  _contract!(R, T1, T2, props, α, β)
  return R
  #end
end

function _contract!(
  CT::DenseTensor{El,NC},
  AT::DenseTensor{El,NA},
  BT::DenseTensor{El,NB},
  props::ContractionProperties,
  α::Number=one(El),
  β::Number=zero(El),
) where {El,NC,NA,NB}
  # TODO: directly use Tensor instead of Array
  C = ReshapedArray(data(storage(CT)), dims(inds(CT)), ())
  A = ReshapedArray(data(storage(AT)), dims(inds(AT)), ())
  B = ReshapedArray(data(storage(BT)), dims(inds(BT)), ())

  tA = 'N'
  if props.permuteA
    pA = NTuple{NA,Int}(props.PA)
    #@timeit_debug timer "_contract!: permutedims A" begin
    @strided Ap = permutedims(A, pA)
    #end # @timeit
    AM = ReshapedArray(Ap, (props.dmid, props.dleft), ())
    tA = 'T'
  else
    #A doesn't have to be permuted
    if Atrans(props)
      AM = ReshapedArray(A.parent, (props.dmid, props.dleft), ())
      tA = 'T'
    else
      AM = ReshapedArray(A.parent, (props.dleft, props.dmid), ())
    end
  end

  tB = 'N'
  if props.permuteB
    pB = NTuple{NB,Int}(props.PB)
    #@timeit_debug timer "_contract!: permutedims B" begin
    @strided Bp = permutedims(B, pB)
    #end # @timeit
    BM = ReshapedArray(Bp, (props.dmid, props.dright), ())
  else
    if Btrans(props)
      BM = ReshapedArray(B.parent, (props.dright, props.dmid), ())
      tB = 'T'
    else
      BM = ReshapedArray(B.parent, (props.dmid, props.dright), ())
    end
  end

  # TODO: this logic may be wrong
  if props.permuteC
    # Need to copy here since we will be permuting
    # into C later
    CM = ReshapedArray(copy(C), (props.dleft, props.dright), ())
  else
    if Ctrans(props)
      CM = ReshapedArray(C.parent, (props.dright, props.dleft), ())
      (AM, BM) = (BM, AM)
      if tA == tB
        tA = tB = (tA == 'T' ? 'N' : 'T')
      end
    else
      CM = ReshapedArray(C.parent, (props.dleft, props.dright), ())
    end
  end

  _gemm!(tA, tB, El(α), AM, BM, El(β), CM)

  if props.permuteC
    pC = NTuple{NC,Int}(props.PC)
    Cr = ReshapedArray(CM.parent, props.newCrange, ())
    # TODO: use invperm(pC) here?
    #@timeit_debug timer "_contract!: permutedims C" begin
    @strided C .= permutedims(Cr, pC)
    #end # @timeit
  end

  return CT
end

"""
    NDTensors.permute_reshape(T::Tensor,pos...)

Takes a permutation that is split up into tuples. Index positions
within the tuples are combined.

For example:

permute_reshape(T,(3,2),1)

First T is permuted as `permutedims(3,2,1)`, then reshaped such
that the original indices 3 and 2 are combined.
"""
function permute_reshape(
  T::DenseTensor{ElT,NT,IndsT}, pos::Vararg{<:Any,N}
) where {ElT,NT,IndsT,N}
  perm = flatten(pos...)

  length(perm) ≠ NT && error("Index positions must add up to order of Tensor ($N)")
  isperm(perm) || error("Index positions must be a permutation")

  dimsT = dims(T)
  indsT = inds(T)
  if !is_trivial_permutation(perm)
    T = permutedims(T, perm)
  end
  if all(p -> length(p) == 1, pos) && N == NT
    return T
  end
  newdims = MVector(ntuple(_ -> eltype(IndsT)(1), Val(N)))
  for i in 1:N
    if length(pos[i]) == 1
      # No reshape needed, just use the
      # original index
      newdims[i] = indsT[pos[i][1]]
    else
      newdim_i = 1
      for p in pos[i]
        newdim_i *= dimsT[p]
      end
      newdims[i] = eltype(IndsT)(newdim_i)
    end
  end
  newinds = similartype(IndsT, Val{N})(Tuple(newdims))
  return reshape(T, newinds)
end

# svd of an order-n tensor according to positions Lpos
# and Rpos
function LinearAlgebra.svd(
  T::DenseTensor{<:Number,N,IndsT}, Lpos::NTuple{NL,Int}, Rpos::NTuple{NR,Int}; kwargs...
) where {N,IndsT,NL,NR}
  M = permute_reshape(T, Lpos, Rpos)
  UM, S, VM, spec = svd(M; kwargs...)
  u = ind(UM, 2)
  v = ind(VM, 2)

  Linds = similartype(IndsT, Val{NL})(ntuple(i -> inds(T)[Lpos[i]], Val(NL)))
  Uinds = push(Linds, u)

  # TODO: do these positions need to be reversed?
  Rinds = similartype(IndsT, Val{NR})(ntuple(i -> inds(T)[Rpos[i]], Val(NR)))
  Vinds = push(Rinds, v)

  U = reshape(UM, Uinds)
  V = reshape(VM, Vinds)

  return U, S, V, spec
end

# qr decomposition of an order-n tensor according to
# positions Lpos and Rpos
function LinearAlgebra.qr(
  T::DenseTensor{<:Number,N,IndsT}, Lpos::NTuple{NL,Int}, Rpos::NTuple{NR,Int}; kwargs...
) where {N,IndsT,NL,NR}
  M = permute_reshape(T, Lpos, Rpos)
  QM, RM = qr(M; kwargs...)
  q = ind(QM, 2)
  r = ind(RM, 1)
  # TODO: simplify this by permuting inds(T) by (Lpos,Rpos)
  # then grab Linds,Rinds
  Linds = similartype(IndsT, Val{NL})(ntuple(i -> inds(T)[Lpos[i]], Val(NL)))
  Qinds = push(Linds, r)
  Q = reshape(QM, Qinds)
  Rinds = similartype(IndsT, Val{NR})(ntuple(i -> inds(T)[Rpos[i]], Val(NR)))
  Rinds = pushfirst(Rinds, r)
  R = reshape(RM, Rinds)
  return Q, R
end

# polar decomposition of an order-n tensor according to positions Lpos
# and Rpos
function polar(
  T::DenseTensor{<:Number,N,IndsT}, Lpos::NTuple{NL,Int}, Rpos::NTuple{NR,Int}
) where {N,IndsT,NL,NR}
  M = permute_reshape(T, Lpos, Rpos)
  UM, PM = polar(M)

  # TODO: turn these into functions
  Linds = similartype(IndsT, Val{NL})(ntuple(i -> inds(T)[Lpos[i]], Val(NL)))
  Rinds = similartype(IndsT, Val{NR})(ntuple(i -> inds(T)[Rpos[i]], Val(NR)))

  # Use sim to create "similar" indices, in case
  # the indices have identifiers. If not this should
  # act as an identity operator
  simRinds = sim(Rinds)
  Uinds = (Linds..., simRinds...)
  Pinds = (simRinds..., Rinds...)

  U = reshape(UM, Uinds)
  P = reshape(PM, Pinds)
  return U, P
end

function LinearAlgebra.exp(
  T::DenseTensor{ElT,N}, Lpos::NTuple{NL,Int}, Rpos::NTuple{NR,Int}; ishermitian::Bool=false
) where {ElT,N,NL,NR}
  M = permute_reshape(T, Lpos, Rpos)
  indsTp = permute(inds(T), (Lpos..., Rpos...))
  if ishermitian
    expM = parent(exp(Hermitian(matrix(M))))
    return tensor(Dense{ElT}(vec(expM)), indsTp)
  else
    expM = exp(M)
    return reshape(expM, indsTp)
  end
end

function HDF5.write(
  parent::Union{HDF5.File,HDF5.Group}, name::String, D::Store
) where {Store<:Dense}
  g = create_group(parent, name)
  attributes(g)["type"] = "Dense{$(eltype(Store))}"
  attributes(g)["version"] = 1
  if eltype(D) != Nothing
    write(g, "data", D.data)
  end
end

function HDF5.read(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{Store}
) where {Store<:Dense}
  g = open_group(parent, name)
  ElT = eltype(Store)
  typestr = "Dense{$ElT}"
  if read(attributes(g)["type"]) != typestr
    error("HDF5 group or file does not contain $typestr data")
  end
  if ElT == Nothing
    return Dense{Nothing}()
  end
  # Attribute __complex__ is attached to the "data" dataset
  # by the h5 library used by C++ version of ITensor:
  if haskey(attributes(g["data"]), "__complex__")
    M = read(g, "data")
    nelt = size(M, 1) * size(M, 2)
    data = Vector(reinterpret(ComplexF64, reshape(M, nelt)))
  else
    data = read(g, "data")
  end
  return Dense{ElT}(data)
end

function show(io::IO, mime::MIME"text/plain", T::DenseTensor)
  summary(io, T)
  return print_tensor(io, T)
end


text:
function integrate!(crr::CorefResolver)
    # top_n matches are taken
    # need to of how many
    space_pair = crr |> main
    crspace = crr.crspace
    crs = crr.crs

    inner_spaces = space_pair |> innerspaces |> collect
    n = space_pair.erspace1.height
    space_fit = get_space_fit(space_pair)
    for k in 1:space_pair.erspace1.width
        name = space_pair.erspace1.entity_types[k]
        for space_ind in 1:n
            for (thisspace, otherspace, lspace) in zip(inner_spaces, reverse(inner_spaces), crr.crspace.latent_match_vars[name])
                i1 = thisspace[space_ind, k]
                i2 = otherspace[space_ind, k]
                proposal = lspace[i1]
                println(i1, i2)
                update!(proposal, i2, crr.crspace.param_log, space_fit)
            end
        end
    end
end

function update!(
    proposal::CorefProposal,
    match::Int,
    steplog::Array{StepParameters, 1},
    space_fit
    )
    proposal.ind = match
end

function get_space_fit(space_pair::SpacePair)
    10
end

ext: jl


text:
module Kuber

using JSON
using Swagger
using HTTP
using Retry, HTTP.IOExtras

include("api/Kubernetes.jl")
using .Kubernetes
import .Kubernetes: getAPIVersions

import Base: convert, get, put!, delete!, show
import Swagger: SwaggerModel

include("typealiases.jl")
include("apialiases.jl")
include("helpers.jl")
include("simpleapi.jl")

export KuberContext, set_server, set_ns, get_server, get_ns, kuber_type, kuber_obj, @K_str
export get, list, put!, update!, delete!, sel, get_logs, list_namespaced_custom_metrics, list_custom_metrics

end # module


text:
function fdist!(f::Vector{Float64},x::Vector{Float64},D::Vector{Float64},C::Float64,
    R0::Float64,R1::Float64,Ws::Float64,W1::Float64,rho::Float64,beta::Float64,
    c0::Float64,a1::Float64,a2::Float64,Pn::Float64,h::Float64)
    f[1] = D[1]*(a1/(a2*Ws^5)*(0.125*(W1-Ws*x[1])+c0)^2 - (W1^5/(8192*Ws^5) -
        3*W1^4*x[1]/(8192*Ws^4) + W1^4*c0/(256*Ws^5) + W1^3*x[1]^2/(4096*Ws^3) -
        W1^3*x[1]*c0/(128*Ws^4) + 3*W1^3*c0^2/(64*Ws^5) + W1^2*x[1]^3/(4096*Ws^2) -
        3*W1^2*x[1]*c0^2/(64*Ws^4) + W1^2*c0^3/(4*Ws^5) - 3*W1*x[1]^4/(8192*Ws) +
        W1*x[1]^3*c0/(128*Ws^2) - 3*W1*x[1]^2*c0^2/(64*Ws^3) + W1*c0^4/(2*Ws^5) +
        x[1]^5/8192 - x[1]^4*c0/(256*Ws) + 3*x[1]^3*c0^2/(64*Ws^2) -
        x[1]^2*c0^3/(4*Ws^3) + x[1]*c0^4/(2*Ws^4)) + Pn/(a2*Ws^5));
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# Brown and Dennis functions
#
#   Source: Problem 16 in
#   J.J. Moré, B.S. Garbow and K.E. Hillstrom,
#   "Testing Unconstrained Optimization Software",
#   ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981
#
#   classification SUR2-AN-4-0

function brownden()

  nlp = Model()

  @variable(nlp, x[1:4])
  setvalue(x, [25.0; 5.0; -5.0; -1.0])

  @NLobjective(
    nlp,
    Min,
    sum( ((x[1] + x[2] * i/5 - exp(i/5))^2 + (x[3] + x[4]*sin(i/5) -
      cos(i/5))^2)^2 for i = 1:20)
  )

  return nlp
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    ArrayShape{T,N} <: AbstractValueShape

Describes the shape of `N`-dimensional arrays of type `T` and a given size.

Constructor:

    ArrayShape{T}(dims::NTuple{N,Integer}) where {T,N}
    ArrayShape{T}(dims::Integer...) where {T}

e.g.

    shape = ArrayShape{Real}(2, 3)

In addition to using the shape as a value constructor

    size(shape(undef)) == (2, 3)
    eltype(shape(undef)) == Float64

(see [`AbstractValueShape`](@ref)), a shape can also be used as an argument of
certains array type constructors to explicitly construct standard `Array`s
or `ElasticArray`s:

    using ElasticArrays

    size(Array(undef, shape)) == (2, 3)
    eltype(Array(undef, shape)) == Float64

    size(ElasticArray(undef, shape)) == (2, 3)

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ArrayShape{T,N} <: AbstractValueShape
    dims::NTuple{N,Int}
end

export ArrayShape

ArrayShape{T}(dims::NTuple{N,Integer}) where {T,N} = ArrayShape{T,N}(map(Int, dims))
ArrayShape{T}(dims::Integer...) where {T} = ArrayShape{T}(dims)

@inline Base.size(shape::ArrayShape) = shape.dims
Base.length(shape::ArrayShape) = prod(size(shape))

@inline Base.eltype(::ArrayShape{T}) where {T} = T

@inline _valshapeoftype(T::Type{<:AbstractArray}) = throw(ArgumentError("Type $T does not have a fixed shape"))

@inline function valshape(x::AbstractArray{T}) where T
    _valshapeoftype(T) # ensure T has a fixed shape
    ArrayShape{T}(size(x))
end

# Possible extension: valshape(x::AbstractArrayOfSimilarArrays{...})

totalndof(shape::ArrayShape{T}) where{T} =
    prod(size(shape)) * totalndof(_valshapeoftype(T))

(shape::ArrayShape{T,N})(::UndefInitializer) where {T,N} = Array{nonabstract_eltype(shape),N}(undef, size(shape)...)

@static if VERSION < v"1.3"
    # Workaround for Julia issue #14919
    @inline (shape::ArrayShape)(data::AbstractVector{<:Real}) =
        _apply_shape_to_data(shape, data)
end

@inline Array{U}(::UndefInitializer, shape::ArrayShape{T}) where {T,U<:T} =
    Array{U}(undef, size(shape)...)

@inline Array(::UndefInitializer, shape::ArrayShape) =
    Array{nonabstract_eltype(shape)}(undef, shape)

@inline ElasticArray{U}(::UndefInitializer, shape::ArrayShape{T}) where {T,U<:T} =
    ElasticArray{U}(undef, size(shape)...)

@inline ElasticArray(::UndefInitializer, shape::ArrayShape) =
    ElasticArray{nonabstract_eltype(shape)}(undef, shape)

const ArrayAccessor{T,N} = ValueAccessor{ArrayShape{T,N}} where {T,N}

Base.@propagate_inbounds Base.getindex(data::AbstractVector{<:Real}, va::ArrayAccessor) = view(data, va)

Base.@propagate_inbounds Base.view(data::AbstractVector{<:Real}, va::ArrayAccessor{T,1}) where {T} =
    view(data, view_idxs(axes(data, 1), va))

Base.@propagate_inbounds Base.view(data::AbstractVector{<:Real}, va::ArrayAccessor{T,N}) where {T,N} =
    reshape(view(data, view_idxs(axes(data, 1), va)), size(va.shape)...)

Base.@propagate_inbounds function _bcasted_getindex(data::AbstractVectorOfSimilarVectors{<:Real}, va::ArrayAccessor{T,1}) where {T,N}
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    fpview = view(flat_data, idxs, :)
    VectorOfSimilarVectors(fpview)
end

Base.@propagate_inbounds function _bcasted_getindex(data::AbstractVectorOfSimilarVectors{<:Real}, va::ArrayAccessor{T,N}) where {T,N}
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    fpview = view(flat_data, idxs, :)
    VectorOfSimilarArrays(reshape(fpview, size(va.shape)..., :))
end

Base.copy(instance::VSBroadcasted2{typeof(getindex),AbstractVectorOfSimilarVectors{<:Real},Ref{<:ArrayAccessor}}) =
    _bcasted_getindex(instance.args[1], instance.args[2][])

# TODO: Add support for StaticArray.

# Possible extension: variable/flexible array shapes?


text:
using Pkg
Pkg.add("Documenter")
using Documenter, DSGE

makedocs(modules = [DSGE],
         clean = false,
         format = Documenter.HTML(),
         sitename = "DSGE.jl",
         authors = "FRBNY-DSGE",
         linkcheck = false,
         strict = false,
         pages = Any[
                     "Home"                                   => "index.md",
                     "Model Design"                           => "model_design.md",
                     "Special Model Types"                    => "special_model_types.md",
                     "Running An Existing Model"              => "running_existing_model.md",
                     "Advanced Usage"                         => "advanced_usage.md",
                     "Input Data"                             => "input_data.md",
                     "FRBNY Model Input Data"                 => "frbny_data.md",
                     "Implementation Details"                 => "implementation_details.md",
                     "Solving the Model"                      => "solving.md",
                     "Estimation"                             => "estimation.md",
                     "Forecasting"                            => "forecast.md",
                     "Impulse Response Functions"             => "irf.md",
                     "Computing Means and Bands"              => "means_bands.md",
                     "Alternative Policies"                   => "altpolicy.md",
                     "Alternative Scenarios"                  => "scenarios.md",
                     "Forecast Decomposition"                 => "forecast_decomposition.md",
                     "Plotting"                               => "plotting.md",
                     "Algorithms"                             => "algorithms.md",
                     "Contributing to DSGE.jl"                => "contributing.md",
                     "MATLAB to Julia Transition: Estimation" => "MatlabToJuliaTransition.md",
                     "MATLAB to Julia Transition: Forecast"   => "julia_forecasting.md",
                     "License"                                => "license.md"
         ],
         doctest = false # for now
)

# if "deploy" in ARGS
#    include("../../fake_travis.jl")
# end
deploydocs(
    repo = "github.com/FRBNY-DSGE/DSGE.jl.git",
    target = "build",
    deps = nothing,
    devbranch = "master",
    branch = "gh-pages",
    # versions = "v#",
    # julia = "0.7",
    # osname = "osx",
    make = nothing
)


text:
"""
"""
function ConformingFESpace(
  reffes::Vector{<:ReferenceFE},
  conf::Conformity,
  model::DiscreteModel,
  face_labeling::FaceLabeling,
  dirichlet_tags,
  dirichlet_components=nothing,
  is_ref=true)

  grid_topology = get_grid_topology(model)

  _dirichlet_components = _convert_dirichlet_components(dirichlet_tags,dirichlet_components)

  cell_dofs, nfree, ndirichlet, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
    reffes,conf,grid_topology,face_labeling,dirichlet_tags,_dirichlet_components)

  ntags = length(dirichlet_tags)

  grid = get_grid(model)
  cell_to_ctype = get_cell_type(grid_topology)
  cell_map = get_cell_map(grid)

  cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map, Val(is_ref))

  UnconstrainedFESpace(
    nfree,
    ndirichlet,
    cell_dofs,
    cell_shapefuns,
    cell_dof_basis,
    dirichlet_dof_tag,
    dirichlet_cells,
    ntags)

end

function ConformingFESpace(
  reffes::Vector{<:ReferenceFE},
  model::DiscreteModel,
  face_labeling::FaceLabeling,
  dirichlet_tags,
  dirichlet_components=nothing,
  is_ref=true)

  conf = get_default_conformity(first(reffes))
  ConformingFESpace(
    reffes,
    conf,
    model,
    face_labeling,
    dirichlet_tags,
    dirichlet_components,
    is_ref)
end

"""
"""
function GradConformingFESpace(
  reffes::Vector{<:ReferenceFE},
  model::DiscreteModel,
  dirichlet_tags,
  dirichlet_components=nothing,
  is_ref=true)

  face_labeling = get_face_labeling(model)

  ConformingFESpace(
    reffes,GradConformity(),model,face_labeling,dirichlet_tags,dirichlet_components,is_ref)

end

function GradConformingFESpace(
  ::Type{T},
  model::DiscreteModel,
  order::Integer,
  dirichlet_tags,
  dirichlet_components=nothing,
  is_ref=true) where T

  face_labeling = get_face_labeling(model)

  GradConformingFESpace(T,model,order,face_labeling,dirichlet_tags,dirichlet_components,is_ref)

end

function GradConformingFESpace(
  ::Type{T},
  model::DiscreteModel,
  order::Integer,
  face_labeling::FaceLabeling,
  dirichlet_tags,
  dirichlet_components=nothing,
  is_ref=true) where T

  grid_topology = get_grid_topology(model)
  polytopes = get_polytopes(grid_topology)

  reffes = [ LagrangianRefFE(T,p,order) for p in polytopes ]

  ConformingFESpace(
    reffes,GradConformity(),model,face_labeling,dirichlet_tags,dirichlet_components,is_ref)

end

function compute_cell_space(reffes, cell_to_ctype, cell_map,ref_style::Val{true})
  compute_cell_space(reffes, cell_to_ctype, cell_map)
end

function compute_cell_space(reffes, cell_to_ctype, cell_map,ref_style::Val{false})
  compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
end

"""
"""
function compute_cell_space(reffes, cell_to_ctype, cell_map)

  dof_basis = map(get_dof_basis,reffes)
  cell_dof_basis = CompressedArray(dof_basis,cell_to_ctype)
  cell_dof_basis = GenericCellDofBasis(Val{true}(),cell_dof_basis)

  shapefuns =  map(get_shapefuns,reffes)
  refshapefuns = CompressedArray(shapefuns,cell_to_ctype)
  cell_shapefuns = attachmap(refshapefuns,cell_map)

  reffe_to_axs = map(reffe->(Base.OneTo(num_dofs(reffe)),),reffes)
  cell_axs = CompressedArray(reffe_to_axs,cell_to_ctype)
  cell_shapefuns = GenericCellField(cell_shapefuns,cell_map,Val{true}(),cell_axs,Val((:,)))

  (cell_shapefuns, cell_dof_basis)
end

"""
It creates the cell-wise DOF basis and shape functions, when the DOFs
are evaluated at the physical space. The DOFs (moments) for the prebasis
are assumed to be computable at a reference FE space.
"""
function compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
# E.g., if one has to implement $\int_F q ϕ_h(x)$ for $q \in P(F)$, we can
# assume that it can be written as $\sum_{x_{gp} \in Q}_{\hat{F}}
# \hat{q}(\tilde{x}_{gp}) ϕ(x_{gp})$ for $\hat{q} \in P(\hat{F})$.

  dof_bases = map(get_dof_basis,reffes)

  cell_dof_basis = _cell_dof_basis_physical_space(dof_bases,cell_to_ctype,cell_map)
  ref_style = Val{false}()
  cell_dof_basis = GenericCellDofBasis(ref_style,cell_dof_basis)

  prebasis =  map(get_prebasis,reffes)
  cell_prebasis = CompressedArray(prebasis,cell_to_ctype)
  reffe_to_axs = map(reffe->(Base.OneTo(num_dofs(reffe)),),reffes)
  cell_axs = CompressedArray(reffe_to_axs,cell_to_ctype)
  cell_prebasis = GenericCellField(cell_prebasis,cell_map,Val{false}(),cell_axs,Val((:,)))
  cell_shapefuns = _cell_shape_functions_physical_space(cell_prebasis,cell_dof_basis,cell_map)

  (cell_shapefuns, cell_dof_basis)
end

function _cell_dof_basis_physical_space(dof_bases,cell_to_ctype,cell_map)
  @notimplemented
end

function _cell_dof_basis_physical_space(
  dof_bases::Vector{MomentBasedDofBasis{T,V}},
  cell_to_ctype,cell_map) where {T,V}

  ctype_to_refnodes = map(get_nodes,dof_bases)
  cell_to_refnodes = CompressedArray(ctype_to_refnodes,cell_to_ctype)
  cell_physnodes = evaluate(cell_map,cell_to_refnodes)

  # Not efficient, create a Kernel
  ct_face_moments = map(get_face_moments,dof_bases)
  c_face_moments = CompressedArray(ct_face_moments,cell_to_ctype)
  ct_face_nodes_dofs = map(get_face_nodes_dofs,dof_bases)
  c_face_nodes_dofs = CompressedArray(ct_face_nodes_dofs,cell_to_ctype)
  cell_dof_basis = apply( (n,m,nd) -> MomentBasedDofBasis(n,m,nd),
                          cell_physnodes, c_face_moments, c_face_nodes_dofs)
end

function _cell_dof_basis_physical_space(
  dof_bases::Vector{LagrangianDofBasis{T,V}},
  cell_to_ctype,cell_map) where {T,V}

  # Create new dof_basis with nodes in the physical space
  # ctype_to_refnodes= map(get_node_coordinates,reffes)
  ctype_to_refnodes= map(get_nodes,dof_bases)
  cell_to_refnodes = CompressedArray(ctype_to_refnodes,cell_to_ctype)
  cell_physnodes = evaluate(cell_map,cell_to_refnodes)
  cell_dof_basis = apply( nodes -> LagrangianDofBasis(V,nodes), cell_physnodes )
end

function _cell_shape_functions_physical_space(cell_prebasis,cell_dof_basis,cell_map)
  cell_matrix = evaluate(cell_dof_basis,cell_prebasis)
  cell_matrix_inv = apply(inv,cell_matrix)
  cell_shapefuns_phys = apply(change_basis,get_array(cell_prebasis),cell_matrix_inv)
  cell_axs = get_cell_axes(cell_prebasis)
  metasize_style = MetaSizeStyle(cell_prebasis)
  ref_style = RefStyle(cell_prebasis)
  GenericCellField(cell_shapefuns_phys,cell_map,ref_style,cell_axs,metasize_style)
   # @santiagobadia : better implementation in the future...
end

"""
  compute_conforming_cell_dofs(
    reffes,
    conf,
    grid_topology,
    face_labeling,
    dirichlet_tags)

  compute_conforming_cell_dofs(
    reffes,
    conf,
    grid_topology,
    face_labeling,
    dirichlet_tags,
    dirichlet_components)

The result is the tuple

    (cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells)

Assumes that the reffes are aligned with the cell type in the grid_topology
and that it is possible to build a conforming space without imposing constraints

If `dirichlet_components`  is given, then `get_dof_to_comp` has to be defined
for the reference elements in `reffes`.
"""
function compute_conforming_cell_dofs(
  reffes,conf,grid_topology,face_labeling,dirichlet_tags,dirichlet_components=nothing)

  D = num_cell_dims(grid_topology)
  n_faces = num_faces(grid_topology)
  cell_to_ctype = get_cell_type(grid_topology)
  d_to_cell_to_dfaces = [ Table(get_faces(grid_topology,D,d)) for d in 0:D]
  d_to_dface_to_cells = [ Table(get_faces(grid_topology,d,D)) for d in 0:D]
  d_to_offset = get_offsets(grid_topology)
  d_to_ctype_to_ldface_to_own_ldofs = [
    [ get_face_own_dofs(reffe,conf,d) for reffe in reffes] for d in 0:D]

  face_to_own_dofs, ntotal, d_to_dface_to_cell, d_to_dface_to_ldface =  _generate_face_to_own_dofs(
     n_faces,
     cell_to_ctype,
     d_to_cell_to_dfaces,
     d_to_dface_to_cells,
     d_to_offset,
     d_to_ctype_to_ldface_to_own_ldofs)

  d_to_dface_to_tag = [ get_face_tag_index(face_labeling,dirichlet_tags,d)  for d in 0:D]
  cell_to_faces = Table(get_cell_faces(grid_topology))

  nfree, ndiri, diri_dof_tag = _split_face_own_dofs_into_free_and_dirichlet!(
    face_to_own_dofs,
    d_to_offset,
    d_to_dface_to_tag,
    d_to_dface_to_cell,
    d_to_dface_to_ldface,
    cell_to_ctype,
    reffes,
    dirichlet_components)

  cell_to_lface_to_pindex = Table(get_cell_permutations(grid_topology))
  ctype_to_lface_to_own_ldofs = [ get_face_own_dofs(reffe,conf) for reffe in reffes]
  ctype_to_num_dofs = [num_dofs(reffe) for reffe in reffes]
  ctype_to_lface_to_pindex_to_pdofs = [get_face_own_dofs_permutations(reffe,conf) for reffe in reffes]

  cell_dofs = CellDofsNonOriented(
    cell_to_faces,
    cell_to_lface_to_pindex,
    cell_to_ctype,
    ctype_to_lface_to_own_ldofs,
    ctype_to_num_dofs,
    face_to_own_dofs,
    ctype_to_lface_to_pindex_to_pdofs)

  diri_cells = _generate_diri_cells(
    d_to_dface_to_tag,
    d_to_cell_to_dfaces)

  (cell_dofs, nfree, ndiri, diri_dof_tag, diri_cells)
end

function _generate_face_to_own_dofs(
  n_faces,
  cell_to_ctype,
  d_to_cell_to_dfaces::Vector{Table{T,Vd,Vp}},
  d_to_dface_to_cells::Vector{Table{T,Vd,Vp}},
  d_to_offset,
  d_to_ctype_to_ldface_to_own_ldofs) where {T,Vd,Vp}

  P=eltype(Vp)
  face_to_own_dofs_ptrs = zeros(P,n_faces+1)

  D = length(d_to_offset)-1

  icell = 1
  d_to_dface_to_cell = [ get_local_item(d_to_dface_to_cells[d+1],icell)  for d in 0:D ]

  d_to_dface_to_ldface = [
    find_local_index(d_to_dface_to_cell[d+1],d_to_cell_to_dfaces[d+1]) for d in 0:D ]

  for d in 0:D
    cell_to_dfaces = d_to_cell_to_dfaces[d+1]
    dface_to_cells = d_to_dface_to_cells[d+1]
    offset = d_to_offset[d+1]
    ctype_to_ldface_to_own_ldofs = d_to_ctype_to_ldface_to_own_ldofs[d+1]
    ctype_to_ldface_to_num_own_ldofs = map( (x) -> length.(x) ,ctype_to_ldface_to_own_ldofs)
    dface_to_cell_owner = d_to_dface_to_cell[d+1]
    dface_to_ldface = d_to_dface_to_ldface[d+1]

    if any( ctype_to_ldface_to_num_own_ldofs .!= 0)
      _generate_face_to_own_dofs_count_d!(
        face_to_own_dofs_ptrs,
        offset,
        cell_to_ctype,
        dface_to_cell_owner,
        dface_to_ldface,
        ctype_to_ldface_to_num_own_ldofs)
    end
  end

  length_to_ptrs!(face_to_own_dofs_ptrs)

  n_dofs = face_to_own_dofs_ptrs[end]-1
  face_to_own_dofs_data = collect(T(1):T(n_dofs))

  face_to_own_dofs = Table(face_to_own_dofs_data,face_to_own_dofs_ptrs)
  (face_to_own_dofs, n_dofs, d_to_dface_to_cell, d_to_dface_to_ldface)
end

function  _generate_face_to_own_dofs_count_d!(
  face_to_own_dofs_ptrs,
  offset,
  cell_to_ctype,
  dface_to_cell_owner,
  dface_to_ldface,
  ctype_to_ldface_to_num_own_ldofs)

  n_dfaces = length(dface_to_ldface)
  for dface in 1:n_dfaces
    cell = dface_to_cell_owner[dface]
    ldface = dface_to_ldface[dface]
    ctype = cell_to_ctype[cell]
    n_own_ldofs = ctype_to_ldface_to_num_own_ldofs[ctype][ldface]
    face = dface + offset
    face_to_own_dofs_ptrs[face+1] = n_own_ldofs
  end
end

function _split_face_own_dofs_into_free_and_dirichlet!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag,
  d_to_dface_to_cell,
  d_to_dface_to_ldface,
  cell_to_ctype,
  reffes,
  dirichlet_components::Nothing)

  _split_face_own_dofs_into_free_and_dirichlet_generic!(
    face_to_own_dofs,
    d_to_offset,
    d_to_dface_to_tag)

end

function _split_face_own_dofs_into_free_and_dirichlet!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag,
  d_to_dface_to_cell,
  d_to_dface_to_ldface,
  cell_to_ctype,
  reffes,
  dirichlet_components)

  D = length(d_to_dface_to_cell)-1

  d_to_ctype_to_ldface_to_own_ldofs = [
    [ get_face_own_dofs(reffe,d) for reffe in reffes] for d in 0:D]

  ctype_to_ldof_to_comp = [get_dof_to_comp(reffe) for reffe in reffes]

  _split_face_own_dofs_into_free_and_dirichlet_with_components!(
    face_to_own_dofs,
    d_to_offset,
    d_to_dface_to_tag,
    d_to_dface_to_cell,
    d_to_dface_to_ldface,
    cell_to_ctype,
    ctype_to_ldof_to_comp,
    d_to_ctype_to_ldface_to_own_ldofs,
    dirichlet_components)
end

function _split_face_own_dofs_into_free_and_dirichlet_generic!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag)

  D = length(d_to_offset)-1
  nfree = 0
  ndiri = 0
  dirichlet_dof_tag = Int8[]
  for d in 0:D
    dface_to_tag = d_to_dface_to_tag[d+1]
    offset = d_to_offset[d+1]
    for (dface, tag) in enumerate(dface_to_tag)
      face = dface + offset
      pini = face_to_own_dofs.ptrs[face]
      pend = face_to_own_dofs.ptrs[face+1]-1
      if tag == UNSET
        for p in pini:pend
          nfree += 1
          face_to_own_dofs.data[p] = nfree
        end
      else
        for p in pini:pend
          ndiri += 1
          face_to_own_dofs.data[p] = -ndiri
          push!(dirichlet_dof_tag,tag)
        end
      end
    end
  end

  (nfree, ndiri, dirichlet_dof_tag)
end

function _split_face_own_dofs_into_free_and_dirichlet_with_components!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag,
  d_to_dface_to_cell,
  d_to_dface_to_ldface,
  cell_to_ctype,
  ctype_to_ldof_to_comp,
  d_to_ctype_to_ldface_to_own_ldofs,
  dirichlet_components)

  D = length(d_to_offset)-1
  nfree = 0
  ndiri = 0
  dirichlet_dof_tag = Int8[]
  for d in 0:D
    dface_to_tag = d_to_dface_to_tag[d+1]
    offset = d_to_offset[d+1]
    dface_to_cell = d_to_dface_to_cell[d+1]
    dface_to_ldface = d_to_dface_to_ldface[d+1]
    ctype_to_ldface_to_own_ldofs = d_to_ctype_to_ldface_to_own_ldofs[d+1]
    ctype_to_ldface_to_num_own_ldofs = map( (x) -> length.(x) ,ctype_to_ldface_to_own_ldofs)
    if all(ctype_to_ldface_to_num_own_ldofs .== 0)
      continue
    end
    for (dface, diritag) in enumerate(dface_to_tag)
      face = dface + offset
      pini = face_to_own_dofs.ptrs[face]
      pend = face_to_own_dofs.ptrs[face+1]-1
      if diritag == UNSET
        for p in pini:pend
          nfree += 1
          face_to_own_dofs.data[p] = nfree
        end
      else
        cell = dface_to_cell[dface]
        ldface = dface_to_ldface[dface]
        ctype = cell_to_ctype[cell]
        ldof_to_comp = ctype_to_ldof_to_comp[ctype]
        own_ldofs = ctype_to_ldface_to_own_ldofs[ctype][ldface]
        comp_to_isdiri = dirichlet_components[diritag]
        for (fldof, p) in enumerate(pini:pend)
          ldof = own_ldofs[fldof]
          comp = ldof_to_comp[ldof]
          isdiri = comp_to_isdiri[comp]
          if isdiri
            ndiri += 1
            face_to_own_dofs.data[p] = -ndiri
            push!(dirichlet_dof_tag,diritag)
          else
            nfree += 1
            face_to_own_dofs.data[p] = nfree
          end
        end
      end
    end
  end

  (nfree, ndiri, dirichlet_dof_tag)
end

function _generate_diri_cells(
  d_to_dface_to_tag,
  d_to_cell_to_dfaces)

  ncells = length(first(d_to_cell_to_dfaces))
  cell_visited = fill(false,ncells)

  diri_cells = Int[]

  D = length(d_to_dface_to_tag)-1
  for d in 0:D
    dface_to_tag = d_to_dface_to_tag[d+1]
    cell_to_dfaces = Table(d_to_cell_to_dfaces[d+1])
    for cell in 1:length(cell_to_dfaces)
      if cell_visited[cell]
        continue
      end
      pini = cell_to_dfaces.ptrs[cell]
      pend = cell_to_dfaces.ptrs[cell+1]-1
      for p in pini:pend
        dface = cell_to_dfaces.data[p]
        tag = dface_to_tag[dface]
        if tag != UNSET
          push!(diri_cells,cell)
          cell_visited[cell] = true
          break
        end
      end
    end
  end

  diri_cells

end

function _convert_dirichlet_components(dirichlet_tags,dirichlet_components)
  dirichlet_components
end

function _convert_dirichlet_components(dirichlet_tags::Integer,dirichlet_components)
  [dirichlet_components,]
end

function _convert_dirichlet_components(dirichlet_tags::Integer,dirichlet_components::Nothing)
  nothing
end

function _convert_dirichlet_components(dirichlet_tags::String,dirichlet_components)
  [dirichlet_components,]
end

function _convert_dirichlet_components(dirichlet_tags::String,dirichlet_components::Nothing)
  nothing
end

struct CellDofsNonOriented <:AbstractVector{Vector{Int}}
  cell_to_faces::Table{Int,Vector{Int},Vector{Int32}}
  cell_to_lface_to_pindex::Table{Int8,Vector{Int8},Vector{Int32}}
  cell_to_ctype::Vector{Int8}
  ctype_to_lface_to_own_ldofs::Vector{Vector{Vector{Int}}}
  ctype_to_num_dofs::Vector{Int}
  face_to_own_dofs::Table{Int,Vector{Int},Vector{Int32}}
  ctype_to_lface_to_pindex_to_pdofs::Vector{Vector{Vector{Vector{Int}}}}
end

Base.size(a::CellDofsNonOriented) = size(a.cell_to_faces)

Base.IndexStyle(::Type{<:CellDofsNonOriented}) = IndexStyle(Table)

function array_cache(a::CellDofsNonOriented)
  n_dofs = testitem(a.ctype_to_num_dofs)
  T = eltype(eltype(a))
  v = zeros(T,n_dofs)
  CachedArray(v)
end

function getindex!(cache,a::CellDofsNonOriented,cell::Integer)
  ctype = a.cell_to_ctype[cell]
  n_dofs = a.ctype_to_num_dofs[ctype]
  setsize!(cache,(n_dofs,))
  dofs = cache.array
  lface_to_own_ldofs = a.ctype_to_lface_to_own_ldofs[ctype]
  p = a.cell_to_faces.ptrs[cell]-1
  for (lface, own_ldofs) in enumerate(lface_to_own_ldofs)
    face = a.cell_to_faces.data[p+lface]
    pindex = a.cell_to_lface_to_pindex.data[p+lface]
    pdofs = a.ctype_to_lface_to_pindex_to_pdofs[ctype][lface][pindex]

    q = a.face_to_own_dofs.ptrs[face]-1
    for (i,ldof) in enumerate(own_ldofs)
      j = pdofs[i]
      dof = a.face_to_own_dofs.data[q+j]
      dofs[ldof] = dof
    end
  end
  dofs
end

function Base.getindex(a::CellDofsNonOriented,cell::Integer)
  cache = array_cache(a)
  getindex!(cache,a,cell)
end

###struct CellDofsOriented{T,P,V<:AbstractVector} <: AbstractVector{Vector{T}}
###  cell_to_faces::Table{T,P}
###  cell_to_ctype::V
###  ctype_to_lface_to_own_ldofs::Vector{Vector{Vector{Int}}}
###  ctype_to_num_dofs::Vector{Int}
###  face_to_own_dofs::Table{T,P}
###end
###
###Base.size(a::CellDofsOriented) = size(a.cell_to_faces)
###
###Base.IndexStyle(::Type{<:CellDofsOriented}) = IndexStyle(Table)
###
###function array_cache(a::CellDofsOriented)
###  n_dofs = testitem(a.ctype_to_num_dofs)
###  T = eltype(eltype(a))
###  v = zeros(T,n_dofs)
###  CachedArray(v)
###end
###
###function getindex!(cache,a::CellDofsOriented,cell::Integer)
###  ctype = a.cell_to_ctype[cell]
###  n_dofs = a.ctype_to_num_dofs[ctype]
###  setsize!(cache,(n_dofs,))
###  dofs = cache.array
###  lface_to_own_ldofs = a.ctype_to_lface_to_own_ldofs[ctype]
###  p = a.cell_to_faces.ptrs[cell]-1
###  for (lface, own_ldofs) in enumerate(lface_to_own_ldofs)
###    face = a.cell_to_faces.data[p+lface]
###    q = a.face_to_own_dofs.ptrs[face]-1
###    for (i,ldof) in enumerate(own_ldofs)
###      dof = a.face_to_own_dofs.data[q+i]
###      dofs[ldof] = dof
###    end
###  end
###  dofs
###end
###
###function Base.getindex(a::CellDofsOriented,cell::Integer)
###  cache = array_cache(a)
###  getindex!(cache,a,cell)
###end
###


text:
struct AbstractSimplex{T} end

export Simplex
struct Simplex{T}  <: AbstractVector{T}
    weights::Vector{T}
end

function Simplex(w)
    N = length(w)
    T = eltype(w)
    Simplex{T}(w)
end

Base.length(s::Simplex{T}) where {T} = length(s.weights)
Base.size(s::Simplex{T}) where {T} = (length(s.weights),)
Base.getindex(s::Simplex{T}, n) where {T} = s.weights[n]

struct Log{T}
    val :: T
end

function Base.:*(x::Log{Tx}, y::Log{Ty}) where {Tx <: Real, Ty <: Real}
    Log(x.val + y.val)
end

Base.log(x::Log{T}) where {T <: Real} = x.val

struct AbstractMix end

# struct MixFor
#     dists :: For
#     weights :: Simplex
# end

struct MixVec{D,W}
    dists :: Vector{D}
    weights :: Simplex{W}
end

export Mix

function Mix(dists :: Vector{D}, weights :: Vector{W}) where {D,W}
    MixVec{D,W}(dists, Simplex(weights))
end

# Mix(w::Vector) = dists -> Mix(dists, log.(w))

function Base.rand(mix::MixVec)
    # This is the "Gumbel max trick"  for categorical sampling
    (j_max,lw_max) = (0,-Inf)

    for j in eachindex(mix.weights)
        lw_gumbel = log(mix.weights[j]) + rand(Gumbel())
        if lw_gumbel > lw_max
            (j_max,lw_max) = (j, lw_gumbel)
        end
    end

    @inbounds mix.dists[j_max] |> rand
end

# function Base.rand(mix::Mix, N::Int)
#     x1 = rand(mix)
#     x = Vector{typeof(x1)}(undef, N)
#     @inbounds x[1] = x1
#     @inbounds for n in 2:N
#         x[n] = rand(mix)
#     end
#     x
# end

xform(mix::MixVec) = xform(mix.dists[1])

function Dists.logpdf(mix::MixVec, x)
    ℓ = 0.0
    @simd for j in eachindex(mix.weights)
        @inbounds ℓ += log(mix.weights[j]) + logdensity(mix.dists[j], x)
    end
    ℓ
end


text:
using GCMAES, Plots

xs = Float32.(collect(-6:0.1:6))
lo = fill(1f0, size(xs))
hi = fill(2f0, size(xs))
trans = BoxLinQuadTransform(lo, hi)
plot(xs, GCMAES.transform(trans, xs))
plot(xs, GCMAES.inverse(trans, xs))


text:
import Base.sprint
export Gol, game, glidergun

mutable struct Gol
    state::Matrix{Int8}
    isperiodic::Bool
    generations::Int
end

function Gol(n::Int,
             m::Int = n;
             isperiodic::Bool = true,
             generations::Int = 100)
    state = rand(Int8[0,1], n, m)

    return Gol(state, isperiodic, generations)
end

"count the number of live neighbors"
function determine_neighbor_state(d::Gol, i::Integer, j::Integer)
    n, m = size(d.state)
    neighbor = [(p,q) for p in -1:1 for q in -1:1 if !(p == 0 && q == 0)]
    total = Int8(0)
    if d.isperiodic # periodic
        for idx in 1:8
            total += d.state[mod1(i + neighbor[idx][1], n), mod1(j + neighbor[idx][2], m)]
        end
    else # non-periodic
        for idx in 1:8
            if !(i + neighbor[idx][1] <= 0 || j + neighbor[idx][2] <= 0 || i + neighbor[idx][1] > n || j + neighbor[idx][2] > m)
                total += d.state[i + neighbor[idx][1], j + neighbor[idx][2]]
            end
        end
    end

    return total
end

"update cells"
function next_generation(d::Gol)
    n,m = size(d.state)
    next_state = zeros(Int8, n, m)
    for col in 1:m
        for row in 1:n
            if d.state[row, col] == one(d.state[row, col]) # live
                if determine_neighbor_state(d, row, col) in Int8[2, 3]
                    next_state[row, col] = Int8(1)
                end
            else # dead
                if determine_neighbor_state(d, row, col) == Int8(3)
                    next_state[row, col] = Int8(1)
                end
            end
        end
    end
    return next_state
end

function sprint(d::Gol)
    n, m = size(d.state)
    sm = Matrix{String}(undef,n,m)
    s = ""
    for i in 1:n
        for j in 1:m
            if Bool(d.state[i,j])
                s *= "X"
            else
                s *= " "
            end
        end
        s *= "\n"
    end
    return s = chomp(s)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
# Spectrum calculator for AtmosLES

struct AtmosLESSpectraDiagnosticsParams <: DiagnosticsGroupParams
    nor::Float64
end

"""
    setup_atmos_spectra_diagnostics(
        ::AtmosLESConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
        nor = Inf,
    )

Create the "AtmosLESSpectra" `DiagnosticsGroup` which contains the following
diagnostic variable:

- spectrum: power spectrum from 3D velocity fields

This variable is output with the `k` dimension (wave number) on an interpolated
grid (`interpol` _must_ be specified) as well as a (unlimited) `time` dimension
at the specified `interval`.
"""
function setup_atmos_spectra_diagnostics(
    ::AtmosLESConfigType,
    interval::String,
    out_prefix::String,
    nor::Float64;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosLESSpectra",
        Diagnostics.atmos_les_spectra_init,
        Diagnostics.atmos_les_spectra_fini,
        Diagnostics.atmos_les_spectra_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        AtmosLESSpectraDiagnosticsParams(nor),
    )
end

function get_spectrum(mpicomm, mpirank, Q, bl, interpol, nor)
    FT = eltype(Q)
    istate = similar(Q.data, interpol.Npl, number_states(bl, Prognostic(), FT))
    interpolate_local!(interpol, Q.data, istate)
    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    if mpirank == 0
        u = all_state_data[:, :, :, 2] ./ all_state_data[:, :, :, 1]
        v = all_state_data[:, :, :, 3] ./ all_state_data[:, :, :, 1]
        w = all_state_data[:, :, :, 4] ./ all_state_data[:, :, :, 1]
        x1 = Array(interpol.x1g)
        d = length(x1)
        s, k = power_spectrum_3d(
            AtmosLESConfigType(),
            u,
            v,
            w,
            x1[d] - x1[1],
            d,
            nor,
        )
        return s, k
    end
    return nothing, nothing
end

function atmos_les_spectra_init(dgngrp, currtime)
    Q = Settings.Q
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    FT = eltype(Q)
    interpol = dgngrp.interpol
    nor = dgngrp.params.nor

    spectrum, wavenumber = get_spectrum(mpicomm, mpirank, Q, bl, interpol, nor)
    if mpirank == 0
        dims = OrderedDict("k" => (wavenumber, Dict()))
        vars = OrderedDict("spectrum" => (("k",), FT, Dict()))

        dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
        dfilename = joinpath(Settings.output_dir, dprefix)
        noov = Settings.no_overwrite
        init_data(dgngrp.writer, dfilename, noov, dims, vars)
    end

    return nothing
end

function atmos_les_spectra_collect(dgngrp, currtime)
    Q = Settings.Q
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    FT = eltype(Q)
    interpol = dgngrp.interpol
    nor = dgngrp.params.nor

    spectrum, _ = get_spectrum(mpicomm, mpirank, Q, bl, interpol, nor)

    if mpirank == 0
        varvals = OrderedDict("spectrum" => spectrum)
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function atmos_les_spectra_fini(dgngrp, currtime) end


text:
"""
    BSplines

Module defining B-spline bases and B-spline related functions.
"""
module BSplines

export
    AbstractBSplineBasis,
    BasisFunction,
    BSplineBasis,
    BSplineOrder,
    boundaries,
    order,
    knots,
    basis,
    evaluate,
    evaluate!,
    nonzero_in_segment,
    support

using ..DifferentialOps

"""
    AbstractBSplineBasis{k,T}

Abstract type defining a B-spline basis, or more generally, a functional basis
defined from B-splines.

The basis is represented by a B-spline order `k` and a knot element type `T`.
"""
abstract type AbstractBSplineBasis{k,T} end

"""
    BSplineOrder(k::Integer)

Specifies the B-spline order `k`.
"""
struct BSplineOrder{k} end

@inline BSplineOrder(k::Integer) = BSplineOrder{Int(k)}()

include("basis.jl")
include("basis_function.jl")
include("knots.jl")

end


text:
module Params

export verbose

const verbose = false

end # module


text:
# FrictionTypes.jl

"""
Abstract type defining a basic friction model.
"""
abstract type FrictionModel end

"""
An abstract heuristic friction model.
"""
abstract type HeuristicFrictionModel <: FrictionModel end

"""
A simple Coulomb friction model the represents friction force proportional
to the normal force via a single friction coefficient.
"""
struct CoulombModel <: FrictionModel
    """
    The Coulomb friction coefficient.
    """
    coefficient::Number
end

"""
The Lu-Gre friction model.
"""
struct LuGreModel <: HeuristicFrictionModel
    static_coefficient::Number
    coulomb_coefficient::Number
    stribeck_velocity::Number
    bristle_stiffness::Number
    bristle_damping::Number
    viscous_damping::Number
end

"""
A friction model based upon hyperbolic functions proposed by Rodriguez et. al.
"""
struct HyperbolicModel <: FrictionModel
    friction_coefficient::Number
    normalization_coefficient::Number
    dissipation_coefficient::Number
    hysteresis_coefficient::Number
    stribeck_velocity::Number
    viscous_damping::Number
end

"""
The elasto-plastic model discussed by Dupont et al.
"""
struct ElastoPlasticModel <: HeuristicFrictionModel
    static_coefficient::Number
    coulomb_coefficient::Number
    stribeck_velocity::Number
    bristle_stiffness::Number
    bristle_damping::Number
    breakaway_displacement::Number
    max_bristle_displacement::Number
    viscous_damping::Number
end


text:
"""
    HashDict{K, V} <: AbstractDict{K, V}

Dict-like behavior, stores key-value pairs by key hash.
"""
mutable struct HashDict{K, V} <: AbstractDict{K, V}
    buckets::Vector{Bucket{K, V}}
    number_entries::Int
end

function HashDict{K, V}(; base_size::Int = 1024) where {K, V}
    buckets = [Bucket{K, V}([]) for _ in 1:base_size]
    return HashDict{K, V}(buckets, 0)
end

function HashDict(items::Vector{Tuple{K, V}}; base_size = 1024) where {K, V}
    d = HashDict{K, V}(base_size = base_size)
    for (k, v) in items
        d[k] = v
    end
    return d
end

function HashDict(items::Vector{Pair{K, V}}; base_size = 1024) where {K, V}
    d = HashDict{K, V}(base_size = base_size)
    for (k, v) in items
        d[k] = v
    end
    return d
end

function HashDict(items::Vararg{Pair{K, V}}; base_size = 1024) where {K, V}
    d = HashDict{K, V}(base_size = base_size)
    for (k, v) in items
        d[k] = v
    end
    return d
end

function Base.setindex!(d::HashDict{K, V}, value::V, key::K) where {K, V}
    bucket = find_bucket(d, key)
    (v, replaced) = setindex!(bucket, value, key)
    if !replaced
        d.number_entries += 1
        if d.number_entries > LOAD_FACTOR * length(d.buckets)
            resize!(d)
        end
    end
    return v
end

function Base.getindex(d::HashDict{K, V}, key::K) where {K, V}
    bucket = find_bucket(d, key)
    return bucket[key]
end

function Base.delete!(d::HashDict{K, V}, key::K) where {K, V}
    bucket = find_bucket(d, key)
    found = delete!(bucket, key)
    if found
        d.number_entries -= 1
    end
    return d
end

"""
    replace_for_delete!(d::HashDict{K, V}, key::K, value::V, former_size::Int)

Delete a (key, value) pair from its former location and re-add it to the `HashDict`.
"""
function replace_for_delete!(d::HashDict{K, V}, key::K, value::V, former_size::Int) where {K, V}
    key_hash = hash(key)
    new_idx = 1 + key_hash % length(d.buckets)
    new_bucket = d.buckets[new_idx]
    old_idx = 1 + key_hash % former_size
    old_bucket = d.buckets[old_idx]
    found = delete!(old_bucket, key)
    if found
        new_bucket[key] = value
    end
    return d
end

function Base.empty!(d::HashDict)
    for bucket in d.buckets
        empty!(bucket)
    end
    d.number_entries = 0
    return d
end

function Base.length(d::HashDict)
    return d.number_entries
end

function Base.show(io::IO, d::HashDict{K, V}) where {K, V}
    show(io, typeof(d))
    write(io, " with $(d.number_entries) entries:\n")
    for (k, v) in d
        Base.write(io, "  ")
        Base.show(io, k)
        Base.write(io, " => ")
        Base.show(io, v)
        Base.write(io, "\n")
    end
    return nothing
end

function Base.iterate(d::HashDict)
    if length(d) == 0
        return nothing
    end
    res = nothing
    bucket_idx = 1
    while res === nothing && bucket_idx <= length(d.buckets)
        res = iterate(d.buckets[bucket_idx])
        if res === nothing
            bucket_idx += 1
            continue
        else
            (value, bucket_state) = res
            return (value, (bucket_idx, bucket_state))
        end
    end
    return nothing
end

function Base.iterate(d::HashDict, state::Tuple{Int, Int})
    (bucket_idx, bucket_state) = state
    current_bucket = d.buckets[bucket_idx]
    res = iterate(d.buckets[bucket_idx], bucket_state)
    while res === nothing
        bucket_idx += 1
        if bucket_idx > length(d.buckets)
            return nothing
        end
        res = iterate(d.buckets[bucket_idx])
    end
    (value, next_bucket_state) = res
    return (value, (bucket_idx, next_bucket_state))
end

function Base.get(d::HashDict{K, V}, key::K, default::V) where {K, V}
    bucket = find_bucket(d, key)
    return get(bucket, key, default)
end

function Base.haskey(d::HashDict{K, V}, key::K) where {K, V}
    bucket = find_bucket(d, key)
    return haskey(bucket, key)
end

"""
Maximum average number of items per bucket.
Resizing is performed when this average is reached.
"""
const LOAD_FACTOR = 3

"""
    Base.resize!(d::HashDict)

Double the size of the data underlying the hash dict.
Called when inserting an item goes above `LOAD_FACTOR`.
"""
function Base.resize!(d::HashDict{K, V}) where {K, V}
    l = length(d.buckets)
    resize!(d.buckets, 2l)
    for idx in l+1:2l
        d.buckets[idx] = Bucket{K, V}([])
    end
    for (k, v) in d
        replace_for_delete!(d, k, v, l)
    end
    return d
end

"""
    Base.resize!(d::HashDict{K, V}, capacity::Integer)

Resize a `HashDict` by adding a given capacity.
Becomes a no-op if `capacity <= 0`.
"""
function Base.resize!(d::HashDict{K, V}, capacity::Integer) where {K, V}
    if capacity <= 0
        return d
    end
    l = length(d.buckets)
    new_size = l + capacity
    resize!(d.buckets, new_size)
    for idx in l+1:new_size
        d.buckets[idx] = Bucket{K, V}([])
    end
    for (k, v) in d
        replace_for_delete!(d, k, v, l)
    end
    return d
end

function find_bucket(d::HashDict{K}, key::K) where {K}
    key_hash = hash(key)
    idx = 1 + key_hash % length(d.buckets)
    return d.buckets[idx]
end


text:
#
# Generic diffusivity functions
#

@inline KU(m, i) = KU₀(m) + Cᴷu(m, i) * onface(m.state.K, i)
@inline KV(m, i) = KV₀(m) + Cᴷv(m, i) * onface(m.state.K, i)
@inline KT(m, i) = KT₀(m) + CᴷT(m, i) * onface(m.state.K, i)
@inline KS(m, i) = KS₀(m) + CᴷS(m, i) * onface(m.state.K, i)
@inline Ke(m, i) = Ke₀(m) + Cᴷe(m, i) * onface(m.state.K, i)

#
# Background diffusivity parameters
#

Base.@kwdef struct BackgroundDiffusivities{T} <: AbstractParameters
    KU₀ :: T = 1e-6 # Background viscosity for momentum [m s⁻²]
    KC₀ :: T = 1e-6 # Background diffusivity for tracers and TKE [m s⁻²]
end

@inline KU₀(m) = m.background_diffusivities.KU₀
@inline KV₀(m) = m.background_diffusivities.KU₀
@inline KC₀(m) = m.background_diffusivities.KC₀
@inline KT₀(m) = m.background_diffusivities.KC₀
@inline KS₀(m) = m.background_diffusivities.KC₀
@inline Ke₀(m) = m.background_diffusivities.KC₀

#
# Eddy diffusivity parameters
#

"""
    struct SinglePrandtlDiffusivities{T} <: AbstractParameters

A `simple' Prandtl number model that assumes the TKE diffusivity is identical
to the tracer diffusivity. The momentum diffusivity parameter is Cᴷu, and the
tracer diffusivities equal to the momentum diffusivity divided by CᴷPr.
"""
Base.@kwdef struct SinglePrandtlDiffusivities{T} <: AbstractParameters
     Cᴷu :: T = 0.1   # Diffusivity parameter for velocity
    CᴷPr :: T = 0.74  # Constant Prandtl number for tracers and TKE
end

const SPD = SinglePrandtlDiffusivities

@inline Cᴷu(m::Model{L, <:SPD}, i) where L = m.eddy_diffusivities.Cᴷu
@inline Cᴷv(m::Model{L, <:SPD}, i) where L = m.eddy_diffusivities.Cᴷu
@inline Cᴷc(m::Model{L, <:SPD}, i) where L = m.eddy_diffusivities.Cᴷu / m.eddy_diffusivities.CᴷPr
@inline CᴷT(m::Model{L, <:SPD}, i) where L = Cᴷc(m, i)
@inline CᴷS(m::Model{L, <:SPD}, i) where L = Cᴷc(m, i)
@inline Cᴷe(m::Model{L, <:SPD}, i) where L = Cᴷc(m, i)

"""
    struct IndependentDiffusivities{T} <: AbstractParameters

A diffusivity model in which momentum, tracers, and TKE
each have their own diffusivity parameter.
"""
Base.@kwdef struct IndependentDiffusivities{T} <: AbstractParameters
     Cᴷu :: T = 0.0274 # Diffusivity parameter for velocity
     Cᴷc :: T = 0.0498 # Diffusivity parameter for tracers
     Cᴷe :: T = 0.0329 # Diffusivity parameter for TKE
end

const IDP = IndependentDiffusivities

@inline Cᴷu(m::Model{L, <:IDP}, i) where L = m.eddy_diffusivities.Cᴷu
@inline Cᴷv(m::Model{L, <:IDP}, i) where L = m.eddy_diffusivities.Cᴷu
@inline Cᴷc(m::Model{L, <:IDP}, i) where L = m.eddy_diffusivities.Cᴷc
@inline CᴷT(m::Model{L, <:IDP}, i) where L = Cᴷc(m, i)
@inline CᴷS(m::Model{L, <:IDP}, i) where L = Cᴷc(m, i)
@inline Cᴷe(m::Model{L, <:IDP}, i) where L = m.eddy_diffusivities.Cᴷe

"""
    struct RiDependentDiffusivities{T} <: AbstractParameters

A diffusivity model in which momentum, tracers, and TKE
each have independent Richardson number dependent diffusivities.

The Richardson number is

    ``Ri = ∂z B / ( (∂z U)² + (∂z V)² )`` ,

where ``B`` is buoyancy and ``∂z`` denotes a vertical derviative.

The Richardson-number dependent diffusivities are multiplied by the stability
function

    ``σ(Ri) = σ⁰ + σᵟ * step(Ri, Riᶜ, Riʷ)``

where ``σ⁰``, ``σᵟ``, ``Riᶜ``, and ``Riʷ`` are free parameters,
and ``step`` is a smooth step function defined by

    ``step(x, c, w) = 1/2 * (1 + tanh((x - c) / w))``.
"""
Base.@kwdef struct RiDependentDiffusivities{T} <: AbstractParameters
     Cᴷu⁰ :: T = 0.0274 # Convection diffusivity parameter for velocity
     Cᴷuᵟ :: T = 0.0    # Shift diffusivity parameter for velocity
     Cᴷuᶜ :: T = 0.0    # Ri inflection diffusivity parameter for velocity
     Cᴷuʷ :: T = 1.0    # Ri width Diffusivity parameter for velocity
     Cᴷc⁰ :: T = 0.0498 # Convection diffusivity parameter for tracers
     Cᴷcᵟ :: T = 0.0    # Shift diffusivity parameter for tracers
     Cᴷcᶜ :: T = 0.0    # Ri inflection diffusivity parameter for tracers
     Cᴷcʷ :: T = 1.0    # Ri width diffusivity parameter for tracers
     Cᴷe⁰ :: T = 0.0329 # Convection diffusivity parameter for TKE
     Cᴷeᵟ :: T = 0.0    # Shift diffusivity parameter for TKE
     Cᴷeᶜ :: T = 0.0    # Ri inflection diffusivity parameter for TKE
     Cᴷeʷ :: T = 1.0    # Ri width diffusivity parameter for TKE
end

const RiD = RiDependentDiffusivities

@inline function Richardson_number(m, i)
    N² = ∂B∂z(m, i)
    return ifelse(N² == 0, 0.0, N² / shear_squared(m, i))
end

@inline step(x, c, w) = 1/2 * (1 + tanh((x - c) / w))

@inline stability_function(Ri, σ⁰, σᵟ, Riᶜ, Riʷ) = σ⁰ + σᵟ * step(Ri, Riᶜ, Riʷ)

@inline Cᴷu(m::Model{L, <:RiD}, i) where L = stability_function(
                                                                Richardson_number(m, i),
                                                                m.eddy_diffusivities.Cᴷu⁰,
                                                                m.eddy_diffusivities.Cᴷuᵟ,
                                                                m.eddy_diffusivities.Cᴷuᶜ,
                                                                m.eddy_diffusivities.Cᴷuʷ,
                                                               )

@inline Cᴷv(m::Model{L, <:RiD}, i) where L = Cᴷu(m, i)

@inline Cᴷc(m::Model{L, <:RiD}, i) where L = stability_function(
                                                                Richardson_number(m, i),
                                                                m.eddy_diffusivities.Cᴷc⁰,
                                                                m.eddy_diffusivities.Cᴷcᵟ,
                                                                m.eddy_diffusivities.Cᴷcᶜ,
                                                                m.eddy_diffusivities.Cᴷcʷ,
                                                               )

@inline CᴷT(m::Model{L, <:RiD}, i) where L = Cᴷc(m, i)
@inline CᴷS(m::Model{L, <:RiD}, i) where L = Cᴷc(m, i)

@inline Cᴷe(m::Model{L, <:RiD}, i) where L = stability_function(
                                                                Richardson_number(m, i),
                                                                m.eddy_diffusivities.Cᴷe⁰,
                                                                m.eddy_diffusivities.Cᴷeᵟ,
                                                                m.eddy_diffusivities.Cᴷeᶜ,
                                                                m.eddy_diffusivities.Cᴷeʷ,
                                                               )


text:
import IntervalArithmetic: diam, isinterior

export branch_and_prune, Bisection, Newton

diam(x::Root) = diam(x.interval)

Base.size(x::Interval) = (1,)

isinterior{N}(X::IntervalBox{N}, Y::IntervalBox{N}) = all(isinterior.(X, Y))

"""
    branch_and_prune(f, X, contractor, tol=1e-3)

Generic branch and prune routine for finding isolated roots of a function
`f:R^n → R^n` in a box.

Inputs:
- `f`: function whose roots will be found
- `X`: `Interval` or `IntervalBox`
- `contractor`: function that, when applied to the function `f`, determines
    the status of a given box `X`. It returns the new box and a symbol indicating
    the status. Current possible values are `Bisection` and `Newton`.

"""
function branch_and_prune(f, X, contractor, tol=1e-3)

    input_dim = length(X)
    output_dim = length(X)

    # @show input_dim
    # @show output_dim

    # if !(input_dim == output_dim)
    #     throw(ArgumentError("Input dimension ($input_dim) and output dimension ($output_dim) must be the same."))
    # end

    contract = contractor(Val{input_dim}, f)

    # main algorithm:

    working = [X]
    outputs = Root{typeof(X)}[]

    sizehint!(outputs, 100)
    sizehint!(working, 1000)

    while !isempty(working)
        # @show working
        X = pop!(working)

        status, output = contract(X)

        if status == :empty
            continue

        elseif status == :unique
            push!(outputs, Root(output, :unique))

        elseif diam(output) < tol
            push!(outputs, Root(output, :unknown))

        else  # branch
            X1, X2 = bisect(X)

            push!(working, X1, X2)
        end

    end

    return outputs
end

branch_and_prune(f, X::Root, contractor, tol=1e-3) =
    branch_and_prune(f, X.interval, contractor, tol)

function branch_and_prune(f, V::Vector{Root{T}}, contractor, tol=1e-3) where {T}
    reduce(append!, Root{T}[], [branch_and_prune(f, X.interval, contractor, tol) for X in V])
end

function branch_and_prune(f, V::Vector{T}, contractor, tol=1e-3) where {T}
    reduce(append!, Root{T}[], [branch_and_prune(f, X, contractor, tol) for X in V])
end

export recursively_branch_and_prune

function recursively_branch_and_prune(h, X, contractor=BisectionContractor, final_tol=1e-14)
    tol = 2
    roots = branch_and_prune(h, X, IntervalRootFinding.BisectionContractor, tol)

    while tol > 1e-14
       tol /= 2
       roots = branch_and_prune(h, roots, IntervalRootFinding.BisectionContractor, tol)
    end

    return roots
end

"""
If the input interval is complex, treat `f` as a complex function, currently of one complex variable `z`.
"""
function branch_and_prune{T}(f, Xc::Complex{Interval{T}}, contractor, tol=1e-3)

    g = realify(f)
    Y = IntervalBox(reim(Xc))

    roots = branch_and_prune(g, Y, contractor, tol)

    # @show roots

    return [Root(Complex(root.interval...), root.status) for root in roots]
end

contains_zero{T}(X::Interval{T}) = zero(T) ∈ X
contains_zero(X::SVector) = all(contains_zero.(X))
contains_zero(X::IntervalBox) = all(contains_zero(X[i]) for i in 1:length(X))

# contractors:

abstract type Contractor{F} end

export Bisection, Newton

struct Bisection{F} <: Contractor{F}
    dimension::Int
    f::F
end

Bisection{n}(::Type{Val{n}}, f) = Bisection(n, f)

function (contractor::Bisection)(X)
    image = contractor.f(X)

    if !(contains_zero(image))
        return :empty, X
    end

    return :unknown, X
end

"""
Generic refine operation for Krawczyk and Newton.
This assumes that it is already known that `X` contains a unique root.
Call using e.g. `op = X -> N(f, f_prime, X)`
"""
function refine(op, X)

    tolerance = 1e-16

    while diam(X) > tolerance  # avoid problem with tiny floating-point numbers if 0 is a root
        NX = op(X) ∩ X
        NX == X && break  # reached limit of precision
        X = NX
    end

    return X
end

struct Newton{F,FP,O} <: Contractor{F}
    dimension::Int
    f::F
    fp::FP
    op::O
end

function Newton(::Type{Val{1}}, f::Function)
    f_prime = x -> ForwardDiff.derivative(f, x)
    Newton(1, f, f_prime, N)
end

function Newton{n}(::Type{Val{n}}, f::Function)
    f_prime = x -> ForwardDiff.jacobian(f, x)
    Newton(n, f, f_prime, N)
end

function (C::Newton)(X)

    # use Bisection contractor for this:
    if !(contains_zero(IntervalBox(C.f(X))))
        return :empty, X
    end

    # given that have the Jacobian, can also do mean value form

    NX = C.op(C.f, C.fp, X) ∩ X

    isempty(NX) && return :empty, X

    if NX ⪽ X  # isinterior; know there's a unique root inside
        NX =  refine(X -> C.op(C.f, C.fp, X), NX)
        return :unique, NX
    end

    return :unknown, NX
end

"""
    roots(f, X, contractor, tol=1e-3)

Generic branch and prune routine for finding isolated roots of a function
`f:R^n → R^n` in a box `X`, or a vector of boxes.

Inputs:
- `f`: function whose roots will be found
- `X`: `Interval` or `IntervalBox`
- `contractor`: function that, when applied to the function `f`, determines
    the status of a given box `X`. It returns the new box and a symbol indicating
    the status. Current possible values are `Bisection` and `Newton`.

"""
roots{C<:Contractor}(f, X, contractor::Type{C}, tol::Float64=1e-3) = branch_and_prune(f, X, contractor, tol)

roots(f, X, tol::Float64=1e-3) = branch_and_prune(f, X, Newton, tol)


text:
@testset "Symmetry" begin
    words = collect(keys(SemanticSimilarity.CACHE.synsets))

    Random.seed!(42)
    for i in 1:1000
        # Similarity is symmetric
        word1 = rand(words)
        word2 = rand(words)
        @test similarity(word1, word2) ≈ similarity(word2, word1)
        @test similarity(WuPalmer(), word1, word2) ≈ similarity(WuPalmer(), word2, word1)
        @test similarity(SimilarityDepth(), word1, word2) ≈ similarity(SimilarityDepth(), word2, word1)
    end
    for i in 1:100
        # No word is similar to a nonexistent word
        word1 = rand(words)
        word2 = "abcdefgh"
        @test similarity(word1, word2) == similarity(word2, word1) == 0
        @test similarity(WuPalmer(), word1, word2) == similarity(WuPalmer(), word2, word1) == 0
        @test similarity(SimilarityDepth(), word1, word2) == similarity(SimilarityDepth(), word2, word1) == 0
    end

    # Nonexistent words are still identitcal to themselves
    @test similarity("abcdefgh", "abcdefgh") == 1.0
    @test similarity(WuPalmer(), "abcdefgh", "abcdefgh") == 1.0
    @test similarity(SimilarityDepth(), "abcdefgh", "abcdefgh") == 1.0

    # Nonexistent words are not similar to other nonexistent words
    @test similarity("abcdefgh", "asdfghj") == 0
end

@testset "Synonyms" begin
    # Members of the same synset should always have similarity == 1
    Random.seed!(1)
    synsets = collect(values(SemanticSimilarity.CACHE.synsets))
    for i in 1:100
        synset = rand(rand(synsets))
        for word1 in words(synset)
            for word2 in words(synset)
                @test similarity(CrypticCrosswords.normalize(word1),
                                 CrypticCrosswords.normalize(word2)) == 1.0
            end
        end
    end
end

@testset "Specific words" begin

≈
end


text:
struct HermiteTensorProductBasis{dim,NF,T} <: AbstractBasis{dim,NF}
    basis::CubicHermiteBasis{T}
    points::Matrix{T}
    function HermiteTensorProductBasis(dim; start = -1.0, stop = 1.0)
        basis1d = CubicHermiteBasis(start = start, stop = stop)
        NF = (number_of_basis_functions(basis1d))^dim
        points = hermite_tensor_product_interpolation_points(dim, start, stop)
        T = type_of_interpolation_points(basis1d)
        new{dim,NF,T}(basis1d, points)
    end
end

function Base.show(io::IO, basis::HermiteTensorProductBasis{dim,NF,T}) where {dim,NF,T}
    str = "HermiteTensorProductBasis\n\tDimension : $dim"
    print(io,str)
end

function order(B::HermiteTensorProductBasis{dim,NF,T}) where {dim,NF,T}
    return order(B.basis)
end

function hermite_tensor_product_interpolation_points(dim, start, stop)
    if dim == 2
        points = tensor_product_points(2, [start, stop]')
        return points
    else
        error("Current support for dim = $dim, got dim = $dim")
    end
end

function hermite_outer_product(u,v)
    return [
        u[1] * v[1],
        u[1] * v[2],
        u[2] * v[1],
        u[2] * v[2],
        u[1] * v[3],
        u[1] * v[4],
        u[2] * v[3],
        u[2] * v[4],
        u[3] * v[1],
        u[3] * v[2],
        u[4] * v[1],
        u[4] * v[2],
        u[3] * v[3],
        u[3] * v[4],
        u[4] * v[3],
        u[4] * v[4],
    ]
end

function (B::HermiteTensorProductBasis{2})(x, y)
    u = B.basis(x)
    v = B.basis(y)
    return hermite_outer_product(u,v)
end

function (B::HermiteTensorProductBasis{2})(x::V) where {V<:AbstractVector}
    @assert length(x) == 2
    return B(x[1],x[2])
end

function gradient(B::HermiteTensorProductBasis{2},dir::Z,x,y) where {Z<:Integer}
    @assert 1 <= dir <= 2
    if dir == 1
        dNx = derivative(B.basis,x)
        Ny = B.basis(y)
        return hermite_outer_product(dNx,Ny)
    else
        Nx = B.basis(x)
        dNy = derivative(B.basis,y)
        return hermite_outer_product(Nx,dNy)
    end
end

function gradient(B::HermiteTensorProductBasis{2},x,y)
    col1 = gradient(B,1,x,y)
    col2 = gradient(B,2,x,y)
    return hcat(col1,col2)
end

function hessian(B::HermiteTensorProductBasis{2},x,y)

    d2Nx = hessian(B.basis,x)
    d1Nx = derivative(B.basis,x)
    d0Nx = B.basis(x)

    d2Ny = hessian(B.basis,y)
    d1Ny = derivative(B.basis,y)
    d0Ny = B.basis(y)

    Nxx = hermite_outer_product(d2Nx,d0Ny)
    Nxy = hermite_outer_product(d1Nx,d1Ny)
    Nyy = hermite_outer_product(d0Nx,d2Ny)

    return hcat(Nxx,Nxy,Nyy)
end

function hessian(B::HermiteTensorProductBasis{2},x::V) where {V<:AbstractVector}
    @assert length(x) == 2
    return hessian(B,x[1],x[2])
end


text:
function main()

    call_your_code()

end

function call_your_code()

    l::Int64 = 3
    arrx = Vector{Float64}([5,5,5])
    arry = Vector{Float64}([2,3,7])
    s::Float64 = 5

    # x = ccall((:__mod_julfort_MOD_dot, "./mod_julfort.so"),
    #           Float64,
    #           (Ref{Int64}, Ptr{Float64}, Ptr{Float64}),
    #           l, arrx, arry)
    # println(x)

    println(arrx)
    ccall((:mult_arr, "./mod_julfort.so"),
          Cvoid,
          (Ref{Int64}, Ptr{Float64}, Ref{Float64}),
          l, arrx, s)
    println(arrx)

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4
coefs = [0.75;0.5]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1,1, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst1 = NewStencilInstance(st_def, m_step=false)

## Input Data size Definition
radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)
data = CreateData(dx,dy,dz)

t_steps = 16
## Compile Functions

gpu_out = ApplyFFTstencil(st_inst1, data, t_steps)
CUDA.cuProfilerStart()
NVTX.@range "r1" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,1) end
NVTX.@range "r2" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,2) end
NVTX.@range "r3" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,4) end
NVTX.@range "r4" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,8) end
NVTX.@range "r5" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,16) end
exit()


text:
function filter_approx!(cur_vector; atol=5e-2, rtol=eps())
  delete_indices = []

  for (cur_index, cur_value) in enumerate(cur_vector)
    is_duplicate = any(
      tmp_value -> loose_isapprox(tmp_value, cur_value, atol, rtol),
      cur_vector[1:cur_index-1]
    )

    is_duplicate || continue
    push!(delete_indices, cur_index)
  end

  deleteat!(cur_vector, delete_indices)

  cur_vector
end

function approx_push!(cur_vector, cur_value; atol=5e-2)
  is_duplicate = any(
    tmp_value -> loose_isapprox(tmp_value, cur_value, atol, rtol),
    cur_vector
  )

  is_duplicate && return cur_vector
  push!(cur_vector, cur_value)
end

function approx_append!(cur_vector, other_vector; atol=5e-2)
  append!(cur_vector, other_vector)
  filter_approx!(cur_vector; atol=atol)

  cur_vector
end


text:
###############################################################################
#
# Temperature Dedenpdency Parameter sets
# Data source: Bernacchi et al. (2001)
# Improved temperature response functions for models of Rubisco‐limited
# photosynthesis
#
###############################################################################
""" [`ArrheniusTD`](@ref) type Kc TD from Bernacchi's data """
KcTDBernacchi(FT) =
    ArrheniusTD{FT}(41.0264925, 79430.0 / GAS_R(), 79430.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Ko TD from Bernacchi's data """
KoTDBernacchi(FT) =
    ArrheniusTD{FT}(  28208.88, 36380.0 / GAS_R(), 36380.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Respiration TD from Bernacchi's data """
RespirationTDBernacchi(FT) =
    ArrheniusTD{FT}(       Inf, 46390.0 / GAS_R(), 46390.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Vcmax TD from Bernacchi's data """
VcmaxTDBernacchi(FT) =
    ArrheniusTD{FT}(       Inf, 65330.0 / GAS_R(), 65330.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Vomax TD from Bernacchi's data """
VomaxTDBernacchi(FT) =
    ArrheniusTD{FT}(       Inf, 60110.0 / GAS_R(), 60110.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type ``Γ^{*}`` TD from Bernacchi's data """
ΓStarTDBernacchi(FT) =
    ArrheniusTD{FT}(4.33164375, 37830.0 / GAS_R(), 37830.0 / RT_25(FT));

###############################################################################
#
# Temperature Dedenpdency Parameter sets
# Data source: Boyd et al. (2001)
# Temperature responses of C4 photosynthesis: biochemical analysis of Rubisco,
# phosphoenolpyruvate carboxylase, and carbonic anhydrase in Setaria viridis
#
###############################################################################
""" [`ArrheniusTD`](@ref) type Kpep TD from Boyd's data """
KpepTDBoyd(FT) =
    ArrheniusTD{FT}(16.0, 36300.0 / GAS_R(), 36300.0 / RT_25(FT));

""" [`ArrheniusPeakTD`](@ref) type Vpmax TD from Boyd's data """
function VpmaxTDBoyd(FT)
    ΔHa_to_RT25::FT = 94800.0 / RT_25(FT);
    ΔHd_to_R::FT = 73300.0 / GAS_R();
    ΔSv_to_R::FT = 250.0   / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

###############################################################################
#
# Temperature Dedenpdency Parameter sets
# Data source: Leuning (2002)
# Temperature dependence of two parameters in a photosynthesis model
#
###############################################################################
""" [`ArrheniusPeakTD`](@ref) type Jmax TD from Leuning's data """
function JmaxTDLeuning(FT)
    ΔHa_to_RT25::FT = 50300.0  / RT_25(FT);
    ΔHd_to_R::FT = 152044.0 / GAS_R();
    ΔSv_to_R::FT = 495.0    / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

""" [`ArrheniusPeakTD`](@ref) type Vcmax TD from Leuning's data """
function VcmaxTDLeuning(FT)
    ΔHa_to_RT25::FT = 73637.0  / RT_25(FT);
    ΔHd_to_R::FT = 149252.0 / GAS_R();
    ΔSv_to_R::FT = 486.0    / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

###############################################################################
#
# Temperature Dedenpdency Parameter sets
# Data source: Lavigne and Ryan (1997)
# Growth and maintenance respiration rates of aspen, blackspruce and jack pine
#     stems at northern and southern BOREAS sites
#
###############################################################################
""" [`Q10TD`](@ref) type Respiration TD for angiosperms per biomass """
Q10TDAngiosperm(FT) = Q10TD{FT}(0.014/8760, 298.15, 1.4)

""" [`Q10TD`](@ref) type Respiration TD for symnosperms per biomass """
Q10TDGymnosperm(FT) = Q10TD{FT}(0.0425/8760, 298.15, 1.7)

###############################################################################
#
# Temperature Dedenpdency Parameter sets
# Data source missing
#
###############################################################################
""" [`ArrheniusTD`](@ref) type Kc TD """
KcTDCLM(FT) =
    ArrheniusTD{FT}(  40.49, 79430.0 / GAS_R(), 79430.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Ko TD """
KoTDCLM(FT) =
    ArrheniusTD{FT}(27840.0, 36380.0 / GAS_R(), 36380.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Kpep TD """
KpepTDCLM(FT) =
    ArrheniusTD{FT}(    8.0, 36000.0 / GAS_R(), 36000.0 / RT_25(FT));

""" [`ArrheniusTD`](@ref) type Γ* TD """
ΓStarTDCLM(FT) =
    ArrheniusTD{FT}(  4.275, 37830.0 / GAS_R(), 37830.0 / RT_25(FT));

""" [`ArrheniusPeakTD`](@ref) type Jmax TD """
function JmaxTDBernacchi(FT)
    ΔHa_to_RT25::FT = 57500.0  / RT_25(FT);
    ΔHd_to_R::FT = 439000.0 / GAS_R();
    ΔSv_to_R::FT = 1400.0   / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

""" [`ArrheniusPeakTD`](@ref) type Jmax TD """
function JmaxTDCLM(FT)
    ΔHa_to_RT25::FT = 43540.0  / RT_25(FT);
    ΔHd_to_R::FT = 150000.0 / GAS_R();
    ΔSv_to_R::FT = 490.0    / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

""" [`ArrheniusPeakTD`](@ref) type Respiration TD """
function RespirationTDCLM(FT)
    ΔHa_to_RT25::FT = 46390.0  / RT_25(FT);
    ΔHd_to_R::FT = 150650.0 / GAS_R();
    ΔSv_to_R::FT = 490.0    / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

""" [`ArrheniusPeakTD`](@ref) type Vcmax TD """
function VcmaxTDCLM(FT)
    ΔHa_to_RT25::FT = 65330.0  / RT_25(FT);
    ΔHd_to_R::FT = 150000.0 / GAS_R();
    ΔSv_to_R::FT = 490.0    / GAS_R();
    C::FT = 1 + exp( ΔSv_to_R - ΔHd_to_R/T_25(FT) );
    return ArrheniusPeakTD{FT}(ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R, C)
end

###############################################################################
#
# Vcmax to R correlation
#
###############################################################################
""" A constant of 0.01 """
VtoRCollatz(FT) = FT(0.010);

""" A constant of 0.015 """
VtoRDefault(FT) = FT(0.015);

###############################################################################
#
# Fluorescence model parameter set
# Data source: van der Tol et al. (2014)
# Models of fluorescence and photosynthesis for interpreting measurements of
#     solar-induced chlorophyll fluorescence
#
###############################################################################
""" [`FluoParaSet`](@ref) type parameter set using all data """
FluorescenceVanDerTol(FT) =
        FluoParaSet{FT}(2.48, 2.83, 0.114);

""" [`FluoParaSet`](@ref) type parameter set using Flexas's data (drought) """
FluorescenceVanDerTolDrought(FT) =
        FluoParaSet{FT}(5.01, 1.93, 10.0);

###############################################################################
#
# Pre-set photosynthesis model parameter sets
#
###############################################################################
""" [`C3ParaSet`](@ref) type C3 photosynthesis using Bernacchi's data """
function C3Bernacchi(FT)
    JT  = JmaxTDBernacchi(FT);
    KcT = KcTDBernacchi(FT);
    KoT = KoTDBernacchi(FT);
    ReT = RespirationTDBernacchi(FT);
    VcT = VcmaxTDBernacchi(FT);
    ΓsT = ΓStarTDBernacchi(FT);
    Flu = FluorescenceVanDerTolDrought(FT);
    VR  = VtoRDefault(FT);
    E1  = FT(4);
    E2  = FT(8);
    return C3ParaSet{FT}(JT, KcT, KoT, ReT, VcT, ΓsT, Flu, VR, E1, E2)
end

""" [`C3ParaSet`](@ref) type C3 photosynthesis using CLM5's data """
function C3CLM(FT)
    JT  = JmaxTDCLM(FT);
    KcT = KcTDCLM(FT);
    KoT = KoTDCLM(FT);
    ReT = RespirationTDCLM(FT);
    VcT = VcmaxTDCLM(FT);
    ΓsT = ΓStarTDCLM(FT);
    Flu = FluorescenceVanDerTolDrought(FT);
    VR  = VtoRDefault(FT);
    E1  = FT(4);
    E2  = FT(8);
    return C3ParaSet{FT}(JT, KcT, KoT, ReT, VcT, ΓsT, Flu, VR, E1, E2)
end

""" [`C4ParaSet`](@ref) type C4 photosynthesis using CLM5's data """
function C4CLM(FT)
    KpT = KpepTDCLM(FT);
    ReT = RespirationTDCLM(FT);
    VcT = VcmaxTDCLM(FT);
    VpT = VpmaxTDBoyd(FT);
    Flu = FluorescenceVanDerTolDrought(FT);
    VR  = VtoRDefault(FT);
    return C4ParaSet{FT}(KpT, ReT, VcT, VpT, Flu, VR)
end


text:
using LiterateTest
@evaltest_throw ":no_throw" begin
end