text:
#src # DGSEM with flux differencing

# This tutorial starts with a presentation of the weak formulation of the discontinuous Galerkin
# spectral element method (DGSEM) in order to fix the notation of the used operators.
# Then, the DGSEM formulation with flux differencing (split form DGSEM) and its implementation in
# [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) is shown.

# We start with the one-dimensional conservation law
# ```math
# u_t + f(u)_x = 0, \qquad t\in \mathbb{R}^+, x\in\Omega
# ```
# with the physical flux $f$.

# We split the domain $\Omega$ into elements $K$ with center $x_K$ and size $\Delta x$. With the
# transformation mapping $x(\xi)=x_K + \frac{\Delta x}{2} \xi$ we can transform the reference element
# $[-1,1]$ to every physical element. So, the equation can be restricted to the reference element using the
# determinant of the Jacobian matrix of the transformation mapping
# $J=\frac{\partial x}{\partial \xi}=\frac{\Delta x}{2}$.
# ```math
# J u_t + f(u)_{\xi} = 0, \qquad t\in \mathbb{R}^+, \xi\in [-1,1]
# ```

# ## The weak form of the DGSEM
# We consider the so-called discontinuous Galerkin spectral element method (DGSEM) with collocation.
# It results from choosing a nodal DG ansatz using $N+1$ Gauss-Lobatto nodes $\xi_i$ in $[-1,1]$
# with matching interpolation weights $w_i$, which are used for numerical integration and interpolation with
# the Lagrange polynomial basis $l_i$ of degree $N$. The Lagrange functions are created with those nodes and
# hence fulfil a Kronecker property at the GL nodes.
# The weak formulation of the DGSEM for one element is
# ```math
# J \underline{\dot{u}}(t) = - M^{-1} B \underline{f}^* + M^{-1} D^T M \underline{f}
# ```
# where $\underline{u}=(u_0, u_1, \dots, u_N)^T\in\mathbb{R}^{N+1}$ is the collected pointwise evaluation
# of $u$ at the discretization nodes and $\dot{u} = \partial u / \partial t = u_t$ is the temporal derivative.
# The nodal values of the flux function $f$ results with collocation in $\underline{f}$, since
# $\underline{f}_j=f(\underline{u}_j)$. Moreover, we got the numerical flux $f^*=f^*(u^-, u^+)$.

# We will now have a short overview over the operators we used.

# The **derivative matrix** $D\in\mathbb{R}^{(N+1)\times (N+1)}$ mimics a spatial derivation on a
# discrete level with $\underline{f}_x \approx D \underline{f}$. It is defined by $D_{ij} = l_j'(\xi_i)$.

# The diagonal **mass matrix** $M$ is defined by $M_{ij}=\langle l_j, l_i\rangle_N$ with the
# numerical scalar product $\langle \cdot, \cdot\rangle_N$ defined for functions $f$ and $g$ by
# ```math
# \langle f, g\rangle_N := \int_{-1, N}^1 f(\xi) g(\xi) d\xi := \sum_{k=0}^N f(\xi_k) g(\xi_k) w_k.
# ```
# The multiplication by $M$ matches a discrete integration
# ```math
#   \int_{-1}^1 f(\xi) \underline{l}(\xi) d\xi \approx M \underline{f},
# ```

# The **boundary matrix** $B=\text{diag}([-1, 0,..., 0, 1])$ represents an evaluation of a
# function at the boundaries $\xi_0=-1$ and $\xi_N=1$.

# For these operators the following property holds:
# ```math
#   M D + (M D)^T = B.
# ```
# This is called the summation-by-parts (SBP) property since it mimics integration by parts on a
# discrete level ([Gassner (2013)](https://doi.org/10.1137/120890144)).

# The explicit definitions of the operators and the contruction of the 1D algorithm can be found
# for instance in the tutorial [introduction to DG methods](@ref scalar_linear_advection_1d)
# or in more detail in [Kopriva (2009)](https://link.springer.com/book/10.1007/978-90-481-2261-5).

# This property shows the equivalence between the weak form and the following strong formulation
# of the DGSEM.
# ```math
# \begin{align*}
# J \underline{\dot{u}}(t)
# &= - M^{-1} B \underline{f}^* + M^{-1} D^T M \underline{f}\\[5pt]
# &= - M^{-1} B \underline{f}^* + M^{-1} (B - MD) \underline{f}\\[5pt]
# &= - M^{-1} B (\underline{f}^* - \underline{f}) - D \underline{f}
# \end{align*}
# ```
# More information about the equivalence you can find in [Kopriva, Gassner (2010)](https://doi.org/10.1007/s10915-010-9372-3).

# ## DGSEM with flux differencing
# When using the diagonal SBP property it is possible to rewrite the application of the derivative
# operator $D$ in the calculation of the volume integral into a subcell based finite volume type
# differencing formulation ([Fisher, Carpenter (2013)](https://doi.org/10.1016/j.jcp.2013.06.014)).
# We replace $D \underline{f}$ in the strong form by $2D \underline{f}_{vol}(u^-, u^+)$ with
# the consistent two-point volume flux $f_{vol}$ and receive the DGSEM formulation with flux differencing
# (split form DGSEM) ([Gassner, Winters, Kopriva (2016)](https://doi.org/10.1016/j.jcp.2016.09.013)).

# ```math
# \begin{align*}
# J \underline{\dot{u}}(t) &= - M^{-1} B (\underline{f}^* - \underline{f}) - 2D \underline{f}_{vol}(u^-, u^+)\\[5pt]
# &= - M^{-1} B (\underline{f}^* - \underline{f}_{vol}(\underline{u}, \underline{u})) - 2D \underline{f}_{vol}(u^-, u^+)\\[5pt]
# &= - M^{-1} B \underline{f}_{sur}^* - (2D - M^{-1} B) \underline{f}_{vol}\\[5pt]
# &= - M^{-1} B \underline{f}_{sur}^* - D_{split} \underline{f}_{vol}
# \end{align*}
# ```
# This formulation is in a weak form type formulation and can be implemented by using the derivative
# split matrix $D_{split}=(2D-M^{-1}B)$ and two different fluxes. We divide between the surface
# flux $f=f_{sur}$ used for the numerical flux $f_{sur}^*$ and the already mentioned volume
# flux $f_{vol}$ especially for this formulation.

# This formulation creates a more stable version of DGSEM, because it fulfils entropy stability.
# Moreover it allows the construction of entropy conserving discretizations without relying on
# exact integration. This is achieved when using a two-point entropy conserving flux function as
# volume flux in the volume flux differencing formulation.
# Then, the numerical surface flux can be used to control the dissipation of the discretization and to
# guarantee decreasing entropy, i.e. entropy stability.

# ## Implementation in Trixi
# Now, we have a look at the implementation of DGSEM with flux differencing with [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).
using OrdinaryDiffEq, Trixi

# We implement a simulation for the compressible Euler equations in 2D
# ```math
# \partial_t \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho e \end{pmatrix}
# + \partial_x \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1 \end{pmatrix}
# + \partial_y \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2 \end{pmatrix}
# = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}
# ```
# for an ideal gas with ratio of specific heats $\gamma=1.4$.
# Here, $\rho$ is the density, $v_1$, $v_2$ the velocities, $e$ the specific total energy and
# ```math
# p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right)
# ```
# the pressure.

gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

# As our initial condition we will use a weak blast wave from [Hennemann, Gassner (2020)](https://arxiv.org/abs/2008.12044).
# The primitive variables are defined by
# ```math
# \begin{pmatrix} \rho \\ v_1 \\ v_2 \\ p \end{pmatrix}
# = \begin{pmatrix} 1.0 \\ 0.0 \\ 0.0 \\ 1.0 \end{pmatrix} \text{if } \|x\|_2 > 0.5,\;
# \text{and } \begin{pmatrix} \rho \\ v_1 \\ v_2 \\ p \end{pmatrix}
# = \begin{pmatrix} 1.1691 \\ 0.1882 * \cos(\phi) \\ 0.1882 * \sin(\phi) \\ 1.245 \end{pmatrix} \text{else}
# ```
# with $\phi = \tan^{-1}(\frac{x_2}{x_1})$.

# This initial condition is implemented in Trixi under the name [`initial_condition_weak_blast_wave`](@ref).
initial_condition = initial_condition_weak_blast_wave

# In Trixi, flux differencing for the volume integral can be implemented with
# [`VolumeIntegralFluxDifferencing`](@ref) using symmetric two-point volume fluxes.
# First, we set up a simulation with the entropy conserving and kinetic energy preserving
# flux [`flux_ranocha`](@ref) by [Hendrik Ranocha (2018)](https://cuvillier.de/en/shop/publications/7743)
# as surface and volume flux.

# We will confirm the entropy conservation property numerically.

volume_flux = flux_ranocha # = f_vol
solver = DGSEM(polydeg=3, surface_flux=volume_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# Now, we implement Trixi's `mesh`, `semi` and `ode` in a simple framework. For more information please
# have a look at the documentation, the basic tutorial [introduction to DG methods](@ref scalar_linear_advection_1d)
# or some basic elixirs.
coordinates_min = (-2.0, -2.0)
coordinates_max = ( 2.0,  2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=5,
                n_cells_max=10_000,
                periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_condition_periodic)

## ODE solvers
tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan);

# To analyse the entropy conservation of the approximation, we will use the analysis calllback
# implemented in Trixi. It provides some information about the approximation including the entropy change.
analysis_callback = AnalysisCallback(semi, interval=100);

# We now run the simulation using `flux_ranocha` for both surface and volume flux.
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
            callback=analysis_callback, save_everystep=false);
# A look at the change in entropy $\sum \partial S/\partial U \cdot U_t$ in the analysis callback
# confirms that the flux is entropy conserving since the change is about machine precision.

# We can plot the approximated solution at the time `t=0.4`.
using Plots
plot(sol)

# Now, we can use for instance the dissipative flux [`flux_lax_friedrichs`](@ref) as surface flux
# to get an entropy stable method.
using OrdinaryDiffEq, Trixi

gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

initial_condition = initial_condition_weak_blast_wave

volume_flux = flux_ranocha # = f_vol
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-2.0, -2.0)
coordinates_max = ( 2.0,  2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=5,
                n_cells_max=10_000,
                periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_condition_periodic)

## ODE solvers
tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan);

analysis_callback = AnalysisCallback(semi, interval=100);

# We now run the simulation using the volume flux `flux_ranocha` and surface flux `flux_lax_friedrichs`.
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
            callback=analysis_callback, save_everystep=false);
# The change in entropy confirms the expected entropy stability.

using Plots
plot(sol)

# Of course, you can use more than these two fluxes in Trixi. Here, we will give a short list
# of possible fluxes for the compressible Euler equations.
# For the volume flux Trixi provides for example [`flux_ranocha`](@ref), [`flux_shima_etal`](@ref),
# [`flux_chandrashekar`](@ref), [`flux_kennedy_gruber`](@ref).
# As surface flux you can use all volume fluxes and additionally for instance [`flux_lax_friedrichs`](@ref),
# [`flux_hll`](@ref), [`flux_hllc`](@ref).


text:
using DelimitedFiles
using DifferentialEquations
using ProgressMeter
using PyCall
include("PSS_function_BIG.jl")
pygui(:qt)
using PyPlot
pygui(true)

Energy=BigFloat("0.125789")
H=(2*Energy)^2
max_hit=BigInt(100000)
t_end=BigFloat("1000")

Q=BigFloat(".0000000001")
P=BigFloat(".0000000001")

Y,Q_PSS, P_PSS=PSS_function_BIG(Q,P,Energy,t_end, max_hit)
plot(P_PSS,Q_PSS, color="k",".",markersize=1, markeredgewidth=.1)

ext: jl


text: using Ion; Ion.comonicon_install()


text:
using DiffEqBayes, BenchmarkTools

