text:
# Linear Transformed Beta distribution

module TransformedBetaDistributions

if VERSION >= v"0.7"
  using Statistics
  using Distributed
  import Statistics: mean, median, maximum, minimum, quantile, std, var, cov, cor
else
  using Compat
  using Compat.Statistics
  using Compat.Distributed
  import Base: mean, median, maximum, minimum, quantile, std, var, cov, cor
end

import Base.length, Base.show
using Distributions
import Distributions.params, Distributions.@check_args
import Distributions.rand
import Distributions.pdf, Distributions.logpdf, Distributions.cdf, Distributions.gradlogpdf
import Distributions.quantile, Distributions.insupport
import Distributions.minimum, Distributions.maximum
import Distributions.mean, Distributions.var
import Distributions.mode, Distributions.modes

export LinearTransformedBeta
export params, rand, pdf, logpdf, cdf, gradlogpdf, quantile, insupport
export minimum, maximum, mean, var, mode, modes

struct LinearTransformedBeta{T<:Real} <: ContinuousUnivariateDistribution
    dist::Beta{T}
    xmin::T
    xmax::T

    function LinearTransformedBeta{T}(α::T, β::T; xmin::T = zero(T), xmax::T = one(T) ) where T
        @check_args(LinearTransformedBeta, α > zero(α) && β > zero(β))
        @check_args(LinearTransformedBeta, -Inf < xmin < xmax < Inf )
        new{T}(Beta(α, β),xmin,xmax)
    end
end
LinearTransformedBeta(α::T, β::T; xmin::T=zero(T), xmax::T = one(T) ) where T = LinearTransformedBeta{T}(α, β, xmin=xmin, xmax=xmax)

#### Conversions
function convert(::Type{LinearTransformedBeta{T}}, α::Real, β::Real, xmin::Real, xmax::Real) where T<:Real
    LinearTransformedBeta(T(α), T(β), T(xmin), T(xmax) )
end
function convert(::Type{LinearTransformedBeta{T}}, d::LinearTransformedBeta{S}) where {T <: Real, S <: Real}
    LinearTransformedBeta(T(d.α), T(d.β), T(d.xmin), T(d.xmax) )
end

### Parameters
params(d::LinearTransformedBeta) = (d.dist.α, d.dist.β,d.xmin,d.xmax)
@inline partype(d::LinearTransformedBeta{T}) where {T<:Real} = T

### Sampling
rand(d::LinearTransformedBeta) = d.xmin+(d.xmax-d.xmin)*rand(d.dist)
#sampler

### Evaluation
pdf(d::LinearTransformedBeta, x::T) where T<:Real = pdf(d.dist, (x.-d.xmin)./(d.xmax-d.xmin) )
logpdf(d::LinearTransformedBeta, x::T) where T<:Real = logpdf(d.dist, (x.-d.xmin)./(d.xmax-d.xmin) )
cdf(d::LinearTransformedBeta, q::T) where T<:Real = cdf(d.dist, (x.-d.xmin)./(d.xmax-d.xmin) )
quantile(d::LinearTransformedBeta, q::T) where T<:Real = quantile(d.dist, (x.-d.xmin)./(d.xmax-d.xmin) )
minimum(d::LinearTransformedBeta) = d.xmin
maximum(d::LinearTransformedBeta) = d.xmax

function insupport(d::LinearTransformedBeta, x::T)  where T<:Real
  d.xmin<d.xmax && insupport(d.dist, (x.-d.xmin)./(d.xmax-d.xmin) )
end

function gradlogpdf(d::LinearTransformedBeta, x::T) where T<:Real
  (d.xmax-d.xmin)*gradlogpdf(d.dist, (x-d.xmin)/(d.xmax-d.xmin) )
end

### Basic statistics
mean(d::LinearTransformedBeta) = d.xmin+(d.xmax-d.xmin)*mean(d.dist)
var(d::LinearTransformedBeta) = (d.xmax-d.xmin)^2*var(d.dist)
mode(d::LinearTransformedBeta) = d.xmin+(d.xmax-d.xmin)*mode(d.dist)
modes(d::LinearTransformedBeta) = d.xmin.+(d.xmax-d.xmin).*modes(d.dist)
# TODO skewness
# TODO kurtosis
# TODO entropy
# TODO mgf
# TODO cf

### Show

distrname(d::LinearTransformedBeta) = "LinearTransformedBeta"
function Base.show(io::IO, d::LinearTransformedBeta)
  show(io,distrname(d) * "(α=" * string(d.dist.α) * ", β=" * string(d.dist.β) * ", min=" * string(d.xmin) * " ,max=" * string(d.xmax) * ")" )
end

end # module


text:
using MVDA, Random, LinearAlgebra, StableRNGs
using StatsBase, Statistics
using ProgressMeter

# Compute an initial solution w/ regularization. Take this estimate as the ground truth.
function __init__(λ, problem, ϵ, gtol)
    MVDA.init!(MMSVD(), problem, ϵ, λ, verbose=false, maxiter=10^6, nesterov_threshold=0, gtol=gtol)
    B0 = copy(problem.coeff.all)
    return B0
end

# Solve for s-sparse solution.
function __solve__(s, problem, ϵ, gtol, dtol)
    MVDA.fit(MMSVD(), problem, ϵ, s, verbose=false, ninner=10^6, nouter=200, nesterov_threshold=0, gtol=gtol, dtol=dtol, rtol=0.0)
    B = problem.proj.all
    return B
end

# Check classification accuracy against X.
function __accuracy__(B, problem, targets, X)
    n = problem.n
    copyto!(problem.proj.all, B)
    labels = MVDA.classify(problem, X)
    return round(100 * (1 - sum(labels .== targets) / n), sigdigits=4)
end

# Check MSE against ground truth B0.
mse(B, B0) = mean( (B - B0) .^ 2 ) / dot(B0, B0)

# Check contingency table by comparing estimate x against ground truth y.
function discovery_metrics(x, y)
    TP = FP = TN = FN = 0
    for (xi, yi) in zip(eachrow(x), eachrow(y))
        nxi, nyi = norm(xi), norm(yi)
        TP += (nxi != 0) && (nyi != 0)
        FP += (nxi != 0) && (nyi == 0)
        TN += (nxi == 0) && (nyi == 0)
        FN += (nxi == 0) && (nyi != 0)
    end
    return (TP, FP, TN, FN)
end

function run(filename, λ)
    # Fix problem parameters.
    c = 10
    p = 50
    rng = StableRNG(1903)

    # Define range of values for various simulation settings.
    separation_levels = (1.0, 3.0)
    noise_levels = (1e0, 5e0, 1e1)
    dataset_sizes = (10, 40, 100)
    sparsity_levels = Float64[1 - k/p for k in p:-1:0]

    # Allocate arrays for interesting metrics.
    m = length(noise_levels)
    train_error = zeros(p+1, m)
    mean_squared_error = zeros(p+1, m)

    @info "Initializing ouptut file: $(filename).dat"
    open("$(filename).dat", "w+") do io
        write(io, "n,p,c,nclass,ncausal,separation,sigma,sparsity,error,MSE,TP,FP,TN,FN\n")
    end

    for samples_per_class in dataset_sizes, separation in separation_levels
        for (j, sigma) in enumerate(noise_levels)
            # Simulate ground truth with the given parameters.
            n = samples_per_class*c
            targets, X, B0 = MVDA.simulate_ground_truth(p, c, samples_per_class, d=separation, rng=rng, sigma=sigma)

            # Standardize ALL data based on the training set.
            F = StatsBase.fit(ZScoreTransform, X, dims=1)
            has_nan = any(isnan, F.scale) || any(isnan, F.mean)
            has_inf = any(isinf, F.scale) || any(isinf, F.mean)
            has_zero = any(iszero, F.scale)
            if has_nan
                error("Detected NaN in z-score.")
            elseif has_inf
                error("Detected Inf in z-score.")
            elseif has_zero
                for idx in eachindex(F.scale)
                    x = F.scale[idx]
                    F.scale[idx] = ifelse(iszero(x), one(x), x)
                end
            end
            StatsBase.transform!(F, X)

            # Figure out how to group samples together and create a problem instance.
            sortidx = sortperm(targets)
            problem = MVDAProblem(targets, X, intercept=false)

            # Fix hyperparameter controlling distance insensitivity.
            ϵ = MVDA.maximal_deadzone(problem)

            # Create closures specific to our problem instance.
            make_data_plot() = __make_data_plot__(targets, n, c, samples_per_class, ccp, colors)
            make_distance_plot(B) = __make_distance_plot__(B, problem, targets, ϵ, sortidx, ccp, colors)
            init() = __init__(λ, problem, ϵ, 1e-6)
            solve(s) = __solve__(s, problem, ϵ, 1e-3, 1e-3)
            accuracy(B) = __accuracy__(B, problem, targets, X)

            @showprogress "$(samples_per_class) per class / d = $(separation) / σ = $(sigma)... " for (i, s) in enumerate(sparsity_levels)
                # Estimate s-sparse coefficients.
                if s != 0
                    B = solve(s)
                else
                    B = init()
                end

                train_error[i,j] = accuracy(B)
                mean_squared_error[i,j] = mse(B, B0)
                TP, FP, TN, FN = discovery_metrics(B[1:p,:], B0[1:p,:])

                open("$(filename).dat", "a+") do io
                    write(io, "$(n),$(p),$(c),$(samples_per_class),$(c),$(separation),$(sigma),$(s),$(train_error[i,j]),$(mean_squared_error[i,j]),$(TP),$(FP),$(TN),$(FN)\n")
                end
            end
        end
    end
end

# Parse input arguments and run.
length(ARGS) < 2 && error("Script requires a filename (no extension required) and a value λ > 0 for regularization.")
filename = ARGS[1]
λ = parse(Float64, ARGS[2])
@info "Output file: $(filename)"
@info "λ value: $(λ)"
run(filename, λ)

ext: jl


text:
using RvSpectMLBase
using Test

@testset "ChunkOfSpectrum" begin
        num_pixels = 100
        num_orders = 10
        lambda = reshape(range(5000.0, stop=6000, length=num_pixels*num_orders), num_pixels,num_orders)
        data = ones(num_pixels,num_orders)
        order = 5
        pixels = 20:30
        @test_nowarn ChunkOfSpectrum(lambda, data, data, (pixels=pixels, order=order) )
        ch = ChunkOfSpectrum(lambda, data, data, (pixels=pixels, order=order) )
        @test RvSpectMLBase.calc_normalization(ch) ≈ 1
        @test RvSpectMLBase.calc_normalization_var_weighted(ch) ≈ 1
        @test RvSpectMLBase.find_orders_with_line(5421.0,lambda)[1] == order
        @test RvSpectMLBase.find_pixels_for_line_in_chunk(ch, 5420, 5422) == 2:3
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
GEOMTYPE = Dict( GEOS_POINT => :Point,
                 GEOS_LINESTRING => :LineString,
                 GEOS_LINEARRING => :LinearRing,
                 GEOS_POLYGON => :Polygon,
                 GEOS_MULTIPOINT => :MultiPoint,
                 GEOS_MULTILINESTRING => :MultiLineString,
                 GEOS_MULTIPOLYGON => :MultiPolygon,
                 GEOS_GEOMETRYCOLLECTION => :GeometryCollection)

for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval writegeom(obj::$geom, wktwriter::WKTWriter, context::GEOSContext = _context) = _writegeom(obj.ptr, wktwriter, context)
    @eval writegeom(obj::$geom, wkbwriter::WKBWriter, context::GEOSContext = _context) = _writegeom(obj.ptr, wkbwriter, context)
    @eval writegeom(obj::$geom, context::GEOSContext = _context) = _writegeom(obj.ptr, context)
end

