text:
module RegressionWrapper

using JavaCall

const AbstractMultipleLinearRegression = @jimport org.hipparchus.stat.regression.AbstractMultipleLinearRegression
const Class = @jimport java.lang.Class
const GLSMultipleLinearRegression = @jimport org.hipparchus.stat.regression.GLSMultipleLinearRegression
const JString = @jimport java.lang.String
const MillerUpdatingRegression = @jimport org.hipparchus.stat.regression.MillerUpdatingRegression
const MultipleLinearRegression = @jimport org.hipparchus.stat.regression.MultipleLinearRegression
const OLSMultipleLinearRegression = @jimport org.hipparchus.stat.regression.OLSMultipleLinearRegression
const Object = @jimport java.lang.Object
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const RegressionResults = @jimport org.hipparchus.stat.regression.RegressionResults
const SimpleRegression = @jimport org.hipparchus.stat.regression.SimpleRegression
const UpdatingMultipleLinearRegression = @jimport org.hipparchus.stat.regression.UpdatingMultipleLinearRegression

export AbstractMultipleLinearRegression
export GLSMultipleLinearRegression
export MillerUpdatingRegression
export MultipleLinearRegression
export OLSMultipleLinearRegression
export RegressionResults
export SimpleRegression
export UpdatingMultipleLinearRegression
export add_data
export add_observation
export add_observations
export append
export calculate_adjusted_r_squared
export calculate_hat
export calculate_r_squared
export calculate_residual_sum_of_squares
export calculate_total_sum_of_squares
export clear
export equals
export estimate_error_variance
export estimate_regressand_variance
export estimate_regression_parameters
export estimate_regression_parameters_standard_errors
export estimate_regression_parameters_variance
export estimate_regression_standard_error
export estimate_residuals
export get_adjusted_r_squared
export get_class
export get_covariance_of_parameters
export get_diagonal_of_hat_matrix
export get_error_sum_squares
export get_intercept
export get_intercept_std_err
export get_mean_square_error
export get_n
export get_number_of_parameters
export get_order_of_regressors
export get_parameter_estimate
export get_parameter_estimates
export get_partial_correlations
export get_r
export get_r_square
export get_r_squared
export get_regression_sum_squares
export get_significance
export get_slope
export get_slope_confidence_interval
export get_slope_std_err
export get_std_error_of_estimate
export get_std_error_of_estimates
export get_sum_of_cross_products
export get_sum_squared_errors
export get_total_sum_squares
export get_x_sum_squares
export has_intercept
export hash_code
export is_no_intercept
export new_sample_data
export notify
export notify_all
export predict
export regress
export remove_data
export set_no_intercept
export to_string
export wait

include("abstract_multiple_linear_regression.jl")
include("gls_multiple_linear_regression.jl")
include("miller_updating_regression.jl")
include("multiple_linear_regression.jl")
include("ols_multiple_linear_regression.jl")
include("regression_results.jl")
include("simple_regression.jl")
include("updating_multiple_linear_regression.jl")

end


text:
function convert_func(op, T1, W1, T2, W2 = W1)
  typ1 = LLVM_TYPES[T1]
  typ2 = LLVM_TYPES[T2]
  vtyp1 = vtype(W1, typ1)
  vtyp2 = vtype(W2, typ2)
  instrs = """
      %res = $op $vtyp2 %0 to $vtyp1
      ret $vtyp1 %res
      """
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W1,$T1}, Tuple{_Vec{$W2,$T2}}, data(v)))
  end
end
# For bitcasting between signed and unsigned integers (LLVM does not draw a distinction, but they're separate in Julia)
function identity_func(W, T1, T2)
  vtyp1 = vtype(W, LLVM_TYPES[T1])
  instrs = """
      ret $vtyp1 %0
      """
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W,$T1}, Tuple{_Vec{$W,$T2}}, data(v)))
  end
end

### `Base.convert(::Type{<:AbstractSIMDVector}, x)` methods
### These are the critical `Base.convert` methods; scalar and `VecUnroll` are implemented with respect to them.
@generated function Base.convert(::Type{Vec{W,F}}, v::Vec{W,T}) where {W,F<:Union{Float32,Float64},T<:IntegerTypesHW}
  convert_func(T <: Signed ? "sitofp" : "uitofp", F, W, T)
end

@generated function Base.convert(::Type{Vec{W,T}}, v::Vec{W,F}) where {W,F<:Union{Float32,Float64},T<:IntegerTypesHW}
  convert_func(T <: Signed ? "fptosi" : "fptoui", T, W, F)
end
@generated function Base.convert(::Type{Vec{W,T1}}, v::Vec{W,T2}) where {W,T1<:IntegerTypesHW,T2<:IntegerTypesHW}
  sz1 = sizeof(T1)::Int; sz2 = sizeof(T2)::Int
  if sz1 < sz2
    convert_func("trunc", T1, W, T2)
  elseif sz1 == sz2
    identity_func(W, T1, T2)
  else
    convert_func(((T1 <: Signed) && (T2 <: Signed)) ? "sext" : "zext", T1, W, T2)
  end
end
@generated function Base.convert(::Type{Vec{W,Float32}}, v::Vec{W,Float64}) where {W}
  convert_func("fptrunc", Float32, W, Float64, W)
end
@generated function Base.convert(::Type{Vec{W,Float64}}, v::Vec{W,Float32}) where {W}
  convert_func("fpext", Float64, W, Float32, W)
end
@inline Base.convert(::Type{<:AbstractMask{W}}, v::Vec{W,Bool}) where {W} = tomask(v)
@inline Base.convert(::Type{M}, v::Vec{W,Bool}) where {W,U,M<:AbstractMask{W,U}} = tomask(v)
@inline Base.convert(::Type{<:AbstractMask{W,U} where U}, v::Vec{W,Bool}) where {W} = tomask(v)
@inline Base.convert(::Type{<:AbstractMask{L,U} where {L,U}}, v::Vec{W,Bool}) where {W} = tomask(v)
# @inline Base.convert(::Type{Mask}, v::Vec{W,Bool}) where {W} = tomask(v)
# @generated function Base.convert(::Type{<:AbstractMask{W}}, v::Vec{W,Bool}) where {W}
#     instrs = String[]
#     push!(instrs, "%m = trunc <$W x i8> %0 to <$W x i1>")
#     zext_mask!(instrs, 'm', W, '0')
#     push!(instrs, "ret i$(max(8,W)) %res.0")
#     U = mask_type_symbol(W);
#     quote
#         $(Expr(:meta,:inline))
#         Mask{$W}($LLVMCALL($(join(instrs, "\n")), $U, Tuple{_Vec{$W,Bool}}, data(v)))
#     end
# end
@inline Base.convert(::Type{Vec{W,Bit}}, v::Vec{W,Bool}) where {W,Bool} = tomask(v)

@inline Base.convert(::Type{Vec{W,T}}, v::Vec{W,T}) where {W,T<:IntegerTypesHW} = v
@inline Base.convert(::Type{Vec{W,T}}, v::Vec{W,T}) where {W,T} = v
@inline Base.convert(::Type{Vec{W,T}}, s::NativeTypes) where {W,T} = vbroadcast(Val{W}(), T(s))
@inline Base.convert(::Type{Vec{W,T}}, s::IntegerTypesHW) where {W,T<:IntegerTypesHW} = _vbroadcast(StaticInt{W}(), s % T, StaticInt{W}() * static_sizeof(T))
@inline Base.convert(::Type{V}, u::VecUnroll) where {V<:AbstractSIMDVector} = VecUnroll(fmap(Base.convert, V, getfield(u, :data)))
@inline Base.convert(::Type{V}, u::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMDVector} = u

@inline Base.convert(::Type{<:AbstractSIMDVector{W,T}}, i::MM{W,X}) where {W,X,T} = vrangeincr(Val{W}(), T(data(i)), Val{0}(), Val{X}())
@inline Base.convert(::Type{MM{W,X,T}}, i::MM{W,X}) where {W,X,T} = MM{W,X}(T(getfield(i, :i)))

@inline Base.convert(::Type{V}, v::AbstractMask{W}) where {W, T <: NativeTypesExceptBit, V <: AbstractSIMDVector{W,T}} = ifelse(v, one(T), zero(T))
@inline Base.convert(::Type{V}, v::AbstractMask{W}) where {W, V <: AbstractSIMDVector{W,Bit}} = v
@inline Base.convert(::Type{V}, v::Vec{W,Bool}) where {W, T <: Base.HWReal, V <: AbstractSIMDVector{W,T}} = ifelse(v, one(T), zero(T))

### `Base.convert(::Type{<:NativeTypes}, x)` methods. These forward to `Base.convert(::Type{Vec{W,T}}, x)`
@inline Base.convert(::Type{T}, v::AbstractSIMD{W,T}) where {T<:NativeTypes,W} = v
@inline Base.convert(::Type{T}, v::AbstractSIMD{W,S}) where {T<:NativeTypes,S,W} = convert(Vec{W,T}, v)

### `Base.convert(::Type{<:VecUnroll}, x)` methods
@inline Base.convert(::Type{VecUnroll{N,W,T,V}}, s::NativeTypes) where {N,W,T,V} = VecUnroll{N}(convert(V, s))
@inline Base.convert(::Type{VecUnroll{N,W,T,V}}, v::AbstractSIMDVector{W}) where {N,W,T,V} = VecUnroll{N}(convert(V, v))
@inline Base.convert(::Type{VecUnroll{N,W,T,V}}, v::VecUnroll{N}) where {N,W,T,V} = VecUnroll(fmap(convert, V, getfield(v, :data)))
@inline Base.convert(::Type{VecUnroll{N,W,T,V}}, v::VecUnroll{N,W,T,V}) where {N,W,T,V} = v
@generated function Base.convert(::Type{VecUnroll{N,1,T,T}}, s::NativeTypes) where {N,T}
  quote
    $(Expr(:meta,:inline))
    x = convert($T, s)
    VecUnroll((Base.Cartesian.@ntuple $(N+1) n -> x))
  end
end

# @inline Base.convert(::Type{T}, v::T) where {T} = v

@generated function splitvectortotuple(::StaticInt{N}, ::StaticInt{W}, v::AbstractMask{L}) where {N,W,L}
  N*W == L || throw(ArgumentError("Can't split a vector of length $L into $N pieces of length $W."))
  t = Expr(:tuple, :(Mask{$W}(u)))
  s = 0
  for n ∈ 2:N
    push!(t.args, :(Mask{$W}(u >>> $(s += W))))
  end
  # This `Base.convert` will dispatch to one of the following two `Base.convert` methods
  Expr(:block, Expr(:meta,:inline), :(u = data(v)), t)
end
@generated function splitvectortotuple(::StaticInt{N}, ::StaticInt{W}, v::AbstractSIMDVector{L}) where {N,W,L}
  N*W == L || throw(ArgumentError("Can't split a vector of length $L into $N pieces of length $W."))
  t = Expr(:tuple);
  j = 0
  for i ∈ 1:N
    val = Expr(:tuple)
    for w ∈ 1:W
      push!(val.args, j)
      j += 1
    end
    push!(t.args, :(shufflevector(v, Val{$val}())))
  end
  Expr(:block, Expr(:meta,:inline), t)
end
@generated function splitvectortotuple(::StaticInt{N}, ::StaticInt{W}, v::LazyMulAdd{M,O}) where {N,W,M,O}
  # LazyMulAdd{M,O}(splitvectortotuple(StaticInt{N}(), StaticInt{W}(), getfield(v, :data)))
  t = Expr(:tuple)
  for n ∈ 1:N
    push!(t.args, :(LazyMulAdd{$M,$O}(splitdata[$n])))
  end
  Expr(:block, Expr(:meta,:inline), :(splitdata = splitvectortotuple(StaticInt{$N}(), StaticInt{$W}(), getfield(v, :data))), t)
end

@generated function Base.convert(::Type{VecUnroll{N, W, T, V}}, v::AbstractSIMDVector{L}) where {N, W, T, V, L}
  if W == L # _vconvert will dispatch to one of the two above
    Expr(:block, Expr(:meta,:inline), :(convert(VecUnroll{$N,$W,$T,$V}, v)))
  else
    Expr(:block, Expr(:meta,:inline), :(convert(VecUnroll{$N,$W,$T,$V}, VecUnroll(splitvectortotuple(StaticInt{$(N+1)}(), StaticInt{$W}(), v)))))
  end
end

@inline Vec{W,T}(v::Vec{W,S}) where {W,T,S} = Base.convert(Vec{W,T}, v)
@inline Vec{W,T}(v::S) where {W,T,S<:NativeTypes} = Base.convert(Vec{W,T}, v)

@inline vsigned(v::AbstractSIMD{W,T}) where {W,T <: Base.BitInteger} = v % signed(T)
@inline vunsigned(v::AbstractSIMD{W,T}) where {W,T <: Base.BitInteger} = v % unsigned(T)

@generated function _vfloat(v::Vec{W,I}, ::StaticInt{RS}) where {W, I <: Integer, RS}
  ex = if 8W ≤ RS
    :(Base.convert(Vec{$W,Float64}, v))
  else
    :(Base.convert(Vec{$W,Float32}, v))
  end
  Expr(:block, Expr(:meta, :inline), ex)
