text:
# Extract out common functions/etc for the Extended Types

for (Name, T) in ((InfExtendedReal, Real), (InfExtendedTime, TimeType))
    @eval begin
        # base.jl
        $Name{T}(x::$Name{T}) where {T<:$T} = x

        Utils.posinf(::Type{T}) where {T<:$Name} = T(PosInf)
        Utils.neginf(::Type{T}) where {T<:$Name} = T(NegInf)

        # arithmetic.jl
        Base.typemin(::Type{T}) where {T<:$Name} = neginf(T)
        Base.typemax(::Type{T}) where {T<:$Name} = posinf(T)

        # io.jl
        function Base.show(io::IO, x::T) where {T<:$Name}
            if get(io, :compact, false)
                print(io, x.val)
            else
                print(io, "$T(")
                show(io, x.val)
                print(io, ")")
            end
        end

        # comparison.jl
        Base.isfinite(x::$Name) = x.flag == FINITE && isfinite(x.finitevalue)
        Base.isinf(x::$Name) = isposinf(x) || isneginf(x)

        Base.:≤(x::$Name, y::$Name) = !(y < x)

        # conversion.jl
        Base.promote_rule(::Type{Infinite}, ::Type{T}) where {T<:$T} = T <: $Name ? T : $Name{T}
        Base.promote_rule(::Type{$Name{T}}, ::Type{$Name{S}}) where {T<:$T, S<:$T} = $Name(promote_type(T, S))
        Base.promote_rule(::Type{$Name{T}}, ::Type{S}) where {T<:$T, S<:$T} = $Name(promote_type(T, S))
        Base.promote_rule(::Type{$Name{T}}, ::Type{Infinite}) where {T<:$T} = $Name{T}

        Base.convert(::Type{T}, x::S) where {T<:$Name, S<:$T} = T(x)
        Base.convert(::Type{T}, x::$Name) where {T<:$Name} = T(x)
        Base.convert(::Type{T}, x::Infinite) where {T<:$Name} = T(x)
    end
end


text:
# nanstd.jl
#
#   Standard deviation, ignoring NaNs
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.4.0
# 09.04.2020
# Last Edit: 09.04.2020

nanstd(x::Array{T}) where T <: Number = std(filter(!isnan,x))

nanstd(x::Array{Union{T,Missing}}) where T <: Number = std(filter(!isnan,filter(!ismissing,x)))

ext: jl


text:
"""
# chebyshevfun
普通递归给出第一类切比雪夫函数的表达式
```math
T_n(x)=2xT_{n-1}(x)-T_{n-2}(x)
```
"""
function chebyshevT(n::Int64,x)
    if n<0
        error("chebyshev(n,x), n>=0 is violated")
    elseif n==0
        return 1
    elseif n==1
        return x
    elseif n==2
        return 2*x^2-1
    elseif n==3
        return 4*x^3-3*x
    elseif n==4
        return 8*x^4-8*x^2+1
    elseif n==5
        return 16*x^5-20*x^3+5*x
    elseif n==6
        return 16*x^5-20*x^3+5*x
    else
        return 2*x*chebyshevT(n-1,x)-chebyshevT(n-2,x)
    end
end

"""
# chebyshevfun
普通递归给出第一类切比雪夫函数的表达式
```math
U_n(x)=2xU_{n-1}(x)-U_{n-2}(x)
```
"""
function chebyshevU(n::Int64,x)
    if n<0
        error("chebyshev(n,x), n>=0 is violated")
    elseif n==0
        return 1
    elseif n==1
        return 2*x
    elseif n==2
        return -1 + 4*x^2
    elseif n==3
        return -4*x + 8*x^3
    elseif n==4
        return 1 - 12*x^2 + 16*x^4
    elseif n==5
        return 6*x - 32*x^3 + 32*x^5
    elseif n==6
        return -1 + 24*x^2 - 80*x^4 + 64*x^6
    else
        return 2*x*chebyshevU(n-1,x)-chebyshevU(n-2,x)
    end
end

# """
# # chebyshevexpand
# 切比雪夫展开函数
# 输入一个函数，返回在指定点处切比雪夫展开后的版本
# """
"""
为切比雪夫T而写的宏
"""
macro chebyT(n,x)
    return esc(:(2*$x*chebyshevT($n-1,$x)-chebyshevT($n-2,$x)))
end

"""
# chebyshevTmeta
为切比雪夫展开的多项式
用元编程进行操作
"""
function chebyshevTmeta(n::Int64,x)
    if n==0
        return 1
    elseif n==1
        return x
    elseif n>=2
        return @chebyT(n,x)
    end
end

# """
# 为切比雪夫U而写的宏
# """
# macro chebyU(n,x)
#     return esc(:(2*x*chebyshevT($n-1,$x)-chebyshevT($n-2,$x)))
# end

# """
# # chebyshevTmeta
# 为切比雪夫展开的多项式
# 用元编程进行操作
# """
# function chebyshevUmeta(n::Int64,x)
#     if n==0
#         return 1
#     elseif n==1
#         return 2*x
#     elseif n>=2
#         return @chebyU(n,x)
#     end
# end


text:
using SparseMatrixDicts
using Test
using SparseArrays

# Local functions:
#   issamematrix, make_test_matrix, convert_test_to_dict
#   generate_range
include("functions.jl")

# two sypes of indices
@testset "Tuple and Vector Index" begin
   m = 7
   n = 5
   A = SparseMatrixDict(m,n)
   B = SparseMatrixDict(m,n)
   M = zeros(Float64,m,n)
   for j=1:n
      for i=1:m
         val = (10.0*n)*i+j
         A[(i,j)] = val
         B[i,j] = val
         M[i,j] = val
      end
   end
   @test nnz(A)==length( myfindnz(Matrix(M))[1] )
   @test nnz(A)==nnz(B)
   for key in collect(keys(A.dict))
      @test B[key]==A[key]
   end
end

# allowance of different types
@testset "Copy and Conversion   " begin
   # to/from dense/sparse matrix
   m = 25
   n = 20
   r = 0.5
   ntests = 10
   for Tx in (Matrix,SparseMatrixCSC)
      for Tv in (Float64,Float32,Float16,Int64,Int32,Int16,Int8,Int)
         for Ti in (Int64,Int32,Int16,Int8,Int,Any)
            for round=1:ntests
               # conversion
               M = make_test_matrix(Tx,Tv,m,n,r)
               A = convert_test_to_dict(M,Ti)
               @test eltype(M)==eltype(A.dict.vals)
               @test issamematrix(A,M)
               M = Matrix{Tv}(A)
               @test eltype(M)==eltype(A.dict.vals)
               @test issamematrix(A,M)
               # copy
               B = copy(A)
               @test eltype(B.dict.vals)==eltype(A.dict.vals)
               @test eltype(B.dict.keys[1])==eltype(A.dict.keys[1])
               @test issamematrix(A,B)
            end
         end
      end
   end
end

# getindex
@testset "GetIndex              " begin
   m = 17
   n = 14
   r = 0.5
   ntests = 20
   Tx = Matrix
   for Tv in (Float64,Float32,Float16,Int64,Int32,Int16,Int8,Int)
      for Ti in (Int64,Int32,Int16,Int8,Int,Any)
         # conversion
         M = make_test_matrix(Tx,Tv,m,n,r)
         A = convert_test_to_dict(M,Ti)
         # indices
         for IDXj in (Number,AbstractRange,Array,Bool,Colon)
            for IDXi in (Number,AbstractRange,Array,Bool,Colon)
               for round=1:ntests
                  J,lenJ = generate_range(IDXj, 1, n)
                  I,lenI = generate_range(IDXi, 1, m)
                  X = A[I,J]
                  Y = M[I,J]
                  if lenI==1 && lenJ>1
                     Y = reshape(Y,1,length(Y))
                  elseif lenI>1 && lenJ==1
                     Y = reshape(Y,length(Y),1)
                  end
                  @test issamematrix(X,Y)
               end
            end
         end
      end
   end
