text:
function thread_blocks_heuristic(len::Integer)
    # TODO better threads default
    threads = clamp(len, 1, 256)
    blocks = max(ceil(Int, len / threads), 1)
    (blocks,), (threads,)
end

ext: jl


text:
# gaussianprocess.jl

"""
    abstract type AbstractGaussianProcess end

Subtypes: GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.
Supertype: Any

Both subtypes represent a Gaussian process expressed as the inner product of a
Gaussian vector of real-valued coefficients and and a vector of functions. The
distribution is represented by a subtype of AbstractMvNormal from the
Distribution package.

See also: GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.
"""
abstract type AbstractGaussianProcess end

"""
    struct GaussianProcess{T<:AbstractMvNormal} <: AbstractGaussianProcess
        basis::Vector{<:Function}
        distribution::T
    end

Implements a Gaussian process defined as an inner product of a Gaussian vector
with distribution 'distribution' of type T<:AbstractMvNormal of the
Distributions package and a vector of functions, the basis of the function
space. For the use of a Gaussian process expanded in Faber-Schauder functions up
to level j, we recommend using FaberSchauderExpansionWithGaussianCoefficients,
as this makes efficient use of the sparsity structure of Faber-Schauder functions.

See also: FaberSchauderExpansionWithGaussianCoefficients, AbstractGaussianProcess

# Example
```julia
using Distributions, PyPlot
α = 0.5
Π = GaussianProcess([fourier(k) for k in 1:10], MvNormalCanon([k^(α+0.5) for k in 1:10]))
X = rand(Π)
clf()
plot(X)
```
"""
struct GaussianProcess{T<:AbstractMvNormal} <: AbstractGaussianProcess
  basis::Vector{<:Function}
  # distribution is the distribution of the coefficients corresponding to the
  # basis. distribution is a AbstractMvNormal subtype of the Distributions
  # package.
  distribution::T
  function GaussianProcess(basis::Vector{<:Function}, distribution::T) where
        T<:AbstractMvNormal
      length(basis)==length(distribution) || error("The basis and distribution are not of equal length.")
      new{T}(basis, distribution)
  end
end

"""
    calculatedependencystructurefaberschauderbasis(higestlevel::Int64)

Higest level should be at least zero. Returns a
2^(higestlevel+1) x 2^(higestlevel+1)-matrix A (say) with A[i,i'] = true when
Faber-Schauder function i=(j,k) (with i=2^j+k) and i'=(j',k') have essentially
non-empty support, and false otherwise. When an entry A[i,i']=false,
\int_0^1ψ_i(X_t)ψ_{i'}(X_t)dX_t=0 and \int_0^1ψ_i(X_t)ψ_{i'}(X_t)dt are zero,
which speeds up calculations in the calculateposterior method.

# Examples
```
calculatedependencystructurefaberschauderbasis(0)
calculatedependencystructurefaberschauderbasis(2)
```
"""
function calculatedependencystructurefaberschauderbasis(higestlevel::Int64)
    numberofbasisfunctions = 2^(higestlevel+1)
    dependentcoefficients = Array{Bool,2}(numberofbasisfunctions,
        numberofbasisfunctions)
    for k in 1:numberofbasisfunctions
        dependentcoefficients[k,k] = true
    end
    for k in 2:numberofbasisfunctions
        dependentcoefficients[1, k] = true
    end
    for k in 3:numberofbasisfunctions
        dependentcoefficients[2, k] = true
    end
    for j in 1:higestlevel
        for k in 1:2^j
            for l in k+1:2^j
                dependentcoefficients[2^j+k, 2^j+l] = false
            end
        end
    end
    for jone in 1:higestlevel
        for kone in 1:2^jone
            for jtwo in jone+1:higestlevel
                jdifference = jtwo - jone
                for ktwo in 1:2^jtwo
                    if (kone-1)*2^jdifference < ktwo ≤ kone*2^jdifference
                        dependentcoefficients[2^jone+kone, 2^jtwo+ktwo] = true
                    else
                        dependentcoefficients[2^jone+kone, 2^jtwo+ktwo] = false
                    end
                end
            end
        end
    end
    return Symmetric(dependentcoefficients)
end

"""
    struct FaberSchauderExpansionWithGaussianCoefficients <: AbstractGaussianProcess
        higestlevel::Int64
        basis::Vector{<:Function}
        distribution::MvNormalCanon{Float64, PDSparseMat{Float64, SparseMatrixCSC{Float64, Int64}}, S} where S<:Union{Vector{Float64}, ZeroVector{Float64}}
        dependentcoefficients::Symmetric{Bool, Array{Bool, 2}}
    end

Implements a Gaussian Process with Faber-Schauder functions and optimally exploits
the sparsity structure of Faber-Schauder functions.

Constructors:
    FaberSchauderExpansionWithGaussianCoefficients(higestlevel, distribution)
    FaberSchauderExpansionWithGaussianCoefficients(inversevariancesperlevel::Vector{Float64})

The user is not allowed to set basis and dependentcoefficients, they are calculated
from the input. The length of distribution should be equal to 2^(higestlevel+1).
The second constructor defines a Gaussian process with Faber-Schauder basis with
independent coefficients with length(inversevariancesperlevel)-1 levels, so
2^(length(inversevariancesperlevel)) number of basis functions. Where the variance
of the coefficients belonging to level k is inversevariancesperlevel[k+1] (we start with level 0).

# Example
```julia
using PyPlot
clf()
α = 0.5
Π = FaberSchauderExpansionWithGaussianCoefficients([2^(2*j*α) for j in 0:5])
f = rand(Π)
x = -1.0:0.01:2.0
plot(x, f.(x))
```
"""
struct FaberSchauderExpansionWithGaussianCoefficients <: AbstractGaussianProcess
    higestlevel::Int64
    basis::Vector{<:Function}
    distribution::MvNormalCanon{Float64, PDSparseMat{Float64, SparseMatrixCSC{Float64, Int64}}, S} where S<:Union{Vector{Float64}, ZeroVector{Float64}}
    dependentcoefficients::Symmetric{Bool, Array{Bool, 2}}
    function FaberSchauderExpansionWithGaussianCoefficients(higestlevel::Int64,
        distribution::MvNormalCanon{Float64, PDSparseMat{Float64,
        SparseMatrixCSC{Float64,Int64}}, T}) where T<:Union{Vector{Float64},ZeroVector{Float64}}
        length(distribution) == 2^(higestlevel+1) || error("The length of the distribution is not equal to 2^(higestlevel+1).")
        dependentcoefficients = calculatedependencystructurefaberschauderbasis(higestlevel)
        basis = vcat(faberschauderone, [faberschauder(j,k) for j in 0:higestlevel for k in 1:2^j])
        new(higestlevel, basis, distribution, dependentcoefficients)
    end
end

# Constructor
function FaberSchauderExpansionWithGaussianCoefficients(inversevariancesperlevel::Vector{Float64})
    lenghtinversevariancesperlevel = length(inversevariancesperlevel)
    lenghtinversevariancesperlevel == 0 && error("inversevariancesperlevel is of zero length.")
    # We start with level zero.
    higestlevel = lenghtinversevariancesperlevel - 1
    # There are two functions of level zero.
    vectorofinversevariances = repmat(inversevariancesperlevel[1:1],2)
    # and 2^k of level k, k=1,2,...
    for k in 1:higestlevel
        vectorofinversevariances = vcat(vectorofinversevariances, repmat(inversevariancesperlevel[k+1:k+1],2^k))
    end
    distribution = MvNormalCanon(PDSparseMat(spdiagm(vectorofinversevariances)))
    return FaberSchauderExpansionWithGaussianCoefficients(higestlevel, distribution)
end

# Extends Base.length to objects of a subtypes of AbstractGaussianProcess.
"""
    length(Π::AbstractGaussianProcess)

Returns the number of basis functions == length of the distribution of the
coefficients. So for a FaberSchauderExpansionWithGaussianCoefficients object it
this is equal to 2^(higestlevel+1).

# Example
```julia
Π = GaussianProcess([sin, cos], MvNormal([1.,1.]))
length(Π)
#
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^j for j in 0:3])
length(Π)
```
"""
length(Π::AbstractGaussianProcess)=length(Π.basis)

"""
    sumoffunctions(vectoroffunctions::Vector{<:Function}, vectorofscalars::Vector{Float64})

Calculates the 'inner product' of the function vector and the scalar vector. In
other words the sum of the functions weigthed by vectorofscalars. Returns a
function.
"""
function sumoffunctions(vectoroffunctions::Vector{<:Function},
    vectorofscalars::Vector{Float64})
  n = length(vectoroffunctions)
  n == length(vectorofscalars) || error("The vector of functions and the vector of scalars should be of equal length")
  return function(x::Float64)
    value = 0.0;
    for i in 1:n
      value += vectorofscalars[i]*vectoroffunctions[i](x)
    end
    return value
  end
end

# Extend Distributions.rand to objects which are subtypes of
# AbstractGaussianProcess. Returns a function.
"""
    rand(Π::AbstractGaussianProcess)

Returns a random function, where the coefficients have distribution Π.distribution
and the basis functions are defined in Π.basis.

# Example
```julia
using Distributions, PyPlot
distribution = MvNormalCanon(collect(1:10))
Π = GaussianProcess([fourier(k) for  k in 1:10], distribution)
f = rand(Π)
x = -1.0:0.01:2.0
clf()
plot(x, f.(x))
```
"""
function rand(Π::AbstractGaussianProcess)
  Z = rand(Π.distribution)
  return sumoffunctions(Π.basis,Z)
end

"""
    rand(Π::FaberSchauderExpansionWithGaussianCoefficients)

Samples a random function from Gaussian process Π.

# Examples

```julia
using PyPlot, Distributions
clf()
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^j for j in 0:3])
f = rand(Π)
x = -1.0:0.01:2.0
plot(x,f.(x))
```
"""
function rand(Π::FaberSchauderExpansionWithGaussianCoefficients)
    h = Π.distribution.h
    if typeof(h) <: ZeroVector
        h = zeros(length(h))
    end
    J = full(Π.distribution.J)
    Z = rand(MvNormalCanon(h,J))
    return sumoffunctions(Π.basis, Z)
end

# Calculates the stochastic Ito integral.
# Is not exported. Used in calculateposterior.
"""
    calculatestochasticintegral(f::Function, X::AbstractSamplePath)

Internal function. Calculates the stochastic (Ito) integral \int_0^T f(X_t)dX_t,
where T=X.timeinterval[end].
"""
function calculatestochasticintegral(f::Function, X::AbstractSamplePath)::Float64
  return sum(map(f, X.samplevalues[1:end-1]).*(X.samplevalues[2:end] -
    X.samplevalues[1:end-1]))
end

# Calculates the Lebesgue integral.
# Is not exported. Used in calculateposterior.
"""
    calculateLebesgueintegral(f::Function, X::SamplePath)

Is an internal function. Calculates the Lebesgue integral \int_0^T f(X_t)dt,
with T=X.timeinterval[end].
"""
function calculateLebesgueintegral(f::Function, X::SamplePath{S})::Float64 where S <: Range{Float64}
  return sum(map(f, X.samplevalues[1:end-1]))*step(X)