end
@inline Base.float(v::Vec{W,I}) where {W, I <: Integer} = _vfloat(v, register_size())
@inline Base.float(v::AbstractSIMD{W,T}) where {W,T <: Union{Float32,Float64}} = v
@inline Base.float(vu::VecUnroll) = VecUnroll(fmap(vfloat, getfield(vu, :data)))
@inline Base.float(x::Union{Float32,Float64}) = x
@inline Base.float(x::UInt64) = Base.uitofp(Float64, x)
@inline Base.float(x::Int64) = Base.sitofp(Float64, x)
@inline Base.float(x::Union{UInt8,UInt16,UInt32}) = Base.uitofp(Float32, x)
@inline Base.float(x::Union{Int8,Int16,Int32}) = Base.sitofp(Float32, x)
# @inline Base.float(v::Vec{W,I}) where {W, I <: Union{UInt64, Int64}} = Vec{W,Float64}(v)

@inline vfloat_fast(v::AbstractSIMDVector{W,T}) where {W,T <: Union{Float32,Float64}} = v
@inline vfloat_fast(vu::VecUnroll{W,T}) where {W,T<:Union{Float32,Float64}} = vu
@inline vfloat_fast(vu::VecUnroll) = VecUnroll(fmap(vfloat_fast, getfield(vu, :data)))

@generated function __vfloat_fast(v::Vec{W,I}, ::StaticInt{RS}) where {W, I <: Integer, RS}
  arg = if (2W*sizeof(I) ≤ RS) || sizeof(I) ≤ 4
    :v
  elseif I <: Signed
    :(v % Int32)
  else
    :(v % UInt32)
  end
  ex = if 8W ≤ RS
    :(Vec{$W,Float64}($arg))
  else
    :(Vec{$W,Float32}($arg))
  end
  Expr(:block, Expr(:meta, :inline), ex)
end
@inline _vfloat_fast(v, ::False) = __vfloat_fast(v, register_size())
@inline _vfloat_fast(v, ::True) = float(v)

@inline Base.reinterpret(::Type{Vec{1,T}}, x::S) where {T,S<:NativeTypes} = reinterpret(T,x)
@inline vrem(x::NativeTypes, ::Type{T}) where {T} = x % T
@generated function Base.reinterpret(::Type{T1}, v::Vec{W2,T2}) where {W2, T1 <: NativeTypes, T2}
  W1 = W2 * sizeof(T2) ÷ sizeof(T1)
  Expr(:block, Expr(:meta, :inline), :(reinterpret(Vec{$W1,$T1}, v)))
end
@inline Base.reinterpret(::Type{Vec{1,T1}}, v::Vec{W,T2}) where {W,T1,T2<:Base.BitInteger} = reinterpret(T1, fuseint(v))
@generated function Base.reinterpret(::Type{Vec{W1,T1}}, v::Vec{W2,T2}) where {W1, W2, T1, T2}
  @assert sizeof(T1) * W1 == W2 * sizeof(T2)
  convert_func("bitcast", T1, W1, T2, W2)
end

@inline vunsafe_trunc(::Type{I}, v::Vec{W,T}) where {W,I,T} = convert(Vec{W,I}, v)
@inline vrem(v::AbstractSIMDVector{W,T}, ::Type{I}) where {W,I,T} = convert(Vec{W,I}, v)
@inline vrem(v::AbstractSIMDVector{W,T}, ::Type{V}) where {W,I,T,V<:AbstractSIMD{W,I}} = convert(V, v)
@inline vrem(r::IntegerTypesHW, ::Type{V}) where {W, I, V <: AbstractSIMD{W,I}} = convert(V, r % I)

ext: jl


text:
# This file includes an example for two dimensional interpolation over triangular interpolation domain.

# TODO: Relative error must decrease with decreasing vertical scaling factor and increasing interpolation point.
# TODO: Plot error levels with colors in 2D surfaces.
# TODO: Increase number of points in the mesh.
# TODO: Relative error calculation for zero. Look at MSE with scaling factor.
# TODO: Construct the interpolation domain while constructing the interpolant.

using FractalTools
using Triangulation
using PyPlot
using Statistics

close("all")

datapath = joinpath(@__DIR__,  "../../data/gmsh_data")
# datapath = joinpath(pathof(FractalTools), "data/gmsh_data")
filename = "t3.msh"

# Read mesh file.
mesh = read_mesh(joinpath(datapath, filename))
scale_factor = 1   # Rescale the mesh.
for node in mesh.nodes
    node.x /= scale_factor
    node.y /= scale_factor
end
meshes = triangular_partition(mesh)
regions = gmsh2matplotlibtri.(meshes)

# Compute z-values on mesh points.
func(x, y) = x.^2 .+ y.^2 .+ 1.
# func(x, y) = -20 * exp.(-0.2 * sqrt.(0.5 * (x.^2 .+ y.^2))) - exp.(0.5 * (cos.(2 * pi * x) + cos.(2 * pi * y))) .+ MathConstants.e .+ 20
# func(x, y) = -(y .+ 47) .* sin.(sqrt.(abs.(x ./ 2 .+ (y .+ 47)))) .- x .* sin.(sqrt.(abs.(x .- (y .+ 47))))
# func(x, y) = 100 * sqrt.(abs.(y - 0.01 * x.^2)) .+ 0.01 * abs.(x .+ 10)
number_of_regions = length(regions)
z = Vector{Vector}(undef, number_of_regions)
for k = 1 : number_of_regions
    z[k] = func(regions[k].x, regions[k].y)
end

# Refine region to evaluate interpolated function
subdiv = 3
x, y = refine(gmsh2matplotlibtri(mesh), subdiv, true)

# Compute initial and real values.
func0(xi, yi) = 0.

# Perform interpolation.
# α = 0.025
# α = 0.151
α = 0.001
interpolant, coeffs = fis(regions, z, α=α, func0=func0, num_iter=10, get_coeffs=true)

# Interpolation and error
real_values = func.(x, y)
interpolated_values = interpolant.(x, y)
absolute_error = abs.(interpolated_values - real_values)
relative_error = absolute_error ./ (abs.(real_values)) * 100
frerr = mean(absolute_error.^2)

# Filtered relative error
threshold = 0.05
mask = x.^2 .+ y.^2 .< threshold
idx = findall(mask)
relative_error_fltr = copy(relative_error)
relative_error_fltr[idx] .= 0.# TODO: NAN value

# # Plot the results.
xd = xcoords(mesh)
yd = ycoords(mesh)
zd = func.(xd, yd)
minx, meanx, maxx = minimum(x), mean(x), maximum(x)
miny, meany, maxy = minimum(y), mean(y), maximum(y)

# Plot the results
figure("Real Values")
plot_trisurf(x, y, real_values)

colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
zticks(range(minimum(real_values), maximum(real_values), length=5))
tight_layout()
display(PyPlot.gcf())

figure("Interpolated Values")
plot_trisurf(x, y, interpolated_values)
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
zticks(range(minimum(interpolated_values), maximum(interpolated_values), length=5))
tight_layout()
display(PyPlot.gcf())

figure("Absolute Values")
plot_trisurf(x, y, absolute_error)
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
zticks(range(minimum(absolute_error), maximum(absolute_error), length=5))
tight_layout()
display(PyPlot.gcf())

figure("Relative Error")
plot_trisurf(x, y, relative_error)
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(relative_error), maximum(relative_error), length=5))
tight_layout()
display(PyPlot.gcf())

figure("Relative Error Filtered")
plot_trisurf(x, y, relative_error_fltr)
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(relative_error_fltr), maximum(relative_error_fltr), length=5))
tight_layout()
display(PyPlot.gcf())

# Plot interpolation domain.
fig, ax = subplots(1)
for region in regions
    ax.triplot(region)
    ax.plot(region.x, region.y, ".", color="orange")
end
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
tight_layout()
display(PyPlot.gcf())


text:
module ExecutionTests

using BenchmarkTools
using Test

seteq(a, b) = length(a) == length(b) == length(intersect(a, b))

#########
# setup #
#########

groups = BenchmarkGroup()
groups["sum"] = BenchmarkGroup(["arithmetic"])
groups["sin"] = BenchmarkGroup(["trig"])
groups["special"] = BenchmarkGroup()

sizes = (5, 10, 20)

for s in sizes
    A = rand(s, s)
    groups["sum"][s] = @benchmarkable sum($A) seconds=3
    groups["sin"][s] = @benchmarkable(sin($s), seconds=1, gctrial=false)
end

groups["special"]["macro"] = @benchmarkable @test(1 == 1)
groups["special"]["nothing"] = @benchmarkable nothing
groups["special"]["block"] = @benchmarkable begin rand(3) end
groups["special"]["comprehension"] = @benchmarkable [s^2 for s in sizes]

function testexpected(received::BenchmarkGroup, expected::BenchmarkGroup)
    @test length(received) == length(expected)
    @test seteq(received.tags, expected.tags)
    @test seteq(keys(received), keys(expected))
    for (k, v) in received
        testexpected(v, expected[k])
    end
end

function testexpected(trial::BenchmarkTools.Trial, args...)
    @test length(trial) > 1
end

testexpected(b::BenchmarkTools.Benchmark, args...) = true

#########
# tune! #
#########

oldgroups = copy(groups)

for id in keys(groups["special"])
    testexpected(tune!(groups["special"][id]))
end

testexpected(tune!(groups["sin"], verbose = true), groups["sin"])
testexpected(tune!(groups, verbose = true), groups)

oldgroupscopy = copy(oldgroups)

loadparams!(oldgroups, params(groups), :evals, :samples)
loadparams!(oldgroups, params(groups))

@test oldgroups == oldgroupscopy == groups

#######
# run #
#######

testexpected(run(groups; verbose = true), groups)
testexpected(run(groups; seconds = 1, verbose = true, gctrial = false), groups)
testexpected(run(groups; verbose = true, seconds = 1, gctrial = false, time_tolerance = 0.10, samples = 2, evals = 2, gcsample = false), groups)

testexpected(run(groups["sin"]; verbose = true), groups["sin"])
testexpected(run(groups["sin"]; seconds = 1, verbose = true, gctrial = false), groups["sin"])
testexpected(run(groups["sin"]; verbose = true, seconds = 1, gctrial = false, time_tolerance = 0.10, samples = 2, evals = 2, gcsample = false), groups["sin"])

testexpected(run(groups["sin"][first(sizes)]))
testexpected(run(groups["sin"][first(sizes)]; seconds = 1, gctrial = false))
testexpected(run(groups["sin"][first(sizes)]; seconds = 1, gctrial = false, time_tolerance = 0.10, samples = 2, evals = 2, gcsample = false))

testexpected(run(groups["sum"][first(sizes)], BenchmarkTools.DEFAULT_PARAMETERS))

###########
# warmup #
###########

p = params(warmup(@benchmarkable sin(1)))

@test p.samples == 1
@test p.evals == 1
@test p.gctrial == false
@test p.gcsample == false

##############
# @benchmark #
##############

mutable struct Foo
    x::Int
end

const foo = Foo(-1)

t = @benchmark sin(foo.x) evals=3 samples=10 setup=(foo.x = 0)

@test foo.x == 0
@test params(t).evals == 3
@test params(t).samples == 10

b = @benchmarkable sin(x) setup=(foo.x = -1; x = foo.x) teardown=(@assert(x == -1); foo.x = 1)
tune!(b)

@test foo.x == 1
@test params(b).evals > 100

foo.x = 0
tune!(b)

@test foo.x == 1
@test params(b).evals > 100

# test variable assignment with `@benchmark args...` form
@benchmark local_var="good" setup=(local_var="bad") teardown=(@test local_var=="good")
@test_throws UndefVarError local_var
@benchmark some_var="whatever" teardown=(@test_throws UndefVarError some_var)
@benchmark foo,bar="good","good" setup=(foo="bad"; bar="bad") teardown=(@test foo=="good" && bar=="good")

# test variable assignment with `@benchmark(args...)` form
@benchmark(local_var="good", setup=(local_var="bad"), teardown=(@test local_var=="good"))
@test_throws UndefVarError local_var
@benchmark(some_var="whatever", teardown=(@test_throws UndefVarError some_var))
@benchmark((foo,bar) = ("good","good"), setup=(foo = "bad"; bar = "bad"), teardown=(@test foo == "good" && bar == "good"))

# test kwargs separated by `,`
@benchmark(output=sin(x), setup=(x=1.0; output=0.0), teardown=(@test output == sin(x)))

########
# misc #
########

# This test is volatile in nonquiescent environments (e.g. Travis)
# BenchmarkTools.DEFAULT_PARAMETERS.overhead = BenchmarkTools.estimate_overhead()
# @test time(minimum(@benchmark nothing)) == 1

@test [:x, :y, :z, :v, :w] == BenchmarkTools.collectvars(quote
           x = 1 + 3
           y = 1 + x
           z = (a = 4; y + a)
           v,w = 1,2
           [u^2 for u in [1,2,3]]
       end)

# this should take < 1 s on any sane machine
@test @belapsed(sin($(foo.x)), evals=3, samples=10, setup=(foo.x = 0)) < 1
@test @belapsed(sin(0)) < 1

@test @ballocated(sin($(foo.x)), evals=3, samples=10, setup=(foo.x = 0)) == 0
@test @ballocated(sin(0)) == 0
@test @ballocated(Ref(1)) == 2*sizeof(Int)  # 1 for the pointer, 1 for content

let fname = tempname()
    try
        ret = open(fname, "w") do f
            redirect_stdout(f) do
                x = 1
                a = nothing
                y = @btime(sin($x))
                @test y == sin(1)
                @test a === nothing
            end
        end
        s = read(fname, String)
        try
            @test occursin(r"[0-9.]+ \w*s \([0-9]* allocations?: [0-9]+ bytes\)", s)
        catch
            println(stderr, "@btime output didn't match ", repr(s))
            rethrow()
        end
    finally
        isfile(fname) && rm(fname)
    end
end

# issue #107
let time = 2
    @benchmark identity(time)
end

end # module