function geomFromGEOS(ptr::GEOSGeom)
    if geomTypeId(ptr) == GEOS_POINT
        return Point(ptr)
    elseif geomTypeId(ptr) == GEOS_LINESTRING
        return LineString(ptr)
    elseif geomTypeId(ptr) == GEOS_LINEARRING
        return LinearRing(ptr)
    elseif geomTypeId(ptr) == GEOS_POLYGON
        return Polygon(ptr)
    elseif geomTypeId(ptr) == GEOS_MULTIPOINT
        return MultiPoint(ptr)
    elseif geomTypeId(ptr) == GEOS_MULTILINESTRING
        return MultiLineString(ptr)
    elseif geomTypeId(ptr) == GEOS_MULTIPOLYGON
        return MultiPolygon(ptr)
    else @assert geomTypeId(ptr) == GEOS_GEOMETRYCOLLECTION
        return GeometryCollection(ptr)
    end
end

readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext=_context) =
    geomFromGEOS(_readgeom(wktstring, wktreader, context))
readgeom(wktstring::String, context::GEOSContext=_context) =
    readgeom(wktstring, WKTReader(context), context)

readgeom(wkbbuffer::Vector{Cuchar}, wkbreader::WKBReader, context::GEOSContext=_context) =
    geomFromGEOS(_readgeom(wkbbuffer, wkbreader, context))
readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext=_context) =
    readgeom(wkbbuffer, WKBReader(context), context)

# -----
# Linear referencing functions -- there are more, but these are probably sufficient for most purposes
# -----
project(line::LineString, point::Point) = project(line.ptr, point.ptr)
projectNormalized(line::LineString, point::Point) = projectprojectNormalized(line.ptr, point.ptr)
interpolate(line::LineString, dist::Real) = Point(interpolate(line.ptr, dist))
interpolateNormalized(line::LineString, dist::Real) = Point(interpolateNormalized(line.ptr, dist))

# # -----
# # Topology operations
# # -----
for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval buffer(obj::$geom, dist::Real, quadsegs::Integer=8) = geomFromGEOS(buffer(obj.ptr, dist, quadsegs))
    @eval envelope(obj::$geom) = geomFromGEOS(envelope(obj.ptr))
    @eval convexhull(obj::$geom) = geomFromGEOS(convexhull(obj.ptr))
    @eval boundary(obj::$geom) = geomFromGEOS(boundary(obj.ptr))
    @eval unaryUnion(obj::$geom) = geomFromGEOS(unaryUnion(obj.ptr))
    @eval pointOnSurface(obj::$geom) = Point(pointOnSurface(obj.ptr))
    @eval centroid(obj::$geom) = Point(centroid(obj.ptr))
    @eval node(obj::$geom) = geomFromGEOS(node(obj.ptr))
end

for g1 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    for g2 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
        @eval intersection(obj1::$g1, obj2::$g2) = geomFromGEOS(intersection(obj1.ptr, obj2.ptr))
        @eval difference(obj1::$g1, obj2::$g2) = geomFromGEOS(difference(obj1.ptr, obj2.ptr))
        @eval symmetricDifference(obj1::$g1, obj2::$g2) = geomFromGEOS(symmetricDifference(obj1.ptr, obj2.ptr))
        @eval union(obj1::$g1, obj2::$g2) = geomFromGEOS(union(obj1.ptr, obj2.ptr))
    end
end

# # all arguments remain ownership of the caller (both Geometries and pointers)
# function polygonize(geoms::Vector{GEOSGeom})
#     result = GEOSPolygonize(pointer(geoms), length(geoms))
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSPolygonize")
#     end
#     result
# end
# # GEOSPolygonizer_getCutEdges
# # GEOSPolygonize_full

# function lineMerge(ptr::GEOSGeom)
#     result = GEOSLineMerge(ptr)
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSLineMerge")
#     end
#     result
# end

for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval simplify(obj::$geom, tol::Real) = geomFromGEOS(simplify(obj.ptr, tol))
    @eval topologyPreserveSimplify(obj::$geom, tol::Real) = geomFromGEOS(topologyPreserveSimplify(obj.ptr, tol))
    @eval uniquePoints(obj::$geom) = MultiPoint(uniquePoints(obj.ptr))
    @eval delaunayTriangulationEdges(obj::$geom, tol::Real=0.0) = MultiLineString(delaunayTriangulation(obj.ptr, tol, true))
    @eval delaunayTriangulation(obj::$geom, tol::Real=0.0) = GeometryCollection(delaunayTriangulation(obj.ptr, tol, false))
end

sharedPaths(obj1::LineString, obj2::LineString) = GeometryCollection(sharedPaths(obj1.ptr, obj2.ptr))

# # Snap first geometry on to second with given tolerance
for g1 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    for g2 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
        @eval snap(obj1::$g1, obj2::$g2, tol::Real) = geomFromGEOS(snap(obj1.ptr, obj2.ptr, tol))
    end
end

# -----
# Binary predicates
# -----
for g1 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    for g2 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
        @eval disjoint(obj1::$g1, obj2::$g2) = disjoint(obj1.ptr, obj2.ptr)
        @eval touches(obj1::$g1, obj2::$g2) = touches(obj1.ptr, obj2.ptr)
        @eval intersects(obj1::$g1, obj2::$g2) = intersects(obj1.ptr, obj2.ptr)
        @eval crosses(obj1::$g1, obj2::$g2) = crosses(obj1.ptr, obj2.ptr)
        @eval within(obj1::$g1, obj2::$g2) = within(obj1.ptr, obj2.ptr)
        @eval contains(obj1::$g1, obj2::$g2) = contains(obj1.ptr, obj2.ptr)
        @eval overlaps(obj1::$g1, obj2::$g2) = overlaps(obj1.ptr, obj2.ptr)
        @eval equals(obj1::$g1, obj2::$g2) = equals(obj1.ptr, obj2.ptr)
        @eval equalsexact(obj1::$g1, obj2::$g2, tol::Real) = equalsexact(obj1.ptr, obj2.ptr, tol)
        @eval covers(obj1::$g1, obj2::$g2) = covers(obj1.ptr, obj2.ptr)
        @eval coveredby(obj1::$g1, obj2::$g2) = coveredby(obj1.ptr, obj2.ptr)
    end
end

# # -----
# # Prepared Geometry Binary predicates - return 2 on exception, 1 on true, 0 on false
# # -----

# # GEOSGeometry ownership is retained by caller
# function prepareGeom(ptr::GEOSGeom)
#     result = GEOSPrepare(ptr)
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSPrepare")
#     end
#     result
# end

# destroyPreparedGeom(ptr::Ptr{GEOSPreparedGeometry}) = GEOSPreparedGeom_destroy(ptr)

# function prepcontains(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedContains(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedContains")
#     end
#     bool(result)
# end

# function prepcontainsproperly(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedContainsProperly(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedContainsProperly")
#     end
#     bool(result)
# end

# function prepcoveredby(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedCoveredBy(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedCoveredBy")
#     end
#     bool(result)
# end

# function prepcovers(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedCovers(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedCovers")
#     end
#     bool(result)
# end

# function prepcrosses(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedCrosses(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedCrosses")
#     end
#     bool(result)
# end

# function prepdisjoint(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedDisjoint(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedDisjoint")
#     end
#     bool(result)
# end

# function prepintersects(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedIntersects(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedIntersects")
#     end
#     bool(result)
# end

# function prepoverlaps(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedOverlaps(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedOverlaps")
#     end
#     bool(result)
# end

# function preptouches(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedTouches(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedTouches")
#     end
#     bool(result)
# end

# function prepwithin(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom)
#     result = GEOSPreparedWithin(g1, g2)
#     if result == 0x02
#         error("LibGEOS: Error in GEOSPreparedWithin")
#     end
#     bool(result)
# end

# # -----
# # STRtree functions
# # -----
# # GEOSSTRtree_create
# # GEOSSTRtree_insert
# # GEOSSTRtree_query
# # GEOSSTRtree_iterate
# # GEOSSTRtree_remove
# # GEOSSTRtree_destroy

# # -----
# # Unary predicate - return 2 on exception, 1 on true, 0 on false
# # -----
for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval isEmpty(obj::$geom) = isEmpty(obj.ptr)
    @eval isSimple(obj::$geom) = isSimple(obj.ptr)
    @eval isRing(obj::$geom) = isRing(obj.ptr)
    @eval isValid(obj::$geom) = isValid(obj.ptr)
    @eval hasZ(obj::$geom) = hasZ(obj.ptr)
end

isClosed(obj::LineString) = isClosed(obj.ptr) # Call only on LINESTRING

# # -----
# # Dimensionally Extended 9 Intersection Model related
# # -----

# # GEOSRelatePattern (return 2 on exception, 1 on true, 0 on false)
# # GEOSRelate (return NULL on exception, a string to GEOSFree otherwise)
# # GEOSRelatePatternMatch (return 2 on exception, 1 on true, 0 on false)
# # GEOSRelateBoundaryNodeRule (return NULL on exception, a string to GEOSFree otherwise)

# # -----
# # Validity checking -- return 2 on exception, 1 on true, 0 on false
# # -----

# # /* These are for use with GEOSisValidDetail (flags param) */
# # enum GEOSValidFlags {
# #     GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE=1
# # };

# # * return NULL on exception, a string to GEOSFree otherwise
# # GEOSisValidReason

# # * Caller has the responsibility to destroy 'reason' (GEOSFree)
# # * and 'location' (GEOSGeom_destroy) params
# # * return 2 on exception, 1 when valid, 0 when invalid
# # GEOSisValidDetail

# # -----
# # Geometry info
# # -----

# # May be called on all geometries in GEOS 3.x, returns -1 on error and 1
# # for non-multi geometries. Older GEOS versions only accept
# # GeometryCollections or Multi* geometries here, and are likely to crash
# # when fed simple geometries, so beware if you need compatibility with
# # old GEOS versions.
# function numGeometries(ptr::GEOSGeom)
#     result = GEOSGetNumGeometries(ptr)
#     if result == -1
#         error("LibGEOS: Error in GEOSGeomTypeId")
#     end
#     result
# end

# # Call only on GEOMETRYCOLLECTION or MULTI*
# # (Return a pointer to the internal Geometry. Return NULL on exception.)
# # Returned object is a pointer to internal storage:
# # it must NOT be destroyed directly.
# # Up to GEOS 3.2.0 the input geometry must be a Collection, in
# # later version it doesn't matter (i.e. getGeometryN(0) for a single will return the input).
# function getGeometry(ptr::GEOSGeom, n::Integer)
#     result = GEOSGetGeometryN(ptr, Int32(n-1))
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSGetGeometryN")
#     end
#     result
# end
# getGeometries(ptr::GEOSGeom) = GEOSGeom[getGeometry(ptr, i) for i=1:numGeometries(ptr)]

# Converts Geometry to normal form (or canonical form).
for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval Base.normalize!(obj::$geom) = normalize!(obj.ptr)
end

# # Return -1 on exception, Geometry must be a Point.
# let out = Array{Float64}(1)
#     global getGeomX
#     function getGeomX(ptr::GEOSGeom)
#         result = GEOSGeomGetX(ptr, pointer(out))
#         if result == -1
#             error("LibGEOS: Error in GEOSGeomGetX")
#         end
#         out[1]
#     end
# end

# let out = Array{Float64}(1)
#     global getGeomY
#     function getGeomY(ptr::GEOSGeom)
#         result = GEOSGeomGetY(ptr, pointer(out))
#         if result == -1
#             error("LibGEOS: Error in GEOSGeomGetY")
#         end
#         out[1]
#     end
# end

interiorRings(obj::Polygon) = map(LinearRing, interiorRings(obj.ptr))
exteriorRing(obj::Polygon) = LinearRing(exteriorRing(obj.ptr))

# # Return -1 on exception
# function numCoordinates(ptr::GEOSGeom)
#     result = GEOSGetNumCoordinates(ptr)
#     if result == -1
#         error("LibGEOS: Error in GEOSGetNumCoordinates")
#     end
#     result
# end

# # Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
# function getCoordSeq(ptr::GEOSGeom)
#     result = GEOSGeom_getCoordSeq(ptr)
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSGeom_getCoordSeq")
#     end
#     result
# end
# # getGeomCoordinates(ptr::GEOSGeom) = getCoordinates(getCoordSeq(ptr))

# # Return 0 on exception (or empty geometry)
# getGeomDimensions(ptr::GEOSGeom) = GEOSGeom_getDimensions(ptr)

# # Return 2 or 3.
# getCoordinateDimension(ptr::GEOSGeom) = int(GEOSGeom_getCoordinateDimension(ptr))

# # Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
# function getPoint(ptr::GEOSGeom, n::Integer)
#     result = GEOSGeomGetPointN(ptr, Int32(n-1))
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSGeomGetPointN")
#     end
#     result
# end

numPoints(obj::LineString) = numPoints(obj.ptr) # Call only on LINESTRING
startPoint(obj::LineString) = Point(startPoint(obj.ptr)) # Call only on LINESTRING
endPoint(obj::LineString) = Point(endPoint(obj.ptr)) # Call only on LINESTRING

# # -----
# # Misc functions
# # -----
for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval area(obj::$geom) = geomArea(obj.ptr)
    @eval geomLength(obj::$geom) = geomLength(obj.ptr)
end

for g1 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    for g2 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
        @eval distance(obj1::$g1, obj2::$g2) = geomDistance(obj1.ptr, obj2.ptr)
        @eval hausdorffdistance(obj1::$g1, obj2::$g2) = hausdorffdistance(obj1.ptr, obj2.ptr)
        @eval hausdorffdistance(obj1::$g1, obj2::$g2, densify::Real) = hausdorffdistance(obj1.ptr, obj2.ptr, densify)
    end
end

# Returns the closest points of the two geometries. The first point comes from g1 geometry and the second point comes from g2.
for g1 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    for g2 in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
        @eval begin
            function nearestPoints(obj1::$g1, obj2::$g2)
                points = nearestPoints(obj1.ptr, obj2.ptr)
                if points == C_NULL
                    return Point[]
                else
                    return Point[Point(getCoordinates(points,1)),Point(getCoordinates(points,2))]
                end
            end
        end
    end
end


text:
using MLDatasets
using CUDA
using Flux
using Flux: Data.DataLoader
using Flux: onehotbatch, onecold, @epochs
using Statistics: mean

x_train, y_train = MNIST.traindata()
x_test, y_test = MNIST.testdata()
x_train = reshape(x_train, 28, 28, 1, :) .|> Float32
x_test = reshape(x_test, 28, 28, 1, :) .|> Float32
# one hot batch
y_train = onehotbatch(y_train, 0:9)
y_test = onehotbatch(y_test, 0:9)
x_train, x_test, y_train, y_test = if CUDA.functional()
    x_train |> gpu, x_test |> gpu, y_train |> gpu, y_test |> gpu
else
    x_train, x_test, y_train, y_test
end
data = DataLoader((x_train, y_train), batchsize=40, shuffle=true)

model = Chain(
  Conv((5, 5), 1=>6, relu),
  MaxPool((2, 2)),
  Conv((5, 5), 6=>16, relu),
  MaxPool((2, 2)),
  flatten,
  Dense(256, 120, relu),
  Dense(120, 84, relu),
  Dense(84, 10),
  softmax,
)
model = if CUDA.functional() model |> gpu else model end

loss(x, y) = Flux.crossentropy(model(x), y)
accurate(x, y) = mean(onecold(model(x)) .== onecold(y))
opt = ADAM()

@epochs 10 begin
  Flux.train!(loss, params(model), data, opt)
  @show(loss(x_test, y_test))
  @show(accurate(x_test, y_test))
end


text:
using ApproxFun, SingularIntegralEquations, Base.Test

## quadratic

a=1+10*im;b=2-6*im
d=Curve(Fun(x->1+a*x+b*x^2))

x=Fun(d)
w=sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x))

@test_approx_eq cauchy(w,2.) (-4.722196879007759+2.347910413861846im)
@test_approx_eq cauchy(w,2.) sum(w/(x-2.))/(2π*im)
@test_approx_eq logkernel(w,2.) linesum(w*log(abs(x-2.)))/π

w=1/(sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x)))

@test_approx_eq cauchy(w,2.) sum(w/(x-2.))/(2π*im)
@test_approx_eq logkernel(w,2.) linesum(w*log(abs(x-2.)))/π
@test_approx_eq (SingularIntegral(0)*w)(fromcanonical(d,0.1)) logkernel(w,fromcanonical(d,0.1))
@test_approx_eq (Hilbert()*w)(fromcanonical(d,0.1)) hilbert(w,fromcanonical(d,0.1))

## cubic

a=1+10*im;b=2-6*im
d=Curve(Fun(x->1+a*x+x^2+b*x^3))

x=Fun(d)
w=sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x))

@test_approx_eq cauchy(w,2.) sum(w/(x-2.))/(2π*im)
@test_approx_eq logkernel(w,2.) linesum(w*log(abs(x-2.)))/π

w=1/(sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x)))

@test_approx_eq cauchy(w,2.) sum(w/(x-2.))/(2π*im)
@test_approx_eq logkernel(w,2.) linesum(w*log(abs(x-2.)))/π
@test_approx_eq (SingularIntegral(0)*w)(fromcanonical(d,0.1)) logkernel(w,fromcanonical(d,0.1))
@test_approx_eq (Hilbert()*w)(fromcanonical(d,0.1)) hilbert(w,fromcanonical(d,0.1))

## quartic

a=1+10*im;b=2-6*im
d=Curve(Fun(x->1+a*x+x^2+x^3+b*x^4))

x=Fun(d)
w=sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x))

@test_approx_eq cauchy(w,2.) sum(w/(x-2.))/(2π*im)
@test_approx_eq logkernel(w,2.) linesum(w*log(abs(x-2.)))/π

w=1/(sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x)))

@test_approx_eq cauchy(w,2.) sum(w/(x-2.))/(2π*im)
@test_approx_eq logkernel(w,2.) linesum(w*log(abs(x-2.)))/π
@test_approx_eq (SingularIntegral(0)*w)(fromcanonical(d,0.1)) logkernel(w,fromcanonical(d,0.1))
@test_approx_eq (Hilbert()*w)(fromcanonical(d,0.1)) hilbert(w,fromcanonical(d,0.1))

## Arc

d=exp(im*Interval(0.1,0.2))
x=Fun(d)
w=sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x))

z=10.;
@test_approx_eq sum(w/(x-z))/(2π*im) cauchy(w,z)
@test_approx_eq sum(w*log(z-x))/(-2π*im) cauchyintegral(w,z)
@test_approx_eq linesum(w*log(abs(z-x)))/π logkernel(w,z)

w=1/(sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x)))
@test_approx_eq sum(w/(x-z))/(2π*im) cauchy(w,z)
@test_approx_eq sum(w*log(z-x))/(-2π*im) cauchyintegral(w,z)
@test_approx_eq linesum(w*log(abs(z-x)))/π logkernel(w,z)
@test_approx_eq (SingularIntegral(0)*w)(fromcanonical(d,0.1)) logkernel(w,fromcanonical(d,0.1))
@test_approx_eq (Hilbert()*w)(fromcanonical(d,0.1)) hilbert(w,fromcanonical(d,0.1))


text:
# Autogenerated wrapper script for libidn_jll for x86_64-apple-darwin
export idn, libidn

