text:
const Kron{T,N,I<:Tuple} = ApplyArray{T,N,typeof(kron),I}

function instantiate(K::Applied{<:Any,typeof(kron)})
    isempty(K.args) && throw(ArgumentError("Cannot take kronecker product of empty vectors"))
    applied(kron, map(instantiate,K.args)...)
end

Kron(A...) = ApplyArray(kron, A...)
Kron{T}(A...) where T = ApplyArray{T}(kron, A...)

_kron_dims() = 0
_kron_dims(A, B...) = max(ndims(A), _kron_dims(B...))

eltype(A::Applied{<:Any,typeof(kron)}) = promote_type(map(eltype,A.args)...)
ndims(A::Applied{<:Any,typeof(kron)}) = _kron_dims(A.args...)

size(K::Kron, j::Int) = prod(size.(K.args, j))
size(a::Kron{<:Any,1}) = (size(a,1),)
size(a::Kron{<:Any,2}) = (size(a,1), size(a,2))
axes(a::Kron{<:Any,1}) = (oneto(size(a,1)),)
axes(a::Kron{<:Any,2}) = (oneto(size(a,1)), oneto(size(a,2)))

issquare(M::AbstractMatrix) = (size(M, 1) == size(M, 2))

adjoint(K::Kron{<:Any, 2}) = Kron(adjoint.(K.args)...)
transpose(K::Kron{<:Any, 2}) = Kron(transpose.(K.args)...)
pinv(K::Kron{<:Number, 2}) = Kron(pinv.(K.args)...)
function inv(K::Kron{<:Number, 2})
    n = checksquare(K)

    # see det below for why the presence of rectangular factors renders
    # the kronecker product singular
    all(issquare, K.args) || throw(LinearAlgebra.SingularException(n))

    return Kron(inv.(K.args)...)
end

function det(K::Kron{T, 2}) where T<:Number
    s = checksquare(K)
    d = one(T)

    for A in K.args
        if issquare(A)
            dA = det(A)
            if iszero(dA)
                return dA
            end
            d *= dA^(s ÷ size(A, 1))
        else
            # The Kronecker Product of rectangular matrices, if it is square, will
            # have determinant zero. This can be shown by using the fact that
            # rank(A ⊗ B) = rank(A)rank(B) and showing that this is strictly less
            # than the number of rows/cols in the resulting Kronecker matrix. Hence,
            # since A ⊗ B does not have full rank, its determinant must be zero.
            return zero(T)
        end
    end
    return d
end

function logabsdet(K::Kron{T, 2}) where T<:Number
    s = checksquare(K)
    d = zero(T)
    sgn = one(T)

    for A in K.args
        if issquare(A)
            ldA, sgnA = logabsdet(A)
            if isinf(ldA)
                return ldA, sgnA
            end
            p = (s ÷ size(A, 1))
            d += p * ldA
            sgn *= sgnA ^ p
        else
            # see definition of det above for explanation
            return float(T)(-Inf), zero(T)
        end
    end
    return d, sgn
end

function tr(K::Kron{<:Any, 2})
    checksquare(K)
    if all(issquare, K.args)
        return prod(tr.(K.args))
    else
        return sum(diag(K))
    end
end

function diag(K::Kron{<:Any, 2})
    if all(issquare, K.args)
        return kron(diag.(K.args)...)
    else
        d = similar(K.args[1], minimum(size(K)))
        @inbounds for i in 1:length(d)
            d[i] = K[i, i]
        end
        return d
    end
end

_compatible_sizes((A, B)) = (size(A, 2) == size(B, 1))

# Implements the mixed-product property for the Kronecker product and matrix
# multiplication
function copy(M::Mul{ApplyLayout{typeof(kron)}, ApplyLayout{typeof(kron)}})
    A, B = M.A, M.B
    MT = promote_type(map(typeof, arguments(A))..., map(typeof, arguments(B))...)
    # Keeping it simple for now, but could potentially make this "alignment"-check
    # more precise. For example, if one factor of A has 6 columns and it's aligned
    # with two factors of B with 2 and 3 rows, then we can still use the
    # mixed-product property after explicitly computing the Kronecker product of
    # the two factors of B (or maybe just deferring to the shuffle algorithm)
    if (length(arguments(A)) == length(arguments(B))) && all(_compatible_sizes, zip(arguments(A), arguments(B)))
        factors = [(A_i * B_i) for (A_i, B_i) in zip(arguments(A), arguments(B))]
        return Kron(factors...)
    else
        algo_type = shuffle_algorithm_type(MT)
        return shuffle_algorithm(algo_type, A, B, eltype(M))
    end
end

# Implements Roth's lemma aka the "vec-trick" for multiplying a 2-factor Kron
# matrix to a vector
function copy(M::Mul{
    LazyArrays.ApplyLayout{typeof(kron)},
    <:Any,
    Kron{T,2,NTuple{2,M}},V
} where {T, M<:AbstractMatrix{T}, V<:AbstractVector})
    K, v = M.A, M.B
    A, B = arguments(K)
    V = reshape(v, size(B, 2), size(A, 2))
    return vec(B * V * transpose(A))
end

abstract type ShuffleAlgorithm end
struct Shuffle <: ShuffleAlgorithm end
struct ModifiedShuffle <: ShuffleAlgorithm end

shuffle_algorithm_type(::Type) = Shuffle()
shuffle_algorithm_type(::Type{<:AbstractArray}) = ModifiedShuffle()
shuffle_algorithm_type(::Type{<:ApplyArray}) = Shuffle()

function copy(M::Mul{ApplyLayout{typeof(kron)}})
    MT = promote_type(map(typeof, arguments(M.A))...)
    algo_type = shuffle_algorithm_type(MT)
    return shuffle_algorithm(algo_type, M.A, M.B, eltype(M))
end

function shuffle_algorithm(
    ::ModifiedShuffle, K::Kron{T,2} where T, p::AbstractVecOrMat, OT::Type{<:Number}
)

    if size(K, 2) != size(p, 1)
        if ndims(p) > 1
            throw(DimensionMismatch("matrix A has dimensions $(size(K)), matrix B has dimensions $(size(p))"))
        else
            throw(DimensionMismatch("matrix A has dimensions $(size(K)), vector B has length $(size(p, 1))"))
        end
    end

    output_sizes = (ndims(p) > 1) ? (size(K, 1), size(p, 2)) : (size(K,1),)
    q = fill!(similar(p, OT, output_sizes), zero(OT))
    if any(iszero, K.args)
        return q
    end

    K_shrunk_factors::Vector{eltype(K.args)} = []
    R_H::Vector{Vector{Int}} = []
    C_H::Vector{Vector{Int}} = []

    is_dense = !any(issparse, K.args)

    # note: the following computation costs are for multiplication against
    #       a single vector.
    nnz_ = [issparse(X_h) ? nnz(X_h) : prod(size(X_h)) for X_h in K.args]
    trad_cost = 2*sum([
        prod(size.(K.args[1:h-1], 2)) * nnz_[h] * prod(size.(K.args[h+1:end], 1))
        for (h, X_h) in enumerate(K.args)
        if X_h != I
    ])

    naive_cost = 2*prod(size(K))

    if is_dense && naive_cost <= trad_cost && naive_cost < 1e6
        return Matrix(K) * p
    end

    for X_h in K.args
        R_h, C_h = Set{Int}(), Set{Int}()

        for i in axes(X_h, 1), j in axes(X_h, 2)
            if X_h[i, j] != 0
                push!(R_h, i)
                push!(C_h, j)
            end
        end

        R_h, C_h = sort(collect(R_h)), sort(collect(C_h))

        is_dense = is_dense && (axes(X_h, 1) == R_h && axes(X_h, 2) == C_h)

        push!(K_shrunk_factors, X_h[R_h, C_h])
        push!(R_H, R_h)
        push!(C_H, C_h)
    end

    if is_dense
        # this means that the component matrices were all dense AND they had no
        # zero rows/cols
        return shuffle_algorithm(Shuffle(), K, p, OT)
    end

    nnz_m = [issparse(X_h) ? nnz(X_h) : prod(size(X_h)) for X_h in K_shrunk_factors]
    modified_cost = 2*sum([
        prod(length.(C_H[1:h-1])) * nnz_m[h] * prod(length.(R_H[h+1:end]))
        for (h, X_h) in enumerate(K.args)
        if X_h != I
    ])

    if trad_cost <= modified_cost
        return shuffle_algorithm(Shuffle(), K, p, OT)
    end

    R_indices = Iterators.product(reverse(R_H)...)
    C_indices = Iterators.product(reverse(C_H)...)

    r_index_v = [1, cumprod([size(m, 1) for m in K.args[end:-1:2]])...]
    c_index_v = [1, cumprod([size(m, 2) for m in K.args[end:-1:2]])...]

    r_indices = zeros(Int, prod(size(R_indices)))
    for (i, r_ind) in enumerate(R_indices)
        r_indices[i] = dot(r_index_v, r_ind .- 1) + 1
    end

    c_indices = zeros(Int, prod(size(C_indices)))
    for (i, c_ind) in enumerate(C_indices)
        c_indices[i] = dot(c_index_v, c_ind .- 1) + 1
    end

    K_ = Kron(K_shrunk_factors...)
    if ndims(p) == 1
        q[r_indices] = shuffle_algorithm(Shuffle(), K_, p[c_indices], OT)
    else
        q[r_indices, :] = shuffle_algorithm(Shuffle(), K_, p[c_indices, :], OT)
    end

    return q
end

function shuffle_algorithm(::Shuffle, K::Kron{T,2} where T, p::AbstractVecOrMat, OT::Type{<:Number})
    q = copy!(similar(p, OT), p)
    r = [size(m, 1) for m in K.args]
    c = [size(m, 2) for m in K.args]

    all_square = all(map(==, r, c))
    is_vector = (ndims(p) == 1)

    i_left, i_right = 1, prod(c)
    H = length(K.args)

    @inbounds for h in 1:H
        X_h = K.args[h]
        r_h, c_h = r[h], c[h]
        i_right ÷= c_h

        if X_h != I
            if all_square
                q′ = q
            else
                size_ = prod(r[1:h]) * prod(c[h+1:H])
                output_sizes = (ndims(p) > 1) ? (size_, size(p, 2)) : (size_,)
                q′ = fill!(similar(p, OT, output_sizes), zero(OT))
            end

            base_i, base_j = 0, 0
            for i_l in 1:i_left
                for i_r in 1:i_right
                    slc_in  = base_i + i_r : i_right : base_i + i_r + (c_h-1)*i_right
                    slc_out = base_j + i_r : i_right : base_j + i_r + (r_h-1)*i_right

                    if is_vector
                        @views q′[slc_out] = X_h * q[slc_in]
                    else
                        @views q′[slc_out, :] = X_h * q[slc_in, :]
                    end
                end

                base_i += i_right*c_h
                base_j += i_right*r_h
            end
            q = q′
        end

        i_left *= r_h
    end

    return q
end

kron_getindex((A,)::Tuple{AbstractVector}, k::Integer) = A[k]
function kron_getindex((A,B)::NTuple{2,AbstractVector}, k::Integer)
    K,κ = divrem(k-1, length(B))
    A[K+1]*B[κ+1]
end
kron_getindex((A,)::Tuple{AbstractMatrix}, k::Integer, j::Integer) = A[k,j]
function kron_getindex((A,B)::NTuple{2,AbstractArray}, k::Integer, j::Integer)
    K,κ = divrem(k-1, size(B,1))
    J,ξ = divrem(j-1, size(B,2))
    A[K+1,J+1]*B[κ+1,ξ+1]
end

kron_getindex(args::Tuple, k::Integer, j::Integer) = kron_getindex(tuple(Kron(args[1:2]...), args[3:end]...), k, j)
kron_getindex(args::Tuple, k::Integer) = kron_getindex(tuple(Kron(args[1:2]...), args[3:end]...), k)

