text:
using TransformVariables
# using Distributions
using OmegaCore:StdNormal, StdUniform
# 1. Space might already be constrained
# 2. Custom proposals?
# 3. Dont modify derived values
# 4.

# function propose(rng, ::StdNormal, x)

# end

# function propose(rng, ::StdUniform, x)
#   x_ = transform(asℝ, x) + σ * randn(rng)
#   inverse()

#   inv_transform(transform(x) + σ * randn(rng))
# end

"Compute a score using the change in prior of the *single* changed site"
function proposalkernel(kernel, x)
  ∇logdensity(x) = x |> transform |> jacobian |> abs |> log
  before = ∇logdensity(x)
  proposed = kernel(x)
  after = ∇logdensity(proposed)
  ratio = after - before
  proposed, ratio
end

# Do the transform (move unconstrained space)

normalkernel(rng, x, σ = 0.1) = proposalkernel(x) do x
        inv_transform(transform(x) + σ * randn(rng))
      end

"Changes a uniformly chosen single site with kernel"
function swapsinglesite(transitionkernel, rng, ω)
  logtranstionp = 0.0
  function updater(x)
    result, logtranstionp = transitionkernel(x)
    result
  end
  update(ω, rand(rng, keys(ω)), updater), logtranstionp
end

function moveproposal(rng, ω; σ = 1.0)
  swapsinglesite(rng, ω) do x
    normalkernel(rng, x, σ)
  end
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text: abstract type Drawable end

ext: jl