JLLWrappers.@generate_wrapper_header("libidn")
JLLWrappers.@declare_executable_product(idn)
JLLWrappers.@declare_library_product(libidn, "@rpath/libidn.12.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        idn,
        "bin/idn",
    )

    JLLWrappers.@init_library_product(
        libidn,
        "lib/libidn.12.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

type IndexedVector{T}
    elts::Vector{T}
    nzidx::Vector{Int}
    nnz::Int
    empty::BitArray{1}
end

IndexedVector{T}(::Type{T},n::Integer) = IndexedVector(zeros(T,n),zeros(Int,n),0,trues(n))

function addelt!{T}(v::IndexedVector{T},i::Integer,val::T)
    if val != zero(T)
        if v.empty[i]  # new index
            v.elts[i] = val
            v.nzidx[v.nnz += 1] = i
            v.empty[i] = false
        else
            v.elts[i] += val
        end
    end
    return nothing
end

function Base.empty!{T}(v::IndexedVector{T})
    elts = v.elts
    nzidx = v.nzidx
    empty = v.empty
    for i in 1:v.nnz
        elts[nzidx[i]] = zero(T)
        empty[nzidx[i]] = true
    end
    v.nnz = 0
end

Base.length(v::IndexedVector) = length(v.elts)
function Base.resize!(v::IndexedVector, n::Integer)
    if n > length(v)
        @assert v.nnz == 0 # only resize empty vector
        resize!(v.elts, n)
        fill!(v.elts,0)
        resize!(v.nzidx, n)
        resize!(v.empty, n)
        fill!(v.empty,true)
    end
end

# lightweight unsafe view for vectors
# it seems that the only way to avoid triggering
# allocations is to have only bitstype fields, so
# we store a pointer.
immutable VectorView{T} <: DenseVector{T}
    offset::Int
    len::Int
    ptr::Ptr{T}
end

Base.getindex(v::VectorView,idx::Integer) = unsafe_load(v.ptr, idx+v.offset)
Base.setindex!(v::VectorView,value,idx::Integer) = unsafe_store!(v.ptr,value,idx+v.offset)
function Base.setindex!{T}(v::VectorView{T},value::T,idx::Vector{Int})
    for i in idx
        v[i] = value
    end
end
Base.length(v::VectorView) = v.len
function Base.fill!{T}(v::VectorView{T},value)
    val = convert(T,value)
    for i in 1:length(v)
        v[i] = val
    end
    nothing
end
function Base.scale!{T<:Number}(v::VectorView{T},value::T)
    for i in 1:length(v)
        v[i] *= value
    end
    nothing
end
function reinterpret_unsafe{T,R}(::Type{T},x::Vector{R})
    # how many T's fit into x?
    @assert isbits(T) && isbits(R)
    len = length(x)*sizeof(R)
    p = reinterpret(Ptr{T},pointer(x))
    return VectorView(0,div(len,sizeof(T)),p)
end


text:
using Test: @testset, @test
using Random: bitrand
using Statistics: mean
using LinearAlgebra: norm
import RestrictedBoltzmannMachines as RBMs
import CenteredRBMs
using CenteredRBMs: center, uncenter

@testset "pcd" begin
    rbm = CenteredRBMs.center(RBMs.BinaryRBM((28,28), 100))
    train_x = bitrand(28,28,1024)

    RBMs.initialize!(rbm, train_x) # fit independent site statistics and center
    @test CenteredRBMs.visible_offset(rbm) ≈ dropdims(mean(train_x; dims=3); dims=3)
    train_h = RBMs.mean_h_from_v(rbm, train_x)
    @test CenteredRBMs.hidden_offset(rbm) ≈ vec(mean(train_h; dims=2))
    @test norm(mean(RBMs.inputs_v_to_h(rbm, train_x); dims=2)) < 1e-6

    RBMs.pcd!(rbm, train_x; batchsize=128)
    @test CenteredRBMs.visible_offset(rbm) ≈ dropdims(mean(train_x; dims=3); dims=3)
    train_h = RBMs.mean_h_from_v(rbm, train_x)
    @test CenteredRBMs.hidden_offset(rbm) ≈ vec(mean(train_h; dims=2))
end


text:
#MDDatasets: meas interface
#Intent: use meas(:MEASTYPE, ...) to minimize namespace pollution.
#-------------------------------------------------------------------------------

meas(mtype::Symbol, args...; kwargs...) =
	meas(DS(mtype), args...; kwargs...)
meas(::DS{T}, args...; kwargs...) where T<:Symbol =
	throw(ArgumentError("Measurement type: meas($mtype, ...) not recognized."))

meas(::DS{:xcross}, args...; kwargs...) = xcross(args...; kwargs...)
meas(::DS{:ycross}, args...; kwargs...) = ycross(args...; kwargs...)
meas(::DS{:xcross1}, args...; kwargs...) = xcross1(args...; kwargs...)
meas(::DS{:ycross1}, args...; kwargs...) = ycross1(args...; kwargs...)

meas(::DS{:delay}, args...; kwargs...) = measdelay(args...; kwargs...)

#Last Line


text:
export
    homogenize,
    homoexponent_vectors,
    collect_homoexponent_vectors,
    pushfirst_sliceofones,
    Collectallrows,
    collect_coefffiltred_exponent_vectors,
    findallpositivevertices,
    findallnegativevertices,
    findallnonpositivevertices,
    findallnonnegativevertices,
    collect_termfiltred_exponent_vectors,
    collectpositivevertices,
    collectnegativevertices,
    collectnonpositivevertices,
    collectnonnegativevertices,
    Findallrows,
    Collectallcols,
    Findallcols,
    texponents,
    tcoeffs,
    realpositiveroots,
    tevalof,
    teval,
    tpointof,
    collect_realpositiveroots

## Code adapted from "enumarate(iter)"
struct Homogenize{I}
    itr::I
end

"""
    homogenize(iter)

Given an iterator `iter` whose values are 1d arrays,
returns an iterator that yields `[1; x]` where `1` is a unit of the
element type of `x`, and `x` is a value of `iter`.

# Examples
```jldoctest; setup = :(using CRNT)
julia> x = [[1, 2], [3, 4], [5, 6]];

julia> collect(homogenize(x))
3-element Array{Array{Int64,1},1}:
 [1, 1, 2]
 [1, 3, 4]
 [1, 5, 6]
```
"""
homogenize(iter) = Homogenize(iter)

Base.length(h::Homogenize) = Base.length(h.itr)
Base.size(h::Homogenize) = Base.size(h.itr)
Base.eltype(::Type{Homogenize{I}}) where {I} = eltype(I)
Base.IteratorSize(::Type{Homogenize{I}}) where {I} = Base.IteratorSize(I)
Base.IteratorEltype(::Type{Homogenize{I}}) where {I} = Base.IteratorEltype(I)

Base.@propagate_inbounds function Base.iterate(h::Homogenize, state)
    n = iterate(h.itr, state...)
    n === nothing && return n
    return pushfirst!(Vector(n[1]), oneunit(Base.eltype(n[1]))), n[2]
end

Base.@propagate_inbounds function Base.iterate(h::Homogenize)
    n = iterate(h.itr)
    n === nothing && return n
    return pushfirst!(Vector(n[1]), oneunit(Base.eltype(n[1]))), n[2]
end

"""
    homoexponent_vectors(p::MPolyElem)

An iterator for the homogenized (with an extra first component set to '1') exponent vectors of the multivariate polynomial `p`.
To retrieve an array, use `collect(homoexponent_vectors(p))`.
# Examples
```jldoctest; setup = :(using CRNT, Nemo)
julia> using Nemo

julia> R, vars = PolynomialRing(ZZ, vcat(["k\$i" for i in 1:5], ["x\$i" for i in 1:4]));

julia> p = vars[1]*vars[2]*vars[6]-vars[8]
k1*k2*x1-x3

julia> collect(exponent_vectors(p))
2-element Array{Array{Int64,1},1}:
 [1, 1, 0, 0, 0, 1, 0, 0, 0]
 [0, 0, 0, 0, 0, 0, 0, 1, 0]

julia> collect(homoexponent_vectors(p))
2-element Array{Any,1}:
 [1, 1, 1, 0, 0, 0, 1, 0, 0, 0]
 [1, 0, 0, 0, 0, 0, 0, 0, 1, 0]
```
"""
homoexponent_vectors(p::MPolyElem) = homogenize(exponent_vectors(p))

# function collect_homoiter(iter)
#     ## First iteration outside the loop to prealocate v
#     (u, state) = iterate(iter)
#     ## Prealocate v
#     v = ones(eltype(u), size(u,1)+1, length(iter))
#     # v = Matrix{eltype(u)}(undef, size(u,1)+1, length(iter))
#     # v[1,:] .= one(eltype(u))
#     v[2:end,1] = u
#     next = iterate(iter, state)
#     ## Start the loop
#     while next !== nothing
#         (u, state) = next
#         v[2:end,state-1] = u
#         next = iterate(iter, state)
#     end
#     return transpose(v)
# end

"""
    CollectExponent_homovectors(p::MPolyElem)

Collects the homogenized (with an extra first component set to 1) exponents of `p` in the rows of a matrix.
The returned matrix can be used to compute the Newton polytope of `p`.
# Examples
```jldoctest; setup = :(using CRNT, Nemo)
julia> using Nemo

julia> R, vars = PolynomialRing(ZZ, vcat(["k\$i" for i in 1:5], ["x\$i" for i in 1:4]));

julia> q = vars[1]*vars[2]*vars[6]-vars[8]
k1*k2*x1-x3

julia> collect(exponent_vectors(q))
2-element Array{Array{Int64,1},1}:
 [1, 1, 0, 0, 0, 1, 0, 0, 0]
 [0, 0, 0, 0, 0, 0, 0, 1, 0]

julia> CollectExponent_homovectors(q)
2×10 Transpose{Int64,Array{Int64,2}}:
 1  1  1  0  0  0  1  0  0  0
 1  0  0  0  0  0  0  0  1  0
```

Recall it is approximately ten times faster than the equivalent
    transpose(reduce(hcat,collect(exponent_homovectors(q))))
"""
function collect_homoexponent_vectors(p)
    # ## Save the iterator (this function just works for iters whose states are Int and start at ZERO)
    # iter = exponent_vectors(p)
    # ## First iteration outside the loop to prealocate v
    # (u, state) = iterate(iter)
    # ## Prealocate v
    # v = ones(eltype(u), size(u, 1) + 1, length(iter))
    # v[2:end,1] = u
    # next = iterate(iter, state)
    # ## Start the loop
    # while next !== nothing
    #     (u, state) = next
    #     v[2:end,state] = u
    #     next = iterate(iter, state)
    # end
    # return v
    return pushfirst_sliceofones(reduce(hcat, exponent_vectors(p)))
end

function pushfirst_sliceofones(M, dims=1)
    indx = collect(size(M))
    indx[dims] = 1
    return cat(ones(eltype(M), Tuple(indx)), M, dims=dims)
end

# R, vars = PolynomialRing(ZZ, vcat(["k$i" for i in 1:5], ["x$i" for i in 1:4]));
# """
#     Collectallrows(predicate, p::MPolyElem, V::AbstractMatrix)

# Returns an `Array` of all the exponents of terms of `p` for which the exponent is a row of `V` and the coefficient satisfies `predicate`.

# # Examples
# ```jldoctest; setup = :(using CRNT)
# julia> using Nemo

# julia> R, vars = PolynomialRing(ZZ, vcat(["k\$i" for i in 1:5], ["x\$i" for i in 1:4]));

# julia> q = vars[1]*vars[2]*vars[6]-vars[3]*vars[8]+vars[5]*vars[9]-2*vars[1]*vars[3]*vars[5]^2
# k1*k2*x1-2*k1*k3*k5^2-k3*x3+k5*x4

# julia> collect(exponent_vectors(q))
# 4-element Array{Array{Int64,1},1}:
#  [1, 1, 0, 0, 0, 1, 0, 0, 0]
#  [1, 0, 1, 0, 2, 0, 0, 0, 0]
#  [0, 0, 1, 0, 0, 0, 0, 1, 0]
#  [0, 0, 0, 0, 1, 0, 0, 0, 1]

# julia> V =  [[1 1 0 0 0 1 0 0 0]; [0 0 1 0 0 0 0 1 0]; [1 0 1 0 2 0 0 0 0]];

# julia> Collectallrows(x->Nemo.isless(x,0), q, V)
# 2-element Array{Array{Int64,1},1}:
#  [1, 0, 1, 0, 2, 0, 0, 0, 0]
#  [0, 0, 1, 0, 0, 0, 0, 1, 0]
# ```
# """
function Collectallrows(predicate, p::MPolyElem, V::AbstractMatrix)
    predisinV(x) = predicate(x[1]) && x[2] in eachrow(V)
    getexponent(x) = x[2]
    return getexponent.(Iterators.filter(predisinV, zip(coeffs(p), exponent_vectors(p))))
end

## Optmized functions to find which rows in matrix vertices correspond to pos/neg/nonpos/nonneg exponents
function collect_coefffiltred_exponent_vectors(predicate, p::MPolyElem)
    return [iterate(exponent_vectors(p),n-1)[1] for n in findallcoeffs(predicate, p)]
end

function findallpositivevertices(p::MPolyElem, vertices::AbstractMatrix)
    return findall(in(collect_coefffiltred_exponent_vectors(ispositive, p)), collect(eachrow(vertices)))
end

function findallnegativevertices(p::MPolyElem, vertices::AbstractMatrix)
    return findall(in(collect_coefffiltred_exponent_vectors(isnegative, p)), collect(eachrow(vertices)))
end

function findallnonpositivevertices(p::MPolyElem, vertices::AbstractMatrix)
    return findall(in(collect_coefffiltred_exponent_vectors(isnonpositive, p)), collect(eachrow(vertices)))
end

function findallnonnegativevertices(p::MPolyElem, vertices::AbstractMatrix)
    return findall(in(collect_coefffiltred_exponent_vectors(isnonnegative, p)), collect(eachrow(vertices)))
end

## Optimized functions to collect rows of matrix vertices correspond to pos/neg/nonpos/nonneg exponents
function collect_termfiltred_exponent_vectors(predicate, p::MPolyElem)
    return getexponent.(filterterms(predicate, p))
end

function isours(predicate, vertices)
    return term -> predicate(getcoeff(term)) && getexponent(term) in eachrow(vertices)
end

function collectpositivevertices(p::MPolyElem, vertices::AbstractMatrix)
    return collect_termfiltred_exponent_vectors(isours(ispositive, vertices), p)
end

function collectnegativevertices(p::MPolyElem, vertices::AbstractMatrix)
    return collect_termfiltred_exponent_vectors(isours(isnegative, vertices), p)
end

function collectnonpositivevertices(p::MPolyElem, vertices::AbstractMatrix)
    return collect_termfiltred_exponent_vectors(isours(isnonpositive, vertices), p)
end

function collectnonnegativevertices(p::MPolyElem, vertices::AbstractMatrix)
    return collect_termfiltred_exponent_vectors(isours(isnonnegative, vertices), p)
end

function collectpositivevertices(p::MPolyElem)
    return collectpositivevertices(p, verticesof(Newtonpolytope(p)))
end

function collectnegativevertices(p::MPolyElem)
    return collectnegativevertices(p, verticesof(Newtonpolytope(p)))
end

function collectnonpositivevertices(p::MPolyElem)
    return collectnonpositivevertices(p, verticesof(Newtonpolytope(p)))
end

function collectnonnegativevertices(p::MPolyElem)
    return collectnonnegativevertices(p, verticesof(Newtonpolytope(p)))
end

function Findallrows(predicate, p::MPolyElem, V::AbstractMatrix)
    predexpposs = findall(predicate, collect(coeffs(p)))
    predexp = [iterate(exponent_vectors(p), n - 1)[1] for n in predexpposs]
    return findall(in(predexp), collect(eachrow(V)))
end

function Collectallcols(predicate, p::MPolyElem, V::AbstractMatrix)
    predisinV(x) = predicate(x[1]) && x[2] in eachcol(V)
    getexponent(x) = x[2]
    return getexponent.(Iterators.filter(predisinV, zip(coeffs(p), exponent_vectors(p))))
end

function Findallcols(predicate, p::MPolyElem, V::AbstractMatrix)
    predexpposs = findall(predicate, collect(coeffs(p)))
    predexp = [iterate(exponent_vectors(p), n - 1)[1] for n in predexpposs]
    return findall(in(predexp), collect(eachcol(V)))
end

# getcoeff(x) = x[1]
# isnegative(x) = x < 0
# iscoeffnegative(x) = isnegative(getcoeff(x))
# negcoeffs_exponent(p) = Iterators.Filter(iscoeffnegative, zip(coeffs(p),exponent_vectors(p)))
# negexp = getexponent.(negcoeffs_exponent(p))

## Method to obtain negexp 4 times faster for the non-commented function
## (you need to define the polynomial q)
# using BenchmarkTools
# @benchmark getexponent.(negcoeffs_exponent(q))
# @benchmark [iterate(exponent_vectors(q),n-1)[1] for n in findall(isnegative, collect(coeffs(q)))]

function texponents(p::MPolyElem, texp)
    ## Compute exponent of t (to know the minimum in case it is negative).
    return [dot(texp, exp) for exp in exponent_vectors(p)]
end

function tcoeffs(p::MPolyElem, texps)
    maxdeg = maximum(texps)
    mindeg = minimum(texps)
    ## If mindeg>0 we divide by t^(mindeg), otherwise we multiply by it.
    ## Both correspond to shift by -mindeg!
    polyn = zeros(maxdeg - mindeg + 1)
    for (c, i) in zip(coeffs(p), texps)
        polyn[i - mindeg + 1] += BigFloat(c)
    end
    return polyn
end

function tcoeffs(coeffs::AbstractArray, texps)
    maxdeg = maximum(texps)
    mindeg = minimum(texps)
    ## If mindeg>0 we divide by t^(mindeg), otherwise we multiply by it.
    ## Both correspond to shift by -mindeg!
    polyn = zeros(maxdeg - mindeg + 1)
    for (c, i) in zip(coeffs, texps)
        polyn[i - mindeg + 1] += BigFloat(c)
    end
    return polyn
end

function realpositiveroots(polyn, rtol::Real=1e-7)
    # poly = evalexponent(p, texp)
    roots = PolynomialRoots.roots(polyn)
    return filter(ispositive, filter_complex(roots, rtol))
end

function tevalof(p::MPolyElem, texps, T=Float64)
    return t -> sum([T(c) * (t^i) for (c, i) in zip(coeffs(p), texps)])
end

function tevalof(polyn::AbstractVector{T}) where {T <: AbstractFloat}
    return t -> sum([T(c) * (t^(i - 1)) for (i, c) in enumerate(polyn)])
end

function teval(p::MPolyElem, point::AbstractVector{T}) where {T <: AbstractFloat}
    return sum([T(c) * prod((point).^(exp)) for (c, exp) in dissect(p)])
end

function tpointof(texp)
    return t -> (t).^(texp)
end

function collect_realpositiveroots(p::MPolyElem, texp, rtol::Real=1e-7)
    texps = texponents(p, texp)
    tpoly = tcoeffs(p, texps)
    troots = realpositiveroots(tpoly, rtol)
    # return tpointof(texp).(troots)
    return [texps, tpoly, troots]
end

function collect_realpositiveroots(exps, coeffs, texp, rtol::Real=1e-7)
    texps = [dot(texp, exp) for exp in exps]
    tpoly = tcoeffs(coeffs, texps)
    troots = realpositiveroots(tpoly, rtol)
    # return tpointof(texp).(troots)
    return [texps, tpoly, troots]
end


text:
"""
The module provides example calculations for the second edition
of the book

> Kriele M. and Wolf, J.
> _Wertorientiertes Risikomanagement von Versicherungsunternehmen_,
>  2nd edition, Springer-Verlag, Berlin Heidelberg,  2016
> (to be published)

It is also intended to use these examples for future editions of
the English  translation of this book,
_Value-Oriented Risk Management of Insurance Companies_.
(The examples in both the first German edition and the first
 English edition are written in R).

The module consists of 4 distinct parts:

  - *SSTLife*: An extremely simplified example of the SST
(Swiss Solvency Test) calculation for life insurance.
The Swiss Solvency Test is the Swiss regulatory capital
requirement.  The resulting monetary requirement is referred
to as the "target capital" `ZK`.
  - *S2Life*: A simplified example of the S2 (Solvency 2)
calculation for non-life insurance. Solvency 2 is the new
regulatory capital requirement in the European Union. The
resulting monetary requirement is referred to as the
"Solvency capital requirement" `SCR`.
  - *S2NonLife*: A simplified example of the S2 calculation
for life insurance
  - *ECModel*: An extremely simplified example of an internal
economic capital model for non-life insurance. This model is
used to illustrate some techniques used in value based
management.

Note that we have simplified and (in part changed for our
exposition) the regulatory requirements for SST and Solvency 2.
Also note that the implementation of Solvency 2 may be slightly
different in different EU countries. For definitive information
about SST or Solvency 2, please consult the original literature
and any guidance issued by the supervisory authorities in the
jurisdiction of interest.
"""
module ValueOrientedRiskManagementInsurance

export es

using Distributions
using DataFrames#, DataArrays
using LinearAlgebra
using Random

# import Base.show, Base.isequal
# import Base.merge!
import Distributions.rand
import LinearAlgebra.⋅

## General functions --------------------------------------------
## Expected shortfall
es(x, α) =
  mean(sort(x, rev = true)[1:ceil(Integer, (1 - α) * length(x))])

# Simplfied Swiss Solvency Test----------------------------------
include("SST/SST__Types.jl")
include("SST/SST_Functions.jl")

# Simplified Life insurer ---------------------------------------
include("Life/Life__Types.jl")
include("Life/Life_Constructors.jl")
include("Life/Life_Functions.jl")

# Simplified Solvency 2 Life ------------------------------------
include("S2Life/S2Life__Types.jl")
include("S2Life/S2Life_Constructors.jl")
include("S2Life/S2Life_Functions.jl")

# Simplified Solvency 2 Non-Life --------------------------------
include("S2NonLife/S2NonLife__Types.jl")
include("S2NonLife/S2NonLife_Functions.jl")

# Simple economic capital model ---------------------------------
include("ECModel/ECModel__Types.jl")
include("ECModel/ECModel_Functions.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ExpressRouteServiceProviderPropertiesFormat <: SwaggerModel
    peeringLocations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: peeringLocations
    bandwidthsOffered::Any # spec type: Union{ Nothing, Vector{ExpressRouteServiceProviderBandwidthsOffered} } # spec name: bandwidthsOffered
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ExpressRouteServiceProviderPropertiesFormat(;peeringLocations=nothing, bandwidthsOffered=nothing, provisioningState=nothing)
        o = new()
        validate_property(ExpressRouteServiceProviderPropertiesFormat, Symbol("peeringLocations"), peeringLocations)
        setfield!(o, Symbol("peeringLocations"), peeringLocations)
        validate_property(ExpressRouteServiceProviderPropertiesFormat, Symbol("bandwidthsOffered"), bandwidthsOffered)
        setfield!(o, Symbol("bandwidthsOffered"), bandwidthsOffered)
        validate_property(ExpressRouteServiceProviderPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ExpressRouteServiceProviderPropertiesFormat

const _property_map_ExpressRouteServiceProviderPropertiesFormat = Dict{Symbol,Symbol}(Symbol("peeringLocations")=>Symbol("peeringLocations"), Symbol("bandwidthsOffered")=>Symbol("bandwidthsOffered"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ExpressRouteServiceProviderPropertiesFormat = Dict{Symbol,String}(Symbol("peeringLocations")=>"Vector{String}", Symbol("bandwidthsOffered")=>"Vector{ExpressRouteServiceProviderBandwidthsOffered}", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ExpressRouteServiceProviderPropertiesFormat }) = collect(keys(_property_map_ExpressRouteServiceProviderPropertiesFormat))
Swagger.property_type(::Type{ ExpressRouteServiceProviderPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ExpressRouteServiceProviderPropertiesFormat[name]))}
Swagger.field_name(::Type{ ExpressRouteServiceProviderPropertiesFormat }, property_name::Symbol) =  _property_map_ExpressRouteServiceProviderPropertiesFormat[property_name]

function check_required(o::ExpressRouteServiceProviderPropertiesFormat)
    true
end

function validate_property(::Type{ ExpressRouteServiceProviderPropertiesFormat }, name::Symbol, val)
end


text:
import QPDAS
#import OSQP

#Type to store n linear equalities + 1 pair of linear equalities
mutable struct SavedPlanes{T}
    A::Array{T,2}
    b::Array{T,1}
    n::Int64
    neq::Int64
    nineq::Int64
end

function projectonnormals!(s::SavedPlanes{T},x,y) where T
    n = length(x)
    m = length(s.b)

    eqidx =               1:(s.neq*(s.n+1))
    ineqidx = (s.neq*(s.n+1)+1):(s.neq+s.nineq)*(s.n+1)

    A = s.A[  eqidx,:]
    b = s.b[eqidx]
    C = s.A[ineqidx,:]
    d = s.b[ineqidx]
    z = -x
    #
    qp = QPDAS.QuadraticProgram(BigFloat.(A), BigFloat.(b), BigFloat.(-C), BigFloat.(-d), BigFloat.(z), I, scaling=true, ϵ=1e-12)
    sol, val = QPDAS.solve!(qp)
    y .= sol

    # model = OSQP.Model()
    # M = [A;-C]
    # u = [b;-d]
    # l = [b;fill(-Inf, length(d))]
    # OSQP.setup!(model; P=SparseMatrixCSC{Float64}(I, n, n), l=l, A=sparse(M), u=u, verbose=false,
    #     eps_abs=0.001*eps(), eps_rel=0.001*eps(),
    #     eps_prim_inf=0.001*eps(), eps_dual_inf=0.001*eps(), max_iter=10000)
    #
    # OSQP.update!(model; q=z)
    # results = OSQP.solve!(model)
    # y .= results.x

    # println(size(s.A))
    # println("A:")
    # println(s.A[  eqidx,:])
    # println("b:")
    # println(s.b[  eqidx])
    # println("C:")
    # println(s.A[ineqidx,:])
    # println("d:")
    # println(s.b[ineqidx])
    # println("x:")
    # println(x)
    return false
end

function SavedPlanes(x::AbstractVector{T}, n::Int, neq, nineq) where T
    total = (n+1)*neq+(n+1)*nineq #Save all eq and last ineq
    SavedPlanes{T}(similar(x,total, length(x)), similar(x,total), n, neq, nineq)
end

#
# #Type to store n linear equalities + 1 pair of linear equalities
# mutable struct SavedPlanes{T}
#     A::Array{T,2}
#     b::Array{T,1}
#     n::Int64
#     neq::Int64
#     nineq::Int64
# end
#
# function projectonnormals!(s::SavedPlanes{T},x,y) where T
#     n = length(x)
#     m = length(s.b)
#     #println(s.A)
#     #println(s.b)
#     AAt = s.A*transpose(s.A)
#     AAtF = try
#         bkfact!(AAt)
#     catch
#         y .= x
#         #println("Non symmetric")
#         return true
#     end
#     tmp = copy(s.b)
#     # tmp = -Ax+b
#     LinearAlgebra.gemv!('N', -one(T), s.A, x, one(T), tmp)
#     # tmp = (AA')⁻¹*(-Ax+b)
#     try
#         ldiv!(AAtF,tmp)
#     catch
#         y .= x
#         println("Noninvertable")
#         return true
#     end
#     # y = A'(AA')⁻¹*(-Ax+b)
#     mul!(y, transpose(s.A), tmp)
#     # y = x + A'(AA')⁻¹*(-Ax+b) = (I-A'(AA')⁻¹*A)x + A'(AA')⁻¹b
#     y .= y .+ x
#     #println("in 2")
#     return false
# end
#
# function SavedPlanes(x::AbstractVector{T}, n::Int, neq, nineq) where T
#     total = (n+1)*neq+nineq #Save all eq and last ineq
#     SavedPlanes{T}(similar(x,total, length(x)), similar(x,total), n, neq, nineq)
# end

#Add planes to a specific location
function addplanesat(v1,b1,i,s::SavedPlanes)
    s.A[i,:] .= v1
    s.b[i]    = b1
end

#Add planes to a random location in 1:n with probability p=0.05
function addplanesrand(v1,b1,s::SavedPlanes, p = 0.05)
    if s.n > 0 && rand() < p
        i = rand(1:s.n)
        addplanesat(v1,b1,i,s)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    slice(object, xmin:xmax, ymin:ymax, ...)

Slice spatial `object` using real coordinate ranges
`xmin:xmax`, `ymin:ymax`, ...
"""
slice(object, ranges::Vararg) =
  inside(object, Rectangle(first.(ranges), last.(ranges)))


text:
function Routes(pairs::Pair...)
    return Routes([pairs...])
end

# Priorities, so that e.g. r".*" doesn't catch absolut matches by e.g a string
pattern_priority(x::Pair) = pattern_priority(x[1])
pattern_priority(x::String) = 1
pattern_priority(x::Tuple) = 2
pattern_priority(x::Regex) = 3

function Base.setindex!(routes::Routes, f, pattern)
    idx = findfirst(pair-> pair[1] == pattern, routes.table)
    if idx !== nothing
        routes.table[idx] = pattern => f
    else
        push!(routes.table, pattern => f)
    end
    # Sort for priority so that exact string matches come first
    sort!(routes.table, by = pattern_priority)
    # return if it was inside already!
    return idx === nothing
end

function route!(application::Server, pattern_f::Pair)
    application.routes[pattern_f[1]] = pattern_f[2]
end

function route!(f, application::Server, pattern)
    route!(application, pattern => f)
end

function websocket_route!(application::Server, pattern_f::Pair)
    application.websocket_routes[pattern_f[1]] = pattern_f[2]
end

function apply_handler(f, args...)
    return f(args...)
end

function apply_handler(chain::Tuple, context, args...)
    f = first(chain)
    result = f(args...)
    return apply_handler(Base.tail(chain), context, result...)
end

function apply_handler(app::App, context)
    server = context.application
    session = insert_session!(server)
    html_dom = Base.invokelatest(app.handler, session, context.request)
    return html(page_html(session, html_dom))
end

function delegate(routes::Routes, application, request::Request, args...)
    for (pattern, f) in routes.table
        match = match_request(pattern, request)
        if match !== nothing
            context = (
                routes = routes,
                application = application,
                request = request,
                match = match
            )
            return apply_handler(f, context, args...)
        end
    end
    # If no route is found we have a classic case of 404!
    # What a classic this response!
    return response_404("Didn't find route for $(request.target)")
end

function match_request(pattern::String, request)
    return request.target == pattern ? pattern : nothing
end

function match_request(pattern::Regex, request)
    return match(pattern, request.target)
end

"""
    local_url(server::Server, url)
The local url to reach the server, on the server
"""
function local_url(server::Server, url)
    return string("http://", server.url, ":", server.port, url)
end

"""
    online_url(server::Server, url)
The url to connect to the server from the internet.
Needs to have `JSSERVE_CONFIGURATION.external_url` set to the IP or dns route of the server
"""
function online_url(server::Server, url)
    base_url = JSSERVE_CONFIGURATION.external_url[]
    if isempty(base_url)
        local_url(server, url)
    else
        base_url * url
    end
end

function websocket_request()
    headers = [
        "Host" => "localhost",
        "User-Agent" => "Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:68.0) Gecko/20100101 Firefox/68.0",
        "Accept" => "*/*",
        "Accept-Encoding" => "gzip, deflate, br",
        "Accept-Language" => "de,en-US;q=0.7,en;q=0.3",
        "Cache-Control" => "no-cache",
        "Connection" => "keep-alive, Upgrade",
        "Dnt" => "1",
        "Origin" => "https://localhost",
        "Pragma" => "no-cache",
        "Sec-Websocket-Extensions" => "permessage-deflate",
        "Sec-Websocket-Key" => "BL3d8I8KC5faPjubRM0riA==",
        "Sec-Websocket-Version" => "13",
        "Upgrade" => "websocket",
    ]
    msg = HTTP.Request(
        "GET",
        "/",
        headers,
        UInt8[],
        parent = nothing,
        version = v"1.1.0"
    )
    return Stream(msg, IOBuffer())
end

"""
warmup(application::Server)

Warms up the application, by sending a couple of request.
"""
function warmup(application::Server)
    yield() # yield to server task to give it a chance to get started
    task = application.server_task[]
    if Base.istaskdone(task)
        error("Webserver doesn't serve! Error: $(fetch(task))")
    end
    # Make a websocket request
    stream = websocket_request()
    try
        @async stream_handler(application, stream)
        write(stream, "blaaa")
    catch e
        # TODO make it not error so we can test this properly
        # This will error, since its not a propper websocket request
        @debug "Error in stream_handler" exception=e
    end
    target = register_local_file(JSServeLib.assets[1].local_path) # http target part
    asset_url = local_url(application, target)
    request = Request("GET", target)

    delegate(application.routes, application, request)

    if Base.istaskdone(task)
        error("Webserver doesn't serve! Error: $(fetch(task))")
    end
    resp = WebSockets.HTTP.get(asset_url, readtimeout=500, retries=1)
    if resp.status != 200
        error("Webserver didn't start succesfully")
    end
    return
end

function stream_handler(application::Server, stream::Stream)
    if WebSockets.is_upgrade(stream)
        try
            WebSockets.upgrade(stream) do request, websocket
                delegate(
                    application.websocket_routes, application, request, websocket
                )
            end
            return
        catch e
            # When browser closes we may get an io error here!
            if !(e isa Base.IOError)
                rethrow(e)
            end
        end
    end
    http_handler = HTTP.RequestHandlerFunction() do request
        delegate(
            application.routes, application, request,
        )
    end
    try
        HTTP.handle(http_handler, stream)
    catch e
        # we expect the IOError to happen, if either the page gets closed
        # or we close the server!
        if !(e isa Base.IOError && e.msg == "stream is closed or unusable")
            rethrow(e)
        end
    end
end

const MATCH_HEX = r"[\da-f]"
const MATCH_UUID4 = MATCH_HEX^8 * r"-" * (MATCH_HEX^4 * r"-")^3 * MATCH_HEX^12
const MATCH_SESSION_ID = MATCH_UUID4 * r"/" * MATCH_HEX^4

"""
Server(
        dom, url::String, port::Int;
        verbose = false
    )

Creates an application that manages the global server state!
"""
function Server(
        app::App, url::String, port::Int;
        verbose = false,
        routes = Routes(
            "/" => app,
            r"/assetserver/" * MATCH_HEX^40 * r"-.*" => file_server,
            r".*" => (context)-> response_404()
        ),
        websocket_routes = Routes(
            r"/" * MATCH_SESSION_ID => websocket_handler
        )
    )

    application = Server(
        url, port, Dict{String, Dict{String, Session}}(),
        Ref{Task}(), Ref{TCPServer}(),
        routes,
        websocket_routes
    )

    try
        start(application; verbose=verbose)
        # warmup server!
        warmup(application)
    catch e
        close(application)
        rethrow(e)
    end

    return application
end

function Server(
        app, url::String, port::Int; kw...
    )
    Server(App(app), url, port; kw...)
end

function isrunning(application::Server)
    return (isassigned(application.server_task) &&
        isassigned(application.server_connection) &&
        !istaskdone(application.server_task[]) &&
        isopen(application.server_connection[]))
end

function Base.close(application::Server)
    # Closing the io connection should shut down the HTTP listen loop
    for (id, session) in application.sessions
        close(session)
    end

    if isassigned(application.server_connection)
        close(application.server_connection[])
        @assert !isopen(application.server_connection[])
    end
    # For good measures, wait until the task finishes!
    if isassigned(application.server_task)
        try
            wait(application.server_task[])
            @assert !Base.isdone(application.server_connection[])
        catch e
            @debug "Server task failed with an (expected) exception on close" exception=e
        end
    end
    # Sometimes, the first request after closing will still go through
    # see: https://github.com/JuliaWeb/HTTP.jl/pull/494
    # We need to make sure, that we are the ones making this request,
    # So that a newly opened connection won't get a faulty response from this server!
    try
        app_url = local_url(application, "/")
        while true
            x = HTTP.get(app_url, readtimeout=1, retries=1)
            x.status != 200 && break # we got a bad request, maining server is closed!
        end
    catch e
        # This is expected to fail!
        @debug "Failed get request successfully after closing server!" exception=e
    end
    @assert !isrunning(application)
end

function start(application::Server; verbose=false)
    isrunning(application) && return
    address = Sockets.InetAddr(Sockets.getaddrinfo(application.url), application.port)
    ioserver = Sockets.listen(address)
    application.server_connection[] = ioserver
    # pass tcp connection to listen, so that we can close the server
    application.server_task[] = @async HTTP.listen(
            application.url, application.port; server=ioserver, verbose=verbose
        ) do stream::Stream
        Base.invokelatest(stream_handler, application, stream)
    end
    return
end

const GLOBAL_SERVER = Ref{Server}()

function get_server()
    if !isassigned(GLOBAL_SERVER) || istaskdone(GLOBAL_SERVER[].server_task[])
        GLOBAL_SERVER[] = Server(
            App("Nothing to see"),
            JSSERVE_CONFIGURATION.listen_url[],
            JSSERVE_CONFIGURATION.listen_port[],
            verbose=JSSERVE_CONFIGURATION.verbose[]
        )
    end
    return GLOBAL_SERVER[]
end

function insert_session!(server::Server, session=Session())
    server.sessions[session.id] = session
    return session
end


text:
import Base: (+), (-), (*), isnan, isfinite

#####################################################
# add basic maths to the point type
#####################################################

function add_maths{geodesy_type}(::Type{geodesy_type}, fields)

    nf = length(fields)  # number of fields and their names for the Geodesy variables

    # allow interaction with these types
    interact_types = [Vec{nf, Float64}, Vector]

    # create an expression to add to
    qb = quote end

    # Allow subtracting the same Geodesy representation
    rhs_expr = :(())
    append!(rhs_expr.args, [:(X.$(field) - Y.$(field)) for field in fields[1:nf]])

    # form the whole expression
    qn = quote
        (-){T <: $(geodesy_type)}(X::T, Y::T) = Vec{$(nf), Float64}($(rhs_expr.args...))
    end
    append!(qb.args, qn.args)

    # Allow Rotations
    if (nf == 3)
        qn = quote
            (*)(R::Mat{3,3,Float64}, X::$(geodesy_type)) = R * Vec(X)
        end
        append!(qb.args, qn.args)  # include them
    end

    # Allow uniform scaling
    # build the rhs epxression
    rhs_expr = :(())
    append!(rhs_expr.args, [:(s * X.$(field)) for field in fields[1:nf]])

    # form the whole expression
    qn = quote
        (*){T <: $(geodesy_type), U <: Real}(X::T, s::U) = T($(rhs_expr.args...))
        (*){T <: $(geodesy_type), U <: Real}(s::U, X::T) = T($(rhs_expr.args...))
    end
    append!(qb.args, qn.args)

    # NaN check
    rhs_nan_expr = Expr(:||, :(isnan(X.$(fields[nf - 1]))), :(isnan(X.$(fields[nf]))))
    for i in nf-2:-1:1
        rhs_nan_expr = Expr(:||,  :(isnan(X.$(fields[i]))) , rhs_nan_expr)
    end

    # finite check
    rhs_finite_expr = Expr(:&&, :(isfinite(X.$(fields[nf - 1]))), :(isfinite(X.$(fields[nf]))))
    for i in nf-2:-1:1
        rhs_finite_expr = Expr(:&&,  :(isfinite(X.$(fields[i]))), rhs_finite_expr)
    end

    #  Fill function for the interactions types
    for iT in interact_types

        # add a vector to it
        rhs_add_expr = :(())
        append!(rhs_add_expr.args, [:(X.$(field) + dX[$(i)]) for (i, field) in enumerate(fields)])

        # subtract a vector from it
        rhs_sub_expr = :(())
        append!(rhs_sub_expr.args, [:(X.$(field) - dX[$(i)]) for (i, field) in enumerate(fields)])

        qn = quote

            # addition
            (+){T <: $(geodesy_type)}(X::T, dX::$(iT)) = T($(rhs_add_expr.args...))
            (+){T <: $(geodesy_type)}(dX::$(iT), X::T) = T($(rhs_add_expr.args...))

            # subtraction
            (-){T <: $(geodesy_type)}(X::T, dX::$(iT)) = T($(rhs_sub_expr.args...))

            # NaN
            isnan(X::$(geodesy_type)) = $(rhs_nan_expr)

            # finite
            isfinite(X::$(geodesy_type)) = $(rhs_finite_expr)

        end
        append!(qb.args, qn.args)  # include them
    end
    return qb

end


text:
module HTLogger

using SerialPorts
using Dates

struct TimeoutError <: Exception end

"""
    run(;path="log", baudrate=9600, port="", interval=5, lines_per_file=135_000, debug=false)

Starts logging the temperature and humidity data from the connected device.

# Keyword Arguments

* `path`: the path where the `.log` files are stored. If the path does not
  exist, it is created.

* `baudrate`: the baudrate the device is configured to

* `port`: the serial port the device is connected to. If `port` is an empty
  string, the program will try to find the device automatically.

* `interval`: the time between the individual measurements in seconds. It is
  implemented in such a way that the main measurement loop is idle for the
  `interval` time. This means that the individual measurements are more than
  `interval` seconds apart because additional time is spent reading the data
  from the device and writing to the file.

* `lines_per_file=135_000`: number of lines in a log file until a new file is
  created. 1MB roughly corresponds to 28,000 entries.

* `debug`: prints messages that help to debug in case there are problems
"""
function run(;path="log", baudrate=57600, port="",interval=5,lines_per_file = 135_000,debug=false)
    # In case an error occurs we will rerun the program recursively indefinetly
    # until the program is interrupted. Collect the keyword arguments that get
    # passed to the run function recursively. Do not specify the port in the
    # rerun because the logger port could have been changed.
    kwargs = (
        path=path,
        baudrate=baudrate,
        interval=interval,
        lines_per_file=lines_per_file,
        debug=debug
    )

    # define the serial port as a constant reference so it is visible in all
    # try-catch blocks
    local s = Ref{SerialPort}()

    try
        # if no port is supplied search automatically
        if port == ""
            port = find_port(; debug, baudrate)
        end

        # connect to the device
        s[] = SerialPort(port, baudrate)

        # wait for board to be ready
        sleep(3)
    catch err
        if err isa InterruptException
            # in case the user interrupts the program
            println("Terminated by user while trying to connect to logging device.")

            return nothing
        else
            # if no port can be found or any other error occurs we will try running again
            @warn "Could not find the device on any port. Make sure it is connected.\nTrying again..."
            sleep(5)
            run(; kwargs...)
        end
    end

    # create a file / get the file to write to
    filepath, nlines = logfile_handling(path, lines_per_file)

    try
        while true
            t = string(Dates.now())
            H, T = get_data(s[])

            # round these values
            H = round(H, digits=3)
            T = round(T, digits=2)

            open(filepath, "a") do io
                newline = "$t\t$T\t$H\n"
                write(io, newline)
            end

            nlines += 1
            if nlines >= lines_per_file
                filepath, nlines = logfile_handling(path, lines_per_file)
            end

            sleep(interval)

        end
    catch err
        if err isa InterruptException
            println("Logging terminated by user")
        else
            @warn "A $err occured. Trying to start over..."
            close(s[])
            # do not specify the port! If the device is disconnected
            run(; kwargs...)
        end
    finally
        # close the serial port
        close(s[])
    end

    nothing
end

function logfile_handling(path, lines_per_file)
    if !isdir(path)
        mkpath(path)
        println("""
        Did not find the specified path at $path, so I
        created that for you!\n
        """)
    end

    filepath = get_logfile_for_writing(path, lines_per_file)
end

function get_logfile_for_writing(path, lines_per_file)::Tuple{String,Int}
    # lines_per_file has to be greater than 0. Otherwise we'll be
    # stuck generating empty files
    lines_per_file > 0 || error("Number of lines per file has to be > 0.")

    files = readdir(path)
    filter!(x -> occursin(".log", x), files)
    sort!(files)

    if length(files) > 0
        println("Did find $(length(files)) logfiles")

        file = last(files)
        filepath = joinpath(path, file)
        nlines = get_num_lines(filepath)

        if nlines < lines_per_file
            println("Found logfile $file with $nlines lines. Writing to that.")
            return (filepath, nlines)
        else
            # Found a logfile but it already has more than the maximum
            # allowed number of lines
            println("""Found logfile $file but it already has $nlines lines.
            The maximum number of lines is set to $lines_per_file.\n""")
        end
    else
        println("Did not find a log file at $path.")
    end
    filename = generate_logfile_name()
    filepath = joinpath(path, filename)
    touch(filepath)

    filepath, 0
end

function generate_logfile_name(prefix="hty939", ext=".log")
    datestring = splitext( string(now()) )[1]
    datestring = replace(datestring, ':' => '_')
    prefix * "_" * datestring * ext
end

function get_num_lines(filepath)::Int
    # open file / read as byte array / close file
    f = open(filepath, "r")
    str = read(f)
    close(f)

    nrows = 0
    for byte in str
        # check for newline (\n)
        byte == 0x0a && (nrows += 1)
    end
    nrows
end

function read_line(s; timeout=10)
    readbuffer = IOBuffer()
    t_start = time()
    byte    = ""
    while byte ≠ "\n"
        t = time()
        if bytesavailable(s) > 0
            byte = read(s, 1)
            write(readbuffer, byte)
        end
        if t - t_start > timeout
            throw(TimeoutError)
        end
    end

    str = String(take!(readbuffer))
end

function query(s, msg; timeout=2)
    write(s, msg)
    read_line(s; timeout=timeout)
end

function read_raw_data(s)
    # clear what ever may be there
    readavailable(s)
    readbuffer = query(s, "m")
    parse.(UInt16, strip.(split(readbuffer, ',')))
end

"""
    convert_temperature(T_raw) -> Float64

Convert the raw temperature data from the sensor to the actual temperature in
degree Celsius.
"""
function convert_temperature(T_raw)
    (165 / 16383) * T_raw - 40
end

"""
    convert_humidity(H_raw) -> Float64

Convert the raw humidity data from the sensor to the actual percent relative
humidity.
"""
function convert_humidity(H_raw)
    100 / 16383 * H_raw
end

function get_data(s)
    H_raw, T_raw = read_raw_data(s)
    H = convert_humidity(H_raw)
    T = convert_temperature(T_raw)
    H, T
end

function find_port(;debug=false, baudrate=57600)
    # get a list with all available serial ports
    ports = SerialPorts.list_serialports()

    print("Searching for logger...")

    for p in ports
        print(".")

        try # try to connect to each port in the list
            s = SerialPort(p, baudrate)
            debug && println("$p is open")

            # wait for the device to initialize
            sleep(2.0)

            try # try to write and read from the port
                # clear everything that is in the buffer
                readavailable(s)
                readbuffer = query(s, "i")
                debug && println("Got '$readbuffer' from $p")

                # If we successfully identified the logger close the connection
                # and return the port. Otherwise continue with the next port in
                # the list.
                if readbuffer == "htlogger\n"
                    println("Found the logger on port $p")
                    close(s)
                    return p
                else
                    debug && println("Could read from $p but got wrong identifier.")
                end
            catch # write / read failed
                debug && println("Could not write to port $p")
            finally
                close(s)
            end
        catch # connection failed
            debug && println("Could not open port $p")
        end
    end

    # make a new line
    println()

    error("Could not find logger.")
end

end


text:
# insert functions

"""
    insert_boundary!(tier, time; split_at = 0)
    insert_boundary!(tg, num, time; split_at = 0)

Insert a boundary at `time` in an interval `tier`, which can also be specified
by its number in a `TextGrid`. This action splits an existing interval and
increases the size of the tier by 1.

The keyword argument `split_at` indicates the starting position in the label of
the interval being split that belongs to the right interval. The remaining part
of the label belongs to the left interval. A value outside any valid index for
the label (i.e. the default value 0, any negative integer, or any positive
integer greater than its length) means that the whole original label belongs to
the left interval.
"""
function insert_boundary!(tier::Tier, time::Real; split_at::Int = 0)
    tier.class == "interval" || error("Cannot insert boundaries in point tiers.")
    index = findfirst(x -> x.xmin ≤ time ≤ x.xmax, tier.contents)
    isnothing(index) && error("Time out of bounds of tier.")
    if tier.contents[index].xmin == time || tier.contents[index].xmax == time
        error("Boundary already exists in Tier $tier at $(time)s.")
    end

    label_length = length(tier.contents[index].label)
    split_at ∈ 1:label_length || (split_at = label_length + 1)
    insert!(tier.contents, index + 1, Interval(index + 1, time, tier.contents[index].xmax, tier.contents[index].label[split_at:end]))
    tier.contents[index].xmax = time
    tier.contents[index].label = tier.contents[index].label[1:(split_at - 1)]
    tier.size += 1
    reindex_intervals!(tier, from = index + 1)
end

insert_boundary!(tg::TextGrid, num::Int, time::Real; split_at::Int = 0) = insert_boundary!(tg[num], time, split_at = split_at)

"""
    insert_interval!(tier, start, stop, label; split_at = 0)
    insert_interval!(tg, num, start, stop, label; split_at = 0)

Insert an interval from `start` to `stop` in an interval `tier`. Inserted
intervals must not straddle an existing boundary, but if either `start` or
`stop` coincides with an existing edge of an interval, then the other boundary
is added to split the original interval.

The keyword argument `split_at` indicates the starting position in the label of
the interval being split that belongs to the interval immediately to the right
of the inserted interval. Note that if either `start` or `stop` coincides with
an existing edge, `split_at` is automatically overridden and the whole original
label goes to the split interval. See also `insert_boundary!`.
"""
function insert_interval!(tier::Tier, start::Real, stop::Real, label::AbstractString; split_at::Int = 0)
    tier.class == "interval" || error("Cannot insert intervals in point tiers.")

    start_index, stop_index = is_interval_insertable(tier, start, stop)

    if start == tier.contents[start_index].xmin
        # If only left edge already exists, insert new right edge and force original label to the right
        insert_boundary!(tier, stop, split_at = 1)
        relabel!(tier, start_index, label)
    elseif stop == tier.contents[stop_index].xmax
        # If only right edge already exists, insert new left edge and force original label to the left
        insert_boundary!(tier, start)
        relabel!(tier, start_index + 1, label)
    else
        insert_boundary!(tier, stop, split_at = split_at)
        insert_boundary!(tier, start)
        relabel!(tier, start_index + 1, label)
    end
    return tier
end

insert_interval!(tg::TextGrid, num::Int, start::Real, stop::Real, label::AbstractString; split_at::Int = 0) = insert_interval!(tg[num], start, stop, label, split_at = split_at)

"""
    copy_interval!(tier, source_tier, index; split_at = 0)
    copy_interval!(tg, num, source_num, index; split_at = 0)

Copy the `index`-th interval from `source_tier` to `tier`.

See also `insert_interval!`.
"""
copy_interval!(tier::Tier, source_tier::Tier, index::Int; split_at::Int = 0) = insert_interval!(tier, source_tier.contents[index].xmin, source_tier.contents[index].xmax, source_tier.contents[index].label, split_at = split_at)
copy_interval!(tg::TextGrid, num::Int, source_num::Int, index::Int; split_at::Int = 0) = copy_interval!(tg[num], tg[source_num], index, split_at = split_at)

"""
Insert multiple intervals to an interval `tier`, with boundaries and labels
defined by `starts`, `stops` and `labels`. The keyword argument `split_at` can
be either a single integer if constant across inserted intervals or a vector of
integers if custom per interval.

Intervals are inserted in the specified order, and intervals that cannot be
inserted are simply discarded rather than throw an error.

See also `insert_interval!`
"""
function insert_intervals!(tier::Tier, starts::Vector{<:Real}, stops::Vector{<:Real}, labels::Vector{<:AbstractString}; split_at::Union{Int, Vector{Int}} = 0)
    num_added = minimum(length, [starts, stops, labels])
    if split_at isa Vector
        length(split_at) < num_added && error("`split_at must not be shorter than `starts`, `stops` and `labels`.")
        for (start, stop, label, split_at_i) in zip(starts, stops, labels, split_at)
            try
                insert_interval!(tier, start, stop, label, split_at = split_at_i)
            catch
                continue
            end
        end
    else
        for (start, stop, label) in zip(starts, stops, labels)
            try
                insert_interval!(tier, start, stop, label, split_at = split_at)
            catch
                continue
            end
        end
    end
    return tier
end

insert_intervals!(tg::TextGrid, num::Int, starts::Vector{<:Real}, stops::Vector{<:Real}, labels::Vector{<:AbstractString}; split_at::Union{Int, Vector{Int}} = 0) = insert_intervals!(tg[num], starts, stops, labels, split_at = split_at)

function is_interval_insertable(tier::Tier, start::Real, stop::Real)
    start ≥ stop && error("Start time must come before stop time.")
    start_index = find_interval(tier, start)
    stop_index = find_low_interval(tier, stop)
    start_index * stop_index == 0 && error("Time out of bounds of tier.")
    start_index != stop_index && error("Cannot insert interval that straddles a boundary.")
    start == tier.contents[start_index].xmin && stop == tier.contents[stop_index].xmax && error("Interval already exists with the same boundaries.")
    return start_index, stop_index
end

# remove functions

"""
    remove_left_boundary!(tier, index; delim = "")
    remove_left_boundary!(tg, num, index; delim = "")

Remove the left boundary of the `index`-th interval in a `tier`.

The keyword argument `delim` defines the string used to join the labels from
the two intervals being combined.
"""
function remove_left_boundary!(tier::Tier, index::Int; delim::AbstractString = "")
    tier.class == "interval" || error("Cannot remove boundaries from point tiers.")
    index == 1 && error("Cannot remove left edge of tier.")

    tier.contents[index - 1].xmax = tier.contents[index].xmax
    relabel!(tier, index - 1, tier.contents[index - 1].label * delim * tier.contents[index].label)
    deleteat!(tier.contents, index)
    tier.size -= 1
    reindex_intervals!(tier, from = index)
end

remove_left_boundary!(tg::TextGrid, num::Int, index::Int; delim::AbstractString = "") = remove_left_boundary!(tg[num], index, delim = delim)

"""
    remove_right_boundary!(tier, index; delim = "")
    remove_right_boundary!(tg, num, index; delim = "")

Remove the right boundary of the `index`-th interval in an interval `tier`.

The keyword argument `delim` defines the string used to join the labels from
the two intervals being combined.
"""
function remove_right_boundary!(tier::Tier, index::Int; delim::AbstractString = "")
    tier.class == "interval" || error("Cannot remove boundaries from point tiers.")
    index == tier.size && error("Cannot remove right edge of tier.")

    remove_left_boundary!(tier, index + 1, delim = delim)
end

remove_right_boundary!(tg::TextGrid, num::Int, index::Int; delim::AbstractString = "") = remove_right_boundary!(tg[num], index, delim = delim)

"""
    remove_boundary!(tier, time; delim = "", tolerance = 0.0001)
    remove_boundary!(tier, num, time; delim = "", tolerance = 0.0001)

Remove the boundary at `time` (± `tolerance`) from an interval `tier`.

The keyword argument `delim` defines the string used to join the labels from
the two intervals being combined.
"""
function remove_boundary!(tier::Tier, time::AbstractFloat; delim::AbstractString = "", tolerance::Real = 0.0001)
    index = findlast(x -> abs(x.xmin - time) ≤ tolerance, tier.contents)
    if !isnothing(index)
        remove_left_boundary!(tier, index, delim)
    end
end

remove_boundary!(tg::TextGrid, num::Int, time::AbstractFloat; delim::AbstractString = "", tolerance::Real = 0.0001) = remove_boundary!(tg[num], time, delim = delim, tolerance = tolerance)

"""
    remove_boundary!(tier, index, edge; delim = "")
    remove_boundary!(tier, num, index, edge; delim = "")

Remove left and/or right boundaries of the `index`-th interval in `tier`.
Depending on value of `edge`, equivalent to:
- `"left"`: `remove_left_boundary!`
- `"right"`: `remove_right_boundary!`
- `"both"`: a combination of the two

The keyword argument `delim` defines the string used to join the labels from
the two intervals being combined.
"""
function remove_boundary!(tier::Tier, index::Int, edge::AbstractString; delim::AbstractString = "")
    edge ∈ ("left", "right", "both") || error("`edge` must be `left`, `right` or `both`.")
    edge != "left" && remove_right_boundary!(tier, index, delim = delim)
    edge != "right" && remove_left_boundary!(tier, index, delim = delim)
    return tier
end

remove_boundary!(tg::TextGrid, num::Int, index::Int, edge::AbstractString; delim::AbstractString = "") = remove_boundary!(tg[num], index, edge, delim = delim)

function remove_interval!(tier::Tier, index::Int)
    tier.xmin == tier.contents[index].xmin && tier.xmax == tier.contents[index].xmax && error("Cannot remove only interval in tier $tier.")

    relabel!(tier, index, "")
    edge = if tier.contents[index].xmin == tier.xmin
        "right"
    elseif tier.contents[index].xmax == tier.xmax
        "left"
    else
        "both"
    end

    remove_boundary!(tier, index, edge)
end

remove_interval!(tg::TextGrid, num::Int, index::Int) = remove_interval!(tg[num], index)

# move functions

"""
    move_left_boundary!(tier, index; by = 0, to = -1)
    move_left_boundary!(tg, num, index; by = 0, to = -1)

Move the left boundary of the `index`-th interval in an interval `tier`, either
`by` a specified shift or `to` a specified destination. Boundary movement is
not allowed to cross or overlap with other boundaries.

See also `move_right_boundary!`, `move_boundary!` and `move_interval!`.
"""
function move_left_boundary!(tier::Tier, index::Int; by::Real = 0, to::Real = -1)
    index == 1 && error("Cannot move left edge of tier.")
    by == 0 && to == -1 && error("Boundary movement not specified.")
    by != 0 && to != -1 && error("Both `by` and `to` specified.")

    by != 0 && (to = tier.contents[index].xmin + by)

    (to < tier.xmin || to > tier.xmax) && error("Cannot move boundary out of bounds.")
    (to ≤ tier.contents[index - 1].xmin || to ≥ tier.contents[index].xmax) && error("Cannot move boundary beyond interval duration.")

    tier.contents[index].xmin = to
    tier.contents[index - 1].xmax = to
    return tier
end

move_left_boundary!(tg::TextGrid, num::Int, index::Int; by::Real = 0, to::Real = -1) = move_left_boundary!(tg[num], index, by = by, to = to)

"""
    move_right_boundary!(tier, index; by = 0, to = -1)
    move_right_boundary!(tg, num, index; by = 0, to = -1)

Move the right boundary of the `index`-th interval in an interval `tier`, either
`by` a specified shift or `to` a specified destination. Boundary movement is
not allowed to cross or overlap with other boundaries.

See also `move_left_boundary!`, `move_boundary!` and `move_interval!`.
"""
function move_right_boundary!(tier::Tier, index::Int; by::Real = 0, to::Real = -1)
    index == tier.size && error("Cannot move right edge of tier.")
    move_left_boundary!(tier, index + 1, by = by, to = to)
end

move_right_boundary!(tg::TextGrid, num::Int, index::Int; by::Real = 0, to::Real = -1) = move_right_boundary!(tg[num], index, by = by, to = to)

"""
    move_interval!(tier, index; by = 0, to = -1)
    move_interval!(tg, num, index; by = 0, to = -1)

Move both left and right boundaries of the `index`-th interval in an interval
`tier`, either `by` a specified shift or `to` a specified destination. `to`
specifies the destination of the left boundary. Boundary movement is not allowed
to cross or overlap with other boundaries.

See also `move_left_boundary!`, `move_right_boundary!` and `move_boundary!`.
"""
function move_interval!(tier::Tier, index::Int; by::Real = 0, to::Real = -1)
    index == 1 && error("Cannot move left edge of tier.")
    index == tier.size && error("Cannot move right edge of tier.")
    by == 0 && to == -1 && error("Boundary movement not specified.")
    by != 0 && to != -1 && error("Both `by` and `to` specified.")

    if by != 0
        left_to = tier.contents[index].xmin + by
        right_to = tier.contents[index].xmax + by
    else
        left_to = to
        right_to = to + tier.contents[index].xmax - tier.contents[index].xmin
    end

    (left_to < tier.xmin || right_to > tier.xmax) && error("Cannot move boundary out of bounds.")
    (left_to ≤ tier.contents[index - 1].xmin || right_to ≥ tier.contents[index + 1].xmax) && error("Cannot move boundary beyond interval duration.")

    tier.contents[index].xmin = left_to
    tier.contents[index - 1].xmax = left_to
    tier.contents[index].xmax = right_to
    tier.contents[index + 1].xmin = right_to

    return tier
end

move_interval!(tg::TextGrid, num::Int, index::Int; by::Real = 0, to::Real = -1) = move_interval!(tg[num], index, by = by, to = to)

"""
    move_boundary!(tier, index, edge; by = 0, to = -1)
    move_boundary!(tg, num, index, edge; by = 0, to = -1)

Move one or both boundaries of the `index`-th interval in an interval `tier`,
either `by` a specified shift or `to` a specified destination. Equivalent to
the following depending on value of `edge`:
- "left": `move_left_boundary!`
- "right": `move_right_boundary!`
- "both": `move_interval!`
"""
function move_boundary!(tier::Tier, index::Int, edge::AbstractString; by::Real = 0, to::Real = -1)
    if edge == "left"
        return move_left_boundary!(tier, index, by = by, to = to)
    elseif edge == "right"
        return move_right_boundary!(tier, index, by = by, to = to)
    elseif edge == "both"
        return move_interval!(tier, index, by = by, to = to)
    else
        error("`edge` must be `left`, `right` or `both`.")
    end
end

move_boundary!(tg::TextGrid, num::Int, index::Int, edge::AbstractString; by::Real = 0, to::Real = -1) = move_boundary!(tg[num], index, edge, by = by, to = to)

"""
    move_boundary!(tier, time; by = 0, to = -1, tolerance = 0.0001)
    move_boundary!(tg, num, time; by = 0, to = -1, tolerance = 0.0001)

Move an interval boundary located at `time` (± `tolerance`), either `by` a
specified shift or `to` a specified destination. Boundary movement is not
allowed to cross or overlap with other boundaries. Nothing happens if no
boundary is found within the tolerance limit.
"""
function move_boundary!(tier::Tier, time::AbstractFloat; by::Real = 0, to::Real = -1, tolerance::Real = 0.0001)
    index = find_interval(tier, time)
    if time - tier.contents[index].xmin ≤ tolerance
        return move_left_boundary!(tier, index, by = by, to = to)
    elseif tier.contents[index].xmax - time ≤ tolerance
        return move_right_boundary!(tier, index, by = by, to = to)
    end
end

move_boundary!(tg::TextGrid, num::Int, time::AbstractFloat; by::Real = 0, to::Real = -1, tolerance::Real = 0.0001) = move_boundary!(tg[num], time, by = by, to = to, tolerance = tolerance)


text:
#
## Complex methods
#

function call(::Type{Sidef_Types_Number_Complex}, a::Sidef_Types_Number_Number, b::Sidef_Types_Number_Number)
    Sidef_Types_Number_Number(Complex(a.value, b.value))
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia