text:
"""
    hopper
    	model inspired by Raibert hopper
            from "Dynamically Stable Legged Locomotion"
		q = (x, z, t, r)
			x - lateral position
			z - vertical position
			t - body orientation
			r - leg length
"""
struct Hopper{T}
    # dimensions
    nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters

    # parameters
    mass_body::T
    mass_foot::T
    inertia_body::T
    body_radius::T
    foot_radius::T
    leg_len_max::T
    leg_len_min::T
    friction_joint::T

    # environment
    friction_body_world::T
    friction_foot_world::T
	gravity::T
end

# kinematics
function kinematics_body(::Hopper, q)
	[q[1],
	 q[2]]
end

function kinematics_body_jacobian(::Hopper, q)
    [1.0 0.0 0.0 0.0;
	 0.0 1.0 0.0 0.0]
end

function kinematics_foot(::Hopper, q)
	[q[1] + q[4] * sin(q[3]),
	 q[2] - q[4] * cos(q[3])]
end

function kinematics_foot_jacobian(::Hopper, q)
    [1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
	 0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]))]
end

# Lagrangian
function lagrangian(model::Hopper, q, q̇)
    L = 0.0

    # body kinetic energy
    v_body = q̇[1:2]
    ω_body = q̇[3]
    L += 0.5 * model.mass_body * dot(v_body, v_body)
    L += 0.5 * model.inertia_body * ω_body^2.0

    # body potential energy
    z_body = q[2]
    L -= model.mass_body * model.gravity * z_body

    # foot kinetic energy
    v_foot = kinematics_foot_jacobian(model, q) * q̇
    L += 0.5 * model.mass_foot * dot(v_foot, v_foot)

    # foot potential energy
    z_foot =  kinematics_foot(model, q)[2]
    L -= model.mass_foot * model.gravity * z_foot

    return L
end

# signed distance
function signed_distance(model::Hopper, q)
    [
     q[2] - model.body_radius,                    # body-world distance
     kinematics_foot(model, q)[2] - model.foot_radius, # foot-world distance
     q[4] - model.leg_len_min,                    # leg length min.
     model.leg_len_max - q[4]                     # leg length max.
    ]
end

function signed_distance(model::Hopper, q, pt1, pt2)
    # error between robot kinematics and point on nearest surface
    e1 = kinematics_body(model, q) - pt1
    e2 = kinematics_foot(model, q) - pt2
    @show pt1
    @show pt2
    [
     sqrt(dot(e1, e1)) - model.body_radius,  # body-world distance
     sqrt(dot(e2, e2)) - model.foot_radius,  # foot-world distance
     q[4] - model.leg_len_min,             # leg length min.
     model.leg_len_max - q[4]              # leg length max.
    ]
end

# input matrix
function input_jacobian(model, q)
	[
     0.0 -sin(q[3]);
     0.0 cos(q[3]);
     1.0 0.0;
     0.0 1.0
    ]
end

# dimensions
nq = 4
nu = 2
nw = 0

# parameters
mass_body = 1.0
mass_foot = 0.1
inertia_body = 0.1
body_radius = 0.1
foot_radius = 0.05
leg_len_max = 1.0
leg_len_min = 0.25
friction_joint = 0.1

# environment
friction_body_world = 0.5
friction_foot_world = 0.5
gravity = 9.81

# model
model = Hopper(nq, nu, nw,
            mass_body, mass_foot, inertia_body,
            body_radius, foot_radius,
            leg_len_max, leg_len_min,
            friction_joint,
            friction_body_world, friction_foot_world,
            gravity)

# code-gen

nq = model.nq
nu = model.nu
nw = model.nw

# variables
@variables q[1:nq] q̇[1:nq]

# Lagrangian
L = lagrangian(model, q, q̇)

ddL = Symbolics.hessian(L, [q; q̇])#, simplify=true)
dLq = Symbolics.gradient(L, q)#, simplify=true)
ddLq̇q = ddL[nq .+ (1:nq), 1:nq]

M = ddL[nq .+ (1:nq), nq .+ (1:nq)]
C = ddLq̇q * q̇ - dLq

mass_matrix_sym = eval(Symbolics.build_function(M, q)[1])
dynamics_bias_sym = eval(Symbolics.build_function(C, q, q̇)[1])

# ϕ = signed_distance(model, q)
# dϕ = Symbolics.jacobian(ϕ, q)
# signed_distance_jacobian = eval(Symbolics.build_function(dϕ, q)[1])

function contact_jacobian(model, q)
    [1.0 0.0 0.0 0.0; # note that body torque due to friction is added in residual
     0.0 1.0 0.0 0.0;
     1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
     0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]));
     0.0 0.0 0.0 1.0;
     0.0 0.0 0.0 -1.0]
end

function lagrangian_derivatives(model, q, v)
	D1L = -1.0 * dynamics_bias_sym(q, v)
    D2L = mass_matrix_sym(q) * v
	return D1L, D2L
end

function dynamics(model, h, q0, q1, u1, w1, λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2
		+ input_jacobian(model, qm2) * u1
        - model.friction_joint * [0.0; 0.0; 0.0; vm2[4]] # joint friction
        + λ1)
end

@variables h[1:1] q0[1:nq] q1[1:nq] u1[1:nu] w1[1:nw] λ1[1:4] q2[1:nq]
d = dynamics(model, h, q0, q1, u1, w1, λ1, q2)

function residual(model, z, θ, η, μ)
    # dimensions
    nq = model.nq
    nu = model.nu
    nw = model.nw

    # unpack data
    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    w1 = θ[2nq + nu .+ (1:nw)]
    friction_body_world = θ[2nq + nu + nw .+ (1:1)]
    friction_foot_world = θ[2nq + nu + nw + 1 .+ (1:1)]
    h = θ[2nq + nu + nw + 2 .+ (1:1)]

    # unpack contact points and directions
    pt_contact_body = η[1:2]
    pt_contact_foot = η[3:4]
    dir_contact_body = η[5:6]
    dir_contact_foot = η[7:8]

    # unpack variables
    q2 = z[1:nq]
    γ1 = z[nq .+ (1:4)]
    sγ1 = z[nq + 4 .+ (1:4)]
    ψ1 = z[nq + 4 + 4 .+ (1:2)]
    b1 = z[nq + 4 + 4 + 2 .+ (1:2)]
    sψ1 = z[nq + 4 + 4 + 2 + 2 .+ (1:2)]
    sb1 = z[nq + 4 + 4 + 2 + 2 + 2 .+ (1:2)]

    rotation_body = rotation(dir_contact_body)
    rotation_foot = rotation(dir_contact_foot)
    # rotation_body = rotation([0.0; 1.0])
    # rotation_foot = rotation([0.0; 1.0])

    # tangential velocity
    v1 = (q2 - q1) / h[1]
    vT_body = (rotation_body * v1[1:2])[1] + model.body_radius * v1[3]
    vT_foot = (rotation_foot * kinematics_foot_jacobian(model, q2) * v1)[1]

    # contact forces
    J = contact_jacobian(model, q2)
    λ1 = transpose(J) * [transpose(rotation_body) * [b1[1]; γ1[1]];
                         transpose(rotation_body) * [b1[2]; γ1[2]];
                         γ1[3:4]]
    λ1[3] += model.body_radius * b1[1] # friction on body creates a moment

	[
     dynamics(model, h, q0, q1, u1, w1, λ1, q2);
	 signed_distance(model, q2, pt_contact_body, pt_contact_foot) - sγ1;
    #  signed_distance(model, q2) - sγ1;
     vT_body - sb1[1];
     vT_foot - sb1[2];
     friction_body_world[1] * γ1[1] - ψ1[1];
     friction_foot_world[1] * γ1[2] - ψ1[2];
     γ1 .* sγ1 .- μ[1];
	 cone_product([ψ1[1]; b1[1]], [sψ1[1]; sb1[1]]) - [μ[1]; 0.0];
     cone_product([ψ1[2]; b1[2]], [sψ1[2]; sb1[2]]) - [μ[1]; 0.0];
    ]
end

nz = nq + 4 + 4 + 2 + 2 + 2 + 2
ny = 4 + 4 + 2 + 2 + 2 + 2
nθ = 2nq + nu + nw + 2 + 1
nη = 2 * (2 + 2)

@variables z[1:nz] θ[1:nθ] η[1:nη] μ[1:1]
r = residual(model, z, θ, η, μ)
rz = Symbolics.jacobian(r, z)
rη = Symbolics.jacobian(r, η)

rθ = Symbolics.jacobian(r, θ)
rz_sp = Symbolics.sparsejacobian(r, z)
rη_sp = Symbolics.sparsejacobian(r, η)
rθ_sp = Symbolics.sparsejacobian(r, θ)

r_hopper! = eval(Symbolics.build_function(r, z, θ, η, μ)[2])
rz_hopper! = eval(Symbolics.build_function(rz, z, θ, η)[2])
rη_hopper! = eval(Symbolics.build_function(rη, z, θ, η)[2])
rθ_hopper! = eval(Symbolics.build_function(rθ, z, θ, η)[2])

# r_ = eval(Symbolics.build_function(r, z, θ, η, μ)[1])
# rz_ = eval(Symbolics.build_function(rz, z, θ, η)[1])
# rη_ = eval(Symbolics.build_function(rη, z, θ, η)[1])
# rθ_ = eval(Symbolics.build_function(rθ, z, θ, η)[1])
# rz_sp! = eval(Symbolics.build_function(rz_sp.nzval, z, θ, η)[2])
# rη_sp! = eval(Symbolics.build_function(rη_sp.nzval, z, θ, η)[2])
# rθ_sp! = eval(Symbolics.build_function(rθ_sp.nzval, z, θ, η)[2])

# r0 = zeros(length(r))
# rz0 = zeros(size(rz))
# rθ0 = zeros(size(rθ))
# rη0 = zeros(size(rη))

# rz0_sp = zeros(size(rz_sp.nzval))
# rθ0_sp = zeros(size(rθ_sp.nzval))
# rη0_sp = zeros(size(rη_sp.nzval))

# z0 = rand(nz)
# θ0 = rand(nθ)
# η0 = rand(nη)

# μ0 = [1.0]

# using BenchmarkTools
# @benchmark r!($r0, $z0, $θ0, $η0, $μ0)
# @benchmark rz!($rz0, $z0, $θ0, $η0)
# @benchmark rη!($rη0, $z0, $θ0, $η0)
# @benchmark rθ!($rθ0, $z0, $θ0, $η0)

# @benchmark $r0 .= r_($z0, $θ0, $η0, $μ0)
# @benchmark $rz0 .= rz_($z0, $θ0, $η0)
# @benchmark $rη0 .= rη_($z0, $θ0, $η0)
# @benchmark $rθ0 .= rθ_($z0, $θ0, $η0)

# @benchmark rz_sp!($rz0_sp, $z0, $θ0, $η0)
# @benchmark rη_sp!($rη0_sp, $z0, $θ0, $η0)
# @benchmark rθ_sp!($rθ0_sp, $z0, $θ0, $η0)

##
# q0 = [0.0, 1.0, 0.0 * π, 0.5]
# z0[1:nq] = copy(q0)
# z0[nq .+ (1:4)] .= 1.0
# z0[nq + 4 .+ (1:4)] .= 1.0
# z0[nq + 4 + 4 .+ (1:2)] .= 1.0
# z0[nq + 4 + 4 + 2 .+ (1:2)] .= 0.1
# z0[nq + 4 + 4 + 2 + 2 .+ (1:2)] .= 1.0
# z0[nq + 4 + 4 + 2 + 2 + 2 .+ (1:2)] .= 0.1

# k_body = kinematics_body(model, q0)
# k_foot = kinematics_foot(model, q0)

# p_body = contact_point(surf, k_body, r=2.0 * model.body_radius)
# p_foot = contact_point(surf, k_foot, r=2.0 * model.body_radius)

# η0 = [p_body[1]; p_foot[1]; p_body[3]; p_foot[3]]
# rz0 = rz_(z0, θ0, η0)

# r1 = rη_(z0, θ0, η0)[:, 1:2] * p_body[6] * kinematics_body_jacobian(model, q0)
# r2 = rη_(z0, θ0, η0)[:, 2 .+ (1:2)] * p_foot[6] * kinematics_foot_jacobian(model, q0)
# r3 = rη_(z0, θ0, η0)[:, 4 .+ (1:2)] * p_body[7] * p_body[6] * kinematics_body_jacobian(model, q0)
# r4 = rη_(z0, θ0, η0)[:, 6 .+ (1:2)] * p_foot[7] * p_foot[6] * kinematics_foot_jacobian(model, q0)

# rη_update = r1 + r2 + r3 + r4

# rz0[:, 1:nq] += rη_update

function res!(r, z, θ, μ)
    q2 = z[1:nq]
    k_body = kinematics_body(model, q2)
    k_foot = kinematics_foot(model, q2)

    p_body = contact_point(surf, k_body, r=2.0 * model.body_radius)
    p_foot = contact_point(surf, k_foot, r=2.0 * model.foot_radius)

    η = [p_body[1]; p_foot[1]; p_body[3]; p_foot[3]]

    r_hopper!(r, z, θ, η, μ)
end

