text:
__precompile__()

module ML

export Kernel, PolynomialKernel, GaussianKernel, ExponentialKernel,
    fit!, grad, ∇, center, KRR, KernelRidgeRegressor, RFKM, RandomFourierKernelMachine,
    KPCA, KernelPCA, negative_log_likelihood, rde, mae, rmse

using Distributions: Normal
using Distances: pairwise, colwise, Euclidean, SqEuclidean, euclidean, sqeuclidean

mae(x,y) = mean(abs(x-y))

rmse(x,y) = sqrt(sumabs2(x-y)/length(x))

abstract Kernel

type PolynomialKernel <: Kernel
    α::Float64
    c::Float64
    ξ::Float64
end

PolynomialKernel(;α=1.,c=0.,ξ=2.) = PolynomialKernel(α,c,ξ)

Base.call(k::PolynomialKernel, x::AbstractVector, y::AbstractVector) = (k.α*vecdot(x,y) + k.c).^k.ξ
Base.call(k::PolynomialKernel, x::AbstractArray, y::AbstractArray) = (k.α*(x'*y) + k.c).^k.ξ
Base.call(k::PolynomialKernel, x::AbstractArray) = k(x,x)

type GaussianKernel <: Kernel
    σ::Float64
end

GaussianKernel(;σ=1.) = GaussianKernel(σ)

Base.call(k::GaussianKernel, x::AbstractVector, y::AbstractVector) = exp(-sqeuclidean(x,y)/(2k.σ^2))
Base.call(k::GaussianKernel, x::AbstractArray, y::AbstractVector) = exp(-colwise(SqEuclidean(),x,y)/(2k.σ^2))
Base.call(k::GaussianKernel, x::AbstractVector, y::AbstractArray) = k(y,x)
Base.call(k::GaussianKernel, x::AbstractArray, y::AbstractArray) = exp(-pairwise(SqEuclidean(),x,y)/(2k.σ^2))
Base.call(k::GaussianKernel, x::AbstractArray) = exp(-pairwise(SqEuclidean(),x)/(2k.σ^2))

function grad(k::GaussianKernel, x::AbstractVector, y::AbstractVector)
    (y-x)*(k(x,y)/k.σ^2)
end

function grad(k::GaussianKernel, x::AbstractArray, y::AbstractVector)
    K = k(x,y)
    n,m = size(x,1),length(K)
    dst = Array{eltype(x)}((n,m))
    for i in 1:m
        dst[:,i] = grad(k,x[:,i],y)
    end
    reshape(dst,(n,size(K)...))
end

function grad(k::GaussianKernel, x::AbstractVector, y::AbstractArray)
    K = k(x,y)
    n,m = size(x,1),length(K)
    dst = Array{eltype(x)}((n,m))
    for i in 1:m
        dst[:,i] = grad(k,x,y[:,i])
    end
    reshape(dst,(n,size(K)...))
end

function grad(k::GaussianKernel, x::AbstractMatrix, y::AbstractMatrix)
    n,m,p = size(x,1),size(x,2),size(y,2)
    K = k(x,y)
    dst = Array{eltype(x)}((n,m,p))
    for j in 1:p
        for i in 1:m
            dst[:,i,j] = grad(k,x[:,i],y[:,j])
        end
    end
    dst
end

# symmetric
# function grad(k::GaussianKernel, x::AbstractMatrix)
#   n,m = size(x,1),size(x,2)
#   K = k(x,x)
#   dst = Array{eltype(x)}((n,m,m))
#   for j in 1:m
#       for i in j:m
#           dst[:,i,j] = (x[:,j]-x[:,i])*(K[i,j]/k.σ^2)
#           dst[:,j,i] = -dst[:,i,j]
#       end
#   end
#   dst
# end

# function grad(k::GaussianKernel, x::AbstractArray, y::AbstractArray)
#   n = size(x,1)
#   m,p = prod(size(x)[2:end]),prod(size(y)[2:end])
#   K = k(x,y)
#   shp = size(K)
#   K = reshape(K,(m,p))
#   dst = Array{eltype(x)}((n,m,p))
#   for j in 1:p
#       for i in 1:m
#           dst[:,i,j] = (y[:,j]-x[:,i])*(K[i,j]/k.σ^2)
#       end
#   end
#   reshape(dst,(n,shp...))
# end

function rand_train_test_ind(n,k)
    ind = randperm(n)
    ind[1:k], ind[k+1:end]
end

function uniform_train_test_ind(n,d)
    train = round(Int,1:n/d:n)
    test = setdiff(1:n,train)
    train,test
end

function perturb(x,i,h)
    y = copy(x)
    y[i] += h
    y
end

function fd_grad(f,x;h=1e-4)
    f0 = f(x)
    eltype(x)[(f(perturb(x,i,h))-f0)/h for i in 1:length(x)]
end

function fd_jacobian(f,x;h=1e-4)
    f0 = f(x)
    hcat([(f(perturb(x,i,h))-f0)/h for i in 1:length(x)]...)'
end

type ExponentialKernel <: Kernel
    σ::Float64
end

ExponentialKernel(;σ=1.) = ExponentialKernel(σ)

Base.call(k::ExponentialKernel, x::AbstractVector, y::AbstractVector) = exp(-euclidean(x,y)/(2k.σ^2))
Base.call(k::ExponentialKernel, x::AbstractArray, y::AbstractVector) = exp(-colwise(Euclidean(),x,y)/(2k.σ^2))
Base.call(k::ExponentialKernel, x::AbstractVector, y::AbstractArray) = k(y,x)
Base.call(k::ExponentialKernel, x::AbstractArray, y::AbstractArray) = exp(-pairwise(Euclidean(),x,y)/(2k.σ^2))
Base.call(k::ExponentialKernel, x::AbstractArray) = exp(-pairwise(Euclidean(),x)/(2k.σ^2))

abstract Predictor

# type KernelRidgeRegressor{T,Tk<:Kernel,Ni,No} <: Predictor
#     x::Array{T,Ni}          # features
#     y::Array{T,No}          # labels
#     k::Tk                   # kernel
#     λ::Float64              # regularization strength
#     w::Array{T,No}          # weights
# end

# typealias KRR KernelRidgeRegressor

# function KernelRidgeRegressor(x  ::AbstractArray,                           # features
#                               y  ::AbstractArray;                           # labels
#                               k  ::Kernel  = GaussianKernel(),              # kernel
#                               λ  ::Float64 = 1e-4)                          # regularization strength
#     KernelRidgeRegressor{eltype(x),typeof(k),ndims(x),ndims(y)}(x,y,k,λ,zeros(y))
# end

# typealias sKRR{T,Tk<:Kernel,Ni} KRR{T,Tk,Ni,1} # single valued

# function fit!(m::sKRR)
#     K = m.k(m.x)                                        # kernel matrix
#     m.w = (K + m.λ*I)\m.y                               # solve for weights
# end

# function fit!(m::KRR)
#     K = m.k(m.x)                                        # kernel matrix
#     m.w = (K + m.λ*I)\m.y'                              # solve for weights
# end

# # fit with new features, labels
# function fit!(m::KRR,
#               x::AbstractArray,                 # features
#               y::AbstractArray)                 # labels
#     m.x,m.y = x,y
#     fit!(m)
# end

# Base.call(m::sKRR, x::AbstractVector) = dot(m.w,m.k(x,m.x))
# Base.call(m::sKRR, x::AbstractArray) = m.k(x,m.x)*m.w

# Base.call(m::KRR, x::AbstractVector) = m.w' * m.k(x,m.x)
# Base.call(m::KRR, x::AbstractArray) = m.w' * m.k(x,m.x)'

# grad(m::KRR, x::AbstractVector) = grad(m.k,x,m.x)*m.w

# function grad(m::KRR, x::AbstractMatrix)
#     n,d = size(m.x)
#     k = size(x,2)
#     reshape(reshape(grad(m.k,x,m.x),(n*k,d)) * m.w, (n,k))
# end

# ∇ = grad

mae(x,y) = mean(abs(x-y))
rmse(x,y) = sqrt(sumabs2(x-y)/length(x))

function L_matrix(n::Int, d::Int)
    d == 1 && return I
    @assert n%d == 0 "n=$n must be divisble by d=$d"
    m = div(n,d)
    L = zeros(n,m)
    for i in 1:m
        L[d*(i-1)+1:d*i,i] = 1
    end
    L
end

L_matrix(x::AbstractMatrix, y::AbstractVector) = L_matrix(size(x,2),div(size(x,2),length(y)))

type KernelRidgeRegressor{T,Tk<:Kernel,Ni,No} <: Predictor
    x::Array{T,Ni}          # features
    y::Array{T,No}          # labels
    k::Tk                   # kernel
    λ::Float64              # regularization strength
    w::Array{T,No}          # weights
end

typealias KRR KernelRidgeRegressor

function KernelRidgeRegressor(x  ::AbstractArray,                           # features
                              y  ::AbstractArray;                           # labels
                              k  ::Kernel  = GaussianKernel(),              # kernel
                              λ  ::Float64 = 1e-4)                          # regularization strength
    KernelRidgeRegressor{eltype(x),typeof(k),ndims(x),ndims(y)}(x,y,k,λ,zeros(y))
end

# 1-dim features
function KernelRidgeRegressor(x  ::Vector{Float64},                         # features
                              y  ::Vector{Float64},                         # labels
                              k  ::Kernel = GaussianKernel(),               # kernel
                              λ  ::Float64 = 1e-4)                          # regularization strength
    KernelRidgeRegressor(reshape(x,(1,length(x))),y,k=k,λ=λ)
end

typealias sKRR{T,Tk<:Kernel,Ni} KRR{T,Tk,Ni,1} # single valued

function fit!(m::sKRR)
    # m.w = (K + m.λ*I)\m.y                             # solve for weights
    K = m.k(m.x)                                        # kernel matrix
    L = L_matrix(m.x,m.y)
    KL = K*L
    Λ = I/m.λ
    KLΛ = KL*Λ
    m.w = (K + KLΛ*KL')\(KLΛ*m.y)
end

# fit with new features, labels
function fit!(m::KRR,
              x::AbstractArray,                 # features
              y::AbstractArray)                 # labels
    m.x,m.y = x,y
    fit!(m)
end

Base.call(m::sKRR, x::AbstractVector) = dot(m.w,m.k(x,m.x))
Base.call(m::sKRR, x::AbstractArray) = m.k(x,m.x)*m.w

Base.call(m::KRR, x::AbstractVector) = m.w' * m.k(x,m.x)
Base.call(m::KRR, x::AbstractArray) = m.w' * m.k(x,m.x)'

grad(m::KRR, x::AbstractVector) = grad(m.k,x,m.x)*m.w

function grad(m::KRR, x::AbstractMatrix)
    n,d = size(m.x)
    k = size(x,2)
    reshape(reshape(grad(m.k,x,m.x),(n*k,d)) * m.w, (n,k))
end

∇ = grad

# # uncentered kernel matrix
# K = kpca._get_kernel(X_train)

# # scaled kPCA eigenvectors
# A = (kpca.alphas_/sqrt(kpca.lambdas_)).T

# # centering matrix
# H = empty((Nt,Nt))
# H.fill(-1.0/Nt)
# for i in range(Nt):
#   H[i,i] += 1

# # projection matrix
# Q = dot(H, dot(dot(A.T, A), H))

# def coeff(x, kernel=None, X=None, Q=None, K=None):
#   n = K.shape[0]
#   kx = squeeze(kernel(x, X))

#   # expansion coefficients of sample over training
#   # samples in uncentered feature space
#   c = dot(Q, kx - sum(K, axis=0)/n) + 1.0/n

#   return c

# def preimage(x0, c, kernel=None, X=None, max_steps=100, tol=1E-15):
#   x = x0.copy()

#   # fixed point iteration
#   for _ in range(max_steps):
#       kx = squeeze(kernel(x, X))
#       dx = dot(c*kx, X) / dot(c, kx) - x

#       x = x + dx

#       if dot(dx, dx) < tol:
#           break

#   return x

# abstract Predictor

# type KernelRidgeRegressor{T<:Kernel} <: Predictor
#   x::Matrix{Float64}      # features
#   y::Vector{Float64}      # labels
#   k::T                    # kernel
#   λ::Float64              # regularization strength
#   β::Float64              # gradient regularization strength
#   w::Vector{Float64}      # weights
# end

# typealias KRR KernelRidgeRegressor

# function KernelRidgeRegressor(x  ::Matrix{Float64},                       # features
#                             y  ::Vector{Float64};                         # labels
#                             k  ::Kernel  = GaussianKernel(),              # kernel
#                             λ  ::Float64 = 1e-4,                          # regularization strength
#                             β  ::Float64 = 1e-4)                          # gradient regularization strength
#   KernelRidgeRegressor{typeof(k)}(x,y,k,λ,β,[])
# end

# # 1-dim features
# # function KernelRidgeRegressor(x  ::Vector{Float64},                         # features
# #                               y  ::Vector{Float64},                         # labels
# #                               # dy ::Vector{Float64};                       # gradients
# #                               k  ::Kernel = GaussianKernel(),               # kernel
# #                               λ  ::Float64 = 1e-4,                          # regularization strength
# #                               β  ::Float64 = 1e-2)                          # gradient regularization strength
# #     KernelRidgeRegressor(reshape(x,(1,length(x))),y,reshape(dy,(1,length(dy))),k=k,λ=λ,β=β)
# # end

# # function loss(m::KernelRidgeRegressor)
# #     n,d = size(m.X)
# #     K = m.k(m.X)
# #     # δK = reshape(grad(m.k,m.X,m.X),(n*d,d))
# #     function loss(w)
# #         Kw = K*w
# #         sumabs2(Kw-m.y) + m.λ*dot(w,Kw) #+ m.β*sumabs2(δK*w - vec(m.dy))
# #         # m.β*sumabs(δK*w - vec(m.dy))
# #     end
# #     loss
# # end

# function fit!(m::KernelRidgeRegressor)
#   K = m.k(m.x)                                        # kernel matrix
#   m.w = (K + m.λ*I)\m.y                               # solve for weights
# end

# # function fit!(m::KernelRidgeRegressor;
# #               method     ::Symbol  = :l_bfgs,
# #               xtol       ::Float64 = 1e-8,
# #               ftol       ::Float64 = 1e-8,
# #               grtol      ::Float64 = 1e-8,
# #               iterations ::Int     = 1000)
# #     res = optimize(loss(m),zeros(size(m.X,2)),method=method,
# #         xtol=xtol,ftol=ftol,grtol=grtol,iterations=iterations)
# #     m.w = res.minimum
# #     res
# # end

# # fit with new features, labels
# function fit!(m::KernelRidgeRegressor,
#             x::Matrix{Float64},                       # features
#             y::Vector{Float64})                       # labels
#   m.x,m.y = x,y
#   fit!(m)
# end

# # function fit!(m::KernelRidgeRegressor,
# #               x::Vector{Float64},                       # features
# #               y::Vector{Float64})                       # labels
# #     fit!(m,reshape(x,(1,length(x))),y)
# # end

# Base.call(m::KernelRidgeRegressor, x::AbstractVector) = dot(m.w,m.k(x,m.x))
# Base.call(m::KernelRidgeRegressor, x::AbstractMatrix) = m.k(x,m.x)*m.w
# # Base.call(m::KernelRidgeRegressor, x::Number) = m([x])

# grad(m::KernelRidgeRegressor, x::AbstractVector) = grad(m.k,x,m.x)*m.w

# function grad(m::KRR, x::AbstractMatrix)
#   n,d = size(m.x)
#   k = size(x,2)
#   reshape(reshape(grad(m.k,x,m.x),(n*k,d)) * m.w, (n,k))
# end

# abstract Predictor

# type KernelRidgeRegressor{T<:Kernel} <: Predictor
#   X::Matrix{Float64}      # features
#   y::Vector{Float64}      # labels
#   k::T                    # kernel
#   λ::Float64              # regularization strength
#   w::Vector{Float64}      # weights
# end

# typealias KRR KernelRidgeRegressor

# function KernelRidgeRegressor(X ::Matrix{Float64},                            # features
#                             y ::Vector{Float64};                          # labels
#                             k ::Kernel  = GaussianKernel(),               # kernel
#                             λ ::Float64 = 1e-4)                           # regularization strength
#   KernelRidgeRegressor{typeof(k)}(X,y,k,λ,[])
# end

# # 1-dim features
# function KernelRidgeRegressor(x ::Vector{Float64},                            # features
#                             y ::Vector{Float64};                          # labels
#                             k ::Kernel = GaussianKernel(),                # kernel
#                             λ ::Float64 = 1e-4)                           # regularization strength
#   KernelRidgeRegressor(reshape(x,(1,length(x))),y,k=k,λ=λ)
# end

# function loss(m::KernelRidgeRegressor)
#   K = m.k(m.X)
#   function loss(w)
#       Kw = K*w
#       sumabs2(Kw-m.y) + m.λ*dot(w,Kw)
#   end
#   loss
# end

# function fit!(m::KernelRidgeRegressor;
#             method     ::Symbol  = :l_bfgs,
#             xtol       ::Float64 = 1e-8,
#             ftol       ::Float64 = 1e-8,
#             grtol      ::Float64 = 1e-8,
#             iterations ::Int     = 1000)
#   res = optimize(loss(m),zeros(size(m.X,2)),method=method,
#       xtol=xtol,ftol=ftol,grtol=grtol,iterations=iterations)
#   m.w = res.minimum
#   res
# end

# # function fit!(m::KernelRidgeRegressor)
# #     K = m.k(m.X)                                        # kernel matrix
# #     K[diagind(K)] += m.λ                                # regularization
# #     m.w = K\m.y                                         # solve for weights
# # end

# # fit with new features, labels
# function fit!(m::KernelRidgeRegressor,
#             X::Matrix{Float64},                       # features
#             y::Vector{Float64})                       # labels
#   m.X,m.y = X,y
#   fit!(m)
# end

# function fit!(m::KernelRidgeRegressor,
#             x::Vector{Float64},                       # features
#             y::Vector{Float64})                       # labels
#   fit!(m,reshape(x,(1,length(x))),y)
# end

# Base.call(m::KernelRidgeRegressor, x::AbstractVector) = dot(m.w,m.k(x,m.X))
# Base.call(m::KernelRidgeRegressor, x::AbstractMatrix) = m.k(x,m.X)*m.w
# Base.call(m::KernelRidgeRegressor, x::Number) = m([x])

# grad(m::KernelRidgeRegressor, x::AbstractVector) = grad(m.k,x,m.X)*m.w

# function grad(m::KRR, x::AbstractMatrix)
#   n,d = size(m.X)
#   k = size(x,2)
#   reshape(reshape(grad(m.k,x,m.X),(n*k,d)) * m.w, (n,k))
# end

type RandomFourierKernelMachine <: Predictor
    X::Matrix{Float64}              # features
    y::Vector{Float64}              # labels
    φ::Matrix{Float64}              # random features
    n::Int                          # input dimensionality
    d::Int                          # number of random features
    σ::Float64                      # kernel width
    λ::Float64                      # regularization strength
    w::Vector{Float64}              # weights
end

typealias RFKM RandomFourierKernelMachine

function RandomFourierKernelMachine(X ::Matrix{Float64},            # features
                                    y ::Vector{Float64};            # labels
                                    d ::Int     = 100,              # number of random features
                                    σ ::Float64 = 1.,               # kernel width
                                    λ ::Float64 = 1e-4)             # regularization strength
  RandomFourierKernelMachine(X,y,Matrix{Float64}(),size(X,1),d,σ,λ,[])
end

function kernel_matrix(m::RFKM, x::AbstractArray)
  z = exp(m.φ*x*im)/sqrt(m.d)
  real(z'*z)
end

function kernel_matrix(m::RFKM, x::AbstractArray, y::AbstractArray)
  zx = exp(m.φ*x*im)/sqrt(m.d)
  zy = exp(m.φ*y*im)/sqrt(m.d)
  real(zx'*zy)
end

function fit!(m::RFKM)
    d,n = m.d,m.n
    dist = Normal(0,1/m.σ)                      # normal distribution
    m.φ = rand(dist, (d,n))                     # random features
    z = cos(m.φ*m.X)/sqrt(m.d)                  # projection of X onto random features
    # z0 = centered ? zeros(n) : squeeze(mean(z,2),2)
    # y0 = centered ? zeros(n) : squeeze(mean(y,2),2)
    # zc = z .- z0
    # yc = y .- y0
    m.w = (z*z' + m.λ*I)\(z*m.y)
end

function Base.call(m::RFKM, x::AbstractVector)
    z = cos(m.φ*x)/sqrt(m.d)
    # y = Φ.α'*(z .- Φ.z0) .+ Φ.y0
    dot(m.w,z)
end

function Base.call(m::RFKM, x::AbstractMatrix)
    z = cos(m.φ*x)/sqrt(m.d)
    z'*m.w
end

function grad(m::RFKM, x::AbstractVector)
  dz = m.φ*im .* exp(m.φ*x*im)/sqrt(m.d)
  real(dz'*m.w)
end

function grad(m::RFKM, x::AbstractArray)
  n,d,k = m.n,m.d,size(x,2)
  c = -sin(m.φ*x)/sqrt(d)
  dz = m.φ .* reshape(c,(d,1,k))
  reshape(m.w' * reshape(dz,(d,n*k)), (n,k))
end

# H = I + fill(-1/n,size(K))
# H*K*H

doc"""
centering kernel matrices

$\varphi = \frac{1}{n} \sum_i \varphi(x_i)$

$\tilde K_{ij} = (\varphi(x_i) - \tilde\varphi)^T (\varphi(x_j) - \tilde\varphi)$

$\tilde K_{ij} = K_{ij} - \frac{1}{n} \sum_j K_{ij} - \frac{1}{n} \sum_i K_{ji} + \frac{1}{n^2} \sum_{ij} K_{ij}$

$\tilde L_{ij} = L_{ij} - \frac{1}{n} \sum_j K_{ij} - \frac{1}{n} \sum_i L_{ji} + \frac{1}{n^2} \sum_{ij} K_{ij}$

$\tilde M_{ij} = M_{ij} - \frac{1}{n} \sum_j L_{ij} - \frac{1}{n} \sum_i L_{ji} + \frac{1}{n^2} \sum_{ij} M_{ij}$

"""
function center(K)
    n = size(K,1)
    A = sum(K,1)/n
    K .- A .- A' + sum(K)/n^2
end

# J = fill(-1/n,(m,n))
# (L + J*K)*H
function center(K,L::Matrix)
    n = size(K,1)
    L .- sum(K,1)/n .- sum(L,2)/n + sum(K)/n^2
end

function center(K,L::Vector)
    n = size(K,1)
    L .- vec(sum(K,1))/n .- sum(L)/n + sum(K)/n^2
end

# M + L*J' + J*L' + J*K*J'
function center(K,L,M)
    n = size(K,1)
    M .- sum(L,2)/n .- sum(L,2)'/n + sum(K)/n^2
end

type KernelPCA{T<:Kernel}
    k::T                        # kernel
    x::Matrix{Float64}          # features
    K::Matrix{Float64}          # kernel matrix
    u::Vector{Float64}          # eigenvalues
    v::Matrix{Float64}          # eigenvectors
    n::Int                      # number of principal components
end

typealias KPCA KernelPCA

function KernelPCA(; k = GaussianKernel(), n::Int = -1)
    KernelPCA{typeof(k)}(k,Matrix{Float64}(),Matrix{Float64}(),[],Matrix{Float64}(),n)
end

function fit!(m::KPCA, x::Matrix{Float64})
    n = m.n
    m.x = x                                                 # features
    m.K = m.k(x)                                            # compute kernel matrix
    K = center(m.K)                                         # center kernel matrix
    u,v = eig(Symmetric(K))                                 # compute principal components
    ind = sortperm(u,rev=true)                              # sort eigenvalues
    u = u[ind]
    v = v[:,ind]
    n = n == -1 ? findlast(u.>0,true) : n                   # ignore zero eigenvalues
    m.u = u[1:n]                                            # select n eigenvalues
    m.v = v[:,1:n]                                          # select n eigenvectors
    m.v ./= sqrt(m.u)'                                      # normalize eigenvectors
    return
end

function Base.call(m::KPCA, x::Vector{Float64})
    kx = center(m.K,m.k(x,m.x))         # compute centered kernel matrix
    m.v' * kx                           # project features onto principal components
end

function Base.call(m::KPCA, x::Matrix{Float64})
    K = center(m.K,m.k(x,m.x))          # compute centered kernel matrix
    (K * m.v)'                          # project features onto principal components
end

function grad(m::KPCA, x::Vector{Float64})
    n,d = size(m.x)
    δK = grad(m.k,x,m.x)
    δK = δK .- sum(δK,2)/d
    δK * m.v
end

function grad(m::KPCA, x::Matrix{Float64})
    n,d = size(m.x)
    l = size(x,2)
    δK = grad(m.k,x,m.x)
    δK = δK .- sum(δK,3)/d
    reshape(reshape(δK,(n*l,d)) * m.v, (n,l,m.n))
end

∇ = grad

function negative_log_likelihood(φ::KPCA, y::AbstractVector)
    n = length(φ.u)
    z = φ.v' * y  # projection of labels onto KPCA eigenvectors
    σ1 = d -> sum(z[1:d].^2) / d
    σ2 = d -> sum(z[d+1:n].^2) / (n-d)
    Float64[d/n * log(σ1(d)) + (n-d)/n * log(σ2(d)) for d in 1:n-1]
end

function rde(φ::KPCA, y::AbstractVector)
    nll = negative_log_likelihood(φ,y)
    d = findfirst(diff(nll).<0)
end

end # module

# abstract Predictor

# abstract Loss

# type L2RegularizedQuadraticLoss <: Loss
#   λ::Float64
# end

# L2RegularizedQuadraticLoss(;λ=1e-8) = L2RegularizedQuadraticLoss(λ)

# type QuadraticLoss <: Loss end

# type L1RegularizedQuadraticLoss <: Loss
#   λ::Float64
# end

# L1RegularizedQuadraticLoss(;λ=1e-8) = L1RegularizedQuadraticLoss(λ)

# type KernelRidgeRegressor{Tk<:Kernel,Tl<:Loss} <: Predictor
#   k::Tk                   # kernel
#   l::Tl                   # loss model
#   X::Matrix{Float64}      # features
#   y::Vector{Float64}      # labels
#   w::Vector{Float64}      # weights
# end

# function KernelRidgeRegressor(X ::Matrix{Float64},                            # features
#                             y ::Vector{Float64};                          # labels
#                             k ::Kernel = GaussianKernel(),                # kernel
#                             l ::Loss   = L2RegularizedQuadraticLoss())    # loss
#   KernelRidgeRegressor{typeof(k),typeof(l)}(k,l,X,y,[])
# end

# # 1-dim features
# function KernelRidgeRegressor(x ::Vector{Float64},                            # features
#                             y ::Vector{Float64};                          # labels
#                             k ::Kernel = GaussianKernel(),                # kernel
#                             l ::Loss   = L2RegularizedQuadraticLoss())    # loss
#   KernelRidgeRegressor(reshape(x,(1,length(x))),y,k=k,l=l)
# end

# function loss{Tk<:Kernel}(m::KernelRidgeRegressor{Tk,L1RegularizedQuadraticLoss})
#   K = m.k(m.X)                                # kernel matrix
#   L = chol(K)                                 # cholesky factorization
#   function loss(w)
#       sumabs2(K*w-m.y) + m.l.λ*sumabs(L*w)
#   end
#   loss
# end

# function loss{Tk<:Kernel}(m::KernelRidgeRegressor{Tk,L2RegularizedQuadraticLoss})
#   K = m.k(m.X)
#   function loss(w)
#       Kw = K*w
#       sumabs2(Kw-m.y) + m.l.λ*dot(w,Kw)
#   end
#   loss
# end

# function loss_gradient!{Tk<:Kernel}(m::KernelRidgeRegressor{Tk,L2RegularizedQuadraticLoss})
#   K = m.k(m.X)
#   function loss_gradient!(w, dst)
#       Kw = K*w
#       dst[:] = 2*K*(Kw-m.y) + 2*m.l.λ*Kw
#   end
#   loss_gradient!
# end

# function loss{Tk<:Kernel}(m::KernelRidgeRegressor{Tk,QuadraticLoss})
#   K = m.k(m.X)
#   loss(w) = sumabs2(K*w-m.y)
#   loss
# end

# function loss_gradient!{Tk<:Kernel}(m::KernelRidgeRegressor{Tk,QuadraticLoss})
#   K = m.k(m.X)
#   function loss_gradient!(w, dst)
#       dst[:] = 2*K*(K*w-m.y)
#   end
#   loss_gradient!
# end

# @generated function fit!(m::KernelRidgeRegressor;
#                        method     ::Symbol  = :l_bfgs,
#                        xtol       ::Float64 = 1e-8,
#                        ftol       ::Float64 = 1e-8,
#                        grtol      ::Float64 = 1e-8,
#                        iterations ::Int     = 1000)
#   if method_exists(:loss_gradient!,(m,))  # if loss gradient is implemented
#       return quote
#           res = optimize(loss(m),loss_gradient!(m),zeros(size(m.X,2)),method=method,
#               xtol=xtol,ftol=ftol,grtol=grtol,iterations=iterations)
#           m.w = res.minimum
#           res
#       end
#   else
#       return quote
#           res = optimize(loss(m),zeros(size(m.X,2)),method=method,
#               xtol=xtol,ftol=ftol,grtol=grtol,iterations=iterations)
#           m.w = res.minimum
#           res
#       end
#   end
# end

# function fit!{Tk<:Kernel}(m::KernelRidgeRegressor{Tk,L2RegularizedQuadraticLoss})
#   K = m.k(m.X)                                        # kernel matrix
#   K[diagind(K)] += m.l.λ                              # regularization
#   m.w = K\m.y                                         # solve for weights
# end

# # fit with new features, labels
# function fit!(m::KernelRidgeRegressor,
#             X::Matrix{Float64},                       # features
#             y::Vector{Float64})                       # labels
#   m.X,m.y = X,y
#   fit!(m)
# end

# function fit!(m::KernelRidgeRegressor,
#             x::Vector{Float64},                       # features
#             y::Vector{Float64})                       # labels
#   fit!(m,reshape(x,(1,length(x))),y)
# end

# predict
# Base.call(m::KernelRidgeRegressor, x::AbstractVector) = dot(m.w,m.k(x,m.X))
# Base.call(m::KernelRidgeRegressor, x::AbstractMatrix) = m.k(x,m.X)*m.w
# Base.call(m::KernelRidgeRegressor, x::Number) = m([x])


text:
const epoch_dt_format3 = Dates.DateFormat("yyyy mm dd HH MM SS.sssssss")

function parse_3epoch_time(ephs::String,toUTC::Millisecond,
    epoch_dt_format3::DateFormat=epoch_dt_format3)

    dtstr = SubString(ephs, 3, 29)
    y = parse_Int_or_Empty(dtstr[1:4])
    m = parse_Int_or_Empty(dtstr[6:7])
    d = parse_Int_or_Empty(dtstr[9:10])
    H = parse_Int_or_Empty(dtstr[12:13])
    M = parse_Int_or_Empty(dtstr[15:16])
    S = parse_Int_or_Empty(dtstr[18:19])
    s = parse_Int_or_Empty(dtstr[21:23])
    this_dt = Dates.DateTime(y, m, d, H, M, S, s)
    # @show ephs
    # @show this_dt
    this_dt += toUTC
    this_dt
end

function add_epoch_line3!(sods::Dict{SatID,SatelliteObsData},
                line::String, this_dt::DateTime,
                systemNumObs::D_SystemNumObs,
                recommend_length::Int64, )

    satid = SatID(SubString(line, 1, 3))
    obsNum, obsIds = get(systemNumObs, satid.ss, (0, Array{ObsID,1}()))

    obsNum == 0 && throw(error("no such system $satid"))
    line = rpad(line, 3 + 16 * obsNum)

    if haskey(sods, satid)
        this_sod = sods[satid]
        this_dts = this_sod.dts
        this_obs = this_sod.obs
    else
        this_dts  = Array{DateTime,1}()
        sizehint!(this_dts, recommend_length)
        this_obs = Array{Array{RinexDatum,1},1}(undef, obsNum)
        for ion in 1:obsNum
            this_obs[ion] = Array{RinexDatum,1}()
            sizehint!(this_obs[ion], recommend_length)
        end
        sods[satid] = SatelliteObsData(satid, obsIds,
                                        this_dts, this_obs)
    end

    push!(this_dts, this_dt)
    si = 4
    for ion in 1:obsNum
        push!(this_obs[ion],
                RinexDatum(SubString(line, si, si + 15)), )
        si += 16
    end

end

function add_this_epoch_comments3!(comments::Dict{DateTime,String},
                dt_before::DateTime, f::IOStream, numSVs::Int64, )
    this_comments = Array{String,1}(undef, numSVs)
    for ii_el in 1:numSVs
        this_comments[ii_el]  = readline(f)
    end
    if haskey(comments, dt_before)
        comments[dt_before] = string(comments[dt_before],"\n",
                                join(this_comments, "\n"))
    else
        comments[dt_before] = join(this_comments, "\n")
    end
end

function read_one_epoch3(f::IOStream,
    epochs::Dict{DateTime,EpochHeader},
    comments::Dict{DateTime,String},
    add_this_epoch_line3::Function,
    toUTC::Millisecond,
    ephs_dr_set::Set{Int64},
    dt_before::DateTime,)::DateTime

    ephs = readline(f)

    if ephs[1] != '>'
        throw(error("epoch doesn't start with > ! \n$ephs"))
    end

    epf = parse(Int64, ephs[32])
    numSVs = parse(Int64, ephs[33:35])

    if epf in ephs_dr_set # 0,1,6 data record
        this_dt = parse_3epoch_time(ephs, toUTC)
        if haskey(epochs, this_dt)
            @warn "duplicate epoch? $ephs"
            delete!(epochs, this_dt)
        end
        len_ephs = length(ephs)
        if len_ephs > 41
            index_end = len_ephs > 55 ? 55 : len_ephs
            rco = parse(Float64, ephs[42:index_end])
        else
            rco = 0.0
        end
        epochs[this_dt] = EpochHeader(this_dt, epf, numSVs, rco)

        for ii_el in 1:numSVs
            line = readline(f)
            add_this_epoch_line3(line, this_dt)
        end
    else
        add_this_epoch_comments3!(comments, dt_before, f, numSVs)
        this_dt = dt_before
    end
    this_dt
end

function read_3obsdata(f::IOStream, sysNumObs::D_SystemNumObs,
    toUTC::Millisecond=Millisecond(0),
    recommend_length::Int64=4800,
    ephs_dr_set::Set{Int64}=Set([0, 1, 6]))

    epochs  = Dict{DateTime,EpochHeader}()
    comments = Dict{DateTime,String}()

    sods = Dict{SatID,SatelliteObsData}()

    add_this_epoch_line3(line::String,
        this_dt::DateTime) = add_epoch_line3!(
                                sods,
                                line, this_dt,
                                sysNumObs, recommend_length, )

    dt_before = DateTime(0, 1, 1, 0, 0, 0)
    while !eof(f)
        dt_before = read_one_epoch3(f,
                epochs, comments,
                add_this_epoch_line3,
                toUTC,
                ephs_dr_set, dt_before)
    end

    RinexObsData(epochs, comments, sods)
end

ext: jl


text:
# NullableArray is dense and allows fast linear indexing.
import Base: LinearFast

@compat Base.IndexStyle(::Type{<:NullableArray}) = IndexLinear()

# resolve ambiguity created by the two definitions that follow.
function Base.getindex{T, N}(X::NullableArray{T, N})
    return X[1]
end

"""
    getindex{T, N}(X::NullableArray{T, N}, I::Int...)

Retrieve a single entry from a `NullableArray`. If the value in the entry
designated by `I` is present, then it will be returned wrapped in a
`Nullable{T}` container. If the value is missing, then this method returns
`Nullable{T}()`.
"""
# Extract a scalar element from a `NullableArray`.
@inline function Base.getindex{T, N}(X::NullableArray{T, N}, I::Int...)
    if isbits(T)
        ifelse(X.isnull[I...], Nullable{T}(), Nullable{T}(X.values[I...]))
    else
        if X.isnull[I...]
            Nullable{T}()
        else
            Nullable{T}(X.values[I...])
        end
    end
end

"""
    getindex{T, N}(X::NullableArray{T, N}, I::Nullable{Int}...)

Just as above, with the additional behavior that this method throws an error if
any component of the index `I` is null.
"""
@inline function Base.getindex{T, N}(X::NullableArray{T, N},
                                     I::Nullable{Int}...)
    any(isnull, I) && throw(NullException())
    values = [ get(i) for i in I ]
    return getindex(X, values...)
end

"""
    setindex!(X::NullableArray, v::Nullable, I::Int...)

Set the entry of `X` at position `I` equal to a `Nullable` value `v`. If
`v` is null, then only `X.isnull` is updated to indicate that the entry at
index `I` is null. If `v` is not null, then `X.isnull` is updated to indicate
that the entry at index `I` is present and `X.values` is updated to store the
value wrapped in `v`.
"""
# Insert a scalar element from a `NullableArray` from a `Nullable` value.
@inline function Base.setindex!(X::NullableArray, v::Nullable, I::Int...)
    if isnull(v)
        X.isnull[I...] = true
    else
        X.isnull[I...] = false
        X.values[I...] = get(v)
    end
    return v
end

"""
    setindex!(X::NullableArray, v::Any, I::Int...)

Set the entry of `X` at position `I` equal to `v`. This method always updates
`X.isnull` to indicate that the entry at index `I` is present and `X.values`
to store `v` at `I`.
"""
# Insert a scalar element from a `NullableArray` from a non-Nullable value.
@inline function Base.setindex!(X::NullableArray, v::Any, I::Int...)
    X.values[I...] = v
    X.isnull[I...] = false
    return v
end

function unsafe_getindex_notnull(X::NullableArray, I::Int...)
    return Nullable(getindex(X.values, I...))
end

function unsafe_getvalue_notnull(X::NullableArray, I::Int...)
    return getindex(X.values, I...)
end

if VERSION >= v"0.5.0-dev+4697"
    function Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, i::Nullable)
        isnull(i) ? throw(NullException()) : checkindex(Bool, inds, get(i))
    end

    function Base.checkindex{N}(::Type{Bool}, inds::AbstractUnitRange, I::NullableArray{Bool, N})
        any(isnull, I) && throw(NullException())
        checkindex(Bool, inds, I.values)
    end

    function Base.checkindex{T<:Real}(::Type{Bool}, inds::AbstractUnitRange, I::NullableArray{T})
        any(isnull, I) && throw(NullException())
        b = true
        for i in 1:length(I)
            @inbounds v = unsafe_getvalue_notnull(I, i)
            b &= checkindex(Bool, inds, v)
        end
        return b
    end
else
    function Base.checkbounds{T<:Real}(::Type{Bool}, sz::Int, x::Nullable{T})
        isnull(x) ? throw(NullException()) : checkbounds(Bool, sz, get(x))
     end

    function Base.checkbounds(::Type{Bool}, sz::Int, I::NullableVector{Bool})
         any(isnull, I) && throw(NullException())
        length(I) == sz
     end

    function Base.checkbounds{T<:Real}(::Type{Bool}, sz::Int, I::NullableArray{T})
        inbounds = true
         any(isnull, I) && throw(NullException())
         for i in 1:length(I)
             @inbounds v = unsafe_getvalue_notnull(I, i)
            inbounds &= checkbounds(Bool, sz, v)
         end
        return inbounds
     end
end

function Base.to_index(X::NullableArray)
    any(isnull, X) && throw(NullException())
    Base.to_index(X.values)
end

"""
    nullify!(X::NullableArray, I...)

This is a convenience method to set the entry of `X` at index `I` to be null
"""
function nullify!(X::NullableArray, I...)
    setindex!(X, Nullable{eltype(X)}(), I...)
end


text:
using TDAC
using TimerOutputs
using MPI

MPI.Init()
mpi_size = MPI.Comm_size(MPI.COMM_WORLD)

# Warmup run
warmup_custom_params = Dict(
    "nobs"=>4,
    "padding"=>0,
    "n_time_step"=>1,
    "n_integration_step"=>1,
    "time_step"=>0.1,
    "nprt"=>mpi_size)
warmup_params = TDAC.get_params(warmup_custom_params)
TDAC.tdac(warmup_params)

# Real run
TimerOutputs.enable_debug_timings(TDAC)
custom_params = Dict(
    "nobs"=>64,
    "verbose"=>true,
    "enable_timers"=>true,
    "padding"=>0,
    "output_filename"=>string("weak_scaling_r",mpi_size,".h5"),
    "nprt"=>mpi_size * 64)
params = TDAC.get_params(custom_params)
TDAC.tdac(params)


text:
"""
    AutoCore
Defines the core Automotive Types
"""
module AutoCore

using Compat
using Vec

export

    # curves
    CurvePt,
    Curve,
    CurveIndex,
    CurveProjection,

    CURVEINDEX_START,

    get_lerp_time,
    index_closest_to_point,
    get_curve_index,
    curveindex_end,

    # roadway
    LaneTag,
    LaneBoundary,
    Lane,
    LaneConnection,
    SpeedLimit,
    RoadSegment,
    Roadway,
    RoadProjection,
    RoadIndex,

    NULL_BOUNDARY,
    NULL_LANETAG,
    NULL_ROADINDEX,

    DEFAULT_SPEED_LIMIT,
    DEFAULT_LANE_WIDTH,

    is_in_exits,
    is_in_entrances,
    is_at_curve_end,
    is_between_segments_lo,
    is_between_segments_hi,
    is_between_segments,
    has_segment,
    has_lanetag,
    has_next,
    has_prev,
    next_lane,
    prev_lane,
    next_lane_point,
    prev_lane_point,
    connect!,
    prev,
    project_to_closest_lane,
    move_along,
    n_lanes_left,
    n_lanes_right,
    get_neighbor_lanetag_left,
    get_neighbor_lanetag_right,
    read_dxf,

    # vehicles
    Frenet,
    VehicleState,
    Vehicle,
    AgentClass,
    VehicleDef,

    DEFAULT_VEHICLE_DEF,
    NULL_VEHICLEDEF,
    NULL_FRENET,

    get_posG,
    get_vel_s,
    get_vel_t,
    get_center,
    get_footpoint,
    get_front_center,
    get_rear_center,

    # trajdata
    TrajdataFrame,
    TrajdataState,
    Trajdata,
    TrajdataVehicleIterator,

    get_ids,
    nframes,
    frame_inbounds,
    carsinframe,
    nth_carid,
    iscarinframe,
    get_first_frame_with_id,
    get_last_frame_with_id,
    get_time,
    get_elapsed_time,
    get_mean_timestep,
    get_vehiclestate,
    get_vehicledef,
    get_vehicle!,
    get_vehicle,

    # scene
    Scene,

    get_by_id,
    get_index_of_first_vehicle_with_id,
    get_first_available_id,
    record_length,
    get_scene,
    push_back_records!,
    update!,

    NeighborLongitudinalResult,
    get_neighbor_fore_along_lane,
    get_neighbor_fore_along_left_lane,
    get_neighbor_fore_along_right_lane,
    get_neighbor_rear_along_lane,
    get_neighbor_rear_along_left_lane,
    get_neighbor_rear_along_right_lane,

    FrenetRelativePosition,
    get_frenet_relative_position,

    # scene record
    SceneRecord,

    pastframe_inbounds

include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "utils.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "splines.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "curves.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "roadways.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "vehicles.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "trajdatas.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "scenes.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "core", "scene_records.jl"))

end # module


text:
using IntervalArithmetic
using Test
using StaticArrays

@testset "rand tests" begin

    X = 3..4
    for i in 1:100
        @test rand(X) ∈ X
    end

    Y = IntervalBox(3..4, 5..6)
    for i in 1:100
        @test rand(Y) ∈ Y
    end

    for T in (Float32, Float64, BigFloat)
        X = Interval{T}(3, 4)
        @test rand(X) isa T

        Y = IntervalBox(X, X)
        @test rand(Y) isa SVector{2,T}
    end

end


text:
## Claudia: commented these tests out because they were not part of runtests.jl
## when I added the test for readNexusTrees. So, I was not sure if they were
## not included for a reason (slow test?)

## @testset "test: read CF data" begin
## # d=readInputData("1.ms"); #tableCF0.txt
## d = (@test_logs PhyloNetworks.readInputData("1.ms",:all,0,["1","2","3","4","5","6"],false,"none",false,false);) #tableCF0.txt
## d = PhyloNetworks.readInputData("1.ms",:all,0,false); # writes summaryTreesQuartets.txt
## d=readInputData("1.ms",:rand,10); #tableCF3.txt
## d=readInputData("1.ms",[1,2,3,4,5]); #tableCF4.txt

## d=readInputData("1.ms","allQuartets.txt"); #tableCF1.txt
## d=readInputData("1.ms","allQuartets.txt",:rand,10); #tableCF2.txt
## d=readInputData("1.ms","allQuartets.txt",true,"try4.txt"); #try4.txt

## descData(d)

## d = readTrees2CF("1.ms");
## d = readTrees2CF("1.ms",filename="try.txt");
## d = readTrees2CF("1.ms","allQuartets.txt",filename="hola.txt");
## d = readTrees2CF("1.ms","allQuartets.txt",whichQ=:rand,numQ=10,filename="hola.txt");

## d= readTableCF("tableCF.txt");
## d2= readTableCF("tableCFbad.txt")
## end

@testset "test: reading nexus file" begin
nexusfile = joinpath(@__DIR__, "..", "examples", "test.nex")
# nexusfile = joinpath(dirname(pathof(PhyloNetworks)), "..","examples","test.nex")
vnet = readNexusTrees(nexusfile);
@test length(vnet) == 10
@test length(vnet[10].edge) == 10
@test vnet[10].edge[7].length ≈ 0.00035
vnet = readNexusTrees(nexusfile, PhyloNetworks.readTopologyUpdate, false, false);
@test length(vnet) == 10
@test length(vnet[10].edge) == 9
@test vnet[10].edge[7].length ≈ 0.00035
end

@testset "test: calculate quartet CF from input gene trees" begin
sixtreestr = ["(E,((A,B),(C,D)),O);","(((A,B),(C,D)),(E,O));","(A,B,((C,D),(E,O)));",
              "(B,((C,D),(E,O)));","((C,D),(A,(B,E)),O);","((C,D),(A,B,E),O);"]
sixtrees = readTopology.(sixtreestr)
df1 = writeTableCF(countquartetsintrees(sixtrees)...)
df2 = writeTableCF(readTrees2CF(sixtrees, writeTab=false, writeSummary=false))
o = [1,2,4,7,11,3,5,8,12,6,9,13,10,14,15]
@test df1 == df2[o,:]
q,t = countquartetsintrees(sixtrees, Dict("A"=>"AB", "B"=>"AB"); showprogressbar=false);
df1 = writeTableCF(q,t)
@test df1[!,:CF12_34] ≈ [0,0,2/3,2/3,1]
@test df1[!,:CF13_24] ≈ [0,0,1/3,1/3,0]
@test df1[!,:CF14_23] ≈ [1.,1,0,0,0]
@test df1[!,:ngenes]  ≈ [6.,6,6,6,6]
# again, but weight each allele
q,t = countquartetsintrees(sixtrees, Dict("A"=>"AB", "B"=>"AB"); weight_byallele=true, showprogressbar=false);
df1 = writeTableCF(q,t)
@test df1[!,:CF12_34] ≈ [0,0,7/11,7/11,1]
@test df1[!,:CF13_24] ≈ [0,0,4/11,4/11,0]
@test df1[!,:CF14_23] ≈ [1.,1,0,0,0]
@test df1[!,:ngenes]  ≈ [11.,11,11,11,6]
# different averaging: first across each set of 4 alleles, then
# across sets of 4 alleles that map to the same set of 4 species
CSV.write("tmp_qCF.csv", df2)
CSV.write("tmp_map.csv", DataFrame(allele = ["A","B"], species = ["AB","AB"]))
df2_byallele = (@test_logs (:warn, r"not all alleles were mapped") mapAllelesCFtable("tmp_map.csv", "tmp_qCF.csv"))
rm.(["tmp_qCF.csv","tmp_map.csv"]);
q = readTableCF!(df2_byallele);
df2 = writeTableCF(q) # 45×8 DataFrames.DataFrame
# df2[7:11,:]
# df12 = innerjoin(df1, df2, on=[:t1,:t2,:t3,:t4], makeunique=true)
# all([df12[:,4+i] ≈ df12[:,8+i] for i in 1:4]) # false: because averaging done differently by the 2 functions
@test df2 == DataFrame(
  t1=["AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","C"],
  t2=["AB__2","AB__2","AB__2","AB__2","AB__2","AB__2","C","C","C","D","D"],
  t3=["C","C","C","D","D","E","D","D","E","E","E"],
  t4=["D","E","O","E","O","O","E","O","O","O","O"],
  CF12_34=[1.0,0.75,1.0,0.75,1.0,0.75,0.0,0.0,(3/5+4/6)/2,(3/5+4/6)/2,1.0],
  CF13_24=[0,0.25,0,0.25,0,0,0,0,(2/5+2/6)/2,(2/5+2/6)/2,0.0],
  CF14_23=[0,0,0,0,0,0.25,1,1,0,0,0],
  ngenes=Union{Missing,Float64}[5,4,5,4,5,4,5.5,5.5,5.5,5.5,6]
)
# basic tests for QuartetT type
nCk = PhyloNetworks.nchoose1234(5)
oneQ = PhyloNetworks.QuartetT(1,3,4,6, [.92,.04,.04, 100], nCk)
@test string(oneQ) == "4-taxon set number 8; taxon numbers: 1,3,4,6\ndata: [0.92, 0.04, 0.04, 100.0]"
end

if false # was used to time `countquartetsintrees` vs `readTrees2CF`
dir = "/Users/ane/Documents/private/concordance/quartetNetwork/multiind/data"
treefile = joinpath(dir, "raxml_1387_sample_5species4alleles.tre")
tree = readMultiTopology(treefile); # 1387 trees
# extrema([t.numTaxa for t in tree]) # 4-16 taxa in each
@time df1 = writeTableCF(countquartetsintrees(tree)...)
# 0.139761 seconds (900.12 k allocations: 52.000 MiB, 11.52% gc time). 3876×8 DataFrames.DataFrame
@time df2 = writeTableCF(readTrees2CF(tree, writeTab=false, writeSummary=false))
# 13.154085 seconds (84.86 M allocations: 10.010 GiB, 7.21% gc time).  3876×8 DataFrames.DataFrame
df12 = innerjoin(df1, df2, on=[:t1,:t2,:t3,:t4], makeunique=true)
@test all([df12[:,4+i] == df12[:,8+i] for i in 1:4])
# using BenchmarkTools
# @benchmark countquartetsintrees(tree)
#=
memory estimate:  49.08 MiB
allocs estimate:  788504
--------------
minimum time:     110.907 ms (10.89% GC)
median time:      121.782 ms (11.20% GC)
mean time:        125.725 ms (14.00% GC)
maximum time:     227.248 ms (54.16% GC)
--------------
samples:          40
evals/sample:     1
=#
# @benchmark readTrees2CF(tree, writeTab=false, writeSummary=false)
#=
BenchmarkTools.Trial:
  memory estimate:  10.01 GiB
  allocs estimate:  84744840
  --------------
  minimum time:     13.547 s (8.63% GC)
  median time:      13.547 s (8.63% GC)
  mean time:        13.547 s (8.63% GC)
  maximum time:     13.547 s (8.63% GC)
  --------------
  samples:          1
  evals/sample:     1
=#
mappingfile = joinpath(dir, "strain2bin_map.csv")
using CSV
taxonmap = DataFrame(CSV.File(mappingfile); copycols=false) # 110×3 DataFrames.DataFrame
taxonmap = Dict(taxonmap[i,:allele] => taxonmap[i,:species] for i in 1:110)
@time df1 = writeTableCF(countquartetsintrees(tree, taxonmap; weight_byallele=true)...)
# 0.119289 seconds (698.57 k allocations: 43.305 MiB, 17.40% gc time). 5×8 DataFrames.DataFrame
## larger examples: 98 to 110 taxa, 1387 trees
tree = readMultiTopology(joinpath(dir, "raxml_1387.tre")) # 1387 trees, 98-110 taxa in each
@time df1 = writeTableCF(countquartetsintrees(tree)...)
# 1219.94 seconds = 20.3 min (298.45 M allocations: 8.509 GiB, 0.79% gc time). 5773185×8 DataFrames.DataFrame
## mid-size example: to be able to run the slower algorithm and compare times
tree = readMultiTopology(joinpath(dir, "raxml_1387_sample_13species4alleles.tre")); # 1387 trees, 19-40 taxa in each
@time df1 = writeTableCF(countquartetsintrees(tree)...)
# 5.639443 seconds (16.84 M allocations: 568.496 MiB, 7.50% gc time). 292825×8 DataFrame
# ~ 600 times faster
@time df2 = writeTableCF(readTrees2CF(tree, writeTab=false, writeSummary=false))
# 3365.783672 seconds = 50.1 min (13.43 G allocations: 2.665 TiB, 21.33% gc time). 292825×8 DataFrame
df12 = innerjoin(df1, df2, on=[:t1,:t2,:t3,:t4], makeunique=true)
@test df12[!,8] ≈ df12[!,12] # number of genes
hasdata = map(iszero, df12[!,8]) # sum: 34 four-taxon sets have data for 0 genes
df12[hasdata,5:12] # countquartetsintrees gives 0s, readTrees2CF gives NaN
@test all([df12[.!hasdata,4+i] ≈ df12[.!hasdata,8+i] for i in 1:3]) # true. yeah!
end


text:
module TestFunctorLens

include("preamble.jl")

using Kaleido: fmap, fcompose

@testset begin
    fst = FLens((f, obj) -> fmap(x -> (x, obj[2:end]...), f(obj[1])))
    @test get((1, 2, 3), fst) == 1
    @test set((1, 2, 3), fst, 100) == (100, 2, 3)
    @test modify(string, (1, 2, 3), fst) == ("1", 2, 3)
end

asflens(lens::Lens) =
    FLens((f, obj) -> fmap(x -> set(obj, lens, x), f(get(obj, lens))))

function test_lenses(obj, lens, flens)
    @test get(obj, lens) === get(obj, flens)
    @test set(obj, lens, 100) === set(obj, flens, 100)
    @test modify(string, obj, lens) === modify(string, obj, flens)
end

@testset for (obj, lens) in [
    ((1, 2), @lens _[1])
    ((a=1, b=2), @lens _.a)
]
    flens = asflens(lens)
    test_lenses(obj, lens, flens)
end

@testset begin
    obj = ((a=1, b=nothing), nothing)
    l1 = asflens(@lens _[1])
    l2 = asflens(@lens _.a)
    lens = l1 ∘ l2
    flens = fcompose(l1, l2)
    test_lenses(obj, lens, flens)
end

end  # module


text:
function lv_turbo(r::Vector{UInt64}, mask::UInt64)
  @turbo for i in 1:length(r)
    r[i] &= mask
  end
  r
end

@testset "Simple Miscellaneous" begin
  @show @__LINE__
  r1 = rand(UInt, 239);
  mask = ~(one(eltype(r1))<<(2))
  r2 = r1 .& mask;
  @test lv_turbo(r1, mask) == r2;
end


text:
order = CATEGORY_DEFAULT # TODO: move inside function

"""Return a dict where keys are a tuple of input parameters (fuel, unit_type) and values are
generator types."""
function get_generator_mapping(filename = nothing)
    if isnothing(filename)
        filename = GENERATOR_MAPPING_FILE
    end
    genmap = open(filename) do file
        YAML.load(file)
    end

    mappings = Dict{NamedTuple, String}()
    for (gen_type, vals) in genmap
        for val in vals
            pm =
                isnothing(val["primemover"]) ? nothing :
                uppercase(string(val["primemover"]))
            key = (fuel = val["fuel"], primemover = pm)
            if haskey(mappings, key)
                error(
                    "duplicate generator mappings: $gen_type $(key.fuel) $(key.primemover)",
                )
            end
            mappings[key] = gen_type
        end
    end

    return mappings
end

"""Return the generator category for this fuel and unit_type."""
function get_generator_category(fuel, primemover, mappings::Dict{NamedTuple, String})
    fuel = isnothing(fuel) ? nothing : uppercase(string(fuel))
    primemover = isnothing(primemover) ? nothing : uppercase(string(primemover))
    generator = nothing

    # Try to match the primemover if it's defined. If it's nothing then just match on fuel.
    for pm in (primemover, nothing), f in (fuel, nothing)
        key = (fuel = f, primemover = pm)
        if haskey(mappings, key)
            generator = mappings[key]
            break
        end
    end

    if isnothing(generator)
        @error "No mapping defined for generator fuel=$fuel primemover=$primemover"
    end

    return generator
end

"""
    generators = make_fuel_dictionary(system::PSY.System, mapping::Dict{NamedTuple, String})

This function makes a dictionary of fuel type and the generators associated.

# Arguments
- `c_sys5_re::PSY.System`: the system that is used to create the results
- `results::IS.Results`: results

# Key Words
- `categories::Dict{String, NamedTuple}`: if stacking by a different category is desired

# Example
results = solve_op_model!(OpModel)
generators = make_fuel_dictionary(c_sys5_re)

"""
function make_fuel_dictionary(sys::PSY.System, mapping::Dict{NamedTuple, String})
    generators = PSY.get_components(PSY.StaticInjection, sys, PSY.get_available)
    gen_categories = Dict()
    for category in unique(values(mapping))
        gen_categories["$category"] = []
    end

    for gen in generators
        gen isa PSY.ElectricLoad && continue
        fuel = hasmethod(PSY.get_fuel, Tuple{typeof(gen)}) ? PSY.get_fuel(gen) : nothing
        pm =
            hasmethod(PSY.get_prime_mover, Tuple{typeof(gen)}) ? PSY.get_prime_mover(gen) :
            nothing
        category = get_generator_category(fuel, pm, mapping)
        push!(gen_categories["$category"], (string(typeof(gen)), (PSY.get_name(gen))))
    end
    [delete!(gen_categories, "$k") for (k, v) in gen_categories if isempty(v)]
    return gen_categories
end

function make_fuel_dictionary(sys::PSY.System; kwargs...)
    mapping = get_generator_mapping(get(kwargs, :generator_mapping_file, nothing))
    return make_fuel_dictionary(sys, mapping)
end


text:
#=
#define FMPR_RND_DOWN  0    RoundingMode{:Down}()
#define FMPR_RND_UP    1
#define FMPR_RND_FLOOR 2
#define FMPR_RND_CEIL  3
#define FMPR_RND_NEAR  4
=#

#=
sigBitsToUse(prec::Int, sig::Int, base::Int) =
    min(prec, ifelse(signbit(sig), -sig, sig) * Float64(nextfloat(Float32(nextfloat(log2(base))))))
=#
sigBitsToUseRadix10(prec::Int, sig::Int) =
    min(prec, ifelse(signbit(sig), -sig, sig) * 3.3219285) # log2(10)

sigBitsToUseRadix2(prec::Int, sig::Int) =
    min(prec, ifelse(signbit(sig), -sig, sig) * 1)         # log2(2)

function sigBitsToUse(prec::Int, sig::Int, base::Int)
    if base==2
        return sigBitsToUseRadix2(prec, sig)
    elseif base==10
        return sigBitsToUseRadix10(prec, sig)
    else
        throw(ErrorException(string("Not Implemented for radix ",base)))
    end
end

function round{P}(x::ArbFloat{P}, sig::Int=P, base::Int=10)
    sig = max(1, abs(sig))
    sigbits = base==2 ? sig : sigBitsToUse(P, sig, base)
    z = ArbFloat{P}()
    ccall(@libarb(arb_set_round), Void,  (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}, Int), &z, &x, sigbits)
    return z
end

function ceil{P}(x::ArbFloat{P}, sig::Int=P, base::Int=10)
    sig = max(1, abs(sig))
    sigbits = sigBitsToUse(P, sig, base)
    z = ArbFloat{P}()
    ccall(@libarb(arb_ceil), Void,  (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}, Int), &z, &x, sigbits)
    return z
end

function floor{P}(x::ArbFloat{P}, sig::Int=P, base::Int=10)
    sig = max(1, abs(sig))
    sigbits = sigBitsToUse(P, sig, base)
    z = ArbFloat{P}()
    ccall(@libarb(arb_floor), Void,  (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}, Int), &z, &x, sigbits)
    return z
end

function trunc{P}(x::ArbFloat{P}, sig::Int=P, base::Int=10)
    sig = max(1, abs(sig))
    sigbits = sigBitsToUse(P, sig, base)
    z = ArbFloat{P}()
    if signbit(x)
        ccall(@libarb(arb_ceil), Void,  (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}, Int), &z, &x, sigbits)
    else
        ccall(@libarb(arb_floor), Void,  (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}, Int), &z, &x, sigbits)
    end
    return z
end

function round{I<:Integer,P}(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=10)
    z = round(x, sig, base)
    return convert(I, z)
end
function ceil{I<:Integer,P}(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=10)
    z = ceil(x, sig, base)
    return convert(I, z)
end
function floor{I<:Integer,P}(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=10)
    z = floor(x, sig, base)
    return convert(I, z)
end
function trunc{I<:Integer,P}(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=10)
    z = trunc(x, sig, base)
    return convert(I, z)
end

fld{P}(x::ArbFloat{P}, y::ArbFloat{P}) = convert(Int, floor(x/y))
cld{P}(x::ArbFloat{P}, y::ArbFloat{P}) = convert(Int, ceil(x/y))
div{P}(x::ArbFloat{P}, y::ArbFloat{P}) = convert(Int, trunc(x/y))

rem{P}(x::ArbFloat{P}, y::ArbFloat{P}) = convert(Int, x - div(x,y)*y)
mod{P}(x::ArbFloat{P}, y::ArbFloat{P}) = convert(Int, x - fld(x,y)*y)

function divrem{P}(x::ArbFloat{P}, y::ArbFloat{P})
   dv = div(x,y)
   r  = x - d*y
   rm = convert(Int, r)
   return dv,rm
end

function fldmod{P}(x::ArbFloat{P}, y::ArbFloat{P})
   fd = fld(x,y)
   m  = x - d*y
   md = convert(Int, m)
   return fd,md
end


text:
module ApFilterTrials

using DocStringExtensions

using ..Configuration
import ..Trials

"""
$(TYPEDSIGNATURES)

Return the elements of a particular hyperparameter's history (identified by `nid`) that
correspond to trials whose losses were above or below the threshold.
"""
function ap_filter_trials(
    nid::Symbol, trials::Vector{Trials.Trial}, config::Config
)::NTuple{2, Vector{T} where T <: Union{Float64, Int}}

    # Splitting is done this way to cope with duplicate loss values.
    # This is the number of below values that are extracted from trials by loss.
    n_below = min(Int(ceil(config.threshold * sqrt(length(trials)))), config.linear_forgetting)

    trials_by_loss = sort(trials; by = t -> t.loss)
    below = [
        trial.vals[nid]
            for trial in trials_by_loss[1:n_below]
                if nid in keys(trial.vals)
    ]
    above = [
        trial.vals[nid]
            for trial in trials_by_loss[n_below + 1:end]
                if nid in keys(trial.vals)
    ]
    if all(isinteger.(below)) && all(isinteger.(above))
        below = convert(Vector{Int}, below)
        above = convert(Vector{Int}, above)
    else
        below = convert(Vector{Float64}, below)
        above = convert(Vector{Float64}, above)
    end
    return below, above
end

end # module ApFilterTrials


text:
function soft_drop(cur_player::AbstractPlayer, is_step::Bool=false)
  cur_piece = cur_player.piece
  cur_grid = cur_player.grid

  did_move = drop!(
    cur_player, cur_piece, cur_grid,
    -1, true, is_step
  )

  did_move
end


text:
module ClosedGroupFunctions

    using ProgressMeter
    using Serialization
    using Bijections

    export
        group_generator_basic,

        label_generators,
        labelled_group_generator_simple,
        labelled_group_generator_shortest,

        calculate_conjugacy_classes,

        is_invariant,
        apply_invariant_to_all_in_class,
        apply_invariant_to_first_in_class,
        apply_invariant_to_first_in_all_classes,

        store_group,
        load_group,

        group_element_calculate

    include("group_generator_basic.jl")
    include("labelled_group_generator.jl")
    include("calculate_conjugacy_classes.jl")
    include("conjugacy_classes_invariants.jl")
    include("storage.jl")
    include("helpers.jl")

end # module


text:
# include("datastructs_module.jl")

# module tr # to be removed
# using ..datastructs
# include("tree_search_helper.jl") # to be removed

# mutable struct PtrTreeNode{T}
#     parent::PtrTreeNode
#     child
# end

mutable struct TreeNode{T}
    parent::TreeNode
    left_child::TreeNode
    right_sibling::TreeNode
    data::T

    function TreeNode{T}(content::T) where {T}
        obj = new()
        obj.data = content
        return obj
    end
    function TreeNode{T}() where {T}
        obj = new()
        return obj
    end
end

TreeNode(content::T) where {T} = TreeNode{T}(content::T)
TreeNode() = TreeNode(undef)

Base.eltype( ::Type{ <: TreeNode{T}} ) where {T} = T
Base.isempty(node::TreeNode) = node.data == undef

function has_parent(node::TreeNode)
    return isdefined(node, :parent)
end

function has_child(node::TreeNode)
    return isdefined(node, :left_child)
end

function has_right_sibling(node::TreeNode)
    return isdefined(node, :right_sibling)
end

function is_root(node::TreeNode)
    return !has_parent(node)
end

function is_leaf(node::TreeNode)
    return !has_child(node)
end

function is_empty(node::TreeNode)
    return !isdefined(node, :data)
end

function set_data(node::TreeNode, data)
    node.data = data
end

function get_data(node::TreeNode)
    return node.data
end

function Base.show(io::IO, node::TreeNode{T}) where {T}
    println("TreeNode{$(T)}")

    print("\tparent: ")
    if has_parent(node)
        isdefined(node.parent, :data) ? show(node.parent.data) : print("<no_data>\n")
    else
        print("not_defined\n")
    end

    print("\tleft_child: ")
    if has_child(node)
        isdefined(node.left_child, :data) ? show(node.left_child.data) : print("<no_data>\n")
    else
        print("not_defined\n")
    end

    print("\right_sibling: ")
    if has_right_sibling(node)
        isdefined(node.right_sibling, :data) ? show(node.right_sibling.data) : print("<no_data>\n")
    else
        print("not_defined\n")
    end

    print("\tdata: ")
    if !is_empty(node)
        show(node.data)
    else
        print("not_defined\n")
    end
end

function clone(node::TreeNode{T}) where {T}
    new_node = TreeNode{T}()

    if has_parent(node)
        new_node.parent         = node.parent
    end
    if has_child(node)
        new_node.left_child     = node.left_child
    end
    if has_right_sibling(node)
        new_node.right_sibling  = node.right_sibling
    end
    if !is_empty(node)
        new_node.data           = node.data
    end

    return new_node
end

function children_vector(node::TreeNode{T}) where {T}

    vec = Vector{TreeNode{T}}(undef, 0)
    push!(vec, node)

    while has_right_sibling(node)
        right = node.right_sibling
        push!(vec, right)
    end

    return vec
end

mutable struct Tree{T}
    count::Int64
    node::TreeNode{T}

    function Tree{T}(content::T) where {T}
        count_ = 1
        node_ = TreeNode{T}(content)
        return new(count_, node_)
    end
    function Tree{T}() where {T}
        count_ = 1
        node_ = TreeNode{T}()
        return new(count_, node_)
    end
end
Tree(content::T) where {T} = Tree{T}(content::T)

function is_empty(tree::Tree)
    return tree.count == 0
end

function count(tree::Tree)
    return tree.count
end

function Base.show(io::IO, tree::Tree{T}) where {T}
    println("Tree{$(T)}")
    println("\tcount: $(tree.count)")
end

function insert_after!(tree::Tree, node::TreeNode, new_node::TreeNode)

    if is_leaf(node)
        node.left_child = new_node
    else
        # has some children, singly linked list
        left_child             = node.left_child
        node.left_child        = new_node
        new_node.right_sibling = left_child
    end

    new_node.parent = node
    tree.count     += 1

    return tree
end

function select(tree::Tree, index_vec)

    node = tree.node

    ii = 1
    L = length(index_vec)
    while ii <= L

        idx = index_vec[ii]
        if idx == 1
            node = node.left_child
        else
            aa = 1
            while has_right_sibling(node)
                aa += 1
                node = node.right_sibling
                if aa >= idx
                    break
                end
            end
            if aa < idx
                error("There is not as many siblings at level: $(ii), index: $(idx)")
            end

        end
        ii += 1
    end

    return node
end

function root(tree::Tree)
    return tree.node
end

function depth_first_search(tree::Tree;
    visited_fcn = do_nothing,
    push_fcn = do_nothing)

    n        = tree.count
    depths   = zeros(Int, n)
    NodeType = typeof(tree.node)
    # nodes    = Vector{NodeType}(undef, n)
    nodes    = Vector{Any}(undef, n)

    init_size    = div(n, 2)
    helper       = TreeSearchHelper(1, false)
    # node_stack   = Stack{NodeType}(init_size)
    node_stack   = Stack{Any}(init_size)

    helper_stack = Stack{TreeSearchHelper}(init_size)
    push!(node_stack,   tree.node)
    push!(helper_stack, helper)

    index = 0
    while !isempty(node_stack)

        node   = pop!(node_stack)
        helper = pop!(helper_stack)

        if helper.visited

            index        += 1
            nodes[index]  = node
            depths[index] = helper.depth

            # call the callback function
            visited_fcn(node, helper.depth, helper.visited)
        else

            # add the same node once again
            helper = TreeSearchHelper(helper.depth, true)
            push!(node_stack,   node)
            push!(helper_stack, helper)

            # add the children
            next_depth = helper.depth + 1
            visited    = false

            if !has_child(node)
                continue
            end

            node = node.left_child
            push!(node_stack,   node)
            push!(helper_stack, TreeSearchHelper(next_depth, visited))
            push_fcn(node, next_depth, visited)

            while has_right_sibling(node)
                node = node.right_sibling
                push!(node_stack,   node)
                push!(helper_stack, TreeSearchHelper(next_depth, visited))
                push_fcn(node, next_depth, visited)
            end

        end

    end

    return (nodes, depths)
end

function breadth_first_search(tree::Tree;
    visited_fcn = do_nothing,
    push_fcn = do_nothing)

    n        = tree.count
    depths   = zeros(Int, n)
    NodeType = typeof(tree.node)
    nodes    = Vector{NodeType}(undef, n)

    init_size    = div(n, 2)
    helper       = TreeSearchHelper(1, false)
    node_queue   = Queue{NodeType}(init_size)
    helper_queue = Queue{TreeSearchHelper}(init_size)
    push!(node_queue,   tree.node)
    push!(helper_queue, helper)

    index = 0
    while !isempty(node_queue)

        node   = pop!(node_queue)
        helper = pop!(helper_queue)

        # println("Popped: ")
        # show(node)

        if helper.visited
            index        += 1
            nodes[index]  = node
            depths[index] = helper.depth

            # call the callback function
            visited_fcn(node, helper.depth, helper.visited)
        else

            # add the same node once again
            helper = TreeSearchHelper(helper.depth, true)
            push!(node_queue,   node)
            push!(helper_queue, helper)

            # add the children
            next_depth = helper.depth + 1
            visited    = false

            if !has_child(node)
                continue
            end

            node = node.left_child
            push!(node_queue,   node)
            push!(helper_queue, TreeSearchHelper(next_depth, visited))
            push_fcn(node, next_depth, visited)

            while has_right_sibling(node)
                node = node.right_sibling
                push!(node_queue,   node)
                push!(helper_queue, TreeSearchHelper(next_depth, visited))
                push_fcn(node, next_depth, visited)
            end
        end

    end

    return (nodes, depths)
end

# function leaves(tree::Tree)
#     # needs e.g. dfs
# end

# function remove_node!(tree::Tree, node::TreeNode)
#     if is_root(node)
#         tree.count = 0
#         return tree
#     end

#     if is_leaf(node)
#         # easy case
#         parent = node.parent

#
#     else
#         # hard case
#     end

#     tree.count -= 1
#     return tree
# end

# function merge!(tree1::TreeNode, tree2::TreeNode)
#     # basically a simultanous breadth first search
#     return tree1
# end

# Creation
#
# clone - shallow copy -> is not implemented because it is dangerous
# split - not implemented yet
# sub tree - not implemented yet
#
# State
# is_empty
# count
# show
#
# Selection
# select == find_indices -> basically the path from the root
# leaves
# root
#
#
# Insertion
# insert_after
# insert_before - not implemented yet
# merge - not implemented yet
#
# Deletion
# remove_node
#
# search from node - bfs, dfs -> apply functions include

# end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module Utilities

using LinearAlgebra # For dot

using MathOptInterface
const MOI = MathOptInterface

const MOIU = MOI.Utilities # used in macro

const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}
const SQF{T} = MOI.ScalarQuadraticFunction{T}
const VQF{T} = MOI.VectorQuadraticFunction{T}

const VI = MOI.VariableIndex
const CI{F,S} = MOI.ConstraintIndex{F,S}

include("functions.jl")
include("sets.jl")
include("constraints.jl")
include("copy.jl")
include("results.jl")

include("model.jl")
include("parser.jl")
include("mockoptimizer.jl")
include("cachingoptimizer.jl")
include("universalfallback.jl")

end # module


text:
# Autogenerated wrapper script for LCIO_jll for powerpc64le-linux-gnu-cxx03
export liblcio, libsio

using Zlib_jll
JLLWrappers.@generate_wrapper_header("LCIO")
JLLWrappers.@declare_library_product(liblcio, "liblcio.so.2.17")
JLLWrappers.@declare_library_product(libsio, "libsio.so.0.0")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        liblcio,
        "lib/liblcio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsio,
        "lib/libsio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using RandomBooleanMatrices
using SparseArrays
using Random
using Test

Random.seed!(1337)

@testset "curveball" begin
    m = sprand(Bool, 8, 6, 0.2)
    m_old = copy(m)
    rsm = sum(m, dims = 1)
    csm = sum(m, dims = 2)
    randomize_matrix!(m, method = curveball)

    @test rsm == sum(m, dims = 1)
    @test csm == sum(m, dims = 2)
    @test m_old != m

    m2 = rand(0:1, 6, 5)
    rsm = sum(m2, dims = 1)
    csm = sum(m2, dims = 2)
    rmg = matrixrandomizer(m2, method = curveball)
    m3 = rand(rmg)

    @test rsm == sum(m3, dims = 1)
    @test csm == sum(m3, dims = 2)

    m4 = rand(rmg)

    @test m3 != m4
end


text:
include("containers.jl")
include("math.jl")

# Enable debug timings `@timeit_debug timer() "name" stuff...`.
# This allows us to disable timings completely by executing
# `TimerOutputs.disable_debug_timings(Trixi)`
# and to enable them again by executing
# `TimerOutputs.enable_debug_timings(Trixi)`
timeit_debug_enabled() = true

# Store main timer for global timing of functions
const main_timer = TimerOutput()

# Always call timer() to hide implementation details
timer() = main_timer

"""
    PerformanceCounter

A `PerformanceCounter` be used to track the runtime performance of some calls.
Add a new runtime measurement via `put!(counter, runtime)` and get the averaged
runtime of all measurements added so far via `take!(counter)`, resetting the
`counter`.
"""
mutable struct PerformanceCounter
  ncalls_since_readout::Int
  runtime::Float64
end

PerformanceCounter() = PerformanceCounter(0, 0.0)

function Base.take!(counter::PerformanceCounter)
  time_per_call = counter.runtime / counter.ncalls_since_readout
  counter.ncalls_since_readout = 0
  counter.runtime = 0.0
  return time_per_call
end

function Base.put!(counter::PerformanceCounter, runtime::Real)
  counter.ncalls_since_readout += 1
  counter.runtime += runtime
end

"""
    examples_dir()

Return the directory where the example files provided with Trixi.jl are located. If Trixi is
installed as a regular package (with `]add Trixi`), these files are read-only and should *not* be
modified. To find out which files are available, use, e.g., `readdir`:

# Examples
```julia
julia> readdir(examples_dir())
5-element Array{String,1}:
 "1d"
 "2d"
 "3d"
 "README.md"
 "paper-self-gravitating-gas-dynamics"
```
"""
examples_dir() = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

"""
    get_examples()

Return a list of all example elixirs that are provided by Trixi.
"""
function get_examples()
  examples = String[]
  for (root, dirs, files) in walkdir(examples_dir())
    for f in files
      if startswith(f, "elixir_") && endswith(f, ".jl")
        push!(examples, joinpath(root, f))
      end
    end
  end

  return examples
end

"""
    default_example()

Return the path to an example elixir that can be used to quickly see Trixi in action.
"""
default_example() = joinpath(examples_dir(), "2d", "elixir_advection_basic.jl")

# Print informative message at startup
function print_startup_message()
  s = """

    ████████╗██████╗ ██╗██╗  ██╗██╗
    ╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
       ██║   ██████╔╝██║ ╚███╔╝ ██║
       ██║   ██╔══██╗██║ ██╔██╗ ██║
       ██║   ██║  ██║██║██╔╝ ██╗██║
       ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝
    """
  mpi_println(s)
end

"""
    get_name(x)

Returns a name of `x` ready for pretty printing.
By default, return `string(y)` if `x isa Val{y}` and return `string(x)` otherwise.

# Examples

```jldoctest
julia> Trixi.get_name("test")
"test"

julia> Trixi.get_name(Val(:test))
"test"
```
"""
get_name(x) = string(x)
get_name(::Val{x}) where x = string(x)

"""
    @threaded for ... end

Basically the same as `Threads.@threads` with an additional check whether only
one thread is used to reduce the overhead of serial execution.

Some discussion can be found at https://discourse.julialang.org/t/overhead-of-threads-threads/53964
"""
macro threaded(expr)
  # esc(quote ... end) as suggested in https://github.com/JuliaLang/julia/issues/23221
  return esc(quote
    if Threads.nthreads() == 1
      $(expr)
    else
      Threads.@threads $(expr)
    end
  end)
end


text:
MAP(fstate::NestedParticleState; kwargs...) = MAP(fstate.model, kwargs...)

function MAP(model::BinomialGridModel; kwargs...)
    MAP(BinomialModel(model), kwargs...)
end

function MAP(model::BinomialModel; marginal::Bool = false)
    if marginal
        return BinomialModel(
            mode(Array(model.N)),
            mode(Array(model.p)),
            mode(Array(model.q)),
            mode(Array(model.σ)),
            mode(Array(model.τ))
        )
    end
    v = hcat(
            model.N,
            model.p,
            model.q,
            model.σ,
            model.τ
        ) |> Array |> eachrow |> mode
    return BinomialModel(v...)
end


text:
module YetAnotherHelperFlagKillBarrier

using ..Ahorn, Maple

@mapdef Entity "YetAnotherHelper/FlagKillBarrier" FlagKillBarrier(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, falseColorHex::String="FFFFFF", trueColorHex::String="FF0000", flag::String="example_flag")

const placements = Ahorn.PlacementDict(
    "Flag Kill Barrier (Yet Another Helper)" => Ahorn.EntityPlacement(
        FlagKillBarrier,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::FlagKillBarrier) = 8, 8
Ahorn.resizable(entity::FlagKillBarrier) = true, true

function Ahorn.selection(entity::FlagKillBarrier)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagKillBarrier, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.79,0.38,0.38, 0.6), (0.79,0.32,0.30, 0.7))
end