end

# setindex simple
@testset "SetIndex Simple       " begin
   m = 18
   n = 21
   r = 0.1
   ntests = 20
   nchanges = 100
   Tx = Matrix
   for Tv in (Float64,Float32,Float16,Int64,Int32,Int16,Int8,Int)
      for Ti in (Int64,Int32,Int16,Int8,Int,Any)
         # conversion
         M = make_test_matrix(Tx,Tv,m,n,r)
         A = convert_test_to_dict(M,Ti)
         for round=1:ntests
            for k=1:nchanges
               i = rand(1:m)
               j = rand(1:m)
               v = rand(Tv)
               A[i,j] = A[i,j] + v
               M[i,j] = M[i,j] + v
            end
            @test issamematrix(A,M)
         end
      end
   end
end

# setindex
@testset "SetIndex              " begin
   m = 17
   n = 21
   r = 0.5
   ntests = 20
   Tx = Matrix
   for Tv in (Float64,Float32,Float16,Int64,Int32,Int16,Int8,Int)
      for Ti in (Int64,Int32,Int16,Int8,Int,Any)
         # conversion
         M = make_test_matrix(Tx,Tv,m,n,r)
         A = convert_test_to_dict(M,Ti)
         N = copy(M)
         B = copy(A)
         # indices
         for IDXj in (Number,AbstractRange,Array,Bool,Colon)
            for IDXi in (Number,AbstractRange,Array,Bool,Colon)
               for round=1:ntests
                  J,lenJ = generate_range(IDXj, 1, n)
                  I,lenI = generate_range(IDXi, 1, m)
                  if lenI==1 && lenJ==1
                     V = rand(Tv)
                  else
                     V = rand(Tv,lenI,lenJ)
                  end
                  # substitution
                  A[I,J] = V
                  if (isa(I,Integer) && isa(J,Integer))
                     # scalar
                     M[I,J] = V
                  elseif (lenI>1 && lenJ>1) || (lenI==1 && lenJ==1)
                     M[I,J] .= V
                  else
                     M[I,J] = V
                  end
                  X = A[I,J]
                  Y = M[I,J]
                  if lenI==1 && lenJ>1
                     Y = reshape(Y,1,length(Y))
                  elseif lenI>1 && lenJ==1
                     Y = reshape(Y,length(Y),1)
                  end
                  @test issamematrix(A,M)
                  @test issamematrix(X,Y)
                  # update
                  B[I,J] = B[I,J] .+ V
                  if length(size(N[I,J]))==1 && (lenI>1 || lenJ>1)
                     N[I,J] = N[I,J] + vec(V)
                  else
                     N[I,J] = N[I,J] .+ V
                  end
                  X = B[I,J]
                  Y = N[I,J]
                  if lenI==1 && lenJ>1
                     Y = reshape(Y,1,length(Y))
                  elseif lenI>1 && lenJ==1
                     Y = reshape(Y,length(Y),1)
                  end
                  @test issamematrix(B,N)
                  @test issamematrix(X,Y)
               end
            end
         end
      end
   end
end


text:
# # titanic
# Titanic is one of many competitions available through the [Kaggle](https://www.kaggle.com/c/titanic) platform.
# Its goal is to train a classifier which recognizes people who survived the Titanic catastrophe based on some known characteristics, such as age, gender, etc.
# This package provides some possible methods which (at least up to some precision) solve this problem.
#
# ## Input data
# Before we describe the used methods, we need to prepare the data.

# First we need to load the package

using titanic

# Then we can load training data from provided datasets. (can be downloaded on [Kaggle](https://www.kaggle.com/c/titanic/data))

traindata, trainlabels = loadtrain("../data/train.csv")

# Note that this function call performs normalization on the data, so for some classificators later (adaboost), we may turn this procedure off with `normalize=false` keyword parameter.

# Now we do almost the same for testing data.
# Here, Kaggle does not provide labels for the testing dataset and the classificator predictions need to be uploaded to their webpage to receive evaluation.

testdata = loadtest("../data/test.csv")

# The datasets were loaded from `.csv` files into inner `Dataset` structures.
# These provide an useful abstraction, so one does not have to operate on matrices and vectors, but it usually suffices to pass these structures.
# As was said before, the constructor of this structure automacially normalizes the data so they are more suitable for most classifiers.
# Also, many of these methods need to add so-called 'bias' term into the feature vector.
# This is conducted in this structure as well by adding $1$ to the beginning of each point in the data matrix.
#
# The same can be said about `trainlabels`, which are stored inside `Labels` structure.
#
# Both of these structures can return their main field when called as functions (or functors).

traindata()
#-
trainlabels()

# ## Logistic regression
# As first and most basic classifier, we use Logistic regression.
#
# This classifier searches for a vector `w` which determines a hyperplane in the feature space that optimally separates the data points.
# However, the data are most likely not linearly separable, so there will be some some misclassified points.
#
# This search occurs in high-level function

w = regression(traindata, trainlabels; epsilon=1e-4, step=1.0)

# The keyword parameters `epsilon` and `step` are optional and default to values in the example.
# `epsilon` indicates precision for the underlying optimizer and `step` determines coefficient of initial step.
#
# As the method for finding the optimal `w` was implemented gradient descent with adaptive stepsize.
# In each iteration, the gradient descent computes so called *negative log-likelihood*
# $$ E(w) = \sum_{x}{\log(1 + e^{-kxw})} $$
# , where $k$ is the class of $x$, and its derivation
# $$ \frac{\partial E(w)}{\partial w} = - \sum_x\frac{kx}{1 + e^{kxw}} $$
# Using this derivation we find the w with minimal *negative log-likelihood* (or *cross entropy*).
#
# The adaptive step size behaves in a way that when the current step size finds better solution, we double the step in the next iteration.
# On the other hand, when the solution is worse (meaning it "jumped" over the optimum), we use half of it in next iteration (this iteration's result is not even used).
# The speed of convergence could be improved by using some more advanced step size adaptation, like the *Armijo conditions*, but the solution would be (almost) the same.

# When we obtain the optimal `w`, we are ready to clasify the test labels.
#
# Before that we can check classification error on the train dataset.

trainprediction = classify(traindata, w)
#-
using titanic: classificationerror
classificationerror(trainlabels, trainprediction)

# We can see that even on the training set, the error is not zero, so we can conclude that the data are indeed linearly unserparable.

# Then we can predict classes for the test dataset.
# For this we cannot get error so straightforwardly, because we have to upload it to Kaggle website.
# So, the package provides a function for exporting `Labels` into `.csv` file.

testprediction = classify(testdata, w)
#-
savepredictions("./logreg.csv", testprediction)

# From the website we obtain slightly worse error than on the training set, but this is expected and the difference is not that huge.
#
# | training error | test error |
# | :------------: | :--------: |
# | 0.19978        | 0.22967    |
#
# These results are not awesome, but for a basic linear classifier not that bad.

# ## Suppor Vector Machines
# Second, we showcase more complicated linear classifier and that is **SVM**.
# We implemented *soft-margin* variant with different kernel possibilities in function

# *Soft-margin* is used because the data are not linearly separable and thus we introduce penalty $C$ for misclassified points.
# $C$ is optional keyword parameter and defaults to $10$.