function res_z!(rz, z, θ)
    q2 = z[1:nq]
    k_body = kinematics_body(model, q2)
    k_foot = kinematics_foot(model, q2)

    p_body = contact_point(surf, k_body, r=2.0 * model.body_radius)
    p_foot = contact_point(surf, k_foot, r=2.0 * model.foot_radius)

    # @show p_body[1]
    # @show p_foot[1]

    η = [p_body[1]; p_foot[1]; p_body[3]; p_foot[3]]

    rz_hopper!(rz, z, θ, η)
    rη = rη_(z, θ, η)

    r1 = rη[:, 1:2] * p_body[6] * kinematics_body_jacobian(model, q2)
    r2 = rη[:, 2 .+ (1:2)] * p_foot[6] * kinematics_foot_jacobian(model, q2)
    r3 = rη[:, 4 .+ (1:2)] * p_body[7] * p_body[6] * kinematics_body_jacobian(model, q2)
    r4 = rη[:, 6 .+ (1:2)] * p_foot[7] * p_foot[6] * kinematics_foot_jacobian(model, q2)

    rη_update = r1 + r2 + r3 + r4

    rz[:, 1:nq] += rη_update
end

# res(z0, θ0, μ0)
# res_z(z0, θ0, μ0)

# res!(r0, z0, θ0, μ0)
# res_z!(rz0, z0, θ0)

# r0


text:
#####
##### `sum`
#####

function frule((_, ẋ), ::typeof(sum), x)
    return sum(x), sum(ẋ)
end

function rrule(::typeof(sum), x::AbstractArray{<:Real}; dims=:)
    y = sum(sum, x; dims=dims)
    function sum_pullback(ȳ)
        # broadcasting the two works out the size no-matter `dims`
        x̄ = @thunk broadcast(x, ȳ) do xi, ȳi
            ȳi
        end
        return (NO_FIELDS, x̄)
    end
    return y, sum_pullback
end

function rrule(::typeof(sum), ::typeof(abs2), x::AbstractArray{<:Real}; dims=:)
    y = sum(abs2, x; dims=dims)
    function sum_abs2_pullback(ȳ)
        return (NO_FIELDS, DoesNotExist(), @thunk(2ȳ .* x))
    end
    return y, sum_abs2_pullback
end

ext: jl


text:
module TestSignalsAPI

  using Test, Dates
  using NumerAPI.Signals

  public_id = ENV["NUMERAI_PUBLIC_ID"]
  secret_key = ENV["NUMERAI_SECRET_KEY"]
  username = ENV["SIGNALS_USERNAME"]
  model_id = ENV["SIGNALS_MODEL_ID"]

  signals_api = SignalsAPI(public_id, secret_key)

  @testset "test get_leaderboard" begin
    data = get_leaderboard(signals_api)

    @test typeof(data) <: Vector
    @test length(data) == 50
  end

  @testset "test upload_predictions" begin

  end

  @testset "test daily_submissions_performances" begin
    data = submission_status(signals_api, model_id = model_id)

    @test typeof(data) <: Vector
  end

  @testset "test public_user_profile" begin
    data = public_user_profile(signals_api, username)

    @test typeof(data) <: Dict
    @test haskey(data, "username")
    @test data["username"] == username
  end

  @testset "test daily_user_performances" begin
    data = daily_user_performances(signals_api, username)

    @test typeof(data) <: Vector{Dict}
  end

  @testset "test ticker_universe" begin
    tickers = ticker_universe(signals_api)

    @test typeof(tickers) == Vector{String}
  end

  @testset "test daily_submissions_performances" begin
    data = daily_submissions_performances(signals_api, username)

    @test typeof(data) <: Vector{Dict}
  end

  @testset "test download_validation_data" begin
    data_path = download_validation_data(signals_api)

    @test data_path === "./numerai_signals_historical.csv"
    @test isfile(data_path)

    rm(data_path, force=true)
  end

  @testset "test stake_get" begin
    data = stake_get(signals_api, username)

    @test typeof(data) <: Real
  end

end # module


text:
using StaticArrays

export delta_tetrahedron
export delta_parallelepiped
export delta_parallelepiped_vertex
export delta_tetrahedron_weights
export delta_tetrahedron_polynomial

# Routines for tetrahedron integration of delta function

# const KLIST = [SVector{3,Int}(-1, -1, -1),
#                SVector{3,Int}(+1, -1, -1),
#                SVector{3,Int}(-1, +1, -1),
#                SVector{3,Int}(+1, +1, -1),
#                SVector{3,Int}(-1, -1, +1),
#                SVector{3,Int}(+1, -1, +1),
#                SVector{3,Int}(-1, +1, +1),
#                SVector{3,Int}(+1, +1, +1)]

@inline function delta_tetrahedron(etarget, e1234)
    """
    Calcultate val = int_{tetrahedron} d^3k delta(etarget - e(k)) / vol
    where e(k) is a linear function, e1, ..., e4 are the values of e(k) at
    the four vertices of the tetrahedron, and vol is the volume of the tetrahedron.

    Ref: Eqs. (7-10) of P. B. Allen, Phys. Stat. Sol. (b) 120, 629 (1983).
    Note in case B of Eq. (8), we use an equivalent formula which is more stable
    when (e2 - e1) is small.

    # Input
    - etarget: energy where the delta function becomes nonzero
    - e1234: Array containing energy at the four vertices of the tetrahedron
    """
    e1, e2, e3, e4 = sort(e1234)
    if e1 < etarget <= e2
        c1 = 3 / ((e2 - e1) * (e3 - e1) * (e4 - e1))
        return (etarget - e1)^2 * c1
    elseif e2 < etarget <= e3
        e31 = e3 - e1
        e41 = e4 - e1
        c2 = 3 / (e41 * e31)
        c3 = 3 * (e4 + e3 - e2 - e1) / (e41 * e31 * (e4 - e2) * (e3 -e2))
        return (2*etarget - e1 - e2) * c2 - (etarget - e2)^2 * c3
    elseif e3 < etarget <= e4
        c4 = 3 / ((e4 - e1) * (e4 - e2) * (e4 - e3))
        return (e4 - etarget)^2 * c4
    else
        # etarget <= e1 or e4 < etarget.
        return zero(etarget)
    end
end

@inline function delta_parallelepiped(etarget::T, e0, v0, L) where {T}
    """
    Calcultate val = int_{parallelepiped} d^3k delta(`etarget` - e(k)) / volume
    where e(k) = `e0` + k ⋅ `v0`. The parallelepiped region is [-L/2, L/2]^3.
    Divide the parallelepiped into six tetrahedra and use tetrahedron integration.
    """
    if norm(v0) < 1E-10
        # Velocity is too small. Single-point tetrahedron method does not work.
        return zero(T)
    end
    L_div_2 = L ./ 2
    v0_L = v0 .* L_div_2
    e1 = e0 - v0_L[1] - v0_L[2] - v0_L[3]
    e2 = e0 + v0_L[1] - v0_L[2] - v0_L[3]
    e3 = e0 - v0_L[1] + v0_L[2] - v0_L[3]
    e4 = e0 + v0_L[1] + v0_L[2] - v0_L[3]
    e5 = e0 - v0_L[1] - v0_L[2] + v0_L[3]
    e6 = e0 + v0_L[1] - v0_L[2] + v0_L[3]
    e7 = e0 - v0_L[1] + v0_L[2] + v0_L[3]
    e8 = e0 + v0_L[1] + v0_L[2] + v0_L[3]
    delta_parallelepiped_vertex(etarget, e1, e2, e3, e4, e5, e6, e7, e8)
end

@inline function delta_parallelepiped_vertex(etarget::FT, e1, e2, e3, e4, e5, e6, e7, e8) where {FT}
    """
    Calcultate val = int_{parallelepiped} d^3k delta(`etarget` - e(k)) / volume.
    `e1`, ..., `e8` are the energies at the eight vertices of the parallelepiped.
    Divide the parallelepiped into six tetrahedra and use tetrahedron integration.
    """
    val = zero(FT)
    # val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e2, e4, e5))
    # val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e3, e4, e5))
    # val += delta_tetrahedron(etarget, SVector{4,FT}(e2, e4, e6, e5))
    # val += delta_tetrahedron(etarget, SVector{4,FT}(e3, e4, e7, e5))
    # val += delta_tetrahedron(etarget, SVector{4,FT}(e4, e7, e8, e5))
    # val += delta_tetrahedron(etarget, SVector{4,FT}(e4, e6, e8, e5))
    val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e8, e6, e2))
    val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e8, e2, e4))
    val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e8, e4, e3))
    val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e8, e3, e7))
    val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e8, e7, e5))
    val += delta_tetrahedron(etarget, SVector{4,FT}(e1, e8, e5, e6))
    return val / 6
end

function delta_tetrahedron_weights(etarget, e1234)
    """
    Calcultate weights for evaluating I = int_{tetrahedron} d^3k f(k) delta(etarget - e(k)) / vol
    where e(k) is a linear function, e1, ..., e4 are the values of e(k) at
    the four vertices of the tetrahedron, and vol is the volume of the tetrahedron.
    With the output `w1234`, `I = sum_{i=1}^{4} f(k_i) w_i` holds.

    Ref: Appendix B of Blöchl et al, PRB 16 232 (1994). Since the reference is for step functions,
    we take the derivative with respect to e_F in all equations in the reference.

    # Input
    - etarget: energy where the delta function becomes nonzero
    - e1234: Array containing energy at the four vertices of the tetrahedron
    # Output
    - w1234: Array containing weights at the four vertices
    """
    inds = sortperm(e1234)
    e1, e2, e3, e4 = e1234[inds]
    allzero = false
    if e1 < etarget <= e2
        c11 = (etarget - e1)^2 / ((e2 - e1) * (e3 - e1) * (e4 - e1))
        c12 = (etarget - e1) * c11
        w2 = c12 / (e2 - e1)
        w3 = c12 / (e3 - e1)
        w4 = c12 / (e4 - e1)
        w1 = 3 * c11 - (w2 + w3 + w4)
    elseif e2 < etarget <= e3
        # local variables
        # eei = etarget - ei
        ee1 = etarget - e1
        ee2 = etarget - e2
        ee3 = etarget - e3
        ee4 = etarget - e4
        # eij = ei - ej
        e41 = e4 - e1
        e31 = e3 - e1
        e42 = e4 - e2
        e32 = e3 - e2
        # Eq. (B11)
        c1  = ee1^2 / (4 * e41 * e31)
        dc1 = ee1 / (2 * e41 * e31)
        # Eq. (B12)
        c2  = - ee1 * ee2 * ee3 / (4 * e41 * e32 * e31 )
        dc2 = ( - ee2 * ee3 - ee1 * ee3 - ee1 * ee2 ) / (4 * e41 * e32 * e31)
        # Eq. (B13)
        c3  = - ee2^2 * ee4 / (4 * e42 * e32 * e41)
        dc3 = - (2 * ee2 * ee4 + ee2^2) / (4 * e42 * e32 * e41)

        # weights: derivatives of Eqs. (B7-B10)
        w1 = dc1 - ((dc1 + dc2) * ee3 + c1 + c2) / e31 - ((dc1 + dc2 + dc3) * ee4 + c1 + c2 + c3) / e41
        w2 = dc1 + dc2 + dc3 - ((dc2 + dc3) * ee3 + c2 + c3) / e32 - (dc3 * ee4 + c3) / e42
        w3 = ((dc1 + dc2 ) * ee1 + c1 + c2) / e31 + ((dc2 + dc3) * ee2 + c2 + c3) / e32
        w4 = ((dc1 + dc2 + dc3) * ee1 + c1 + c2 + c3) / e41 + (dc3 * ee2 + c3) / e42
    elseif e3 < etarget <= e4
        c31 = (e4 - etarget)^2 / ((e4 - e1) * (e4 - e2) * (e4 - e3))
        c32 = (e4 - etarget) * c31
        w1 = c32 / (e4 - e1)
        w2 = c32 / (e4 - e2)
        w3 = c32 / (e4 - e3)
        w4 = 3 * c31 - (w1 + w2 + w3)
    else
        w1, w2, w3, w4 = zero(etarget), zero(etarget), zero(etarget), zero(etarget)
        allzero = true
    end

    # Blöchl correction Eq. (22)
    esum = sum(e1234)
    if e1 < etarget <= e2
        dw = 6 * (etarget - e1) / ((e2 - e1) * (e3 - e1) * (e4 - e1))
    elseif e2 < etarget <= e3
        dw = 6 * (1 - (etarget - e2) * (e4 + e3 - e2 - e1) / ((e4 - e2) * (e3 -e2))) / (e41 * e31)
    elseif e3 < etarget <= e4
        dw = 6 * (etarget - e4) / ((e4 - e1) * (e4 - e2) * (e4 - e3))
    end
    if ! allzero
        w1 += (esum - 4 * e1) * dw / 40
        w2 += (esum - 4 * e2) * dw / 40
        w3 += (esum - 4 * e3) * dw / 40
        w4 += (esum - 4 * e4) * dw / 40
    end

    # w1234 are the weights in the sorted index. Need to invert the indices to original index.
    w1234 = SVector(w1, w2, w3, w4)
    w1234[sortperm(inds)]
end

