text:
### This File contains all compression routines for the HssMatrices.jl package
#
# Direct compression
# as described in
# Xia, J., Chandrasekaran, S., Gu, M., & Li, X. S. (2010). Fast algorithms for hierarchically semiseparable matrices.
# Numerical Linear Algebra with Applications, 17(6), 953–976. https://doi.org/10.1002/nla.691
#
# Recompression
# as described in
# Xia, J., Chandrasekaran, S., Gu, M., & Li, X. S. (2010). Fast algorithms for hierarchically semiseparable matrices.
# Numerical Linear Algebra with Applications, 17(6), 953–976. https://doi.org/10.1002/nla.691
#
# Randomized compression
# as described in
# Martinsson, P. G. (2011). A Fast Randomized Algorithm for Computing a Hierarchically Semiseparable Representation of a Matrix.
# SIAM Journal on Matrix Analysis and Applications, 32(4), 1251–1274. https://doi.org/10.1137/100786617
#
# Re-Written by Boris Bonev, Jan. 2021

# Utility routine to provide access to pivoted rank-revealing qr
function _compress_block!(A::AbstractMatrix{T}, atol::Float64, rtol::Float64) where T
  #Q, R, p = prrqr!(A, atol, rtol)
  #rk = min(size(R)...)
  #return Q[:,1:rk], R[1:rk, invperm(p)]
  # temporarily using prrqr of LowRankApprox.jl - may be replaced in the future to reduce dependencies
  F = pqrfact(A; atol = atol, rtol = rtol, sketch=:none, pqrfact_retval = "qr")
  #rk = min(size(F.R)...)
  return F.Q, F.R[:, invperm(F.p)]
end

"""
  compress(A, rcl, ccl; args...)

Direct HSS compression.

# Arguments
* `A`: Matrix to be compressed
* `rcl`: row-cluster tree
* `ccl`: column-cluster tree

# Examples
```julia
julia> hssA = compress(A, rcl, ccl)
```
"""
function compress(A::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, opts::HssOptions=HssOptions(T); args...) where T
  opts = copy(opts; args...)
  chkopts!(opts)
  compatible(rcl, ccl) || throw(ArgumentError("row and column clusters are not compatible"))
  m = length(rcl.data); n = length(ccl.data)
  if size(A) != (m,n) throw(ArgumentError("size of row- and column-cluster-trees must match")) end
  Brow = Array{T}(undef, m, 0)
  Bcol = Array{T}(undef, 0, n)
  if isleaf(rcl) && isleaf(ccl)
    hssA, _, _ = _compress!(A, Brow, Bcol, rcl.data, ccl.data, opts.atol, opts.rtol; rootnode=true)
  elseif isbranch(rcl) && isbranch(ccl)
    hssA, _, _ = _compress!(A, Brow, Bcol, rcl, ccl, opts.atol, opts.rtol; rootnode=true)
  else
    throw(ArgumentError("row and column clusters are not compatible"))
  end
  return hssA
end

