text:
# using Revise, Plots, Test
using BifurcationKit, LinearAlgebra, Setfield, SparseArrays
const BK = BifurcationKit

k = 2
N = 1
F = (x, p) -> p[1] .* x .+ x.^(k+1)/(k+1) .+ 0.01
Jac_m = (x, p) -> diagm(0 => p[1] .+ x.^k)

####################################################################################################
# test creation of specific scalar product
_dt = BK.DotTheta()
# tests for the predictors
BK.emptypredictor!(nothing)
BK.mergefromuser(1., (a=1,))
BK.mergefromuser(rand(2), (a=1,))
BK.mergefromuser((1,2), (a=1,))
####################################################################################################

normInf = x -> norm(x, Inf)

opts = ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds=0.001, maxSteps = 140, pMin = -3., saveSolEveryStep = 0, newtonOptions = NewtonPar(tol = 1e-8, verbose = false), saveEigenvectors = false, detectBifurcation = 0)
x0 = 0.01 * ones(N)

opts = @set opts.doArcLengthScaling = true
br0, = @time continuation(F, Jac_m, x0, -1.5, (@lens _), opts) #(16.12 k allocations: 772.250 KiB)
BK.getfirstusertype(br0)
BK.propertynames(br0)
BK.computeEigenvalues(opts)
BK.computeEigenvectors(opts)

# test with callbacks
br0, = continuation(F,Jac_m,x0, -1.5, (@lens _), (@set opts.maxSteps = 3), callbackN = (x, f, J, res, iteration, itlinear, optionsN; kwargs...)->(true))

###### Used to check type stability of the methods
# using RecursiveArrayTools
iter = ContIterable(F, Jac_m, x0, -1.5, (@lens _), opts)
state = iterate(iter)[1]
# test copy, copyto!
state1 = copy(state);copyto!(state1, state)
contRes = ContResult(iter, state)
continuation!(iter, state, contRes)
eltype(iter)
length(iter)
#
typeof(contRes)

# state = iterate(iter)[1]
# 	 contRes = BK.ContResult(iter, state)
# 	 @code_warntype continuation!(iter, state, contRes)
#####

opts = ContinuationPar(opts;detectBifurcation = 1,saveEigenvectors=true)
br1, sol, _ = continuation(F,Jac_m,x0,-1.5, (@lens _),opts) #(14.28 k allocations: 1001.500 KiB)
show(br1)
length(br1)
br1[1]
BK.eigenvals(br1,20)
BK.eigenvec(br1,20,1)
BK.haseigenvector(br1)
BK.getvectortype(br1)
BK.getvectoreltype(br1)
br1.param
br1.params

br2, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, printSolution = (x,p) -> norm(x,2))

# test for different norms
br3, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, normC = normInf)

# test for linesearch in Newton method
opts = @set opts.newtonOptions.linesearch = true
br4, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, normC = normInf) # (15.61 k allocations: 1.020 MiB)

# test for different ways to solve the bordered linear system arising during the continuation step
opts = @set opts.newtonOptions.linesearch = false
br5, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, normC = normInf, linearAlgo = BorderingBLS())

br5, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, normC = normInf, linearAlgo = MatrixBLS())

# test for stopping continuation based on user defined function
finaliseSolution = (z, tau, step, contResult; k...) -> (step < 20)
br5a, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, finaliseSolution = finaliseSolution)
@test length(br5a.branch) == 21

# test for different predictors
br6, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, tangentAlgo = SecantPred())

optsnat = setproperties(opts; ds = 0.001, dsmax = 0.1, dsmin = 0.0001)
br7, = continuation(F,Jac_m,x0,-1.5, (@lens _),optsnat, tangentAlgo = NaturalPred(),printSolution = (x,p)->x[1])

# tangent prediction with Bordered predictor
br8, sol, _ = continuation(F,Jac_m,x0,-1.5, (@lens _),opts, tangentAlgo = BorderedPred(),printSolution = (x,p)->x[1])

# tangent prediction with Multiple predictor
opts9 = (@set opts.newtonOptions.verbose=false)
	opts9 = ContinuationPar(opts9; maxSteps = 48, ds = 0.015, dsmin = 1e-5, dsmax = 0.05)
	br9, = continuation(F,Jac_m,x0,-1.5, (@lens _),opts9,
	printSolution = (x,p)->x[1],
	tangentAlgo = MultiplePred(copy(x0), 0.01,13)
	)
	BK.emptypredictor!(BK.MultiplePred(copy(x0), 0.01,13))
	# plot(br9, title = "$(length(br9))",marker=:d,vars=(:p,:sol),plotfold=false)

# tangent prediction with Polynomial predictor
polpred = PolynomialPred(BorderedPred(),2,6,x0)
	opts9 = (@set opts.newtonOptions.verbose=false)
	opts9 = ContinuationPar(opts9; maxSteps = 76, ds = 0.005, dsmin = 1e-4, dsmax = 0.02, plotEveryStep = 3,)
	br10, = continuation(F, Jac_m, x0, -1.5, (@lens _), opts9,
	tangentAlgo = polpred, plot=false,
	printSolution = (x,p)->x[1],
	)
	# plot(br10) |> display
	polpred(0.1)
	BK.emptypredictor!(polpred)
	# plot(br10, title = "$(length(br10))",marker=:dplot,fold=false)
	# plot!(br9)

# polpred(0.0)
# 	for _ds in LinRange(-0.51,.1,161)
# 		_x,_p = polpred(_ds)
# 		# @show _ds,_x[1], _p
# 		scatter!([_p],[_x[1]],label="",color=:blue, markersize = 1)
# 	end
# 	scatter!(polpred.parameters,map(x->x[1],polpred.solutions), color=:green)
# 	scatter!([polpred(0.01)[2]],[polpred(0.01)[1][1]],marker=:cross)
# 	# title!("",ylims=(-0.1,0.22))
# # BK.isready(polpred)
# # polpred.coeffsPar

# test for polynomial predictor interpolation
# polpred = BK.PolynomialPred(4,9,x0)
# 	for (ii,v) in enumerate(LinRange(-5,1.,10))
# 		if length(polpred.arclengths)==0
# 			push!(polpred.arclengths, 0.1)
# 		else
# 			push!(polpred.arclengths, polpred.arclengths[end]+0.1)
# 		end
# 		push!(polpred.solutions, [v])
# 		push!(polpred.parameters, 1-v^2+0.001v^4)
# 	end
# 	BK.updatePred!(polpred)
# 	polpred(-0.5)
#
# plot()
# 	scatter(polpred.parameters, reduce(vcat,polpred.solutions))
# 	for _ds in LinRange(-1.5,.75,121)
# 		_x,_p = polpred(_ds)
# 		scatter!([_p],[_x[1]],label="", color=:blue, markersize = 1)
# 	end
# 	scatter!([polpred(0.1)[2]],[polpred(0.1)[1][1]],marker=:cross)
# 	title!("",)

# further testing with sparse Jacobian operator
Jac_sp_simple = (x, p) -> SparseArrays.spdiagm(0 => p  .+ x.^k)
brsp, = continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts)
brsp, = continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts, printSolution = (x,p) -> norm(x,2))
brsp, = continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts,linearAlgo = BK.BorderingBLS())
brsp, = continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts,linearAlgo = BK.MatrixBLS())
# plotBranch(br1,marker=:d);title!("")
# plotBranch!(br8,marker=:d);title!("")
####################################################################################################
# check bounds for all predictors / correctors
for talgo in [BorderedPred(), SecantPred()]
	brbd,  = continuation(F,Jac_m,x0,-3.0, (@lens _), opts; tangentAlgo = talgo, verbosity = 0)
	@test length(brbd) > 2
end
brbd,  = continuation(F,Jac_m,ones(N)*3,-3.0, (@lens _), ContinuationPar(opts, pMax = -2))
brbd,  = continuation(F,Jac_m,ones(N)*3,-3.2, (@lens _), ContinuationPar(opts, pMax = -2), verbosity = 3)
@test isnothing(brbd)
####################################################################################################
# testing when starting with 2 points on the branch
opts = BK.ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds=0.001, maxSteps = 140, pMin = -3., saveSolEveryStep = 0, newtonOptions = NewtonPar(verbose = false), detectBifurcation = 3)
x0 = 0.01 * ones(2)

x0, = newton(F,Jac_m,x0, -1.5, opts.newtonOptions)
x1, = newton(F,Jac_m,x0, -1.45, opts.newtonOptions)

br0, = continuation(F,Jac_m, x0, -1.5, (@lens _), opts, verbosity=0)
BK.getEigenelements(br0, br0.specialpoint[1])

br1, = continuation(F,Jac_m, x1, -1.45, x0, -1.5, (@lens _), ContinuationPar(opts; ds = -0.001))