end

max_stars_repo_name: YetAnotherHelper/YetAnotherHelper

max_issues_repo_name: YetAnotherHelper/YetAnotherHelper

max_forks_repo_name: YetAnotherHelper/YetAnotherHelper


text:
include("../../ApproximateCOSMO.jl/src/COSMO.jl")
using Main.COSMO, JLD2
include("../src/sdpa.jl")
using Printf
using DataFrames, CSV
using Glob
using Base.Filesystem
# include("../utils/plot_utils.jl")

@show LinearAlgebra.BLAS.openblas_get_config()

filepaths = glob("../data/sdpa/*.dat-s")

df = DataFrame(Name = String[], MaxPsdDim = Int[], Rank = Int[], Time = Float64[], TimeLOBCPG = Float64[],
    TimePrj = Float64[], TimePrjLOBPCG = Float64[],
    Iter = Int[], IterLOBPCG = Int[],
    Objective = Float64[], ObjectiveLobpcg = Float64[],
    OptimalObjective = Float64[], Status = [], StatusLobpcg = []
)

filepaths = [filepaths[1]; filepaths]
# filepaths = filepaths[20:end]
for filepath in filepaths[1:end]
    # if basename(filepath)[1:end-6] == "maxG55" #&& !in(basename(filepath)[1:end-6], ["maxG60" "maxG55" "qpG51" "maxG32"])
    c, F, A, b, optimal_objective = load_sdpa_file(filepath)

    _, _, objective_lobpcg, status_lobpcg, solver_lobpcg = solve_sdpa_jump_dual(c, F, A, b, COSMO.Optimizer,
        psd_projector = COSMO.with_options(COSMO.PsdConeTriangleLOBPCG, verbosity=0, max_iter=10),
        eps_abs = 1e-4, eps_rel = 1e-4,
        verbose = true,
        adaptive_rho_tolerance = 10.0,
        # max_iter = 5000000, check_termination = 100000,
        # adaptive_rho = false, rho=1.2
    )
    max_psd_dim = -1
    max_rank = -1
    for set in solver_lobpcg.p.C.sets
        if isa(set, COSMO.PsdConeTriangleLOBPCG)
            max_psd_dim = max(max_psd_dim, set.sqrt_dim)
            max_rank = max(max_rank, size(set.U, 2))
        end
    end
    t_lobpcg = solver_lobpcg.times.solver_time
    t_proj_lobpcg = solver_lobpcg.times.proj_time
    iterations_lobpcg = solver_lobpcg.iterations

    _, _, objective, status, solver = solve_sdpa_jump_dual(c, F, A, b, COSMO.Optimizer,
        eps_abs = 1e-4, eps_rel = 1e-4,
        verbose = true,
        adaptive_rho_tolerance = 10.0,
        # max_iter = 5000000, check_termination = 100000,
        # adaptive_rho = false, rho=1.2
    )
    t = solver.times.solver_time
    t_proj = solver.times.proj_time
    iterations = solver.iterations
    #=
    status = NaN
    objective = NaN
    t = NaN
    t_proj = NaN
    iterations = -1
    =#

    push!(df, [basename(filepath)[1:end-6],
    max_psd_dim, max_rank,
    t, t_lobpcg,
    t_proj, t_proj_lobpcg,
    iterations, iterations_lobpcg,
     objective, objective_lobpcg, optimal_objective,
     string(status), string(status_lobpcg)])
    df |> CSV.write(string("results.csv"))
    flush(stdout)
    # end
end

#=
using MosekTools
solution, solution_objective = solve_sdpa_jump_dual(c, F, A, b, Mosek.Optimizer)

using ProxSDP
for i = 1:2
    solution, solution_objective = solve_sdpa_jump_dual(c, F, A, b, ProxSDP.Optimizer, log_verbose=true,
        tol_primal = 1e-4, tol_dual = 1e-4)
end
=#


text:
sniffslope(df::DataFrameRate, σ::Quantity) = slope(smooth(df.sniff, σ; fs=df.fs))

function sniffgrams(sniff, tds::AbstractVector{TrialData}; fs=default_fs, sniffwindow=0s..2s, Δf=1/width(sniffwindow))
    sgs = Matrix{eltype(sniff)}[]
    freq = time = nothing
    n = ceil(Int, convert(AbstractFloat, fs/Δf))
    for td in tds
        ostart = minimum(td.odorinterval)
        tr = first(td.trialrange) .+ iv2r(ostart + minimum(sniffwindow) .. ostart + maximum(sniffwindow), fs)
        if length(tr) < n
            tr = first(tr) : first(tr) + n -1
        end
        if last(tr) <= lastindex(sniff)
            sg = Periodograms.spectrogram(view(sniff, tr), n; fs=fs/(1Hz))
            freq, time = sg.freq, sg.time
            push!(sgs, sg.power)
        end
    end
    nt = length(time)
    Δt = width(sniffwindow)/nt
    return range(minimum(sniffwindow)+Δt/2, stop=maximum(sniffwindow)-Δt/2, length=nt), freq*Hz, sgs