# Return the piecewise polynomial function for the delta function integration result.
@inline function delta_tetrahedron_polynomial(e1234)
    T = eltype(e1234)
    e1, e2, e3, e4 = sort(e1234)
    # e1 < etarget <= e2
    c1 = 3 / ((e2 - e1) * (e3 - e1) * (e4 - e1))
    # e2 < etarget <= e3
    e31 = e3 - e1
    e41 = e4 - e1
    c2 = 3 / (e41 * e31)
    c3 = 3 * (e4 + e3 - e2 - e1) / (e41 * e31 * (e4 - e2) * (e3 -e2))
    # e3 < etarget <= e4
    c4 = 3 / ((e4 - e1) * (e4 - e2) * (e4 - e3))
    ((e1, e2, c1 .* (e1^2, -2*e1, 1)),
     (e2, e3, c2 .* (-e1-e2, 2, 0) .- c3 .* (e2^2, -2*e2, 1)),
     (e3, e4, c4 .* (e4^2, -2*e4, 1))
    )
end

function _evaluate_piecewise_polynomial(x, polys)
    for (e1, e2, coeffs) in polys
        e1 < x <= e2 && return evalpoly(x, coeffs)
    end
    return zero(x)
end

# TODO: Blöchl correction (Eq. (22))
# TODO: Add tests (e.g. sum of weights == total weight)
# TODO: Which diagonal to cut tetrahedron

max_stars_repo_name: jaemolihm/TetrahedronIntegration

max_issues_repo_name: jaemolihm/TetrahedronIntegration

max_forks_repo_name: jaemolihm/TetrahedronIntegration


text:
struct LayerNotFoundException <: Exception
    var::String
end

function Base.showerror(io::IO, e::LayerNotFoundException)
    println(io, typeof(e), ": ", e.var)
end


text:
module Core

export callinstancemethod, callstaticmethod, callconstructor

using JavaCall.JNI
using JavaCall.Signatures
using JavaCall.CodeGeneration

dispatch_rettype(::Type{jboolean}) = jboolean
dispatch_rettype(::Type{jbyte}) = jbyte
dispatch_rettype(::Type{jchar}) = jchar
dispatch_rettype(::Type{jshort}) = jshort
dispatch_rettype(::Type{jint}) = jint
dispatch_rettype(::Type{jlong}) = jlong
dispatch_rettype(::Type{jfloat}) = jfloat
dispatch_rettype(::Type{jdouble}) = jdouble
dispatch_rettype(::Type{jvoid}) = jvoid
# Default return type for all types is jobject when calling a method
dispatch_rettype(::T) where {T} = jobject

callinstancemethod(receiver::jobject, methodname::Symbol, rettype::Any, argtypes::Vector{Any}, args::Vararg{Any, N}) where N =
    callinstancemethod(Val(dispatch_rettype(rettype)), receiver, methodname, rettype, argtypes, args...)

for (type, method) in [
        (:jboolean, :call_boolean_method_a),
        (:jbyte, :call_byte_method_a),
        (:jchar, :call_char_method_a),
        (:jshort, :call_short_method_a),
        (:jint, :call_int_method_a),
        (:jlong, :call_long_method_a),
        (:jfloat, :call_float_method_a),
        (:jdouble, :call_double_method_a),
        (:jobject, :call_object_method_a),
        (:jvoid, :call_void_method_a)
    ]
    params = [
        :(::Val{$type}),
        :(receiver::Any),
        :(methodname::Symbol),
        :(rettype::Any),
        :(argtypes::Vector{Any}),
        :(args::Vararg{Any, N})
    ]
    body = quote
        @assert length(args) == length(argtypes)

        receiver_class = JNI.get_object_class(receiver)
        methodsignature = MethodSignature(rettype, argtypes)
        method = JNI.get_method_id(receiver_class, string(methodname), signature(methodsignature))
        JNI.$method(receiver, method, jvalue[args...])
    end
    eval(generatemethod(:callinstancemethod, params, body, :N))
end

callstaticmethod(class::jclass, methodname::Symbol, rettype::Any, argtypes::Vector{Any}, args::Vararg{Any, N}) where N =
    callstaticmethod(Val(dispatch_rettype(rettype)), class, methodname, rettype, argtypes, args...)

for (type, method) in [
        (:jboolean, :call_static_boolean_method_a),
        (:jbyte, :call_static_byte_method_a),
        (:jchar, :call_static_char_method_a),
        (:jshort, :call_static_short_method_a),
        (:jint, :call_static_int_method_a),
        (:jlong, :call_static_long_method_a),
        (:jfloat, :call_static_float_method_a),
        (:jdouble, :call_static_double_method_a),
        (:jobject, :call_static_object_method_a),
        (:jvoid, :call_static_void_method_a)
    ]
    params = [
        :(::Val{$type}),
        :(class::jclass),
        :(methodname::Symbol),
        :(rettype::Any),
        :(argtypes::Vector{Any}),
        :(args::Vararg{Any, N})
    ]
    body = quote
        @assert length(args) == length(argtypes)

        methodsignature = MethodSignature(rettype, argtypes)
        method = JNI.get_static_method_id(class, string(methodname), signature(methodsignature))
        JNI.$method(class, method, jvalue[args...])
    end
    eval(generatemethod(:callstaticmethod, params, body, :N))
end

callinstancemethod(receiver::jobject, methodname::Symbol, rettype::Any, signature::String, args::Vararg{Any, N}) where N =
    callinstancemethod(Val(dispatch_rettype(rettype)), receiver, methodname, signature, args...)

for (type, method) in [
        (:jboolean, :call_boolean_method_a),
        (:jbyte, :call_byte_method_a),
        (:jchar, :call_char_method_a),
        (:jshort, :call_short_method_a),
        (:jint, :call_int_method_a),
        (:jlong, :call_long_method_a),
        (:jfloat, :call_float_method_a),
        (:jdouble, :call_double_method_a),
        (:jobject, :call_object_method_a),
        (:jvoid, :call_void_method_a)
    ]
    params = [
        :(::Val{$type}),
        :(receiver::Any),
        :(methodname::Symbol),
        :(signature::String),
        :(args::Vararg{Any, N})
    ]
    body = quote
        receiver_class = JNI.get_object_class(receiver)
        method = JNI.get_method_id(receiver_class, string(methodname), signature)
        JNI.$method(receiver, method, jvalue[args...])
    end
    eval(generatemethod(:callinstancemethod, params, body, :N))
end

callstaticmethod(class::jclass, methodname::Symbol, rettype::Any, signature::String, args::Vararg{Any, N}) where N =
    callstaticmethod(Val(dispatch_rettype(rettype)), class, methodname, signature, args...)

for (type, method) in [
    (:jboolean, :call_static_boolean_method_a),
    (:jbyte, :call_static_byte_method_a),
    (:jchar, :call_static_char_method_a),
    (:jshort, :call_static_short_method_a),
    (:jint, :call_static_int_method_a),
    (:jlong, :call_static_long_method_a),
    (:jfloat, :call_static_float_method_a),
    (:jdouble, :call_static_double_method_a),
    (:jobject, :call_static_object_method_a),
    (:jvoid, :call_static_void_method_a)
]
    params = [
        :(::Val{$type}),
        :(class::jclass),
        :(methodname::Symbol),
        :(signature::String),
        :(args::Vararg{Any, N})
    ]
    body = quote
        method = JNI.get_static_method_id(class, string(methodname), signature)
        JNI.$method(class, method, jvalue[args...])
    end
    eval(generatemethod(:callstaticmethod, params, body, :N))
end

function callconstructor(class::jclass, signature::String, args::Vararg{Any, N}) where N
    # Construct name is always <init> as specified by the JNI reference
    constructor = JNI.get_method_id(class, "<init>", signature)
    JNI.new_object_a(class, constructor, jvalue[args...])
end

end


text:
"""
Abstract idea which defines the external boundary of object.
"""
abstract type Shape{T<:AbstractFloat,Dim} end

Symmetry(s::Shape{T,Dim}) where {T,Dim} = WithoutSymmetry{Dim}()

"""
    origin(shape::Shape)::SVector

Origin of shape, typically the center
"""
origin(shape::Shape) = shape.origin

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

True if shapes are the same but in different positions (origins), standard mathematical definition.
"""
iscongruent(s1::Shape, s2::Shape) = false # false by default, overload in specific examples

# Define synonym for iscongruent ≅, and add documentation
≅(s1::Shape, s2::Shape) = iscongruent(s1, s2)
@doc (@doc iscongruent(::Shape, ::Shape)) (≅(::Shape, ::Shape))

"""
    congruent(s::Shape, x)::Shape

Create shape congruent to `s` but with origin at `x`
"""
function congruent end

"Generic helper function which tests if boundary coordinate is between 0 and 1"
function check_boundary_coord_range(t)
    if t < 0 || t > 1
        throw(DomainError("Boundary coordinate must be between 0 and 1"))
    end
end

# Concrete shapes
# include("rectangle.jl")
include("box.jl")
include("sphere.jl")
include("halfspace.jl")
include("plate.jl")
include("time_of_flight.jl")
include("time_of_flight_from_point.jl")
include("empty_shape.jl")

"""
    points_in_shape(Shape; res = 20, xres = res, yres = res,
             exclude_region = EmptyShape(region), kws...)

returns `(x_vec, region_inds)` where `x_vec` is a vector of points that cover a box which bounds `Shape`, and `region_inds` is an array of linear indices such that `x_vec[region_inds]` are points contained `Shape`. For 3D we use `zres` instead of `yres`.

"""
function points_in_shape(region::Shape{T,2};
        res::Number = 20, xres::Number = res, yres::Number = res,
        exclude_region::Shape = EmptyShape(region),
        kws...) where T

    rect = bounding_box(region)

    #Size of the step in x and y direction
    x_vec_step = rect.dimensions ./ [xres, yres]
    bl = bottomleft(rect)
    x_vec = [SVector{2}(bl + x_vec_step .* [i,j]) for i=0:xres, j=0:yres][:]
    region_inds = findall(x -> !(x ∈ exclude_region) && x ∈ region, x_vec)

    return x_vec, region_inds
end

function points_in_shape(region::Shape{T,3};
        y::T = zero(T),
        res::Number = 20, xres::Number = res, zres::Number = res,
        exclude_region::Shape = EmptyShape(region)) where T

    box = bounding_box(region)

    #Size of the step in x and z direction
    x_vec_step = [box.dimensions[1] / xres, zero(T), box.dimensions[3] / zres]

    bl = corners(box)[1]
    x_vec = [SVector{3}(bl + x_vec_step .* [i,y,j]) for i = 0:xres, j = 0:zres][:]
    region_inds = findall(x -> !(x ∈ exclude_region) && x ∈ region, x_vec)

    return x_vec, region_inds
end

"Returns a set of points on the boundary of a 2D shape."
function boundary_points(shape2D::Shape{T,2}, num_points::Int = 4; dr = zero(T)) where T
    x, y = boundary_functions(shape2D)
    v(τ) = SVector(x(τ),y(τ)) + dr * (SVector(x(τ),y(τ)) - origin(shape2D))
    return [ v(τ) for τ in LinRange(zero(T),one(T),num_points+1)[1:end-1] ]
end

"Returns a set of points on the boundary of a 3D shape."
function boundary_points(shape3D::Shape{T,3}, num_points::Int = 4; dr = zero(T)) where T
    x, y, z = boundary_functions(shape3D)
    v(τ,s) = SVector(x(τ,s),y(τ,s),z(τ,s)) + dr * (SVector(x(τ,s),y(τ,s),z(τ,s)) - origin(shape3D))
    mesh = LinRange(zero(T),one(T),num_points+1)[1:end-1]
    return [v(τ,s) for τ in mesh, s in mesh]
end

"Returns box which completely encloses the shapes"
bounding_box(shape1::Shape, shape2::Shape) = bounding_box([shape1, shape2])

# Create a box which bounds an array of shapes
function bounding_box(shapes::Vector{S}) where S<:Shape
    boxes = bounding_box.(shapes)
    corners_mat = hcat(vcat((corners.(boxes))...)...)

    maxdims = maximum(corners_mat, dims=2)
    mindims = minimum(corners_mat, dims=2)

    c = (maxdims + mindims) ./ 2
    dimensions = maxdims - mindims

    return Box(c[:], dimensions[:])
end

# Docstrings
"""
    name(shape::Shape)::String

Name of a shape
"""
function name end

"""
    outer_radius(shape::Shape{T})::T

The radius of a circle which completely contains the shape
"""
function outer_radius end

"""
    volume(shape::Shape{T})::T

Volume of a shape
"""
function volume end

"""
    volume(shape::Shape)::NTuple{Function,Dim)

Returns Tuple of Dim Functions which define outer boundary of shape when given boundary coordinate t∈[0,1]
"""
function boundary_functions end

"""
    issubset(shape1, shape2)::Bool

Returns true if shape1 is entirely contained within shape2, false otherwise (also works with particles).
"""
function issubset(s1::Shape,s2::Shape) throw(MethodError(issubset, (s1,s2))) end

"""
    in(vector, shape)::Bool

Returns true if vector is in interior of shape, false otherwise.
"""
function in(::AbstractVector,::Shape) end


text:
# ---------------------------------------------------------------------------
# automation.jl: Various functions for synthetic data generation and analysis
# ---------------------------------------------------------------------------

"""
Randomly sample parameters for tumour evolution simulation.

Input:
    - None

