text:
using RawArray
using MAT
using HDF5
using BenchmarkTools

nfiles = 10000
n = 10
path = mktempdir() * "/"

println("Temporary files will be placed in $path.")
println("If code crashes, or is terminated early, but sure to delete that directory.")
function cleanup(n, nfiles)
	for i in 1:nfiles
		rm("$path/$i.ra", force=true)
		rm("$path/$i.h5", force=true)
		rm("$path/$i.mat", force=true)
	end
	rm("$path/big.ra", force=true)
	rm("$path/big.h5", force=true)
	rm("$path/small.h5", force=true)
	rm("$path/big.mat", force=true)
	rm("$path/small.mat", force=true)
end

function rasmalltest(x, n, nfiles)
	for i in 1:nfiles
		rawrite(x, "$path/$i.ra")
	end
	for i in 1:nfiles
		x = raread("$path/$i.ra")
	end
	return x[1]
end

function ralargetest(x, n, nfiles)
	rawrite(x, path*"big.ra")
	x = raread("$path/big.ra")
	return x[1]
end

function hdf5largetest(x, n, nfiles)
	h5write("$path/big.h5", "test/x", x)
	x = h5read("$path/big.h5", "test/x")
	return x[1]
end

function hdf5smalltestone(x, n, nfiles)
	h5open("$path/small.h5", "w") do file
		for i in 1:nfiles
			write(file, "a$i", x)  # alternatively, say "@write file A"
		end
	end
	c = h5open("$path/small.h5", "r") do file
		for i in 1:nfiles
			x = read(file, "a$i")
		end
	end
	return x[1]
end

function hdf5smalltestmany(x, n, nfiles)
	for i in 1:nfiles
		h5write("$path/$i.h5", "test/x", x)
	end
	for i in 1:nfiles
		x = h5read("$path/$i.h5", "test/x")
	end
	return x[1]
end

function matlargetest(x, n, nfiles)
	matwrite("$path/big.mat", Dict("x"=>x))
	x = matread("$path/big.mat")
	return x["x"]
end

function matsmalltestmany(x, n, nfiles)
	for i in 1:nfiles
		matwrite("$path/$i.mat", Dict("x"=> x))
	end
	for i in 1:nfiles
		x = matread("$path/$i.mat")
	end
	return x["x"]
end

xsmall = rand(n)
xbig = rand(nfiles, n, n, n)

print("MAT small separate:      ")
@time matsmalltestmany(xsmall, n, nfiles)
print("HDF5 small separate:     ")
@time hdf5smalltestmany(xsmall, n, nfiles)
print("HDF5 small together:     ")
@time hdf5smalltestone(xsmall, n, nfiles)
print("RawArray small separate: ")
@time rasmalltest(xsmall, n, nfiles)
print("MAT large:               ")
@time matlargetest(xbig, n, nfiles)
print("HDF5 large:              ")
@time hdf5largetest(xbig, n, nfiles)
print("RawArray large:          ")
@time ralargetest(xbig, n, nfiles)
cleanup(n, nfiles)

max_stars_repo_name: davidssmith/ra

max_issues_repo_name: davidssmith/ra

max_forks_repo_name: davidssmith/ra


text:
isdefined(Base, :__precompile__) && __precompile__()

module JSON

using Compat

export json # returns a compact (or indented) JSON representation as a string

include("Common.jl")
include("Parser.jl")

using .Common
import .Parser.parse

# These are temporary ways to bypass excess memory allocation
# They can be removed once types define their own serialization behaviour again
"Internal JSON.jl implementation detail; do not depend on this type."
immutable AssociativeWrapper{T} <: Associative{Symbol, Any}
    wrapped::T
    fns::Vector{Symbol}
end
AssociativeWrapper(x) = AssociativeWrapper(x, fieldnames(x))

const JSONPrimitive = Union{
        Associative, Tuple, AbstractArray, AbstractString, Integer,
        AbstractFloat, Bool, Void}

Base.getindex(w::AssociativeWrapper, s::Symbol) = getfield(w.wrapped, s)
Base.keys(w::AssociativeWrapper) = w.fns
Base.length(w::AssociativeWrapper) = length(w.fns)

"""
Return a value of a JSON-encodable primitive type that `x` should be lowered
into before encoding as JSON. Supported types are: `Associative` to JSON
objects, `Tuple` and `AbstractVector` to JSON arrays, `AbstractArray` to nested
JSON arrays, `AbstractString` to JSON string, `Integer` and `AbstractFloat` to
JSON number, `Bool` to JSON boolean, and `Void` to JSON null.

Extensions of this method should preserve the property that the return value is
one of the aforementioned types. If first lowering to some intermediate type is
required, then extensions should call `lower` before returning a value.

Note that the return value need not be *recursively* lowered—this function may
for instance return an `AbstractArray{Any, 1}` whose elements are not JSON
primitives.
"""
lower(a) = AssociativeWrapper(a)
lower(a::JSONPrimitive) = a

if isdefined(Base, :Dates)
    lower(s::Base.Dates.TimeType) = string(s)
end

lower(s::Symbol) = string(s)

if VERSION < v"0.5.0-dev+2396"
    lower(f::Function) = "function at $(f.fptr)"
end

lower(c::Char) = string(c)
lower(d::Type) = string(d)
lower(m::Module) = throw(ArgumentError("cannot serialize Module $m as JSON"))
lower(x::Real) = Float64(x)

const INDENT=true
const NOINDENT=false

@compat type State{I}
    indentstep::Int
    indentlen::Int
    prefix::AbstractString
    otype::Vector{Bool}
    (::Type{State{I}}){I}(indentstep::Int) = new{I}(indentstep,
                                                    0,
                                                    "",
                                                    Bool[])
end
State(indentstep::Int=0) = State{indentstep>0}(indentstep)

function set_state(state::State{INDENT}, operate::Int)
    state.indentlen += state.indentstep * operate
    state.prefix = " "^state.indentlen
end

set_state(state::State{NOINDENT}, operate::Int) = nothing

suffix(::State{INDENT})   = "\n"
suffix(::State{NOINDENT}) = ""

prefix(s::State{INDENT})  = s.prefix
prefix(::State{NOINDENT}) = ""

separator(::State{INDENT})   = ": "
separator(::State{NOINDENT}) = ":"

# short hand for printing suffix then prefix
printsp(io::IO, state::State{INDENT}) = Base.print(io, suffix(state), prefix(state))
printsp(io::IO, state::State{NOINDENT}) = nothing

function start_object(io::IO, state::State{INDENT}, is_dict::Bool)
    push!(state.otype, is_dict)
    Base.print(io, is_dict ? "{": "[", suffix(state))
    set_state(state, 1)
end

function start_object(io::IO, state::State{NOINDENT}, is_dict::Bool)
    Base.print(io, is_dict ? "{": "[")
end

function end_object(io::IO, state::State{INDENT}, is_dict::Bool)
    set_state(state, -1)
    pop!(state.otype)
    printsp(io, state)
    Base.print(io, is_dict ? "}": "]")
end

function end_object(io::IO, state::State{NOINDENT}, is_dict::Bool)
    Base.print(io, is_dict ? "}": "]")
end

function print_escaped(io::IO, s::AbstractString)
    @inbounds for c in s
        c <= '\x7f' ? Base.write(io, ESCAPED_ARRAY[@compat UInt8(c) + 0x01]) :
                      Base.print(io, c) #JSON is UTF8 encoded
    end
end

function print_escaped(io::IO, s::Compat.UTF8String)
    @inbounds for c in Vector{UInt8}(s)
        Base.write(io, ESCAPED_ARRAY[c + 0x01])
    end
end

function _writejson(io::IO, state::State, s::AbstractString)
    Base.print(io, '"')
    JSON.print_escaped(io, s)
    Base.print(io, '"')
end

# workaround for issue in Julia 0.5.x where Float32 values are printed as
# 3.4f-5 instead of 3.4e-5
if v"0.5-" <= VERSION < v"0.6.0-dev.788"
    _writejson(io::IO, state::State, s::Float32) = _writejson(io, state, Float64(s))
end

function _writejson(io::IO, state::State, s::Union{Integer, AbstractFloat})
    if isnan(s) || isinf(s)
        Base.print(io, "null")
    else
        Base.print(io, s)
    end
end

function _writejson(io::IO, state::State, n::Void)
    Base.print(io, "null")
end

function _writejson(io::IO, state::State, a::Nullable)
    if isnull(a)
        Base.print(io, "null")
    else
        _writejson(io, state, get(a))
    end
end

function _writejson(io::IO, state::State, a::Associative)
    if length(a) == 0
        Base.print(io, "{}")
        return
    end
    start_object(io, state, true)
    first = true
    for key in keys(a)
        first ? (first = false) : Base.print(io, ",", suffix(state))
        Base.print(io, prefix(state))
        _writejson(io, state, string(key))
        Base.print(io, separator(state))
        _writejson(io, state, a[key])
    end
    end_object(io, state, true)
end

function _writejson(io::IO, state::State, a::Union{AbstractVector,Tuple})
    if length(a) == 0
        Base.print(io, "[]")
        return
    end
    start_object(io, state, false)
    Base.print(io, prefix(state))
    i = start(a)
    !done(a,i) && ((x, i) = next(a, i); _writejson(io, state, x); )

    while !done(a,i)
        (x, i) = next(a, i)
        Base.print(io, ",")
        printsp(io, state)
        _writejson(io, state, x)
    end
    end_object(io, state, false)
end

function _writejson(io::IO, state::State, a)
    # FIXME: This fallback is harming performance substantially.
    # Remove this fallback when _print removed.
    if applicable(_print, io, state, a)
        Base.depwarn(
            "Overloads to `_print` are deprecated; extend `lower` instead.",
            :_print)
        _print(io, state, a)
    else
        _writejson(io, state, lower(a))
    end
end

# Note: Arrays are printed in COLUMN MAJOR format.
# i.e. json([1 2 3; 4 5 6]) == "[[1,4],[2,5],[3,6]]"
function _writejson{T, N}(io::IO, state::State, a::AbstractArray{T, N})
    lengthN = size(a, N)
    if lengthN > 0
        start_object(io, state, false)
        if VERSION <= v"0.3"
            newdims = ntuple(N - 1, i -> 1:size(a, i))
        else
            newdims = ntuple(i -> 1:size(a, i), N - 1)
        end
        Base.print(io, prefix(state))
        _writejson(io, state, Compat.view(a, newdims..., 1))

        for j in 2:lengthN
            Base.print(io, ",")
            printsp(io, state)
            _writejson(io, state, Compat.view(a, newdims..., j))
        end
        end_object(io, state, false)
    else
        Base.print(io, "[]")
    end
end

# this is _print() instead of _print because we need to support v0.3
# FIXME: drop the parentheses when v0.3 support dropped
"Deprecated way to overload JSON printing behaviour. Use `lower` instead."
function _print(io::IO, s::State, a::JSONPrimitive)
    Base.depwarn(
        "Do not call internal function `JSON._print`; use `JSON.print`",
        :_print)
    _writejson(io, s, a)
end

function print(io::IO, a, indent=0)
    _writejson(io, State(indent), a)
    if indent > 0
        Base.print(io, "\n")
    end
end

function print(a, indent=0)
    _writejson(STDOUT, State(indent), a)
    if indent > 0
        println()
    end
end

json(a, indent=0) = sprint(JSON.print, a, indent)

function parsefile{T<:Associative}(filename::AbstractString; dicttype::Type{T}=Dict{Compat.UTF8String, Any}, use_mmap=true)
    sz = filesize(filename)
    open(filename) do io
        s = use_mmap ? Compat.UTF8String(Mmap.mmap(io, Vector{UInt8}, sz)) : readstring(io)
        JSON.parse(s; dicttype=dicttype)
    end