# The support vector machine in its base form again searches for some optimal vector `w` which defines the hyperplane.
# However, for optimization it is more convenient to optimize its dual problem
# \begin{align*}
# \alpha = \text{argmax}\quad&\sum_{i = 1}^{N}\alpha_i - \frac{1}{2}\sum_{i,j = 1}^{N} \alpha_i \alpha_j y_i y_j K(x_i, x_j) \\
# \text{subject to}\quad&\sum_{i=1}^N \alpha_i y_i = 0 \\
# & 0 \leq \alpha_i \leq C
# \end{align*}

alpha = svm(traindata, trainlabels; C=10.0)

# Kernels on the other hand are supposed to handle the non-separability as they modify the feature space and possibly making the data more separable.
# They usually take some form of a dot product, as in vanilla svm.
# More generally, each symetric positive semidefinite matrix is a kernel function.
#
# Used kernels:
# - Linear kernel: $$ K(x_i, x_j) = x_i \cdot x_j $$
# - Polynomial kernel of power $p$: $$ K(x_i, x_j) = (1 + x_i \cdot x_j)^p $$
# - Gaussian kernel with variance $\sigma^2$: $$ K(x_i, x_j) = e^{\frac{-||x_i - x_j||_2^2}{2\sigma^2}} $$

# Kernels can be changed with setting the kernel function to one of these functions:

titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.linearkernel(xi, xj)
linearalpha = svm(traindata, trainlabels)
#-
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.polynomialkernel(xi, xj; degree=2)
polynomialalpha = svm(traindata, trainlabels)
#-
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.gaussiankernel(xi, xj; variance=1.0)
gaussianalpha = svm(traindata, trainlabels)

# When we obtain alpha from the dual problem, the classification then goes as follows:
# $$ \text{class}(x) = \text{sign}(\sum_{i=1}^N \alpha_i y_i K(x_i, x)) $$

# Again, we can compute errors on the training set

titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.linearkernel(xi, xj)
classificationerror(trainlabels, classify(traindata, traindata, trainlabels, linearalpha))
#-
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.polynomialkernel(xi, xj; degree=2)
classificationerror(trainlabels, classify(traindata, traindata, trainlabels, polynomialalpha))
#-
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.gaussiankernel(xi, xj; variance=1.0)
classificationerror(trainlabels, classify(traindata, traindata, trainlabels, gaussianalpha))

# As in logistic regression, for getting test error on test data, we need to upload predictions to the Kaggle page.
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.linearkernel(xi, xj)
linearpredictions = classify(testdata, traindata, trainlabels, linearalpha)
savepredictions("./linearsvm.csv", linearpredictions)
#-
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.polynomialkernel(xi, xj; degree=2)
polynomialpredictions = classify(testdata, traindata, trainlabels, polynomialalpha)
savepredictions("./polynomialsvm.csv", polynomialpredictions)
#-
titanic.kernel(xi::Vector{<:Real}, xj::Vector{<:Real}) = titanic.gaussiankernel(xi, xj; variance=1.0)
gaussianpredictions = classify(testdata, traindata, trainlabels, gaussianalpha)
savepredictions("./gaussiansvm.csv", gaussianpredictions)

# Results on training and test data are summarized in table below
#
# | kernel     | training error | test error |
# | :--------- | :------------: | :--------: |
# | linear     | 0.21324        | 0.23445    |
# | polynomial | 0.16498        | 0.24642    |
# | gaussian   | 0.12233        | 0.23924    |
#
# Interestingly, we can observe that while polynomial and gaussian kernels have significantly smaller error on the training set, they perform worse than the linear kernel on the test set.

# ## Adaboost
#
# Last implemented classifier is adaboost.
# Adaboost is actually a multitude of *weak (simple) classifiers* united into one so-called *strong classifier*.
# Each of these weak classifiers classifies a point based on some criterion and returns its decision to the strong classifier, which decides based on weighted sum of all these decisions.
# Weights are computed from classification errors of the weak classifiers.
# Simultaneously, it assigns different weights to data samples based on how hard it is for the weak classifiers to classify them correctly.
# Also, thanks to partial errors of the weak classifiers we can compute upper bound of the classification error, which can be convenient.
#
# Pseudocode for adaboost follows
#
# 1. initialize data weights $$D_1(i) = \frac{1}{N} \qquad \forall i = 1, \dots, N$$

# 2. for t in $1,\dots,T$
#
#     3. find the weak classifier $h_t$ with lowest error $$ h_t = \text{argmin }\epsilon(h) \qquad \epsilon(h) = \sum_{i=1}^N D_t(i)\llbracket y_i \neq  h_t(x_i)\rrbracket $$
#     4. if $\epsilon(h_t) \geq 0.5$ then stop
#     5. compute classifer weight $$\alpha_t = \frac{1}{2}\log(\frac{1 - \epsilon_t}{\epsilon_t}) $$
#     6. adjust data weights $$D_{t+1}(i) = \frac{1}{Z_t}D_t(i) e^{-\alpha_t y_i h_t(x_i)} \quad \forall i = 1, \dots, N \qquad Z_t = \sum_{i=1}^N D_t(i) e^{-\alpha_t y_i h_t(x_i)} $$
#
# $Z_t$ normalizes the weights so they sum to $1$ and also is the already mentioned upper bound on classification error.

# This choice seemed sensible as the features used in data samples are usually not related and adaboost gives the option to classify based on some feature only.
# With this thought, we designed class of weak classifiers that selects one data feature and finds threshold that misclassifies the least amount of samples based only on this feature.
# For each pair feature and threshold there are two weak classifiers, one that classifies all featuers $<$ than threhsold as $1$ and others as $-1$, and second that classifies all features $\geq$ than threshold as $1$ and others as $-1$.

# As this classifier operates on individual features, the data normalization is not needed.
# Thus, the data can be loaded like this

traindata, trainlabels = loadtrain("../data/train.csv"; normalize=false)
#-
testdata = loadtest("../data/test.csv"; normalize=false)

# Training of our classifier is performed by

H = boost(traindata, trainlabels; limit=100)

# The keyword parameter `limit` corresponds to $T$ in the pseudocode and limits the number of weak classifiers constructing the strong one.

# Classification of a point then proceeds as follows
# $$ H(x) = \text{sign}(\sum_{t=1}^T \alpha_t h_t(x)) $$
# Again, we compare errors on train and test datasets.

trainprediction = classify(traindata, H)
classificationerror(trainlabels, trainprediction)
#-
testprediction = classify(testdata, H)
savepredictions("./adaboost.csv", testprediction)

# Again, the training error is quite low, but the test error is very similar to results from other classifiers.
# It is surprising, that it did not work as expected.
# It is possible, that some features are present more in the training set and are more discriminative than in the test set.
# Perhaps, some experiments with what features use to decide would find an answer to this possibility and the algorithm could be improved.
#
# | training error | test error |
# | :------------: | :--------: |
# | 0.16835        | 0.23924    |

# The implementation of the algorithm is sped up by presorting feature spaces, so it is not necessary to iterate over all values in range, but only until the error starts increasing again
# If this feature was removed, we would get the same results but it would be a lot slower.
#
# Another speedup could be achieved by decreasing the `limit`.
# When we reduce it to $10$, the learning is approximately $10$ times faster, but the error on training set increased only by $0.2$.
# If speed was prefered to accuracy, this is another area to take advantage of.

