text:
module Extrema_Test

using Compat
using Compat.Test
using EAGOSmoothMcCormickGrad
using IntervalArithmetic
using StaticArrays

function about(calc,val,tol)
    return (val - tol <= calc <= val + tol)
end

EAGOSmoothMcCormickGrad.set_diff_relax(0)

a = seed_g(Float64,1,2)
xIBox = [Interval(-3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = min(3,X)
println("out 1: $out")
@test about(out.cc,3.0,1E-1)
@test about(out.cv,1.0909090909090908,1E-1)
@test about(out.cc_grad[1],0.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],0.545455,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,-3,1E-4)
@test about(out.Intv.hi,8,1E-4)

a = seed_g(Float64,1,2)
xIBox = [Interval(-3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = max(5,X)
println("out 2: $out")
@test about(out.cc,7.045454545454545,1E-1)
@test about(out.cv,5.0,1E-1)
@test about(out.cc_grad[1],0.272727,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],0.0,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,-3,1E-4)
@test about(out.Intv.hi,8,1E-4)

EAGOSmoothMcCormickGrad.set_diff_relax(1)

a = seed_g(Float64,1,2)
xIBox = [Interval(-3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = min(3,X)
println("out 3: $out")
@test about(out.cc,3.0,1E-1)
@test about(out.cv,1.0909090909090908,1E-1)
@test about(out.cc_grad[1],0.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],0.545455,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,-3,1E-4)
@test about(out.Intv.hi,8,1E-4)

a = seed_g(Float64,1,2)
xIBox = [Interval(-3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = max(5,X)
println("out 4: $out")
@test about(out.cc,7.045454545454545,1E-1)
@test about(out.cv,5.0,1E-1)
@test about(out.cc_grad[1],0.272727,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],0.0,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,-3,1E-4)
@test about(out.Intv.hi,8,1E-4)

end


text:
#
# driverdriver.jl::
#
# trivial coupled system mock up
#
# very simple top level driver that show proof of concept implementation for
# simple couplng API sufficient for v0 atmos-land-ocean DG tools.
# Components
# _A - atmos
# _L - land
# _O - ocean
# _C - coupler
#
# driverdriver multiplexes land ocean atmos and coupler components.
#
# all exist on aligned horizontal grids decomposed on ranks in aligned way.
#
# components in this example run sequentially on shared set/subset of
# underlying resources without spinwait contention.
#
# components have their own MPI communicators to support MPI rank disaggregation onto
# distinct resources or to interleave on same resources with some lightweight
# resource contention coordination.
#
# disaggregated ranks require coupler to handle asynchronous send/receive posting
#
#
abstract type CLIMAcomponent end
abstract type  DGModelComponent <: CLIMAcomponent end
struct CouplerComponent         <: CLIMAcomponent end
struct AtmosModelComponent      <: DGModelComponent end
struct LandModelComponent       <: DGModelComponent end
struct OceanModelComponent      <: DGModelComponent end

comp_config(args...; kargs...) = nothing;
#
using MPI
using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes

 #####
 # Basic initialization
 #####
 ClimateMachine.init()
 ArrayType = ClimateMachine.array_type()
 mpicomm = MPI.COMM_WORLD
 FT = Float64

 # In general we might dup communicator and/or create intercommunicators/
 # intracommicators depending on resource use and parallelism mode. Here we
 # assume sequential, all componenta on all same ranks, aligned. So each component
 # just uses default mpicomm - (also not sure if CLiMA as written always uses
 # component communicator). We carry different names so we can make them different
 # in the future.
 mpicomm_A=mpicomm
 mpicomm_O=mpicomm
 mpicomm_L=mpicomm
 mpicomm_C=mpicomm

 #####
 # Create mesh
 #####

 # Horiz mesh
 # Same horizontal mesh size, extents and polynomical order for all components
 # We still create distinct grids so that they can use different communicators when
 # needed.
 const N = 4
 const Nˣ = 10
 const Nʸ = 10
 const Lˣ = 4e6  # m
 const Lʸ = 4e6  # m
 xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
 yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
 brickrange_2D = (xrange, yrange)
 topol_2D( comm ) = (
   BrickTopology( comm, brickrange_2D, periodicity = (false, false) ) );
 grid_2D( topol ) = ( DiscontinuousSpectralElementGrid(
   topol,
   FloatType = FT,
   DeviceArray = ArrayType,
   polynomialorder = N,
 )
 )
 topo_2D_A  = topol_2D( mpicomm_A )
 grid_2D_A  = grid_2D(  topo_2D_A )
 topo_2D_O  = topol_2D( mpicomm_O )
 grid_2D_O  = grid_2D(  topo_2D_O )
 topo_2D_L  = topol_2D( mpicomm_L )
 grid_2D_L  = grid_2D(  topo_2D_L )
 topo_2D_C  = topol_2D( mpicomm_C )
 grid_2D_C  = grid_2D(  topo_2D_C )

 # Vert mesh (stacked topology)
 # Different components can have different vertical extents and element counts
 const Nᶻ = 20
 # Atmos vert
 const Nᶻ_A = 20
 # Ocean vert
 const Nᶻ_O = 50
 # Land vert - lets makes this 2d as an example
 const Nᶻ_L =  0
 # Coupler vert - lets makes this 2d too for now
 const Nᶻ_C =  0

 const H_A = 100000  # m
 const H_O = 5000    # m
 zrange_A  = range(FT(0   ); length = Nᶻ_A + 1, stop = H_A)
 zrange_O  = range(FT(-H_O); length = Nᶻ_O + 1, stop = 0  )
 brickrange_3D_A = (xrange, yrange, zrange_A)
 brickrange_3D_O = (xrange, yrange, zrange_O)
 topol_3D( comm, br ) = (
  StackedBrickTopology(
    mpicomm,
    br;
    periodicity = (false, false, false),
    boundary = ((1, 1), (1, 1), (2, 3)),
  )
 );
 topol_3D_A = topol_3D( mpicomm_A, brickrange_3D_A  )
 topol_3D_O = topol_3D( mpicomm_A, brickrange_3D_O  )
 grid_3D( topol ) = (
  DiscontinuousSpectralElementGrid(
   topol,
   FloatType = FT,
   DeviceArray = ArrayType,
   polynomialorder = N,
  )
 );
 grid_3D_A=grid_3D( topol_3D_A )
 grid_3D_O=grid_3D( topol_3D_O )

 # Config and init coupler
 comp_C=comp_config( CouplerComponent(), ( (), ) )

 # Config (including mask) and init drivers of individual components
 comp_A=comp_config( AtmosModelComponent() ; grids=( (grid=grid_3D_A, label="g3da",), ) )
 comp_L=comp_config( LandModelComponent()  ; grid2D=grid_2D_L )
 comp_O=comp_config( OceanModelComponent() ; grid3D=grid_3D_O )

ext: jl

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
import Printf

export heritability,
  logpdf,
  gradient!, gradient,
  fisher_Σ!, fisher_Σ, fisher_B!, fisher_B,
  mle_fs!, mle_mm!,
  fit_mle!, fit_reml!,
  update_meanparam!

#---------------------------------------------------------------------------#
# Evaluate log-pdf
#---------------------------------------------------------------------------#

"""
    logpdf(vcmrot, vcobsrot)

Calculate log-pdf of a [`TwoVarCompVariateRotate`](@ref) instance under a
[`TwoVarCompModelRotate`](@ref).
"""
function logpdf(
  vcmrot::T1,
  vcobsrot::T2,
  vcaux::T3 = VarianceComponentAuxData(vcobsrot)
  ) where {
      T1 <: TwoVarCompModelRotate,
      T2 <: TwoVarCompVariateRotate,
      T3 <: VarianceComponentAuxData}

  n, d = size(vcobsrot.Yrot, 1), size(vcobsrot.Yrot, 2)
  T = eltype(vcmrot)
  zeroT, oneT = zero(T), one(T)
  residual!(vcaux.res, vcmrot, vcobsrot)
  # evaluate 2(log-likehood)
  objval = convert(T, - n * d * log(2π) - d * vcobsrot.logdetV2 - n * vcmrot.logdetΣ2)
  tmp, λj = zeroT, zeroT
  @inbounds for j in 1:d
    λj = vcmrot.eigval[j]
    @simd for i in 1:n
      tmp = oneT / (vcobsrot.eigval[i] * λj + oneT)
      objval += log(tmp) - tmp * vcaux.res[i, j]^2
    end
  end
  objval /= 2
end

"""
    logpdf(vcm, vcobsrot)

Calculate log-pdf of a [`TwoVarCompVariateRotate`](@ref) instance under a
[`VarianceComponentModel`](@ref).
"""
function logpdf(
  vcm::T1,
  vcobs::T2,
  vcaux::T3 = VarianceComponentAuxData(vcobs)
  ) where {
    T1 <: VarianceComponentModel,
    T2 <: TwoVarCompVariateRotate,
    T3 <: VarianceComponentAuxData}

  logpdf(TwoVarCompModelRotate(vcm), vcobs, vcaux)
end

function logpdf(
  vcm::T1,
  vcobs::T2,
  vcaux::T3 = VarianceComponentAuxData(vcobs)
  ) where {
    T1 <: VarianceComponentModel,
    T2 <: VarianceComponentVariate,
    T3 <: VarianceComponentAuxData}

  logpdf(TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobs), vcaux)
end

"""
Calculate log-pdf of an array of variance component instances.
"""
function logpdf(
  vcm::T1,
  vcdata::Array{T2},
  vcaux::Array{T3} = VarianceComponentAuxData(vcdata)
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate},
    T3 <: VarianceComponentAuxData}

  mapreduce(x -> logpdf(vcm, x...), +, zip(vcdata, vcaux))
end

#---------------------------------------------------------------------------#
# Evaluate gradient
#---------------------------------------------------------------------------#

"""
    gradient!(∇, vcmrot, vcobsrot)

Evaluate gradient of `Σ` at `vcmrot.Σ` and overwrite `∇`.

# Input
- `∇::Vector`: gradient vector.
- `vcmrot`: *rotated* two variance component model [`TwoVarCompModelRotate`](@ref).
- `vcobsrot`: *rotated* two variance component data [`TwoVarCompVariateRotate`](@ref).

# Output
- `∇`: gradient vector at `Σ = (Σ[1], Σ[2])`.

# TODO
- optimize computation
"""
function gradient!(
  ∇::AbstractVector{T},
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobsrot)
  ) where {T <: AbstractFloat}

  n, d = size(vcobsrot.Yrot, 1), size(vcobsrot.Yrot, 2)
  zeroT, oneT = zero(T), one(T)
  residual!(vcaux.res, vcmrot, vcobsrot)
  # gradient wrt Σ[2]
  tmp, λj, dg = zeroT, zeroT, zeros(T, d)
  @inbounds for j in 1:d
    λj = vcmrot.eigval[j]
    @simd for i in 1:n
      tmp = oneT / (vcobsrot.eigval[i] * λj + oneT)
      vcaux.res[i, j] *= tmp
      dg[j] += tmp
    end
  end
  N = transpose(vcaux.res) * vcaux.res
  @inbounds for j in 1:d
    N[j, j] -= dg[j]
  end
  #A_mul_Bt!(N, N, vcmrot.eigvec), A_mul_B!(N, vcmrot.eigvec, N)
  N = vcmrot.eigvec * N * vcmrot.eigvec'
  copyto!(∇, d^2 + 1, N, 1, d^2)
  # gradient wrt Σ[1]
  @inbounds for j in 1:d
    λj = vcmrot.eigval[j]
    dg[j] = zeroT
    @simd for i in 1:n
      tmp = vcobsrot.eigval[i]
      vcaux.res[i, j] *= √tmp
      dg[j] += tmp / (tmp * λj + oneT)
    end
  end
  mul!(N, transpose(vcaux.res), vcaux.res)
  @inbounds for j in 1:d
    N[j, j] -= dg[j]
  end
  #A_mul_Bt!(N, N, vcmrot.eigvec), A_mul_B!(N, vcmrot.eigvec, N)
  N = vcmrot.eigvec * N * vcmrot.eigvec'
  copyto!(∇, N)
  # scale by 0.5
  rmul!(∇, 1//2)
end # function gradient!

"""

    gradient(vcmrot, vcobsrot)

Evaluate gradient of `Σ` at `vcmrot.Σ` and overwrite `∇`.

# Input
- `vcmrot`: *rotated* two variance component model [`TwoVarCompModelRotate`](@ref).
- `vcobsrot`: *rotated* two variance component data [`TwoVarCompVariateRotate`](@ref).

# Output
- `∇`: gradient vector at `Σ = (Σ[1], Σ[2])`.

# TODO
- optimize computation
"""
function gradient(
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobsrot)
  ) where {T <: AbstractFloat}

  d = length(vcmrot)
  ∇ = zeros(T, 2d^2)
  gradient!(∇, vcmrot, vcobsrot, vcaux)
end

"""
    gradient!(∇, vcm, vcobsrot)

Evaluate gradient of `Σ` at `vcmrot.Σ` and overwrite `∇`.

# Input
- `∇::Vector`: gradient vector.
- `vcmrot`: two variance component model [`VarianceComponentModel`](@ref).
- `vcobsrot`: *rotated* two variance component data [`TwoVarCompVariateRotate`](@ref).

# Output
- `∇`: gradient vector at `Σ = (Σ[1], Σ[2])`.

# TODO
- optimize computation
"""
function gradient!(
  ∇::AbstractVector{T},
  vcm::VarianceComponentModel{T, 2},
  vcobs::TwoVarCompVariateRotate{T},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobs)
  ) where {T <: AbstractFloat}

  gradient!(∇, TwoVarCompModelRotate(vcm), vcobs, vcaux)
end

function gradient!(
  ∇::AbstractVector{T},
  vcm::VarianceComponentModel{T, 2},
  vcobs::VarianceComponentVariate{T, 2},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobs)
  ) where {T <: AbstractFloat}

  gradient!(∇, TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobs), vcaux)
end

"""
    gradient(vcm, vcobsrot)

Evaluate gradient of `Σ` at `vcmrot.Σ` and overwrite `∇`.

# Input
- `vcmrot`: two variance component model [`VarianceComponentModel`](@ref).
- `vcobsrot`: *rotated* two variance component data [`TwoVarCompVariateRotate`](@ref).

# Output
- `∇`: gradient vector at `Σ = (Σ[1], Σ[2])`.

# TODO
- optimize computation
"""
function gradient(
  vcm::VarianceComponentModel{T, 2},
  vcobs::TwoVarCompVariateRotate,
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobs)
  ) where {T <: AbstractFloat}

  gradient(TwoVarCompModelRotate(vcm), vcobs, vcaux)
end

function gradient(
  vcm::VarianceComponentModel{T, 2},
  vcobs::VarianceComponentVariate{T, 2},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobs)
  ) where {T <: AbstractFloat}

  gradient(TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobs), vcaux)
end