text:
function show()::Vector
    s = []
    try
        throw(Exception("foo"))
    catch exn
        let e = exn
            if e isa Exception
                push!(s, "foo")
            end
        end
    finally
        push!(s, "Finally")
    end
    try
        3 / 0
    catch exn
        if exn isa ZeroDivisionError
            push!(s, "ZeroDivisionError")
        end
    end
    try
        throw(Exception("foo"))
    catch exn
        push!(s, "foo_2")
    end
    return s
end

if abspath(PROGRAM_FILE) == @__FILE__
    @assert(show() == ["foo", "Finally", "ZeroDivisionError", "foo_2"])
end


text:
## We start with definitions of instruction mnemonics, broken down by category:
instruction_categories = Dict(
    # These are mnemonics that we purposefully aren't going to classify.
    # Also, this initializes an "unknown" category for us, which is nice.
    "unknown" => String[
        "hlt",
    ],
    # We might want to know if this binary does CPUID detection
    "cpuid" => String[
        "cpuid",
    ],
    "nop" => String[
        "nop",
    ],
    "call" => String[
        "call",
        "callq",
    ],
    "486" => String[
        "bswap",
        "cmpxchg",
        "invd",
        "invlpg",
        "wbinvd",
        "xadd",
    ],
    "586" => String[
        "cmpxchg8b",
        "rdmsr",
        "rdtsc",
        "wrmsr",
    ],
    "686" => String[
        "cmova",
        "cmovae",
        "cmovb",
        "cmovbe",
        "cmovc",
        "cmove",
        "cmovg",
        "cmovge",
        "cmovl",
        "cmovle",
        "cmovna",
        "cmovnae",
        "cmovnb",
        "cmovnbe",
        "cmovnc",
        "cmovne",
        "cmovng",
        "cmovnge",
        "cmovnl",
        "cmovnle",
        "cmovno",
        "cmovnp",
        "cmovns",
        "cmovnz",
        "cmovo",
        "cmovp",
        "cmovs",
        "cmovz",
        "fcmova",
        "fcmovae",
        "fcmovb",
        "fcmovbe",
        "fcmove",
        "fcmovna",
        "fcmovnae",
        "fcmovnb",
        "fcmovnbe",
        "fcmovne",
        "fcmovnu",
        "fcmovu",
        "fcomi",
        "fcomip",
        "fcompi",
        "fucomi",
        "fucomip",
        "fucompi",
        "fxrstor",
        "fxsave",
        "rdpmc",
        "sysenter",
        "sysexit",
        "ud2",
        "ud2a",
        "ud2b",
    ],
    "mmx" => String[
        "emms",
        "movd",
        "movq",
        "packssdw",
        "packsswb",
        "packuswb",
        "paddb",
        "paddd",
        "paddsb",
        "paddsw",
        "paddusb",
        "paddusw",
        "paddw",
        "pand",
        "pandn",
        "pcmpeqb",
        "pcmpeqd",
        "pcmpeqw",
        "pcmpgtb",
        "pcmpgtd",
        "pcmpgtw",
        "pmaddwd",
        "pmulhw",
        "pmullw",
        "por",
        "pslld",
        "psllq",
        "psllw",
        "psrad",
        "psraw",
        "psrld",
        "psrlq",
        "psrlw",
        "psubb",
        "psubd",
        "psubsb",
        "psubsw",
        "psubusb",
        "psubusw",
        "psubw",
        "punpckhbw",
        "punpckhdq",
        "punpckhwd",
        "punpcklbw",
        "punpckldq",
        "punpcklwd",
        "pxor",
    ],
    "sse" => String[
        "addps",
        "addss",
        "andnps",
        "andps",
        "cmpeqps",
        "cmpeqss",
        "cmpleps",
        "cmpless",
        "cmpltps",
        "cmpltss",
        "cmpneqps",
        "cmpneqss",
        "cmpnleps",
        "cmpnless",
        "cmpnltps",
        "cmpnltss",
        "cmpordps",
        "cmpordss",
        "cmpps",
        "cmpss",
        "cmpunordps",
        "cmpunordss",
        "comiss",
        "cvtpi2ps",
        "cvtps2pi",
        "cvtsi2ss",
        "cvtss2si",
        "cvttps2pi",
        "cvttss2si",
        "divps",
        "divss",
        "ldmxcsr",
        "maskmovq",
        "maxps",
        "maxss",
        "minps",
        "minss",
        "movaps",
        "movhlps",
        "movhps",
        "movlhps",
        "movlps",
        "movmskps",
        "movntps",
        "movntq",
        "movss",
        "movups",
        "mulps",
        "mulss",
        "orps",
        "pavgb",
        "pavgw",
        "pextrw",
        "pinsrw",
        "pmaxsw",
        "pmaxub",
        "pminsw",
        "pminub",
        "pmovmskb",
        "pmulhuw",
        "prefetchnta",
        "prefetcht0",
        "prefetcht1",
        "prefetcht2",
        "psadbw",
        "pshufw",
        "rcpps",
        "rcpss",
        "rsqrtps",
        "rsqrtss",
        "sfence",
        "shufps",
        "sqrtps",
        "sqrtss",
        "stmxcsr",
        "subps",
        "subss",
        "ucomiss",
        "unpckhps",
        "unpcklps",
        "xorps",
    ],
    "sse2" => String[
        "addpd",
        "addsd",
        "andnpd",
        "andpd",
        "clflush",
        "cmpeqpd",
        "cmpeqsd",
        "cmplepd",
        "cmplesd",
        "cmpltpd",
        "cmpltsd",
        "cmpneqpd",
        "cmpneqsd",
        "cmpnlepd",
        "cmpnlesd",
        "cmpnltpd",
        "cmpnltsd",
        "cmpordpd",
        "cmpordsd",
        "cmppd",
        "cmpsd",
        "cmpunordpd",
        "cmpunordsd",
        "comisd",
        "cvtdq2pd",
        "cvtdq2ps",
        "cvtpd2dq",
        "cvtpd2pi",
        "cvtpd2ps",
        "cvtpi2pd",
        "cvtps2dq",
        "cvtps2pd",
        "cvtsd2si",
        "cvtsd2ss",
        "cvtsi2sd",
        "cvtss2sd",
        "cvttpd2dq",
        "cvttpd2pi",
        "cvttps2dq",
        "cvttsd2si",
        "divpd",
        "divsd",
        "lfence",
        "maskmovdqu",
        "maxpd",
        "maxsd",
        "mfence",
        "minpd",
        "minsd",
        "movapd",
        "movd",
        "movdq2q",
        "movdqa",
        "movdqu",
        "movhpd",
        "movlpd",
        "movmskpd",
        "movntdq",
        "movnti",
        "movntpd",
        "movq",
        "movq2dq",
        "movsd",
        "movupd",
        "mulpd",
        "mulsd",
        "orpd",
        "packssdw",
        "packsswb",
        "packuswb",
        "paddb",
        "paddd",
        "paddq",
        "paddsb",
        "paddsw",
        "paddusb",
        "paddusw",
        "paddw",
        "pand",
        "pandn",
        "pause",
        "pavgb",
        "pavgw",
        "pcmpeqb",
        "pcmpeqd",
        "pcmpeqw",
        "pcmpgtb",
        "pcmpgtd",
        "pcmpgtw",
        "pextrw",
        "pinsrw",
        "pmaddwd",
        "pmaxsw",
        "pmaxub",
        "pminsw",
        "pminub",
        "pmovmskb",
        "pmulhuw",
        "pmulhw",
        "pmullw",
        "pmuludq",
        "por",
        "psadbw",
        "pshufd",
        "pshufhw",
        "pshuflw",
        "pslld",
        "pslldq",
        "psllq",
        "psllw",
        "psrad",
        "psraw",
        "psrld",
        "psrldq",
        "psrlq",
        "psrlw",
        "psubb",
        "psubd",
        "psubq",
        "psubsb",
        "psubsw",
        "psubusb",
        "psubusw",
        "psubw",
        "punpckhbw",
        "punpckhdq",
        "punpckhqdq",
        "punpckhwd",
        "punpcklbw",
        "punpckldq",
    ],
    "sse3" => String[
        "addsubpd",
        "addsubps",
        "fisttp",
        "fisttpl",
        "fisttpll",
        "haddpd",
        "haddps",
        "hsubpd",
        "hsubps",
        "lddqu",
        "monitor",
        "movddup",
        "movshdup",
        "movsldup",
        "mwait",
    ],
    "sse4.1" => String[
        "blendpd",
        "blendps",
        "blendvpd",
        "blendvps",
        "dppd",
        "dpps",
        "extractps",
        "insertps",
        "movntdqa",
        "mpsadbw",
        "packusdw",
        "pblendvb",
        "pblendw",
        "pcmpeqq",
        "pextrb",
        "pextrd",
        "pextrq",
        "pextrw",
        "phminposuw",
        "pinsrb",
        "pinsrd",
        "pinsrq",
        "pmaxsb",
        "pmaxsd",
        "pmaxud",
        "pmaxuw",
        "pminsb",
        "pminsd",
        "pminud",
        "pminuw",
        "pmovsxbd",
        "pmovsxbq",
        "pmovsxbw",
        "pmovsxdq",
        "pmovsxwd",
        "pmovsxwq",
        "pmovzxbd",
        "pmovzxbq",
        "pmovzxbw",
        "pmovzxdq",
        "pmovzxwd",
        "pmovzxwq",
        "pmuldq",
        "pmulld",
        "ptest",
        "roundpd",
        "roundps",
        "roundsd",
        "roundss",
    ],
    "sse4.2" => String[
        "crc32",
        "pcmpestri",
        "pcmpestrm",
        "pcmpgtq",
        "pcmpistri",
        "pcmpistrm",
        "popcnt",
    ],
    "sse4a" => String[
        "extrq",
        "insertq",
        "lzcnt",
        "movntsd",
        "movntss",
        "popcnt",
    ],
    "ssse3" => String[
        "pabsb",
        "pabsw",
        "palignr",
        "phaddd",
        "phaddsw",
        "phaddw",
        "phsubd",
        "phsubsw",
        "phsubw",
        "pmaddubsw",
        "pmulhrsw",
        "pshufb",
        "psignb",
        "psignd",
        "psignw",
    ],
    "3dnow" => String[
        "femms",
        "pavgusb",
        "pf2id",
        "pf2iw",
        "pfacc",
        "pfadd",
        "pfcmpeq",
        "pfcmpge",
        "pfcmpgt",
        "pfmax",
        "pfmin",
        "pfmul",
        "pfrcp",
        "pfrcpit1",
        "pfrcpit2",
        "pfrsqit1",
        "pfrsqrt",
        "pfsub",
        "pfsubr",
        "pi2fd",
        "pi2fw",
        "pmulhrw",
        "prefetch",
        "pswapw",
    ],
    "3dnowext" => String[
        "maskmovq",
        "movntq",
        "pavgb",
        "pavgw",
        "pextrw",
        "pf2iw",
        "pfnacc",
        "pfpnacc",
        "pi2fw",
        "pinsrw",
        "pmaxsw",
        "pmaxub",
        "pminsw",
        "pminub",
        "pmovmskb",
        "pmulhuw",
        "prefetchnta",
        "prefetcht0",
        "prefetcht1",
        "prefetcht2",
        "psadbw",
        "pshufw",
        "pswapd",
        "sfence",
    ],
    "aes" => String[
        "aesdec",
        "aesdeclast",
        "aesenc",
        "aesenclast",
        "aesimc",
        "aeskeygenassist",
    ],
    "pclmul" => String[
        "pclmulqdq",
    ],
    "avx" => String[
        "vaddpd",
        "vaddps",
        "vaddsd",
        "vaddss",
        "vaddsubpd",
        "vaddsubps",
        "vaesdec",
        "vaesdeclast",
        "vaesenc",
        "vaesenclast",
        "vaesimc",
        "vaeskeygenassist",
        "vandnpd",
        "vandnps",
        "vandpd",
        "vandps",
        "vblendpd",
        "vblendps",
        "vblendvpd",
        "vblendvps",
        "vbroadcastf128",
        "vbroadcastsd",
        "vbroadcastss",
        "vcmpeq_osps",
        "vcmpeq_uqps",
        "vcmpeq_usps",
        "vcmpeqpd",
        "vcmpeqps",
        "vcmpeqsd",
        "vcmpeqss",
        "vcmpfalse_osps",
        "vcmpfalseps",
        "vcmpge_oqps",
        "vcmpgeps",
        "vcmpgt_oqps",
        "vcmpgtps",
        "vcmple_oqps",
        "vcmplepd",
        "vcmpleps",
        "vcmplesd",
        "vcmpless",
        "vcmplt_oqps",
        "vcmpltpd",
        "vcmpltps",
        "vcmpltsd",
        "vcmpltss",
        "vcmpneq_oqps",
        "vcmpneq_osps",
        "vcmpneq_usps",
        "vcmpneqpd",
        "vcmpneqps",
        "vcmpneqsd",
        "vcmpneqss",
        "vcmpnge_uqps",
        "vcmpngeps",
        "vcmpngt_uqps",
        "vcmpngtps",
        "vcmpnle_uqps",
        "vcmpnlepd",
        "vcmpnleps",
        "vcmpnlesd",
        "vcmpnless",
        "vcmpnlt_uqps",
        "vcmpnltpd",
        "vcmpnltps",
        "vcmpnltsd",
        "vcmpnltss",
        "vcmpord_sps",
        "vcmpordpd",
        "vcmpordps",
        "vcmpordsd",
        "vcmpordss",
        "vcmppd",
        "vcmpps",
        "vcmpsd",
        "vcmpss",
        "vcmptrue_usps",
        "vcmptrueps",
        "vcmpunord_sps",
        "vcmpunordpd",
        "vcmpunordps",
        "vcmpunordsd",
        "vcmpunordss",
        "vcomisd",
        "vcomiss",
        "vcvtdq2pd",
        "vcvtdq2ps",
        "vcvtpd2dq",
        "vcvtpd2dqx",
        "vcvtpd2dqy",
        "vcvtpd2ps",
        "vcvtpd2psx",
        "vcvtpd2psy",
        "vcvtps2dq",
        "vcvtps2pd",
        "vcvtsd2si",
        "vcvtsd2sil",
        "vcvtsd2ss",
        "vcvtsi2sd",
        "vcvtsi2sdl",
        "vcvtsi2ss",
        "vcvtsi2ssl",
        "vcvtss2sd",
        "vcvtss2si",
        "vcvtss2sil",
        "vcvttpd2dq",
        "vcvttpd2dqx",
        "vcvttpd2dqy",
        "vcvttps2dq",
        "vcvttsd2si",
        "vcvttss2si",
        "vdivpd",
        "vdivps",
        "vdivsd",
        "vdivss",
        "vdppd",
        "vdpps",
        "vextractf128",
        "vextractps",
        "vhaddpd",
        "vhaddps",
        "vhsubpd",
        "vhsubps",
        "vinsertf128",
        "vinsertps",
        "vlddqu",
        "vldmxcsr",
        "vmaskmovdqu",
        "vmaskmovpd",
        "vmaskmovps",
        "vmaxpd",
        "vmaxps",
        "vmaxsd",
        "vmaxss",
        "vminpd",
        "vminps",
        "vminsd",
        "vminss",
        "vmovapd",
        "vmovaps",
        "vmovd",
        "vmovddup",
        "vmovdqa",
        "vmovdqu",
        "vmovhlps",
        "vmovhpd",
        "vmovhps",
        "vmovlhps",
        "vmovlpd",
        "vmovlps",
        "vmovmskpd",
        "vmovmskps",
        "vmovntdq",
        "vmovntdqa",
        "vmovntpd",
        "vmovntps",
        "vmovq",
        "vmovsd",
        "vmovshdup",
        "vmovsldup",
        "vmovss",
        "vmovupd",
        "vmovups",
        "vmpsadbw",
        "vmulpd",
        "vmulps",
        "vmulsd",
        "vmulss",
        "vorpd",
        "vorps",
        "vpabsb",
        "vpabsd",
        "vpabsw",
        "vpackssdw",
        "vpacksswb",
        "vpackusdw",
        "vpackuswb",
        "vpaddb",
        "vpaddd",
        "vpaddq",
        "vpaddsb",
        "vpaddsw",
        "vpaddusb",
        "vpaddusw",
        "vpaddw",
        "vpalignr",
        "vpand",
        "vpandn",
        "vpavgb",
        "vpavgw",
        "vpblendvb",
        "vpblendw",
        "vpclmulhqhqdq",
        "vpclmulhqlqdq",
        "vpclmullqhqdq",
        "vpclmullqlqdq",
        "vpclmulqdq",
        "vpcmpeqb",
        "vpcmpeqd",
        "vpcmpeqq",
        "vpcmpeqw",
        "vpcmpestri",
        "vpcmpestrm",
        "vpcmpgtb",
        "vpcmpgtd",
        "vpcmpgtq",
        "vpcmpgtw",
        "vpcmpistri",
        "vpcmpistrm",
        "vperm2f128",
        "vpermilpd",
        "vpermilps",
        "vpextrb",
        "vpextrd",
        "vpextrw",
        "vphaddd",
        "vphaddsw",
        "vphaddw",
        "vphminposuw",
        "vphsubd",
        "vphsubsw",
        "vphsubw",
        "vpinsrb",
        "vpinsrd",
        "vpinsrw",
        "vpmaddubsw",
        "vpmaxsb",
        "vpmaxsd",
        "vpmaxsw",
        "vpmaxub",
        "vpmaxud",
        "vpmaxuw",
        "vpminsb",
        "vpminsd",
        "vpminsw",
        "vpminub",
        "vpminud",
        "vpminuw",
        "vpmovmskb",
        "vpmovsxbd",
        "vpmovsxbq",
        "vpmovsxbw",
        "vpmovsxdq",
        "vpmovsxwd",
        "vpmovsxwq",
        "vpmovzxbd",
        "vpmovzxbq",
        "vpmovzxbw",
        "vpmovzxdq",
        "vpmovzxwd",
        "vpmovzxwq",
        "vpmuldq",
        "vpmulhrsw",
        "vpmulhuw",
        "vpmulhw",
        "vpmulld",
        "vpmullw",
        "vpmuludq",
        "vpor",
        "vpsadbw",
        "vpshufb",
        "vpshufd",
        "vpshufhw",
        "vpshuflw",
        "vpsignb",
        "vpsignd",
        "vpsignw",
        "vpslld",
        "vpslldq",
        "vpsllq",
        "vpsllw",
        "vpsrad",
        "vpsraw",
        "vpsrld",
        "vpsrldq",
        "vpsrlq",
        "vpsrlw",
        "vpsubb",
        "vpsubd",
        "vpsubq",
        "vpsubsb",
        "vpsubsw",
        "vpsubusb",
        "vpsubusw",
        "vpsubw",
        "vptest",
        "vpunpckhbw",
        "vpunpckhdq",
        "vpunpckhqdq",
        "vpunpckhwd",
        "vpunpcklbw",
        "vpunpckldq",
        "vpunpcklqdq",
        "vpunpcklwd",
        "vpxor",
        "vrcpps",
        "vrcpss",
        "vroundpd",
        "vroundps",
        "vroundsd",
        "vroundss",
        "vrsqrtps",
        "vrsqrtss",
        "vshufpd",
        "vshufps",
        "vsqrtpd",
        "vsqrtps",
        "vsqrtsd",
        "vsqrtss",
        "vstmxcsr",
        "vsubpd",
        "vsubps",
        "vsubsd",
        "vsubss",
        "vtestpd",
        "vtestps",
        "vucomisd",
        "vucomiss",
        "vunpckhpd",
        "vunpckhps",
        "vunpcklpd",
        "vunpcklps",
        "vxorpd",
        "vxorps",
        "vzeroall",
        "vzeroupper",
    ],
    "avx2" => [
        "vblendmpd",
        "vblendmps",
        "vpblendmd",
        "vpblendmq",
        "vpblendmb",
        "vpblendmw",
        "vpcmpd",
        "vpcmpud",
        "vpcmpq",
        "vpcmpuq",
        "vpcmpb",
        "vpcmpub",
        "vpcmpw",
        "vpcmpuw",
        "vptestmq",
        "vptestmd",
        "vptestnmd",
        "vptestnmq",
        "vptestmb",
        "vptestmw",
        "vptestnmb",
        "vptestnmw",
        "vcompresspd",
        "vcompressps",
        "vpcompressd",
        "vpcompressq",
        "vexpandpd",
        "vexpandps",
        "vpexpandd",
        "vpexpandq",
        "vpermb",
        "vpermw",
        "vpermt2b",
        "vpermt2w",
        "vpermi2pd",
        "vpermi2ps",
        "vpermi2d",
        "vpermi2q",
        "vpermi2b",
        "vpermi2w",
        "vpermt2ps",
        "vpermt2pd",
        "vpermt2d",
        "vpermt2q",
        "vshuff32x4",
        "vshuff64x2",
        "vshuffi32x4",
        "vshuffi64x2",
        "vpternlogd",
        "vpternlogq",
        "vpmovqd",
        "vpmovsqd",
        "vpmovusqd",
        "vpmovqw",
        "vpmovsqw",
        "vpmovusqw",
        "vpmovqb",
        "vpmovsqb",
        "vpmovusqb",
        "vpmovdw",
        "vpmovsdw",
        "vpmovusdw",
        "vpmovdb",
        "vpmovsdb",
        "vpmovusdb",
        "vpmovwb",
        "vpmovswb",
        "vpmovuswb",
        "vcvtps2udq",
        "vcvtpd2udq",
        "vcvttps2udq",
        "vcvttpd2udq",
        "vcvtss2usi ",
        "vcvtsd2usi ",
        "vcvttss2usi ",
        "vcvttsd2usi",
        "vcvtps2qq",
        "vcvtpd2qq",
        "vcvtps2uqq",
        "vcvtpd2uqq",
        "vcvttps2qq",
        "vcvttpd2qq",
        "vcvttps2uqq",
        "vcvttpd2uqq",
        "vcvtudq2ps ",
        "vcvtudq2pd",
        "vcvtusi2ps ",
        "vcvtusi2pd",
        "vcvtusi2sd",
        "vcvtusi2ss",
        "vcvtuqq2ps",
        "vcvtuqq2pd",
        "vcvtqq2pd",
        "vcvtqq2ps",
        "vgetexppd",
        "vgetexpps",
        "vgetexpsd",
        "vgetexpss",
        "vgetmantpd",
        "vgetmantps",
        "vgetmantsd",
        "vgetmantss",
        "vfixupimmpd",
        "vfixupimmps",
        "vfixupimmsd",
        "vfixupimmss",
        "vrcp14pd",
        "vrcp14ps",
        "vrcp14sd",
        "vrcp14ss",
        "vrndscaleps",
        "vrndscalepd",
        "vrndscaless",
        "vrndscalesd",
        "vrsqrt14pd",
        "vrsqrt14ps",
        "vrsqrt14sd",
        "vrsqrt14ss",
        "vscalefps",
        "vscalefpd",
        "vscalefss",
        "vscalefsd",
        "valignd",
        "valignq",
        "vdbpsadbw",
        "vpbroadcastd",
        "vpbroadcastq",
        "vpbroadcastb",
        "vpbroadcastw",
        "vbroadcastf32x4",
        "vbroadcastf64x4",
        "vbroadcasti32x4",
        "vbroadcasti64x4",
        "vbroadcastf32x2",
        "vbroadcastf64x2",
        "vbroadcastf32x8",
        "vbroadcasti32x2",
        "vbroadcasti64x2",
        "vbroadcasti32x8",
        "vextractf32x4",
        "vextractf64x4",
        "vextracti32x4",
        "vextracti64x4",
        "vextractf64x2",
        "vextractf32x8",
        "vextracti64x2",
        "vextracti32x8",
        "vinsertf32x4",
        "vinsertf64x4",
        "vinserti32x4",
        "vinserti64x4",
        "vinsertf64x2",
        "vinsertf32x8",
        "vinserti64x2",
        "vinserti32x8",
        "vpabsq",
        "vpmaxsq",
        "vpmaxuq",
        "vpminsq",
        "vpminuq",
        "vpmultishiftqb",
        "vprold",
        "vprolvd",
        "vprolq",
        "vprolvq",
        "vprord",
        "vprorvd",
        "vprorq",
        "vprorvq",
        "vpscatterdd",
        "vpscatterdq",
        "vpscatterqd",
        "vpscatterqq",
        "vscatterdps",
        "vscatterdpd",
        "vscatterqps",
        "vscatterqpd",
        "vmovdqa32",
        "vmovdqa64",
        "vmovdqu8",
        "vmovdqu16",
        "vmovdqu32",
        "vmovdqu64",
        "vpxord",
        "vpxorq",
        "vpconflictd",
        "vpconflictq",
        "vplzcntd",
        "vplzcntq",
        "vpbroadcastmb2q",
        "pbroadcastmw2d",
        "vexp2pd",
        "vexp2ps",
        "vrcp28pd",
        "vrcp28ps",
        "vrcp28sd",
        "vrcp28ss",
        "vrsqrt28pd",
        "vrsqrt28ps",
        "vrsqrt28sd",
        "vrsqrt28ss",
        "vgatherpf0dps",
        "vgatherpf0qps",
        "vgatherpf0dpd",
        "vgatherpf0qpd",
        "vgatherpf1dps",
        "vgatherpf1qps",
        "vgatherpf1dpd",
        "vgatherpf1qpd",
        "vscatterpf0dps",
        "vscatterpf0qps",
        "vscatterpf0dpd",
        "vscatterpf0qpd",
        "vscatterpf1dps",
        "vscatterpf1qps",
        "vscatterpf1dpd",
        "vscatterpf1qpd",
        "vfpclassps",
        "vfpclasspd",
        "vfpclassss",
        "vfpclasssd",
        "vrangeps",
        "vrangepd",
        "vrangess",
        "vrangesd",
        "vreduceps",
        "vreducepd",
        "vreducess",
        "vreducesd",
        "vpmovm2d",
        "vpmovm2q",
        "vpmovm2b",
        "vpmovm2w",
        "vpmovd2m",
        "vpmovq2m",
        "vpmovb2m",
        "vpmovw2m",
        "vpmullq",
        "vpmadd52luq",
        "vpmadd52huq",
        "v4fmaddps",
        "v4fmaddss",
        "v4fnmaddps",
        "v4fnmaddss",
        "vp4dpwssd",
        "vp4dpwssds",
        "vpopcntd",
        "vpopcntq",
    ],

)