end

# Calculates the Lebesgue integral for SamplePath objects, but would also work
# for SamplePath objects.
function calculateLebesgueintegral(f::Function, X::SamplePath)::Float64
  return sum(map(f, X.samplevalues[1:end-1]).*(X.timeinterval[2:end]-
    X.timeinterval[1:end-1]))
end

# plot(X::SamplePath) extends PyPlot.plot.
# TO DO: needs to be extended, so that other PyPlot options may be used.
"""
    plot(X::AbstractSamplePath)

Equivalent to plot(X.timeinterval, X.samplevalues).
"""
plot(X::AbstractSamplePath) = plot(X.timeinterval, X.samplevalues)

function dividefunctions(f::Function, g::Function)
    return function(x) return f(x)/g(x) end
end

function dividefunctions(f::Function, g::Number)
    return function(x) return f(x)/g end
end

dividefunctions(f::Number, g::Number) = f*g

function multiplyfunctions(f::Function, g::Function)
    return function(x) return f(x)*g(x) end
end

function multiplyfunctions(f::Function, g::Number)
    return function(x) return g*f(x) end
end

multiplyfunctions(f::Number, g::Function) = multiplyfunctions(g, f)

multiplyfunctions(f::Number, g::Float64) = f*g

function calculateGirsanovmatrix(Π::AbstractGaussianProcess,
        X::AbstractSamplePath, M::SDEModel)
    d = length(Π)
    girsanovmatrix = Array{Float64}(d, d)
    for k in 1:d
        for l in k:d
            girsanovmatrix[k,l] = calculateLebesgueintegral(dividefunctions(multiplyfunctions(Π.basis[k], Π.basis[l]),multiplyfunctions(M.σ, M.σ)), X)
        end
    end
    return Symmetric(girsanovmatrix)
end

function calculateGirsanovvector(Π::AbstractGaussianProcess, X::AbstractSamplePath, M::SDEModel)
    return [calculatestochasticintegral(dividefunctions(Π.basis[k], multiplyfunctions(M.σ,M.σ)),X) for k in 1:length(Π)]
end
# The following function calculates the posterior, with prior Π and data X. Π is
# a Gaussian process and X is data satisfying an unknown stochastic differential
# equation dX_t=θ(X_t)dt+dt.
"""
    calculateposterior(Π::GaussianProcess, X::AbstractSamplePath, M::SDEModel{Float64})
    calculateposterior(Π::GaussianProcess, X::AbstractSamplePath, M::SDEModel{<:Function})
    calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients, X::AbstractSamplePath, M::SDEModel{Float64})
    calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients, X::AbstractSamplePath, M::SDEModel{<:Function})

Calculates the posterior distribution Π(⋅∣X), and returns it as a
GaussianProcess object. Uses M to determine the right likelihood.

# Example

```julia
using PyPlot
clf()
sde = SDEWithConstantVariance(x->sinpi(2*x), 1.0, 0.0, 10000.0, 0.01)
X = rand(sde)
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^j for j in 0:5])
M = SDEModel(1.0, 0.0)
postΠ = calculateposterior(Π, X, M)
x=0.0:0.01:1.0
for k in 1:10
    f = rand(postΠ)
    plot(x, f.(x))
end
```
"""
function calculateposterior(Π::GaussianProcess,
        X::AbstractSamplePath, M::SDEModel)::GaussianProcess
    girsanovvector = calculateGirsanovvector(Π, X, M)
    girsanovmatrix = calculateGirsanovmatrix(Π, X, M)
    precisionmatrixposterior = girsanovmatrix + invcov(Π.distribution)
    potentialposterior = girsanovvector + mean(Π.distribution)
    posteriordistribution = MvNormalCanon(potentialposterior,
        precisionmatrixposterior)
    return GaussianProcess(Π.basis, posteriordistribution)
end

function calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients,
        X::AbstractSamplePath, M::SDEModel)::FaberSchauderExpansionWithGaussianCoefficients
    girsanovvector = calculateGirsanovvector(Π, X, M)
    girsanovmatrix = calculateGirsanovmatrix(Π, X, M)
    girsanovmatrix = sparse(girsanovmatrix)
    girsanovmatrix = PDSparseMat(girsanovmatrix)
    precisionmatrixposterior = girsanovmatrix + Π.distribution.J
    potentialposterior = girsanovvector + mean(Π.distribution)
    meanposterior = precisionmatrixposterior \ potentialposterior
    posteriordistribution = MvNormalCanon(meanposterior, potentialposterior,
        precisionmatrixposterior)
    return FaberSchauderExpansionWithGaussianCoefficients(Π.higestlevel,
        posteriordistribution)
end


text:
const LP = MathOptFormat.LP

const LP_TEST_FILE = "test.lp"

@test sprint(show, LP.Model()) == "A .LP-file model"

@testset "write_to_file" begin
    @testset "comprehensive write" begin
        model = LP.Model()
        MOIU.loadfromstring!(model, """
        variables: a, x, y, z
        minobjective: x
        c1: x >= -1.0
        c2: x <= 2.0
        c3: y == 3.0
        c4: z in Interval(4.0, 5.0)
        c5: 1.1x + 0.0 <= 5.1
        c6: 1.3x + -1.4 >= -0.1
        c7: 1.5a + 1.6 == 0.2
        c8: 1.7a + 1.8 in Interval(0.3, 0.4)
        c9: x in ZeroOne()
        c10: y in Integer()
        """)
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
            "minimize\n" *
            "obj: x\n" *
            "subject to\n" *
            "c5: 1.1 x <= 5.1\n" *
            "c6: -1.4 + 1.3 x >= -0.1\n" *
            "c7: 1.6 + 1.5 a = 0.2\n" *
            "c8: 0.3 <= 1.8 + 1.7 a <= 0.4\n" *
            "Bounds\n" *
            "x <= 2\n" *
            "x >= -1\n" *
            "y = 3\n" *
            "4 <= z <= 5\n" *
            "General\n" *
            "y\n" *
            "Binary\n" *
            "x\n" *
            "End\n"

        @test !MOI.is_empty(model)
        MOI.empty!(model)
        @test MOI.is_empty(model)
    end

    @testset "Name sanitisation" begin
        @testset "sanitized_name" begin
            max_length = 15
            o = LP.Options(max_length, true, false, Set{Char}(), Set{Char}())

            @test LP.sanitized_name("x", o) == "x"
            @test LP.sanitized_name(repeat("x", max_length), o) == repeat("x", max_length)

            too_long = repeat("x", max_length + 1)
            @test (@test_logs (:warn, "Name $(too_long) too long (length: $(length(too_long)); maximum: $(max_length)). Truncating.") LP.sanitized_name(too_long, o)) == repeat("x", max_length)

            @test (@test_logs (:warn, "Name .x cannot start with a period, a number, e, or E. Prepending an underscore to name.") LP.sanitized_name(".x", o)) == "_.x"
            @test (@test_logs (:warn, "Name 0x cannot start with a period, a number, e, or E. Prepending an underscore to name.") LP.sanitized_name("0x", o)) == "_0x"
            @test (@test_logs (:warn, "Name Ex cannot start with a period, a number, e, or E. Prepending an underscore to name.") LP.sanitized_name("Ex", o)) == "_Ex"
            @test (@test_logs (:warn, "Name ex cannot start with a period, a number, e, or E. Prepending an underscore to name.") LP.sanitized_name("ex", o)) == "_ex"

            @test (@test_logs (:warn, "Name x*ds contains an illegal character: \"*\". Removing the offending character from name.") LP.sanitized_name("x*ds", o)) == "x_ds"
            @test (@test_logs (:warn, "Name x^ contains an illegal character: \"^\". Removing the offending character from name.") LP.sanitized_name("x^", o)) == "x_"
            @test (@test_logs (:warn, "Name x*ds[1] contains an illegal character: \"*\". Removing the offending character from name.") LP.sanitized_name("x*ds[1]", o)) == "x_ds_1_"
            @test (@test_logs(
                    (:warn, "Name Ex*ds[1] cannot start with a period, a number, e, or E. Prepending an underscore to name."),
                    (:warn, "Name _Ex*ds[1] contains an illegal character: \"*\". Removing the offending character from name."),
                    LP.sanitized_name("Ex*ds[1]", o))
                ) == "_Ex_ds_1_"
        end

        @testset "Whole chain" begin
            model = LP.Model(warn=true)
            MOIU.loadfromstring!(model, """
            variables: a
            minobjective: a
            c1: a >= -1.0
            """)
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[1], "a[]")

            @test_logs (:warn, "Name a[] contains an illegal character: \"[\". Removing the offending character from name.") MOI.write_to_file(model, LP_TEST_FILE)
        end

        @testset "Warn multiple times" begin
            model = LP.Model(warn=true)
            MOIU.loadfromstring!(model, """
            variables: a, b, c
            minobjective: a + b + c
            c1: a + b + c >= -1.0
            """)
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[1], "a[]")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[2], "b[]")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[3], "c[]")

            @test_logs(
                (:warn, "Name a[] contains an illegal character: \"[\". Removing the offending character from name."),
                (:warn, "Name b[] contains an illegal character: \"[\". Removing the offending character from name."),
                (:warn, "Name c[] contains an illegal character: \"[\". Removing the offending character from name."),
                MOI.write_to_file(model, LP_TEST_FILE))
        end

        @testset "Warn once" begin
            model = LP.Model(warn_once=true)
            MOIU.loadfromstring!(model, """
            variables: a, b, c
            minobjective: a + b + c
            c1: a + b + c >= -1.0
            """)
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[1], "a[]")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[2], "b[]")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[3], "c[]")

            @test_logs (:warn, "Name a[] contains an illegal character: \"[\". Removing the offending character from name.") MOI.write_to_file(model, LP_TEST_FILE)
        end

        @testset "Duplicate names after sanitization" begin
            model = LP.Model(warn=true)
            MOIU.loadfromstring!(model, """
            variables: a, b, c, d
            minobjective: a + b + c + d
            c1: a + b + c + d >= -1.0
            """)
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[1], "a[")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[2], "a]")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[3], "a*")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[4], "a^")

            @test_logs (:warn, "Name a[ contains an illegal character: \"[\". Removing the offending character from name.") (:warn, "Name a] contains an illegal character: \"]\". Removing the offending character from name.") (:warn, "Name a* contains an illegal character: \"*\". Removing the offending character from name.") (:warn, "Name a^ contains an illegal character: \"^\". Removing the offending character from name.") MOI.write_to_file(model, LP_TEST_FILE)
            @test read(LP_TEST_FILE, String) ==
                "minimize\n" *
                "obj: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3\n" *
                "subject to\n" *
                "c1: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3 >= -1\n" *
                "Bounds\n" *
                "End\n"
        end

        @testset "Too long duplicate names after sanitization" begin
            max_length = 15

            model = LP.Model(maximum_length=max_length, warn=true)
            MOIU.loadfromstring!(model, """
            variables: a, b, c, d
            minobjective: a + b + c + d
            c1: a + b + c + d >= -1.0
            """)
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[1], "a[")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[2], "a]")
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[3], "a*" * repeat("x", max_length + 1))
            MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[4], "a^" * repeat("x", max_length + 1))

            @test_logs(
                (:warn, "Name a[ contains an illegal character: \"[\". Removing the offending character from name."),
                (:warn, "Name a] contains an illegal character: \"]\". Removing the offending character from name."),
                (:warn, "Name a*$(repeat("x", max_length + 1)) contains an illegal character: \"*\". Removing the offending character from name."),
                (:warn, "Name a_$(repeat("x", max_length + 1)) too long (length: $(2 + max_length + 1); maximum: $(max_length)). Truncating."),
                (:warn, "Name a^$(repeat("x", max_length + 1)) contains an illegal character: \"^\". Removing the offending character from name."),
                (:warn, "Name a_$(repeat("x", max_length + 1)) too long (length: $(2 + max_length + 1); maximum: $(max_length)). Truncating."),
                MOI.write_to_file(model, LP_TEST_FILE))
            @test read(LP_TEST_FILE, String) ==
                "minimize\n" *
                "obj: 1 a_ + 1 a__1 + 1 a_$(repeat("x", max_length - 2)) + 1 a_$(repeat("x", max_length - 4))_1\n" *
                "subject to\n" *
                "c1: 1 a_ + 1 a__1 + 1 a_$(repeat("x", max_length - 2)) + 1 a_$(repeat("x", max_length - 4))_1 >= -1\n" *
                "Bounds\n" *
                "End\n"
        end

        @testset "Many too long duplicate names after sanitization" begin
            max_length = 6

            model = LP.Model(maximum_length=max_length, warn=true)
            MOIU.loadfromstring!(model, """
            variables: a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r
            minobjective: a + b + c + d + e + f + g + h + i + j + k + l + m + n + o + p + q + r
            c1: a + b + c + d + e + f + g + h + i + j + k + l + m + n + o + p + q + r >= -1.0
            """)
            for i in 1:18
                MOI.set(model, MOI.VariableName(), MOI.get(model, MOI.ListOfVariableIndices())[i], repeat("x", max_length + 1))
            end

            @test_logs(
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(2). Renamed to xxxxxxx_1."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(3). Renamed to xxxxxxx_2."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(4). Renamed to xxxxxxx_3."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(5). Renamed to xxxxxxx_4."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(6). Renamed to xxxxxxx_5."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(7). Renamed to xxxxxxx_6."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(8). Renamed to xxxxxxx_7."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(9). Renamed to xxxxxxx_8."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(10). Renamed to xxxxxxx_9."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(11). Renamed to xxxxxxx_10."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(12). Renamed to xxxxxxx_11."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(13). Renamed to xxxxxxx_12."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(14). Renamed to xxxxxxx_13."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(15). Renamed to xxxxxxx_14."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(16). Renamed to xxxxxxx_15."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(17). Renamed to xxxxxxx_16."),
                (:warn, "Duplicate name xxxxxxx detected for variable MathOptInterface.VariableIndex(18). Renamed to xxxxxxx_17."),
                (:warn, "Name $(repeat("x", max_length + 1)) too long (length: $(1 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_1 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_2 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_3 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_4 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_5 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_6 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_7 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_8 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_9 too long (length: $(3 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_10 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_11 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_12 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_13 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_14 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_15 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_16 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                (:warn, "Name $(repeat("x", max_length + 1))_17 too long (length: $(4 + max_length); maximum: $(max_length)). Truncating."),
                MOI.write_to_file(model, LP_TEST_FILE))

            expr = "1 " * repeat("x", max_length) * " + " *
                join(["1 " * repeat("x", max_length - 2) * "_" * string(i) for i in 1:9], " + ") * " + " *
                join(["1 " * repeat("x", max_length - 3) * "_" * string(i) for i in 10:17], " + ")
            @test read(LP_TEST_FILE, String) ==
                "minimize\n" *
                "obj: $(expr)\n" *
                "subject to\n" *
                "c1: $(expr) >= -1\n" *
                "Bounds\n" *
                "End\n"
        end
    end

    @testset "other features" begin
        model = LP.Model()
        MOIU.loadfromstring!(model, """
        variables: x
        maxobjective: 2.0 * x + -1.0
        """)
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
            "maximize\n" *
            "obj: -1 + 2 x\n" *
            "subject to\n" *
            "Bounds\n" *
            "End\n"
    end
