text:
# T , TSSM_HANDLE
println("including tssm_schroedinger_rotating.jl for type ", T)
for (METHOD, SUF, COMPLEX_METHOD, DIM) in (
                 (:SchroedingerRotating2D, :_schroedinger_rotating_2d, true, 2 ),
                 (:SchroedingerRotating3D, :_schroedinger_rotating_3d, true, 3 ),
                 (:SchroedingerRotatingReal2D, :_schroedinger_rotating_real_2d, false, 2 ),
                 (:SchroedingerRotatingReal3D, :_schroedinger_rotating_real_3d, false, 3 ),
                )
println("    ", METHOD)
if T == :Float128
    SUF = Symbol(SUF, "_wf128")
end
WF = Symbol(:Wf,METHOD)

if DIM==2
@eval begin
    function ($METHOD)(T::Type{($T)}, nx::Integer, xmin::Real, xmax::Real,
                                      ny::Integer, ymin::Real, ymax::Real;
                       Omega::Real=0.0,
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_2D,
                       potential_t::Function=none_3D,
                       potential_t_derivative::Function=none_3D,
                       cubic_coupling::Real=0.0)
        with_potential = potential!=none_2D
        with_potential_t = potential_t!=none_3D
        with_potential_t_derivative = potential_t_derivative!=none_3D
        V_c = cfunction_check_return_type(potential, ($T), (($T),($T)))
        V_t_c = cfunction_check_return_type(potential_t, ($T), (($T),($T),($T)))
        V_t_derivative_c = cfunction_check_return_type(potential_t_derivative, ($T), (($T),($T),($T)))
        c = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))), Ptr{Nothing},
                   (Int32, ($T), ($T), Int32, ($T), ($T),
                   ($T), ($T), ($T), Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, ($T)),
                   nx, xmin, xmax, ny, ymin, ymax,
                   Omega, hbar, mass, V_c, with_potential,
                   V_t_c, with_potential_t, V_t_derivative_c, with_potential_t_derivative,
                   cubic_coupling)
        m = ($METHOD){$T}(c)
        finalizer(m, x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                        $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m) )
        m
    end
end # eval
if T == :Float64
@eval begin
    function ($METHOD)(nx::Integer, xmin::Real, xmax::Real,
                       ny::Integer, ymin::Real, ymax::Real;
                       Omega::Real=0.0,
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_2D,
                       potential_t::Function=none_3D,
                       potential_t_derivative::Function=none_3D,
                       cubic_coupling::Real=0.0)
        ($METHOD)(($T), nx, xmin, xmax,  ny, ymin, ymax,
                       Omega=Omega, hbar=hbar, mass=mass, potential=potential,
                       potential_t=potential_t, potential_t_derivative=potential_t_derivative,
                       cubic_coupling=cubic_coupling)
    end
end # eval
end #if
elseif DIM==3
@eval begin
    function ($METHOD)(T::Type{($T)}, nx::Integer, xmin::Real, xmax::Real,
                                      ny::Integer, ymin::Real, ymax::Real,
                                      nz::Integer, zmin::Real, zmax::Real;
                       Omega::Real=0.0,
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_3D,
                       potential_t::Function=none_4D,
                       potential_t_derivative::Function=none_4D,
                       cubic_coupling::Real=0.0)
        with_potential = potential!=none_3D
        with_potential_t = potential_t!=none_4D
        with_potential_t_derivative = potential_t_derivative!=none_4D
        V_c = cfunction_check_return_type(potential, ($T), (($T),($T),($T)))
        V_t_c = cfunction_check_return_type(potential_t, ($T), (($T),($T),($T),($T)))
        V_t_derivative_c = cfunction_check_return_type(potential_t_derivative, ($T), (($T),($T),($T),($T)))
        c = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))), Ptr{Nothing},
                   (Int32, ($T), ($T), Int32, ($T), ($T), Int32, ($T), ($T),
                   ($T), ($T), ($T), Ptr{Nothing}, Bool,  Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, ($T)),
                   nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax,
                   Omega, hbar, mass, V_c, with_potential,
                   V_t_c, with_potential_t, V_t_derivative_c, with_potential_t_derivative,
                   cubic_coupling)
        m = ($METHOD){$T}(c)
        finalizer(m, x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                        $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m) )
        m
    end
end # eval
if T == :Float64
@eval begin
    function ($METHOD)(nx::Integer, xmin::Real, xmax::Real,
                       ny::Integer, ymin::Real, ymax::Real,
                       nz::Integer, zmin::Real, zmax::Real;
                       Omega::Real=0.0,
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_3D,
                       potential_t::Function=none_4D,
                       potential_t_derivative::Function=none_4D,
                       cubic_coupling::Real=0.0)
        ($METHOD)(($T), nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax,
                       Omega=Omega, hbar=hbar, mass=mass, potential=potential,
                       potential_t=potential_t, potential_t_derivative=potential_t_derivative,
                       cubic_coupling=cubic_coupling)
    end
end # eval
end #if

end # if

@eval begin

    function get_Omega(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_Omega",SUF))), ($T),
             (Ptr{Nothing}, ), m.m )
    end

    function propagate_C!(psi::($WF){$T}, dt::Number)
        ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_C_wf",SUF))), Nothing,
                (Ptr{Nothing}, Complex{$T},), psi.p, dt)
    end

    function propagate_C_derivative!(this::($WF){$T}, other::($WF){$T},
                         dt::Number)
       if this.m ≠ other.m
           error("this and other must belong to the same method")
       end
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_C_derivative_wf",SUF))), Nothing,
                (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                 this.p, other.p, dt)
    end

    function add_apply_C!(this::($WF){$T}, other::($WF){$T},
                         coefficient::Number=1.0)
       if this.m ≠ other.m
           error("this and other must belong to the same method")
       end
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_apply_C_wf",SUF))), Nothing,
                (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                 this.p, other.p, coefficient)
    end

    function imaginary_time_propagate_C!(psi::($WF){$T}, dt::Number)
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"imaginary_time_propagate_C_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Complex{($T)},), psi.p, dt)
    end

    function get_nx(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nx",SUF))), Int32,
             (Ptr{Nothing}, ), m.m )
    end

    function get_xmin(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_xmin",SUF))), ($T),
             (Ptr{Nothing}, ), m.m )
    end

    function get_xmax(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_xmax",SUF))), $T,
             (Ptr{Nothing}, ), m.m )
    end

    function is_real_space_x(psi::($WF){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_x_wf",SUF))), Int32,
           (Ptr{Nothing},), psi.p) == 1
    end

    function is_frequency_space_x(psi::($WF){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_x_wf",SUF))), Int32,
           (Ptr{Nothing},), psi.p) != 1
    end

    function to_real_space_x!(psi::($WF){$T})
         ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_real_space_x_wf",SUF))), Nothing,
                (Ptr{Nothing},), psi.p)
    end

    function to_frequency_space_x!(psi::($WF){$T})
         ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_frequency_space_x_wf",SUF))), Nothing,
                (Ptr{Nothing},), psi.p)
    end
end # eval

if DIM>=2
    @eval begin

        function get_ny(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ny",SUF))), Int32,
                 (Ptr{Nothing}, ), m.m )
        end

        function get_ymin(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ymin",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function get_ymax(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ymax",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end
        function is_real_space_y(psi::($WF){$T})
        ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_y_wf",SUF))), Int32,
            (Ptr{Nothing},), psi.p) == 1
        end

        function is_frequency_space_y(psi::($WF){$T})
        ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_y_wf",SUF))), Int32,
            (Ptr{Nothing},), psi.p) != 1
        end

        function to_real_space_y!(psi::($WF){$T})
            ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_real_space_y_wf",SUF))), Nothing,
                    (Ptr{Nothing},), psi.p)
        end

        function to_frequency_space_y!(psi::($WF){$T})
            ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_frequency_space_y_wf",SUF))), Nothing,
                    (Ptr{Nothing},), psi.p)
        end
    end # eval
end

if DIM>=3
    @eval begin
        function get_nz(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nz",SUF))), Int32,
                 (Ptr{Nothing}, ), m.m )
        end

        function get_zmin(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_zmin",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function get_zmax(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_zmax",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function is_real_space_z(psi::($WF){$T})
        ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_z_wf",SUF))), Int32,
            (Ptr{Nothing},), psi.p) == 1
        end

        function is_frequency_space_z(psi::($WF){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_z_wf",SUF))), Int32,
            (Ptr{Nothing},), psi.p) != 1
        end

        function to_real_space_z!(psi::($WF){$T})
            ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_real_space_z_wf",SUF))), Nothing,
                    (Ptr{Nothing},), psi.p)
        end

        function to_frequency_space_z!(psi::($WF){$T})
            ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_frequency_space_z_wf",SUF))), Nothing,
                    (Ptr{Nothing},), psi.p)
        end

    end # eval
end

end # for


text:
using NonlinearBandits, Test

X = [3 1 2]
x = [1]
@test nearest_neighbours(x, X) == ([2, 3, 1], [0.0, 1.0, 2.0])

ext: jl


text:
"""
    Chain(layers...)

Chain multiple layers / functions together, so that they are called in sequence
on a given input.

```julia
m = Chain(x -> x^2, x -> x+1)
m(5) == 26

m = Chain(Dense(10, 5), Dense(5, 2))
x = rand(10)
m(x) == m[2](m[1](x))
```

`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.
"""
type Chain
  layers::Vector{Any}
  Chain(xs...) = new([xs...])
end

@forward Chain.layers Base.getindex, Base.first, Base.last, Base.endof, Base.push!
@forward Chain.layers Base.start, Base.next, Base.done

children(c::Chain) = c.layers
mapchildren(f, c::Chain) = Chain(f.(c.layers)...)
adapt(T, c::Chain) = Chain(map(x -> adapt(T, x), c.layers)...)

(c::Chain)(x) = foldl((x, m) -> m(x), x, c.layers)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)

function Base.show(io::IO, c::Chain)
  print(io, "Chain(")
  join(io, c.layers, ", ")
  print(io, ")")
end

"""
    Dense(in::Integer, out::Integer, σ = identity)

Creates a traditional `Dense` layer with parameters `W` and `b`.

    y = σ.(W * x .+ b)

The input `x` must be a vector of length `in`, or a batch of vectors represented
as an `in × N` matrix. The out `y` will be a vector or batch of length `out`.

```julia
julia> d = Dense(5, 2)
Dense(5, 2)

julia> d(rand(5))
Tracked 2-element Array{Float64,1}:
  0.00257447
  -0.00449443
```
"""
struct Dense{F,S,T}
  W::S
  b::T
  σ::F
end

Dense(W, b) = Dense(W, b, identity)

function Dense(in::Integer, out::Integer, σ = identity;
               initW = glorot_uniform, initb = zeros)
  return Dense(param(initW(out, in)), param(initb(out)), σ)
end

treelike(Dense)

function (a::Dense)(x)
  W, b, σ = a.W, a.b, a.σ
  @fix σ.(W*x .+ b)
end

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

"""
    Diagonal(in::Integer)

Creates an element-wise linear transformation layer with learnable
vectors `α` and `β`:

    y = α .* x .+ β

The input `x` must be a array where `size(x, 1) == in`.
"""
struct Diagonal{T}
  α::T
  β::T
end

Diagonal(in::Integer; initα = ones, initβ = zeros) =
  Diagonal(param(initα(in)), param(initβ(in)))

treelike(Diagonal)

function (a::Diagonal)(x)
  α, β = a.α, a.β
  α.*x .+ β
end

function Base.show(io::IO, l::Diagonal)
  print(io, "Diagonal(", length(l.α), ")")
end


text:
using ContinuousWavelets, Wavelets, Interpolations, LinearAlgebra
using Test
using FFTW
using Logging
@testset "ContinuousWavelets.jl" begin
    # Write your tests here.
    include("basicTypesAndNumber.jl")
    include("deltaSpikes.jl")
    include("defaultProperties.jl")
    include("inversionTests.jl")
end
# TODO:
#       test averaging types
#            various extra dimensions
#            inverse is actually functional


text:
export projSD, projSDhis

mutable struct  projSDhis
	Jc::Array
	dJ::Array
	F::Array
	Dc::Array
	Rc::Array
	alphas::Array
	Active::Array
	stepNorm::Array
	lsIter::Array
	timeMisfit::Array
	timeReg::Array
	timeGradMisfit::Array
end

function getProjSDhis(maxIter)
	Jc             = zeros(maxIter+1)
	dJ             = zeros(maxIter+1)
	F              = zeros(maxIter+1)
	Dc             = []
	Rc             = zeros(maxIter+1)
	alphas         = zeros(maxIter+1)
	Active         = zeros(maxIter+1)
	stepNorm       = zeros(maxIter+1)
	lsIter         = zeros(Int,maxIter+1)
	timeMisfit     = zeros(maxIter+1,4)
	timeReg        = zeros(maxIter+1)
	timeGradMisfit = zeros(maxIter+1,2)

	return projSDhis(Jc,dJ,F,Dc,Rc,alphas,Active,stepNorm,lsIter,timeMisfit,timeReg,timeGradMisfit)
end

function updateHis!(iter::Int64,His::projSDhis,Jc::Real,dJ::Real,Fc,Dc,Rc::Real,alpha::Real,
					nActive::Int64,stepNorm::Real,lsIter::Int,timeMisfit::Vector,timeReg::Real)
	His.Jc[iter+1]            = Jc
	His.F[iter+1]             = Fc
	push!(His.Dc,Dc)
	His.Rc[iter+1]            = Rc
	His.alphas[iter+1]        = alpha
	His.Active[iter+1]        = nActive
	His.stepNorm[iter+1]      = stepNorm
	His.timeMisfit[iter+1,:] += timeMisfit
	His.timeReg[iter+1]      += timeReg
end

"""
	mc,Dc,outerFlag = projSD(mc,pInv::InverseParam,pMis, indFor = [], dumpResults::Function = dummy)

	(Projected) Steepest Descent method for solving

		min_x misfit(x) + regularizer(x) subject to  x in C

	where C is a convex set and a projection operator proj(x) needs to be provided.

	Input:

		mc::Vector          - intial guess for model
		pInv::InverseParam  - parameter for inversion
		pMis                - misfit terms
		indCredit           - indices of forward problems to work on
		dumpResults			- A function pointer for saving the results throughout the iterations.
							- We assume that dumpResults is dumpResults(mc,Dc,iter,pInv,pMis),
							- where mc is the recovered model, Dc is the predicted data.
							- If dumpResults is not given, nothing is done (dummy() is called).
		out::Int            - flag for output (-1: no output, 1: final status, 2: residual norm at each iteration)

	Output:
		mc                  - final model
		Dc                  - data
		outerFlag           - flag for convergence
		His                 - iteration history

"""
function  projSD(mc,pInv::InverseParam,pMis; proj=x->min.(max.(x,pInv.boundsLow),pInv.boundsHigh),
	indCredit=[], dumpResults::Function = dummy,out::Int=2)

	maxIter     = pInv.maxIter      #  Max. no. iterations.
	stepTol     = pInv.minUpdate    #  Step norm stopping tol.
	maxStep     = pInv.maxStep
	alpha       = pInv.alpha
	low         = pInv.boundsLow
	high        = pInv.boundsHigh

	His = getProjSDhis(maxIter)
	#---------------------------------------------------------------------------
	#  Initialization.
	#---------------------------------------------------------------------------
	mc = proj(mc)

	Active = (mc .<=low) .| (mc.>=high)  # Compute active set

	## evaluate function and derivatives
	sig,dsig = pInv.modelfun(mc)
	if isempty(indCredit)
		Dc,F,dF,d2F,pMis,tMis = computeMisfit(sig,pMis,true)
	else
		Dc,F,dF,d2F,pMis,tMis,indDebit = computeMisfit(sig,pMis,true,indCredit=indCredit)
	end
	dF = dsig'*dF

	# compute regularizer
	tReg = @elapsed R,dR,d2R = computeRegularizer(pInv.regularizer,mc,pInv.mref,pInv.MInv,alpha)

	# objective function
	Jc  = F  + R
	gc  = dF + dR

	F0 = F; J0 = Jc
	############################################################################
	##  Outer iteration.                                                        #
	############################################################################
	iter = 0
	outerFlag = -1; stepNorm=0.0

	outStr = @sprintf("%4s\t%08s\t%08s\t%08s\t%08s\t%08s\n",
					  	"i.LS", "F", "R","alpha[1]","Jc/J0","#Active")
	updateHis!(0,His,Jc,norm(projGrad(gc,mc,low,high)),F,Dc,R,alpha[1],count(!iszero, Active),0.0,-1,tMis,tReg)

	if out>=2; print(outStr); end
	f = open("projSD.out", "w")
	write(f, outStr)
	close(f)

	while outerFlag == -1

		iter += 1
		outStr = @sprintf("%3d.0\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3d\n",
		         iter, F, R,alpha[1],Jc/J0,count(!iszero, Active))
		if out>=2; print(outStr); end
		f = open("jInv.out", "a")
		write(f, outStr)
		close(f)

		# scale step
		if maximum(abs.(gc)) > maxStep; gc = gc./maximum(abs.(gc))*maxStep; end

		## Begin projected Armijo line search
		muLS = 1; lsIter = 1; mt = zeros(size(mc)); Jt = Jc
		while true
			mt = proj(mc - muLS*gc)
			## evaluate function
			sigt, = pInv.modelfun(mt)
			if isempty(indCredit)
				Dc,F,dF,d2F,pMis,tMis = computeMisfit(sigt,pMis,false)
			else
				Dc,F,dF,d2F,pMis,tMis,indDebit = computeMisfit(sigt,false,indCredit=indCredit)
			end
			His.timeMisfit[iter+1,:]+=tMis

			His.timeReg[iter+1] += @elapsed R,dR,d2R = computeRegularizer(pInv.regularizer,mt,pInv.mref,pInv.MInv,alpha)

			# objective function
			Jt  = F  + R
			if out>=2;
				println(@sprintf( "   .%d\t%3.2e\t%3.2e\t\t\t%3.2e",
			           lsIter, F,       R,       Jt/J0))
			end

			if Jt < Jc
			    break
			end
			muLS /=2; lsIter += 1
			if lsIter > 6
			    outerFlag = -2
				break
			end
		end
		## End Line search

		## Check for termination
		stepNorm = norm(mt-mc,Inf)
		mc = mt
		Jc = Jt

		sig, dsig = pInv.modelfun(mc)

		Active = (mc .<=low) .| (mc.>=high)  # Compute active set

		#  Check stopping criteria for outer iteration.
		updateHis!(iter,His,Jc,-1,F,Dc,R,alpha[1],count(!iszero, Active),stepNorm,lsIter,tMis,tReg)

		dumpResults(mc,Dc,iter,pInv,pMis);
		if stepNorm < stepTol
			outerFlag = 1
			break
		elseif iter >= maxIter
			break
		end
		# Evaluate gradient
		t = time_ns();
		if isempty(indCredit)
			dF = computeGradMisfit(sig,Dc,pMis)
		else
			dF = computeGradMisfit(sig,Dcp,pMis,indDebit)
		end
		His.timeGradMisfit[iter+1]+=(time_ns() - t)/1e+9;

		dF = dsig'*dF
		gc = dF + dR
		His.dJ[iter+1] = norm(projGrad(gc,mc,low,high))
	end # while outer_flag == 0

	if out>=1
		if outerFlag==-1
			println("projSD iterated maxIter=$maxIter times but reached only stepNorm of $(stepNorm) instead $(stepTol)." )
		elseif outerFlag==-2
			println("projSD stopped at iteration $iter with a line search fail.")
		elseif outerFlag==1
			println("projSD reached desired accuracy at iteration $iter.")
		end
	end

	return mc,Dc,outerFlag,His
end  # Optimization code


text:
module LinearResponse

    using Symbolics: variables
    using LinearAlgebra
    using Printf
    using Symbolics
    using ..HarmonicBalance
    using ..HC_wrapper
    using PyPlot
    using ProgressMeter
    using DocStringExtensions
    using Latexify

    import Base: *, show; export *, show

    include("LinearResponse/types.jl")
    include("LinearResponse/perturbations.jl")
    include("LinearResponse/jacobian_spectrum.jl")
    include("LinearResponse/linear_response.jl")

end


text:
dict = Dict(
    "first" => ["single"],
    "second" => ["dual1", "dual2"]
)

single_dict = Dict("first" => ["single"])
single_str = LightOSM.concatenate_exclusions(single_dict)
@test startswith(single_str, '[')
@test endswith(single_str, ']')
args_only = chop(single_str, head=1, tail=1)
args = split(args_only, "!~")
@test args[1] == "\"first\""
@test args[2] == "\"single\""

dual_dict = Dict("second" => ["dual1", "dual2"])
dual_str = LightOSM.concatenate_exclusions(dual_dict)
@test startswith(dual_str, '[')
@test endswith(dual_str, ']')
args_only = chop(dual_str, head=1, tail=1)
args = split(args_only, "!~")
@test args[1] == "\"second\""
@test args[2] == "\"dual1|dual2\""

combined_dict = dict = Dict(
    "first" => ["single"],
    "second" => ["dual1", "dual2"]
)
str = LightOSM.concatenate_exclusions(combined_dict)
@test str == single_str * dual_str || str == dual_str * single_str # either order


text:
println("Dm) Test of the RadiativeAuger module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-RadiativeAuger.sum")
setDefaults("method: continuum, Galerkin")           ## setDefaults("method: continuum, Galerkin")  "method: continuum, asymptotic Coulomb"
                                                     ## setDefaults("method: normalization, Ong-Russek")
setDefaults("method: normalization, pure sine")      ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")
setDefaults("unit: rate", "a.u.")

grid = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 1.0e-2, rbox = 10.0)

if  false
    # Green function for radiative Auger decay of 1s photoionized neon
    name             = "radiative Auger decay of 1s photoionized neon"
    refConfigs       = [Configuration("1s 2s^2 2p^6")]
    levelSymmetries  = [LevelSymmetry(1//2, Basics.minus), LevelSymmetry(1//2, Basics.plus)]
    greenSettings    = GreenSettings(3, [0, 1], 0.01, true, LevelSelection())
    greenRep         = Representation(name, Nuclear.Model(10.), Radial.Grid(true), refConfigs,
                                      GreenExpansion( AtomicState.DampedSpaceCI(), Basics.DeExciteSingleElectron(), levelSymmetries, 3, greenSettings) )
    greenOut         = generate(greenRep, output=true)
    raGreen          = greenOut["Green channels"]

elseif true
    # Radiative Auger decay of 1s photoionized neon
    raSettings = RadiativeAuger.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], raGreen, 2, 3, true, CoulombInteraction(),
                                         LineSelection(true, indexPairs=[(1,1), (1,2)]))

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(10.),
                            initialConfigs=[Configuration("1s 2s^2 2p^6")],
                            finalConfigs  =[Configuration("1s^2 2s^2 2p^4"), Configuration("1s^2 2s 2p^5"), Configuration("1s^2 2p^6") ],
                            process = RAuger(), processSettings = raSettings )

    wb = perform(wa)
end
setDefaults("print summary: close", "")


text:
function BoxCox(x, lambda, IFinv = 0)
    if IFinv == 0
        if lambda == 0
            return [log(x), max(float(x).^(-1), 1e5)]
        else
            return [(expm1(lambda * log(x)))/lambda, max(float(x).^(lambda .-1), 1e5)]
        end
    else # compute the inverse function and inverse derivative
        if lambda == 0
            return [exp(x), exp(x)]
        else
            return [(lambda * x + 1)^(1/lambda), (lambda * x + 1)^(1/lambda-1)/lambda]
        end
    end