# Turn instructions "inside out", so e.g. we have "vzeroall" => "avx"
mnemonics_by_category = Dict(
    inst => cat for (cat, insts) in instruction_categories for inst in insts
)

"""
    instruction_mnemonics(path::AbstractString, platform::Platform)

Dump a binary object with `objdump`, returning a list of instruction mnemonics
for further analysis with `analyze_instruction_set()`.

Note that this function only really makes sense for x86/x64 binaries.  Don't
run this on armv7l, aarch64, ppc64le etc... binaries and expect it to work.

This function returns the list of mnemonics as well as the counts of each,
binned by the mapping defined within `instruction_categories`.
"""
function instruction_mnemonics(path::AbstractString, platform::Platform)
    # The outputs we are calculating
    counts = Dict(k => 0 for k in keys(instruction_categories))
    mnemonics = String[]

    ur = preferred_runner()(
        abspath(dirname(path));
        cwd="/workspace/",
        platform=platform,
        verbose=false,
    )
    output = IOBuffer()

    # Run objdump to disassemble the input binary
    objdump_cmd = "\$OBJDUMP -d $(basename(path))"
    run_interactive(ur, `/bin/bash -c "$(objdump_cmd)"`; stdout=output)
    seekstart(output)

    for line in eachline(output)
        # First, ensure that this line of output is 3 fields long at least
        fields = split(line, '\t')
        if length(fields) < 3
            continue
        end

        # Grab the mnemonic for this line as the first word of the 3rd field
        m = split(fields[3])[1]
        push!(mnemonics, m)

        # For each mnemonic, find it in mnemonics_by_category, if we can, and
        # increment the appropriate `counts` member:
        if haskey(mnemonics_by_category, m)
            counts[mnemonics_by_category[m]] += 1
        else
            counts["unknown"] += 1
        end
    end

    # Return both the list of mnemonics as well as the binned counts
    return mnemonics, counts