using OrdinaryDiffEq, RecursiveArrayTools, Distributions, ParameterizedFunctions, CmdStan, DynamicHMC
using Plots

gr(fmt=:png)

fitz = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

prob_ode_fitzhughnagumo = ODEProblem(fitz,[1.0,1.0],(0.0,10.0),[0.7,0.8,1/12.5,0.5])
sol = solve(prob_ode_fitzhughnagumo, Tsit5())

t = collect(range(1,stop=10,length=10))
sig = 0.20
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))

scatter(t, data[1,:])
scatter!(t, data[2,:])
plot!(sol)

priors = [Truncated(Normal(1.0,0.5),0,1.5),Truncated(Normal(1.0,0.5),0,1.5),Truncated(Normal(0.0,0.5),0.0,0.5),Truncated(Normal(0.5,0.5),0,1)]

@btime bayesian_result_stan = stan_inference(prob_ode_fitzhughnagumo,t,data,priors;num_samples = 10_000,printsummary=false)

@btime bayesian_result_turing = turing_inference(prob_ode_fitzhughnagumo,Tsit5(),t,data,priors;num_samples = 10_000)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using PyPlot, Random, Revise, Statistics, LinearAlgebra,
      Distributed, DelimitedFiles, HiQGA.transD_GP
## 1D functions
Random.seed!(200)
x = readdlm("func2.txt", ',', Float64, '\n')[:,1]
y = readdlm("func2.txt", ',', Float64, '\n')[:,2]
σ, fractrain = 0.275, 0.5
ntrain = round(Int, fractrain*length(y))
ynoisy = similar(y) .+ NaN
linidx = randperm(length(y))[1:ntrain]
ynoisy[linidx] = y[linidx] + σ*randn(ntrain)
line = transD_GP.Line(ynoisy;useML=false, σ=σ)
figure(figsize=(4,3))
plot(x[:], y)
plot(x[:], ynoisy, ".m", alpha=0.5)
xlabel("x")
ylabel("f(x)")
plt.tight_layout()
savefig("jump1D.png", dpi=300)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
#=move_user_drop.jl - > drops the piece at the coordinates
given. accepts 4 command line argument,<filename> => database <piece> => piece
<xtarget> => xTarget <ytarget> => yTarget
=#
#include("square.jl")
include("dParse.jl")

module move_user_drop
#using ST
using SQLite
function drop(database,pieceToDrop,xTarget,yTarget)
  #database = ARGS[1] #/path/to/database/file {string}
  #pieceToDrop = chomp(ARGS[2])[1] #piece {char}
  #xTarget = parse(Int,chomp(ARGS[3])) #x coordinate to place piece {Int}
  #yTarget = parse(Int,chomp(ARGS[4])) #y coordinate to place piece {Int}
  board = ST.loadBoard()
  db = SQLite.DB(database) #Opens the database gamefile

    #= ---- Determines the move_number ---- =#

    res = SQLite.query(db,"""SELECT MAX("move_number") FROM moves;""") #Finds the last played move (maximum move_number)
    lastMoveIDNullable = res[1,1] #SQL query with max move_number (POSSIBLY "NULL" if no moves have been made)

    if (!isnull(lastMoveIDNullable)) #Checks that lastMoveID was not NULL
      lastMoveID = get(res[1,1]) #Parses the last move move_number as an Int

    else #lastMoveID is NULL
      lastMoveID = 0 #move_number "0" is unsused. This implies no moves have been made

    end
    move_number = lastMoveID+1 #Current move number

  #=-----UPDATE DATABASE W/MOVE-----=#
  #Option will be dropped pieces abv
    SQLite.query(db,"INSERT INTO moves (move_number, move_type, targetx, targety, option)
    VALUES ($(move_number),'drop',$xTarget, $yTarget, '$pieceToDrop');")
    board[xTarget,yTarget].piece = pieceToDrop
    if iseven(move_number)
      board[xTarget,yTarget].team = 'w'
    else
      board[xTarget,yTarget].team = 'b'
    end

  ST.saveBoard(board)
end
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using Plots
using jInv.Mesh
using Flux
using Printf
using Statistics
using LinearAlgebra
using JLD
using MFGnet
using MAT

include("../ROLNWF2019/viewers.jl")
include("../ROLNWF2019/runOMThelpers.jl")

d   = 2

mu0 = [0.0;  3.0]
mu1 = [0.0; -3.0]
sig0 = [0.3;0.3];
sig1 = [0.3;0.3];
Qheight = 50;
sigQ = [1.0;0.5];
domain = [-3.0 3 -4.5 4.5];

rho0 = Gaussian(d,sig0,mu0)
rho1 = Gaussian(d,sig1,mu1)
Q    = Gaussian(2,sigQ,zeros(R,2),R(Qheight))

M      = getRegularMesh(domain,[512, 512])
X0     = Matrix(getCellCenteredGrid(M)')

# setup validation
rho0x = rho0(X0)
rho1x = rho1(X0)
Qx    = Q(X0)

matwrite("ObstacleProblem2D.mat", Dict(
	"Qx" => Qx,
    "rho0x" => rho0x,
    "rho1x" => rho1x,
    "domain" => M.domain,
    "n" => M.n))

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using Documenter, NCEI

makedocs(modules = [NCEI],
         sitename = "NCEI.jl",
         pages = ["index.md",
                  "GettingStarted.md",
                  "Walkthrough.md"
                  ])
deploydocs(repo = "github.com/Nosferican/NCEI.jl.git")


text:
using Documenter
using QuantumOpticsBase

pages = [
        "index.md",
        "api.md"
    ]

makedocs(
    sitename = "QuantumOpticsBase.jl",
    modules = [QuantumOpticsBase],
    pages = pages,
    checkdocs=:exports
    )

deploydocs(
    repo = "github.com/qojulia/QuantumOpticsBase.jl.git",
    )


text:
using Documenter, JSON3, StructTypes

makedocs(;
    modules=[JSON3],
    format=Documenter.HTML(edit_link="main"),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/quinnj/JSON3.jl/blob/{commit}{path}#L{line}",
    sitename="JSON3.jl",
    authors="Jacob Quinn",
    strict=true,
)

deploydocs(;
    repo="github.com/quinnj/JSON3.jl",
    devbranch = "main",
)


text:
function report_progress(callback::Function, name::AbstractString)
    global btask_report
    put!(btask_report, (name, callback))
    return nothing
end

function progress_bar(progress::Float64, bar_width::Int)
    progress = max(0,min(1,progress))
    num_hashes = floor(Int,progress*bar_width)
    return "#"^num_hashes * " "^(bar_width - num_hashes)
end

type Cycle
    t::Union{Void,Timer}
    name::AbstractString
    progress::Float64
    speed::Float64

    function Cycle(name::AbstractString)
        # Generate random speed
        speed = rand()*.01

        # Create Cycle object (have to create dummy Timer to avoid Nullable)
        cycle = new(nothing, name, 0, speed)

        # Create timer to increment progress and wrap it around
        cycle.t = Timer(t -> begin
            cycle.progress += cycle.speed
            while cycle.progress >= 1.0
                cycle.progress -= 1.0
            end
        end, 0, .01)

        # Don't forget to clean ourselves up!
        finalizer(cycle, x -> begin
            if x.t != nothing
                close(x.t)
            end
            x.progress = 1.0
        end)
        return cycle
    end
end

function stop(c::Cycle)
    close(c.t)
    c.progress = 1.0
end

# Return the progress line for a particular cycle
function cycle_progress_line(c::Cycle, width)
    if c.progress >= 1.0
        return nothing
    end
    bar_width = max(0, width - strwidth(c.name)) - 32
    return @sprintf("%s: [%s] (%02.1f%%)", c.name, progress_bar(c.progress, bar_width), c.progress*100)
end

cycles = Dict{AbstractString,Cycle}()
function start_cycling(name::AbstractString)
    global cycles

    # Cycle object creates a Timer object that runs itself ragged for us
    cycle = Cycle(name)
    cycles[name] = cycle

    # Tell btask how to ask for our status line of width `width`
    report_progress(name) do width
        return cycle_progress_line(cycle, width)
    end
    return nothing
end

function stop_cycling(name::AbstractString)
    global cycles
    if haskey(cycles, name)
        stop(cycles[name])
        delete!(cycles, name)

        # Kill this name
        report_progress(name) do width
            return nothing
        end
    end
    return nothing
end

function stop_cycling()
    global cycles
    cycles = []
end

function cursor_position(state)
    const term = TTYTerminal("",STDIN,STDOUT,STDERR)
    const term_width = LineEdit.width(term)
    buffstr = takebuf_string(copy(LineEdit.buffer(state)))

    # Absolute offset of cursor
    pos = position(LineEdit.buffer(state))

    # Calculate which line that is on
    splitbuff = split(buffstr[1:pos], "\n")
    line = sum(Int[div(max(0,strwidth(line)-1),term_width)+1 for line in splitbuff]) - 1

    # The length of the last line is the offset into that line of the cursor
    offset = strwidth(splitbuff[end])

    return line, offset
end

function prompt_lines(state)
    const term = TTYTerminal("",STDIN,STDOUT,STDERR)
    const term_width = LineEdit.width(term)
    buffstr = takebuf_string(copy(LineEdit.buffer(state)))

    # Calculate number of lines in the full buffer
    return sum(Int[div(max(0,strwidth(line)-1),term_width)+1 for line in split(buffstr, "\n")])
end


text:
@inbounds function residual(ts, y, p, u0, Tign, idt)
    ng = length(ts)
    pred = reshape(y, :, ng)
    F = similar(pred) * p[1]
    _ts = ts ./ idt
    i = 1
    @view(F[1:end - 1, i]) .= @views((pred[1:end - 1, i] .- u0))
    @view(F[end, i]) .= @views((pred[end, i + 1] - pred[end, i]))
    du = similar(pred[1:end - 1, 1]) * p[1]
    for i = 2:ng - 1
        @view(F[1:end - 1, i]) .=
            @views((pred[1:end - 1, i] .- pred[1:end - 1, i - 1]) ./ (_ts[i] - _ts[i - 1])) .-
            dudt!(du, @view(pred[1:end - 1, i]), p, 0.0) .* pred[end, i]
        @view(F[end, i]) .=
            @views(pred[end, i + 1] - pred[end, i])
    end
    i = ng
    @view(F[1:end - 1, i]) .=
        @views((pred[1:end - 1, i] .- pred[1:end - 1, i - 1]) ./ (_ts[i] - _ts[i - 1])) .-
        dudt!(du, @view(pred[1:end - 1, i]), p, 0.0) .* pred[end, i]
    @view(F[end, i]) .= @views((pred[end - 1, i] - Tign))
    return vcat(F...)
end

function downsampling(ts, pred; dT=0.1, n_max=100, verbose=false)
    ind_sample = [1]
    Ts = pred[end, :]
    _T = Ts[1]
    _t = ts[1]
    for i = 2:length(ts) - 1
        if (abs(Ts[i] - _T) > dT) & (ts[i] - _t > ts[end] / n_max)
            _T = Ts[i]
            _t = ts[i]
            push!(ind_sample, i)
        end
    end
    push!(ind_sample, length(ts))

    if verbose
        println("\n original sample size $(length(ts))
                 downsampled size $(length(ind_sample)) \n")
    end

    ts = ts[ind_sample]
    pred = pred[:, ind_sample]

    return ts, pred
end

@inbounds function sensBVP(ts, pred, p)
    local ng = length(ts)
    Fp_ = zeros(ng * nu, np)
    Fy_ = BandedMatrix(Zeros(ng * nu, ng * nu), (nu, nu))
    return sensBVP!(Fy_, Fp_, ts, pred, p)
end

@inbounds function sensBVP!(Fy, Fp, ts, pred, p)
    local idt = ts[end]
    local Tign = pred[end, end]
    local ng = length(ts)
    Fy .*= 0.0
    Fp .*= 0.0
    i = 1
    i_F = 1 + (i - 1) * nu:i * nu - 1
    @view(Fy[i_F, i_F])[ind_diag] .= ones_nu
    Fy[i * nu, i * nu] = -1.0
    Fy[i * nu, (i + 1) * nu] = 1.0
    du = similar(@view(pred[:, i]))

    dts = @views(ts[2:end] .- ts[1:end - 1]) ./ idt
    for i = 2:ng
        u = @view(pred[:, i])
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fy[i_F, i_F]) .= jacobian((du, x) -> dudt!(du, x, p, 0.0),
                                         du, u)::Array{Float64,2} .* (-idt)
        @view(Fy[i_F, i * nu]) .= - dudt!(du, u, p, 0.0)
        @view(Fy[i_F, i_F])[ind_diag] .+= ones_nu ./ (dts[i - 1])
        @view(Fy[i_F, i_F .- nu])[ind_diag] .+= ones_nu ./ (-dts[i - 1])
        if i < ng
            Fy[i * nu, i * nu] = -1.0
            Fy[i * nu, (i + 1) * nu] = 1.0
        else
            Fy[i * nu, i * nu - 1] = 1.0
        end
        @view(Fp[i_F, :]) .= jacobian((du, x) -> dudt!(du, u, x, 0.0),
                                       du, p)::Array{Float64,2} .* (-idt)
    end
    dydp = - @views(Fy[1:ng * nu, 1:ng * nu] \ Fp[1:ng * nu, :])
    return @view(dydp[end, :]) ./ idt
end

function sensBVP_mthread(ts, pred, p)
    idt = ts[end]
    Tign = pred[end, end]
    ng = length(ts)
    Fy = BandedMatrix(Zeros(ng * nu, ng * nu), (nu, nu));
    Fp = zeros(ng * nu, np)
    i = 1
    i_F = 1 + (i - 1) * nu:i * nu - 1
    @view(Fy[i_F, i_F])[ind_diag] .= ones_nu
    Fy[i * nu, i * nu] = -1.0
    Fy[i * nu, (i + 1) * nu] = 1.0

    dts = @views(ts[2:end] .- ts[1:end - 1]) ./ idt
    @threads for i = 2:ng
        u = @view(pred[:, i])
        du = similar(u)
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fp[i_F, :]) .= jacobian((du, x) -> dudt!(du, u, x, 0.0),
                                    du, p)::Array{Float64,2} .* (-idt)
    end
    @threads for i = 2:ng
        u = @view(pred[:, i])
        du = similar(u)
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fy[i_F, i_F]) .= jacobian((du, x) -> dudt!(du, x, p, 0.0),
                                    du, u)::Array{Float64,2} .* (-idt)
    end
    for i = 2:ng
        u = @view(pred[:, i])
        du = similar(u)
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fy[i_F, i * nu]) .= - dudt!(du, u, p, 0.0)
        @view(Fy[i_F, i_F])[ind_diag] .+= ones_nu ./ (dts[i - 1])
        @view(Fy[i_F, i_F .- nu])[ind_diag] .+= ones_nu ./ (-dts[i - 1])
        if i < ng
            Fy[i * nu, i * nu] = -1.0
            Fy[i * nu, (i + 1) * nu] = 1.0
        else
            Fy[i * nu, i * nu - 1] = 1.0
        end
    end
    dydp = - Fy \ sparse(Fp)
    return @view(dydp[end, :]) ./ idt