end

end # module

ext: jl


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function rhs!(du, u, t,
              mesh::StructuredMesh{3}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @trixi_timeit timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @trixi_timeit timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Calculate interface fluxes
  @trixi_timeit timer() "interface flux" calc_interface_flux!(
    cache, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.surface_integral, dg)

  # Calculate boundary fluxes
  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(
    cache, u, t, boundary_conditions, mesh, equations, dg.surface_integral, dg)

  # Calculate surface integrals
  @trixi_timeit timer() "surface integral" calc_surface_integral!(
    du, u, mesh, equations, dg.surface_integral, dg, cache)

  # Apply Jacobian from mapping to reference element
  @trixi_timeit timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @trixi_timeit timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  return nothing
end

function calc_volume_integral!(du, u,
                               mesh::Union{StructuredMesh{3}, P4estMesh{3}},
                               nonconservative_terms::Val{false}, equations,
                               volume_integral::VolumeIntegralWeakForm,
                               dg::DGSEM, cache)
  @unpack derivative_dhat = dg.basis
  @unpack contravariant_vectors = cache.elements

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, k, element)

      flux1 = flux(u_node, 1, equations)
      flux2 = flux(u_node, 2, equations)
      flux3 = flux(u_node, 3, equations)

      # Compute the contravariant flux by taking the scalar product of the
      # first contravariant vector Ja^1 and the flux vector
      Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
      contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3
      for ii in eachnode(dg)
        multiply_add_to_node_vars!(du, derivative_dhat[ii, i], contravariant_flux1, equations, dg, ii, j, k, element)
      end

      # Compute the contravariant flux by taking the scalar product of the
      # second contravariant vector Ja^2 and the flux vector
      Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
      contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3
      for jj in eachnode(dg)
        multiply_add_to_node_vars!(du, derivative_dhat[jj, j], contravariant_flux2, equations, dg, i, jj, k, element)
      end

      # Compute the contravariant flux by taking the scalar product of the
      # third contravariant vector Ja^3 and the flux vector
      Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)
      contravariant_flux3 = Ja31 * flux1 + Ja32 * flux2 + Ja33 * flux3
      for kk in eachnode(dg)
        multiply_add_to_node_vars!(du, derivative_dhat[kk, k], contravariant_flux3, equations, dg, i, j, kk, element)
      end
    end
  end

  return nothing
end

# flux differencing volume integral on curvilinear hexahedral elements. Averaging of the
# mapping terms, stored in `contravariant_vectors`, is peeled apart from the evaluation of
# the physical fluxes in each Cartesian direction
@inline function split_form_kernel!(du::AbstractArray{<:Any,5}, u,
                                    element, mesh::Union{StructuredMesh{3}, P4estMesh{3}},
                                    nonconservative_terms::Val{false}, equations,
                                    volume_flux, dg::DGSEM, cache, alpha=true)
  # true * [some floating point value] == [exactly the same floating point value]
  # This can (hopefully) be optimized away due to constant propagation.
  @unpack derivative_split = dg.basis
  @unpack contravariant_vectors = cache.elements

  # Calculate volume integral in one element
  for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, j, k, element)

    # pull the contravariant vectors in each coordinate direction
    Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
    Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
    Ja3_node = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)

    # All diagonal entries of `derivative_split` are zero. Thus, we can skip
    # the computation of the diagonal terms. In addition, we use the symmetry
    # of the `volume_flux` to save half of the possible two-point flux
    # computations.

    # x direction
    for ii in (i+1):nnodes(dg)
      u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
      # pull the contravariant vectors and compute the average
      Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors,
                                             ii, j, k, element)
      Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii)
      # compute the contravariant sharp flux in the direction of the
      # averaged contravariant vector
      fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations)
      multiply_add_to_node_vars!(du, alpha * derivative_split[i, ii], fluxtilde1, equations, dg, i,  j, k, element)
      multiply_add_to_node_vars!(du, alpha * derivative_split[ii, i], fluxtilde1, equations, dg, ii, j, k, element)
    end

    # y direction
    for jj in (j+1):nnodes(dg)
      u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
      # pull the contravariant vectors and compute the average
      Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors,
                                             i, jj, k, element)
      Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj)
      # compute the contravariant sharp flux in the direction of the
      # averaged contravariant vector
      fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations)
      multiply_add_to_node_vars!(du, alpha * derivative_split[j, jj], fluxtilde2, equations, dg, i, j,  k, element)
      multiply_add_to_node_vars!(du, alpha * derivative_split[jj, j], fluxtilde2, equations, dg, i, jj, k, element)
    end

    # z direction
    for kk in (k+1):nnodes(dg)
      u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
      # pull the contravariant vectors and compute the average
      Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors,
                                             i, j, kk, element)
      Ja3_avg = 0.5 * (Ja3_node + Ja3_node_kk)
      # compute the contravariant sharp flux in the direction of the
      # averaged contravariant vector
      fluxtilde3 = volume_flux(u_node, u_node_kk, Ja3_avg, equations)
      multiply_add_to_node_vars!(du, alpha * derivative_split[k, kk], fluxtilde3, equations, dg, i, j, k,  element)
      multiply_add_to_node_vars!(du, alpha * derivative_split[kk, k], fluxtilde3, equations, dg, i, j, kk, element)
    end
  end
end

@inline function split_form_kernel!(du::AbstractArray{<:Any,5}, u,
                                    element, mesh::StructuredMesh{3},
                                    nonconservative_terms::Val{true}, equations,
                                    volume_flux, dg::DGSEM, cache, alpha=true)
  @unpack derivative_split = dg.basis
  @unpack contravariant_vectors = cache.elements
  symmetric_flux, nonconservative_flux = volume_flux

  # Apply the symmetric flux as usual
  split_form_kernel!(du, u, element, mesh, Val(false), equations, symmetric_flux, dg, cache, alpha)

  # Calculate the remaining volume terms using the nonsymmetric generalized flux
  for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, j, k, element)

    # pull the contravariant vectors in each coordinate direction
    Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
    Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
    Ja3_node = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)

    # The diagonal terms are zero since the diagonal of `derivative_split`
    # is zero. We ignore this for now.
    # In general, nonconservative fluxes can depend on both the contravariant
    # vectors (normal direction) at the current node and the averaged ones.
    # Thus, we need to pass both to the nonconservative flux.

    # x direction
    integral_contribution = zero(u_node)
    for ii in eachnode(dg)
      u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
      # pull the contravariant vectors and compute the average
      Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, k, element)
      Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii)
      # compute the contravariant nonconservative flux in the direction of the
      # averaged contravariant vector
      fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg, equations)
      integral_contribution = integral_contribution + derivative_split[i, ii] * fluxtilde1
    end

    # y direction
    for jj in eachnode(dg)
      u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
      # pull the contravariant vectors and compute the average
      Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, k, element)
      Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj)
      # compute the contravariant nonconservative flux in the direction of the
      # averaged contravariant vector
      fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg, equations)
      integral_contribution = integral_contribution + derivative_split[j, jj] * fluxtilde2
    end

    # z direction
    for kk in eachnode(dg)
      u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
      # pull the contravariant vectors and compute the average
      Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, i, j, kk, element)
      Ja3_avg = 0.5 * (Ja3_node + Ja3_node_kk)
      # compute the contravariant nonconservative flux in the direction of the
      # averaged contravariant vector
      fluxtilde3 = nonconservative_flux(u_node, u_node_kk, Ja3_node, Ja3_avg, equations)
      integral_contribution = integral_contribution + derivative_split[k, kk] * fluxtilde3
    end

    # The factor 0.5 cancels the factor 2 in the flux differencing form
    multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, dg, i, j, k, element)
  end
end

function calc_interface_flux!(cache, u, mesh::StructuredMesh{3},
                              nonconservative_terms, # can be Val{true}/Val{false}
                              equations, surface_integral, dg::DG)
  @unpack elements = cache

  @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2, 3)"

    # Interfaces in x-direction (`orientation` = 1)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[1, element],
                         element, 1, u, mesh,
                         nonconservative_terms, equations,
                         surface_integral, dg, cache)

    # Interfaces in y-direction (`orientation` = 2)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[2, element],
                         element, 2, u, mesh,
                         nonconservative_terms, equations,
                         surface_integral, dg, cache)

    # Interfaces in z-direction (`orientation` = 3)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[3, element],
                         element, 3, u, mesh,
                         nonconservative_terms, equations,
                         surface_integral, dg, cache)
  end

  return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
                                      orientation, u,
                                      mesh::StructuredMesh{3},
                                      nonconservative_terms::Val{false}, equations,
                                      surface_integral, dg::DG, cache)
  # This is slow for LSA, but for some reason faster for Euler (see #519)
  if left_element <= 0 # left_element = 0 at boundaries
    return surface_flux_values
  end

  @unpack surface_flux = surface_integral
  @unpack contravariant_vectors, inverse_jacobian = cache.elements

  right_direction = 2 * orientation
  left_direction = right_direction - 1

  for j in eachnode(dg), i in eachnode(dg)
    if orientation == 1
      u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, j, left_element)
      u_rr = get_node_vars(u, equations, dg, 1,          i, j, right_element)

      # If the mapping is orientation-reversing, the contravariant vectors' orientation
      # is reversed as well. The normal vector must be oriented in the direction
      # from `left_element` to `right_element`, or the numerical flux will be computed
      # incorrectly (downwind direction).
      sign_jacobian = sign(inverse_jacobian[1, i, j, right_element])

      # First contravariant vector Ja^1 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(1, contravariant_vectors,
                                                                  1, i, j, right_element)
    elseif orientation == 2
      u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), j, left_element)
      u_rr = get_node_vars(u, equations, dg, i, 1,          j, right_element)

      # See above
      sign_jacobian = sign(inverse_jacobian[i, 1, j, right_element])

      # Second contravariant vector Ja^2 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(2, contravariant_vectors,
                                                                  i, 1, j, right_element)
    else # orientation == 3
      u_ll = get_node_vars(u, equations, dg, i, j, nnodes(dg), left_element)
      u_rr = get_node_vars(u, equations, dg, i, j, 1,          right_element)

      # See above
      sign_jacobian = sign(inverse_jacobian[i, j, 1, right_element])

      # Third contravariant vector Ja^3 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(3, contravariant_vectors,
                                                                  i, j, 1, right_element)
    end

    # If the mapping is orientation-reversing, the normal vector will be reversed (see above).
    # However, the flux now has the wrong sign, since we need the physical flux in normal direction.
    flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations)

    for v in eachvariable(equations)
      surface_flux_values[v, i, j, right_direction, left_element] = flux[v]
      surface_flux_values[v, i, j, left_direction, right_element] = flux[v]
    end
  end

  return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
                                      orientation, u,
                                      mesh::StructuredMesh{3},
                                      nonconservative_terms::Val{true}, equations,
                                      surface_integral, dg::DG, cache)
  # See comment on `calc_interface_flux!` with `nonconservative_terms::Val{false}`
  if left_element <= 0 # left_element = 0 at boundaries
    return surface_flux_values
  end

  surface_flux, nonconservative_flux = surface_integral.surface_flux
  @unpack contravariant_vectors, inverse_jacobian = cache.elements

  right_direction = 2 * orientation
  left_direction = right_direction - 1

  for j in eachnode(dg), i in eachnode(dg)
    if orientation == 1
      u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, j, left_element)
      u_rr = get_node_vars(u, equations, dg, 1,          i, j, right_element)

      # If the mapping is orientation-reversing, the contravariant vectors' orientation
      # is reversed as well. The normal vector must be oriented in the direction
      # from `left_element` to `right_element`, or the numerical flux will be computed
      # incorrectly (downwind direction).
      sign_jacobian = sign(inverse_jacobian[1, i, j, right_element])

      # First contravariant vector Ja^1 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(1, contravariant_vectors,
                                                                  1, i, j, right_element)
    elseif orientation == 2
      u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), j, left_element)
      u_rr = get_node_vars(u, equations, dg, i, 1,          j, right_element)

      # See above
      sign_jacobian = sign(inverse_jacobian[i, 1, j, right_element])

      # Second contravariant vector Ja^2 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(2, contravariant_vectors,
                                                                  i, 1, j, right_element)
    else # orientation == 3
      u_ll = get_node_vars(u, equations, dg, i, j, nnodes(dg), left_element)
      u_rr = get_node_vars(u, equations, dg, i, j, 1,          right_element)

      # See above
      sign_jacobian = sign(inverse_jacobian[i, j, 1, right_element])

      # Third contravariant vector Ja^3 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(3, contravariant_vectors,
                                                                  i, j, 1, right_element)
    end

    # If the mapping is orientation-reversing, the normal vector will be reversed (see above).
    # However, the flux now has the wrong sign, since we need the physical flux in normal direction.
    flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations)

    # Compute both nonconservative fluxes
    # In general, nonconservative fluxes can depend on both the contravariant
    # vectors (normal direction) at the current node and the averaged ones.
    # However, both are the same at watertight interfaces, so we pass the
    # `normal_direction` twice.
    # Scale with sign_jacobian to ensure that the normal_direction matches that
    # from the flux above
    noncons_left  = sign_jacobian * nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations)
    noncons_right = sign_jacobian * nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations)

    for v in eachvariable(equations)
      # Note the factor 0.5 necessary for the nonconservative fluxes based on
      # the interpretation of global SBP operators coupled discontinuously via
      # central fluxes/SATs
      surface_flux_values[v, i, j, right_direction, left_element] = flux[v] + 0.5 * noncons_left[v]
      surface_flux_values[v, i, j, left_direction, right_element] = flux[v] + 0.5 * noncons_right[v]
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::StructuredMesh{3}, equations, surface_integral, dg::DG)
  @assert isperiodic(mesh)