end

"""
    minimum_instruction_set(counts::Dict, is_64bit::Bool)

This function returns the minimum instruction set required, depending on
whether the object file being pointed to is a 32-bit or 64-bit one:

* For 32-bit object files, this returns one of [:pentium4, :prescott]

* For 64-bit object files, this returns one of [:core2, :sandybridge, :haswell]
"""
function minimum_instruction_set(counts::Dict, is_64bit::Bool)
    if is_64bit
        if counts["avx2"] > 0
            return :haswell
        end
        if counts["avx"] > 0
            return :sandybridge
        end
        return :core2
    else
        if counts["sse3"] > 0
            return :prescott
        end
        return :pentium4
    end
end

"""
    analyze_instruction_set(oh::ObjectHandle, platform::Platform; verbose::Bool = false)

Analyze the instructions within the binary located at the given path for which
minimum instruction set it requires, taking note of groups of instruction sets
used such as `avx`, `sse4.2`, `i486`, etc....

Some binary files (such as libopenblas) contain multiple versions of functions,
internally determining which version to call by using the `cpuid` instruction
to determine processor support.  In an effort to detect this, we make note of
any usage of the `cpuid` instruction, disabling our minimum instruction set
calculations if such an instruction is found, and notifying the user of this
if `verbose` is set to `true`.

Note that this function only really makes sense for x86/x64 binaries.  Don't
run this on armv7l, aarch64, ppc64le etc... binaries and expect it to work.
"""
function analyze_instruction_set(oh::ObjectHandle, platform::Platform; verbose::Bool = false, io::IO = stdout)
    # Get list of mnemonics
    mnemonics, counts = instruction_mnemonics(path(oh), platform)

    # Analyze for minimum instruction set
    min_isa = minimum_instruction_set(counts, is64bit(oh))

    # If the binary uses `cpuid`, we can't know what it's doing, so just
    # return the most conservative ISA and warn the user if `verbose` is set.
    if counts["cpuid"] > 0
        new_min_isa = is64bit(oh) ? :core2 : :pentium4
        if verbose && new_min_isa != min_isa
            msg = replace("""
            $(basename(path(oh))) contains a `cpuid` instruction; refusing to
            analyze for minimum instruction set, as it may dynamically select
            the proper instruction set internally.  Would have chosen
            $(min_isa), instead choosing $(new_min_isa).
            """, '\n' => ' ')
            warn(io, strip(msg))
        end
        return new_min_isa
    end

    # Otherwise, return `min_isa` and let 'em know!
    return min_isa
end


text:
##########################################################################
###
### Fletcher-Powell
###
##########################################################################

function fletcher_powell(x::Vector)
    function theta(x::Vector)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end

    return 100.0 * (x[3] - 10.0 * theta(x[1], x[2]))^2 +
            (sqrt(x[1]^2 + x[2]^2) - 1.0)^2 + x[3]^2
end

function fletcher_powell_gradient(x::Vector, storage::Vector)
    error("Not yet defined")
end

function fletcher_powell_hessian(x::Vector, storage::Matrix)
    error("Not yet defined")
end

##########################################################################
###
### Himmelbrau's Function
###
##########################################################################