getindex(K::Kron{<:Any,1}, k::Integer) = kron_getindex(K.args, k)
getindex(K::Kron{<:Any,2}, k::Integer, j::Integer) = kron_getindex(K.args, k, j)
getindex(K::Applied{DefaultArrayApplyStyle,typeof(kron)}, k::Integer) = kron_getindex(K.args, k)
getindex(K::Applied{DefaultArrayApplyStyle,typeof(kron)}, k::Integer, j::Integer) = kron_getindex(K.args, k, j)

## Adapted from julia/stdlib/LinearAlgebra/src/dense.jl kron definition
function _kron!(R, a, b)
    size(R) == size(Kron(a,b)) || throw(DimensionMismatch("Matrices have wrong dimensions"))
    require_one_based_indexing(a, b)
    m = 1
    @inbounds for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1)
        aij = a[i,j]
        for k = 1:size(b,1)
            R[m] = aij*b[k,l]
            m += 1
        end
    end
    R
end

_kron!(R, a) = copyto!(R, a)

_kron!(R, a, b, c, d...) = _copyto!(UnknownLayout(), UnknownLayout(), R, Kron(a,b,c,d...))

_copyto!(_, ::ApplyLayout{typeof(kron)}, R::AbstractMatrix, K::AbstractMatrix) =
    _kron!(R, arguments(K)...)
_copyto!(_, ::ApplyLayout{typeof(kron)}, R::AbstractVector, K::AbstractVector) =
    _kron!(R, arguments(K)...)

struct Diff{T, N, Arr} <: LazyArray{T, N}
    v::Arr
    dims::Int
end

Diff(v::AbstractVector{T}) where T = Diff{T,1,typeof(v)}(v, 1)

function Diff(A::AbstractMatrix{T}; dims::Integer) where T
    dims == 1 || dims == 2 || throw(ArgumentError("dimension must be 1 or 2, got $dims"))
    Diff{T,2,typeof(A)}(A, dims)
end

IndexStyle(::Type{<:Diff{<:Any,1}}) = IndexLinear()
IndexStyle(::Type{<:Diff{<:Any,2}}) = IndexCartesian()

size(D::Diff{<:Any,1}) = (length(D.v)-1,)
function size(D::Diff{<:Any,2})
    if D.dims == 1
        (size(D.v,1)-1,size(D.v,2))
    else #dims == 2
        (size(D.v,1),size(D.v,2)-1)
    end
end

getindex(D::Diff{<:Any, 1}, k::Integer) = D.v[k+1] - D.v[k]
function getindex(D::Diff, k::Integer, j::Integer)
    if D.dims == 1
        D.v[k+1,j] - D.v[k,j]
    else # dims == 2
        D.v[k,j+1] - D.v[k,j]
    end
end

mutable struct Accumulate{T, N, Op, DM<:AbstractArray{T,N}, Arr<:AbstractArray{T,N}} <: AbstractCachedArray{T, N}
    op::Op
    data::DM
    v::Arr
    dims::Int
    datasize::NTuple{N,Int}
end

Accumulate(op, data::AbstractArray, v::AbstractArray, d::Int) = Accumulate(op, data, v, d, size(data))
function Accumulate(op, v::AbstractVector, d::Int)
    @assert d == 1
    Accumulate(op, [v[1]], v, d)
end

function Accumulate(op, v::AbstractMatrix, d::Int)
    @assert 1 ≤ d ≤ 2
    data = d == 1 ? v[1:1,:] : v[:,1:1]
    Accumulate(op, data, v, d)
end
Accumulate(op, v::AbstractVector{T}; dims::Integer=1) where T = Accumulate(op, v, dims)

function Accumulate(op, A::AbstractMatrix{T}; dims::Integer=1) where T
    dims == 1 || dims == 2 || throw(ArgumentError("dimension must be 1 or 2, got $dims"))
    Accumulate(op, A, dims)
end

const Cumsum{T,N,Arr} = Accumulate{T,N,typeof(+),Array{T,N},Arr}
Cumsum(v::AbstractArray; dims::Integer=1) = Accumulate(+, v; dims=dims)

IndexStyle(::Type{<:Accumulate{<:Any,1}}) = IndexLinear()
IndexStyle(::Type{<:Accumulate{<:Any,2}}) = IndexCartesian()

size(Q::Accumulate) = size(Q.v)

copy(Q::Accumulate) = Accumulate(Q.op, copy(Q.data), copy(Q.v), Q.dims, Q.datasize)

==(a::Cumsum{<:Any,1}, b::Cumsum{<:Any,1}) = a.v == b.v

copyto!(x::AbstractArray{<:Any,N}, C::Cumsum{<:Any,N}) where N = cumsum!(x, C.v)

# How we populate the data
function cache_filldata!(K::Accumulate{<:Any,1}, inds)
    copyto!(view(K.data,inds), view(K.v,inds))
    @inbounds for k in inds
        K.data[k] = K.op(K.data[k-1], K.data[k])
    end
end

function cache_filldata!(K::Accumulate{<:Any,2}, kr, jr)
    if K.dims == 1
        @inbounds for j in jr, k in kr
            K.data[k,j] = K.op(K.data[k-1,j], K.v[k,j])
        end
    else # K.dims == 2
        @inbounds for j in jr, k in kr
            K.data[k,j] = K.op(K.data[k,j-1], K.v[k,j])
        end
    end
end

# keep lazy
cumsum(a::LazyArray; kwds...) = Cumsum(a; kwds...)
accumulate(op, a::LazyArray; kwds...) = Accumulate(op, a; kwds...)

const AccumulateAbstractVector{T,Op} = Accumulate{T,1,Op,Vector{T},AbstractVector{T}}
AccumulateAbstractVector(op, A::AbstractVector{T}) where T = AccumulateAbstractVector{T,typeof(op)}(op, [A[1]], A, 1, (1,))

## Rotations

for op in (:rot180, :rotl90, :rotr90)
    @eval begin
        ndims(::Applied{<:Any,typeof($op)}) = 2
        eltype(A::Applied{<:Any,typeof($op)}) = eltype(A.args...)
    end
end
size(A::Applied{<:Any,typeof(rot180)}) = size(A.args...)
axes(A::Applied{<:Any,typeof(rot180)}) = axes(A.args...)
size(A::Applied{<:Any,typeof(rotl90)}) = reverse(size(A.args...))
size(A::Applied{<:Any,typeof(rotr90)}) = reverse(size(A.args...))

getindex(A::Applied{<:Any,typeof(rot180)}, k::Int, j::Int) = A.args[1][end-k+1,end-j+1]
getindex(A::Applied{<:Any,typeof(rotl90)}, k::Int, j::Int) = A.args[1][j,end-k+1]
getindex(A::Applied{<:Any,typeof(rotr90)}, k::Int, j::Int) = A.args[1][end-j+1,k]

applylayout(::Type{typeof(rot180)}, ::AbstractStridedLayout) = StridedLayout()
function strides(A::ApplyMatrix{<:Any,typeof(rot180)})
    a,b = strides(A.args...)
    -a,-b
end
unsafe_convert(::Type{Ptr{T}}, A::ApplyMatrix{T,typeof(rot180)}) where T =
    pointer(A.args..., length(A))

applylayout(::Type{typeof(rot180)}, ::ApplyLayout{typeof(*)}) = ApplyLayout{typeof(*)}()
call(::ApplyLayout{typeof(*)}, A::ApplyMatrix{<:Any,typeof(rot180)}) = *
arguments(::ApplyLayout{typeof(*)}, A::ApplyMatrix{<:Any,typeof(rot180)}) = ApplyMatrix.(rot180, arguments(A.args...))


text:
using SafeTestsets

@safetestset "My f test" begin include("my_test.jl") end

ext: jl


text:
r, c = parse.(Int, split(readline()))
table = zeros(Int, r+1, c+1)
for row in 1:r
    table[row, 1:end-1] = parse.(Int, split(readline()))
end

for i in 1:r
    table[i,end] = sum(table[i,1:end-1])
end

for i in 1:c
    table[end,i] = sum(table[1:end-1, i])
end

table[end,end] = (sum(table[end,:]) + sum(table[:,end])) / 2

for i in 1:r+1
    for j in 1:c+1
        if j == 1
            print(table[i,j])
        else
            print(" ", table[i,j])
        end
    end
    println()
end


text:
# Autogenerated wrapper script for Nettle_jll for aarch64-linux-gnu
export libhogweed, libnettle, nettle_hash

using GMP_jll
JLLWrappers.@generate_wrapper_header("Nettle")
JLLWrappers.@declare_library_product(libhogweed, "libhogweed.so.6")
JLLWrappers.@declare_library_product(libnettle, "libnettle.so.8")
JLLWrappers.@declare_executable_product(nettle_hash)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libhogweed,
        "lib/libhogweed.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@init_executable_product(
        nettle_hash,
        "bin/nettle-hash",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module HW1

using Obfuscatee
using Nettle: hexdigest
import JSON

import DMUStudent

export
    fx,
    fy

# mostly to get binaries to recompile
version = v"1.0.11"

@binclude(".bin/hw1_4")

"""
    evaluate(f, [email]; fname="results.json")

Evaluate homework 1 programming assignment.

# Arguments
- `f::Function`: submitted function
- `email::String`: email of submitting student - must match gradescope email!

# Keyword Arguments
- `fname::String`: submission output file name
"""
evaluate = @binclude(".bin/hw1_eval")

end


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

@time myfunc()
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 1
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdata.csv"
    df = CSV.read(HTTP.get(url).body, DataFrame)
    #show(df, allcols=true)

    # create bus id variable
    df = @transform(df, bus_id = 1:size(df,1))

    #---------------------------------------------------
    # reshape from wide to long (must do this twice be-
    # cause DataFrames.stack() requires doing it one
    # variable at a time)
    #---------------------------------------------------
    # first reshape the decision variable
    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))

    # next reshape the odometer 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))

    # join reshaped df's back together
    df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
    sort!(df_long,[:bus_id,:time])

    flogit = glm(@formula(Y~Odometer + Odometer^2 + RouteUsage + RouteUsage^2+ Branded + time + time^2 +
                Odometer*RouteUsage + Odometer*RouteUsage^2 + Odometer*Branded + Odometer*time +
            Odometer*time^2 + Odometer^2*RouteUsage + Odometer^2*RouteUsage^2 + Odometer^2*Branded +
            Odometer^2*time + Odometer^2*time^2 + RouteUsage*Branded + RouteUsage*time + RouteUsage*time^2 +
            RouteUsage^2*Branded + RouteUsage^2*time + RouteUsage^2*time^2 + Branded*time + Branded*time^2 +
            Odometer*RouteUsage * Branded + Odometer*RouteUsage * time + Odometer*RouteUsage * time^2 +
            Odometer*RouteUsage^2 * Branded + Odometer*RouteUsage^2 * time + Odometer*RouteUsage^2 * time^2 +
            Odometer*Branded * time + Odometer*Branded * time^2 + Odometer^2*RouteUsage * Branded +
            Odometer^2*RouteUsage * time + Odometer^2*RouteUsage * time^2 + Odometer^2*RouteUsage^2 * Branded +
            Odometer^2*RouteUsage^2 * time + Odometer^2*RouteUsage^2 * time^2 + Odometer^2*Branded * time +
            Odometer^2*Branded * time^2 + RouteUsage*Branded * time + RouteUsage*Branded * time^2 +
            RouteUsage^2*Branded * time + RouteUsage^2*Branded * time^2 +
            Odometer*RouteUsage * Branded * time + Odometer*RouteUsage * Branded * time^2 +
            Odometer*RouteUsage^2 * Branded * time + Odometer*RouteUsage^2 * Branded * time^2 +
            Odometer^2*RouteUsage * Branded * time + Odometer^2*RouteUsage * Branded * time^2 +
            Odometer^2*RouteUsage^2 * Branded * time + Odometer^2*RouteUsage^2 * Branded * time^2),df_long, Binomial(), LogitLink())
    println(flogit)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 3
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    include("create_grids.jl")
    Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])

    # [ convert other data frame columns to matrices]
    O = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
    X = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
    Z = Matrix(df[:,[:Zst]])
    B = Matrix(df[:,[:Branded]])

    zval,zbin,xval,xbin,xtran = create_grids()
    T = 20
    FV = zeros((size(xtran,1),2,T+1))
    beta = 0.9
    theta = [0.1,0.2,0.3]
    sdf = DataFrame(Odometer = kron(ones(zbin),xval), RouteUsage = kron(ones(xbin),zval),Branded = zeros(xbin*size(zval,1)),time = zeros(xbin*size(zval,1)))

    function future_val(sdf,flogit,X, Z, B, xtran)
        for t = 2:T
            for i =0:1
                sdf.time = t * ones(size(sdf,1))
                sdf.Branded = i * ones(size(sdf,1))
                p0 = predict(flogit,sdf)
                FV[:,i+1,t] = -beta*log.(p0)
            end
        end

        FVT1 = zeros(size(X,1),T)
        for t = 2:T
            for i = 1:size(X,1)
                t = (T+1)-t
                z = Z[i]
                x = X[i,t]
                b = B[i]
                row1 = x + (z-1)*xbin
                row0 = 1+(z-1)*xbin
                FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,b+1,t+1]
            end
        end

        return FVT1'[:]

    end

    fvt1 = future_val(sdf,flogit,X, Z, B, xtran)
    df_long = @transform(df_long,fv = fvt1)
    theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),df_long, Binomial(), LogitLink(),offset=df_long.fv)
    println(theta_hat_ccp_glm)