Returns:
    - integer: Mutation rate (per genome per doubling)
    - float: Probability of subclone/driver event (per division probability)
    - integer: Number of clonal mutations at t = 0
    - integer: Number of potential subclone/driver events during tumour growth
    - float: Mean of exponential selection coefficient distribution
    - integer: Mean sequencing depth
    - float: Sequencing overdispersion parameter for beta-binomial distribution
"""
function sampleParameters()
    u = rand(DiscreteUniform(1, 500), 1)[1]
    rho = 10^-rand(Uniform(2, 3), 1)[1]
    phi_b = rand(Uniform(0.001, 0.1), 1)[1] / u # Normalize by u so time of driver occurences are similar across different mutation rates
    nclonal = rand(DiscreteUniform(0, 5000), 1)[1]
    ndrivers = rand(DiscreteUniform(1, 3), 1)[1]
    lambda = rand(Uniform(1, 3), 1)[1]
    depth = rand(DiscreteUniform(50, 250), 1)[1]
    return u, phi_b, nclonal, ndrivers, lambda, depth, rho
end

"""
Generates synthetic sequenced tumour data for a population subject to positive selection and a population evolving neutrally.

Input:
    - b (float): Birth rate
    - d (float): Death rate
    - Nfinal (integer): Final tumour population size (default = 1000)
    - noise (string): Sequencing noise model (default = "betabinomial")
    - lower_cutoff (float): Lower frequency (range 0 - 0.5) bound for detecting subclones (default = 0.09)
    - upper_cutoff (float): Upper frequency (range 0 - 0.5) bound for detecting subclones (default = 0.41)
    - nsubclone_min (integer): Minimum number of subclones within (lower_cutoff, upper_cutoff) to accept for a positive selection simulation
    - nsubclone_max (integer): Maximum number of subclones within (lower_cutoff, upper_cutoff) to accept for a positive selection simulation
    - alt_reads (integer): Minimum number of alternate reads to call a mutation during virtual sequencing

Returns:
    - nested array: Positive selection simulation output
        - string: Simulation id
        - array: VAF, reads, depth for each mutation where each entry is an array with a length equal to number of detectable mutations
        - array: Subclone/driver fitness, timing, frequency information
        - array: Mutation count information
        - array: Simulation parameters
    - nested array: Neutral evolution simulation output
        - string: Simulation id
        - array: VAF, reads, depth for each mutation where each entry is an array with a length equal to number of detectable mutations
        - array: Subclone/driver fitness, timing, frequency information
        - array: Mutation count information
        - array: Simulation parameters
"""
function autoSimulation(b, d; Nfinal::Int64 = 1000, noise::String = "betabinomial", lower_cutoff::Float64 = 0.1, upper_cutoff::Float64 = 0.4, nsubclone_min::Int64 = 1, nsubclone_max::Int64 = 2, alt_reads::Int64 = 2, trim = true)

    # Run positive selection simulation
    u, phi_b, nclonal, ndrivers, lambda, depth, rho = sampleParameters()
    println("{A} Starting positive simulation")
    println("--------------------------------")
    println("")
    n_detectable, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters = 0, [], [], [], []
    counter = 0
    while true
        counter += 1
        if counter == 10
            u, phi_b, nclonal, ndrivers, lambda, depth = sampleParameters()
            counter = 0
            println("Parameters failed to generate subclone number of interest, resampling parameters...")
        end
        n_detectable, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters = simulation(b, d, u, Nfinal = Nfinal, phi_b = phi_b, phi_d = 0.0, nclonal = nclonal, depth = depth, noise = noise, rho = rho, lower_cutoff = lower_cutoff, upper_cutoff = upper_cutoff, ndrivers = ndrivers, lambda = lambda, alt_reads = alt_reads)
        if nsubclone_max >= n_detectable >= nsubclone_min
            println("")
            println("Correct number of subclones generated for synthetic tumour with positive selection...")
            println("")
            break
        end
        println("Non-specified number of subclones generated, re-simulating positive selection...")
    end

    # Generate a synthetic neutral simulation using matched mutational burden parameters from positive selection simulation
    println("{B} Generating parameter paired fully synthetic neutral data")
    println("------------------------------")
    println("")

    scale_pos, scale_clonal = rand(Uniform(0.8, 1.2)), rand(Uniform(0.8, 1.2)) # Add variability to clonal and subclonal mutation counts in neutral synthetic tumours
    if trim == true
        trim_tail = ifelse(rand(Uniform(0, 1), 1)[1] < 0.05, 1, 0) # Randomly remove tail from neutral distributions to inject more sparsity in neutral samples
    else
        trim_tail = 0
    end
    n_detectable, neutral_VAF, neutral_haplotypes, neutral_mutations, neutral_parameters = full_synthetic_neutral(depth = depth, Nfinal = Nfinal, nclonal = Int(round(nclonal * scale_clonal)), match_pos_muts = Int(round(positive_mutations[5] * scale_pos)), noise = noise, rho = rho, alt_reads = alt_reads, trim_tail = trim_tail)
    neutral_parameters[3] = u

    println("Saving synthetic data...")
    simid = randstring(10)
    positive_id, neutral_id = join(["P_", simid]), join(["N_", simid])

    return [positive_id, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters], [neutral_id, neutral_VAF, neutral_haplotypes, neutral_mutations, neutral_parameters]
end

"""
Generates a histogram of the mutation frequency information for a given tumour evolution simulation.

Input:
    - CanEvolve simulation output: A single simulation object e.g. p from p, n = autoSimulation(1, 0.1)
    - bins (integer): Number of bins in histogram

Returns:
    - Plots.Plot{Plots.GRBackend}: A histogram of the VAF information from a simulation. If sim_object is positive selection, blue lines indicate subclone frequency (cellular fraction / 2)
"""
function plotVAF(sim_object; bins::Int64 = 100)
    VAF = sim_object[2][1]
    scaler = fit(Histogram, VAF, nbins = bins)
    scaler = normalize(scaler, mode = :pdf)
    plot(scaler, fillcolor = "#008EA0", linecolor = "#008EA0", linealpha = 0.5, legend = false, xlims = (0,1), ytickfont = font(8, "Helvetica"), xtickfont = font(8, "Helvetica"), xtickfontsize=8,ytickfontsize=8)
    if sum(sim_object[3][1]) != 0
        vline!(sim_object[3][1], color = RGB(118/255,140/255,195/255), alpha = 1, linewidth = 3)
    end
    xlabel!("Variant allele frequency (VAF)")
    ylabel!("Density")
    density!(VAF, linecolor = RGBA(249/255,101/255,103/255,0), linewidth = 5, linealpha = 0.95)
end

"""
Computes frequency threshold under a binomial sequencing noise, where variance = np(1-p), model where the lower cutoff is variable to account for additional variability in empirical data

Input:
    - depth (integer/float): Mean sequencing depth
    - alt_reads (integer): Minimum number of reads covering mutation (default = 2)

Returns:
    - float: The lower frequency cutoff for VAFs that determines the limit of detection prior to adding sequencing noise.
    - float: The upper frequency cutoff that is fixed at 3 binomial standard deviations from 0.5. ## Note, this is generally not used in this framework
"""
function frequencyThresholds(depth; alt_reads::Int64 = 2)
    alt_reads = rand(Uniform(2, 4), 1)[1]
    tail_shift = rand(Uniform(1, 4), 1)[1]
    f_min = (alt_reads/depth) + ( ( tail_shift*sqrt(alt_reads*(1-(alt_reads/depth))) ) / depth)
    f_max = 0.5 - ( ( 3.0*sqrt((0.5*depth)*(1-0.5)) ) / depth)
    return f_min, f_max
end

"""
Generates a histogram representation of the VAF distribution.

Input:
    - array: VAF information for each mutation
    - depth (integer): Mean sequencing depth (default = 100)
    - alt_reads (integer): Minimum number of reads covering mutation (default = 2)
    - k (integer): Number of bins to segment the VAF distribution
    - range (array): An array indicating lower and upper bound to constrain histogram range (default = [0.02, 0.5])
    - cutoff (integer): A 0 or 1 indicating if the lower bound should be based on range[1] or on 2 binomial standard devations above alt_reads/depth

Retuns:
    - array: An array of length k where each entry indicates the number of mutations in bins of width (range[2] - range[1])/k
"""
function uniformdensity(vaf; depth::Int64 = 100, alt_reads::Int64 = 2, k::Int64 = 100, range = [0.02, 0.5], cutoff::Int64 = 1)
    if cutoff == 1
        f_min = (alt_reads/depth) + ( ( 2*sqrt(alt_reads*(1-(alt_reads/depth))) ) / depth)
        h = np.histogram(vaf[(range[2] .> vaf .> f_min)], range = range, bins = k)
        nd = h[1]
    else
        h = np.histogram(vaf[(range[2] .> vaf .> range[1])], range = range, bins = k)
        nd = h[1]
    end
    return nd
end

"""
Normalizes a vector by dividing all entries by the maximum value in the vector.
"""
function rescale!(arr)
    mx = maximum(arr)
    arr = arr ./ mx
    return arr
end

"""
Transformation of the VAF information into normalized feature vectors.

Input:
    - CanEvolve simulation output: A single simulation object e.g. p from p, n = autoSimulation(1, 0.1)
    - lower_cutoff (float): Lower frequency (range 0 - 0.5) bound for detecting subclones (default = 0.09)
    - upper_cutoff (float): Upper frequency (range 0 - 0.5) bound for detecting subclones (default = 0.41)
    - k (array): Number of bins in histogram feature vector. Multiple entries of different sizes will generate multiple feature vectors (default = [100])

Returns:
    - array: Features
    - array: Labels (mode, nsubclones, subclone frequency(ies), subclone emergence time, subclone fitness, subclone fitness relative to population average, mutation rate, mean depth, sequencing overdispersion)
"""
function engineer(sim_object; lower_cutoff::Float64 = 0.09, upper_cutoff::Float64 = 0.41, k::Array{Int64} = [100])

    VAF, depth, alt_reads = sim_object[2][1], Int64(round(mean(sim_object[2][3]))), sim_object[5][end]

    features = []
    for bin_number in k
        ud1 = rescale!(uniformdensity(VAF, k = bin_number, alt_reads = alt_reads, range = [0.02, 0.5]))
        push!(features, ud1)
    end

    mode = ifelse(sum(sim_object[3][1]) == 0.0, 0, 1)
    if mode == 1
        detectable = intersect(findall(2*upper_cutoff .> sim_object[3][2] .> 2*lower_cutoff), findall(upper_cutoff .> sim_object[3][1] .> lower_cutoff))
        labels = [1, length(sim_object[3][1][detectable]), sim_object[3][1][detectable], sim_object[3][2][detectable], sim_object[3][4][detectable], sim_object[3][5][detectable], sim_object[5][3], depth, sim_object[5][10]]
    elseif mode == 0
        labels = [0, 0, 0, 0, 0, 0, sim_object[5][3], depth, sim_object[5][10]]
    else
        println("Did not identify proper mode.")
    end

    return features, labels
end

"""
Transformation of the VAF information into nonnormalized feature vectors.

Input:
    - CanEvolve simulation output: A single simulation object e.g. p from p, n = autoSimulation(1, 0.1)
    - lower_cutoff (float): Lower frequency (range 0 - 0.5) bound for detecting subclones (default = 0.09)
    - upper_cutoff (float): Upper frequency (range 0 - 0.5) bound for detecting subclones (default = 0.41)
    - k (array): Number of bins in histogram feature vector. Multiple entries of different sizes will generate multiple feature vectors (default = [100])

Returns:
    - array: Features
    - array: Labels (mode, nsubclones, subclone frequency(ies), subclone emergence time, subclone fitness, subclone fitness relative to population average, mutation rate, mean depth, sequencing overdispersion)
"""
function engineer_un(sim_object; lower_cutoff::Float64 = 0.09, upper_cutoff::Float64 = 0.41, k::Array{Int64} = [100])

    VAF, depth, alt_reads = sim_object[2][1], Int64(round(mean(sim_object[2][3]))), sim_object[5][end]
    dp = sim_object[5][8]
    features = []
    for bin_number in k
        ud1 = uniformdensity(VAF, k = bin_number, depth = dp, alt_reads = alt_reads, range = [0.02, 0.5], cutoff = 1)
        if length(k) == 1
            append!(features, ud1)
        else
            push!(features, ud1)
        end
    end

    mode = ifelse(sum(sim_object[3][1]) == 0.0, 0, 1)
    if mode == 1
        detectable = intersect(findall(2*upper_cutoff .> sim_object[3][2] .> 2*lower_cutoff), findall(upper_cutoff .> sim_object[3][1] .> lower_cutoff))
        labels = [1, length(sim_object[3][1][detectable]), sim_object[3][1][detectable], sim_object[3][3][detectable], sim_object[3][4][detectable], sim_object[3][5][detectable], sim_object[5][3], depth, sim_object[5][10]]
    elseif mode == 0
        labels = [0, 0, 0, 0, 0, 0, sim_object[5][3], depth, sim_object[5][10]]
    else
        println("Did not identify proper mode.")
    end

    return features, labels
end

"""
A replicate of autoSimulation but with fixed mutation rate sequencing depth and overdispersion (rho) parameters.