function himmelbrau(x::Vector)
    return (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
end

function himmelbrau_gradient!(x::Vector, storage::Vector)
    storage[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] +
                  44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    storage[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
                  4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
end

function himmelbrau_hessian!(x::Vector, storage::Matrix)
    error("Not yet defined")
end

##########################################################################
###
### Parabola
###
##########################################################################

function parabola(x::Vector)
    return (1.0 - x[1])^2 + (2.0 - x[2])^2 + (3.0 - x[3])^2 +
            (5.0 - x[4])^2 + (8.0 - x[5])^2
end

function parabola_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (1.0 - x[1])
    storage[2] = -2.0 * (2.0 - x[2])
    storage[3] = -2.0 * (3.0 - x[3])
    storage[4] = -2.0 * (5.0 - x[4])
    storage[5] = -2.0 * (8.0 - x[5])
end

function parabola_hessian!(x::Vector, storage::Matrix)
    for i in 1:5
        for j in 1:5
            if i == j
                storage[i, j] = 2.0
            else
                storage[i, j] = 0.0
            end
        end
    end
end

##########################################################################
###
### Powell
###
##########################################################################

function powell(x::Vector)
    return (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
            (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4
end

function powell_gradient!(x::Vector, storage::Vector)
    storage[1] = 2.0 * (x[1] + 10.0 * x[2]) + 40.0 * (x[1] - x[4])^3
    storage[2] = 20.0 * (x[1] + 10.0 * x[2]) + 4.0 * (x[2] - 2.0 * x[3])^3
    storage[3] = 10.0 * (x[3] - x[4]) - 8.0 * (x[2] - 2.0 * x[3])^3
    storage[4] = -10.0 * (x[3] - x[4]) - 40.0 * (x[1] - x[4])^3
end

function powell_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0 + 120.0 * (x[1] - x[4])^2
    storage[1, 2] = 20.0
    storage[1, 3] = 0.0
    storage[1, 4] = 0.0
    storage[2, 1] = 20.0
    storage[2, 2] = 200.0 + 12.0 * (x[2] - 2.0 * x[3])^2
    storage[2, 3] = -24.0 * (x[2] - 2.0 * x[3])^2
    storage[2, 4] = 0.0
    storage[3, 1] = 0.0
    storage[3, 2] = -24.0 * (x[2] - 2.0 * x[3])^2
    storage[3, 3] = 10.0 + 48.0 * (x[2] - 2.0 * x[3])^2
    storage[3, 4] = -10.0
    storage[4, 1] = 0.0
    storage[4, 2] = 0.0
    storage[4, 3] = -10.0
    storage[4, 4] = 10.0 + 120.0 * (x[1] - x[4])^2
end

##########################################################################
###
### Rosenbrock
###
##########################################################################

function rosenbrock(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

function rosenbrock_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    storage[2] = 200.0 * (x[2] - x[1]^2)
end

function rosenbrock_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
    storage[1, 2] = -400.0 * x[1]
    storage[2, 1] = -400.0 * x[1]
    storage[2, 2] = 200.0
end

##########################################################################
###
### Simple 4th-Degree Polynomial Example
###
##########################################################################

function polynomial(x::Vector)
    return (10.0 - x[1])^2 + (7.0 - x[2])^4 + (108.0 - x[3])^4
end

function polynomial_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (10.0 - x[1])
    storage[2] = -4.0 * (7.0 - x[2])^3
    storage[3] = -4.0 * (108.0 - x[3])^3
end

function polynomial_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0
    storage[1, 2] = 0.0
    storage[1, 3] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = 12.0 * (7.0 - x[2])^2
    storage[2, 3] = 0.0
    storage[3, 1] = 0.0
    storage[3, 2] = 0.0
    storage[3, 3] = 12.0 * (108.0 - x[3])^2
end

##########################################################################
###
### Simple Exponential Function Example
###
##########################################################################

function exponential(x::Vector)
    return exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
end

function exponential_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
    storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
end

function exponential_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    storage[1, 2] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = 2.0 * exp((3.0 - x[1])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
end

##########################################################################
###
### Large-Scale Quadratic
###
##########################################################################

function large_polynomial(x::Vector)
    res = 0.0
    for i in 1:250
        res += (i - x[i])^2
    end
    return res
end

function large_polynomial_gradient!(x::Vector, storage::Vector)
    for i in 1:250
        storage[i] = -2.0 * (i - x[i])
    end
end

function large_polynomial_hessian!(x::Vector, storage::Matrix)
    for i in 1:250
        for j in 1:250
            if i == j
                storage[i, j] = 2.0
            else
                storage[i, j] = 0.0
            end
        end
    end
end


text:
using Test
using PressureFieldContact
using StaticArrays
using Rotations
using LinearAlgebra
using RigidBodyDynamics
using PressureFieldContact.Clip
using PressureFieldContact.MathKernel
using PressureFieldContact.Binary_BB_Trees
using PressureFieldContact.Radau
using PressureFieldContact.Geometry

set_zero_subnormals(true)
BLAS.set_num_threads(1)  # NOTE: comment out this line if using IntelMKL

include("test_friction.jl")
include(joinpath("test_radau", "runtests.jl"))
include(joinpath("test_geometry", "runtests.jl"))
include(joinpath("test_obb", "runtests.jl"))
include(joinpath("test_math_kernel", "runtests.jl"))
include(joinpath("test_clip", "runtests.jl"))
include("test_exports.jl")
include("test_normal.jl")
include("test_mesh_body_utility.jl")
include("test_vol_vol.jl")
include("boxes.jl")


text:
using Distances
using HDF5
using Shrike
using NearestNeighbors
using NPZ

#DATAFOLDER = "/pine/scr/d/j/djpassey/rpdata/"
#NAMES = ["fashionmnist", "mnist", "sift", "gist"]
#FNAMES = ["fashion-mnist-784-euclidean.hdf5", "mnist-784-euclidean.hdf5", "sift-128-euclidean.hdf5", "gist-960-euclidean.hdf5"]
DATAFOLDER = "/Users/djpassey/Downloads/"
NAMES = ["mist"]
FNAMES = ["mnist-784-euclidean.hdf5"]

println("Loading Datasets")
FILES = map(x -> string(DATAFOLDER, x), FNAMES)
SHRIKE_NTREES = [2]#, 5, 8, 10, 20, 100, 200]
SHRIKE_DEPTH = [4]#, 5, 7, 8, 9, 11, 12, 13]
SHRIKE_NVOTES = [1]#, 2, 3, 10, 20]
KDTREE_LEAFSIZE = [100]#, 300, 500, 1000, 2000, 5000, 10000]

struct NNData
    name::String
    fname::String
    nfeatures::Int
    npoints::Int
    ntest::Int
    k::Int
end

function NNData(name::String, fname::String)
    data = read(h5open(fname))
    nfeatures, npoints = size(data["train"])
    k = size(data["neighbors"], 1)
    ntest, _ = size(data["test"])
    return NNData(name, fname, nfeatures, npoints, ntest, k)
end

train(d::NNData) = read(h5open(d.fname))["train"]
test(d::NNData) = read(h5open(d.fname))["test"]
neighbors(d::NNData) = read(h5open(d.fname))["neighbors"]

brutenn(x::AbstractArray, X::AbstractArray, k::Int) =
    sortperm(map(i -> sqeuclidean(x, X[:,i]), 1:size(X,2)), alg=PartialQuickSort(k))[1:k]

allNNDatasets() = map(NNData, NAMES, FILES)

println("Datasets Loaded")
for nndata in allNNDatasets()
    println("Starting Dataset: $(nndata.name)")
    X = train(nndata)
    Xts = test(nndata)
    println("Data shape: $(size(X))")
    rpt_results::Array{Array{Float64, 1}, 1} = [[],[]]
    kdt_results::Array{Array{Float64, 1}, 1} = [[],[]]
    maxdepth = floor(log(2, nndata.npoints) - log(2, nndata.k))
    # Shrike Algorithm
    for ntrees in SHRIKE_NTREES
        for depth in SHRIKE_DEPTH[SHRIKE_DEPTH .< maxdepth]
            for nvotes in SHRIKE_NVOTES
                index = ShrikeIndex(X, nndata.k, depth=depth, ntrees=ntrees)
                tot_time = 0.0
                tot_recall = 0.0
                for i in 1:nndata.ntest
                    q = Xts[:, i]
                    tnn = brutenn(q, X, nndata.k)
                    td = @timed approx_nn = ann(index, q, nndata.k, vote_cutoff=nvotes)
                    tot_time += td.time
                    tot_recall += (length(intersect(td.value, tnn)) / nndata.k)
                end
                push!(rpt_results[1], nndata.ntest / tot_time)
                push!(rpt_results[2], tot_recall / nndata.ntest)
            end
        end
    end
    npzwrite("rpt-$(nndata.name).npy", hcat(rpt_results...))
    println("Shrike Complete")

    for leafsize in KDTREE_LEAFSIZE
        index = KDTree(X, leafsize=leafsize)
        tot_time = 0.0
        tot_recall = 0.0
        for i in 1:nndata.ntest
            q = Xts[:, i]
            tnn =  brutenn(q, X, nndata.k)
            td = @timed approxnn = knn(index, q, nndata.k)
            tot_time += td.time
            tot_recall += (length(intersect(td.value[1], tnn)) / nndata.k)
        end
        push!(kdt_results[1], nndata.ntest / tot_time)
        push!(kdt_results[2], tot_recall / nndata.ntest)
    end
    npzwrite("kdt-$(nndata.name).npy", hcat(kdt_results...))
    println("NearestNeighbors Complete")
end


text:
#=
 Copyright 2013 - 2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015 Oliver Heimlich (oheim@posteo.de)

 Original author: Marco Nehmeier (unit tests in libieeep1788)
 Converted into portable ITL format by Oliver Heimlich with minor corrections.

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
=#
#Language imports

#Test library imports
if VERSION < v"0.7.0-DEV.2004"
    using Base.Test
else
    using Test
end

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal_inf_test" begin
    @test inf(∅) == Inf
    @test inf(Interval(-Inf, Inf)) == -Inf
    @test inf(Interval(1.0, 2.0)) == 1.0
    @test inf(Interval(-3.0, -2.0)) == -3.0
    @test inf(Interval(-Inf, 2.0)) == -Inf
    @test inf(Interval(-Inf, 0.0)) == -Inf
    @test inf(Interval(-Inf, -0.0)) == -Inf
    @test inf(Interval(-2.0, Inf)) == -2.0
    @test inf(Interval(0.0, Inf)) == -0.0
    @test inf(Interval(-0.0, Inf)) == -0.0
    @test inf(Interval(-0.0, 0.0)) == -0.0
    @test inf(Interval(0.0, -0.0)) == -0.0
    @test inf(Interval(0.0, 0.0)) == -0.0
    @test inf(Interval(-0.0, -0.0)) == -0.0
end

@testset "minimal_inf_dec_test" begin
    @test inf(DecoratedInterval(∅, trv)) == Inf
    @test inf(DecoratedInterval(Interval(-Inf, Inf), def)) == -Inf
    @test inf(DecoratedInterval(Interval(1.0, 2.0), com)) == 1.0
    @test inf(DecoratedInterval(Interval(-3.0, -2.0), trv)) == -3.0
    @test inf(DecoratedInterval(Interval(-Inf, 2.0), dac)) == -Inf
    @test inf(DecoratedInterval(Interval(-Inf, 0.0), def)) == -Inf
    @test inf(DecoratedInterval(Interval(-Inf, -0.0), trv)) == -Inf
    @test inf(DecoratedInterval(Interval(-2.0, Inf), trv)) == -2.0
    @test inf(DecoratedInterval(Interval(0.0, Inf), def)) == -0.0
    @test inf(DecoratedInterval(Interval(-0.0, Inf), trv)) == -0.0
    @test inf(DecoratedInterval(Interval(-0.0, 0.0), dac)) == -0.0
    @test inf(DecoratedInterval(Interval(0.0, -0.0), trv)) == -0.0
    @test inf(DecoratedInterval(Interval(0.0, 0.0), trv)) == -0.0
    @test inf(DecoratedInterval(Interval(-0.0, -0.0), trv)) == -0.0
end

@testset "minimal_sup_test" begin
    @test sup(∅) == -Inf
    @test sup(Interval(-Inf, Inf)) == Inf
    @test sup(Interval(1.0, 2.0)) == 2.0
    @test sup(Interval(-3.0, -2.0)) == -2.0
    @test sup(Interval(-Inf, 2.0)) == 2.0
    @test sup(Interval(-Inf, 0.0)) == 0.0
    @test sup(Interval(-Inf, -0.0)) == 0.0
    @test sup(Interval(-2.0, Inf)) == Inf
    @test sup(Interval(0.0, Inf)) == Inf
    @test sup(Interval(-0.0, Inf)) == Inf
    @test sup(Interval(-0.0, 0.0)) == 0.0
    @test sup(Interval(0.0, -0.0)) == 0.0
    @test sup(Interval(0.0, 0.0)) == 0.0
    @test sup(Interval(-0.0, -0.0)) == 0.0
end

@testset "minimal_sup_dec_test" begin
    @test sup(DecoratedInterval(∅, trv)) == -Inf
    @test sup(DecoratedInterval(Interval(-Inf, Inf), def)) == Inf
    @test sup(DecoratedInterval(Interval(1.0, 2.0), com)) == 2.0
    @test sup(DecoratedInterval(Interval(-3.0, -2.0), trv)) == -2.0
    @test sup(DecoratedInterval(Interval(-Inf, 2.0), dac)) == 2.0
    @test sup(DecoratedInterval(Interval(-Inf, 0.0), def)) == 0.0
    @test sup(DecoratedInterval(Interval(-Inf, -0.0), trv)) == 0.0
    @test sup(DecoratedInterval(Interval(-2.0, Inf), trv)) == Inf
    @test sup(DecoratedInterval(Interval(0.0, Inf), def)) == Inf
    @test sup(DecoratedInterval(Interval(-0.0, Inf), trv)) == Inf
    @test sup(DecoratedInterval(Interval(-0.0, 0.0), dac)) == +0.0
    @test sup(DecoratedInterval(Interval(0.0, -0.0), trv)) == +0.0
    @test sup(DecoratedInterval(Interval(0.0, 0.0), trv)) == +0.0
    @test sup(DecoratedInterval(Interval(-0.0, -0.0), trv)) == +0.0
end

@testset "minimal_mid_test" begin
    @test mid(Interval(-Inf, Inf)) == 0.0
    @test mid(Interval(-0x1.fffffffffffffp1023, +0x1.fffffffffffffp1023)) == 0.0
    @test mid(Interval(0.0, 2.0)) == 1.0
    @test mid(Interval(2.0, 2.0)) == 2.0
    @test mid(Interval(-2.0, 2.0)) == 0.0
    @test mid(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022)) == 0.0
    @test mid(Interval(-0x0.0000000000001p-1022, 0x0.0000000000002p-1022)) == 0.0
end

@testset "minimal_mid_dec_test" begin
    @test mid(DecoratedInterval(Interval(-Inf, Inf), def)) == 0.0
    @test mid(DecoratedInterval(Interval(-0x1.fffffffffffffp1023, +0x1.fffffffffffffp1023), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(0.0, 2.0), com)) == 1.0
    @test mid(DecoratedInterval(Interval(2.0, 2.0), dac)) == 2.0
    @test mid(DecoratedInterval(Interval(-2.0, 2.0), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(-0x0.0000000000001p-1022, 0x0.0000000000002p-1022), trv)) == 0.0
end

@testset "minimal_rad_test" begin
    @test radius(Interval(0.0, 2.0)) == 1.0
    @test radius(Interval(2.0, 2.0)) == 0.0
    @test radius(Interval(-Inf, Inf)) == Inf
    @test radius(Interval(0.0, Inf)) == Inf
    @test radius(Interval(-Inf, 1.2)) == Inf
end

@testset "minimal_rad_dec_test" begin
    @test radius(DecoratedInterval(Interval(0.0, 2.0), trv)) == 1.0
    @test radius(DecoratedInterval(Interval(2.0, 2.0), com)) == 0.0
    @test radius(DecoratedInterval(Interval(-Inf, Inf), trv)) == Inf
    @test radius(DecoratedInterval(Interval(0.0, Inf), def)) == Inf
    @test radius(DecoratedInterval(Interval(-Inf, 1.2), trv)) == Inf
end

@testset "minimal_wid_test" begin
    @test diam(Interval(2.0, 2.0)) == 0.0
    @test diam(Interval(1.0, 2.0)) == 1.0
    @test diam(Interval(1.0, Inf)) == Inf
    @test diam(Interval(-Inf, 2.0)) == Inf
    @test diam(Interval(-Inf, Inf)) == Inf
end

@testset "minimal_wid_dec_test" begin
    @test diam(DecoratedInterval(Interval(2.0, 2.0), com)) == 0.0
    @test diam(DecoratedInterval(Interval(1.0, 2.0), trv)) == 1.0
    @test diam(DecoratedInterval(Interval(1.0, Inf), trv)) == Inf
    @test diam(DecoratedInterval(Interval(-Inf, 2.0), def)) == Inf
    @test diam(DecoratedInterval(Interval(-Inf, Inf), trv)) == Inf
end

@testset "minimal_mag_test" begin
    @test mag(Interval(1.0, 2.0)) == 2.0
    @test mag(Interval(-4.0, 2.0)) == 4.0
    @test mag(Interval(-Inf, 2.0)) == Inf
    @test mag(Interval(1.0, Inf)) == Inf
    @test mag(Interval(-Inf, Inf)) == Inf
    @test mag(Interval(-0.0, 0.0)) == 0.0
    @test mag(Interval(-0.0, -0.0)) == 0.0
end

@testset "minimal_mag_dec_test" begin
    @test mag(DecoratedInterval(Interval(1.0, 2.0), com)) == 2.0
    @test mag(DecoratedInterval(Interval(-4.0, 2.0), trv)) == 4.0
    @test mag(DecoratedInterval(Interval(-Inf, 2.0), trv)) == Inf
    @test mag(DecoratedInterval(Interval(1.0, Inf), def)) == Inf
    @test mag(DecoratedInterval(Interval(-Inf, Inf), trv)) == Inf
    @test mag(DecoratedInterval(Interval(-0.0, 0.0), trv)) == 0.0
    @test mag(DecoratedInterval(Interval(-0.0, -0.0), trv)) == 0.0
end

@testset "minimal_mig_test" begin
    @test mig(Interval(1.0, 2.0)) == 1.0
    @test mig(Interval(-4.0, 2.0)) == 0.0
    @test mig(Interval(-4.0, -2.0)) == 2.0
    @test mig(Interval(-Inf, 2.0)) == 0.0
    @test mig(Interval(-Inf, -2.0)) == 2.0
    @test mig(Interval(-1.0, Inf)) == 0.0
    @test mig(Interval(1.0, Inf)) == 1.0
    @test mig(Interval(-Inf, Inf)) == 0.0
    @test mig(Interval(-0.0, 0.0)) == 0.0
    @test mig(Interval(-0.0, -0.0)) == 0.0
end

@testset "minimal_mig_dec_test" begin
    @test mig(DecoratedInterval(Interval(1.0, 2.0), com)) == 1.0
    @test mig(DecoratedInterval(Interval(-4.0, 2.0), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(-4.0, -2.0), trv)) == 2.0
    @test mig(DecoratedInterval(Interval(-Inf, 2.0), def)) == 0.0
    @test mig(DecoratedInterval(Interval(-Inf, -2.0), trv)) == 2.0
    @test mig(DecoratedInterval(Interval(-1.0, Inf), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(1.0, Inf), trv)) == 1.0
    @test mig(DecoratedInterval(Interval(-Inf, Inf), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(-0.0, 0.0), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(-0.0, -0.0), trv)) == 0.0
end
# FactCheck.exitstatus()


text:
function main(window)
    push!(window.assets, "codemirror")
    inp = Signal(Dict())

    s = sampler() # A thing that lets you watch widgets/behaviors upon updates to other behaviors

    editor = watch!(s, :code, codemirror("some code"))
    code_cell = trigger!(s, :submit, keypress("ctrl+enter shift+enter", editor))

    vbox(
        intent(s, code_cell) >>> inp,
        inp
    )
end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 3.1.
#

fl = Float64

# define the step
δ = fl(2^(-52))

# we can modify where we want to start “counting”
starting_point = fl(0)
x = fl(1)

# print out the consecutive numbers
for i in range(starting_point, stop=2^52)
    println(
        bitstring(
            fl(x + fl(δ * i))
        ),
        " (",
        fl(x + fl(δ * i)),
        ")"
    )
    sleep(1)
end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
using RobustAdaptiveMetropolisSampler
using Distributions
using Statistics
using LinearAlgebra
using Random
using Test

@testset "RobustAdaptiveMetropolisSampler" begin

    Random.seed!(10)

    result = RAM_sample(p -> logpdf(Normal(3., 2.), p[1]), [0.], 0.5, 100_000, show_progress=false, output_log_probability_x=true)

    @test length(result.chain) == 100_000
    @test result.acceptance_rate ≈ 0.234 atol = 0.01
    @test mean(result.chain[:,1]) ≈ 3. atol = 0.1
    @test std(result.chain[:,1]) ≈ 2. atol = 0.01
    @test result.log_probabilities_x[1] == logpdf(Normal(3., 2.), result.chain[1])
    @test result.log_probabilities_x[length(result.chain)] == logpdf(Normal(3., 2.), result.chain[length(result.chain)])

    result = RAM_sample(p -> logpdf(Normal(3., 2.), p[1]), [0.], [0.5], 1_000, show_progress=false, output_log_probability_x=false)
    @test length(result.chain) == 1_000

    result = RAM_sample(p -> logpdf(Normal(3., 2.), p[1]), [0.], Diagonal([0.5]), 1_000, show_progress=false, output_log_probability_x=false)
    @test length(result.chain) == 1_000

    result = RAM_sample(p -> logpdf(Normal(3., 2.), p[1]), [0.], fill(0.5, 1, 1), 1_000, show_progress=false, output_log_probability_x=false)
    @test length(result.chain) == 1_000

end


text:
module Camera

using LinearAlgebra
using GeometryTypes
using Rotations
using ..LavaCore
using ..Utils

include("CameraBase.jl")
include("GenericCamera.jl")

end


text:
using Oceananigans.Grids: total_extent, halo_size

#####
##### Regular rectilinear grids
#####

function test_regular_rectilinear_correct_size(FT)
    grid = RegularRectilinearGrid(FT, size=(4, 6, 8), extent=(2π, 4π, 9π))

    @test grid.Nx == 4
    @test grid.Ny == 6
    @test grid.Nz == 8

    # Checking ≈ as the grid could be storing Float32 values.
    @test grid.Lx ≈ 2π
    @test grid.Ly ≈ 4π
    @test grid.Lz ≈ 9π

    return nothing
end

function test_regular_rectilinear_correct_extent(FT)
    grid = RegularRectilinearGrid(FT, size=(4, 6, 8), x=(1, 2), y=(π, 3π), z=(0, 4))

    @test grid.Lx ≈ 1
    @test grid.Ly ≈ 2π
    @test grid.Lz ≈ 4

    return nothing
end

function test_regular_rectilinear_correct_coordinate_lengths(FT)
    grid = RegularRectilinearGrid(FT, size=(2, 3, 4), extent=(1, 1, 1), halo=(1, 1, 1),
                                  topology=(Periodic, Bounded, Bounded))

    Nx, Ny, Nz = size(grid)
    Hx, Hy, Hz = halo_size(grid)

    @test length(grid.xC) == Nx + 2Hx
    @test length(grid.yC) == Ny + 2Hy
    @test length(grid.zC) == Nz + 2Hz
    @test length(grid.xF) == Nx + 2Hx
    @test length(grid.yF) == Ny + 2Hy + 1
    @test length(grid.zF) == Nz + 2Hz + 1

    return nothing
end

function test_regular_rectilinear_correct_halo_size(FT)
    grid = RegularRectilinearGrid(FT, size=(4, 6, 8), extent=(2π, 4π, 9π), halo=(1, 2, 3))

    @test grid.Hx == 1
    @test grid.Hy == 2
    @test grid.Hz == 3

    return nothing
end

function test_regular_rectilinear_correct_halo_faces(FT)
    N = 4
    H = 1
    L = 2.0
    Δ = L / N

    topo = (Periodic, Bounded, Bounded)
    grid = RegularRectilinearGrid(FT, topology=topo, size=(N, N, N), x=(0, L), y=(0, L), z=(0, L), halo=(H, H, H))

    @test grid.xF[0] == - H * Δ
    @test grid.yF[0] == - H * Δ
    @test grid.zF[0] == - H * Δ

    @test grid.xF[N+1] == L  # Periodic
    @test grid.yF[N+2] == L + H * Δ
    @test grid.zF[N+2] == L + H * Δ

    return nothing
end

function test_regular_rectilinear_correct_first_cells(FT)
    N = 4
    H = 1
    L = 4.0
    Δ = L / N

    grid = RegularRectilinearGrid(FT, size=(N, N, N), x=(0, L), y=(0, L), z=(0, L), halo=(H, H, H))

    @test grid.xC[1] == Δ/2
    @test grid.yC[1] == Δ/2
    @test grid.zC[1] == Δ/2

    return nothing
end

function test_regular_rectilinear_correct_end_faces(FT)
    N = 4
    L = 2.0
    Δ = L / N

    grid = RegularRectilinearGrid(FT, size=(N, N, N), x=(0, L), y=(0, L), z=(0, L), halo=(1, 1, 1),
                                  topology=(Periodic, Bounded, Bounded))

    @test grid.xF[N+1] == L
    @test grid.yF[N+2] == L + Δ
    @test grid.zF[N+2] == L + Δ

    return nothing
end

function test_regular_rectilinear_ranges_have_correct_length(FT)
    Nx, Ny, Nz = 8, 9, 10
    Hx, Hy, Hz = 1, 2, 1

    grid = RegularRectilinearGrid(FT, size=(Nx, Ny, Nz), extent=(1, 1, 1), halo=(Hx, Hy, Hz),
                                  topology=(Bounded, Bounded, Bounded))

    @test length(grid.xC) == Nx + 2Hx
    @test length(grid.yC) == Ny + 2Hy
    @test length(grid.zC) == Nz + 2Hz
    @test length(grid.xF) == Nx + 1 + 2Hx
    @test length(grid.yF) == Ny + 1 + 2Hy
    @test length(grid.zF) == Nz + 1 + 2Hz

    return nothing
end

# See: https://github.com/climate-machine/Oceananigans.jl/issues/480
function test_regular_rectilinear_no_roundoff_error_in_ranges(FT)
    Nx = Ny = 1
    Nz = 64
    Hz = 1

    grid = RegularRectilinearGrid(FT, size=(Nx, Ny, Nz), extent=(1, 1, π/2), halo=(1, 1, Hz))

    @test length(grid.zC) == Nz + 2Hz
    @test length(grid.zF) == Nz + 2Hz + 1

    return nothing
end

function test_regular_rectilinear_grid_properties_are_same_type(FT)
    grid = RegularRectilinearGrid(FT, size=(10, 10, 10), extent=(1, 1//7, 2π))

    @test grid.Lx isa FT
    @test grid.Ly isa FT
    @test grid.Lz isa FT
    @test grid.Δx isa FT
    @test grid.Δy isa FT
    @test grid.Δz isa FT

    @test eltype(grid.xF) == FT
    @test eltype(grid.yF) == FT
    @test eltype(grid.zF) == FT
    @test eltype(grid.xC) == FT
    @test eltype(grid.yC) == FT
    @test eltype(grid.zC) == FT

    return nothing
end

function test_xnode_ynode_znode_are_correct(FT)
    N = 3
    grid = RegularRectilinearGrid(FT, size=(N, N, N), x=(0, π), y=(0, π), z=(0, π),
                                  topology=(Periodic, Periodic, Bounded))

    @test xnode(Center(), 2, grid) ≈ FT(π/2)
    @test ynode(Center(), 2, grid) ≈ FT(π/2)
    @test znode(Center(), 2, grid) ≈ FT(π/2)

    @test xnode(Face(), 2, grid) ≈ FT(π/3)
    @test ynode(Face(), 2, grid) ≈ FT(π/3)
    @test znode(Face(), 2, grid) ≈ FT(π/3)

    return nothing
end

function test_regular_rectilinear_constructor_errors(FT)
    @test isbitstype(typeof(RegularRectilinearGrid(FT, size=(16, 16, 16), extent=(1, 1, 1))))

    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32,), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 64), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 32, 32, 16), extent=(1, 1, 1))

    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 32, 32.0), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(20.1, 32, 32), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, nothing, 32), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, "32", 32), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 32, 32), extent=(1, nothing, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 32, 32), extent=(1, "1", 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 32, 32), extent=(1, 1, 1), halo=(1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(32, 32, 32), extent=(1, 1, 1), halo=(1.0, 1, 1))

    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=2)
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), y=[1, 2])
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), z=(-π, π))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=1, y=2, z=3)
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=(0, 1), y=(0, 2), z=4)
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=("0", "1"))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=(1, 2, 3))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=(1, 0), y=(1//7, 5//7), z=(1, 2))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), extent=(1, 2, 3), x=(0, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), extent=(1, 2, 3), x=(0, 1), y=(1, 5), z=(-π, π))

    @test_throws ArgumentError RegularRectilinearGrid(FT, size=(16, 16, 16), extent=(1, 1, 1), topology=(Periodic, Periodic, Flux))

    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Flat, Periodic, Periodic), size=(16, 16, 16), extent=1)
    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Periodic, Flat, Periodic), size=(16, 16, 16), extent=(1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Periodic, Periodic, Flat), size=(16, 16, 16), extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Periodic, Periodic, Flat), size=(16, 16),     extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Periodic, Periodic, Flat), size=16,           extent=(1, 1, 1))

    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Periodic, Flat, Flat), size=16, extent=(1, 1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Flat, Periodic, Flat), size=16, extent=(1, 1))
    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Flat, Flat, Periodic), size=(16, 16), extent=1)

    @test_throws ArgumentError RegularRectilinearGrid(FT, topology=(Flat, Flat, Flat), size=16, extent=1)

    return nothing