br2, = continuation(F,Jac_m,x0, -1.5, x1, -1.45, (@lens _), opts; tangentAlgo = BorderedPred())
####################################################################################################
# test for computing both sides
br3, = continuation(F,Jac_m,x0, -1.5, (@lens _), opts; tangentAlgo = BorderedPred(), bothside = true)
####################################################################################################
# test for deflated continuation
brdc, = continuation(F,Jac_m, 0.5, (@lens _),
	ContinuationPar(opts, ds = -0.001, maxSteps = 800, newtonOptions = NewtonPar(verbose = false, maxIter = 6), plotEveryStep = 40),
	DeflationOperator(2, .001, [[0.]]); showplot=false, verbosity = 0,
	perturbSolution = (x,p,id) -> (x  .+ 0.1 .* rand(length(x))),
	callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> res <1e3)


text:
export @codegen
macro codegen()
  str="""
	function print(ex...)
		println(ex...)
	end"""
  eval(parse(str))
end

ext: jl


text:
module TestPaperSelfgravitatingGasDynamics

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "paper-self-gravitating-gas-dynamics")

# Numerical examples from the Euler-gravity paper
@testset "paper-self-gravitating-gas-dynamics" begin
  @testset "parameters_eoc_test_euler.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_euler.toml"),
            l2   = [0.00017409779099463607, 0.0003369287450282371, 0.00033692874502819616, 0.0006099035183426747],
            linf = [0.0010793454782482836, 0.0018836374478419238, 0.0018836374478410356, 0.003971446179607874])
  end

  @testset "parameters_eoc_test_euler.toml with polydeg=4" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_euler.toml"),
            l2   = [1.7187032983384504e-5, 2.6780178144541376e-5, 2.678017814469407e-5, 4.952410417693103e-5],
            linf = [0.00015018092862240096, 0.00016548331714294484, 0.00016548331714405506, 0.00043726245511699346],
            polydeg=4)
  end

  @testset "parameters_eoc_test_hyperbolic_diffusion.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_hyperbolic_diffusion.toml"),
            l2   = [0.00315402168051244, 0.012394424055283394, 0.021859728673870843],
            linf = [0.017332075119072865, 0.07843510773347322, 0.11325788389718668])
  end

  @testset "parameters_eoc_test_hyperbolic_diffusion.toml with polydeg=4" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_hyperbolic_diffusion.toml"),
            l2   = [0.00025112830138292663, 0.0008808243851096586, 0.0016313343234903468],
            linf = [0.001719090967553516, 0.0031291844657076145, 0.00994609342322228],
            polydeg=4)
  end

  @testset "parameters_eoc_test_coupled_euler_gravity.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_coupled_euler_gravity.toml"),
            l2   = [0.0002487158370511598, 0.0003370291440916084, 0.00033702914409161063, 0.0007231934514459757, 0.00013852406160669235, 0.0007541252869723029, 0.0007541252869723299],
            linf = [0.001581173125044355, 0.002049389755695241, 0.0020493897556961294, 0.004793721268126383, 0.0009549587622960237, 0.0030981236291827237, 0.003098123629182964],
            t_end=0.1)
  end

  @testset "parameters_eoc_test_coupled_euler_gravity.toml with polydeg=4" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_coupled_euler_gravity.toml"),
            l2   = [1.9536732064098693e-5, 2.756381055173374e-5, 2.7563810551703437e-5, 5.688705902953846e-5, 1.0684325470325204e-5, 5.829033623593028e-5, 5.829033623591347e-5],
            linf = [0.00012335977351507488, 0.00020086338378089152, 0.00020086338378044744, 0.0004962132679873221, 8.5358666522109e-5, 0.0002927883863423353, 0.00029278838634330673],
            t_end=0.1, polydeg=4)
  end

  @testset "parameters_eoc_test_coupled_euler_gravity.toml with update_gravity_once_per_stage=false" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_coupled_euler_gravity.toml"),
            l2   = [0.00039754984075255105, 0.0011317710289437735, 0.00113177102894379, 0.002302567979915388, 0.0002449228820755184, 0.0009838245219995854, 0.0009838245219995676],
            linf = [0.0013753628419428399, 0.0031706120730756737, 0.0031706120730756737, 0.0069469754604232214, 0.0008251489152860739, 0.0030597255494218545, 0.0030597255494215977],
            t_end=0.1, update_gravity_once_per_stage=false)
  end

  @testset "parameters_eoc_test_coupled_euler_gravity.toml with 1st order RK3S*" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_coupled_euler_gravity.toml"),
            l2   = [0.00024871583705119436, 0.0003370291440915927, 0.0003370291440916112, 0.0007231934514459859, 0.00013852406160669225, 0.0007541252869723031, 0.0007541252869723208],
            linf = [0.001581173125044355, 0.002049389755695241, 0.0020493897556961294, 0.004793721268126383, 0.0009549587622960237, 0.0030981236291827237, 0.003098123629182964],
            t_end=0.1, time_integration_scheme_gravity="timestep_gravity_erk51_3Sstar!")
  end

  @testset "parameters_eoc_test_coupled_euler_gravity.toml with 3rd order RK3S*" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_coupled_euler_gravity.toml"),
            l2   = [0.000248715837047647, 0.0003370291440257414, 0.00033702914402587556, 0.0007231934513057375, 0.0001385240616293125, 0.0007541252869544295, 0.0007541252869544261],
            linf = [0.00158117312532835, 0.0020493897540796446, 0.0020493897540800887, 0.0047937212650124295, 0.000954958762033685, 0.0030981236303003834, 0.003098123630300921],
            t_end=0.1, time_integration_scheme_gravity="timestep_gravity_erk53_3Sstar!")
  end

  @testset "parameters_jeans_instability.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_jeans_instability.toml"),
            l2   = [10734.053145619571, 13357.12732342943, 1.689152281530251e-6, 26835.125221793754],
            linf = [15194.889945842326, 18881.97158566768, 8.787452968447721e-6, 37974.480818059295],
            t_end=0.1)
  end

  @testset "parameters_jeans_instability.toml with update_gravity_once_per_stage=false" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_jeans_instability.toml"),
            l2   = [10723.612237729343, 13336.21827245215, 2.602694645664843e-6, 26809.022944817934],
            linf = [15180.113784657791, 18852.413654859887, 9.130954929048019e-6, 37937.54015407339],
            t_end=0.1, update_gravity_once_per_stage=false)
  end

  @testset "parameters_jeans_instability.toml with RK3S*" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_jeans_instability.toml"),
            l2   = [10734.450531008228, 13358.15805096243, 2.0427912643500984e-6, 26836.118685698464],
            linf = [15195.451979616657, 18883.428240704106, 1.0980410384074066e-5, 37975.88591177389],
            t_end=0.1, time_integration_scheme_gravity="timestep_gravity_erk52_3Sstar!",
            cfl_gravity=1.2)
  end

  @testset "parameters_sedov_self_gravity.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_sedov_self_gravity.toml"),
            l2   = [0.04630745182870653, 0.06507397069667138, 0.06507397069667123, 0.48971269294890085],
            linf = [2.3861430058270847, 4.083635578775231, 4.083635578775232, 16.246070713311475],
            t_end=0.05)
  end

  @testset "parameters_sedov_self_gravity.toml with amr_interval=0 and ref-level=8" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_sedov_self_gravity.toml"),
            l2   = [0.0028922121586238323, 0.013724796675028317, 0.013724796675028307, 0.05822941648860658],
            linf = [0.26747911779347966, 1.385822018653034, 1.385822018653034, 4.071204772447614],
            t_end=0.005, amr_interval=0, initial_refinement_level=8)
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Word2Vec_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Word2Vec")
JLLWrappers.@generate_main_file("Word2Vec", UUID("9fbe4022-c126-5389-b4b2-756cc9f654d0"))
end  # module Word2Vec_jll


text:
## Preformance Evaluation
using BSON: @load

"""
confusion_matrix(k::Integer, gts, preds):
Generates confusion matrix between gts and pred.
"""
function confusion_matrix(k::Integer, gts, preds)
    n = length(gts)
    length(preds) == n || throw(DimensionMismatch("Inconsistent lengths."))
    R = zeros(Int, k, k)
    for i = 1:n
        @inbounds g = gts[i]
        @inbounds p = preds[i]
        R[g, p] += 1
    end
    return R
end

"""
function error_stats(gts, pred, err):
Calculates error for each data point and overall mean error and std. Type of error is specified by function passed to err.
"""
function error_stats(gts, pred, err)
    error = err.(gts, pred)
    μ = round(mean(error), digits=6)
    σ = round(std(error), digits=6)
    return error, μ, σ
end

"""
rel_error(gts,pred):
Returns the relative error between gts and pred
"""
rel_error(gts,pred) = abs((gts-pred)/gts)

"""
calculate_leadingpiece(X):
Calculates leading order term of X. Depends on global parameter ONTOLOGY.
"""
function calculate_leadingpiece(X)
    leading = []
    for i ∈ 1:size(X,2)
        v = X[:,i]
        for x ∈ v[collect(length(v):-1:1)]
            (x == 0.0) && pop!(v)
            (x != 0.0) && break
        end
        while length(v) > 4; pop!(v);end
        pts = v[2]
        dualpts = v[4]
        if ONTOLOGY == "h11"
            push!(leading, pts - (DIM + 1))
        elseif ONTOLOGY == "h21"
            push!(leading, dualpts - (DIM + 1))
        else ONTOLOGY == "euler"
            push!(leading, 2*(pts - dualpts))
        end
    end
    return leading