See autoSimulation for a full description of input parameters.
"""
function autoSimulation_fixed(b, d; u::Int64 = 20, depth::Int64 = 100, Nfinal::Int64 = 1000, noise::String = "betabinomial", rho::Float64 = 0.0, lower_cutoff::Float64 = 0.1, upper_cutoff::Float64 = 0.4, nsubclone_min::Int64 = 1, nsubclone_max::Int64 = 2, alt_reads::Int64 = 2)

    # Run simulation
    _, phi_b, nclonal, ndrivers, lambda, _ = sampleParameters()
    println("{A} Starting positive simulation")
    println("--------------------------------")
    println("")
    n_detectable, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters = 0, [], [], [], []
    counter = 0
    while true
        counter += 1
        if counter == 10
            _, phi_b, nclonal, ndrivers, lambda, _ = sampleParameters()
            counter = 0
            println("Parameters failed to generate subclone number of interest, resampling parameters...")
        end
        n_detectable, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters = simulation(b, d, u, Nfinal = Nfinal, phi_b = phi_b, phi_d = 0.0, nclonal = nclonal, depth = depth, noise = noise, rho = rho, lower_cutoff = lower_cutoff, upper_cutoff = upper_cutoff, ndrivers = ndrivers, lambda = lambda, alt_reads = alt_reads)
        if nsubclone_max >= n_detectable >= nsubclone_min
            println("")
            println("Correct number of subclones generated for synthetic tumour with positive selection...")
            println("")
            break
        end
        println("Non-specified number of subclones generated, re-simulating positive selection...")
    end

    println("{B} Generating parameter paired fully synthetic neutral data")
    println("------------------------------")
    println("")

    scale_pos, scale_clonal = rand(Uniform(1, 1.5)), rand(Uniform(0.5, 1.5)) # Add variability to clonal and subclonal mutation counts in neutral synthetic tumours
    trim_tail = ifelse(rand(Uniform(0, 1), 1)[1] < 0.05, 1, 0) # Randomly remove tail from neutral distributions to inject more sparsity in neutral samples
    n_detectable, neutral_VAF, neutral_haplotypes, neutral_mutations, neutral_parameters = full_synthetic_neutral(depth = depth, Nfinal = Nfinal, nclonal = Int(round(nclonal * scale_clonal)), match_pos_muts = Int(round(positive_mutations[5] * scale_pos)), noise = noise, rho = rho, alt_reads = alt_reads, trim_tail = trim_tail)

    println("Saving synthetic data...")
    # Add identifiers for future mapping if needed
    simid = randstring(10)
    positive_id, neutral_id = join(["P_", simid]), join(["N_", simid])

    return [positive_id, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters], [neutral_id, neutral_VAF, neutral_haplotypes, neutral_mutations, neutral_parameters]
end

"""
A replicate of autoSimulation but with fixed mutation rate and sequencing depth parameters.

See autoSimulation for a full description of input parameters.
"""
function autoSimulation_fixed_depth(b, d; depth::Int64 = 100, Nfinal::Int64 = 1000, noise::String = "betabinomial", lower_cutoff::Float64 = 0.1, upper_cutoff::Float64 = 0.4, nsubclone_min::Int64 = 1, nsubclone_max::Int64 = 2, alt_reads::Int64 = 2)

    # Run simulation
    u, phi_b, nclonal, ndrivers, lambda, _, rho = sampleParameters()
    println("{A} Starting positive simulation")
    println("--------------------------------")
    println("")
    n_detectable, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters = 0, [], [], [], []
    counter = 0
    while true
        counter += 1
        if counter == 10
            _, phi_b, nclonal, ndrivers, lambda, _ = sampleParameters()
            counter = 0
            println("Parameters failed to generate subclone number of interest, resampling parameters...")
        end
        n_detectable, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters = simulation(b, d, u, Nfinal = Nfinal, phi_b = phi_b, phi_d = 0.0, nclonal = nclonal, depth = depth, noise = noise, rho = rho, lower_cutoff = lower_cutoff, upper_cutoff = upper_cutoff, ndrivers = ndrivers, lambda = lambda, alt_reads = alt_reads)
        if nsubclone_max >= n_detectable >= nsubclone_min
            println("")
            println("Correct number of subclones generated for synthetic tumour with positive selection...")
            println("")
            break
        end
        println("Non-specified number of subclones generated, re-simulating positive selection...")
    end

    println("{B} Generating parameter paired fully synthetic neutral data")
    println("------------------------------")
    println("")

    scale_pos, scale_clonal = rand(Uniform(1, 1.5)), rand(Uniform(0.5, 1.5)) # Add variability to clonal and subclonal mutation counts in neutral synthetic tumours
    trim_tail = ifelse(rand(Uniform(0, 1), 1)[1] < 0.05, 1, 0) # Randomly remove tail from neutral distributions to inject more sparsity in neutral samples
    n_detectable, neutral_VAF, neutral_haplotypes, neutral_mutations, neutral_parameters = full_synthetic_neutral(depth = depth, Nfinal = Nfinal, nclonal = Int(round(nclonal * scale_clonal)), match_pos_muts = Int(round(positive_mutations[5] * scale_pos)), noise = noise, rho = rho, alt_reads = alt_reads, trim_tail = trim_tail)

    println("Saving synthetic data...")
    # Add identifiers for future mapping if needed
    simid = randstring(10)
    positive_id, neutral_id = join(["P_", simid]), join(["N_", simid])

    return [positive_id, positive_VAF, positive_haplotypes, positive_mutations, positive_parameters], [neutral_id, neutral_VAF, neutral_haplotypes, neutral_mutations, neutral_parameters]
end


text:
#=
    this splitting criterion considers the pointwise ratio = numerator./denominator
    we minimize the variance(ratios in left child) + variance(ratios in right child) =#

# calculateSplitValue(::DecisionTrees.msePointwiseSplit, ::Symbol, ::Int64, ::Array{UInt8,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64, ::DecisionTrees.MyIncreasingSubsets, ::Array{OnlineStats.Series,1})
function calculateSplitValue(a::msePointwiseSplit, fname::Symbol, number_of_char_features::Int, labellist::Vector{T}, sumnumerator::Array{Float64,1}, sumdenominator::Array{Float64,1}, sumweight::Array{Float64,1}, countlistfloat::Array{Float64,1}, minweight::Float64, subs::DTSubsets, moments_per_pdaclass) where T <: Unsigned
# here randomweight==0
# for subsets, exhaustive search with flipping members (gray code) or "increasing" subset search ({1}, {1,2}, {1,2,3}, .... {1,2,3, ....., n-1,2})
# all input lists (labellist,sumnumerator,sumdenominator,sumweight,countlistfloat) need to be sorted in the same manner
# convention: in the beginning everything is on the right side
    elementsInLeftChildBV = BitVector(undef, length(labellist));fill!(elementsInLeftChildBV, false) # indicates which classes belong to right child
    chosen_subset_bitarray = BitVector(undef, 0)

    totalMeanAndVariance = copy(emptyMeanVarSeries)
    for idx = 1:length(moments_per_pdaclass)
        @inbounds  mi = moments_per_pdaclass[idx]
        @inbounds statsI = mi.stats[1]
        if statsI.n > 0 # this is needed because of https://github.com/joshday/OnlineStats.jl/issues/125
            @inbounds merge!(totalMeanAndVariance, mi)
        end
    end

    weighttot = sum(sumweight)
    weighttot_minw = weighttot - minweight

    momentsl = copy(emptyMeanVarSeries) # Series(OnlineStats.Mean(),OnlineStats.Variance())
# momentsr needs to be the 'total' at the start of this.
    momentsr = totalMeanAndVariance # copy(emptyMeanVarSeries) #Series(OnlineStats.Mean(),OnlineStats.Variance())

    sumwl = 0.0
    varl = varr = val = -Inf
    chosen_sumwl = NaN
    for i in subs
  # @show "vv=$(i)" #should be the element which switches
  # i is an index, it indicates which element of labellist will flip sides for the next calculation
        @inbounds switching_class = moments_per_pdaclass[i]
        @inbounds if elementsInLeftChildBV[i]
    # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds  elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true) # updating needs to occur here before we copy the array (in case it is a better split than what we have seen so far)
    # move class from left to right side
            unmerge!(momentsl, switching_class) # switching_class is 'removed' from the left variance/mean statistcs momentsl
            merge!(momentsr, switching_class) # switching_class is 'added' to the right variance/mean statistcs momentsl
    # sumnl-=sumnumerator[i]
    # sumdl-=sumdenominator[i]
            @inbounds  sumwl -= sumweight[i]
        else
      # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds  elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true)
      # move class from right to left side
            merge!(momentsl, switching_class)
            unmerge!(momentsr, switching_class)
      # sumnl+=sumnumerator[i]
    	# sumdl+=sumdenominator[i]
    	@inbounds  sumwl += sumweight[i]
        end
        if (sumwl > minweight) && (weighttot_minw > sumwl) # do we have enough exposure? is the split valid?
	  # vold=valnew
            @inbounds varl = momentsl.stats[2].σ2
            @inbounds varr = momentsr.stats[2].σ2
	  # @show fname,sumwl,weighttot-sumwl
            valnew = -(varl + varr)
      # if vold==valnew;@show i,valnew,sumnl,sumdl;end;
	  if valnew > val
                val = valnew
                chosen_subset_bitarray = copy(elementsInLeftChildBV)
                chosen_sumwl = sumwl
            end
        end

    end
    if isfinite(val)
  # warning: if the labellist is not sorted 1:n we need to be careful here!
        chosen_subset = labellist[chosen_subset_bitarray]
    else
        chosen_subset = Array{UInt}(undef, 0)
    end
    return val, chosen_subset, chosen_sumwl, weighttot - chosen_sumwl
end

function calculateSplitValue(a::msePointwiseSplit, number_of_char_features::Int, labellist::Array{UInt8,1}, sumnumerator::Array{Float64,1}, sumdenominator::Array{Float64,1}, sumweight::Array{Float64,1}, countlistfloat::Array{Float64,1}, minweight::Float64, subs::DTSubsets, feature_column_id::Int, moments_per_pdaclass)
    @warn("221 this will not work yet.....")
# here randomweight>0
# for subsets, exhaustive search with flipping members (gray code) or "increasing" subset search ({1}, {1,2}, {1,2,3}, .... {1,2,3, ....., n-1,2})
# all input lists (labellist,sumnumerator,sumdenominator,sumweight,countlistfloat) need to be sorted in the same manner
# convention, in the beginning everything is on the right side
    elementsInLeftChildBV = BitVector(undef, length(labellist));fill!(elementsInLeftChildBV, false) # indicates which classes belong to right child

    totalMeanAndVariance = copy(emptyMeanVarSeries)
    for idx = 1:length(moments_per_pdaclass)
        @inbounds  mi = moments_per_pdaclass[idx]
        @inbounds statsI = mi.stats[1]
        if statsI.n > 0 # this is needed because of https://github.com/joshday/OnlineStats.jl/issues/125
            @inbounds merge!(totalMeanAndVariance, mi)
        end
    end

    weighttot = sum(sumweight)
    weighttot_minw = weighttot - minweight

    momentsl = copy(emptyMeanVarSeries) # Series(OnlineStats.Mean(),OnlineStats.Variance())
# momentsr needs to be the 'total' at the start of this.
    momentsr = totalMeanAndVariance # copy(emptyMeanVarSeries) #Series(OnlineStats.Mean(),OnlineStats.Variance())

    sumwl = 0.0
    this_splitlist = Array{Splitdef}(undef, 0)

    for i in subs
    # i is an index, it indicates which element of labellist will flip sides for the next calculation
        if elementsInLeftChildBV[i]
    # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true) # updating needs to occur here before we copy the array (in case it is a better split than what we have seen so far)
    # move class from left to right side
            unmerge!(momentsl, switching_class)
            merge!(momentsr, switching_class)
	@inbounds   sumwl -= sumweight[i]
        else
    # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true)
    # move class from right to left side
            merge!(momentsl, switching_class)
            unmerge!(momentsr, switching_class)
	@inbounds   sumwl += sumweight[i]
        end
        if (sumwl > minweight) && (weighttot_minw > sumwl) # do we have enough exposure? is the split valid?
            @inbounds varl = momentsl.stats[2].σ2
            @inbounds varr = momentsr.stats[2].σ2
	  	valnew = -(varl + varr)
            @inbounds push!(this_splitlist, Splitdef{T}(feature_column_id, feature_column_id2, fname, labellist[elementsInLeftChildBV], valnew, sumwl, weighttot - sumwl))
        end
    end
    return this_splitlist
end

function aggregate_data_msePointwise(f, numerator::Array{Float64,1}, denominator::Array{Float64,1}, weight::Array{Float64,1})
    # this is the core function of the modelling process
    # besides copying of the data, the vast majority of time is spent in here!
    # most of the time is spent here, if we can improve the for loop below, that would improve performance greatly!
    # one possibility would be to introduce parallelization here (which is not straightforward, I think....)
    lo = one(eltype(f.parent.refs))
    hi = convert(eltype(f.parent.refs), length(f.parent.pool))
    ooo = one(lo) - lo
    vecsize = hi + ooo
    cnt = zeros(Int, vecsize)
    sumnumerator = zeros(Float64, vecsize)
    sumdenominator = zeros(Float64, vecsize)
    sumweight = zeros(Float64, vecsize)
    moments_per_pdaclass = Array{OnlineStats.Series}(undef, vecsize)
    for i = 1:vecsize
        @inbounds moments_per_pdaclass[i] = copy(emptyMeanVarSeries)
    end
    @inbounds for count in f.indices[1]
        @inbounds idx = f.parent.refs[count] + ooo
        @inbounds cnt[idx] += 1
        @inbounds this_num = numerator[count]
        @inbounds sumnumerator[idx] += this_num
        @inbounds this_denom = denominator[count]
        @inbounds sumdenominator[idx] += this_denom
        @inbounds sumweight[idx] += weight[count]
        @inbounds OnlineStats.fit!(moments_per_pdaclass[idx], this_num / this_denom)
    end
    return cnt, sumnumerator, sumdenominator, sumweight, moments_per_pdaclass
end

"""
function unmerge!(s1::T, s2::T) where {T <: OnlineStats.Series}
this function updates the mean and variance of s1
The assumption is that the underlying data of s2 is a subset of s1.
If this is not the case, the result of this function is most likely meaningless
That is to say "subtract" the series s2 from the series s1 and update the mean and variance of s1 accordingly
both s1 and s2 must be of the 'form' Series(OnlineStats.Mean(),OnlineStats.Variance())
Mean must be the fist statistic and Variance the second
Example
vec1=[1.2,3.4,99.1];
vec2=[11.1,3.2];
vecCombined=vcat(vec1,vec2);
#@show isapprox(StatsBase.var(vecCombined,correcte=true),thisVar)