end

# sensitivity calculated by BFSA method
# d\tau/dp = d\tau/dT * dT/dp
sensealg = ForwardDiffSensitivity()
function sensBFSA(phi, P, T0, p; dT=200, dTabort=600, tfinal=1.0)
    ts, pred = get_Tcurve(phi, P, T0, p;
                          dT=dT, dTabort=dTabort, tfinal=tfinal);
    idt = interpx(ts, pred[end,:], pred[end,1] + dT);
    dTdidt = (pred[end,end] - pred[end,end - 1]) / (ts[end] - ts[end - 1]);

    prob = make_prob(phi, P, T0, p; tfinal=idt)
    function fsol_T(p, idt)
        sol = solve(prob, Trapezoid(), p=p, saveat=[0,idt],
                    reltol=1e-6, abstol=1e-9, sensealg=sensealg);
        return sol[end,end];
    end

    dTdp = ForwardDiff.gradient(x -> fsol_T(x, idt), p);

    return dTdp ./ dTdidt ./ idt
end

# sensitivity calculated by Brute-force method
function sensBF_mthread(phi, P, T0, p; dT=200, dTabort=600, pdiff=5e-3)
    idt = get_idt(phi, P, T0, p; dT=dT, dTabort=dTabort)
    np = length(p);
    sens = zeros(np);
    @threads for i = 1:np
        pp = deepcopy(p);
        pp[i] += pdiff;
        pidt = get_idt(phi, P, T0, pp; dT=dT, dTabort=dTabort)
        sens[i] = (pidt - idt) / pdiff
    end
    return sens ./ idt;
end

max_stars_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_issues_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_forks_repo_name: DENG-MIT/ArrheniusActiveSubspace


text:
"""
Abstract type for sampling with Base.rand(ProbabilisticGraphicalModel, Sampler, nsamples)
                                Base.rand!(Assignment, ProbabilisticGraphicalModel, Sampler)
                                Base.rand(ProbabilisticGraphicalModel, Sampler)
"""
abstract type Sampler end

"""
Overwrites Assignment with a sample from the PGM using the given Sampler
"""
Base.rand!(a::Assignment, pgm::ProbabilisticGraphicalModel, sampler::Sampler) = error("rand! not implemented for $(typeof(Sampler))")

"""
Returns a new Assignment sampled from the PGM using the provided sampler
"""
Base.rand(pgm::ProbabilisticGraphicalModel, sampler::Sampler) = rand!(Assignment(), bn, sampler)

"""
Generates a DataFrame containing a dataset of variable assignments.
Always return a DataFrame with `nsamples` rows.
"""
function Base.rand(pgm::ProbabilisticGraphicalModel, sampler::Sampler, nsamples::Integer)

    a = rand(pgm, sampler)
    df = DataFrame()
    nodenames = names(pgm)
    for nodename in nodenames
        df[nodename] = Array{typeof(a[nodename])}(nsamples)
    end

    for i in 1:nsamples
        rand!(a, pgm, sampler)
        for name in nodenames
            df[i, name] = a[name]
        end
    end

    df
end


text:
# what is the difference between schur direct and schur????

mutable struct Schur_KKT_solver_direct <: abstract_schur_solver
    # abstract_KKT_system_solver
    ls_solver::abstract_linear_system_solver
    factor_it::Class_iterate
    delta_x_vec::Array{Float64,1}
    delta_s_vec::Array{Float64,1}
    rhs::System_rhs
    dir::Class_point
    kkt_err_norm::Class_kkt_error
    rhs_norm::Float64
    pars::Class_parameters

    ready::Symbol #

    # Schur_KKT_solver only
    true_diag::Array{Float64,1}
    M::SparseMatrixCSC{Float64,Int64}
    K::SparseMatrixCSC{Float64,Int64}
    current_it::Class_iterate
    reduct_factors::Class_reduction_factors

    function Schur_KKT_solver_direct()
      this = new()
      this.ready = :not_ready

      return this
    end
end

function compute_direction_implementation!(kkt_solver::Schur_KKT_solver_direct, timer::class_advanced_timer)
    start_advanced_timer(timer, "SCHUR")

    it = kkt_solver.current_it
    reduct = kkt_solver.reduct_factors
    pars = kkt_solver.pars

    #∇a = get_jac(it);
    y = get_y(it);
    s = get_s(it);

    rhs = kkt_solver.rhs

    #r1 + ∇a_org' * (( r3 + (r2 .* y_org) ) ./ s_org)
    start_advanced_timer(timer, "SCHUR/rhs");
    symmetric_primal_rhs = rhs.primal_r + rhs.comp_r ./ y
    Σ_vec = ( y ./ s )
    y_ = rhs.primal_r .* Σ_vec + rhs.comp_r ./ s
    schur_rhs = rhs.dual_r + eval_jac_T_prod(it, y_)
    pause_advanced_timer(timer, "SCHUR/rhs");

    dir = kkt_solver.dir
    dir.x = solver_schur_rhs(schur_rhs, kkt_solver, timer)
    dir.y = -(eval_jac_prod(it, dir.x) - symmetric_primal_rhs) .* Σ_vec
    dir.s = ( rhs.comp_r - dir.y .* s ) ./ y
    #dir.s = get_jac(it) * dir.x - rhs.primal_r
    #dir.y = ( rhs.comp_r - dir.s .* get_y(it) ) ./ get_s(it) # (mu_target - s_cur .* y_cur)

    check_for_nan(dir)

    start_advanced_timer(timer, "SCHUR/kkt_err");
    update_kkt_error!(kkt_solver, Inf, timer)
    pause_advanced_timer(timer, "SCHUR/kkt_err");
    pause_advanced_timer(timer, "SCHUR")
end


text:
mutable struct LinearBucket <: AbstractRRModel
    K::Float64
    S::Float64

    function LinearBucket(K::Float64, S::Float64)
        0.0 <= K <= 1.0 || throw(DomainError(K, "K parameter should have a value between 0 and 1"))
        S >= 0.0 || throw(DomainError(S, "S initial storage should have a non-negative value"))

        new(K, S)
    end
end

function LinearBucket()
    LinearBucket(0.1, 0.0)
end

function LinearBucket(model::LinearBucket)
    LinearBucket(model.K, model.S)
end

function LinearBucket(d::Dict)
    LinearBucket(d[:K], d[:S])
end

function ==(a::LinearBucket, b::LinearBucket)
    return a.K == b.K && a.S == b.S
end

function show(io::IO, model::LinearBucket)
    println("Linear Bucket Model:")
    println(io, "K = $(model.K)")
    println(io, "S = $(model.S)")
    return
end

function Base.setproperty!(obj::LinearBucket, field::Symbol, value)
    if field == :K
        0.0 <= value <= 1.0 || throw(DomainError(value, "K parameter should have a value between 0 and 1"))
    elseif field == :S
        value >= 0.0 || throw(DomainError(value, "S initial storage should have a non-negative value"))
    end
    setfield!(obj, field, value)
end

function stash(model::LinearBucket)::Dict
    return Dict(:K=>model.K, :S=>model.S)
end

function simulate(model::LinearBucket, P::Array{Float64,1}, PET::Array{Float64,1})
    if size(P) != size(PET)
        throw(DimensionMismath("P and PET have different sizes"))
    end

    Q = zeros(Float64, size(P))
    AET = zeros(Float64, size(P))
    for i in 1:length(P)
        Q[i], AET[i] = simulate(model, P[i], PET[i])
    end
    return Q, AET
end

function simulate(model::LinearBucket, P::Float64, PET::Float64)
    Q = 0.0
    AET = 0.0
    for i in 1:length(P)
        model.S += P
        AET = model.S >= PET ? PET : model.S
        model.S -= AET
        Q = model.K * model.S
        model.S -= Q
    end
    return Q, AET
end


text:
using Metatheory, Metatheory.EGraphs
using Test

dbgproof(n::ENode) = println("$(n.proof_src) ⩜ $(n.proof_trg)")

function prove(g::EGraph, t::Vector{<:Rule}, exprs...;
    params=SaturationParams())
    # @info "Checking equality for " exprs
    n = length(exprs)
    if n == 1; return true end
    # rebuild!(G)

    ids = Vector{EClassId}(undef, n)
    nodes = Vector{ENode}(undef, n)
    for i ∈ 1:n
        ec, node = addexpr!(g, exprs[i])
        ids[i] = ec.id
        nodes[i] = node
    end

    goal = EqualityGoal(collect(exprs), ids)

    # params.goal = goal
    report = saturate!(g, t, params; mod=mod)

    # display(g.classes); println()
    if !(report.reason isa EGraphs.Saturated) && !reached(g, goal)
        return missing # failed to prove
    end

    # @show reached(g, goal)

    for (id, ec) in g.classes
        for n in ec
            # println(id => n)
            # dbgproof(n)
        end
    end

    for i in 1:n
        node = nodes[i]
        if haskey(g.memo, node)
            # TODO really override the proof step here?
            eclass = geteclass(g, g.memo[node])
            for nn in eclass
                if node == nn
                    # dbgproof(node)
                    # dbgproof(nn)
                    # println("$node == $nn")
                    nodes[i] = nn
                end
            end
        end
    end
    # println("========================================")
    # for i in 1:n
    #     nn = nodes[i]
    #     dbgproof(nn)
    # end
    # println("========================================")
    @show reached(g, goal)
    proof_bfs(g, nodes[1], nodes[2])
end

mutable struct ProofNode
    state::ENode
    why::Union{Nothing,Rule}
    when::Int
    cost::Number
    parent::Union{Nothing, ProofNode}
end

struct Proof
    g::EGraph
    head::ProofNode
end

"""
A closured cost function that considers the age of the enode and the ast size.
"""
function oldestatage(age::Int)
    return (n::ENode, g::EGraph, an::Type{<:AbstractAnalysis}) -> begin
        # cost = 0
        # println("current age is $age")
        # println("enode $n age is $(n.age)")
        # cost = n.age - age
        cost = n.age - age
        # println("cost is $cost")
        return cost
    end
end

function oldest_node_extract(g::EGraph, n::ENode, age::Int)
    costfun = oldestatage(age)
    ex = EGraphs.extractnode(g, n, ExtractionAnalysis{costfun})
    # println("extracted $ex aged $(n.age) at age $age")
    return ex
end

function Base.show(io::IO, mime::MIME"text/plain", proof::Proof)
    lines = []
    curr = proof.head
    while !isnothing(curr.parent)
        ex = oldest_node_extract(proof.g, curr.state, curr.when)
        pushfirst!(lines, "$ex")
        pushfirst!(lines, "from $(repr("text/plain", curr.why))")
        curr = curr.parent
    end
    ex = oldest_node_extract(proof.g, curr.state, curr.when)
    pushfirst!(lines, "given $ex")

    for line in lines
        println(io, line)
    end
end

# TODO go through each path in the proof. do a BFS?
using DataStructures
function proof_bfs(g::EGraph, src::ENode, trg::ENode)
    root = ProofNode(src, nothing, src.age, 0, nothing)
    if src == trg
        return Proof(g, root)
    end
    frontier = ProofNode[]
    explored = Set{ENode}()
    push!(frontier, root)
    while !isempty(frontier)
        node = popfirst!(frontier)
        # println("exploring $node")
        push!(explored, node.state)
        # todo take rules in account
        for (rule, child_enode, age) in unique(node.state.proof_trg) #∪ node.state.proof_src
            child = ProofNode(child_enode, rule, age, node.cost+1, node)
            if child_enode ∉ explored && child ∉ frontier
                # goal test
                if child_enode == trg
                    return Proof(g, child)
                end
                push!(frontier, child)
            end
        end
    end
    error("proof not found!")
end

t = @theory begin
    a * b == b * a
    a * 2 == a + a
end

# g = EGraph()
# addexpr!(g, :(2 * a))
# addexpr!(g, :(a + a))
# saturate!(g, t)
# proof = prove(g, t, :(2 * a), :(a + a))

# Base.print(proof)

prove(EGraph(), t, :(2 * x), :(x + x))

include("../logic/prop_logic_theory.jl")

@metatheory_init

ex = :((x => (y => z)) => ((x => y) => (x => z)))
proof = prove(EGraph(), t, ex, true)

proof = prove(EGraph(), t, :((x => (x ∨ x))), :((¬(x) ∧ y) => y))

# TODO introduce a mechanism of enode age and egraph age to be able to extract
# precisely which enode was there at that moment
# to be printed, if rule was applied at time `t`, then extract the enodes that was
# applied earliest but before (or after???) `t`
ex = :(((x ∨ y) ∨ ¬(z ∧ a)) ∨ a)
proof = prove(EGraph(), t, ex, true)

Metatheory.options.verbose = false

# chat with oliver
# introduce a new type of enode at the proof generation stage
# set of unknown variables

# keep proof at the expr level
# do not store src and trg in enode
# store holes
#
# keep another unionfind for the same holes


text:
using Base
using Test
using LinearAlgebra
using Surrogates

@testset "Radial" begin

    @testset "1D" begin

        #WITHOUT ADD_POINT
        x = [1.0,2.0,3.0]
        y = [4.0,5.0,6.0]
        a = 0
        b = 4
        my_rad = RadialBasis(x,y,a,b,z->norm(z),1)
        est = my_rad(3.0)
        @test est ≈ 7.875

        #WITH ADD_POINT, adding singleton
        x = [1.0,2.0,3.0]
        y = [4.0,5.0,6.0]
        a = 0
        b = 4
        my_rad = RadialBasis(x,y,a,b,z->norm(z),1)
        add_point!(my_rad,4.0,10.0)
        est = my_rad(3.0)
        @test est ≈ 6.499999999999991

        #WITH ADD_POINT, adding more

        x = [1.0,2.0,3.0]
        y = [4.0,5.0,6.0]
        a = 0
        b = 4
        my_rad = RadialBasis(x,y,a,b,z->norm(z),1)
        est_rad = my_rad(3.0)
        add_point!(my_rad,[3.2,3.3,3.4],[8.0,9.0,10.0])
        est_rad2 = my_rad(3.0)
        @test est_rad2 ≈ 6.49416593818781

    end
    @testset "ND" begin
        #WITHOUT ADD_POINT

        x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
        y = [4.0,5.0,6.0]
        bounds = [[0.0,4.0],[3.0,7.0],[6.0,10.0]]
        my_rad = RadialBasis(x,y,bounds,z->norm(z),1)
        est = my_rad((1.0,2.0,3.0))
        @test est ≈ 4.0

        #WITH ADD_POINT, adding singleton
        x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
        y = [4.0,5.0,6.0]
        bounds = [[0.0,4.0],[3.0,7.0],[6.0,10.0]]
        my_rad = RadialBasis(x,y,bounds,z->norm(z),1)
        add_point!(my_rad,(9.0,10.0,11.0),10.0)
        est = my_rad((1.0,2.0,3.0))
        @test est ≈ 4.0

        #WITH ADD_POINT, adding more

        x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
        y = [4.0,5.0,6.0]
        bounds = [[0.0,4.0],[3.0,7.0],[6.0,10.0]]
        my_rad = RadialBasis(x,y,bounds,z->norm(z),1)
        add_point!(my_rad,[(9.0,10.0,11.0),(12.0,13.0,14.0)],[10.0,11.0])
        est = my_rad((1.0,2.0,3.0))
        @test est ≈ 4.0

    end
end


text:
export BlockIndex,
	outerindex,
	innerindex

############
# Overview #
############

# We collect some methods and definitions having to do with various kinds of
# indexing. We make some assumptions on all index sets:
# - an index set is mathematically a set, i.e., it has no duplicates
# - in addition to being a set, the elements can be ordered and the ordering
#   is defined by the iterator of the set.
#
# The functionality in this file includes:
# - bounds checking: see `checkbounds` function
# - conversion between indices of different types
# - efficient iterators

"""
The `type IndexList` implements a map from linear indices to another family
of indices, and vice-versa.

It is implemented as an abstract vector. Hence, the functionality of vectors
is available. This also means that the dimension of an index list is always a
vector, whose length equals that of the dictionary it corresponds to. Still,
in several cases one may index the list with other kinds of indices.

Concrete subtypes should implement:
- `getindex(l::MyIndexList, idx::Int)` -> this is the map from linear indices
  to native indices
- `getindex(l::MyIndexList, idxn::MyIndex)` -> this is the inverse map
- `size(l::MyIndexList)`
"""
abstract type IndexList{T} <: AbstractVector{T}
end

# The concrete subtype should implement size (and perhaps length)

# Assume linear indexing, override if this is not appropriate
Base.IndexStyle(list::IndexList) = Base.IndexLinear()

##################
# Native indices
##################

# Dictionaries can be indexed in various ways. We assume that the semantics
# of the index is determined by its type and, moreover, that linear indices
# are always Int's. This means that no other index can have type Int.
const LinearIndex = Int

"""
An `AbstractIntegerIndex` represents an integer that is being used as an index
in the BasisFunctions package.

The type implements basic functionality of integers such that one can do
computations with the indices or construct ranges.
"""
abstract type AbstractIntegerIndex <: Integer end

# Convert to or back from other integers
convert(I::Type{<:AbstractIntegerIndex}, value::Int) = I(value)
convert(::Type{T}, idx::AbstractIntegerIndex) where {T<:Integer} = convert(T, value(idx))

# Some functionality of numbers that one might encounter in other conversions:
# - These lines support conversion to BigFloat
Base.div(idx::AbstractIntegerIndex, x::Integer, r::RoundingMode) = div(value(idx), x, r)
Base.div(x::Integer, idx::AbstractIntegerIndex, r::RoundingMode) = div(x, value(idx), r)
Base.div(i::AbstractIntegerIndex, j::AbstractIntegerIndex, r::RoundingMode) = div(value(i), value(j), r)
# - These lines support conversion to Double64 (among other things)
convert(::Type{Float32}, idx::AbstractIntegerIndex) = convert(Float32, value(idx))
convert(::Type{Float64}, idx::AbstractIntegerIndex) = convert(Float64, value(idx))
# Float32(idx::AbstractIntegerIndex) = Float32(value(idx))
Float64(idx::AbstractIntegerIndex) = Float64(value(idx))

# convert(::Type{T}, idx::AbstractIntegerIndex) where {T} = convert(T, value(idx))
# Resolve an ambiguity with mpfr code...
# convert(::Type{BigFloat}, idx::AbstractIntegerIndex) = convert(BigFloat, value(idx))
# convert(::Type{Double64}, idx::AbstractIntegerIndex) = convert(BigFloat, value(idx))

# With this line we inherit binary operations involving integers and native indices
Base.promote_rule(I::Type{<:AbstractIntegerIndex}, ::Type{Int}) = I
# For floating points, we choose to convert to the numeric value
Base.promote_rule(I::Type{<:AbstractIntegerIndex}, F::Type{<:AbstractFloat}) = F

for op in (:+, :-, :*)
    @eval $op(a::I, b::I) where {I<:AbstractIntegerIndex} = I($op(value(a),value(b)))
end

for op in (:<, :<=, :>, :>=)
    @eval $op(a::I, b::I) where {I<:AbstractIntegerIndex} = $op(value(a),value(b))
end

(-)(a::AbstractIntegerIndex) = typeof(a)(-value(a))

# Convenience: make a vector indexable using native indices, if possible
# It is possible whenever the size and element type of the vector completely
# determine the index map from native to linear indices
getindex(v::Array, idxn::AbstractIntegerIndex) =
    getindex(v, linear_index(idxn, size(v), eltype(v)))

setindex!(v::Array, val, idxn::AbstractIntegerIndex) =
    setindex!(v, val, linear_index(idxn, size(v), eltype(v)))

"""
A native index is distinguishable from a linear index by its type, but otherwise
it simply wraps an integer and it acts like an integer.
"""
struct NativeIndex{S} <: AbstractIntegerIndex
    value   ::  Int
end

NativeIndex{S}(a::NativeIndex{S}) where {S} = a

value(idx::NativeIndex) = idx.value

"""
The native index list is an identity map from integers to themselves, wrapped
in a `NativeIndex{S}` for some symbol `S`.
"""
struct NativeIndexList{S} <: IndexList{NativeIndex{S}}
    n		::	Int
end

size(list::NativeIndexList) = (list.n,)

