text:
using CovarianceMatrices
using Base.Test

############################################################
## HAC
############################################################

X = randn(100, 5);

@time vcov(X, TruncatedKernel(2.))
@time vcov(X, BartlettKernel(2.))
@time vcov(X, ParzenKernel(2.))
@time vcov(X, QuadraticSpectralKernel(2.))
@time vcov(X, TruncatedKernel())
@time vcov(X, BartlettKernel())
@time vcov(X, ParzenKernel())
@time vcov(X, QuadraticSpectralKernel())

############################################################
## HC
############################################################

# A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
clotting = DataFrame(
    u    = log([5,10,15,20,30,40,60,80,100]),
    lot1 = [118,58,42,35,27,25,21,19,18],
    lot2 = [69,35,26,21,18,16,13,12,12],
    w    = [1/8, 1/9, 1/25, 1/6, 1/14, 1/25, 1/15, 1/13, 0.3022039]
)

## Unweighted OLS though GLM interface
OLS = fit(GeneralizedLinearModel, lot1~u,clotting, Normal(), IdentityLink())
S0 = vcov(OLS, HC0())
S1 = vcov(OLS, HC1())
S2 = vcov(OLS, HC2())
S3 = vcov(OLS, HC3())
S4 = vcov(OLS, HC4())
S4m = vcov(OLS, HC4m())
S5 = vcov(OLS, HC5())

St0 = [720.6213 -190.0645; -190.0645 51.16333]
St1 = [926.5131 -244.3687; -244.3687 65.78143]
St2 = [1300.896 -343.3307; -343.3307 91.99719]
St3 = [2384.504 -628.975 ; -628.975 167.7898]
St4 = [2538.746 -667.9597; -667.9597 177.2631]
St4m= [3221.095 -849.648 ; -849.648 226.1705]
St5 = St4

@test abs(maximum(S0 .- St0)) < 1e-04
@test abs(maximum(S1 .- St1)) < 1e-04
@test abs(maximum(S2 .- St2)) < 1e-04
@test abs(maximum(S3 .- St3)) < 1e-04
@test abs(maximum(S4 .- St4)) < 1e-03
@test abs(maximum(S4m .- St4m)) < 1e-03
@test abs(maximum(S5 .- St5)) < 1e-03

M0 = meat(OLS.model, HC0())
M1 = meat(OLS.model, HC1())
M2 = meat(OLS.model, HC2())
M3 = meat(OLS.model, HC3())
M4 = meat(OLS.model, HC4())
M4m = meat(OLS.model, HC4m())
M5 = meat(OLS.model, HC5())

Mt0 = [206.6103 518.7871; 518.7871 1531.173]
Mt1 = [265.6418 667.012;  667.012 1968.651]
Mt2 = [323.993 763.0424;  763.0424 2149.767]
Mt3 = [531.478 1172.751;  1172.751 3122.79]
Mt4 = [531.1047 1110.762; 1110.762 2783.269]
Mt4m= [669.8647 1412.227; 1412.227 3603.247]
Mt5 = Mt4

@test abs(maximum(M0 .- Mt0)) < 1e-04
@test abs(maximum(M1 .- Mt1)) < 1e-04
@test abs(maximum(M2 .- Mt2)) < 1e-04
@test abs(maximum(M3 .- Mt3)) < 1e-03
@test abs(maximum(M4 .- Mt4)) < 1e-03
@test abs(maximum(M4m .- Mt4m)) < 1e-03
@test abs(maximum(M5 .- Mt5)) < 1e-03

## Unweighted
OLS = glm(lot1~u,clotting, Normal(), IdentityLink())
S0 = vcov(OLS, HC0())
S1 = vcov(OLS, HC1())
S2 = vcov(OLS, HC2())
S3 = vcov(OLS, HC3())
S4 = vcov(OLS, HC4())

## Weighted OLS though GLM interface
wOLS = fit(GeneralizedLinearModel, lot1~u, clotting, Normal(),
           IdentityLink(), wts = Vector{Float64}(clotting[:w]))
S0 = vcov(wOLS, HC0())
S1 = vcov(wOLS, HC1())
S2 = vcov(wOLS, HC2())
S3 = vcov(wOLS, HC3())
S4 = vcov(wOLS, HC4())
S4m= vcov(wOLS, HC4m())
S5 = vcov(wOLS, HC5())

St0 = [717.7362 -178.4043; -178.4043 45.82273]
St1 = [922.8037 -229.3769; -229.3769 58.91494]
St2 = [1412.94 -361.33; -361.33 95.91252]
St3 = [2869.531 -756.2976; -756.2976 208.2344]
St4 = [3969.913 -1131.358; -1131.358 342.2859]
St4m= [4111.626 -1103.174; -1103.174   310.194]
St5 = St4

@test abs(maximum(S0 .- St0)) < 1e-04
@test abs(maximum(S1 .- St1)) < 1e-04
@test abs(maximum(S2 .- St2)) < 1e-03
@test abs(maximum(S3 .- St3)) < 1e-03
@test abs(maximum(S4 .- St4)) < 1e-03
@test abs(maximum(S4m .- St4m)) < 1e-03
@test abs(maximum(S5 .- St5)) < 1e-03

## Unweighted GLM - Gamma
GAMMA = glm(lot1~u, clotting, Gamma(),InverseLink())

S0 = vcov(GAMMA, HC0())
S1 = vcov(GAMMA, HC1())
S2 = vcov(GAMMA, HC2())
S3 = vcov(GAMMA, HC3())
S4 = vcov(GAMMA, HC4())
S4m = vcov(GAMMA, HC4m())
S5 = vcov(GAMMA, HC5())

St0 = [4.504287921232951e-07 -1.700020601541489e-07;
       -1.700020601541490e-07  8.203697048568913e-08]

St1 = [5.791227327299548e-07 -2.185740773410504e-07;
       -2.185740773410510e-07  1.054761049101728e-07]

St2 = [3.192633083111232e-06 -9.942484630848573e-07;
       -9.942484630848578e-07  3.329973305723091e-07]

St3 = [2.982697811926944e-05 -8.948137019946751e-06;
       -8.948137019946738e-06  2.712024459305714e-06]

St4 = [0.002840158946368653 -0.0008474436578800430;
       -0.000847443657880045  0.0002528819761961959]

St4m= [8.49306e-05 -2.43618e-05; -2.43618e-05  7.04210e-06]

St5 = St4

@test abs(maximum(S0 .- St0)) < 1e-06
@test abs(maximum(S1 .- St1)) < 1e-06
@test abs(maximum(S2 .- St2)) < 1e-06
@test abs(maximum(S3 .- St3)) < 1e-06
@test abs(maximum(S4 .- St4)) < 1e-06
@test abs(maximum(S4m .- St4m)) < 1e-05
@test abs(maximum(S5 .- St5)) < 1e-05

## Weighted Gamma

GAMMA = glm(lot1~u, clotting, Gamma(),InverseLink(), wts = convert(Array, clotting[:w]))

S0 = vcov(GAMMA, HC0())
S1 = vcov(GAMMA, HC1())
S2 = vcov(GAMMA, HC2())
S3 = vcov(GAMMA, HC3())
S4 = vcov(GAMMA, HC4())
S4m = vcov(GAMMA, HC4m())
S5 = vcov(GAMMA, HC5())

St0 = [4.015104e-07 -1.615094e-07;
       -1.615094e-07  8.378363e-08]

St1 = [5.162277e-07 -2.076549e-07;
       -2.076549e-07  1.077218e-07]

St2 = [2.720127e-06 -8.490977e-07;
       -8.490977e-07  2.963563e-07]

St3 = [2.638128e-05 -7.639883e-06;
       -7.639883e-06  2.259590e-06]

St4 = [0.0029025754 -0.0008275858;
       -0.0008275858  0.0002360053]

St4m = [8.493064e-05 -2.436180e-05; -2.436180e-05  7.042101e-06]

St5 = St4

@test abs(maximum(S0 .- St0)) < 1e-06
@test abs(maximum(S1 .- St1)) < 1e-06
@test abs(maximum(S2 .- St2)) < 1e-06
@test abs(maximum(S3 .- St3)) < 1e-06
@test abs(maximum(S4 .- St4)) < 1e-06
@test abs(maximum(S4m .- St4m)) < 1e-05
@test abs(maximum(S5 .- St5)) < 1e-05

### Cluster basic interface

## Construct Fake ols data

## srand(1)

## df = DataFrame( Y = randn(500),
##                X1 = randn(500),
##                X2 = randn(500),
##                X3 = randn(500),
##                X4 = randn(500),
##                X5 = randn(500),
##                w  = rand(500),
##                cl = repmat(collect(1:25), 20))

df = readtable("wols_test.csv")

OLS = fit(GeneralizedLinearModel, Y~X1+X2+X3+X4+X5, df,
          Normal(), IdentityLink())

S0 = vcov(OLS, HC0())
S1 = vcov(OLS, HC1())
S2 = vcov(OLS, HC2())
S3 = vcov(OLS, HC3())
S4 = vcov(OLS, HC4())
S5 = vcov(OLS, HC5())

cl = convert(Array, df[:cl])
S0 = stderr(OLS, CRHC0(cl))
S1 = stderr(OLS, CRHC1(cl))
S2 = stderr(OLS, CRHC2(cl))
S3 = stderr(OLS, CRHC3(cl))

## STATA
St1 = [.0374668, .0497666, .0472636, .0437952, .0513613, .0435369]

@test maximum(abs(S0 .- St1)) < 1e-02
@test maximum(abs(S1 .- St1)) < 1e-04
@test maximum(abs(S2 .- St1)) < 1e-02
@test maximum(abs(S3 .- St1)) < 1e-02

wOLS = fit(GeneralizedLinearModel, Y~X1+X2+X3+X4+X5, df,
          Normal(), IdentityLink(), wts = convert(Array, df[:w]))

S0 = stderr(wOLS, CRHC0(cl))
S1 = stderr(wOLS, CRHC1(cl))
S2 = stderr(wOLS, CRHC2(cl))
S3 = stderr(wOLS, CRHC3(cl))

St1 = [0.042839848169137905,0.04927285387211425,
       0.05229519531359171,0.041417170723876025,
       0.04748115282615204,0.04758615959662984]

@test maximum(abs(S1 .- St1)) < 1e-10

############################################################
## Linear Model
############################################################
srand(1)
y = randn(100);
x = randn(100, 5);

lm1 = lm(x, y)
@test_approx_eq_eps stderr(lm1, HC0()) [0.0941998, 0.0946132, 0.0961678, 0.0960445, 0.101651] 1e-06
@test_approx_eq_eps diag(vcov(lm1, HC0())) [0.0941998, 0.0946132, 0.0961678, 0.0960445, 0.101651].^2 1e-06

############################################################
## HAC
############################################################

include("ols_hac.jl")


text:
__precompile__()

module PandaRobot

using RigidBodyDynamics
using StaticArrays

export Panda,
  packagepath,
  urdfpath

packagepath() = joinpath(@__DIR__, "..", "deps")
urdfpath() = joinpath(packagepath(), "Panda", "panda.urdf")

mutable struct Panda{T}
  mechanism::Mechanism{T}

  function Panda{T}() where T
    mechanism = RigidBodyDynamics.parse_urdf(T,urdfpath())
    new{T}(mechanism)
  end
end
Panda(::Type{T} = Float64; kwargs...) where {T} = Panda{T}(; kwargs...)

function __init__()
  if !isfile(urdfpath())
    error("Could not find $(urdfpath()). Please run `import Pkg; Pkg.build(\"Panda\")`.")
  end
end

end # module

ext: jl


text:
@testset "Constants" begin
    @test BibInternal.maxfieldlength == 13
    @info "Collections: " BibInternal.fields BibInternal.entries
    @test BibInternal.Fields == Dict{String,String}
end


text:
module GroupRings

using AbstractAlgebra
import AbstractAlgebra: Group, NCRing, NCRingElem, parent, elem_type, parent_type, addeq!, mul!

using SparseArrays
import LinearAlgebra

import Base: convert, show, hash, ==, +, -, *, ^, //, /, length, getindex, setindex!, eltype, one, zero

GroupOrNCRing = Union{AbstractAlgebra.Group, AbstractAlgebra.NCRing}
GroupOrNCRingElem = Union{AbstractAlgebra.GroupElem, AbstractAlgebra.NCRingElem}

###############################################################################
#
#   GroupRings / GroupRingsElem
#
###############################################################################