end
sniffgrams(dfr::DataFrameRate, tds::AbstractVector{TrialData}; kwargs...) = sniffgrams(dfr.df.sniff, tds; fs=dfr.fs, kwargs...)

function sniffgrams_power(freq, sgs; finterval=4Hz .. 10Hz)
    idx = searchsortedfirst(freq, minimum(finterval)) : searchsortedfirst(freq, maximum(finterval))
    return [sum(sg[idx,:]) for sg in sgs]
end


text:
function build_despot(p::DESPOTPlanner, b_0)
    D = DESPOT(p, b_0)
    b = 1
    trial = 1
    start = CPUtime_us()

    while D.mu[1]-D.l[1] > p.sol.epsilon_0 &&
          CPUtime_us()-start < p.sol.T_max*1e6 &&
          trial <= p.sol.max_trials
        b = explore!(D, 1, p)
        backup!(D, b, p)
        trial += 1
    end

    return D
end

function explore!(D::DESPOT, b::Int, p::DESPOTPlanner)
    while D.Delta[b] <= p.sol.D &&
        excess_uncertainty(D, b, p) > 0.0 &&
        !prune!(D, b, p)

        if isempty(D.children[b]) # a leaf
            expand!(D, b, p)
        end
        b = next_best(D, b, p)
    end
    if D.Delta[b] > p.sol.D
        make_default!(D, b)
    end
    return b