end

function ArandaOrdaz(x, lambda)
    # return function values and deivatives
    @assert 0 <= x <= 1
    if lambda == 0
        return [log(-log(1-x)),  1/(log(1-x)*(x-1))]
    else
        return [log(((1-x)^(-lambda) - 1)/lambda), (lambda * (1-x)^(-1)) / (1 - (1-x)^(lambda))]
    end
end

function DetJ(z, dg)
    n = size(z, 1)
    J = 1
    for i in 1:n
        J = J * abs(dg(z[i]))
    end
    return J
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
export get_guild_audit_log

"""
    get_guild_audit_log(c::Client, guild::Integer; kwargs...) -> AuditLog

Get a [`Guild`](@ref)'s [`AuditLog`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/audit-log#get-guild-audit-log).
"""
get_guild_audit_log(c::Client, guild::Integer; kwargs...) = Response{AuditLog}(c, :GET, "/guilds/$guild/audit-logs"; kwargs...)


text:
export angle,
       approximate_angle,
       approximate_dihedral_angle,
       are_ordered_along_line,
       are_strictly_ordered_along_line,
       area,
       barycenter,
       bisector,
       centroid,
       circumcenter,
       collinear,
       collinear_are_ordered_along_line,
       collinear_are_strictly_ordered_along_line,
       compare_dihedral_angle,
       compare_distance_to_point,
       compare_lexicographically,
       compare_signed_distance_to_line,
       compare_signed_distance_to_plane,
       compare_slope,
       compare_squared_distance,
       compare_squared_radius,
       compare_x,
       compare_x_at_y,
       compare_xy,
       compare_xyz,
       compare_y,
       compare_y_at_x,
       compare_yx,
       coplanar,
       coplanar_orientation,
       coplanar_side_of_bounded_circle,
       cross_product,
       determinant,
       do_intersect,
       do_overlap,
       has_larger_distance_to_point,
       has_larger_signed_distance_to_line,
       has_larger_signed_distance_to_plane,
       has_smaller_distance_to_point,
       has_smaller_signed_distance_to_line,
       has_smaller_signed_distance_to_plane,
       intersection,
       l_infinity_distance,
       left_turn,
       lexicographically_xy_larger,
       lexicographically_xy_larger_or_equal,
       lexicographically_xy_smaller,
       lexicographically_xy_smaller_or_equal,
       lexicographically_xyz_smaller,
       lexicographically_xyz_smaller_or_equal,
       max_vertex,
       midpoint,
       min_vertex,
       normal,
       orientation,
       orthogonal_vector,
       parallel,
       radical_line,
       radical_plane,
       rational_rotation_approximation,
       right_turn,
       scalar_product,
       side_of_bounded_circle,
       side_of_bounded_sphere,
       side_of_oriented_circle,
       side_of_oriented_sphere,
       squared_area,
       squared_distance,
       squared_radius,
       unit_normal,
       volume,
       x_equal,
       y_equal,
       z_equal

"""
    angle(u::Vector23, v::Vector23)

Returns [`OBTUSE`](@ref), [`RIGHT`](@ref) or [`ACUTE`](@ref) depending on the
angle formed by the two vectors `u` and `v`.
"""
angle(u::Vector23, v::Vector23)

"""
    angle(p::Point23, q::Point23, r::Point23)

Returns [`OBTUSE`](@ref), [`RIGHT`](@ref) or [`ACUTE`](@ref) depending on the
angle formed by the three points `p`, `q`, `r` (`q` being the vertex of the
angle).

The returned value is the same as `angle(p - q, r - q)`.
"""
angle(p::Point23, q::Point23, r::Point23)

"""
    angle(p::Point23, q::Point23, r::Point23, s::Point23)

Returns [`OBTUSE`](@ref), [`RIGHT`](@ref) or [`ACUTE`](@ref) depending on the
angle formed by the two vectors `pq`, `rs`.

The returned value is the same as `angle(q - p, s - r)`.
"""
angle(p::Point23, q::Point23, r::Point23, s::Point23)

"""
    angle(p::Point3, q::Point3, r::Point3, v::Vector3)

Returns [`OBTUSE`](@ref), [`RIGHT`](@ref) or [`ACUTE`](@ref) depending on the
angle formed by the normal of the triangle `pqr` and the vector `v`.
"""
angle(p::Point3, q::Point3, r::Point3, v::Vector3)

"""
    approximate_angle(p::Point3, q::Point3, r::Point3)

Returns an approximation of the angle between `p - q` and `r - q`.

The angle is given in degrees.

!!! info "Precondition"

    `p` and `r` are not equal to `q`.
"""
approximate_angle(p::Point3, q::Point3, r::Point3)

"""
    approximate_angle(u::Vector3, v::Vector3)

Returns an approximation of the angle between `u` and `v`.

The angle is given in degrees.

!!! info "Precondition"

    `u` and `v` are not equal to the null vector.
"""
approximate_angle(u::Vector3, v::Vector3)

"""
    approximate_dihedral_angle(p::Point3, q::Point3, r::Point3, s::Point3)

Returns an approximation of the signed dihedral angle in the tetrahedron `pqrs`
of edge `pq`.

The sign is negative if `orientation(p,q,r,s)` is `NEGATIVE` and positive
otherwise.  The angle is given in degrees.

!!! info "Precondition"
    `p`,`q`,`r` and `p`,`q`,`s` are not collinear.
"""
approximate_dihedral_angle(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    area(p::Point2, q::Point2, r::Point2)

Returns the signed area of the triangle defined by the points `p`, `q` and `r`.
"""
area(p::Point2, q::Point2, r::Point2)

"""
    are_ordered_along_line(p::Point23, q::Point23, r::Point23)

Returns `true`, iff the three points are collinear and `q` lies between `p` and
`r`.

Note that `true` is returned, if `q==p` or `q==r`.
"""
are_ordered_along_line(p::Point23, q::Point23, r::Point23)

"""
    are_strictly_ordered_along_line(p::Point23, q::Point23, r::Point23)

Returns `true`, iff the three points are collinear and `q` lies strictly between
`p` and `r`.

Note that `false` is returned, if `q==p` or `q==r`.
"""
are_strictly_ordered_along_line(p::Point23, q::Point23, r::Point23)

"""
    barycenter(p₁::Point23, w₁::Real, p₂::Point23)

Compute the barycenter of the points `p₁` and `p₂` with corresponding weights
`w₁` and `1-w₁`.
"""
barycenter(p₁::Point23, w₁::Real, p₂::Point23)

"""
    barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real)

Compute the barycenter of the points `p₁` and `p₂` with corresponding weights
`w₁` and `w₂`.
"""
barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real)

"""
    barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real, p₃::Point23)

Compute the barycenter of the points `p₁`,`p₂` and `p₃` with corresponding
weights `w₁,w₂` and `1-w₁-w₂`.
"""
barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real, p₃::Point23)

"""
    barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real,
               p₃::Point23, w₃::Real)

Compute the barycenter of the points `p₁`,`p₂` and `p₃` with corresponding
weights `w₁,w₂` and `w₃`.
"""
barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real, p₃::Point23, w₃::Real)

"""
    barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real,
               p₃::Point23, w₃::Real, p₄::Point23)

Compute the barycenter of the points `p₁`,`p₂`, `p₃` and `p₄` with corresponding
weights `w₁,w₂,w₃` and `1-w₁-w₂-w₃`.
"""
barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real, p₃::Point23, w₃::Real,
           p₄::Point23)

"""
    barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real,
               p₃::Point23, w₃::Real, p4::Point23, w4::Real)

Compute the barycenter of the points `p₁`,`p₂`, `p₃` and `p₄` with corresponding
weights `w₁,w₂,w₃` and `w₄`.
"""
barycenter(p₁::Point23, w₁::Real, p₂::Point23, w₂::Real, p₃::Point23, w₃::Real,
           p₄::Point23, w₄::Real)

for P ∈ (:Point2, :Point3)
    @eval begin
        @cxxdereference barycenter(p₁::$P, w₁::Real, p₂::$P) =
            barycenter(p₁, convert(FT, w₁), p₂)
        @cxxdereference barycenter(p₁::$P, w₁::Real, p₂::$P, w₂::Real) =
            barycenter(p₁, convert(FT, w₁), p₂, convert(FT, w₂))
        @cxxdereference barycenter(p₁::$P, w₁::Real, p₂::$P, w₂::Real, p₃::$P) =
            barycenter(p₁, convert(FT, w₁), p₂, convert(FT, w₂), p₃)
        @cxxdereference barycenter(p₁::$P, w₁::Real, p₂::$P, w₂::Real, p₃::$P, w₃::Real) =
            barycenter(p₁, convert(FT, w₁), p₂, convert(FT, w₂), p₃, convert(FT, w₃))
        @cxxdereference barycenter(p₁::$P, w₁::Real, p₂::$P, w₂::Real, p₃::$P, w₃::Real, p₄::$P) =
            barycenter(p₁, convert(FT, w₁), p₂, convert(FT, w₂), p₃, convert(FT, w₃), p₄)
        @cxxdereference barycenter(p₁::$P, w₁::Real, p₂::$P, w₂::Real, p₃::$P, w₃::Real, p₄::$P, w₄::Real) =
            barycenter(p₁, convert(FT, w₁), p₂, convert(FT, w₂), p₃, convert(FT, w₃), p₄, convert(FT, w₄))
    end
end

"""
    bisector(p::Point2, q::Point2)

Constructs the bisector line of two points `p` and `q`.

The bisector is oriented in such a way that `p` lies on its positive side.

!!! info "Precondition"

    `p` and `q` are not equal.
"""
bisector(p::Point2, q::Point2)

"""
    bisector(l₁::Line2, l₂::Line2)

Constructs the bisector of the two lines `l₁` and `l₂`.

In the general case, the bisector has the direction of the vector which is the
sum of the normalized directions of the two lines, and which passes through the
intersection of `l₁` and `l₂`.  If `l₁` and `l₂` are parallel, then the bisector
is defined as the line which has the same direction as `l₁`, and which is at the
same distance from `l₁` and `l₂`.  This function requires that [`RT`](@ref)
supports the [`sqrt()`](@ref) operation.
"""
bisector(l₁::Line2, l₂::Line2)

"""
    bisector(p::Point3, q::Point3)

Constructs the bisector plane of the two points `p` and `q`.

The bisector is oriented in such a way that `p` lies on its positive side.

!!! info "Precondition"

    `p` and `q` are not equal.
"""
bisector(p::Point3, q::Point3)

"""
    bisector(h₁::Plane3, h₂::Plane3)

Constructs the bisector of the two planes `h₁` and `h₂`.

In the general case, the bisector has a normal vector which has the same
direction as the sum of the normalized normal vectors of the two planes, and
passes through the intersection of `h₁` and `h₂`.  If `h₁` and `h₂` are parallel,
then the bisector is defined as the plane which has the same oriented normal
vector as `l₁`, and which is at the same distance from `h₁` and `h₂`.  This
function requires that [`RT`](@ref) supports the [`sqrt()`](@ref) operation.
"""
bisector(h₁::Plane3, h₂::Plane3)

"""
    centroid(p::Point23, q::Point23, r::Point23)

Compute the centroid of the points `p`, `q`, and `r`.
"""
centroid(p::Point23, q::Point23, r::Point23)

"""
    centroid(p::Point23, q::Point23, r::Point23, s::Point23)

Compute the centroid of the points `p`, `q`, `r`, and `s`.
"""
centroid(p::Point23, q::Point23, r::Point23, s::Point23)

"""
    centroid(t::Triangle23)