# ## Conclusion
#
# | method              | training error | test error |
# | :------------------ | :------------: | :--------: |
# | logistic regression | 0.19978        | 0.22967    |
# | kernel svm          | 0.21324        | 0.23445    |
# | adaboost            | 0.16835        | 0.23924    |
#
# From the results we can see, that logistic regression is the best of the algorithms on test dataset, even though it does not have the lowest error on training data.
# On the other hand, adaboost has the lowest training error, but is the worst on testing data.
# This could mean, that it is overfitted or that the test data are more different from the training set than expected.
# The kernel svm generalizes quite well as the training error and test error are similar.
# It is disappointing that the kernels with lower training error performed worse on the test data.
# Possibly, it could be improved by finding better exponent or variance.
#
# In terms of speed, the worst is probably adaboost, because it does not have an easily optimizable form as gradient descent or quadratic optimization problem.
# We need to search for the best classifier over all features and that makes it slower.
#
# In conclusion, the results cannot be considered as a success as all the methods classified with quite high errors on both training and testing datasets.
# Some more powerful classifiers, such as neural nets, would probably work better or some smarter data preprocessing could be investigated.
# For example, some features could be omitted as unimportant and better results could be achieved.


text:
"""Load the stream graph from figure 1."""
function load_stream_graph_1()
	s=StreamGraph("stream-graph 1", Intervals([(0.0,10.0)]))
	record!(s,0.0,10.0,"a")
	record!(s,0.0,4.0,"b")
	record!(s,5.0,10.0,"b")
	record!(s,4.0,9.0,"c")
	record!(s,1.0,3.0,"d")
	record!(s,1.0,3.0,"a","b")
	record!(s,7.0,8.0,"a","b")
	record!(s,4.5,7.5,"a","c")
	record!(s,6.0,9.0,"b","c")
	record!(s,2.0,3.0,"b","d")
	s
end

"""Load the stream graph from figure 2a."""
function load_stream_graph_2a()
	s=StreamGraph("stream-graph 2a", Intervals([(0.0,3.0)]))
	record!(s,0.0,1.0,"a")
	record!(s,0.0,3.0,"b")
	record!(s,0.0,1.0,"c")
	record!(s,2.0,3.0,"c")
	record!(s,0.0,1.0,"a","b")
	record!(s,0.0,1.0,"b","c")
	record!(s,2.0,3.0,"b","c")
	s
end

"""Load the stream graph from figure 2b."""
function load_stream_graph_2b()
	s=StreamGraph("stream-graph 2b", Intervals([(0.0,3.0)]))
	record!(s,0.0,1.0,"a")
	record!(s,0.0,3.0,"b")
	record!(s,1.0,3.0,"c")
	record!(s,0.0,1.0,"a","b")
	record!(s,1.0,3.0,"b","c")
	s
end

"""Load the stream graph from figure 6a."""
function load_stream_graph_6a()
	s=StreamGraph("stream-graph 6a", Intervals([(0.0,4.0)]))
	record!(s,0.0,2.0,"a")
	record!(s,0.0,4.0,"b")
	record!(s,2.0,4.0,"c")
	record!(s,0.0,2.0,"a","b")
	record!(s,2.0,4.0,"b","c")
	s
end

"""Load the stream graph from figure 6b."""
function load_stream_graph_6b()
	s=StreamGraph("stream-graph 6b", Intervals([(0.0,4.0)]))
	record!(s,0.0,0.0,"a")
	record!(s,0.0,4.0,"b")
	record!(s,0.0,4.0,"c")
	record!(s,0.0,4.0,"b","c")
	s
end

"""Load the stream graph from figure 7."""
function load_stream_graph_7()
	s=StreamGraph("stream-graph 7", Intervals([(0.0,10.0)]))
	record!(s,0.0,10.0,"a")
	record!(s,0.0,10.0,"b")
	record!(s,0.0,5.0,"c")
	record!(s,0.0,5.0,"b","c")
	record!(s,5.0,10.0,"a","b")
	s
end

# StreamGraph tests
s=load_stream_graph_1()
@test length(s.V)==4
@test duration(s)==10.0
@test times(s,"a")==Intervals([(0.0,10.0)])
@test times(s,"b")==Intervals([(0.0,4.0),(5.0,10.0)])
@test times(s,"c")==Intervals([(4.0,9.0)])
@test times(s,"d")==Intervals([(1.0,3.0)])
@test times(s,"a","b")==Intervals([(1.0,3.0),(7.0,8.0)])
@test times(s,"b","a")==Intervals([(1.0,3.0),(7.0,8.0)])
@test times(s,"a","c")==Intervals([(4.5,7.5)])
@test times(s,"c","a")==Intervals([(4.5,7.5)])
@test times(s,"b","c")==Intervals([(6.0,9.0)])
@test times(s,"c","b")==Intervals([(6.0,9.0)])
@test times(s,"b","d")==Intervals([(2.0,3.0)])
@test times(s,"d","b")==Intervals([(2.0,3.0)])
@test times(s,"a","d")==Intervals([])
@test Set(nodes(s,0.0))==Set(["a","b"])
@test Set(nodes(s,2.0))==Set(["a","b","d"])
@test Set(nodes(s,4.0))==Set(["a","b","c"])
@test Set(nodes(s,4.00001))==Set(["a","c"])
@test Set(nodes(s,1.5,2.5))==Set(["a","b","d"])
@test Set(nodes(s,0.0,4.0))==Set(["a","b"])
@test Set(nodes(s,1.5,9.0))==Set(["a"])
@test Set(nodes(s,5.0,10.0))==Set(["a","b"])
@test Set(nodes(s,0.0,10.0))==Set(["a"])
@test Set(nodes(s,0.0,0.0))==Set(["a","b"])
#@test Set(links(s,0.0,3.0))==Set([Link("a to b",Intervals([(1.0,3.0)]),"a","b",1.0),Link("b to d",Intervals([(2.0,3.0)]),"b","d",1.0)])
@test coverage(s)==0.65
@test number_of_nodes(s)==2.6
@test number_of_links(s)==1.0
@test node_duration(s)==6.5
@test isapprox(link_duration(s),1.6666666;atol=0.00001)
@test uniformity(s)==0.39285714285714285
@test density(s)==0.45454545454545453
@test density(s,2.0)==2.0/3.0
@test density(s,"a","b")==1.0/3.0
@test density(s,"b","a")==1.0/3.0
@test density(s,"b","d")==1.0/2.0
@test density(s,"d","b")==1.0/2.0
@test density(s,"d")==0.25
Na=neighborhood(s,"a")
@test length(Na)==2
@test Na["b"].presence==Intervals([(1.0,3.0),(7.0,8.0)])
@test Na["c"].presence==Intervals([(4.5,7.5)])
@test degree(s,"a")==0.6
@test neighborhood(s,"a",0.0)==[]
@test neighborhood(s,"a",1.0)==["b"]
@test Set(neighborhood(s,"a",7.4))==Set(["b","c"])
@test Set(neighborhood(s,"b",7.4))==Set(["a","c"])
@test Set(neighborhood(s,"d",2.4))==Set(["b"])
@test degree(s,"a",0.0)==0
@test degree(s,"a",1.0)==1
@test degree(s,"a",7.4)==2
@test degree(s,"b",7.4)==2
@test degree(s,"d",2.4)==1
@test degree(s,0.0)==0.0
@test degree(s,1.0)==0.5
@test degree(s,2.0)==1.0
@test degree(s,6.0)==1.0
@test degree(s,7.1)==6.0/4.0
@test degree(s)==0.5
@test times(s)==[0.0,1.0,2.0,3.0,4.0,4.5,5.0,6.0,7.0,7.5,8.0,9.0,10.0]
@test time_clustering(s)==0.11538461538461539
@test clustering(s)==0.037500000000000006