end

function flat_size_regular_rectilinear_grid(FT; topology, size, extent)
    grid = RegularRectilinearGrid(FT; size=size, topology=topology, extent=extent)
    return grid.Nx, grid.Ny, grid.Nz
end

function flat_halo_regular_rectilinear_grid(FT; topology, size, halo, extent)
    grid = RegularRectilinearGrid(FT; size=size, halo=halo, topology=topology, extent=extent)
    return grid.Hx, grid.Hy, grid.Hz
end

function flat_extent_regular_rectilinear_grid(FT; topology, size, extent)
    grid = RegularRectilinearGrid(FT; size=size, topology=topology, extent=extent)
    return grid.Lx, grid.Ly, grid.Lz
end

function test_flat_size_regular_rectilinear_grid(FT)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Periodic), size=(2, 3), extent=(1, 1)) === (1, 2, 3)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Bounded),  size=(2, 3), extent=(1, 1)) === (2, 1, 3)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Periodic, Bounded, Flat),  size=(2, 3), extent=(1, 1)) === (2, 3, 1)

    @test flat_size_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Periodic), size=(2, 3), extent=(1, 1)) === (1, 2, 3)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Bounded),  size=(2, 3), extent=(1, 1)) === (2, 1, 3)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Periodic, Bounded, Flat),  size=(2, 3), extent=(1, 1)) === (2, 3, 1)

    @test flat_size_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Flat), size=2, extent=1) === (2, 1, 1)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Flat), size=2, extent=1) === (1, 2, 1)
    @test flat_size_regular_rectilinear_grid(FT, topology=(Flat, Flat, Bounded),  size=2, extent=1) === (1, 1, 2)

    @test flat_size_regular_rectilinear_grid(FT, topology=(Flat, Flat, Flat), size=(), extent=()) === (1, 1, 1)

    @test flat_halo_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Periodic), size=(1, 1), extent=(1, 1), halo=nothing) === (0, 1, 1)
    @test flat_halo_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Bounded),  size=(1, 1), extent=(1, 1), halo=nothing) === (1, 0, 1)
    @test flat_halo_regular_rectilinear_grid(FT, topology=(Periodic, Bounded, Flat),  size=(1, 1), extent=(1, 1), halo=nothing) === (1, 1, 0)

    @test flat_halo_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Periodic), size=(1, 1), extent=(1, 1), halo=(2, 3)) === (0, 2, 3)
    @test flat_halo_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Bounded),  size=(1, 1), extent=(1, 1), halo=(2, 3)) === (2, 0, 3)
    @test flat_halo_regular_rectilinear_grid(FT, topology=(Periodic, Bounded, Flat),  size=(1, 1), extent=(1, 1), halo=(2, 3)) === (2, 3, 0)

    @test flat_halo_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Flat), size=1, extent=1, halo=2) === (2, 0, 0)
    @test flat_halo_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Flat), size=1, extent=1, halo=2) === (0, 2, 0)
    @test flat_halo_regular_rectilinear_grid(FT, topology=(Flat, Flat, Bounded),  size=1, extent=1, halo=2) === (0, 0, 2)

    @test flat_halo_regular_rectilinear_grid(FT, topology=(Flat, Flat, Flat), size=(), extent=(), halo=()) === (0, 0, 0)

    @test flat_extent_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Periodic), size=(2, 3), extent=(1, 1)) == (0, 1, 1)
    @test flat_extent_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Periodic), size=(2, 3), extent=(1, 1)) == (1, 0, 1)
    @test flat_extent_regular_rectilinear_grid(FT, topology=(Periodic, Periodic, Flat), size=(2, 3), extent=(1, 1)) == (1, 1, 0)

    @test flat_extent_regular_rectilinear_grid(FT, topology=(Periodic, Flat, Flat), size=2, extent=1) == (1, 0, 0)
    @test flat_extent_regular_rectilinear_grid(FT, topology=(Flat, Periodic, Flat), size=2, extent=1) == (0, 1, 0)
    @test flat_extent_regular_rectilinear_grid(FT, topology=(Flat, Flat, Periodic), size=2, extent=1) == (0, 0, 1)

    @test flat_extent_regular_rectilinear_grid(FT, topology=(Flat, Flat, Flat), size=(), extent=()) == (0, 0, 0)

    return nothing
end

#####
##### Vertically stretched grids
#####

function test_vertically_stretched_grid_properties_are_same_type(FT, arch)
    grid = VerticallyStretchedRectilinearGrid(FT, architecture=arch, size=(1, 1, 16), x=(0,1), y=(0,1), z_faces=collect(0:16))

    @test grid.Lx isa FT
    @test grid.Ly isa FT
    @test grid.Lz isa FT
    @test grid.Δx isa FT
    @test grid.Δy isa FT

    @test eltype(grid.xᶠᵃᵃ) == FT
    @test eltype(grid.xᶜᵃᵃ) == FT
    @test eltype(grid.yᵃᶠᵃ) == FT
    @test eltype(grid.yᵃᶜᵃ) == FT
    @test eltype(grid.zᵃᵃᶠ) == FT
    @test eltype(grid.zᵃᵃᶜ) == FT

    @test eltype(grid.Δzᵃᵃᶜ) == FT
    @test eltype(grid.Δzᵃᵃᶠ) == FT

    return nothing
end

function test_architecturally_correct_stretched_grid(FT, arch, zF)
    grid = VerticallyStretchedRectilinearGrid(FT, architecture=arch, size=(1, 1, length(zF)-1), x=(0, 1), y=(0, 1), z_faces=zF)

    ArrayType = array_type(arch)
    @test grid.zᵃᵃᶠ  isa OffsetArray{FT, 1, <:ArrayType}
    @test grid.zᵃᵃᶜ  isa OffsetArray{FT, 1, <:ArrayType}
    @test grid.Δzᵃᵃᶠ isa OffsetArray{FT, 1, <:ArrayType}
    @test grid.Δzᵃᵃᶜ isa OffsetArray{FT, 1, <:ArrayType}

    return nothing
end