getindex(list::NativeIndexList{S}, idx::Int) where {S} = NativeIndex{S}(idx)
getindex(list::NativeIndexList{S}, idxn::NativeIndex{S}) where {S} = value(idxn)

# In this case, we can compute the linear index without reference to
# any dictionary or list:
linear_index(idx::NativeIndex) = value(idx)
# The line below may be called with the size and eltype of an array when indexing
# into an array using a native index.
linear_index(idx::NativeIndex, size, T) = linear_index(idx)

"""
If a dictionary does not have a native index, we use the `DefaultNativeIndex`
type to distinguish between linear and native indices. However, the default
native index simply wraps the integer linear index.
"""
const DefaultNativeIndex = NativeIndex{:default}
const DefaultIndexList = NativeIndexList{:default}

"""
A `ShiftedIndex{S}` is a linear index shifted by `S`. Thus, its value starts at
`1-S` and ranges up to `n-S`. A typical case is `S=1`.
"""
abstract type AbstractShiftedIndex{S} <: AbstractIntegerIndex end

shift(idx::AbstractShiftedIndex{S}) where {S} = S
value(idx::AbstractShiftedIndex) = idx.value

# Here too, as with the DefaultIndex above, we can compute the linear index
# without reference to any dictionary or list
linear_index(idx::AbstractShiftedIndex) = value(idx) + shift(idx)
linear_index(idx::AbstractShiftedIndex, size, T) = linear_index(idx)

struct ShiftedIndex{S} <: AbstractShiftedIndex{S}
    value   ::  Int
end

ShiftedIndex{S}(a::ShiftedIndex{S}) where {S} = a

struct ShiftedIndexList{S,T<:AbstractShiftedIndex{S}} <: IndexList{T}
    n       ::  Int
end

ShiftedIndexList(n, T::Type{<:AbstractShiftedIndex{S}}) where {S} =
	ShiftedIndexList{S,T}(n)

size(list::ShiftedIndexList) = (list.n,)
shift(list::ShiftedIndexList{S,T}) where {S,T} = S

getindex(list::ShiftedIndexList{S,T}, idx::Int) where {S,T} = T(idx-S)
getindex(list::ShiftedIndexList{S,T}, idxn::AbstractShiftedIndex) where {S,T} = value(idxn)+S

########################
# Product indices
########################

# We reuse the CartesianIndex defined in Base
const ProductIndex{N} = CartesianIndex{N}

# Convenience functions
index(idx::ProductIndex) = value(idx)
# Return the index as a tuple
indextuple(idx::ProductIndex) = idx.I

"""
`ProductIndexList` is a list of product indices, suitable for use as the ordering
of a dictionary with tensor product structure.
"""
struct ProductIndexList{N} <: IndexList{ProductIndex{N}}
    size    ::  NTuple{N,Int}
end

# There is something dodgy about defining the index list to be an abstract vector,
# yet having its size be that of a tensor. What matters is that the length of
# the list is correct, and that the list can be indexed with integers.
# All code which assumes that length(list) == size(list)[1] will be fooled...
size(list::ProductIndexList) = list.size

# We use Cartesian indexing, because that is more efficient than linear indexing
# when iterating over product sets: we have to convert int's to cartesian indices
# anyway
Base.IndexStyle(list::ProductIndexList) = Base.IndexCartesian()

# We have to override eachindex, because the default eachindex in Base returns
# a linear index for a vector (because any IndexList is an AbstractVector), and
# the most efficient iteration over product dictionaries is using cartesian indices.
eachindex(list::ProductIndexList) = CartesianIndices(axes(list))

# We convert between integers and product indices using ind2sub and sub2ind
product_native_index(size, idx::Int) = ProductIndex(CartesianIndices(size)[idx])
product_linear_index(size, idxn::ProductIndex) = LinearIndices(size)[indextuple(idxn)...]

# We also know how to convert a tuple
product_native_index(size::NTuple{N,Int}, idx::NTuple{N,Int}) where {N} = ProductIndex(idx)
product_native_index(size::NTuple{N,Int}, idx::ProductIndex{N}) where {N} = idx

getindex(list::ProductIndexList, idx::Int) =
    product_native_index(size(list), idx)
getindex(list::ProductIndexList, idxn::ProductIndex) =
    product_linear_index(size(list), idxn)

# Convert a tuple of int's or a list of int's to a linear index
getindex(list::ProductIndexList{N}, idx::NTuple{N,Int}) where {N} =
    getindex(list, ProductIndex(idx))
getindex(list::ProductIndexList{N}, idx1::Int, idx2::Int, idx::Int...) where {N} =
    getindex(list, ProductIndex{N}(idx1, idx2, idx...))

## We implement a promotion system for two indices.
# The goal is to convert both indices to the Cartesian{N} type, if possible.
# We leave the indices unchanged if it is not possible.

# This is an exhaustive list of types we know how to convert
ProductIndices = Union{NTuple,ProductIndex}

# - Both indices are fine:
promote_product_indices(size::NTuple{N,Int}, idx1::ProductIndex{N}, idx2::ProductIndex{N}) where N = (idx1,idx2)
# - One or both of them are not fine, but they are both in PTI. The formulation
#   includes the case above, since they could both be ProductIndex{N}. By leaving
#   out the type of size, we make sure that the line above is more specific and
#   will be chosen by the compiler, thus avoiding an infinite loop.
promote_product_indices(size, idx1::ProductIndices, idx2::ProductIndices) =
    (product_native_index(size, idx1), product_native_index(size, idx2))
# - At least one of the indices does not have a suitable type, we don't know what to do
promote_product_indices(size, idx1, idx2) = (idx1,idx2)

########################
# Composite indices
########################

const MultilinearIndex = BlockIndex{1}

const MultilinearIndices = Union{MultilinearIndex,Tuple{Int,Any}}

# At the time of writing, BlockIndex has two fields I and α.
# For BlockIndex{1} they are a tuple with one entry. We provide
# outerindex and innerindex as generic access routines for the values they contain.
outerindex(idx::MultilinearIndex) = idx.I[1]
innerindex(idx::MultilinearIndex) = idx.α[1]

# For compatibility, we do the same for tuples with two elements.
outerindex(idx::Tuple{Int,Any}) = idx[1]
innerindex(idx::Tuple{Int,Any}) = idx[2]

"""
`MultiLinearIndexList` is a list of indices suitable for use as the ordering
of a dictionary with composite structure.
"""
struct MultilinearIndexList{O} <: IndexList{MultilinearIndex}
    offsets ::  O
end

const MLIndexList = MultilinearIndexList

size(list::MLIndexList) = (list.offsets[end],)

# Convert a linear index into a multilinear index using the offsets information
function offsets_multilinear_index(offsets, idx::Int)
    i = 0
    while idx > offsets[i+1]
        i += 1
    end
    BlockIndex(i, idx-offsets[i])
end

# and vice-versa
offsets_linear_index(offsets, idx::MultilinearIndices) = offsets[outerindex(idx)] + innerindex(idx)

getindex(list::MLIndexList, idx::Int) = offsets_multilinear_index(list.offsets, idx)
getindex(list::MLIndexList, idx::MultilinearIndices) = offsets_linear_index(list.offsets, idx)

"""
A `MultilinearIndexIterator` iterates over a sequence of linear indices. Each index
is a tuple, where the first entry refers to the linear index, and the second
entry is the linear index.

This resembles the `Flatten` iterator in `Base`, but the latter would yield only the
second entry. For example, the composite indices of `(1:5,1:3)` are:
```
(1,1), (1,2), (1,3), (1,4), (1,5), (2,1), (2,2), (2,3).
```
In contrast, the `Flatten` iterator in Base would yield:
```
1, 2, 3, 4, 5, 1, 2, 3
```
"""
struct MultilinearIndexIterator{L}
    lengths ::  L
end

length(it::MultilinearIndexIterator) = sum(it.lengths)

iterate(it::MultilinearIndexIterator) = BlockIndex(1,1), (1,1)

function iterate(it::MultilinearIndexIterator, state)
    i, j = state
    if j == it.lengths[i]
        next = (i+1,1)
    else
        next = (i,j+1)
    end
    if next[1] <= length(it.lengths)
        BlockIndex(next...), next
    end
end

Base.eltype(::Type{MultilinearIndexIterator}) = Tuple{Vararg{Int}}

last(iter::MultilinearIndexIterator) = BlockIndex(length(iter.lengths), iter.lengths[end])


text:
"""
    PS.jl

This is a Julia port of a MATLAB implementation of batch and incremental
Partition Separation (PS) Cluster Validity Index.

Authors:
MATLAB implementation: Leonardo Enzo Brito da Silva
Julia port: Sasha Petrenko <sap625@mst.edu>

REFERENCES
[1] Miin-Shen Yang and Kuo-Lung Wu, "A new validity index for fuzzy clustering,"
10th IEEE International Conference on Fuzzy Systems. (Cat. No.01CH37297), Melbourne,
Victoria, Australia, 2001, pp. 89-92, vol.1.
[2] E. Lughofer, "Extensions of vector quantization for incremental clustering," Pattern
Recognit., vol. 41, no. 3, pp. 995–1011, 2008.
"""

"""
    PS

The stateful information of the Partition Separation (PS) Cluster Validity Index.

# References
1. Miin-Shen Yang and Kuo-Lung Wu, "A new validity index for fuzzy clustering," 10th IEEE International Conference on Fuzzy Systems. (Cat. No.01CH37297), Melbourne, Victoria, Australia, 2001, pp. 89-92, vol.1.
2. E. Lughofer, "Extensions of vector quantization for incremental clustering," Pattern Recognit., vol. 41, no. 3, pp. 995–1011, 2008.
"""
mutable struct PS <: CVI
    label_map::LabelMap
    dim::Integer
    n_samples::Integer
    n::IntegerVector        # dim
    v::RealMatrix           # dim x n_clusters
    D::RealMatrix           # n_clusters x n_clusters
    v_bar::RealVector       # dim
    beta_t::Float
    PS_i::RealVector        # n_clusters
    n_clusters::Integer
    criterion_value::Float
end # PS <: CVI

"""
    PS()

Default constructor for the Partition Separation (PS) Cluster Validity Index.
"""
function PS()
    PS(
        LabelMap(),                     # label_map
        0,                              # dim
        0,                              # n_samples
        Array{Integer, 1}(undef, 0),    # n
        Array{Float, 2}(undef, 0, 0),   # v
        Array{Float, 2}(undef, 0, 0),   # D
        Array{Float, 1}(undef, 0),      # v_bar
        0.0,                            # beta_t
        Array{Float, 1}(undef, 0),      # PS_i
        0,                              # n_clusters
        0.0                             # criterion_value
    )
end # PS()

"""
    setup!(cvi::PS, sample::Vector{T}) where {T<:RealFP}
"""
function setup!(cvi::PS, sample::Vector{T}) where {T<:RealFP}
    # Get the feature dimension
    cvi.dim = length(sample)
    # Initialize the 2-D arrays with the correct feature dimension
    cvi.v = Array{T, 2}(undef, cvi.dim, 0)
end # setup!(cvi::PS, sample::Vector{T}) where {T<:RealFP}