end


text:
#= ######################################################
fct for struct
###################################################### =#

function dot(v1::Point2D{T}, v2::Point2D{T}) where T<:AbstractFloat
    d = v1.x * v2.x + v1.y * v2.y
    return d
end

#= ######################################################
inside functions
###################################################### =#

function get_angle(v1::Point2D{T},v2::Point2D{T}) where T<:AbstractFloat
    # calculate angle between two vectors
    EPSILON = 1E-8
    cosval = dot(v1,v2) / (norm(v1) * norm(v2))
    if cosval < -1 || cosval > 1
        if cosval > (-1 - EPSILON) && cosval < -1
            cosval = -1
        elseif cosval < (1 + EPSILON) && cosval > 1
            cosval = 1
        else
            error("custom error message: cosval in fct: get_angle")
        end
    end
    phir = acos(cosval)
    #phid = phir * 180 / pi
    return phir
end

function line_segment_intersection(p1::Point2D{T}, p2::Point2D{T}, p3::Point2D{T}, p4::Point2D{T}) where T<:AbstractFloat
    # check if two 2D line segments are intersecting
    hit = false
    x1 = p1.x
    x2 = p2.x
    x3 = p3.x
    x4 = p4.x
    y1 = p1.y
    y2 = p2.y
    y3 = p3.y
    y4 = p4.y
    denom = ((y2-y1)*(x3-x4)-(x2-x1)*(y3-y4))
    if denom != 0
        t = ((x1-x3)*(y3-y4)-(y1-y3)*(x3-x4))/denom
        u = ((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/denom
        if (t >= 0.0 && t <= 1.0) && (u >= 0.0 && u <= 1.0)
            hit = true
        end
    end
    return hit
end

function line_segment_intersection_test(p1::Point2D{T}, p2::Point2D{T}, p3::Point2D{T}, p4::Point2D{T}) where T<:AbstractFloat
    # check if two 2D line segments are intersecting
    #print("intersection check")
    x1 = p1.x
    x2 = p2.x
    x3 = p3.x
    x4 = p4.x
    y1 = p1.y
    y2 = p2.y
    y3 = p3.y
    y4 = p4.y
    denom = ((y2-y1)*(x3-x4)-(x2-x1)*(y3-y4))
    #print("denominator = ",denom)
    if denom == 0
        #print("no intersection -> parallel vectors")
        return false
    else
        t = ((x1-x3)*(y3-y4)-(y1-y3)*(x3-x4))/denom
        u = ((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/denom
        #print("t = ",t,"; u = ",u)
        if t >= 0.0 && t <= 1.0
            # pointX = x1+t*(x2-x1)
            # pointY = y1+t*(y2-y1)
            if u >= 0.0 && u <= 1.0
                # pointX = x3+u*(x4-x3)
                # pointY = y3+u*(y4-y3)
                # print("intersection at ",point)
                #plt.plot((x1,x2), (y1,y2))
                #plt.plot((x3,x4), (y3,y4))
                #plt.plot(point[0],point[1],'r+')
                return true
            else
                return false
            end
        else
            return false
        end
    end
end

function is_inside_tile_BB(n1::Point2D{T}, n2::Point2D{T}, tmin_x::T, tmax_x::T, tmin_y::T, tmax_y::T) where T<:AbstractFloat
    # check if element e is inside tile with bounding box
    box1 = @SMatrix [min(n1.x, n2.x) max(n1.x, n2.x); min(n1.y, n2.y) max(n1.y, n2.y)]
    box2 = @SMatrix [tmin_x tmax_x; tmin_y tmax_y]
    return isOverlapping2D(box1, box2)
end

function isOverlapping2D(box1,box2)
    # box1 = (x:(xmin1,xmax1),y:(ymin1,ymax1))
    # box2 = (x:(xmin2,xmax2),y:(ymin2,ymax2))
    # isOverlapping2D(box1,box2)
    # = isOverlapping1D(box1.x, box2.x) and
    #   isOverlapping1D(box1.y, box2.y)
    box1x = @view box1[1,:]
    box1y = @view box1[2,:]
    box2x = @view box2[1,:]
    box2y = @view box2[2,:]
    inside = false
    if isOverlapping1D(box1x, box2x) && isOverlapping1D(box1y, box2y)
        inside = true
    end
    return inside
end

function isOverlapping1D(box1,box2)
    # box1 = (xmin1, xmax1)
    # box2 = (xmin2, xmax2)
    # isOverlapping1D(box1,box2)
    # = xmax1 >= xmin2 and xmax2 >= xmin1
    xmin1 = box1[1]
    xmax1 = box1[2]
    xmin2 = box2[1]
    xmax2 = box2[2]
    inside = false
    if xmax1 >= xmin2 && xmax2 >= xmin1
        inside = true
    end
    return inside
end

function get_area_of_part(m, p)
    # returns the area of part p
    e1 = m.elem2par[p].first
    e2 = m.elem2par[p].last
    # println("calculate area between elements ", e1, " and ", e2)
    a = 0.0
    for i = e1:e2
        a += m.elem[i].area
    end
    return a
end

#= ######################################################
final functions
###################################################### =#

function existing_vf!(m::Model, mat)
    # check if elements are faced towards each other
    # angle between nvec and connector < 90°
    for i1 = 1:m.nelem, i2 = (i1+1):m.nelem
        con = m.elem[i2].com - m.elem[i1].com
        phir1 = get_angle(m.elem[i1].nvec,con)
        phir2 = get_angle(m.elem[i2].nvec,con*(-1))
        tol = 1E-8
        if phir1 < (pi/2 - tol) && phir2 < (pi/2 - tol)
            mat[i1,i2] = 1
            mat[i2,i1] = 1
        end
    end
end

function blocking_vf!(m::Model, mat)
    # check if element pairs are blocked by others with brute force method
    for i1 = 1:m.nelem, i2 = (i1+1):m.nelem
        if mat[i1,i2] == 1
            for ib in 1:m.nelem
                if ib != i1 && ib != i2
                    p1 = m.elem[i1].com
                    p2 = m.elem[i2].com
                    p3 = m.nodes[m.elem[ib].node1]
                    p4 = m.nodes[m.elem[ib].node2]
                    blocked = line_segment_intersection(p1, p2, p3, p4)
                    if blocked
                        # println(i1," and ",i2," intersected by ",ib)
                        mat[i1,i2] = 0
                        mat[i2,i1] = 0
                        break
                    end
                end
            end
        end
    end
end

function create_tiles(m::Model, n::Integer)::Tuple{Float64,Float64}
    # create tiles based on model nodes
    # get min and max x and y of nodes
    nmin, nmax = get_nodes_min_and_max(m)
    dx = (nmax.x - nmin.x) / n
    dy = (nmax.y - nmin.y) / n
    return dx, dy
end

function check_tile_occupation(m::Model, dx::T1, dy::T1, n::T2; blockparts = 1:m.npar)::Matrix{Union{Vector{T2},Missing}} where {T1<:AbstractFloat, T2<:Integer}
    # tile sorting
    t_occ = Matrix{Union{Vector{T2},Missing}}(missing,n,n)
    e_in_tile = Vector{T2}(undef,m.nelem) # empty vec for appending
    for t1 = 1:n
        for t2 = 1:n
            tmin_x = dx*(t1-1)
            tmax_x = dx*t1
            tmin_y = dy*(t2-1)
            tmax_y = dy*t2
            hit = 0
            for p in blockparts
                e1 = m.elem2par[p].first
                e2 = m.elem2par[p].last
                for i = e1:e2
                    node1 = m.nodes[m.elem[i].node1]
                    node2 = m.nodes[m.elem[i].node2]
                    if is_inside_tile_BB(node1, node2, tmin_x, tmax_x, tmin_y, tmax_y)
                        # println("tile (",t1,",",t2,") --> ",i)
                        hit += 1
                        e_in_tile[hit] = i
                    end
                end
            end
            if hit > 0
                t_occ[t1,t2] = e_in_tile[1:hit]
            end
        end
    end
    return t_occ
end

function blocking_vf_with_tiles!(m::Model, mat, dx::T1, dy::T1, n::T2, t_occ::Matrix{Union{Vector{T2},Missing}}) where {T1<:AbstractFloat, T2<:Integer}
    # check if element pairs are blocked by others with tiles
    max_steps = get_max_steps(n)
    tiles = Vector{Index2D{T2}}(undef,max_steps)
    @inbounds for i1 = 1:m.nelem, i2 = (i1+1):m.nelem
        # println("--------> checking ",i1," and ",i2)
        p1 = m.elem[i1].com
        p2 = m.elem[i2].com
        if mat[i1,i2] == 1
            # first tilewalk and get all tiles between i1 and i2
            ntiles = tilewalk_with_return!(tiles, p1, p2, dx, dy, n)
            hitten = false
            @inbounds for i = 1:ntiles
                t = tiles[i]
                # println("next tile number: ", t.x, ", ", t.y)
                if !ismissing(t_occ[t.x,t.y])
                    # println("------> occupied ")
                    for it in t_occ[t.x,t.y]
                        if it != i1 && it != i2
                            # println("------> ",i1," and ",i2, " checking with ",it)
                            p3 = m.nodes[m.elem[it].node1]
                            p4 = m.nodes[m.elem[it].node2]
                            blocked = line_segment_intersection(p1, p2, p3, p4)
                            if blocked
                                # println("----------------> ",i1," and ",i2," intersected by ",it)
                                mat[i1,i2] = 0
                                mat[i2,i1] = 0
                                hitten = true
                                break
                            end
                        end
                    end
                end
                if hitten
                    break
                end
            end
        end
    end
end

function blocking_vf_with_tiles_simplified!(m::Model, mat, dx::T1, dy::T1, n::T2, t_occ::Matrix{Union{Vector{T2},Missing}}; elem_in_t = 12::T2, skipped_t::T2 = 2) where {T1<:AbstractFloat, T2<:Integer}
    # check if element pairs are blocked by others with tiles
    max_steps = get_max_steps(n)
    tiles = Vector{Index2D{T2}}(undef,max_steps)
    count_simp = 0
    count_all = 0
    @inbounds for i1 = 1:m.nelem, i2 = (i1+1):m.nelem
        # println("--------> checking ",i1," and ",i2)
        p1 = m.elem[i1].com
        p2 = m.elem[i2].com
        if mat[i1,i2] == 1
            # first tilewalk and get all tiles between i1 and i2
            ntiles = tilewalk_with_return!(tiles, p1, p2, dx, dy, n)
            hitten = false
            @inbounds for i = 1:ntiles
                t = tiles[i]
                # println("next tile number: ", t.x, ", ", t.y)
                if !ismissing(t_occ[t.x,t.y])
                    # println("------> occupied ")
                    count_all += 1
                    if size(t_occ[t.x,t.y],1) <= elem_in_t || i <= 1+skipped_t || i >= ntiles-skipped_t
                        for it in t_occ[t.x,t.y]
                            if it != i1 && it != i2
                                # println("------> ",i1," and ",i2, " checking with ",it)
                                p3 = m.nodes[m.elem[it].node1]
                                p4 = m.nodes[m.elem[it].node2]
                                blocked = line_segment_intersection(p1, p2, p3, p4)
                                if blocked
                                    # println("----------------> ",i1," and ",i2," intersected by ",it)
                                    mat[i1,i2] = 0
                                    mat[i2,i1] = 0
                                    hitten = true
                                    break
                                end
                            end
                        end
                    else
                        # println("----------------> SIMPLIFIED -> ",i1," and ",i2," intersected by ",it)
                        mat[i1,i2] = 0
                        mat[i2,i1] = 0
                        count_simp += 1
                    end
                end
                if hitten
                    break
                end
            end
        end
    end
    # println("counter simp/all: ", count_simp, " / ", count_all)
    println("simplification rate: ", round(count_simp/count_all*100, digits = 2), " %")
end

function calculating_vf!(m, mat; normit = false)
    # calculate viewfactors if vf is existing
    for i1 = 1:m.nelem, i2 = (i1+1):m.nelem
        if mat[i1,i2] == 1
            # viewfactor existing
            a = m.nodes[m.elem[i1].node1]
            b = m.nodes[m.elem[i1].node2]
            c = m.nodes[m.elem[i2].node1]
            d = m.nodes[m.elem[i2].node2]
            # punkte müssen richtig sortiert sein
            # nur aus python rauskopiert - nicht mehr kontrolliert
            # hier umgangen mit betrag
            # 4 Strecken berechnen -> Abstand zwischen 2 Punkten
            ac = norm(c - a)
            bd = norm(d - b)
            bc = norm(c - b)
            ad = norm(d - a)
            # println("AC: ", ac," BD: ", bd," BC: ", bc," AD: ", ad)
            # cross strings formula - Fläche berechnen und vf bestimmen
            A1 = norm(b - a)
            A2 = norm(d - c)
            mat[i1,i2] = abs((0.5 * ((ac+bd)-(bc+ad))) / A1)
            mat[i2,i1] = abs((0.5 * ((ac+bd)-(bc+ad))) / A2)
            #print(i1," -> ",i2,":",vf12,vf21)
        end
    end
    if normit
        control = sum(mat, dims=2)
        mat .= mat ./ control
    end
end

#= ######################################################
blocking test with 2 elements
###################################################### =#

function existing_vf_2elem(m, i1, i2)
    # check if elements are faced towards each other
    # angle between nvec and connector < 90°
    exists = false
    con = m.elem[i2].com - m.elem[i1].com
    phir1 = get_angle(m.elem[i1].nvec,con)
    phir2 = get_angle(m.elem[i2].nvec,con*(-1))
    println("Angles: ", phir1, " and ", phir2)
    tol = 1E-8
    if phir1 < (pi/2 - tol) && phir2 < (pi/2 - tol)
        println("is existing between ", i1, " and ", i2)
        exists = true
    else
        println("is not existing between ", i1, " and ", i2)
    end
    return exists
end

function blocking_vf_with_tiles_2elem(fig, ax, i1, i2, m::Model, dx, dy, n::T2, t_occ) where T2<:Integer
    # check if element pair is blocked by others with tiles
    p1 = m.elem[i1].com
    p2 = m.elem[i2].com
    max_steps = get_max_steps(n)
    tiles = Vector{Index2D{T2}}(undef,max_steps)
    # first tilewalk and get all tiles between i1 and i2
    ntiles = tilewalk_with_return!(tiles, p1, p2, dx, dy, n)
    println(tiles[1:ntiles])
    hitten = false
    for i = 1:ntiles
        t = tiles[i]
        println("next tile number: ", t.x, ", ", t.y)
        # print tile in color
        poly!(ax, Rect((t.x-1) * dx, (t.y-1) * dy, dx, dy), color = (:grey, 0.6))
        if !ismissing(t_occ[t.x,t.y])
            println("------> occupied ")
            for it in t_occ[t.x,t.y]
                if it != i1 && it != i2
                    println("------> ",i1," and ",i2, " checking with ",it)
                    p3 = m.nodes[m.elem[it].node1]
                    p4 = m.nodes[m.elem[it].node2]
                    blocked = line_segment_intersection(p1, p2, p3, p4)
                    if blocked
                        println("----------------> ",i1," and ",i2," intersected by ",it)
                        hitten = true
                        break
                    end
                end
            end
        end
        if hitten
            break
        end
    end
end

function blocking_vf_with_tiles_2elem_tiles(fig, ax, i1, i2, m::Model, dx, dy, n::T2, t_occ) where T2<:Integer
    # check if element pair is blocked by others with tiles
    p1 = m.elem[i1].com
    p2 = m.elem[i2].com
    max_steps = get_max_steps(n)
    tiles = Vector{Index2D{T2}}(undef,max_steps)
    # first tilewalk and get all tiles between i1 and i2
    ntiles = tilewalk_with_return!(tiles, p1, p2, dx, dy, n)
    # println(tiles[1:ntiles])
    for i = 1:ntiles
        t = tiles[i]
        println("next tile number: ", t.x, ", ", t.y)
        # print tile in color
        poly!(ax, Rect((t.x-1) * dx, (t.y-1) * dy, dx, dy), color = (:grey, 0.6))
    end
end

#= ######################################################
post processing of viewfactors
###################################################### =#

function vfmat_to_parts(m, mat; normit = false)
    # get viewfactor matrix parts out of viewfactor matrix elements
    matp = zeros(Float64, m.npar, m.npar)
    areap = [get_area_of_part(m, i) for i = 1:m.npar]
    for p1 = 1:m.npar, p2 = 1:m.npar
        vfsplit = mat[m.elem2par[p1].first:m.elem2par[p1].last, m.elem2par[p2].first:m.elem2par[p2].last]
        acol = [m.elem[i].area for i = m.elem2par[p1].first:m.elem2par[p1].last]
        asplit = repeat(acol, 1, (m.elem2par[p2].last-m.elem2par[p2].first+1))
        vfsplit = vfsplit .* asplit
        matp[p1,p2] = sum(vfsplit) / areap[p1]
    end
    if normit
        controlp = sum(matp, dims=2)
        matp .= matp ./ controlp
    end
    return matp
end

function save_vfmat_part_to_csv(mat, path, name; control = false)
    # save vfmat parts to parts as csv
    # option of saving with control sum
    if control == true
        control = sum(mat, dims=2)
        mat = hcat(mat, control)
        nparts = size(mat,1)
        header = Matrix{String}(undef, 1, nparts+1)
        for i = 1:nparts
            header[1,i] = "part"*string(i)
        end
        header[1,end] = "control"
    else
        header = []
    end
    filename = path*"/"*name*".txt"
    open(filename, "w") do io
        writedlm(io, header)
        writedlm(io, mat)
    end
end

function get_number_of_elements_in_tiles(t_occ)
    nx = size(t_occ,1)
    ny = size(t_occ,2)
    t_num = zeros(Int64, nx, ny)
    for i = 1:nx, j = 1:ny
        if ismissing(t_occ[i,j])
            t_num[i,j] = 0
        else
            t_num[i,j] = size(t_occ[i,j],1)
        end
    end
    return t_num
end

function tile_occ_analysis(t_occ; printit = true)
    # analysis of tile occupation matrix
    t_num = get_number_of_elements_in_tiles(t_occ)
    t_num_v = vec(t_num)
    t_num_v_occ = t_num_v[t_num_v[:] .> 0]
    t_min = minimum(t_num_v_occ)
    t_max = maximum(t_num_v_occ)
    n_all = size(t_num_v,1)
    n_occ = size(t_num_v_occ,1)
    t_mean = sum(t_num_v_occ) / n_all
    t_mean = round(t_mean, digits=1)
    t_mean_occ = sum(t_num_v_occ) / n_occ
    t_mean_occ = round(t_mean_occ, digits=1)
    if printit
        println("Tile occupation with elements:")
        println("    Tiles occupied: ", n_occ, " / ", n_all)
        println("    Max: ", t_max, " // Min: ", t_min)
        println("    Mean: ", t_mean, " // Mean_occ: ", t_mean_occ)
    else
        return t_max, t_min, t_mean_occ, n_occ/n_all
    end
end


text:
struct TaskStats{T}
    rwd::T
    dpl::T
end

function no_reward(prob::TaskStats)
    p, g = prob.rwd, prob.dpl
    return (1-p)*g / (1-(1-p)*(1-g))
end

function no_reward(prob::TaskStats, n)
    @assert n >= 0
    p, g = prob.rwd, prob.dpl
    return n == 0 ? 1.0 : (1-p)*(g+(1-g)*no_reward(prob, n-1))
end

function no_more_reward(prob::TaskStats, args...)
    g = prob.dpl
    g+(1-g)*no_reward(prob, args...)
end

function high_sequence_probability(prob::TaskStats, res::AbstractArray{Bool})
    p, g = prob.rwd, prob.dpl
    n = length(res)
    r = count(res)
    o = n-r
    return p^r*(1-p)^o*(1-g)^(n-1)
end

function sequence_probability(prob::TaskStats, res::AbstractArray{Bool})
    p, g = prob.rwd, prob.dpl
    ls = findlast(res)
    ls === nothing && return no_reward(prob)
    ph = high_sequence_probability(prob, res[1:ls])
    return ph*no_more_reward(prob, length(res)-ls)
end

one_reward(prob::TaskStats) = prob.rwd*no_more_reward(prob)


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

function arrow_buttons_with_repeater(store::AbstractStore, get_state=Redux.get_state)
    if AppDemo.Repeaters.Repeater(store, s->get_state(s).repeaters["basic_repeater"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
    end
end


text:
#=
This is a Julia port of tiny_sha3/sha3.c, 19-Nov-11  Dr. Markku-Juhani O. Saarinen <mjos@iki.fi>
https://github.com/mjosaarinen/tiny_sha3, copyright (c) 2015 Markku-Juhani O. Saarinen, MIT Licensed

This Julia code is Copyright (c) 2019, William R. Burdick Jr, MIT Licensed

The MIT License (MIT)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

const keccakf_rndc = UInt64[0x0000000000000001, 0x0000000000008082, 0x800000000000808a, 0x8000000080008000, 0x000000000000808b, 0x0000000080000001, 0x8000000080008081, 0x8000000000008009, 0x000000000000008a, 0x0000000000000088, 0x0000000080008009, 0x000000008000000a, 0x000000008000808b, 0x800000000000008b, 0x8000000000008089, 0x8000000000008003, 0x8000000000008002, 0x8000000000000080, 0x000000000000800a, 0x800000008000000a, 0x8000000080008081, 0x8000000000008080, 0x0000000080000001, 0x8000000080008008]
const keccakf_rotc = UInt32[1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 2, 14, 27, 41, 56, 8, 25, 43, 62, 18, 39, 61, 20, 44]
const keccakf_piln = map((x)-> x + 1, UInt32[10, 7, 11, 17, 18, 3, 5, 16, 8, 21, 24, 4, 15, 23, 19, 13, 12, 2, 20, 14, 22, 9, 6, 1])
const KECCAKF_ROUNDS = 24

@inline rotl64(x, y) = (x << y) | (x >> (64 - y))

mutable struct Ctx
    q::Vector{UInt64}
    b::AbstractArray{UInt8}
    pt::UInt64 # location in b, starts at 1
    rsiz::UInt64
    mdlen::UInt64
    function Ctx(mdlen)
        q = zeros(UInt64, 25)
        new(q, reinterpret(UInt8, q), 1, 200 - 2 * mdlen, mdlen)
    end
end

macro reorderfrom(ctx)
    if reinterpret(UInt8, UInt64[3])[1] == 3
        :()
    else
        :(fromendian($(esc(ctx))))
    end
end

function fromendian(ctx::Ctx)
    offset = 1
    q = ctx.q
    b = ctx.b
    for i in 1:length(ctx.q)
        q[i] = b[offset] | b[offset + 1] << 8 |
          b[offset + 2] << 16 | b[offset + 3] << 24 |
          b[offset + 4] << 32 | b[offset + 5] << 40 |
          b[offset + 6] << 48 | b[offset + 7] << 56
    end
end

macro reorderto(ctx)
    if reinterpret(UInt8, UInt64[3])[1] == 3
        :()
    else
        :(toendian($(esc(ctx))))
    end
end

function toendian(ctx::Ctx)
    q = ctx.q
    b = ctx.b
    for i in 1:length(a)
        offset = (i - 1) * 8 + 1
        t = q[i]
        b[offset] = t & 0xFF
        b[offset + 1] = (t >> 8) & 0xFF
        b[offset + 2] = (t >> 16) & 0xFF
        b[offset + 3] = (t >> 24) & 0xFF
        b[offset + 4] = (t >> 32) & 0xFF
        b[offset + 5] = (t >> 40) & 0xFF
        b[offset + 6] = (t >> 48) & 0xFF
        b[offset + 7] = (t >> 56) & 0xFF
    end
end

function sha3_keccakf(ctx::Ctx)
    i = j = 0
    t::UInt64 = 0
    bc = UInt64[0, 0, 0, 0, 0]
    st = ctx.q

    @reorderfrom(ctx)
    for r in 1:KECCAKF_ROUNDS
        # Theta
        for i in 1:5
            bc[i] = st[i] ⊻ st[i + 5] ⊻ st[i + 10] ⊻ st[i + 15] ⊻ st[i + 20]
        end
        for i in 1:5
            t = bc[(i + 3) % 5 + 1] ⊻ rotl64(bc[i % 5 + 1], 1)
            for j in 1:5:25
                st[j + i - 1] ⊻= t
            end
        end
        # Rho Pi
        t = st[2]
        for i in 1:24
            j = keccakf_piln[i]
            bc[1] = st[j]
            st[j] = rotl64(t, keccakf_rotc[i])
            t = bc[1]
        end
        # Chi
        for j in 1:5:25
            for i in 1:5
                bc[i] = st[j + i - 1]
            end
            for i in 1:5
                st[j + i - 1] ⊻= ~bc[i % 5 + 1] & bc[(i + 1) % 5 + 1]
            end
        end
        # Iota
        st[1] ⊻= keccakf_rndc[r]
    end
    @reorderto(ctx)
end

# Initialize the context for SHA3
const sha3_init = Ctx # alias for sha3_init(mdlen)

# update state with more data
function sha3_update(c::Ctx, data::Vector{UInt8})
    j = c.pt
    for i in 1:length(data)
        c.b[j] ⊻= data[i]
        j += 1
        if j > c.rsiz
            sha3_keccakf(c)
            j = 1
        end
    end
    c.pt = j
end

# finalize and output a hash
function sha3_final(c::Ctx)
    c.b[c.pt] ⊻= 0x06
    c.b[c.rsiz] ⊻= 0x80
    sha3_keccakf(c)
    c.b[1:c.mdlen] # return first mdlen bytes of c.b
end

# compute a SHA-3 hash (md) of given byte length from "in"
function sha3(in::Vector{UInt8}, mdlen::Int64)
    sha3 = Ctx(mdlen)
    sha3_update(sha3, in)
    sha3_final(sha3)
end

# SHAKE128 and SHAKE256 extensible-output functionality

function shake_xof(c::Ctx)
    c.b[c.pt] ⊻= 0x1F
    c.b[c.rsiz] ⊻= 0x80
    sha3_keccakf(c)
    c.pt = 1
end

function shake_out(c::Ctx, out::Vector{UInt8})
    j = c.pt
    for i in 1:length(out)
        if j > c.rsiz
            sha3_keccakf(c)
            j = 1
        end
        out[i] = c.b[j]
        j += 1
    end
    c.pt = j
end

function shake(in::Vector{UInt8}, mdlen::Int64)
    buf = zeros(UInt8, 32)
    c = Ctx(mdlen)

    sha3_update(c, in)
    shake_xof(c)
    for j in 1:32:512
        shake_out(c, buf)
    end
    buf
end

# finalize and output a hash
function keccak_final(c::Ctx)
    c.b[c.pt] ⊻= 0x01
    c.b[c.rsiz] ⊻= 0x80
    sha3_keccakf(c)
    c.b[1:c.mdlen] # return first mdlen bytes of c.b
end

# compute a keccak hash (md) of given byte length from "in"
function keccak(in::Vector{UInt8}, mdlen::Int64)
    sha3 = Ctx(mdlen)
    sha3_update(sha3, in)
    keccak_final(sha3)
end

macro defsha3(bits)
    :($(Symbol("sha3_" * string(bits) * "_init"))() = sha3_init($(UInt(trunc(bits / 8)))))
end

@defsha3(224)
@defsha3(256)
@defsha3(384)
@defsha3(512)

sha3_224(in::Vector{UInt8}) = sha3(in, 28)
sha3_256(in::Vector{UInt8}) = sha3(in, 32)
sha3_384(in::Vector{UInt8}) = sha3(in, 48)
sha3_512(in::Vector{UInt8}) = sha3(in, 64)

keccak224(in::Vector{UInt8}) = keccak(in, 28)
keccak256(in::Vector{UInt8}) = keccak(in, 32)
keccak384(in::Vector{UInt8}) = keccak(in, 48)
keccak512(in::Vector{UInt8}) = keccak(in, 64)

function test()
    println("sha3_224 zero\n  expect: 6B4E03423667DBB73B6E15454F0EB1ABD4597F9A1B078E3F5B5A6BC7\n"
      * "  got:    " * uppercase(bytes2hex(sha3(UInt8[], Int(224/8)))))
    println("sha3_256 zero\n  expect: A7FFC6F8BF1ED76651C14756A061D662F580FF4DE43B49FA82D80A4B80F8434A\n"
      * "  got:    " * uppercase(bytes2hex(sha3(UInt8[], Int(256/8)))))
    println("shake128 zero\n  expect: 43E41B45A653F2A5C4492C1ADD544512DDA2529833462B71A41A45BE97290B6F\n"
      * "  got:    " * uppercase(bytes2hex(shake(UInt8[], 16))))
    println("shake256 zero\n  expect: AB0BAE316339894304E35877B0C28A9B1FD166C796B9CC258A064A8F57E27F2A\n"
      * "  got:    " * uppercase(bytes2hex(shake(UInt8[], 32))))
    println("shake128 1600-bit\n  expect: 44C9FB359FD56AC0A9A75A743CFF6862F17D7259AB075216C0699511643B6439\n"
      * "  got:    " * uppercase(bytes2hex(shake(fill(0xA3, UInt(1600 / 8)), 16))))
    println("shake256 1600-bit\n  expect: 6A1A9D7846436E4DCA5728B6F760EEF0CA92BF0BE5615E96959D767197A0BEEB\n"
      * "  got:    " * uppercase(bytes2hex(shake(fill(0xA3, UInt(1600 / 8)), 32))))
    println("keccak256(hello)\n  expect: 1c8aff950685c2ed4bc3174f3472287b56d9517b9c948127319a09a7a36deac8\n"
      * "  got:    " * lowercase(bytes2hex(keccak(collect(UInt8, "hello"), 32))))
    println("keccak256 zero\n  expect: c5d2460186f7233c927e7db2dcc703c0e500b653ca82273b7bfad8045d85a470\n"
      * "  got:    " * lowercase(bytes2hex(keccak256(UInt8[]))))
    println("keccak256 2000-bit\n  expect: ???\n"
      * "  got:    " * lowercase(bytes2hex(keccak256(fill(UInt8(0xEF), 2000)))))
end


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "Real SHIP Implementation" begin

##
using SymPy
using ACE, JuLIP, BenchmarkTools, LinearAlgebra, Test, Random, StaticArrays
using JuLIP
using JuLIP: evaluate, evaluate_d, evaluate!, evaluate_d!
using JuLIP.Testing

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N), 0
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

##

trans = PolyTransform(3, 1.0)
fcut = PolyCutoff2s(2, 0.5, 3.0)
B2 = SHIPBasis(SparseSHIP(2, 5, wL=1.0), trans, fcut)
B3 = SHIPBasis(SparseSHIP(3, 5, wL=1.0), trans, fcut)
B4 = SHIPBasis(SparseSHIP(4, 5, wL=1.0), trans, fcut)
B5 = SHIPBasis(SparseSHIP(5, 5, wL=1.0), trans, fcut)
BB = [B2, B3, B4, B5]

##
@info("Testing Correctness of the C->R Ship conversion")
for B in BB
   coeffs = randcoeffs(B)
   ship = SHIP(B, coeffs)
   @info("bo = $(bodyorder(B)); converting to RSHIP ...")
   rship = ACE.convertc2r(ship)

   tmp = ACE.alloc_temp(ship, 10)
   rtmp = ACE.alloc_temp(rship, 10)

   for nsamples = 1:30
      Rs, Zs, z0 = randR(10)
      Es = evaluate!(tmp, ship, Rs, Zs, z0)
      rEs = evaluate!(rtmp, rship, Rs, Zs, z0)
      print_tf(@test Es ≈ rEs)
   end
   println()
end

##

@info("Check Correctness of RSHIP gradients")
for B in BB
   @info("   body-order = $(ACE.bodyorder(B))")
   coeffs = randcoeffs(B)
   ship = SHIP(B, coeffs)
   @info("   converting to RSHIP ...")
   rship = ACE.convertc2r(ship)

   Rs, Zs, z0 = randR(10)
   tmp = ACE.alloc_temp_d(rship, length(Rs))
   dEs = zeros(JVecF, length(Rs))
   evaluate_d!(dEs, tmp, rship, Rs, Zs, z0)
   Es = evaluate!(tmp, rship, Rs, Zs, z0)
   println(@test Es ≈ evaluate(rship, Rs, Zs, z0))
   println(@test dEs ≈ evaluate_d(rship, Rs, Zs, z0))
   @info("      Correctness of directional derivatives")
   for ndir = 1:20
      U = [rand(JVecF) .- 0.5 for _=1:length(Rs)]
      errs = Float64[]
      for p = 2:10
         h = 0.1^p
         dEs_U = dot(dEs, U)
         dEs_h = (evaluate!(tmp, rship, Rs + h * U, Zs, z0) - Es) / h
         push!(errs, abs(dEs_h - dEs_U))
      end
      success = (minimum(errs[2:end]) < 1e-3 * maximum(errs[1:3])) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end
   println()
end

##

# @info("Check Correctness of SHIP calculators")
#
# naive_energy(ship::SHIP, at) =
#       sum( evaluate(ship, R, at.Z[j], at.Z[i])
#             for (i, j, R) in sites(at, cutoff(ship)) )
#
# for B in BB
#    @info("   body-order = $(ACE.bodyorder(B))")
#    coeffs = randcoeffs(B)
#    ship = SHIP(B, coeffs)
#    at = bulk(:Si) * 3
#    at.Z[:] .= 0
#    rattle!(at, 0.1)
#    print("     energy: ")
#    println(@test energy(ship, at) ≈ naive_energy(ship, at) )
#    print("site-energy: ")
#    println(@test energy(ship, at) ≈ sum( site_energy(ship, at, n)
#                                          for n = 1:length(at) ) )
#    println("forces: ")
#    println(@test JuLIP.Testing.fdtest(ship, at))
#    println("site-forces: ")
#    println(@test JuLIP.Testing.fdtest( x -> site_energy(ship, set_dofs!(at, x), 3),
#                                        x -> mat(site_energy_d(ship, set_dofs!(at, x), 3))[:],
#                                        dofs(at) ) )
# end
#

end


text:
push!(LOAD_PATH,"../src/")

using Pkg
Pkg.activate(".")

_required = ["Weave", "Plots"]

Pkg.add.(_required)

using Mangal
using Weave

_path_elements = ["docs", "src", "vignettes"]
_jmd_files = filter(x -> endswith(x, ".Jmd"), readdir(joinpath(_path_elements...)))
_files_to_compile = [joinpath(_path_elements..., f) for f in _jmd_files]

for _file in _files_to_compile
    weave(_file, doctype="github")
end

Pkg.rm.(_required)


text:
const totalname = "Total"
const cellsname = "Cells"
const errorname = "Error"
const remaindername = "Remainder"

"""
    anova(observations::Array{Union{Number, Vector{Number}}}, factortypes = FactorType[]; factornames = String[], hasreplicates = true)
    anova(observations::Vector{Number}, factorassignments::Vector{Vector{Any}}, factortypes = FactorType[]; factornames = String[], hasreplicates = true)
    anova(df::DataFrame, observationscolumn::Symbol, factorcolumns::Vector{Symbol}, factortypes = FactorType[]; factornames = String[])

Performs an Analysis of Variance (ANOVA) calculation.

Operates on fixed or random factors, subject/block factors, and nested factors, with or without replicates, on balanced data.

Limitations:
- If any factors are `random` type, limited to 3-way
- Repeated measures ANOVA (with subject/block factor) limited to 3 `fixed` factors which may be partitioned as within or among subjects

Operates on up to 3 crossed factors (fixed or random) and arbitrarily many random nested factors, with or without
replicates, on balanced data.

# Arguments
- `observations`: Array containing the values to test. For the array, each dimension is a factor level, such that observations[2,5,3] indicates the 2nd level of the first factor, the 5th level of the second factor, and the 3rd level of the third factor. May contain values or vectors of values, where the vector contains replicates. Factors should be ordered with least significant first. For the vector, must provide `factorassignments` to specify factor levels.
- `factorassignments`: Vector of vectors of integers specifying how each observation is assigned to a factor level. Provide this when `observations` is given as a vector. Factor levels do not have to be consecutive or ordered. Nested factors must reuse factor levels currently.
- `factortypes`: Vector indicating the `FactorType` for each factor. Factors must be ordered with `nested` first, then `random`/`measured` or `fixed`/`manipulated` in any order. For repeated measures, `subject` or `block` must appear after within-subject `fixed` and before among-subject `fixed`. If too few values are provided, remaining are assumed to be `fixed`.
- `factornames`: Vector of names for each factor, excluding the replicate factor. If empty, will be automatically populated alphabetically.
- `hasreplicates`: Boolean to specify if the first level should be considered

Notes: The last index will be the top factor in the table.

Output: `AnovaData` structure containing the test results for each factor.

# Examples
```julia
anova(observations)                           # N-way fixed-effects ANOVA with replicates (vectors or first dimension)
anova(observations, hasreplicates = false)    # N-way fixed-effects ANOVA without replicates (first dimension)
anova(observations, [random])                 # N-way ANOVA with lower random factor and 1 or 2 upper fixed factors
anova(observations, [random])                 # N-way ANOVA with lower random factor and 1 or 2 upper fixed factors
anova(observations, [fixed, random])          # N-way ANOVA with 1 lower fixed factor, 1 random factor, and 0 or 1 upper fixed factor
anova(observations, [nested, random])         # N-way fixed-effects ANOVA with 1 random nested factor, 1 random factor, and 1-2 fixed factors
anova(observations, [fixed, subject])         # N-way repeated measures ANOVA with 1 within-subjects fixed factor
anova(observations, [fixed, block])           # N-way repeated measures ANOVA with 1 within-block fixed factor
anova(observations, [nested, fixed, subject]) # N-way repeated measures ANOVA with 1 within-subjects fixed factor and 1 nested factor
```

# Glossary
- observation: The dependent variable.
- factor: An independent variable.
- factor level: A value of a factor.
- balanced: All combinations of factor levels have the same number of observations.
- crossed factor: A factor with levels that combine with the levels of all other crossed factors.
- fixed/manipulated factor: A factor with fixed effects (e.g. treatment, concentration, exposure time).
- random/measured factor: A factor with random effects (e.g. location, individual).
- nested factor: A random factor where the levels are unique to a combination of crossed factor levels (e.g. replicate).
- subject/block factor: A nested factor that is subjected to multiple levels of another factor.
- sum of squares (SS): A measure of variance that is dependent on sample size. Also called "sum of squared deviations."
- degrees of freedom (DF, ν): The number of bins in which the values could have been moved, if random.
- mean square (MS): SS / DF. Corrects for the larger variance expected if random values can be assigned to more bins. Also called "mean squared error" or "mean squared deviation."
- F-statistic: The division of MS values produce a result belonging to the "F distribution", the shape of which depends on the DF of the numerator and error. The location of this value on the distribution provides the p-value.
- p-value: The probability that, if all measurements had been drawn from the same population, you would obtain data at least as extreme as contained in your observations.
- effect size (ω²): The standardized difference in the measurement caused by the factor. This is the generalized version which is largely stable with respect to study design.
"""
function anova(observations::AbstractVector{T}, factorassignments::AbstractVector{<:AbstractVector}, factortypes::Vector{FactorType} = FactorType[]; factornames::Vector{<:AbstractString} = String[]) where {T <: Number}
    # convert vector arguments into a multidimensional matrix

    length(observations) > 0 || return
    nfactors = length(factorassignments)
    N = length(observations)
    all(length.(factorassignments) .== N) || error("Each observation must have an assignment for each factor.")

    factorassignments = categorical.(factorassignments)
    droplevels!.(factorassignments)

    factorlevels = levels.(factorassignments)
    nfactorlevels = length.(factorlevels)
    N % prod(nfactorlevels) == 0 || error("Design is unbalanced.")

    factorlevelcounts = [[count(l -> l == factorlevels[i][j], factorassignments[i]) for j ∈ 1:nfactorlevels[i]] for i ∈ 1:nfactors]
    nperfactorlevel = factorlevelcounts .|> unique
    all(nperfactorlevel .|> length .== 1) || error("Design is unbalanced.")
    nperfactorlevel = nperfactorlevel .|> first

    nreplicates = Int(N / prod(nfactorlevels))
    hasreplicates = nreplicates > 1
    if hasreplicates
        nlevels = [nreplicates; nfactorlevels]
        sortorder = sortperm(repeat(1:nreplicates, Int(N / nreplicates)) .+
                             sum([levelcode.(factorassignments[i]) .* prod(nlevels[1:i]) for i ∈ 1:nfactors]))
    else
        nlevels = nfactorlevels
        sortorder = sortperm(sum([levelcode.(factorassignments[i]) .* prod(nlevels[1:i - 1]) for i ∈ 1:nfactors]))
    end

    observationsmatrix = reshape(observations[sortorder], nlevels...)
    anova(observationsmatrix, factortypes, factornames = factornames, hasreplicates = nreplicates > 1)
end

function anova(observations::AbstractArray{T}, factortypes::Vector{FactorType} = FactorType[]; factornames::Vector{<:AbstractString} = String[], hasreplicates = true) where {T <: Union{Number, AbstractVector{<:Number}}}
    length(observations) > 0 || return

    isrepeatedmeasures = subject ∈ factortypes

    eltype(observations) <: Number ||
        all(c -> length(c) == length(observations[1]), observations) ||
        error("All cells must have the same number of replicates.")

    if eltype(observations) <: AbstractVector && length(observations[1]) == 1
        observations = first.(observations)
        hasreplicates = false
    else
        observations = observations |> upcat
    end

    observations = observations .|> Float64

    nfactors = ndims(observations) - (hasreplicates ? 1 : 0)

    # defaults to assuming all unspecified factors are fixed
    if length(factortypes) < nfactors
        nremaining = nfactors - length(factortypes)
        append!(factortypes, fill(fixed, nremaining))
    end
    replace!(factortypes, block => subject)
    replace!(factortypes, manipulated => fixed)
    replace!(factortypes, measured => random)

    # automatically assigns alphabetical names if not provided.
    if isempty(factornames)
        factornames = string.('A':'Z')[1:min(nfactors, 26)]
        reverse!(factornames)
    end

    validate(observations, factortypes, factornames, nfactors, hasreplicates)

    anovakernel(observations, factornames, factortypes, hasreplicates)
end

function validate(observations, factortypes::Vector{FactorType}, factornames::Vector{<:AbstractString}, nfactors, hasreplicates)
    length(factortypes) == nfactors ||
            error("`factortypes` must have an entry for each factor.")

    nested ∉ factortypes ||
        factortypes[1:count(isnested, factortypes)] |> unique |> length == 1 ||
        error("Nested factors must come before any other factors.")

    if subject ∈ factortypes
        count(issubject, factortypes) == 1 ||
            error("Maximum of one subject/block factor.")

        notnestedfactortypes = filter(f -> !isnested(f), factortypes)
        subject ∈ notnestedfactortypes[2:min(4, end)] ||
            error("Subject/block factor must be in the second through fourth factors.")

        length(notnestedfactortypes) < 5 ||
            error("Maximum of 3 within-subjects or among-subjects factors.")

       # random ∉ notnestedfactortypes ||
        #    error("Random factors are not supported with subject/block factor.") # change to only filter out random within factors
    end

    crossedfactortypes = filter(f -> f ∈ [fixed, random], factortypes)
    length(crossedfactortypes) < 4 ||
        all(isfixed, crossedfactortypes) ||
        error("ANOVA with 4 or more crossed factors is not supported if any are random.")

    length(factornames) == nfactors || error("factornames must have an entry for each factor.")
end

#=
function anova(data::AnovaData, crossedfactors::Vector{Int}, )
    # performs a subtest of the specified crossed factors within level of the remaing crossed factors, using the original errors
end

# possible bug: function anova(data::AnovaData, crossedfactors::Vector{Int}, ) with normal functions ==> hang when called?
=#

function anovakernel(observations::AbstractArray{<:Number}, factornames, factortypes, hasreplicates)
    isrepeatedmeasures = subject ∈ factortypes

    N = length(observations)
    totaldf = N - 1
    totalvar = anovavalue(totalname, var(observations), totaldf)

    if hasreplicates
        nreplicates = size(observations, 1)
        cellmeans = meanfirstdim(observations)
    else
        nreplicates = 1
        cellmeans = observations
    end

    cellsdf = length(cellmeans) - 1
    cellsvar = anovavalue(cellsname, var(cellmeans) * nreplicates, cellsdf)

    errorvar = AnovaFactor(errorname, totalvar - cellsvar)

    nnested = count(isnested, factortypes)
    nestedfactornames = @view factornames[1:nnested]
    nnestedfactorlevels = size(cellmeans)[1:nnested]
    notnestedfactornames = @view factornames[(nnested + 1):end]
    notnestedfactortypes = @view factortypes[(nnested + 1):end]

    nestedvars, nestederrorvars, cellmeans, nnestedreplicates = nestedfactors(cellmeans, nreplicates, nestedfactornames, errorvar)

    # I need to check myself, but I'm 75% sure that without any within-subject factors, "subjects" is just like a nested factor
    # and thus this can be further simplified.
    if isrepeatedmeasures
        si = findfirst(x -> x == subject, notnestedfactortypes)
        factorvars, factorvarsdict = anovafactors(dropdims(mean(cellmeans, dims = si), dims = si), nnestedreplicates * size(cellmeans, si), notnestedfactornames[Not(si)])
        subjectinteractionvarsdict = subjectinteractions(cellmeans, si, nnestedreplicates, notnestedfactornames)
        # adjust keys of factorvarsdict to adjust for missing subject factor
        nfactors = ndims(cellmeans)
        si = invertfactorindex(si, nfactors)
        dictkeys = collect(keys(factorvarsdict))
        for k ∈ dictkeys
            knew = ntuple(i -> k[i] + (k[i] < si ? 0 : 1), length(k))
            factorvarsdict[knew] = factorvarsdict[k]
        end
        merge!(factorvarsdict, subjectinteractionvarsdict)
        factorerrorvars = anovasubjecterrors(factorvarsdict, notnestedfactortypes)
    else
        factorvars, factorvarsdict = anovafactors(cellmeans, nnestedreplicates, notnestedfactornames)
        crossedfactorerrorvar = nnested > 0 ?   nestedvars[1] :
                                hasreplicates ? errorvar :
                                                factorvars[end]
        factorerrorvars = anovaerrors(factorvarsdict, notnestedfactortypes, crossedfactorerrorvar)
    end

    allfactorvars = [factorvars; nestedvars]
    allfactorerrorvars = [factorerrorvars; nestederrorvars]

    if !hasreplicates
        if isrepeatedmeasures
            remaindervar = factorerrorvars[end]
        else
            remaindervar = pop!(allfactorvars)
            pop!(allfactorerrorvars)
            if nnested == 0
                pop!(factorvars)
                pop!(factorerrorvars)
            else
                pop!(nestedvars)
                pop!(nestederrorvars)
            end
        end
        errorvar = AnovaFactor(remaindername, remaindervar)
    end

    mainfactortypes = reverse(notnestedfactortypes[.!issubject.(notnestedfactortypes)])
    subjectfactors = isrepeatedmeasures ? collect(values(subjectinteractionvarsdict)) : AnovaEffect[]
    filter!(x -> x.ms != errorvar.ms, subjectfactors)
    ω² = effectsizeω²(factorvars, factorerrorvars, nestedvars, nestederrorvars, N, subjectfactors, errorvar, mainfactortypes)

    factorresults = AnovaResult.(ftest.(allfactorvars, allfactorerrorvars), ω²)

    allresults = [totalvar; factorresults; subjectfactors; errorvar]

    nfactorlevels = Int[size(cellmeans)...]
    return AnovaData(allresults, factorerrorvars, nfactorlevels, nreplicates, cellmeans)
end

function anovafactors(cellmeans, nreplicates, factornames)
    # calculate all simple factors and all interactions

    N = length(cellmeans) * nreplicates

    nfactors = ndims(cellmeans)
    factors = 1:nfactors

    allfactors = []
    allfactorvars = AnovaFactor[]
    allfactorsdict = Dict{Tuple{Vararg{Int}}, AnovaEffect}()
    for i ∈ factors
        ifactors = collect(combinations(reverse(factors), i))
        iotherfactors = [factors[Not(i...)] for i ∈ ifactors]
        iupperfactorvars = [allfactorvars[findall(x -> x ⊆ i, allfactors)] for i ∈ ifactors]

        ifactornames = [makefactorname(factornames[i]) for i ∈ ifactors]
        ifactorss = [var(mean(cellmeans, dims = iotherfactors[i]), corrected = false) * N - sum(iupperfactorvars[i]).ss for i ∈ eachindex(iotherfactors)]
        ifactordf = isempty(iupperfactorvars[1]) ? [size(cellmeans, invertfactorindex(i, nfactors)) - 1 for i ∈ factors] :
                                                   [prod(f.df for f ∈ iupperfactorvars[j][1:i]) for j ∈ eachindex(iotherfactors)]
        ifactorvars = AnovaFactor.(ifactornames, ifactorss, ifactordf)

        append!(allfactors, ifactors)
        append!(allfactorvars, ifactorvars)
        merge!(allfactorsdict, Dict(invertfactorindex.((tuple.(i...) for i ∈ ifactors), nfactors) .=> ifactorvars))
    end

    return allfactorvars, allfactorsdict
end

function anovaerrors(factorvarsdict, factortypes, errorvar)
    # assign proper error terms for each factor

    length(factortypes) == 1 && return [errorvar]
    all(isfixed, factortypes) && return fill(errorvar, length(factorvarsdict))

    factortypes = reverse(factortypes)

    if length(factortypes) == 2
        if factortypes[1] == factortypes[2]
            factor1error = factor2error = factorvarsdict[(1,2)]
        else
            if factortypes[1] == fixed
                factor1error = factorvarsdict[(1,2)]
                factor2error = errorvar
            else
                factor1error = errorvar
                factor2error = factorvarsdict[(1,2)]
            end
        end
        factorerrorvars = [factor1error; factor2error; errorvar]

    elseif length(factortypes) == 3
        if factortypes[1] == factortypes[2] == factortypes[3]
            factor1error = threeway_random_error(factorvarsdict[(1,2)], factorvarsdict[(1,3)], factorvarsdict[(1,2,3)])
            factor2error = threeway_random_error(factorvarsdict[(1,2)], factorvarsdict[(2,3)], factorvarsdict[(1,2,3)])
            factor3error = threeway_random_error(factorvarsdict[(1,3)], factorvarsdict[(2,3)], factorvarsdict[(1,2,3)])
            interaction12error = interaction13error = interaction23error = factorvarsdict[(1,2,3)]

        elseif factortypes[1] == factortypes[2]
            if factortypes[1] == fixed
                factor1error = factorvarsdict[(1,3)]
                factor2error = factorvarsdict[(2,3)]
                factor3error = errorvar
                interaction12error = factorvarsdict[(1,2,3)]
                interaction13error = interaction23error = errorvar

            else
                factor1error = factor2error = factorvarsdict[(1,2)]
                factor3error = threeway_random_error(factorvarsdict[(1,3)], factorvarsdict[(2,3)], factorvarsdict[(1,2,3)])
                interaction12error = errorvar
                interaction13error = interaction23error = factorvarsdict[(1,2,3)]

            end
        elseif factortypes[1] == factortypes[3]
            if factortypes[1] == fixed
                factor1error = factorvarsdict[(1,2)]
                factor2error = errorvar
                factor3error = factorvarsdict[(2,3)]
                interaction12error = interaction23error = errorvar
                interaction13error = factorvarsdict[(1,2,3)]

            else
                factor1error = factor3error = factorvarsdict[(1,3)]
                factor2error = threeway_random_error(factorvarsdict[(1,2)], factorvarsdict[(2,3)], factorvarsdict[(1,2,3)])
                interaction12error = interaction23error = factorvarsdict[(1,2,3)]
                interaction13error = errorvar

            end
        else
            if factortypes[2] == fixed
                factor1error = errorvar
                factor2error = factorvarsdict[(1,2)]
                factor3error = factorvarsdict[(1,3)]
                interaction12error = interaction13error = errorvar
                interaction23error = factorvarsdict[(1,2,3)]

            else
                factor1error = threeway_random_error(factorvarsdict[(1,2)], factorvarsdict[(1,3)], factorvarsdict[(1,2,3)])
                factor2error = factor3error = factorvarsdict[(2,3)]
                interaction12error = interaction13error = factorvarsdict[(1,2,3)]
                interaction23error = errorvar

            end
        end
        factorerrorvars = [factor1error; factor2error; factor3error; interaction12error; interaction13error; interaction23error; errorvar]

    else
        error("More than 3 factors with any random are not supported.")
    end

    return factorerrorvars
end

function anovasubjecterrors(factorvarsdict, factortypes)
    factortypes = reverse(factortypes)
    si = findfirst(x -> x == subject, factortypes)

    nfactors = length(factortypes)

    if nfactors == 2
        if si == 2
            # one among-subject factor
            factorerrorvars = [factorvarsdict[(1, si)]]
        else
            # one within-subject factor
            factorerrorvars = [factorvarsdict[(si, 2)]]
        end

    elseif nfactors == 3
        if si == 3
            # two among-subject factors
            factorerrorvars = [factorvarsdict[(1, si)], factorerrorvars[(2, si)],
                               factorerrorvars[(1, 2, si)]]
        elseif si == 2
            # one among-subject factor and one within-subject factor
            factorerrorvars = [factorvarsdict[(1, si)]; factorvarsdict[(1, si, 3)];
                               factorvarsdict[(1, si, 3)]]
        else
            # two within-subject factors
            factorerrorvars = [factorvarsdict[(si, 2)]; factorvarsdict[(si, 3)];
                               factorvarsdict[(si, 2, 3)]]
        end
    elseif nfactors == 4
        if si == 4
            # three among factors
            factorerrorvars = [factorvarsdict[(1, 2, 3, si)]; factorvarsdict[(1, 2, 3, si)]; factorvarsdict[(1, 2, 3, si)];
                               factorvarsdict[(1, 2, 3, si)]; factorvarsdict[(1, 2, 3, si)]; factorvarsdict[(1, 2, 3, si)];
                               factorvarsdict[(1, 2, 3, si)]]
        elseif si == 3
            # two among factors and one within-subject factor
            factorerrorvars = [factorvarsdict[(1, 2, si)]; factorvarsdict[(1, 2, si)];
                               factorvarsdict[(1, 2, si)];
                               factorvarsdict[(1, 2, si, 4)]; factorvarsdict[(1, 2, si, 4)]; factorvarsdict[(1, 2, si, 4)];
                               factorvarsdict[(1, 2, si, 4)]]
        elseif si == 2
            # one among factor and two within-subject factors
            factorerrorvars = [factorvarsdict[(1, si)];
                               factorvarsdict[(1, si, 3)]; factorvarsdict[(1, si, 4)];
                               factorvarsdict[(1, si, 3)]; factorvarsdict[(1, si, 4)]; factorvarsdict[(1, si, 3, 4)];
                               factorvarsdict[(1, si, 3, 4)]]
        else
            # three within-subject factors
            factorerrorvars = [factorvarsdict[(si, 2)]; factorvarsdict[(si, 3)]; factorvarsdict[(si, 4)];
                               factorvarsdict[(si, 2, 3)]; factorvarsdict[(si, 2, 4)]; factorvarsdict[(si, 3, 4)];
                               factorvarsdict[(si, 2, 3, 4)]]
        end
    else
        error("More than 3 non-subject factors are not supported.")
    end

    return AnovaFactor.(factorerrorvars)
end

function nestedfactors(cellmeans, nreplicates, factornames, errorvar)
    # compute nested factors and collapse nested levels

    nnested = length(factornames)

    nestedvars = AnovaFactor[]
    for i ∈ 1:nnested
        nestedmeans = mean(cellmeans, dims = 1)
        ss = sum((cellmeans .- nestedmeans) .^ 2) * nreplicates
        df = (size(cellmeans, 1) - 1) * prod(size(nestedmeans))
        push!(nestedvars, AnovaFactor(factornames[i], ss, df))
        nreplicates *= size(cellmeans, 1)
        cellmeans = dropdims(nestedmeans, dims = 1)
    end
    nestederrorvars = nnested > 0 ? [errorvar; nestedvars[1:(end - 1)]] :
                                    []
    return reverse!(nestedvars), reverse!(nestederrorvars), cellmeans, nreplicates
end

function subjectinteractions(cellmeans, si, nreplicates, factornames)
    cellmeans = reshape(cellmeans, (size(cellmeans)[1:si]..., prod(size(cellmeans)[(si + 1):end])))
    namongfactors = length(factornames) - si

    factorvarsdict = Dict{Tuple{Vararg{Int}}, AnovaEffect}()
    for i ∈ axes(cellmeans, si + 1)
        _, ifactorvarsdict = @views anovafactors(selectdim(cellmeans, si + 1, i), nreplicates, factornames[1:si])

        if isempty(factorvarsdict)
            for k ∈ keys(ifactorvarsdict)
                1 ∈ k || continue

                newk = ((1:namongfactors)..., (k .+ namongfactors)...)
                ifactor = ifactorvarsdict[k]
                factorvarsdict[newk] = AnovaFactor(makefactorname(factornames[1:si][invertfactorindex.(collect(k[2:end]), si)], factornames[si], factornames[(si + 1):end]), ifactor)
            end
        else
            for k ∈ keys(ifactorvarsdict)
                1 ∈ k || continue

                newk = ((1:namongfactors)..., (k .+ namongfactors)...)
                factorvarsdict[newk] += ifactorvarsdict[k]
            end
        end
    end
    return factorvarsdict
end

anovavalue(name, variance, df) = AnovaValue(name, variance * df, df)

meanfirstdim(observations::AbstractArray{<:Number}) = dropdims(mean(observations, dims = 1), dims = 1)
meanfirstdim(observations::AbstractArray{<:Array}) = dropdims(mean(observations, dims = 1), dims = 1)

upcat(x::AbstractArray) = x
upcat(x::AbstractArray{<:AbstractVector}) = reshape(vcat(x...), (size(x[1])..., size(x)...))

invertfactorindex(index, nfactors) = nfactors .- index .+ 1

makefactorname(factorname::AbstractString) = factorname

function makefactorname(factornames::AbstractVector{<:AbstractString})
    length(factornames) > 0 || return ""

    factorname = factornames[1]

    for i ∈ 2:length(factornames)
        factorname *= " × " * factornames[i]
    end

    return factorname
end

function makefactorname(withinfactornames::AbstractVector{<:AbstractString}, subjectfactorname, amongfactornames::AbstractVector{<:AbstractString})
    subjectname = subjectfactorname

    if !isempty(amongfactornames)
        subjectname *= "[" * makefactorname(amongfactornames) * "]"
    end

    makefactorname([withinfactornames; subjectname])
end

function threeway_random_error(interaction_ab, interaction_bc, interaction_abc)
    reducedmeansquare(factor::AnovaFactor) = factor.ms ^ 2 / factor.df
    ms = interaction_ab.ms + interaction_bc.ms - interaction_abc.ms
    df = ms ^ 2 / (reducedmeansquare(interaction_ab) + reducedmeansquare(interaction_bc) + reducedmeansquare(interaction_abc))
    AnovaFactor("", ms * df, df, ms)
end

function ftest(x, y)
    f = x.ms / y.ms
    fdist = FDist(x.df, y.df)
    p = ccdf(fdist, f)
    AnovaResult(x, f, p)
end

#=
Sources:
Olejnik and Algina, Generalized Eta and Omega Squared Statistics: Measures of Effect Size for Some Common Research Designs
Howell, Statistical methods of psychology (8th edition)
=#
function effectsizeω²(factorvars, factorerrorvars, nestedvars, nestederrorvars, N, subjectvars, errorvar, factortypes)
    σ²factorvars = σ²component(factorvars, factorerrorvars)
    σ²nestedvars = σ²component(nestedvars, nestederrorvars)

    containsrandomfactor = random .∈ combinations(factortypes)
    length(containsrandomfactor) == length(factorvars) || pop!(containsrandomfactor)
    σ²random = sum(σ²factorvars[containsrandomfactor])

    subjecterrorvars = (errorvar, subjectvars...)
    σ²nested = !isempty(nestedvars) ? sum(σ²nestedvars) : 0
    σ²error = N * sum(s.ss for s ∈ subjecterrorvars) / sum(s.df for s ∈ subjecterrorvars) # pool error from subject factors

    σ²id = σ²random + σ²nested + σ²error # sum of all random-effects variance
    ω² = σ²factorvars ./ (.!containsrandomfactor .* σ²factorvars .+ σ²id) # factor only needs to appear in denominator if a fixed-effect
    ω²nested = σ²nestedvars ./ σ²id

    return [ω²; ω²nested]
end

σ²component(factorvars, factorerrorvars) = [factorvars[i].ss - factorvars[i].df * factorerrorvars[i].ms for i ∈ eachindex(factorvars)]


text:
"""
    Cube Geoemtry

The cube goes from [-kL/2, kL/2] (with homogeneous distribution)
"""
function Atom(geometry::Cube, N::Int64, kL::Union{Real, Integer}; createFunction=ftn_AtomsOnCube::Function)
    @debug "start: Shape - Cube"

    dimensions = 3
    ρ = N / kL^3
    rₘᵢₙ = get_rₘᵢₙ(ρ)
    if rₘᵢₙ ≥ kL / 10
        rₘᵢₙ = kL / 100
    end

    r = get_atoms(dimensions, N, rₘᵢₙ; createFunction, kL)

    @debug "end  : Shape - Cube"
    return Atom(Cube(), r, N, Float64(kL))
end

function Atom(geometry::Cube, r::Matrix, kL::Union{Real, Integer})
    N = size(r, 2) # remember to use each collum as a atom position
    return Atom(Cube(), r, N, Float64(kL))
end

function ftn_AtomsOnCube(;kwargs...)
    kL = kwargs[:kL]

    x = -kL * rand() + (kL / 2)
    y = -kL * rand() + (kL / 2)
    z = -kL * rand() + (kL / 2)
    return [x, y, z]
end


text:
module CommunalHelperStationBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper

const themes = ["Normal", "Moon"]
const behaviors = ["Pulling", "Pushing"]

@mapdef Entity "CommunalHelper/StationBlock" StationBlock(
    x::Integer,
    y::Integer,
    width::Integer=16,
    height::Integer=16,
    theme::String="Normal",
    behavior::String="Pulling",
    customBlockPath::String="",
    customArrowPath::String="",
    customTrackPath::String="",
    speedFactor::Number=1.0,
    allowWavedash::Bool=false,
    wavedashButtonColor::String="5BF75B",
    wavedashButtonPressedColor::String="F25EFF",
    dashCornerCorrection::Bool=true,
)

const placements = Ahorn.PlacementDict(
    "Station Block ($theme, $behavior) (Communal Helper)" => Ahorn.EntityPlacement(
        StationBlock,
        "rectangle",
        Dict{String,Any}(
            "theme" => theme,
            "behavior" => behavior,
        ),
    ) for theme in themes, behavior in behaviors
)

Ahorn.editingOptions(entity::StationBlock) = Dict{String,Any}(
    "theme" => themes,
    "behavior" => behaviors,
)

Ahorn.minimumSize(entity::StationBlock) = 16, 16
Ahorn.resizable(entity::StationBlock) = true, true

function Ahorn.selection(entity::StationBlock)
    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

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::StationBlock) = renderStationBlock(ctx, entity)

function getSprites(blockSize::Integer, entity::StationBlock)
    theme = lowercase(get(entity, "theme", "normal"))
    behavior = lowercase(get(entity, "behavior", "pulling"))

    path = "objects/CommunalHelper/stationBlock/"
    if theme == "normal"
        arrowDirectory, block = (behavior == "pulling") ? ("arrow/", "/block") : ("altArrow/", "/alt_block")
    else
        arrowDirectory, block = (behavior == "pulling") ? ("moonArrow/", "/moon_block") : ("altMoonArrow/", "/alt_moon_block")
    end

    arrow = (blockSize <= 16) ? "small00" : ((blockSize <= 24) ? "med00" : "big00")
    button = (Bool(get(entity.data, "allowWavedash", false)) ? "_button" : "")

    return (
        (path * arrowDirectory * arrow),
        (path * "blocks" * block * button),
        (path * "button"),
        (path * "button_outline"),
    )
end

function renderStationBlock(ctx::Ahorn.Cairo.CairoContext, entity::StationBlock)
    x, y = Ahorn.position(entity)

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

    arrow, block, button, buttonOutline = getSprites(Integer(min(width, height)), entity)
    arrowSprite = Ahorn.getSprite(arrow, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    # Button
    if Bool(get(entity.data, "allowWavedash", false))
        hexColor = String(get(entity.data, "wavedashButtonColor", "FFFFFF"))
        color = (length(hexColor) == 6) ? hexToRGBA(hexColor) : (1.0, 1.0, 1.0, 1.0)

        for i in 1:tilesWidth
            tx = (i == 1) ? 0 : ((i == tilesWidth) ? 16 : 8)

            Ahorn.drawImage(ctx, buttonOutline, x + (i - 1) * 8, y - 4, tx, 0, 8, 8, tint=color)
            Ahorn.drawImage(ctx, button, x + (i - 1) * 8, y - 4, tx, 0, 8, 8, tint=color)
        end
    end

    # Block
    for i in 1:tilesWidth, j in 1:tilesHeight
        tx = (i == 1) ? 0 : ((i == tilesWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tilesHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, block, x + (i - 1) * 8, y + (j - 1) * 8, tx, ty, 8, 8)
    end

    Ahorn.drawImage(ctx, arrowSprite, x + floor(Int, (width - arrowSprite.width) / 2), y + floor(Int, (height - arrowSprite.height) / 2))
end

end


text:
using ADI
using Documenter
using HCIToolbox
using Literate: markdown

# weave in examples using literate
examplesdir(args...) = joinpath(@__DIR__, "..", "examples", args...)
outdir = joinpath(@__DIR__, "src", "examples")
markdown(examplesdir("betapictoris.jl"), outdir)
markdown(examplesdir("geometries.jl"), outdir)

# now do Documenter commands

setup = quote
    using ADI
    using Random
    Random.seed!(8799)
end

DocMeta.setdocmeta!(ADI, :DocTestSetup, setup; recursive = true)

doctest(ADI)

makedocs(;
    modules = [ADI, HCIToolbox],
    authors = "Miles Lucas <mdlucas@hawaii.edu>",
    repo = "https://github.com/juliahci/ADI.jl/blob/{commit}{path}#L{line}",
    sitename = "ADI.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliahci.github.io/ADI.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Introduction to HCI" => "introduction.md",
        "Getting Started with ADI.jl" => "gettingstarted.md",
        "Algorithms" => [
            "algorithms/classic.md",
            "algorithms/pca.md",
            "algorithms/nmf.md",
            "algorithms/greeds.md",
            "algorithms/api.md"
        ],
        "Framewise" => "framewise.md",
        "SDI" => "sdi.md",
        "Metrics" => "metrics.md",
        "Examples" => [
            "examples/betapictoris.md",
            "examples/geometries.md"
        ],
        "Benchmarks" => "benchmarks.md",
        "Index" => "api.md"
    ],
    doctest = false
)

deploydocs(;
    repo = "github.com/JuliaHCI/ADI.jl",
    push_preview = true
)


text:
@testset "exponential" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    @testset "SqExponentialKernel" begin
        k = SqExponentialKernel()
        @test kappa(k, x) ≈ exp(-x / 2)
        @test k(v1, v2) ≈ exp(-norm(v1 - v2)^2 / 2)
        @test kappa(SqExponentialKernel(), x) == kappa(k, x)
        @test metric(SqExponentialKernel()) == SqEuclidean()
        @test RBFKernel == SqExponentialKernel
        @test GaussianKernel == SqExponentialKernel
        @test SEKernel == SqExponentialKernel
        @test repr(k) == "Squared Exponential Kernel (metric = Euclidean(0.0))"
        @test KernelFunctions.iskroncompatible(k) == true

        k2 = SqExponentialKernel(; metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k)
        test_ADs(SEKernel)
    end
    @testset "ExponentialKernel" begin
        k = ExponentialKernel()
        @test kappa(k, x) ≈ exp(-x)
        @test k(v1, v2) ≈ exp(-norm(v1 - v2))
        @test kappa(ExponentialKernel(), x) == kappa(k, x)
        @test metric(ExponentialKernel()) == Euclidean()
        @test repr(k) == "Exponential Kernel (metric = Euclidean(0.0))"
        @test LaplacianKernel == ExponentialKernel
        @test KernelFunctions.iskroncompatible(k) == true

        k2 = ExponentialKernel(; metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k)
        test_ADs(ExponentialKernel)
    end
    @testset "GammaExponentialKernel" begin
        γ = 1.0
        k = GammaExponentialKernel(; γ=γ)
        @test k(v1, v2) ≈ exp(-norm(v1 - v2)^γ)
        @test kappa(GammaExponentialKernel(), x) == kappa(k, x)
        @test GammaExponentialKernel(; gamma=γ).γ == [γ]
        @test metric(GammaExponentialKernel()) == Euclidean()
        @test metric(GammaExponentialKernel(; γ=2.0)) == Euclidean()
        @test repr(k) == "Gamma Exponential Kernel (γ = $(γ), metric = Euclidean(0.0))"
        @test KernelFunctions.iskroncompatible(k) == true

        k2 = GammaExponentialKernel(; γ=γ, metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        test_ADs(γ -> GammaExponentialKernel(; gamma=only(γ)), [1 + 0.5 * rand()])
        test_params(k, ([γ],))
        TestUtils.test_interface(GammaExponentialKernel(; γ=1.36))

        #Coherence :
        @test isapprox(
            GammaExponentialKernel(; γ=2.0)(sqrt(0.5) * v1, sqrt(0.5) * v2),
            SqExponentialKernel()(v1, v2),
        )
        @test GammaExponentialKernel()(v1, v2) ≈ ExponentialKernel()(v1, v2)
    end
end


text:
using Meshing
using GeometryTypes
using LinearAlgebra: dot, norm
using FileIO

# generate an SDF of a sphere
sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
    sqrt(sum(dot(v,v))) - 1 # sphere
end

m = GLNormalMesh(sdf_sphere, MarchingCubes())

save("sphere.ply",m)