text:
# Autogenerated wrapper script for Libxc_jll for x86_64-apple-darwin-libgfortran4
export libxc

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Libxc")
JLLWrappers.@declare_library_product(libxc, "@rpath/libxc.9.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libxc,
        "lib/libxc.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Manifolds, Manopt, Test, ManifoldsBase

@testset "Differentials" begin
    p = [1.0, 0.0, 0.0]
    q = [0.0, 1.0, 0.0]
    M = Sphere(2)
    X = log(M, p, q)
    Y = similar(X)
    @testset "Differentials on Sn(2)" begin
        @test differential_log_basepoint(M, p, p, X) == -X
        differential_log_basepoint!(M, Y, p, p, X)
        @test Y == -X
        @test differential_log_basepoint(M, p, q, X) == -X
        differential_log_basepoint!(M, Y, p, q, X) == -X
        @test Y == -X
        @test differential_log_argument(M, p, p, X) == X
        differential_log_argument!(M, Y, p, p, X) == X
        @test Y == X
        @test differential_log_argument(M, p, q, X) == zero_tangent_vector(M, q)
        differential_log_argument!(M, Y, p, q, X)
        @test Y == zero_tangent_vector(M, q)
        @test differential_exp_basepoint(M, p, zero_tangent_vector(M, p), X) == X
        differential_exp_basepoint!(M, Y, p, zero_tangent_vector(M, p), X)
        @test Y == X
        @test norm(M, q, differential_exp_basepoint(M, p, X, X) - [-π / 2, 0.0, 0.0]) ≈ 0 atol =
            6 * 10^(-16)
        differential_exp_basepoint!(M, Y, p, X, X)
        @test norm(M, q, Y - [-π / 2, 0.0, 0.0]) ≈ 0 atol = 6 * 10^(-16)
        @test differential_exp_argument(M, p, zero_tangent_vector(M, p), X) == X
        differential_exp_argument!(M, Y, p, zero_tangent_vector(M, p), X) == X
        @test Y == X
        @test norm(M, q, differential_exp_argument(M, p, X, zero_tangent_vector(M, p))) ≈ 0
        differential_exp_argument!(M, Y, p, X, zero_tangent_vector(M, p))
        @test norm(M, q, Y) ≈ 0
        for t in [0, 0.15, 0.33, 0.66, 0.9]
            @test differential_geodesic_startpoint(M, p, p, t, X) == (1 - t) * X
            differential_geodesic_startpoint!(M, Y, p, p, t, X)
            @test Y == (1 - t) * X
            @test norm(M, p, differential_geodesic_endpoint(M, p, p, t, X) - t * X) ≈ 0 atol =
                10.0^(-16)
            differential_geodesic_endpoint!(M, Y, p, p, t, X)
            @test norm(M, p, Y - t * X) ≈ 0 atol = 10.0^(-16)
        end
    end
    @testset "Differentials on Power of Sn(2)" begin
        N = PowerManifold(M, NestedPowerRepresentation(), 3)
        x = [p, q, p]
        y = [p, p, q]
        V = [X, zero_tangent_vector(M, p), -X]
        W = similar.(V)
        @test norm(
            N,
            x,
            differential_forward_logs(N, x, V) -
            [-X, [π / 2, 0.0, 0.0], zero_tangent_vector(M, p)],
        ) ≈ 0 atol = 8 * 10.0^(-16)
        differential_forward_logs!(N, W, x, V)
        @test norm(N, x, W - [-X, [π / 2, 0.0, 0.0], zero_tangent_vector(M, p)]) ≈ 0 atol =
            8 * 10.0^(-16)
        @test differential_log_argument(N, x, y, V) == [V[1], V[2], V[2]]
        differential_log_argument!(N, W, x, y, V)
        @test W == [V[1], V[2], V[2]]
    end
    @testset "Differentials on SPD(2)" begin
        #
        # Single differentials on Hn
        M2 = SymmetricPositiveDefinite(2)
        p2 = [1.0 0.0; 0.0 1.0]
        X2 = [0.5 1.0; 1.0 0.5]
        q2 = exp(M2, p2, X2)
        # Text differentials (1) Dx of Log_xy
        @test norm(M2, p2, differential_log_basepoint(M2, p2, p2, X2) + X2) ≈ 0 atol =
            4 * 10^(-16)
        @test norm(
            M2, q2, differential_log_argument(M2, p2, q2, zero_tangent_vector(M2, p2))
        ) ≈ 0 atol = 4 * 10^(-16)
        @test norm(
            M2, p2, differential_exp_basepoint(M2, p2, zero_tangent_vector(M2, p2), X2) - X2
        ) ≈ 0 atol = 4 * 10^(-16)
        @test norm(
            M2, p2, differential_exp_argument(M2, p2, zero_tangent_vector(M2, p2), X2) - X2
        ) ≈ 0 atol = 4 * 10^(-16)
        for t in [0, 0.15, 0.33, 0.66, 0.9]
            @test norm(
                M2, p2, differential_geodesic_startpoint(M2, p2, p2, t, X2) - (1 - t) * X2
            ) ≈ 0 atol = 4 * 10^(-16)
            @test norm(M2, p2, differential_geodesic_endpoint(M2, p2, p2, t, X2) - t * X2) ≈
                  0 atol = 4 * 10.0^(-16)
        end
        @test norm(M2, q2, differential_geodesic_startpoint(M2, p2, q2, 1.0, X2)) ≈ 0 atol =
            4 * 10.0^(-16)
        @test norm(
            M2, q2, differential_exp_basepoint(M2, p2, X2, zero_tangent_vector(M2, p2))
        ) ≈ 0 atol = 4 * 10.0^(-16)
        @test norm(
            M2, q2, differential_exp_argument(M2, p2, X2, zero_tangent_vector(M2, p2))
        ) ≈ 0 atol = 4 * 10.0^(-16)
        # test coeff of log_basepoint, since it is not always expicitly used.
        @test βdifferential_log_basepoint(-1.0, 1.0, 2.0) ≈ -2 * cosh(2.0) / sinh(2.0)
    end
    @testset "Differentials on Euclidean(2)" begin
        M3 = Euclidean(2)
        x3 = [1.0, 2.0]
        ξ3 = [1.0, 0.0]
        @test norm(M3, x3, differential_exp_basepoint(M3, x3, ξ3, ξ3) - ξ3) ≈ 0 atol =
            4 * 10.0^(-16)
    end
    @testset "Differentials on the Circle" begin
        M = Circle()
        p = 0
        q = π / 4
        X = π / 8
        @test differential_log_argument(M, p, q, X) == X
    end
    @testset "forward logs on a multivariate power manifold" begin
        S = Sphere(2)
        M = PowerManifold(S, NestedPowerRepresentation(), 2, 2)
        p = [zeros(3) for i in [1, 2], j in [1, 2]]
        p[1, 1] = [1.0, 0.0, 0.0]
        p[1, 2] = 1 / sqrt(2) .* [1.0, 1.0, 0.0]
        p[2, 1] = 1 / sqrt(2) .* [1.0, 0.0, 1.0]
        p[2, 2] = [0.0, 1.0, 0.0]
        t1 = forward_logs(M, p)
        @test t1[1, 1, 1] ≈ log(S, p[1, 1], p[2, 1])
        @test t1[1, 1, 2] ≈ log(S, p[1, 1], p[1, 2])
        @test t1[1, 2, 1] ≈ log(S, p[1, 2], p[2, 2])
        @test t1[1, 2, 2] ≈ log(S, p[1, 2], p[1, 2]) atol = 1e-15
        @test t1[2, 1, 1] ≈ log(S, p[2, 1], p[2, 1]) atol = 1e-15
        @test t1[2, 1, 2] ≈ log(S, p[2, 1], p[2, 2])
        @test t1[2, 2, 1] ≈ log(S, p[2, 2], p[2, 2])
        @test t1[2, 2, 2] ≈ log(S, p[2, 2], p[2, 2])
        t1a = zero.(t1)
        forward_logs!(M, t1a, p)
        @test all(t1 .== t1a)
        X = zero_tangent_vector(M, p)
        X[1, 1] .= [0.0, 0.5, 0.5]
        t2 = differential_forward_logs(M, p, X)
        a =
            differential_log_basepoint(S, p[1, 1], p[2, 1], X[1, 1]) +
            differential_log_argument(S, p[1, 1], p[2, 1], X[2, 1])
        @test t2[1, 1, 1] ≈ a
        @test t2[1, 2, 1] ≈ zero_tangent_vector(S, p[1, 2]) atol = 1e-17
        @test t2[2, 1, 1] ≈ zero_tangent_vector(S, p[2, 1]) atol = 1e-17
        @test t2[2, 2, 1] ≈ zero_tangent_vector(S, p[2, 2]) atol = 1e-17
        b =
            differential_log_basepoint(S, p[1, 1], p[1, 2], X[1, 1]) +
            differential_log_argument(S, p[1, 1], p[1, 2], X[1, 2])
        @test t2[1, 1, 2] ≈ b
        @test t2[1, 2, 2] ≈ zero_tangent_vector(S, p[1, 2]) atol = 1e-17
        @test t2[2, 1, 2] ≈ zero_tangent_vector(S, p[2, 1]) atol = 1e-17
        @test t2[2, 2, 2] ≈ zero_tangent_vector(S, p[2, 2]) atol = 1e-17
    end
end


text:
module PProf

export pprof, @pprof

using Profile
using ProtoBuf
using OrderedCollections
import pprof_jll

using Profile: clear

"""
    PProf.clear()

Alias for `Profile.clear()`
"""
clear

include(joinpath("..", "lib", "perftools.jl"))

import .perftools.profiles: ValueType, Sample, Function,
                            Location, Line
const PProfile = perftools.profiles.Profile

const proc = Ref{Union{Base.Process, Nothing}}(nothing)

"""
    _enter!(dict::OrderedDict{T, Int64}, key::T) where T

Resolves from `key` to the index (zero-based) in the dict.
Useful for the Strings table

NOTE: We must use Int64 throughout this package (regardless of system word-size) b/c the
proto file specifies 64-bit integers.
"""
function _enter!(dict::OrderedDict{T, Int64}, key::T) where T
    if haskey(dict, key)
        return dict[key]
    else
        l = Int64(length(dict))
        dict[key] = l
        return l
    end
end

using Base.StackTraces: StackFrame

# TODO:
# - Mappings

"""
    pprof([data, [lidict]];
            web = true, webhost = "localhost", webport = 57599,
            out = "profile.pb.gz", from_c = true, full_signatures = true, drop_frames = "",
            keep_frames = "", ui_relative_percentages = true, sampling_delay = nothing,
         )
    pprof(FlameGraphs.flamegraph(); kwargs...)

Fetches the collected `Profile` data, exports to the `pprof` format, and (optionally) opens
a `pprof` web-server for interactively viewing the results.

If `web=true`, the web-server is opened in the background. Re-running `pprof()` will refresh
the web-server to use the new output.

If you manually edit the output file, `PProf.refresh()` will refresh the server without
overwriting the output file. `PProf.kill()` will kill the server.

You can also use `PProf.refresh(file="...")` to open a new file in the server.

# Arguments:
- `data::Vector{UInt}`: The data provided by `Profile.retrieve()` [optional].
- `lidict::Dict`: The lookup dictionary provided by `Profile.retrieve()` [optional].
    - Note that you need both the `data` and the `lidict` returned from
      `Profile.retrieve()` in order to export profiling data between julia processes.
- `flamegraph`: PProf also accepts profile data passed as a `FlameGraphs.jl` graph object.

# Keyword Arguments
- `sampling_delay::UInt64`: The period between samples in nanoseconds [optional].
- `web::Bool`: Whether to launch the `go tool pprof` interactive webserver for viewing results.
- `webhost::AbstractString`: If using `web`, which host to launch the webserver on.
- `webport::Integer`: If using `web`, which port to launch the webserver on.
- `out::String`: Filename for output.
- `from_c::Bool`: If `false`, exclude frames that come from from_c. Defaults to `true`.
- `full_signatures::Bool`: If `true`, methods are printed as signatures with full argument
                           types. If `false`, as only names. E.g. `eval(::Module, ::Any)` vs `eval`.
- `drop_frames`: frames with function_name fully matching regexp string will be dropped from the samples,
                 along with their successors.
- `keep_frames`: frames with function_name fully matching regexp string will be kept, even if it matches drop_functions.
- `ui_relative_percentages`: Passes `-relative_percentages` to pprof. Causes nodes
  ignored/hidden through the web UI to be ignored from totals when computing percentages.
"""
function pprof(data::Union{Nothing, Vector{UInt}} = nothing,
               lidict::Union{Nothing, Dict} = nothing;
               sampling_delay::Union{Nothing, UInt64} = nothing,
               web::Bool = true,
               webhost::AbstractString = "localhost",
               webport::Integer = 57599,
               out::AbstractString = "profile.pb.gz",
               from_c::Bool = true,
               full_signatures::Bool = true,
               drop_frames::Union{Nothing, AbstractString} = nothing,
               keep_frames::Union{Nothing, AbstractString} = nothing,
               ui_relative_percentages::Bool = true,
            )
    if data === nothing
        data = copy(Profile.fetch())
    end
    lookup = lidict
    if lookup === nothing
        lookup = Profile.getdict(data)
    end
    if sampling_delay === nothing
        sampling_delay = ccall(:jl_profile_delay_nsec, UInt64, ())
    end
    @assert !isempty(basename(out)) "`out=` must specify a file path to write to. Got unexpected: '$out'"
    if !endswith(out, ".pb.gz")
        out = "$out.pb.gz"
        @info "Writing output to $out"
    end

    string_table = OrderedDict{AbstractString, Int64}()
    enter!(string) = _enter!(string_table, string)
    enter!(::Nothing) = _enter!(string_table, "nothing")
    ValueType!(_type, unit) = ValueType(_type = enter!(_type), unit = enter!(unit))

    # Setup:
    enter!("")  # NOTE: pprof requires first entry to be ""
    # Functions need a uid, we'll use the pointer for the method instance
    seen_funcs = Set{UInt64}()
    funcs = Dict{UInt64, Function}()

    seen_locs = Set{UInt64}()
    locs  = Dict{UInt64, Location}()
    locs_from_c  = Dict{UInt64, Bool}()

    sample_type = [
        ValueType!("events",      "count"), # Mandatory
        ValueType!("stack_depth", "count")
    ]

    prof = PProfile(
        sample = [], location = [], _function = [],
        mapping = [], string_table = [],
        sample_type = sample_type, default_sample_type = 1, # events
        period = sampling_delay, period_type = ValueType!("cpu", "nanoseconds")
    )

    if drop_frames !== nothing
        prof.drop_frames = enter!(drop_frames)
    end
    if keep_frames !== nothing
        prof.keep_frames = enter!(keep_frames)
    end

    # start decoding backtraces
    location_id = Vector{eltype(data)}()
    lastwaszero = true

    for ip in data
        # ip == 0x0 is the sentinel value for finishing a backtrace, therefore finising a sample
        if ip == 0
            # Avoid creating empty samples
            if lastwaszero
                @assert length(location_id) == 0
                continue
            end

            # End of sample
            value = [
                1,                   # events
                length(location_id), # stack_depth
            ]
            push!(prof.sample, Sample(;location_id = location_id, value = value))
            location_id = Vector{eltype(data)}()
            lastwaszero = true
            continue
        end
        lastwaszero = false

        # A backtrace consists of a set of IP (Instruction Pointers), each IP points
        # a single line of code and `litrace` has the necessary information to decode
        # that IP to a specific frame (or set of frames, if inlining occured).

        # if we have already seen this IP avoid decoding it again
        if ip in seen_locs
            # Only keep C frames if from_c=true
            if (from_c || !locs_from_c[ip])
                push!(location_id, ip)
            end
            continue
        end
        push!(seen_locs, ip)

        # Decode the IP into information about this stack frame (or frames given inlining)
        location = Location(;id = ip, address = ip, line=[])
        location_from_c = true
        # Will have multiple frames if frames were inlined (the last frame is the "real
        # function", the inlinee)
        for frame in lookup[ip]
            # ip 0 is reserved
            frame.pointer == 0 && continue

            # if any of the frames is not from_c the entire location is not from_c
            location_from_c &= frame.from_c

            # `func_id` - Uniquely identifies this function (a method instance in julia, and
            # a function in C/C++).
            # Note that this should be unique even for several different functions all
            # inlined into the same frame.
            func_id = if frame.linfo !== nothing
                hash(frame.linfo)
            else
                hash((frame.func, frame.file, frame.line))
            end
            push!(location.line, Line(function_id = func_id, line = frame.line))

            # Known function
            func_id in seen_funcs && continue
            push!(seen_funcs, func_id)

            # Store the function in our functions dict
            funcProto = Function()
            funcProto.id = func_id
            file = nothing
            simple_name = _escape_name_for_pprof(frame.func)
            local full_name_with_args
            if frame.linfo !== nothing && frame.linfo isa Core.MethodInstance
                linfo = frame.linfo::Core.MethodInstance
                meth = linfo.def
                file = string(meth.file)
                io = IOBuffer()
                Base.show_tuple_as_call(io, meth.name, linfo.specTypes)
                name = String(take!(io))
                full_name_with_args = _escape_name_for_pprof(name)
                funcProto.start_line = convert(Int64, meth.line)
            else
                # frame.linfo either nothing or CodeInfo, either way fallback
                file = string(frame.file)
                full_name_with_args = _escape_name_for_pprof(string(frame.func))
                funcProto.start_line = convert(Int64, frame.line) # TODO: Get start_line properly
            end
            # WEIRD TRICK: By entering a separate copy of the string (with a
            # different string id) for the name and system_name, pprof will use
            # the supplied `name` *verbatim*, without pruning off the arguments.
            # So even when full_signatures == false, we want to generate two `enter!` ids.
            funcProto.system_name = enter!(simple_name)
            if full_signatures
                funcProto.name = enter!(full_name_with_args)
            else
                funcProto.name = enter!(simple_name)
            end
            file = Base.find_source_file(file)
            funcProto.filename   = enter!(file)
            # Only keep C functions if from_c=true
            if (from_c || !frame.from_c)
                funcs[func_id] = funcProto
            end
        end
        locs_from_c[ip] = location_from_c
        # Only keep C frames if from_c=true
        if (from_c || !location_from_c)
            locs[ip] = location
            push!(location_id, ip)
        end
    end

    # Build Profile
    prof.string_table = collect(keys(string_table))
    # If from_c=false funcs and locs should NOT contain C functions
    prof._function = collect(values(funcs))
    prof.location  = collect(values(locs))

    # Write to disk
    open(out, "w") do io
        writeproto(io, prof)
    end

    if web
        refresh(webhost = webhost, webport = webport, file = out,
            ui_relative_percentages = ui_relative_percentages)
    end

    out
end

function _escape_name_for_pprof(name)
    # HACK: Apparently proto doesn't escape func names with `"` in them ... >.<
    # TODO: Remove this hack after https://github.com/google/pprof/pull/564
    quoted = repr(string(name))
    quoted = quoted[2:thisind(quoted, end-1)]
    return quoted
end

"""
    refresh(; webhost = "localhost", webport = 57599, file = "profile.pb.gz",
            ui_relative_percentages = true)

Start or restart the go pprof webserver.

- `webhost::AbstractString`: Which host to launch the webserver on.
- `webport::Integer`: Which port to launch the webserver on.
- `file::String`: Profile file to open.
- `ui_relative_percentages::Bool`: Passes `-relative_percentages` to pprof. Causes nodes
  ignored/hidden through the web UI to be ignored from totals when computing percentages.
"""
function refresh(; webhost::AbstractString = "localhost",
                   webport::Integer = 57599,
                   file::AbstractString = "profile.pb.gz",
                   ui_relative_percentages::Bool = true,
                )

    if proc[] === nothing
        # The first time, register an atexit hook to kill the web server.
        atexit(PProf.kill)
    else
        # On subsequent calls, restart the pprof web server.
        Base.kill(proc[])
    end

    relative_percentages_flag = ui_relative_percentages ? "-relative_percentages" : ""

    proc[] = pprof_jll.pprof() do pprof_path
        open(pipeline(`$pprof_path -http=$webhost:$webport $relative_percentages_flag $file`))
    end
end

"""
    pprof_kill()

Kills the pprof server if running.
"""
function kill()
    if proc[] !== nothing
        Base.kill(proc[])
        proc[] = nothing
    end
end

"""
    @pprof ex

Profiles the expression using `@profile` and starts or restarts the `pprof()` web UI with
default arguments.
"""
macro pprof(ex)
    esc(quote
        $Profile.@profile $ex
        $(@__MODULE__).pprof()
    end)
end

include("flamegraphs.jl")

if isdefined(Profile, :Allocs)  # PR https://github.com/JuliaLang/julia/pull/42768
    include("Allocs.jl")
end

# Precompile as much as possible, so that profiling doesn't end up measuring our own
# compilation.
function __init__()
    precompile(pprof, ()) || error("precompilation of package functions is not supposed to fail")
    precompile(kill, ()) || error("precompilation of package functions is not supposed to fail")
    precompile(refresh, ()) || error("precompilation of package functions is not supposed to fail")
end
end # module


text:
function gradient_replicates!(gstore, dfns, parms, pts, datav)
  fill!(gstore, zero(eltype(gstore)))
  @assert length(gstore) == length(dfns) "Length of gstore and number of derivative functions don't match."
  _S = Symmetric([matern(x, y, parms) for x in pts, y in pts])
  S  = cholesky(_S)
  for j in eachindex(gstore)
    dfj = dfns[j]
    Sj  = Symmetric([dfj(x, y, parms) for x in pts, y in pts])
    gstore[j] = (tr(S\Sj)*length(datav) - sum(z->dot(z, S\(Sj*(S\z))), datav))/2
  end
  gstore
end


text:
# p36.m - Laplace eq. on [-1,1]x[-1,1] with nonzero BCs

# Set up grid and 2D Laplacian, boundary points included:
N = 24; (D,x) = cheb(N); y = x;
xx = repmat(x',N+1,1); yy = repmat(y,1,N+1);
D2 = D^2; I = eye(N+1); L = kron(I,D2) + kron(D2,I);

# Impose boundary conditions by replacing appropriate rows of L:
b = @. (abs(xx[:])==1) | (abs(yy[:])==1);            # boundary pts
L[b,:] = 0; L[b,b] = eye(4*N);
rhs = zeros((N+1)^2);
rhs[b] = @. (yy[b]==1)*(xx[b]<0)*sin(pi*xx[b])^4 + .2*(xx[b]==1)*sin(3*pi*yy[b]);

# Solve Laplace equation, reshape to 2D, and plot:
u = L\rhs; uu = reshape(u,N+1,N+1);
xxx = yyy = -1:.04:1;
s = interpolate((x[end:-1:1],y[end:-1:1]),reduce(flipdim,uu,1:2),Gridded(Linear()));
clf();
surf(xxx,yyy,s[xxx,yyy]);
xlim(-1,1); ylim(-1,1); zlim(-.2,1);
gca()[:view_init](45,-110);
umid = uu[Int(N/2)+1,Int(N/2)+1];
text3D(0,.8,.4,"u(0,0) = $(signif(umid,9))");


text:
using Plots
include_ddp()

# Model
include_model("rocket3D")

function fd(model::Rocket3D{Midpoint, FixedTime}, x, u, w, h, t)
	return view(x, 1:model.n) + h * f(model, view(x, 1:model.n) + 0.5 * h * f(model, view(x, 1:model.n), view(u, 1:model.m), w), view(u, 1:model.m), w)
end

n = model.n
m = model.m

# Time
T = 201
h = 0.01

# Initial conditions, controls, disturbances
x1 = zeros(model.n)
x1[1] = 2.5
x1[2] = 0.0
x1[3] = 10.0
mrp = MRP(RotY(-0.5 * π) * RotX(0.0 * π))
x1[4:6] = [mrp.x; mrp.y; mrp.z]
x1[9] = -5.0

# visualize!(vis, model, [x1], Δt = h)

# visualize!(vis, model, [x1], Δt = h)

xT = zeros(model.n)
# xT[1] = 2.5
# xT[2] = 0.0
xT[3] = model.length

u_ref = [0.0; 0.0; 0.0]#model.mass * 9.81]
ū = [u_ref + [1.0e-2; 1.0e-2; 1.0e-2] .* randn(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)
# x̄ = linear_interpolation(x1, xT, T)
# plot(hcat(x̄...)')

# Objective
Q = h * [(t < T ? 1.0 * Diagonal([1.0e-1 * ones(3); 0.0 * ones(3); 1.0e-1 * ones(3); 1000.0 * ones(3)])
        : 0.0 * Diagonal(0.0 * ones(model.n))) for t = 1:T]
q = h * [-2.0 * Q[t] * xT for t = 1:T]

R = h * [Diagonal([10000.0; 10000.0; 100.0]) for t = 1:T-1]
r = h * [-2.0 * R[t] * u_ref  for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => [-5.0; -5.0; 0.0], :uu => [5.0; 5.0; 15.0], :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		# c .= [ul - u; u - uu]
	elseif t == T
		xT = cons.con[T].info[:xT]
		c .= x - xT
	else
		nothing
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	con_tol = 1.0e-3,
	ρ_init = 1.0, ρ_scale = 10.0)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

x̄_nominal = x̄
ū_nominal = ū

@save "/home/taylor/Research/motion_planning/examples/differential_dynamic_programming/implicit_dynamics/rocket.jld2" x̄_nominal ū_nominal
@load "/home/taylor/Research/motion_planning/examples/differential_dynamic_programming/implicit_dynamics/rocket.jld2"

# Trajectories
plot(hcat(ū...)', linetype = :steppost)
plot(hcat(x̄...)[1:3, :]', linetype = :steppost)

# t = 1
# idx = [3; 1; 2]
# second_order_cone_projection(ū[t][idx])
# ū[t]
# # Visualize
# include(joinpath(pwd(), "models/visualize.jl"))
# vis = Visualizer()
# render(vis)
# # open(vis)
# visualize!(vis, model, x̄, Δt = h)
#
# # ū_fixed_time = ū
#
# # Visualize
obj_rocket = joinpath(pwd(), "models/starship/Starship.obj")
mtl_rocket = joinpath(pwd(), "models/starship/Starship.mtl")
ctm = ModifiedMeshFileObject(obj_rocket, mtl_rocket, scale=1.0)
# setobject!(vis["rocket"]["starship"], ctm)
#
# settransform!(vis["rocket"]["starship"],
# 	compose(Translation(0.0, 0.0, -model.length),
# 		LinearMap(0.25 * RotY(0.0) * RotZ(0.5 * π) * RotX(0.5 * π))))
#
# default_background!(vis)
# settransform!(vis["rocket"],
# 	compose(Translation(0.0, 0.0, 0.0),
# 	LinearMap(RotY(0.0))))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Test

@testset "PREM" begin
include("PREM_tests.jl")
end

@testset "Oscillation" begin
include("oscillation_tests.jl")
end

@testset "Path" begin
include("path_tests.jl")
end

@testset "Matter" begin
include("matter_tests.jl")
end


text:
function testGrid()
    grid = Grid2d(7,7)
    @assert frequencies(1, 1, grid) == ( 0,0)
    @assert frequencies(2, 3, grid) == ( 1,2)
    @assert frequencies(5, 4, grid) == (-3,3)
end

testGrid()


text:
# This file is a part of UnsafeArrays.jl, licensed under the MIT License (MIT).

import Test

Test.@testset "Package UnsafeArrays" begin
    include("util.jl")
    include("unsafe_array.jl")
    include("uview.jl")
end


text:
using Revise
using Dice
using Dice: num_flips, num_nodes, ifelse

# Number of nodes BWH needs to model a distribution on n bits
bwh_num_nodes(n) = 2^(n + 1) - n

function generate_code_bwh(p::Vector{Float64})
    @dice begin
        @assert ispow2(length(p)) "Distribution length should be power of 2"
        prefix_sums = accumulate(+, p)
        @assert last(prefix_sums) ≈ 1 "Distribution should sum to 1"
        function helper(i, j)
            if i == j
                DistInt(dicecontext(), i - 1)
            else
                first_half_end = (i + j - 1) ÷ 2
                first_half_p = prefix_sums[first_half_end] - prefix_sums[i] + p[i]
                total_p = prefix_sums[j] - prefix_sums[i] + p[i]
                ifelse(flip(first_half_p/total_p),
                    helper(i, first_half_end),
                    helper(first_half_end+1, j))
            end
        end
        helper(1, length(p))
    end
end

function test_bwh_num_nodes(p::Vector{Float64})
    code = generate_code_bwh(p)
    bdd = compile(code)
    @assert infer(code, :bdd) ≈ p  # Verify correctness

    bdd_nodes = num_nodes(bdd)
    num_bits = round(Int, log2(length(p)))
    @assert bdd_nodes == bwh_num_nodes(num_bits)  # Verify num bits

    println("BWH for a $(num_bits)-bit distribution: $(bdd_nodes) nodes, $(num_flips(bdd)) flips")
end

dist_2 = [0.5, 0.5]
dist_4 = [0.1, 0.2, 0.3, 0.4]
dist_8 = [0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.44]
dist_16 = [0.005, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.001, 0.11, 0.12, 0.13, 0.14, 0.022, 0.032]

for dist in [dist_2, dist_4, dist_8, dist_16]
    test_bwh_num_nodes(dist)
end


text:
function eigsolve(A, x₀, howmany::Int, which::Selector, alg::Lanczos)
    krylovdim = alg.krylovdim
    maxiter = alg.maxiter
    howmany > krylovdim &&
        error("krylov dimension $(krylovdim) too small to compute $howmany eigenvalues")

    ## FIRST ITERATION: setting up

    # Initialize Lanczos factorization
    iter = LanczosIterator(A, x₀, alg.orth)
    fact = initialize(iter; verbosity = alg.verbosity - 2)
    numops = 1
    numiter = 1
    sizehint!(fact, krylovdim)
    β = normres(fact)
    tol::typeof(β) = alg.tol

    # allocate storage
    HH = fill(zero(eltype(fact)), krylovdim + 1, krylovdim)
    UU = fill(zero(eltype(fact)), krylovdim, krylovdim)

    converged = 0
    local D, U, f
    while true
        β = normres(fact)
        K = length(fact)

        # diagonalize Krylov factorization
        if β <= tol
            if K < howmany
                @warn "Invariant subspace of dimension $K (up to requested tolerance `tol = $tol`), which is smaller than the number of requested eigenvalues (i.e. `howmany == $howmany`); setting `howmany = $K`."
                howmany = K
            end
        end
        if K == krylovdim || β <= tol || (alg.eager && K >= howmany)
            U = copyto!(view(UU, 1:K, 1:K), I)
            f = view(HH, K + 1, 1:K)
            T = rayleighquotient(fact) # symtridiagonal

            # compute eigenvalues
            if K == 1
                D = [T[1, 1]]
                f[1] = β
                converged = Int(β <= tol)
            else
                if K < krylovdim
                    T = deepcopy(T)
                end
                D, U = tridiageigh!(T, U)
                by, rev = eigsort(which)
                p = sortperm(D, by = by, rev = rev)
                D, U = permuteeig!(D, U, p)
                mul!(f, view(U, K, :), β)
                converged = 0
                while converged < K && abs(f[converged+1]) <= tol
                    converged += 1
                end
            end

            if converged >= howmany
                break
            elseif alg.verbosity > 1
                msg = "Lanczos eigsolve in iter $numiter, krylovdim = $K: "
                msg *= "$converged values converged, normres = ("
                msg *= @sprintf("%.2e", abs(f[1]))
                for i in 2:howmany
                    msg *= ", "
                    msg *= @sprintf("%.2e", abs(f[i]))
                end
                msg *= ")"
                @info msg
            end
        end

        if K < krylovdim# expand Krylov factorization
            fact = expand!(iter, fact; verbosity = alg.verbosity - 2)
            numops += 1
        else ## shrink and restart
            if numiter == maxiter
                break
            end

            # Determine how many to keep
            keep = div(3 * krylovdim + 2 * converged, 5) # strictly smaller than krylovdim since converged < howmany <= krylovdim, at least equal to converged

            # Restore Lanczos form in the first keep columns
            H = fill!(view(HH, 1:keep+1, 1:keep), zero(eltype(HH)))
            @inbounds for j in 1:keep
                H[j, j] = D[j]
                H[keep+1, j] = f[j]
            end
            @inbounds for j in keep:-1:1
                h, ν = householder(H, j + 1, 1:j, j)
                H[j+1, j] = ν
                H[j+1, 1:j-1] .= zero(eltype(H))
                lmul!(h, H)
                rmul!(view(H, 1:j, :), h')
                rmul!(U, h')
            end
            @inbounds for j in 1:keep
                fact.αs[j] = H[j, j]
                fact.βs[j] = H[j+1, j]
            end

            # Update B by applying U using Householder reflections
            B = basis(fact)
            basistransform!(B, view(U, :, 1:keep))
            r = residual(fact)
            B[keep+1] = rmul!(r, 1 / β)

            # Shrink Lanczos factorization
            fact = shrink!(fact, keep)
            numiter += 1
        end
    end

    if converged > howmany
        howmany = converged
    end
    values = D[1:howmany]

    # Compute eigenvectors
    V = view(U, :, 1:howmany)

    # Compute convergence information
    vectors = let B = basis(fact)
        [B * v for v in cols(V)]
    end
    residuals = let r = residual(fact)
        [last(v) * r for v in cols(V)]
    end
    normresiduals = let f = f
        map(i -> abs(f[i]), 1:howmany)
    end

    if alg.verbosity > 0
        if converged < howmany
            @warn """Lanczos eigsolve finished without convergence after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        else
            @info """Lanczos eigsolve finished after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        end
    end

    return values,
    vectors,
    ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
end


text:
"""
    Controlled(C, P)

"Controlled" iterator `P` with the state `t` of `C`
calling
```
(t => x) = evolve(P, t=>x, tᵒ)
```

# Example
```
# Apply collatz twice each step using the default for Evolve

collatz(n) = n % 2 == 0 ? n÷2 : 3n + 1
collectfrom(Controlled(1:2:20, Evolve(collatz)), (1,14))
```
"""
struct Controlled{T,S} <: DynamicIterator
        C::S
        P::T
end

function dyniterate(M::Controlled, start::Union{Start,Nothing})
    tᵒ, c = @returnnothing dyniterate(M.C, nothing)
    u, p = @returnnothing dyniterate(M.P, Control(start),  tᵒ)
    u, (c, p)
end
function dyniterate(M::Controlled, (c, p)::Tuple)
    tᵒ, c = @returnnothing iterate(M.C, c)
    u, p = @returnnothing dyniterate(M.P, p, tᵒ)
    u, (c, p)
end

struct Attach{F,T} <: DynamicIterator
        f::F
        P::T
end
attach(f, P) = Attach(f, P)

function dyniterate(M::Attach, start::Union{Start,Nothing}, args...)
    dyniterate(M.P, M.f(start),  args...)
end
function dyniterate(M::Attach, state, args...)
    dyniterate(M.P, state, args...)
end


text:
module IndividualRewrite
export Rule, sim
using Catlab.CategoricalAlgebra, Catlab.Present, Catlab.Theories
using Random: shuffle

# Data structures
#################

mutable struct Rule
  L::ACSetTransformation
  R::ACSetTransformation
  delay::Int
  frequency::Int
end

mutable struct RuleCount
  r::Rule
  count::Int
end

RuleCount(r::Rule) = RuleCount(r, 0)

mutable struct MatchCount
  name::Symbol
  rule::Rule
  match::ACSetTransformation
  count::Int
end

MatchCount(name::Symbol, r::Rule, m::ACSetTransformation) =
  MatchCount(name, r, m, r.delay)

# Simulation
############
"""
Run simulation for `n` timesteps.
"""
function sim(init_state::StructACSet, rules::Dict{Symbol,Rule}, n::Int;
             verbose=false)
  rule_counts = Dict(k=>RuleCount(v) for (k,v) in collect(rules))
  matches, states = [], [init_state]
  for step_i in 1:n
    if verbose
      println("\nStep $step_i ($(length(matches)) rewrites queued)")
    end
    # Look for matches if frequency count has reached zero
    for (r_name,r) in collect(rule_counts)
      if r.count > 0
        r.count -= 1
      else
        r.count = r.r.frequency
        rmatches = shuffle(homomorphisms(codom(r.r.L), states[end]))
        if !isempty(rmatches)
          if verbose println("\tAdding match for rewrite $r_name") end
          push!(matches, MatchCount(r_name, r.r, first(rmatches)))
        end
      end
    end
    # Execute matches if their delay count has reached zero
    for (i,m) in enumerate(matches)
      if isnothing(m)
        continue
      elseif m.count > 0
        m.count -= 1
      else
        matches[i] = nothing
        _, kg,_,kh = rewrite_match_maps(m.rule.L, m.rule.R, m.match)
        push!(states, codom(kh))

        if verbose
          println("\tnew state after firing $(m.name)")
          print("\t"); show(stdout, "text/plain", states[end])
        end
        for (j, m_) in enumerate(matches )
          matches[j] = postcompose_partial(kg, kh, m_)
        end
      end
    end
    matches = filter(x->!isnothing(x), matches) # remove executed/invalid match
  end
  return states
end

"""
Convert a match L->G to a match L->H using a partial morphism G->H, if possible.
       L ========= L
     m ↓           ↓ m'
       G <-- K --> H
"""
function postcompose_partial(
    kg::ACSetTransformation,kh::ACSetTransformation, mc::MatchCount)
  println("")
  m = mc.match
  d = Dict()
  for (k,vs) in pairs(components(m))
    vs_ = Int[]
    for v in collect(vs)
      kv = findfirst(==(v), collect(kg[k]))
      if isnothing(kv)
        mc = nothing
        return nothing
      else
        push!(vs_, kh[k](kv))
      end
    end
    d[k] = vs_
  end
  MatchCount(mc.name,mc.rule,ACSetTransformation(dom(m), codom(kh); d...), mc.count)
end

postcompose_partial(::ACSetTransformation,::ACSetTransformation,::Nothing) =
  nothing

end # module


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

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

Performs an Analysis of Variance (ANOVA) calculation.

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

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

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

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

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

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

    firstlevelreplicates = eltype(observations) <: Number ? hasreplicates : false
    nfactors = ndims(observations) - (firstlevelreplicates ? 1 : 0)

    # defaults to assuming all unspecified factors are fixed
    if length(factortypes) < nfactors
        nremaining = nfactors - length(factortypes)
        append!(factortypes, repeat([fixed], nremaining))
    end

    validate(factortypes, factornames, nfactors)

    # automatically assigns alphabetical names if not provided.
    if isempty(factornames)
        factornames = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"][1:nfactors]
        reverse!(factornames)
    end

    nnestedfactors = count(f -> f == nested, factortypes)
    crossedfactortypes = filter(f -> f ∈ [fixed, random], factortypes)
    reverse!(crossedfactortypes)
    ncrossedfactors = length(crossedfactortypes)
    ncrossedfactors < 4 || error("ANOVA with 4 or more crossed factors is not supported.")
    nreplicates = firstlevelreplicates ? size(observations, 1) : length(observations[1])
    firstlevelreplicates || all(c -> length(c) == nreplicates, observations) || throw(ErrorException("All cells must have the same number of replicates."))
    ncells = Int.(length(observations) / (firstlevelreplicates ? nreplicates : 1))
    nfactorlevels = firstlevelreplicates ? [size(observations)...][Not(1)] : [size(observations)...]

    crossedfactornames = factornames[factortypes .≠ nested]
    nestedfactornames = factornames[factortypes .== nested]

    #10kb allocated before this point
    anovakernel(observations, nreplicates, ncells, nnestedfactors, ncrossedfactors, nfactorlevels, crossedfactortypes, crossedfactornames, nestedfactornames)
end

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

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

function anova(observations::AbstractVector{T}, factorassignments::AbstractVector{<:AbstractVector}, factortypes::Vector{FactorType} = FactorType[]; factornames::Vector{<:AbstractString} = String[]) where {T <: Number}
    length(observations) > 0 || return
    nfactors = length(factorassignments)
    N = length(observations)
    N % nfactors == 0 || error("Design is unbalanced.")
    all(length.(factorassignments) .== N) || error("Each observation must have an assignment for each factor.")

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

    if !(isa(factorassignments, Number)) || any(maximum.(factorlevels) .> nfactorlevels)
        compressedfactorlevels = [1:i for i ∈ nfactorlevels]
        factorlevelremapping = [factorlevels[i] .=> compressedfactorlevels[i] for i ∈ 1:nfactors]
        factorassignments = [replace(factorassignments[i], factorlevelremapping[i]...) for i ∈ 1:nfactors]
    end

    nreplicates = Int(N / prod(nfactorlevels))

    nlevels = [nreplicates; nfactorlevels]
    sortorder = sortperm(repeat(1:nreplicates, Int(N / nreplicates)) .+
                         sum([factorassignments[i] .* prod(nlevels[1:i]) for i ∈ 1:nfactors]))
    observationsmatrix = reshape(observations[sortorder], nlevels...)
    anova(observationsmatrix, factortypes, factornames = factornames, hasreplicates = nreplicates > 1)
end

function validate(factortypes::Vector{FactorType}, factornames::Vector{<:AbstractString}, nfactors)
    if !isempty(factortypes)
        length(factortypes) == nfactors || error("factortypes must have an entry for each factor.")
        nested ∉ factortypes || factortypes[1:count(t -> t == nested, factortypes)] |> unique |> length == 1 || error("nested entries must come before crossed factors")
    end
    if !isempty(factornames)
        nfactors ≤ 26 || error("Can only automatically name up to 26 factors. Provide names explicitly.")
        length(factornames) == nfactors || error("factornames must have an entry for each factor.")
    end
end

function anovakernel(observations, nreplicates, ncells, nnestedfactors, ncrossedfactors, nfactorlevels, crossedfactortypes, crossedfactornames, nestedfactornames)
    N = ncells * nreplicates
    nfactors = nnestedfactors + ncrossedfactors

    # collapse replicate dimension
    cellsums = eltype(observations) <: Number && nreplicates == 1 ? observations : sumfirstdim(observations)
    constant = sum(cellsums) ^ 2 / N
    total = totalcalc(observations, N, constant)
    amongallnested, crossedcellsums, ncrossedfactorlevels, nnestedfactorlevels = amongnestedfactorscalc(cellsums, nfactorlevels, nnestedfactors, nreplicates, constant)
    cells = cellscalc(cellsums, nreplicates, ncells, constant)

    npercrossedcell = nreplicates * prod(nnestedfactorlevels)

    crossedfactors = factorscalc(crossedcellsums, ncrossedfactors, ncrossedfactorlevels, N, constant, crossedfactornames) # 3kb allocated here, possibly can't be avoided
    interactions = interactionscalc(cells, crossedcellsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, constant, crossedfactornames, npercrossedcell) # 4 kb allocated here, down from 9 kb
    nestedfactors = nestedfactorscalc(amongallnested, nnestedfactors, crossedfactors, interactions, nestedfactornames)

    reverse!(crossedfactors)
    reverse!(ncrossedfactorlevels)
    if ncrossedfactors == 3 && nreplicates > 1
        interactions[1:3] = reverse(interactions[1:3])
    end

    error = errorcalc(total, amongallnested, cells, [crossedfactors; interactions[1:end-1]], nnestedfactors, nreplicates) # could probably condense the otherfactor here
    numerators = getnumerators(crossedfactors, ncrossedfactors, nnestedfactors, nestedfactors, interactions)
    denominators = getdenominators(interactions, nestedfactors, error, crossedfactortypes, ncrossedfactors, nnestedfactors) # 2 kb allocated here

    # drop least significant test if nreplicates == 1; lowest nesting level. Lowest interaction previously dropped if no nested factors. Should look to do same for nesting.
    if nnestedfactors > 0 && nreplicates == 1
        droppedfactor = pop!(numerators)
        pop!(denominators)
    end

    # perform test
    results = ftest.(numerators, denominators)

    crossedcellmeans = crossedcellsums ./ npercrossedcell
    results = effectsizescalc(results, denominators, total, ncrossedfactors, npercrossedcell, ncrossedfactorlevels, crossedfactortypes, nnestedfactors, nnestedfactorlevels, nreplicates) # 12 kb allocated
    data = AnovaData([total; results], total, ncrossedfactors, ncrossedfactorlevels, npercrossedcell, crossedfactors, denominators[1:ncrossedfactors], crossedcellmeans) # 1 kb allocated here
    nnestedfactors > 0 && nreplicates == 1 && push!(data.effects, droppedfactor)
    error.df > 0 && push!(data.effects, error)

    return data
end

function sumfirstdim(observations::T) where {T <: AbstractArray{<:AbstractVector}}
    map(sumfirstdim, observations)
end

function sumfirstdim(observations::T) where {T <: AbstractArray{<:Number}}
    dropdims(sum(observations, dims = 1), dims = 1)
end

function sumfirstdim(observations::T) where {T <: AbstractVector{<:Number}}
    sum(observations, dims = 1) |> first
end

function totalcalc(observations, N, C)
    ss = sum(c -> sum(c.^2), observations) - C
    df = N - 1
    AnovaValue(totalname, ss, df)
end

function factorscalc(cellsums, nfactors, nfactorlevels, N, C, factornames)
    factorindices = 1:nfactors
    ss = map(i -> sum(sum(cellsums, dims = factorindices[Not(i)]) .^ 2) / (N / nfactorlevels[i]), factorindices) .- C
    df = nfactorlevels .- 1
    AnovaFactor.(factornames, ss, df)
end

function cellscalc(cellsums, nreplicates, ncells, C)
    ss = sum(cellsums .^ 2) / nreplicates - C
    df = ncells - 1
    AnovaValue(cellsname, ss, df)
end

function errorcalc(total, amongallnested, cells, otherfactors, nnestedfactors, nreplicates)
    if nnestedfactors > 0
        otherfactor = amongallnested[1]
        name = errorname
    elseif nreplicates > 1
        otherfactor = cells
        name = errorname
    else
        otherfactor = AnovaValue("", sum(f -> f.ss, otherfactors), sum(f -> f.df, otherfactors))
        name = remaindername
    end

    ss = total.ss - otherfactor.ss
    df = total.df - otherfactor.df
    AnovaFactor(name, ss, df)
end

function amongnestedfactorscalc(cellsums, nfactorlevels, nnestedfactors, nreplicates, C)
    nestedsums = cellsums
    nlowerfactorlevels = 1
    nupperfactorlevels = nfactorlevels
    amongallnested = Vector{AnovaValue}(undef, nnestedfactors)
    for i ∈ 1:nnestedfactors
        # collapse each nested factor
        ss = sum(nestedsums .^ 2 ./ (nreplicates * prod(nlowerfactorlevels))) - C
        df = prod(nupperfactorlevels) - 1
        amongallnested[i] = AnovaValue(ss, df)
        nlowerfactorlevels = nfactorlevels[1:i]
        nupperfactorlevels = nfactorlevels[(i+1):end]
        nestedsums = sumfirstdim(nestedsums)
    end
    amongallnested, nestedsums, nupperfactorlevels, nlowerfactorlevels
end

function interactionscalc(cells, nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactorlabels, npercrossedcell)
    interactions = AnovaFactor[]
    if (ncrossedfactors == 2 && npercrossedcell > 1) || ncrossedfactors == 3
        pairwise = pairwisecalc(nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactorlabels)
        push!(interactions, pairwise[1,2])
        if ncrossedfactors == 3
            push!(interactions, pairwise[1,3])
            push!(interactions, pairwise[2,3])
            if npercrossedcell > 1
                threewise = threewisecalc(cells, crossedfactors, pairwise, crossedfactorlabels)
                push!(interactions, threewise)
            end
        end
    end
    interactions
end

function pairwisecalc(nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactorlabels)
    pairwise = Array{AnovaFactor,2}(undef, ncrossedfactors, ncrossedfactors)
    factorindexes = 1:ncrossedfactors
    for i ∈ factorindexes
        for j ∈ (i+1):ncrossedfactors
            otherfactorindexes = factorindexes[Not([i, j])]
            ss = sum(sum(nestedsums, dims = otherfactorindexes) .^ 2 ./ (prod(ncrossedfactorlevels[otherfactorindexes]) * prod(nnestedfactorlevels) * nreplicates)) - C - crossedfactors[i].ss - crossedfactors[j].ss
            df = crossedfactors[i].df * crossedfactors[j].df
            pairwise[i,j] = pairwise[j,i] = AnovaFactor("$(crossedfactorlabels[j]) × $(crossedfactorlabels[i])", ss, df)
        end
    end
    pairwise
end

function threewisecalc(cells, factors, pairwise, crossedfactorlabels)
    ss = cells.ss - sum(f -> f.ss, factors) - pairwise[1,2].ss - pairwise[1,3].ss - pairwise[2,3].ss
    df = prod(f -> f.df, factors)
    AnovaFactor("$(crossedfactorlabels[3]) × $(crossedfactorlabels[2]) × $(crossedfactorlabels[1])", ss, df)
end

function nestedfactorscalc(amongallnested, nnestedfactors, crossedfactors, interactions, nestedfactorlabels)
    nestedfactors = []
    if nnestedfactors > 0
        otherfactors = [crossedfactors; interactions]
        otherss = sum(f -> f.ss, otherfactors)
        otherdf = sum(f -> f.df, otherfactors)
        for i ∈ nnestedfactors:-1:1
            ss = amongallnested[i].ss - otherss
            df = amongallnested[i].df - otherdf
            otherss += ss
            otherdf += df
            push!(nestedfactors, AnovaFactor(nestedfactorlabels[i], ss, df))
        end
    end
    nestedfactors
end

function getnumerators(crossedfactors, ncrossedfactors, nnestedfactors, nestedfactors, interactions)
    numerators = copy(crossedfactors)
    if ncrossedfactors > 1 && length(interactions) > 0
        push!(numerators, interactions[1])
        ncrossedfactors > 2 && push!(numerators, interactions[2], interactions[3], interactions[4])
    end
    nnestedfactors > 0 && append!(numerators, nestedfactors)
    numerators
end

function getdenominators(interactions, nestedfactors, error, crossedfactortypes, ncrossedfactors, nnestedfactors)
    nesteddenominators = Vector{AnovaFactor}(undef, nnestedfactors)

    if nnestedfactors > 0
        crossedbasedenominator = nestedfactors[1]
        nesteddenominators[1:end-1] = nestedfactors[2:end]
        nesteddenominators[end] = error
    else
        crossedbasedenominator = error
    end

    if ncrossedfactors == 1
        denominators = [crossedbasedenominator]
    elseif length(interactions) == 0
        denominators = repeat([crossedbasedenominator], ncrossedfactors)
    else # 2 or 3 crossed factors
        ninteractions = ncrossedfactors == 2 ? 1 : 4
        crosseddenominators = Vector{AnovaFactor}(undef, ncrossedfactors)
        pairinteractiondenominators = Vector{AnovaFactor}(undef, ninteractions - 1)

        if all(f -> f == fixed, crossedfactortypes)
            crosseddenominators = repeat([crossedbasedenominator], ncrossedfactors)
            pairinteractiondenominators = repeat([crossedbasedenominator], ninteractions - 1)

        elseif all(f -> f == random, crossedfactortypes)
            if ncrossedfactors == 2
                crosseddenominators = repeat(interactions, ncrossedfactors)

            else # 3 crossed factors
                crosseddenominators = Vector{AnovaFactor}(undef, 3)
                for i ∈ 1:3
                    j1, j2 = (1:3)[Not(i)]
                    crosseddenominators[i] = threewayinteraction(interactions[i + j1 - 2], interactions[i + j2 - 2], interactions[4])
                end
                pairinteractiondenominators = repeat([interactions[4]], ninteractions - 1)

            end

        else # has both random and fixed factors
            if ncrossedfactors == 2
                crosseddenominators = map(f -> f == fixed ? interactions[1] : crossedbasedenominator, crossedfactortypes)

            else # 3 crossed factors
                randomcount = count(f -> f == random, crossedfactortypes)
                if randomcount == 1
                    randomi = findfirst(f -> f == random, crossedfactortypes)
                    fixedi = (1:3)[Not(randomi)]
                    crosseddenominators = Vector{AnovaFactor}(undef, 3)
                    crosseddenominators[fixedi] = [interactions[i + randomi - 2] for i ∈ fixedi]
                    crosseddenominators[randomi] = crossedbasedenominator

                    fixedinteractioni = sum(fixedi) - 2
                    pairinteractiondenominators = Vector{AnovaFactor}(undef, 3)
                    pairinteractiondenominators[fixedinteractioni] = interactions[4]
                    pairinteractiondenominators[(1:3)[Not(fixedinteractioni)]] .= crossedbasedenominator

                else # 2 random factors
                    fixedi = findfirst(f -> f == fixed, crossedfactortypes)
                    randomi = (1:3)[Not(fixedi)]
                    crosseddenominators = Vector{AnovaFactor}(undef, 3)
                    crosseddenominators[fixedi] = threewayinteraction(interactions[fixedi + randomi[1] - 2], interactions[fixedi + randomi[2] - 2], interactions[4])
                    crosseddenominators[randomi] .= interactions[randomi[1] + randomi[2] - 2]

                    randominteractioni = sum(randomi) - 2
                    pairinteractiondenominators = Vector{AnovaFactor}(undef, 3)
                    pairinteractiondenominators[(1:3)[Not(randominteractioni)]] .= interactions[4]
                    pairinteractiondenominators[randominteractioni] = crossedbasedenominator

                end
            end
        end
        denominators = [crosseddenominators; pairinteractiondenominators; crossedbasedenominator]
    end
    append!(denominators, nesteddenominators)
    denominators
end

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

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

function effectsizescalc(results, denominators, total, ncrossedfactors, npercrossedcell, ncrossedfactorlevels, crossedfactortypes, nnestedfactors, nnestedfactorlevels, nreplicates)
    differences = [results[i].ms - denominators[i].ms for i ∈ eachindex(results)] # 1 kb between this line and next
    crossedfactordfs = [r.df for r ∈ results[1:ncrossedfactors]]

    if nreplicates == 1 && nnestedfactors > 0
        nnestedfactors -= 1
        nnestedfactorlevels = nnestedfactorlevels[1:(end-1)]
    end

    if ncrossedfactors == 1
        if nnestedfactors == 0
            ω² = [(results[1].ss - results[1].df * denominators[1].ms) / (total.ss + denominators[1].ms)]
        else
            effectdenominators = repeat([nreplicates], nnestedfactors + 1)
            nfactorlevels = [ncrossedfactorlevels; nnestedfactorlevels]
            effectdenominators[1] *= prod(nfactorlevels)
            factors = ones(Int, nnestedfactors + 1)
            factors[1] = crossedfactordfs[1]
            for i ∈ 2:nnestedfactors
                effectdenominators[2:(end - i + 1)] .*= nfactorlevels[end - i + 2]
            end
            σ² = factors .* differences ./ effectdenominators
            σ²total = sum(σ²) + denominators[end].ms
            ω² = σ² ./ σ²total
        end
    else
        if ncrossedfactors == 2 # this whole block not quite 1 kb
            if npercrossedcell > 1
                interactionindexes = [[1,2]]
                imax = 3
            else
                interactionindexes = []
                imax = 2
            end
        else
            if npercrossedcell > 1
                interactionindexes = [[1,2], [1,3], [2,3], [1,2,3]]
                imax = 7
            else
                interactionindexes = [[1,2], [1,3], [2,3]]
                imax = 6
            end
        end

        icrossed = 1:ncrossedfactors # this whole block 1 kb
        iother = ncrossedfactors < imax ? ((ncrossedfactors + 1):imax) : []
        factors = Vector{Int}(undef, imax)
        factors[icrossed] = [crossedfactortypes[i] == fixed ? crossedfactordfs[i] : 1 for i ∈ icrossed]
        factors[iother] = [prod(factors[x]) for x ∈ interactionindexes]

        effectsdenominators = repeat([npercrossedcell], imax)
        israndom = [x == random for x ∈ crossedfactortypes] # Originally used broadcasted equality (.==) but causes high allocations as of 1.3.0-rc3
        isfixed = [x == fixed for x ∈ crossedfactortypes]
        crossedeffectsdenominators = effectsdenominators[icrossed]
        crossedeffectsdenominators[isfixed] .*= prod(ncrossedfactorlevels)
        crossedeffectsdenominators[israndom] .*= [prod(ncrossedfactorlevels[Not(i)]) for i ∈ icrossed[israndom]]
        effectsdenominators[icrossed] = crossedeffectsdenominators
        effectsdenominators[iother] .*= [prod(ncrossedfactorlevels[Not(icrossed[israndom] ∩ x)]) for x ∈ interactionindexes] # 7kb - set intersection is 1kb, has to be done for each interaction

        σ² = factors .* differences[1:imax] ./ effectsdenominators

        if nnestedfactors > 0
            nestedrange = (length(results) .- nnestedfactors .+ 1):length(results)
            nestedeffectdenominators = repeat([nreplicates], nnestedfactors)
            for i ∈ 1:(nnestedfactors - 1)
                nestedeffectdenominators[1:(end - i + 1)] .*= nnestedfactorlevels[end - i + 2]
            end
            σ²nested = differences[nestedrange] ./ nestedeffectdenominators
            σ² = [σ²; σ²nested]
        end

        σ²total = sum(σ²) + denominators[end].ms
        ω² = σ² ./ σ²total
    end
    AnovaResult.(results, ω²)
end


text:
struct CSV_Foo
end

struct CSVString
    s::String
end

Base.parse(::Type{CSVString}, s::String) = CSVString(s)
Base.zero(::Type{CSVString}) = CSVString("")

struct Dec64
    x::Float64
end

Base.parse(::Type{Dec64}, s::String) = Dec64(parse(Float64, s))
Base.zero(::Type{Dec64}) = Dec64(0.0)

@testset "CSV.File basics" begin

#test on non-existent file
@test_throws ArgumentError CSV.File("");

#test where datarow > headerrow
@test_throws ArgumentError CSV.File(joinpath(dir, "test_no_header.csv"); datarow=1, header=2);

#test bad types
@test_throws CSV.Error CSV.File(joinpath(dir, "test_float_in_int_column.csv"); types=[Int, Int, Int], strict=true)

# Integer overflow; #100
@test_throws CSV.Error CSV.File(joinpath(dir, "int64_overflow.csv"); types=[Int8], strict=true)

# #172
@test_throws ArgumentError CSV.File(joinpath(dir, "test_newline_line_endings.csv"), types=Dict(1=>Integer))

# #289
tmp = CSV.File(IOBuffer(" \"a, b\", \"c\" "), datarow=1) |> columntable
@test length(tmp) == 2
@test length(tmp[1]) == 1
@test tmp.Column1[1] == "a, b"
@test tmp.Column2[1] == "c"

tmp = CSV.File(IOBuffer(" \"2018-01-01\", \"1\" ,1,2,3"), datarow=1) |> columntable
@test length(tmp) == 5
@test length(tmp[1]) == 1
@test tmp.Column1[1] == Date(2018, 1, 1)
@test tmp.Column2[1] == 1
@test tmp.Column3[1] == 1
@test tmp.Column4[1] == 2
@test tmp.Column5[1] == 3

# #329
f = CSV.File(joinpath(dir, "test_types.csv"), types=Dict(:string=>Union{Missing,DateTime}), silencewarnings=true)
@test f.string[1] === missing

# #352
@test_throws KeyError first(CSV.File(joinpath(dir, "test_types.csv"))).a

@test_throws ArgumentError CSV.File(IOBuffer("a\0b\n1\02\n"); delim='\0')
@test_throws ArgumentError CSV.File(IOBuffer("a\0b\n1\02\n"); delim="\0")

f = CSV.File(IOBuffer("a,b\n1,2\n"); header=3)
@test isempty(f.types)
@test f.rows == 0

f = CSV.File(IOBuffer("a,b\n1,2\n"); datarow=3)
@test f.names == [:a, :b]
@test f.rows == 0

f = CSV.File(IOBuffer("a,b\n1,2\n"); skipto=3)
@test f.names == [:a, :b]
@test f.rows == 0

f = CSV.File(IOBuffer("a,b\n1,2\n"); limit=0)
@test f.names == [:a, :b]
@test f.rows == 0

# delim same as quotechar
f = CSV.File(IOBuffer("a\"b\n1\"2\n"); delim='"')
@test f.rows == 0
@test f.cols == 2
@test f.types == [Missing, Missing]

# delim same as quotechar w/ quoted field
f = CSV.File(IOBuffer("a\"b\n1\"\"2\"\n"); delim='"')
@test f.rows == 0
@test f.cols == 2
@test f.types == [Missing, Missing]

# 387
f = CSV.File(IOBuffer("a,b\n1,1\n1,2\n1,3\n"))

for rows in Iterators.partition(f, 2)
    for row in rows
        # get last property
        getproperty(row, :b)
        @test row.a == 1
    end
end

# 391
rows = collect(f)
@test rows[1].a == 1
@test rows[1].a == 1

# 459
rows = collect(CSV.File(joinpath(dir, "time.csv"); dateformat="H:M:S"))
@test rows[1].time == Time(0)
@test rows[2].time == Time(0, 10)

# 388
f = CSV.File(joinpath(dir, "GSM2230757_human1_umifm_counts.csv"); threaded=false)
@test length(f.names) == 20128
@test length(f) == 3

# 386
f = CSV.File(IOBuffer("fullVisitorId,PredictedLogRevenue\n18966949534117,0\n39738481224681,0\n"), limit=3)
@test f.rows == 2
rows = collect(f)
@test rows[1].fullVisitorId == 18966949534117
@test rows[2].fullVisitorId == 39738481224681

f = CSV.File(IOBuffer("col1,col2\n1.0,hi"), limit=3)
@test (length(f), length(f.names)) == (1, 2)

# 288
f = CSV.File(IOBuffer("x\n\",\"\n\",\""))
@test (length(f), length(f.names)) == (2, 1)

# delimiter auto-detection
f = CSV.File(IOBuffer("x\n1\n2\n"))
@test (length(f), length(f.names)) == (2, 1)

f = CSV.File(IOBuffer("x,y\n1,\n2,\n"))
@test (length(f), length(f.names)) == (2, 2)

f = CSV.File(IOBuffer("x y\n1 \n2 \n"))
@test (length(f), length(f.names)) == (2, 2)

f = CSV.File(IOBuffer("x\ty\n1\t\n2\t\n"))
@test (length(f), length(f.names)) == (2, 2)

f = CSV.File(IOBuffer("x|y\n1|\n2|\n"))
@test (length(f), length(f.names)) == (2, 2)

# type promotion
# int => float
f = CSV.File(IOBuffer("x\n1\n3.14"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === 1.0
@test f.x[2] === 3.14

# int => missing
f = CSV.File(IOBuffer("x\n1\n\n"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] == 1
@test f.x[2] === missing

# missing => int
f = CSV.File(IOBuffer("x\n\n1\n"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === missing
@test f.x[2] == 1

# missing => int => float
f = CSV.File(IOBuffer("x\n\n1\n3.14\n"))
@test (length(f), length(f.names)) == (3, 1)
@test f.x[1] === missing
@test f.x[2] === 1.0
@test f.x[3] === 3.14

# int => missing => float
f = CSV.File(IOBuffer("x\n1\n\n3.14\n"))
@test (length(f), length(f.names)) == (3, 1)
@test f.x[1] === 1.0
@test f.x[2] === missing
@test f.x[3] === 3.14

# int => float => missing
f = CSV.File(IOBuffer("x\n1\n3.14\n\n"))
@test (length(f), length(f.names)) == (3, 1)
@test f.x[1] === 1.0
@test f.x[2] === 3.14
@test f.x[3] === missing

# int => string
f = CSV.File(IOBuffer("x\n1\nabc"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] == "1"
@test f.x[2] == "abc"

# float => string
f = CSV.File(IOBuffer("x\n3.14\nabc"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] == "3.14"
@test f.x[2] == "abc"

# int => float => string
f = CSV.File(IOBuffer("x\n1\n3.14\nabc"))
@test (length(f), length(f.names)) == (3, 1)
@test f.x[1] == "1"
@test f.x[2] == "3.14"
@test f.x[3] == "abc"

# missing => int => float => string
f = CSV.File(IOBuffer("x\n\n1\n3.14\nabc"))
@test (length(f), length(f.names)) == (4, 1)
@test f.x[1] === missing
@test f.x[2] == "1"
@test f.x[3] == "3.14"
@test f.x[4] == "abc"

# missing => catg
f = CSV.File(IOBuffer("x\n\na\n"), pool=true)
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === missing
@test f.x[2] == "a"

# catg => missing
f = CSV.File(IOBuffer("x\na\n\n"), pool=true)
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] == "a"
@test f.x[2] === missing

# catg => string
f = CSV.File(IOBuffer("x\na\nb\na\nb\na\nb\na\nc\nd\ne\nf\ng\nh\ni\nj\nk\nl\nn\nm\no\np\nq\nr\n"), pool=0.5)
@test typeof(f.x) == Vector{String}

# a few corner cases for escape strings
f = CSV.File(IOBuffer("\"column name with \"\" escape character inside\"\n1\n"))
@test keys(f[1])[1] == Symbol("column name with \" escape character inside")

f = CSV.File(IOBuffer("\"column name with \"\" escape character inside\",1\n,2"), transpose=true)
@test keys(f[1])[1] == Symbol("column name with \" escape character inside")
@test keys(f[1])[2] == :Column2

f = CSV.File(IOBuffer("x\na\nb\n\"quoted field with \"\" escape character inside\"\n"), pool=true)
@test f.x[1] == "a"
@test f.x[3] == "quoted field with \" escape character inside"

f = CSV.File(IOBuffer("x\na\nb\n\"quoted field with \"\" escape character inside\"\n"), pool=true)
@test f.x[1] == "a"
@test f.x[3] == "quoted field with \" escape character inside"

# invalid quoted field is fatal error
@test_throws CSV.Error CSV.File(IOBuffer("x\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.File(IOBuffer("x\nhey\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.File(IOBuffer("x\n\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.File(IOBuffer("x\n1\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.File(IOBuffer("x\n1.0\n\"quoted field that never ends"))
@test_throws CSV.Error CSV.File(IOBuffer("x\na\n\"quoted field that never ends"), pool=true)

# invalid integer
f = CSV.File(IOBuffer("x\nabc\n"), type=Int)
@test (length(f), length(f.names)) == (1, 1)
@test f.x[1] === missing

@test_throws CSV.Error CSV.File(IOBuffer("x\nabc\n"), type=Int, strict=true)

# transpose corner cases
f = CSV.File(IOBuffer("x,y,1\nx2,y2,2\n"), transpose=true, header=2)
@test f.names == [:y, :y2]
@test (length(f), length(f.names)) == (1, 2)
@test f.y[1] == 1
@test f.y2[1] == 2

f = CSV.File(IOBuffer("x,y,1\nx2,y2,2\n"), transpose=true, header=1, datarow=3)
@test f.names == [:x, :x2]
@test (length(f), length(f.names)) == (1, 2)
@test f.x[1] == 1
@test f.x2[1] == 2

f = CSV.File(IOBuffer("x,y,1\nx2,y2,2\n"), transpose=true, header=false, datarow=3)
@test f.names == [:Column1, :Column2]
@test (length(f), length(f.names)) == (1, 2)
@test f.Column1[1] == 1
@test f.Column2[1] == 2

f = CSV.File(IOBuffer(""), transpose=true, header=false)
@test (length(f), length(f.names)) == (0, 0)

f = CSV.File(IOBuffer(""), transpose=true, header=Symbol[])
@test (length(f), length(f.names)) == (0, 0)

# providing empty header vector
f = CSV.File(IOBuffer("x\nabc\n"), header=Symbol[])
@test f.names == [:Column1]

# Union{Bool, Missing}
f = CSV.File(IOBuffer("x\ntrue\n\n"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === true
@test f.x[2] === missing

# Union{Date, Missing}
f = CSV.File(IOBuffer("x\n2019-01-01\n\n"))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === Date(2019, 1, 1)
@test f.x[2] === missing

# use_mmap=false
f = @test_deprecated CSV.File(IOBuffer("x\n2019-01-01\n\n"), use_mmap=false)
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === Date(2019, 1, 1)
@test f.x[2] === missing

# types is Dict{String, Type}
f = CSV.File(IOBuffer("x\n2019-01-01\n\n"), types=Dict("x"=>Date))
@test (length(f), length(f.names)) == (2, 1)
@test f.x[1] === Date(2019, 1, 1)
@test f.x[2] === missing

# various CSV.File/CSV.Row properties
f = CSV.File(IOBuffer("int,float,date,datetime,bool,null,str,catg,int_float\n1,3.14,2019-01-01,2019-01-01T01:02:03,true,,hey,abc,2\n2,NaN,2019-01-02,2019-01-03T01:02:03,false,,there,abc,3.14\n"), pool=0.3)
@test Tables.istable(f)
@test Tables.rowaccess(typeof(f))
@test Tables.columnaccess(typeof(f))
@test Tables.schema(f) == Tables.Schema([:int, :float, :date, :datetime, :bool, :null, :str, :catg, :int_float], [Int64, Float64, Date, DateTime, Bool, Missing, String, String, Float64])
@test Tables.rows(f) === f
@test eltype(f) <: CSV.Row
row = first(f)
@test propertynames(row) == [:int, :float, :date, :datetime, :bool, :null, :str, :catg, :int_float]
@test row.int == 1
@test row.float == 3.14
@test row.date == Date(2019, 1, 1)
@test row.datetime == DateTime(2019, 1, 1, 1, 2, 3)
@test row.bool === true
@test row.null === missing
@test row.str == "hey"
@test row.catg == "abc"
@test typeof(row.catg) == String
@test row.int_float === 2.0
row = iterate(f, 2)[1]
@test row.int_float === 3.14

# 448
@test_throws ArgumentError CSV.File(IOBuffer("x\n1\n2\n3\n#4"), ignorerepeated=true)

# reported by oxinabox on slack; issue w/ manually specified pool column type and 0 rows
f = CSV.File(IOBuffer("x\n"), types=[CSV.PooledString])
@test (length(f), length(f.names)) == (0, 1)

f = CSV.File(IOBuffer("x\n"), types=[Union{CSV.PooledString, Missing}])
@test (length(f), length(f.names)) == (0, 1)

f = CSV.File(IOBuffer("x\n1\n2\n3\n#4"), comment="#")
@test length(f.x) == 3
@test f.x[end] == 3

# 453
@test_throws ArgumentError CSV.File(IOBuffer("x\n1\n2\n3\n#4"), types=[CSV_Foo])
@test_throws ArgumentError CSV.File(IOBuffer("x\n1\n2\n3\n#4"), types=Dict(:x=>CSV_Foo))

# 447
f = CSV.File(IOBuffer("a,b,c\n1,2,3\n\n"), ignoreemptylines=true)
@test (length(f), length(f.names)) == (1, 3)

f = CSV.File(IOBuffer("zip\n11111-1111\n"), dateformat = "y-m-dTH:M:S.s")
@test (length(f), length(f.names)) == (1, 1)

# Supporting commands across multiple platforms cribbed from julia/test/spawn.jl
catcmd = `cat`
havebb = false
if Sys.iswindows()
    busybox = download("https://frippery.org/files/busybox/busybox.exe", joinpath(tempdir(), "busybox.exe"))
    havebb = try # use busybox-w32 on windows, if available
        success(`$busybox`)
        true
    catch
        false
    end
    if havebb
        catcmd = `$busybox cat`
    end
end
f = @test_deprecated CSV.File(`$(catcmd) $(joinpath(dir, "test_basic.csv"))`)
@test columntable(f) == columntable(CSV.File(joinpath(dir, "test_basic.csv")))

#476
f = CSV.File(transcode(GzipDecompressor, Mmap.mmap(joinpath(dir, "randoms.csv.gz"))))
@test (length(f), length(f.names)) == (70000, 7)

f = CSV.File(IOBuffer("thistime\n10:00:00.0\n12:00:00.0"))
@test typeof(f.thistime[1]) <: Dates.Time
@test f.thistime[1] === Time(10)

# 530
f = CSV.File(IOBuffer(",column2\nNA,2\n2,3"), missingstrings=["NA"])
@test f.names == [:Column1, :column2]

# reported on slack from Kevin Bonham
f = CSV.File(IOBuffer("x\n01:02:03\n\n04:05:06\n"), delim=',')
@test isequal(f.x, [Dates.Time(1,2,3), missing, Dates.Time(4,5,6)])

# 566
f = CSV.File(IOBuffer("x\r\n1\r\n2\r\n3\r\n4\r\n5\r\n"), footerskip=3)
@test length(f) == 2
@test f[1][1] == 1

# 578
f = CSV.File(IOBuffer("h1234567890123456\t"^2262 * "lasthdr\r\n" *"dummy dummy dummy\r\n"* ("1.23\t"^2262 * "2.46\r\n")^10), datarow=3, threaded=false)
@test (length(f), length(f.names)) == (10, 2263)
@test all(x -> eltype(x) == Float64, f.columns)

# Date constructor throws which breaks CSV parsing error handling flow
f = CSV.File(IOBuffer("date\n2020-05-05\n2020-05-32"))

# multiple dateformats
f = CSV.File(IOBuffer("time,date,datetime\n10:00:00.0,04/16/2020,2020-04-16 23:14:00\n"), dateformats=Dict(2=>"mm/dd/yyyy", 3=>"yyyy-mm-dd HH:MM:SS"))
@test length(f) == 1
@test f[1].time == Dates.Time(10)
@test f[1].date == Dates.Date(2020, 4, 16)
@test f[1].datetime == Dates.DateTime(2020, 4, 16, 23, 14)

# custom types
f = CSV.File(
    IOBuffer("int8,uint32,bigint,bigfloat,dec64,csvstring\n1,2,170141183460469231731687303715884105727,3.14,1.02,hey there sailor\n2,,,,,\n");
    types=[Int8, UInt32, BigInt, BigFloat, Dec64, CSVString]
)
@test f.int8 isa Vector{Int8}
@test f.int8 == Int8[1, 2]
@test f.uint32 isa Vector{Union{UInt32, Missing}}
@test isequal(f.uint32, [0x00000002, missing])
@test f.bigint isa CSV.SVec2{BigInt}
@test isequal(f.bigint, [BigInt(170141183460469231731687303715884105727), missing])
@test f.bigfloat isa CSV.SVec2{BigFloat}
@test isequal(f.bigfloat, [parse(BigFloat, "3.14"), missing])
@test f.dec64 isa CSV.SVec{Dec64}
@test isequal(f.dec64, [Dec64(1.02), missing])
@test f.csvstring isa CSV.SVec2{CSVString}
@test isequal(f.csvstring, [CSVString("hey there sailor"), missing])

f = CSV.File(transcode(GzipDecompressor, Mmap.mmap(joinpath(dir, "randoms.csv.gz"))); types=[Int32, CSVString, String, Float64, Dec64, Date, DateTime])
@test f.id isa AbstractVector{Int32}
@test f.first isa AbstractVector{CSVString}
@test f.wage isa AbstractVector{Union{Missing, Dec64}}

f = CSV.File(joinpath(dir, "promotions.csv"); lazystrings=true)
@test eltype.(f.columns) == [Float64, Union{Missing, Int64}, Union{Missing, Float64}, String, Union{Missing, String}, String, String, Union{Missing, Int64}]
@test f.int_string isa CSV.LazyStringVector

f = CSV.File(joinpath(dir, "promotions.csv"); limit=7500, threaded=true)
@test length(f) == 7500

f = CSV.File(IOBuffer("1,2\r\n3,4\r\n\r\n5,6\r\n"); header=["col1", "col2"], ignoreemptylines=true)
@test f.col1 == [1, 3, 5]

f = CSV.File(joinpath(dir, "escape_row_starts.csv"); tasks=2)
@test length(f) == 10000
@test eltype(f.col1) == String
@test eltype(f.col2) == Int64

f = CSV.File(IOBuffer("col1\nhey\nthere\nsailor"); lazystrings=true)
@test f.col1 isa CSV.LazyStringVector
@test Tables.columnnames(f) == [:col1]
@test propertynames(f) == [:col1]
@test occursin("IOBuffer", CSV.getname(f))
@test CSV.getcols(f) == 1
@test Base.IndexStyle(f) == Base.IndexLinear()
@test f.col1 === Tables.getcolumn(f, 1)
@test columntable(f) == columntable(collect(f))
show(f)

f = CSV.File(joinpath(dir, "big_types.csv"); lazystrings=true, pool=false)
@test eltype(f.time) == Dates.Time
@test eltype(f.bool) == Bool
@test f.lazy isa CSV.LazyStringVector
@test eltype(f.lazy) == String
@test eltype(f.lazy_missing) == Union{String, Missing}

r = CSV.Rows(joinpath(dir, "big_types.csv"); lazystrings=false, types=[Dates.Time, Bool, String, Union{String, Missing}])
row = first(r)
@test row.time == Dates.Time(12)
@test row.bool
@test row.lazy == "hey"
@test row.lazy_missing === missing

@test CSV.File(IOBuffer("col1\n1")).col1 == [1]

rows = 0
chunks = CSV.Chunks(joinpath(dir, "promotions.csv"); lazystrings=true)
for chunk in chunks
    rows += length(chunk)
end
@test rows == 10000

# 668
buf = IOBuffer("""
       garbage
       a,b
       1,2
       """)
readline(buf)
f = CSV.File(buf, header=["A", "B"])
@test length(f) == 2
@test f.names == [:A, :B]
@test (f[1].A, f[1].B) == ("a", "b")
@test (f[2].A, f[2].B) == ("1", "2")

# 680: ensure typemap works with custom types
f = CSV.File(IOBuffer("a\n1\n2\n3"); typemap=Dict(Int64=>Int32))
@test f.a isa Vector{Int32}
f = CSV.File(IOBuffer("a\n1\n2\n3"); typemap=Dict(Int64=>String))
@test f.a isa Vector{String}

# 678
f = CSV.File(IOBuffer("""x,y
                                       a,b
                                       a,b
                                       a,b
                                       a,b
                                       a,b
                                       a,b
                                       a,b
                                       a,b
                                       a,b

                                       """))
@test f.x[end] === missing
@test f.y[end] === missing

# 679
f = CSV.File(IOBuffer("a,b,c\n1,2,3\n4,5,6\n"); select=["a"], types=Dict(2=>Int8))
@test f.a == [1, 4]

f = CSV.File(transcode(GzipDecompressor, Mmap.mmap(joinpath(dir, "randoms.csv.gz"))); types=Dict(:id=>Int32), select=["first"])
@test length(f) == 70000
@test eltype(f.first) == String

# 706
f = CSV.File(IOBuffer("a,b\n1,2"); types = Dict{Symbol,Type}(
               :a => CategoricalValue{String,UInt32},
               :b => Union{Missing, CategoricalValue{String,UInt32}},
       ))
@test f.a[1] == "1"
@test f.b[1] == "2"

end


text:
data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
    datadep"Mutagenesis"
end

@testset "Mutagenesis" begin
    train_x, train_y = Mutagenesis.traindata()
    test_x, test_y = Mutagenesis.testdata()
    val_x, val_y = Mutagenesis.valdata()

    @test length(train_x) == length(train_y) == 100
    @test length(test_x) == length(test_y) == 44
    @test length(val_x) == length(val_y) == 44
    # test that label is not contained in features
    @test !any(haskey.(train_x, :mutagenic))
    @test !any(haskey.(test_x, :mutagenic))
    @test !any(haskey.(val_x, :mutagenic))
    # test data is materialized
    @test train_x isa Vector{<:Dict}
    @test test_x isa Vector{<:Dict}
    @test val_x isa Vector{<:Dict}
end


text:
## test the main entry point in Celeste: the `infer` function
import JLD

import Celeste: Config
using Celeste.SDSSIO

"""
test infer with a single (run, camcol, field).
This is basically just to make sure it runs at all.
"""
function test_infer_single()
    # very small patch of sky that turns out to have 4 sources.
    # We checked that this patch is in the given field.
    box = ParallelRun.BoundingBox(164.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]
    result = ParallelRun.one_node_infer(field_triplets, datadir; box=box)
end

function test_load_active_pixels()
    ea, vp, catalog = gen_sample_star_dataset()

    # these images have 20 * 23 * 5 = 2300 pixels in total.
    # the star is bright but it doesn't cover the whole image.
    # it's hard to say exactly how many pixels should be active,
    # but not all of them, and not none of them.
    config = Config()
    config.min_radius_pix = 0
    Infer.load_active_pixels!(config, ea.images, ea.patches)

    # most star light (>90%) should be recorded by the active pixels
    num_active_photons = 0.0
    num_active_pixels = 0
    total_pixels = 0
    for n in 1:ea.N
        img = ea.images[n]
        p = ea.patches[1, n]
        H2, W2 = size(p.active_pixel_bitmap)
        total_pixels += length(img.pixels)
        for w2 in 1:W2, h2 in 1:H2
            # (h2, w2) index the local patch, while (h, w) index the image
            h = p.bitmap_offset[1] + h2
            w = p.bitmap_offset[2] + w2
            num_active_photons += img.pixels[h, w] - img.sky[h, w]
            num_active_pixels += 1
        end
    end

    @test 100 < num_active_pixels < total_pixels

    total_photons = 0.0
    for img in ea.images
        for w in 1:img.W, h in 1:img.H
            total_photons += img.pixels[h, w] - img.sky[h, w]
        end
    end

    @test num_active_photons <= total_photons  # sanity check
    @test num_active_photons > 0.9 * total_photons

    # super dim images
    for img in ea.images
        for w in 1:img.W, h in 1:img.H
            img.pixels[h, w] = img.sky[h, w]
        end
    end

    # only 2 pixels per image are within 0.6 pixels of the
    # source's center (10.9, 11.5)
    config = Config()
    config.min_radius_pix = 0.6
    Infer.load_active_pixels!(config, ea.images, ea.patches)

    for n in 1:ea.N
#  FIXME: is load active pixels off by (0.5, 0.5)?
#        @test sum(ea.patches[1,n].active_pixel_bitmap) == 2
    end
end

function test_patch_pixel_selection()
    ea, vp, catalog = gen_two_body_dataset();
    patches = Infer.get_sky_patches(ea.images, catalog; radius_override_pix=5);
    config = Config()
    config.min_radius_pix = 5
    Infer.load_active_pixels!(config, ea.images, patches, noise_fraction=Inf)

    for n in 1:ea.N
        # Make sure, for testing purposes, that the whole bitmap isn't full.
        for s in 1:ea.S
            @test sum(patches[s, n].active_pixel_bitmap) <
                  prod(size(patches[s, n].active_pixel_bitmap))
        end

        function patch_in_whole_image(p::SkyPatch)
            patch_image = zeros(size(ea.images[n].pixels))
            for h in 1:ea.images[n].H, w in 1:ea.images[n].W
                if Infer.is_pixel_in_patch(h, w, p)
                    patch_image[h, w] += 1
                end
            end
            return patch_image
        end
        patch_images = [ patch_in_whole_image(patches[s, 3]) for s in 1:ea.S ];
        for s in 1:ea.S
            @test sum(patch_images[s]) == sum(patches[s, n].active_pixel_bitmap)
        end

        H_min, W_min, H_max, W_max =
            Infer.get_active_pixel_range(patches, collect(1:ea.S), n);
        patch_image = zeros(H_max - H_min + 1, W_max - W_min + 1);

        for h in H_min:H_max, w in W_min:W_max, s in 1:ea.S
            p = patches[s, n]
            if Infer.is_pixel_in_patch(h, w, p)
                patch_image[h - H_min + 1, w - W_min + 1] += 1
            end
        end
        @test all(patch_image .== sum(patch_images)[H_min:H_max, W_min:W_max])
    end
end

@testset "test that we don't select a patch that is way too big" begin
    wd = pwd()
    cd(datadir)
    run(`make RUN=4114 CAMCOL=3 FIELD=127`)
    run(`make RUN=4114 CAMCOL=4 FIELD=127`)
    cd(wd)

    rcfs = [RunCamcolField(4114, 3, 127), RunCamcolField(4114, 4, 127)]
    strategy = PlainFITSStrategy(datadir)
    images = SDSSIO.load_field_images(strategy, rcfs)
    catalog = SDSSIO.read_photoobj_files(strategy, rcfs)

    entry_id = 429  # star at RA, Dec = (309.49754066435867, 45.54976572870953)

    neighbors = Infer.find_neighbors([entry_id,], catalog, images)[1]

    # there's a lot near this star, but not a lot that overlaps with it, see
    # http://skyserver.sdss.org/dr10/en/tools/explore/summary.aspx?id=0x112d1012607f050a
    @test length(neighbors) < 5
end

test_patch_pixel_selection()
test_load_active_pixels()
test_infer_single()


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents, Test, SafeTestsets, .Threads

# test_basics.jl has to be evaluated in global scope
@testset     "Basics"       begin include("test_basics.jl") end

@safetestset "Events"       begin include("test_events.jl") end
@safetestset "Clock"        begin include("test_clock.jl")  end
@safetestset "Units"        begin include("test_units.jl")  end
@safetestset "Threads"      begin include("test_threads.jl") end
@safetestset "Channel 1"    begin include("test_channels1.jl") end
@safetestset "Processes"    begin include("test_process.jl") end
@safetestset "Resources"    begin include("test_resources.jl") end
@safetestset "Utilities"    begin include("test_utils.jl") end
@safetestset "Macros"       begin include("test_macros.jl") end

if (VERSION ≥ v"1.3") && (nthreads() > 1)
    @safetestset "timer.jl" begin include("test_timer.jl") end
else
    println("... timer requires Julia ≥ 1.3, nthreads() > 1")
end

@safetestset "examples"     begin include("test_examples.jl") end

println(".... finished testing DiscreteEvents.jl ....")


text:
# BPF helpers

const ptr_sk_buff = API.pointertype(API.sk_buff)
const ptr_task_struct = API.pointertype(API.task_struct)

bpfconvert(x) = x
bpfconvert(x::AbstractBuffer) = pointer(x)

@inline function refwrap(f, x::T) where T
    x_ref = Ref{T}(x)
    GC.@preserve x_ref begin
        f(Base.unsafe_convert(Ptr{T}, x_ref))
    end
end

@inline probe_read(buf::AbstractSizedBuffer, addr::BufPtr) =
    bpfcall(API.probe_read, Clong, Tuple{BufPtr, UInt32, BufPtr}, pointer(buf), length(buf), addr)

@inline ktime_get_ns() = bpfcall(API.ktime_get_ns, UInt64)

@inline trace_printk(fmt::AbstractSizedBuffer, x...) = # TODO: Allow only 3 trailing arguments
    bpfcall(API.trace_printk, Clong, Tuple{BufPtr, UInt32, typeof(bpfconvert.(x))...}, pointer(fmt), length(fmt), map(bpfconvert, x)...)

@inline get_prandom_u32() = bpfcall(API.get_prandom_u32, UInt32)

@inline get_smp_processor_id() = bpfcall(API.get_smp_processor_id, UInt32)

@inline skb_store_bytes(skb::ptr_sk_buff, offset::UInt32, from::Ptr{Cvoid}, len::UInt32, flags::UInt64) =
    bpfcall(API.skb_store_bytes, Clong, Tuple{ptr_sk_buff, UInt32, Ptr{Cvoid}, UInt32, UInt64}, skb, offset, from, len, flags)
@inline skb_store_bytes(skb::ptr_sk_buff, offset, from, len, flags) = refwrap(from) do from_ref
    skb_store_bytes(skb, unsafe_trunc(UInt32, offset), reinterpret(Ptr{Cvoid}, from_ref), unsafe_trunc(UInt32, len), unsafe_trunc(UInt64, flags))
end
@inline l3_csum_replace(skb::ptr_sk_buff, offset::UInt32, from::UInt64, to::UInt64, size::UInt64) =
    bpfcall(API.l3_csum_replace, Clong, Tuple{ptr_sk_buff, UInt32, UInt64, UInt64, UInt64}, skb, offset, from, to, size)
@inline l3_csum_replace(skb::ptr_sk_buff, offset, from, to, size) =
    l3_csum_replace(skb, unsafe_trunc(UInt32, offset), unsafe_trunc(UInt64, from), unsafe_trunc(UInt64, to), unsafe_trunc(UInt64, size))
@inline l4_csum_replace(skb::ptr_sk_buff, offset::UInt32, from::UInt64, to::UInt64, flags::UInt64) =
    bpfcall(API.l4_csum_replace, Clong, Tuple{ptr_sk_buff, UInt32, UInt64, UInt64, UInt64}, skb, offset, from, to, flags)
@inline l4_csum_replace(skb::ptr_sk_buff, offset, from, to, flags) =
    l4_csum_replace(skb, unsafe_trunc(UInt32, offset), unsafe_trunc(UInt64, from), unsafe_trunc(UInt64, to), unsafe_trunc(UInt64, flags))

# TODO: tail_call

@inline clone_redirect(skb::ptr_sk_buff, ifindex::UInt32, flags::UInt64) =
    bpfcall(API.clone_redirect, Clong, Tuple{ptr_sk_buff, UInt32, UInt64}, skb, ifindex, flags)

function split_u64_u32(x::UInt64)
    lower = Base.unsafe_trunc(UInt32, x)
    upper = Base.unsafe_trunc(UInt32, x >> 32)
    return lower, upper
end

@inline get_current_pid_tgid() = split_u64_u32(bpfcall(API.get_current_pid_tgid, UInt64))
@inline get_current_uid_gid() = split_u64_u32(bpfcall(API.get_current_uid_gid, UInt64))
@inline get_current_comm(buf::AbstractSizedBuffer) =
    bpfcall(API.get_current_comm, Clong, Tuple{BufPtr, UInt32}, pointer(buf), length(buf))

@inline perf_event_output(ctx::T, map::M, flags, buf::AbstractSizedBuffer) where {T, M<:RTMap} =
    bpfcall(API.perf_event_output, Clong, Tuple{T, M, UInt64, BufPtr, UInt64}, ctx, map, unsafe_trunc(UInt64, flags), pointer(buf), length(buf))

@inline get_stackid(ctx::T, map::M, flags::Integer) where {T,M<:RTMap} =
    unsafe_trunc(Int32, bpfcall(API.get_stackid, Clong, Tuple{T, M, UInt64}, ctx, map, unsafe_trunc(UInt64,flags)))
@inline get_current_task() = bpfcall(API.get_current_task, ptr_task_struct)
@inline get_stack(ctx::T, buf::AbstractSizedBuffer, flags::UInt64) where {T} =
    bpfcall(API.get_stack, Clong, Tuple{T, BufPtr, UInt32, UInt64}, ctx, pointer(buf), length(buf), flags)

@inline probe_read_user(dest::T, size::Integer, unsafe_ptr::T) where {T<:Ptr} =
    bpfcall(API.probe_read_user, Clong, Tuple{T, UInt32, T}, dest, unsafe_trunc(UInt32, size), unsafe_ptr)
@inline probe_read_kernel(dest::T, size::Integer, unsafe_ptr::T) where {T<:Ptr} =
    bpfcall(API.probe_read_kernel, Clong, Tuple{T, UInt32, T}, dest, unsafe_trunc(UInt32, size), unsafe_ptr)

@inline get_task_stack(ctx::ptr_task_struct, buf::AbstractSizedBuffer, flags::UInt64) where {T} =
    bpfcall(API.get_task_stack, Clong, Tuple{ptr_task_struct, BufPtr, UInt32, UInt64}, ctx, pointer(buf), length(buf), flags)

# TODO: The rest!


text:
using Remarkable
using Test

@testset "Remarkable.jl" begin
   client = RemarkableClient(path_to_token = @__DIR__)
   @test Remarkable.token(client) != String(read(joinpath(@__DIR__, ".token")))
   list = list_items(client)
   @test list isa Remarkable.Collection
   @test list[1] isa Remarkable.Collection
   @test list[2] isa Remarkable.Document
   print_tree(list)
end


text:
using CausalityTools
using CausalityToolsBase
using TimeseriesSurrogates
using PerronFrobenius
using PyCall, Conda
using CrossMappings
using HypothesisTests
using Documenter, DocumenterMarkdown
#Conda.add("scipy")
using Plots
using DynamicalSystems
using SimpleDiffEq
using DynamicalSystems
using Distributions
using StaticArrays
using Statistics
using StatsBase
using TransferEntropy
using UncertainData

ENV["GKSwstype"] = "100"

PAGES = [
    #"index.md",
    #"FAQ.md",
    "overview.md",
    "Syntax overview" => "syntax_overview.md",
    "CHANGELOG.md",
    "System models" => [
        "SystemModels/SystemModels_overview.md",
        "SystemModels/ContinuousModels.md",
        "SystemModels/randomising.md",
        "SystemModels/abstract_types.md"
    ],

    "Causality tests" => [

        # Causal analysis
        "causalitytests/CausalAnalysis.md",

        # From uncertain data
        "causalitytests/causality_from_uncertain_data_naive.md",
        "causalitytests/causality_from_uncertain_data_binneddatacausalitytest.md",
        "causalitytests/causality_from_uncertain_data_strictlyincreasing_interpolated.md",
        "causalitytests/causality_from_uncertain_data_naive_constrained.md",
        "causalitytests/causality_from_uncertain_data_InterpolateBinTest.md",
        "causalitytests/causality_from_uncertain_data_RandomSequencesTest.md",
        "causalitytests/causality_from_uncertain_data.md",

        # From dynamical systems
        "causalitytests/causality_from_dynamical_systems.md",

        # Tests
        "causalitytests/causality_tests.md",
        "causalitytests/abstract_test_types.md",
        "causalitytests/causality_tests_overview.md",
        "causalitytests/ConvergentCrossMappingTest.md",
        "causalitytests/CrossMappingTest.md",
        "causalitytests/JointDistanceDistributionTest.md",
        "causalitytests/SMeasureTest.md",

        "causalitytests/NearestNeighbourMITest.md",
        "causalitytests/TransferOperatorGridTest.md",
        "causalitytests/VisitationFrequencyTest.md",
        "causalitytests/ApproximateSimplexIntersectionTest.md",
        "causalitytests/ExactSimplexIntersectionTest.md",

        "causalitytests/PredictiveAsymmetryTest.md",
        "causalitytests/NormalisedPredictiveAsymmetryTest.md"
    ],
    "CausalityToolsBase" => [
        "Discretization" => "causalitytoolsbase/discretization.md",
        "Delay reconstructions" => "causalitytoolsbase/delay_reconstructions.md"
    ],
    "Transfer operator estimation" => [
        "Transfer operator" => "perronfrobenius/transferoperator.md",
        "Invariant measure" => "perronfrobenius/invariantmeasure.md"
    ],
    "Transfer entropy" => [
        "Overview" => "transferentropy/overview_te.md",
        "Estimators" => "transferentropy/TE_estimators.md",
        "Estimating TE" => "transferentropy/estimating_TE.md",
        "Assigning marginals" => "transferentropy/assigning_marginals.md",
        "Generalised delay embedding" => "transferentropy/generalised_delay_reconstructions.md",
        "Convenience methods" => "transferentropy/convenience_methods_te.md",
        "Effect of discretization scheme" => "transferentropy/examples_TE_different_partitionings.md",
    ],
    "PredictiveAsymmetry" => [
        "Predictive asymmetry" => "PredictiveAsymmetry/predictive_asymmetry.md"
    ],
    "Distance based statistics" => [
        "CCM" => [
            "Overview" => "crossmappings/ccm/overview.md",
            "Cross mapping" => "crossmappings/ccm/crossmapping.md",
            "Converent cross mapping" => "crossmappings/ccm/convergentcrossmapping.md"
        ]
    ],

    "Joint distance distribution" => [
        "JointDistanceDistribution/joint_distance_distribution.md"
    ],

    "SMeasure" => [
        "S-measure" => "SMeasure/s_measure.md"
    ],
    "Worked example" => [
        "worked_examples/worked_example_transferentropy.md"
    ],

    "Simplices" => [
        "simplices/Simplices.md"
    ],

    "Example systems" => [
        "example_systems/noise.md",
        "example_systems/example_systems_overview.md",
        "example_systems/example_systems_discrete.md",
        "example_systems/example_systems_continuous.md"
    ],

    "Surrogate data" => [
        "surrogates/iaaft_docs.md",
        "surrogates/aaft_docs.md",
        "surrogates/randomphases_docs.md",
        "surrogates/randomamplitudes_docs.md",
        "surrogates/randomshuffle_docs.md",
        "surrogates/surrogates_overview.md"
    ],
    "UncertainData" => [
        "uncertaindata/BinnedDataCausalityTest.md"
    ],
    "Tutorials" => [
        "causality" => [
            "tutorials/list_of_tutorials.md",
            "tutorials/causality/Tutorial_BinnedResampling_RandomSequencesTest_PredictiveAsymmetryTest_ar1_unidir.md",
            "tutorials/causality/binned_uncertain_data/tutorial_BinnedDataCausalityTest_PredictiveAsymmetryTest_BinnedResampling.md",
            "tutorials/causality/binned_uncertain_data/tutorial_BinnedDataCausalityTest_PredictiveAsymmetryTest_BinnedMeanResampling.md"
        ],
        "tutorials/causality/tutorial_SMeasureTest.md"
    ]
]

# Some things need to be expanded before index.md for headers to work properly.
expandfirst = [
    "causalitytoolsbase/discretization.md",
    "causalitytoolsbase/delay_reconstructions.md",
    "simplices/Simplices.md",
    "causalitytests/causality_tests_overview.md",
    "perronfrobenius/transferoperator.md",
    "perronfrobenius/invariantmeasure.md",
    "causalitytests/NearestNeighbourMITest.md",
    "crossmappings/ccm/overview.md"
]

makedocs(
    sitename = "CausalityTools.jl documentation",
    modules = [CausalityTools, TransferEntropy, PerronFrobenius, CrossMappings, CausalityToolsBase, UncertainData, Simplices],
    format = DocumenterMarkdown.Markdown(),
    linkcheck = false,
    clean = true,
    expandfirst = expandfirst,
    pages = PAGES,
    highlightsig = true
)

if !Sys.iswindows()
    deploydocs(
        deps   = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4",
        "python-markdown-math", "pygments", "pymdown-extensions"),
        repo   = "github.com/kahaaga/CausalityTools.jl.git",
        target = "site",
        make = () -> run(`mkdocs build`)
    )
end