end

@testset "read_from_file" begin
    model = LP.Model()
    exception = ErrorException("Read from file is not implemented for LP files.")
    @test_throws exception MOI.read_from_file(model, LP_TEST_FILE)
    @test_throws exception MathOptFormat.read_from_file(LP_TEST_FILE)
end

# Clean up
sleep(1.0)  # Allow time for unlink to happen.
rm(LP_TEST_FILE, force = true)


text: getrealfloor(n::Int)::Int = n - (n > 0) - (n > 13)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
"""
    typetree( s::Symbol )

Returns a subgraph type tree plot of a given type

This is literally a wrapper for basic GraphRecipes functionality.

"""
function typetree(d::Detective, s::Symbol)
    type_in = getfield(d.moduleinst, s)
    return plot(
        type_in;
        method=:tree,
        color=:black,
        nodesize=0.07,
        nodecolor=:lightgray,
        fontsize=10,
        nodeshape=:rect,
    )
end

"""
    magnify( d::Detective, s::Symbol)

Returns a subgraph plot of the first order interactions of a
 given object in a `Detective` instance.

"""
function magnify(d::Detective, s::Symbol, style::Symbol)
    if style == :MIMO
        return MIMO(d, s)
    elseif style == Symbol("Type Tree")
        return typetree(d, s)
    else #This shouldn't happen... ever...
        @error("Plot style($style) requested, not recognized.")
    end
end

"""
    MIMO( d::Detective, s::Symbol)

Returns a subgraph plot of the first order interactions of a
 given object in a `Detective` instance.

"""
function MIMO(d::Detective, s::Symbol)
    if haskey(d.lookup, s)
        ndx = d.lookup[s]
        innb, outnb = inneighbors(d.graph, ndx), outneighbors(d.graph, ndx)
        if (length(innb) + length(outnb)) > 0
            subgraph = SimpleDiGraph(length(unique(vcat(innb, outnb))) + 1)

            for (nb, v) in enumerate(innb)
                add_edge!(subgraph, nb, length(innb) + 1)
            end
            for (nb, v) in enumerate(outnb)
                add_edge!(subgraph, length(innb) + 1, length(innb) + 1 + nb)
            end
            names = vcat([d.tag[i] for i in innb], s, [d.tag[o] for o in outnb])
            x = vcat([0.0 for i in 1:length(innb)], 0.5, [1.0 for o in 1:length(outnb)])
            y = vcat(
                [(i - 1) / (length(innb) - 1) for i in 1:length(innb)],
                0.5,
                [(o - 1) / (length(outnb) - 1) for o in 1:length(outnb)],
            )

            return graphplot(
                subgraph;
                x=x,
                y=y,
                curvature_scalar=0.05,
                nodesize=0.08,
                names=names,
                nodecolor=:lightgray,
                color=:black,
                nodeshape=:rect,
                fontsize=10,
            )
        else
            return "Selected item has no found connections..."
        end
    else
        return "Key $s not found..."
    end
end

"""
    sherlockplot(d::Detective)

Returns a LightGraphs plot recipe of a `Detective` instances internal graph.

"""
function sherlockplot(d::Detective)
    try
        l = @layout [a{0.05h}; b]

        p1 = scatter(
            [0, 1, 2, 3, 4],
            [0, 0, 0, 0, 0];
            markershape=:rect,
            xlim=[-1.65, 5.65],
            ylim=[-0.01, 0.01],
            legend=false,
            axis=nothing,
            ticks=nothing,
            markerstrokecolor=:white,
            border=:none,
            markersize=45,
            size=(999, 70),
            markercolor=[:lightblue, :lightgreen, :tomato, :orchid1, :grey],
            title="$(string(d.modulename)) Overview",
        )
        annotate!(p1, [0.0], [0.00], "Function")
        annotate!(p1, [1.0], [0.00], "Type")
        annotate!(p1, [2.0], [0.00], "Abstract")
        annotate!(p1, [3.0], [0.00], "Dynamic")
        annotate!(p1, [4.0], [0.00], "Unknown")

        colormap = Dict(
            is_function => :lightblue,
            is_type => :lightgreen,
            is_abstract_type => :tomato,
            is_untyped => :orchid1,
            not_found => :grey,
        )

        colors = [colormap[inquire(d, v)] for (k, v) in d.tag]
        colororder = sortperm([k for (k, _) in d.tag])
        p2 = graphplot(
            d.graph;
            size=(999, 999),
            markersize=0.111,
            nodeshape=:rect,
            markercolor=colors[colororder],
            names=[d.tag[i] for i in 1:(d.nv)],
            fontsize=10,
            linecolor=:black,
        )
        return plot(p1, p2; layout=l)
    catch
        return "No edges/connections found..."
    end
end

"""
    sherlock_UI()

This is a convenience function to create a Blink window which hosts an Interact
UI. The intention of the UI is to explore a module (typed in by a user), and display
some of the knowledge obtained about it graphically.

"""
function sherlock_UI()
    module_lbl = "Module Name: "
    available_modules = Observable{String}("Sherlock")
    module_txt = Widgets.textbox(available_modules[])
    module_btn = Widgets.button("Inspect")
    types_to_functions = Widgets.toggle(true; label="Types -> Functions")
    functions_to_functions = Widgets.toggle(true; label="Functions -> Functions")
    views = vbox(types_to_functions, functions_to_functions)
    topload = hbox(pad(1em, module_lbl), pad(1em, module_txt), views, pad(1em, module_btn))

    graphdisplay = Observable{Any}(
        plot(
            [0.0],
            [0.0];
            color=:white,
            border=:none,
            legend=false,
            axis=nothing,
            ticks=nothing,
            title="Please Inspect a Module that is Loaded into Scope. \n Inspection may require 30s for the first run.",
        ),
    )
    mainwindow = vbox(topload, graphdisplay)
    sherlock = Observable(Detective(Sherlock))

    function make_graph(d, selected_module::Symbol, types_to_fns, fns_to_fns)
        try
            if (types_to_fns || fns_to_fns)
                d = Detective(getfield(Main, selected_module))
                if types_to_fns
                    typetype_edges(d)
                end
                if fns_to_fns
                    functiontype_edges(d)
                end
                focus_lbl = "Focus On: "
                focus_options = radiobuttons(["MIMO", "Type Tree"])
                focus_btn = Widgets.button("Focus")
                #get types and functions which have connections
                ts = vcat(d.types, d.abstracttypes)
                fs = d.functions
                nbs = [LightGraphs.neighbors(d.graph, d.lookup[f]) for f in fs]
                fs = fs[length.(nbs) .> 0]
                focus_txt = Widgets.dropdown(vcat(ts, fs))
                focus_frame = vbox(
                    hbox(
                        pad(1em, focus_lbl),
                        pad(1em, focus_txt),
                        pad(1em, focus_options),
                        pad(1em, focus_btn),
                    ),
                )
                map!(
                    x -> vbox(
                        Interact.hline(),
                        focus_frame,
                        magnify(d, Symbol(focus_txt[]), Symbol(focus_options[])),
                    ),
                    graphdisplay,
                    focus_btn,
                )
                return vbox(Interact.hline(), focus_frame, sherlockplot(d))
            else
                return "Please Choose a Module or a View..."
            end
        catch
            return "Module not found, does not exist, or serious error occurred!"
        end
    end
    map!(
        x -> make_graph(
            sherlock[],
            Symbol(module_txt[]),
            types_to_functions[],
            functions_to_functions[],
        ),
        graphdisplay,
        module_btn,
    )

    w = Window(; async=true)
    body!(w, mainwindow; async=true)
    return w