Consider also the merge method: merge!{T <: Series}(s1::T, s2::T, method::Symbol = :append) of OnlineStats
https://github.com/joshday/OnlineStats.jl
"""
 function unmerge!(s1::T, s2::T) where {T <: OnlineStats.Series}
    n2 = OnlineStats.nobs(s2)
    n2 == 0 && return s1
    # OnlineStats.updatecounter!(s1, -n2)
    s1.stats[1].n -= n2
    s1.stats[2].n -= n2

    sizetot = OnlineStats.nobs(s1) + n2
    sizewo = sizetot - n2
    gamma = n2 / sizetot
    # we currently access the statistics (i.e. variance and mean) by index which is suboptimal
    # this will throw a typeassert error if the index is not right or the series consists of different types
     o = s1.stats[2]::OnlineStats.Variance # WARNING if this is not of type variance, then all this will fail miserably
     o2 = s2.stats[2]::OnlineStats.Variance
    # update variance
        one_minus_gamma = 1.0 - gamma
        one_minus_gamma_inverted = 1 / one_minus_gamma
        muwo = o.μ = 1 / one_minus_gamma * (o.μ - gamma * o2.μ)
        o.σ2 = one_minus_gamma_inverted * (o.σ2 - gamma * (o2.σ2 + (o2.μ - muwo)^2 * one_minus_gamma))
        o.μ = muwo # note the variance "struct" has a mu parameter as well (in addition to the mean of the series)
    # update mean
    s1.stats[1].μ = muwo
    return s1
end


text:
# example 6: commands & subtables

using ArgParse

function main(args)

    s = ArgParseSettings("Example 6 for argparse.jl: " *
                         "commands and their associated sub-tables.")

    @add_arg_table s begin
        "run"
            action = :command        # adds a command which will be read from an argument
            help = "start running mode"
        "jump"
            action = :command
            help = "start jumping mode"
    end

    @add_arg_table s["run"] begin    # add command arg_table: same as usual, but invoked on s["cmd"]
        "--speed"
            arg_type = Float64
            default = 10.
            help = "running speed, in Å/month"
    end

    s["jump"].description = "Jump mode for example 6"  # this is how settings are tweaked
                                                       # for commands
    s["jump"].commands_are_required = false            # this makes the sub-commands optional
    s["jump"].autofix_names = true                     # this uses dashes in long options, underscores
                                                       # in auto-generated dest_names

    @add_arg_table s["jump"] begin
        "--higher"
            action = :store_true
            help = "enhance jumping"
        "--somersault"
            action = :command        # this adds a sub-command (passed via an option instead)
            dest_name = "som"        # flag commands can set a dest_name
            help = "somersault jumping mode"
        "--clap-feet"                # dest_name will be "clap_feet" (see the "autofix_names" settings")
            action = :command
            help = "clap feet jumping mode"
    end

    s["jump"]["som"].description = "Somersault jump " *  # this is how settings are tweaked
                                   "mode for example 6"  # for sub-commands

    s["jump"]["clap_feet"].description = "Clap-feet jump " *  # notice the underscore in the name
                                         "mode for example 6"

    parsed_args = parse_args(args, s)
    println("Parsed args:")
    for (key,val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end
    println()

    # parsed_args will have a special field "%COMMAND%"
    # which will hold the executed command name (or 'nothing')
    println("Command: ", parsed_args["%COMMAND%"])

    # thus, the command args are in parsed_args[parsed_args["%COMMAND%]]
    println("Parsed command args:")
    command_args = parsed_args[parsed_args["%COMMAND%"]]
    for (key,val) in command_args
        println("  $key  =>  $(repr(val))")
    end
end

main(ARGS)


text:
#
# Empty storage
#

struct Empty{ElT, StoreT <: TensorStorage} <: TensorStorage{ElT}
end

# Get the Empty version of the TensorStorage
function empty(::Type{StoreT}) where {StoreT <: TensorStorage{ElT}} where {ElT}
  return Empty{ElT, StoreT}
end

# Defaults to Dense
function Empty(::Type{ElT}) where {ElT}
  return empty(Dense{ElT, Vector{ElT}})()
end

Empty() = Empty(Float64)

Base.copy(S::Empty) = S

Base.isempty(::Empty) = true

nnzblocks(::Empty) = 0

nnz(::Empty) = 0

function Base.complex(::Type{<: Empty{ElT, StoreT}}) where {ElT,
                                                            StoreT}
  return Empty{complex(ElT), complex(StoreT)}
end

function Base.complex(S::Empty)
  return complex(typeof(S))()
end

#Base.size(::Empty) = 0

function Base.show(io::IO,
                   mime::MIME"text/plain",
                   S::Empty)
  println(io, typeof(S))
end

#
# EmptyTensor (Tensor using Empty storage)
#

const EmptyTensor{ElT,
                  N,
                  StoreT,
                  IndsT} = Tensor{ElT,
                                  N,
                                  StoreT,
                                  IndsT} where {StoreT <: Empty}

Base.isempty(::EmptyTensor) = true

function EmptyTensor(::Type{ElT}, inds) where {ElT <: Number}
  return tensor(Empty(ElT), inds)
end

function EmptyTensor(::Type{StoreT}, inds) where {StoreT <: TensorStorage}
  return tensor(empty(StoreT)(), inds)
end

function EmptyBlockSparseTensor(::Type{ElT}, inds) where {ElT <: Number}
  StoreT = BlockSparse{ElT, Vector{ElT}, length(inds)}
  return EmptyTensor(StoreT, inds)
end

# From an EmptyTensor, return the closest Tensor type
function Base.fill(::Type{<: Tensor{ElT,
                                    N,
                                    EStoreT,
                                    IndsT}}) where {ElT <: Number,
                                                    N,
                                                    EStoreT <: Empty{ElT,
                                                                     StoreT},
                                                    IndsT} where {StoreT}
  return Tensor{ElT, N, StoreT, IndsT}
end

function Base.zeros(T::TensorT) where {TensorT <: EmptyTensor}
  TensorR = fill(TensorT)
  return zeros(TensorR, inds(T))
end

function addblock(T::EmptyTensor{<: Number, N},
                  block::Block{N}) where {N}
  R = zeros(T)
  addblock!(R, block)
  return R
end

addblock!!(T::EmptyTensor{<: Number, N},
           block::Block{N}) where {N} = addblock(T, block)

Base.@propagate_inbounds function Base.setindex(T::EmptyTensor{<: Number, N},
                                                x::Number,
                                                I::Vararg{Int, N}) where {N}
  R = zeros(T)
  R[I...] = x
  return R
end

function Base.setindex(T::EmptyTensor{<: Number, Any},
                                      x::Number,
                                      I::Int...)
  error("Setting element of EmptyTensor with Any number of dimensions not defined")
end

setindex!!(T::EmptyTensor,
           x::Number,
           I::Int...) = setindex(T, x, I...)

function Base.show(io::IO,
                   mime::MIME"text/plain",
                   T::EmptyTensor)
  summary(io, T)
  println(io)
end


text:
import StatsBase

"""
    VisitationFrequencyTest(k::Int = 1, l::Int = 1, m::Int = 1, n::Int = 1,
        τ::Int = 1, b = 2, estimator::VisitationFrequency = VisitationFrequency(),
        binning_summary_statistic::Function = StatsBase.mean,
        binning::RectangularBinning, ηs)

The parameters for a transfer entropy test using the `VisitationFrequency` estimator.
This is the original transfer entropy estimator from [1], as implemented in [2].

## Mandatory keyword arguments

- **`binning::RectangularBinning`**: An instance of a [`RectangularBinning`](@ref) that dictates
    how the delay embedding is discretized.

- **`ηs`**: The prediction lags (that go into the ``T_{f}`` component of the embedding).

## Optional keyword arguments

- **`k::Int`**: The dimension of the ``T_{f}`` component of the embedding.

- **`l::Int`**: The dimension of the ``T_{pp}`` component of the embedding.

- **`m::Int`**: The dimension of the ``S_{pp}`` component of the embedding.

- **`n::Int`**: The dimension of the ``C_{pp}`` component of the embedding.

- **`τ::Int`**: The embedding lag. Default is `τ = 1`.

- **`b`**: Base of the logarithm. The default (`b = 2`) gives the TE in bits.

- **`estimator::VisitationFrequency`**: A `VisitationFrequency` estimator instance.

- **`binning_summary_statistic::Function`**: A summary statistic to summarise the
    transfer entropy values if multiple binnings are provided.

## Estimation of the invariant measure

With a ` VisitationFrequencyTest`, the invariant distribution from which transfer
entropy is computed is estimated using simple counting (hence, visitation frequency).
Counting is done over the partition elements of a [discretization](@ref discretization)
of an appropriate [delay reconstruction](@ref custom_delay_reconstruction) of the time
series to be analysed.

## About the delay reconstruction for transfer entropy analysis

Denote the time series for the source process ``S`` as ``S(t)``, and the time series for
the target process ``T`` as ``T(t)``, and ``C_i(t)`` as the time series for any conditional
processes ``C_i`` that also may influence ``T``. To compute (conditional) TE, we need a
generalised embedding [3, 4] incorporating all of these processes.

For convenience, define the state vectors

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1))\\}, \\label{eq:Tf} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\\\
S_{pp}^{(m)} &= \\{(S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1}))\\},\\\\
C_{pp}^{(n)} &= \\{ (C_1(t), C_1(t-\\tau_1), \\ldots,  C_2(t), C_2(t-\\tau_1) \\},
\\end{align}
```

where the state vectors ``T_f^{(k)}`` contain ``k`` future values of the target
variable, ``T_{pp}^{(l)}`` contain ``l`` present and past values of the target
variable, ``S_{pp}^{(m)}`` contain ``m`` present and past values of the source
variable, ``C_{pp}^{(n)}`` contain a total of ``n`` present and past values of any
conditional variable(s).  Combining all variables, we have the generalised embedding

```math
\\begin{align}
\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)}, C_{pp}^{(n)})
\\end{align}
```

with a total embedding dimension of ``k + l + m + n``.
Here, ``\\tau`` indicates the [embedding](@ref custom_delay_reconstruction) lag
(in the `VisitationFrequencyTest`, we set ``\\tau_1 = \\tau_2 = \\tau_3 = \\ldots``), and
``\\eta`` indicates the prediction lag (the lag of the influence the source
has on the target).

Hence, in the generalised embedding, only ``T_f`` depends on the prediction lag ``\\eta``,
which is to be determined by the analyst. For transfer entropy analysis, ``\\eta`` is chosen
to be some positive integer, while for
[predictive asymmetry analysis](@ref predictive_asymmetry),
symmetric negative and positive ``\\eta``s are used for computing ``\\mathbb{A}``.

In terms of this generalised embedding, transfer entropy from a source variable ``S`` to a
target variable ``T`` with conditioning on variable(s) ``C`` is thus defined as

```math
\\begin{align}
TE_{S \\rightarrow T|C} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}, C_{pp}) \\log_{2}{\\left( \\frac{P(T_f | T_{pp}, S_{pp}, C_{pp})}{P(T_f | T_{pp}, C_{pp})}\\right)}
\\end{align}
```

Without conditioning, we have

```math
\\begin{align}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}) \\log_{2}{\\left(\\frac{P(T_f | T_{pp}, S_{pp})}{P(T_f | T_{pp})}\\right)}
\\end{align}
```

## Low-level estimator

This test uses the `VisitationFrequency` estimator on the following low-level method
under the hood.

- [`transferentropy(::Any, ::TEVars, ::RectangularBinning, ::TransferEntropyEstimator)`](@ref)

In this estimator, the mapping between variables of the
[generalised embedding](@ref custom_delay_reconstruction) and the marginals during
transfer entropy computation is controlled using a [`TEVars`](@ref)
instance. It is *highly* recommended that you check the documentation for this
method, because it describes the transfer entropy estimation procedure in detail.

## Notes:

- Use `causality(source, target, params::VisitationFrequencyTest)` for regular
    transfer entropy analysis. This method uses only the `k`, `l`, `m` and ignores `n`
    when constructing the delay reconstruction.

- Use `causality(source, target, cond, params::VisitationFrequencyTest)` for conditional
    transfer entropy analysis. This method uses the `k`, `l`, `m` *and* `n` when constructing
    the delay reconstruction.

## Example

```julia
# Prediction lags
ηs = 1:10
binning = RectangularBinning(10)