end

function prune!(D::DESPOT, b::Int, p::DESPOTPlanner)
    x = b
    blocked = false
    while x != 1
        n = find_blocker(D, x, p)
        if n > 0
            make_default!(D, x)
            backup!(D, x, p)
            blocked = true
        else
            break
        end
        x = D.parent_b[x]
    end
    return blocked
end

function find_blocker(D::DESPOT, b::Int, p::DESPOTPlanner)
    len = 1
    bp = D.parent_b[b] # Note: unlike the normal use of bp, here bp is a parent following equation (12)
    while bp != 1
        left_side_eq_12 = length(D.scenarios[bp])/p.sol.K*discount(p.pomdp)^D.Delta[bp]*D.U[bp] - D.l_0[bp]
        if left_side_eq_12 <= p.sol.lambda * len
            return bp
        else
            bp = D.parent_b[bp]
            len += 1
        end
    end
    return 0 # no blocker
end

function make_default!(D::DESPOT, b::Int)
    l_0 = D.l_0[b]
    D.mu[b] = l_0
    D.l[b] = l_0
end

function backup!(D::DESPOT, b::Int, p::DESPOTPlanner)
    # Note: maybe this could be sped up by just calculating the change in the one mu and l corresponding to bp, rather than summing up over all bp
    while b != 1
        ba = D.parent[b]
        b = D.parent_b[b]

        # https://github.com/JuliaLang/julia/issues/19398
        #=
        D.ba_mu[ba] = D.ba_rho[ba] + sum(D.mu[bp] for bp in D.ba_children[ba])
        =#
        sum_mu = 0.0
        for bp in D.ba_children[ba]
            sum_mu += D.mu[bp]
        end
        D.ba_mu[ba] = D.ba_rho[ba] + sum_mu

        #=
        max_mu = maximum(D.ba_rho[ba] + sum(D.mu[bp] for bp in D.ba_children[ba]) for ba in D.children[b])
        max_l = maximum(D.ba_rho[ba] + sum(D.l[bp] for bp in D.ba_children[ba]) for ba in D.children[b])
        =#
        max_U = -Inf
        max_mu = -Inf
        max_l = -Inf
        for ba in D.children[b]
            weighted_sum_U = 0.0
            sum_mu = 0.0
            sum_l = 0.0
            for bp in D.ba_children[ba]
                weighted_sum_U += length(D.scenarios[bp]) * D.U[bp]
                sum_mu += D.mu[bp]
                sum_l += D.l[bp]
            end
            new_U = (D.ba_Rsum[ba] + discount(p.pomdp) * weighted_sum_U)/length(D.scenarios[b])
            new_mu = D.ba_rho[ba] + sum_mu
            new_l = D.ba_rho[ba] + sum_l
            max_U = max(max_U, new_U)
            max_mu = max(max_mu, new_mu)
            max_l = max(max_l, new_l)
        end

        l_0 = D.l_0[b]
        D.U[b] = max_U
        D.mu[b] = max(l_0, max_mu)
        D.l[b] = max(l_0, max_l)
    end