mutable struct GroupRing{Gr<:GroupOrNCRing, T<:GroupOrNCRingElem} <: NCRing
   group::Gr
   basis::Vector{T}
   basis_dict::Dict{T, Int}
   pm::Matrix{Int}

   function GroupRing(G::Gr, basis::Vector{T};
         cachedmul::Bool=false) where {Gr, T}
      RG = new{Gr, T}(G, basis, reverse_dict(basis))
      cachedmul && initializepm!(RG)
      return RG
   end

   function GroupRing(G::Gr, b::Vector{T}, b_d::Dict{T, Int}, pm::Array{Int,2}) where {Gr,T}
      return new{Gr, T}(G, b, b_d, pm)
   end

   function GroupRing(G::Gr, pm::Array{Int,2}) where {Gr}
      RG = new{Gr, elem_type(G)}(G)
      RG.pm = pm
      return RG
   end
end

mutable struct GroupRingElem{T, A<:AbstractVector, GR<:GroupRing} <: NCRingElem
   coeffs::A
   parent::GR

   function GroupRingElem{T, A, GR}(c::AbstractVector{T}, RG::GR, check=true) where {T, A, GR}
      if check
         if hasbasis(RG)
            length(c) == length(RG.basis) || throw(
            "Can't create GroupRingElem -- lengths differ: length(c) =
            $(length(c)) != $(length(RG.basis)) = length(RG.basis)")
         else
            # @warn("Basis of the GroupRing is not defined.")
         end
      end
      return new{T, A, GR}(c, RG)
   end
end

export GroupRing, GroupRingElem, complete!, create_pm, star, aug, supp

###############################################################################
#
#   GroupRing / GroupRingElem constructors
#
###############################################################################

function GroupRingElem(c::AbstractVector, RG::GroupRing)
   return GroupRingElem{eltype(c), typeof(c), typeof(RG)}(c, RG)
end

function GroupRing(G::Generic.SymmetricGroup; cachedmul::Bool=false)
   return GroupRing(G, vec(collect(G)), cachedmul=cachedmul)
end

function GroupRing(G::Group, basis::AbstractVector, pm::AbstractMatrix{<:Integer})
   size(pm,1) == size(pm,2) || throw("pm must be square, got $(size(pm))")
   eltype(basis) == elem_type(G) || throw("Basis must consist of elements of $G")
   return GroupRing(G, basis, reverse_dict(basis), pm)
end

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

elem_type(::Type{GroupRing}) = GroupRingElem

eltype(::Type{GroupRingElem{T, A, Gr}}) where {T, A, Gr} = T

parent(g::GroupRingElem) = g.parent

parent_type(X::GroupRingElem) = typeof(parent(X))

Base.promote_rule(::Type{GroupRingElem{T}}, ::Type{GroupRingElem{S}}) where {T,S} =
   GroupRingElem{promote_type(T,S)}

function convert(::Type{T}, X::GroupRingElem) where T<:Number
   return GroupRingElem(Vector{T}(X.coeffs), parent(X))
end

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

# sparse storage:

zero(RG::GroupRing, T::Type=Int) = RG(T)
one(RG::GroupRing, T::Type=Int) = RG(one(RG.group), T)
one(RG::GroupRing{<:AbstractAlgebra.NCRing}, T::Type=Int) = RG(one(RG.group), T)

function (RG::GroupRing)(T::Type=Int)
   hasbasis(RG) || throw("Can not coerce without basis of GroupRing")
   return GroupRingElem(spzeros(T,length(RG.basis)), RG)
end

function (RG::GroupRing)(i::Int, T::Type=Int)
   elt = RG(T)
   elt[one(RG.group)] = i
   return elt
end

function (RG::GroupRing{<:AbstractAlgebra.NCRing})(i::Int, T::Type=Int)
   elt = RG(T)
   elt[one(RG.group)] = i
   return elt
end

function (RG::GroupRing)(g::GroupOrNCRingElem, T::Type=Int)
   result = RG(T)
   result[RG.group(g)] = one(T)
   return result
end

function (RG::GroupRing{Gr,T})(V::Vector{T}, S::Type=Int) where {Gr, T}
   res = RG(S)
   for g in V
      res[g] += one(S)
    end
    return res
end

function (RG::GroupRing)(f::Function, X::GroupRingElem{T}) where T
   hasbasis(RG) || throw("Can not coerce without basis of GroupRing")
   res = RG(T)
   for g in supp(X)
      res[f(g)] = X[g]
   end
   return res
end

# keep storage type

function (RG::GroupRing)(x::AbstractVector{T}) where T
   hasbasis(RG) || throw("Basis of GroupRing not defined. For advanced use the direct constructor of GroupRingElem is provided.")
   length(x) == length(RG.basis) || throw("Can not coerce to $RG: lengths differ")
   return GroupRingElem(x, RG)
end

function (RG::GroupRing)(X::GroupRingElem)
   RG == parent(X) || throw("Can not coerce!")
   return RG(X.coeffs)
end

###############################################################################
#
#   Basic manipulation && Array protocol
#
###############################################################################

function hash(X::GroupRingElem, h::UInt)
   return hash(X.coeffs, hash(parent(X), hash(GroupRingElem, h)))
end

function getindex(X::GroupRingElem, n::Int)
   return X.coeffs[n]
end

function getindex(X::GroupRingElem, g::GroupOrNCRingElem)
   return X.coeffs[parent(X).basis_dict[g]]
end

function setindex!(X::GroupRingElem, value, n::Int)
   X.coeffs[n] = value
end

function setindex!(X::GroupRingElem, value, g::GroupOrNCRingElem)
   RG = parent(X)
   if !(g in keys(RG.basis_dict))
      g = (RG.group)(g)
   end
   X.coeffs[RG.basis_dict[g]] = value
end

Base.size(X::GroupRingElem) = size(X.coeffs)
Base.IndexStyle(::Type{GroupRingElem}) = Base.LinearFast()

dense(X::GroupRingElem{T, A}) where {T, A<:DenseVector} = X

function dense(X::GroupRingElem{T, Sp}) where {T, Sp<:SparseVector}
   return GroupRingElem(Vector(X.coeffs), parent(X))
end

SparseArrays.sparse(X::GroupRingElem{T, Sp}) where {T, Sp<:SparseVector} = X

function SparseArrays.sparse(X::GroupRingElem{T, A}) where {T, A<:Vector}
   return GroupRingElem(sparse(X.coeffs), parent(X))
end

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, A::GroupRing)
   print(io, "Group Ring of $(A.group)")
end

function show(io::IO, X::GroupRingElem)
   RG = parent(X)
   T = eltype(X.coeffs)
   if X.coeffs == zero(X.coeffs)
      print(io, "$(zero(T))*$(one(RG.group))")
   elseif hasbasis(RG)
      non_zeros = ((X.coeffs[i], RG.basis[i]) for i in findall(!iszero, X.coeffs))
      elts = String[]
      for (c,g) in non_zeros
      	  sgn = (sign(c)>=0 ? " + " : " - ")
	  if c == T(1)
	      coeff = ""
	  else
	      coeff = "$(abs(c))"
	  end
	  push!(elts, sgn*coeff*"$(g)")
      end
      str = join(elts, "")
      if sign(first(non_zeros)[1]) > 0
         str = str[4:end]
      end
      print(io, str)
   else
      @warn("Basis of the parent Group is not defined, showing coeffs")
      show(io, MIME("text/plain"), X.coeffs)
   end
end

###############################################################################
#
#   Comparison
#
###############################################################################

function (==)(X::GroupRingElem, Y::GroupRingElem)
   suppX = supp(X)
   suppX == supp(Y) || return false

   for g in suppX
      X[g] == Y[g] || return false
   end

   return true
end

function (==)(A::GroupRing, B::GroupRing)
   A.group == B.group || return false
   if hasbasis(A) && hasbasis(B)
      A.basis == B.basis || return false
   elseif isdefined(A, :pm) && isdefined(B, :pm)
      A.pm == B.pm || return false
   end
   return true
end

hasbasis(A::GroupRing) = isdefined(A, :basis)

Base.deepcopy_internal(x::GroupRingElem, dict::IdDict) =
   GroupRingElem(deepcopy(x.coeffs), parent(x))

###############################################################################
#
#   Scalar operators
#
###############################################################################

(-)(X::GroupRingElem) = GroupRingElem(-X.coeffs, parent(X))

function mul!(a::T, X::GroupRingElem{T}) where T
   X.coeffs .*= a
   return X
end

_mul(a::Number, X::GroupRingElem) = GroupRingElem(a.*X.coeffs, parent(X))

Base.:*(a::Number, X::GroupRingElem) = _mul(a, X)
Base.:*(X::GroupRingElem, a::Number) = _mul(a, X)

# disallow Rings to hijack *(::, ::GroupRingElem)
*(a::Union{AbstractFloat, Integer, RingElem, Rational}, X::GroupRingElem) = _mul(a, X)

(/)(X::GroupRingElem, a) = 1/a*X
(//)(X::GroupRingElem, a::Union{Integer, Rational}) = 1//a*X

(^)(X::GroupRingElem, n::Integer) = Base.power_by_squaring(X, n)

###############################################################################
#
#   Binary operators
#
###############################################################################

function addeq!(X::GroupRingElem, Y::GroupRingElem)
   X.coeffs += Y.coeffs
   return X
end

Base.:+(X::GroupRingElem, Y::GroupRingElem) = addeq!(deepcopy(X), Y)
-(X::GroupRingElem, Y::GroupRingElem) where T = addeq!((-Y), X)

"""
    fmac!(result::AbstractVector{T},
              X::AbstractVector,
              Y::AbstractVector,
             pm::Array{Int,2}) where T
> Fused multiply-add for group ring coeffs using multiplication table `pm`.
> The result of X*Y in GroupRing is added in-place to `result`.
> Notes:
> * this method will silently produce false results if `X[k]` is non-zero for
> `k > size(pm,1)`.
> * This method will fail if any zeros (i.e. uninitialised entries) are present
> in `pm`.
> Use with extreme care!
"""

function fmac!(result::AbstractVector{T},
                   X::AbstractVector,
                   Y::AbstractVector,
                  pm::Array{Int,2}) where T
   z = zero(T)
   s1 = size(pm,1)
   s2 = size(pm,2)

   @inbounds for j in 1:s2
      if Y[j] != z
         for i in 1:s1
            if X[i] != z
               result[pm[i,j]] += X[i]*Y[j]
            end
         end
      end
   end
   return result
end

"""
    GRmul!(result::AbstractVector{T}, X::AbstractVector, Y::AbstractVector,
           pm::Matrix{<:Integer}) where T
> The most specialised multiplication for `X` and `Y` (intended for `coeffs` of
> `GroupRingElems`), using multiplication table `pm`.
> Notes:
> * this method will silently produce false results if `X[k]` is non-zero for
> `k > size(pm,1)`.
> * This method will fail if any zeros (i.e. uninitialised entries) are present
> in `pm`.
> Use with extreme care!
"""
function GRmul!(result::AbstractVector{T},
                   X::AbstractVector,
                   Y::AbstractVector,
                  pm::AbstractMatrix{<:Integer}) where T
   z = zero(T)
   result .= z

   return fmac!(result, X, Y, pm)
end

"""
    mul!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem)
> In-place multiplication for `GroupRingElem`s `X` and `Y`.
> `mul!` will make use the initialised entries of `pm` attribute of
> `parent(X)::GroupRing` (if available), and will compute and store in `pm` the
> remaining products necessary to perform the multiplication.
> The method will fail with `KeyError` if product `X*Y` is not supported on
> `parent(X).basis`.
"""
function mul!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem)
   if result === X
      result = deepcopy(result)
   end
   T = eltype(result.coeffs)
   z = zero(T)
   result.coeffs .= z

   RG = parent(X)

   lX = length(X.coeffs)
   lY = length(Y.coeffs)

   if isdefined(RG, :pm)
      s = size(RG.pm)
      k = findprev(!iszero, X.coeffs, lX)
      (k == nothing ? 0 : k) <= s[1] || throw("Element in X outside of support of parents product: $k $(RG.basis[k])")
      k = findprev(!iszero, Y.coeffs, lY)
      (k == nothing ? 0 : k) <= s[2] || throw("Element in Y outside of support of parents product")

      for j in 1:lY
         if Y.coeffs[j] != z
            for i in 1:lX
               if X.coeffs[i] != z
                  if RG.pm[i,j] == 0
                     RG.pm[i,j] = RG.basis_dict[RG.basis[i]*RG.basis[j]]
                  end
                  result.coeffs[RG.pm[i,j]] += X[i]*Y[j]
               end
            end
         end
      end
   else
      for j in 1:lY
         if Y.coeffs[j] != z
            for i in 1:lX
               if X.coeffs[i] != z
                  result[RG.basis[i]*RG.basis[j]] += X[i]*Y[j]
               end
            end
         end
      end
   end
   return result
end

function *(X::GroupRingElem{T}, Y::GroupRingElem{S}, check::Bool=true) where {T,S}
   if check
      parent(X) === parent(Y) || throw("Elements don't seem to belong to the same Group Ring!")
   end

   TT = promote_type(T,S)

   if hasbasis(parent(X))
      result = parent(X)(similar(X.coeffs, TT))
      result = mul!(result, X, Y)
   else
      result = similar(X.coeffs, TT)
      result = GRmul!(result, X.coeffs, Y.coeffs, parent(X).pm)
      result = GroupRingElem(result, parent(X))
   end
   return result