# leaf node function for compression
function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rows::UnitRange{Int}, cols::UnitRange{Int}, atol::Float64, rtol::Float64; rootnode=false) where T
  if rootnode
    return HssMatrix(A[rows, cols]), Brow, Bcol
  else
    U, Brow = _compress_block!(Brow, atol, rtol)
    V, Bcol = _compress_block!(Bcol', atol, rtol) #TODO: write code that is better at dealing with Julia's lazy transpose
    return HssMatrix(A[rows, cols], U, V), Brow, copy(Bcol')
  end
end

# branch node function for compression
function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, atol::Float64, rtol::Float64; rootnode=false) where T
  m1 = length(rcl.left.data); m2 = length(rcl.right.data)
  n1 = length(ccl.left.data); n2 = length(ccl.right.data)

  rows1 = rcl.left.data; rows2 = rcl.right.data
  cols1 = ccl.left.data; cols2 = ccl.right.data

  # left node
  Brow1 = [A[rows1, cols2]  Brow[1:m1, :]]
  Bcol1 = [A[rows2, cols1]; Bcol[:, 1:n1]]
  if isleaf(rcl.left) && isleaf(ccl.left)
    A11, Brow1, Bcol1 = _compress!(A, Brow1, Bcol1, rows1, cols1, atol, rtol)
  elseif isbranch(rcl.left) && isbranch(ccl.left)
    A11, Brow1, Bcol1 = _compress!(A, Brow1, Bcol1, rcl.left, ccl.left, atol, rtol)
  else
    throw(ArgumentError("row and column clusters are not compatible"))
  end

  # right node
  Brow2 = [Bcol1[1:m2, :]  Brow[m1+1:end, :]]
  Bcol2 = [Brow1[:, 1:n2]; Bcol[:, n1+1:end]]
  if isleaf(rcl.right) && isleaf(ccl.right)
    A22, Brow2, Bcol2 = _compress!(A, Brow2, Bcol2, rows2, cols2, atol, rtol)
  elseif isbranch(rcl.right) && isbranch(ccl.right)
    A22, Brow2, Bcol2 = _compress!(A, Brow2, Bcol2, rcl.right, ccl.right, atol, rtol)
  else
    throw(ArgumentError("row and column clusters are not compatible"))
  end

  # figure out compressed off-diagonal blocks
  rm1 = size(Brow1,1); rn1 = size(Bcol1, 2)
  rm2 = size(Brow2,1); rn2 = size(Bcol2, 2)
  B12 = Bcol2[1:rm1, :]
  B21 = Brow2[:, 1:rn1]

  # clean up stuff from the front and form the composed HSS block row/col for compression
  Brow = [Brow1[:, n2+1:end]; Brow2[:, rn1+1:end]]
  Bcol = [Bcol1[m2+1:end, :]  Bcol2[rm1+1:end, :]]

  if !rootnode
    # do the actual compression and disentangle blocks of the translation operators
    R, Brow = _compress_block!(Brow, atol, rtol)
    R1 = R[1:rm1, :]
    R2 = R[rm1+1:end, :]

    X = copy(Bcol')
    W, Bcol = _compress_block!(copy(Bcol'), atol, rtol);
    Bcol = copy(Bcol')
    W1 = W[1:rn1, :]
    W2 = W[rn1+1:end, :]

    hssA = HssMatrix(A11, A22, B12, B21, R1, W1, R2, W2, rootnode)
  else
    hssA = HssMatrix(A11, A22, B12, B21, rootnode)
  end
  return hssA, Brow, Bcol
end

"""
  recompress!(hssA; args...)

Recompression of an HSS matrix

# Arguments
* `hssA`: HSS matrix for recompression

# Examples
```julia
julia> hssA = recompress!(hssA, atol=1e-3, rtol=1e-3)
```
"""
function recompress!(hssA::HssMatrix{T}, opts::HssOptions=HssOptions(T); args...) where T
  opts = copy(opts; args...)
  chkopts!(opts)
  rtol = opts.rtol; atol = opts.atol;

  if isleaf(hssA); return hssA; end

  # a prerequisite for this algorithm to work is that generators are orthonormal
  orthonormalize_generators!(hssA)

  # compress B12, B21 via something that resembles the SVD
  P1, S2 = _compress_block!(hssA.B12, atol, rtol)
  Q2, T1 = _compress_block!(copy(hssA.B12'), atol, rtol)
  P2, S1 = _compress_block!(hssA.B21, atol, rtol)
  Q1, T2 = _compress_block!(copy(hssA.B21'), atol, rtol)

  hssA.B12 = S2*Q2
  hssA.B21 = S1*Q1

  # fix dimensions of ghost-translators in rootnode
  hssA.R1 = hssA.R1[1:size(hssA.B12, 1),:]
  hssA.W1 = hssA.W1[1:size(hssA.B21, 2),:]
  hssA.R2 = hssA.R2[1:size(hssA.B21, 1),:]
  hssA.W2 = hssA.W2[1:size(hssA.B12, 2),:]

  # pass information to children and proceed recursively
  if isbranch(hssA.A11)
    hssA.A11.R1 = hssA.A11.R1*P1
    hssA.A11.R2 = hssA.A11.R2*P1
    hssA.A11.W1 = hssA.A11.W1*Q1
    hssA.A11.W2 = hssA.A11.W2*Q1
  elseif isleaf(hssA.A11)
    hssA.A11.U = hssA.A11.U*P1
    hssA.A11.V = hssA.A11.V*Q1
  end
  # repeat for A22
  if isbranch(hssA.A22)
    hssA.A22.R1 = hssA.A22.R1*P2
    hssA.A22.R2 = hssA.A22.R2*P2
    hssA.A22.W1 = hssA.A22.W1*Q2
    hssA.A22.W2 = hssA.A22.W2*Q2
  elseif isleaf(hssA.A22)
    hssA.A22.U = hssA.A22.U*P2
    hssA.A22.V = hssA.A22.V*Q2
  end

  # call recompression recursively
  if isbranch(hssA.A11)
    _recompress!(hssA.A11, hssA.B12, copy(hssA.B21'), atol, rtol)
  end
  if isbranch(hssA.A22)
    _recompress!(hssA.A22, hssA.B21, copy(hssA.B12'), atol, rtol)
  end

  return hssA
end

function _recompress!(hssA::HssMatrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, atol, rtol) where T
  # compress B12
  Brow1 = [hssA.B12  hssA.R1*Brow]
  Bcol2 = [hssA.B12' hssA.W2*Bcol]
  P1, S2 = _compress_block!(Brow1, atol, rtol)
  Q2, T1 = _compress_block!(Bcol2, atol, rtol)
  # get the original number of columns in B12
  rm1, rn2 = size(hssA.B12)
  hssA.B12 = S2[:,1:rn2]*Q2
  hssA.R1 = P1'*hssA.R1
  hssA.W2 = Q2'*hssA.W2
  # compress B21
  Brow2 = [hssA.B21  hssA.R2*Brow]
  Bcol1 = [hssA.B21' hssA.W1*Bcol]
  P2, S1 = _compress_block!(Brow2, atol, rtol)
  Q1, T2 = _compress_block!(Bcol1, atol, rtol)
  # get the original number of columns in B21
  rm2, rn1 = size(hssA.B21)
  hssA.B21 = S1[:,1:rn1]*Q1
  hssA.R2 = P2'*hssA.R2
  hssA.W1 = Q1'*hssA.W1
  # update generators of A11
  if isbranch(hssA.A11)
    hssA.A11.R1 = hssA.A11.R1*P1
    hssA.A11.R2 = hssA.A11.R2*P1
    hssA.A11.W1 = hssA.A11.W1*Q1
    hssA.A11.W2 = hssA.A11.W2*Q1
  elseif isleaf(hssA.A11)
    hssA.A11.U = hssA.A11.U*P1
    hssA.A11.V = hssA.A11.V*Q1
  end
  # update generators in A22
  if isbranch(hssA.A22)
    hssA.A22.R1 = hssA.A22.R1*P2
    hssA.A22.R2 = hssA.A22.R2*P2
    hssA.A22.W1 = hssA.A22.W1*Q2
    hssA.A22.W2 = hssA.A22.W2*Q2
  else isleaf(hssA.A22)
    hssA.A22.U = hssA.A22.U*P2
    hssA.A22.V = hssA.A22.V*Q2
  end

  # call recompression recursively
  if isbranch(hssA.A11)
    Brow1 = hcat(hssA.B12, S2[:,rn2+1:end])
    Bcol1 = hcat(hssA.B21', T2[:,rm2+1:end])
    _recompress!(hssA.A11, Brow1, Bcol1, atol, rtol)
  end
  if isbranch(hssA.A22)
    Brow2 = hcat(hssA.B21, S1[:,rn1+1:end])
    Bcol2 = hcat(hssA.B12', T1[:,rm1+1:end])
    _recompress!(hssA.A22, Brow2, Bcol2, atol, rtol)
  end

  return hssA
end

"""
  randcompress(A, rcl, ccl, kest; args...)

Randomized HSS compression.

# Arguments
* `A`: Either a Matrix or a LinearMap type, to be compressed
* `rcl`: row-cluster tree
* `ccl`: column-cluster tree
* `kest`: rank estimate determines the size of the sampling matrix, as well as maximum off-diagonal rank

# Examples
```julia
julia> hssA = randcompress(A, rcl, ccl, kest)
```
"""
function randcompress(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree, kest::Int, opts::HssOptions=HssOptions(T), args...) where T
  opts = copy(opts; args...)
  chkopts!(opts)
  m, n = size(A)
  compatible(rcl, ccl) || throw(ArgumentError("row and column clusters are not compatible"))
  #if typeof(A) <: AbstractMatrix A = LinOp(A) end
  hssA = hss_blkdiag(A, rcl, ccl)

  # compute initial sampling
  k = kest; r = opts.noversampling;
  Ωcol = randn(n, k+r)
  Ωrow = randn(m, k+r)
  Scol = A*Ωcol # this should invoke the magic of the linearoperator.jl type
  Srow = A'*Ωrow
  hssA, _, _, _, _, _, _, _, _ = _randcompress!(hssA, A, Scol, Srow, Ωcol, Ωrow, 0, 0, opts.atol, opts.rtol; rootnode=true)
  return hssA
end

"""
  randcompress_adaptive(A, rcl, ccl; kest=10 args...)

Randomized HSS compression.

# Arguments
* `A`: Either a Matrix or a LinearMap type, to be compressed
* `rcl`: row-cluster tree
* `ccl`: column-cluster tree
* `kest`: rank estimate determines the size of the sampling matrix, as well as maximum off-diagonal rank

# Examples
```julia
julia> hssA = randcompress_adaptive(A, rcl, ccl, kest=20)
```
"""
function randcompress_adaptive(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree, opts=HssOptions(T); kest=10, args...) where T
  opts = copy(opts; args...)
  chkopts!(opts)
  m, n = size(A)
  maxrank = n
  compatible(rcl, ccl) || throw(ArgumentError("row and column clusters are not compatible"))
  #if typeof(A) <: AbstractMatrix A = LinMap(A) end
  hssA = hss_blkdiag(A, rcl, ccl)

  # compute initial sampling
  k = kest; r = opts.noversampling; bs = opts.stepsize
  #bs = Integer(ceil(n*0.01)) # this should probably be better estimated
  Ωcol = randn(n, k+r)
  Ωrow = randn(m, k+r)
  Scol = A*Ωcol # this should invoke the magic of the linearoperator.jl type
  Srow = A'*Ωrow
  failed = true

  while failed && k < maxrank
    # TODO: In further versions we might want to re-use the information gained during previous attempts
    hssA, _, _, _, _, _, _, _, _  = _randcompress!(hssA, A, Scol, Srow, Ωcol, Ωrow, 0, 0, opts.atol, opts.rtol; rootnode=true)

    Ωcol_test = randn(n, bs)
    Ωrow_test = randn(m, bs)
    Scol_test = A*Ωcol_test
    Srow_test = A'*Ωrow_test

    #@infiltrate
    nrm = sqrt(1/bs)*norm(Scol_test)
    nrm_est = sqrt(1/bs)*norm(Scol_test - hssA*Ωcol_test)
    failed = nrm_est > opts.atol && nrm_est > opts.rtol*nrm

    if failed
      opts.verbose && println("Enlarging sampling space to ", k+bs)
      Ωcol = [Ωcol Ωcol_test]
      Ωrow = [Ωrow Ωrow_test]
      Scol = [Scol Scol_test]
      Srow = [Srow Srow_test]
      k = k+bs
    end
  end

  return hssA
end

# this function compresses given the sampling matrix of rank k
function _randcompress!(hssA::HssMatrix, A::AbstractMatOrLinOp, Scol::Matrix, Srow::Matrix, Ωcol::Matrix, Ωrow::Matrix, ro::Int, co::Int, atol::Float64, rtol::Float64; rootnode=false)
  if isleaf(hssA)
    Scol .= Scol .- hssA.D * Ωcol
    Srow .= Srow .- hssA.D' * Ωrow

    #if rootnode return hssA, Scol, Srow, Ωcol, Ωrow, Jcol, Jrow, hssA.U, hssA.V  end

    # take care of column-space
    Xcol, Jcol = _interpolate(Scol', atol, rtol)
    hssA.U = Xcol'
    Scol = Scol[Jcol, :]
    U = Xcol'
    Jcol .= ro .+ Jcol

    # same for the row-space
    Xrow, Jrow = _interpolate(Srow', atol, rtol)
    hssA.V = Xrow'
    Srow = Srow[Jrow, :]
    V = Xrow'
    Jrow .= co .+ Jrow
  else
    m1, n1 = hssA.sz1; m2, n2 = hssA.sz2
    hssA.A11, Scol1, Srow1, Ωcol1, Ωrow1, Jcol1, Jrow1, U1, V1 = _randcompress!(hssA.A11, A, Scol[1:m1, :], Srow[1:n1, :], Ωcol[1:n1, :], Ωrow[1:m1, :], ro, co, atol, rtol)
    hssA.A22, Scol2, Srow2, Ωcol2, Ωrow2, Jcol2, Jrow2, U2, V2 = _randcompress!(hssA.A22, A, Scol[m1+1:end, :], Srow[n1+1:end, :], Ωcol[n1+1:end, :], Ωrow[m1+1:end, :], ro+m1, co+n1, atol, rtol)

    # update the sampling matrix based on the extracted generators
    Ωcol2 = V2' * Ωcol2
    Ωcol1 = V1' * Ωcol1
    Ωrow2 = U2' * Ωrow2
    Ωrow1 = U1' * Ωrow1
    # step 1 in the algorithm: look only at extracted rows/cols
    Jcol = [Jcol1; Jcol2]
    Jrow = [Jrow1; Jrow2]
    Ωcol = [Ωcol1; Ωcol2]
    Ωrow = [Ωrow1; Ωrow2]
    # extract the correct generator blocks
    hssA.B12 = A[Jcol1, Jrow2]
    hssA.B21 = A[Jcol2, Jrow1]
    # subtract the diagonal block
    Scol = [Scol1 .- hssA.B12  * Ωcol2;  Scol2 .- hssA.B21  * Ωcol1 ];
    Srow = [Srow1 .- hssA.B21' * Ωrow2;  Srow2 .- hssA.B12' * Ωrow1 ];

    if rootnode
      rk1, wk1 = gensize(hssA.A11)
      rk2, wk2 = gensize(hssA.A22)
      hssA.R1 = Matrix{eltype(hssA)}(undef, rk1, 0)
      hssA.R2 = Matrix{eltype(hssA)}(undef, rk2, 0)
      hssA.W1 = Matrix{eltype(hssA)}(undef, wk1, 0)
      hssA.W2 = Matrix{eltype(hssA)}(undef, wk2, 0)
      hssA.rootnode = true

      U = similar(U1, rk1+rk2, 0)
      V = similar(V1, wk1+wk2, 0)
    else
      # take care of the columns
      Xcol, Jcolloc = _interpolate(Scol', atol, rtol)
      hssA.R1 = Xcol[:, 1:size(Scol1, 1)]'
      hssA.R2 = Xcol[:, size(Scol1, 1)+1:end]'
      Scol = Scol[Jcolloc, :]
      Jcol = Jcol[Jcolloc]
      U = [hssA.R1; hssA.R2]

      # take care of the rows
      Xrow, Jrowloc = _interpolate(Srow', atol, rtol)
      hssA.W1 = Xrow[:, 1:size(Srow1, 1)]'
      hssA.W2 = Xrow[:, size(Srow1, 1)+1:end]'
      Srow = Srow[Jrowloc, :]
      Jrow = Jrow[Jrowloc]
      V = [hssA.W1; hssA.W2]
    end
  end
  return hssA, Scol, Srow, Ωcol, Ωrow, Jcol, Jrow, U, V
end

# interpolative decomposition
# still gotta figure out which qr decomposition to use
function _interpolate(A::AbstractMatrix{T}, atol::Float64, rtol::Float64) where T
  size(A,2) == 0 && return Matrix{eltype(A)}(undef, 0, 0), Vector{Int}()
  #_, R, p = prrqr(A, tol; reltol=reltol)
  #rk = min(size(R)...)
  _, R, p  = qr(A, Val(true))
  tol = min( atol, rtol * abs(R[1,1]) )
  rk = sum(abs.(diag(R)) .> tol)
  J = p[1:rk]
  X = R[1:rk, 1:rk]\R[1:rk,invperm(p)]
  return X, J
end

# extracts the block-diagonal of a matrix as HSS matrix of rank 0
function hss_blkdiag(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree; rootnode=true) where T
  m = length(rcl.data); n = length(ccl.data)
  if isleaf(rcl) # only check row cluster as we have already checked cluster equality

    D = convert(Matrix{T}, A[rcl.data, ccl.data])
    if rootnode
      hssA = HssMatrix(D)
    else
      hssA = HssMatrix(D, Matrix{T}(undef,m,0), Matrix{T}(undef,n,0))
    end
  elseif isbranch(rcl)
    A11 = hss_blkdiag(A, rcl.left, ccl.left; rootnode=false)
    A22 = hss_blkdiag(A, rcl.right, ccl.right; rootnode=false)
    B12 = Matrix{T}(undef,0,0)
    B21 = Matrix{T}(undef,0,0)
    if rootnode
      hssA = HssMatrix(A11, A22, B12, B21, true)
    else
      R1 = Matrix{T}(undef,0,0)
      W1 = Matrix{T}(undef,0,0)
      R2 = Matrix{T}(undef,0,0)
      W2 = Matrix{T}(undef,0,0)
      hssA = HssMatrix(A11, A22, B12, B21, R1, W1, R2, W2, )
    end
  else
    throw(ErrorException("Encountered node with only one child. Make sure that each node in the binary tree has either two children or is a leaf."))
  end
  return hssA
end


text:
### Abstract Interface

Base.length(sol::DESolution) = length(sol.u) # Must be on u for the test solutions!
Base.size(sol::DESolution,n::Int) = n==1 ? length(sol.u) : error("Only dimension 1 has a well-defined size.")
Base.endof(sol::DESolution) = length(sol)
Base.getindex(sol::DESolution,i::Int) = sol.u[i]
Base.getindex(sol::DESolution,c::AbstractArray) = [sol.u[j] for j in c]
Base.getindex(sol::DESolution,i::Int,I::Int...) = sol.u[i][I...]
Base.getindex(sol::DESolution,i::Int,::Colon,I::Int) = sol.u[i][:,I]
Base.getindex(sol::DESolution,i::Int,::Colon,I::Int...) = sol.u[i][:,I...]
Base.getindex(sol::DESolution,::Colon) = sol.u
Base.getindex(sol::DESolution,::Colon,i::Int...) = [sol.u[j][i...] for j in eachindex(sol)]
Base.getindex(sol::DESolution,c::AbstractArray,i::Int...) =  [sol.u[j][i...] for j in c]
eachindex(sol::DESolution) = eachindex(sol.t)
tuples(sol::DESolution) = tuple.(sol.t,sol.u)

function start(sol::DESolution)
  sol.tslocation = 1
  1
end

function next(sol::DESolution,state)
  state += 1
  sol.tslocation = state
  (sol,state)
end

function done(sol::DESolution,state)
  state >= length(sol)
end

function eltype(sol::DESolution)
  if typeof(sol[1]) <: AbstractArray
    return typeof(sol[1][1])
  else
    return typeof(sol[1])
  end
end

const DEFAULT_PLOT_FUNC = (x...) -> (x...)

@recipe function f(sol::DESolution;
                   plot_analytic=false,
                   denseplot = (sol.dense || typeof(sol.prob) <: AbstractDiscreteProblem) && !(typeof(sol) <: AbstractRODESolution),
                   plotdensity = sol.tslocation==0 ? (typeof(sol.prob) <: AbstractDiscreteProblem ? 100*length(sol) : 10*length(sol)) : 100*sol.tslocation,
                   vars=nothing)

  int_vars = interpret_vars(vars,sol)
  if sol.tslocation == 0
    end_idx = length(sol)
  else
    end_idx = sol.tslocation
  end

  if denseplot
    # Generate the points from the plot from dense function
    plott = collect(Ranges.linspace(sol.t[1],sol.t[end_idx],plotdensity))
    plot_timeseries = sol(plott)
    if plot_analytic
      plot_analytic_timeseries = [sol.prob.f(Val{:analytic},t,sol.prob.u0) for t in plott]
    else
      plot_analytic_timeseries = nothing
    end
  else
    # Plot for sparse output: use the timeseries itself
    if sol.tslocation == 0
      plott = sol.t
      plot_timeseries = sol.u
      if plot_analytic
        plot_analytic_timeseries = sol.u_analytic
      else
        plot_analytic_timeseries = nothing
      end
    else
      plott = sol.t[1:end_idx]
      plot_timeseries = sol.u[1:end_idx]
      if plot_analytic
        plot_analytic_timeseries = sol.u_analytic[1:end_idx]
      else
        plot_analytic_timeseries = nothing
      end
    end
  end

  dims = length(int_vars[1])
  for var in int_vars
    @assert length(var) == dims
  end
  # Should check that all have the same dims!
  plot_vecs,labels = solplot_vecs_and_labels(dims,int_vars,plot_timeseries,plott,sol,plot_analytic,plot_analytic_timeseries)

  tdir = sign(sol.t[end]-sol.t[1])
  xflip --> tdir < 0
  seriestype --> :path
  linewidth --> 3

  # Special case labels when vars = (:x,:y,:z) or (:x) or [:x,:y] ...
  if typeof(vars) <: Tuple && (typeof(vars[1]) == Symbol && typeof(vars[2]) == Symbol)
    xlabel --> vars[1]
    ylabel --> vars[2]
    if length(vars) > 2
      zlabel --> vars[3]
    end
  end
  if getindex.(int_vars,1) == zeros(length(int_vars)) || getindex.(int_vars,2) == zeros(length(int_vars))
    xlabel --> "t"
  end
  label --> reshape(labels,1,length(labels))
  (plot_vecs...)
end

function interpret_vars(vars,sol)
  if vars != nothing && has_syms(sol.prob.f)
    # Do syms conversion
    tmp_vars = []
    for var in vars
      if typeof(var) <: Symbol
        var_int = findfirst(sol.prob.f.syms,var)
      elseif typeof(var) <: Union{Tuple,AbstractArray} #eltype(var) <: Symbol # Some kind of iterable
        tmp = []
        for x in var
          if typeof(x) <: Symbol
            push!(tmp,findfirst(sol.prob.f.syms,x))
          else
            push!(tmp,x)
          end
        end
        if typeof(var) <: Tuple
          var_int = tuple(tmp...)
        else
          var_int = tmp
        end
      else
        var_int = var
      end
      push!(tmp_vars,var_int)
    end
    if typeof(vars) <: Tuple
      vars = tuple(tmp_vars...)
    else
      vars = tmp_vars
    end
  end

  if vars == nothing
    # Default: plot all timeseries
    if typeof(sol[1]) <: Union{Tuple,AbstractArray}
      vars = collect((DEFAULT_PLOT_FUNC,0, i) for i in plot_indices(sol[1]))
    else
      vars = [(DEFAULT_PLOT_FUNC,0, 1)]
    end
  end

  if typeof(vars) <: Integer
    vars = [(DEFAULT_PLOT_FUNC,0, vars)]
  end

  if typeof(vars) <: AbstractArray
    # If list given, its elements should be tuples, or we assume x = time
    tmp = Tuple[]
    for x in vars
      if typeof(x) <: Tuple
        if typeof(x[1]) <: Int
          push!(tmp,tuple(DEFAULT_PLOT_FUNC,x...))
        else
          push!(tmp,x)
        end
      else
        push!(tmp,(DEFAULT_PLOT_FUNC,0, x))
      end
    end
    vars = tmp
  end

  if typeof(vars) <: Tuple
    # If tuple given...
    if typeof(vars[end-1]) <: AbstractArray
      if typeof(vars[end]) <: AbstractArray
        # If both axes are lists we zip (will fail if different lengths)
        vars = collect(zip([DEFAULT_PLOT_FUNC for i in eachindex(vars[end-1])],vars[end-1], vars[end]))
      else
        # Just the x axis is a list
        vars = [(DEFAULT_PLOT_FUNC,x, vars[end]) for x in vars[end-1]]
      end
    else
      if typeof(vars[2]) <: AbstractArray
        # Just the y axis is a list
        vars = [(DEFAULT_PLOT_FUNC,vars[end-1], y) for y in vars[end]]
      else
        # Both axes are numbers
        if typeof(vars[1]) <: Int
          vars = [tuple(DEFAULT_PLOT_FUNC,vars...)]
        else
          vars = [vars]
        end
      end
    end
  end

  # Here `vars` should be a list of tuples (x, y).
  assert(typeof(vars) <: AbstractArray)
  assert(eltype(vars) <: Tuple)
  vars
end

function add_labels!(labels,x,dims,sol)
  lys = []
  for j in 3:dims
    if x[j] == 0
      push!(lys,"t,")
    else
      if has_syms(sol.prob.f)
        push!(lys,"$(sol.prob.f.syms[x[j]]),")
      else
        push!(lys,"u$(x[j]),")
      end
    end
  end
  lys[end] = lys[end][1:end-1] # Take off the last comma
  if x[2] == 0
    tmp_lab = "$(lys...)(t)"
  else
    if has_syms(sol.prob.f)
      tmp = sol.prob.f.syms[x[2]]
      tmp_lab = "($tmp,$(lys...))"
    else
      tmp_lab = "(u$(x[2]),$(lys...))"
    end
  end
  if x[1] != DEFAULT_PLOT_FUNC
    push!(labels,"$(x[1])$(tmp_lab)")
  else
    push!(labels,tmp_lab)
  end
  labels
end

function add_analytic_labels!(labels,x,dims,sol)
  lys = []
  for j in 3:dims
    if x[j] == 0
      push!(lys,"t,")
    else
      if has_syms(sol.prob.f)
        push!(lys,string("True ",sol.prob.f.syms[x[j]],","))
      else
        push!(lys,"True u$(x[j]),")
      end
    end
  end
  lys[end] = lys[end][1:end-1] # Take off the last comma
  if x[2] == 0
    tmp_lab = "$(lys...)(t)"
  else
    if has_syms(sol.prob.f)
      tmp = string("True ",sol.prob.f.syms[x[2]])
      tmp_lab = "($tmp,$(lys...))"
    else
      tmp_lab = "(True u$(x[2]),$(lys...))"
    end
  end
  if x[1] != DEFAULT_PLOT_FUNC
    push!(labels,"$(x[1])$(tmp_lab)")
  else
    push!(labels,tmp_lab)
  end
end

function u_n(timeseries::AbstractArray, n::Int,sol,plott,plot_timeseries)
  # Returns the nth variable from the timeseries, t if n == 0
  if n == 0
    return plott
  elseif n == 1 && !(typeof(sol[1]) <: Union{AbstractArray,ArrayPartition})
    return timeseries
  else
    tmp = Vector{eltype(sol[1])}(length(plot_timeseries))
    for j in 1:length(plot_timeseries)
      tmp[j] = plot_timeseries[j][n]
    end
    return tmp
  end
end

function solplot_vecs_and_labels(dims,vars,plot_timeseries,plott,sol,plot_analytic,plot_analytic_timeseries)
  plot_vecs = []
  for i in 2:dims
    push!(plot_vecs,[])
  end
  labels = String[]
  for x in vars
    tmp = []
    for j in 2:dims
      push!(tmp, u_n(plot_timeseries, x[j],sol,plott,plot_timeseries))
    end
    f = x[1]
    tmp = f.(tmp...)

    tmp = tuple((getindex.(tmp,i) for i in eachindex(tmp[1]))...)
    for i in eachindex(tmp)
      push!(plot_vecs[i],tmp[i])
    end
    add_labels!(labels,x,dims,sol)
  end

  if plot_analytic
    analytic_plot_vecs = []
    for i in 2:dims
      push!(analytic_plot_vecs,[])
    end
    for x in vars
      tmp = []
      for j in 2:dims
        push!(tmp, u_n(plot_analytic_timeseries, x[j],sol,plott,plot_analytic_timeseries))
      end
      f = x[1]
      tmp = f.(tmp...)
      tmp = tuple((getindex.(tmp,i) for i in eachindex(tmp[1]))...)
      for i in eachindex(tmp)
        push!(plot_vecs[i],tmp[i])
      end
      add_analytic_labels!(labels,x,dims,sol)
    end
  end
  plot_vecs = [hcat(x...) for x in plot_vecs]
  plot_vecs,labels
end

plot_indices(A::AbstractArray) = eachindex(A)
plot_indices(A::ArrayPartition) = eachindex(A)

ext: jl


text:
"""
A DiskProjMPO computes and stores the projection of an
MPO into a basis defined by an MPS, leaving a
certain number of site indices of the MPO unprojected.
Which sites are unprojected can be shifted by calling
the `position!` method.

Drawing of the network represented by a ProjMPO `P(H)`,
showing the case of `nsite(P)==2` and `position!(P,psi,4)`
for an MPS `psi`:

```
o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o-      -o--o--o--o--o--o |psi>
```

The environment tensors are stored on disk, which is helpful
for large bond dimensions if they cannot fit in memory.
"""
mutable struct DiskProjMPO <: AbstractProjMPO
  lpos::Int
  rpos::Int
  nsite::Int
  H::MPO
  LR::DiskVector{ITensor}
  Lcache::Union{ITensor,OneITensor}
  lposcache::Union{Int,Nothing}
  Rcache::Union{ITensor,OneITensor}
  rposcache::Union{Int,Nothing}
end

function DiskProjMPO(H::MPO)
  return new(
    0,
    length(H) + 1,
    2,
    H,
    disk(Vector{ITensor}(undef, length(H))),
    OneITensor,
    nothing,
    OneITensor,
    nothing,
  )
end

function disk(pm::ProjMPO)
  return DiskProjMPO(
    pm.lpos, pm.rpos, pm.nsite, pm.H, disk(pm.LR), lproj(pm), pm.lpos, rproj(pm), pm.rpos
  )
end
disk(pm::DiskProjMPO) = pm

# Special overload of lproj which uses the cached
# version of the left projected MPO, and if the
# cache doesn't exist it loads it from disk.
function lproj(P::DiskProjMPO)::Union{ITensor,OneITensor}
  (P.lpos <= 0) && return OneITensor()
  if (P.lpos ≠ P.lposcache) || (P.lpos == 1)
    # Need to update the cache
    P.Lcache = P.LR[P.lpos]
    P.lposcache = P.lpos
  end
  return P.Lcache
end

# Special overload of rproj which uses the cached
# version of the right projected MPO, and if the
# cache doesn't exist it loads it from disk.
function rproj(P::DiskProjMPO)::Union{ITensor,OneITensor}
  (P.rpos >= length(P) + 1) && return OneITensor()
  if (P.rpos ≠ P.rposcache) || (P.rpos == length(P))
    # Need to update the cache
    P.Rcache = P.LR[P.rpos]
    P.rposcache = P.rpos
  end
  return P.Rcache
end

function makeL!(P::DiskProjMPO, psi::MPS, k::Int)
  L = _makeL!(P, psi, k)
  if !isnothing(L)
    # Cache the result
    P.Lcache = L
    P.lposcache = P.lpos
  end
  return P
end

function makeR!(P::DiskProjMPO, psi::MPS, k::Int)
  R = _makeR!(P, psi, k)
  if !isnothing(R)
    # Cache the result
    P.Rcache = R
    P.rposcache = P.rpos
  end
  return P
end


text:
#=
DLC import Methods
=#

#=
Read data table and return julia arrays for x position, y position, and liklihood for each whisker point
=#
function dlc_hd5_to_array(path,l_thres,pole=true)

    #deeplabcut stores points as pandas dataframe
    #Assuming a whisker is labeled as a discrete set of points
    file=h5open(path)
    mytable=read(file,"df_with_missing")["table"];

    w_ind_length=length(mytable[1].data[2])
    if pole
        w_inds=collect(1:3:w_ind_length-3)
    else
        w_inds=collect(1:3:w_ind_length)
    end

    t1=1
    t2=length(mytable)

    #Each datapoint in an image of the pandas array is represented by values in 3 columns:
    #x,y and liklihood.
    #w_inds specifices the first index of each data point of interest
    xx=zeros(Float64,length(w_inds),length(mytable))
    yy=zeros(Float64,length(w_inds),length(mytable))
    ll=falses(length(w_inds),length(mytable))

    v_x=zeros(Float64,5)
    v_y=zeros(Float64,5)
    l=zeros(Float64,5)

    w_x=zeros(Float64,5)
    w_y=zeros(Float64,5)

    for i=t1:t2
        mydata=mytable[i].data[2]

        for j=1:length(w_inds)

            this_x=mydata[w_inds[j]]
            this_y=mydata[w_inds[j]+1]

            #We can set a liklihood threshold
            if mydata[w_inds[j]+2]>l_thres

                xx[j,i]=this_x
                yy[j,i]=this_y

                ll[j,i]=true
            end
        end
    end

    close(file)

    (xx,yy,ll)
end

#=
Updates liklihood matrix to false if points jump past
distance threshold from smoothed temporal traces
=#
function dlc_smooth_liklihood(xx::Array{Float64,2},yy::Array{Float64,2},kernel_size,ll,dist_thres)

    #outlier removal
    #smooth with gaussian kernel
    #calculate residual distance threshold and set Threshold
    x_smooth=zeros(Float64,size(xx))
    y_smooth=zeros(Float64,size(yy))

    for i=1:size(xx,1)
        x_smooth[i,:]=smooth(xx[i,:],kernel_size)
        y_smooth[i,:]=smooth(yy[i,:],kernel_size)
    end

    for i=1:size(xx,2)
        for j=1:size(xx,1)
            mydist=sqrt((xx[j,i]-x_smooth[j,i])^2+(yy[j,i]-y_smooth[j,i])^2)

            if mydist>dist_thres
                ll[j,i]=false
            end

        end
    end

    nothing
end

function convert_whisker_points_to_janelia(xx,yy,tracked)

    woi=[WhiskerTracking.Whisker1() for i=1:length(xx)]

    #Assuming is stored such that first data point centered on the follicle,
    #But Janelia assumes that the last point in the array is on the follicle, so we flip here
    #To take on the Janelia convention
    #Should probably be a flag for this
    for i=1:length(xx)
        num_points=length(xx[i])

        w=WhiskerTracking.Whisker1(0,i,num_points,xx[i],yy[i],ones(Float64,num_points),ones(Float64,num_points))
        woi[i]=w
    end

    woi
end

function convert_discrete_to_janelia(preds::Array,conf_thres::Float64,pad_pos::Tuple)
    woi=[WhiskerTracking.Whisker1() for i=1:size(preds,3)]
    tracked=trues(size(preds,3))

    myinds=falses(size(preds,1))
    for i=1:size(preds,3)
        num_points=0
        myinds[:]=falses(size(preds,1))
        for j=1:size(preds,1)
            if preds[j,3,i]>conf_thres
                myinds[j]=true
                num_points+=1
            end
        end

        woi[i]=WhiskerTracking.Whisker1(0,i,num_points,preds[myinds,1,i],preds[myinds,2,i],ones(Float64,num_points),ones(Float64,num_points))

        if num_points<5
            tracked[i]=false
        end

        #Points need to have a nearest neighbor at least 20 pixels away, otherwise they are just assumed
        #to be stray noise.
        for j=2:length(woi[i].x)
            if sqrt((woi[i].x[j]-woi[i].x[j-1])^2+(woi[i].y[j]-woi[i].y[j-1])^2) > 20.0
                tracked[i]=false
                break
            end
        end

        if tracked[i]
            WhiskerTracking.WT_reorder_whisker(woi[i],pad_pos)
        end
    end

    (woi,tracked)
end

function read_whisker_hdf5(path;l_thres=0.5,pole=true)

    (xx,yy,ll)=dlc_hd5_to_array(path,l_thres,pole)

    dlc_smooth_liklihood(xx,yy,15,ll,50.0)

    convert_whisker_points_to_janelia(xx,yy,ll)
end

function read_pole_hdf5(path)

    file=h5open(path)
    mytable=read(file,"df_with_missing")["table"];

    #Assumple pole is at the end of the array
    col_pos=length(mytable[1].data[2]) - 2

    p=zeros(Float64,length(mytable),2)

    for i=1:length(mytable)
        mydata=mytable[i].data[2]

        if mydata[col_pos+2]>0.1
            p[i,1]=mydata[col_pos]
            p[i,2]=mydata[col_pos+1]
        else
           p[i,1]=NaN
           p[i,2]=NaN
        end
    end

    close(file)

    p
end

function read_pole_and_whisker_hdf5(path,l_thres_in=0.5)

    woi=read_whisker_hdf5(path,l_thres=l_thres_in)
    p=read_pole_hdf5(path,16)

    (woi,p)
end

#=
Snapshot
=#
function save_single_image(han::Tracker_Handles,path=han.paths.images,use_contrast=false)
    save_single_image(han,han.current_frame,han.displayed_frame,path,use_contrast)
end

function save_single_image(han::Tracker_Handles,img::AbstractArray,num::Int,path=han.paths.images,use_contrast=false)

    img_name = name_img(han,"img",num)

    save_img_with_dir_change(han,img_name,img,path,use_contrast)

    nothing
end

function save_label_image(han::Tracker_Handles,path=han.paths.images,interp=false)

    img_name = name_img(han,"w",han.displayed_frame)

    img = create_label_image(han,0,interp)

    save_img_with_dir_change(han,img_name,img,path)

    nothing
end

function save_follicle_image(han::Tracker_Handles,path=han.paths.images)

    img_name = name_img(han,"w",han.displayed_frame)

    img = create_follicle_image(han)

    save_img_with_dir_change(han,img_name,img,path)

    nothing
end

function save_img_with_dir_change(han::Tracker_Handles,img_name::String,img::AbstractArray,path=han.paths.images,use_contrast=false)
    my_wd=pwd()
    cd(path)
    if use_contrast
        img2 = deepcopy(img)
        adjust_contrast(img2,han.contrast_min,han.contrast_max)
        Images.save(img_name,img2)
    else
        Images.save(img_name, img)
    end
    cd(my_wd)
end

function name_img(han::Tracker_Handles,name_prefix,num)

    #Deeplabcut needs left padded zeros
    max_dig=length(digits(han.max_frames))
    img_name = string(name_prefix,lpad(num,max_dig,"0"),".png")
end

function create_label_image(han::Tracker_Handles,rad=0,interp=false)

    img = zeros(UInt8,size(han.current_frame))

    w1=han.woi[han.displayed_frame]

    #Janelia whiskers are 0 in
    new_x = w1.x .+ 1
    new_y = w1.y .+ 1

    if interp
        (new_x, new_y) = interpolate_whisker(new_x,new_y,0.2)
    end

    for i=1:length(new_x)
        x=round(Int64,new_x[i])
        y=round(Int64,new_y[i])

        img[y,x] = 255 #Backbone

        for xx=-1*rad:rad
            for yy=-1*rad:rad
                if !(((y+yy)<1)|(y+yy>size(img,1))|((x+xx)<1)|(x+xx>size(img,2)))
                     img[y+yy,x+xx]=255
                end
            end
        end

    end
    img
end

#=
function create_follicle_image(han::Tracker_Handles,rad=5)

    img = zeros(UInt8,size(han.current_frame))

    x0 = han.woi[han.displayed_frame].x[1]]
    y0 = han.woi[han.displayed_frame].y[1]]

end
=#

function load_whisker_into_gui(han,whiskers::Array{Float64,2})

    han.tracked_whiskers = whiskers

    nothing
end


text:
# Functions to load complete test cases from files hosted in the repository

include("create_profile.jl")

"""
    load_cigre_mv_eu(<keyword arguments>)

Load an extended version of CIGRE MV benchmark network (EU configuration).

Source: <https://e-cigre.org/publication/ELT_273_8-benchmark-systems-for-network-integration-of-renewable-and-distributed-energy-resources>, chapter 6.2

Extensions:
- storage at bus 14 (in addition to storage at buses 5 and 10, already present);
- candidate storage at buses 5, 10, and 14;
- time series (8760 hours) for loads and RES generators.

# Arguments
- `flex_load::Bool = false`: toggles flexibility of loads.
- `ne_storage::Bool = false`: toggles candidate storage.
- `scale_gen::Float64 = 1.0`: scale factor of all generators, wind included.
- `scale_wind::Float64 = 1.0`: further scaling factor of wind generator.
- `scale_load::Float64 = 1.0`: scale factor of loads.
- `number_of_hours::Int = 8760`: number of hourly optimization periods.
- `start_period::Int = 1`: first period of time series to use.
- `year_scale_factor::Int = 10`: how many years a representative year should represent [years].
- `energy_cost::Float64 = 50.0`: cost of energy exchanged with transmission network [€/MWh].
- `cost_scale_factor::Float64 = 1.0`: scale factor for all costs.
"""
function load_cigre_mv_eu(;
        flex_load::Bool = false,
        ne_storage::Bool = false,
        scale_gen::Float64 = 1.0,
        scale_wind::Float64 = 1.0,
        scale_load::Float64 = 1.0,
        number_of_hours::Int = 8760,
        start_period::Int = 1,
        year_scale_factor::Int = 10, # years
        energy_cost::Float64 = 50.0, # €/MWh
        cost_scale_factor::Float64 = 1.0,
    )

    grid_file = normpath(@__DIR__,"..","data","cigre_mv_eu","cigre_mv_eu_more_storage.m")
    sn_data = _FP.parse_file(grid_file)
    _FP.add_dimension!(sn_data, :hour, number_of_hours)
    _FP.add_dimension!(sn_data, :scenario, Dict(1 => Dict{String,Any}("probability"=>1)), metadata = Dict{String,Any}("mc"=>true))
    _FP.add_dimension!(sn_data, :year, 1; metadata = Dict{String,Any}("scale_factor"=>year_scale_factor))
    _FP.add_dimension!(sn_data, :sub_nw, 1)

    # Set cost of energy exchanged with transmission network
    sn_data["gen"]["14"]["ncost"] = 2
    sn_data["gen"]["14"]["cost"] = [energy_cost, 0.0]

    # Scale wind generation
    sn_data["gen"]["6"]["pmin"] *= scale_wind
    sn_data["gen"]["6"]["pmax"] *= scale_wind
    sn_data["gen"]["6"]["qmin"] *= scale_wind
    sn_data["gen"]["6"]["qmax"] *= scale_wind

    # Toggle flexible demand
    for load in values(sn_data["load"])
        load["flex"] = flex_load ? 1 : 0
    end

    # Toggle candidate storage
    if !ne_storage
        sn_data["ne_storage"] = Dict{String,Any}()
    end

    _FP.scale_data!(sn_data; cost_scale_factor)
    _FP.add_td_coupling_data!(sn_data; sub_nw = 1)
    d_time_series = create_profile_data_cigre(sn_data, number_of_hours; start_period, scale_load, scale_gen, file_profiles_pu=normpath(@__DIR__,"..","data","cigre_mv_eu","time_series","CIGRE_profiles_per_unit_Italy.csv"))
    d_mn_data = _FP.make_multinetwork(sn_data, d_time_series)

    return d_mn_data
end


text:
import Distributions.sample

#=
    randmvn(S)

Sampling singular semidefinite multivariate normal distributions ``x ~ N(0, A)``

```
S = cholfact(A, :L, Val{true})
x = randmvn(S)
```

=#
randmvn(S) = ipermute!(S[:L]*permute!(randn(size(S,1)),S[:p]), S[:p])

function randmvn(S::Base.LinAlg.CholeskyPivoted, n)
    d = size(S,1)
    Z = zeros(d, n)
    for i in 1:n
        Z[:,i] = randmvn(S)
    end
    Z
end

"""
```
Y, X = sample(n, m, M::LinearHomogSystem)
```

Sample observations generated by `M`.

- n -- number of observations per simulation
- m -- number of independent simulations
- M -- linear dynamic system and observation model

- Y -- observations (d2xnxm array)
- X -- simulated processes (dxnxm array)

"""
function sample(n, m, M::LinearHomogSystem{Vector{T}}) where {T}

    H, Phi, b = M.H, M.Phi, M.b

    d2, d = size(M.H)

    X = zeros(d, n, m)
    Y = zeros(d2, n, m)

    Nv = MvNormal(zeros(d2), M.R)
    S = cholfact(M.Q, :L, Val{true}) # allow singular Q
    Sx0 = cholfact(M.P0, :L, Val{true}) #
    #Nw = MvNormal(zeros(d), M.Q)
    #Nx0 = MvNormal(M.x0, M.P0)

    for j in 1:m
        x = M.x0 + randmvn(Sx0)
        V = rand(Nv, n)
        W = randmvn(S, n)
        X[.., 1, j] = Phi*x + b + W[:, 1]
        Y[.., 1, j] = H*x + V[:, 1]
        for i in 2:n
            x[:] = Phi*x + b + W[:, i]
            X[.., i, j] = x
            Y[.., i, j] = H*x + V[:, i]
        end

    end

    return Y, X
end

function sample(n, m, M::LinearHomogSystem{Tx, TP, Ty}) where {Tx, TP, Ty}

    H, Phi, b = M.H, M.Phi, M.b

    X = zeros(Tx, n, m)
    Y = zeros(Ty, n, m)

    Nx0 = Gaussian(M.x0, M.P0)
    Nw = Gaussian(zero(Tx), M.Q)
    Nv = Gaussian(zero(Ty), M.R)

    for j in 1:m
        x = rand(Nx0)
        V = rand(Nv, (n,))
        W = rand(Nv, (n,))
        X[1, j] = Phi*x + b + W[1]
        Y[1, j] = H*x + V[1]
        for i in 2:n
            x = Phi*x + b + W[i]
            X[i, j] = x
            Y[i, j] = H*x + V[i]
        end

    end

    return Y, X
end

function sample(n, M)
    Y, X = sample(n, 1, M)
    squeeze(Y, ndims(Y)), squeeze(X, ndims(X))
end


text:
module Place

	export buildmodel, placebo, predict, freerun, description, vembed, BasisFunc, PlaceModel

	using ToeplitzMatrices, Random, StatsBase, Distributions, LinearAlgebra

	include("BasisFunctionTypes.jl")

	struct BasisFunc
    	funct
    	embed
    	radius
    	centre
    	params
	end

	struct PlaceModel
		rbf
		λ
		params
		basis
		embedv
	end

	function predict(model,y)
		#use the model to predict y
		#
		#embed
		X,yt =embedding(y,model.embedv)
		#compute ϕ
		ϕ, offset = placebo(X, model.rbf)
		#compute prediction
		yp=ϕ[:,model.basis]*model.λ
		ep=yp-yt
		#
		return yp,yt,ep
	end

	function freerun(model,y,npts ::Int64=0,i ::Int64=1)
		#use the model to predict y
		#
		#embed
		X,yt =embedding(y,model.embedv)
		dX,nx =size(X)
		if npts<1
			npts=length(yt)
		end
		Xi=X[:,i]
		yt=yt[i:end]
		Xp=Array{Float64,2}(undef,dX,npts)
		for j in 1:npts
			ϕ, offset = placebo(Xi[:,:], model.rbf) #Xi needs to be 2-D!
			yp=ϕ[:,model.basis]*model.λ
			Xi=[yp; Xi[1:end-1]]     #future is at the TOP!!!!
			Xp[:,j]=Xi
		end
		yp=Xp[1,:]
		if length(yt)>npts
			yt=yt[1:npts]
		end
		#
		return yp, yt
	end

	function buildmodel(y,options=[])
		#build a time series model of y splitting at point td between build and test

		#sort out all the options
		if "nbasis"∈keys(options)
			nbasis=options["nbasis"]
		else
			nbasis=500 #default 500
		end
		if "testdatum"∈keys(options)
			td=options["testdatum"]
		else
			td=[]
		end
		if "stopstep"∈keys(options)
			stopstep=options["stopstep"]
		else
			stopstep=6 #default 6
		end
		if "embedding"∈keys(options)
			v=options["embedding"]
		else
			v=[0, 1, 2, 3]
		end
		if v isa Tuple
			nv=length(v)
		else
			nv=1
		end
		if "penalty"∈keys(options)
			penalty=options["penalty"]
		else
			penalty = :(nx*log(mss)+2*nk)  	#Akaike Information Criterion
		#	penalty = :(nx*log(mss)+nk*log(nx))	#Schwarz Information Crtierion
			penalty = :description(mss,λ,δ,nx) #Rissanen desciption length
		end
		if "functions"∈keys(options)
			functype=options["functions"]
		else
			functype=(gaussian, tophat)
		end
		if "nneighbours"∈keys(options)
			nneighbours=options["nneighbours"]
		else
			nneighbours=1
		end

		#check for NaNs and any other potential trouble makers
		if any(x->isnan(x),y)
			println("NaN detected in input to Place.buildmodel")
		end

		#embed the data
		X,z,Xp,zp = embedding(y,v,td)

		#generate some basis functions
		rbfset = getbasis(X,z,functype,v, nbasis, nneighbours)

		#choose basis functions
		rbf,λ,basis,mdlv = topdown(X,z,rbfset,penalty,stopstep)

		#build and output model structure
		model = PlaceModel(rbf,λ,params,basis,v)
		#
		#done
		return model, X, z, mdlv
	end

	function embedding(y,v,td=0)
		#embed y with strategy adequate for v
		#remember - future is at the TOP!
		de=maximum([maximum(vv) for vv in v])+1 #maximum embedding lag required
		x0=y[1:(de+1)] #initial condition
		z=y[(de+1):end] #predicted target
		X=ToeplitzMatrices.Hankel(x0,z) #sparse embedded matrix
		X=X[end-1:-1:1,:] #state matrix flip it
		(de,nX)=size(X)
		if isempty(td)
			td=Int(floor(nX/2))
		end
		if td>nX
			td=Int(floor(nX*0.9))
		end
		if td==0
			Xp=[]
			zp=[]
			td=nX+1
		else
			Xp=X[:, td:end]
			zp=z[td:end] #Xp and zp are the test data
		end
		X=X[:, 1:(td-1)]
		z=z[1:(td-1)]#X and z are used to build the model
		#
		return X, z, Xp, zp
	end

	function vembed(v::Vector{T}) where T
		#compute all possible variable embeddings from v
		#i.e. all possible nonempty subsets of v as a Tuple
    	result = Vector{T}[[]]
    	for elem in v, j in eachindex(result)
        	push!(result, [result[j] ; elem])
    	end
	    return tuple(result[2:end]...)
	end

	function nearestofnneighbours(X::Array{Float64,2},randi::Array{Int64,1},y,n::Int64=1)
		#distances to the closest of n randomly choosen neighbours - multiple reference points
		dx,nx =size(X)
		ny=length(y)
		nbasis=length(randi)
		dd=Array{Float64,1}(undef,nbasis)
		dd .= Inf
		for i in 1:n
			randomneighbour=sample(1:ny, Weights(abs.(y)),nbasis)
			randr=X[:,randi]-X[:,randomneighbour]
			randr=normn(randr)
			dd=minimum([dd'; randr], dims=1)'
		end
		return dd
	end

	function nearestofnneighbours(X::Array{Float64,2},randi::Int64,y,n::Int64=1)
		#distances to the closest of n randomly choosen neighbours - single reference point
		dx = length(X)
		ny = length(y)
		randomneighbour=sample(1:ny, Weights(abs.(y)),n)
		randr=X[:,randi]*ones(1,n) - X[:,randomneighbour]
		randr=normn(randr)
		dd = minimum(randr)
		return dd
	end

	function normn(X,n:: Int64=2)
		#n-norm over columns of a square array X
		#
		return (sum(X.^n, dims=1)).^(1/n)
	end

	function getbasis(X :: Array{Float64,2}, y :: Array{Float64, 1}, functions, v, nbasis, nnearest) :: Array{BasisFunc,1}
		#generate a random set of basis functions, distributed over the data X according to pdf (X,y)
		de,nx=size(X)
		rbfs=Array{BasisFunc,1}(undef,nbasis)
		randfunct=Array{Any,1}(undef,nbasis)
		if functions isa Tuple
			randfunct=rand(functions,nbasis)
		else
			randfunct.=functions
		end
		randi=sample(1:nx, Weights(abs.(y)),nbasis)
	#	randomneighbour=sample(1:nx, Weights(abs.(y)),nbasis)
		randr=Array{Float64,1}(undef,nbasis)
		if v isa Tuple
			randv=rand(v,nbasis)
			for i in 1:nbasis
				#randr[i]=normn(X[randv[i].+1,randi[i]]-X[randv[i].+1,randomneighbour[i]])[1]
				randr[i]=nearestofnneighbours(X[randv[i].+1,:],randi[i],y,nnearest)
			end
		else
			randv=[]
			#randr=X[v.+1,randi]-X[v.+1,randomneighbour]
			#randr=normn(randr)
			randr=nearestofnneighbours(X[v.+1,:],randi,y,nnearest)
		end
		randr=abs.(randr.*randn(nbasis))
		for i in 1:nbasis
			#
			thisfunct=randfunct[i]
			theseparams=[]
			tfp=thisfunct()
			if ~isempty(tfp)
				np,dp=size(tfp)
				theseparams=Array{Float64,1}(undef,np)
				for npi in 1:np #choose option parameters, if any
					theseparams[npi] = rand(Uniform(tfp[npi,1],tfp[npi,2]))
				end
			end
			if ~isempty(randv)
				thisv=randv[i]
			else
				thisv=v
			end
			rbf=BasisFunc(thisfunct, thisv, randr[i], X[thisv.+1,randi], theseparams)
			rbfs[i]=rbf
			#
		end
		#
		#
		return rbfs
	end

	function finddeltas(Q,δ)
		#compute solution of Qδ=1/δ from initial guess δ
		#via Newtons method
	#	print(size(Q))
	#	println(" in finddeltas")
		nq,mq=size(Q) # should be the same and the same as length(δ)
		if nq==1 #solveable for dimension 1
			δ= 1 ./ sqrt(abs.(Q))
		elseif nq==2 #still solevable but messy
			δ=find2deltas(Q,δ)
		else #otherwise...
			#Newton method
			#
			#set up parameters
			max_dd= 1e-4
			max_err=1e-4
			min_delta=1e-6
			max_count=1e2
			now_count =0
			# Initialise
			δ = abs.(δ) #should be not necessary
			err = Q*δ - 1 ./ δ
			dd = 2*max_dd
			bestδ = δ
			while any(abs.(dd).>max_dd) || any(err.>max_err) || now_count>max_count
				dQ=(Q + diagm(1 ./ δ.^2))
  				if false #CHECK rcond(SS)<eps;
  					#deal with ill -conditioning
 				else
     				dd = -(Q + diagm(1 ./ δ.^2))\err
  				end
  				δ = δ + dd
  				#push negative \deltas back
  				δ = map(x->max(x,min_delta),δ)
  				err = Q*δ - 1 ./ δ
  				now_count += 1
			end
		end
	#	println(size(Q))
		return δ
	end

	function find2deltas(Q::Array{Float64,2},δ)
		#compute solution of Qδ=1/δ for the case of 2-by-2 Q
		nq,mq=size(Q) # should be the same and the same as length(δ)
		if mq!=2 || nq!=2
			println("WARNING Q not 2-by-2 in find2deltas")
		end
		q11=Q[1,1] #just easier to do this way
		q12=Q[1,2]
		q21=Q[2,1]
		q22=Q[2,2]
		a = q22*q11^2/q12^2-q21*q11/q12
		b = q21/q12-2*q11*q22^2/q12-1
		c = q22/q12^2
		Δ = (b^2-4*a*c)
		if Δ<0
			println("WARNING Qδ=1/δ has no exact solution in 2-D case - giving up")
			return δ
		end
		δ1 = [ (-b+sqrt(Δ))/(2*a), (-b-sqrt(Δ))/(2*a) ]
		δ1 = filter(x-> x>0 ,δ1) #positive roots only, because we next...
		δ1 = sqrt(δ1)
		if isempty(δ1)
			println("WARNING: no positive δ, giving up and guessing")
			δ1 = -[ (-b+sqrt(Δ))/(2*a), (-b-sqrt(Δ))/(2*a) ]
			δ1 = min(δ1)
		end
		δ2 = (1 ./ δ1 .- q11*δ1)/q12
		δ=[δ1, δ2]
		return δ
	end

	function lstar(λ,δ) # λ::Float64,δ::Float64 ?
		#Compute universal prior encoding of positive float λ to precision δ
        λ=abs.(λ)
		δ=abs.(δ)
	#	println("In lstar: λ=$λ, δ=$δ")
		bigl1 = lstar(2*ceil.(λ ./ δ))
		bigl2 = lstar(ceil.(abs.(2*log.(map(x->max(x,1/x),λ)))))
	#	println("Output: $bigl1+$bigl2=$(bigl1+bigl2)")
		return bigl1+bigl2
	end

	function lstar(λ) # λ::Int64 ?
		#Compute universal prior encoding of positive integers λ
        λ=abs.(λ)
		bigl=log(2.865)
		λ=1 .+ ceil.(λ)
		while any(x-> x>1, λ)
			λ=log.(map(x->max(x,1),λ))
			bigl=bigl .+ λ
		end
		return bigl
	end

	function description(mss,λ,δ,nx)
		#Compute description length encoding (ala Rissanen) of model with
		#error mss, parameters λ to precision δ and nx data
		global needδ=true
		return 0.5*nx*(1+log(2*π*mss)) + sum(lstar(λ,δ))
	end

	function topdown(X,y,allrbfs:: Array{BasisFunc,1},penalty::Expr,stopstep::Int64)
		#select an optimal set of basis functions to fit X to y.
		global needδ, mss, nx, nk, λ, δ   #needed for penalty function evaluation
		ϕ,offset=placebo(X,allrbfs)
		nb=length(allrbfs)
		modelbasis=[]
		ϕQ=Array{Float64,2}(undef,0,0)
		ϕR=Array{Float64,2}(undef,0,0)
		λ=Array{Float64,1}(undef,0)
		δ=Array{Float64,1}(undef,0)
		err=y
		mss=mean(y'*y)
		ny=length(y)
		dx,nx=size(X)
		dϕ,nϕ=size(ϕ)
		println("dx=$dx, nx=$nx, nϕ=$nϕ, offset=$offset")
		nk=0
		mdl= Inf
		mdlv=[]
		notimproved=0
		bestmodelbasis=modelbasis
		needδ=false
		eval(penalty) #just to set needδ=true if necessary...
		println("needδ=$needδ")
		while nk<nϕ && notimproved<stopstep
			#
			μ = -ϕ'*err #compute sensitivity
			#first try an expand
			~,ind = findmax(abs.(μ))
			append!(modelbasis, ind)
			ϕQ,ϕR = qrappend(ϕQ,ϕR,ϕ[:,ind])
			λ = ϕR\(ϕQ'*y) #new model weights
			if needδ
				if isempty(δ)
					δ=sqrt(abs(mean(λ)))
				else
					δ=δ[:] #needs a vector not an 1-by-n matrix - damn you Julia
					append!(δ, mean(δ))
				end
			end

			#now delete
			~,ind = findmin(abs.(λ))
			if ind==(nk+1) #last added is worst so keep it
				nk += 1
			else
				#last added is not worst, so delete new worst
				deleteat!(modelbasis,ind)
				ϕQ,ϕR = qrdelete(ϕQ,ϕR,ind)
				if needδ
					δ=δ[:] #needs a vector not an 1-by-n matrix - damn you Julia
					deleteat!(δ,ind)
				end
			end
			#
			if needδ
				Q = ϕQ'*ϕQ/mss
				δ = finddeltas(Q,δ)
			end
			#
			λ = ϕR\(ϕQ'*y) #new model weights
			err = y - ϕ[:,modelbasis]*λ
         	mss = (err'*err/ny)
         	dl=eval(penalty)
         	if dl<mdl #is this model best?
         		mdl=dl
         		bestmodelbasis=deepcopy(modelbasis)
         		notimproved = 0
         		print("*")
         	else
         		notimproved += 1
         	end
         	println("MSS=$mss DL=$dl size=$nk")
         	if nk>length(mdlv)
         		append!(mdlv,dl)
         	else
         		mdlv[nk]=dl
         	end
		end
		#λ = ϕR\(ϕQ'*y)
		modelbasis=bestmodelbasis
		rbfs = allrbfs[filter(x -> x>offset, modelbasis).-offset]
		mbi = [filter(x->x<=offset, modelbasis); (1:count(modelbasis.>offset)) .+ offset]
		ϕ, = placebo(X,rbfs)
		λ = ϕ[:,mbi]\y
		#
		return rbfs, λ, mbi, mdlv
 	end

	function placebo(X :: Array{Float64,2}, rbfs ::Array{BasisFunc,1}, constant::Bool=true, linear::Bool=true)
		#evaluate the basis functions rbfs at the points of X
		#include constant and linear terms in the model candidates
		(de,nx)=size(X)
		nb=length(rbfs)
		offset=Int64(constant)+Int64(linear)*de
		ϕ=Array{Float64,2}(undef,nx,nb+offset)
		if constant
			ϕ[:,1] .= 1
		end
		if linear
			ϕ[:,2:(de+1)]=X'
		end
		for (i,rbf) in enumerate(rbfs)
			#
			ϕX = X[rbf.embed.+1, : ] - rbf.centre[rbf.embed.+1]*ones(1,nx)
			ϕX = normn(ϕX) ./ rbf.radius
			if isempty(rbf.params)
				ϕi = rbf.funct(ϕX)
			else
				ϕi = rbf.funct(ϕX, rbf.params)
			end
			ϕ[:,i+offset] = ϕi
			#
		end
		#
		if any(x->isnan(x),ϕ) #something is screwed up to cause this to arise in the first place...
			println("NaN generated in placebo")
			ϕ[isnan.(ϕ)] .= 0
		end
		#
		return ϕ, offset
	end

	function qrappend(Q,R,x)
		# [Q,R] = qrappend(Q,R,x) append a column to a QR factorization.
		# If [Q,R] = qr(A) is the original QR factorization of A,
		# then [Q,R] = qrappend(Q,R,x) changes Q and R to be the
		# factorization of the matrix obtained by appending an extra
		# column, x, to A.
		# ported from MATLAB code
		#
		if isempty(Q)
  			fact = qr(x)
  			Q=Matrix(fact.Q)
  			R=fact.R
		else
			~,m = size(Q)
			m += 1
			r = Q'*x			 # best fit of x by Q
			R = [R r] 		 	 # add coeff to R
		  	q= x - Q*r 		 	 # q is orthogonal part of x
			f= norm(q)
		  	R = [R; zeros(1,m);] # update R for q
  			R[m,m] = f[1] 		 # f is coeff of q when normalized
			Q = [Q q/f] 		 # extend basis by normalized q
		end
		#
		#
		return Q,R
	end

	function qrdelete(Q,R,j)
		# qrdelete(Q,R,j) delete the j-th column of A=QR from the QR factorisation
		# That is, Q1*R1 is the matrix Q*R with it's j-th column removed
		# adapted from MATLAB implementation of the same thing
		#
		mq,nq = size(Q)
		R=R[:, 1:end .!=j]
		m,n = size(R)
		#
		for k in j:min(n,m-1)
			p=k:(k+1)
			G,R[p,k] = planerot(R[p,k])
			if k<n
				R[p,k+1:n] = G*R[p,k+1:n]
			end
			Q[:,p]=Q[:,p]*G'
		end
		Q=Q[:,1:n]      #additional clean-up
		R=R[1:n,1:n]	#bog-standard qrdelete misses this
		#
		#
		return Q,R
	end

	function planerot(x)
		#planerot(x) should do what Givens does - Givens rotation
		#this code, I did just plunder from github
	    if length(x) != 2
	        error()
	    end
	    da, db = x
	    roe = db
	    if abs(da) > abs(db)
	       roe = da
	    end
	    scale = abs(da) + abs(db)
	    if scale == 0
	        c = 1.
	        s = 0.
	        r = 0.
	    else
	        r = scale * sqrt((da/scale)^2 + (db/scale)^2)
	        r = sign(roe) * r
	        c = da/r
	        s = db/r
	    end
	    #
	    #
	    return [c s; -s c], [r, 0.]
	end

end # of module

max_stars_repo_name: m-small/place

max_issues_repo_name: m-small/place

max_forks_repo_name: m-small/place


text:
"""
# Usage

    solve(prob::FODESystem, h, NonLinearAlg())

Nonlinear algorithm for nonlinear fractional differential equations.

### References

Dingyu Xue, Northeastern University, China ISBN:9787030543981
"""
struct NonLinearAlg <: FractionalDiffEqAlgorithm end

function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
    @unpack f, α, u0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    n = length(u0)
    m = round(Int, (T-t0)/h)+1
    g = genfun(1)
    g = g[:]
    u0 = u0[:]
    ha = h.^α
    z = zeros(n, m)
    x1 = copy(u0) # Here we pass the value of x0 to x1. Honestly, I kept finding this bug for almost a whole night😅

    # All of the min(m, L0+1) is to set the memory effect.
    SetMemoryEffect = Int64(min(m, L0+1))
    W = zeros(n, SetMemoryEffect) #Initializing W a n*m matrix

    @fastmath @inbounds @simd for i = 1:n
        W[i, :] = getvec(α[i], SetMemoryEffect, g)
    end

    du = zeros(n)
    @fastmath @inbounds @simd for k = 2:m
        tk = t0+(k-1)*h
        L = min(Int64(k-1), Int64(L0))
        f(du, x1, nothing, tk)

        @fastmath @inbounds @simd for i = 1:n
            x1[i] = du[i]*ha[i] - W[i, 2:L+1]'*z[i, k-1:-1:k-L] + u0[i]
        end
        z[:, k] = x1 - u0
    end

    result = (z + repeat(u0, 1, m))'

    return result
end

"""
P-th precision polynomial generate function

```math
g_p(z)=\\sum_{k=1}^p \\frac{1}{k}(1-z)^k
```
"""
function genfun(p)
    a = collect(1:p+1)
    A = Vandermonde(a)'
    return (1 .-a')*inv(A')
end

function getvec(α, n, g)
    p = length(g)-1
    b = 1 + α
    g0 = g[1]
    w = Float64[]
    push!(w, g[1]^α)

    @fastmath @inbounds @simd for m = 2:p
        M = m-1
        dA = b/M
        temp = (-(g[2:m] .*collect((1-dA):-dA:(1-b))))' *w[M:-1:1]/g0
        push!(w, temp)
    end

    @fastmath @inbounds @simd for k = p+1:n
        M = k-1
        dA = b/M
        temp = (-(g[2:(p+1)] .*collect((1-dA):-dA:(1-p*dA))))' *w[M:-1:(k-p)]/g0
        push!(w, temp)
    end
    return w
end


text:
using HTbase
using Test

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


text:
const synset_test_line = string(
    "13134947 20 n 01 fruit 0 037 @ 11675842 n 0000 + 02397119 a 0101 + 101",
    "13997 n 0101 + 13135692 n 0101 + 01652731 v 0101 + 01652895 v 0101 + 0",
    "0506672 v 0101 + 00056188 v 0101 ~ 07705931 n 0000 ~ 11636835 n 0000 ~",
    " 11700279 n 0000 ~ 12036067 n 0000 ~ 12158031 n 0000 ~ 12162758 n 0000",
    " ~ 12193334 n 0000 ~ 12267677 n 0000 ~ 12301445 n 0000 ~ 12620546 n 00",
    "00 ~ 12642090 n 0000 ~ 12644283 n 0000 ~ 12647787 n 0000 ~ 12650805 n ",
    "0000 ~ 12658481 n 0000 ~ 12737251 n 0000 ~ 13133613 n 0000 ~ 13135692 ",
    "n 0000 ~ 13135832 n 0000 ~ 13137409 n 0000 ~ 13137672 n 0000 ~ 1313830",
    "8 n 0000 ~ 13138842 n 0000 ~ 13139055 n 0000 ~ 13139482 n 0000 ~ 13140",
    "367 n 0000 ~ 13141415 n 0000 ~ 13150378 n 0000 ~ 13150592 n 0000 | the",
    " ripened reproductive body of a seed plant"
)

const synset_test_line_2 = string(
    "00074790 04 n 0b blunder 0 blooper 0 bloomer 0 bungle 0 pratfall 0 fou",
    "l-up 0 fuckup 0 flub 0 botch 0 boner 0 boo-boo 0 019 @ 00070965 n 0000",
    " + 02229000 a 0901 + 02527651 v 0901 + 02527651 v 0808 + 02527651 v 07",
    "18 + 02527651 v 0616 + 02527651 v 040d + 00013172 v 0401 + 02566227 v ",
    "0103 ~ 00071864 n 0000 ~ 00075283 n 0000 ~ 00075471 n 0000 ~ 00075790 ",
    "n 0000 ~ 00075912 n 0000 ~ 00076072 n 0000 ~ 00076196 n 0000 ~ 0007632",
    "3 n 0000 ~ 00076393 n 0000 ~ 00076563 n 0000 | an embarrassing mistake"
)

const synset_test_line_3 = string(
    "10951697 18 n 03 Edward a Edward_VIII 0 Duke_of_Windsor 0 002 @i 10233",
    "445 n 0000 #m 08159740 n 0000 | King of England and Ireland in 1936; h",
    "is marriage to Wallis Warfield Simpson created a constitutional crisis"
)

facts("Synset") do
    ss = Synset(synset_test_line, 'n')

    context("constructor") do
        @fact ss.synset_type --> 'n'
        @fact ss.gloss --> "the ripened reproductive body of a seed plant"
        @fact collect(words(ss)) --> ["fruit"]
    end

    context("parses n_words from a hex string") do
        hex_ss = Synset(synset_test_line_2, 'n')
        @fact hex_ss.gloss --> "an embarrassing mistake"
    end

    context("parses the lex_id as a hex string") do
        hex_ss = Synset(synset_test_line_3, 'n')
    end

    context("words") do
        @fact collect(words(ss)) --> ["fruit"]
    end

    context("word_count") do
        @fact word_count(ss) --> 1
    end

    context("show") do
        expected = "(n) fruit (the ripened reproductive body of a seed plant)"
        @fact io_to_string(io -> show(io, ss)) --> expected
    end

    context("is_nothing") do
        @fact WordNet.is_nothing(WordNet.∅) --> true
    end
end

facts("helpers") do
    context("parse_int_hex") do
        @fact WordNet.parse_int_hex("a") --> 10
    end
end


text:
"""
    AbstractClassFunction{T}
Abstract type representing functions constant on conjugacy classes of a group.

The following functionality is required for an `AbstractClassFunction`:
 * `parent(χ)`: the underlying group
 * `conjugacy_classes(χ)`: iterator over conjugacy classes of `χ`.
 * `values(χ)`: a _typed_ iterator over values
 * `getindex(χ, i::Integer)`: the value on `i`-th conjugacy class.
 Note: Indexing with negative integers should return values on the class which
 contains inverses of the `i`-th class.
 * `χ(g::GroupElement)`: value of `χ` on a group element (only for julia < 1.3)
"""
abstract type AbstractClassFunction{T} end # <:AbstractVector{T}?

Base.eltype(::AbstractClassFunction{T}) where T = T

function LinearAlgebra.dot(χ::AbstractClassFunction, ψ::AbstractClassFunction)
    @assert conjugacy_classes(χ) === conjugacy_classes(ψ)
    val = sum(
        length(cc) * χ[i] * ψ[-i] for (i, cc) in enumerate(conjugacy_classes(χ))
    )
    orderG = sum(length, conjugacy_classes(χ))
    val = _div(val, orderG)
    return val
end

function (χ::AbstractClassFunction)(g::GroupElement)
    for (i, cc) in enumerate(conjugacy_classes(χ))
        g ∈ cc && return χ[i]
    end
    throw(
        DomainError(
            g,
            "element does not belong to conjugacy classes of $χ",
        ),
    )
end

_div(val, orderG) = div(val, orderG)
_div(val::AbstractFloat, orderG) = val / orderG
_div(val::Complex{<:AbstractFloat}, orderG) = val / orderG

## Arbitrary ClassFunctions without decomposition into irreps
struct ClassFunction{T,CCl} <: AbstractClassFunction{T}
    vals::Vector{T}
    conjugacy_classes::CCl
    inv_of::Vector{Int}
end

ClassFunction(vals, cclasses) = ClassFunction(vals, cclasses, _inv_of(cclasses))

## AbstractClassFunction api
Base.parent(χ::ClassFunction) = parent(first(first(conjugacy_classes(χ))))
conjugacy_classes(χ::ClassFunction) = χ.conjugacy_classes
Base.values(χ::ClassFunction) = χ.vals

Base.@propagate_inbounds function Base.getindex(χ::ClassFunction, i::Integer)
    i = i < 0 ? χ.inv_of[-i] : i
    @boundscheck 1 ≤ i ≤ length(conjugacy_classes(χ))
    return values(χ)[i]
end

# TODO: move to AbstractClassFunction when we drop julia-1.0
# adding call methods to abstract types was not supported before julia-1.3
function (χ::ClassFunction)(g::GroupElement)
    for (i, cc) in enumerate(conjugacy_classes(χ))
        g ∈ cc && return χ[i]
    end
    throw(
        DomainError(
            g,
            "element does not belong to conjugacy classes of $χ",
        ),
    )
end

"""
    Character <: AbstractClassFunction
Struct representing (possibly virtual) character of a group.

Characters are backed by `table(χ)::CharacterTable` which actually stores the
character values. The constituents (decomposition into the irreducible summands)
of a given character can be obtained by calling `constituents(χ)` which returns a
vector of coefficients of `χ` in the basis of `irreducible_characters(table(χ))`.

It is assumed that equal class functions on the same group will have **identical**
(ie. `===`) character tables.
"""
struct Character{T, S, ChT<:CharacterTable} <: AbstractClassFunction{T}
    table::ChT
    constituents::Vector{S}
end

function Character{R}(
    chtbl::CharacterTable{Gr, T},
    constituents::AbstractVector{S}
) where {R, Gr, T, S}
    return Character{R, S, typeof(chtbl)}(chtbl, constituents)
end

function Character(chtbl::CharacterTable, constituents::AbstractVector)
    R = Base._return_type(*, Tuple{eltype(chtbl), eltype(constituents)})
    @assert R ≠ Any
    return Character{R}(chtbl, constituents)
end

Character(chtbl::CharacterTable, i::Integer) = Character{eltype(chtbl)}(chtbl, i)

function Character{T}(chtbl::CharacterTable, i::Integer) where T
    v = zeros(Int, nconjugacy_classes(chtbl))
    v[i] = 1
    return Character{T, Int, typeof(chtbl)}(chtbl, v)
end

function Character{T}(χ::Character) where T
    S = eltype(constituents(χ))
    ChT = typeof(table(χ))
    return Character{T, S, ChT}(table(χ), constituents(χ))
end

## Accessors
table(χ::Character) = χ.table
constituents(χ::Character) = χ.constituents

## AbstractClassFunction api
Base.parent(χ::Character) = parent(table(χ))
conjugacy_classes(χ::Character) = conjugacy_classes(table(χ))
Base.values(χ::Character{T}) where T = T[χ[i] for i in 1:nconjugacy_classes(table(χ))]

Base.@propagate_inbounds function Base.getindex(
    χ::Character{T},
    i::Integer,
) where {T}
    i = i < 0 ? table(χ).inv_of[abs(i)] : i
    @boundscheck 1 ≤ i ≤ nconjugacy_classes(table(χ))

    if isirreducible(χ)
        k = findfirst(!iszero, constituents(χ))::Int
        return convert(T, table(χ)[k, i])
    else
        return convert(
            T,
            sum(
                c * table(χ)[idx, i] for
                (idx, c) in enumerate(constituents(χ)) if !iszero(c)
            ),
        )
    end
end

## Basic functionality

Base.:(==)(χ::Character, ψ::Character) =
    table(χ) === table(ψ) && constituents(χ) == constituents(ψ)
Base.hash(χ::Character, h::UInt) = hash(table(χ), hash(constituents(χ), h))

Base.deepcopy_internal(χ::Character, ::IdDict) =
    Character(table(χ), copy(constituents(χ)))

## Character arithmetic

for f in (:+, :-)
    @eval begin
        function Base.$f(χ::Character, ψ::Character)
            @assert table(χ) === table(ψ)
            return Character(table(χ), $f(constituents(χ), constituents(ψ)))
        end
    end
end

Base.:*(χ::Character, c::Number) = Character(table(χ), c .* constituents(χ))
Base.:*(c::Number, χ::Character) = χ * c
Base.:/(χ::Character, c::Number) = Character(table(χ), constituents(χ) ./ c)

## Group-theoretic functions:

PermutationGroups.degree(χ::Character) = Int(χ(one(parent(χ))))
PermutationGroups.degree(
    χ::Character{T,CCl},
) where {T,CCl <: AbstractOrbit{<:AbstractMatrix}} = Int(χ[1])

function Base.conj(χ::Character{T, S}) where {T, S}
    vals = collect(values(χ))
    all(isreal, vals) && return Character{T}(χ)
    tbl = table(χ)
    ψ = ClassFunction(vals[tbl.inv_of], conjugacy_classes(tbl), tbl.inv_of)
    constituents = S[dot(ψ, χ) for χ in irreducible_characters(tbl)]
    return Character{T, eltype(constituents), typeof(tbl)}(tbl, constituents)
end

isvirtual(χ::Character) =
    any(<(0), constituents(χ)) || any(!isinteger, constituents(χ))

function isirreducible(χ::Character)
    C = constituents(χ)
    k = findfirst(!iszero, C)
    k !== nothing || return false # χ is zero
    isone(C[k]) || return false # muliplicity is ≠ 1
    kn = findnext(!iszero, C, k+1)
    kn === nothing && return true # there is only one ≠ 0 entry
    return false
end

"""
    affordable_real!(χ::Character)
Return either `χ` or `2re(χ)` depending whether `χ` is afforded by a real
representation, modifying `χ` in place.
"""
function affordable_real!(χ::Character)
    ι = frobenius_schur(χ)
    if ι <= 0 # i.e. χ is complex or quaternionic
        χ.constituents .+= constituents(conj(χ))
    end
    return χ
end

"""
    frobenius_schur(χ::AbstractClassFunction[, pmap::PowerMap])
Return Frobenius-Schur indicator of `χ`, i.e. `Σχ(g²)` where sum is taken over
the whole group.

If χ is an irreducible `Character`, Frobenius-Schur indicator takes values in
`{1, 0, -1}` which correspond to the following situations:
 1. `χ` is real-valued and is afforded by an irreducible real representation,
 2. `χ` is a complex character which is not afforded by a real representation, and
 3. `χ` is quaternionic character, i.e. it is real valued, but is not afforded by a
 real representation.

In cases 2. and 3. `2re(χ) = χ + conj(χ)` corresponds to an irreducible character
afforded by a real representation.
"""
function frobenius_schur(χ::Character)
    @assert isirreducible(χ)

    pmap = powermap(table(χ))
    ι = sum(
        length(c) * χ[pmap[i, 2]] for (i, c) in enumerate(conjugacy_classes(χ))
    )

    ι_int = Int(ι)
    ordG = sum(length, conjugacy_classes(χ))
    d, r = divrem(ι_int, ordG)
    @assert r == 0 "Non integral Frobenius Schur Indicator: $(ι_int) = $d * $ordG + $r"
    return d
end

Base.isreal(χ::Character) = frobenius_schur(χ) > 0

function Base.show(io::IO, ::MIME"text/plain", χ::Character)
    println(io, "Character over ", eltype(χ))
    _print_char(io, χ)
end

Base.show(io::IO, χ::Character) = _print_char(io, χ)

function _print_char(io::IO, χ::Character)
    first = true
    for (i, c) in enumerate(constituents(χ))
        iszero(c) && continue
        first || print(io, " ")
        print(io, ((c < 0 || first) ? "" : '+'))
        !isone(c) && print(io, c, '·')
        print(io, 'χ', FiniteFields.subscriptify(i))
        first = false
    end
end


text:
	function get_keys(dict)
		return sort([i for i=keys(dict)])
	end


text:
using Markov
using Catlab
using Catlab.Theories
using Test

X = Ob(FreeCartesianCategory, :Float64)
u = Hom(x->x*rand(), X, X)

u₀ = Hom(()->rand(), munit(FreeCartesianCategory.Ob), X)
crand(u₀⊗u₀)

meantest(f::Function, μ::Real, n::Int, ϵ::Real) = begin
    μ̂ = sum(map(f, 1:n))/n
    @test μ - ϵ < μ̂
    @test μ̂ < μ + ϵ
end

meantest(1/1, 10000, 1e-1) do i
    crand(((u⊗u)⋅braid(X,X)), 1,2)[1]
end
meantest(1/2, 10000, 1e-1) do i
    crand(((u⊗u)⋅braid(X,X)), 1,2)[2]
end
meantest(1/1, 10000, 1e1) do i
    crand(mcopy(X),2)[1]
end
meantest(1/1, 10000, 1e1) do i
    crand(mcopy(X),2)[2]
end

@test crand(delete(X), 1) == ()
@test crand(u⋅delete(X), 1) == ()
meantest(1/2, 1000, 1e-1) do i
    crand(u⋅mcopy(X)⋅(id(X)⊗delete(X)), 1)[1][1]
end

@testset "Uniforms" begin
  u₀ = Hom(()->rand(), munit(FreeCartesianCategory.Ob), X)
  @test 0 <= crand(u, 1.0) <= 1.0
  @test 0 <= crand(u, 2) <= 2
  @test compose(u,u).args |> length == 2
  @test dom(compose(u,u)) == X
  @test codom(compose(u,u)) == X
  @test 0 <= crand(compose(u,u), 0.5) <= 1/2
  @test 0 <= crand(u₀) <= 1
  @test crand(u₀⊗u₀) |> length == 2
  @test crand(u⊗u₀, 1.0) |> length == 2
  @test crand(u₀⊗u, 1.0) |> length == 2
  @test crand(otimes(u₀,u₀,u), 2.0) |> length == 3
  @test crand(otimes(u₀,u,u₀), 2.0) |> length == 3
  @test crand(otimes(u₀,u,u), 1.0, 2.0) |> length == 3
  @test crand(otimes(u₀,u,u)⋅otimes(u,u,u), 1.0, 2.0) |> length == 3
  @test crand(otimes(u,u,u), 1.0, 2.0, 3.0) |> length == 3
end

@test crand(compose(id(X), id(X), id(X)), 1)[1] == 1


text:
using Test
using StatsModelComparisons, StanSample
using StatsFuns, CSV, Random

@testset "Cars"  begin

    ProjDir = @__DIR__

    df = CSV.read(joinpath(ProjDir, "..", "data", "cars.csv"), DataFrame)

    cars_stan = "
    data {
        int N;
        vector[N] speed;
        vector[N] dist;
    }
    parameters {
        real a;
        real b;
        real sigma;
    }
    transformed parameters{
        vector[N] mu;
        mu = a + b * speed;
    }
    model {
        a ~ normal(0, 100);
        b ~ normal(0, 10);
        sigma ~ exponential(1);
        dist ~ normal(mu, sigma)    ;
    }
    generated quantities {
        vector[N] log_lik;
        for (i in 1:N)
            log_lik[i] = normal_lpdf(dist[i] | mu[i], sigma);
    }
    "

    Random.seed!(1)
    cars_stan_model = SampleModel("cars.model", cars_stan)
    data = (N = size(df, 1), speed = df.Speed, dist = df.Dist)
    rc = stan_sample(cars_stan_model; data)

    if success(rc)
        nt_cars = read_samples(cars_stan_model);
    end

    log_lik = nt_cars.log_lik'
    n_sam, n_obs = size(log_lik)
    lppd = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs)

    pwaic = [var(log_lik[:, i]) for i in 1:n_obs]
    @test -2(sum(lppd) - sum(pwaic)) ≈ 421.0 atol=2.0
    loo, loos, pk = psisloo(log_lik)
    @test -2loo ≈ 421.0 atol=2.0

end


text:
module json

using Base: peek # not exported until Julia 1.5

function parse(str::AbstractString)
    try
        tokens = Iterators.Stateful(tokenize(str))
        value = parse_value!(tokens)
        isempty(tokens) || error()
        return value
    catch
        error("invalid JSON: $str")
    end
end

function parse_value!(tokens::Iterators.Stateful)
    token = peek(tokens)
    token == '[' ? parse_array!(tokens) :
    token == '{' ? parse_dict!(tokens) :
    token isa Char ? error() : popfirst!(tokens)
end

function pop_comma!(tokens::Iterators.Stateful)
    peek(tokens) == ',' || return false
    popfirst!(tokens)
    return true
end

function parse_array!(tokens::Iterators.Stateful)
    array = []
    popfirst!(tokens) == '[' || error()
    while peek(tokens) != ']'
        value = parse_value!(tokens)
        push!(array, value)
        pop_comma!(tokens) || break
    end
    popfirst!(tokens) == ']' || error()
    return array
end

function parse_dict!(tokens::Iterators.Stateful)
    dict = Dict{String,Any}()
    popfirst!(tokens) == '{' || error()
    while peek(tokens) != '}'
        key = popfirst!(tokens)
        key isa String || error()
        popfirst!(tokens) == ':' || error()
        value = parse_value!(tokens)
        dict[key] = value
        pop_comma!(tokens) || break
    end
    popfirst!(tokens) == '}' || error()
    return dict
end

function tokenize(str::AbstractString)
    tokens = []
    chars = Iterators.Stateful(str)
    isspecial(c::Char) = isspace(c) || c in "{}[],:\""
    for c in chars
        if isspace(c)
            continue
        elseif c == '"'
            s = sprint() do io
                for c in chars
                    if c == '"'
                        break
                    elseif c == '\\'
                        c = popfirst!(chars)
                        if c == 'u'
                            s = popfirst!(chars) * popfirst!(chars) *
                                popfirst!(chars) * popfirst!(chars)
                            print(io, Char(parse(Int, s, base=16)))
                        else
                            c = c == 'b' ? '\b' :
                                c == 'f' ? '\f' :
                                c == 'n' ? '\n' :
                                c == 'r' ? '\r' :
                                c == 't' ? '\t' :
                                c == 'v' ? '\v' :
                                c in "\\\"'" ? c : error()
                            print(io, c)
                        end
                    else
                        print(io, c)
                    end
                end
            end
            push!(tokens, s)
        elseif isspecial(c)
            push!(tokens, c)
        else
            word = sprint() do io
                print(io, c)
                while !isempty(chars) && !isspecial(peek(chars))
                    print(io, popfirst!(chars))
                end
            end
            token = word == "true"  ? true :
                    word == "false" ? false :
                    word == "null"  ? nothing :
                    word[1] in "0123456789.+-" ? parse(Float64, word) : error()
            push!(tokens, token)
        end
    end
    return tokens
end

end # module


text:
fm1 = lmm(Yield ~ 1 + (1|Batch), ds)
fm1w = lmm(Yield ~ 1 + (1|Batch), ds; wt = ones(size(ds, 1)))

@test size(fm1.A) == (3, 3)
@test size(fm1.wttrms) == (3,)
@test size(fm1.R) == (3, 3)
@test size(fm1.Λ) == (1,)
@test lowerbd(fm1) == zeros(1)
@test fm1[:θ] == ones(1)

fm1[:θ] = [0.713]
@test objective(fm1) ≈ 327.34216280955366
MixedModels.describeblocks(IOBuffer(), fm1)

fit!(fm1);
@test_approx_eq_eps objective(fm1) 327.3270598811428 1e-3
@test_approx_eq_eps fm1[:θ] [0.752580] 1.e-5
@test_approx_eq_eps deviance(fm1) 327.32705988 1.e-3
@test_approx_eq_eps AIC(fm1) 333.3270598811394 1.e-3
@test_approx_eq_eps BIC(fm1) 337.5306520261259 1.e-3
@test fixef(fm1) ≈ [1527.5]
@test MixedModels.fixef!(zeros(1),fm1) ≈ [1527.5]
@test coef(fm1) ≈ [1527.5]
@test cond(fm1) == ones(1)
cm = coeftable(fm1)
@test size(cm.mat) == (1,3)
@test MixedModels.fnames(fm1) == [:Batch]
@test model_response(fm1) == convert(Vector,ds[:Yield])
@test abs(sum(ranef(fm1,true)[1])) < 1.e-5

@test_approx_eq_eps logdet(fm1) 8.06014522999825 1.e-3
@test_approx_eq_eps varest(fm1) 2451.2501089607676 1.e-3
@test_approx_eq_eps pwrss(fm1) 73537.50326882303 1.e-1
@test_approx_eq_eps stderr(fm1) [17.69455188898009] 1.e-4

vc = VarCorr(fm1)
show(IOBuffer(), vc)
@test vc.s == sdest(fm1)
srand(1234321)
simulate!(fm1)   # changes the fit
@test_approx_eq_eps deviance(fm1) 339.0218639362958 1.e-3
bootstrap(fm1,1,(i,x) -> fixef(x)) # bootstrap restores the fit
@test_approx_eq_eps deviance(fm1) 339.0218639362958 1.e-3

fm2 = lmm(Yield ~ 1 + (1|Batch),ds2)
@test lowerbd(fm2) == zeros(1)
#redirect_stdout(IOBuffer())
fit!(fm2,true)
#redirect_stdout()
show(IOBuffer(), fm2)
@test abs(fm2[:θ][1]) < 1.0e-9
@test objective(fm2) ≈ 162.87303665382575
@test abs(std(fm2)[1][1]) < 1.0e-9
@test std(fm2)[2] ≈ [3.653231351374652]
@test stderr(fm2) ≈ [0.6669857396443261]
@test coef(fm2) ≈ [5.6656]
@test logdet(fm2) ≈ 0.0
refit!(fm2,convert(Vector,ds[:Yield]))

fm3 = lmm(Reaction ~ 1+Days + (1+Days|Subject),slp)
@test lowerbd(fm3) == [0.,-Inf,0.]
@test isa(fm3.A[1,1],MixedModels.HBlkDiag{Float64})
@test size(fm3.A[1,1]) == (36,36)
@test fm3.A[1,1][1,1] == 10.
@test fm3.A[1,1][6,1] == 0.
@test_throws BoundsError size(fm3.A[1,1],0)
@test size(fm3.A[1,1],1) == 36
@test full(fm3.A[1,1])[1:2,1:2] == reshape([10.,45,45,285],(2,2))

fit!(fm3)

@test objective(fm3) ≈ 1751.9393444647046
@test_approx_eq_eps fm3[:θ] [0.9292213074888169,0.01816838485113137,0.22264487095998978] 1.e-6
@test pwrss(fm3) ≈ 117889.46144025437
@test_approx_eq_eps logdet(fm3) 73.90322021999222 1e-3
@test_approx_eq_eps stderr(fm3) [6.632257721914501,1.5022354739749826]  1.e-4
@test coef(fm3) ≈ [251.40510484848477,10.4672859595959]
@test fixef(fm3) ≈ [251.40510484848477,10.4672859595959]
@test_approx_eq_eps stderr(fm3) [6.632246393963571,1.502190605041084] 1.e-2
@test_approx_eq_eps std(fm3)[1] [23.780468100188497,5.716827903196682] 1.e-2
@test_approx_eq_eps logdet(fm3) 73.90337187545992 1.e-3
@test_approx_eq diag(cor(fm3)[1]) ones(2)
@test_approx_eq_eps cond(fm3) [4.1752507630514915] 1.e-4
@test loglikelihood(fm3) ≈ -875.9696722323523
@test eltype(fm3.wttrms[1]) === Float64
const u3 = ranef(fm3,true)
@test length(u3) == 1
@test size(u3[1]) == (2,18)
@test_approx_eq_eps u3[1][1,1] 3.030300122575336 1.e-3
const b3 = ranef(fm3)
@test length(b3) == 1
@test size(b3[1]) == (2,18)
@test_approx_eq_eps b3[1][1,1] 2.815819441982976 1.e-3

fm4 = lmm(Reaction ~ Days + (1|Subject) + (0+Days|Subject), slp);
@test size(fm4) == (180,2,36,2)
@test fm4[:θ] == ones(2)
@test lowerbd(fm4) == zeros(2)

fit!(fm4);

@test_approx_eq_eps deviance(fm4) 1752.0032551398835 1.e-3
@test_approx_eq_eps objective(fm4) 1752.0032551398835 1.e-3
@test_approx_eq coef(fm4) [251.40510484848585,10.467285959595715]
@test_approx_eq fixef(fm4) [251.40510484848585,10.467285959595715]
@test_approx_eq_eps stderr(fm4) [6.707710260366577,1.5193083237479683] 1.e-3
@test_approx_eq_eps fm4[:θ] [0.9458106880922268,0.22692826607677266] 1.e-4
@test std(fm4)[1] ≈ [24.171449463289047]
@test std(fm4)[2] ≈ [5.799379721123582]
@test std(fm4)[3] ≈ [25.556130034081047]
@test_approx_eq_eps logdet(fm4) 74.46952585564611 1.e-3

tbl = MixedModels.lrt(fm4,fm3)

@test_approx_eq_eps tbl[:Deviance] [1752.0032551398835,1751.9393444636157] 1e-3
@test tbl[:Df] == [5,6]

simulate!(fm3)  # to test one of the unscaledre methods

fm5 = lmm(Diameter ~ (1|Plate) + (1|Sample), pen);
@test size(fm5) == (144,1,30,2)
@test fm5[:θ] == ones(2)
@test lowerbd(fm5) == zeros(2)

fit!(fm5)

@test_approx_eq_eps objective(fm5) 332.18834867227616 1.e-3
@test_approx_eq_eps coef(fm5) [22.97222222222222] 1.e-3
@test_approx_eq_eps fixef(fm5) [22.97222222222222] 1.e-3
@test coef(fm5)[1] ≈ mean(pen[:Diameter])
@test_approx_eq_eps stderr(fm5) [0.7445960346851368] 1.e-4
@test_approx_eq_eps fm5[:θ] [1.5375772376554968,3.219751321180035] 1.e-3
@test_approx_eq_eps std(fm5)[1] [0.8455645948223015] 1.e-4
@test_approx_eq_eps std(fm5)[2] [1.770647779277388] 1.e-4
@test_approx_eq_eps varest(fm5) 0.3024263987592062 1.e-4
@test_approx_eq_eps logdet(fm5) 95.74614821367786 1.e-3
@test length(ranef(fm5)) == 2

fm6 = lmm(Strength ~ (1|Sample) + (1|Batch), psts);
@test size(fm6) == (60,1,40,2)
@test fm6[:θ] == ones(2)
@test lowerbd(fm6) == zeros(2)

fit!(fm6);

@test_approx_eq_eps objective(fm6) 247.99446586289676 1.e-3
@test_approx_eq_eps coef(fm6) [60.05333333333329] 1.e-3
@test_approx_eq_eps fixef(fm6) [60.05333333333329] 1.e-3
@test_approx_eq_eps stderr(fm6) [0.6421359883527029] 1.e-4
@test_approx_eq_eps fm6[:θ] [3.5268858714382905,1.3299230213750168] 1.e-3
@test_approx_eq_eps std(fm6)[1] [2.904069002535747] 1.e-3
@test_approx_eq_eps std(fm6)[2] [1.095070371687089] 1.e-4
@test_approx_eq_eps std(fm6)[3] [0.8234088395243269] 1.e-4
@test_approx_eq_eps varest(fm6) 0.6780020742644107 1.e-4
@test_approx_eq_eps logdet(fm6) 101.0381339953986 1.e-3

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using BenchmarkTools
using StaticArrays

using ComputedGeodesicEquations

randpos() = (rand(0.0:1000.0), rand(1.1:10.0), rand(0.1:π/1.1), rand(-π/3:π/3))

us = [ SVector{4, Float64}([randpos()...]) for i in 1:1000 ]
vs = [ SVector{4, Float64}([randpos()...]) for i in 1:1000 ]

function runner(us, vs, m)
    for (u, v) in zip(us, vs)
        geodesic_eq(m, u, v)
    end
end

suite = BenchmarkGroup()
suite["regular"] = BenchmarkGroup()
suite["jacobian"] = BenchmarkGroup()

suite["regular"]["boyer-lindquist"] = @benchmarkable runner($us, $vs, m) setup=(m = BoyerLindquist())
suite["jacobian"]["boyer-lindquist"] = @benchmarkable runner($us, $vs, m) setup=(m = BoyerLindquistJac())

suite["regular"]["johannsen-psaltis"] = @benchmarkable runner($us, $vs, m) setup=(m = JohannsenPsaltis())
suite["jacobian"]["johannsen-psaltis"] = @benchmarkable runner($us, $vs, m) setup=(m = JohannsenPsaltisJac())

suite["regular"]["eddington-finkelstein"] = @benchmarkable runner($us, $vs, m) setup=(m = EddingtonFinkelstein())
suite["jacobian"]["eddington-finkelstein"] = @benchmarkable runner($us, $vs, m) setup=(m = EddingtonFinkelsteinJac())

suite["regular"]["morris-thorne"] = @benchmarkable runner($us, $vs, m) setup=(m = MorrisThorne())
suite["jacobian"]["morris-thorne"] = @benchmarkable runner($us, $vs, m) setup=(m = MorrisThorneJac())

tune!(suite)
results = run(suite, verbose = true)


text:
module NoiseFourier

include("FourierBasis.jl")
include("Fourier1DMatrix.jl")

using .Fourier1D
export assemble_matrix, noise_matrix, truncation_error, rigorous_norm, assemble_matrix_high_precision
export evaluate_trig_polynomial, abscissas

end


text:
function syrkd!(C::StridedMatrix{T}, A::SparseMatrixCSC{T,BlasInt}, op::Char, α::Number, β::Number) where {T}
    return syrkd!(C, csrptr(A'), op, α, β)
end