end

function next_best(D::DESPOT, b::Int, p::DESPOTPlanner)
    max_mu = -Inf
    best_ba = first(D.children[b])
    for ba in D.children[b]
        mu = D.ba_mu[ba]
        if mu > max_mu
            max_mu = mu
            best_ba = ba
        end
    end

    max_eu = -Inf
    best_bp = first(D.ba_children[best_ba])
    for bp in D.ba_children[best_ba]
        eu = excess_uncertainty(D, bp, p)
        if eu > max_eu
            max_eu = eu
            best_bp = bp
        end
    end

    return best_bp

    # ai = indmax(D.ba_mu[ba] for ba in D.children[b])
    # ba = D.children[b][ai]
    # zi = indmax(excess_uncertainty(D, bp, p) for bp in D.ba_children[ba])
    # return D.ba_children[ba][zi]
end

function excess_uncertainty(D::DESPOT, b::Int, p::DESPOTPlanner)
    return D.mu[b]-D.l[b] - length(D.scenarios[b])/p.sol.K * p.sol.xi * (D.mu[1]-D.l[1])
end


text:
struct Spectrum

    amps
    fs
    power
    signal::Signal

    framerate::Number
    #start # start time of the signal
    d # length of the converted wave
end

struct Integrated_Spectrum

    cumulated
    fs

end

function Integrated_Spectrum(spec::Spectrum)

    cs = cumsum(spec.power)

    cs /= cs[end] #normalizing, so total power = 1

    return Integrated_Spectrum(cs, spec.fs)
end

function filter(spec::Spectrum, type::Symbol, args...)

    amps = deepcopy(spec.amps)

    if type == :band
        @assert length(args) == 2
        @assert args[1] < args[2]
        amps[spec.fs .< args[1] ] .= 0
        amps[spec.fs .> args[2] ] .= 0
    end

    if type == :low
        @assert length(args) == 1
        amps[spec.fs .> args[1] ] .= 0
    end

    if type == :high
        @assert length(args) == 1
        amps[spec.fs .< args[1] ] .= 0
    end

    return Spectrum(amps, spec.fs, abs2.(amps), spec.signal, spec.framerate, spec.d)

end

@recipe function plot(spec::Spectrum;noise=false)

    #=
    if haskey(w.signal.info, "name")
        label --> (w.signal.info["name"]*" framerate $(w.framerate)")
    else
        label --> "Framerate $(w.framerate)"
    end
    =#
    xguide --> "Frequeny in Hz"
    legend --> :none
    if noise
        xscale --> :log10
        yscale --> :log10
        xlims --> [0.01,spec.fs[end]*10]
        ylims --> [(spec.power |> median) * 1e-5, (spec.power |> maximum) * 10 ]
        yguide --> "Power"
        spec.fs, abs.(spec.amps).^2
    else
        yguide --> "Amplitude"
        spec.fs, abs.(spec.amps)
    end

end