end

"""
relevance(W, b, X):
Outputs vector R containing neuron-by-neuron relevance for each layer.
"""
function relevance(W, b, X)
    a = get_activations(W, b, X, relu)
    R = [zeros(length(a[i])) for i ∈ 1:length(a)]
    R₀ = [0.0 for i ∈ 1:length(a)]
    R[end][argmax(softmax(a[end]))] = 1
    # ϵ = [0.0, 0.25, 0.0]
    # γ = [0.25, 0.0, 0.0]
    # lrp_layers = [floor((length(W) - 1) / 3) for i=1:3]
    # remainder = length(W) - sum(lrp_layers)
    # remainder > 0 && (lrp_layers[3] = lrp_layers[3] + remainder)
    # lrp_layers = cumsum(lrp_layers)
    # for m ∈ length(W):-1:1
    #     ms = m*ones(3)
    #     α = findall(x -> 0 ≤ x ≤ m + 1, lrp_layers - ms)
    #     for j ∈ 1:length(R[m])
    #         for k ∈ 1:length(R[m + 1])
    #             z = sum([a[m][l]*(W[m][k,l] + (γ[α] * relu(W[m][k,l]) )[1] ) for l ∈ 1:length(a[m])]) + b[m][k] + ϵ[α][1]
    #             s = R[m + 1][k] / z
    #             c1 = W[m][k,j] * s
    #             j == 1 && (R₀[m] += (b[m][k] * s))
    #             R[m][j] += a[m][j] * c1
    #         end
    #     end
    # end
    for m ∈ length(W):-1:1
        for j ∈ 1:length(R[m])
            for k ∈ 1:length(R[m + 1])
                z = sum([a[m][l]*W[m][k,l] for l ∈ 1:length(a[m])]) + b[m][k]
                s = R[m + 1][k] / z
                c1 = W[m][k,j] * s
                j == 1 && (R₀[m] += (b[m][k] * s))
                R[m][j] += a[m][j] * c1
            end
        end
    end
    return (R, R₀)
end

"""
get_activations(W, b, X, f):
Feeds X forward through the model defined by (W, b) with activation function f.
Returns list a containing the values of each neuron for every layer.
"""
function get_activations(W, b, X, f)
    a = [X]
    for m ∈ 1:length(W)
        push!(a, f.(W[m]*a[m] + b[m]))
    end
    return a
end

"""
Layer-wise Relevance Propogation.
LRP(model_dir, data) : (model, X) ⟶ R = [r₁, … ,rₙ] where rᵢ is the relevance of layer i
"""
function LRP(model_dir, data; normalize=false)
    @load model_dir model
    P = params(model)
    W = [P[i] for i ∈ 1:2:length(P)]
    b = [P[i] for i ∈ 2:2:length(P)]
    R = zeros(length(data[1]), length(data))
    for (i,d) ∈ enumerate(data)
        r = relevance(W, b, d)
        R[:,i] = r[1][1]
    end
    if normalize == true
        R_max = map(maximum, eachcol(abs.(R)))
        for i = 1:10
            R[:,i] = abs.(R[:,i] / R_max[i])
        end
    end
    return R
end

function LRP(model::Chain, data; normalize=false)
    P = params(model)
    W = [P[i] for i ∈ 1:2:length(P)]
    b = [P[i] for i ∈ 2:2:length(P)]
    R = zeros(length(data[1]), length(data))
    R₀ = []
    for (i,d) ∈ enumerate(data)
        r = relevance(W, b, d)
        R[:,i] = r[1][1]
        push!(R₀ ,r[2])
    end
    if normalize == true
        R_max = map(maximum, eachcol(abs.(R)))
        for i = 1:10
            R[:,i] = abs.(R[:,i] / R_max[i])
        end
    end
    return R, R₀
end


text:
"""
$(TYPEDEF)

Paired sample test for proportion

Constructors
------------
* `McNemarProp(p1::Real, p2::Real)`

Arguments
---------
* `p1`: Proportion for cell of no event - event

* `p2`: Proportion for cell of event - no event

Fields
------
$(FIELDS)
"""
type McNemarProp <: TrialTest
    p1::Real
    p2::Real

    # Validator
    function McNemarProp(p1, p2)

        if !((0 < p1 < 1) & (0 < p2 < 1))

            error("Proportion values must be in (0, 1)")

        end # end if
        new(p1, p2)

    end # end function

end # function

# Two sample test for proportion
function hypotheses{T <: McNemarProp}(test::T, n::Real, std::Union{Real, Tuple{Real, Real}, Void}, alpha::Real, side::String)

    pdisc = test.p1 + test.p2
    pdiff = test.p2 - test.p1
    if side == "two"

        p = cdf(ZDIST, (pdiff * sqrt(n) - quantile(ZDIST, 1 - alpha / 2) * sqrt(pdisc)) / sqrt(pdisc - pdiff ^ 2)) +
            cdf(ZDIST, (-pdiff * sqrt(n) - quantile(ZDIST, 1 - alpha / 2) * sqrt(pdisc)) / sqrt(pdisc - pdiff ^ 2))

    else

        p = cdf(ZDIST, (abs(pdiff) * sqrt(n) - quantile(ZDIST, 1 - alpha) * sqrt(pdisc)) / sqrt(pdisc - pdiff ^ 2))

    end # end if
    return p

end # end function


text:
module CorrelationWrapper

using JavaCall

const Class = @jimport java.lang.Class
const Covariance = @jimport org.hipparchus.stat.correlation.Covariance
const JString = @jimport java.lang.String
const KendallsCorrelation = @jimport org.hipparchus.stat.correlation.KendallsCorrelation
const Object = @jimport java.lang.Object
const PearsonsCorrelation = @jimport org.hipparchus.stat.correlation.PearsonsCorrelation
const RankingAlgorithm = @jimport org.hipparchus.stat.ranking.RankingAlgorithm
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const SpearmansCorrelation = @jimport org.hipparchus.stat.correlation.SpearmansCorrelation
const StorelessBivariateCovariance = @jimport org.hipparchus.stat.correlation.StorelessBivariateCovariance
const StorelessCovariance = @jimport org.hipparchus.stat.correlation.StorelessCovariance

export Covariance
export KendallsCorrelation
export PearsonsCorrelation
export SpearmansCorrelation
export StorelessBivariateCovariance
export StorelessCovariance
export append
export compute_correlation_matrix
export correlation
export covariance
export covariance_to_correlation
export equals
export get_class
export get_correlation_matrix
export get_correlation_p_values
export get_correlation_standard_errors
export get_covariance
export get_covariance_matrix
export get_data
export get_n
export get_rank_correlation
export get_result
export hash_code
export increment
export notify
export notify_all
export to_string
export wait

include("covariance.jl")
include("kendalls_correlation.jl")
include("pearsons_correlation.jl")
include("spearmans_correlation.jl")
include("storeless_bivariate_covariance.jl")
include("storeless_covariance.jl")

end


text:
"""
test_args(args...)

# Description

Testing to see if I can create documentation from a different folder

# Arguments
- `args...`: (arbitrary as long as it works with println).

# Return
- Nothing

# Additional details
"""
function test_args(args...)
    for i in args
        println(i)
    end
    return nothing
end

"""
test_args_kw(args..., keyword = "keyword")

# Description

Testing to see if I can create documentation from a different folder

# Arguments
- `args...`: (arbitrary as long as it works with println).

# Keyword Arguments
- `keyword`: default = "keyword"

# Return
- Nothing

# Additional details
Just another test
"""
function test_args_kw(args...; keyword = "keyword")
    println("The function arguments are")
    for i in args
        println(i)
    end
    println("The keyword argument is ")
    println(keyword)
    return nothing
end

max_stars_repo_name: sandreza/TrialDocs

max_issues_repo_name: sandreza/TrialDocs

max_forks_repo_name: sandreza/TrialDocs


text:
module TwoBodyUnitTests

using Test
using UnitfulAstrodynamics
using UnitfulAstrodynamics.TwoBody.Systems: Earth

@testset "Transforms" begin

    rᵢ = [0.0, 11681.0, 0.0] * u"km"
    vᵢ = [5.134, 4.226, 2.787] * u"km/s"
    orbit = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s") |> KeplerianOrbit

    @test semimajor_axis(orbit) == 24509.265399338536 * u"km"
    @test eccentricity(orbit) == 0.723452708202361
    @test inclination(orbit) == 151.50460766373865 * u"°"
    @test true_anomoly(orbit) == 89.99652573907436  * u"°"

    @test orbit ≈ CartesianOrbit(orbit)

    e      =  0.3
    a      =  semimajor_axis(orbit)
    i      =  10.      * u"°"
    Ω      =  0.       * u"°"
    ω      =  10.      * u"°"
    ν      =  0.       * u"°"
    orbit  =  KeplerianOrbit(e, a, i, Ω, ω, ν, Earth, 0.0u"s")

    @test isapprox(orbit, CartesianOrbit(orbit); atol=1e-6)

end