end

function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple,
                             mesh::StructuredMesh{3}, equations, surface_integral, dg::DG)
  @unpack surface_flux_values = cache.elements
  linear_indices = LinearIndices(size(mesh))

  for cell_z in axes(mesh, 3), cell_y in axes(mesh, 2)
    # Negative x-direction
    direction = 1
    element = linear_indices[begin, cell_y, cell_z]

    for k in eachnode(dg), j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (1, j, k), (j, k), element)
    end

    # Positive x-direction
    direction = 2
    element = linear_indices[end, cell_y, cell_z]

    for k in eachnode(dg), j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (nnodes(dg), j, k), (j, k), element)
    end
  end

  for cell_z in axes(mesh, 3), cell_x in axes(mesh, 1)
    # Negative y-direction
    direction = 3
    element = linear_indices[cell_x, begin, cell_z]

    for k in eachnode(dg), i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (i, 1, k), (i, k), element)
    end

    # Positive y-direction
    direction = 4
    element = linear_indices[cell_x, end, cell_z]

    for k in eachnode(dg), i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (i, nnodes(dg), k), (i, k), element)
    end
  end

  for cell_y in axes(mesh, 2), cell_x in axes(mesh, 1)
    # Negative z-direction
    direction = 5
    element = linear_indices[cell_x, cell_y, begin]

    for j in eachnode(dg), i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 3,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (i, j, 1), (i, j), element)
    end

    # Positive z-direction
    direction = 6
    element = linear_indices[cell_x, cell_y, end]

    for j in eachnode(dg), i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 3,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (i, j, nnodes(dg)), (i, j), element)
    end
  end
end

function apply_jacobian!(du,
                         mesh::Union{StructuredMesh{3}, P4estMesh{3}},
                         equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      factor = -cache.elements.inverse_jacobian[i, j, k, element]

      for v in eachvariable(equations)
        du[v, i, j, k, element] *= factor
      end
    end
  end

  return nothing
end

end # @muladd


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_OnOffRenderer" begin
     sphere = Sphere(translation(Vec(2, 0, 0)) * scaling(Vec(0.2, 0.2, 0.2)),
                    Material(DiffuseBRDF(UniformPigment(WHITE))))
     image = HDRimage(3, 3)
     camera = OrthogonalCamera()
     tracer = ImageTracer(image, camera, 0) # without antialiasing
     world = World()
     add_shape!(world, sphere)
     renderer = OnOffRenderer(world)
     fire_all_rays!(tracer, renderer)

     @test Raytracing.get_pixel(image, 0, 0) ≈ BLACK
     @test Raytracing.get_pixel(image, 1, 0) ≈ BLACK
     @test Raytracing.get_pixel(image, 2, 0) ≈ BLACK

     @test Raytracing.get_pixel(image, 0, 1) ≈ BLACK
     @test Raytracing.get_pixel(image, 1, 1) ≈ WHITE # gives error with antialiasing
     @test Raytracing.get_pixel(image, 2, 1) ≈ BLACK

     @test Raytracing.get_pixel(image, 0, 2) ≈ BLACK
     @test Raytracing.get_pixel(image, 1, 2) ≈ BLACK
     @test Raytracing.get_pixel(image, 2, 2) ≈ BLACK
end

@testset "test_FlatRenderer" begin
     sphere_color = RGB{Float32}(5.0, 6.0, 7.0)
     α = 0.2
     sphere = Sphere(translation(Vec(2, 0, 0)) * scaling(Vec(α, α, α)),
                    Material(DiffuseBRDF(UniformPigment(sphere_color))))

     @test sphere.Material.brdf.pigment == UniformPigment(sphere_color)

     image = HDRimage(3, 3)
     camera = OrthogonalCamera()
     tracer = ImageTracer(image, camera, 0) # without antialiasing

     world = World()
     add_shape!(world, sphere)
     renderer = FlatRenderer(world)
     fire_all_rays!(tracer, renderer)

     @test Raytracing.get_pixel(image, 0, 0) ≈ BLACK
     @test Raytracing.get_pixel(image, 1, 0) ≈ BLACK
     @test Raytracing.get_pixel(image, 2, 0) ≈ BLACK

     @test Raytracing.get_pixel(image, 0, 1) ≈ BLACK
     @test Raytracing.get_pixel(image, 1, 1) ≈ sphere_color # gives error with antialiasing
     @test Raytracing.get_pixel(image, 2, 1) ≈ BLACK

     @test Raytracing.get_pixel(image, 0, 2) ≈ BLACK
     @test Raytracing.get_pixel(image, 1, 2) ≈ BLACK
     @test Raytracing.get_pixel(image, 2, 2) ≈ BLACK
end

@testset "test_PathTracer" begin
     # Here is impemented the Furnace test.
     # It is runned several times using random values
     # for the emitted radiance and reflectance.

     pcg = PCG()

     for i in 1:10
          world = World()
          emitted_radiance = random(pcg)
          reflectance = random(pcg)
          enclosure_material =
               Material(
                    DiffuseBRDF(UniformPigment(RGB{Float32}(1., 1., 1.) * reflectance)),
                    UniformPigment(RGB{Float32}(1., 1., 1.) * emitted_radiance),
               )

          add_shape!(world, Sphere(enclosure_material))

          path_tracer = PathTracer(world, BLACK, pcg, 1, 100, 101)

          ray = Ray(Point(0., 0., 0.), Vec(1., 0., 0.))
          color = path_tracer(ray)

          expected = emitted_radiance / (1.0 - reflectance)

          err=1e-3
          @test are_close(expected, color.r, err)
          @test are_close(expected, color.g, err)
          @test are_close(expected, color.b, err)

     end
end

@testset "test_PointLightTracer" begin
     @testset "quick_ray_intersection_sphere" begin
          world = World()

          sphere1 = Sphere(translation(Vec(2.0, 0.0, 0.0)))
          sphere2 = Sphere(translation(Vec(8.0, 0.0, 0.0)))
          add_shape!(world, sphere1)
          add_shape!(world, sphere2)

          @test is_point_visible(world, Point(-10.0, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test !is_point_visible(world, Point(10.0, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test !is_point_visible(world, Point(5.0, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test is_point_visible(world, Point(5.0, 0.0, 0.0), Point(4.0, 0.0, 0.0))
          @test is_point_visible(world, Point(0.5, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test is_point_visible(world, Point(0.0, 10.0, 0.0), Point(0.0, 0.0, 0.0))
          @test is_point_visible(world, Point(0.0, 0.0, 10.0), Point(0.0, 0.0, 0.0))
     end

     @testset "quick_ray_intersection_plane" begin
          world = World()

          plane1 =
               Plane(
                    translation(Vec(2.0, 0.0, 0.0)) * rotation_y(π/2.)
               )
          plane2 =
               Plane(
                    translation(Vec(8.0, 0.0, 0.0)) * rotation_y(π/2.)
               )
          add_shape!(world, plane1)
          add_shape!(world, plane2)

          @test is_point_visible(world, Point(-10.0, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test !is_point_visible(world, Point(10.0, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test !is_point_visible(world, Point(5.0, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test is_point_visible(world, Point(5.0, 0.0, 0.0), Point(4.0, 0.0, 0.0))
          @test is_point_visible(world, Point(0.5, 0.0, 0.0), Point(0.0, 0.0, 0.0))
          @test is_point_visible(world, Point(0.0, 10.0, 0.0), Point(0.0, 0.0, 0.0))
          @test is_point_visible(world, Point(0.0, 0.0, 10.0), Point(0.0, 0.0, 0.0))
     end
end


text:
# ---
# title: 1404. Number of Steps to Reduce a Number in Binary Representation to One
# id: problem1404
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Bit Manipulation
# link: <https://leetcode.com/problems/number-of-steps-to-reduce-a-number-in-binary-representation-to-one/description/>
# hidden: true
# ---
#
# Given a number `s` in their binary representation. Return the number of steps
# to reduce it to 1 under the following rules:
#
#   * If the current number is even, you have to divide it by 2.
#
#   * If the current number is odd, you have to add 1 to it.
#
# It's guaranteed that you can always reach to one for all testcases.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "1101"
#     Output: 6
#     Explanation: "1101" corressponds to number 13 in their decimal representation.
#     Step 1) 13 is odd, add 1 and obtain 14.
#     Step 2) 14 is even, divide by 2 and obtain 7.
#     Step 3) 7 is odd, add 1 and obtain 8.
#     Step 4) 8 is even, divide by 2 and obtain 4.
#     Step 5) 4 is even, divide by 2 and obtain 2.
#     Step 6) 2 is even, divide by 2 and obtain 1.
#
#
# **Example 2:**
#
#
#
#     Input: s = "10"
#     Output: 1
#     Explanation: "10" corressponds to number 2 in their decimal representation.
#     Step 1) 2 is even, divide by 2 and obtain 1.
#
#
# **Example 3:**
#
#
#
#     Input: s = "1"
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 500`
#   * `s` consists of characters '0' or '1'
#   * `s[0] == '1'`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# High-level interface

"""
    throttle(callback, Event, freq = 1)
    throttle(callback, Event, seconds = 1)

Throttle `Event` type for `callback` so that it is triggered either only every
`freq`'th time  or every `seconds` seconds.

## Examples

If you want to only sporadically log metrics ([`LogMetrics`](#)) or images
([`LogVisualization`](#)), `throttle` can be used as follows.

Every 10 steps:

```julia
callback = throttle(LogMetrics(TensorBoardBackend()), StepEnd, freq = 10)
learner = Learner(<args>, callback)
```

Or every 5 seconds:

```julia
callback = throttle(LogMetrics(TensorBoardBackend()), StepEnd, seconds = 5)
learner = Learner(<args>, callback)
```
"""
function throttle(callback, event::Type{<:Event}; freq = nothing, seconds = nothing)
    xor(isnothing(freq), isnothing(seconds)) || error("Pass either `every` OR `seconds`.")
    if !isnothing(freq)
        return ConditionalCallback(callback, FrequencyThrottle(freq, event))
    elseif !isnothing(seconds)
        return ConditionalCallback(callback, TimeThrottle(seconds, event))
    end
end

# Implementation

"""
    abstract type CallbackCondition

Supertype for conditions to use with [`ConditionalCallback`](#).
To implement a `CallbackCondition`, implement
[`shouldrun`](#)`(::MyCondition, event, phase)`.

See [`FrequencyThrottle`](#), [`TimeThrottle`](#) and [`throttle`](#).
"""
abstract type CallbackCondition end

"""
    ConditionalCallback(callback, condition) <: Callback

Wrapper callback that only forwards events to the wrapped callback
if [`CallbackCondition`](#) `condition` is met. See [`throttle`](#).
"""
struct ConditionalCallback <: Callback
    callback::Callback
    condition::CallbackCondition
end

stateaccess(cc::ConditionalCallback) = stateaccess(cc.callback)
init!(cc::ConditionalCallback, learner) = init!(cc.callback, learner)

function on(event::Event, phase::Phase, cb::ConditionalCallback, learner)
    if shouldrun(cb.condition, event, phase)
        on(event, phase, cb.callback, learner)
    end
end

mutable struct FrequencyThrottle <: CallbackCondition
    freq
    event
    counter
    FrequencyThrottle(f, e) = new(f, e, 1)
end

function shouldrun(c::FrequencyThrottle, event, phase)
    if typeof(event) == c.event
        if c.counter == 1
            c.counter = c.freq
            return true
        else
            c.counter -= 1
            return false
        end
    else
        return true
    end
end

mutable struct TimeThrottle <: CallbackCondition
    seconds
    event
    timer
    TimeThrottle(s, e, t = nothing) = new(s, e, t)
end

function shouldrun(c::TimeThrottle, event, phase)
    if typeof(event) == c.event
        if isnothing(c.timer) || (time() - c.timer) > c.seconds
            c.timer = time()
            return true
        else
            return false
        end
    else
        return true
    end
end

@testset "throttle" begin
    function train(cb)
        learner = testlearner(cb)
        epoch!(learner, TrainingPhase())
        return learner.cbstate.history[TrainingPhase()]
    end

    @test train(Recorder()).steps == 16
    @testset "freq" begin
        @test train(throttle(Recorder(), StepEnd, freq = 2)).steps == 8
        @test train(throttle(Recorder(), StepEnd, freq = 16)).steps == 1
    end

    if !Base.Sys.iswindows()
        @testset "seconds" begin
            @test train(throttle(Recorder(), StepEnd, seconds = 0)).steps == 16
            @test train(throttle(Recorder(), StepEnd, seconds = 10)).steps == 1
        end
    end
end


text:
import StaticArrays
import Polyhedra

struct Library <: Polyhedra.Library
    solver::Polyhedra.SolverOrNot
end
Library() = Library(nothing)
_isone(d::Int) = isone(d)
_isone(d::StaticArrays.Size{(1,)}) = true
_isone(d::StaticArrays.Size) where T = false
Polyhedra.similar_library(l::Library, d::Polyhedra.FullDim, T::Type) = Polyhedra.default_library(d, T)
function Polyhedra.similar_library(l::Library, d::Polyhedra.FullDim, ::Type{Float64})
    if _isone(d)
        return Polyhedra.default_library(d, Float64)
    else
        return Library(l.solver)
    end
end

mutable struct Polyhedron <: Polyhedra.Polyhedron{Float64}
    ine::Union{Nothing, Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}}
    ext::Union{Nothing, Polyhedra.MixedMatVRep{Float64, Matrix{Float64}}}
    nohredundancy::Bool
    novredundancy::Bool
    solver::Polyhedra.SolverOrNot
    area::Union{Nothing, Float64}
    volume::Union{Nothing, Float64}

    function Polyhedron(ine, ext, nri::Bool, nrg::Bool, solver)
        new(ine, ext, nri, nrg, solver, nothing, nothing)
    end
    function Polyhedron(ine::Polyhedra.HRepresentation{Float64}, solver)
        new(ine, nothing, false, false, solver, nothing, nothing)
    end
    function Polyhedron(ext::Polyhedra.VRepresentation{Float64}, solver)
        new(nothing, ext, false, false, solver, nothing, nothing)
    end
end

Polyhedra.FullDim(p::Polyhedron) = Polyhedra.FullDim_rep(p.ine, p.ext)
Polyhedra.library(p::Polyhedron) = Library(p.solver)
Polyhedra.similar_type(::Type{<:Polyhedron}, d::Polyhedra.FullDim, T::Type) = Polyhedra.default_type(d, T)
function Polyhedra.similar_type(::Type{<:Polyhedron}, d::Polyhedra.FullDim, ::Type{Float64})
    if _isone(d)
        return Polyhedra.default_type(d, Float64)
    else
        Polyhedron
    end
end

Polyhedra.default_solver(p::Polyhedron; T=nothing) = p.solver
Polyhedra.supportssolver(::Type{<:Polyhedron}) = true

Polyhedra.hvectortype(::Union{Polyhedron, Type{Polyhedron}}) = Polyhedra.hvectortype(Polyhedra.MixedMatHRep{Float64, Matrix{Float64}})
Polyhedra.vvectortype(::Union{Polyhedron, Type{Polyhedron}}) = Polyhedra.vvectortype(Polyhedra.MixedMatVRep{Float64, Matrix{Float64}})

# Helpers
epsz = 1e-8

function qhull(p::Polyhedron, rep=:Auto)
    if rep == :V || (rep == :Auto && (p.ext !== nothing))
        p.ext, ine, p.area, p.volume = qhull(getext(p), p.solver)
        if p.ine === nothing
            # Otherwise, it is not interesting as it may have redundancy
            p.ine = ine
        end
    else
        @assert rep == :H || rep == :Auto
        p.ine, ext = qhull(getine(p), p.solver)
        if p.ext === nothing
            p.ext = ext
        end
    end
    p.nohredundancy = true
    p.novredundancy = true
end

function qhull(h::Polyhedra.MixedMatVRep{T}, solver=nothing) where T
    if Polyhedra.hasrays(h)
        error("Rays are not supported.")
    end
    V = h.V
    ch = chull(V)
    V = ch.points[ch.vertices, :]
    vnored = Polyhedra.vrep(V)
    N = Polyhedra.fulldim(h)
    A = ch.facets[:, 1:N]
    b = -ch.facets[:, N+1]
    h = Polyhedra.hrep(A, b)
    # The same facet may be listed several times in case the facet is not
    # simplicial, e.g. a quad facet is listed twice as it is made of two
    # triangles. After `removeduplicates`, there should be no redundancy.
    hnored = Polyhedra.removeduplicates(h)
    vnored, hnored, ch.area, ch.volume
end

function qhull(h::Polyhedra.MixedMatHRep{T}, solver=nothing) where {T<:Real}
    linset = h.linset
    if !isempty(linset)
        error("Equalities are not supported.")
    end
    N = Polyhedra.fulldim(h)
    containorigin = Polyhedra.ininterior(zeros(N), h)
    if !containorigin
        chebycenter, chebyradius = Polyhedra.chebyshevcenter(h, solver)
        h = Polyhedra.translate(h, -chebycenter)
    end

    A = h.A
    b = h.b
    m = size(A, 1)
    B = Matrix{T}(undef, m, N)
    for i in 1:m
        @assert !(i in linset)
        if b[i] < epsz
            error("The origin should be in the interior of the polytope but the $(i)th inequality is not safisfied at the origin.")
        end
        B[i,:] = (@view A[i,:]) / b[i]
    end
    ch = chull(B)
    C = ch.points[ch.vertices, :]
    hnored = Polyhedra.hrep(C, ones(size(C, 1)))
    Vlift = ch.facets
    nvreps = size(Vlift, 1)
    irays = BitSet()
    ipoints = BitSet()
    for i in 1:nvreps
        if Vlift[i, N+1] > epsz
            error("The origin should be in the interior of the convex hull")
        end
        if Vlift[i, N+1] > -epsz
            push!(irays, i)
        else
            push!(ipoints, i)
        end
    end
    rays = collect(irays)
    points = collect(ipoints)
    R = Matrix{T}(undef, length(rays), N)
    V = Matrix{T}(undef, length(points), N)
    nr = nv = 0
    for i in 1:nvreps
        if i in irays
            nr += 1
            R[nr, :] = -@view Vlift[i, 1:N]
        else
            nv += 1
            V[nv, :] = -(@view Vlift[i, 1:N]) / Vlift[i, N+1]
        end
    end
    v = Polyhedra.vrep(V, R)
    if !containorigin
        hnored = Polyhedra.translate(hnored, chebycenter)
    end
    if !containorigin
        v = Polyhedra.translate(v, chebycenter)
    end
    # See comment in `qhull(::MixedMatVRep)`.
    vnored = Polyhedra.removeduplicates(v)
    # `ch.area`, `ch.volume` give area and volume of the polar of the shifted
    # polytope, not the area and volume of the polytope with reps `h` and `v`.
    hnored, vnored
end

function getine(p::Polyhedron)
    if p.ine === nothing
        qhull(p)
    end
    p.ine
end
function getext(p::Polyhedron)
    if p.ext === nothing
        qhull(p)
    end
    p.ext
end

function clearfield!(p::Polyhedron)
    p.ine = nothing
    p.ext = nothing
    p.nohredundancy = false
    p.novredundancy = false
end

# Implementation of Polyhedron's mandatory interface
function Polyhedra.polyhedron(repit::Polyhedra.Representation, lib::Library)
    Polyhedron(repit, lib.solver)
end

Polyhedron(h::Polyhedra.HRepresentation, solver) = Polyhedron(convert(Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}, h), solver)
Polyhedron(v::Polyhedra.VRepresentation, solver) = Polyhedron(convert(Polyhedra.MixedMatVRep{Float64, Matrix{Float64}}, v), solver)

Polyhedron(d::Polyhedra.FullDim, hps::Polyhedra.HyperPlaneIt, hss::Polyhedra.HalfSpaceIt; solver=nothing) = Polyhedron(Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}(d, hps, hss), solver)
Polyhedron(d::Polyhedra.FullDim, ps::Polyhedra.PointIt, ls::Polyhedra.LineIt, rs::Polyhedra.RayIt; solver=nothing) = Polyhedron(Polyhedra.MixedMatVRep{Float64, Matrix{Float64}}(d, ps, ls, rs), solver)

function Base.copy(p::Polyhedron)
    ine = nothing
    if p.ine !== nothing
        ine = copy(p.ine)
    end
    ext = nothing
    if p.ext !== nothing
        ext = copy(p.ext)
    end
    Polyhedron(ine, ext, p.nohredundancy, p.novredundancy, p.solver)
end
function Polyhedra.hrepiscomputed(p::Polyhedron)
    p.ine !== nothing
end
function Polyhedra.hrep(p::Polyhedron)
    getine(p)
end
function Polyhedra.vrepiscomputed(p::Polyhedron)
    p.ext !== nothing
end
function Polyhedra.vrep(p::Polyhedron)
    getext(p)
end
function Polyhedra.sethrep!(p::Polyhedron, h::Polyhedra.HRepresentation)
    p.ine = h
end
function Polyhedra.setvrep!(p::Polyhedron, v::Polyhedra.VRepresentation)
    p.ext = v
end
function Polyhedra.resethrep!(p::Polyhedron, h::Polyhedra.HRepresentation)
    clearfield!(p)
    p.ine = h
end
function Polyhedra.resetvrep!(p::Polyhedron, v::Polyhedra.VRepresentation)
    clearfield!(p)
    p.ext = v
end
function Polyhedra.removehredundancy!(p::Polyhedron)
    if !p.nohredundancy
        qhull(p, :H)
    end
end
function Polyhedra.removevredundancy!(p::Polyhedron)
    if !p.novredundancy
        qhull(p, :V)
    end
end

function Polyhedra.volume(p::Polyhedron)
    if p.volume === nothing
        qhull(p, :V)
    end
    return p.volume
end


text:
"""Euclidean norm of segment"""
length{D, T}(line::Segment{D, T}) = norm(line.B - line.A)

"""Iterator over edges of the curve as Segment entities"""
function iter_segments{D, T}(curve::Curve{D, T})
    (Segment(curve.points[i0], curve.points[i1]) for (i0, i1) in curve.segments)
end

"""Total length of the curve"""
length{D, T}(curve::Curve{D, T}) = sum(imap(length, iter_segments(curve)))

"""Smallest segment of the curve"""
min_segment{D, T}(curve::Curve{D, T}) = minimum(imap(length, iter_segments(curve)))

"""Longest segment of the curve"""
max_segment{D, T}(curve::Curve{D, T}) = maximum(imap(length, iter_segments(curve)))

"""Length of branch refered to by segments indices in curve.segments"""
function branch_length(curve::Curve, segments::Vector{Int})
    l = 0.
    for segment in segments
        # index -> indices
        v0, v1 = curve.segments[segment]
        # Vertices as physical points
        l += norm(curve.points[v1]-curve.points[v0])
    end
    l
end

"""Produce an array of branch lengths in the order they appear in bmap"""
function branch_lengths(curve::Curve,  bmap::Dict{Tuple{Int, Int}, Vector{Vector{Int}}})
    lengths = Vector{Real}(num_branches(curve))
    for key in keys(bmap)
        for branch in bmap[key]
            push!(lengths, branch_length(branch))
        end
    end
    lengths
end

"""Length of the longest branch"""
function min_branch(curve::Curve,  bmap::Dict{Tuple{Int, Int}, Vector{Vector{Int}}})
    minimum(branch_lengths(curve, bmap))
end

"""Length of the shortest branch"""
function max_branch(curve::Curve,  bmap::Dict{Tuple{Int, Int}, Vector{Vector{Int}}})
    maximum(branch_lengths(curve, bmap))
end

# Generate those which compute the map
for foo in (:branch_lengths, :min_branch, :max_branch)
    @eval $(foo)(curve::Curve) = $(foo)(curve, branch_lengths(curve))
end

######################################################################
# Box counting and fractional dimension
######################################################################
include("fractal_dim.jl")

"""
Scale the curve (in place) such as the distances between points are
changed by factor a.
"""
function scale!{D, T}(c::Curve{D, T}, a::Number)
    o = BoundingBox(c).low

    points = c.points
    for i in 1:length(c.points)
        points[i] = a*(points[i] - o) + o
    end
end

"""
The distance vector between points is scaled by diagm(a). This allows for
sheer.
"""
function scale!{D, T}(c::Curve{D, T}, a::NTuple{D, T})
    o = BoundingBox(c).low

    X = diagm(SVector(a...))
    points = c.points
    for i in 1:length(c.points)
        points[i] = X*(points[i] - o) + o
    end
end

"""Translate! the curve by vector b"""
function translate!{D, T}(c::Curve{D, T}, b::SVector{D, T})
    points = c.points
    for i in 1:length(c.points)
        points[i] += b
    end
end

"""Scale! curve to rectangle b, i.e. the bbox of c should be b."""
function scale!{D, T}(c::Curve{D, T}, b::BoundingBox{D, T})
    target_low, target_high = b.low, b.high

    cframe = BoundingBox(c)
    c_low, c_high = cframe.low, cframe.high
    # Figure out the relative scaling
    scale = (target_high - target_low)./(c_high - c_low)

    X = diagm(scale)
    points = c.points
    for i in 1:length(c.points)
        # Scale and shift to next origin
        points[i] = X*(points[i] - c_low) + target_low
    end
end

# IDEA: fd in parallel
# TODO: tortuosity


text:
# This file was generated, do not modify it. # hide
DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree
forest = EnsembleModel(model=DecisionTreeRegressor())

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
gen_scenario(length::Int) = Float64[2 + sin(2π * t / 52) + rand() for t in 1:length]
gen_dates(length::Int, year::Int) = DateTime[DateTime(year, October, 1) + Week(i) for i in 0:length - 1]
gen_ts(length::Int, year::Int) = TimeArray(gen_dates(length, year), gen_scenario(length), ["Discharge (m³/s)"], Week(1))
gen_tses(length::Int, year::Int, n::Int) = [gen_ts(length, year + i - 1) for i in 1:n]
gen_river(name::AbstractString, length::Int, year::Int, n::Int) =
  NaturalRiver(name=name, scenarios=gen_tses(length, year, n))
gen_river(name::AbstractString, length::Int, year::Int, n::Int, envFlow::Discharge, maxFlow::Discharge) =
  DivertedRiver(name=name, scenarios=gen_tses(length, year, n), environmental_flow=envFlow, maximum_flow=maxFlow)

function compare_rivers(r1::NaturalRiver, r2::NaturalRiver)
  @fact getName(r1) --> getName(r2)
  @fact countScenarios(r1) --> countScenarios(r2)
  @fact countTimeSteps(r1) --> countTimeSteps(r2)
  @fact getPeriod(r1) --> getPeriod(r2)

  for s in 1:countScenarios(r1)
    @fact timestamp(getScenarioTS(r1, s)) --> timestamp(getScenarioTS(r2, s))
    @fact getScenario(r1, s) --> getScenario(r2, s)
    @fact colnames(getScenarioTS(r1, s)) --> colnames(getScenarioTS(r2, s))
    @fact meta(getScenarioTS(r1, s)) --> meta(getScenarioTS(r2, s))
  end
end


text:
using Test
using StaticArrays
using RobotDynamics: DynamicsJacobian

n,m = 4,1
a = SizedMatrix{4,5}(rand(n,n+m))
F = DynamicsJacobian(a)
F[1,1] = 2
@test F[1,1] == 2
@test F[1,2] == a[1,2]
@test F.data === a

@test size(F.A) == (n,n)
@test size(F.B) == (n,m)
@test F.A == a[:,1:n]
@test F.B == a[:,n .+ (1:m)]

F.A[1,2] = 3
@test F[1,2] == 3
F.B[3,1] = 4
@test F[3,n+1] == 4

@test RobotDynamics.get_A(F) ≈ F.A
@test RobotDynamics.get_B(F) ≈ F.B
@test RobotDynamics.get_static_A(F) ≈ F.A
@test RobotDynamics.get_static_B(F) ≈ F.B
@test RobotDynamics.get_static_A(F) isa SMatrix{n,n}
@test RobotDynamics.get_static_B(F) isa SMatrix{n,m}

F_ = SMatrix(F)
@test DynamicsJacobian{n,n+m}(Tuple(F_)) isa DynamicsJacobian{n,n+m,Float64}
@test DynamicsJacobian{n,n+m,Float32}(Tuple(F_)) isa DynamicsJacobian{n,n+m,Float32}
@test Tuple(F) == Tuple(F_)
@test DynamicsJacobian(F_) isa DynamicsJacobian
@test DynamicsJacobian(F_) ≈ F

@test RobotDynamics.get_data(F_) === F_
@test RobotDynamics.get_data(F) === F.data.data
@test RobotDynamics.get_data(F.data) === F.data.data
@test RobotDynamics.get_data(F) isa Matrix
@test RobotDynamics.get_data(F.data) isa Matrix

## Test automatic diff methods
model = Cartpole()
RD.diffmethod(::Cartpole) = RD.ForwardAD()  # reset it
@test RD.diffmethod(model) == RD.ForwardAD()
F1 = DynamicsJacobian(model)
F2 = copy(F1)
@test F2 isa DynamicsJacobian
x,u = rand(model)
z = KnotPoint(x,u,0.1)

jacobian!(F1, model, z)
RD._jacobian!(RD.FiniteDifference(), F2, model, z)
err = norm(F1-F2)
@test 1e-12 < err < 1e-6   # shouldn't be the same, but close

# pass in the cache
cache = FiniteDiff.JacobianCache(model, Val(:central), Float32)
@test eltype(cache.x1) == Float32
cache = FiniteDiff.JacobianCache(model, Val(:central))
RD._jacobian!(RD.FiniteDifference(), F2, model, z, cache)
@test eps(Float64) < norm(F1-F2) < err
@test_nowarn jacobian!(F2, model, z, cache)
@test_nowarn RD._jacobian!(RD.ForwardAD(), F2, model, z, cache)

# test discrete jacobian
discrete_jacobian!(RK4, F1, model, z)
@test (@allocated discrete_jacobian!(RK4, F1, model, z)) == 0
RD._discrete_jacobian!(RD.FiniteDifference(), RK4, F2, model, z)
err = norm(F1-F2)
@test eps(Float64) < err < 1e-6

# test finite differencing
RD.diffmethod(::Cartpole) = RD.FiniteDifference()
discrete_jacobian!(RK4, F2, model, z, cache)
@test eps(Float64) < norm(F1-F2) < err

# check that finite diff doesn't allocate
cache = FiniteDiff.JacobianCache(model)
discrete_jacobian!(RK4, F2, model, z, cache)
@test (@allocated discrete_jacobian!(RK4, F2, model, z, cache)) == 0  # it's 0 for b

RD.diffmethod(::Cartpole) = RD.ForwardAD()  # reset it

# Sparse Jacobians
using Rotations
model = Quadrotor()
# model = Cartpole()
RD.diffmethod(::Quadrotor) = RD.FiniteDifference()
cache0 = RD.gen_cache(model)
@test cache0 isa FiniteDiff.JacobianCache
@test cache0.colorvec == 1:sum(size(model))
F1 = RD.DynamicsJacobian(model)
x,u = rand(model)
# @test norm(x[4:7]) ≈ 1
z = StaticKnotPoint(x,u,0.1)
discrete_jacobian!(RK4, F1, model, z, cache0)

# Detect sparsity and generated a colored cache
sparsity = RD.detect_sparsity(RK4, model)
cache_sparse = FiniteDiff.JacobianCache(model, colored=true, sparsity=sparsity)
@test maximum(cache_sparse.colorvec) < sum(size(model))
F2 = copy(F1)
discrete_jacobian!(RK4, F2, model, z, cache_sparse)
@test norm(F1 - F2) ≈ 0

############################################################################################
#                                    JACOBIAN-VECTOR PRODUCTS
############################################################################################
model = Cartpole()
λ = @SVector rand(state_dim(model))
F = RD.DynamicsJacobian(model)
x,u = rand(model)
z = StaticKnotPoint(x,u,0.1)
discrete_jacobian!(RK4, F, model, z)
dg0 = F'λ
jacobian!(F, model, z)
g0 = F'λ

RD.diffmethod(::Cartpole) = RD.ForwardAD()
g = zeros(sum(size(model)))
RD.discrete_jvp!(RK4, g, model, z, λ)
@test g ≈ dg0
@test (@allocated RD.discrete_jvp!(RK4, g, model, z, λ)) == 0

RD.jvp!(g, model, z, λ)
@test g ≈ g0
@test (@allocated RD.jvp!(g, model, z, λ)) == 0

# Finite difference
RD.diffmethod(::Cartpole) = RD.FiniteDifference()
cache = RD.gen_grad_cache(model)
@test cache isa FiniteDiff.GradientCache
g .= 0
RD.discrete_jvp!(RK4, g, model, z, λ, cache)
@test g ≈ dg0 atol=1e-6
@test (@allocated RD.discrete_jvp!(RK4, g, model, z, λ, cache)) == 0

RD.jvp!(g, model, z, λ, cache)
@test g ≈ g0 atol=1e-6
@test (@allocated RD.jvp!(g, model, z, λ, cache)) == 0

############################################################################################
#                               DYNAMICS HESSIANS
############################################################################################
n,m = size(model)
H = zeros(n+m,n+m)
RD.∇jacobian!(H, model, z, λ)
RD.∇discrete_jacobian!(RK4, H, model, z, λ)


text:
module BoardGameGeek

using HTTP
using EzXML
using Dates: Date
using ProgressMeter: @showprogress
import Base: show

# Using the BGG XML API2
# Reference: https://boardgamegeek.com/wiki/page/BGG_XML_API2
const XMLAPI2 = "https://boardgamegeek.com/xmlapi2"

include("compat.jl")
include("utils.jl")
include("game.jl")
include("user.jl")

export BGGUser
export userinfo, userreviews, buddies, collection
export gameinfo, gamereviews

end # module


text:
"""
    convert_WOAunits_to_unitful(s)

Converts the string from WOA's `units` attribute to a Unitful.jl unit.
"""
convert_to_Unitful(v::String) = @match v begin
    "micromoles_per_liter"             => u"μmol/l"
    "micromoles_per_kilogram"          => u"μmol/kg"
    "percent"                          => u"percent"
    "meters"                           => u"m"
    "degrees_north"                    => u"°"
    "degrees_east"                     => u"°"
    "degrees_celsius"                  => u"°C"
    _ => nothing
end


text:
include(Pkg.dir()*"/Processing/src/Processing2D.jl")

Processing2D.coordSystem(-1, 1, 1, -1) # user-defined coordinate system (x-axis min = -1,
						  # x-axis max = 1, y-axis min = -1, y-axis max = 1,
						  # xcent = 0, ycent = 0)
						  #
						  # I can just think easier this way...
						  #
						  # default coordinate system is in units of pixels with
						  # top-left corner of drawing surface labeled (0,0),
						  # just like Processing's default coordinate system
Processing2D.smooth()

Processing2D.ellipse(0,0,0.3,0.5)
Processing2D.ellipse(-.5,.5,0.2,0.2)
Processing2D.triangle(.3, .75, .58, .20, .86, .75)
Processing2D.triangle(.3, .6, .58, .40, .86, .7)
Processing2D.point(0, 0)
Processing2D.fill(1,0,0,1)
Processing2D.quad(-.3, -.75, -.58, -.20, -.86, -.75, -.2, -.4)
Processing2D.noFill()
Processing2D.rect(-.6, -.6, .2, .5)
Processing2D.line(.5, -.4, .7, -.5)
Processing2D.line(.7, -.4, .5, -.5)
Processing2D.strokeWeight(1)
Processing2D.stroke(0.9, 0, 0, 1)
Processing2D.arc(-.6, .6, .3, .3, 0.4pi, pi, CHORD)
Processing2D.animate()

x = 0
while true
	Processing2D.fill(sin(x)/2+0.5,cos(x)/2+0.5,(sin(x)*cos(x))/2+0.5,1)
	Processing2D.ellipse(0,0,0.2,0.2)
	Processing2D.point(0, 0)
	Processing2D.animate()
	x = x + 0.01
	if Processing2D.mousePressed()
		println("mouse pressed")
	end
end


text:
const prime = Symbol("'")

"""
    isindex(ex)

Check if the input `ex` is a valid index for a tensor. It should be a `Symbol` or an `Int`,
and prime of them are also allowed.
"""
function isindex(ex)
    if isa(ex, Symbol) || isa(ex, Int)
        return true
    elseif isa(ex, Expr) && ex.head == prime && length(ex.args) == 1
        return isindex(ex.args[1])
    else
        return false
    end
end

"""
    istensor(ex::Expr)

Check if the input `ex` is a tensor object with valid indices. The possible forms of the
tensor could be:
1. `A[]` with head `:ref`
2. `A[a, b; c, d]` with head `:ref`
3. `A[a, b, c]` with head `:ref`
4. `A[a b c]` with head `:typed_hcat`
5. `A[a b; c d]` with head `:typed_vcat`
6. `A[(a,b);(c,d)]` with head `:typed_vcat`
where `A` can itself be an expression, such as `A[1:2:end, :][a, b, c]`.
"""
istensor(ex) = false
function istensor(ex::Expr)
    if ex.head == :ref || ex.head == :typed_hcat
        if length(ex.args) == 1
            # A[]
            return true
        elseif isa(ex.args[2], Expr) && ex.args[2].head == :parameters
            # A[a,b,...;c,d,...]
            return all(isindex, ex.args[2].args) || all(isindex, ex.args[3:end])
        else
            # A[a,b,c,...] or A[a b c ...]
            return all(isindex, ex.args[2:end])
        end
    elseif ex.head == :typed_vcat && length(ex.args) == 3
        length(ex.args) == 3 || return false
        if isa(ex.args[2], Expr) && (ex.args[2].head == :row || ex.args[2].head == :tuple)
            # row: A[a b ...; ...]; or tuple: A[(a,b,...); ...]
            all(isindex, ex.args[2].args) || return false
        else
            # A[a; ...]
            isindex(ex.args[2]) || return false
        end
        if isa(ex.args[3], Expr) && (ex.args[3].head == :row || ex.args[3].head == :tuple)
            # row: A[...; b c ...]; or tuple: A[...;(b,c,...)]
            all(isindex, ex.args[3].args) || return false
        else
            # A[...; c]
            isindex(ex.args[3]) || return false
        end
        return true
    end
    return false
end

"""
    isgeneraltensor(ex::Expr)

Check if the input `ex` is a general tensor in one of the following forms:
1. `ex` itself is a tensor
2. `+ tensor` or `- tensor`
3. `conj(tensor)`, or `adjoint(tensor)`, or `tensor'`, or `transpose(tensor)`
4. `α * tensor` or `tensor * α`, where `α` is a scalar
5. `tensor/α` or `α\tensor`, where `α` is a scalar
"""
isgeneraltensor(ex) = false
function isgeneraltensor(ex::Expr)
    if istensor(ex)
        return true
    elseif ex.head == :call && ex.args[1] == :+ && length(ex.args) == 2
        # unary plus
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :- && length(ex.args) == 2
        # unary minus
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :conj && length(ex.args) == 2
        # conjugation
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :adjoint && length(ex.args) == 2
        # adjoint
        return isgeneraltensor(ex.args[2])
    elseif ex.head == prime && length(ex.args) == 1
        # adjoint
        return isgeneraltensor(ex.args[1])
    elseif ex.head == :call && ex.args[1] == :transpose && length(ex.args) == 2
        # conjugation
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :*
        # scalar multiplication
        count = 0
        for i = 2:length(ex.args)
            if isgeneraltensor(ex.args[i])
                count += 1
            elseif !isscalarexpr(ex.args[i])
                return false
            end
        end
        return count == 1
    elseif ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3
        # scalar multiplication
        return (isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2]))
    elseif ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3
        # scalar multiplication
        return (isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3]))
    end
    return false
end

"""
    isscalarexpr(ex::Expr)

Check if the input `ex` is a scalar expression with no indices.
"""
isscalarexpr(ex::Symbol) = true
isscalarexpr(ex::Number) = true
isscalarexpr(ex) = true
function isscalarexpr(ex::Expr)
    if ex.head == :call && ex.args[1] == :scalar
        return istensorexpr(ex.args[2])
    elseif ex.head in (:ref, :typed_vcat, :typed_hcat)
        return false
    else
        return all(isscalarexpr, ex.args)
    end
end

"""
    istensorexpr(ex)

Check if the input `ex` is a tensor expression, i.e. something that can be evaluated to a
tensor.
"""
function istensorexpr(ex)
    if isgeneraltensor(ex)
        return true
    elseif isa(ex, Expr) && ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-)
        # linear combination of general tensors
        return all(istensorexpr, ex.args[2:end])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :*
        # multiplications of several general tensors
        count = 0
        for i = 2:length(ex.args)
            if istensorexpr(ex.args[i])
                count += 1
            elseif !isscalarexpr(ex.args[i])
                return false
            end
        end
        return count > 0
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3
        return istensorexpr(ex.args[2]) && isscalarexpr(ex.args[3])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3
        return istensorexpr(ex.args[3]) && isscalarexpr(ex.args[2])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :conj && length(ex.args) == 2
        # conj all general tensors
        return istensorexpr(ex.args[2])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :adjoint && length(ex.args) == 2
        # adjoint all general tensors
        return istensorexpr(ex.args[2])
    elseif isa(ex, Expr) && ex.head == prime
        # adjoint all general tensors
        return istensorexpr(ex.args[1])
    end
    return false
end

"""
    hastraceindices(ex)

Check if the input tensor `ex` has indices that needed to be contracted within the tensor,
i.e., whether any trace operation is applied on the tensor.
"""
function hastraceindices(ex)
    obj, leftind, rightind, = decomposegeneraltensor(ex)
    allind = vcat(leftind, rightind)
    return length(allind) != length(unique(allind))
end

"""
    istensorcontraction(ex)

Check if the input `ex` is an expression for a tensor contraction operation in the form such
as `A[a b]*B[c a]*...`.
"""
function istensorcontraction(ex)
    if isa(ex, Expr) && ex.head == :call && ex.args[1] == :*
        return count(istensorexpr, ex.args[2:end]) >= 2
    end
    return false
end

"""
    isassignment(ex)

Check if the input `ex` is an assignment to an existing tensor.
"""
isassignment(ex) = false
isassignment(ex::Expr) = (ex.head == :(=) || ex.head == :(+=) || ex.head == :(-=))

"""
    isdefinition(ex)

Check if the input `ex` is an expression that create a new tensor by definition.
"""
isdefinition(ex) = false
isdefinition(ex::Expr) = (ex.head == :(:=) || ex.head == :(≔))


text:
module WeatherStation

using CSV
using DataFrames
using Dates
using Dash
using Glob
using Statistics

include("server.jl")
include("data.jl")

end # module


text:
"""
    @unimpl

Implement a function from a signature which raises an error.

# Example

```jldoctest; setup = :(import Neet.@unimpl)
julia> abstract type A end

julia> struct B <: A end

julia> @unimpl foo(a::A)
foo (generic function with 1 method)

julia> foo(B())
ERROR: `foo(a::A)` is not implemented
[...]

julia> foo(b::B) = "works"
foo (generic function with 2 methods)

julia> foo(B())
"works"
```
"""
macro unimpl(func)
    if func.head != :call && func.head != :where
        error("illformed function identifier")
    end
    expr = Expr(:string, "`$(func)` is not implemented")
    :($(esc(func)) = $(esc(error))($(esc(expr))))
end


text:
export inner

"""
    inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
         ) where {N}

Compute the inner product ``xᵀ A y``.

### Input

- `x` -- vector on the left
- `A` -- matrix
- `y` -- vector on the right

### Output

The (scalar) result of the multiplication.
"""
function inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
              ) where {N}
    dot(x, A * y)
end

"""
    _vector_type(T)

Return a corresponding vector type with respect to type `T`.

### Input

- `T` -- vector or matrix type

### Output

A vector type that corresponds in some sense (see Notes below) to `T`.

### Notes

- If `T` is a sparse vector or a sparse matrix, the corresponding type is also
  a sparse vector.
- If `T` is a regular vector (i.e. `Vector`) or a regular matrix (i.e. `Matrix`),
  the corresponding type is also a regular vector.
- Otherwise, the corresponding type is a regular vector.
"""
function _vector_type end

"""
    _matrix_type(T)

Return a corresponding matrix type with respect to type `T`.

### Input

- `T` -- vector or matrix type

### Output

A matrix type that corresponds in some sense (see Notes below) to `T`.

### Notes

- If `T` is a sparse vector or a sparse matrix, the corresponding type is also
  a sparse matrix.
- If `T` is a regular vector (i.e. `Vector`) or a regular matrix (i.e. `Matrix`),
  the corresponding type is also a regular matrix.
- Otherwise, the corresponding type is a regular matrix.
"""
function _matrix_type end

_vector_type(::Type{<:AbstractSparseArray{T}}) where T = SparseVector{T, Int}
_vector_type(VT::Type{<:AbstractVector{T}}) where T = VT
_vector_type(::Type{<:AbstractMatrix{T}}) where T = Vector{T}

_matrix_type(::Type{<:AbstractVector{T}}) where T = Matrix{T}
_matrix_type(MT::Type{<:AbstractMatrix{T}}) where T = MT
_matrix_type(::Type{<:AbstractSparseVector{T}}) where T = SparseMatrixCSC{T, Int}


text:
eval(quote
    using .MultivariatePolynomials: AbstractVariable, AbstractMonomialLike,
                                    exponents, variables, powers
end)

eval(load_kron_multivariate())


text:
struct TextReader
    file::String
    vocab::Vocab
end

Base.IteratorSize(::Type{TextReader}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{TextReader}) = Base.HasEltype()
Base.eltype(::Type{TextReader}) = Vector{Int}

function Base.iterate(r::TextReader, s=nothing)
    s === nothing && (s = open(r.file))
    if eof(s); close(s); return nothing; end
    v = r.vocab
    [ get(v.w2i, w, v.unk) for w in v.tokenizer(readline(s)) ], s
end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
import Base: size, deepcopy

abstract type Constraint end

struct GridCSP
    variables::Array{Int, 2}
    constraints::Array{Constraint}
    domains::Array{BitArray, 2}
end

function GridCSP(constraints::Array{Constraint, 1}, k::Int)
    vars = -ones(Int, k, k)
    domains = Array{BitArray, 2}(undef, k, k)
    for i = 1:k
        for j = 1:k
            domains[i, j] = BitArray(undef, k)
            domains[i, j] .+= 1
        end
    end
    return GridCSP(vars, constraints, domains)
end

size(csp::GridCSP) = size(csp.domains, 1)
deepcopy(csp::GridCSP) = GridCSP(
    deepcopy(csp.variables),
    deepcopy(csp.constraints),
    deepcopy(csp.domains)
)

struct RowConstraint <: Constraint
    idx::Int
end

struct ColumnConstraint <: Constraint
    idx::Int
end

struct CellConstraint <: Constraint
    value::Int
    coord::Tuple{Int, Int}
end

function set_domain!(i, j, value, domains)
    domains[i, j] .= 0
    domains[i, j][value] = 1
end

function set_cell!(i, j, value, variables)
    variables[i, j] = value
end

"""
Fills in values that have been completely determined.
"""
function fill_in!(csp::GridCSP)
    k = size(csp)
    for i = 1:k
        for j = 1:k
            dom = csp.domains[i, j]
            if sum(dom) == 1 && csp.variables[i, j] < 0
                v = argmax(dom)
                csp.variables[i, j] = v
            end
        end
    end
end

function _eliminate_rc!(doms::Array{BitArray}, vars::Array{Int})
    k = length(doms)
    # for each value in the row or column
    for (i, v) ∈ enumerate(vars)
        # if the value has been determined
        if v > 0
            for j = 1:k
                # remove the value as a candidate from every other cell in the
                # row/column
                if j ≠ i
                    doms[j][v] = 0
                end
            end
        end
    end
end

"""
Eliminates candidate values for a row constraint.
"""
function eliminate!(constr::RowConstraint, csp::GridCSP)
    return _eliminate_rc!(csp.domains[constr.idx, :], csp.variables[constr.idx, :])
end

"""
Eliminates candidate values for a column constraint.
"""
function eliminate!(constr::ColumnConstraint, csp::GridCSP)
    return _eliminate_rc!(csp.domains[:, constr.idx], csp.variables[:, constr.idx])
end

"""
Eliminates values from the domains of the Kenken's cells based on the constraints.
"""
function eliminate!(csp::GridCSP)
    # To complete an elimination pass, we eliminate for each constraint and fill
    # in values that are completely determined by each step
    for constr ∈ csp.constraints
        eliminate!(constr, csp)
        fill_in!(csp)
    end
end

"""
A Kenken is consistent if there are no domains that have no more viable candidates.
"""
function consistent(csp::GridCSP)
    return all(sum.(csp.domains) .> 0)
end

"""
A Kenken is complete if there are no values that have not been completely determined.
"""
function complete(csp::GridCSP)
    return all(csp.variables .> 0)
end

"""
Find a cell whose value has not been determined which contains the minimum number of
viable candidates.
"""
function find_next_unset(csp::GridCSP)
    k = size(csp)
    min_cands = k
    imin, jmin = nothing, nothing
    for i = 1:k
        for j = 1:k
            s = sum(csp.domains[i, j])
            if csp.variables[i, j] < 0 && s ≤ min_cands
                min_cands = s
                imin, jmin = i, j
                if s == 2
                    return imin, jmin
                end
            end
        end
    end
    return imin, jmin
end

"""
Find solutions to arbitrary GridCSP by backtracking.
"""
function backtrack(csp::GridCSP, n_solutions::Int; randomize::Bool = false)
    solutions = GridCSP[]
    backtrack!(csp, solutions, n_solutions, randomize)
    return n_solutions == 1 ? solutions[1] : solutions
end

function backtrack!(csp::GridCSP, solutions::Array{GridCSP, 1}, n_solutions::Int, randomize::Bool)
    eliminate!(csp)
    cons = consistent(csp)
    if !cons
        return
    end
    if complete(csp) && cons
        push!(solutions, csp)
        return
    end
    k = size(csp)
    i, j = find_next_unset(csp)
    dom = csp.domains[i, j]
    cands = 1:k
    cands = cands[dom[cands] .== 1]
    if randomize
        cands = cands[randperm(length(cands))]
    end
    for l = cands
        cspcpy = deepcopy(csp)
        set_domain!(i, j, l, cspcpy.domains)
        set_cell!(i, j, l, cspcpy.variables)
        backtrack!(cspcpy, solutions, n_solutions, randomize)
        # if we have more solutions than we want, terminate
        if length(solutions) ≥ n_solutions
            break
        end
    end
    return solutions[1:min(length(solutions), n_solutions)]
end

max_stars_repo_name: gindij/CSPSolver

max_issues_repo_name: gindij/CSPSolver

max_forks_repo_name: gindij/CSPSolver


text:
# SteinTensorizedKernel
#
# These represent a family of kernels of the form
#
# K(x, y) = \prod_i ki(x_i, y_i)
#
# where ki(x,y) is a kernel in one dimension. In general, ki
# will be the same for each dimension i; the exception being
# when the support is not the same amongst all dimensions.

abstract type SteinTensorizedKernel <: SteinKernel end

### METHODS TO IMPLEMENT
function ki(ker::SteinTensorizedKernel, x::Float64, y::Float64)
    error("Must implement the method ki")
end

function gradxki(ker::SteinTensorizedKernel, x::Float64, y::Float64)
    error("Must implement the method gradxki")
end

function gradxyki(ker::SteinTensorizedKernel, x::Float64, y::Float64)
    error("Must implement the method gradxyki")
end

### INDUCED METHODS
function ki(ker::SteinTensorizedKernel, x::Float64, y::Float64, i::Int)
    # if not overriden, we assume each dimension is the same
    ki(ker, x, y)
end

function gradxki(ker::SteinTensorizedKernel, x::Float64, y::Float64, i::Int)
    gradxki(ker, x, y)
end

function gradxyki(ker::SteinTensorizedKernel, x::Float64, y::Float64, i::Int)
    gradxyki(ker, x, y)
end

# this is a helpful utility for computing \prod_{j != i} kj(x_j, y_j)
function k_minusi(ker::SteinTensorizedKernel,
                  x::Array{Float64, 1},
                  y::Array{Float64, 1},
                  i::Int)
    d = length(x)
    noti = filter(x -> x != i, 1:d)
    prod([ki(ker, x[j], y[j], j) for j in noti])
end

function k(ker::SteinTensorizedKernel, x::Array{Float64,1}, y::Array{Float64,1})
    d = length(x)
    prod([ki(ker, x[i], y[i], i) for i=1:d])
end

function gradxk(ker::SteinTensorizedKernel, x::Array{Float64,1}, y::Array{Float64,1})
    d = length(x)
    [gradxki(ker, x[i], y[i], i) * k_minusi(ker, x, y, i) for i=1:d]
end

function gradxyk(ker::SteinTensorizedKernel, x::Array{Float64,1}, y::Array{Float64,1})
    d = length(x)
    sum([gradxyki(ker, x[i], y[i], i) * k_minusi(ker, x, y, i) for i=1:d])
end


text:
function testdata()
    X = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]
    y = [0.5, 1.0, 1.5, 2.0]
    f = Constant()
    Fx = covariate(f, X)
    g = Identity()
    k = Gaussian()
    return X, y, Fx, f, k, g
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
function modes_x(y,byp,H,Np,B,m,omega)

# Computes the longitudinal component of the magnetic
# perturbation using a solution for the meridional
# component
#
# input
#  y     - grid points in meridional coordinate
#  byp   - meridional magnetic perturbation
#  H     - thickness of layer (km)
#  Np    - dimensionless buoyancy frequency  (N/Omega)
#  B     - radial magnetic field  (T)
#  m     - angular order
#  omega0 - initial estimate of frequency

# a few constants   (must be consistent with ModelParameters)
Omega = 0.7292e-4;
R = 3.48e6;
N = Np * Omega;
k = pi/(1000.0*H);

# compute model parameters
C,M = ModelParameters(H,Np,B,m,omega);

# evaluate derivative matrices
d1,d2 = derivatives(y);

#   right-hand side of Eq. (17) for bxp
n = length(y);
y2 = y.*y;
yc2 = -(y2 .- 1.0);
a1 = -(m^2*ones(n) + M * yc2);

#   left-hand side
a2 = -((C)im * y .* byp);
a3 = (m * d1*byp)im;

bxp = (a2 + a3)./a1;

return bxp

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
using MLStyle
using Random: AbstractRNG
export @parameterized

# A fold over ASTs. Example usage in `replace`
function foldast(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, args...) => branch(head, map(go, args); kwargs...)
            x                   => leaf(x; kwargs...)
        end
    end

    return go
end

function replace(p, f, expr)
    leaf(s) = p(s) ? f(s) : s
    branch(head, newargs) = Expr(head, newargs...)
    foldast(leaf, branch)(expr)
end

# from https://thautwarm.github.io/MLStyle.jl/latest/tutorials/capture.html
function capture(template, ex, action)
    let template = Expr(:quote, template)
        quote
            @match $ex begin
                $template => $action
                _         => nothing
            end
        end
    end
end

macro capture(template, ex, action)
    capture(template, ex, action) |> esc
end

function _parameterized(__module__, expr)
    @capture ($op($μ($(p...)), $base)) expr begin
        @assert op ∈ [:<<, :≪, :≃]
        μbase = Symbol(:__, μ, :_base)

        μ = esc(μ)
        base = esc(base)

        q = quote
            struct $μ{N,T} <: MeasureTheory.ParameterizedMeasure{N}
                par :: NamedTuple{N,T}

            end

            MeasureTheory.basemeasure(::$μ) = $base
        end

        if !isempty(p)
            # e.g. Normal(μ,σ) = Normal((μ=μ, σ=σ))
            pnames = QuoteNode.(p)
            push!(q.args, :($μ($(p...)) = $μ(NamedTuple{($(pnames...),)}(($(p...),)))))
        end

        return q
    end

    @capture $μ($(p...)) expr begin
        μ = esc(μ)

        q = quote
            struct $μ{N,T} <: MeasureTheory.ParameterizedMeasure{N}
                par :: NamedTuple{N,T}
            end
        end

        if !isempty(p)
            # e.g. Normal(μ,σ) = Normal((μ=μ, σ=σ))
            pnames = QuoteNode.(p)
            push!(q.args, :($μ($(p...)) = $μ(NamedTuple{($(pnames...),)}(($(p...),)))))
        end

        return q
    end
end

"""
    @parameterized <declaration>

The <declaration> gives a measure and its default parameters, and specifies
its relation to its base measure. For example,

    @parameterized Normal(μ,σ)

declares the `Normal` is a measure with default parameters `μ and σ`. The result is equivalent to
```
struct Normal{N,T} <: ParameterizedMeasure{N}
    par :: NamedTuple{N,T}
end

KeywordCalls.@kwstruct Normal(μ,σ)

Normal(μ,σ) = Normal((μ=μ, σ=σ))
```

See [KeywordCalls.jl](https://github.com/cscherrer/KeywordCalls.jl) for details on `@kwstruct`.
"""
macro parameterized(expr)
    _parameterized(__module__, expr)
end

# (@macroexpand @parameterized Normal(μ,σ) ≃ (1/sqrt2π) * Lebesgue(X)) |> MacroTools.prettify

using MLStyle

macro μσ_methods(ex)
    _μσ_methods(__module__, ex)
end

function _μσ_methods(__module__, ex)
    @match ex begin
        :($dist($(args...))) => begin
            argnames = QuoteNode.(args)

            d_args = (:(d.$arg) for arg in args)

            method_μσ = KeywordCalls._kwstruct(__module__, :($dist($(args...), μ, σ)))
            method_μ  = KeywordCalls._kwstruct(__module__, :($dist($(args...), μ)))
            method_σ  = KeywordCalls._kwstruct(__module__, :($dist($(args...), σ)))

            q = quote

                 $method_μσ
                 $method_μ
                 $method_σ

                function Base.rand(rng::AbstractRNG, T::Type, d::$dist{($(argnames...), :μ, :σ)})
                    d.σ * rand(rng, T, $dist($(d_args...))) + d.μ
                end

                function MeasureTheory.logdensity(d::$dist{($(argnames...), :μ, :σ)}, x)
                    z = (x - d.μ) / d.σ
                    return logdensity($dist($(d_args...)), z) - log(d.σ)
                end

                function Base.rand(rng::AbstractRNG, T::Type, d::$dist{($(argnames...), :σ)})
                    d.σ * rand(rng, T, $dist($(d_args...)))
                end

                function MeasureTheory.logdensity(d::$dist{($(argnames...), :σ)}, x)
                    z = x / d.σ
                    return logdensity($dist($(d_args...)), z) - log(d.σ)
                end

                function Base.rand(rng::AbstractRNG, T::Type, d::$dist{($(argnames...), :μ)})
                    rand(rng, T, $dist($(d_args...))) + d.μ
                end

                function MeasureTheory.logdensity(d::$dist{($(argnames...), :μ)}, x)
                    z = x - d.μ
                    return logdensity($dist($(d_args...)), z)
                end
            end

            return q
        end
    end
end

macro σ_methods(ex)
    _σ_methods(__module__, ex)
end

function _σ_methods(__module__, ex)
    @match ex begin
        :($dist($(args...))) => begin
            argnames = QuoteNode.(args)

            d_args = (:(d.$arg) for arg in args)

            method_σ  = KeywordCalls._kwstruct(__module__, :($dist($(args...), σ)))

            q = quote
                $method_σ

                function Base.rand(rng::AbstractRNG, T::Type, d::$dist{($(argnames...), :σ)})
                    d.σ * rand(rng, T, $dist($(d_args...)))
                end

                function MeasureTheory.logdensity(d::$dist{($(argnames...), :σ)}, x)
                    z = x / d.σ
                    return logdensity($dist($(d_args...)), z) - log(d.σ)
                end
            end

            return q
        end
    end
end

"""
    @half dist([paramnames])

Starting from a symmetric univariate measure `dist ≪ Lebesgue(ℝ)`, create a new
measure `Halfdist ≪ Lebesgue(ℝ₊)`. For example,

    @half Normal()

creates `HalfNormal()`, and

    @half StudentT(ν)

creates `HalfStudentT(ν)`.
"""
macro half(ex)
    esc(_half(ex))
end

function _half(ex)
    @match ex begin
        :($dist($(args...))) => begin
            halfdist = Symbol(:Half, dist)

            quote
                struct $halfdist{N,T} <: ParameterizedMeasure{N}
                    par :: NamedTuple{N,T}
                end

                unhalf(μ::$halfdist) = $dist(getfield(μ, :par))

                function MeasureTheory.basemeasure(μ::$halfdist)
                    b = basemeasure(unhalf(μ))
                    @assert basemeasure(b) == Lebesgue(ℝ)
                    lw = b.logweight
                    return WeightedMeasure(logtwo + lw, Lebesgue(ℝ₊))
                end

                function MeasureTheory.logdensity(μ::$halfdist, x)
                    return logdensity(unhalf(μ), x)
                end

                function Base.rand(rng::AbstractRNG, T::Type, μ::$halfdist)
                    return abs(rand(rng, T, unhalf(μ)))
                end

                (::$halfdist ≪ ::Lebesgue{ℝ₊}) = true
            end
        end
    end
end


text:
using Treesitter

parser = Parser()
set_language(parser, tree_sitter_c())
codestring = read(joinpath(@__DIR__, "corpus/easystruct.h"), String)

tree =
    ts_parser_parse_string(parser.ptr, C_NULL, codestring, length(codestring))
root_node = node_rootof(tree)

tc = cursor_new(root_node)

nodecurr = cursor_current_node(tc)

node_type(nodecurr)

crs = cursor_goto_parent!(crs)

node_type(crs)


text:
module WildBootTests
export BootTestResult, wildboottest, teststat, stattype, p, padj, reps, repsfeas, nbootclust, dof, dof_r, plotpoints, peak, ci, dist, statnumer, statvar, auxweights

using LinearAlgebra, Random, Distributions, SortingAlgorithms, Printf

include("StrBootTest.jl")
include("utilities.jl")
include("estimators.jl")
include("init.jl")
include("WRE.jl")
include("nonWRE.jl")
include("plot-CI.jl")
include("interface.jl")

# top-level computation routine for OLS/arubin (and score BS on IV/2SLS); split off to reduce latency when just doing WRE
function boottestOLSARubin!(o::StrBootTest{T}) where T
	if !o.initialized
		Init!(o)
	elseif !o.null
		NoNullUpdate!(o)
		return
	end

	if o.Nw > 1  # if more than one weight group to save memory, make on every call to boottest(), not just once in Init!()
		Random.seed!(o.rng,o.seed)
		MakeWildWeights!(o, last(o.WeightGrp[1])-1, first=true)
	end

	MakeInterpolables!(o)  # make stuff that depends linearly on r, possibly by interpolating, for first weight group

	MakeNonWREStats!(o, 1)
	for w ∈ 2:o.Nw  # do group 1 first because it includes col 1, which is all that might need updating in constructing CI in WCU
		MakeWildWeights!(o, length(o.WeightGrp[w]), first=false)
		MakeNonWREStats!(o, w)
	end
	!o.bootstrapt && UpdateBootstrapcDenom!(o)

	o.BFeas = isnan(o.dist[1]) ? 0 : sum(.!(isnan.(o.dist) .| isinf.(o.dist))) - 1
	o.distCDR = zeros(T,0,0)
	nothing
end

# top-level computation routine for non-arubin WRE; split off to reduce latency when just doing other tests
function boottestWRE!(o::StrBootTest{T}) where T
	if !o.initialized
		Init!(o)
	elseif !o.null
		NoNullUpdate!(o)
		return
	end

	if o.Nw > 1  # if more than one weight group to save memory, make on every call to boottest(), not just once in Init!()
		Random.seed!(o.rng,o.seed)
		MakeWildWeights!(o, last(o.WeightGrp[1])-1, first=true)
	end

	PrepWRE!(o)

	MakeWREStats!(o, 1)
	for w ∈ 2:o.Nw  # do group 1 first because it includes col 1, which is all that might need updating in constructing CI in WCU
		MakeWildWeights!(o, length(o.WeightGrp[w]), first=false)
		MakeWREStats!(o, w)
	end
	!o.bootstrapt && UpdateBootstrapcDenom!(o)

	o.BFeas = isnan(o.dist[1]) ? 0 : sum(.!(isnan.(o.dist) .| isinf.(o.dist))) - 1
	o.distCDR = zeros(T,0,0)
	nothing
end

# if not imposing null and we have returned to boottest!(), then dof=1 or 2; we're plotting or finding CI, and only test stat, not distribution, changes with r
function NoNullUpdate!(o::StrBootTest{T} where T)
	if o.WREnonARubin
		o.numer[:,1] = o.R * o.DGP.Rpar * o.β̈s[1] - o.r
	elseif o.arubin
		EstimateARubin!(o.DGP, o, o.r)
		o.numer[:,1] = o.v_sd * @view o.DGP.β̈[o.kX₁+1:end,:]  # coefficients on excluded instruments in arubin OLS
	else
		o.numer[:,1] = o.v_sd * (o.R * (o.ml ? o.β̈ : iszero(o.κ) ? o.M.β̈ : o.M.Rpar * o.M.β̈) - o.r)  # Analytical Wald numerator; if imposing null then numer[:,1] already equals this. If not, then it's 0 before this
	end
	o.dist[1] = isone(o.dof) ? o.numer[1] / sqrtNaN(o.statDenom[1]) : o.numer[:,1]'invsym(o.statDenom)*o.numer[:,1]
	nothing
end

# compute bootstrap-c denominator from all bootstrap numerators
function UpdateBootstrapcDenom!(o::StrBootTest{T} where T)
	numer1 = o.numer[:,1]
	o.numersum = rowsum(o.numer) - numer1
	o.statDenom = (o.numer * o.numer' - numer1 * numer1' - o.numersum * o.numersum' / o.B) / o.B
	if o.sqrt
		o.dist .= o.numer ./ sqrtNaN.(o.statDenom)
	else
		negcolquadform!(o, o.dist, -invsym(o.statDenom), o.numer)  # to reduce latency by minimizing #=@tturbo=# instances, work with negative of colquadform in order to fuse code with colquadformminus!
	end
	nothing
end

include("precompile_WildBootTests.jl")  # https://timholy.github.io/SnoopCompile.jl/stable/snoopi_deep_parcel/#SnoopCompile.write
_precompile_()

end


text:
# # Default recipes
# Includes stuff from Base/stdlib.
# -------------------------------------------------
# ## Dates & Times

dateformatter(dt) = string(Date(Dates.UTD(dt)))
datetimeformatter(dt) = string(DateTime(Dates.UTM(dt)))
timeformatter(t) = string(Dates.Time(Dates.Nanosecond(round(t))))

@recipe f(::Type{Date}, dt::Date) = (dt -> Dates.value(dt), dateformatter)
@recipe f(::Type{DateTime}, dt::DateTime) =
    (dt -> Dates.value(dt), datetimeformatter)
@recipe f(::Type{Dates.Time}, t::Dates.Time) = (t -> Dates.value(t), timeformatter)
@recipe f(::Type{P}, t::P) where {P<:Dates.Period} =
    (t -> Dates.value(t), t -> string(P(round(t))))

# -------------------------------------------------
# ## Characters

@recipe f(::Type{<:AbstractChar}, ::AbstractChar) = (string, string)