end

###############################################################################
#
#   *-involution
#
###############################################################################

function star(X::GroupRingElem{T}) where T
   RG = parent(X)
   hasbasis(RG) || throw("*-involution without basis is not possible")
   result = RG(T)
   for (i,c) in enumerate(X.coeffs)
      if c != zero(T)
         g = star(RG.basis[i])
         result[g] = c
      end
   end
   return result
end

star(g::GroupElem) = inv(g)
star(r::NCRingElem) = inv(r)

###############################################################################
#
#   Misc
#
###############################################################################

length(X::GroupRingElem) = count(!iszero, X.coeffs)

LinearAlgebra.norm(X::GroupRingElem, p::Int=2) = LinearAlgebra.norm(X.coeffs, p)

aug(X::GroupRingElem) = sum(X.coeffs)

supp(X::GroupRingElem) = parent(X).basis[findall(!iszero, X.coeffs)]

function reverse_dict(::Type{I}, iter) where I<:Integer
   length(iter) > typemax(I) && error("Can not produce reverse dict: $(length(iter)) is too large for $T")
   return Dict{eltype(iter), I}(x => i for (i,x) in enumerate(iter))
end

reverse_dict(iter) = reverse_dict(Int, iter)

function create_pm(basis::AbstractVector{T}, basis_dict::Dict{T, Int},
   limit::Int=length(basis); twisted::Bool=false, check=true) where T
   product_matrix = zeros(Int, (limit,limit))
   Threads.@threads for i in 1:limit
      x = basis[i]
      if twisted
         x = star(x)
      end
      for j in 1:limit
         product_matrix[i,j] = get(basis_dict, x*basis[j], 0)
      end
   end

   check && check_pm(product_matrix, basis, twisted)

   return product_matrix
end

create_pm(b::AbstractVector{<:GroupOrNCRingElem}) = create_pm(b, reverse_dict(b))

function check_pm(product_matrix, basis, twisted=false)
   idx = findfirst(product_matrix' .== 0)
   if idx != nothing
      @warn("Product is not supported on basis")
      i,j = Tuple(idx)
      x = (twisted ? star(basis[i]) : basis[i])
      throw(KeyError(x*basis[j]))
   end
   return true
end

function complete!(RG::GroupRing, twisted::Bool=false)
   hasbasis(RG) || throw(ArgumentError("Provide basis for completion first!"))
   if !isdefined(RG, :pm)
      initializepm!(RG, fill=false)
      return RG
   end

   warning = false
   for idx in findall(iszero, RG.pm)
      i,j = Tuple(idx)
      g = (twisted ? star(RG.basis[i]) : RG.basis[i])*RG.basis[j]
      if haskey(RG.basis_dict, g)
          RG.pm[i,j] = RG.basis_dict[g]
      else
         warning = true
      end
   end
   warning && @warn("Some products were not supported on basis")
   return RG
end

function initializepm!(RG::GroupRing; fill::Bool=false)
   hasbasis(RG) || throw("For baseless Group Rings You need to provide pm.")
   isdefined(RG, :pm) && return RG
   if fill
      RG.pm = try
         create_pm(RG.basis, RG.basis_dict)
      catch err
         isa(err, KeyError) && throw("Product is not supported on basis, $err.")
         throw(err)
      end
   else
      RG.pm = zeros(Int, length(RG.basis), length(RG.basis))
   end
   return RG
end

end # of module GroupRings


text:
using LowRankModels
using Base.Test

TOL = 1e-5

# QuadConstraint
r = QuadConstraint(7)
@test evaluate(r, ones(7)) == 0
@test evaluate(r, ones(100)) == Inf
@test_approx_eq prox(r, ones(100), 1) ones(100)/sqrt(100)*7

r = KSparseConstraint(3)
u = [-1,2,-3,4,-5.]
@test evaluate(r, u) == Inf
@test prox!(r, u, 1) == [0,0.,-3,4,-5.]


text:
"""
to_df(root::GeneralNode)::Tuple{Array{Float64}, Vector{String}}

This function returns a matrix representation of the tree structure and a vector with the column names.
The entry `mat[i,j]` is the length of the edge connecting node `i` with node `j`.
Returns Tuple containing the matrix and a vector of names.

* `root` : root of tree used to create matrix represenation.
"""
function to_df(root::GeneralNode)::Tuple{Array{Float64},Vector{String}}

    post_order_iteration = post_order(root)

    name_list = [i.name for i in post_order_iteration]
    num_list = [i.num for i in post_order_iteration]
    temp_ar = zeros(Float64, (length(post_order_iteration), length(post_order_iteration)))
    for i in post_order_iteration
        if i.nchild != 0
            ind = findfirst(isequal(i.num), num_list)
            for j in i.children
                ind2 = findfirst(isequal(j.num), num_list)
                temp_ar[ind[1], ind2[1]] = j.inc_length
            end # end for
        end # end if
    end # end for

    return temp_ar, name_list
end # end function to_df

#################### Covariance wrapper ####################

function to_covariance(tree::N) where {N<:AbstractNode}
    blv = get_branchlength_vector(tree)
    to_covariance(tree, blv)
end # end to_covariance

#################### To Matrix convertes ####################
"""
    to_distance_matrix(tree::T)::Array{Float64,2} where T <:AbstractNode

Calculate the distance matrix over the set of leaves.

Returns an Array of Floats.

* `tree` : root node of tree used to perform caclulcation.
"""
function to_distance_matrix(tree::T)::Array{Float64,2} where {T<:AbstractNode}
    leaves::Vector{T} = get_leaves(tree)
    ll = length(leaves)
    distance_mat = zeros(Float64, ll, ll)
    for i = 1:ll
        for j = 1:ll
            if i > j
                d = node_distance(tree, leaves[i], leaves[j])
                distance_mat[i, j] = d
                distance_mat[j, i] = d
            end # if
        end # for
    end #for
    distance_mat
end # function to_distance_matrix

"""
    to_covariance(tree::N, blv::Array{T})::Array{T,2} where {N<:AbstractNode,T<: Real}

Calcualte the variance-covariance matrix from `tree`. An entry (i,j) of the matrix
is defined as the length of the path connecting the latest common ancestor
of i and j with the root of the tree.

Returns an Array of Real numbers.

* `tree` : Node in tree of interest.

* `blv` : branchlength vector of tree.

"""
function to_covariance(tree::N, blv::Vector{T})::Array{T,2} where {N<:AbstractNode,T<:Real}
    leaves::Vector{N} = sort(get_leaves(tree), by = x->x.num)
    ll = length(leaves)
    covmat = zeros(T, ll, ll)

    @inbounds for ind = 1:ll, jnd = 1:ind
        itm = leaves[ind]
        if ind == jnd
            covmat[ind,jnd] =  reduce(+, @view blv[get_path(tree, itm)])
        else
            lca = find_lca(tree, itm, leaves[jnd])
            if !lca.root
                tmp = reduce(+, @view blv[get_path(tree, lca)])
                covmat[ind,jnd] = covmat[jnd,ind] = tmp
            end # if
        end # if

    end # for

    covmat
end# function to_covariance

"""
    function leave_incidence_matrix(root::G)::Matrix{Float64} where {G<:AbstractNode}

Calculate the incidence matrix of the tree whos root node is `root`
For a tree with ``m`` leaves and ``n`` vertecies this function returns an ``m \\times n`` matrix ``L``,
where ``L_{ij} = 1`` if vertex ``j`` is on the path from leave ``i`` to the root of the tree and ``0`` otherwise.

Returns leave incidence matrix.

* `root` : Root node of the tree
"""
function leave_incidence_matrix(root::G)::Matrix{Float64} where {G<:AbstractNode}
    n = length(get_branchlength_vector(root))
    leave_incidence_matrix(root, n)
end

function leave_incidence_matrix(root::G, n::Int)::Matrix{Float64} where {G<:AbstractNode}
    leaves::Vector{G} = get_leaves(root)
    out = zeros(length(leaves), n)
    for (i, leave) in enumerate(leaves)
        mother = leave
        while !mother.root
            out[i, mother.num] = 1.0
            mother = get_mother(mother)
        end
    end
    out
end


text:
#########
# Tests #
#########
__precompile__()
module tests

# include("catalan_test.jl")

# include("GeneralizedChartparsing\\src\\GeneralizedChartparsing.jl")
# using .GeneralizedChartparsing
# using .GeneralizedChartparsing: run_chartparser
# using LogProbs

include("FragmentGrammars.jl")
using .FragmentGrammars
using .FragmentGrammars: ApproxRule, update_approx_probs!, categorical_sample, LogProb, clone, boolean_dependency_dict, TreeNode, Tree, tree, logscore, initialize_fragments, init_mcmc, data, run_mcmc!

# TODO:
# -For Constituent conditioning (span-wise conditioning), modify parser_methods line 166 run_chartparser(input::Vector, grammar, dependency_matrix::AbstractMatrix{Bool}, parsing_method=:full; epsilon=missing) so that the dependency matrix is something like a Dict where we index by category and span.

# Test FG with base grammar:
# S -> S T | T S | T
# T -> a

println("-----start-----")
@time fg = FragmentGrammar(["S"], ["S"], [BaseRule("S", ("S", "T")), BaseRule("S", ("T", "S")), BaseRule("S", ("T",))], ["a"], [BaseRule("T", ("a",))], 0.5, 0.5)
# @time fg = FragmentGrammar([:S], [:S], [BaseRule(:S, (:S, :T)), BaseRule(:S, (:T, :S)), BaseRule(:S, (:T,))], [:a], [BaseRule(:T, (:a,))], 0.5, 0.5)
@show test_trees = Tree[tree("[S[S[S[S[T[a]]][T[a]]][T[a]]][T[a]]]") for i in 1:10]
# @show collect(fg.startstate)
# collect(fg.startstate)
@show logscore(fg)
@show logscore(fg).log
# @show initialize_fragments(test_tree, fg)
# @show init_mcmc(fg, Tree[test_tree])[test_tree][2] |> typeof
# update_approx_probs!(fg)

analyses = run_mcmc!(fg, test_trees, 100)
for i in keys(analyses)
    @show "Tree $i"
    for t in analyses[i][2]
        @show t.data[1]
    end
end
@show logscore(fg)
@show logscore(fg).log
# @show anal = Analysis(sample(fg, "S")...)
# add_obs!(fg, anal)
#
# # TO CHECK IF SUMS TO 1 -- WORKS FOR MCMC TOO
# for i in 1:1
#     marg = zero(LogProb)
#     num = zero(Int)
#     @show fg.startstate
#     for state in fg.startstate
#         for comp in state.comp
#             if comp[1] == "S"
#                 marg += @show comp[2].prob # sum(comp[2].probs)
#                 num += @show length(comp[2].probs)
#             end
#         end
#     end
#     @show marg
#     @show num
# end

# ---------- Parsing ------------
# for x in 1:1
#     for j in 1:1
#         for i in 1:1
#            # @time
#            anal = Analysis(sample(fg, "S")...)
#            # @time
#            add_obs!(fg, anal)
#            # rm_obs!(fg, anal)
#            # @time rm_obs!(fg, anal)
#         end
#         @show j
#     end
# end
end
    # @time update_approx_probs!(fg)
    # # @time update_approx_probs!(fg)
    # # @time forest = run_chartparser(["a" for i in 1:5], fg)
    # @time forest1 = run_chartparser(test_tree, fg)
    # @time sampled_approx_tree1 = sample_tree(forest1)
    # @time forest2 = run_chartparser(["a", "a", "a", "a"], fg)
    # @time sampled_approx_tree2 = sample_tree(forest2)
    # # @time sampled_approx_tree = sample_tree(forest)
    # # @show typeof(sampled_approx_tree)
    # @show 1
    # @time sampled_analysis_tree1 = sample(sampled_approx_tree2, fg)
    # @show sampled_analysis_tree1
    # for tree in sampled_approx_tree1
    #     println(tree.data[1], tree.data[2])
    # end
    # @show 2
    # @time sampled_analysis_tree2 = sample(sampled_approx_tree2, fg)
    # @show sampled_analysis_tree2
    # for tree in sampled_approx_tree2
    #     println(tree.data[1], tree.data[2])
    # end
    # anal = Analysis(sample(sampled_approx_tree, fg)...) # TODO: actually implement this.
    # @show anal.pointer

    # @time approx_sampled_rule = sample(sampled_approx_tree.data[2]) # sample from ApproxRule. Seems like this method is sometimes extremely slow for some reason.
    # @time approx_sampled_rule = sample(sampled_approx_tree.data[2])
    #
    # marg = zero(LogProb)
    # @show marg
    # num = zero(Int)
    # for state in fg.startstate
    #     for comp in state.comp
    #         if comp[1] === :S
    #             marg += @show sum(comp[2].probs)
    #             num += @show length(comp[2].probs)
    #         end
    #     end
    # end
    # @show marg
    # @show num
    # @show fg.CRP[:S].num_tables
    # @show fg.CRP[:S].num_customers
    # @show fg.DM
    # @show fg.BB

# @show fg
# @time for i in 1:1
#     anal = Analysis(sample(fg, "S")...)
#     add_obs!(fg, anal)
#     # println(fg.CRP)
#     println("----before----")
#     for state in fg.startstate
#         println("----new state----")
#         for comp in state.comp
#             println("$(comp[2].rules)")
#         end
#     end
#     rm_obs!(fg, anal)
#     # println(fg.CRP)
#     println("----after----")
#     for state in fg.startstate
#         println("----new state----")
#         for comp in state.comp
#             println("$(comp[2].rules)")
#         end
#     end
# end
# This is specifically for the parser + translating from approx PCFG to an analysis
# It should always be called before doing anything ApproxRule-related.
# for state in fg.startstate
#     for comp in state.comp
#         println(comp[2].rhs, length(comp[2].rules))
#     end
# end

# @time update_approx_probs!(fg)

# @time sampled_approx_tree = (run_chartparser(["a" for i in 1:10], fg) |> sample_tree)
# for tr in sampled_approx_tree
#     if tr.data[2] isa ApproxRule
#         for i in 1:10
#         @show sample(tr.data[2]) # sample from ApproxRule
#         end
#         break
#     end
# end

# --------- Observations --------
# for i in 1:5
#     println("-----New obs-----")
#     anal = Analysis(sample(fg, "S")...)
#     # println(anal)
#     add_obs!(fg, anal)
#     rm_obs!(fg, anal)
#     println(fg.startstate)
# end

# --------- Time tests ----------
# @time for i in 1:1 add_obs!(fg, Analysis(sample(fg, "S")...)) end

# println(Analysis(sample(fg,"S")...))

# @time for i in 1:100 add_obs!(fg, Analysis(sample(fg, "S")...)) end
#
# println(Analysis(sample(fg,"S")...))

# @time for i in 1:10000 add_obs!(fg, Analysis(sample(fg, "S")...)) end
#
# println(Analysis(sample(fg,"S")...))

# # for tr in fg.startstate.trans
# #     println(tr)
# # end
# #
# println("-----analysis-----")
# anal = Analysis(sample(fg, "S")...) # lol
# println(anal)
# println(anal)
#
# println("-----add-----")
# add_obs!(fg, anal)
#
# println("-----rm-----")
# rm_obs!(fg, anal)

# NOTE: The grammar below doesn't work. Always have unique terminal rules.
# Fragments must end with preterminals.
# g = Grammar(
# """0.5 S S a""",
# """0.5 S a""",
# ["S"], LogProb)

# NOTE This one works, for example.
# g = Grammar(
# """0.5 S S T
# 0.5 S T""",
# """1.0 T a""",
# ["S"], LogProb)
#
# @show run_chartparser(["a" for i in 1:4], g) |> best_tree |> leafs

# @show ContextFreeRule("S", ("S", "T")) == ContextFreeRule("S", ("S", "T"))
# @show ContextFreeRule("S", ("S", "T")) === ContextFreeRule("S", ("S", "T"))

# a = ContextFreeRule("S", ("S", "T"))
# @show a("test") |> typeof

# using Profile # Lets you see how many backtraces (function calls) come from a specific code block during runtime. Good for finding bottlenecks.


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Antic_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Antic")
JLLWrappers.@generate_main_file("Antic", UUID("e21ec000-9f72-519e-ba6d-10061e575a27"))
end  # module Antic_jll


text:
@testset "ResultConsumptionTestCase" begin
    tc = setup()

    @testset "test_can_consume_result_immediately" begin
        function w(tx, args, kwargs)
            result = run(tx, "UNWIND range(1, 3) AS n RETURN n")
            @test [value(record, 1) for record in result] == [1, 2, 3]
        end

        session(tc.driver) do sess
            read_transaction(sess, w)
        end
    end

    @testset "test_can_consume_result_from_buffer" begin
        function w(tx, args, kwargs)
            result = run(tx, "UNWIND range(1, 3) AS n RETURN n")
            Neo4jBolt.detach(result)
            @test [value(record, 1) for record in result] == [1, 2, 3]
        end

        session(tc.driver) do sess
            read_transaction(sess, w)
        end
    end

    @testset "test_can_consume_result_after_commit" begin
        session(tc.driver) do sess
            tx = begin_transaction(sess)
            result = run(tx, "UNWIND range(1, 3) AS n RETURN n")
            commit(tx)
            @test [value(record, 1) for record in result] == [1, 2, 3]
        end
    end

    @testset "test_can_consume_result_after_rollback" begin
        session(tc.driver) do sess
            tx = begin_transaction(sess)
            result = run(tx, "UNWIND range(1, 3) AS n RETURN n")
            rollback(tx)
            @test [value(record, 1) for record in result] == [1, 2, 3]
        end
    end

    @testset "test_can_consume_result_after_session_close" begin
        result = nothing
        session(tc.driver) do sess
            tx = begin_transaction(sess)
            result = run(tx, "UNWIND range(1, 3) AS n RETURN n")
            commit(tx)
        end
        @test [value(record, 1) for record in result] == [1, 2, 3]
    end

    @testset "test_can_consume_result_after_session_reuse" begin
        sess = session(tc.driver)
        tx = begin_transaction(sess)
        result_a = run(tx, "UNWIND range(1, 3) AS n RETURN n")
        commit(tx)
        close(sess)
        sess = session(tc.driver)
        tx = begin_transaction(sess)
        result_b = run(tx, "UNWIND range(4, 6) AS n RETURN n")
        commit(tx)
        close(sess)
        @test [value(record, 1) for record in result_a] == [1, 2, 3]
        @test [value(record, 1) for record in result_b] == [4, 5, 6]
    end

    @testset "test_can_consume_results_after_harsh_session_death" begin
        sess = session(tc.driver)
        result_a = run(sess, "UNWIND range(1, 3) AS n RETURN n")
        close(sess)
        sess = nothing
        sess = session(tc.driver)
        result_b = run(sess, "UNWIND range(4, 6) AS n RETURN n")
        close(sess)
        sess = nothing
        @test [value(record, 1) for record in result_a] == [1, 2, 3]
        @test [value(record, 1) for record in result_b] == [4, 5, 6]
    end

    @testset "test_can_consume_result_after_session_with_error" begin
        sess = session(tc.driver)
        @test_throws Neo4jBolt.JuliaBolt.CypherSyntaxError consume(run(sess, "X"))
        close(sess)

        sess = session(tc.driver)
        tx = begin_transaction(sess)
        result = run(tx, "UNWIND range(1, 3) AS n RETURN n")
        commit(tx)
        close(sess)

        @test [value(record, 1) for record in result] == [1, 2, 3]
    end

    @testset "test_single_with_exactly_one_record" begin
        sess = session(tc.driver)
        result = run(sess, "UNWIND range(1, 3) AS n RETURN n")
        record = single(result)
        values(record) == [1]
    end

    teardown(tc.driver)
end


text:
"""
Object we can scatter waves off

Subtypes will contain information about shape and material properties. Most
crucially, they will implement the [`t_matrix`](@ref) function
"""
abstract type AbstractParticle{T,Dim} end

"""
    Particle(medium::PhysicalMedium, shape::Shape)

Create particle with inner medium and shape (types and dimension must agree).
"""
struct Particle{T<:AbstractFloat,Dim,P<:PhysicalMedium,S<:Shape} <: AbstractParticle{T,Dim}
    medium::P
    shape::S
    # Enforce that the Dims and Types are all the same
    function Particle{T,Dim,P,S}(medium::P,shape::S) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim},S<:Shape{T,Dim}}
        new{T,Dim,P,S}(medium,shape)
    end