@recipe function plot(ispec::Integrated_Spectrum)

    #=
    if haskey(w.signal.info, "name")
        label --> (w.signal.info["name"]*" framerate $(w.framerate)")
    else
        label --> "Framerate $(w.framerate)"
    end
    =#
    grid --> true
    linewidth --> 3
    label --> :none
    yguide --> "Cumulative Fraction of Power"
    xguide --> "Frequeny in Hz"
    ispec.fs, ispec.cumulated

end


text:
using DESPOTAlpha
using Test

using POMDPs
using POMDPModels
using POMDPSimulators
using Random
using POMDPModelTools
using ParticleFilters
using BeliefUpdaters
using POMDPPolicies
using Plots
theme(:mute)

# include("baby_sanity_check.jl")

pomdp = BabyPOMDP()
pomdp.discount = 1.0
p = solve(DESPOTAlphaSolver(tree_in_info=true), pomdp)
Random.seed!(p, 1)
b0 = initialstate(pomdp)
a, info = action_info(p, b0)
info_analysis(info)

# Light Dark Test
# pomdp = LightDark1D()
# random = solve(RandomSolver(), pomdp)
# bds = IndependentBounds(FORollout(random), pomdp.correct_r)
# solver = DESPOTAlphaSolver(bounds=bds,
#                     K=100,
#                     C=20,
#                     tree_in_info=true
#                     )
# planner = solve(solver, pomdp)
# a, info = action_info(planner, initialstate(pomdp))
# info_analysis(info)
# hr = HistoryRecorder(max_steps=50)
# @time hist = simulate(hr, pomdp, planner)
# # hist_analysis(hist)
# println("Discounted reward is $(discounted_reward(hist))")

# BabyPOMDP Test
# Type stability
pomdp = BabyPOMDP()
bds = (pomdp, b)->(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = DESPOTAlphaSolver(bounds=bds,
                      rng=MersenneTwister(4),
                      K=100,
                      tree_in_info=true
                     )
p = solve(solver, pomdp)

b0 = initialstate(pomdp)
D, Depth = @inferred DESPOTAlpha.build_tree(p, b0)
@inferred action_info(p, b0)
a, info = action_info(p, b0)
info_analysis(info)
@inferred DESPOTAlpha.explore!(D, 1, Bool[], p)
@inferred DESPOTAlpha.expand!(D, D.b-1, DESPOTAlpha.retrieve_as(D, D.b-1), p)
@inferred DESPOTAlpha.backup!(D, D.b-1, DESPOTAlpha.retrieve_as(D, D.b-1), p)
@inferred DESPOTAlpha.next_best(D, 1, p)
@inferred DESPOTAlpha.excess_uncertainty(D, 1, DESPOTAlpha.solver(p).xi, p)
@inferred action(p, b0)

pomdp = BabyPOMDP()

# constant bounds
bds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = DESPOTAlphaSolver(bounds=bds, K=200, tree_in_info=true, num_b=10_000)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=20)
@time hist = simulate(hr, pomdp, planner)
# hist_analysis(hist)
println("Discounted reward is $(discounted_reward(hist))")

# FO policy lower bound
bds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), FORollout(FeedWhenCrying()))
solver = DESPOTAlphaSolver(bounds=bds, K=200, tree_in_info=true, num_b=10_000)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=20)
@time hist = simulate(hr, pomdp, planner)
# hist_analysis(hist)
println("Discounted reward is $(discounted_reward(hist))")

# from README:
using POMDPs, POMDPModels, POMDPSimulators, DESPOTAlpha

pomdp = TigerPOMDP()

solver = DESPOTAlphaSolver(bounds=IndependentBounds(-20.0, 0.0))
planner = solve(solver, pomdp)

for (s, a, o) in stepthrough(pomdp, planner, "s,a,o", max_steps=10)
    println("State was $s,")
    println("action $a was taken,")
    println("and observation $o was received.\n")
end


text:
# PS6 Dayuan Xie

using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using DataFramesMeta
using CSV

function q3e()
# Question 1
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body)

df = @transform(df, bus_id = 1:size(df,1))

dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

# Question 2
θ̂_glm = glm(@formula(Y ~ Odometer*RouteUsage*Branded*time + Odometer^2*RouteUsage^2*Branded*time^2), df_long, Binomial(), LogitLink())
println(θ̂_glm)

# Question 3a
function create_grids()

    function xgrid(theta,xval)
        N      = length(xval)
        xub    = vcat(xval[2:N],Inf)
        xtran1  = zeros(N,N)
        xtran1c = zeros(N,N)
        lcdf   = zeros(N)
        for i=1:length(xval)
            xtran1[:,i]   = (xub[i] .>= xval).*(1 .- exp.(-theta*(xub[i] .- xval)) .- lcdf)
            lcdf        .+= xtran1[:,i]
            xtran1c[:,i] .+= lcdf
        end
        return xtran1, xtran1c
    end

    zval   = collect([.25:.01:1.25]...)
    zbin   = length(zval)
    xval   = collect([0:.125:25]...)
    xbin   = length(xval)
    tbin   = xbin*zbin
    xtran  = zeros(tbin,xbin)
    xtranc = zeros(xbin,xbin,xbin)
    for z=1:zbin
        xtran[(z-1)*xbin+1:z*xbin,:],xtranc[:,:,z] = xgrid(zval[z],xval)
    end

    return zval,zbin,xval,xbin,xtran
end

zval,zbin,xval,xbin,xtran = create_grids()

#Question 3b
Odometer=kron(ones(zbin),xval)
RouteUsage= kron(ones(xbin),zval)
Branded=zeros(size(Odometer))
time=zeros(size(Branded))

df1 = DataFrame(Odometer=Odometer,RouteUsage=RouteUsage,Branded=Branded,time=time)

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body)
Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
Z = Vector(df[:,:RouteUsage])
B = Vector(df[:,:Branded])
N = size(Y,1)
T = size(Y,2)
Xstate = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Zstate = Vector(df[:,:Zst])

FV=zeros(zbin*xbin,2,T+1)

for t= 2:T
    for b= 0:1
        for i=1:size(Branded)
        df1.time[i]=t
        df1.Branded[i]=b
        p0=predict(θ̂_glm,df1)
        row = x + (z-1)*xbin
        FV[row,b+1,t] =-β*log(p0)
        end
    end
end

FVT1=zeros(N,T)
for i= 1:N
    for t= 1:T
    row0 = (Zstate[i]-1)*xbin+1
    row1  = Xstate[i,t] + (Zstate[i]-1)*xbin
    FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,B[i]+1,t+1]
    end
end

fvt1= FVT1'[:]

# Question 3c
df_long = @transform(df_long,fv = fvt1)
theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),df_long, Binomial(), LogitLink(),offset=df_long.fv)

# Question 3d

end

# Question 3e
@time q3e()

max_stars_repo_path: ProblemSets/PS6-ccp/PS6_ Dayuan Xie.jl

max_issues_repo_path: ProblemSets/PS6-ccp/PS6_ Dayuan Xie.jl

max_forks_repo_path: ProblemSets/PS6-ccp/PS6_ Dayuan Xie.jl


text:
##############################################################################
##
## The following function follows the command ranktest (called by ivreg2)
## RANKTEST: Stata module to test the rank of a matrix using the Kleibergen-Paap rk statistic
## Authors: Frank Kleibergen, Mark E Schaffer
## IVREG2: Stata module for extended instrumental variables/2SLS and GMM estimation
## Authors: Christopher F Baum, Mark E Schaffer, Steven Stillman
## More precisely, it corresponds to the Stata command:  ranktest  (X) (Z), wald full
##############################################################################

function ranktest!(X::Matrix{Float64},
                    Z::Matrix{Float64},
                    Pi::Matrix{Float64},
                    vcov_method::AbstractVcovMethod,
                    df_small::Int,
                    df_absorb::Int)

    K = size(X, 2)
    L = size(Z, 2)

    crossz = cholesky!(Symmetric(Z' * Z))
    crossx = cholesky!(Symmetric(X' * X))

    Fmatrix = crossz.U
    Gmatrix = crossx.U
    theta = Fmatrix * (Gmatrix' \ Pi')'

    svddecomposition = svd(theta, full = true)
    u = svddecomposition.U
    vt = svddecomposition.Vt

    # compute lambda
    if K == 1
        a_qq = sqrt(u * u')
        b_qq = sqrt(vt * vt')
    else
        u_12 = u[1:(K-1),(K:L)]
        v_12 = vt[1:(K-1),K]
        u_22 = u[(K:L),(K:L)]
        v_22 = vt[K,K]
        a_qq = vcat(u_12, u_22) * (u_22 \ sqrt(u_22 * u_22'))
        b_qq = sqrt(v_22 * v_22') * (v_22' \ vcat(v_12, v_22)')
    end
    kronv = kron(b_qq, a_qq')
    lambda = kronv * vec(theta)

    # compute vhat
    if typeof(vcov_method) == VcovSimpleMethod
        vhat= Matrix(I/size(X, 1), L * K, L * K)
    else
        temp1 = convert(Matrix{eltype(Gmatrix)}, Gmatrix)
        temp2 = convert(Matrix{eltype(Fmatrix)}, Fmatrix)
        k = kron(temp1, temp2)'
        vcovmodel = VcovData(Z, k, X, size(Z, 1) - df_small - df_absorb)
        matrix_vcov2 = shat!(vcov_method, vcovmodel)
        vhat = k \ (k \ matrix_vcov2)'
    end

    # return statistics
    # why do I need to add Hermitian? (since 0.5)
    vlab = cholesky!(Hermitian(kronv * vhat * kronv'))
    r_kp = lambda' * (vlab \ lambda)
    p_kp = ccdf(Chisq((L-K+1 )), r_kp[1])
    F_kp = r_kp[1] / size(Z, 2)
    return F_kp, p_kp
end


text:
# 2D Objects
struct Circle <: Terminal
  r::Real
  center::Bool
end
struct Square <: Terminal
    pos::AbstractArray{Real, 1}
    center::Bool
end
struct Polygon
    fields
end

toString(c::Circle) = "circle(r = $(c.r), center = $(c.center));"
toString(s::Square) = "square(size = [$(s.pos[1]), $(s.pos[2])], center = $(s.center));"


text:
using GeoData, Test

A = [missing 7.0f0; 2.0f0 missing]
ga = GeoArray(A, (X, Y); missingval=missing)
ga99 = replace_missing(ga, -9999)
@test all(ga99 .=== [-9999.0f0 7.0f0; 2.0f0 -9999.0f0])
gaNaN = replace_missing(ga, NaN32)
@test all(gaNaN .=== [NaN32 7.0f0; 2.0f0 NaN32])
gaMi = replace_missing(ga, missing)
@test all(gaMi .=== ga)

@testset "boolmask" begin
    @test boolmask(A) == [false true; true false]
    @test boolmask(ga) == [false true; true false]
    @test boolmask(ga99) == [false true; true false]
    @test boolmask(gaNaN) == [false true; true false]
    @test dims(boolmask(ga)) == (X(Base.OneTo(2), mode=NoIndex()), Y(Base.OneTo(2), mode=NoIndex()))
end

@testset "missingmask" begin
    @test all(missingmask(A) .=== [missing true; true missing])
    @test all(missingmask(ga) .=== [missing true; true missing])
    @test all(missingmask(ga99) .=== [missing true; true missing])
    @test all(missingmask(gaNaN) .=== [missing true; true missing])
    @test dims(missingmask(ga)) == (X(Base.OneTo(2), mode=NoIndex()), Y(Base.OneTo(2), mode=NoIndex()))
end


text:
module mcssf1
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNonlinearModule: estimatesoundspeed
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using BenchmarkTools
using Test

function test()
    mr = DeforModelRed3D
    mass_density, E, nu = 8000*phun("kg/m^3"), 200000.0*phun("MPa"), 0.3
    m = FinEtoolsDeforNonlinear.MatDeforNeohookeanModule.MatDeforNeohookean(mr, mass_density, E, nu)
    ss = estimatesoundspeed(m)
    @test abs(ss - 5801.193511153214) / 5801.193511153214 < 1.0e-6
    true
end
end
using .mcssf1
mcssf1.test()

module mcmatopt1
using FinEtools
using FinEtoolsDeforNonlinear
using LinearAlgebra: Transpose, I, diagm, mul!, norm
using BenchmarkTools
using Test

function test()
	F = rand(3, 3)
	b = fill(0.0, 3, 3)
	mul!(b, F, Transpose(F))
	@test norm(b - F*F') < 1.0e-6
	true
end
end
using .mcmatopt1
mcmatopt1.test()

module mneotest1
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNonlinearModule: estimatesoundspeed
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.MatDeforNeohookeanNaiveModule: MatDeforNeohookeanNaive
using LinearAlgebra
using BenchmarkTools
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3

    m1 = MatDeforNeohookean(mr, E, nu)
    m2 = MatDeforNeohookeanNaive(mr, E, nu)
    # @show m
    update! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.update!

    stress = fill(0.0, 6)
    D = fill(0.0, 6, 6)
    output = FFlt[]
    statev = FFlt[]
    tn = 0.0
    dtn = 0.0
    loc = [0.0 0.0 0.0]
    label = 0
    quantity=:nothing
    Fn1 = [1.1 0 0; 0 1.0 0.02; 0.1 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]

    update!(m1, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    stress1 = deepcopy(stress)
    # @show stress
    update!(m2, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress - stress1) / norm(stress) <= 1.0e-12
    true
end
end
using .mneotest1
mneotest1.test()