function param_inc!(cvi::PS, sample::RealVector, label::Integer)
    # Get the internal label
    i_label = get_internal_label!(cvi.label_map, label)

    n_samples_new = cvi.n_samples + 1
    if isempty(cvi.v)
        setup!(cvi, sample)
    end

    if i_label > cvi.n_clusters
        n_new = 1
        v_new = sample
        if cvi.n_clusters == 0
            D_new = zeros(1, 1)
        else
            D_new = zeros(cvi.n_clusters + 1, cvi.n_clusters + 1)
            D_new[1:cvi.n_clusters, 1:cvi.n_clusters] = cvi.D
            d_column_new = zeros(cvi.n_clusters + 1)
            for jx = 1:cvi.n_clusters
                d_column_new[jx] = sum((v_new - cvi.v[:, jx]).^2)
            end
            D_new[:, i_label] = d_column_new
            D_new[i_label, :] = transpose(d_column_new)
        end
        # Update 1-D parameters with a push
        cvi.n_clusters += 1
        push!(cvi.n, n_new)
        # Update 2-D parameters with appending and reassignment
        cvi.v = [cvi.v v_new]
        cvi.D = D_new
    else
        n_new = cvi.n[i_label] + 1
        v_new = (1 - 1/n_new) .* cvi.v[:, i_label] + (1/n_new) .* sample
        d_column_new = zeros(cvi.n_clusters)
        for jx = 1:cvi.n_clusters
            if jx == i_label
                continue
            end
            d_column_new[jx] = sum((v_new - cvi.v[:, jx]).^2)
        end
        # Update parameters
        cvi.n[i_label] = n_new
        cvi.v[:, i_label] = v_new
        cvi.D[:, i_label] = d_column_new
        cvi.D[i_label, :] = transpose(d_column_new)
    end
    cvi.n_samples = n_samples_new
end # param_inc!(cvi::PS, sample::RealVector, label::Integer)

function param_batch!(cvi::PS, data::RealMatrix, labels::IntegerVector)
    cvi.dim, cvi.n_samples = size(data)
    # Take the average across all samples, but cast to 1-D vector
    u = unique(labels)
    cvi.n_clusters = length(u)
    cvi.n = zeros(Integer, cvi.n_clusters)
    cvi.v = zeros(cvi.dim, cvi.n_clusters)
    cvi.D = zeros(cvi.n_clusters, cvi.n_clusters)
    for ix = 1:cvi.n_clusters
        subset = data[:, findall(x->x==u[ix], labels)]
        cvi.n[ix] = size(subset, 2)
        cvi.v[1:cvi.dim, ix] = mean(subset, dims=2)
    end
    for ix = 1 : (cvi.n_clusters - 1)
        for jx = ix + 1 : cvi.n_clusters
            cvi.D[jx, ix] = sum((cvi.v[:, ix] - cvi.v[:, jx]).^2)
        end
    end
    cvi.D = cvi.D + transpose(cvi.D)
end # param_batch!(cvi::PS, data::RealMatrix, labels::IntegerVector)

function evaluate!(cvi::PS)
    if cvi.n_clusters > 1
        cvi.v_bar = vec(mean(cvi.v, dims=2))
        cvi.beta_t = 0.0
        cvi.PS_i = zeros(cvi.n_clusters)
        for ix = 1:cvi.n_clusters
            delta_v = cvi.v[:, ix] - cvi.v_bar
            cvi.beta_t = cvi.beta_t + transpose(delta_v) * delta_v
        end
        cvi.beta_t /= cvi.n_clusters
        n_max = maximum(cvi.n)
        for ix = 1:cvi.n_clusters
            d = cvi.D[:, ix]
            deleteat!(d, ix)
            cvi.PS_i[ix] = (cvi.n[ix] / n_max) - exp(-minimum(d) / cvi.beta_t)
        end
        cvi.criterion_value = sum(cvi.PS_i)
    end
end # evaluate(cvi::PS)


text:
using SymbolicML

using Test
using Statistics
using StatsBase

include("functions/runtests_summary.jl")


text:
import GeometryTypes
const gt = GeometryTypes

@testset "GeometryUtils" begin

    @testset "rotate2d_left90 Tests" begin
        @test sb.rotate2d_left90([1, 1]) == [-1, 1]
        @test sb.rotate2d_left90([0, 0]) == [0, 0]
        @test sb.rotate2d_left90([1, 0]) == [0, 1]
        @test sb.rotate2d_left90([0, 1]) == [-1, 0]
        @test_throws AssertionError sb.rotate2d_left90([0, 1, 3])
    end;

    @testset "rotate2d" begin
        @test sb.rotate2d([1, 0], π/2) ≈ [0, 1]
        @test sb.rotate2d([1, 0], 3π/2) ≈ [0, -1]
        @test sb.rotate2d([1, 1], -π/4) ≈ [sqrt(2), 0]
        @test_throws AssertionError sb.rotate2d([0, 1, 3], π/2)
    end

    @testset "normal2d Tests" begin
        @test sb.normal2d([1, 1]) == [-1, 1]
        @test sb.normal2d([0, 0]) == [0, 0]
        @test sb.normal2d([1, 0]) == [0, 1]
        @test sb.normal2d([0, 1]) == [-1, 0]
        @test_throws AssertionError sb.normal2d([0, 1, 3])
    end;

    @testset "vector Tests" begin
        @test sb.vector([2], [5]) == [3]
        @test sb.vector([2, 1], [5, -1]) == [3, -2]
    end;

    @testset "distance Tests" begin
        @test sb.distance(gt.Point2(0), gt.Point2(0)) == 0
        @test sb.distance(gt.Point2(1.3), gt.Point2(1.3)) == 0
        @test sb.distance(gt.Point2(0,1), gt.Point2(1,1)) == 1
        @test sb.distance(gt.Point2(0,-1), gt.Point2(-1,-1)) == 1
        @test sb.distance(gt.Point3(0,-1,0), gt.Point3(-1,-1,0)) == 1
        @test sb.distance(gt.Point{5}([0,-1,0,5,5]), gt.Point{5}([-1,-1,0,5,5])) == 1
        @test sb.distance([0,-1], [-1,-1]) == 1
        @test sb.distance([0,-1,0,5,5], [-1,-1,0,5,5]) == 1
    end;

    @testset "contained Tests" begin
        @test sb.contained(gt.Circle(gt.Point2(0, 0), 0), gt.Circle(gt.Point2(1, 1),1)) == false # point outside of the gt.Circle
        @test sb.contained(gt.Circle(gt.Point2(0, 0), 0), gt.Circle(gt.Point2(1, 1),2)) == true # point insinde of the gt.Circle
        @test sb.contained(gt.Circle(gt.Point2(1, 0), 0), gt.Circle(gt.Point2(1, 1),1)) == true # point on the gt.Circle ring
        @test sb.contained(gt.Circle(gt.Point2(3.5, -0.5), 0.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == false # gt.Circle with center on the gt.Circle ring and overlapping
        @test sb.contained(gt.Circle(gt.Point2(4, -0.5), 0.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == true # gt.Circle insinde other gt.Circle, touching on gt.Circle ring
        @test sb.contained(gt.Circle(gt.Point2(4.5, -0.5), 0.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == true # gt.Circle insinde other gt.Circle
        @test sb.contained(gt.Circle(gt.Point2(5, -0.5), 1.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == true # same gt.Circles
        @test sb.contained(gt.Circle(gt.Point2(2.5, -0.5), 1.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == false # disjunct gt.Circles
        @test sb.contained(gt.Circle(gt.Point2(3, -0.5), 1.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == false # touching in one point
    end;

    @testset "outer_tangent Tests" begin
        @test string(sb.outer_tangent(gt.Point2(0), gt.Circle(gt.Point(3,0),1))) == "GeometryTypes.Point{2,Float64}[[0.0, 0.0], [2.66667, 0.942809]]"
        @test string(sb.outer_tangent(gt.Circle(gt.Point(3,0),1), gt.Point2(0))) == "GeometryTypes.Point{2,Float64}[[2.66667, -0.942809], [0.0, 0.0]]"
        @test string(sb.outer_tangent(gt.Circle(gt.Point2(0), 1), gt.Circle(gt.Point(0, 3), 1))) == string(gt.LineSegment(gt.Point2(-1.0, 1.22465e-16), gt.Point(-1.0, 3.0)))
        @test string(sb.outer_tangent(gt.Circle(gt.Point2(0), 1), gt.Circle(gt.Point(0, 3), 1), true)) == string(gt.LineSegment(gt.Point(-1.0, 3.0), gt.Point2(-1.0, 1.22465e-16)))
        @test string(sb.outer_tangent(gt.Circle(gt.Point(0, 3), 1), gt.Circle(gt.Point2(0), 1))) == string(gt.LineSegment(gt.Point2(1.0, 3.0), gt.Point(1.0, 0.0)))
        @test string(sb.outer_tangent(gt.Circle(gt.Point2(0.5, 0.5), 1.5), gt.Circle(gt.Point(5.5, -0.5), 2.0))) == string(gt.LineSegment(gt.Point2(0.648526, 1.99263), gt.Point(5.69803, 1.49017)))
        @test string(sb.outer_tangent(gt.Circle(gt.Point(5.5, -0.5), 2.0), gt.Circle(gt.Point2(0.5, 0.5), 1.5))) == string(gt.LineSegment(gt.Point2(4.91735, -2.41325), gt.Point(0.0630127, -0.934936)))
        @test sb.outer_tangent(gt.Circle(gt.Point2(4, -0.5), 0.5), gt.Circle(gt.Point2(5, -0.5), 1.5)) == gt.LineSegment(gt.Point2(0.0), gt.Point2(0.0))
        @test sb.outer_tangent(gt.Circle(gt.Point2(5, -0.5), 1.5), gt.Circle(gt.Point2(4, -0.5), 0.5)) == gt.LineSegment(gt.Point2(0.0), gt.Point2(0.0))
    end;

    @testset "intersect gt.LineSegment/gt.LineSegment Tests" begin
        @test sb.intersection(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.LineSegment(gt.Point2(3.0, 2.0), gt.Point2(2.0, 3.0))) == (true, gt.Point2(2.5, 2.5))
        @test sb.intersection(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.LineSegment(gt.Point2(3.0, 2.0), gt.Point2(2.0, 3.0)), true) == (true, gt.Point2(2.5, 2.5))
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 1.0)), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(0.0, 4.0))) == (true, gt.Point2(0.0, 1.0))
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 1.0)), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(0.0, 4.0)), true) == (false, gt.Point2(0))
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 4.0)), gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 4.0))) == (false, gt.Point2(0))
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 4.0)), gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 4.0))) == (false, gt.Point2(0))
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point2(4.0, 4.0)), gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(5.0, 5.0))) == (false, gt.Point2(0))
    end;

    @testset "intersect gt.LineSegment/Rectangel Tests" begin
        @test sb.intersection(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.HyperRectangle{2,Float64}([0.0, 0.0], [7.0, 4.0])) == (true, gt.Point2(4.0, 4.0))
        @test sb.intersection(gt.LineSegment(gt.Point2(3.0, 3.0), gt.Point2(2.0, 2.0)), gt.HyperRectangle{2,Float64}([0.0, 0.0], [7.0, 4.0])) == (true, gt.Point2(0.0, 0.0))
        @test sb.intersection(gt.LineSegment(gt.Point2(3.0, 2.0), gt.Point2(3.0, 3.0)), gt.HyperRectangle{2,Float64}([0.0, 0.0], [7.0, 4.0])) == (true, gt.Point2(3.0, 4.0))
        @test sb.intersection(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.HyperRectangle{2,Float64}([-1.0, -2.0], [7.0, 4.0])) == (true, gt.Point2(2.0, 2.0))
        @test sb.intersection(gt.LineSegment(gt.Point2(8.0, 0.0), gt.Point2(9.0, 0.0)), gt.HyperRectangle{2,Float64}([0.0,0.0], [7.0, 4.0])) == (false, gt.Point2(0.0, 0.0))
    end;

    @testset "intersect gt.LineSegment/gt.Circle Tests" begin
        @test string(sb.intersection(gt.LineSegment(gt.Point2(0.0, 0.5), gt.Point2(3.0, 0.5)), gt.Circle(gt.Point2(1.5, 1.5), 1.5))) == string((true, [0.381966, 0.5], true, [2.61803, 0.5]))
        @test string(sb.intersection(gt.LineSegment(gt.Point2(0.0, 2.5), gt.Point2(3.0, 2.5)), gt.Circle(gt.Point2(1.5, 1.5), 1.5))) == string((true, [0.381966, 2.5], true, [2.61803, 2.5]))
        @test sb.intersection(gt.LineSegment(gt.Point2(-1.0, 1.0), gt.Point(1.0, 1.0)), gt.Circle(gt.Point2(0.0, 0.0), 1.0)) == (true, [0.0, 1.0], true, [0.0, 1.0])
        @test sb.intersection(gt.LineSegment(gt.Point2(-1.5, 1.5), gt.Point(1.0, 1.0)), gt.Circle(gt.Point2(0.0, 0.0), 1.0)) == (false, [0.0, 0.0],false,  [0.0, 0.0])
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point(2.0, 1.0)), gt.Circle(gt.Point2(0.0, 0.0), 1.0)) == (true, [0.0, 1.0], true, [0.0, 1.0])
        @test sb.intersection(gt.LineSegment(gt.Point2(1.0, 1.0), gt.Point(2.0, 1.0)), gt.Circle(gt.Point2(0.0, 0.0), 1.0), only_in_lineSegment=true) == (false, [0.0, 1.0],false,  [0.0, 1.0])
    end;

    @testset "touching Tests" begin
        @test sb.is_touching(gt.Circle(gt.Point2(0, 0), 0), gt.Circle(gt.Point2(1, 1), 1)) == false # point outside of the gt.Circle
        @test sb.is_touching(gt.Circle(gt.Point2(0, 0), 1), gt.Circle(gt.Point2(2, 2), 1)) == false # distinct gt.Circles
        @test sb.is_touching(gt.Circle(gt.Point2(1, 0), 0), gt.Circle(gt.Point2(1, 1), 1)) == true # point touching gt.Circle
        @test sb.is_touching(gt.Circle(gt.Point2(1, 1), 0), gt.Circle(gt.Point2(1, 1), 1)) == true # point inside gt.Circle
        @test sb.is_touching(gt.Circle(gt.Point2(0, 0), 1), gt.Circle(gt.Point2(0, 2), 1)) == true # touchinggt.Circle
        @test sb.is_touching(gt.Circle(gt.Point2(0, 0), 1), gt.Circle(gt.Point2(1, 1), 1)) == true # overlapping gt.Circles
        @test sb.is_touching(gt.Circle(gt.Point2(0, 0), 1), gt.Circle(gt.Point2(0, 0), 1)) == true # same gt.Circles
        @test sb.is_touching(gt.Point2(0, 0), 1, gt.Point2(0, 0), 1) == true # same gt.Circles
        @test sb.is_touching([0, 0], 1, [0, 0], 1) == true # same gt.Circles
        @test sb.is_touching([0, 0, 0, 0], 1, [0, 0, 0, 0], 1) == true # same gt.Circles
        @test_throws AssertionError sb.is_touching([0, 0, 0, 0], 1, [0, 0, 0], 1)
    end;

    @testset "does_intersect gt.LineSegment/gt.Circle Tests" begin
        @test sb.does_intersect(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.Circle(gt.Point2(0.0, 0.0), 1.0)) == true
        @test sb.does_intersect(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.Circle(gt.Point2(2.5, 2.5), 1.0)) == true
        @test sb.does_intersect(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.Circle(gt.Point2(10.0, 0.0), 1.0)) == false
        @test sb.does_intersect(gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0)), gt.Circle(gt.Point2(1.0, 1.0), 1.0)) == false
        @test sb.does_intersect(gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0)), gt.Circle(gt.Point2(1.0, 1.0), 1.0), r_is_outside = false) == true
        @test sb.does_intersect(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.Circle(gt.Point2(0.0, 0.0), 1.0), only_in_lineSegment = true) == false
        @test sb.does_intersect(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.Circle(gt.Point2(2.5, 2.5), 1.0), only_in_lineSegment = true) == true
        @test sb.does_intersect(gt.LineSegment(gt.Point2(2.0, 2.0), gt.Point2(3.0, 3.0)), gt.Circle(gt.Point2(10.0, 0.0), 1.0), only_in_lineSegment = true) == false
        @test sb.does_intersect(gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0)), gt.Circle(gt.Point2(1.0, 1.0), 1.0), only_in_lineSegment = true) == false
        @test sb.does_intersect(gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0)), gt.Circle(gt.Point2(1.0, 1.0), 1.0), only_in_lineSegment = true, r_is_outside = false) == true
    end;

    @testset "is_inside Point/gt.Circle Tests" begin
        @test sb.is_inside(gt.Point2(2.0, 2.0), gt.Circle(gt.Point2(0.0, 0.0), 1.0)) == false
        @test sb.is_inside(gt.Point2(1.0, 0.0), gt.Circle(gt.Point2(0.0, 0.0), 1.0)) == true
        @test sb.is_inside(gt.Point2(1.0, 0.0), gt.Circle(gt.Point2(0.0, 0.0), 1.0), r_is_outside = true) == false
    end;

    @testset "polygon_area" begin
        @test sb.polygon_area(Vector{Float64}[]) ≈ 0.0
        @test sb.polygon_area([[0, 0]]) ≈ 0.0
        @test sb.polygon_area([[0, 0], [0, 0]]) ≈ 0.0
        @test sb.polygon_area([[0, 0], [1, 0]]) ≈ 0.0
        @test sb.polygon_area([[0, 0], [1, 0], [0, 0]]) ≈ 0.0
        @test sb.polygon_area([[0, 0], [1, 0], [1, 1], [0, 1]]) ≈ 1.0
        @test sb.polygon_area([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]) ≈ 1.0
        @test sb.polygon_area([[20, 20], [21, 20], [21, 21], [20, 21]]) ≈ 1.0
        @test sb.polygon_area([[20, 20], [21, 20], [21, 21], [20, 21], [20, 20]]) ≈ 1.0
        @test sb.polygon_area([[0, 0], [2, 0], [2, 1], [1, 1], [1, 2], [0, 2]]) ≈ 3.0
        @test sb.polygon_area([[0, 0], [2, 0], [2, 1], [1, 1], [1, 2], [0, 2], [0, 0]]) ≈ 3.0
        @test sb.polygon_area([[20, 20], [22, 20], [22, 21], [21, 21], [21, 22], [20, 22], [20, 20]]) ≈ 3.0
    end;

    @testset "get_all_polygon_sides" begin
        @test sb.get_all_polygon_sides([gt.Point2(0, 0), gt.Point2(1, 0), gt.Point2(1, 1), gt.Point2(0, 1)]) == [gt.LineSegment(gt.Point2(0, 0), gt.Point2(1, 0)), gt.LineSegment(gt.Point2(1, 0), gt.Point2(1, 1)), gt.LineSegment(gt.Point2(1, 1), gt.Point2(0, 1)), gt.LineSegment(gt.Point2(0, 1), gt.Point2(0, 0))]
        @test sb.get_all_polygon_sides([gt.Point2(0, 0), gt.Point2(0, 1)]) == [gt.LineSegment(gt.Point2(0, 0), gt.Point2(0, 1)), gt.LineSegment(gt.Point2(0, 1), gt.Point2(0, 0))]
        @test sb.get_all_polygon_sides([gt.Point2(0, 0)]) == [gt.LineSegment(gt.Point2(0, 0), gt.Point2(0, 0))]
        @test sb.get_all_polygon_sides(gt.Point2{Float64}[]) == gt.LineSegment{gt.Point2{Float64}}[]
    end;

    @testset "is_inside point in polygon" begin
        # point fully inside, list of of points polygon
        @test sb.is_inside(gt.Point2(0.5, 0.5), [gt.Point2(0, 0), gt.Point2(1, 0), gt.Point2(1, 1), gt.Point2(0, 1)]) == true
        # point fully outside, list of of points polygon
        @test sb.is_inside(gt.Point2(1.5, 1.5), [gt.Point2(0, 0), gt.Point2(1, 0), gt.Point2(1, 1), gt.Point2(0, 1)]) == false
        # point on edge, edge is outside (default)
        @test sb.is_inside(gt.Point2(0.0, 0.0), [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)]) == true
        polygon = [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0)]
        for r in 1:4
            for p in [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(3.0, 0.0), gt.Point2(3.0, 1.0), gt.Point2(3.0, 2.0), gt.Point2(3.0, 3.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 3.0), gt.Point2(0.0, 3.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 1.0)]
                @test sb.is_inside(p, polygon) == false
            end
            for p in [gt.Point2(1.5, 1.5), gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0), gt.Point2(1.5, 2.0), gt.Point2(1.5, 1.0), gt.Point2(1.0, 1.5), gt.Point2(2.0, 1.5)]
                @test sb.is_inside(p, polygon) == true
            end
            polygon = circshift(polygon, 1)
        end
        @test sb.is_inside(gt.Point2(0.0, 0.0), [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 2.0)]) == false

        @test sb.is_inside(gt.Point2(2.5, 1.5), [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 3.0)]) == false
        @test sb.is_inside(gt.Point2(2.6, 1.4), [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 3.0)]) == false
        @test sb.is_inside(gt.Point2(2.4, 1.6), [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 3.0)]) == false

        @test sb.is_inside(gt.Point2(1.0, 1.0), [gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)]) == false
        @test sb.is_inside(gt.Point2(3.0, 3.0), [gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)]) == true
        @test sb.is_inside(gt.Point2(1.0, 5.0), [gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)]) == false

        @test sb.is_inside(gt.Point2(2.0, 3.0), [gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)]) == true
        @test sb.is_inside(gt.Point2(3.0, 2.0), [gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)]) == true

        @test sb.is_inside(gt.Point2(2.0, 3.0), [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 5.0)]) == true
        @test sb.is_inside(gt.Point2(4.0, 1.0), [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 5.0)]) == false
        @test sb.is_inside(gt.Point2(4.0, 5.0), [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 5.0)]) == false

        @test sb.is_inside(gt.Point2(1.0, 0.0), [gt.Point2(0.0, 0.0), gt.Point2(-1.0, 0.0), gt.Point2(-1.0, -1.0), gt.Point2(0.0, -1.0)]) == false

        @test sb.is_inside(gt.Point2(0.0), [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == false

        @test sb.is_inside(gt.Point2(3.0, 4.0), [gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0)]) == false
    end;

    @testset "is_inside point in rectangle" begin
        # point fully inside, list of of points polygon
        @test sb.is_inside(gt.Point2(1.0, 1.0), gt.HyperRectangle([gt.Point2(0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(2.0, 0.0)])) == true
        @test sb.is_inside(gt.Point2(5.0, 5.0), gt.HyperRectangle([gt.Point2(0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(2.0, 0.0)])) == false
        @test sb.is_inside(gt.Point2(0.0, 0.0), gt.HyperRectangle([gt.Point2(0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(2.0, 0.0)])) == true
        @test sb.is_inside(gt.Point2(1.0, 0.0), gt.HyperRectangle([gt.Point2(0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(2.0, 0.0)])) == true
        @test sb.is_inside(gt.Point2(7.0, 7.0), gt.HyperRectangle([gt.Point2(6.0, 6.0), gt.Point2(8.0, 6.0), gt.Point2(8.0, 8.0), gt.Point2(6.0, 8.0)])) == true
        @test sb.is_inside(gt.Point2(8.0, 6.0), gt.HyperRectangle([gt.Point2(6.0, 6.0), gt.Point2(8.0, 6.0), gt.Point2(8.0, 8.0), gt.Point2(6.0, 8.0)])) == true
        @test sb.is_inside(gt.Point2(1.0, 1.0), gt.HyperRectangle([gt.Point2(6.0, 6.0), gt.Point2(8.0, 6.0), gt.Point2(8.0, 8.0), gt.Point2(6.0, 8.0)])) == false
    end;

    @testset "polygon_intersection" begin

        function check_all_combination(polygonA::Vector{<:gt.Point2}, polygonB::Vector{<:gt.Point2}, solution::Vector{<:gt.Point2})
            sol_length = length(solution)
            for itB in 1:length(polygonB)
                for itA in 1:length(polygonA)
                    res = sb.polygon_intersection(polygonA, polygonB)
                    if sol_length == 0
                        @test res == solution
                        continue
                    else
                        res = res[1]
                    end
                    correct = any(map(it -> circshift(solution, it) ≈ res, 1:length(solution)))
                    if !correct
                        println(polygonA)
                        println(polygonB)
                        println(res)
                        println(solution)
                    end
                    # @test length(res) == sol_length && all(map(e -> e in solution, res))
                    @test correct
                    polygonA = circshift(polygonA, 1)
                end
                polygonB = circshift(polygonB, 1)
            end
        end
        # same polygons
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0), gt.Point2(0.0, 0.0)])

        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)])

        check_all_combination([gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0)])

        check_all_combination([gt.Point2(-1.0, -1.0), gt.Point2(2.0, -1.0), gt.Point2(2.0, 2.0), gt.Point2(-1.0, 2.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)])
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(-1.0, -1.0), gt.Point2(2.0, -1.0), gt.Point2(2.0, 2.0), gt.Point2(-1.0, 2.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)])

        check_all_combination([gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 5.0)], [gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)], [gt.Point2(2.5,2.5), gt.Point2(3.0, 3.0), gt.Point2(2.5, 3.5), gt.Point2(2.0, 3.0)])
        check_all_combination([gt.Point2(2.0, 3.0), gt.Point2(4.0, 1.0), gt.Point2(4.0, 5.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 5.0)], [gt.Point2(2.0, 3.0), gt.Point2(2.5,2.5), gt.Point2(3.0, 3.0), gt.Point2(2.5, 3.5)])

        # C - correct
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(0.0, 4.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0)], [gt.Point2(1.5, 1.0), gt.Point2(2.0, 2.0), gt.Point2(2.5, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 1.0)])
        check_all_combination([gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(0.0, 4.0)], [gt.Point2(1.0, 1.0), gt.Point2(1.5, 1.0), gt.Point2(2.0, 2.0), gt.Point2(2.5, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0)])

        # D - correct
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 3.0), gt.Point2(0.0, 3.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0)], [gt.Point2(1.5, 1.0), gt.Point2(2.0, 2.0), gt.Point2(2.5, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 1.0)])
        check_all_combination([gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 2.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 3.0), gt.Point2(0.0, 3.0)], [gt.Point2(1.0, 1.0), gt.Point2(1.5, 1.0), gt.Point2(2.0, 2.0), gt.Point2(2.5, 1.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 3.0)])

        # K - infinite loop for two combinations
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0), gt.Point2(0.0, 0.0)])
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)])

        # L - correct
        check_all_combination([gt.Point2(0.0, 2.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 2.0)], [gt.Point2(0.0, 1.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 2.0)], [gt.Point2(0.4, 1.2), gt.Point2(0.6667, 0.6667), gt.Point2(1.2, 0.4), gt.Point2(2.0, 2.0)])
        check_all_combination([gt.Point2(0.0, 1.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 2.0)], [gt.Point2(0.0, 2.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 2.0)], [gt.Point2(0.4, 1.2), gt.Point2(0.6667, 0.6667), gt.Point2(1.2, 0.4), gt.Point2(2.0, 2.0)]) # [gt.Point2(2.0, 2.0)]

        # correct - Arrrgh - Error
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(-1.0, 0.0), gt.Point2(-1.0, -1.0), gt.Point2(0.0, -1.0)], gt.Point2[])
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(-1.0, 0.0), gt.Point2(-1.0, -1.0), gt.Point2(0.0, -1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], gt.Point2[])

        # A - correct
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(3.0, 4.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 4.0), gt.Point2(0.0, 4.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(0.0, 4.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(3.0, 4.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 4.0), gt.Point2(0.0, 4.0)])
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(0.0, 4.0)], [gt.Point2(0.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(3.0, 4.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 4.0), gt.Point2(0.0, 4.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(4.0, 4.0), gt.Point2(3.0, 4.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 4.0), gt.Point2(0.0, 4.0)])

        # B correct
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(2.0, 3.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(2.0, 2.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(2.0, 2.0)])
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(2.0, 2.0)], [gt.Point2(0.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(2.0, 3.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(3.0, 0.0), gt.Point2(4.0, 0.0), gt.Point2(2.0, 2.0)])

        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)])
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)])

        # H
        @test sb.polygon_intersection([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(2.0, 0.0), gt.Point2(3.0, 0.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 1.0)]) == Vector[]
        @test sb.polygon_intersection([gt.Point2(2.0, 0.0), gt.Point2(3.0, 0.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)]) == Vector[]

        # I - infinite loop
        check_all_combination([gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], [gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], gt.Point2[])
        check_all_combination([gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(0.0, 1.0)], gt.Point2[])

        # M
        check_all_combination([gt.Point2(0.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)], [gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 3.0)], [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0)])
        check_all_combination([gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 3.0)], [gt.Point2(0.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 2.0), gt.Point2(1.0, 2.0)])

        # complicated tests

        # F - seperatist - correct
        @test sb.polygon_intersection([gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]

        # no intersection
        @test sb.polygon_intersection(gt.Point2{Float64}[], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)]) == Vector{gt.Point2{Float64}}[]
        @test sb.polygon_intersection([gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)], gt.Point2{Float64}[]) == Vector{gt.Point2{Float64}}[]
        @test sb.polygon_intersection(gt.Point2{Float64}[], gt.Point2{Float64}[]) == Vector{gt.Point2{Float64}}[]

        @test sb.polygon_intersection([gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]

        @test sb.polygon_intersection([gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(3.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]

        # G - seperatist
        @test sb.polygon_intersection([gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]

        @test sb.polygon_intersection([gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(0.0, 2.0), gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0)], [gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]

        @test sb.polygon_intersection([gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(5.0, 1.0), gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
        @test sb.polygon_intersection([gt.Point2(3.0, 4.0), gt.Point2(1.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(5.0, 1.0)], [gt.Point2(0.0, 0.0), gt.Point2(6.0, 0.0), gt.Point2(6.0, 2.0), gt.Point2(0.0, 2.0)]) == [[gt.Point2(2.0, 2.0), gt.Point2(1.6667, 2.0), gt.Point2(1.0, 1.0)], [gt.Point2(4.3333, 2.0), gt.Point2(4.0, 2.0), gt.Point2(5.0, 1.0)]]
    end;

    @testset "polygon_union" begin
        @test sb.polygon_union([gt.Point2{Float64}[], [gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0)]]) == [[gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0)]]
        @test sb.polygon_union([[gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0)], gt.Point2{Float64}[]]) == [[gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0)]]

        @test sb.polygon_union([gt.Point2{Float64}[]]) == [gt.Point2{Float64}[]]
        @test sb.polygon_union(Vector{gt.Point2{Float64}}[]) == [gt.Point2{Float64}[]]
        @test sb.polygon_union([[gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0)]]) == [[gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0), gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0)]]

        @test sb.polygon_union([[gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0)]]) == [[gt.Point2(0.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(0.0, 1.0)]]
        @test sb.polygon_union([[gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0), gt.Point2(0.0, 0.0)], [gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0)]]) == [[gt.Point2(0.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(0.0, 1.0)]]
        @test sb.polygon_union([[gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)], [gt.Point2(1.0, 0.0), gt.Point2(2.0, 0.0), gt.Point2(2.0, 1.0), gt.Point2(1.0, 1.0)], [gt.Point2(2.0, 0.0), gt.Point2(3.0, 0.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 1.0)]]) == [[gt.Point2(0.0, 0.0), gt.Point2(3.0, 0.0), gt.Point2(3.0, 1.0), gt.Point2(0.0, 1.0)]]
        @test sb.polygon_union([[gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0), gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0)],  [gt.Point2(2.0, 0.0), gt.Point2(3.0, 0.0), gt.Point2(3.0, 1.0), gt.Point2(2.0, 1.0)]]) == [[gt.Point2(1.0, 1.0), gt.Point2(0.0, 1.0), gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0)],  [gt.Point2(3.0, 1.0), gt.Point2(2.0, 1.0), gt.Point2(2.0, 0.0), gt.Point2(3.0, 0.0)]]
        @test sb.polygon_union([[gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 3.0)],  [gt.Point2(2.0, 2.0), gt.Point2(4.0, 2.0), gt.Point2(4.0, 4.0), gt.Point2(2.0, 4.0)]]) == [[gt.Point2(3.0, 2.0), gt.Point2(4.0, 2.0), gt.Point2(4.0, 4.0), gt.Point2(2.0, 4.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 3.0), gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0)]]
        @test sb.polygon_union([[gt.Point2(1.0, 1.0), gt.Point2(3.0, 1.0), gt.Point2(3.0, 3.0), gt.Point2(1.0, 3.0)],  [gt.Point2(2.0, 2.0), gt.Point2(4.0, 2.0), gt.Point2(4.0, 4.0), gt.Point2(2.0, 4.0)], [gt.Point2(3.0, 1.0), gt.Point2(5.0, 1.0), gt.Point2(5.0, 3.0), gt.Point2(3.0, 3.0)]]) == [[gt.Point2(5.0, 3.0), gt.Point2(4.0, 3.0), gt.Point2(4.0, 4.0), gt.Point2(2.0, 4.0), gt.Point2(2.0, 3.0), gt.Point2(1.0, 3.0), gt.Point2(1.0, 1.0), gt.Point2(5.0, 1.0)]]
    end

    @testset "is_on_line_segment" begin
        @test sb.is_on_LineSegment(gt.Point2(0.0, 0.0), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0))) == true
        @test sb.is_on_LineSegment(gt.Point2(0.5, 0.0), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0))) == true
        @test sb.is_on_LineSegment(gt.Point2(1.0, 0.0), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0))) == true
        @test sb.is_on_LineSegment(gt.Point2(2.0, 0.0), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0))) == false
        @test sb.is_on_LineSegment(gt.Point2(2.0, 2.0), gt.LineSegment(gt.Point2(0.0, 0.0), gt.Point2(1.0, 0.0))) == false
    end;

end


text:
module XSBenchThreads

# struct
export Inputs, SimulationData
# functions
export read_CLI, grid_init_do_not_profile

mutable struct Inputs
    nthreads::Int32
    n_isotopes::Int64
    n_gridpoints::Int64
    lookups::Int32
    HM::String
    grid_type::String
    hash_bins::Int32
    particles::Int32
    simulation_method::String
    binary_mode::String
    kernel_id::Int32

    Inputs() = new(
        Threads.nthreads(),
        Int64(355),
        Int64(11303),
        Int32(34),
        "large",
        "unionized",
        Int32(10_000),
        Int32(500_000),
        "history_based",
        "none",
        Int32(0),
    )
end

mutable struct NuclideGridPoint
    energy::Float64
    total_xs::Float64
    elastic_xs::Float64
    absorption_xs::Float64
    fission_xs::Float64
    nu_fission::Float64

    NuclideGridPoint() = new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
end

mutable struct SimulationData
    num_nucs::Array{Int32,1}
    length_nuclide_grid::Int32
    nuclide_grid::Array{NuclideGridPoint,1}

    SimulationData() = new(Array{Int32,1}[], Int32(0), Array{NuclideGridPoint,1}[])
end

include("io.jl")
include("GridInit.jl")

end


text:
module Constraint

using LinearAlgebra

using MutableArithmetics
const MA = MutableArithmetics

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

import MultivariatePolynomials
const MP = MultivariatePolynomials
import MultivariateBases
import SemialgebraicSets
import MultivariateMoments
import PolyJuMP
import SumOfSquares
const SOS = SumOfSquares
const Certificate = SOS.Certificate

# Symmetric PSD matrix bridges
include("empty.jl")
include("psd2x2.jl")
include("diagonally_dominant.jl")
include("scaled_diagonally_dominant.jl")

# SOS polynomial bridges
include("utilities.jl")
include("sos_polynomial.jl")
include("sos_polynomial_in_semialgebraic_set.jl")

end