@testset "Kepler" begin

    rᵢ = [0.0, 11681.0, 0.0]u"km"
    vᵢ = [5.134, 4.226, 2.787]u"km/s"
    orbit = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s")

    @test kepler(orbit) ≈ orbit

end

@testset "Lambert" begin

    rᵢ = [0.0, 11681.0, 0.0]u"km"
    vᵢ = [5.134, 4.226, 2.787]u"km/s"
    initial = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s")

    Δt = 1000u"s"
    final = kepler(initial, Δt; tol=1e-12)

    v₁, v₂ = lambert(position_vector(initial), position_vector(final), mass_parameter(Earth), Δt, :short; tol=1e-6, max_iter=1000)

    @test isapprox(ustrip.(u"km/s", v₁), ustrip.(u"km/s", velocity_vector(initial)); atol=1e-6)
    @test isapprox(ustrip.(u"km/s", v₂), ustrip.(u"km/s", velocity_vector(final)); atol=1e-6)

end

end


text:
mutable struct UpdateState{LT, GT, PST, OPT}
    loss::LT
    grads::GT
    params::PST
    opt::OPT
end

UpdateState() = UpdateState(nothing, nothing, nothing, nothing)


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 912825ec-1e5f-11ec-13b6-f7222876a7d5
using Dates

# ╔═╡ b7a923d4-f637-4f9e-8576-5371e9d72c88
isleapyear(1970)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═912825ec-1e5f-11ec-13b6-f7222876a7d5
# ╠═b7a923d4-f637-4f9e-8576-5371e9d72c88
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
include("../../src/ToyModelsMD.jl")
using .ToyModelsMD

function cbfunc(xstore::Array,rstore::Array,system::AbstractSystem,args...)
  # calculate the stuff
  push!(xstore,copy(system.x))
  push!(rstore,copy(system.thermostat.rands))
end

# define system
Mx = 21
My = 21
Mt = 21
shape = (Mx*My*Mt,2)
mux = rand(Float64,21)
muy = rand(Float64,21)
mut = rand(Float64,21)
sigx = rand(Float64,21)
sigy = rand(Float64,21)
sigt = rand(Float64,21)
cs = zeros(Float64,shape)
gm = GaussianModel2D(Mx,My,Mt,mux,muy,mut,sigx,sigy,sigt,cs)
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(gm,langevin)

# define callback function
xstore = []
rstore = []
cbf(system,args...) = cbfunc(xstore,rstore,system,args)
cb = Callback(1,cbf)

dt = 0.0001
tf = 0.0015
nsteps = convert(Int64, round(tf/dt))

system.x = [0.63, 0.03]
integrator = StochasticEuler(system,gm)
runtraj!(nsteps,dt,integrator,cb)
println(xstore)
println(dt*rstore)


text:
using Exam
using Test
using CSV
using Random
using DataFrames

function generate_data()
    open("exam_results.csv", "w") do io
        write(io, "name,id,bonus,1,2,3\n")
        for i = 1:20
           write(io, randstring('a':'z', 6), ",", string(rand(1:Int(1e12))), ",",string(rand((false, true))), ",",
                     string(rand(0:5)), ",", string(rand(0:2)), ",", string(rand(0:10)), "\n")
        end
    end
end

@testset  begin
    correct_results = CSV.read("simple_test_output.csv", DataFrame)
    results = evaluate("simple_results.csv", "simple_exam.yml", histogram=false)
    @test results == correct_results

    correct_results = CSV.read("simple_test_output_anonymous.csv", DataFrame)
    results = evaluate("simple_results.csv", "simple_exam.yml", anonymize=true, histogram=false)
    @test results == correct_results

end


text:
module IncrementalSVD

using ProgressMeter, InfoZIP, Random, LinearAlgebra

export Rating, RatingSet, RatingsModel
export train
export split_ratings, rmse, truncate_model!, cosine_similarity
export items, users, item_features, user_features, show_items_by_features, item_search
export similar_items, similar_users, user_ratings, get_predicted_rating
export load_book_crossing_dataset, load_small_movielens_dataset, load_large_movielens_dataset

include("types.jl")
include("train.jl")
include("util.jl")

include("book_crossing_dataset.jl")
include("movielens_dataset.jl")

end


text:
module NativeSVG

using NodeJS

export Drawing
export str, cdata, latex

const PRIMITIVES = Dict(
    :svg => true,
    :g => true,
    :defs => true,
    :use => false,
    :desc => true,
    :title => true,
    :style => true,
    :path => false,
    :line => false,
    :rect => false,
    :polygon => false,
    :polyline => false,
    :circle => false,
    :ellipse => false,
    :text => true,
    :tspan => true,
    :textPath => true,
    :linearGradient => true,
    :radialGradient => true,
    :stop => false,
    :pattern => true,
    :marker => true,
    :clipPath => true,
    :style => true,
    :foreignObject => true,
)

for primitive in keys(PRIMITIVES)
    eval(quote
        export $primitive
    end)
end

export Turtle

const COMMANDS = (:forward, :turn, :penup, :pendown, :pencolor,)

for command in COMMANDS
    eval(quote
        export $command
    end)
end

include("svg.jl")
include("turtles.jl")

end

max_stars_repo_name: BenLauwens/NativeSVG

max_issues_repo_name: BenLauwens/NativeSVG

max_forks_repo_name: BenLauwens/NativeSVG


text:
#=
    The File contains tests for find_basis Algorithms
=#

#------------------------------------------------------------------------------

include("../src/wiedemannchik.jl")

#------------------------------------------------------------------------------

if !isdefined(Main, :testset)
    using Test
    using TestSetExtensions
end

#------------------------------------------------------------------------------

function test_minimal_polynomial(minpoly)
    #=
        Only annihilating property is tested
        Minimality property is ommited
    =#

    S, x = QQ["x"]

    # hardcoded tests
    A = from_COO(1, 1, [(1, 1, 1)], QQ)
    f = minpoly(A, S)
    @test !iszero(f) && iszero(evaluate(f, A))

    A = from_dense([1 2; 3 4;], QQ)
    f = minpoly(A, S)
    @test !iszero(f) && iszero(evaluate(f, A))

    A = from_dense([1 0 0; 0 1 0; 0 0 0], QQ)
    f = minpoly(A, S)
    @test !iszero(f) && iszero(evaluate(f, A))

    A = from_dense([3 -3 2; -1 5 -2; -1 3 0], QQ)
    f = minpoly(A, S)
    @test !iszero(f) && iszero(evaluate(f, A))

    A = from_dense([0 0; 0 0;], QQ)
    @test_throws DomainError minpoly(A, S)

    # some random tests
    # - dense case
    sizes = [3, 5, 10]
    densities = [0.2, 0.5, 1.0]
    ZZ = GF(2^31 - 1)
    S, _ = ZZ["x"]
    for n in sizes
        for λ in densities
            A = random_sparsik((n, n), ZZ, density=λ)
            try
                f = minpoly(A, S)
                @test !iszero(f) && iszero(evaluate(f, A))
            catch DomainError
                !iszero(A) && rethrow()
            end
        end
    end

    # - highly singular case
    n = 5
    S, _ = QQ["x"]
    for _ in 1 : n^2
        i, j, k, r = rand(1:n^2, 4)
        # sum of four random unit matrices
        A = sum(unit_sparsik((n, n), pos, QQ) for pos in (i, j, k, r))
        f = minpoly(A, S)
        @test !iszero(f) && iszero(evaluate(f, A))
    end

end

#------------------------------------------------------------------------------

@testset "minpolynomial, basic deterministic" begin

    test_minimal_polynomial(__deterministic_simple_minpoly)

end

@info "OK"


text:
## All of these methods are defined in
## Huber, M., 2017. Optimal linear Bernoulli factories for small mean problems.
## Methodology and Computing in Applied Probability, 19(2), pp.631-645.

function _logistic(f::F, C::Float64, rng::RNG) where {F<:Function,
  RNG<:AbstractRNG}
  A::Float64 = randexp(rng)
  T::Float64 = randexp(rng)/C
  flips::Int64 = 0
  while T < A
    B::Bool = f()
    flips += 1
    if B
      return true, flips
    else
      T += randexp(rng)/C
    end
  end
  return false, flips
end

function _algorithmA(f::F, m::Int64, C::Float64, rng::RNG) where {F<:Function,
  RNG<:AbstractRNG}
  s::Int64 = 1
  flips::Int64 = 0
  while 0 < s < m
    B, newFlips = _logistic(f, C, rng)
    flips += newFlips
    s += 1-2*B
  end
  return s == 0, flips
end

function _highPowerLogistic(f::F, m::Int64, β::Float64, C::Float64,
  rng::RNG) where {F<:Function, RNG<:AbstractRNG}
  s::Int64 = 1
  flips::Int64 = 0
  α::Float64 = β*C
  while 0 < s <= m
    B, newFlips = _logistic(f, α, rng)
    flips += newFlips
    s += 2*B-1
  end
  return s == m+1, flips
end

function _algorithmB(f::F, ϵ::Float64, m::Int64, β::Float64, C::Float64,
  rng::RNG) where {F<:Function, RNG<:AbstractRNG}

  γ::Float64 = 1.0-(1.0-ϵ)*β
  α::Float64 = β*C
  flips::Int64 = 0
  while true
    B1, newFlips = _huber2017(f, α, γ, rng)
    flips += newFlips
    if !B1
      return false, flips
    end
    B2, newFlips = _highPowerLogistic(f, m-2, β, C, rng)
    flips += newFlips
    if B2
      return true, flips
    end
    m -= 1
  end
end

function _huber2017(f::F, C::Float64, ϵ::Float64, rng::RNG) where {F<:Function,
  RNG<:AbstractRNG}
  m::Int64 = ceil(Int64, 4.5/ϵ) + 1
  β::Float64 = 1.0 + 1.0/(m-1.0)
  flips::Int64 = 0
  B1, newFlips = _algorithmA(f, m, β*C, rng)
  flips += newFlips
  if !B1 return false, flips end
  if rand(rng) < 1.0/β
    return true, flips
  else
    B, newFlips = _algorithmB(f, ϵ, m, β, C, rng)
    flips += newFlips
    return B, flips
  end
end

function _huber2017Small(f::F, C::Float64, ϵ::Float64, rng::RNG) where {F<:Function,
  RNG<:AbstractRNG}
  @assert ϵ > 0.5
  M::Float64 = 1.0 - ϵ
  β::Float64 = 1.0/(1.0-2*M)
  Y, flips = _logistic(f, β*C, rng)
  if !Y return false, flips end
  B::Bool = rand(rng) < 1/β
  if B return true, flips end
  ϵ = 1.0 - β/(β-1.0)*M
  X, newFlips = _huber2017(f, C*β/(β-1.0), ϵ, rng)
  return X, flips + newFlips
end


text:
mutable struct JacobiPreconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
    A::ExtendableSparseMatrix{Tv,Ti}
    invdiag::Array{Tv,1}
    function JacobiPreconditioner{Tv,Ti}() where {Tv,Ti}
        p=new()
        p.invdiag=zeros(Tv,0)
        p
    end
end

"""
```
JacobiPreconditioner(;valuetype=Float64, indextype=Int64)
JacobiPreconditioner(matrix)
```

Jacobi preconditioner
"""
JacobiPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64)=JacobiPreconditioner{valuetype,indextype}()

function update!(precon::JacobiPreconditioner{Tv,Ti}) where {Tv,Ti}
    flush!(precon.A)
    cscmatrix=precon.A.cscmatrix
    if length(precon.invdiag)==0
        precon.invdiag=Array{Tv,1}(undef,cscmatrix.n)
    end
    invdiag=precon.invdiag
    n=cscmatrix.n
    @inbounds for i=1:n
        invdiag[i]=1.0/cscmatrix[i,i]
    end
    precon
end

function  LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::JacobiPreconditioner, v::AbstractArray{T,1} where T)
    invdiag=precon.invdiag
    n=length(invdiag)
    @inbounds @simd for i=1:n
        u[i]=invdiag[i]*v[i]
    end
end

function LinearAlgebra.ldiv!(precon::JacobiPreconditioner, v::AbstractArray{T,1} where T)
    ldiv!(v, precon, v)
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

# Parameters
ex = [1.;0.;0.]

l1 = 1.0
x, y = .1, .1

vert11 = [0.;0.;l1 / 2]
vert12 = -vert11
verts = [[vert11];[vert12]]

# Initial orientation
offset1 = pi / 8
offset2 = pi / 4
phi1 = pi / 8
q1 = UnitQuaternion(RotX(phi1))
qoff1 = UnitQuaternion(RotX(offset1))
qoff2 = UnitQuaternion(RotX(offset2))

# Links
origin = Origin{Float64}()
links = [Box(x, y, l1, l1, color = RGBA(1., 1., 0.)) for i = 1:4]

function fourbar(links, vertices, axis)
    j1 = EqualityConstraint(Revolute(links[1], links[2], axis; p1=vertices[1], p2=vertices[2]))
    j2 = EqualityConstraint(Revolute(links[2], links[3], axis; p1=vertices[3], p2=vertices[2]), Cylindrical(links[2], links[4], axis; p1=vertices[2], p2=vertices[2]))
    j3 = EqualityConstraint(Revolute(links[4], links[5], axis; p1=vertices[3], p2=vertices[2]))
    j4 = EqualityConstraint(Revolute(links[3], links[5], axis; p1=vertices[3], p2=vertices[3]))

    return j1, j2, j3, j4
end

function initfourbar!(mechanism, links, vertices, Δq1, Δq2)
    setPosition!(links[1], links[2], p1 = vertices[1], p2 = vertices[2], Δq = Δq1)
    setPosition!(links[2], links[3], p1 = vertices[3], p2 = vertices[2], Δq = inv(Δq2) * inv(Δq2))
    setPosition!(links[2], links[4], p1 = vertices[2], p2 = vertices[2], Δq = inv(Δq2) * inv(Δq2))
    setPosition!(links[4], links[5], p1 = vertices[3], p2 = vertices[2], Δq = Δq2 * Δq2)

    return
end

# Constraints
constraints = [fourbar([origin;links[1:4]], [[zeros(3)];verts], ex)...]

mech = Mechanism(origin, links, constraints)

initfourbar!(mech, [origin;links[1:4]], [[zeros(3)];verts], q1 * qoff2, q1)

function control!(mechanism,k)
    if k==75
        deactivate!(mechanism,constraints[3].id)
    end
end

storage = simulate!(mech, 10., control!, record = true)
visualize(mech, storage)


text:
"""
	This File is used to calculate the benchmarks.
    1, Influence Maximization
    2, Max Degree
    3, Random Sampling
"""

include("../src/Graph.jl")
include("../src/Tools.jl")
include("../src/Algorithm.jl")
include("../src/Sampling.jl")
include("../src/Methods.jl")

using StatsBase
using JSON
using JLD2

"""
Set Seed Nodes & change innate opinions
Set parameters.
k = 50
epsilon : s + epsilon
delta : edge probability parameter
beta: threshold
json struct js
"""

epsilon = 0.1
delta = 0.5
beta = 0.5
tags = []

jsdict = Dict{String, Any}()

# Read Graph
buf = split(ARGS[1], ',')
fileName = string("../data/all/", buf[1], ".jld2")

(name, n, m, V, E) = load(fileName, "name", "n", "m", "V", "E")
G = Graph(n, m, V, E)

merge!(jsdict, Dict("name"=>name, "n"=>n, "m"=>m))

if ARGS[2] == "1"
    #Weighted Cascade
    push!(tags, "WC")
    E_p = load(fileName, "E_WC")
elseif ARGS[2] == "2"
    #Trivalency
    push!(tags, "TM") # Trivalency Model
    E_p = load(fileName, "E_TM")
end
G2 = Graph(n, m, V, E_p)

# Before Making changes
if ARGS[3] == "1"
    # Uniform distribution
    push!(tags, "Uniform")
    s = load(fileName, "s_Uni")
elseif ARGS[3] == "2"
    # Exponential distribution
    push!(tags, "Exponential")
    s = load(fileName, "s_Exp")
elseif ARGS[3] == "3"
    # Power-law distribution
    push!(tags, "Power-law")
    s = load(fileName, "s_Pow")
end

# Different parameters to select k
if ARGS[4] == "1"
    push!(tags, "0.005")
    k = load(fileName, "k1")
elseif ARGS[4] == "2"
    push!(tags, "0.015")
    k = load(fileName, "k2")
elseif ARGS[4] == "3"
    push!(tags, "0.02")
    k = load(fileName, "k3")
end

merge!(jsdict, Dict("k"=>k))
merge!(jsdict, Dict("model"=>tags[1], "s_dist"=>tags[2], "k_par"=>tags[3]))

name = ["z_sum", "aci", "ad", "ap", "aidc"]

##################### Above are configs
#####################
# The only difference for marketing setting and backfire setting is that Δs is different

if ARGS[5] == "1"
    Δs = delta_s_marketing(s, epsilon)
elseif ARGS[5] == "2"
    Δs = delta_s_backfire(s, epsilon, beta)
end

#######################
####################### Baselines

base_labels = ["max_influence_", "high_degree_", "random_"]
base_time = []
base_p = []

# rr-algorithm
(_, seeds, _, T) = ris_influence_max(G2, k, delta) ## RR to select seed nodes
(_, p_rr) = cascadeIC(G2, seeds, round(Int, 2*(10*10)*log(2*n)), delta) ## direct calculation, just to check
push!(base_time, T)
push!(base_p, p_rr)

# high-degree-algorithm
(T, seeds_high_degree) = choose_high_degree(G2, k)
(_, p_degree) = cascadeIC(G2, seeds_high_degree, round(Int, 2*(10*10)*log(2*n)), delta) ## direct calculation, just to check
push!(base_time, T)
push!(base_p, p_degree)

# random
(T, p_random) = random_spread(G2, k, delta)
push!(base_time, T)
push!(base_p, p_random)

for i in 1:length(base_labels)
    nn = base_labels[i].*name
    base_t = base_labels[i]*"time"

    s2 = Δs .* base_p[i] .+ s

    _, aci, ad, ap, aidc = Approx(G, s2)
    z_sum = sum(s2)

    #println("z_sum gain is ", z_sum - sum(s))

    result = [z_sum, aci, ad, ap, aidc]

    merge!(jsdict, Dict(nn[1]=>result[1], nn[2]=>result[2], nn[3]=>result[3], nn[4]=>result[4], nn[5]=>result[5]))
    merge!(jsdict, Dict(base_t=>base_time[i]))
end

js = JSON.json(jsdict)

println(jsdict)

if ARGS[5] == "1"
    open("out/log_benchmark_marketing.json", "a+") do ff
        write(ff, js)
        write(ff, "\n")
    end
elseif ARGS[5] == "2"
    open("out/log_benchmark_backfire.json", "a+") do ff
        write(ff, js)
        write(ff, "\n")
    end
end


text:
# extend hashing tests from Julia Base to Float128

hashtypes = Any[
    Bool,
    Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64, Float128,
    Rational{Int8}, Rational{UInt8}, Rational{Int16}, Rational{UInt16},
    Rational{Int32}, Rational{UInt32}, Rational{Int64}, Rational{UInt64}
]
hashvals = vcat(
    typemin(Int64),
    -Int64(maxintfloat(Float64)) .+ Int64[-4:1;],
    typemin(Int32),
    -Integer(maxintfloat(Float32)) .+ (-4:1),
    -2:2,
    Integer(maxintfloat(Float32)) .+ (-1:4),
    typemax(Int32),
    Int64(maxintfloat(Float64)) .+ Int64[-1:4;],
    typemax(Int64),
)

function coerce(T::Type, x)
    if T<:Rational
        convert(T, coerce(typeof(numerator(zero(T))), x))
    elseif !(T<:Integer)
        convert(T, x)
    else
        x % T
    end
end

@testset "hash 0,1" begin
  for x = hashvals,
      a = coerce(Float128, x)
      @test hash(a,zero(UInt)) == invoke(hash, Tuple{Real, UInt}, a, zero(UInt))
      @test hash(a,one(UInt)) == invoke(hash, Tuple{Real, UInt}, a, one(UInt))
  end
end

@testset "hashing" begin
  for T = hashtypes,
    S = hashtypes,
    x = hashvals,
    a = coerce(T, x),
    b = coerce(S, x)
    T == Float128 || S == Float128 || continue # other types are not our problem
    # equality must imply same hash; converse is not guaranteed
    @test !(isequal(a,b) && (hash(a) != hash(b)))
  end
end


text:
# Original Java version Copyright 2005 Google Inc.
# Julia translation Copyright 2022 Ted Dunning
#
# 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.