Compute the centroid of the triangle `t`.
"""
centroid(t::Triangle23)

"""
    centroid(t::Tetrahedron3)

Compute the centroid of the tetrahedron `t`.
"""
centroid(t::Tetrahedron3)

"""
    circumcenter(p::Point23, q::Point23)

Compute the center of the smallest circle passing through the points `p` and
`q`.

Note: this is the same as `midpoint(p, q)` but is provided for homogeneity.
"""
circumcenter(p::Point23, q::Point23)

"""
    circumcenter(p::Point23, q::Point23, s::Point23)

Compute the center of the circle passing through the points `p`, `q`, and `r`.

!!! info "Precondition"

    `p`, `q`, and `r` are not collinear.
"""
circumcenter(p::Point23, q::Point23, r::Point23)

"""
    circumcenter(t::Triangle23)

Compute the center of the circle passing through the vertices of `t`.

!!! info "Precondition"

    `t` is not degenerate.
"""
circumcenter(t::Triangle23)

"""
    circumcenter(t::Tetrahedron3)

Compute the center of the sphere passing through the vertices of `t`.
"""
circumcenter(t::Tetrahedron3)

"""
    collinear_are_ordered_along_line(p::Point23, q::Point23, r::Point23)

Returns `true`, iff `q` lies between `p` and `r`.

!!! info "Precondition"

    `p`, `q` and `r` are collinear.
"""
collinear_are_ordered_along_line(p::Point23, q::Point23)

"""
    collinear_are_strictly_ordered_along_line(p::Point23, q::Point23, r::Point23)

Returns `true`, iff `q` lies strictly between `p` and `r`.

!!! info "Precondition"

    `p`, `q` and `r` are collinear.
"""
collinear_are_strictly_ordered_along_line(p::Point23, q::Point23)

"""
    collinear(p::Point23, q::Point23, r::Point23)

Returns `true`, iff `p`, `q`, and `r` are collinear.
"""
collinear(p::Point23, q::Point23, r::Point23)

@doc raw"""
    compare_dihedral_angle(a₁::Point3, b₁::Point3, c₁::Point3, d₁::Point3, cosine::Real)

Compares the dihedral angles ``θ_1`` and ``θ_2``, where ``θ_1`` is the dihedral
angle, in ``[0, π]``, of the tetrahedron `(a₁, b₁, c₁, d₁)` at the edge `(a₁,
b₁)`, and ``θ_2`` is the angle, in ``[0, π]``, such that ``cos(θ_2) = cosine``.

!!! info "Precondition"

    `a₁`, `b₁`, `c₁` are not collinear, and `a₁`, `b₁`, `d₁` are not collinear.
"""
compare_dihedral_angle(a₁::Point3, b₁::Point3, c₁::Point3, d₁::Point3, cosine::Real)

@cxxdereference compare_dihedral_angle(a₁::Point3, b₁::Point3, c₁::Point3, d₁::Point3, cosine::Real) =
    compare_dihedral_angle(a₁, b₁, c₁, d₁, convert(FT, cosine))

@doc raw"""
    compare_dihedral_angle(a₁::Point3, b₁::Point3, c₁::Point3, d₁::Point3,
                           a₂::Point3, b₂::Point3, c₂::Point3, d₂::Point3)

Compares the dihedral angles ``θ_1`` and ``θ_2``, where ``θ_i`` is the dihedral
angle in the tetrahedron `(aᵢ, bᵢ, cᵢ, dᵢ)` at the edge `(aᵢ, bᵢ)`.

These two angles are computed in ``[0, π]``.  The result is the same as
`compare_dihedral_angle(b₁-a₁, c₁-a₁, d₁-a₁, b₂-a₂, c₂-a₂, d₂-a₂)`.

!!! info "Precondition"

    For ``i ∈ {1, 2}``, `aᵢ`, `bᵢ`, `cᵢ` are not collinear, and `aᵢ`,
    `bᵢ`, `dᵢ` are not collinear.
"""
compare_dihedral_angle(a₁::Point3, b₁::Point3, c₁::Point3, d₁::Point3,
                       a₂::Point3, b₂::Point3, c₂::Point3, d₂::Point3)

@doc raw"""
    compare_dihedral_angle(u₁::Vector3, v₁::Vector3, w₁::Vector3, cosine::Real)