end

import Base.show
function show(io::IO, p::Particle)
    # Particle paramaters can be determined entirely from the medium and shape so we do not need to print them
    write(io, "Particle(")
    show(io, p.medium)
    write(io, ", ")
    show(io, p.shape)
    write(io, ")")
    return
end

"""
    CapsuleParticle(outer::Particle, inner::Particle)

A particle within another particle, both with the same shape type and origin.
"""
struct CapsuleParticle{T<:AbstractFloat,Dim,P<:PhysicalMedium,S<:Shape} <: AbstractParticle{T,Dim}
    outer::Particle{T,Dim,P,S}
    inner::Particle{T,Dim,P,S}
    # Enforce that particles are concentric
    function CapsuleParticle{T,Dim,P,S}(p2::Particle{T,Dim,P,S},p1::Particle{T,Dim,P,S}) where {T,Dim,P<:PhysicalMedium{T,Dim},S<:Shape}
        if origin(p1) != origin(p2) error("outer and inner particles should share the same origin") end
        if outer_radius(p1) >= outer_radius(p2)
            new{T,Dim,P,S}(p1,p2)
        else
            new{T,Dim,P,S}(p2,p1)
        end
    end
end

# Shorthand for all Vectors of particles
AbstractParticles{T<:AbstractFloat,Dim} = Vector{Pt} where Pt<:AbstractParticle{T,Dim}

# Convenience constructor which does not require explicit types/parameters
function Particle(medium::P,s::S) where {Dim,T,P<:PhysicalMedium{T,Dim},S<:Shape{T,Dim}}
    Particle{T,Dim,P,S}(medium,s)
end

"""
    Particle(medium, radius)

Returns a particle shaped like a sphere or circle, when the particle shape is not given and with the specified `radius`.
"""
function Particle(medium::P, radius::T) where {T, Dim, P <: PhysicalMedium{T,Dim}}
    Particle{T,Dim,P,Sphere{T,Dim}}(medium,Sphere(Dim,radius))
end

function CapsuleParticle(p1::Particle{T,Dim,P,S},p2::Particle{T,Dim,P,S}) where {T,Dim,S<:Shape,P<:PhysicalMedium}
    CapsuleParticle{T,Dim,P,S}(p1,p2)
end

shape(p::Particle) = p.shape
shape(p::CapsuleParticle) = p.outer.shape

# Shape holds infomation about origin of Particle
origin(p::AbstractParticle) = origin(shape(p))

boundary_points(p::AbstractParticle, num_points::Int = 3; kws...) = boundary_points(shape(p),num_points; kws...)

CircleParticle{T,P} = Particle{T,2,P,Sphere{T,2}}

outer_radius(p::AbstractParticle) = outer_radius(shape(p))
volume(p::AbstractParticle) = volume(shape(p))

bounding_box(p::AbstractParticle) = bounding_box(shape(p))
bounding_box(ps::AbstractParticles) = bounding_box([shape(p) for p in ps])

function volume(particles::AbstractParticles)
    mapreduce(volume, +, particles)
end

import Base.(==)
function ==(p1::Particle, p2::Particle)
    p1.medium == p2.medium &&
    p1.shape == p2.shape
end

function ==(p1::CapsuleParticle, p2::CapsuleParticle)
    p1.outer == p2.outer &&
    p1.inner == p2.inner
end

import Base.isequal
function isequal(p1::Particle, p2::Particle)
    isequal(p1.medium, p2.medium) &&
    isequal(p1.shape, p2.shape)
end

"""
    iscongruent(p1::AbstractParticle, p2::AbstractParticle)::Bool
    ≅(p1::AbstractParticle, p2::AbstractParticle)::Bool

Returns true if medium and shape of particles are the same, ignoring origin, false otherwise.
"""
iscongruent(p1::AbstractParticle, p2::AbstractParticle) = false # false by default, overload in specific examples

# Define synonym for iscongruent ≅, and add documentation
≅(p1::AbstractParticle, p2::AbstractParticle) = iscongruent(p1, p2)
@doc (@doc iscongruent(::AbstractParticle, ::AbstractParticle)) (≅(::AbstractParticle, ::AbstractParticle))

function iscongruent(p1::Particle, p2::Particle)
    p1.medium == p2.medium &&
    iscongruent(p1.shape, p2.shape)
end

function iscongruent(p1::CapsuleParticle, p2::CapsuleParticle)
    iscongruent(p1.inner, p2.inner) && iscongruent(p1.outer, p2.outer)
end

import Base.in
"""
    in(vector, particle)::Bool

Returns true if vector is in interior of shape of particle, false otherwise.
"""
in(x::AbstractVector, particle::AbstractParticle) = in(x, shape(particle))

import Base.issubset
issubset(s::Shape, particle::AbstractParticle) = issubset(s, shape(particle))
issubset(particle::AbstractParticle, s::Shape) = issubset(shape(particle), s)
issubset(particle1::AbstractParticle, particle2::AbstractParticle) = issubset(shape(particle1), shape(particle2))