function gradient!(
  ∇::AbstractVector,
  vcm::T1,
  vcdata::Array{T2},
  vcaux::Array{T3} = VarianceComponentAuxData(vcdata)# map(x -> VarianceComponentAuxData(x), vcdata)
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate},
    T3 <: VarianceComponentAuxData}

  copyto!(∇, gradient(vcm, vcdata, vcaux))
end

function gradient(
  vcm::T1,
  vcdata::Array{T2},
  vcaux::Array{T3} = VarianceComponentAuxData(vcdata)
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate},
    T3 <: VarianceComponentAuxData}

  mapreduce(x -> gradient(vcm, x...), +, zip(vcdata, vcaux))
end

#---------------------------------------------------------------------------#
# Evaluate Fisher information matrix with respect to Σ
#---------------------------------------------------------------------------#

"""
    fisher_Σ!(H, vcmrot, vcobsrot)

Calculate Fisher information matrix at `Σ = (Σ[1], Σ[2])` and overwrite `H`
based on *rotated* two variance component model `vcmrot` and *rotated* two
variance component data `vcobsrot`.

# Input
- `H`: Hessian matrix.
- `vcmrot::TwoVarCompModelRotate`: *rotated* two variance component model.
- `vcobsrot::TwoVarCompVariateRotate`: *rotated* two variance component data.

# Output
- `H`: Fisher information matrix at `Σ = (Σ[1], Σ[2])`.
"""
function fisher_Σ!(
  H::AbstractMatrix{T},
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T}
  ) where {T <: AbstractFloat}

  n, d = size(vcobsrot.Yrot, 1), size(vcobsrot.Yrot, 2)
  zeroT, oneT = zero(T), one(T)
  fill!(H, zeroT)
  # evaluate Hessian with respect to Σ[1], Σ[2]
  C = zeros(T, d, d)
  Φ2 = kron(vcmrot.eigvec, vcmrot.eigvec)
  # (1, 1) block
  λi, λj, deval = zeroT, zeroT, zeroT
  @inbounds for j in 1:d, i in j:d
    λi, λj = vcmrot.eigval[i], vcmrot.eigval[j]
    @simd for obs in 1:n
      deval = vcobsrot.eigval[obs]
      C[i, j] += deval * deval / (λi * deval + oneT) / (λj * deval + oneT)
    end
  end
  LinearAlgebra.copytri!(C, 'L')
  mul!(view(H, 1:d^2, 1:d^2), Φ2 * Diagonal(vec(C)), transpose(Φ2))
  # (2, 1) block
  fill!(C, zeroT)
  @inbounds for j in 1:d, i in j:d
    λi, λj = vcmrot.eigval[i], vcmrot.eigval[j]
    @simd for obs in 1:n
      deval = vcobsrot.eigval[obs]
      C[i, j] += deval / (λi * deval + oneT) / (λj * deval + oneT)
    end
  end
  LinearAlgebra.copytri!(C, 'L')
  mul!(view(H, (d^2+1):(2d^2), 1:d^2), Φ2 * Diagonal(vec(C)), transpose(Φ2))
  # d-by-d (2, 2) block
  fill!(C, zeroT)
  @inbounds for j in 1:d, i in j:d
    λi, λj = vcmrot.eigval[i], vcmrot.eigval[j]
    @simd for obs in 1:n
      deval = vcobsrot.eigval[obs]
      C[i, j] += oneT / (λi * deval + oneT) / (λj * deval + oneT)
    end
  end
  LinearAlgebra.copytri!(C, 'L')
  mul!(view(H, (d^2+1):(2d^2), (d^2+1):(2d^2)), Φ2 * Diagonal(vec(C)), transpose(Φ2))
  # copy to upper triangular part
  LinearAlgebra.copytri!(H, 'L')
  rmul!(H, convert(T, 0.5))
  # output
  return H
end # function fisher!

function fisher_Σ(
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T}
  ) where {T <: AbstractFloat}

  d = length(vcmrot.eigval)
  H = zeros(T, 2d^2, 2d^2)
  fisher_Σ!(H, vcmrot, vcobsrot)
end

function fisher_Σ!(
  H::AbstractMatrix{T},
  vcm::VarianceComponentModel{T, 2},
  vcobsrot::TwoVarCompVariateRotate{T}
  ) where {T <: AbstractFloat}

  fisher_Σ!(H, TwoVarCompModelRotate(vcm), vcobsrot)
end

function fisher_Σ!(
  H::AbstractMatrix{T},
  vcm::VarianceComponentModel{T, 2},
  vcobsrot::VarianceComponentVariate{T, 2}
  ) where {T <: AbstractFloat}

  fisher_Σ!(H, TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobsrot))
end

function fisher_Σ(
  vcm::VarianceComponentModel{T, 2},
  vcobsrot::VarianceComponentVariate{T, 2}
  ) where {T <: AbstractFloat}

  fisher_Σ(TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobsrot))
end

function fisher_Σ(
  vcm::VarianceComponentModel{T, 2},
  vcobsrot::TwoVarCompVariateRotate{T}
  ) where {T <: AbstractFloat}

  fisher_Σ(TwoVarCompModelRotate(vcm), vcobsrot)
end

function fisher_Σ!(
  H::AbstractMatrix,
  vcm::T1,
  vcobs::Array{T2}
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate}}

  copyto!(H, fisher_Σ(vcm, vcobs))
end

function fisher_Σ(
  vcm::T1,
  vcobs::Array{T2}
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate}}

  mapreduce(x -> fisher_Σ(vcm, x), +, vcobs)
end

#---------------------------------------------------------------------------#
# Evaluate Fisher information matrix with respect to B
#---------------------------------------------------------------------------#

"""

    fisher_B!(H, vcmrot, vcobsrot)

Calculate Fisher information matrix of `B` and overwrite `H`
based on two variance component model `vcmrot` and *rotated* two
variance component data `vcobsrot`.

# Input
- `H`: Hessian matrix.
- `vcmrot::VarianceComponentModel{T, 2}`: two variance component model.
- `vcobsrot::TwoVarCompVariateRotate`: *rotated* two variance component data.

# Output
- `H`: Fisher information matrix at `B`.
"""
function fisher_B!(
  H::AbstractMatrix{T},
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobsrot)
  ) where {T <: AbstractFloat}

  zeroT, oneT = zero(T), one(T)
  # working X
  fill!(vcaux.Xwork, zeroT)
  kronaxpy!(vcmrot.eigvec', vcobsrot.Xrot, vcaux.Xwork)
  # working weights
  fill!(vcaux.obswt, zeroT)
  kronaxpy!(vcobsrot.eigval, vcmrot.eigval, vcaux.obswt)
  @inbounds @simd for i in eachindex(vcaux.obswt)
    vcaux.obswt[i] = oneT / √(vcaux.obswt[i] + oneT)
  end
  # weighted least squares
  lmul!(Diagonal(vcaux.obswt), vcaux.Xwork)
  mul!(H, transpose(vcaux.Xwork), vcaux.Xwork)
end

function fisher_B(
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T}
  ) where {T <: AbstractFloat}

  H = zeros(T, nmeanparams(vcmrot), nmeanparams(vcmrot))
  fisher_B!(H, vcmrot, vcobsrot)
end

function fisher_B!(
  H::AbstractMatrix{T},
  vcm::Union{VarianceComponentModel{T, 2}, TwoVarCompVariateRotate{T}},
  vcobs::Union{TwoVarCompVariateRotate{T}, VarianceComponentVariate{T, 2}},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobs)
  ) where {T <: AbstractFloat}

  fisher_B!(H, TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobs), vcaux)
end

function fisher_B(
  vcm::T1,
  vcobs::T2,
  vcaux::T3 = VarianceComponentAuxData(vcobs)
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate},
    T3 <: VarianceComponentAuxData}

  H = zeros(eltype(vcm), nmeanparams(vcm), nmeanparams(vcm))
  fisher_B!(H, TwoVarCompModelRotate(vcm), TwoVarCompVariateRotate(vcobs), vcaux)
end

function fisher_B!(
  H::AbstractMatrix,
  vcm::T1,
  vcdata::Array{T2},
  vcaux::Array{T3} = VarianceComponentAuxData(vcdata)
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate},
    T3 <: VarianceComponentAuxData}

  copyto!(H, fisher_B(vcm, vcdata, vcaux))
end

function fisher_B(
  vcm::T1,
  vcdata::Array{T2},
  vcaux::Array{T3} = VarianceComponentAuxData(vcdata)
  ) where {
    T1 <: Union{VarianceComponentModel, TwoVarCompModelRotate},
    T2 <: Union{VarianceComponentVariate, TwoVarCompVariateRotate},
    T3 <: VarianceComponentAuxData}

  mapreduce(x -> fisher_B(vcm, x...), +, zip(vcdata, vcaux))
end

#---------------------------------------------------------------------------#
# Update mean parameters by quadratic programming
#---------------------------------------------------------------------------#