s2a=load_stream_graph_2a()
s2b=load_stream_graph_2b()
@test isapprox(number_of_nodes(s2a),2.0;atol=0.00001)
@test isapprox(number_of_nodes(s2b),2.0;atol=0.00001)
@test number_of_links(s2a)==1
@test number_of_links(s2b)==1
@test density(s2a)==0.75
@test density(s2b)==1.0
@test times(s2a)==[0.0,1.0,2.0,3.0]
@test times(s2b)==[0.0,1.0,3.0]

s6a=load_stream_graph_6a()
s6b=load_stream_graph_6b()
@test density(s6a)==1.0
@test density(s6b)==1.0
@test degree(s6a,"a")==0.5
@test degree(s6a,"b")==1.0
@test degree(s6a,"c")==0.5
@test average_node_degree(s6a)==0.75

s7=load_stream_graph_7()
@test number_of_nodes(s7)==2.5
@test degree(s7,"a")==0.5
@test degree(s7,"b")==1.0
@test degree(s7,"c")==0.5
@test isapprox(average_node_degree(s7),0.7;atol=0.000001)
@test node_duration(s7)==25.0/3.0
@test times(s7)==[0.0,5.0,10.0]
@test node_contribution(s7,2.5)==1.0
@test degree(s7,2.5)==2.0/3.0
@test node_contribution(s7,7.5)==2.0/3.0
@test degree(s7,7.5)==2.0/3.0
@test average_time_degree(s7)==2.0/3.0
@test time_clustering(s7,2.5)==0.0
@test time_clustering(s7,5.0)==0.0
@test time_clustering(s7)==0.0
@test clustering(s7)==0.0


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -37.92418043260588 <= q <= 19.457860872935036)

                     add_NL_constraint(m, :(log(1 + exp(-0.5948345021955626 + 0.7909419235855366*log(1 + exp(0.16572281076834727 + 0.7701313689304001*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.3960198610694281*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.006655112331988011*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.5604852483433853*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.5611932928423773*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.9350053029495258*log(1 + exp(0.6292157071886879 + 0.7447123643459479*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.8071073342022652*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.7192221993944883*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1911667575027427*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.30728485701247665*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.6696024550400765*log(1 + exp(-0.5341808090129341 + 0.4032770914369115*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.6061337759397625*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.5534397384266061*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.16208042832895408*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.7376214833764658*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.9733717899660395*log(1 + exp(0.43846326345255404 + -0.02862602361460853*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.5029397949622734*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.20686891235628035*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.9777195341190259*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.2597965376833278*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.07261217388131413*log(1 + exp(0.02496868545483233 + -0.04485855107115411*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.47774406442351136*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.5903487696420049*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1719483424118926*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.8023455576856944*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))))) + log(1 + exp(-0.7044745597894595 + -0.85619817917634*log(1 + exp(0.16572281076834727 + 0.7701313689304001*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.3960198610694281*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.006655112331988011*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.5604852483433853*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.5611932928423773*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.37830231556965854*log(1 + exp(0.6292157071886879 + 0.7447123643459479*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.8071073342022652*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.7192221993944883*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1911667575027427*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.30728485701247665*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.12156626196118481*log(1 + exp(-0.5341808090129341 + 0.4032770914369115*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.6061337759397625*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.5534397384266061*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.16208042832895408*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.7376214833764658*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.3744591738356289*log(1 + exp(0.43846326345255404 + -0.02862602361460853*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.5029397949622734*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.20686891235628035*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.9777195341190259*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.2597965376833278*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.4985780257939494*log(1 + exp(0.02496868545483233 + -0.04485855107115411*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.47774406442351136*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.5903487696420049*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1719483424118926*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.8023455576856944*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))))) + log(1 + exp(-0.9683275825894269 + -0.7637278723892034*log(1 + exp(0.16572281076834727 + 0.7701313689304001*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.3960198610694281*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.006655112331988011*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.5604852483433853*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.5611932928423773*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.9143661573707975*log(1 + exp(0.6292157071886879 + 0.7447123643459479*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.8071073342022652*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.7192221993944883*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1911667575027427*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.30728485701247665*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.7824888723979511*log(1 + exp(-0.5341808090129341 + 0.4032770914369115*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.6061337759397625*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.5534397384266061*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.16208042832895408*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.7376214833764658*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.6866950641423624*log(1 + exp(0.43846326345255404 + -0.02862602361460853*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.5029397949622734*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.20686891235628035*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.9777195341190259*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.2597965376833278*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.09633196432882363*log(1 + exp(0.02496868545483233 + -0.04485855107115411*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.47774406442351136*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.5903487696420049*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1719483424118926*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.8023455576856944*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))))) + log(1 + exp(-0.5143076265809325 + -0.11312975790997548*log(1 + exp(0.16572281076834727 + 0.7701313689304001*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.3960198610694281*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.006655112331988011*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.5604852483433853*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.5611932928423773*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.7215302240627146*log(1 + exp(0.6292157071886879 + 0.7447123643459479*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.8071073342022652*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.7192221993944883*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1911667575027427*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.30728485701247665*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.8436266175918532*log(1 + exp(-0.5341808090129341 + 0.4032770914369115*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.6061337759397625*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.5534397384266061*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.16208042832895408*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.7376214833764658*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.8029468027318871*log(1 + exp(0.43846326345255404 + -0.02862602361460853*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.5029397949622734*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.20686891235628035*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.9777195341190259*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.2597965376833278*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.07363029133667531*log(1 + exp(0.02496868545483233 + -0.04485855107115411*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.47774406442351136*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.5903487696420049*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1719483424118926*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.8023455576856944*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))))) + log(1 + exp(-0.8626839052591899 + -0.6726505975051698*log(1 + exp(0.16572281076834727 + 0.7701313689304001*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.3960198610694281*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.006655112331988011*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.5604852483433853*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.5611932928423773*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + -0.7431517048783043*log(1 + exp(0.6292157071886879 + 0.7447123643459479*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.8071073342022652*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.7192221993944883*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1911667575027427*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.30728485701247665*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.03587478002549904*log(1 + exp(-0.5341808090129341 + 0.4032770914369115*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + -0.6061337759397625*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + 0.5534397384266061*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.16208042832895408*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.7376214833764658*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.7126331133171049*log(1 + exp(0.43846326345255404 + -0.02862602361460853*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.5029397949622734*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.20686891235628035*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.9777195341190259*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.2597965376833278*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))) + 0.7945114361960401*log(1 + exp(0.02496868545483233 + -0.04485855107115411*log(1 + exp(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2]))) + 0.47774406442351136*log(1 + exp(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2]))) + -0.5903487696420049*log(1 + exp(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2]))) + -0.1719483424118926*log(1 + exp(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2]))) + -0.8023455576856944*log(1 + exp(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))))))) - $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:
# SPDX-License-Identifier: X11
# 2020-12-09
# Day 5, Part 2

const TOTAL = 128 * 8

function main()
  passes = readlines()
  seats = falses(TOTAL)

  for p ∈ passes
    seats[seatid(p)] = true
  end

  for i ∈ 2:TOTAL - 1
    seats[i - 1] && !seats[i] && seats[i + 1] && (println(i); return)
  end
end

function seatid(str::AbstractString)
  row = -1
  s = 0
  e = 127
  for (i, c) ∈ enumerate(str[1:7])
    c == 'F' ? e -= 2^(7 - i) : s += 2^(7 - i)
  end
  @assert s == e
  row = s

  col = -1
  s = 0
  e = 7
  for (i, c) ∈ enumerate(str[8:10])
    c == 'L' ? e -= 2^(3 - i) : s += 2^(3 - i)
  end
  @assert s == e
  col = s

  return row * 8 + col
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
#
# LatticeRule32
#
struct LatticeRule32{s} <: AbstractLatticeRule{s}
    z::Vector{UInt32} # generating vector
    n::Int64 # max number of points in the lattice rule
end

# default lattice rule type
const LatticeRule = LatticeRule32

# access max number of points in the lattice
Base.length(lattice_rule::LatticeRule32) = lattice_rule.n

# uinttype
uinttype(::LatticeRule32) = UInt32

"""
    LatticeRule32(z, s, n)
    LatticeRule32(z, s)
    LatticeRule32(z)

Returns a rank-1 lattice rule in `s` dimensions with generating vector `z` and at most `n` points.

When no maximum number of points `n` is provided, we assume `n = typemax(UInt32) = 2^32 - 1`. When no number of dimensions `s` is provided, we assume `s = length(z)`.

!!! info

    Technically, we return an extensible lattice sequence where the `k`-th point is transformed using the gray coded radical inverse function. This has the advantage that we can add points to the lattice without changing the already computed points.

More generating vectors can be found online [here](https://web.maths.unsw.edu.au/~fkuo/lattice/index.html) or [here](https://people.cs.kuleuven.be/~dirk.nuyens/qmc-generators/).

# Examples
```jldoctest; setup = :(using LatticeRules; import Random; Random.seed!(1))
julia> lattice_rule = LatticeRule32([UInt32(1), UInt32(5)], 2, 8) # Fibonacci lattice
LatticeRule32{2}

julia> getpoint(lattice_rule, 2)
2-element Array{Float64,1}:
 0.25
 0.25

```
See also: [`getpoint`](@ref), [`ShiftedLatticeRule32`](@ref)
"""
LatticeRule32(z::Vector{UInt32}) = LatticeRule32(z, length(z)) # specify generating vector

# specify generating vector and number of dimensions
LatticeRule32(z::Vector{UInt32}, s::Integer) = LatticeRule32(z, s, typemax(UInt32) + 1)

# specify generating vector, number of dimensions and maximum number of points
function LatticeRule32(z::Vector{UInt32}, s::Integer, n::Integer)
    s > 0 || throw(ArgumentError("number of dimensions s must be larger than 0"))
    s ≤ length(z) || throw(ArgumentError("number of dimensions s must be less than or equal to the length of the generating vector z"))
    n > 0 || throw(ArgumentError("maximum number of points n must be larger than 0"))
    n ≤ typemax(UInt32) + 1 || throw(ArgumentError("maximum number of points n must be less than or equal to 2^32, consider implementing a LatticeRule64 type"))
    LatticeRule32{s}(view(z, 1:s), n)
end

"""
    LatticeRule32(file, s, n)
    LatticeRule32(file, s)
    LatticeRule32(file)

Returns a rank-1 lattice rule in `s` dimensions with generating vector `z` read from the file `file` and with at most `n` points.

# Examples
```jldoctest; setup = :(using LatticeRules)
julia> z_file = K_3600_32_file;

julia> lattice_rule = LatticeRule32(z_file, 16)
LatticeRule32{16}

julia> getpoint(lattice_rule, 123)
16-element Array{Float64,1}:
 0.8671875
 0.9609375
 0.6015625
 0.8984375
 0.6484375
 0.6328125
 0.3203125
 0.2890625
 0.0234375
 0.1015625
 0.7890625
 0.0703125
 0.6953125
 0.0234375
 0.1171875
 0.0859375

```
See also: [`getpoint`](@ref), [`ShiftedLatticeRule32`](@ref)
"""
LatticeRule32(file::AbstractString) = LatticeRule32(read32(file)) # specify file containing generating vector

# specify file containting generating vector and number of dimensions
LatticeRule32(file::AbstractString, s::Integer) = LatticeRule32(read32(file), s)

# specify file containing generating vector, number of dimensions and maximum number of points
LatticeRule32(file::AbstractString, s::Integer, n::Integer) = LatticeRule32(read32(file), s, n)

"""
    LatticeRule32(s)

Returns a rank-1 lattice rule in `s` dimensions that uses a default generating vector with order-2 weights.

# Examples
```jldoctest; setup = :(using LatticeRules)
julia> lattice_rule = LatticeRule32(16)
LatticeRule32{16}

julia> getpoint(lattice_rule, 123)
16-element Array{Float64,1}:
 0.8671875
 0.5390625
 0.6015625
 0.3671875
 0.6796875
 0.8203125
 0.3046875
 0.8515625
 0.7109375
 0.6328125
 0.5703125
 0.2578125
 0.6953125
 0.0390625
 0.2421875
 0.4453125

```
See also: [`getpoint`](@ref), [`ShiftedLatticeRule32`](@ref)
"""
function LatticeRule32(s::Integer) # specify number of dimensions only
    s ≤ 3600 || throw(ArgumentError("number of dimensions s must be less than or equal to 3600, please supply your own generating vector z"))
    s ≤ 250 ? LatticeRule32(CKN_250_20, s, 2^20) : LatticeRule32(K_3600_32, s)
end

# in-place version of unsafe_getpoint (with 0 memory allocations)
@inline function unsafe_getpoint!(x::Vector{<:AbstractFloat}, lattice_rule::LatticeRule32, k::UInt32)
    ϕ_k = reversebits(k) * 2.0^(-32) # gray coded radical inverse function in base 2
    @inbounds for i in 1:length(x)
        x[i] = ϕ_k * lattice_rule.z[i]
        x[i] -= floor(x[i]) # mod 1
    end
    x
end

# fancy printing
Base.show(io::IO, lattice_rule::LatticeRule32{s}) where s = print(io, string("LatticeRule32{", s, "}"))


text:
using GLMNet

function logrange(min, max, num)
  shift = 1.0 - min
  exp(linrange(log(min + shift), log(max + shift), num)) - shift
end

macro validate_and_init_for_stability_selection()
  esc(quote
    @assert 0.0 <= alpha <= 1.0
    X = float64(X)
    y = float64(y)
    N, P = size(X)
    halfsize = ifloor(N / 2)
    counts = zeros(Int64, P, nlambda)
    updatecounts = zeros(Int64, nlambda)
    Xn, xcolnorms = normalize_columns(X)
    Xnw = zeros(Float64, N, P)
    indices = collect(1:N)
  end)
end

function normalize_columns(A)
  acolnorms = map((j) -> norm(A[:,j]), 1:size(A, 2))
  An = broadcast(/, A, acolnorms')
  return (An, acolnorms)
end

macro update_counts(counts, updatecounts, r)
  quote
    for li in 1:length($(r).lambda)
      $(updatecounts)[li] += 1
      for bi in 1:P
        if $(r).betas[bi, li] != 0.0
          $(counts)[bi, li] += 1
        end
      end
    end
  end
end

# Stability selection in the spirit of Meinshausen&Buhlman but with the Complementary Pairs
# variant described in
#   Rajen Shah and Richard J. Samworth, "Variable selection with error control: Another
#   look at Stability Selection", arxiv.org, 2011, http://arxiv.org/pdf/1105.5578.pdf
#
#   selfreqs = stability_selection(X, y; options...)
#
# Inputs:
#  X is the N*P design matrix
#  y is the N*1 variable to predict
#
# Output:
#  selfreqs is the probability that a feature was selected during the first
#             maxit iterations of the GLMNet LASSO algorithm. The features are
#             randomly reweighted in [alpha, 1.0] and the probability is estimated
#             by numbootstraps samples.
#
# Options:
#  maxit is the maximum number of iterations of the LASSO
#  alpha is the minimum weight of a feature
#  numbootstraps is the number of bootstrap samples
#  parallel is true if it should run the
function stability_selection(X::AbstractMatrix, y::AbstractVector;
                             numbootstraps = 50,
                             maxit = 500,
                             alpha = 0.20,
                             nlambda = 25,
                             lambdas = [],
                             options...)

  @validate_and_init_for_stability_selection

  # Run one glmnet fit with only a few lambdas so we can extract their min and max.
  if length(lambdas) == 0
    res = glmnet(Xn, y; nlambda = 3, options...)
    lambdamax, lambdamin = res.lambda[1], res.lambda[end]
    lambdas = reverse(logrange(lambdamin, lambdamax, nlambda))
  end

  # Now take numbootstraps (parallel) samples
  map(1:numbootstraps) do i

    # Reweight each covariate in each round. We can do this via the penalty_factors
    # option instead?!
    weights = alpha + (1.0 - alpha) * rand(P)
    for j in 1:P
      for i in 1:N
        Xnw[i, j] = Xn[i, j] / weights[j]
      end
    end

    # Randomly split sample in two sets
    shuffle!(indices)
    idxs1 = indices[1:halfsize]
    idxs2 = indices[(halfsize+1):(2*halfsize)] # Can't take to end since they need to be same size

    # Run GLMNet on each of the samples and count which variables are selected
    r1 = glmnet(Xnw[idxs1,:], y[idxs1]; lambda = lambdas, standardize = false, intercept = true, maxit = maxit)
    r2 = glmnet(Xnw[idxs2,:], y[idxs2]; lambda = lambdas, standardize = false, intercept = true, maxit = maxit)

    @update_counts(counts, updatecounts, r1)
    @update_counts(counts, updatecounts, r2)
  end

  # Complementary Pairs SS just sums all counts per feature
  #freqs = sum(counts, 2) / sum(updatecounts)

  # Original muhleinbein paper used the max??
  freqs = broadcast(*, counts, (1 ./ updatecounts)')
  # Return the max frequency achieved per feature
  freqs = maximum(freqs, 2)

  freqs[:]
end

# This implements the treshold (tau) selection procedure as described
# at the end of page 10 in the Rajan&Shamworth paper.
#
# Input:
#  P is the number of original features
#  l is max number of low probability features selected, on average
#  B is number of bootstraps, i.e. 2*B is the number of subsets used
function select_treshold(P, l = 0.50, B = 50)
  q = sqrt(0.8*l*P)
  minD(q/P, B)
end

function minD(theta, B, rs = [-0.5, -0.25])
  minimum(vcat(ones(B), tailprobabilities(theta^2, B, rs[1]), tailprobabilities(theta, B, rs[2])))
end

# Calculate the tail probabilities (based on r-concave funcs) for each
# relevant tau. Based on Shah and Shamworth paper and their R code:
#   http://www.statslab.cam.ac.uk/~rds37/papers/r_concave_tail.R
function tailprobabilities(eta, B, r)
  MAXa = 100000
  MINa = 0.00001
  s = -1/r
  etaB = eta * B
  k_start = iceil(2 * etaB) + 1
  output = ones(B)
  if k_start > B
    return(output)
  end

  a_vec = MAXa * ones(B)

end

function indover(v, treshold = 0.70)
  collect(1:length(v))[v .> treshold]
end

N, P, k = (50, 200, 5)
coefs = shuffle(vcat(100*rand(k), zeros(P-k)))
X = randn(N, P)
errors = 1*randn(N)
y = X * coefs .+ errors

tic()
selfreqs = stability_selection(X, y; maxit = 500, alpha = 0.20, nlambda = 20)
t = toq()
reverse(sort(coefs))[1:k]
reverse(sort(selfreqs))[1:(2k)]
indover(coefs, 0.0)
indover(selfreqs, min(N/P, 1.0)) # Shah&Samworth recommends N/P as the treshold, I modified to min(N/P, 1.0) since frequencies <= 1.0


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 04-2019
# Define operations for 1D spaces
#--------------------------------------------------------------------

using LinearAlgebra
import Base: maximum, minimum, range, size, ndims
export order, identity, incomingboundary, outgoingboundary

ndims(space::S)   where {S<:Cardinal{Tag, N, T}} where {Tag, N, T} = 1
minimum(space::S) where {S<:Cardinal{Tag, N, T}} where {Tag, N, T} = space.min
maximum(space::S) where {S<:Cardinal{Tag, N, T}} where {Tag, N, T} = space.max
order(space::S)   where {S<:Cardinal{Tag, N, T}} where {Tag, N, T} = N-1
size(space::S)    where {S<:Cardinal{Tag, N, T}} where {Tag, N, T} = N
range(space::S)   where {S<:Cardinal{Tag, N, T}} where {Tag, N, T} = (minimum(space), maximum(space))

function Base. identity(space::S)::Operator{S} where {S<:Cardinal{Tag, N, T}} where {Tag, N, T}
    return Operator(space, Diagonal(ones(T, size(space))))
end

function outgoingboundary(space::S)::Operator{S} where {S<:Cardinal{Tag, N, T}} where {Tag, N, T}
    return Operator(space, Diagonal([1, zeros(T, order(space))...]))
end

function incomingboundary(space::S)::Operator{S} where {S<:Cardinal{Tag, N, T}} where {Tag, N, T}
    return Operator(space, Diagonal([zeros(T, order(space))..., 1]))
end

function derivative(space::S)::Operator{S} where {S<:Cardinal{Tag, N, T}} where {Tag, N, T}
    return Operator(space, [derivative(space, i, j) for i in 1:size(space), j in 1:size(space)])
end

function integral(space::S)::Operator{S} where {S<:Cardinal{Tag, N, T}} where {Tag, N, T}
    return Operator(space, Diagonal([integral(space, i) for i in 1:size(space)]))
end

function Field(space::S, umap::Function)::Field{S} where {S<:Cardinal{Tag, N, T}} where {Tag, N, T}
    value = zeros(T, size(space))
    for index in CartesianIndices(value)
        value[index] = umap(collocation(space, index.I[1]))
    end
    return Field(space, value)
end


text:
using Optim, StatsBase

"""
    ols_cdf(parameters, experimental, experimental_ecdf, ν0, I)

Compute ordinary least squares comparing the experimental ecdf with the
theoretical cdf, calculated with the nmr parameters, the spin (I), and the
Larmor frequency (ν0) at each x-value in the experimental data

"""
function ols_cdf(
    parameters::Quadrupolar,
    exp::Array{Float64, 2},
    exp_ecdf::Array{Float64, 1},
    ν0::Float64;
    I::Int64 = 3,
    samples::Int64 = 1_000_000,
    transitions::UnitRange{Int64} = 1:(2*I),
    range::Tuple{Float64,Float64} = (exp[1, 1], exp[end, 1]),
)
    #range = (findfirst(x -> range[1] < x, exp[:, 1]),
    #    findlast(x -> range[2] > x, exp[:, 1]) - 1)
    th_ecdf = ecdf(estimate_powder_pattern(parameters, samples, ν0, I,
        transitions = transitions)).(exp[:, 1])
    th_ecdf .-= th_ecdf[1]
    th_ecdf ./= th_ecdf[end]
    return sum((exp_ecdf .- th_ecdf) .^ 2)
end

function ols_cdf(
    parameters::Missing,
    exp::Array{Float64, 2},
    exp_ecdf::Array{Float64, 1},
    ν0::Float64;
    I::Int64 = 3,
    samples::Int64 = 1_000_000,
    transitions::UnitRange{Int64} = 1:(2*I),
)
    return size(exp)[1] / 10
end

function ols_cdf(
    parameters::ChemicalShift,
    exp::Array{Float64, 2},
    exp_ecdf::Array{Float64, 1};
    samples::Int64 = 1_000_000,
)
    th_ecdf = ecdf(estimate_powder_pattern(parameters, samples)).(exp[:, 1])
    th_ecdf .-= th_ecdf[1]
    th_ecdf ./= th_ecdf[end]
    return sum((exp_ecdf .- th_ecdf) .^ 2)
end

function ols_cdf(
    parameters::Missing,
    exp::Array{Float64, 2},
    exp_ecdf::Array{Float64, 1};
    samples::Int64 = 1_000_000,
)
    return size(exp)[1] / 10
end

function transform_params(x::Array{Float64}, T::DataType)
    if T == Quadrupolar
        x = x .^ 2
    elseif T == ChemicalShift
        for i = 1:length(x) ÷ 7
            x[7 * (i - 1) + 2] ^= 2
            x[7 * (i - 1) + 4] ^= 2
            x[7 * (i - 1) + 5] ^= 2
            x[7 * (i - 1) + 6] ^= 2
            x[7 * (i - 1) + 7] ^= 2
        end
    end
    return x
end

"""
    fit_quadrupolar(experimental; sites, iters, options, method)

Fit the NMR parameters, assuming a normal distribution and using the specified
optimization method (currently implementer are Nelder-Mead and Simulated
Annealing) for with Optim options (options) for iterations, timeout, etc, with
the loss function of ordinary least squares of the CDFs, assuming a specified
number of sites (sites, defaults to 1)

"""
function fit_quadrupolar(
    experimental::Array{Float64, 2};
    sites::Int64 = 1,
    iters::Int64 = 1000,
    options = Optim.Options(iterations = iters),
    method::String = "NelderMead",
    I::Int64 = 3,
    samples::Int64 = 1_000_000,
    starting_values = get_quadrupolar_starting_values(sites),
    transitions::UnitRange{Int64} = 1:(2*I),
    range::Tuple{Float64,Float64} = (experimental[1, 1], experimental[end, 1]),
)
    experimental_ecdf = get_experimental_ecdf(experimental)
    ν0 =  get_ν0(experimental, experimental_ecdf)

    if method == "SimulatedAnnealing"
        upper_bounds, lower_bounds = zeros(5 * sites), zeros(5 * sites)
        upper_bounds[1:5:end] .= 9  # Qcc
        upper_bounds[2:5:end] .= 1  # σQcc
        upper_bounds[3:5:end] .= 1  # η
        upper_bounds[4:5:end] .= 1  # ση
        upper_bounds[5:5:end] .= 1  # weights
        result = optimize(
            x -> ols_cdf(  # objective function
                Quadrupolar(x),
                experimental,
                experimental_ecdf,
                ν0,
                I = I,
                samples = samples,
                transitions = transitions,
                range = range,
            ),
            lower_bounds,
            upper_bounds,
            starting_values,
            SAMIN(),
            options,
        )
    elseif method == "ParticleSwarm"
        upper_bounds, lower_bounds = zeros(5 * sites), zeros(5 * sites)
        upper_bounds[1:5:end] .= 9  # Qcc
        upper_bounds[2:5:end] .= 1  # σQcc
        upper_bounds[3:5:end] .= 1  # η
        upper_bounds[4:5:end] .= 1  # ση
        upper_bounds[5:5:end] .= 1  # weights
        result = optimize(
            x -> ols_cdf(  # objective function
                Quadrupolar(x),
                experimental,
                experimental_ecdf,
                ν0,
                I = I,
                samples = samples,
                transitions = transitions,
            ),
            lower_bounds,
            upper_bounds,
            starting_values,
            ParticleSwarm(),
            options,
        )
    elseif method == "NelderMead"
        result = optimize(
            x -> ols_cdf(
                Quadrupolar(transform_params(x, Quadrupolar)),
                experimental,
                experimental_ecdf,
                ν0,
                I = I,
                samples = samples,
                transitions = transitions,
                range = range,
            ),
            starting_values,
            options,
        )
    else
        println("Unknown method $(method). Try NelderMead, SimulatedAnnealing, or ParticleSwarm")
    end
    return result
end

function fit_chemical_shift(
    experimental::Array{Float64, 2};
    sites::Int64 = 1,
    iters::Int64 = 1000,
    options = Optim.Options(iterations = iters),
    method::String = "NelderMead",
    samples::Int64 = 1_000_000,
    starting_values = get_chemical_shift_starting_values(sites),
    range::Tuple{Float64,Float64} = (experimental[1, 1], experimental[end, 1]),
)
    experimental_ecdf = get_experimental_ecdf(experimental)

    if method == "SimulatedAnnealing"
        upper_bounds, lower_bounds = zeros(7 * sites), zeros(7 * sites)
        upper_bounds[1:7:end] .= 4000
        upper_bounds[2:7:end] .= 800
        upper_bounds[3:7:end] .= 4000
        upper_bounds[4:7:end] .= 400
        upper_bounds[5:7:end] .= 1
        upper_bounds[6:7:end] .= 1
        upper_bounds[7:7:end] .= 1
        lower_bounds[1:7:end] .= -4000
        lower_bounds[2:7:end] .= 0.0000001
        lower_bounds[3:7:end] .= -4000
        lower_bounds[4:7:end] .= 0.0000001
        lower_bounds[5:7:end] .= 0.0000001
        lower_bounds[6:7:end] .= 0.0000001
        lower_bounds[7:7:end] .= 0.0000001
        result = optimize(
            x -> ols_cdf(  # objective function
                ChemicalShift(x),
                experimental,
                experimental_ecdf,
                samples = samples,
            ),
            lower_bounds,
            upper_bounds,
            starting_values,
            SAMIN(),
            options,
        )
    elseif method == "ParticleSwarm"
        upper_bounds, lower_bounds = zeros(7 * sites), zeros(7 * sites)
        upper_bounds[1:7:end] .= 4000
        upper_bounds[2:7:end] .= 800
        upper_bounds[3:7:end] .= 4000
        upper_bounds[4:7:end] .= 400
        upper_bounds[5:7:end] .= 1
        upper_bounds[6:7:end] .= 1
        upper_bounds[7:7:end] .= 1
        lower_bounds[1:7:end] .= -4000
        lower_bounds[2:7:end] .= 0.0000001
        lower_bounds[3:7:end] .= -4000
        lower_bounds[4:7:end] .= 0.0000001
        lower_bounds[5:7:end] .= 0.0000001
        lower_bounds[6:7:end] .= 0.0000001
        lower_bounds[7:7:end] .= 0.0000001
        result = optimize(
            x -> ols_cdf(  # objective function
                ChemicalShift(x),
                experimental,
                experimental_ecdf,
                samples = samples,
            ),
            lower_bounds,
            upper_bounds,
            starting_values,
            ParticleSwarm(),
            options,
        )
    elseif method == "NelderMead"
        result = optimize(
            x -> ols_cdf(
                ChemicalShift([Normal(x[1], x[2] ^ 2)], [Normal(x[3], x[4] ^ 2)],
                    [truncated(Normal(x[5] ^ 2, x[6] ^ 2), 0, 1)], [1.0]),
                experimental,
                experimental_ecdf,
                samples = samples,
            ),
            starting_values,
            options,
        )
    end
    return result
end