Compares the dihedral angles ``θ_1`` and ``θ_2``, where ``θ_1`` is the dihedral
angle, in ``[0, π]``, between the vectorial planes defined by `(u₁, v₁)` and
`(u₁, w₁)`, and ``θ_2`` is the angle, in ``[0, π]``, such that ``cos(θ_2) =
cosine`.

!!! info "Precondition"

    `u₁` and `v₁` are not collinear, and `u₁` and `w₁` are not collinear.
"""
compare_dihedral_angle(u₁::Vector3, v₁::Vector3, w₁::Vector3, cosine::Real)

@cxxdereference compare_dihedral_angle(u₁::Vector3, v₁::Vector3, w₁::Vector3, cosine) =
    compare_dihedral_angle(u₁, v₁, w₁, convert(FT, cosine))

@doc raw"""
    compare_dihedral_angle(u₁::Vector3, v₁::Vector3, w₁::Vector3,
                           u₂::Vector3, v₂::Vector3, w₂::Vector3)

Compares the dihedral angles ``θ_1`` and ``θ_2``, where ``θ_i`` is the dihedral
angle between the vectorial planes defined by `(uᵢ, vᵢ)` and `(uᵢ, wᵢ)`.

These two angles are computed in ``[0, π]``.

!!! info "Precondition"

    For ``i ∈ {1, 2}``, `uᵢ` and `vᵢ` are not collinear, and `uᵢ` and `wᵢ` are
    not collinear.
"""
compare_dihedral_angle(u₁::Vector3, v₁::Vector3, w₁::Vector3,
                       u₂::Vector3, v₂::Vector3, w₂::Vector3)

"""
    compare_distance_to_point(p::Point23, q::Point23, r::Point23)

Compares the distances of points `q` and `r` to `p`.

Returns [`SMALLER`](@ref), iff `q` is closer to `p` than `r`, [`LARGER`](@ref)
iff `r` is closer to `p` than `q`, and [`EQUAL`](@ref) otherwise.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point`()](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
compare_distance_to_point(p::Point23, q::Point23, r::Point23)

"""
    compare_lexicographically(p::Point2, q::Point2)

Compares the Cartesian coordinates of points `p` and `q` lexicographically in
``xy`` order: first ``x``-coordinates are compared, if they are equal,
``y``-coordinates are compared.

This is the same function as `compare_xy`.
"""
compare_lexicographically(p::Point2, q::Point2)

"""
    compare_lexicographically(p::Point3, q::Point3)

Compares the Cartesian coordinates of points `p` and `q` lexicographically in
``xyz`` order: first ``x``-coordinates are compared, if they are equal,
``y``-coordinates are compared, and if both ``x``- and ``y``- coordinate are
equal, ``z``-coordinates are compared.

This is the same function as `compare_xyz`.
"""
compare_lexicographically(p::Point3, q::Point3)

"""
    compare_signed_distance_to_line(l::Line2, p::Point2, q::Point2)

Returns [`LARGER`](@ref) iff the signed distance of `p` and `l` is larger than
the signed distance of `q` and `l`, [`SMALLER`](@ref), iff it is smaller, and
[`EQUAL`](@ref) iff both are equal.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
compare_signed_distance_to_line(l::Line2, p::Point2, q::Point2)

"""
    compare_signed_distance_to_line(p::Point2, q::Point2, r::Point2, s::Point2)

Returns [`LARGER`](@ref) iff the signed distance of `r` and `l` is larger than
the signed distance of `s` and `l`, [`SMALLER`](@ref), iff it is smaller, and
[`EQUAL`](@ref) iff both are equal, where `l` is the directed line through `p`
and `q`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
compare_signed_distance_to_line(p::Point2, q::Point2, r::Point2, s::Point2)

"""
    compare_signed_distance_to_plane(h::Plane3, p::Point3, q::Point3)

Returns [`LARGER`](@ref) iff the signed distance of `p` and `h` is larger than
the signed distance of `q` and `h`, [`SMALLER`](@ref), iff it is smaller, and
[`EQUAL`](@ref) iff both are equal.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
compare_signed_distance_to_plane(h::Plane3, p::Point3, q::Point3)

"""
    compare_signed_distance_to_plane(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

Returns [`LARGER`](@ref) iff the signed distance of `s` and `h` is larger than
the signed distance of `t` and `h`, [`SMALLER`](@ref), iff it is smaller, and
[`EQUAL`](@ref) iff both are equal, where `h` is the oriented plane through `p`,
`q` and `r`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
compare_signed_distance_to_plane(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

"""
    compare_slope(l₁::Line2, l₂::Line2)

Compares the slopes of the lines `l₁` and `l₂`.
"""
compare_slope(l₁::Line2, l₂::Line2)

"""
    compare_slope(s₁::Segment2, s₂::Segment2)

Compares the slopes of the segments `s₁` and `s₂`, where the slope is the
variation of the ``y``-coordinate from the left to the right endpoint of the
segments.
"""
compare_slope(s₁::Segment2, s₂::Segment2)

"""
    compare_slope(p::Point3, q::Point3, r::Point3, s::Point3)

Compares the slopes of the segments `(p, q)` and `(r, s)`, where the slope is
the variation of the ``z``-coordinate from the first to the second Point23 of the
segment divided by the length of the segment.
"""
compare_slope(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    compare_squared_distance(p::Point23, q::Point23, d²::Real)

Compares the squared distance of points `p` and `q` to `d²`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
compare_squared_distance(p::Point23, q::Point23, d²::Real)

compare_squared_distance(p, q, d²::Real) =
    compare_squared_distance(p, q, convert(FT, d²))

"""
    compare_squared_radius(p::Point3, r²::Real)

Compares the squared radius of the sphere of radius 0 center at `p` to `r²`.

This returns the opposite sign of `r²`.
"""
compare_squared_radius(p::Point3, r²::Real)

"""
    compare_squared_radius(p::Point3, q::Point3, r²::Real)

Compares the squared radius of the sphere defined by the points `p` and `q` to
`r²`.
"""
compare_squared_radius(p::Point3, q::Point3, r²::Real)

"""
    compare_squared_radius(p::Point3, q::Point3, r::Point3, r²::Real)

Compares the squared radius of the sphere defined by the points `p`, `q`, and
`r` to `r²`.
"""
compare_squared_radius(p::Point3, q::Point3, r::Point3, r²::Real)

"""
    compare_squared_radius(p::Point3, q::Point3, r::Point3, s::Point3, r²::Real)

Compares the squared radius of the sphere defined by the points `p`, `q`, `r`,
and `s` to `r²`.
"""
compare_squared_radius(p::Point3, q::Point3, r::Point3, s::Point3, r²::Real)

@cxxdereference compare_squared_radius(p::Point3, r²::Real) =
    compare_squared_radius(p, convert(FT, r²))
@cxxdereference compare_squared_radius(p::Point3, q::Point3, r²::Real) =
    compare_squared_radius(p, q, convert(FT, r²))
@cxxdereference compare_squared_radius(p::Point3, q::Point3, r::Point3, r²::Real) =
    compare_squared_radius(p, q, r, convert(FT, r²))
@cxxdereference compare_squared_radius(p::Point3, q::Point3, r::Point3, s::Point3, r²::Real) =
    compare_squared_radius(p, q, r, s, convert(FT, r²))

"""
    compare_x(p::Point23, q::Point23)

Compares the ``x``-coordinates of `p` and `q`.

See also: [`compare_xy()`](@ref), [`compare_xyz()`](@ref),
[`compare_x_at_y()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x(p::Point23, q::Point23)

"""
    compare_x(p::Point2, l₁::Line2, l₂::Line2)

Compares the ``x``-coordinates of `p` and the intersection of lines `l₁` and
`l₂`.

See also: [`compare_xy()`](@ref), [`compare_xyz()`](@ref),
[`compare_x_at_y()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x(p::Point2, l₁::Line2, l₂::Line2)

"""
    compare_x(l::Line2, h₁::Line2, h₂::Line2)

Compares the ``x``-coordinates of the intersection of line `l` with of line `h₁`
and with line `h₂`.

See also: [`compare_xy()`](@ref), [`compare_xyz()`](@ref),
[`compare_x_at_y()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x(l::Line2, h₁::Line2, h₂::Line2)

"""
    compare_x(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

Compares the ``x``-coordinates of the intersection of lines `l₁` and `l₂` and
the intersection of lines `h₁` and `h₂`.

See also: [`compare_xy()`](@ref), [`compare_xyz()`](@ref),
[`compare_x_at_y()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

"""
    compare_xy(p::Point23, q::Point23)

Compares the Cartesian coordinates of points `p` and `q` lexicographically in
``xy`` order: first ``x``-coordinates are compared, if they are equal,
``y``-coordinates are compared.

See also: [`compare_x()`](@ref), [`compare_xyz()`](@ref),
[`compare_x_at_y()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_xy(p::Point23, q::Point23)

"""
    compare_x_at_y(p::Point2, h::Line2)

Compares the ``x``-coordinates of `p` and the horizontal projection of `p` on
`h`.

!!! info "Precondition"

    `h` is not horizontal.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x_at_y(p::Point2, h::Line2)

"""
    compare_x_at_y(p::Point2, h₁::Line2, h₂::Line2)

This function compares the ``x``-coordinates of the horizontal projection of `p`
on `h₁` and on `h₂`.

!!! info "Precondition"

    `h₁` and `h₂` are not horizontal.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x_at_y(p::Point2, h₁::Line2, h₂::Line2)

"""
    compare_x_at_y(l₁::Line2, l₂::Line2, h::Line2)

Let `p` be the intersection of lines `l₁` and `l₂`.

This function compares the ``x``-coordinates of `p` and the horizontal
projection of `p` on `h`.

!!! info "Precondition"

    `l₁` and `l₂` intersect and are not horizontal; `h` is not horizontal.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x_at_y(l₁::Line2, l₂::Line2, h::Line2)

"""
    compare_x_at_y(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

Let `p` be the intersection of lines `l₁` and `l₂`.

This function compares the ``x`-coordinates of the horizontal projects of `p`
on `h₁` and on `h₂`.

!!! info "Precondition"

    `l₁` and `l₂` intersect and are not horizontal; `h₁` and `h₂` are not
    horizontal.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_x_at_y(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

"""
    compare_y_at_x(p::Point2, h::Line2)

Compares the ``y``-coordinates of `p` and the vertical projection of `p` on `h`.

!!! info "Precondition"

    `h` is not vertical.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_z()`](@ref).
"""
compare_y_at_x(p::Point2, h::Line2)

"""
    compare_x_at_y(p::Point2, h₁::Line2, h₂::Line2)

This function compares the ``y``-coordinates of the vertical projection of `p`
on `h₁` and on `h₂`.

!!! info "Precondition"

    `h₁` and `h₂` are not vertical.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_z()`](@ref).
"""
compare_y_at_x(p::Point2, h₁::Line2, h₂::Line2)

"""
    compare_y_at_x(l₁::Line2, l₂::Line2, h::Line2)

Let `p` be the intersection of lines `l₁` and `l₂`.

This function compares the ``y``-coordinates of `p` and the vertical projection
of `p` on `h`.

!!! info "Precondition"

    `l₁`, `l₂` intersect and `h` is not vertical.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_z()`](@ref).
"""
compare_y_at_x(l₁::Line2, l₂::Line2, h::Line2)

"""
    compare_y_at_x(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

Let `p` be the intersection of lines `l₁` and `l₂`.

This function compares the ``y`-coordinates of the vertical projects of `p` on
`h₁` and on `h₂`.

!!! info "Precondition"

    `l₁` and `l₂` intersect; `h₁` and `h₂` are not vertical.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_z()`](@ref).
"""
compare_y_at_x(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

"""
    compare_y_at_x(p::Point2, s::Segment2)

Compares the ``y``-coordinates of `p` and the vertical projection of `p` on `s`.

If `s` is vertical, then return [`EQUAL`](@ref) when `p` lies on `s`,
[`SMALLER`](@ref) when `p` lies under `s`, and [`LARGER`](@ref) otherwise.

!!! info "Precondition"

    `p` is withing the x range of `s`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_z()`](@ref).
"""
compare_y_at_x(p::Point2, s::Segment2)

"""
    compare_y_at_x(p::Point2, s₁::Segment2, s₂::Segment2)

Compares the ``y``-coordinates of `p` and the vertical projection of `p` on
`s₁` and on `s₂`.

If `s₁` or `s₂` is vertical, then return [`EQUAL`](@ref) if they intersect,
otherwise return [`SMALLER`](@ref) if `s₁` lies below `s₂`, and return
[`LARGER`](@ref) otherwise.

!!! info "Precondition"

    `p` is withing the x range of `s₁` and `s₂`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_z()`](@ref).
"""
compare_y_at_x(p::Point2, s₁::Segment2, s₂::Segment2)

"""
    compare_y(p::Point23, q::Point23)

Compares the ``y``-coordinates of `p` and `q`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_y(p::Point23, q::Point23)

"""
    compare_y(p::Point2, l₁::Line2, l₂::Line2)

Compares the ``y``-coordinates of `p` and the intersection of lines `l₁` and
`l₂`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_y(p::Point2, l₁::Line2, l₂::Line2)

"""
    compare_y(l::Line2, h₁::Line2, h₂::Line2)

Compares the ``y``-coordinates of the intersection of line `l` with of line `h₁`
and with line `h₂`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_y(l::Line2, h₁::Line2, h₂::Line2)

"""
    compare_y(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

Compares the ``y``-coordinates of the intersection of lines `l₁` and `l₂` and
the intersection of lines `h₁` and `h₂`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_y(l₁::Line2, l₂::Line2, h₁::Line2, h₂::Line2)

"""
    compare_xyz(p::Point3, q::Point3)

Compares the Cartesian coordinates of points `p` and `q` lexicographically in
``xyz`` order: first ``x``-coordinates are compared, if they are equal,
``y``-coordinates are compared, and if both ``x``- and ``y``- coordinates are
equal, ``z``-coordinates are compared.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_x_at_y()`](@ref), [`compare_y()`](@ref), [`compare_yx()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_xyz(p::Point3, q::Point3)

"""
    compare_z(p::Point3, q::Point3)

Compares the ``z``-coordinates of `p` and `q`.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_yx()`](@ref), [`compare_y_at_x()`](@ref).
"""
compare_z(p::Point3, q::Point3)

"""
    compare_yx(p::Point2, q::Point2)

Compares the Cartesian coordinates of points `p` and `q` lexicographically in
``yx`` order: first ``y``-coordinates are compared, if they are equal,
``x``-coordinates are compared.

See also: [`compare_x()`](@ref), [`compare_xy()`](@ref),
[`compare_xyz()`](@ref), [`compare_x_at_y()`](@ref), [`compare_y()`](@ref),
[`compare_y_at_x()`](@ref), [`compare_z()`](@ref).
"""
compare_yx(p::Point2, q::Point2)

"""
    coplanar(p::Point3, q::Point3, r::Point3, s::Point3)

Returns `true`, if `p`, `q`, `r`, and `s` are coplanar.
"""
coplanar(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    coplanar_orientation(p::Point3, q::Point3, r::Point3, s::Point3)

Let `P` be the plane defined by the points `p`, `q`, and `r`.

Note that the order defined the orientation of `P`.  The function computes the
orientation of points `p`, `q`, and `s` in `P`: Iff `p`, `q`, `s` are collinear,
[`COLLINEAR`](@ref) is returned.  Iff `P` and the plane defined by `p`, `q`, and
`s` have the same orientation, [`POSITIVE`](@ref) is returned; otherwise
[`NEGATIVE`](@ref) is returned.

!!! info "Precondition"

    `p`, `q`, `r`, and `s` are coplanar and `p`, `q`, and `r` are not collinear.
"""
coplanar_orientation(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    coplanar_orientation(p::Point3, q::Point3, r::Point3)

If `p`, `q`, `r` are collinear, then [`COLLINEAR`](@ref) is returned.

If not, then `p`, `q`, `r` define a plane `P`.  The return value, in this case,
is either [`POSITIVE`](@ref) or [`NEGATIVE`](@ref), but we don't specify it
explicitly.  However, we guarantee that all calls to this predicate over 3 points
in `P` will return a coherent orientation if considered a 2D orientation in `P`.
"""
coplanar_orientation(p::Point3, q::Point3, r::Point3)

"""
    coplanar_side_of_bounded_circle(p::Point3, q::Point3, r::Point3, s::Point3)

Returns the bounded side of the circle defined by `p`, `q`, and `r` on which `s`
lies.

!!! info "Precondition"

    `p`, `q`, `r`, and `s` are coplanar and `p`, `q`, and `r` are not collinear.
"""
coplanar_side_of_bounded_circle(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    cross_product(v::Vector3, u::Vector3)

Returns the cross product of `u` and `v`.
"""
cross_product(v::Vector3, u::Vector3)

"""
    determinant(v::Vector2, w::Vector2)

Returns the determinant of `v` and `w`.
"""
determinant(v::Vector2, w::Vector2)

"""
    determinant(u::Vector3, v::Vector3, w::Vector3)

Returns the determinant of `u`, `v`, and `w`.
"""
determinant(u::Vector3, v::Vector3, w::Vector3)

"""
    do_intersect(obj₁::Type1, obj₂::Type2)

Checks whether `obj₁` and `obj₂` intersect.

Two objects `obj₁` and `obj₂` intersect if there is a Point23 `p` that is
part of both `obj₁` and `obj₂`.  The intersection region of those two objects
is defined as the set of all points `p` that are part of both `obj₁` and
`obj₂`.  Note that for objects like triangles and polygons that enclose a
bounded region, this region is part of the object.

The types `Type1` and `Type2` can be any of the following:

- [`Point2`](@ref)
- [`Line2`](@ref)
- [`Ray2`](@ref)
- [`Segment2`](@ref)
- [`Triangle2`](@ref)
- [`IsoRectangle2`](@ref)

Also, `Type1` and `Type2` can be respectively of types

- [`Circle2`](@ref) and [`CircularArc2`](@ref) (or the contrary)
- [`Circle2`](@ref) and [`Line2`](@ref) (or the contrary)
- [`Circle2`](@ref) and [`Segment2`](@ref) (or the contrary)
- [`CircularArc2`](@ref) and [`Line2`](@ref) (or the contrary)
- [`CircularArc2`](@ref) and [`Segment2`](@ref) (or the contrary)

In three-dimensional space, the types `Type1` and `Type2` can be any of the
following:

- [`Bbox3`](@ref)
- [`IsoCuboid3`](@ref)
- [`Line3`](@ref)
- [`Plane3`](@ref)
- [`Point3`](@ref)
- [`Ray3`](@ref)
- [`Segment3`](@ref)
- [`Sphere3`](@ref)
- [`Tetrahedron3`](@ref)
- [`Triangle3`](@ref)

Also, `Type1` and `Type2` can be both of type

- [`Circle3`](@ref)
- [`CircularArc3`](@ref)

Furthermore, `Type1` and `Type2` can be interchangeably of types

- [`Circle3`](@ref) and [`Line3`](@ref)
- [`Circle3`](@ref) and [`Plane3`](@ref)
- [`Circle3`](@ref) and [`Sphere3`](@ref)
- [`CircularArc3`](@ref) and [`Plane3`](@ref)

See also: [`intersection`](@ref).
"""
do_intersect(obj₁, obj₂)

"""
    do_intersect(obj₁::Type1, obj₂::Type2, obj₃::Type3)

Checks whether `obj₁`, `obj₂`, and `obj₃` intersect.

The types `Type1`, `Type2` and `Type3` can be respectively of types

- [`Sphere3`](@ref), [`Sphere3`](@ref), and [`Plane3`](@ref)
- [`Plane3`](@ref), [`Sphere3`](@ref), and [`Sphere3`](@ref)
- [`Plane3`](@ref), [`Plane3`](@ref), and [`Sphere3`](@ref)
- [`Sphere3`](@ref), [`Plane3`](@ref), and [`Plane3`](@ref)
"""
do_intersect(obj₁, obj₂, obj₃)

"""
    has_larger_distance_to_point(p::Point23, q::Point23, r::Point23)

Returns `true` iff the distance between `q` and `p` is larger than the distance
between `r` and `p`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_larger_distance_to_point(p::Point23, q::Point23, r::Point23)

"""
    has_larger_signed_distance_to_line(l::Line2, p::Point2, q::Point2)

Returns `true` iff the signed distance of `p` and `l` is larger than the signed
distance of `q` and `l`

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_larger_signed_distance_to_line(l::Line2, p::Point2, q::Point2)

"""
    has_larger_signed_distance_to_line(p::Point2, q::Point2, r::Point2, s::Point2)

Returns `true` iff the signed distance of `r` and `l` is larger than the signed
distance of `s` and `l`, where `l` is the directed line through points `p` and
`q`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_larger_signed_distance_to_line(p::Point2, q::Point2, r::Point2, s::Point2)

"""
    has_larger_signed_distance_to_plane(h::Plane3, p::Point3, q::Point3)

Returns `true` iff the signed distance of `p` and `h` is larger than
the signed distance of `q` and `h`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_larger_signed_distance_to_plane(h::Plane3, p::Point3, q::Point3)

"""
    has_larger_signed_distance_to_plane(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

Returns `true` iff the signed distance of `s` and `h` is larger than the signed
distance of `t` and `h`, where `h` is the oriented plane through `p`, `q`, and
`r`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_larger_signed_distance_to_plane(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

"""
    has_smaller_distance_to_point(p::Point23, q::Point23, r::Point23)

Returns `true` iff the distance between `q` and `p` is smaller than the distance
between `r` and `p`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_smaller_distance_to_point(p::Point23, q::Point23, r::Point23)

"""
    has_smaller_signed_distance_to_line(l::Line2, p::Point2, q::Point2)

Returns `true` iff the signed distance of `p` and `l` is larger than the signed
distance of `q` and `l`

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_smaller_signed_distance_to_line(l::Line2, p::Point2, q::Point2)

"""
    has_smaller_signed_distance_to_line(p::Point2, q::Point2, r::Point2, s::Point2)

Returns `true` iff the signed distance of `r` and `l` is smaller than the signed
distance of `s` and `l`, where `l` is the directed line through points `p` and
`q`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_smaller_signed_distance_to_line(p::Point2, q::Point2, r::Point2, s::Point2)

"""
    has_smaller_signed_distance_to_plane(h::Plane3, p::Point3, q::Point3)

Returns `true` iff the signed distance of `p` and `h` is smaller than
the signed distance of `q` and `h`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_smaller_signed_distance_to_plane(h::Plane3, p::Point3, q::Point3)

"""
    has_smaller_signed_distance_to_plane(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

Returns `true` iff the signed distance of `s` and `h` is smaller than the signed
distance of `t` and `h`, where `h` is the oriented plane through `p`, `q`, and
`r`.

See also: [`compare_squared_distance()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
has_smaller_signed_distance_to_plane(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

"""
    intersection(obj₁::Type1, obj₂::Type2)

Two objects `obj₁` and `obj₂` intersect if there is a Point23 `p` that is
part of both `obj₁` and `obj₂`.

The intersection region of those two objects is defined as the set of all points
`p` that are part of both `obj₁` and `obj₂`.  Note that for objects like
triangles and polygons that enclose a bounded region, this region is considered
part of the object.  If a segment lies completely inside a triangle, then those
two objects intersect and the intersection region is the complete segment.

The return type is `Union{T...,Nothing}`, the last column in the table providing
the template parameter pack.

The following tables give the possible values for `Type1` and `Type2`.

### 2D Intersections

|          Type1          |          Type2          |                                       Return Type: T...                                      |
|:------------------------|:------------------------|:---------------------------------------------------------------------------------------------|
| [`Circle2`](@ref)       | [`Circle2`](@ref)       | [`Point2`](@ref), [`Vector{Point2}`](@ref), or [`Circle2`](@ref)                             |
| [`Circle2`](@ref)       | [`CircularArc2`](@ref)  | [`Point2`](@ref), [`Vector{Point2}`](@ref), or [`CircularArc2`](@ref)                        |
| [`Circle2`](@ref)       | [`Line2`](@ref)         | [`Point2`](@ref), or [`Vector{Point2}`](@ref)                                                |
| [`Circle2`](@ref)       | [`Segment2`](@ref)      | [`Point2`](@ref), or [`Vector{Point2}`](@ref)                                                |
| [`CircularArc2`](@ref)  | [`CircularArc2`](@ref)  | [`Point2`](@ref), [`Vector{Point2}`](@ref), or [`CircularArc2`](@ref)                        |
| [`CircularArc2`](@ref)  | [`Line2`](@ref)         | [`Point2`](@ref), or [`Vector{Point2}`](@ref)                                                |
| [`CircularArc2`](@ref)  | [`Segment2`](@ref)      | [`Point2`](@ref), or [`Vector{Point2}`](@ref)                                                |
| [`IsoRectangle2`](@ref) | [`IsoRectangle2`](@ref) | [`IsoRectangle2`](@ref)                                                                      |
| [`IsoRectangle2`](@ref) | [`Line2`](@ref)         | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`IsoRectangle2`](@ref) | [`Ray2`](@ref)          | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`IsoRectangle2`](@ref) | [`Segment2`](@ref)      | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`IsoRectangle2`](@ref) | [`Triangle2`](@ref)     | [`Point2`](@ref), or [`Segment2`](@ref), or [`Triangle2`](@ref), or [`Vector{Point2}`](@ref) |
| [`Line2`](@ref)         | [`Line2`](@ref)         | [`Point2`](@ref), or [`Line2`](@ref)                                                         |
| [`Line2`](@ref)         | [`Ray2`](@ref)          | [`Point2`](@ref), or [`Ray2`](@ref)                                                          |
| [`Line2`](@ref)         | [`Segment2`](@ref)      | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`Line2`](@ref)         | [`Triangle2`](@ref)     | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`Ray2`](@ref)          | [`Ray2`](@ref)          | [`Point2`](@ref), or [`Segment2`](@ref), or [`Ray2`](@ref)                                   |
| [`Ray2`](@ref)          | [`Segment2`](@ref)      | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`Ray2`](@ref)          | [`Triangle2`](@ref)     | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`Segment2`](@ref)      | [`Segment2`](@ref)      | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`Segment2`](@ref)      | [`Triangle2`](@ref)     | [`Point2`](@ref), or [`Segment2`](@ref)                                                      |
| [`Triangle2`](@ref)     | [`Triangle2`](@ref)     | [`Point2`](@ref), or [`Segment2`](@ref), or [`Triangle2`](@ref), or [`Vector{Point2}`](@ref) |

**3D Intersections**

|         Type1          |         Type2          |                                      Return Type: T...                                       |
|:-----------------------|:-----------------------|:---------------------------------------------------------------------------------------------|
| [`Circle3`](@ref)      | [`Circle3`](@ref)      | [`Point3`](@ref), [`Circle3`](@ref), or [`Vector{Point3}`](@ref)                             |
| [`Circle3`](@ref)      | [`Line3`](@ref)        | [`Point3`](@ref), [`Vector{Point3}`](@ref)                                                   |
| [`Circle3`](@ref)      | [`Plane3`](@ref)       | [`Point3`](@ref), [`Circle3`](@ref), or [`Vector{Point3}`](@ref)                             |
| [`Circle3`](@ref)      | [`Sphere3`](@ref)      | [`Point3`](@ref), [`Circle3`](@ref), or [`Vector{Point3}`](@ref)                             |
| [`CircularArc3`](@ref) | [`CircularArc3`](@ref) | [`Point3`](@ref), [`Circle3`](@ref), or [`CircularArc3`](@ref), or [`Vector{Point3}`](@ref)  |
| [`CircularArc3`](@ref) | [`Plane3`](@ref)       | [`Point3`](@ref), [`CircularArc3`](@ref), or [`Vector{Point3}`](@ref)                        |
| [`Line3`](@ref)        | [`Line3`](@ref)        | [`Point3`](@ref), or [`Line3`](@ref)                                                         |
| [`Line3`](@ref)        | [`Line3`](@ref)        | [`Point3`](@ref), or [`Line3`](@ref)                                                         |
| [`Line3`](@ref)        | [`Plane3`](@ref)       | [`Point3`](@ref), or [`Line3`](@ref)                                                         |
| [`Line3`](@ref)        | [`Ray3`](@ref)         | [`Point3`](@ref), or [`Ray3`](@ref)                                                          |
| [`Line3`](@ref)        | [`Segment3`](@ref)     | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Line3`](@ref)        | [`Triangle3`](@ref)    | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Plane3`](@ref)       | [`Plane3`](@ref)       | [`Line3`](@ref), or [`Segment3`](@ref)                                                       |
| [`Plane3`](@ref)       | [`Ray3`](@ref)         | [`Point3`](@ref), or [`Ray3`](@ref)                                                          |
| [`Plane3`](@ref)       | [`Segment3`](@ref)     | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Plane3`](@ref)       | [`Sphere3`](@ref)      | [`Point3`](@ref), or [`Circle3`](@ref)                                                       |
| [`Plane3`](@ref)       | [`Triangle3`](@ref)    | [`Point3`](@ref), or [`Segment3`](@ref), or [`Triangle3`](@ref)                              |
| [`Ray3`](@ref)         | [`Ray3`](@ref)         | [`Point3`](@ref), or [`Ray3`](@ref), or [`Segment3`](@ref)                                   |
| [`Ray3`](@ref)         | [`Segment3`](@ref)     | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Ray3`](@ref)         | [`Triangle3`](@ref)    | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Segment3`](@ref)     | [`Segment3`](@ref)     | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Segment3`](@ref)     | [`Triangle3`](@ref)    | [`Point3`](@ref), or [`Segment3`](@ref)                                                      |
| [`Segment3`](@ref)     | [`Sphere3`](@ref)      | [`Point3`](@ref), or [`Circle3`](@ref), or [`Sphere3`](@ref)                                 |
| [`Triangle3`](@ref)    | [`Triangle3`](@ref)    | [`Point3`](@ref), or [`Segment3`](@ref), or [`Triangle3`](@ref), or [`Vector{Point3}`](@ref) |

# Examples

The following examples demonstrate the most common use of
[`intersection()`](@ref) functions.

In the first two examples we intersect a segment and a line.

```julia-repl
julia> seg = Segment2(Point2(0,0), Point2(2,2))
Segment_2(PointC2(0, 0), PointC2(2, 2))

julia> lin = Line2(1,-1,0)
Line_2(0, 0, 0)

julia> result = intersection(seg, lin)
Segment_2(PointC2(0, 0), PointC2(2, 2))
```
"""
intersection(obj₁, obj₂)

"""
    intersection(obj₁::Type1, obj₂::Type2, obj₃::Type3)

Computes the intersection of the three given objects.

The return type is `Union{Circle3,Plane3,Sphere3,Vector{Point3},Nothing}`.

The following table gives the possible types for `Type1`, `Type2`, and `Type3`.

|         Type1          |         Type2          |         Type3          |
|:-----------------------|:-----------------------|:-----------------------|
| [`Sphere3`](@ref)      | [`Sphere3`](@ref)      | [`Sphere3`](@ref)      |
| [`Sphere3`](@ref)      | [`Sphere3`](@ref)      | [`Plane3`](@ref)       |
| [`Plane3`](@ref)       | [`Sphere3`](@ref)      | [`Sphere3`](@ref)      |
| [`Plane3`](@ref)       | [`Plane3`](@ref)       | [`Sphere3`](@ref)      |
| [`Sphere3`](@ref)      | [`Plane3`](@ref)       | [`Plane3`](@ref)       |
"""
intersection(obj₁, obj₂, obj₃)

"""
    l_infinity_distance(p::Point23, q::Point23)

Returns the distance between `p` and `q` in the L-infinity metric.
"""
l_infinity_distance(p::Point23, q::Point23)

"""
    left_turn(p::Point2, q::Point2, r::Point2)

Returns `true` iff `p`, `q`, and `r` form a left turn.

See also: [`collinear()`](@ref), [`orientation()`](@ref), [`right_turn`](@ref).
"""
left_turn(p::Point2, q::Point2, r::Point2)

"""
    lexicographically_xy_larger(p::Point2, q::Point2)

Returns `true` iff `p` is lexicographically larger than `q` with respect to
``xy`` order.

See also: [`compare_xy()`](@ref),
[`lexicographically_xy_larger_or_equal()`](@ref),
[`lexicographically_xy_smaller()`](@ref),
[`lexicographically_xy_smaller_or_equal()`](@ref).
"""
lexicographically_xy_larger(p::Point2, q::Point2)

"""
    lexicographically_xy_larger_or_equal(p::Point2, q::Point2)

Returns `true` iff `p` is lexicographically not smaller than `q` with respect to
``xy`` order.

See also: [`compare_xy()`](@ref),
[`lexicographically_xy_larger()`](@ref),
[`lexicographically_xy_smaller()`](@ref),
[`lexicographically_xy_smaller_or_equal()`](@ref).
"""
lexicographically_xy_larger_or_equal(p::Point2, q::Point2)

"""
    lexicographically_xy_smaller(p::Point2, q::Point2)

Returns `true` iff `p` is lexicographically smaller than `q` with respect to
``xy`` order.

See also: [`compare_xy()`](@ref),
[`lexicographically_xy_larger()`](@ref),
[`lexicographically_xy_larger_or_equal()`](@ref),
[`lexicographically_xy_smaller_or_equal()`](@ref).
"""
lexicographically_xy_smaller(p::Point2, q::Point2)

"""
    lexicographically_xy_smaller_or_equal(p::Point2, q::Point2)

Returns `true` iff `p` is lexicographically not larger than `q` with respect to
``xy`` order.

See also: [`compare_xy()`](@ref),
[`lexicographically_xy_larger()`](@ref),
[`lexicographically_xy_larger_or_equal()`](@ref),
[`lexicographically_xy_smaller()`](@ref).
"""
lexicographically_xy_smaller_or_equal(p::Point2, q::Point2)

"""
    lexicographically_xyz_smaller(p::Point3, q::Point3)

Returns `true` iff `p` is lexicographically smaller than `q` with respect to
``xyz`` order.

See also: [`compare_xyz()`](@ref),
[`lexicographically_xyz_smaller_or_equal()`](@ref).
"""
lexicographically_xyz_smaller(p::Point3, q::Point3)

"""
    lexicographically_xyz_smaller_or_equal(p::Point3, q::Point3)

Returns `true` iff `p` is lexicographically not larger than `q` with respect to
``xyz`` order.

See also: [`compare_xyz()`](@ref),
[`lexicographically_xyz_smaller()`](@ref).
"""
lexicographically_xyz_smaller_or_equal(p::Point3, q::Point3)

"""
    max_vertex(ir::IsoRectangle2)

Computes the vertex with the lexicographically largest coordinates of the iso
rectangle `ir`.
"""
max_vertex(ir::IsoRectangle2)

"""
    max_vertex(ic::IsoCuboid3)

Computes the vertex with the lexicographically largest coordinates of the iso
cuboid `ic`.
"""
max_vertex(ic::IsoCuboid3)

"""
    midpoint(p::Point23, q::Point23)

Computes the midpoint of the segment `pq`.
"""
midpoint(p::Point23, q::Point23)

"""
    min_vertex(ir::IsoRectangle2)

Computes the vertex with the lexicographically smallest coordinates of the iso
rectangle `ir`.
"""
min_vertex(ir::IsoRectangle2)

"""
    min_vertex(ic::IsoCuboid3)

Computes the vertex with the lexicographically smallest coordinates of the iso
cuboid `ic`.
"""
min_vertex(ic::IsoCuboid3)

"""
    normal(p::Point3, q::Point3, r::Point3)

Computes the normal vector for the vectors `q-p` and `r-p`.

!!! info "Precondition"

    The points `p`, `q`, and `r` must not be collinear.
"""
normal(p::Point3, q::Point3, r::Point3)

"""
    orientation(p::Point2, q::Point2, r::Point2)

Returns [`LEFT_TURN`](@ref), if `r` lies to the left of the oriented line `l`
defined by `p` and `q`, returns [`RIGHT_TURN`](@ref) if `r` lies to the right
of `l`, and returns [`COLLINEAR`](@ref) if `r` lies on `l`.

See also: [`collinear()`](@ref), [`left_turn()`](@ref), [`right_turn()`](@ref).
"""
orientation(p::Point2, q::Point2, r::Point2)

"""
    orientation(u::Vector2, v::Vector2, w::Vector2)

Returns [`LEFT_TURN`](@ref), if `u` and `v` form a left turn, returns
[`RIGHT_TURN`](@ref) if `u` and `v` form a right turn, and returns
[`COLLINEAR`](@ref) if `u` and `v` are collinear.

See also: [`collinear()`](@ref), [`left_turn()`](@ref), [`right_turn()`](@ref).
"""
orientation(u::Vector2, v::Vector2, w::Vector2)

"""
    orientation(p::Point3, q::Point3, r::Point3, s::Point3)

Returns [`POSITIVE`](@ref), if `s` lies on the positive side of the oriented
plane `h` defined by `p`, `q`, and `r`, returns [`NEGATIVE`](@ref) if if `s`
lies on the negative side of `h`, and returns [`COPLANAR`](@ref) if `s` lies on
`h`.

See also: [`collinear()`](@ref), [`left_turn()`](@ref), [`right_turn()`](@ref).
"""
orientation(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    orientation(u::Vector3, v::Vector3, w::Vector3)

Returns [`NEGATIVE`](@ref) if `u`, `v` and `w` are negatively oriented,
[`POSITIVE`](@ref) if `u`, `v` and `w` are positively oriented, and
[`COPLANAR`](@ref) if `u`, `v` and `w` are coplanar.

See also: [`collinear()`](@ref), [`left_turn()`](@ref), [`right_turn()`](@ref).
"""
orientation(u::Vector3, v::Vector3, w::Vector3)

"""
    orthogonal_vector(p::Point3, q::Point3, r::Point3)

Computes an orthogonal vector of the plane defined by `p`, `q` and `r`,  which
is directed to the positive side of this plane.
"""
orthogonal_vector(p::Point3, q::Point3, r::Point3)

for (T, t) ∈ ((:Line23, :l)
            , (:Ray23, :r)
            , (:Segment23, :s)
            , (:Plane3, :h))
    t1, t2 = string.(t, ("₁", "₂"))
    args = join(string.((t1, t2), "::", T), ", ")
    @eval begin
        """
            parallel($($args))

        Returns `true`, if `$($t1)` and `$($t2)` are parallel or if one of those
        (or both) is degenerate.
        """
        parallel($t1::$T, $t2::$T)
    end
end

"""
    radical_line(c₁::Circle2, c₂::Circle2)

Returns the radical line of the two circles.

!!! info "Precondition"

    `c₁` and `c₂` are not concentric.
"""
radical_line(c₁::Circle2, c₂::Circle2)

"""
    radical_plane(s₁::Sphere3, s₂::Sphere3)

Returns the radical plane of the two spheres.

!!! info "Precondition"

    `s₁` and `s₂` are not concentric.
"""
radical_plane(s₁::Sphere3, s₂::Sphere3)

"""
    rational_rotation_approximation(dirx::Real, diry::Real
                                   ;num_eps::Real = eps(Float64),
                                    den_eps::Real = 1)

Returns integers `sin_num`, `cos_num` and `denom`, such that `sin_num`/`denom`
approximates the sine of direction `(dirx,diry)`.

The difference between the sine and the approximating rational is bounded by
`num_eps`/`den_eps`.

!!! info "Precondition"

    `num_eps` ``≠ 0``.

**Implementation**

The approximation is based on Farey sequences as described in the rational
rotation method presented by Canny and Ressler at the 8th SoCG 1992. The authors
of the implementation use a slower version which needs no division operation in
the approximation.

See also: [`AffTransformation2`](@ref)
"""
rational_rotation_approximation(dirx::Real, diry::Real
                               ;num_eps::Real=eps(Float64), den_eps::Real=1) =
    iszero(num_eps) ? error("num_eps = 0") :
    let (sin_num, cos_num, denom) = Ref{RT}.((0, 0, 0)),
        res = rational_rotation_approximation(
                  convert(RT, dirx), convert(RT, diry),
                  sin_num, cos_num, denom,
                  convert(RT, num_eps), convert(RT, den_eps))

        sin_num[], cos_num[], denom[]
    end

"""
    right_turn(p::Point2, q::Point2, r::Point2)

Returns `true` iff `p`, `q`, and `r` form a right turn.

See also: [`collinear()`](@ref), [`left_turn()`](@ref), [`orientation()`](@ref).
"""
right_turn(p::Point2, q::Point2, r::Point2)

"""
    scalar_product(u::Vector23, v::Vector23)

Returns the scalar product of `u` and `v`.
"""
scalar_product(u::Vector23, v::Vector23)

"""
    side_of_bounded_circle(p::Point2, q::Point2, r::Point2, t::Point2)

Returns the relative position of point `t` to the circle defined by `p`, `q` and
`r`.

The order of the points `p`, `q`, and `r` does not matter.

!!! info "Precondition"

    `p`, `q`, and `r` are not collinear.

See also: [`coplanar_side_of_bounded_circle()`](@ref),
[`side_of_oriented_circle()`](@ref).
"""
side_of_bounded_circle(p::Point2, q::Point2, r::Point2, t::Point2)

"""
    side_of_bounded_circle(p::Point2, q::Point2, t::Point2)

Returns the position of the point `t` relative to the circle that has `pq` as
its diameter.

See also: [`coplanar_side_of_bounded_circle()`](@ref),
[`side_of_oriented_circle()`](@ref).
"""
side_of_bounded_circle(p::Point2, q::Point2, t::Point2)

"""
    side_of_bounded_sphere(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

Returns the relative position of point `t` to the sphere defined by `p`, `q`,
`r`, and `s`.

The order of the points `p`, `q`, `r`, and `s` does not matter.

!!! info "Precondition"

    `p`, `q`, `r` and `s` are not coplanar.

See also: [`side_of_oriented_sphere()`](@ref).
"""
side_of_bounded_sphere(p::Point3, q::Point3, r::Point3, s::Point3, t::Point3)

"""
    side_of_bounded_sphere(p::Point3, q::Point3, r::Point3, t::Point3)

Returns the position of the point `t` relative to the sphere passing through
`p`, `q`, and `r` and whose center is in the plane defined by these three
points.

See also: [`side_of_oriented_sphere()`](@ref).
"""
side_of_bounded_sphere(p::Point3, q::Point3, r::Point3, t::Point3)

"""
    side_of_bounded_sphere(p::Point3, q::Point3, t::Point3)

Returns the position of the point `t` relative to the sphere that has `pq` as
its diameter.

See also: [`side_of_oriented_sphere()`](@ref).
"""
side_of_bounded_sphere(p::Point3, q::Point3, t::Point3)

"""
    side_of_oriented_circle(p::Point2, q::Point2, r::Point2, test::Point2)

Returns the relative position of point `test` to the oriented circle defined by
`p`, `q` and `r`.

The order of the points `p`, `q` and `r` is important, since it determines the
orientation of the implicitly constructed circle.

If `p`, `q` and `r` are collinear, the circle degenerates in a line.
[`ON_ORIENTED_BOUNDARY`](@ref) is returned if `test` is also collinear or if two
points are identical, otherwise, `side_of_oriented_circle(r, q, test, p)` is
returned.

See also: [`side_of_bounded_circle()`](@ref)
"""
side_of_oriented_circle(p::Point2, q::Point2, r::Point2, test::Point2)

"""
    side_of_oriented_sphere(p::Point3, q::Point3, r::Point3, s::Point3, test::Point3)

Returns the relative position of point `test` to the oriented circle defined by
`p`, `q`, `r` and `s`.

The order of the points `p`, `q`, `r` and `s` is important, since it determines
the orientation of the implicitly constructed sphere. If the points `p`, `q`,
`r` and `s` are positively oriented, positive side is the bounded interior of
the sphere.

In case of degeneracies, [`ON_ORIENTED_BOUNDARY`](@ref) is returned if all
points are coplanar. Otherwise, there is a cyclic permutation of the five points
that puts four non coplanar points first, it is used to answer the predicate:
e.g. `side_of_oriented_sphere(q, r, s, test, p)` is returned if `q`, `r`, `s`
and `test` are non coplanar.

See also: [`side_of_bounded_sphere()`](@ref)
"""
side_of_oriented_sphere(p::Point3, q::Point3, r::Point3, s::Point3, test::Point3)

"""
    squared_area(p::Point3, q::Point3, r::Point3)

Returns the squared area of the triangle defined by the points `p`, `q` and `r`.
"""
squared_area(p::Point3, q::Point3, r::Point3)

"""
    squared_distance(obj₁::Type1, obj₂::Type2)

Computes the square of the Euclidean distance between two geometric objects.

For arbitrary geometric objects `obj₁` and `obj₂`, the squared distance is
defined as the minimal `squared_distance(p₁, p₂)`, where `p₁` is a point of
`obj₁` and `p₂` is a point of `obj₂`. Note that for objects that have an inside
(a bounded region), this inside is part of the object. So, the squared distance
from a point inside is zero, not the squared distance to the closest point on
the boundary.

In 2D, the types `Type2` and `Type2` can be any of the following:

- [`Point2`](@ref)
- [`Line2`](@ref)
- [`Ray2`](@ref)
- [`Segment2`](@ref)
- [`Triangle2`](@ref).

In 3D, the types `Type1` and `Type2` can be any of the following:

- [`Point3`](@ref)
- [`Line3`](@ref)
- [`Ray3`](@ref)
- [`Segment3`](@ref)
- [`Plane3`](@ref).

See also: [`compare_distance_to_point()`](@ref),
[`compare_signed_distance_to_line()`](@ref),
[`compare_signed_distance_to_plane()`](@ref),
[`has_larger_distance_to_point()`](@ref),
[`has_larger_signed_distance_to_line()`](@ref),
[`has_larger_signed_distance_to_plane()`](@ref),
[`has_smaller_distance_to_point()`](@ref),
[`has_smaller_signed_distance_to_line()`](@ref),
[`has_smaller_signed_distance_to_plane()`](@ref).
"""
squared_distance(obj₁, obj₂)

"""
    squared_radius(p::Point2, q::Point2, r::Point2)

Compute the squared radius of the circle passing through the points `p`, `q`,
and `r`.

!!! info "Precondition"

    `p`, `q`, and `r` are not collinear.
"""
squared_radius(p::Point2, q::Point2, r::Point2)

"""
    squared_radius(p::Point23, q::Point23)

Compute the squared radius of the smallest circle passing through `p`, and `q`,
i.e. one fourth of the squared distance between `p` and `q`.
"""
squared_radius(p::Point23, q::Point23)

"""
    squared_radius(p::Point23)

Computes the squared radius of the smallest circle passing through `p`, i.e.
``0``.
"""
squared_radius(p::Point23)

"""
    squared_radius(p::Point3, q::Point3, r::Point3, s::Point3)

Computes the squared radius of the sphere passing through the points `p`, `q`,
`r` and `s`.

!!! info "Precondition"

    `p`, `q`, `r`, and `s` are not coplanar.
"""
squared_radius(p::Point3, q::Point3, r::Point3, s::Point3)

"""
    squared_radius(p::Point3, q::Point3, r::Point3)

Computes the squared radius of the sphere passing through the points `p`, `q`,
and `r` and whose center is in the same plane as those three points.
"""
squared_radius(p::Point3, q::Point3, r::Point3)

"""
    unit_normal(p::Point3, q::Point3, r::Point3)

Computes the unit normal vector for the vectors `q-p` and `r-p`.

!!! info "Precondition"

    The points `p`, `q`, and `r` must not be collinear.
"""
unit_normal(p::Point3, q::Point3, r::Point3)

"""
    volume(p₀::Point3, p₁::Point3, p₂::Point3, p₃::Point3)

Computes the signed volume of the tetrahedron defined by the four points `p₀`,
`p₁`, `p₂` and `p₃`.
"""
volume(p₀::Point3, p₁::Point3, p₂::Point3, p₃::Point3)

"""
    x_equal(p::Point23, q::Point23)

Returns `true`, iff `p` and `q` have the same x-coordinate.

See also: [`compare_x()`](@ref), [`y_equal()`](@ref), [`z_equal()`](@ref).
"""
x_equal(p::Point23, q::Point23)

"""
    y_equal(p::Point23, q::Point23)

Returns `true`, iff `p` and `q` have the same y-coordinate.

See also: [`compare_y()`](@ref), [`x_equal()`](@ref), [`z_equal()`](@ref).
"""
y_equal(p::Point23, q::Point23)

"""
    z_equal(p::Point3, q::Point3)

Returns `true`, iff `p` and `q` have the same z-coordinate.

See also: [`compare_z()`](@ref), [`x_equal()`](@ref), [`y_equal()`](@ref).
"""
z_equal(p::Point3, q::Point3)

"""
    +(p::Point23, v::Vector23)

Returns the point obtained by translating `p` by the vector `v`.
"""
+(p::Point23, v::Vector23)

"""
    -(p::Point23, v::Vector23)

Returns the point obtained by translating `p` by the vector -`v`.
"""
-(p::Point23, v::Vector23)

"""
    *(v::Vector23, s::Real)

Multiplication with a scalar from the right.
"""
*(v::Vector23, s::Real)

"""
    *(s::Real, v::Vector23)

Multiplication with a scalar from the left.
"""
*(s::Real, v::Vector23)

for V ∈ (:Vector2, :Vector3)
    @eval begin
        @cxxdereference Base.:*(v::$V, s::Real) = v * convert(RT, s)
        @cxxdereference Base.:*(s::Real, v::$V) = convert(RT, s) * v
    end
end

"""
    do_overlap(bb₁::Bbox23, bb₂::Bbox23)

Returns `true` iff `bb₁` and `bb₂` overlap, i.e., iff their intersection is
non-empty.
"""
do_overlap(bb₁::Bbox23, bb₂::Bbox23)


text:
module ProximalBase

export
  # types
  DifferentiableFunction,
  ProximableFunction,
  SparseIterate, SymmetricSparseIterate, AtomIterate, numCoordinates,

  # smooth functions
  QuadraticFunction, L2Loss, LeastSquaresLoss,
  value, value_and_gradient!, gradient!,

  # proximal functions
  ProxZero,
  ProxL1,
  ProxL2, ProxL2Sq, ProxNuclear,
  ProxL1Fused,
  ProxGaussLikelihood,
  prox!, prox, cdprox!,

  # utils
  shrink, shrinkL2!,
  proxL1Fused,
  A_mul_B_row, At_mul_B_row,
  A_mul_X_mul_B, A_mul_X_mul_B_rc,
  A_mul_UUt_mul_B, A_mul_UUt_mul_B_rc,
  ind2subLowerTriangular,sub2indLowerTriangular,
  norm_diff,

  # function for SparseIterate
  nnz

using SparseArrays
using LinearAlgebra

import Base: *

include("sparseIterate.jl")

include("utils.jl")
# DifferentiableFunctions
include("differentiable_functions.jl")

# ProximableFunctions
include("proximal_functions.jl")

end


text:
function ComputeCentroids(X,idx)
    n=size(X,2)
    k=size(unique(idx),1)
    C=zeros(k,n)

    for i=1:k
        C[i,:]=mean(X[idx.==i,:],dims=1)
    end

    return C
end


text:
# https://www.schemecolor.com/stamping-it-out.php
palette_StampingItOut = [
                         [42, 55, 52],
                         [34, 99, 102],
                         [213, 176, 138],
                         [251, 204, 153],
                         [250, 231, 181],
                        ]

# https://www.schemecolor.com/quietly-ethusiastic.php
palette_QuietlyEnthusiastic = [
                               [37, 92, 156],
                               [47, 64, 115],
                               [128, 183, 191],
                               [187, 209, 211],
                               [117, 126, 148],
                               [79, 88, 108],
                              ]


text:
using Omega
using Test

function testpw()
  x = normal(0.0, 1.0)
  f(x::Real, y::Real) = x * y
  y = pw() do
    f(x, x)
  end
  x_, y_ = rand((x, y))
  @test y_ == x_ * x_
end


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

using BAT
using Test

using DensityInterface
using Distributions, Statistics, StatsBase, IntervalSets, ValueShapes

@testset "truncated_density" begin
    parent_dist = NamedTupleDist(a = Normal(), b = Weibull())
    vs = varshape(parent_dist)
    logrenormf = 4.2
    parent_density = convert(AbstractDensity, parent_dist)

    @test @inferred(BAT.renormalize_density(parent_density, logrenormf)) isa BAT.RenormalizedDensity
    density = renormalize_density(parent_density, logrenormf)

    @test @inferred(parent(density)) === parent_density
    @test @inferred(BAT.var_bounds(density)) == BAT.var_bounds(parent_density)
    @test @inferred(varshape(density)) == varshape(parent_density)
    @test @inferred(unshaped(density)) == renormalize_density(unshaped(parent_density), logrenormf)
    @test @inferred(vs(unshaped(density))) == renormalize_density(vs(unshaped(parent_density)), logrenormf)

    v = rand(parent_dist)
    @test @inferred(BAT.checked_logdensityof(density, v)) == BAT.checked_logdensityof(parent_density, v) + logrenormf
    @test @inferred(DensityInterface.logdensityof(density, v)) == DensityInterface.logdensityof(parent_density, v) + logrenormf
    @test @inferred(logdensityof(density, v)) == logdensityof(parent_density, v) + logrenormf

    rng = bat_rng()
    @test @inferred(rand(deepcopy(rng), BAT.sampler(density), 10)) == rand(deepcopy(rng), BAT.sampler(parent_density), 10)
    @test @inferred(rand(deepcopy(rng), BAT.bat_sampler(density), 10)) == rand(deepcopy(rng), BAT.bat_sampler(parent_density), 10)

    @test cov(unshaped(density)) == cov(unshaped(parent_density))

    @test @inferred(renormalize_density(density, 7.9)) == renormalize_density(parent_density, logrenormf + 7.9)
end