begin
    # Create some standard intervals to use in the tests.  These include the
    # empty and full intervals, intervals containing a single point, and
    # intervals spanning one or more "quadrants" which are numbered as follows:
    #    quad1 == [0, π/2]
    #    quad2 == [π/2, π]
    #    quad3 == [-π, -π/2]
    #    quad4 == [-π/2, 0]
    emptyInterval = empty(S1Interval);
    fullInterval = full(S1Interval);

    # Single-point intervals:
    zero =  S1Interval(0, 0);
    pi2 =  S1Interval(π/2, π/2);
    pi = S1Interval(π, π);
    mipi = S1Interval(-π, -π) # Same as "pi" after normalization.
    mipi2 = S1Interval(-π/2, -π/2);

    # Single quadrants:
    quad1 = S1Interval(0, π/2);
    quad2 = S1Interval(π/2, -π);
    quad3 = S1Interval(π, -π/2);
    quad4 = S1Interval(-π/2, 0);
    # Quadrant pairs:
    quad12 = S1Interval(0, -π);
    quad23 = S1Interval(π/2, -π/2);
    quad34 = S1Interval(-π, 0);
    quad41 = S1Interval(-π/2, π/2);
    # Quadrant triples:
    quad123 = S1Interval(0, -π/2);
    quad234 = S1Interval(π/2, 0);
    quad341 = S1Interval(π, π/2);
    quad412 = S1Interval(-π/2, -π);
    # Small intervals around the midpoints between quadrants, such that
    # the center of each interval is offset slightly CCW from the midpoint.
    mid12 = S1Interval(π/2 - 0.01, π/2 + 0.02);
    mid23 = S1Interval(π - 0.01, -π + 0.02);
    mid34 = S1Interval(-π/2 - 0.01, -π/2 + 0.02);
    mid41 = S1Interval(-0.01, 0.02);

    # Spot-check the constructors and accessors.
    @test 0.0 == quad12.lo
    @test Float64(π) == quad12.hi
    @test Float64(π) == pi.lo
    @test Float64(π) == pi.hi

    # Check that [-π, -π] is normalized to [π, π].
    @test S2.S2_M_PI == mipi.lo
    @test S2.S2_M_PI == mipi.hi
    @test Float64(S2.S2_M_PI/2) == quad23.lo
    @test Float64(-S2.S2_M_PI/2) == quad23.hi

    # is_valid(), is_empty(), is_full(), is_inverted()
    @test isvalid(zero) && !isempty(zero) && !isfull(zero)
    @test isvalid(emptyInterval) && isempty(emptyInterval) && !isfull(emptyInterval)
    @test isinverted(emptyInterval)
    @test isvalid(fullInterval) && !isempty(fullInterval) && isfull(fullInterval)
    @test !isempty(quad12) && !isfull(quad12) && !isinverted(quad12)
    @test !isempty(quad23) && !isfull(quad23) && isinverted(quad23)
    @test isvalid(pi) && !isempty(pi) && !isinverted(pi)
    @test isvalid(mipi) && !isempty(mipi) && !isinverted(mipi)

    @test S2.S2_M_PI/2 == center(quad12)
    @test 3.0 - S2.S2_M_PI == center(S1Interval(3.1, 2.9))
    @test S2.S2_M_PI - 3.0 == center(S1Interval(-2.9, -3.1))
    @test S2.S2_M_PI == center(S1Interval(2.1, -2.1))
    @test S2.S2_M_PI == center(pi)
    @test S2.S2_M_PI == center(mipi)
    @test S2.S2_M_PI == abs(center(quad23))
    @test 0.75 * S2.S2_M_PI == center(quad123)

    @test S2.S2_M_PI == length(quad12)
    @test 0.0 == length(pi)
    @test 0.0 == length(mipi)
    @test 1.5 * S2.S2_M_PI == length(quad123)
    @test S2.S2_M_PI == abs(length(quad23))
    @test 2 * S2.S2_M_PI == length(fullInterval)
    @test length(emptyInterval) < 0

    @test isfull(complement(emptyInterval))
    @test isempty(complement(fullInterval))
    @test isfull(complement(pi))
    @test isfull(complement(mipi))
    @test isfull(complement(zero))
    @test complement(quad12) ≈ quad34
    @test complement(quad34) ≈ quad12
    @test complement(quad123) ≈ quad4

    @test !contains(emptyInterval, 0) && !contains(emptyInterval, S2.S2_M_PI) && !contains(emptyInterval, -S2.S2_M_PI)
    @test !interiorcontains(emptyInterval, S2.S2_M_PI) && !interiorcontains(emptyInterval, -S2.S2_M_PI)
    @test contains(fullInterval, 0) && contains(fullInterval, S2.S2_M_PI) && contains(fullInterval, -S2.S2_M_PI)
    @test interiorcontains(fullInterval, S2.S2_M_PI) && interiorcontains(fullInterval, -S2.S2_M_PI)
    @test contains(quad12, 0) && contains(quad12, S2.S2_M_PI) && contains(quad12, -S2.S2_M_PI)
    @test interiorcontains(quad12, S2.S2_M_PI/2) && !interiorcontains(quad12, 0)
    @test !interiorcontains(quad12, S2.S2_M_PI) && !interiorcontains(quad12, -S2.S2_M_PI)
    @test contains(quad23, S2.S2_M_PI/2) && contains(quad23, -S2.S2_M_PI/2)
    @test contains(quad23, S2.S2_M_PI) && contains(quad23, -S2.S2_M_PI)
    @test !contains(quad23, 0)
    @test !interiorcontains(quad23, S2.S2_M_PI/2) && !interiorcontains(quad23, -S2.S2_M_PI/2)
    @test interiorcontains(quad23, S2.S2_M_PI) && interiorcontains(quad23, -S2.S2_M_PI)
    @test !interiorcontains(quad23, 0)
    @test contains(pi, S2.S2_M_PI) && contains(pi, -S2.S2_M_PI) && !contains(pi, 0)
    @test !interiorcontains(pi, S2.S2_M_PI) && !interiorcontains(pi, -S2.S2_M_PI)
    @test contains(mipi, S2.S2_M_PI) && contains(mipi, -S2.S2_M_PI) && !contains(mipi, 0)
    @test !interiorcontains(mipi, S2.S2_M_PI) && !interiorcontains(mipi, -S2.S2_M_PI)
    @test contains(zero, 0) && !interiorcontains(zero, 0)

    function testIntervalOps(x::S1Interval, y::S1Interval,
                             expectedRelation::String,
                             expectedUnion::S1Interval, expectedIntersection::S1Interval)
        # Test all of the interval operations on the given pair of intervals.
        # "expectedRelation" is a sequence of "T" and "F" characters corresponding
        # to the expected results of Contains(), InteriorContains(), Intersects(),
        # and InteriorIntersects() respectively.
        @test contains(x, y) == (expectedRelation[1] == 'T')
        @test interiorcontains(x, y) == (expectedRelation[2] == 'T')
        @test intersects(x, y) == (expectedRelation[3] == 'T')
        @test interiorintersects(x, y) == (expectedRelation[4] == 'T')

        # bounds() returns a const reference to a member variable, so we need to
        # make a copy when invoking it on a temporary object.
        @test expectedUnion == union(x, y)
        @test expectedIntersection == intersection(x, y)

        @test contains(x, y) == (union(x, y) == x)
        @test intersects(x, y) == !isempty(intersection(x, y))

        if y.lo == y.hi
            r = addPoint(x, y.lo);
            @test expectedUnion == r
        end
    end

    begin
        # Contains(S1Interval), InteriorContains(S1Interval),
        # Intersects(), InteriorIntersects(), Union(), Intersection()
        testIntervalOps(emptyInterval, emptyInterval, "TTFF", emptyInterval, emptyInterval);
        testIntervalOps(emptyInterval, fullInterval, "FFFF", fullInterval, emptyInterval);
        testIntervalOps(emptyInterval, zero, "FFFF", zero, emptyInterval);
        testIntervalOps(emptyInterval, pi, "FFFF", pi, emptyInterval);
        testIntervalOps(emptyInterval, mipi, "FFFF", mipi, emptyInterval);

        testIntervalOps(fullInterval, emptyInterval, "TTFF", fullInterval, emptyInterval);
        testIntervalOps(fullInterval, fullInterval, "TTTT", fullInterval, fullInterval);
        testIntervalOps(fullInterval, zero, "TTTT", fullInterval, zero);
        testIntervalOps(fullInterval, pi, "TTTT", fullInterval, pi);
        testIntervalOps(fullInterval, mipi, "TTTT", fullInterval, mipi);
        testIntervalOps(fullInterval, quad12, "TTTT", fullInterval, quad12);
        testIntervalOps(fullInterval, quad23, "TTTT", fullInterval, quad23);

        testIntervalOps(zero, emptyInterval, "TTFF", zero, emptyInterval);
        testIntervalOps(zero, fullInterval, "FFTF", fullInterval, zero);
        testIntervalOps(zero, zero, "TFTF", zero, zero);
        testIntervalOps(zero, pi, "FFFF", S1Interval(0, S2.S2_M_PI), emptyInterval);
        testIntervalOps(zero, pi2, "FFFF", quad1, emptyInterval);
        testIntervalOps(zero, mipi, "FFFF", quad12, emptyInterval);
        testIntervalOps(zero, mipi2, "FFFF", quad4, emptyInterval);
        testIntervalOps(zero, quad12, "FFTF", quad12, zero);
        testIntervalOps(zero, quad23, "FFFF", quad123, emptyInterval);

        testIntervalOps(pi2, emptyInterval, "TTFF", pi2, emptyInterval);
        testIntervalOps(pi2, fullInterval, "FFTF", fullInterval, pi2);
        testIntervalOps(pi2, zero, "FFFF", quad1, emptyInterval);
        testIntervalOps(pi2, pi, "FFFF", S1Interval(S2.S2_M_PI/2, S2.S2_M_PI), emptyInterval);
        testIntervalOps(pi2, pi2, "TFTF", pi2, pi2);
        testIntervalOps(pi2, mipi, "FFFF", quad2, emptyInterval);
        testIntervalOps(pi2, mipi2, "FFFF", quad23, emptyInterval);
        testIntervalOps(pi2, quad12, "FFTF", quad12, pi2);
        testIntervalOps(pi2, quad23, "FFTF", quad23, pi2);

        testIntervalOps(pi, emptyInterval, "TTFF", pi, emptyInterval);
        testIntervalOps(pi, fullInterval, "FFTF", fullInterval, pi);
        testIntervalOps(pi, zero, "FFFF", S1Interval(S2.S2_M_PI, 0), emptyInterval);
        testIntervalOps(pi, pi, "TFTF", pi, pi);
        testIntervalOps(pi, pi2, "FFFF", S1Interval(S2.S2_M_PI/2, S2.S2_M_PI), emptyInterval);
        testIntervalOps(pi, mipi, "TFTF", pi, pi);
        testIntervalOps(pi, mipi2, "FFFF", quad3, emptyInterval);
        testIntervalOps(pi, quad12, "FFTF", S1Interval(0, S2.S2_M_PI), pi);
        testIntervalOps(pi, quad23, "FFTF", quad23, pi);

        testIntervalOps(mipi, emptyInterval, "TTFF", mipi, emptyInterval);
        testIntervalOps(mipi, fullInterval, "FFTF", fullInterval, mipi);
        testIntervalOps(mipi, zero, "FFFF", quad34, emptyInterval);
        testIntervalOps(mipi, pi, "TFTF", mipi, mipi);
        testIntervalOps(mipi, pi2, "FFFF", quad2, emptyInterval);
        testIntervalOps(mipi, mipi, "TFTF", mipi, mipi);
        testIntervalOps(mipi, mipi2, "FFFF", S1Interval(-S2.S2_M_PI, -S2.S2_M_PI/2), emptyInterval);
        testIntervalOps(mipi, quad12, "FFTF", quad12, mipi);
        testIntervalOps(mipi, quad23, "FFTF", quad23, mipi);

        testIntervalOps(quad12, emptyInterval, "TTFF", quad12, emptyInterval);
        testIntervalOps(quad12, fullInterval, "FFTT", fullInterval, quad12);
        testIntervalOps(quad12, zero, "TFTF", quad12, zero);
        testIntervalOps(quad12, pi, "TFTF", quad12, pi);
        testIntervalOps(quad12, mipi, "TFTF", quad12, mipi);
        testIntervalOps(quad12, quad12, "TFTT", quad12, quad12);
        testIntervalOps(quad12, quad23, "FFTT", quad123, quad2);
        testIntervalOps(quad12, quad34, "FFTF", fullInterval, quad12);

        testIntervalOps(quad23, emptyInterval, "TTFF", quad23, emptyInterval);
        testIntervalOps(quad23, fullInterval, "FFTT", fullInterval, quad23);
        testIntervalOps(quad23, zero, "FFFF", quad234, emptyInterval);
        testIntervalOps(quad23, pi, "TTTT", quad23, pi);
        testIntervalOps(quad23, mipi, "TTTT", quad23, mipi);
        testIntervalOps(quad23, quad12, "FFTT", quad123, quad2);
        testIntervalOps(quad23, quad23, "TFTT", quad23, quad23);
        testIntervalOps(quad23, quad34, "FFTT", quad234, S1Interval(-S2.S2_M_PI, -S2.S2_M_PI/2));

        testIntervalOps(quad1, quad23, "FFTF", quad123, S1Interval(S2.S2_M_PI/2, S2.S2_M_PI/2));
        testIntervalOps(quad2, quad3, "FFTF", quad23, mipi);
        testIntervalOps(quad3, quad2, "FFTF", quad23, pi);
        testIntervalOps(quad2, pi, "TFTF", quad2, pi);
        testIntervalOps(quad2, mipi, "TFTF", quad2, mipi);
        testIntervalOps(quad3, pi, "TFTF", quad3, pi);
        testIntervalOps(quad3, mipi, "TFTF", quad3, mipi);

        testIntervalOps(quad12, mid12, "TTTT", quad12, mid12);
        testIntervalOps(mid12, quad12, "FFTT", quad12, mid12);

        quad12eps = S1Interval(quad12.lo, mid23.hi);
        quad2hi = S1Interval(mid23.lo, quad12.hi);
        testIntervalOps(quad12, mid23, "FFTT", quad12eps, quad2hi);
        testIntervalOps(mid23, quad12, "FFTT", quad12eps, quad2hi);

        # This test checks that the union of two disjoint intervals is the smallest
        # interval that contains both of them.  Note that the center of "mid34"
        # slightly CCW of -S2.S2_M_PI/2 so that there is no ambiguity about the result.
        quad412eps = S1Interval(mid34.lo, quad12.hi)
        testIntervalOps(quad12, mid34, "FFFF", quad412eps, emptyInterval);
        testIntervalOps(mid34, quad12, "FFFF", quad412eps, emptyInterval);

        quadeps12 = S1Interval(mid41.lo, quad12.hi);
        quad1lo = S1Interval(quad12.lo, mid41.hi);
        testIntervalOps(quad12, mid41, "FFTT", quadeps12, quad1lo);
        testIntervalOps(mid41, quad12, "FFTT", quadeps12, quad1lo);

        quad2lo = S1Interval(quad23.lo, mid12.hi);
        quad3hi = S1Interval(mid34.lo, quad23.hi);
        quadeps23 = S1Interval(mid12.lo, quad23.hi);
        quad23eps = S1Interval(quad23.lo, mid34.hi);
        quadeps123 = S1Interval(mid41.lo, quad23.hi);
        testIntervalOps(quad23, mid12, "FFTT", quadeps23, quad2lo);
        testIntervalOps(mid12, quad23, "FFTT", quadeps23, quad2lo);
        testIntervalOps(quad23, mid23, "TTTT", quad23, mid23);
        testIntervalOps(mid23, quad23, "FFTT", quad23, mid23);
        testIntervalOps(quad23, mid34, "FFTT", quad23eps, quad3hi);
        testIntervalOps(mid34, quad23, "FFTT", quad23eps, quad3hi);
        testIntervalOps(quad23, mid41, "FFFF", quadeps123, emptyInterval);
        testIntervalOps(mid41, quad23, "FFFF", quadeps123, emptyInterval);
    end

    begin
    @test addPoint(emptyInterval, 0) == zero
    @test addPoint(emptyInterval, S2.S2_M_PI) == pi
    @test addPoint(emptyInterval, -S2.S2_M_PI) == mipi
    @test addPoint(addPoint(emptyInterval, S2.S2_M_PI), -S2.S2_M_PI) == pi
    @test addPoint(addPoint(emptyInterval, -S2.S2_M_PI), S2.S2_M_PI) == pi
    @test addPoint(addPoint(emptyInterval, mid12.lo), mid12.hi) == mid12
    @test addPoint(addPoint(emptyInterval, mid23.lo), mid23.hi) == mid23
    @test addPoint(addPoint(quad1, -0.9 * S2.S2_M_PI), -S2.S2_M_PI/2) == quad123
    @test isfull(addPoint(fullInterval, 0))
    @test isfull(addPoint(fullInterval, S2.S2_M_PI))
    @test isfull(addPoint(fullInterval, -S2.S2_M_PI))
    end

    begin
        r = S1Interval(-S2.S2_M_PI, -S2.S2_M_PI);
        @test S2.S2_M_PI == clampPoint(r, -S2.S2_M_PI)
        @test S2.S2_M_PI == clampPoint(r, 0)
        r = S1Interval(0, S2.S2_M_PI);
        @test 0.1 == clampPoint(r, 0.1)
        @test 0.0 == clampPoint(r, -S2.S2_M_PI/2 + 1e-15)
        @test S2.S2_M_PI == clampPoint(r, -S2.S2_M_PI/2 - 1e-15)
        r = S1Interval(S2.S2_M_PI - 0.1, -S2.S2_M_PI + 0.1);
        @test S2.S2_M_PI == clampPoint(r, S2.S2_M_PI)
        @test S2.S2_M_PI - 0.1 == clampPoint(r, 1e-15)
        @test -S2.S2_M_PI + 0.1 == clampPoint(r, -1e-15)
        @test 0.0 == clampPoint(full(S1Interval), 0)
        @test S2.S2_M_PI == clampPoint(full(S1Interval), S2.S2_M_PI)
        @test S2.S2_M_PI == clampPoint(full(S1Interval), -S2.S2_M_PI)
    end

    @test fromPointPair(-S2.S2_M_PI, S2.S2_M_PI) == pi
    @test fromPointPair(S2.S2_M_PI, -S2.S2_M_PI) == pi
    @test fromPointPair(mid34.hi, mid34.lo) == mid34
    @test fromPointPair(mid23.lo, mid23.hi) == mid23

    @test expanded(emptyInterval, 1) == emptyInterval
    @test expanded(fullInterval, 1) == fullInterval
    @test expanded(zero, 1) == S1Interval(-1, 1)
    @test expanded(mipi, 0.01) == S1Interval(S2.S2_M_PI - 0.01, -S2.S2_M_PI + 0.01)
    @test expanded(pi, 27) == fullInterval
    @test expanded(pi, S2.S2_M_PI/2) == quad23
    @test expanded(pi2, S2.S2_M_PI/2) == quad12
    @test expanded(mipi2, S2.S2_M_PI/2) == quad34

    @test expanded(emptyInterval, -1) == emptyInterval
    @test expanded(fullInterval, -1) == fullInterval
    @test expanded(quad123, -27) == emptyInterval
    @test expanded(quad234, -27) == emptyInterval
    @test expanded(quad123, -S2.S2_M_PI/2) == quad2
    @test expanded(quad341, -S2.S2_M_PI/2) == quad4
    @test expanded(quad412, -S2.S2_M_PI/2) == quad1

    begin
        # Choose two values kLo and kHi such that it's okay to shift an endpoint by
        # kLo (i.e., the resulting interval is equivalent) but not by kHi.
        kLo = 4 * eps(Float64) # < maxError default
        kHi = 6 * eps(Float64) # > maxError default

        # Empty intervals.
        @test emptyInterval ≈ emptyInterval
        @test zero ≈ emptyInterval && emptyInterval ≈ zero
        @test pi ≈ emptyInterval && emptyInterval ≈ pi
        @test mipi ≈ emptyInterval && emptyInterval ≈ mipi
        @test !(emptyInterval ≈ fullInterval)
        @test emptyInterval ≈ S1Interval(1, 1 + 2 * kLo)
        @test !(emptyInterval ≈ S1Interval(1, 1 + 2 * kHi))
        @test S1Interval(S2.S2_M_PI - kLo, -S2.S2_M_PI + kLo) ≈ emptyInterval

        # Full intervals.
        @test fullInterval ≈ fullInterval
        @test !(fullInterval ≈ emptyInterval)
        @test !(fullInterval ≈ zero)
        @test !(fullInterval ≈ pi)
        @test fullInterval ≈ S1Interval(kLo, -kLo)
        @test !(fullInterval ≈ S1Interval(2 * kHi, 0))
        @test S1Interval(-S2.S2_M_PI + kLo, S2.S2_M_PI - kLo) ≈ fullInterval
        @test !(S1Interval(-S2.S2_M_PI, S2.S2_M_PI - 2 * kHi) ≈ fullInterval)

        # Singleton intervals.
        @test pi ≈ pi && mipi ≈ pi
        @test pi ≈ S1Interval(S2.S2_M_PI - kLo, S2.S2_M_PI - kLo)
        @test !(pi ≈ S1Interval(S2.S2_M_PI - kHi, S2.S2_M_PI - kHi))
        @test pi ≈ S1Interval(S2.S2_M_PI - kLo, -S2.S2_M_PI + kLo)
        @test !(pi ≈ S1Interval(S2.S2_M_PI - kHi, -S2.S2_M_PI))
        @test !(zero ≈ pi)
        @test union(union(pi, mid12), zero) ≈ quad12
        @test intersection(quad2, quad3) ≈ pi
        @test intersection(quad3, quad2) ≈ pi

        # Intervals whose corresponding endpoints are nearly the same but where the
        # endpoints are in opposite order (i.e., inverted intervals).
        @test !(S1Interval(0, kLo) ≈ S1Interval(kLo, 0))
        @test !(S1Interval(S2.S2_M_PI - 0.5 * kLo, -S2.S2_M_PI + 0.5 * kLo) ≈ S1Interval(-S2.S2_M_PI + 0.5 * kLo, S2.S2_M_PI - 0.5 * kLo))

        # Other intervals.
        @test S1Interval(1 - kLo, 2 + kLo) ≈ S1Interval(1, 2)
        @test S1Interval(1 + kLo, 2 - kLo) ≈ S1Interval(1, 2)
        @test S1Interval(2 - kLo, 1 + kLo) ≈ S1Interval(2, 1)
        @test S1Interval(2 + kLo, 1 - kLo) ≈ S1Interval(2, 1)
        @test !(S1Interval(1 - kHi, 2 + kLo) ≈ S1Interval(1, 2))
        @test !(S1Interval(1 + kHi, 2 - kLo) ≈ S1Interval(1, 2))
        @test !(S1Interval(2 - kHi, 1 + kLo) ≈ S1Interval(2, 1))
        @test !(S1Interval(2 + kHi, 1 - kLo) ≈ S1Interval(2, 1))
        @test !(S1Interval(1 - kLo, 2 + kHi) ≈ S1Interval(1, 2))
        @test !(S1Interval(1 + kLo, 2 - kHi) ≈ S1Interval(1, 2))
        @test !(S1Interval(2 - kLo, 1 + kHi) ≈ S1Interval(2, 1))
        @test !(S1Interval(2 + kLo, 1 - kHi) ≈ S1Interval(2, 1))
    end

    begin
        @test 0.0 == getDirectedHausdorffDistance(emptyInterval, emptyInterval)
        @test 0.0 == getDirectedHausdorffDistance(emptyInterval, mid12)
        @test S2.S2_M_PI == getDirectedHausdorffDistance(mid12, emptyInterval)

        @test 0.0 == getDirectedHausdorffDistance(quad12, quad123)
        in = S1Interval(3.0, -3.0) # an interval whose complement center is 0.
        @test 3.0 == getDirectedHausdorffDistance(S1Interval(-0.1, 0.2), in)
        @test 3.0 - 0.1 == getDirectedHausdorffDistance(S1Interval(0.1, 0.2), in)
        @test 3.0 - 0.1 == getDirectedHausdorffDistance(S1Interval(-0.2, -0.1), in)
    end
end