"""
    estimate_regular_basis_order(medium::P, ω::Number, radius::Number; tol = 1e-6)
"""
function estimate_outgoing_basisorder(medium::PhysicalMedium, p::Particle, ω::Number; tol = 1e-6)

    k = ω / real(medium.c)

    # A large initial guess
    L = Int(round(4 * abs(k*outer_radius(p)))) + 1

    ts = [
        norm(diag(t_matrix(p, medium, ω, l)))
    for l = 1:L]
    ts = ts ./ ts[1];

    l = findfirst(diff(ts) .< tol)

    return l
end

# Add generic documentation from shape
@doc (@doc issubset(::Shape, ::Shape)) issubset(::Shape, ::AbstractParticle)
@doc (@doc issubset(::Shape, ::Shape)) issubset(::AbstractParticle, ::Shape)
@doc (@doc issubset(::Shape, ::Shape)) issubset(::AbstractParticle, ::AbstractParticle)


text:
"""
    coefplot(x::MixedModel;
             conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...)
    coefplot(x::MixedModelBootstrap;
             conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...)

Create a coefficient plot of the fixed-effects and associated confidence intervals.

!!! note
    This functionality is implemented using [Makie recipes](https://makie.juliaplots.org/v0.15.0/recipes.html)
    and thus there are also additional auto-generated methods for `coefplot` and `coefplot!` that may be useful
    when constructing more complex figures.
"""
function coefplot(
    x::Union{MixedModel,MixedModelBootstrap};
    conf_level=0.95,
    vline_at_zero=true,
    show_intercept=true,
    attributes...,
)
    # need to guarantee a min height of 150
    fig = Figure(; resolution=(640, max(150, 75 * _npreds(x; show_intercept))))
    ax = Axis(fig[1, 1])
    pl = coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, attributes...)
    return Makie.FigureAxisPlot(fig, ax, pl)
end

@recipe(CoefPlot, x) do scene
    return Attributes(; conf_level=0.95, vline_at_zero=true, show_intercept=true)
end

function Makie.plot!(ax::Axis, P::Type{<:CoefPlot}, allattrs::Makie.Attributes, x)
    plot = Makie.plot!(ax.scene, P, allattrs, x)

    if haskey(allattrs, :title)
        ax.title = allattrs.title[]
    end
    if haskey(allattrs, :xlabel)
        ax.xlabel = allattrs.xlabel[]
    else
        ax.xlabel = @sprintf "Estimate and %g%% confidence interval" (
            allattrs.conf_level[] * 100
        )
    end
    if haskey(allattrs, :ylabel)
        ax.ylabel = allattrs.ylabel[]
    end
    reset_limits!(ax)
    show_intercept = allattrs.show_intercept[]
    cn = _coefnames(x; show_intercept)
    nticks = _npreds(x; show_intercept)
    ax.yticks = (nticks:-1:1, cn)
    ylims!(ax, 0, nticks + 1)
    allattrs.vline_at_zero[] && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash)
    return plot
end

function Makie.plot!(plot::CoefPlot{<:Tuple{Union{MixedModel,MixedModelBootstrap}}})
    model_or_boot, conf_level = plot[1][], plot.conf_level[]
    ci = confint_table(model_or_boot, conf_level)
    plot.show_intercept[] || filter!(:coefname => !=("(Intercept)"), ci)
    y = nrow(ci):-1:1
    xvals = ci.estimate
    scatter!(plot, xvals, y)
    errorbars!(plot, xvals, y, xvals .- ci.lower, ci.upper .- xvals; direction=:x)

    return plot
end


text:
function check_model_form(m)

  if typeof(m) <:Blanchard_Kahn_Form
    assessment = check_blanchard_kahn_form(m)
  elseif typeof(m) <:Klein_Form
    assessment = check_klein_form(m)
  elseif typeof(m) <:Structural_Form
    assessment = check_structural_form(m)
  elseif typeof(m) <:Sims_Form
    assessment = check_sims_form(m)
  elseif typeof(m) <:Gomme_Klein_Form
    assessment = check_gomme_klein_form(m)
  elseif typeof(m) <:Lombardo_Sutherland_Form
    assessment = check_lombardo_sutherland_form(m)
  end

  return assessment

end

function check_blanchard_kahn_form(m::Blanchard_Kahn_Form)

  assessment = "Pass"

  if m.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif m.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif m.nx+m.ny != size(m.a,1)
    assessment = "Fail: A is not conformable with nx and ny"
  elseif size(m.a,1) != size(m.a,2)
    assessment = "Fail: A is not square"
  elseif size(m.a,1) != size(m.c,1)
    assessment = "Fail: A and C must have the same number of rows"
  elseif size(m.c,2) > m.nx
    assessment = "Fail: There are more innovations than predetermined variables"
    elseif size(m.sigma,1) != size(m.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(m.c,2) != size(m.sigma,1)
    assessment = "Fail: C and SIGMA are not conformable"
  end

  return assessment

end

function check_klein_form(m::Klein_Form)

  assessment = "Pass"

  if m.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif m.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif m.nx+m.ny != size(m.a,1)
    assessment = "Fail: A is not conformable with nx and ny"
  elseif size(m.a,1) != size(m.a,2)
    assessment = "Fail: A is not square"
  elseif size(m.b,1) != size(m.b,2)
    assessment = "Fail: B is not square"
  elseif size(m.a) != size(m.b)
    assessment = "Fail: A and B are not conformable"
  elseif size(m.a,1) != size(m.c,1)
    assessment = "Fail: A and C must have the same number of rows"
  elseif size(m.c,2) > m.nx
    assessment = "Fail: There are more innovations than predetermined variables"
    elseif size(m.sigma,1) != size(m.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(m.c,2) != size(m.sigma,1)
    assessment = "Fail: C and SIGMA are not conformable"
  end

  return assessment

end

function check_structural_form(m::Structural_Form)

  assessment = "Pass"

  if size(m.a,1) != size(m.a,2)
    assessment = "Fail: A is not square"
  elseif size(m.a) != size(m.a1)
    assessment = "Fail: A and A1 are not conformable"
  elseif size(m.a) != size(m.b)
    assessment = "Fail: A and B are not conformable"
  elseif size(m.a,1) != size(m.c,1)
    assessment = "Fail: A and C must have the same number of rows"
  elseif size(m.sigma,1) != size(m.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(m.c,2) != size(m.sigma,1)
    assessment = "Fail: C and SIGMA are not conformable"
  end

  return assessment

end

function check_sims_form(m::Sims_Form)

  assessment = "Pass"

  if size(m.gamma0,1) != size(m.gamma0,2)
    assessment = "Fail: GAMMA0 is not square"
  elseif size(m.gamma1,1) != size(m.gamma1,2)
    assessment = "Fail: GAMMA1 is not square"
  elseif size(m.gamma0) != size(m.gamma1)
    assessment = "Fail: GAMMA0 and GAMMA1 are not conformable"
  elseif size(m.gamma0,1) != size(m.c,1)
    assessment = "Fail: GAMMA0 and C must have the same number of rows"
  elseif size(m.c,2) != 1
    assessment = "Fail: C must be a vector"
  elseif size(m.gamma0,1) != size(m.psi,1)
    assessment = "Fail: GAMMA0 and PSI must have the same number of rows"
  elseif size(m.gamma0,1) > size(m.pi,1)
    assessment = "Fail: GAMMA0 and PI must have the same number of rows"
  elseif size(m.sigma,1) != size(m.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(m.psi,2) != size(m.sigma,1)
    assessment = "Fail: PSI and SIGMA are not conformable"
  end

  return assessment

end

function check_gomme_klein_form(m::Gomme_Klein_Form)

  assessment = "Pass"

  if m.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif m.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif m.nx+m.ny != size(m.derivs,1)
    assessment = "Fail: Need first derivative for all model equations"
  elseif 2*size(m.derivs,1) != size(m.derivs,2)
    assessment = "Fail: The dimensions of the first-derivative matrix are incorrect"
  elseif size(m.derivs,2) != size(m.hessians,2)
    assessment = "Fail: The derivative matrix and the hessian matrix are not conformable"
  elseif size(m.hessians,1) != (m.nx+m.ny)*size(m.hessians,2)
    assessment = "Fail: The dimensions of the hessian matrix are incorrect"
  elseif size(m.eta,1) != m.nx
    assessment = "Fail: Eta has too few rows"
  elseif size(m.eta,2) != size(m.sigma,1)
    assessment = "Fail: Eta has too few columns"
  elseif size(m.sigma,1) != size(m.sigma,2)
    assessment = "Fail: SIGMA must be square"
  end

  return assessment

end

function check_lombardo_sutherland_form(m::Lombardo_Sutherland_Form)

  assessment = "Pass"

  if m.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif m.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif m.nx+m.ny != size(m.derivs,1)
    assessment = "Fail: Need first derivative for all model equations"
  elseif 2*size(m.derivs,1) != size(m.derivs,2)
    assessment = "Fail: The dimensions of the first-derivative matrix are incorrect"
  elseif size(m.derivs,2) != size(m.hessians,2)
    assessment = "Fail: The derivative matrix and the hessian matrix are not conformable"
  elseif size(m.hessians,1) != (m.nx+m.ny)*size(m.hessians,2)
    assessment = "Fail: The dimensions of the hessian matrix are incorrect"
  elseif size(m.eta,1) != m.nx
    assessment = "Fail: Eta has too few rows"
  elseif size(m.eta,2) != size(m.sigma,1)
    assessment = "Fail: Eta has too few columns"
  elseif size(m.sigma,1) != size(m.sigma,2)
    assessment = "Fail: SIGMA must be square"
  end

  return assessment

end

max_stars_repo_name: dingliumath/SolveDSGE

max_issues_repo_name: dingliumath/SolveDSGE

max_forks_repo_name: dingliumath/SolveDSGE


text:
"""
does the same thing as `x[..., :-1]` in pytorch.
"""
function head(x::T) where {T <: AbstractArray}
  return x[[axes(x, i) for i in 1:ndims(x)-1]..., 1:end-1]
end

"""
returns: x[[axes(x, i) for i in 1:ndims(x)-1]..., l]
"""
function slice(x::T, l::Integer) where {T <: AbstractArray}
  return x[[axes(x, i) for i in 1:ndims(x)-1]..., l]
end

"""
returns: x[[axes(x, i) for i in 1:ndims(x)-1]..., end]
"""
function tail(x::T) where {T <: AbstractArray}
  return x[[axes(x, i) for i in 1:ndims(x)-1]..., end]
end

"""
logsumexp
"""
function logsumexp(logx::T; dims::Integer=1) where {T <: AbstractArray}
  d = dims < 0 ? ndims(logx) + 1 + dims : dims
  mx = maximum(logx, dims=d)
  return mx .+ log.(sum(exp.(logx .- mx), dims=d))
end

"""
return: dropdims(logsumexp(logx, dims=dims), dims=dims)
"""
function logsumexpdd(logx::T; dims::Integer=1) where {T <: AbstractArray}
  d = dims < 0 ? ndims(logx) + 1 + dims : dims
  return dropdims(logsumexp(logx, dims=d), dims=d)
end

"""
return: dropdims(sum(x, dims=dims), dims=dims)
"""
function sumdd(x::T; dims::Integer=1) where {T <: AbstractArray}
  d = dims < 0 ? ndims(x) + 1 + dims : dims
  return dropdims(sum(x, dims=d), dims=d)
end

"""
return: dropdims(prod(x, dims=dims), dims=dims)
"""
function proddd(x::T; dims::Integer=1) where {T <: AbstractArray}
  d = dims < 0 ? ndims(x) + 1 + dims : dims
  return dropdims(prod(x, dims=d), dims=d)
end

function lpdf_normal(x::X, m::M, s::S) where {X <: Real, M <: Real, S<:Real}
  z = (x - m) / s
  # return oftype(z, -0.5 * log(2*pi) - z^2 * 0.5 - log(s))
  return -0.5 * log(2*pi) - z^2 * 0.5 - log(s)
end

"""
log pdf of gaussian mixture model
"""
function lpdf_gmm(x::TX, m::TM, s::TS, w::TW; dims::Integer, dropdim::Bool=false) where {TX, TW, TM, TS}
  if dropdim
    return logsumexpdd(log.(w) .+ lpdf_normal.(x, m, s), dims=dims)
  else
    return logsumexp(log.(w) .+ lpdf_normal.(x, m, s), dims=dims)
  end
end

"""
cat(a, b, dims=ndims(A) + 1)
"""
function stack(a::A, b::B) where {A, B}
  @assert size(a) == size(b)
  return cat(a, b, dims=ndims(a) + 1)
end

function lpdf_dirichlet(p::AbstractArray, a::AbstractArray)
  if size(a) != size(p)
    a = reshape(a, (one.(size(p)[1:end-1])..., length(a)))
  end
  K = ndims(p)
  return Tracker.lgamma.(sumdd(a, dims=K)) .- sumdd(Tracker.lgamma.(a), dims=K) .+ sumdd((a .- 1) .* log.(p), dims=K)
end

function lpdf_beta(p::P, a::A, b::B) where {P <: Real, A <: Real, B <: Real}
  if 0 < p < 1
    return Tracker.lgamma(a + b) - Tracker.lgamma(a) - Tracker.lgamma(b) + (a-1)*log(p) + (b-1)*log1p(-p)
  else
    return -Inf
  end
end

function lpdf_gamma(x::X, shape::A, scale::B) where {X <: Real, A <: Real, B <: Real}
  if x <= 0
    return -Inf
  else
    return -(Tracker.lgamma(shape) + shape*log(scale)) + (shape - 1)*log(x) - x / scale
  end
end

function lpdf_lognormal(x::X, m::A, s::B) where {X <: Real, A <: Real, B <: Real}
  if x <= 0
    return -Inf
  else
    lx = log(x)
    return lpdf_normal(lx, m, s) - lx
  end
end

function lpdf_uniform(x::X, a::A, b::B) where {X <: Real, A <: Real, B <: Real}
  @assert b > a
  if a <= x <= b
    return -log(b - a)
  else
    return -Inf
  end
end

compute_lpdf(dist::Dirichlet, x::X) where X = lpdf_dirichlet(x, dist.alpha)
compute_lpdf(dist::Normal, x::X) where X = lpdf_normal.(x, dist.μ, dist.σ)
compute_lpdf(dist::Beta, x::X) where X = lpdf_beta.(x, dist.α, dist.β)
compute_lpdf(dist::Gamma, x::X) where X = lpdf_gamma.(x, dist.α, dist.θ)
compute_lpdf(dist::LogNormal, x::X) where X = lpdf_lognormal.(x, dist.μ, dist.σ)
compute_lpdf(dist::Uniform, x::X) where X = lpdf_uniform.(x, dist.a, dist.b)

function sigmoid_safe(x::T; eps::Float64=1e-8) where {T <: Real}
  return T(clamp(sigmoid(x), eps, 1 - eps))
end

function logit_safe(p::T; eps::Float64=1e-8) where {T <: Real}
  p = T(clamp(p, eps , 1 - eps))
  return log(p) - log1p(-p)
end


text:
# clique state machine for tree based initialization and inference

# newer exports
# export towardUpOrDwnSolve_StateMachine, maybeNeedDwnMsg_StateMachine
# export prepInitUp_StateMachine, doCliqUpSolveInitialized_StateMachine
# export rmUpLikeliSaveSubFg_StateMachine
# export blockCliqSiblingsParentChildrenNeedDown_StateMachine

## ============================================================================================
## Only fetch related functions after this, all to be deprecated
## ============================================================================================

# XXX only fetch
"""
    $SIGNATURES

Function to iterate through while initializing various child cliques that start off `needdownmsg`.

Notes
- State machine function nr. 7e
"""
function slowWhileInit_StateMachine(csmc::CliqStateMachineContainer)

  if doAnyChildrenNeedDwnMsg(csmc.tree, csmc.cliq)
    infocsm(csmc, "7e, slowWhileInit_StateMachine, must wait for new child messages.")

    # wait for THIS clique to be notified (PUSH NOTIFICATION FROM CHILDREN at `prepPutCliqueStatusMsgUp!`)
    wait(getSolveCondition(csmc.cliq))
  end

  # go to 8f
  return prepInitUp_StateMachine
end

# XXX only fetch
"""
    $SIGNATURES

Delay loop if waiting on upsolves to complete.

Notes
- State machine 7b
- Differs from 4e in that here children must be "upsolved" or equivalent to continue.
  - Also waits on condition, so less succeptible to stale messages from children

DevNotes
- CONSIDER "recursion", to come back to this function to make sure that child clique updates are in fact upsolved.
  - other steps in CSM require each CSM to progress until cascading init problem is solved
"""
function slowIfChildrenNotUpSolved_StateMachine(csmc::CliqStateMachineContainer)

  # childs = getChildren(csmc.tree, csmc.cliq)
  # len = length(childs)
  for chld in getChildren(csmc.tree, csmc.cliq)
    chst = getCliqueStatus(chld)
    if !(chst in [:upsolved;:uprecycled;:marginalized;])
      infocsm(csmc, "7b, slowIfChildrenNotUpSolved_StateMachine, wait $(chst), cliq=$(chld.index), ch_lbl=$(getCliqFrontalVarIds(chld)[1]).")

      # wait for child clique status/msg to be updated
      wait(getSolveCondition(chld))
      # tsk = @async
      # timedwait(()->tsk.state==:done, 10)

      # check again and reroute if :needdownmsg
      # chst = getCliqueStatus(chld)
      # if chst == :needdownmsg
      #   return prntPrepDwnInitMsg_StateMachine
      # end
    end
  end

  # go to 4b
  return trafficRedirectConsolidate459_StateMachine
end

# XXX only fetch
"""
    $SIGNATURES

Block until all children have a csm status, using `fetch---Condition` model.

Notes
- State machine function nr.4e
- Blocking call if there are no status messages available
  - Will not block on stale message
"""
function blockUntilChildrenHaveStatus_StateMachine(csmc::CliqStateMachineContainer)
  #must happen before if :null

  notsolved = true

  while notsolved
    notsolved = false
    infocsm(csmc, "4e, blockUntilChildrenHaveStatus_StateMachine, get new status msgs.")

    stdict = fetchChildrenStatusUp(csmc.tree, csmc.cliq, csmc.logger)
    for (cid, st) in stdict
      infocsm(csmc, "4e, blockUntilChildrenHaveStatus_StateMachine, maybe wait cliq=$(cid), child status=$(st).")

      if st in [:null;]
        infocsm(csmc, "4e, blockUntilChildrenHaveStatus_StateMachine, waiting cliq=$(cid), child status=$(st).")
        wait(getSolveCondition(csmc.tree.cliques[cid]))
        notsolved = true
        break
      end
    end
  end

  # go to 4b
  return trafficRedirectConsolidate459_StateMachine
end

# XXX only fetch
"""
    $SIGNATURES

Decide whether to pursue and upward or downward solve with present state.

Notes
- State machine function nr. 7c
"""
function towardUpOrDwnSolve_StateMachine(csmc::CliqStateMachineContainer)

  sleep(0.1) # FIXME remove after #459 resolved
  # return doCliqInferAttempt_StateMachine
  cliqst = getCliqueStatus(csmc.cliq)
  infocsm(csmc, "7c, status=$(cliqst), before picking direction")
  # d1,d2,cliqst = doCliqInitUpOrDown!(csmc.cliqSubFg, csmc.tree, csmc.cliq, isprntnddw)
  if cliqst == :needdownmsg && !isCliqParentNeedDownMsg(csmc.tree, csmc.cliq, csmc.logger)
      # FIXME, end 459, was collectDwnInitMsgFromParent_StateMachine but 674 requires pull model
      # go to 8c
      # notifyCSMCondition(csmc.cliq)
      return waitChangeOnParentCondition_StateMachine
    # # go to 8a
    # return collectDwnInitMsgFromParent_StateMachine
  # HALF DUPLICATED IN STEP 4
  elseif cliqst == :marginalized
    # go to 1
    return isCliqUpSolved_StateMachine
  end

  # go to 8b
  return attemptCliqInitUp_StateMachine
end

# """
#     $SIGNATURES

# Blocking case when all siblings and parent :needdownmsg.

# Notes
# - State machine function nr. 6c

# DevNotes
# - FIXME understand if this should be consolidated with 4b. `trafficRedirectConsolidate459_StateMachine`?
# - FIXME understand if this should be consolidated with 7c. `towardUpOrDwnSolve_StateMachine`
# """
# function doesParentNeedDwn_StateMachine(csmc::CliqStateMachineContainer)

#   infocsm(csmc, "6c, check/block sibl&prnt :needdownmsg")
#   prnt = getParent(csmc.tree, csmc.cliq)
#   if 0 == length(prnt) || getCliqueStatus(prnt[1]) != :needdownmsg
#     infocsm(csmc, "6c, prnt $(getCliqueStatus(prnt[1]))")
#     # go to 7
#     return determineCliqNeedDownMsg_StateMachine
#   # elseif
#   end

#   # go to 6d
#   return doAllSiblingsNeedDwn_StateMachine
# end

# XXX only fetch
"""
    $SIGNATURES

Redirect CSM traffic in various directions

Notes
- State machine function nr.4b
- Was refactored during #459 dwnMsg effort.

DevNotes
- Consolidate with 7?
"""
function trafficRedirectConsolidate459_StateMachine(csmc::CliqStateMachineContainer)

  cliqst = getCliqueStatus(csmc.cliq)
  infocsm(csmc, "4b, trafficRedirectConsolidate459_StateMachine, cliqst=$cliqst")

  # if no parent or parent will not update

  # for recycle computed clique values case
  if csmc.incremental && cliqst == :downsolved
    csmc.incremental = false
    # might be able to recycle the previous clique solve, go to 0b
    return checkChildrenAllUpRecycled_StateMachine
  end

  # Some traffic direction
  if cliqst == :null
    # go to 4d
    return maybeNeedDwnMsg_StateMachine
  end

  # go to 6c
  # return doesParentNeedDwn_StateMachine
  prnt = getParent(csmc.tree, csmc.cliq)
  if 0 < length(prnt) && cliqst == :needdownmsg
    if getCliqueStatus(prnt[1]) == :needdownmsg
      # go to 8c
      return waitChangeOnParentCondition_StateMachine
    end
    # go to 6d
    return doAllSiblingsNeedDwn_StateMachine
  end

  # go to 7
  return determineCliqNeedDownMsg_StateMachine
end

## ============================================================================================
# START of dwnmsg consolidation bonanza
## ============================================================================================
#
# blockCliqSiblingsParentChildrenNeedDown_  # Blocking case when all siblings and parent :needdownmsg.
# doAllSiblingsNeedDwn_StateMachine   # Trying to figure out when to block on siblings for cascade down init.
# maybeNeedDwnMsg_StateMachine         # If all children (then also escalate to) :needdownmsgs and block until sibling status.
# determineCliqNeedDownMsg_StateMachine     # Try decide whether this `csmc.cliq` needs a downward initialization message.
# doAnyChildrenNeedDwn_StateMachine         # Determine if any one of the children :needdownmsg.
# downInitRequirement_StateMachine          # Place fake up msg and notify down init status if any children :needdownmsg

# XXX only fetch
"""
$SIGNATURES

If all children (then also escalate to) :needdownmsgs and block until sibling status.

Notes
- State machine function nr.4d

DevNotes
- TODO consolidate with 6d?????
- TODO Any overlap with nr.4c??
"""
function maybeNeedDwnMsg_StateMachine(csmc::CliqStateMachineContainer)

  # fetch (should not block)
  stdict = fetchChildrenStatusUp(csmc.tree, csmc.cliq, csmc.logger)
  chstatus = collect(values(stdict))
  infocsm(csmc,"fetched all, keys=$(keys(stdict)), values=$(chstatus).")
  len = length(chstatus)

  # if all children needdownmsg
  if 0 < len && sum(chstatus .== :needdownmsg) == len
    # Can this cliq init with local information?
    # get initial estimate of solvable dims in this cliq
    sdims = getCliqVariableMoreInitDims(csmc.cliqSubFg, csmc.cliq)
    infocsm(csmc, "4d, maybeNeedDwnMsg_StateMachine, sdims=$(sdims)")
    updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)
    if 0 < sum(collect(values(sdims)))
      # try initialize
      # go to 8f
      return prepInitUp_StateMachine
    end

    # TODO maybe can happen where some children need more information?
    infocsm(csmc, "4d, maybeNeedDwnMsg_StateMachine, escalating to :needdownmsg since all children :needdownmsg")
    # NOTE, trying consolidation with prepPutUp for #459 effort
    setCliqueDrawColor!(csmc.cliq, "orchid1")
    prepPutCliqueStatusMsgUp!(csmc, :needdownmsg)

    # debuggin #459 transition
    infocsm(csmc, "4d, maybeNeedDwnMsg_StateMachine -- finishing before going to  blockSiblingStatus_StateMachine")

    # go to 5
    return blockSiblingStatus_StateMachine
  end

  if doAnyChildrenNeedDwnMsg(csmc.tree, csmc.cliq)
    # go to 7e (#754)
    return slowWhileInit_StateMachine
  end

  # go to 7
  return determineCliqNeedDownMsg_StateMachine
end

# XXX only fetch
"""
    $SIGNATURES

Try decide whether this `csmc.cliq` needs a downward initialization message.

Notes
- State machine function nr. 7

DevNotes
- Consolidate with 4b?
- found by trail and error, TODO review and consolidate with rest of CSM after major #459 and PCSM consolidation work is done.
"""
function determineCliqNeedDownMsg_StateMachine(csmc::CliqStateMachineContainer)

  # fetch children status
  stdict = fetchChildrenStatusUp(csmc.tree, csmc.cliq, csmc.logger)
  infocsm(csmc,"fetched all, keys=$(keys(stdict)).")

  # hard assumption here on upsolve from leaves to root
  childst = collect(values(stdict))
  # are child cliques sufficiently solved
  resolveinit = (filter(x-> x in [:upsolved;:marginalized;:downsolved;:uprecycled], childst) |> length) == length(childst)
  chldupandinit = sum(childst .|> x-> (x in [:initialized;:upsolved;:marginalized;:downsolved;:uprecycled])) == length(childst)
  allneeddwn = (filter(x-> x == :needdownmsg, childst) |> length) == length(childst) && 0 < length(childst)
  chldneeddwn = :needdownmsg in childst
  chldnull = :null in childst

  cliqst = getCliqueStatus(csmc.cliq)

  infocsm(csmc, "7, determineCliqNeedDownMsg_StateMachine, childst=$childst")
  infocsm(csmc, "7, determineCliqNeedDownMsg_StateMachine, cliqst=$cliqst, resolveinit=$resolveinit, allneeddwn=$allneeddwn, chldneeddwn=$chldneeddwn, chldupandinit=$chldupandinit, chldnull=$chldnull")

  # merged in from 4c here into 7, part of dwnMsg #459
  if cliqst == :needdownmsg
    if allneeddwn || length(stdict) == 0
      infocsm(csmc, "7, determineCliqNeedDownMsg_StateMachine, at least some children :needdownmsg")
      # # go to 8j
      return dwnInitSiblingWaitOrder_StateMachine
    end
    # FIXME FIXME FIXME hex 2.16 & 3.21 should go for downinitfrom parent
    # perhaps should check if all siblings also needdownmsg???

    if chldneeddwn
      # go to 5
      return blockSiblingStatus_StateMachine
      # # go to 7b
      # return slowIfChildrenNotUpSolved_StateMachine
    end
  end

  # includes case of no children
  if resolveinit && !chldneeddwn
    # go to 8j (dwnMsg #459 WIP 9)
    # return dwnInitSiblingWaitOrder_StateMachine
    # go to 7c
    return towardUpOrDwnSolve_StateMachine
  end

  if chldnull || cliqst == :null && !chldupandinit
    # go to 4e
    return blockUntilChildrenHaveStatus_StateMachine
  elseif chldneeddwn || chldupandinit
    # go to 7b
    return slowIfChildrenNotUpSolved_StateMachine
  end

  # go to 6d
  return doAllSiblingsNeedDwn_StateMachine
end

# XXX only fetch
"""
    $SIGNATURES

Detect block on clique's parent.

Notes
- State machine function nr. 5

DevNotes
- FIXME refactor this for when prnt==:null, cliq==:needdownmsg/init
"""
function blockSiblingStatus_StateMachine(csmc::CliqStateMachineContainer)
    # infocsm(csmc, "5, blocking on parent until all sibling cliques have valid status")
    # setCliqueDrawColor!(csmc.cliq, "blueviolet")

    cliqst = getCliqueStatus(csmc.cliq)
    infocsm(csmc, "5, block on siblings")
    prnt = getParent(csmc.tree, csmc.cliq)
    if cliqst == :needdownmsg && 0 < length(prnt) && getCliqueStatus(prnt[1]) in [:null;:needdownmsg]
      # go to 8c
      return waitChangeOnParentCondition_StateMachine
    end

    # infocsm(csmc, "5, has parent clique=$(prnt[1].index)")
    # ret = fetchChildrenStatusUp(csmc.tree, prnt[1], csmc.logger)
    # infocsm(csmc,"prnt $(prnt[1].index), fetched all, keys=$(keys(ret)).")

  infocsm(csmc, "5, finishing")
  # go to 6c
  # return doesParentNeedDwn_StateMachine

  # go to 6d
  return doAllSiblingsNeedDwn_StateMachine
end

# XXX only fetch
"""
$SIGNATURES

Trying to figure out when to block on siblings for cascade down init.

Notes
- State machine function nr.6d
- Assume there must be a parent.
- Part of #459 dwnMsg consolidation work.
- used for regulating long need down message chains.
- exit strategy is parent becomes status `:initialized`.

DevNotes
- Consolidation with work with similar likely required.
"""
function doAllSiblingsNeedDwn_StateMachine(csmc::CliqStateMachineContainer)

  # go to 6c
  # return doesParentNeedDwn_StateMachine
  infocsm(csmc, "7, check/block sibl&prnt :needdownmsg")
  prnt = getParent(csmc.tree, csmc.cliq)
  if 0 == length(prnt) # || getCliqueStatus(prnt[1]) != :needdownmsg
    # infocsm(csmc, "4d, prnt $(getCliqueStatus(prnt[1]))")
    # go to 7
    return determineCliqNeedDownMsg_StateMachine
  # elseif
  end

  prnt = getParent(csmc.tree, csmc.cliq)
  prnt_ = prnt[1]

  stdict = fetchChildrenStatusUp(csmc.tree, prnt_, csmc.logger)
  # remove this clique's data from dict
  delete!(stdict, csmc.cliq.index)
  # hard assumption here on upsolve from leaves to root
  siblst = values(collect(stdict))
  # are child cliques sufficiently solved
  allneeddwn = (filter(x-> x == :needdownmsg, siblst) |> length) == length(siblst) && 0 < length(siblst)

  # FIXME, understand why is there another status event from parent msg here... How to consolidate this with CSM 8a
  if allneeddwn
    # # go to 6e
    # return slowOnPrntAsChildrNeedDwn_StateMachine
    # go to 8j
    return dwnInitSiblingWaitOrder_StateMachine
  end

  # go to 7
  return determineCliqNeedDownMsg_StateMachine
end

# XXX only fetch
"""
$SIGNATURES

Place up msg and notify down init status if any children :needdownmsg.

Notes
- StateMachine function nr. 8d
- TODO figure out if this function is duplication of other needdwnmsg functionality?
  - Consolidate with 7???
  - Consolidate with 4c???
"""
function decideUpMsgOrInit_StateMachine(csmc::CliqStateMachineContainer)
  #
  infocsm(csmc, "8d, downInitRequirement_StateMachine., start")

  # children = getChildren(csmc.tree, csmc.cliq)
  # if doAnyChildrenNeedDwnMsg(children)
  someChildrenNeedDwn = false
  # for ch in getChildren(csmc.tree, csmc.cliq)
  for (clid, chst) in fetchChildrenStatusUp(csmc.tree, csmc.cliq, csmc.logger)
    # if getCliqueStatus(ch) == :needdownmsg
    if chst == :needdownmsg # NOTE was :needdowninit
      someChildrenNeedDwn = true
      break
    end
  end

  if someChildrenNeedDwn
    # send a down init message
    # TODO down message is sent more than once?
    # here and in sendCurrentUpMsg_StateMachine
    prepPutCliqueStatusMsgDwn!(csmc)
    # go to 8k
    return sendCurrentUpMsg_StateMachine
  end

  # go to 8b
  return attemptCliqInitUp_StateMachine
end

# XXX only fetch
"""
    $SIGNATURES

Determine if up initialization calculations should be attempted.

Notes
- State machine function nr. 8b
"""
function attemptCliqInitUp_StateMachine(csmc::CliqStateMachineContainer)

  # should calculations be avoided.
  notChildNeedDwn = !doAnyChildrenNeedDwnMsg(csmc.tree, csmc.cliq)
  infocsm(csmc, "8b, attemptCliqInitUp, !doAnyChildrenNeedDwnMsg()=$(notChildNeedDwn)" )
  if getCliqueStatus(csmc.cliq) in [:initialized; :null; :needdownmsg] && notChildNeedDwn
    # go to 8f.
    return prepInitUp_StateMachine
  end

  if !notChildNeedDwn
    # go to 8m
    return tryUpInitCliq_StateMachine
  end

  # go to 9
  return checkUpsolveFinished_StateMachine
end

## ============================================================================================
# End of dwnmsg consolidation bonanza
## ============================================================================================

#


text:
using AWSTools
using FilePathsBase
using UUIDs

using AWSTools.CloudFormation: stack_output
using AWSTools.S3: sync, upload
using AWSS3: AWSS3, S3Path, s3_create_bucket, s3_put, s3_delete_bucket

# Enables the running of the "batch" online tests. e.g ONLINE=batch
const ONLINE = strip.(split(get(ENV, "ONLINE", ""), r"\s*,\s*"))

# Run the online S3 tests on the bucket specified
const TEST_BUCKET_AND_PREFIX = let
    p = if haskey(ENV, "TEST_BUCKET_AND_PREFIX")
        ENV["TEST_BUCKET_AND_PREFIX"]
    elseif haskey(ENV, "AWS_STACKNAME")
        output = stack_output(ENV["AWS_STACKNAME"])

        if haskey(output, "TestBucketAndPrefix")
            "s3://" * output["TestBucketAndPrefix"] * "/AWSTools.jl"
        else
            nothing
        end
    else
        nothing
    end

    p !== nothing ? rstrip(p, '/') : p
end

function compare(src_file::AbstractPath, dest_file::AbstractPath)
    if isdir(src_file)
        @test isdir(dest_file)
        @test basename(src_file) == basename(dest_file)
    else
        @test isfile(src_file)
        @test isfile(dest_file)
        @test basename(dest_file) == basename(src_file)
        @test size(dest_file) == size(src_file)
        @test modified(dest_file) >= modified(src_file)

        # Test file contents are equal
        @test read(dest_file) == read(src_file)
    end
end

function compare_dir(src_dir::AbstractPath, dest_dir::AbstractPath)
    @test isdir(dest_dir)

    src_contents = collect(walkpath(src_dir))
    dest_contents = collect(walkpath(dest_dir))

    for (src, dst) in zip(src_contents, dest_contents)
        compare(src, dst)
    end
end

function bucket_and_key(s3_path::AbstractString)
    m = match(r"^s3://(?<bucket>[^/]++)/?+(?<key>.*)$", s3_path)

    if m !== nothing
        return m[:bucket], m[:key]
    else
        throw(ArgumentError("String is not an S3 path: \"$s3_path\""))
    end
end

@testset "S3" begin
    @testset "presign" begin
        @testset "minimum period" begin
            @test_throws InexactError AWSTools.S3.presign(
                p"s3://bucket/file",
                Dates.Millisecond(999),
            )
        end
    end

    @testset "Syncing" begin
        @testset "Sync two local directories" begin
            # Create files to sync
            src = mktmpdir()
            src_file = join(src, "file1")
            write(src_file, "Hello World!")

            src_folder = join(src, "folder1")
            mkdir(src_folder)
            src_folder_file = join(src_folder, "file2")
            write(src_folder_file, "") # empty file

            src_folder2 = join(src_folder, "folder2") # nested folders
            mkdir(src_folder2)
            src_folder2_file = join(src_folder2, "file3")
            write(src_folder2_file, "Test")

            # Sync files
            dest = mktmpdir()
            @test isdir(dest)
            sync(src, dest)

            # Test directories are the same
            compare_dir(src, dest)

            # Get paths of new files
            dest_file = join(dest, "file1")
            dest_file_mtime = modified(dest_file)
            dest_folder = join(dest, basename(src_folder))
            dest_folder_file = join(dest_folder, "file2")
            dest_folder2 = join(dest_folder, basename(src_folder2))
            dest_folder2_file = join(dest_folder2, "file3")

            # Test that contents get copied over and size is equal
            compare(src_file, dest_file)
            compare(src_folder_file, dest_folder_file)
            compare(src_folder2_file, dest_folder2_file)

            compare_dir(src_folder, dest_folder)
            compare_dir(src_folder2, dest_folder2)

            @testset "Sync modified dest file" begin
                # Modify a file in dest
                write(dest_folder_file, "Modified in dest.")

                # Syncing overwrites the newer file in dest because it is of different size
                sync(src, dest)
                compare(src_folder_file, dest_folder_file)
            end

            @testset "Sync modified src file" begin
                 # Modify a file in src
                write(src_folder_file, "Modified in src.")

                # Test that syncing overwrites the modified file in dest
                sync(src, dest)
                compare(src_folder_file, dest_folder_file)

                # Test other files weren't updated
                @test modified(dest_file) == dest_file_mtime
            end

            @testset "Sync newer dest file" begin
                # This is the case because a newer file of the same size is usually the
                # result of an uploaded file always having a newer last_modified time.

                # Modify a file in dest
                write(dest_folder_file, "Modified in dest")

                # Test that syncing doesn't overwrite the newer file in dest
                sync(src, dest)
                @test read(dest_folder_file) != read(src_folder_file)
            end

            @testset "Sync incompatible types" begin
                @test_throws ArgumentError sync(src, dest_folder_file)
                @test_throws ArgumentError sync(src_folder_file, dest)
            end

            rm(src_file)

            @testset "Sync deleted file with no delete flag" begin
                # Syncing should not delete the file in the destination
                sync(src, dest)
                @test isfile(dest_file)
            end

            @testset "Sync deleted files with delete flag" begin
                # Test that syncing deletes the file in dest
                sync(src, dest; delete=true)
                @test !isfile(dest_file)

                rm(Path(src_folder2); recursive=true)

                @test isfile(dest_folder2_file)
                @test isdir(dest_folder2)

                sync(src, dest; delete=true)

                @test !isfile(dest_folder2_file)
                @test !isdir(dest_folder2)
                @test isdir(dest_folder)
            end

            @testset "Sync files" begin
                @test !isfile(dest_file)

                write(src_file, "Test")

                sync(src_file, dest_file)

                @test isfile(dest_file)
                compare(src_file, dest_file)
            end

            @testset "Sync empty directory" begin
                rm(src; recursive=true)

                @test_throws ArgumentError sync(src, dest, delete=true)
            end

            @testset "Sync non existent directories" begin
                isdir(src) && rm(src; recursive=true)
                isdir(dest) && rm(dest; recursive=true)

                # Test syncing creates non existent local directories
                @test_throws ArgumentError sync(src, dest)
            end
        end

        if "S3" in ONLINE
            @testset "Online" begin
                @info "Running ONLINE S3 tests"

                # Create bucket for tests
                test_run_id = string(uuid4())
                s3_prefix = if TEST_BUCKET_AND_PREFIX === nothing
                    bucket = string("awstools-s3-test-temp-", uuid4())
                    @info "Creating S3 bucket $bucket"
                    s3_create_bucket(bucket)
                    "s3://$bucket/$test_run_id"
                else
                    "$TEST_BUCKET_AND_PREFIX/$test_run_id"
                end

                try
                    @testset "Upload to S3" begin
                        dest = S3Path("$s3_prefix/folder3/testfile")

                        try
                            mktemp() do src, stream
                                write(stream, "Local file src")
                                close(stream)

                                @test !exists(dest)

                                uploaded_file = upload(Path(src), dest)
                                @test isa(uploaded_file, S3Path)

                                @test exists(dest)
                                @test isfile(dest)
                                @test isdir(parent(dest))
                                @test read(dest, String) == "Local file src"
                            end
                        finally
                            rm(dest; recursive=true)
                        end
                    end

                    @testset "Download from S3" begin
                        src = S3Path("$s3_prefix/folder4/testfile")

                        try
                            s3_put(src.bucket, src.key, "Remote content")

                            @testset "Download to a directory" begin
                                mktempdir() do dest_dir
                                    dest = Path(dest_dir)

                                    dest_file = download(src, dest)
                                    @test isa(dest_file, AbstractPath)

                                    @test exists(Path(dest_file))
                                    @test read(dest_file, String) == "Remote content"
                                end
                            end

                            @testset "Download to a local file" begin
                                mktemp() do dest_file, stream
                                    dest = Path(dest_file)
                                    close(stream)

                                    dest_file = download(src, dest)
                                    @test isa(dest_file, AbstractPath)

                                    @test dest_file == dest
                                    @test read(dest, String) == "Remote content"
                                end
                            end

                        finally
                            rm(src; recursive=true)
                        end
                    end

                    @testset "Download via presign" begin
                        src = S3Path("$s3_prefix/presign/file")
                        content = "presigned content"
                        s3_put(src.bucket, src.key, content)

                        @testset "file" begin
                            url = AWSTools.S3.presign(src, Dates.Minute(1))
                            r = HTTP.get(url)
                            @test String(r.body) == content
                        end

                        @testset "directory" begin
                            url = AWSTools.S3.presign(parent(src), Dates.Minute(1))
                            r = HTTP.get(url, status_exception=false)
                            @test r.status == 404
                            @test occursin("The specified key does not exist.", String(r.body))
                        end

                        @testset "expired" begin
                            url = AWSTools.S3.presign(src, Dates.Second(1))
                            sleep(2)
                            r = HTTP.get(url, status_exception=false)
                            @test r.status == 403
                            @test occursin("Request has expired", String(r.body))
                        end
                    end

                    @testset "Sync S3 directories" begin
                        bucket, key_prefix = bucket_and_key(s3_prefix)

                        src_dir = S3Path("$s3_prefix/folder1/")
                        dest_dir = S3Path("$s3_prefix/folder2/")

                        # Directories should be empty, but just in case
                        # delete any pre-existing objects in the s3 bucket directories
                        rm(src_dir; recursive=true)
                        rm(dest_dir; recursive=true)

                        # Make the src S3 directory
                        mkdir(src_dir; recursive=true)
                        mkdir(dest_dir; recursive=true)

                        # Note: using `lstrip` for when `key_prefix` is empty
                        # We also include the "folder" objects in this list.
                        s3_objects = [
                            Dict(
                                "Bucket" => bucket,
                                "Key" => lstrip("$key_prefix/folder1/", '/'),
                                "Content" => "",
                            ),
                            Dict(
                                "Bucket" => bucket,
                                "Key" => lstrip("$key_prefix/folder1/folder/", '/'),
                                "Content" => "",
                            ),
                            Dict(
                                "Bucket" => bucket,
                                "Key" => lstrip("$key_prefix/folder1/file1", '/'),
                                "Content" => "Hello World!",
                            ),
                            Dict(
                                "Bucket" => bucket,
                                "Key" => lstrip("$key_prefix/folder1/file2", '/'),
                                "Content" => "",
                            ),
                            Dict(
                                "Bucket" => bucket,
                                "Key" => lstrip("$key_prefix/folder1/folder/file3", '/'),
                                "Content" => "Test",
                            ),
                        ]

                        # Set up the source s3 directory
                        for object in s3_objects
                            s3_put(object["Bucket"], object["Key"], object["Content"])
                        end

                        # Sync files passing in directories as strings
                        sync(string(src_dir), string(dest_dir))

                        src_files = collect(walkpath(src_dir))
                        dest_files = collect(walkpath(dest_dir))

                        # Test destination directory has expected files
                        @test !isempty(dest_files)
                        @test length(dest_files) == length(src_files)

                        # Test directories are the same
                        compare_dir(src_dir, dest_dir)

                        # Test readdir only lists files and "dirs" within this S3 "dir"
                        @test readdir(dest_dir) == ["file1", "file2", "folder/"]
                        @test readdir(S3Path("$s3_prefix/")) == [
                            "folder1/", "folder2/", "presign/"
                        ]
                        # Not including the ending `/` means this refers to an object and
                        # not a directory prefix in S3
                        @test_throws ArgumentError readdir(S3Path(s3_prefix))

                        @testset "Sync modified dest file" begin
                            # Modify a file in dest
                            s3_put(
                                dest_files[1].bucket,
                                dest_files[1].key,
                                "Modified in dest.",
                            )

                            # Syncing overwrites the newer file in dest because it is of
                            # a different size
                            sync(src_dir, dest_dir)

                            # Test directories are the same
                            compare_dir(src_dir, dest_dir)

                            @test read(dest_files[1], String) == s3_objects[3]["Content"]
                        end

                        @testset "Sync modified src file" begin
                            # Modify a file in src
                            s3_objects[3]["Content"] = "Modified in src."
                            s3_put(
                                s3_objects[3]["Bucket"],
                                s3_objects[3]["Key"],
                                s3_objects[3]["Content"],
                            )

                            # Test that syncing overwrites the modified file in dest
                            sync(src_dir, dest_dir)

                            # Test directories are the same
                            compare_dir(src_dir, dest_dir)

                            @test read(dest_files[1], String) == s3_objects[3]["Content"]
                        end

                        @testset "Sync newer file in dest" begin
                            # Modify a file in dest
                            s3_put(
                                dest_files[1].bucket,
                                dest_files[1].key,
                                "Modified in dest",
                            )

                            # Test that syncing doesn't overwrite the newer file in dest
                            # This is the case because a newer file of the same size is
                            # usually the result of an uploaded file always having a newer
                            # last_modified time.
                            sync(src_dir, dest_dir)

                            file_contents = read(dest_files[1], String)
                            @test file_contents != s3_objects[3]["Content"]
                            @test file_contents ==  "Modified in dest"
                        end

                        @testset "Sync s3 bucket with object prefix" begin
                            obj = s3_objects[3]
                            file = "s3://" * join([obj["Bucket"], obj["Key"]], '/')

                            @test startswith(file, "s3://")
                            @test_throws ArgumentError sync(S3Path(file), dest_dir)
                        end

                        rm(src_files[1])

                        @testset "Sync deleted file with no delete flag" begin
                            sync(src_dir, dest_dir)

                            # We collect the walkpath result because we want it to complete
                            # before checking the length
                            src_files = collect(walkpath(src_dir))
                            dest_files = collect(walkpath(dest_dir))

                            @test length(dest_files) == length(src_files) + 1
                        end

                        @testset "Sync deleted files with delete flag" begin
                            # Test that syncing deletes the file in dest
                            sync(src_dir, dest_dir; delete=true)

                            # We collect the walkpath result because we want it to complete
                            # before checking the length
                            src_files = collect(walkpath(src_dir))
                            dest_files = collect(walkpath(dest_dir))
                            @test length(dest_files) == length(src_files)
                        end

                        @testset "Sync files" begin
                            src_file = S3Path("$s3_prefix/folder1/file")
                            dest_file = S3Path("$s3_prefix/folder2/file")

                            # Modify file in src dir
                            write(src_file, "Test modified file in source")

                            # Test syncing individual files instead of directories
                            sync(src_file, dest_file)

                            # Test directories are the same
                            compare_dir(src_dir, dest_dir)

                            # Test destination file was modifies
                            @test read(dest_file, String) == "Test modified file in source"
                        end

                        @testset "Sync non-existent directory" begin
                            rm(src_dir; recursive=true)

                            # Error because src folder doesn't exist
                            @test_throws ArgumentError sync(src_dir, dest_dir; delete=true)

                            src_files = collect(walkpath(src_dir))
                            @test isempty(src_files)

                            rm(dest_dir; recursive=true)
                            dest_files = collect(walkpath(dest_dir))
                            @test isempty(dest_files)
                        end
                    end

                finally
                    # Clean up any files left in the test directory
                    rm(S3Path("$s3_prefix/"); recursive=true)

                    # Delete bucket if it was explicitly created
                    if TEST_BUCKET_AND_PREFIX === nothing
                        s3_bucket_dir = replace(s3_prefix, r"^(s3://[^/]+).*$" => s"\1")
                        bucket, key = bucket_and_key(s3_bucket_dir)
                        @info "Deleting S3 bucket $bucket"
                        rm(S3Path("s3://$bucket/"); recursive=true)
                        s3_delete_bucket(bucket)
                    end
                end
            end
        else
            @warn """
                Skipping AWSTools.S3 ONLINE tests. Set `ENV["ONLINE"] = "S3"` to run."
                Can also optionally specify a test bucket name `ENV["TEST_BUCKET_AND_PREFIX"] = "s3://bucket/dir/"`.
                If `TEST_BUCKET_AND_PREFIX` is not specified, a temporary bucket will be created, used, and then deleted.
                """
        end
    end
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLinSysSolveRoutine

using Test, RLinearAlgebra

@testset "Linear System Solver Routine Abstractions" begin

    # Verify parent types
    @test supertype(LinSysVecRowProjection) == LinSysSolveRoutine
    @test supertype(LinSysVecColProjection) == LinSysSolveRoutine
    @test supertype(LinSysBlkRowProjection) == LinSysSolveRoutine
    @test supertype(LinSysBlkColProjection) == LinSysSolveRoutine
    @test supertype(LinSysPreconKrylov) == LinSysSolveRoutine

end

end # End module


text:
using Test;

@testset "BusinessModel Tests" begin
   @testset "add_users $i" for i in 1:3
       u = 2;
       r = 3;
       m = CyberSecurity.BusinessModel(u, r, .02);
       CyberSecurity.add_users(m, i);
       @test CyberSecurity.count_users(m) == u + i;
   end
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
@testset "New Function Tests" begin

  @test isdefined(Julz, :new) == true

  cd("tmp") do

    package_name = "NewTestPackage.jl"

    rm(package_name, force=true, recursive=true)

    originalSTDERR = STDERR

    (errRead, errWrite) = redirect_stderr()

    @test_throws Pkg.PkgError Julz.new(package_name)

    license = "MIT"

    Julz.new(package_name, license)

    close(errWrite)

    close(errRead)

    redirect_stderr(originalSTDERR)

    cd(package_name) do

      @test isfile("src/$package_name")

    end

    rm(package_name, force=true, recursive=true)

  end

end