function test_correct_constant_grid_spacings(FT, Nz)
    grid = VerticallyStretchedRectilinearGrid(FT, size=(1, 1, Nz), x=(0, 1), y=(0, 1), z_faces=collect(0:Nz))

    @test all(grid.Δzᵃᵃᶜ .== 1)
    @test all(grid.Δzᵃᵃᶠ .== 1)

    return nothing
end

function test_correct_quadratic_grid_spacings(FT, Nz)
    grid = VerticallyStretchedRectilinearGrid(FT, size=(1, 1, Nz), x=(0, 1), y=(0, 1), z_faces=collect(0:Nz).^2)

     zF(k) = (k-1)^2
     zC(k) = (k^2 + (k-1)^2) / 2
    ΔzF(k) = k^2 - (k-1)^2
    ΔzC(k) = zC(k+1) - zC(k)

    @test all(isapprox.(  grid.zᵃᵃᶠ[1:Nz+1],  zF.(1:Nz+1) ))
    @test all(isapprox.(  grid.zᵃᵃᶜ[1:Nz],    zC.(1:Nz)   ))
    @test all(isapprox.( grid.Δzᵃᵃᶜ[1:Nz],   ΔzF.(1:Nz)   ))

    # Note that Δzᵃᵃᶠ[1] involves a halo point, which is not directly determined by
    # the user-supplied zF
    @test all(isapprox.( grid.Δzᵃᵃᶠ[2:Nz-1], ΔzC.(2:Nz-1) ))

    return nothing
end

function test_correct_tanh_grid_spacings(FT, Nz)
    S = 3  # Stretching factor
    zF(k) = tanh(S * (2 * (k - 1) / Nz - 1)) / tanh(S)

    grid = VerticallyStretchedRectilinearGrid(FT, size=(1, 1, Nz), x=(0, 1), y=(0, 1), z_faces=zF)

     zC(k) = (zF(k) + zF(k+1)) / 2
    ΔzF(k) = zF(k+1) - zF(k)
    ΔzC(k) = zC(k+1) - zC(k)

    @test all(isapprox.(  grid.zᵃᵃᶠ[1:Nz+1],  zF.(1:Nz+1) ))
    @test all(isapprox.(  grid.zᵃᵃᶜ[1:Nz],    zC.(1:Nz)   ))
    @test all(isapprox.( grid.Δzᵃᵃᶜ[1:Nz],   ΔzF.(1:Nz)   ))

    # Note that Δzᵃᵃᶠ[1] involves a halo point, which is not directly determined by
    # the user-supplied zF
    @test all(isapprox.( grid.Δzᵃᵃᶠ[2:Nz-1], ΔzC.(2:Nz-1) ))

   return nothing
end

#####
##### Regular latitude-longitude grid tests
#####

function test_basic_lat_lon_bounded_domain(FT)
    Nλ = Nφ = 18
    Hλ = Hφ = 1

    grid = RegularLatitudeLongitudeGrid(FT, size=(Nλ, Nφ, 1), longitude=(-90, 90), latitude=(-45, 45), z=(0, 1), halo=(Hλ, Hφ, 1))

    @test topology(grid) == (Bounded, Bounded, Bounded)

    @test grid.Nx == Nλ
    @test grid.Ny == Nφ
    @test grid.Nz == 1

    @test grid.Lx == 180
    @test grid.Ly == 90
    @test grid.Lz == 1

    @test grid.Δλ == 10
    @test grid.Δφ == 5
    @test grid.Δz == 1

    @test length(grid.λᶠᵃᵃ) == Nλ + 2Hλ + 1
    @test length(grid.λᶜᵃᵃ) == Nλ + 2Hλ

    @test length(grid.φᵃᶠᵃ) == Nφ + 2Hφ + 1
    @test length(grid.φᵃᶜᵃ) == Nφ + 2Hφ

    @test grid.λᶠᵃᵃ[1] == -90
    @test grid.λᶠᵃᵃ[Nλ+1] == 90

    @test grid.φᵃᶠᵃ[1] == -45
    @test grid.φᵃᶠᵃ[Nφ+1] == 45

    @test grid.λᶠᵃᵃ[0] == -90 - grid.Δλ
    @test grid.λᶠᵃᵃ[Nλ+2] == 90 + grid.Δλ

    @test grid.φᵃᶠᵃ[0] == -45 - grid.Δφ
    @test grid.φᵃᶠᵃ[Nφ+2] == 45 + grid.Δφ

    @test all(diff(grid.λᶠᵃᵃ.parent) .== grid.Δλ)
    @test all(diff(grid.λᶜᵃᵃ.parent) .== grid.Δλ)

    @test all(diff(grid.φᵃᶠᵃ.parent) .== grid.Δφ)
    @test all(diff(grid.φᵃᶜᵃ.parent) .== grid.Δφ)

    return nothing
end

function test_basic_lat_lon_periodic_domain(FT)
    Nλ = 36
    Nφ = 32
    Hλ = Hφ = 1

    grid = RegularLatitudeLongitudeGrid(FT, size=(Nλ, Nφ, 1), longitude=(-180, 180), latitude=(-80, 80), z=(0, 1), halo=(Hλ, Hφ, 1))

    @test topology(grid) == (Periodic, Bounded, Bounded)

    @test grid.Nx == Nλ
    @test grid.Ny == Nφ
    @test grid.Nz == 1

    @test grid.Lx == 360
    @test grid.Ly == 160
    @test grid.Lz == 1

    @test grid.Δλ == 10
    @test grid.Δφ == 5
    @test grid.Δz == 1

    @test length(grid.λᶠᵃᵃ) == Nλ + 2Hλ
    @test length(grid.λᶜᵃᵃ) == Nλ + 2Hλ

    @test length(grid.φᵃᶠᵃ) == Nφ + 2Hφ + 1
    @test length(grid.φᵃᶜᵃ) == Nφ + 2Hφ

    @test grid.λᶠᵃᵃ[1] == -180
    @test grid.λᶠᵃᵃ[Nλ] == 180 - grid.Δλ

    @test grid.φᵃᶠᵃ[1] == -80
    @test grid.φᵃᶠᵃ[Nφ+1] == 80

    @test grid.λᶠᵃᵃ[0] == -180 - grid.Δλ
    @test grid.λᶠᵃᵃ[Nλ+1] == 180

    @test grid.φᵃᶠᵃ[0] == -80 - grid.Δφ
    @test grid.φᵃᶠᵃ[Nφ+2] == 80 + grid.Δφ

    @test all(diff(grid.λᶠᵃᵃ.parent) .== grid.Δλ)
    @test all(diff(grid.λᶜᵃᵃ.parent) .== grid.Δλ)

    @test all(diff(grid.φᵃᶠᵃ.parent) .== grid.Δφ)
    @test all(diff(grid.φᵃᶜᵃ.parent) .== grid.Δφ)

    return nothing
end

#####
##### Conformal cubed sphere face grid
#####

function test_cubed_sphere_face_array_size(FT)
    grid = ConformalCubedSphereFaceGrid(FT, size=(10, 10, 1), z=(0, 1))

    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz

    @test grid.λᶜᶜᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.λᶠᶜᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.λᶜᶠᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.λᶠᶠᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.φᶜᶜᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.φᶠᶜᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.φᶜᶠᵃ isa OffsetArray{FT, 2, <:Array}
    @test grid.φᶠᶠᵃ isa OffsetArray{FT, 2, <:Array}

    @test size(grid.λᶜᶜᵃ) == (Nx + 2Hx,     Ny + 2Hy    )
    @test size(grid.λᶠᶜᵃ) == (Nx + 2Hx + 1, Ny + 2Hy    )
    @test size(grid.λᶜᶠᵃ) == (Nx + 2Hx,     Ny + 2Hy + 1)
    @test size(grid.λᶠᶠᵃ) == (Nx + 2Hx + 1, Ny + 2Hy + 1)

    @test size(grid.φᶜᶜᵃ) == (Nx + 2Hx,     Ny + 2Hy    )
    @test size(grid.φᶠᶜᵃ) == (Nx + 2Hx + 1, Ny + 2Hy    )
    @test size(grid.φᶜᶠᵃ) == (Nx + 2Hx,     Ny + 2Hy + 1)
    @test size(grid.φᶠᶠᵃ) == (Nx + 2Hx + 1, Ny + 2Hy + 1)

    return nothing
end

#####
##### Test the tests
#####

@testset "Grids" begin
    @info "Testing grids..."

    @testset "Grid utils" begin
        @info "  Testing grid utilities..."

        @test total_extent(Periodic, 1, 0.2, 1.0) == 1.2
        @test total_extent(Bounded, 1, 0.2, 1.0) == 1.4
    end

    @testset "Regular rectilinear grid" begin
        @info "  Testing regular rectilinear grid..."

        @testset "Grid initialization" begin
            @info "    Testing grid initialization..."

            for FT in float_types
                test_regular_rectilinear_correct_size(FT)
                test_regular_rectilinear_correct_extent(FT)
                test_regular_rectilinear_correct_coordinate_lengths(FT)
                test_regular_rectilinear_correct_halo_size(FT)
                test_regular_rectilinear_correct_halo_faces(FT)
                test_regular_rectilinear_correct_first_cells(FT)
                test_regular_rectilinear_correct_end_faces(FT)
                test_regular_rectilinear_ranges_have_correct_length(FT)
                test_regular_rectilinear_no_roundoff_error_in_ranges(FT)
                test_regular_rectilinear_grid_properties_are_same_type(FT)
                test_xnode_ynode_znode_are_correct(FT)
            end
        end

        @testset "Grid dimensions" begin
            @info "    Testing grid constructor errors..."

            for FT in float_types
                test_regular_rectilinear_constructor_errors(FT)
            end
        end

        @testset "Grids with flat dimensions" begin
            @info "    Testing construction of grids with Flat dimensions..."

            for FT in float_types
                test_flat_size_regular_rectilinear_grid(FT)
            end
        end

        # Testing show function
        topo = (Periodic, Periodic, Periodic)

        grid = RegularRectilinearGrid(topology=topo, size=(3, 7, 9), x=(0, 1), y=(-π, π), z=(0, 2π))

        @test try
            CUDA.allowscalar(false)
            show(grid); println()
            CUDA.allowscalar(true)
            true
        catch err
            println("error in show(::RegularRectilinearGrid)")
            println(sprint(showerror, err))
            false
        end

        @test grid isa RegularRectilinearGrid
    end

    @testset "Vertically stretched rectilinear grid" begin
        @info "  Testing vertically stretched rectilinear grid..."

        for arch in archs, FT in float_types
            @testset "Vertically stretched rectilinear grid construction [$(typeof(arch)), $FT]" begin
                @info "    Testing vertically stretched rectilinear grid construction [$(typeof(arch)), $FT]..."

                test_vertically_stretched_grid_properties_are_same_type(FT, arch)

                zF1 = collect(0:10).^2
                zF2 = [1, 3, 5, 10, 15, 33, 50]
                for zF in [zF1, zF2]
                    test_architecturally_correct_stretched_grid(FT, arch, zF)
                end
            end

            @testset "Vertically stretched rectilinear grid spacings [$(typeof(arch)), $FT]" begin
                @info "    Testing vertically stretched rectilinear grid spacings [$(typeof(arch)), $FT]..."
                for Nz in [16, 17]
                    test_correct_constant_grid_spacings(FT, Nz)
                    test_correct_quadratic_grid_spacings(FT, Nz)
                    test_correct_tanh_grid_spacings(FT, Nz)
                end
            end

            # Testing show function
            Nz = 20
            grid = VerticallyStretchedRectilinearGrid(architecture=arch, size=(1, 1, Nz-1), x=(0, 1), y=(0, 1), z_faces=collect(0:Nz).^2)

            @test try
            CUDA.allowscalar(false)
            show(grid); println()
            CUDA.allowscalar(true)
                true
            catch err
                println("error in show(::VerticallyStretchedRectilinearGrid)")
                println(sprint(showerror, err))
                false
            end

            @test grid isa VerticallyStretchedRectilinearGrid
        end
    end

    @testset "Regular latitude-longitude grid" begin
        @info "  Testing regular latitude-longitude grid..."

        for FT in float_types
            test_basic_lat_lon_bounded_domain(FT)
            test_basic_lat_lon_periodic_domain(FT)
        end

        # Testing show function
        grid = RegularLatitudeLongitudeGrid(size=(36, 32, 1), longitude=(-180, 180), latitude=(-80, 80), z=(0, 1))

        @test try
            CUDA.allowscalar(false)
            show(grid); println()
            CUDA.allowscalar(true)
            true
        catch err
            println("error in show(::RegularLatitudeLongitudeGrid)")
            println(sprint(showerror, err))
            false
        end

        @test grid isa RegularLatitudeLongitudeGrid
    end

    @testset "Conformal cubed sphere face grid" begin
        @info "  Testing conformal cubed sphere face grid..."

        for FT in float_types
            test_cubed_sphere_face_array_size(Float64)
        end

        # Testing show function
        grid = ConformalCubedSphereFaceGrid(size=(10, 10, 1), z=(0, 1))

        @test try
            CUDA.allowscalar(false)
            show(grid); println()
            CUDA.allowscalar(true)
            true
        catch err
            println("error in show(::ConformalCubedSphereFaceGrid)")
            println(sprint(showerror, err))
            false
        end

        @test grid isa ConformalCubedSphereFaceGrid
    end

    @testset "Conformal cubed sphere face grid from file" begin
        @info "  Testing conformal cubed sphere face grid construction from file..."

        cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"

        for face in 1:6
            grid = ConformalCubedSphereFaceGrid(cs32_filepath, face=face, Nz=1, z=(-1, 0))
            @test grid isa ConformalCubedSphereFaceGrid
        end
    end
end