end


text:
# test evolving graph
g = EvolvingGraph{Node{Int}, String}()
add_edge!(g, 1, 2, "t1")
add_edge!(g, 2, 3, "t2")
add_edge!(g, 4, 2, "t2")
add_edge!(g, 4, 2, "t1")
add_edge!(g, 2, 1, "t3")

@test issorted(g) == false

g1 = sort_timestamps(g)

@test issorted(g1) == true
@test num_timestamps(g1) == 3

slice_timestamps!(g1, "t1", "t2")

@test num_timestamps(g1) == 2

p = sortperm(g.timestamps)

for e in edges(g1)
    @test (e in edges(g)) == true
end


text:
# Driver script for all unit tests

# using Distributed
# @everywhere using Pkg

# if abspath(PROGRAM_FILE) == @__FILE__
#     # when run as a script, activate the test environment of this package so that test dependencies are available
#     @everywhere Pkg.activate(dirname(@__FILE__))
#     @everywhere Pkg.instantiate()
# end

# using Jive

# runtests(@__DIR__, skip=["runtests.jl", "helper"])

include("marginalization.jl")


text:
# #############################################################################
# #############################################################################
#Cavity boundary shape

#This file contains functions for specifying the shape of the cavity boundary, by subtypes of the abstract type Boundary. The following are the available constructors, which can take both positional or keyword arguments:

#   CircularBnd()
#   constructs a circular boundary of radius 1.0

#   EllipticalBnd(eccentricity::Float64=0.5)
#   constructs an elliptical boundary of semimajor axis 1.0 and a specified eccentricity (default 0.5).

#   PolygonalBnd(sides::Int64=3,roundness::Float64=0.0)
#   constructs a regular polygonal boundary of a specified number of sides (default 3) and amount of deformation into a circle (0.0 to 1.0, default 0.0).

#   FourierBnd(coeff_c::Array{Float64,1}=[0.0,0.05], coeff_s::Array{Float64,1}=Float64[])
#   constructs a general boundary based on the specified Fourier expansion coefficients of the boundary radial position in polar coordinates, with a constant component of 1.0 (default is r = 1.0 + 0.05*cos(2*theta)).

#Information about the boundary in polar coordinates can be extracted by the following functions, which dispatch according to the boundary subtype to the appropriate calculations for the shape:

#   r = rfunc(bnd::Boundary,theta)
#   returns the radial position of the boundary <bnd> at the angular position <theta>. <theta> can be a Float64 or Array{Float64,N}, and the result will be of the same type and dimensions as <theta>.

#   r,dtheta_r = rderiv(bnd::Boundary,theta)
#   returns a tuple (<r>,<dtheta_r>) of the radial position <r> and its derivative <dtheta_r> with respect to angular position <theta>, for the boundary <bnd> at the angular position <theta>. <theta> can be a Float64 or Array{Float64,N}, and each entry in the tuple result will be of the same type and dimensions as <theta>. There are no functions to compute the derivatives alone since the calculations involved almost always requires the radial position to also be computed.

#   r,dtheta_r,dtheta2_r = rderiv2(bnd::Boundary,theta)
#   this is like rderiv but also includes the second derivative <dtheta2_r> in the resulting tuple.

#   alpha = normang(theta,r,dtheta_r)
#   returns the angle <alpha> of the normal vector at angular position <theta>. <theta> can be a Float64 or Array{Float64,N}, but <r> and <dtheta_r> must follow its type and dimensions, and the result <alpha> will also follow its type and dimensions.

#   r,alpha = rsys(bnd::Boundary,theta)
#   returns the tuple (<r>,<alpha>) of the radial position r and the angle alpha of the normal vector, for the boundary bnd at the angular position theta. Theta can be a Float64 or Array{Float64,N}, and each entry in the tuple result will be of the same type and dimensions as theta. This function is provided because the radial position and normal vector angle are the particular quantities of interest.

#   rsys!(bnd:Boundary,theta::Float64,result::Array{Float64,1})
#   This is the in-place version of rsys storing the results as the first two entries of the array <result>.

# #############################################################################
# #############################################################################
#Initialize

#Dependencies
#require("util.jl")
import Base: ==, hash

#General boundary type
abstract Boundary

# #############################################################################
# #############################################################################
#Hashing and equality

#Extend base hash function with method to hash equivalent Boundary objects identically.
hash(bnd::Boundary) = fieldhashing(bnd)
hash(bnd::Boundary,h::UInt64) = hash(hash(bnd),h)

#Equality is equality of type and field values
==(bnd1::Boundary,bnd2::Boundary) = fieldequality(bnd1,bnd2)

# #############################################################################
# #############################################################################
#Functions for generic boundary

#These functions are vectorized versions of rfunc, rderiv and rderiv2 for a general Boundary, as well as (vectorized/in-place versions of) rsys for use in solving for the ray trajectories. The implementations for specific Boundary subtypes, if present, are slightly more efficient.

#Radial position --------------------------------------------------------------
function rfunc{N}(bnd::Boundary,theta::Array{Float64,N})
    #Initialize array
    r::Array{Float64,N} = similar(theta)
    #Set values to rfunc output
    for i=1:length(theta)
        r[i] = rfunc(bnd,theta[i])
    end
    return r::Array{Float64,N}
end

#Radial position and derivatives ----------------------------------------------
function rderiv{N}(bnd::Boundary,theta::Array{Float64,N})
    #Initialize arrays
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    #Set values to rderiv output
    for i=1:length(theta)
        r[i],dtheta_r[i] = rderiv(bnd,theta[i])
    end
    return (r,dtheta_r)
end

#Radial position and derivatives2 ---------------------------------------------
function rderiv2{N}(bnd::Boundary,theta::Array{Float64,N})
    #Initialize arrays
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    dtheta2_r::Array{Float64,N} = similar(theta)
    #Set values to rderiv2 output
    for i=1:length(theta)
        r[i],dtheta_r[i],dtheta2_r[i] = rderiv2(bnd,theta[i])
    end
    return (r,dtheta_r,dtheta2_r)
end

#Angular position of normal vector --------------------------------------------
function normang(theta::Float64,r::Float64,dtheta_r::Float64)
    return (theta - atan2(dtheta_r,r))::Float64
end

function normang{N}(theta::Array{Float64,N},r::Array{Float64,N},dtheta_r::Array{Float64,N})
    #Initialize arrays
    alpha::Array{Float64,N}  = similar(theta)
    #Set values to normang output
    for i=1:length(theta)
        alpha[i] = normang(theta[i],r[i],dtheta_r[i])
    end
    return alpha
end

#Radial position and normal vector angle --------------------------------------
function rsys(bnd::Boundary,theta::Float64)
    r::Float64, dtheta_r::Float64 = rderiv(bnd,theta)
    alpha::Float64 = normang(theta,r,dtheta_r)
    return (r,alpha)
end

function rsys{N}(bnd::Boundary,theta::Array{Float64,N})
    #Initialize arrays
    r::Array{Float64,N} = similar(theta)
    alpha::Array{Float64,N} = similar(theta)
    #Set values to rsys output
    for i=1:length(theta)
        r[i],alpha[i] = rsys(bnd,theta[i])
    end
    return (r,alpha)
end

#In-place storing version of rsys
function rsys!(bnd::Boundary,theta::Float64,result::Array{Float64,1})
    result[1],dtheta_r::Float64 = rderiv(bnd,theta)
    result[2] = normang(theta,result[1],dtheta_r)
    return nothing
end

# #############################################################################
# #############################################################################
#List of Boundaries

#Code is replicated a lot here to reduce the overhead of function calls, because these calculations are repeated a lot, so speed is fairly important.

# #############################################################################
#Circle

#Circular Boundary of radius 1.0
type CircularBnd <: Boundary
    #Contains no fields
end

#Radial position --------------------------------------------------------------
rfunc(bnd::CircularBnd,theta::Float64) = 1.0::Float64
rfunc{N}(bnd::CircularBnd,theta::Array{Float64,N}) = ones(Float64,size(theta))::Array{Float64,N}

#Radial position and derivatives ----------------------------------------------
rderiv(bnd::CircularBnd,theta::Float64) = (1.0,0.0)
function rderiv{N}(bnd::CircularBnd,theta::Array{Float64,N})
    return (ones(Float64,size(theta)),zeros(Float64,size(theta)))
end

#Radial position and derivatives2 ---------------------------------------------
rderiv2(bnd::CircularBnd,theta::Float64) = (1.0,0.0,0.0)
function rderiv2{N}(bnd::CircularBnd,theta::Array{Float64,N})
    return (ones(Float64,size(theta)),zeros(Float64,size(theta)),zeros(Float64,size(theta)))
end

# #############################################################################
#Ellipse

#Elliptical Boundary with semimajor axis of length 1.0
type EllipticalBnd <: Boundary
    #Eccentricity of the ellipse
    eccentricity::Float64

    #Construct by fieldname
    EllipticalBnd(;eccentricity::Float64=0.5) = new(eccentricity)
    #Construct by position
    EllipticalBnd(eccentricity::Float64=0.5) = new(eccentricity)
end

#Radial position --------------------------------------------------------------
function rfunc(bnd::EllipticalBnd,theta::Float64)
    r::Float64 = 1.0/sqrt(sq(cos(theta))+sq(sin(theta))/(1-sq(bnd.eccentricity)))
    return r::Float64
end