"""
Compute the quadratic and linear parts of weighted least squares criterion.
"""
function suffstats_for_B(
  vcmrot::TwoVarCompModelRotate{T},
  vcobsrot::TwoVarCompVariateRotate{T},
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobsrot)
  ) where {T <: AbstractFloat}

  zeroT, oneT = zero(T), one(T)
  # working weights
  fill!(vcaux.obswt, zeroT)
  kronaxpy!(vcmrot.eigval, vcobsrot.eigval, vcaux.obswt)
  @inbounds @simd for i in eachindex(vcaux.obswt)
    vcaux.obswt[i] = oneT / √(vcaux.obswt[i] + oneT)
  end
  # working X
  fill!(vcaux.Xwork, zeroT)
  kronaxpy!(vcmrot.eigvec', vcobsrot.Xrot, vcaux.Xwork)
  lmul!(Diagonal(vcaux.obswt), vcaux.Xwork)
  # working y
  mul!(vcaux.res, vcobsrot.Yrot, vcmrot.eigvec)
  @inbounds @simd for i in eachindex(vcaux.ywork)
    vcaux.ywork[i] = vcaux.res[i] * vcaux.obswt[i]
  end
  # output
  return transpose(vcaux.Xwork) * vcaux.Xwork, transpose(vcaux.Xwork) * vcaux.ywork
end

function suffstats_for_B(
  vcmrot::T1,
  vcdatarot::Array{T2},
  vcaux::Array{T3}
  ) where {
    T1 <: TwoVarCompModelRotate,
    T2 <: TwoVarCompVariateRotate,
    T3 <: VarianceComponentAuxData}

  mapreduce(x -> suffstats_for_B(vcmrot, x...), +,
    zip(vcdatarot, vcaux))
end

function +(
  x::Tuple{AbstractMatrix, AbstractVector},
  y::Tuple{AbstractMatrix, AbstractVector}
  )

  x[1] + y[1], x[2] + y[2]
end

"""
Update mean parameters `vcm.B` by quadratic programming.

## Input
- `vcm::VarianceComponentModel{T, 2}`: two variance component model
- `vcdatarot::TwoVarCompVariateRotate`: *rotated* two variance component data
- `qpsolver`: QP solver
- `Q`: quadratic part in the QP
- `c`: linear part in the QP

## Output
- `vcm.B`: updated mean parameters
"""
function update_meanparam!(
  vcm::T1,
  vcdatarot::Union{T2, Array{T2}},
  qpsolver::MathProgBase.SolverInterface.AbstractMathProgSolver
    = MathProgBase.defaultQPsolver,
  vcaux::Union{T3, Array{T3}} = VarianceComponentAuxData(vcdatarot)
  ) where {
    T1 <: VarianceComponentModel,
    T2 <: TwoVarCompVariateRotate,
    T3 <: VarianceComponentAuxData}

  # quick return if there is no mean parameters
  isempty(vcm.B) && return vcm.B
  # rotate the model
  vcmrot = TwoVarCompModelRotate(vcm)
  # accumlate the quadratic and linear parts of QP
  Q, c = suffstats_for_B(vcmrot, vcdatarot, vcaux)
  # quadratic programming
  qpsol = quadprog(-c, Q, vcm.A, vcm.sense, vcm.b, vcm.lb, vcm.ub, qpsolver)
  if qpsol.status ≠ :Optimal
    println("Error in quadratic programming $(qpsol.status)")
  end
  copyto!(vcm.B, qpsol.sol)
  return vcm.B
end

#---------------------------------------------------------------------------#
# Fisher scoring algorithm
#---------------------------------------------------------------------------#

"""
`TwoVarCompOptProb` holds data and model for solving two variance component
problem

# Fields
- `vcmodel`: [`VarianceComponentModel`](@ref)
- `qcdatarot`: [`TwoVarCompVariateRotate`](@ref)
- `qpsolver`: `MathProgBase.SolverInterface.AbstractMathProgSolver`
- `L`: working variable, a tuple of two Cholesky factors `(L[1], L[2])`
- `∇Σ`: working variable, gradient with respect to `(Σ[1], Σ[2])`
- `HΣ`: working variable, Hessian with respect to `(Σ[1], Σ[2])`
- `HL`: working variable, Hessian with respect to `(L[1], L[2])`
- `Q`: working variable, quadratic part of QP for updating `B`
- `c`: working variable, linear part of QP for updating `B`
"""
mutable struct TwoVarCompOptProb{
  T1 <: VarianceComponentModel,
  T2 <: Union{TwoVarCompVariateRotate, AbstractArray},
  T3 <: AbstractMatrix,
  T4 <: AbstractVector,
  T5 <: Union{VarianceComponentAuxData, AbstractArray}} <: MathProgBase.AbstractNLPEvaluator
  # variance component model and data
  vcmodel::T1
  vcdatarot::T2
  # QP solver for updating B given Σ
  qpsolver::MathProgBase.SolverInterface.AbstractMathProgSolver
  # intermediate variables
  L::NTuple{2, T3} # Cholesky factors
  ∇Σ::T4           # graident wrt (Σ1, Σ2)
  HΣ::T3           # Hessian wrt (Σ1, Σ2)
  HL::T3           # Hessian wrt (L1, L2)
  vcaux::T5
end

"""
    TwoVarCompOptProb(vcm, vcobsrot)
    TwoVarCompOptProb(vcm, vcobsrot, qpsolver)

Constructor of [`TwoVarCompOptProb`](@ref) from [`VarianceComponentModel`](@ref) `vcm`,
[`TwoVarCompVariateRotate`](@ref) `vcobsrot`, and input quadratic programming
sovler `qpsolver`.
"""
function TwoVarCompOptProb(
  vcm::T1,
  vcdatarot::Union{T2, Array{T2}},
  qpsolver::MathProgBase.SolverInterface.AbstractMathProgSolver
    = MathProgBase.defaultQPsolver
  ) where {
    T1 <: VarianceComponentModel,
    T2 <: TwoVarCompVariateRotate}

  T = eltype(vcm)
  d, pd = length(vcm), nmeanparams(vcm)
  # number of optimization parameters in variance
  nvar = nvarparams(vcm)
  # allocate intermediate variables
  L  = (zeros(T, d, d), zeros(T, d, d))
  ∇Σ = zeros(T, 2d^2)       # graident wrt (Σ1, Σ2)
  HΣ = zeros(T, 2d^2, 2d^2) # Hessian wrt (Σ1, Σ2)
  HL = zeros(T, nvar, nvar) # Hessian wrt Ls
  vcaux = VarianceComponentAuxData(vcdatarot)

  # constructor
  TwoVarCompOptProb{typeof(vcm), typeof(vcdatarot), typeof(HΣ),
    typeof(∇Σ), typeof(vcaux)}(vcm, vcdatarot, qpsolver, L, ∇Σ, HΣ, HL, vcaux)
end

"""
    optimparam_to_vcparam(dd, x)

Translate optimization variables `x` to variance component parmeters
`dd.vcmodel.Σ` in [`TwoVarCompOptProb`](@ref).
"""
function optimparam_to_vcparam!(dd::TwoVarCompOptProb, x::Vector{T}) where {T}

  d        = size(dd.L[1], 1)
  nvar     = nvarparams(dd.vcmodel)
  nvarhalf = div(nvar, 2)
  # variance parameter
  dd.L[1][trilind(d)] = x[1:nvarhalf]
  dd.L[2][trilind(d)] = x[(nvarhalf+1):end]
  # exponentiate diagonal entries
  @inbounds @simd for i in 1:d
    dd.L[1][i, i] = exp(dd.L[1][i, i])
    dd.L[2][i, i] = exp(dd.L[2][i, i])
  end
  # Σi = Li Li'
  mul!(dd.vcmodel.Σ[1], dd.L[1], transpose(dd.L[1]))
  mul!(dd.vcmodel.Σ[2], dd.L[2], transpose(dd.L[2]))
  # make sure the last variance component is pos. def.
  ϵ = convert(T, 1e-8)
  clamp_diagonal!(dd.vcmodel.Σ[2], ϵ, T(Inf))
end

function MathProgBase.initialize(dd::TwoVarCompOptProb,
  requested_features::Vector{Symbol})

  for feat in requested_features
    if !(feat in [:Grad, :Jac, :Hess])
      error("Unsupported feature $feat")
    end
  end
end # function MathProgBase.initialize

MathProgBase.features_available(dd::TwoVarCompOptProb) = [:Grad, :Jac, :Hess]
MathProgBase.eval_g(dd::TwoVarCompOptProb, g, x) = nothing
MathProgBase.jac_structure(dd::TwoVarCompOptProb) = Int[], Int[]
MathProgBase.eval_jac_g(dd::TwoVarCompOptProb, J, x) = nothing

function MathProgBase.eval_f(dd::TwoVarCompOptProb, x::Vector{T}) where {T}

  # update variance parameter from optim variable x
  optimparam_to_vcparam!(dd, x)
  # update mean parameters
  update_meanparam!(dd.vcmodel, dd.vcdatarot, dd.qpsolver, dd.vcaux)
  # evaluate profile log-pdf
  logpdf(dd.vcmodel, dd.vcdatarot, dd.vcaux)
end # function MathProgBase.eval_f

function MathProgBase.eval_grad_f(
  dd::TwoVarCompOptProb,
  grad_f::Vector{T},
  x::Vector{T}
  ) where {T}

  d        = size(dd.L[1], 1)
  nvar     = nvarparams(dd.vcmodel)
  nvarhalf = div(nvar, 2)
  # update variance parameter from optim variable x
  optimparam_to_vcparam!(dd, x)
  # update mean parameters
  update_meanparam!(dd.vcmodel, dd.vcdatarot, dd.qpsolver, dd.vcaux)
  # gradient wrt (Σ[1], Σ[2])
  gradient!(dd.∇Σ, dd.vcmodel, dd.vcdatarot, dd.vcaux)
  # chain rule for gradient wrt Cholesky factor
  chol_gradient!(view(grad_f, 1:nvarhalf), dd.∇Σ[1:d^2], dd.L[1])
  chol_gradient!(view(grad_f, (nvarhalf+1):nvar), dd.∇Σ[(d^2+1):end], dd.L[2])
  # chain rule for exponential of diagonal entries
  @inbounds for j in 1:d
    # linear index of diagonal entries of L1
    idx = 1 + (j - 1) * d - div((j - 1) * (j - 2), 2)
    grad_f[idx] *= dd.L[1][j, j]
    # linear index of diagonal entries of L2
    idx += nvarhalf
    grad_f[idx] *= dd.L[2][j, j]
  end
  return grad_f
end # function MathProgBase.eval_grad_f

function MathProgBase.hesslag_structure(dd::TwoVarCompOptProb)

  nvar = nvarparams(dd.vcmodel)
  # linear indices for variance parameters
  ## ind2sub((nvar, nvar), trilind(nvar))
  #Tuple(CartesianIndices((nvar, nvar))[trilind(nvar)])
  arr1, arr2 = Int64[], Int64[]
  for i in 1:nvar
    arr1 = vcat(arr1, i:nvar)
    arr2 = vcat(arr2, fill(i, nvar-i+1))
  end
  return (arr1, arr2)
end # function MathProgBase.hesslag_structure

function MathProgBase.eval_hesslag(dd::TwoVarCompOptProb, H::Vector{T},
  x::Vector{T}, σ::T, μ::Vector{T}) where {T}

  d        = size(dd.L[1], 1)
  nvar     = nvarparams(dd.vcmodel)
  nvarhalf = div(nvar, 2)
  # update variance parameter from optim variable x
  optimparam_to_vcparam!(dd, x)
  # update mean parameters
  update_meanparam!(dd.vcmodel, dd.vcdatarot, dd.qpsolver, dd.vcaux)
  # Hessian wrt (Σ1, Σ2)
  fisher_Σ!(dd.HΣ, dd.vcmodel, dd.vcdatarot)
  # chain rule for Hessian wrt Cholesky factor
  # only the lower left triangle
  # (1, 1) block
  chol_gradient!(view(dd.HL, 1:nvarhalf, 1:nvarhalf),
    chol_gradient(dd.HΣ[1:d^2, 1:d^2], dd.L[1])',
    dd.L[1])
  # (2, 1) block
  chol_gradient!(view(dd.HL, nvarhalf+1:nvar, 1:nvarhalf),
    chol_gradient(dd.HΣ[(d^2+1):(2d^2), 1:d^2], dd.L[1])',
    dd.L[2])
  # (2, 2) block
  chol_gradient!(view(dd.HL, nvarhalf+1:nvar, nvarhalf+1:nvar),
    chol_gradient(dd.HΣ[(d^2+1):(2d^2), (d^2+1):(2d^2)], dd.L[2])',
    dd.L[2])
  # chain rule for exponential of diagonal entries
  @inbounds for j in 1:d
    # linear index of diagonal entries of L1
    idx = 1 + (j - 1) * d - div((j - 1) * (j - 2), 2)
    for i in 1:nvar
      dd.HL[i, idx] *= dd.L[1][j, j]
      dd.HL[idx, i] *= dd.L[1][j, j]
    end
    # linear index of diagonal entries of L2
    idx += nvarhalf
    for i in 1:nvar
      dd.HL[i, idx] *= dd.L[2][j, j]
      dd.HL[idx, i] *= dd.L[2][j, j]
    end
  end
  # output
  copyto!(H, vech(dd.HL))
  rmul!(H, -σ)

end

"""
    mle_fs!(vcmodel, vcdatarot; maxiter, solver, qpsolver, verbose)

Find MLE by Fisher scoring algorithm.

# Input
- `vcmodel`: two variane component model [`VarianceComponentModel`](@ref), with
 `vcmodel.B` and `vcmodel.Σ` used as starting point
- `vcdatarot`: rotated two varianec component data [`TwoVarCompVariateRotate`](@ref)

# Keyword
- `maxiter::Int`: maximum number of iterations, default is 1000
- `solver::Symbol`: backend nonlinear programming solver, `:Ipopt` (default) or `:Knitro`
- `qpsolver::Symbol`: backend quadratic programming solver, `:Ipopt` (default) or `:Gurobi` or `Mosek`
- `verbose::Bool`: display information

# Output
- `maxlogl`: log-likelihood at solution
- `vcmodel`: [`VarianceComponentModel`](@ref) with updated model parameters
- `Σse=(Σse[1],Σse[2])`: standard errors of estimate `Σ=(Σ[1],Σ[2])`
- `Σcov`: covariance matrix of estimate `Σ=(Σ[1],Σ[2])`
- `Bse`: standard errors of estimate `B`
- `Bcov`: covariance of estimate `B`
"""
function mle_fs!(
  vcmodel::T1,
  vcdatarot::Union{T2, Array{T2}};
  maxiter::Integer = 1000,
  solver::Symbol = :Ipopt,
  qpsolver::Symbol = :Ipopt,
  verbose::Bool = true,
  ) where {T1<:VarianceComponentModel, T2<:TwoVarCompVariateRotate}

  T = eltype(vcmodel)
  d = length(vcmodel)
  Ltrilind = trilind(d, d)
  # number of mean parameters
  nmean = nmeanparams(vcmodel)
  # number of optimization parameters in Cholesky factors
  nvar = nvarparams(vcmodel)
  nvarhalf = div(nvar, 2)
  # pre-allocate variables for optimization
  zeroT = convert(T, 0)
  # data for the optimization problem
  if qpsolver == :Ipopt
    qs = IpoptSolver(print_level = 0)
  elseif qpsolver == :Gurobi
    qs = GurobiSolver(OutputFlag = 0)
  elseif qpsolver == :Mosek
    qs = MosekSolver(MSK_IPAR_LOG = 0)
  end
  dd = TwoVarCompOptProb(vcmodel, vcdatarot, qs)
  # set up MathProgBase interface
  if solver == :Ipopt
    # see http://www.coin-or.org/Ipopt/documentation/documentation.html for IPOPT
    solver = IpoptSolver(
      hessian_approximation = "exact",
      tol = 1.0e-8, # default is 1.0e-8
      acceptable_tol = 1.0e-5, # default is 1.0e-6
      max_iter = maxiter, # default is 3000
      print_frequency_iter = 5, # default is 1
      print_level = verbose ? 5 : 0,
      print_info_string = "yes",
      #nlp_scaling_method = "none",
      #derivative_test = "second-order",
      #linear_solver = "mumps",
      #linear_solver = "pardiso",
      mehrotra_algorithm = "yes",
      )
  elseif solver == :Mosek
    # see http://docs.mosek.com/7.0/capi/Parameters.html for Mosek options
    solver = MosekSolver(
      MSK_IPAR_INTPNT_MAX_ITERATIONS = maxiter,
      MSK_DPAR_INTPNT_NL_TOL_REL_GAP = 1.0e-8,
      MSK_IPAR_LOG = verbose ? 10 : 0, # deafult value is 10
      #MSK_IPAR_OPTIMIZER = MSK_OPTIMIZER_NONCONVEX,
      #MSK_IPAR_LOG_NONCONVEX = 20,
      #MSK_IPAR_NONCONVEX_MAX_ITERATIONS = 100,
      #MSK_DPAR_INTPNT_NL_TOL_NEAR_REL = 1e8,
      #MSK_IPAR_LOG_CHECK_CONVEXITY = 1,
      #MSK_IPAR_INFEAS_PREFER_PRIMAL = MSK_OFF
      )
  elseif solver == :Knitro
    # see https://www.artelys.com/tools/knitro_doc/3_referenceManual/userOptions.html for Mosek options
    solver = KnitroSolver(
      KTR_PARAM_ALG = 0,
      KTR_PARAM_OUTLEV = verbose ? 2 : 0,
      #KTR_PARAM_MAXCGIT = 5,
      #KTR_PARAM_SCALE = 0,
      #KTR_PARAM_HONORBNDS = 1,
      #KTR_PARAM_GRADOPT = 1,
      #KTR_PARAM_HESSOPT = 1,
      #KTR_PARAM_DERIVCHECK = 3,
      #KTR_PARAM_TUNER = 1,
      #KTR_PARAM_MAXTIMECPU = 5.0,
      #KTR_PARAM_BAR_MURULE = 4,
      #KTR_PARAM_BAR_FEASIBLE = 3,
      #KTR_PARAM_BAR_FEASMODETOL = 1.0e-8,
      #KTR_PARAM_BAR_SWITCHRULE = 3,
      #KTR_PARAM_BAR_PENCONS = 2,
      #KTR_PARAM_BAR_DIRECTINTERVAL = 0,
      #KTR_PARAM_BAR_MAXBACKTRACK = 10,
      #KTR_PARAM_BAR_MAXCROSSIT = 5,
      #KTR_PARAM_BAR_MAXREFACTOR = 5,
      #KTR_PARAM_BAR_REFINEMENT = 1,
      #KTR_PARAM_BAR_WATCHDOG = 1,
      )
  end
  m = MathProgBase.NonlinearModel(solver)
  # lower and upper bounds for variance parameter
  lb = zeros(nvar); fill!(lb, -T(Inf))
  ub = zeros(nvar); fill!(ub,  T(Inf))
  MathProgBase.loadproblem!(m, nvar, 0, lb, ub, T[], T[], :Max, dd)
  # start point
  copyto!(dd.L[1], cholesky(vcmodel.Σ[1], Val(true)).L.data)
  copyto!(dd.L[2], cholesky(vcmodel.Σ[2], Val(true)).L.data)
  # reparameterize diagonal entries to exponential
  @inbounds @simd for i in 1:d
    dd.L[1][i, i] = log(dd.L[1][i, i] + convert(T, 1e-8))
    dd.L[2][i, i] = log(dd.L[2][i, i] + convert(T, 1e-8))
  end
  x0 = [vech(dd.L[1]); vech(dd.L[2])]
  MathProgBase.setwarmstart!(m, x0)
  # optimize
  MathProgBase.optimize!(m)
  stat = MathProgBase.status(m)
  x = MathProgBase.getsolution(m)
  # retrieve variance component parameters
  optimparam_to_vcparam!(dd, x)
  # update mean parameters
  update_meanparam!(dd.vcmodel, dd.vcdatarot, dd.qpsolver, dd.vcaux)
  # update final objective value
  maxlogl = logpdf(vcmodel, vcdatarot, dd.vcaux)

  # standard errors
  Bcov = zeros(T, nmean, nmean)
  Bcov = pinv(fisher_B(vcmodel, vcdatarot, dd.vcaux))
  Bse = similar(vcmodel.B)
  copyto!(Bse, sqrt.(diag(Bcov)))
  Σcov = pinv(fisher_Σ(vcmodel, vcdatarot))
  Σse = (zeros(T, d, d), zeros(T, d, d))
  copyto!(Σse[1], sqrt.(diag(view(Σcov, 1:d^2, 1:d^2))))
  copyto!(Σse[2], sqrt.(diag(view(Σcov, d^2+1:2d^2, d^2+1:2d^2))))

  # output
  maxlogl, vcmodel, Σse, Σcov, Bse, Bcov
end # function mle_fs

#---------------------------------------------------------------------------#
# MM algorithm
#---------------------------------------------------------------------------#

function suffstats_for_Σ!(
  vcmrot::TwoVarCompModelRotate,
  vcobsrot::TwoVarCompVariateRotate,
  vcaux::VarianceComponentAuxData = VarianceComponentAuxData(vcobsrot)
  )

  T = eltype(vcmrot)
  n, d = size(vcobsrot.Yrot, 1), size(vcobsrot.Yrot, 2)
  A1, b1 = zeros(T, d, d), zeros(T, d)
  A2, b2 = zeros(T, d, d), zeros(T, d)
  zeroT, oneT = zero(T), one(T)
  λj, tmp     = zeroT, zeroT
  # (rotated) residual
  residual!(vcaux.res, vcmrot, vcobsrot)
  # sufficient statistics for Σ2
  @inbounds for j in 1:d
    λj = vcmrot.eigval[j]
    for i in 1:n
      tmp = oneT / (vcobsrot.eigval[i] * λj + oneT)
      vcaux.res[i, j] *= tmp
      b2[j] += tmp
    end
  end
  mul!(A2, transpose(vcaux.res), vcaux.res)
  # sufficient statistics for Σ1
  @inbounds for j in 1:d
    λj = vcmrot.eigval[j]
    for i in 1:n
      tmp = vcobsrot.eigval[i]
      vcaux.res[i, j] *= √tmp * λj
      b1[j] += tmp / (tmp * λj + oneT)
    end
  end
  mul!(A1, transpose(vcaux.res), vcaux.res)
  # output
  return A1, b1, A2, b2
end

function suffstats_for_Σ!(
  vcmrot::T1,
  vcdatarot::Array{T2},
  vcaux::Array{T3} = VarianceComponentAuxData(vcdatarot)
  ) where {
    T1 <: TwoVarCompModelRotate,
    T2 <: TwoVarCompVariateRotate,
    T3 <: VarianceComponentAuxData}

  mapreduce(x -> suffstats_for_Σ!(vcmrot, x...), +, zip(vcdatarot, vcaux))
end

function +(
  x::Tuple{AbstractMatrix, AbstractVector, AbstractMatrix, AbstractVector},
  y::Tuple{AbstractMatrix, AbstractVector, AbstractMatrix, AbstractVector}
  )

  x[1] + y[1], x[2] + y[2], x[3] + y[3], x[4] + y[4]
end

function mm_update_Σ!(
  vcm::T1,
  vcdatarot::Union{T2, Array{T2}},
  vcaux::Union{T3, Array{T3}} = VarianceComponentAuxData(vcdatarot)
  ) where {
    T1 <: VarianceComponentModel,
    T2 <: TwoVarCompVariateRotate,
    T3 <: VarianceComponentAuxData}

  T, d = eltype(vcm), length(vcm)
  zeroT, oneT = zero(T), one(T)
  # eigen-decomposition of (vcm.Σ[1], vcm.Σ[2])
  vcmrot = TwoVarCompModelRotate(vcm)
  # sufficient statistics for updating Σ1, Σ2
  A1, b1, A2, b2 = suffstats_for_Σ!(vcmrot, vcdatarot, vcaux)
  @inbounds for j in 1:d
    b1[j] = √b1[j]
    b2[j] = √b2[j]
  end
  Φinv = inv(vcmrot.eigvec)
  # update Σ1
  lmul!(Diagonal(b1), A1), rmul!(A1, Diagonal(b1))
  storage = eigen!(Symmetric(A1))
  @inbounds for i in 1:d
    storage.values[i] = storage.values[i] > zeroT ? √√storage.values[i] : zeroT
  end
  rmul!(storage.vectors, Diagonal(storage.values)) #scale!(storage.vectors, storage.values)
  lmul!(Diagonal(oneT ./ b1), storage.vectors) #scale!(storage.vectors, storage.values)
  mul!(vcm.Σ[1], transpose(Φinv), storage.vectors)
  copyto!(vcm.Σ[1], vcm.Σ[1] * transpose(vcm.Σ[1]))
  # update Σ2
  lmul!(Diagonal(b2), A2), rmul!(A2, Diagonal(b2))
  storage = eigen!(Symmetric(A2))
  @inbounds for i in 1:d
    storage.values[i] = storage.values[i] > zeroT ? √√storage.values[i] : zeroT
  end
  rmul!(storage.vectors, Diagonal(storage.values)) #scale!(storage.vectors, storage.values)
  lmul!(Diagonal(oneT ./ b2), storage.vectors) #scale!(storage.vectors, storage.values)
  mul!(vcm.Σ[2], transpose(Φinv), storage.vectors)
  copyto!(vcm.Σ[2], vcm.Σ[2] * transpose(vcm.Σ[2]))
end

"""
    mle_mm!(vcmodel, vcdatarot; maxiter, qpsolver, verbose)

Find MLE by minorization-maximization (MM) algorithm.

# Input
- `vcmodel`: two variane component model [`VarianceComponentModel`](@ref), with
 `vcmodel.B` and `vcmodel.Σ` used as starting point
- `vcdatarot`: rotated two varianec component data [`TwoVarCompVariateRotate`](@ref)

# Keyword
- `maxiter::Int`: maximum number of iterations, default is 1000
- `qpsolver::Symbol`: backend quadratic programming solver, `:Ipopt` (default) or `:Gurobi` or `Mosek`
- `verbose::Bool`: display information

# Output
- `maxlogl`: log-likelihood at solution
- `vcmodel`: [`VarianceComponentModel`](@ref) with updated model parameters
- `Σse=(Σse[1],Σse[2])`: standard errors of estimate `Σ=(Σ[1],Σ[2])`
- `Σcov`: covariance matrix of estimate `Σ=(Σ[1],Σ[2])`
- `Bse`: standard errors of estimate `B`
- `Bcov`: covariance of estimate `B`

# Reference
- H. Zhou, L. Hu, J. Zhou, and K. Lange (2015)
  MM algorithms for variance components models.
  [http://arxiv.org/abs/1509.07426](http://arxiv.org/abs/1509.07426)
"""
function mle_mm!(
  vcm::T1,
  vcdatarot::Union{T2, Array{T2}};
  maxiter::Integer = 10000,
  funtol::Real = convert(eltype(vcm), 1e-8),
  qpsolver::Symbol = :Ipopt,
  verbose::Bool = true,
  ) where {T1 <: VarianceComponentModel, T2 <: TwoVarCompVariateRotate}

  T, d, pd = eltype(vcm), length(vcm), nmeanparams(vcm)
  # set up QP solver
  if qpsolver == :Ipopt
    qs = IpoptSolver(print_level = 0)
  elseif qpsolver == :Gurobi
    qs = GurobiSolver(OutputFlag = 0)
  elseif qpsolver == :Mosek
    qs = MosekSolver(MSK_IPAR_LOG = 0)
  end
  # allocate intermediate variables
  vcaux = VarianceComponentAuxData(vcdatarot)
  # initial log-likelihood
  logl::T = logpdf(vcm, vcdatarot, vcaux)
  if verbose
    println()
    println("     MM Algorithm")
    println("  Iter      Objective  ")
    println("--------  -------------")
    Printf.@printf("%8.d  %13.e\n", 0, logl)
  end

  # MM loop
  for iter in 1:maxiter

    # update Σ
    mm_update_Σ!(vcm, vcdatarot, vcaux)
    # make sure the last variance component is pos. def.
    ϵ = convert(T, 1e-8)
    clamp_diagonal!(vcm.Σ[2], ϵ, T(Inf))

    # update mean parameters
    if pd > 0; update_meanparam!(vcm, vcdatarot, qs, vcaux); end

    # check convergence
    loglold = logl
    logl    = logpdf(vcm, vcdatarot, vcaux)
    if verbose
      if (iter <= 10) || (iter > 10 && iter % 10 == 0)
        Printf.@printf("%8.d  %13.e\n", iter, logl)
      end
    end
    if abs(logl - loglold) < funtol * (abs(logl) + one(T))
      break
    end
  end
  if verbose; println(); end

  # standard errors
  Bcov = pinv(fisher_B(vcm, vcdatarot, vcaux))
  Bse = similar(vcm.B)
  copyto!(Bse, sqrt.(diag(Bcov)))
  Σcov = pinv(fisher_Σ(vcm, vcdatarot))
  Σse = (zeros(T, d, d), zeros(T, d, d))
  copyto!(Σse[1], sqrt.(diag(view(Σcov, 1:d^2, 1:d^2))))
  copyto!(Σse[2], sqrt.(diag(view(Σcov, d^2+1:2d^2, d^2+1:2d^2))))

  # output
  logl, vcm, Σse, Σcov, Bse, Bcov
end # function mle_mm

#---------------------------------------------------------------------------#
# Estimation gateway
#---------------------------------------------------------------------------#

"""
    fit_mle!(vcmodel, vcdata; algo)

Find MLE of variane component model.

# Input
- `vcmodel`: two variane component model [`VarianceComponentModel`](@ref), with
 `vcmodel.B` and `vcmodel.Σ` used as starting point
- `vcdata`: two varianec component data [`VarianceComponentVariate`](@ref)

# Keyword
- `algo::Symbol`: algorithm, `:FS` (Fisher scoring) for `:MM`
(minorization-maximization algorithm)

# Output
- `maxlogl`: log-likelihood at solution
- `vcmodel`: [`VarianceComponentModel`](@ref) with updated model parameters
- `Σse=(Σse[1],Σse[2])`: standard errors of estimate `Σ=(Σ[1],Σ[2])`
- `Σcov`: covariance matrix of estimate `Σ=(Σ[1],Σ[2])`
- `Bse`: standard errors of estimate `B`
- `Bcov`: covariance of estimate `B`
"""
function fit_mle!(
  vcmodel::T1,
  vcdata::Union{T2, Array{T2}};
  algo::Symbol = :FS
  ) where {T1 <: VarianceComponentModel, T2 <: VarianceComponentVariate}

  # generalized-eigendecomposition and rotate data
  vcdatarot = TwoVarCompVariateRotate(vcdata)
  if algo == :FS
    return mle_fs!(vcmodel, vcdatarot)
  elseif algo == :MM
    return mle_mm!(vcmodel, vcdatarot)
  end
end

"""
    fit_reml!(vcmodel, vcdata; algo)

Find restricted MLE (REML) of variane component model.

# Input
- `vcmodel`: two variane component model [`VarianceComponentModel`](@ref), with
 `vcmodel.B` and `vcmodel.Σ` used as starting point
- `vcdata`: two varianec component data [`VarianceComponentVariate`](@ref)

# Keyword
- `algo::Symbol`: algorithm, `:FS` (Fisher scoring) for `:MM`
(minorization-maximization algorithm)

# Output
- `maxlogl`: log-likelihood at solution
- `vcmodel`: [`VarianceComponentModel`](@ref) with updated model parameters
- `Σse=(Σse[1],Σse[2])`: standard errors of estimate `Σ=(Σ[1],Σ[2])`
- `Σcov`: covariance matrix of estimate `Σ=(Σ[1],Σ[2])`
- `Bse`: standard errors of estimate `B`
- `Bcov`: covariance of estimate `B`
"""
function fit_reml!(
  vcmodel::T1,
  vcdata::Union{T2, Array{T2}};
  algo::Symbol = :FS,
  qpsolver::Symbol = :Ipopt
  ) where {
    T1 <: VarianceComponentModel,
    T2 <: VarianceComponentVariate}

  T, d = eltype(vcmodel), length(vcmodel)
  vcaux = VarianceComponentAuxData(vcdata)
  # solve regular least squares problem
  if qpsolver == :Ipopt
    qs = IpoptSolver(print_level = 0)
  elseif qpsolver == :Gurobi
    qs = GurobiSolver(OutputFlag = 0)
  elseif qpsolver == :Mosek
    qs = MosekSolver(MSK_IPAR_LOG = 0)
  end
  if typeof(vcdata) <: AbstractArray
    Q = kron(Matrix{T}(I, d, d), mapreduce(x -> transpose(x.X) * x.X, +, vcdata))
    c = vec(mapreduce(x -> transpose(x.X) * x.Y, +, vcdata))
  else
    Q = kron(Matrix{T}(I, d, d), transpose(vcdata.X) * vcdata.X)
    c = vec(transpose(vcdata.X) * vcdata.Y)
  end
  qpsol = quadprog(-c, Q, vcmodel.A, vcmodel.sense, vcmodel.b,
    vcmodel.lb, vcmodel.ub, qs)
  if qpsol.status ≠ :Optimal
    println("Error in quadratic programming $(qpsol.status)")
  end
  copyto!(vcmodel.B, qpsol.sol)
  # use residuals as responses
  resdata = deepcopy(vcdata)
  if typeof(vcdata) <: AbstractArray
    for i in eachindex(resdata)
      residual!(resdata[i].Y, vcmodel, vcdata[i])
      resdata[i].X = zeros(T, size(resdata[i].Y, 1), 0)
    end
  else
    residual!(resdata.Y, vcmodel, vcdata)
    resdata.X = zeros(T, size(resdata.Y, 1), 0)
  end
  resdatarot = TwoVarCompVariateRotate(resdata)
  resmodel = deepcopy(vcmodel)
  resmodel.B = zeros(T, 0, d)
  if algo == :FS
    _, _, Σse, Σcov, = mle_fs!(resmodel, resdatarot)
  elseif algo == :MM
    _, _, Σse, Σcov, = mle_mm!(resmodel, resdatarot)
  end

  # estimate mean parameters from covariance estimate
  copyto!(vcmodel.Σ[1], resmodel.Σ[1]), copyto!(vcmodel.Σ[2], resmodel.Σ[2])
  vcdatarot = TwoVarCompVariateRotate(vcdata)
  update_meanparam!(vcmodel, vcdatarot, qs, vcaux)

  # standard errors and covariance of mean parameters
  Bcov = inv(fisher_B(vcmodel, vcdatarot, vcaux))
  Bse = similar(vcmodel.B)
  copyto!(Bse, sqrt.(diag(Bcov)))

  ## output
  logpdf(vcmodel, vcdatarot), vcmodel, (Σse[1], Σse[2]), Σcov, Bse, Bcov
end

#---------------------------------------------------------------------------#
# Heritability estimation
#---------------------------------------------------------------------------#

"""
    heritability(Σ, Σcov, which=1)

Calcualte the heritability of each trait and its standard error from variance
component estimates and their covariance.

# Input
- `Σ=(Σ[1], ..., Σ[m])`: estimate of `m` `d x d` variance components.
- `Σcov`: `md^2 x md^2` covariance matrix of the variance component estiamte.
- `which`: indicator which is the additive genetic variance component.

# Output
- `h`: estimated heritability of `d` traits.
- `hse`: standard errors of the estiamted heritability.
"""
function heritability(
  Σ::NTuple{N, AbstractMatrix{T}},
  Σcov::AbstractMatrix{T},
  which::Integer = 1
  ) where {T, N}

  d = size(Σ[1], 1) # number of traits
  m = length(Σ)     # number of variance components
  @assert size(Σcov, 1) == m * d^2 "Dimensions don't match"
  h, hse = zeros(T, d), zeros(T, d)
  hgrad = zeros(T, m)
  for trait in 1:d
    σ2trait = [Σ[i][trait, trait] for i = 1:m]
    totalσ2trait = sum(σ2trait)
    # heritability
    h[trait] = σ2trait[which] / totalσ2trait
    # standard error of heritability by Delta method
    for i in 1:m
      if i == which
        hgrad[i] = (totalσ2trait - σ2trait[i]) / totalσ2trait^2
      else
        hgrad[i] = - σ2trait[i] / totalσ2trait^2
      end
    end
    dgidx = diagind(Σ[1])[trait]
    traitidx = dgidx:d^2:((m - 1) * d^2 + dgidx)
    hse[trait] = sqrt(dot(hgrad, view(Σcov, traitidx, traitidx) * hgrad))
  end
  # output
  h, hse
end # function heritability

max_stars_repo_name: Hua-Zhou/VarianceComponentModels

max_issues_repo_name: Hua-Zhou/VarianceComponentModels

max_forks_repo_name: Hua-Zhou/VarianceComponentModels


text:
abstract type GreaterElement <: Element end

abstract type GreaterBlock <: GreaterElement end

mutable struct CenterBlock <: GreaterBlock
    parameters::Union{AbstractString, Nothing}
    contents::Vector{Element}
end
mutable struct QuoteBlock <: GreaterBlock
    parameters::Union{AbstractString, Nothing}
    contents::Vector{Element}
end
mutable struct SpecialBlock <: GreaterBlock
    name::AbstractString
    parameters::Union{AbstractString, Nothing}
    contents::Vector{Element}
end

mutable struct Drawer <: GreaterElement
    name::AbstractString
    contents::Vector{Element}
end

mutable struct DynamicBlock <: GreaterElement
    name::AbstractString
    parameters::Union{AbstractString, Nothing}
    contents::Vector{Element}
end

mutable struct FootnoteDefinition <: GreaterElement
    label::AbstractString
    definition::Vector{Element}
end

# See elements.jl for InlineTask for load order reasons

abstract type List <: GreaterElement end

mutable struct Item <: GreaterElement
    bullet::AbstractString
    counterset::Union{AbstractString, Nothing}
    checkbox::Union{Char, Nothing}
    tag::Union{Vector{OrgComponent}, Nothing}
    contents::Vector{OrgComponent}
end

mutable struct UnorderedList <: List
    items::Vector{Item}
end
mutable struct OrderedList <: List
    items::Vector{Item}
end
mutable struct DescriptiveList <: List
    items::Vector{Item}
end

mutable struct PropertyDrawer <: GreaterElement
    contents::Vector{NodeProperty}
end

mutable struct Table <: GreaterElement
    rows::Vector{Union{TableRow, TableHrule}}
    formulas::Vector{AbstractString}
end


text:
#!/usr/bin/env julia

function main()
    random_test_filenames = Dict(
        "A"  => joinpath(pwd(), "test1.random.ucsc.nochrm.bed"),
        "B"  => joinpath(pwd(), "test2.random.ucsc.nochrm.bed"),
        "B_" => joinpath(pwd(), "test3.random.ucsc.bed"))

    random_methods = Dict(
        "coitrees"      => "./coitrees-cov",
        "cgranges"      => "./cgranges-cov",
        "CITree"        => "./CITree-cov",
    )

    sorted_test_filenames = Dict(
        "A"  => joinpath(pwd(), "test1.ucsc.nochrm.bed"),
        "B"  => joinpath(pwd(), "test2.ucsc.nochrm.bed"),
        "B_" => joinpath(pwd(), "test3.ucsc.bed"))

    sorted_methods = copy(random_methods)
    # sorted_methods["coitrees (`--sorted`)"] = "./coitrees-sorted"
    # sorted_methods["bedtools (`-sorted`)"] = "./bedtools-sorted"

    tests = [("A", "B"), ("B", "A"), ("A", "A"), ("B_", "B_")]

    test_sets = [
        ("sorted-cov-bench.csv", sorted_test_filenames, sorted_methods),
        ("random-cov-bench.csv", random_test_filenames, random_methods)
    ]

    # run testts
    for (output_filename, test_filenames, methods) in test_sets
        @show output_filename
        output = open(output_filename, "w")
        println(output, "method,test,time,mem")
        for (method_name, method) in methods
            for (test1, test2) in tests
                println((method_name, test1, test2))
                secs, kilos = runtest(method, test_filenames[test1], test_filenames[test2])
                println(
                    output, method_name, ",", test1, "vs", test2, ",",
                    round(secs, digits=1), ",", round(kilos/1000, digits=1))
                flush(output)
            end
        end
    end
end

function runtest(method, filename1, filename2)
    time_output_path, time_io = mktemp()
    run(pipeline(
        `time -f "%e,%M" -o $(time_output_path) $(method) $(filename1) $(filename2)`,
        stdout=open("/dev/null", "w")))
    time_str = String(read(time_io))
    return parse.(Float64, split(time_str, ','))
end

main()

max_stars_repo_name: dcjones/coitrees

max_issues_repo_name: dcjones/coitrees

max_forks_repo_name: dcjones/coitrees


text:
mutable struct Forest
    T :: Int
    trees :: Vector{Node}
end

Flux.@functor Forest
trainable(F :: Forest) = (trainable.(F.trees)...,)

mutable struct ForestSlice
    F :: Forest
    slice :: Vector{Int}
end

Flux.@functor ForestSlice
trainable(FS :: ForestSlice) = (trainable.(FS.F.trees[FS.slice])...,)

function Forest(num_trees :: Int, depth :: Int, hidden :: Int, num_feats :: Int, ratio :: Float64, threads :: Int = Threads.nthreads())
    tree_feats = Vector{Dict{Int, Vector{Int}}}()
    n = Int(ceil(ratio * num_feats))
    for i ∈ 1:num_trees
        ctr = 1
        feats = Dict{Int, Vector{Int}}()
        for j ∈ 1:(2^depth-1)
            rg = Array(1:num_feats)
            shuffle!(rg)
            feats[ctr] = rg[1:n]
            ctr += 1
        end
        push!(tree_feats, feats)
    end

    trees = Vector{Node}()
    for i ∈ 1:num_trees
        push!(trees, Node(tree_feats[i], hidden ,depth, 1))
    end

    Forest(threads, trees)
end

function precompute!(FS :: ForestSlice, x :: Array{Float32, 2}, y :: Vector{Int})
    Threads.@threads for n ∈ FS.F.trees[FS.slice]
        precompute!(n, x, y)
    end
end

function predict(FS :: ForestSlice, x :: Array{Float32, 2})
    preds = ones(Int, length(FS.F.trees[FS.slice]), size(x, 2))
    Threads.@threads for i ∈ 1:length(FS.F.trees[FS.slice])
        preds[i, :] = predict(FS.F.trees[FS.slice][i], x)
    end

    final_preds = rand(Int, size(x, 2))
    Threads.@threads for i ∈ 1:size(x, 2)
        if !isempty(preds[:, i])
            final_preds[i] = mode(preds[:, i])
        end
    end

    return final_preds
end

function loss!(FS :: ForestSlice, x :: Array{Float32, 2}, y :: Vector{Int})
    loss = 0
    for n ∈ FS.F.trees[FS.slice]
        loss += loss!(n, x, y)
    end
    return loss
end

function precompute!(F :: Forest, x :: Array{Float32, 2}, y :: Vector{Int})
    FS = ForestSlice(F, 1:length(F.trees))
    precompute!(FS, x, y)
end

function predict(F :: Forest, x ::Array{Float32, 2})
    FS = ForestSlice(F, 1:length(F.trees))
    predict(FS, x)
end

function loss!(F :: Forest, x :: Array{Float32, 2}, y :: Vector{Int})
    FS = ForestSlice(F, 1:length(F.trees))
    loss!(FS, x, y)
end

function importance(F :: Forest, x :: Array{Float32, 2})
    scores = Dict{Int, Float32}()
    pbar = Progress(length(F.trees))
    Threads.@threads for n ∈ F.trees
        tree_imp = importance(n, x)
        for (k, v) ∈ tree_imp
            if k ∈ keys(scores)
                scores[k] += v
            else
                scores[k] = v
            end
        end
        next!(pbar)
    end
    ρ = sum(values(scores))
    for (k, v) ∈ scores
        scores[k] = v / ρ
    end
    return scores
end

function forest_train!(epochs :: Int, F :: Forest, x :: Array{Float32, 2}, y :: Vector{Int})
    slices = Vector{ForestSlice}()
    pbar = Progress(epochs * F.T)
    per_thread = Int(ceil(length(F.trees) / F.T))
    for i ∈ 1:F.T
        FS = ForestSlice(F, Array(((i-1)*per_thread+1):min(length(F.trees), (i*per_thread))))
        push!(slices, FS)
    end

    Threads.@threads for i ∈ 1:F.T
        FS = slices[i]
        ps = params(FS)
        opt = NADAM(0.00025)

        for epoch ∈ 1:epochs
            local training_loss

            DeepForest.precompute!(FS, x, y)
            gs = gradient(ps) do
                training_loss = loss!(FS, x, y)
                return training_loss
            end
            next!(pbar; showvalues = [(:loss, training_loss/length(FS.slice)), (:acc, mean(predict(FS, x) .== y))])
            Flux.update!(opt, ps, gs)
        end
    end
end


text:
## karhunen_loeve.jl : approximate Gaussian random field generator using a Karhunen-Lo\`eve decomposition

"""
    KarhunenLoeve{n} <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) using a Karhunen-Lo\u00e8ve (KL) expansion with `n` terms.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

julia> plot_eigenvalues(grf) # plot the eigenvalue decay
[...]

julia> plot_eigenfunction(grf,4) # plots the fourth eigenfunction
[...]

```
The more terms are retained in the expansion, the better the approximation will be.
```jldoctest
julia> nterms = [1 2 5 10 20 50 100 200 500 1000]
1×10 Array{Int64,2}:
 1  2  5  10  20  50  100  200  500  1000

julia> for n in nterms
       grf = GaussianRandomField(c,KarhunenLoeve(n),pts1,pts2)
       @show rel_error(grf)
       end
rel_error(grf) = 0.7499982529722711
rel_error(grf) = 0.49999825591379987
rel_error(grf) = 0.4425751861338164
rel_error(grf) = 0.35789475278408045
rel_error(grf) = 0.16805079842673853
rel_error(grf) = 0.11187098338277579
rel_error(grf) = 0.05130466343704787
rel_error(grf) = 0.017343327476498027
rel_error(grf) = 0.0034278579378175245
rel_error(grf) = 0.0007216777927243623

```
The KL expansion is computed using the Nystrom method. Optional argument are the quadrature rule `quad` and the number of quadrature points in each direction `nq`. Possible values for `quad` are `GaussLegendre()`, `Midpoint()`, `Simpson()`, `Trapezoidal` and `EOLE()` (default). The total number of quadrature points must be larger than or equal to the requested number of terms. Default is `ceil(n^(1/d))`. This value is best left untouched.

```jldoctest
julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2,quad=GaussLegendre())
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2,quad=Simpson(),nq=20)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

```
See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`CirculantEmbedding`](@ref)
"""
struct KarhunenLoeve{n} <: GaussianRandomFieldGenerator end

KarhunenLoeve(n::Integer) = n > 0 ? KarhunenLoeve{n}() : throw(ArgumentError("in KarhunenLoeve(n), number of terms n must be positive"))

const KarhunenLoeveGRF{n} = GaussianRandomField{C,KarhunenLoeve{n}} where {n,C}
const KL{n} = KarhunenLoeve{n}

function _GaussianRandomField(mean,cov::CovarianceFunction{d},method::KL{n},pts...;
                              nq::T=ceil(typeof(n),n^(1/d)), quad::QuadratureRule=EOLE()) where {d,n,T<:Integer}
    # check if number of terms and number of quadrature points are compatible
    nq = nq > 0 ? nq .* tuple(ones(T,d)...) : length.(pts)
    nq = prod(nq) == n ? nq.+1 : nq # adjustment for ARPACK error when looking for ALL eigenvalues
    prod(nq) < n && throw(ArgumentError("too many terms requested, increase nq or lower n"))

    # determine bounding box when irregular mesh
    is_irr = typeof(pts[1]) <: AbstractMatrix
    a = is_irr ? Tuple(minimum(pts[1],dims=2)) : minimum.(pts)
    b = is_irr ? Tuple(maximum(pts[1],dims=2)) : maximum.(pts)

    # compute quadrature nodes and weights
    struc = get_nodes_and_weights.(nq,a,b,Ref(quad))
    nodes = first.(struc)
    weights = last.(struc)

    # eigenvalue problem
    C = apply(cov,nodes,nodes)
    W = d == 1 ? diagm(0 => weights...) : diagm(0 => kron(weights...))
    Wsqrt = sqrt.(W)
    B = Symmetric(Wsqrt*C*Wsqrt) # should be symmetric and positive semi-definite
    isposdef(B) || warn("equivalent eigenvalue problem is not SPD, results may be wrong or inaccurate")

    # solve
    (eigenval,eigenfunc) = eigs(B,nev=n,ritzvec=true,which=:LM,v0=randn(size(B,1)))
 	N = find_last_positive(eigenval)
	N == length(eigenval) || warn("negative eigenvalue $(eigenval[N+1]) detected, Gaussian random field will be approximated (ignoring all negative eigenvalues)")

    # compute eigenfunctions in nodes
    K = apply(cov,pts,nodes)
    Λ = diagm(0 => 1 ./ eigenval)
    eigenfunc = K*Wsqrt*eigenfunc*Λ

    # store eigenvalues and eigenfunctions
    data = SpectralData(sqrt.(eigenval[1:N]),eigenfunc[:,1:N]) # note: store sqrt of eigenval for more efficient sampling

    GaussianRandomField{typeof(cov),KL{N},typeof(pts)}(mean,cov,pts,data)
end

# returns the required dimension of the random points
randdim(grf::KarhunenLoeveGRF{n}) where {n} = n

# relative error in the KL approximation
function rel_error(grf::KarhunenLoeveGRF)
    Leb = prod(maximum.(grf.pts).-minimum.(grf.pts))
	return (std(grf.cov)^2*Leb - sum(grf.data.eigenval.^2))/(std(grf.cov)^2*Leb)
end

# sample function for both Spectral() and KarhunenLoeve(n) type
function _sample(grf::Union{SpectralGRF,KarhunenLoeveGRF}, xi)
    grf.mean + std(grf.cov)*reshape(( grf.data.eigenfunc*diagm(0 => grf.data.eigenval) )*xi,size(grf.mean))
end

show(io::IO,::KarhunenLoeve{n}) where {n} = print(io,"KL expansion with $(n) terms")


text:
# SECTION 2: data & expr prepare (CONTRIBUTION RATIO INCREASE)
# -----------------------------------
# 0.2 short presentation for solving a model
expr_RunModel = :(begin;
   Guess = ( r = 0.08, L = 0.2 );EasySearch.SteadyState!( 1, Guess, Dt, Dst, Pt, Ps, Pc, env, atol = 1E-8, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 ); Guess = ( r = 0.12, L = 0.75 );
   EasySearch.SteadyState!( env.T, Guess, Dt, Dst, Pt, Ps, Pc, env, atol = 1E-6, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 );
   PerfLog = EasySearch.Transition!( Dt, Dst, Pt, Ps, Pc, env, atol = 1E-3, MaxIter = 500, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.7, ReturnLog = true );
   EasySearch.ProcAfterTransition!( Dt, Dst, Pt, Ps, Pc, env );
   end;
)
# 0.3 short expr for extracting current demographic data
expr_GetDemogDat = :(
   Dict( "TotalPopu" => Dt[:N],
      "LaborPopu" => sum(Ps[:N][:,1:env.Sr],dims=2)[:],
      "AgingPopu_65plus" => sum(Ps[:N][:, (65 - env.START_AGE + 1):env.S],dims=2)[:],
      "AgingPopuRatio" => sum(Ps[:N][:, (65 - env.START_AGE + 1):env.S],dims=2)[:] ./ sum(Ps[:N][:,1:env.Sr],dims=2)[:]
   )
)
# 0.4 refresh demography to benchmark
f_resetEnv( demogcsvpath::String ) = (
   T = 400, S = 80, Sr = 40, START_AGE = 20, START_YEAR = 1945,
   PATH_DEMOGCSV = demogcsvpath,
   PATH_WAGEPROFILE = "./data/WageProfileCoef.csv", PATH_MA2MB = "./data/MA2MBCoef.csv",
   PATH_M2C = "./data/M2C.csv", PATH_TFPGROW = "./data/tfpGrowthProfile.csv",
)
env = f_resetEnv( "./data/Demography_base.csv" )
# 0.5 year range to plot
idx_year2plot = 2010:2110
idx_plot = idx_year2plot .- env.START_YEAR .+ 1
# 0.6 reforming macro
macro reform( ReformYear, ReformLev )
   return :( Pt[:ζ][ ( $ReformYear - env.START_YEAR + 1 ) : end] .= $ReformLev )
end # end macro

# SECTION 2: REFORM - contribution rate rising (Rising of firm’s contribution rate to UE-BMI)
# --------------------------------
# 2.1 prepare a dict to store result vars
# NOTE: each member of res_ContriRatRise is a list of lists,
# i.e. typeof(res_ContriRatRise) == Dict{String,Vector{Vector{Float64}}}
res_ContriRatRise = Dict(
   "gap/exp" => [],
   "gap/in" => [],
   "gap/gdp" => [],
   "gap/taxrev" => [],
)
# 2.2 initialize
reform_zeta_levs = [ 0.06, 0.08, 0.10, 0.12 ] # bench lev: 0.06
reform_zeta_stys = [ "-", "-.b", "--r", ":m" ] # style
# 2.3 run & save
for tmpzeta in reform_zeta_levs
   println("Current zeta: ",tmpzeta); println("-"^20);
   # ------ data loading
   include("$(pwd())\\src\\proc_VarsDeclare.jl"); # use pwd() to avoid possible ambiguity of working directory
   include("$(pwd())\\src\\proc_InitPars.jl");
   # ------ reform
   Pt[:ζ][ ( 2020 - env.START_YEAR + 1 ) : end ] .= tmpzeta
   # ------ evaluate model & save
   eval(expr_RunModel)
   push!( res_ContriRatRise["gap/exp"], 100 .* (Dt[:LI] ./ Dt[:AggPoolExp])[idx_plot] )
   push!( res_ContriRatRise["gap/in"] , 100 .* (Dt[:LI] ./ Dt[:AggPoolIn])[idx_plot] )
   push!( res_ContriRatRise["gap/gdp"], 100 .* (Dt[:LI] ./ Dt[:Y])[idx_plot] )
   push!( res_ContriRatRise["gap/taxrev"], 100 .* (Dt[:LI] ./ (Dt[:TRc] .+ Dt[:TRw]))[idx_plot] )
end # end loops

#


text:
using DataFrames
using Crayons.Box
using NumericalDataManipulation.Types
using NumericalDataManipulation.Interpolation
using NumericalDataManipulation.Strings

abstract type Table end
struct NumericalTable <: Table
    argument::Vector{Float64}
    data::DataFrame
    functions::Vector{SplineFunction}
    function NumericalTable(data::DataFrame)
        @assert(!isempty(data), "The data must not be empty.")
        @assert(all(type -> isasupertype(type, Number), eltypes(data)), "All column types must by of the type Number.")
        @assert(size(data, 1) ≥ 3, "The data must have at least 3 rows.")
        @assert(size(data, 2) ≥ 2, "The data must have at least 2 columns.")

        # We assume that the first column is an argument grid.
        functions = Vector{SplineFunction}()
        X = data[:, 1]
        for j = 2:size(data, 2)
            Y = data[:, j]
            spline_function = SplineFunction(X, Y)
            push!(functions, spline_function)
        end
        numerical_table = new(X, data, functions)
        return numerical_table
    end
end
function Base.show(io::IO, table::NumericalTable)
    compact = get(io, :compact, false)

    X = table.argument
    data = table.data
    functions = table.functions
    N_L = size(data, 1)
    N_C = size(data, 2)

    sb = IOBuffer()
    for (i, func) ∈ enumerate(functions)
        println(sb, "    ⚫ $(LIGHT_GRAY_FG(supindex(i)))$func")
    end
    functions_description = String(take!(sb))
    if compact
        print(io, LIGHT_GRAY_FG("NumericalTable(data = X[$(X[1]), $(X[end])]($N_L), Y($N_L×$(N_C - 1)):\n$functions_description"))
    else
        print(io, LIGHT_GRAY_FG("NumericalTable(data = X[$(X[1]), $(X[end])]($N_L), Y($N_L×$(N_C - 1)):\n$functions_description"))
    end

end
function Base.show(io::IO, ::MIME"text/plain", table::NumericalTable)
    show(io, table)
end


text:
struct RunTimeData
    df::DataFrame
    loggers::Dict{String, Any}
end

struct ExperimentData
    name::String
    df::DataFrame
end


text:
# setup model
# struct to hold values
struct Inputs
    parameters::Dict{Symbol,Float64}
    initial_conditions::Dict{Symbol,Float64}
    torque_parameters::Dict{Symbol,Dict{Symbol,Float64}}
    activation_parameters::Dict{Symbol,Vector{Float64}}
    swing_data::Dict{Symbol,Vector{Float64}}
    hat_data::Dict{Symbol,Vector{Float64}}
    matching_data::Dict{Symbol,Vector{Float64}}
end

function load_inputs(;
    parameters="data/parameters.csv",
    initial_conditions="data/initial_conditions.csv",
    torque_generators="data/torque_generator_parameters.csv",
    activations="data/activation_parameters.csv",
    swing="data/matching_swing.csv",
    hat="data/HAT.csv",
    matching="data/matchingData.csv")

    # load parameters
    input_p, headers_p = readdlm(parameters, ',', header=true)
    input_p = Vector{Float64}(input_p[2, :]) # remove units row
    p = Dict(Symbol(h) => v for (h, v) in zip(headers_p, input_p))

    # load initial conditions
    input_u, headers_u = readdlm(initial_conditions, ',', header=true)
    input_u = Vector{Float64}(input_u[2, :]) # remove units row
    u = Dict(Symbol(h) => v for (h, v) in zip(headers_u, input_u))

    # load torque generator parameters
    input_tq, headers_tq = readdlm(torque_generators, ',', header=true)
    tq = Dict(
        Symbol(x[1]) => Dict(
            Symbol(h) => Float64(v) for (h, v) in zip(headers_tq[2:end], x[2:end])
        ) for x in eachrow(input_tq)
    )

    # load activation profiles
    input_act, headers_act = readdlm(activations, ',', header=true)
    α = Dict(Symbol(x[1]) => Vector{Float64}(x[2:end]) for x in eachrow(input_act))

    # load swing leg and HAT data
    input_swing, headers_swing = readdlm(swing, ',', header=true)
    input_swing = Matrix{Float64}(input_swing[2:end, :]) # remove units row
    swing_data = Dict(Symbol(h) => ts for (h, ts) in zip(headers_swing, eachcol(input_swing)))
    input_hat, headers_hat = readdlm(hat, ',', header=true)
    input_hat = Matrix{Float64}(input_hat[2:end, :]) # remove units row
    hat_data = Dict(Symbol(h) => ts for (h, ts) in zip(headers_hat, eachcol(input_hat)))

    # load matching data
    input_matching, headers_matching = readdlm(matching, ',', header=true)
    input_matching = Matrix{Float64}(input_matching[2:end, :]) # remove units row
    matching_data = Dict(Symbol(h) => ts for (h, ts) in zip(headers_matching, eachcol(input_matching)))

    # return in struct
    return Inputs(p, u, tq, α, swing_data, hat_data, matching_data)
end

# returns the parameters for a given speed from a csv file
function read_from_results(fname, speed)
    p, header = readdlm(fname, ',', Float64, header=true)
    row = p[p[:, 1].==speed, :]

    return row[3:end]
end

# returns u₀ and p created from results.csv
function load_from_results(inputs, fname, speed)
    p = read_from_results(fname, speed)

    # create new inputs
    _inputs = deepcopy(inputs)
    _inputs.initial_conditions[:vcmx] = speed
    # _inputs.activation_parameters[:, 2:end] = reshape(p, :, 6) |> permutedims # reshape to 6xn matrix, row major
    n = length(p)
    m = n ÷ 6
    [_inputs.activation_parameters[k] = p[i:i+m-1] for (k, i) in zip([:he, :ke, :ae, :hf, :kf, :af], collect(1:m:n))]

    p, u₀ = set_values(_inputs)
    return p, u₀
end

function set_values(inputs)
    @unpack parameters, initial_conditions, torque_parameters, activation_parameters, swing_data, hat_data = inputs
    parameters = NamedTuple(parameters)
    initial_conditions = NamedTuple(initial_conditions)

    # set parameters
    @inbounds (; footang, g, ina, inb, inc, ind, ine, inf, ing, inh, ini, k1, k2, k3, k4, k5, k6, k7, k8, l1, l10, l11, l2, l3, l4, l5, l6, l7, l8, l9, ma, mb, mc, md, me, mf, mg, mh, mi, mtpb, mtpk, pop1xi, pop2xi) = parameters
    mt = ma + mb + mc + md + me + mf + mg
    u8 = u9 = u10 = u11 = 0.0
    footang = deg2rad(footang)

    # set intial conditions
    @inbounds (; mtp_x, mtp_y, q3, mtp_angle, ankle_angle, knee_angle, hip_angle, vcmx, vcmy, u3, mtp_angular_velocity, ankle_angular_velocity, knee_angular_velocity, hip_angular_velocity) = initial_conditions

    # set torque parameters
    tq_p = Dict(k => [CCParameters(v), SECParameters(v)] for (k, v) in torque_parameters)

    # load activation profiles
    α_p = Dict(k => ActivationProfile(v) for (k, v) in activation_parameters)

    # set swing and hat splines
    # time = swing_data[:time]
    # hip_spl = Spline1D(time, swing_data[:hip_angle], k = 5)
    # knee_spl = Spline1D(time, swing_data[:knee_angle], k = 5)
    # ankle_spl = Spline1D(time, swing_data[:ankle_angle], k = 5)
    # mtp_spl = Spline1D(time, swing_data[:mtp_angle], k = 5)
    # # swing thigh angle
    # ea(t) = evaluate(hip_spl, t) |> deg2rad
    # eap(t) = derivative(hip_spl, t, 1) |> deg2rad
    # eapp(t) = derivative(hip_spl, t, 2) |> deg2rad
    # # swing shank angle
    # fa(t) = evaluate(knee_spl, t) |> deg2rad
    # fap(t) = derivative(knee_spl, t, 1) |> deg2rad
    # fapp(t) = derivative(knee_spl, t, 2) |> deg2rad
    # # HAT CoM location
    # hat_spl = Spline1D(hat_data[:time], hat_data[:hat_com_distance], k = 5)
    # gs(t) = evaluate(hat_spl, t)
    # gsp(t) = derivative(hat_spl, t, 1)
    # gspp(t) = derivative(hat_spl, t, 2)
    # # swing ankle angle
    # ha(t) = evaluate(ankle_spl, t) |> deg2rad
    # hap(t) = derivative(ankle_spl, t, 1) |> deg2rad
    # happ(t) = derivative(ankle_spl, t, 2) |> deg2rad
    # # swing mtp angle
    # ia(t) = evaluate(mtp_spl, t) |> deg2rad
    # iap(t) = derivative(mtp_spl, t, 1) |> deg2rad
    # iapp(t) = derivative(mtp_spl, t, 2) |> deg2rad

    # fit splines and evaluate derivatives
    time = swing_data[:time]
    hip_spl = Spline1D(time, swing_data[:hip_angle], k=5)
    knee_spl = Spline1D(time, swing_data[:knee_angle], k=5)
    ankle_spl = Spline1D(time, swing_data[:ankle_angle], k=5)
    mtp_spl = Spline1D(time, swing_data[:mtp_angle], k=5)
    hat_spl = Spline1D(hat_data[:time], hat_data[:hat_com_distance])
    # get data for 0,1,2 derivatives
    swing = [swing_data[:hip_angle] swing_data[:knee_angle] swing_data[:ankle_angle] swing_data[:mtp_angle]]
    swingp = [derivative(hip_spl, time, 1) derivative(knee_spl, time, 1) derivative(ankle_spl, time, 1) derivative(mtp_spl, time, 1)]
    swingpp = [derivative(hip_spl, time, 2) derivative(knee_spl, time, 2) derivative(ankle_spl, time, 2) derivative(mtp_spl, time, 2)]
    hat = [hat_data[:time] hat_data[:hat_com_distance]]
    hatp = [hat_data[:time] derivative(hat_spl, hat_data[:time], 1)]
    hatpp = [hat_data[:time] derivative(hat_spl, hat_data[:time], 2)]

    # interpolate 0th derivatives
    dt = time[2] - time[1]
    t = time[1]:dt:time[end] # need x values as a range
    hip_spl = CubicSplineInterpolation(t, swing[:, 1])
    knee_spl = CubicSplineInterpolation(t, swing[:, 2])
    ankle_spl = CubicSplineInterpolation(t, swing[:, 3])
    mtp_spl = CubicSplineInterpolation(t, swing[:, 4])
    hat_spl = CubicSplineInterpolation(t, hat[:, 2])
    # 1st derivatives
    hipp_spl = CubicSplineInterpolation(t, swingp[:, 1])
    kneep_spl = CubicSplineInterpolation(t, swingp[:, 2])
    anklep_spl = CubicSplineInterpolation(t, swingp[:, 3])
    mtpp_spl = CubicSplineInterpolation(t, swingp[:, 4])
    hatp_spl = CubicSplineInterpolation(t, hatp[:, 2])
    # 2nd dervatives
    hippp_spl = CubicSplineInterpolation(t, swingpp[:, 1])
    kneepp_spl = CubicSplineInterpolation(t, swingpp[:, 2])
    anklepp_spl = CubicSplineInterpolation(t, swingpp[:, 3])
    mtppp_spl = CubicSplineInterpolation(t, swingpp[:, 4])
    hatpp_spl = CubicSplineInterpolation(t, hatpp[:, 2])

    # evaluation functions
    ea(t)::Float64 = hip_spl(t) |> deg2rad
    fa(t)::Float64 = knee_spl(t) |> deg2rad
    gs(t)::Float64 = hat_spl(t) |> deg2rad
    ha(t)::Float64 = ankle_spl(t) |> deg2rad
    ia(t)::Float64 = mtp_spl(t) |> deg2rad
    eap(t)::Float64 = hipp_spl(t) |> deg2rad
    fap(t)::Float64 = kneep_spl(t) |> deg2rad
    gsp(t)::Float64 = hatp_spl(t) |> deg2rad
    hap(t)::Float64 = anklep_spl(t) |> deg2rad
    iap(t)::Float64 = mtpp_spl(t) |> deg2rad
    eapp(t)::Float64 = hippp_spl(t) |> deg2rad
    fapp(t)::Float64 = kneepp_spl(t) |> deg2rad
    gspp(t)::Float64 = hatpp_spl(t) |> deg2rad
    happ(t)::Float64 = anklepp_spl(t) |> deg2rad
    iapp(t)::Float64 = mtppp_spl(t) |> deg2rad

    # convert initial conditions into generalised coordinates and speeds
    q3 = q3 |> deg2rad              # q3 = q3
    q4 = mtp_angle |> deg2rad              # mang = q4
    q5 = (ankle_angle - 180) |> deg2rad      # aang = 180 + q5
    q6 = (180 - knee_angle) |> deg2rad      # kang = 180 - q6
    q7 = (hip_angle - 180) |> deg2rad      # hang = 180 + q7
    q1 = mtp_x + l2 * cos(q3 - q4 - q5 - q6 - q7)   # pop2x = q1 - l2*cos(q3-q4-q5-q6-q7)  (MTP contacts first)
    q2 = mtp_y + l2 * sin(q3 - q4 - q5 - q6 - q7)   # pop2y = q2 - l2*sin(q3-q4-q5-q6-q7)
    u3 = u3 |> deg2rad              # u3 = u3
    u4 = mtp_angular_velocity |> deg2rad              # mangvel = u4
    u5 = ankle_angular_velocity |> deg2rad              # aangvel = u5
    u6 = -knee_angular_velocity |> deg2rad             # kangvel = -u6
    u7 = hip_angular_velocity |> deg2rad              # hangvel = u7

    # initial angle driven values
    ea₀ = ea(0.0)
    fa₀ = fa(0.0)
    gs₀ = gs(0.0)
    ha₀ = ha(0.0)
    ia₀ = ia(0.0)
    eap₀ = eap(0.0)
    fap₀ = fap(0.0)
    gsp₀ = gsp(0.0)
    hap₀ = hap(0.0)
    iap₀ = iap(0.0)

    u1_fun(ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap) = vcmx - 0.5 * (2 * mg * gsp * cos(q3) + 2 * (l8 * mh + l8 * mi - mf * (l7 - l8)) * sin(ea - fa - q3) * (eap - fap - u3 - u8 - u9) + 2 * (l1 * ma + l2 * mb + l2 * mc + l2 * md + l2 * me + l2 * mf + l2 * mg + l2 * mh + l2 * mi) * sin(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7) - 2 * mg * gs * sin(q3) * u3 - l3 * mb * sin(footang + q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7) - 2 * (l10 * me + l10 * mf + l10 * mg + l10 * mh + l10 * mi + l9 * md) * sin(q3 - q7) * (u3 - u7) - 2 * (l10 * mf + l10 * mh + l10 * mi + me * (l10 - l9)) * sin(ea - q3) * (eap - u3 - u8) - 2 * (l7 * mc + l8 * md + l8 * me + l8 * mf + l8 * mg + l8 * mh + l8 * mi) * sin(q3 - q6 - q7) * (u3 - u6 - u7) - mh * (l3 - l5) * sin(ea + ha - footang - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) - 2 * mi * (l1 - l2) * sin(ea + ha + ia - fa - q3) * (eap + hap + iap - fap - u10 - u11 - u3 - u8 - u9) - (2 * l6 * mi - mh * (l4 - 2 * l6)) * sin(ea + ha - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) - (l4 * mb + 2 * l6 * mc + 2 * l6 * md + 2 * l6 * me + 2 * l6 * mf + 2 * l6 * mg + 2 * l6 * mh + 2 * l6 * mi) * sin(q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7)) / (ma + mb + mc + md + me + mf + mg + mh + mi)
    u2_fun(ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap) = vcmy + 0.5 * (2 * (l10 * mf + l10 * mh + l10 * mi + me * (l10 - l9)) * cos(ea - q3) * (eap - u3 - u8) + mh * (l3 - l5) * cos(ea + ha - footang - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) + (2 * l6 * mi - mh * (l4 - 2 * l6)) * cos(ea + ha - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) + 2 * mi * (l1 - l2) * cos(ea + ha + ia - fa - q3) * (eap + hap + iap - fap - u10 - u11 - u3 - u8 - u9) + 2 * (l1 * ma + l2 * mb + l2 * mc + l2 * md + l2 * me + l2 * mf + l2 * mg + l2 * mh + l2 * mi) * cos(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7) - 2 * mg * gsp * sin(q3) - 2 * mg * gs * cos(q3) * u3 - l3 * mb * cos(footang + q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7) - 2 * (l10 * me + l10 * mf + l10 * mg + l10 * mh + l10 * mi + l9 * md) * cos(q3 - q7) * (u3 - u7) - 2 * (l7 * mc + l8 * md + l8 * me + l8 * mf + l8 * mg + l8 * mh + l8 * mi) * cos(q3 - q6 - q7) * (u3 - u6 - u7) - 2 * (l8 * mh + l8 * mi - mf * (l7 - l8)) * cos(ea - fa - q3) * (eap - fap - u3 - u8 - u9) - (l4 * mb + 2 * l6 * mc + 2 * l6 * md + 2 * l6 * me + 2 * l6 * mf + 2 * l6 * mg + 2 * l6 * mh + 2 * l6 * mi) * cos(q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7)) / (ma + mb + mc + md + me + mf + mg + mh + mi)
    u1 = u1_fun(ea₀, fa₀, gs₀, ha₀, ia₀, eap₀, fap₀, gsp₀, hap₀, iap₀)
    u2 = u2_fun(ea₀, fa₀, gs₀, ha₀, ia₀, eap₀, fap₀, gsp₀, hap₀, iap₀)

    # set up torque generators
    he = TorqueGenerator(2π - deg2rad(hip_angle), deg2rad(-hip_angular_velocity), tq_p[:he][1], tq_p[:he][2], α_p[:he])
    ke = TorqueGenerator(2π - deg2rad(knee_angle), deg2rad(-knee_angular_velocity), tq_p[:ke][1], tq_p[:ke][2], α_p[:ke])
    ae = TorqueGenerator(2π - deg2rad(ankle_angle), deg2rad(-ankle_angular_velocity), tq_p[:ae][1], tq_p[:ae][2], α_p[:ae])
    hf = TorqueGenerator(deg2rad(hip_angle), deg2rad(hip_angular_velocity), tq_p[:hf][1], tq_p[:hf][2], α_p[:hf])
    kf = TorqueGenerator(deg2rad(knee_angle), deg2rad(knee_angular_velocity), tq_p[:kf][1], tq_p[:kf][2], α_p[:kf])
    af = TorqueGenerator(deg2rad(ankle_angle), deg2rad(ankle_angular_velocity), tq_p[:af][1], tq_p[:af][2], α_p[:af])

    # intial CC angles
    θcc₀ = map(x -> x.cc.θ, [he, ke, ae, hf, kf, af])

    # parameters
    p = Params(ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap, eapp, fapp, gspp, happ, iapp, ae, af, footang, g, he, hf, ina, inb, inc, ind, ine, inf, ing, inh, ini, k1, k2, k3, k4, k5, k6, k7, k8, ke, kf, l1, l10, l11, l2, l3, l4, l5, l6, l7, l8, l9, ma, mb, mc, md, me, mf, mg, mh, mi, mtpb, mtpk, pop1xi, pop2xi)

    # intial conditions
    u₀ = SVector(q1, q2, q3, q4, q5, q6, q7, u1, u2, u3, u4, u5, u6, u7, θcc₀...)

    return p, u₀

end

function setup(;
    parameters="data/parameters.csv",
    initial_conditions="data/initial_conditions.csv",
    torque_generator_parameters="data/torque_generator_parameters.csv",
    activation_parameters="data/activation_parameters.csv",
    swing="data/matching_swing.csv",
    hat="data/HAT.csv",
    matching_data="data/matchingData.csv")

    # load parameters
    input_p, headers_p = readdlm(parameters, ',', Float64, header=true)
    @inbounds footang, g, ina, inb, inc, ind, ine, inf, ing, k1, k2, k3, k4, k5, k6, k7, k8, l1, l10, l11, l12, l2, l3, l4, l5, l6, l7, l8, l9, ma, mb, mc, md, me, mf, mg, mtpb, mtpk, pop1xi, pop2xi = input_p
    mt = ma + mb + mc + md + me + mf + mg
    u8 = u9 = 0.0
    footang = deg2rad(footang)

    # load initial conditions
    input_u, headers_u = readdlm(initial_conditions, ',', Float64, header=true)
    @inbounds Xmtp, Ymtp, q3, θm, θa, θk, θh, vcmx, vcmy, u3, ωm, ωa, ωk, ωh = input_u

    # load torque generator parameters
    tq_p, headers_tq = readdlm(torque_generator_parameters, ',', header=true)
    tq_p = map(eachrow(tq_p)) do x
        x = x[2:end] .|> Float64    # remove row names, convert to float
        [CCParameters(x[1:end-1]...), SECParameters(x[end])]
    end

    # load activation profiles
    α_p, headers_α = readdlm(activation_parameters, ',', header=true)
    α_p = map(eachrow(α_p)) do x
        x = x[2:end] .|> Float64    # remove row names, convert to float
        ActivationProfile(x...)
    end

    # load swing leg and HAT data
    swing_data = readdlm(swing, ',', Float64, skipstart=1)
    hip_spl = Spline1D(swing_data[:, 1], swing_data[:, 2], k=5)
    knee_spl = Spline1D(swing_data[:, 1], swing_data[:, 3], k=5)
    # set swing and hat splines
    time = swing_data[:time]
    hip_spl = Spline1D(time, swing_data[:hip_angle], k=5)
    knee_spl = Spline1D(time, swing_data[:knee_angle], k=5)
    ankle_spl = Spline1D(time, swing_data[:ankle_angle], k=5)
    mtp_spl = Spline1D(time, swing_data[:mtp_angle], k=5)
    # swing thigh angle
    ea(t) = evaluate(hip_spl, t) |> deg2rad
    eap(t) = derivative(hip_spl, t, 1) |> deg2rad
    eapp(t) = derivative(hip_spl, t, 2) |> deg2rad
    # swing shank angle
    fa(t) = evaluate(knee_spl, t) |> deg2rad
    fap(t) = derivative(knee_spl, t, 1) |> deg2rad
    fapp(t) = derivative(knee_spl, t, 2) |> deg2rad
    # HAT CoM location
    hat_spl = Spline1D(hat_data[:time], hat_data[:hat_com_distance], k=5)
    gs(t) = evaluate(hat_spl, t)
    gsp(t) = derivative(hat_spl, t, 1)
    gspp(t) = derivative(hat_spl, t, 2)
    # swing ankle angle
    ha(t) = evaluate(ankle_spl, t) |> deg2rad
    hap(t) = derivative(ankle_spl, t, 1) |> deg2rad
    happ(t) = derivative(ankle_spl, t, 2) |> deg2rad
    # swing mtp angle
    ia(t) = evaluate(mtp_spl, t) |> deg2rad
    iap(t) = derivative(mtp_spl, t, 1) |> deg2rad
    iapp(t) = derivative(mtp_spl, t, 2) |> deg2rad

    # load matching data
    matching_data, headers_matching = readdlm(matching_data, ',', Float64, skipstart=1, header=true)
    matching_data = matching_data[1:2:end-2] # only use left side data

    # convert initial conditions into generalised coordinates and speeds
    q3 = q3 |> deg2rad              # q3 = q3
    q4 = θm |> deg2rad              # mang = q4
    q5 = (θa - 180) |> deg2rad      # aang = 180 + q5
    q6 = (180 - θk) |> deg2rad      # kang = 180 - q6
    q7 = (θh - 180) |> deg2rad      # hang = 180 + q7
    q1 = Xmtp + l2 * cos(q3 - q4 - q5 - q6 - q7)   # pop2x = q1 - l2*cos(q3-q4-q5-q6-q7)  (MTP contacts first)
    q2 = Ymtp + l2 * sin(q3 - q4 - q5 - q6 - q7)   # pop2y = q2 - l2*sin(q3-q4-q5-q6-q7)
    u3 = u3 |> deg2rad              # u3 = u3
    u4 = ωm |> deg2rad              # mangvel = u4
    u5 = ωa |> deg2rad              # aangvel = u5
    u6 = -ωk |> deg2rad             # kangvel = -u6
    u7 = ωh |> deg2rad              # hangevl = u7

    ea₀ = ea(0.0)
    fa₀ = fa(0.0)
    gs₀ = gs(0.0)
    ha₀ = ha(0.0)
    ia₀ = ia(0.0)
    eap₀ = eap(0.0)
    fap₀ = fap(0.0)
    gsp₀ = gsp(0.0)
    hap₀ = hap(0.0)
    iap₀ = iap(0.0)

    # u1 = 1.78872
    # u2 = -3.877485
    u1_fun(ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap) = vcmx - 0.5 * (2 * mg * gsp * cos(q3) + 2 * (l8 * mh + l8 * mi - mf * (l7 - l8)) * sin(ea - fa - q3) * (eap - fap - u3 - u8 - u9) + 2 * (l1 * ma + l2 * mb + l2 * mc + l2 * md + l2 * me + l2 * mf + l2 * mg + l2 * mh + l2 * mi) * sin(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7) - 2 * mg * gs * sin(q3) * u3 - l3 * mb * sin(footang + q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7) - 2 * (l10 * me + l10 * mf + l10 * mg + l10 * mh + l10 * mi + l9 * md) * sin(q3 - q7) * (u3 - u7) - 2 * (l10 * mf + l10 * mh + l10 * mi + me * (l10 - l9)) * sin(ea - q3) * (eap - u3 - u8) - 2 * (l7 * mc + l8 * md + l8 * me + l8 * mf + l8 * mg + l8 * mh + l8 * mi) * sin(q3 - q6 - q7) * (u3 - u6 - u7) - mh * (l3 - l5) * sin(ea + ha - footang - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) - 2 * mi * (l1 - l2) * sin(ea + ha + ia - fa - q3) * (eap + hap + iap - fap - u10 - u11 - u3 - u8 - u9) - (2 * l6 * mi - mh * (l4 - 2 * l6)) * sin(ea + ha - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) - (l4 * mb + 2 * l6 * mc + 2 * l6 * md + 2 * l6 * me + 2 * l6 * mf + 2 * l6 * mg + 2 * l6 * mh + 2 * l6 * mi) * sin(q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7)) / (ma + mb + mc + md + me + mf + mg + mh + mi)
    u2_fun(ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap) = vcmy + 0.5 * (2 * (l10 * mf + l10 * mh + l10 * mi + me * (l10 - l9)) * cos(ea - q3) * (eap - u3 - u8) + mh * (l3 - l5) * cos(ea + ha - footang - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) + (2 * l6 * mi - mh * (l4 - 2 * l6)) * cos(ea + ha - fa - q3) * (eap + hap - fap - u10 - u3 - u8 - u9) + 2 * mi * (l1 - l2) * cos(ea + ha + ia - fa - q3) * (eap + hap + iap - fap - u10 - u11 - u3 - u8 - u9) + 2 * (l1 * ma + l2 * mb + l2 * mc + l2 * md + l2 * me + l2 * mf + l2 * mg + l2 * mh + l2 * mi) * cos(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7) - 2 * mg * gsp * sin(q3) - 2 * mg * gs * cos(q3) * u3 - l3 * mb * cos(footang + q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7) - 2 * (l10 * me + l10 * mf + l10 * mg + l10 * mh + l10 * mi + l9 * md) * cos(q3 - q7) * (u3 - u7) - 2 * (l7 * mc + l8 * md + l8 * me + l8 * mf + l8 * mg + l8 * mh + l8 * mi) * cos(q3 - q6 - q7) * (u3 - u6 - u7) - 2 * (l8 * mh + l8 * mi - mf * (l7 - l8)) * cos(ea - fa - q3) * (eap - fap - u3 - u8 - u9) - (l4 * mb + 2 * l6 * mc + 2 * l6 * md + 2 * l6 * me + 2 * l6 * mf + 2 * l6 * mg + 2 * l6 * mh + 2 * l6 * mi) * cos(q3 - q5 - q6 - q7) * (u3 - u5 - u6 - u7)) / (ma + mb + mc + md + me + mf + mg + mh + mi)
    u1 = u1_fun(ea₀, fa₀, gs₀, ha₀, ia₀, eap₀, fap₀, gsp₀, hap₀, iap₀)
    u2 = u2_fun(ea₀, fa₀, gs₀, ha₀, ia₀, eap₀, fap₀, gsp₀, hap₀, iap₀)

    # set up torque generators
    he = TorqueGenerator(2π - deg2rad(θh), deg2rad(-ωh), tq_p[1][1], tq_p[1][2], α_p[1])
    ke = TorqueGenerator(2π - deg2rad(θk), deg2rad(-ωk), tq_p[2][1], tq_p[2][2], α_p[2])
    ae = TorqueGenerator(2π - deg2rad(θa), deg2rad(-ωa), tq_p[3][1], tq_p[3][2], α_p[3])
    hf = TorqueGenerator(deg2rad(θh), deg2rad(ωh), tq_p[4][1], tq_p[4][2], α_p[4])
    kf = TorqueGenerator(deg2rad(θk), deg2rad(ωk), tq_p[5][1], tq_p[5][2], α_p[5])
    af = TorqueGenerator(deg2rad(θa), deg2rad(ωa), tq_p[6][1], tq_p[6][2], α_p[6])
    θcc₀ = map(x -> x.cc.θ, [he, ke, ae, hf, kf, af])

    # parameters
    p = Params(ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap, eapp, fapp, gspp, happ, iapp, ae, af, footang, g, he, hf, ina, inb, inc, ind, ine, inf, ing, inh, ini, k1, k2, k3, k4, k5, k6, k7, k8, ke, kf, l1, l10, l11, l2, l3, l4, l5, l6, l7, l8, l9, ma, mb, mc, md, me, mf, mg, mh, mi, mtpb, mtpk, pop1xi, pop2xi)

    # intial conditions
    u₀ = SVector(q1, q2, q3, q4, q5, q6, q7, u1, u2, u3, u4, u5, u6, u7, θcc₀...)

    return p, u₀, matching_data
end


text:
function Base.summary(curs::Cursor)
    colnames = get_colnames(curs.cols)
    res = "Cursor with columns $colnames"
end

function get_colnames(colnames::Vector{UTF8String})
    n = length(colnames)
    res = "("
    for i in eachindex(colnames)
        i == n ? res = res * "$(colnames[i]))" :
                 res = res * "$(colnames[i]), "
    end
    return res
end

function Base.show(io::IO, curs::Cursor)
    println(io, summary(curs))
    curs.status == 100 ? print(io, "CURRENT: $(readline(curs))") :
    curs.status == 101 && println(io, "No more results to show")
    info("status: ", curs.status)
end

function Base.summary(tab::Table)
    dims = Base.dims2string(size(tab))
    colnames = get_colnames(tab.colnames)
    return "$dims Table with columns $colnames"
end

function Base.show(io::IO, tab::Table)
    println(io, summary(table))
    Base.showlimited(io, tab)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dugc_sidebar

"""
    dugc_sidebar(;kwargs...)
    dugc_sidebar(children::Any;kwargs...)
    dugc_sidebar(children_maker::Function;kwargs...)

A Sidebar component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Custom content in the "no results found" panel
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `urls` (Array; optional): URLs
- `depth` (Real; optional): depth
"""
function dugc_sidebar(; kwargs...)
        available_props = Symbol[:children, :id, :urls, :depth]
        wild_props = Symbol[]
        return Component("dugc_sidebar", "Sidebar", "dash_user_guide_components", available_props, wild_props; kwargs...)
end

dugc_sidebar(children::Any; kwargs...) = dugc_sidebar(;kwargs..., children = children)
dugc_sidebar(children_maker::Function; kwargs...) = dugc_sidebar(children_maker(); kwargs...)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
@testset "787.cheapest-flights-within-k-stops.jl" begin
    @test find_cheapest_price(3, [[1,2,100],[2,3,100],[1,3,500]], 1, 3, 1) == 200
    @test find_cheapest_price(3, [[1,2,100],[2,3,100],[1,3,500]], 1, 3, 0) == 500
end


text:
__precompile__()

module MetaICVI

# -----------------------------------------------------------------------------
# PREPROCESSING
# -----------------------------------------------------------------------------

# Load the scikitlearn definitions here for runtime loading
# using ScikitLearn: fit!, score, predict, predict_proba, @sk_import

# Load the PyObject definition at runtime because the compiled pointer will be different
# function __init__()
#     # @eval @sk_import linear_model:RidgeClassifier
#     @eval @sk_import linear_model:SGDClassifier
    # @eval using PyCallJLD
    # @eval @sk_import linear_model as lm
# end

# -----------------------------------------------------------------------------
# DEFINITIONS
# -----------------------------------------------------------------------------

# Common types and method
include("common.jl")

# Meta-icvi module definition
include("meta-icvi/meta-icvi.jl")

# -----------------------------------------------------------------------------
# EXPORTED NAMES
# -----------------------------------------------------------------------------

export

    # Structures
    MetaICVIModule,
    MetaICVIOpts,

    # Methods
    get_metaicvi,           # Compute features and "performance"
    get_features,           # Convenience function for just features
    get_training_features,  # Load data and process features for classifier training
    train_and_save

end


text:
"""
	$(TYPEDEF)
Struct to store tag id with translations from body to landmark
"""
struct TagkTl
	tagID::Symbol
	ktl::Translation
	kQl::SVector{4}#Quat #Quat does not serialize as wanted out of the box so use SVector instead
end

# ============================================================
# --------------------TagOnPose--------------------
# ============================================================
import Amphitheatre: visualize!

struct TagOnPose <: AbstractAmphitheatre
	robotId::String
	sessionId::String

	tags::Dict{Symbol,Vector{TagkTl}}

	isdrawn::Dict{Symbol, Bool}
	#tag on pose drawing propeties
	tagProp::plDrawProp
end

TagOnPose(robotId::String, sessionId::String, tagsOnPoses::Dict{Symbol,Vector{TagkTl}};
								tagProp::plDrawProp = plDrawProp(0.2, 0.1, RGBA(1,0,1,0.5))) =
			TagOnPose(robotId, sessionId, tagsOnPoses, Dict{Symbol, Bool}(), tagProp)

"""
    $(SIGNATURES)
Basic visualizer object visualize! function.
"""
function visualize!(vis::Visualizer, top::TagOnPose)::Nothing
	#TODO maybe improve this function to lower memmory allocations

	robotId = top.robotId
	sessionId = top.sessionId

	tagsOnPoses = top.tags
	isdrawn = top.isdrawn
	tagProp = top.tagProp

	for poseKey in keys(tagsOnPoses)

		if !get(isdrawn, poseKey, false)
			haskey(isdrawn, poseKey) ? isdrawn[poseKey] = true : push!(isdrawn, poseKey=>true)

			for tag = tagsOnPoses[poseKey]

				tform = tag.ktl ∘ LinearMap(UnitQuaternion(tag.kQl...))

				visPose!(vis[robotId][sessionId][:poses][poseKey][tag.tagID], tform,
							scale = tagProp.scale,
							sphereScale = tagProp.sphereScale,
							color = tagProp.color)
			end
			# sleep(0.1)
		end
	end

    return nothing
end

max_stars_repo_name: Affie/Amphitheatre

max_issues_repo_name: Affie/Amphitheatre


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1NodeAffinity <: SwaggerModel
    preferredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PreferredSchedulingTerm} } # spec name: preferredDuringSchedulingIgnoredDuringExecution
    requiredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeSelector } # spec name: requiredDuringSchedulingIgnoredDuringExecution

    function IoK8sApiCoreV1NodeAffinity(;preferredDuringSchedulingIgnoredDuringExecution=nothing, requiredDuringSchedulingIgnoredDuringExecution=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeAffinity, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        validate_property(IoK8sApiCoreV1NodeAffinity, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        o
    end
end # type IoK8sApiCoreV1NodeAffinity

const _property_map_IoK8sApiCoreV1NodeAffinity = Dict{Symbol,Symbol}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>Symbol("preferredDuringSchedulingIgnoredDuringExecution"), Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>Symbol("requiredDuringSchedulingIgnoredDuringExecution"))
const _property_types_IoK8sApiCoreV1NodeAffinity = Dict{Symbol,String}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1PreferredSchedulingTerm}", Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>"IoK8sApiCoreV1NodeSelector")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeAffinity }) = collect(keys(_property_map_IoK8sApiCoreV1NodeAffinity))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeAffinity }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1NodeAffinity[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeAffinity }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeAffinity[property_name]

function check_required(o::IoK8sApiCoreV1NodeAffinity)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NodeAffinity }, name::Symbol, val)
end


text:
#=ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (без внутренних перегородок)

РЕЗУЛЬТАТ: Робот - в исходном положении, в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядке=#
function chess(r)
    to_Sud=moves!(r,Sud)
    to_West=moves!(r,West)
    chess_parametr=to_Sud + to_West #сумма шагов до Ю.З. угла

    if iseven(chess_parametr)==true #если число чётное-угол маркеруется(исходная точка будет всегда промаркирована)
        putmarker!(r)
    end
    side = Ost
    putmarkers!(r,side)
    while isborder(r,Nord)==false
       if ismarker(r)==false
        move!(r,Nord)
        putmarker!(r)
       else
        move!(r,Nord)
       end
        side=HorizonSide(mod(Int(side)+2, 4))
        putmarkers!(r,side)
    end
    back(r,Sud)
    back(r,West)
    #принцип решения как в 3 задаче(некоторые функции изменены)

    moves!(r,Nord,to_Sud)
    moves!(r,Ost,to_West)
end
function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end
function putmarkers!(r::Robot, side::HorizonSide) #функция-ставим маркеры в линию с чередованием(через клетку)
    while isborder(r,side)==false
        check_of_marker(r,side)
    end
end
function check_of_marker(r::Robot,side::HorizonSide)#если робот на клетке с маркером-следующая без маркера(и наоборот)
    if ismarker(r)==true
       move!(r,side)
    else
       move!(r,side)
       putmarker!(r)
    end
end
function back(r::Robot,side::HorizonSide) #возвращаемся к границам сторон
    while isborder(r,side)==false
        move!(r,side)
    end
end
function moves!(r::Robot, side::HorizonSide, num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
export KozachenkoLeonenko, genentropy

"""
    KozachenkoLeonenko(; k::Int = 1, w::Int = 0) <: EntropyEstimator

Entropy estimator based on nearest neighbors. This implementation is based on Kozachenko
& Leonenko (1987)[^KozachenkoLeonenko1987],
as described in Charzyńska and Gambin (2016)[^Charzyńska2016].

`w` is the Theiler window (defaults to `0`, meaning that only the point itself is excluded
when searching for neighbours).

!!! info
    This estimator is only available for entropy estimation. Probabilities
    cannot be obtained directly.

[^Charzyńska2016]: Charzyńska, A., & Gambin, A. (2016). Improvement of the k-NN entropy estimator with applications in systems biology. Entropy, 18(1), 13.
[^KozachenkoLeonenko1987]: Kozachenko, L. F., & Leonenko, N. N. (1987). Sample estimate of the entropy of a random vector. Problemy Peredachi Informatsii, 23(2), 9-16.
"""
Base.@kwdef struct KozachenkoLeonenko <: NearestNeighborEntropyEstimator
    k::Int = 1
    w::Int = 0
end

function genentropy(x::AbstractDataset{D, T}, est::KozachenkoLeonenko; base::Real = MathConstants.e) where {D, T}
    N = length(x)
    ρs = maximum_neighbor_distances(x, est)
    h = D/N*sum(log.(base, ρs)) + log(base, ball_volume(D)) +
        MathConstants.eulergamma + log(base, N - 1)
    return h
end


text:
module PackageCompilerConvenienceTools

import PackageCompilerX
import Pkg
import Random
import UUIDs

const NATIVE_CPU_TARGET = PackageCompilerX.NATIVE_CPU_TARGET
const APP_CPU_TARGET = PackageCompilerX.APP_CPU_TARGET

const audit_app = PackageCompilerX.audit_app
const create_app = PackageCompilerX.create_app
const create_sysimage = PackageCompilerX.create_sysimage
const restore_default_sysimage = PackageCompilerX.restore_default_sysimage

include("public.jl")

include("dependencies.jl")
include("pkg_dir.jl")
include("generate_precompile_execution_file.jl")
include("utils.jl")

end # module