# Use defaults, binning and prediction lags are required.
# Note that `binning` and `ηs` are *mandatory* keyword arguments.
VisitationFrequencyTest(binning = binning, ηs = ηs)

# The other keywords can also be adjusted
VisitationFrequencyTest(k = 1, l = 2, binning = binning, ηs = ηs)
```

## References

1. Schreiber, Thomas. "Measuring information transfer." Physical review letters 85.2 (2000): 461.
    [https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.85.461](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.85.461)

2. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
    [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.99.042212](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.99.042212)

3.  Sauer, Tim, James A. Yorke, and Martin Casdagli. "Embedology." Journal of statistical Physics 65.3-4 (1991): 579-616.

4.  Deyle, Ethan R., and George Sugihara. "Generalized theorems for nonlinear state space reconstruction." PLoS One 6.3 (2011): e18295.
"""
Base.@kwdef struct VisitationFrequencyTest{N} <: TransferEntropyCausalityTest{N}

    """ The delay reconstruction parameter k (controls dimension of ``T_{f}`` component of embedding). """
    k::Int = 1

    """ The delay reconstruction parameter l (controls dimension of ``T_{pp}`` component of embedding). """
    l::Int = 1

    """ The delay reconstruction parameter m (controls dimension of ``S_{pp}`` component of embedding). """
    m::Int = 1

    """ The delay reconstruction parameter n (controls dimension of ``C_{pp}`` component of embedding). """
    n::Int = 1

    """ The delay reconstruction lag for the ``T_{pp}`` component of the embedding. """
    τ::Int = 1

    """ The base of the logarithm for computing TE. """
    b = 2

    """ The transfer entropy estimator. """
    estimator::VisitationFrequency = VisitationFrequency()

    """
    If there are several binnings provided, what is the statistic used to summarise the
    transfer entropy values to a single value?
    """
    binning_summary_statistic::Function = StatsBase.mean

    """ The binning scheme. """
    binning::Union{RectangularBinning, Vector{RectangularBinning}}

    """ The prediction lags"""
    ηs

    function VisitationFrequencyTest(k::Int, l::Int, m::Int, n::Int, τ::Int, b,
            estimator::VisitationFrequency,
            binning_summary_statistic::Function,
            binning::Union{RectangularBinning, Vector{RectangularBinning}},
            ηs)

        N = length(ηs) # length of return vector when used with `causality`
        return new{N}(k, l, m, n, τ, b, estimator, binning_summary_statistic, binning, ηs)
    end
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, p::VisitationFrequencyTest) where {T <: Real}
    [p.binning_summary_statistic(
        transferentropy(source, target, p.binning,
            p.k, p.l, p.m, η = η, τ = p.τ,
            estimator = p.estimator)) for η in p.ηs]
end

export VisitationFrequencyTest


text:
# debug script

pos = particles[2].x

for particle in particles
    vi = vinduced(particle,pos)
    println(vi)
end
# passed

for particle in particles
    vi = netv(particles, particle.x) .+ Uinf(particle.x)
    println(vi)
end
# passed

for particle in particles
    println(netv(particles, particle.x) .+ Uinf(particle.x))
end
# passed

timestep = 0.01
for particle in particles
    newx = particle.x + particle.v * timestep
    println(newx)
end
# passed

for particle in particles
    particle.v .= netv(particles, particle.x) .+ Uinf(particle.x)
end
# DID NOT PASS

function advance(particles::Array{Particle,1},timestep::Float64,Uinf)
    # update velocity of each particle
    for particle in particles
        particle.v .= netv(particles, particle.x) .+ Uinf(particle.x)
    end

    # move each particle
    for particle in particles
        particle.x .= particle.x + particle.v * timestep
    end
end


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

mutable struct IoK8sApiAuthorizationV1beta1ResourceAttributes <: SwaggerModel
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace
    resource::Any # spec type: Union{ Nothing, String } # spec name: resource
    subresource::Any # spec type: Union{ Nothing, String } # spec name: subresource
    verb::Any # spec type: Union{ Nothing, String } # spec name: verb
    version::Any # spec type: Union{ Nothing, String } # spec name: version

    function IoK8sApiAuthorizationV1beta1ResourceAttributes(;group=nothing, name=nothing, namespace=nothing, resource=nothing, subresource=nothing, verb=nothing, version=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("resource"), resource)
        setfield!(o, Symbol("resource"), resource)
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("subresource"), subresource)
        setfield!(o, Symbol("subresource"), subresource)
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("verb"), verb)
        setfield!(o, Symbol("verb"), verb)
        validate_property(IoK8sApiAuthorizationV1beta1ResourceAttributes, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        o
    end
end # type IoK8sApiAuthorizationV1beta1ResourceAttributes

const _property_map_IoK8sApiAuthorizationV1beta1ResourceAttributes = Dict{Symbol,Symbol}(Symbol("group")=>Symbol("group"), Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"), Symbol("resource")=>Symbol("resource"), Symbol("subresource")=>Symbol("subresource"), Symbol("verb")=>Symbol("verb"), Symbol("version")=>Symbol("version"))
const _property_types_IoK8sApiAuthorizationV1beta1ResourceAttributes = Dict{Symbol,String}(Symbol("group")=>"String", Symbol("name")=>"String", Symbol("namespace")=>"String", Symbol("resource")=>"String", Symbol("subresource")=>"String", Symbol("verb")=>"String", Symbol("version")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1ResourceAttributes }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1ResourceAttributes))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1ResourceAttributes }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1ResourceAttributes[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1ResourceAttributes }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1ResourceAttributes[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1ResourceAttributes)
    true
end

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


text:
export regularize

"derivative of regularizar"
function regularize(W::Array{Float64})
  reg=W # L1:sign(W)
end


text:
module Vlasiator

using Requires
using StaticArrays: SVector, @SVector, SMatrix, @SMatrix
using Printf: @sprintf
using LinearAlgebra: ×, dot, ⋅, norm, normalize, normalize!
using Statistics: mean
using EzXML
using Mmap: mmap
using WriteVTK
using LazyGrids: ndgrid
using LaTeXStrings
using Dates

include("utility/rotation.jl")
include("utility/log.jl")
include("utility/monitor.jl")
include("utility/array.jl")
include("utility/fluxfunction.jl")
include("vlsv/vlsvreader.jl")
include("vlsv/vlsvutility.jl")
include("utility/plot.jl")

export
   # vlsvreader
   MetaVLSV, VarInfo,
   load, readvariable, readparameter, readvariablemeta, readvcells,
   hasvariable, hasparameter, hasname,
   # vlsvutility
   getcell, getslicecell, getlevel, refineslice, getcellcoordinates, getvcellcoordinates,
   getchildren, getparent, isparent, getsiblings,
   getcellinline, getnearestcellwithvdf, getcellwithvdf,
   getdensity, getvelocity, getpressure, getmaxwellianity, write_vtk, write_vlsv, issame,
   # plot helper
   SI, EARTH, Log, Linear, SymLog,
   # log
   readlog,
   # physical parameter monitor
   check_plasma_characteristics,
   # fluxfunction
   compute_flux_function, find_reconnection_points

precompile(load, (String,))
precompile(readvariable, (MetaVLSV, String, Bool))

function __init__()
   @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
      include("plot/pyplot.jl")
   end
   @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
      include("plot/plots.jl")
   end
end

end


text:
type PageSize
	xscale::Float64
	yscale::Float64
	units::String
end

type PlotFrame
	wi::Int
	name::String
	size::PageSize
	title::String
	xlabel::String
	ylabel::String
	xlim::Vector{Float64}
	ylim::Vector{Float64}
	legendInside::Bool
	showLegend::Bool
	legendPos::String
	topAxis::Bool
	leftAxis::Bool
	bottomAxis::Bool
	rightAxis::Bool
	wpx::Float64
	hpx::Float64
	abxb::Float64
	abxe::Float64
	aby::Float64
	oryb::Float64
	orye::Float64
	orx::Float64
	plen::Float64
	phei::Float64
end

function PlotFrame(name::AbstractString, title, xlabel, ylabel, width, height)
	w = p2d.screen(width, height)
	p2d.createFont(w, "Arial", 13*(width/500))
	p2d.background(w, 1.0, 1.0, 1.0, 1.0)

	abxb = 0.18*width; abxe = 0.96*width
	orx = 0.16*width; aby = 0.84*height
	orye = 0.055*height

	plen = abxe - abxb
	phei = aby - orye

	drawAxisLabels(w, title, xlabel, ylabel, width, height, abxb, abxe, aby, orx, orye, plen, phei)

	return PlotFrame(w, name, PageSize(9, 9, "cm"), title, xlabel, ylabel, [NaN, NaN], [NaN, NaN], true, true, "tl", false, true, true, false, width, height, abxb, abxe, aby, aby, orye, orx, plen, phei)
end

PlotFrame(width, height) = PlotFrame("sparrow", "", "x", "y", width, height)
PlotFrame(s) = PlotFrame("sparrow", "", "x", "y", s, s)
PlotFrame() = PlotFrame("sparrow", "", "x", "y", 400, 400)


text:
using GeometryTypes, Visualize, FileIO
using ModernGL, StaticArrays
using GLAbstraction, GLWindow, ModernGL, GeometryTypes
using GLAbstraction: compile_shader
using Visualize: Window, WindowEvents, Mouse, Keyboard, PerspectiveCamera, add!
using Visualize: Area, Projection, View, EyePosition, LookAt, UpVector, Translation
using Visualize: Rotation, ProjectionView, TranslationSpeed, RotationSpeed
import GLFW

function read_ifs(filename)
    open(filename) do io
        function str()
            n = read(io, UInt32)
            String(read(io, UInt8, n))
        end
        ifs = str()
        zero = read(io, Float32)
        name = str()
        vertices = str()
        nverts = read(io, UInt32)
        verts = Vector{Point3f0}(nverts)
        for i = 1:nverts
            verts[i] = Point3f0(
                read(io, Float32),
                read(io, Float32),
                read(io, Float32)
            )
        end
        tris = str()
        nfaces = read(io, UInt32)
        faces = Vector{GLTriangle}(nfaces)
        for i = 1:nfaces
            faces[i] = GLTriangle(
                read(io, UInt32) + UInt32(1),
                read(io, UInt32) + UInt32(1),
                read(io, UInt32) + UInt32(1)
            )
        end
        nv = normals(verts, faces)
        view(collect(zip(verts, nv)), faces)
    end
end

vert = """
#version 450
#extension GL_ARB_shader_draw_parameters : enable

struct VertexArgument{
    vec4 color;
    mat4 model;
};

layout (location = 0) in vec3 position;
layout (location = 1) in vec3 normal;

layout (std140) uniform Scene{
    vec4 lightposition;
    mat4 proj;
    mat4 view;
    mat4 projview;
    vec2 resolution;
} scene;

layout (std140) uniform VertexArguments{
    VertexArgument[2048] args;
} vertex_arguments;

out VertexOut{
    vec3 vertex;
    vec3 normal;
    vec3 lightdir;
    vec4 color;
} vertex_out;

void main(){
    VertexArgument arg = vertex_arguments.args[gl_DrawIDARB];
    vec4 position_camspace = scene.view * arg.model * vec4(position,  1.0);
    gl_Position = scene.proj * position_camspace;
    vertex_out.lightdir = normalize(vec3(-10) - position.xyz);
    vertex_out.vertex = -position_camspace.xyz;
    vertex_out.normal = normal;
    vertex_out.color = arg.color;
}
"""

frag = """
#version 450

vec3 blinnphong(vec3 V, vec3 N, vec3 L, vec3 color){

    float diff_coeff = max(dot(L,N), 0.0);

    // specular coefficient
    vec3 H = normalize(L+V);

    float spec_coeff = pow(max(dot(H,N), 0.0), 8.0);
    if (diff_coeff <= 0.0)
        spec_coeff = 0.0;

    // final lighting model
    return vec3(
        vec3(0.1) * vec3(0.3)  +
        vec3(0.9) * color * diff_coeff +
        vec3(0.3) * spec_coeff
    );
}

in VertexOut{
    vec3 vertex;
    vec3 normal;
    vec3 lightdir;
    vec4 color;
} vertex_in;

layout (location = 0) out vec4 frag_color;

void main(){
    vec3 L = normalize(vertex_in.lightdir);
    vec3 N = normalize(vertex_in.normal);
    vec3 light1 = blinnphong(vertex_in.vertex, N, L, vertex_in.color.rgb);
    vec3 light2 = blinnphong(vertex_in.vertex, N, -L, vertex_in.color.rgb);
    frag_color = vec4(light1 + light2, 1.0);
}
"""

window = create_glcontext(
    major = 4, minor = 5, debugging = false,
    windowhints = [
        (GLFW.SAMPLES,      0),
        (GLFW.DEPTH_BITS,   32),

        (GLFW.ALPHA_BITS,   8),
        (GLFW.RED_BITS,     8),
        (GLFW.GREEN_BITS,   8),
        (GLFW.BLUE_BITS,    8),

        (GLFW.STENCIL_BITS, 0),
        (GLFW.AUX_BUFFERS,  0)
    ]
)

events = WindowEvents(Window => window)

cam = PerspectiveCamera(
    TranslationSpeed => 1f0,
    LookAt => Vec3f0(0),
    EyePosition => Vec3f0(6, 6, 8),
    Rotation => Vec3f0(0),
    Area => events[Area],
    RotationSpeed => 0.1f0
)

for elem in (Mouse.Position, Mouse.Buttons, Mouse.Drag, Area, Keyboard.Buttons, Mouse.Scroll)
    add!(events, elem)
end

add!(cam, Translation, events,
    Mouse.right, Set([Keyboard.left_control])
)
add!(cam, Rotation, events,
    Mouse.left, Set([Keyboard.left_control])
)
add!(cam, Visualize.Projection, View)
add!(cam, EyePosition, LookAt)

cam[LookAt] = Vec3f0(0)
cam[EyePosition] = Vec3f0(3)
vertshader = compile_shader(Vector{UInt8}(vert), GL_VERTEX_SHADER, :vertexshader)
fragshader = compile_shader(Vector{UInt8}(frag), GL_FRAGMENT_SHADER, :fragshader)

program = compile_program(vertshader, fragshader)

scene = (
    Vec4f0(10),
    cam[Projection],
    cam[View],
    cam[ProjectionView],
    Vec2f0(widths(cam[Area]))
)

scene_buff = UniformBuffer(scene)
FieldTraits.on(cam, ProjectionView) do projview
    scene_buff[1] = (
        Vec4f0(10),
        cam[Visualize.Projection],
        cam[View],
        projview,
        Vec2f0(widths(cam[Area]))
    )
end

# Create empty Buffer
uniform_array = UniformBuffer(Tuple{Vec4f0, Mat4f0}, 1024 * 2)

function loadmeshes(folder)
    meshpaths = filter(x-> endswith(x, ".ifs"), readdir(folder))[1:1024]
    faces = GLTriangle[]
    vertices = Tuple{Point3f0, Normal{3, Float32}}[]
    v0 = (Point3f0(typemax(Float32)), Point3f0(typemin(Float32)))
    fidx = 0; vidx = 0;
    drawcommands = Vector{Command}(length(meshpaths))
    for (i, meshpath) in enumerate(meshpaths)
        mesh = read_ifs(joinpath(folder, meshpath))
        fs, vs = mesh.indexes[1], mesh.parent
        append!(faces, fs)
        ns = normals(vs, fs)
        append!(vertices, zip(vs, ns))
        mini, maxi = reduce((x, v)-> (min.(x[1], v), max.(x[2], v)), v0, mesh.parent)
        x, y = ind2sub((32, 32), i)
        trans = translationmatrix(Vec3f0(x, y, 0f0))
        s = maximum(maxi .- mini)
        scale = scalematrix(Vec3f0(1f0 ./ s))
        push!(uniform_array, (
            Vec4f0(rand(Vec3f0)..., 1f0),
            trans * scale * translationmatrix(-Vec3f0(mini))
        ))
        drawcommands[i] = Command(length(fs) * 3, 1, fidx, vidx, 0)
        fidx += length(fs) * 3; vidx += length(vs)
    end
    vbo = VertexArray(view(vertices, faces))
    ibuff = GLBuffer(drawcommands, buffertype = GL_DRAW_INDIRECT_BUFFER)
    vbo, ibuff
end

vbo, commandbuff = loadmeshes(homedir() * "/3dstuff/models")

sceneidx = glGetUniformBlockIndex(program, "Scene")
vertex_arts_idx = glGetUniformBlockIndex(program, "VertexArguments")

glUniformBlockBinding(program, sceneidx, 0)
glUniformBlockBinding(program, vertex_arts_idx, 1)
glBindBufferBase(GL_UNIFORM_BUFFER, 0, scene_buff.buffer.id)
glBindBufferBase(GL_UNIFORM_BUFFER, 1, uniform_array.buffer.id)

function rendloop(window, N, frame_times, commandbuff)
    glUseProgram(program)
    glEnable(GL_DEPTH_TEST)
    glEnable(GL_DEPTH_TEST)
    glDepthMask(GL_TRUE)
    glDepthFunc(GL_LEQUAL)
    glDisable(GL_CULL_FACE)
    glEnable(GL_BLEND)
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
    glClearColor(1, 1, 1, 1)
    GLAbstraction.bind(commandbuff)
    n = 0
    glBindVertexArray(vbo.id)
    while isopen(window) && n <= N
        tic()
        GLWindow.poll_glfw()
        GLAbstraction.@gputime begin
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

            glMultiDrawElementsIndirect(
                GL_TRIANGLES,
                GL_UNSIGNED_INT,
                C_NULL, length(commandbuff), 0
            )
            # glBindVertexArray(0)
            GLWindow.swapbuffers(window)
        end
        push!(frame_times, toq())
        n += 1
    end
    frame_times
end
times = Float64[]
rendloop(window, 60, times, commandbuff)
median(times) * 1000

GLFW.DestroyWindow(window)
window.handle = C_NULL

Profile.clear()
@profile rendloop(window, 2000, timsies)
Profile.print(noisefloor = 2.0)

#timsies = NTuple{3, Float64}[]
timsies = Float64[]
rendloop(window, timsies)
maximum(timsies) * 1000

mean(first.(timsies)) * 1000
mean(map(x->x[3], timsies)) * 1000

open("test.csv", "w") do io
    for t in timsies
        println(io, t)
    end
end
println(pwd())

for i = 1:100
    buzz = i % 5 == 0 ? "Buzz" : ""
    fizz = i % 3 == 0 ? "Fizz" : ""
    println(fizz*buzz)
    if buzz || fizz
        if fizz && buzz
            println("FizzBuzz")
        else
            fizz && println("Buzz")
            buzz && println("Buzz")
        end
    else
        println(i)
    end
end


text:
module Utilities

using DynamicPPL: AbstractSampler, Sampler
using DynamicPPL: init, inittrans, reconstruct, reconstruct!, vectorize
using Distributions, Bijectors
using StatsFuns, SpecialFunctions
import Distributions: sample

export  vectorize,
        reconstruct,
        reconstruct!,
        Sample,
        Chain,
        init,
        set_resume!,
        FlattenIterator

include("helper.jl")

end # module


text:
function content_headers!(headers::AbstractDict, ::Val{:json})
    headers["Accept"] = "application/json"
    headers["Content-Type"] = "application/json"
    return headers
end

function content_headers!(headers::AbstractDict, ::Val{:jpeg})
    headers["Accept"] = "multipart/related; type = \"image/jpeg\""
    headers["Content-Type"] = "application/json"
    return headers
end


text:
using FixedPointNumbers, Test

if VERSION >= v"1.6.0-DEV.816" # JuliaLang/julia #36962 # FIXME
    @test isempty(detect_ambiguities(FixedPointNumbers))
else
    @test isempty(detect_ambiguities(FixedPointNumbers, Base, Core))
end

using Documenter
doctest(FixedPointNumbers, manual = false)

@testset "normed" begin
    include("normed.jl")
end
@testset "fixed" begin
    include("fixed.jl")
end

@testset "traits" begin
    include("traits.jl")
end


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

mutable struct ApplicationGatewayPrivateLinkIpConfiguration <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function ApplicationGatewayPrivateLinkIpConfiguration(;id=nothing)
        o = new()
        validate_property(ApplicationGatewayPrivateLinkIpConfiguration, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type ApplicationGatewayPrivateLinkIpConfiguration

const _property_map_ApplicationGatewayPrivateLinkIpConfiguration = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_ApplicationGatewayPrivateLinkIpConfiguration = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ ApplicationGatewayPrivateLinkIpConfiguration }) = collect(keys(_property_map_ApplicationGatewayPrivateLinkIpConfiguration))
Swagger.property_type(::Type{ ApplicationGatewayPrivateLinkIpConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayPrivateLinkIpConfiguration[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPrivateLinkIpConfiguration }, property_name::Symbol) =  _property_map_ApplicationGatewayPrivateLinkIpConfiguration[property_name]

function check_required(o::ApplicationGatewayPrivateLinkIpConfiguration)
    true
end

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


text:
@testset "DTM: TextHashFunction" begin
    sample_text = "This is sample text"
    sample_doc = StringDocument(sample_text)
    thf = TextHashFunction()
    @test thf isa TextHashFunction
    @test thf.cardinality == StringAnalysis.DEFAULT_CARDINALITY
    @test thf.hash_function == hash
    @test length(hash_dtv(sample_doc, thf)) ==
            StringAnalysis.DEFAULT_CARDINALITY
    c = 15
    thf = TextHashFunction(c)
    @test c == cardinality(thf) == length(hash_dtv(sample_doc, thf))
end

@testset "DTM: Generics" begin
    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    fd = FileDocument(sample_file)
    sd = StringDocument(text(fd))

    docs = [text(fd), text(sd)]
    crps = Corpus([fd, sd])
    for ngram_complexity in [1, 2]
        m = DocumentTermMatrix(crps, ngram_complexity=ngram_complexity)
        @test m isa DocumentTermMatrix
        @test dtm(m) isa SparseMatrixCSC

        update_lexicon!(crps, ngram_complexity)

        m = DocumentTermMatrix(crps, ngram_complexity=ngram_complexity)
        dtm(m)
        @test dtm(DocumentTermMatrix(docs; ngram_complexity=ngram_complexity)) == dtm(m)

        tf_idf(dtm(m))

        doc_idx = 1
        dtv(crps[doc_idx], lexicon(crps), ngram_complexity=ngram_complexity) ==
            dtv(crps, doc_idx, ngram_complexity=ngram_complexity)
        @test_throws ErrorException dtv(crps[1])  # test failure

        @test hash_dtv(text(crps[1]), TextHashFunction(), ngram_complexity=ngram_complexity) ==
            hash_dtv(crps[1], TextHashFunction(), ngram_complexity=ngram_complexity)
        v = hash_dtv(text(crps[1]), cardinality=25, ngram_complexity=ngram_complexity)
        @test v == hash_dtv(crps[1], cardinality=25, ngram_complexity=ngram_complexity)
        @test v isa SparseVector{StringAnalysis.DEFAULT_DTM_TYPE}
        @test length(v) == 25

    end

    dtm1 = dtm(crps)
    dtm1sp = sparse(dtm(crps))
    hash_dtm(crps)

    # Regex dtv
    doc = "a..b"
    lex = OrderedDict("aaa"=>1, "aaab"=>2, "accb"=>3, "bbb"=>4)
    v = dtv_regex(doc, lex, Float32)
    v2 = dtv_regex(NGramDocument(doc), lex, Float32)
    @test v == v2 == Float32[0, 1, 1, 0]

    for ngram_complexity in [1, 2]
        # construct a DocumentTermMatrix from a crps and a custom terms vector
        terms = ["And", "notincrps"]
        m = DocumentTermMatrix(crps, terms, ngram_complexity=ngram_complexity)
        @test size(dtm(m),2) == length(terms)
        @test terms == m.terms
        @test size(dtm(m),1) == length(crps)

        # construct a DocumentTermMatrix from a crps and a custom lexicon
        lex = OrderedDict("And"=>1, "notincrps"=>4)
        m = DocumentTermMatrix(crps, lex, ngram_complexity=ngram_complexity)
        @test size(dtm(m), 2) == length(keys(lex))
        @test size(dtm(m), 2) == length(m.terms)
        @test size(dtm(m), 1) == length(crps)

        # construct a DocumentTermMatrix from a dtm and terms vector
        terms = m.terms
        m2 = DocumentTermMatrix(dtm1, terms)
        @test m.row_indices == m2.row_indices
        m2 = DocumentTermMatrix(dtm1sp, terms)
        @test m.row_indices == m2.row_indices
    end
end

@testset "DTM: Iteration, Indexing" begin
    txt = "This is a Document."
    txt2 = "This is yet another document."
    txt3 = "Its full of documents !"
    txt4 = "This is yet another document dude"
    docs = map(StringDocument, [txt, txt2, txt3, txt4])
    crps = Corpus(docs)
    for ngram_complexity in [1, 2]
        m = DocumentTermMatrix(crps, ngram_complexity=ngram_complexity)
        # Iteration iterface tests
        T = Int8
        for (i,v) in enumerate(each_dtv(crps, eltype=T, ngram_complexity=ngram_complexity))
            @test v == m.dtm[1:end, i]
            i==1 && @test v isa SparseVector{T}
        end
        for (i,v) in enumerate(each_hash_dtv(crps, eltype=T, ngram_complexity=ngram_complexity))
            @test v == hash_dtv(crps[i], ngram_complexity=ngram_complexity)
            i==1 && @test v isa SparseVector{T}
        end
        # Indexing into the DTM
        word = "This"
        @test m[word] == m.dtm[m.row_indices[word], :]
        i = 1; j = 2; ii = 1:2
        @test m[i] == m.dtm[i]
        @test m[i, j] == m.dtm[i, j]
        @test m[ii] == m.dtm[ii]
    end
end