function rfunc{N}(bnd::EllipticalBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    for i=1:length(theta)
        r[i] = 1.0/sqrt(sq(cos(theta[i]))+sq(sin(theta[i]))/(1-sq(bnd.eccentricity)))
    end
    return r::Array{Float64,N}
end

#Radial position and derivatives ----------------------------------------------
function rderiv(bnd::EllipticalBnd,theta::Float64)
    r::Float64 = 1.0/sqrt(sq(cos(theta))+sq(sin(theta))/(1-sq(bnd.eccentricity)))
    dtheta_r::Float64 = -0.5*sin(2*theta)*sq(bnd.eccentricity)*cb(r)/(1-sq(bnd.eccentricity))
    return (r,dtheta_r)
end

function rderiv{N}(bnd::EllipticalBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    for i=1:length(theta)
        r[i] = 1.0/sqrt(sq(cos(theta[i]))+sq(sin(theta[i]))/(1-sq(bnd.eccentricity)))
        dtheta_r[i] = -0.5*sin(2*theta[i])*sq(bnd.eccentricity)*cb(r[i])/(1-sq(bnd.eccentricity))
    end
    return (r,dtheta_r)
end

#Radial position and derivatives2 ---------------------------------------------
function rderiv2(bnd::EllipticalBnd,theta::Float64)
    r::Float64 = 1.0/sqrt(sq(cos(theta))+sq(sin(theta))/(1-sq(bnd.eccentricity)))
    dtheta_r::Float64 = -0.5*sin(2*theta)*sq(bnd.eccentricity)*cb(r)/(1-sq(bnd.eccentricity))
    dtheta2_r::Float64 = -(cos(2*theta)*r + 1.5*sin(2*theta)*dtheta_r)*sq(r)*sq(bnd.eccentricity)/(1-sq(bnd.eccentricity))
    return (r,dtheta_r,dtheta2_r)
end

function rderiv2{N}(bnd::EllipticalBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    dtheta2_r::Array{Float64,N} = similar(theta)
    for i=1:length(theta)
        r[i] = 1.0/sqrt(sq(cos(theta[i]))+sq(sin(theta[i]))/(1-sq(bnd.eccentricity)))
        dtheta_r[i] = -0.5*sin(2*theta[i])*sq(bnd.eccentricity)*cb(r[i])/(1-sq(bnd.eccentricity))
        dtheta2_r[i] = -(cos(2*theta[i])*r[i] + 1.5*sin(2*theta[i])*dtheta_r[i])*sq(r[i])*sq(bnd.eccentricity)/(1-sq(bnd.eccentricity))
    end
    return (r,dtheta_r,dtheta2_r)
end

# #############################################################################
#Polygon

#Regular Polygon with center to vertex distance of 1.0
type PolygonalBnd <: Boundary
    #Number of sides of the regular polygon
    sides::Int64
    #Roundness (the degree of deformation into a circle) taking a value of 0.0 for a polygon to 1.0 for a circle.
    roundedness::Float64;

    #Construct by fieldname
    PolygonalBnd(;sides::Int64=3,roundness::Float64=0.0) = new(sides,roundness)
    #Construct by position
    PolygonalBnd(sides::Int64=3,roundness::Float64=0.0) = new(sides,roundness)
end

#Radial position --------------------------------------------------------------
function rfunc(bnd::PolygonalBnd,theta::Float64)
    angdiv::Float64 = 2*pi/bnd.sides
    thetap::Float64 = mod(theta,angdiv)-0.5*angdiv
    r::Float64 = (1-bnd.roundedness)/cos(thetap) + bnd.roundedness/cos(0.5*angdiv)
    return r::Float64
end

function rfunc{N}(bnd::PolygonalBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    angdiv::Float64 = 2*pi/bnd.sides
    for i=1:length(theta)
        thetap::Float64 = mod(theta[i],angdiv)-0.5*angdiv
        r[i] = (1-bnd.roundedness)./cos(thetap) + bnd.roundedness./cos(0.5*angdiv)
    end
    return r::Array{Float64,N}
end

#Radial position and derivatives ----------------------------------------------
function rderiv(bnd::PolygonalBnd,theta::Float64)
    angdiv::Float64 = 2*pi/bnd.sides
    thetap::Float64 = mod(theta,angdiv)-0.5*angdiv
    r::Float64 = (1-bnd.roundedness)/cos(thetap) + bnd.roundedness/cos(0.5*angdiv)
    if mod(theta,angdiv) != 0.0
        dtheta_r::Float64 = (1-bnd.roundedness)*sec(thetap)*tan(thetap)
    else
        dtheta_r = 0.0
    end
    return (r,dtheta_r)
end

function rderiv{N}(bnd::PolygonalBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    angdiv::Float64 = 2*pi/bnd.sides
    for i=1:length(theta)
        thetap::Float64 = mod(theta[i],angdiv)-0.5*angdiv
        r[i] = (1-bnd.roundedness)./cos(thetap) + bnd.roundedness./cos(0.5*angdiv)
        if mod(theta[i],angdiv) != 0.0
            dtheta_r[i] = (1-bnd.roundedness)*sec(thetap)*tan(thetap)
        else
            dtheta_r[i] = 0.0
        end
    end
    return (r,dtheta_r)
end

#Radial position and derivatives2 ---------------------------------------------
function rderiv2(bnd::PolygonalBnd,theta::Float64)
    angdiv::Float64 = 2*pi/bnd.sides
    thetap::Float64 = mod(theta,angdiv)-0.5*angdiv
    r::Float64 = (1-bnd.roundedness)/cos(thetap) + bnd.roundedness/cos(0.5*angdiv)
    if mod(theta,angdiv) != 0.0
        dtheta_r::Float64 = (1-bnd.roundedness)*sec(thetap)*tan(thetap)
        dtheta2_r::Float64 = (1-bnd.roundedness)*sec(thetap)*(1+2*sq(tan(thetap)))
    else
        dtheta_r = 0.0
        dtheta2_r = 0.0
    end
    return (r,dtheta_r,dtheta2_r)
end

function rderiv2{N}(bnd::PolygonalBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    dtheta2_r::Array{Float64,N} = similar(theta)
    angdiv::Float64 = 2*pi/bnd.sides
    for i=1:length(theta)
        thetap::Float64 = mod(theta[i],angdiv)-0.5*angdiv
        r[i] = (1-bnd.roundedness)./cos(thetap) + bnd.roundedness./cos(0.5*angdiv)
        if mod(theta[i],angdiv) != 0.0
            dtheta_r[i] = (1-bnd.roundedness)*sec(thetap)*tan(thetap)
            dtheta2_r[i] = (1-bnd.roundedness)*sec(thetap)*(1+2*sq(tan(thetap)))
        else
            dtheta_r[i] = 0.0
            dtheta2_r[i] = 0.0
        end
    end
    return (r,dtheta_r,dtheta2_r)
end

# #############################################################################
#Fourier

#General boundary generated by Fourier expansion, normalized to have a constant component of 1.0
type FourierBnd <: Boundary
    #coeff_c and coeff_s are arrays containing respectively the Fourier coefficients for the cosine and sine series expansion of the radius function. Specifically, the nth value of coeff_c is the coefficient of the cos(n*theta) term, etc. They need only be long enough to include the last nonzero coefficient of each series.
    coeff_c::Array{Float64,1}
    coeff_s::Array{Float64,1}

    #Construct by fieldname
    FourierBnd(;coeff_c::Array{Float64,1}=[0.0,0.05],coeff_s::Array{Float64,1}=Float64[]) = new(coeffmatrix(coeff_c),coeffmatrix(coeff_s))
    #Construct by position
    FourierBnd(coeff_c::Array{Float64,1}=[0.0,0.05],coeff_s::Array{Float64,1}=Float64[]) = new(coeffmatrix(coeff_c),coeffmatrix(coeff_s))
end

#Radial position --------------------------------------------------------------
function rfunc(bnd::FourierBnd,theta::Float64)
    r::Float64 = 1.0
    for i = 1:length(bnd.coeff_c)
        r += bnd.coeff_c[i]*cos(i*theta)
    end
    for i = 1:length(bnd.coeff_s)
        r += bnd.coeff_s[i]*sin(i*theta)
    end
    return r::Float64
end

function rfunc{N}(bnd::FourierBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    for i=1:length(theta)
        r[i] = 1.0
        for j = 1:length(bnd.coeff_c)
            r[i] += bnd.coeff_c[j]*cos(j*theta[i])
        end
        for j = 1:length(bnd.coeff_s)
            r[i] += bnd.coeff_s[j]*sin(j*theta[i])
        end
    end
    return r::Array{Float64,N}
end

#Radial position and derivatives ----------------------------------------------
function rderiv(bnd::FourierBnd,theta::Float64)
    r::Float64 = 1.0
    dtheta_r::Float64 = 0.0
    for i = 1:length(bnd.coeff_c)
        r += bnd.coeff_c[i]*cos(i*theta)
        dtheta_r -= bnd.coeff_c[i]*i*sin(i*theta)
    end
    for i = 1:length(bnd.coeff_s)
        r += bnd.coeff_s[i]*sin(i*theta)
        dtheta_r += bnd.coeff_s[i]*i*cos(i*theta)
    end
    return (r,dtheta_r)
end

function rderiv{N}(bnd::FourierBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    for i=1:length(theta)
        r[i] = 1.0
        dtheta_r[i] = 0.0
        for j = 1:length(bnd.coeff_c)
            r[i] += bnd.coeff_c[j]*cos(j*theta[i])
            dtheta_r[i] -= bnd.coeff_c[j]*j*sin(j*theta[i])
        end
        for j = 1:length(bnd.coeff_s)
            r[i] += bnd.coeff_s[j]*sin(j*theta[i])
            dtheta_r[i] += bnd.coeff_s[j]*j*cos(j*theta[i])
        end
    end
    return (r,dtheta_r)
end

#Radial position and derivatives2 ---------------------------------------------
function rderiv2(bnd::FourierBnd,theta::Float64)
    r::Float64 = 1.0
    dtheta_r::Float64 = 0.0
    dtheta2_r::Float64 = 0.0
    for i = 1:length(bnd.coeff_c)
        r += bnd.coeff_c[i]*cos(i*theta)
        dtheta_r -= bnd.coeff_c[i]*i*sin(i*theta)
        dtheta2_r -= bnd.coeff_c[i]*sq(i)*cos(i*theta)
    end
    for i = 1:length(bnd.coeff_s)
        r += bnd.coeff_s[i]*sin(i*theta)
        dtheta_r += bnd.coeff_s[i]*i*cos(i*theta)
        dtheta2_r -= bnd.coeff_s[i]*sq(i)*sin(i*theta)
    end
    return (r,dtheta_r,dtheta2_r)
end

function rderiv2{N}(bnd::FourierBnd,theta::Array{Float64,N})
    r::Array{Float64,N} = similar(theta)
    dtheta_r::Array{Float64,N} = similar(theta)
    dtheta2_r::Array{Float64,N} = similar(theta)
    for i=1:length(theta)
        r[i] = 1.0
        dtheta_r[i] = 0.0
        dtheta2_r[i] = 0.0
        for j = 1:length(bnd.coeff_c)
            r[i] += bnd.coeff_c[j]*cos(j*theta[i])
            dtheta_r[i] -= bnd.coeff_c[j]*j*sin(j*theta[i])
            dtheta2_r[i] -= bnd.coeff_c[j]*sq(j)*cos(j*theta[i])
        end
        for j = 1:length(bnd.coeff_s)
            r[i] += bnd.coeff_s[j]*sin(j*theta[i])
            dtheta_r[i] += bnd.coeff_s[j]*j*cos(j*theta[i])
            dtheta2_r[i] -= bnd.coeff_s[j]*sq(j)*sin(j*theta[i])
        end
    end
    return (r,dtheta_r,dtheta2_r)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -32.35079088597315 <= q <= 33.04707617725232)

                     add_NL_constraint(m, :(gelu(0.7516814527569986 + 0.07211062597860263*gelu(0.800956367833014 + 0.1575422515944389*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.28099609060437203*gelu(-0.8241959750660866 + 0.48752067922981546*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.19401527858092393*gelu(-0.014178711859067938 + -0.6496362145686132*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.43991506890783594*gelu(-0.38559964579408224 + -0.6882459841770032*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + gelu(0.5442086655639207 + -0.859033717055464*gelu(0.800956367833014 + 0.1575422515944389*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.6994107397450926*gelu(-0.8241959750660866 + 0.48752067922981546*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.21447765318835765*gelu(-0.014178711859067938 + -0.6496362145686132*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9247182813449877*gelu(-0.38559964579408224 + -0.6882459841770032*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + gelu(-0.33457802196328057 + 0.8511359109817849*gelu(0.800956367833014 + 0.1575422515944389*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9056512987027685*gelu(-0.8241959750660866 + 0.48752067922981546*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.901983268277478*gelu(-0.014178711859067938 + -0.6496362145686132*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.02854306019224717*gelu(-0.38559964579408224 + -0.6882459841770032*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + gelu(0.85223750045835 + -0.7468516096711739*gelu(0.800956367833014 + 0.1575422515944389*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.1269236474608202*gelu(-0.8241959750660866 + 0.48752067922981546*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.5407492362268407*gelu(-0.014178711859067938 + -0.6496362145686132*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.8670343803826994*gelu(-0.38559964579408224 + -0.6882459841770032*gelu(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*gelu(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*gelu(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*gelu(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using CSV
using DataFrames
using DataFramesMeta
using DataStructures
using Dates
using Statistics
using TimeZones

function read_tw_data(filename)

    tw_data = CSV.read(filename,DataFrame, ntasks = 1)

    @transform!(tw_data, :date = ZonedDateTime.(:created_at,TW_dt))

    return tw_data

end

mentreg = r"((?<=@)\w+)"
hashreg = r"((?<=#)\w+)"
TW_dt = dateformat"yyyy-mm-ddTHH:MM:SS.sssz"

function get_mentions(tweet;regsel = mentreg)
    [lowercase(capture.match) for capture in eachmatch(regsel,lowercase(tweet))]
end

function get_all_username(dataset)
    everyone = Set(dataset.user_username) ∪ Set(vcat(dataset.mentions...))
    return everyone
end

most_common(c::Accumulator) = most_common(c, length(c))
most_common(c::Accumulator, k) = sort(collect(c), by=kv->kv[2], rev=true)[1:k]

get_mention_count(tweet_data) = counter(vcat(get_mentions.(tweet_data...)))


text:
# PETSc error codes and error handling
export chk, PetscError, PetscErrorMessage
#typealias PetscErrorCode Cint # from petscsys.h

# wrap this around all PETSc ccalls, in order to catch exceptions
function chk(errnum)
  if errnum != 0
    throw(PetscError(errnum))
  end
  errnum
end

struct PetscError <: Exception
  errnum::PetscErrorCode
  msg::AbstractString

  function PetscError(n::Integer, msg::AbstractString="")
    if isempty(msg)
      new(PetscErrorCode(n), PetscErrorMessage(n))
    else
      new(PetscErrorMessage(n), msg)
    end
  end
end

# return error message string corresponding to errnum
function PetscErrorMessage(errnum)
  msg_ptr = Ref{Ptr{UInt8}}(C_NULL)
  # use the first petsc library for all cases
  println(stderr, "errnum = ", errnum)
  C.PetscErrorMessage(C.petsc_type[1], PetscErrorCode(errnum), msg_ptr, Ref{Ptr{UInt8}}(C_NULL))
  # error num strings are stored in a constant table so
  # the pointer does not have to be free'd here
  if msg_ptr[] == C_NULL
    return "Petsc did not supply an error message"
  else
    return unsafe_string(msg_ptr[])
  end
end


text:
const DISTMSG = "Right-hand side of a ~ must be subtype of Distribution or a vector of " *
    "Distributions."

const INTERNALNAMES = (:_model, :_sampler, :_context, :_varinfo, :_rng)

"""
    isassumption(expr)

Return an expression that can be evaluated to check if `expr` is an assumption in the
model.

Let `expr` be `:(x[1])`. It is an assumption in the following cases:
    1. `x` is not among the input data to the model,
    2. `x` is among the input data to the model but with a value `missing`, or
    3. `x` is among the input data to the model with a value other than missing,
       but `x[1] === missing`.

When `expr` is not an expression or symbol (i.e., a literal), this expands to `false`.
"""
function isassumption(expr::Union{Symbol, Expr})
    vn = gensym(:vn)

    return quote
        let $vn = $(varname(expr))
            # This branch should compile nicely in all cases except for partial missing data
            # For example, when `expr` is `:(x[i])` and `x isa Vector{Union{Missing, Float64}}`
            if !$(DynamicPPL.inargnames)($vn, _model) || $(DynamicPPL.inmissings)($vn, _model)
                true
            else
                # Evaluate the LHS
                $expr === missing
            end
        end
    end
end

# failsafe: a literal is never an assumption
isassumption(expr) = :(false)

#################
# Main Compiler #
#################

"""
    @model(expr[, warn = true])

Macro to specify a probabilistic model.

If `warn` is `true`, a warning is displayed if internal variable names are used in the model
definition.

# Examples

Model definition:

```julia
@model function model(x, y = 42)
    ...
end
```

To generate a `Model`, call `model(xvalue)` or `model(xvalue, yvalue)`.
"""
macro model(expr, warn=true)
    # include `LineNumberNode` with information about the call site in the
    # generated function for easier debugging and interpretation of error messages
    esc(model(expr, __source__, warn))
end

function model(expr, linenumbernode, warn)
    modelinfo = build_model_info(expr)

    # Generate main body
    modelinfo[:body] = generate_mainbody(
        modelinfo[:modeldef][:body], modelinfo[:allargs_syms], warn
    )

    return build_output(modelinfo, linenumbernode)
end

"""
    build_model_info(input_expr)

Builds the `model_info` dictionary from the model's expression.
"""
function build_model_info(input_expr)
    # Break up the model definition and extract its name, arguments, and function body
    modeldef = MacroTools.splitdef(input_expr)

    # Print a warning if function body of the model is empty
    warn_empty(modeldef[:body])

    ## Construct model_info dictionary

    # Shortcut if the model does not have any arguments
    if !haskey(modeldef, :args) && !haskey(modeldef, :kwargs)
        modelinfo = Dict(
            :allargs_exprs => [],
            :allargs_syms => [],
            :allargs_namedtuple => NamedTuple(),
            :defaults_namedtuple => NamedTuple(),
            :modeldef => modeldef,
        )
        return modelinfo
    end

    # Extract the positional and keyword arguments from the model definition.
    allargs = vcat(modeldef[:args], modeldef[:kwargs])

    # Split the argument expressions and the default values.
    allargs_exprs_defaults = map(allargs) do arg
        MacroTools.@match arg begin
            (x_ = val_) => (x, val)
            x_ => (x, NO_DEFAULT)
        end
    end

    # Extract the expressions of the arguments, without default values.
    allargs_exprs = first.(allargs_exprs_defaults)

    # Extract the names of the arguments.
    allargs_syms = map(allargs_exprs) do arg
        MacroTools.@match arg begin
            (::Type{T_}) | (name_::Type{T_}) => T
            name_::T_ => name
            x_ => x
        end
    end

    # Build named tuple expression of the argument symbols and variables of the same name.
    allargs_namedtuple = to_namedtuple_expr(allargs_syms)

    # Extract default values of the positional and keyword arguments.
    default_syms = []
    default_vals = []
    for (sym, (expr, val)) in zip(allargs_syms, allargs_exprs_defaults)
        if val !== NO_DEFAULT
            push!(default_syms, sym)
            push!(default_vals, val)
        end
    end

    # Build named tuple expression of the argument symbols with default values.
    defaults_namedtuple = to_namedtuple_expr(default_syms, default_vals)

    modelinfo = Dict(
        :allargs_exprs => allargs_exprs,
        :allargs_syms => allargs_syms,
        :allargs_namedtuple => allargs_namedtuple,
        :defaults_namedtuple => defaults_namedtuple,
        :modeldef => modeldef,
    )

    return modelinfo
end

"""
    generate_mainbody(expr, args, warn)

Generate the body of the main evaluation function from expression `expr` and arguments
`args`.

If `warn` is true, a warning is displayed if internal variables are used in the model
definition.
"""
generate_mainbody(expr, args, warn) = generate_mainbody!(Symbol[], expr, args, warn)

generate_mainbody!(found, x, args, warn) = x
function generate_mainbody!(found, sym::Symbol, args, warn)
    if warn && sym in INTERNALNAMES && sym ∉ found
        @warn "you are using the internal variable `$(sym)`"
        push!(found, sym)
    end
    return sym
end
function generate_mainbody!(found, expr::Expr, args, warn)
    # Do not touch interpolated expressions
    expr.head === :$ && return expr.args[1]

    # Apply the `@.` macro first.
    if Meta.isexpr(expr, :macrocall) && length(expr.args) > 1 &&
        expr.args[1] === Symbol("@__dot__")
        return generate_mainbody!(found, Base.Broadcast.__dot__(expr.args[end]), args, warn)
    end

    # Modify dotted tilde operators.
    args_dottilde = getargs_dottilde(expr)
    if args_dottilde !== nothing
        L, R = args_dottilde
        return generate_dot_tilde(generate_mainbody!(found, L, args, warn),
                                  generate_mainbody!(found, R, args, warn),
                                  args) |> Base.remove_linenums!
    end

    # Modify tilde operators.
    args_tilde = getargs_tilde(expr)
    if args_tilde !== nothing
        L, R = args_tilde
        return generate_tilde(generate_mainbody!(found, L, args, warn),
                              generate_mainbody!(found, R, args, warn),
                              args) |> Base.remove_linenums!
    end

    return Expr(expr.head, map(x -> generate_mainbody!(found, x, args, warn), expr.args)...)
end

"""
    generate_tilde(left, right, args)

Generate an `observe` expression for data variables and `assume` expression for parameter
variables.
"""
function generate_tilde(left, right, args)
    @gensym tmpright
    top = [:($tmpright = $right),
           :($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
             || throw(ArgumentError($DISTMSG)))]

    if left isa Symbol || left isa Expr
        @gensym out vn inds
        push!(top, :($vn = $(varname(left))), :($inds = $(vinds(left))))

        # It can only be an observation if the LHS is an argument of the model
        if vsym(left) in args
            @gensym isassumption
            return quote
                $(top...)
                $isassumption = $(DynamicPPL.isassumption(left))
                if $isassumption
                    $left = $(DynamicPPL.tilde_assume)(
                        _rng, _context, _sampler, $tmpright, $vn, $inds, _varinfo)
                else
                    $(DynamicPPL.tilde_observe)(
                        _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
                end
            end
        end

        return quote
            $(top...)
            $left = $(DynamicPPL.tilde_assume)(_rng, _context, _sampler, $tmpright, $vn,
                                               $inds, _varinfo)
        end
    end

    # If the LHS is a literal, it is always an observation
    return quote
        $(top...)
        $(DynamicPPL.tilde_observe)(_context, _sampler, $tmpright, $left, _varinfo)
    end
end

"""
    generate_dot_tilde(left, right, args)

Generate the expression that replaces `left .~ right` in the model body.
"""
function generate_dot_tilde(left, right, args)
    @gensym tmpright
    top = [:($tmpright = $right),
           :($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
             || throw(ArgumentError($DISTMSG)))]

    if left isa Symbol || left isa Expr
        @gensym out vn inds
        push!(top, :($vn = $(varname(left))), :($inds = $(vinds(left))))

        # It can only be an observation if the LHS is an argument of the model
        if vsym(left) in args
            @gensym isassumption
            return quote
                $(top...)
                $isassumption = $(DynamicPPL.isassumption(left))
                if $isassumption
                    $left .= $(DynamicPPL.dot_tilde_assume)(
                        _rng, _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
                else
                    $(DynamicPPL.dot_tilde_observe)(
                        _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
                end
            end
        end

        return quote
            $(top...)
            $left .= $(DynamicPPL.dot_tilde_assume)(
                _rng, _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
        end
    end

    # If the LHS is a literal, it is always an observation
    return quote
        $(top...)
        $(DynamicPPL.dot_tilde_observe)(_context, _sampler, $tmpright, $left, _varinfo)
    end
end

const FloatOrArrayType = Type{<:Union{AbstractFloat, AbstractArray}}
hasmissing(T::Type{<:AbstractArray{TA}}) where {TA <: AbstractArray} = hasmissing(TA)
hasmissing(T::Type{<:AbstractArray{>:Missing}}) = true
hasmissing(T::Type) = false

"""
    build_output(modelinfo, linenumbernode)

Builds the output expression.
"""
function build_output(modelinfo, linenumbernode)
    ## Build the anonymous evaluator from the user-provided model definition.

    # Remove the name.
    evaluatordef = deepcopy(modelinfo[:modeldef])
    delete!(evaluatordef, :name)

    # Add the internal arguments to the user-specified arguments (positional + keywords).
    evaluatordef[:args] = vcat(
        [
            :(_rng::$(Random.AbstractRNG)),
            :(_model::$(DynamicPPL.Model)),
            :(_varinfo::$(DynamicPPL.AbstractVarInfo)),
            :(_sampler::$(DynamicPPL.AbstractSampler)),
            :(_context::$(DynamicPPL.AbstractContext)),
        ],
        modelinfo[:allargs_exprs],
    )

    # Delete the keyword arguments.
    evaluatordef[:kwargs] = []

    # Replace the user-provided function body with the version created by DynamicPPL.
    evaluatordef[:body] = modelinfo[:body]

    ## Build the model function.

    # Extract the named tuple expression of all arguments and the default values.
    allargs_namedtuple = modelinfo[:allargs_namedtuple]
    defaults_namedtuple = modelinfo[:defaults_namedtuple]

    # Update the function body of the user-specified model.
    # We use a name for the anonymous evaluator that does not conflict with other variables.
    modeldef = modelinfo[:modeldef]
    @gensym evaluator
    # We use `MacroTools.@q begin ... end` instead of regular `quote ... end` to ensure
    # that no new `LineNumberNode`s are added apart from the reference `linenumbernode`
    # to the call site
    modeldef[:body] = MacroTools.@q begin
        $(linenumbernode)
        $evaluator = $(MacroTools.combinedef(evaluatordef))
        return $(DynamicPPL.Model)(
            $(QuoteNode(modeldef[:name])),
            $evaluator,
            $allargs_namedtuple,
            $defaults_namedtuple,
        )
    end

    return :($(Base).@__doc__ $(MacroTools.combinedef(modeldef)))
end

function warn_empty(body)
    if all(l -> isa(l, LineNumberNode), body.args)
        @warn("Model definition seems empty, still continue.")
    end
    return
end

"""
    matchingvalue(sampler, vi, value)

Convert the `value` to the correct type for the `sampler` and the `vi` object.
"""
function matchingvalue(sampler, vi, value)
    T = typeof(value)
    if hasmissing(T)
        _value = convert(get_matching_type(sampler, vi, T), value)
        if _value === value
            return deepcopy(_value)
        else
            return _value
        end
    else
        return value
    end
end
matchingvalue(sampler, vi, value::FloatOrArrayType) = get_matching_type(sampler, vi, value)

"""
    get_matching_type(spl::AbstractSampler, vi, ::Type{T}) where {T}

Get the specialized version of type `T` for sampler `spl`.

For example, if `T === Float64` and `spl::Hamiltonian`, the matching type is
`eltype(vi[spl])`.
"""
get_matching_type(spl::AbstractSampler, vi, ::Type{T}) where {T} = T
function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{<:Union{Missing, AbstractFloat}},
)
    return Union{Missing, floatof(eltype(vi, spl))}
end
function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T,N}}) where {T,N}
    return Array{get_matching_type(spl, vi, T), N}
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T}}) where T
    return Array{get_matching_type(spl, vi, T)}
end

floatof(::Type{T}) where {T <: Real} = typeof(one(T)/one(T))
floatof(::Type) = Real # fallback if type inference failed


text:
# assert file to test polynomial implementation
using Base.Test
using Polynomials

pNULL = Poly(Float32[])
p0 = Poly([0])
p1 = Poly([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
p2 = Poly([1,1,0,0])
p3 = Poly([1,2,1,0,0,0,0])
p4 = Poly([1,3,3,1,0,0])
p5 = Poly([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
pN = Poly([276,3,87,15,24,0])
pR = Poly([3//4, -2//1, 1//1])
p1000 = Poly(randn(1000))

@test length(pNULL) == 0
@test length(p1000) == 1000
sprint(show, p1000)
sprint(show, pNULL)

@test p3 == Poly([1,2,1])
@test pN*10 == Poly([2760, 30, 870, 150, 240])
@test pN/10 == Poly([27.6, 0.3, 8.7, 1.5, 2.4])
@test 10*pNULL + pN == pN
@test 10*p0 + pN == pN
@test p5 + 2*p1 == Poly([3,4,6,4,1])
@test 10*pNULL - pN == -pN
@test p0 - pN == -pN
@test p5 - 2*p1 == Poly([-1,4,6,4,1])
@test p2*p2*p2 == p4
@test p2^4 == p5
@test pNULL^3 == pNULL
@test pNULL*pNULL == pNULL

@test polyval(pN, -.125) == 276.9609375
@test polyval(pNULL, 10) == 0
@test polyval(p0, -10) == 0
@test polyval(poly([1//2, 3//2]), 1//2) == 0//1
@test polyder(polyint(pN)) == pN
@test polyder(pR) == Poly([-2//1,2//1])
@test polyint(pNULL,1) == p1
@test polyint(Poly(Rational[1,2,3])) == Poly(Rational[0, 1, 1, 1])
@test polyder(p3) == Poly([2,2])
@test polyder(p1) == polyder(p0) == polyder(pNULL) == pNULL

@test poly([-1,-1]) == p3
@test roots(p0)==roots(p1)==roots(pNULL)==[]
@test roots(p2) == [-1]
a_roots = copy(pN.a)
@test all(abs(sort(roots(poly(a_roots))) - sort(a_roots)) .< 1e6)
@test length(roots(p5)) == 4
@test roots(pNULL) == []
@test roots(pR) == [1//2, 3//2]

@test pNULL + 2 == p0 + 2 == 2 + p0 == Poly([2])
@test p2 - 2 == -2 + p2 == Poly([-1,1])
@test 2 - p2 == Poly([1,-1])

p0 = Poly([0])
p1 = Poly([1])
p2 = Poly([5, 6, -3, 2 ,4])
p3 = Poly([7, -3, 2, 6])
p4 = p2 * p3
@test divrem(p4, p2) == (p3, zero(p3))
@test p3%p2 == p3
@test all((abs((p2/p3 - Poly([1/9,2/3])).a)) .< eps())
@test divrem(p0,p1) == (p0,p0)
@test divrem(p1,p1) == (p1,p0)
@test divrem(p2,p2) == (p1,p0)
@test divrem(pR, pR) == (one(pR), zero(pR))
@test_throws p1/p0
@test_throws divrem(p0,p0)

#Tests for multivariable support
pX = Poly([1, 2, 3, 4, 5])
pS1 = Poly([1, 2, 3, 4, 5], "s")
pS2 = Poly([1, 2, 3, 4, 5], 's')
pS3 = Poly([1, 2, 3, 4, 5], :s)
@test pX != pS1
@test pS1 == pS2
@test pS1 == pS3
@test_throws pS1 + pX
@test_throws pS1 - pX
@test_throws pS1 * pX
@test_throws pS1 / pX
@test_throws pS1 % pX


text:
module GeoJSON

import GeoInterface
import JSON3

export dict2geo, geo2dict

"""
    read(input::Union{AbstractString, IO, AbstractVector{UInt8}})

Read a GeoJSON string or IO stream into a GeoInterface object.

To read a file, use `GeoJSON.read(read(path))`.

!!! note
    GeoJSON.jl loads features into the GeoInterface.jl format and that this differs from GeoJSON in the following ways:

    - Julia Geometries do not provide a `bbox` and `crs` method. If you wish to provide a `bbox` or `crs` attribute, wrap the geometry into a `Feature` or `FeatureCollection`.
    - Features do not have special fields for `id`, `bbox`, and `crs`. These are to be provided (or found) in the `properties` field, under the keys `featureid`, `bbox`, and `crs` respectively (if they exist).

    When saving GeoJSON, these transformations will be reversed: if `properties` has a key `featureid`, that will be removed from `properties` and a matching member `id` will be added to the Feature; similarly for `crs` and `bbox`.

# Examples
```julia
julia> GeoJSON.read("{\"type\": \"Point\", \"coordinates\": [30, 10]}")
GeoInterface.Point([30.0, 10.0])
```
"""
read(input) = dict2geo(JSON3.read(input))

"""
    write(obj)

Create a GeoJSON string from an object that implements the GeoInterface, either
`AbstractGeometry`, `AbstractFeature` or `AbstractFeatureCollection`.

!!! note
    GeoJSON.jl loads features into the GeoInterface.jl format and that this differs from GeoJSON in the following ways:

    - Julia Geometries do not provide a `bbox` and `crs` method. If you wish to provide a `bbox` or `crs` attribute, wrap the geometry into a `Feature` or `FeatureCollection`.
    - Features do not have special fields for `id`, `bbox`, and `crs`. These are to be provided (or found) in the `properties` field, under the keys `featureid`, `bbox`, and `crs` respectively (if they exist).

    When saving GeoJSON, these transformations will be reversed: if `properties` has a key `featureid`, that will be removed from `properties` and a matching member `id` will be added to the Feature; similarly for `crs` and `bbox`.

# Examples
```julia
julia> GeoJSON.write(Point([30.0, 10.0]))
\"{\"coordinates\":[30.0,10.0],\"type\":\"Point\"}\"
```
"""
function write end

for geom in (:AbstractFeatureCollection, :AbstractGeometryCollection, :AbstractFeature,
        :AbstractMultiPolygon, :AbstractPolygon, :AbstractMultiLineString,
        :AbstractLineString, :AbstractMultiPoint, :AbstractPoint)
    @eval write(obj::GeoInterface.$geom) = JSON3.write(geo2dict(obj))
end

"""
    dict2geo(obj::AbstractDict{<:Union{Symbol, String}, Any})

Transform a JSON dictionary to a GeoInterface object.

See also: [`geo2dict`](@ref)

# Examples
```julia
julia> dict2geo(Dict("type" => "Point", "coordinates" => [30.0, 10.0]))
Point([30.0, 10.0])
```
"""
function dict2geo(obj::AbstractDict{<:Union{Symbol, String}, Any})
    t = Symbol(obj["type"])
    if t == :FeatureCollection
        return parseFeatureCollection(obj)
    elseif t == :Feature
        return parseFeature(obj)
    elseif t == :GeometryCollection
        return parseGeometryCollection(obj)
    elseif t == :MultiPolygon
        return GeoInterface.MultiPolygon(obj["coordinates"])
    elseif t == :Polygon
        return GeoInterface.Polygon(obj["coordinates"])
    elseif t == :MultiLineString
        return GeoInterface.MultiLineString(obj["coordinates"])
    elseif t == :LineString
        return GeoInterface.LineString(obj["coordinates"])
    elseif t == :MultiPoint
        return GeoInterface.MultiPoint(obj["coordinates"])
    elseif t == :Point
        return GeoInterface.Point(obj["coordinates"])
    end
end

dict2geo(obj::Nothing) = nothing

parseGeometryCollection(obj::AbstractDict{<:Union{Symbol, String}, Any}) =
    GeoInterface.GeometryCollection(dict2geo.(obj["geometries"]))

function parseFeature(obj::AbstractDict{<:Union{Symbol, String}, Any})
    properties = Dict{String, Any}(String(k) => v for (k, v) in obj["properties"])
    feature = GeoInterface.Feature(dict2geo(obj["geometry"]), properties)
    if haskey(obj, "id")
        feature.properties["featureid"] = obj["id"]
    end
    if haskey(obj, "bbox")
        feature.properties["bbox"] = GeoInterface.BBox(obj["bbox"])
    end
    if haskey(obj, "crs")
        feature.properties["crs"] = obj["crs"]
    end
    feature
end

function parseFeatureCollection(obj::AbstractDict{<:Union{Symbol, String}, Any})
    features = parseFeature.(obj["features"])
    featurecollection = GeoInterface.FeatureCollection(features)
    if haskey(obj, "bbox")
        featurecollection.bbox = GeoInterface.BBox(obj["bbox"])
    end
    if haskey(obj, "crs")
        featurecollection.crs = Dict{String, Any}(String(k) => v for (k, v) in obj["crs"])
    end
    featurecollection
end

"""
    geo2dict(obj)

Transform a GeoInterface object to a JSON dictionary.

See also: [`dict2geo`](@ref)

# Examples
```julia
julia> geo2dict(Point([30.0, 10.0]))
Dict{String,Any} with 2 entries:
  "coordinates" => [30.0, 10.0]
  "type"        => "Point"
```
"""
function geo2dict end

function geo2dict(obj::GeoInterface.AbstractGeometry)
    Dict("type" => string(GeoInterface.geotype(obj)),
        "coordinates" => GeoInterface.coordinates(obj))
end

function geo2dict(obj::GeoInterface.AbstractGeometryCollection)
    Dict("type" => string(GeoInterface.geotype(obj)),
        "geometries" => map(geo2dict, GeoInterface.geometries(obj)))
end

function geo2dict(obj::GeoInterface.AbstractFeature)
    result = Dict("type" => string(GeoInterface.geotype(obj)),
        "geometry" => geo2dict(GeoInterface.geometry(obj)),
        "properties" => copy(GeoInterface.properties(obj)))
    if haskey(result["properties"], "bbox")
        result["bbox"] = result["properties"]["bbox"]
        delete!(result["properties"], "bbox")
    end
    if haskey(result["properties"], "crs")
        result["crs"] = result["properties"]["crs"]
        delete!(result["properties"], "crs")
    end
    if haskey(result["properties"], "featureid")
        result["id"] = result["properties"]["featureid"]
        delete!(result["properties"], "featureid")
    end
    result
end

function geo2dict(obj::GeoInterface.AbstractFeatureCollection)
    result = Dict("type" => string(GeoInterface.geotype(obj)),
        "features" => map(geo2dict, GeoInterface.features(obj)))
    if GeoInterface.bbox(obj) != nothing
        result["bbox"] = GeoInterface.bbox(obj)
    end
    if GeoInterface.crs(obj) != nothing
        result["crs"] = GeoInterface.crs(obj)
    end
    result
end

geo2dict(obj::Nothing) = nothing

include("deprecations.jl")

end  # module


text:
"""
**Logistic Likelihood**

Bernoulli likelihood with a logistic link for the Bernoulli likelihood
    ``p(y|f) = \\sigma(yf) = \\frac{1}{1+\\exp(-yf)}``, (for more info see : [wiki page](https://en.wikipedia.org/wiki/Logistic_function))

```julia
LogisticLikelihood()
```

---

For the analytic version the likelihood, it is augmented via:
```math
p(y|f,\\omega) = \\exp\\left(\\frac{1}{2}\\left(yf - (yf)^2 \\omega\\right)\\right)
```
where ``\\omega \\sim \\text{PG}(\\omega\\mid 1, 0)``, and PG is the Polya-Gamma distribution
See paper : [Efficient Gaussian Process Classification Using Polya-Gamma Data Augmentation](https://arxiv.org/abs/1802.06383)
"""
struct LogisticLikelihood{T<:Real} <: ClassificationLikelihood{T}
    c::AbstractVector{AbstractVector{T}}
    θ::AbstractVector{AbstractVector{T}}
    function LogisticLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function LogisticLikelihood{T}(c::AbstractVector{<:AbstractVector{<:Real}},θ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(c,θ)
    end
end

function LogisticLikelihood()
    LogisticLikelihood{Float64}()
end

function init_likelihood(likelihood::LogisticLikelihood{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Integer) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        LogisticLikelihood{T}([abs.(rand(T,nSamplesUsed)) for _ in 1:nLatent],[zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        LogisticLikelihood{T}()
    end
end

function pdf(l::LogisticLikelihood,y::Real,f::Real)
    logistic(y*f)
end

function Base.show(io::IO,model::LogisticLikelihood{T}) where T
    print(io,"Bernoulli Likelihood with Logistic Link")
end

function compute_proba(l::LogisticLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    for i in 1:N
        if σ²[i] <= 0.0
            pred[i] = logistic(μ[i])
        else
            pred[i] =  expectation(logistic,Normal(μ[i],sqrt(σ²[i])))
        end
    end
    return pred
end

### Local Updates Section ###

function local_updates!(model::VGP{<:LogisticLikelihood,<:AnalyticVI})
    model.likelihood.c .= broadcast((μ,Σ)->sqrt.(Σ+abs2.(μ)),model.μ,diag.(model.Σ))
    model.likelihood.θ .= broadcast(c->0.5*tanh.(0.5*c)./c,model.likelihood.c)
end

function local_updates!(model::SVGP{<:LogisticLikelihood,<:AnalyticVI})
    model.likelihood.c .= broadcast((μ,Σ,K̃,κ)->sqrt.(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)),model.μ,model.Σ,model.K̃,model.κ)
    model.likelihood.θ .= broadcast(c->0.5*tanh.(0.5*c)./c,model.likelihood.c)
end

function sample_local!(model::VGP{<:LogisticLikelihood,<:GibbsSampling})
    pg = PolyaGammaDist()
    model.likelihood.θ .= broadcast((μ::AbstractVector{<:Real})->draw.([pg],[1.0],μ),model.μ)
    return nothing
end

### Natural Gradient Section ###

function cond_mean(model::VGP{<:LogisticLikelihood,<:AnalyticVI},index::Integer)
    return 0.5*model.y[index]
end

function ∇μ(model::VGP{<:LogisticLikelihood})
    return 0.5*model.y
end

function cond_mean(model::SVGP{<:LogisticLikelihood,<:AnalyticVI},index::Integer)
    return 0.5.*model.y[index][model.inference.MBIndices]
end

function ∇μ(model::SVGP{<:LogisticLikelihood})
    return 0.5.*getindex.(model.y,[model.inference.MBIndices])
end

function ∇Σ(model::AbstractGP{<:LogisticLikelihood})
    return model.likelihood.θ
end

### ELBO Section ###

function ELBO(model::AbstractGP{<:LogisticLikelihood,<:AnalyticVI})
    return expecLogLikelihood(model) - GaussianKL(model) - PolyaGammaKL(model)
end

function expecLogLikelihood(model::VGP{<:LogisticLikelihood,<:AnalyticVI})
    tot = -model.nLatent*(0.5*model.nSample*logtwo)
    tot += sum(broadcast((μ,y,θ,Σ)->0.5.*(sum(μ.*y)-dot(θ,Σ+abs2.(μ))),
                        model.μ,model.y,model.likelihood.θ,diag.(model.Σ)))
    return tot
end

function expecLogLikelihood(model::SVGP{<:LogisticLikelihood,<:AnalyticVI})
    tot = -model.nLatent*(0.5*model.inference.nSamplesUsed*logtwo)
    tot += sum(broadcast((κμ,y,θ,κΣκ,K̃)->0.5.*(sum(κμ.*y[model.inference.MBIndices])-dot(θ,K̃+κΣκ+abs2.(κμ))),
                        model.κ.*model.μ,model.y,model.likelihood.θ,opt_diag.(model.κ.*model.Σ,model.κ),model.K̃))
    return model.inference.ρ*tot
end

function PolyaGammaKL(model::AbstractGP{<:LogisticLikelihood})
    model.inference.ρ*sum(broadcast(PolyaGammaKL,[ones(length(model.likelihood.c[1]))],model.likelihood.c,model.likelihood.θ))
end

### Gradient Section ###

function gradpdf(::LogisticLikelihood,y::Int,f::T) where {T<:Real}
    σ=y*f
    σ*(one(T)-σ)
end

function hessiandiagpdf(::LogisticLikelihood,y::Int,f::T) where {T<:Real}
    σ=y*f
    σ*(one(T)-2σ + abs2(σ))
end