text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 6a41e9d0-c58c-11eb-078e-c36bf22a38d2
md"""
# Support Vector Machines

[Reference 1](https://scikit-learn.org/stable/modules/svm.html)\
[Reference 2](https://en.wikipedia.org/wiki/Support_vector_machine)\
[Reference 3](http://www.robots.ox.ac.uk/~az/lectures/ml/lect3.pdf)
"""

# ╔═╡ 9cfc2328-0e4f-4a86-944c-31c0c133fe9a
md"""
The SVM classification problem (directly taken from sklearn website):\

Given training vectors $x_i\in \mathbb{R}^P$, $i=1,...,n$, and a vector $y\in\{1,-1\}^n$\
The goal is to find $w\in \mathbb{R}^P$ and $b\in \mathbb{R}$ such that prediction given by $\text{sign}(w^T\phi(x)+b)$ is accurate for most samples\

Solve the following __primal problem__:\
$$\begin{align}
\min_{w,b,\zeta}\frac{1}{2}w^Tw+C\sum^n_{i=1}\zeta_i\\
\text{subject to }y_i(w^T\phi(x_i)+b)\geq1-\zeta_i,\\
\zeta_i\geq0,i=1,...,n
\end{align}$$\

Intuition is we are maximizing the margin (by minimizing $\|w\|^2=w^Tw$), while penalize when a sample is misclassified\

The __dual problem__ to primal is:\
$$\begin{align}
\min_\alpha\frac{1}{2}\alpha^TQ\alpha-e^T\alpha\\
\text{subject to }y^T\alpha=0,\\
0\leq\alpha_i\leq C,i=1,...,n
\end{align}$$\

 $e$ is vector of all ones, $Q$ is nxn positive semidefinite matrix, $Q_{ij}\equiv y_iy_jK(x_i,x_j)$, and $K(x_i,x_j)=\phi(x_i)^T\phi(x_j)$ is the kernel
$\alpha_i$ are dual coefficients, upper-bounded by $C$\
Highlighting the fact that training vectors are implicitly mapped to a higher (or infinite) dimensional space by function $\phi$\

After optimization problem solved, the __output__ for a given sample $x$ is:
$$\sum_{i\in SV}y_i\alpha_iK(x_i,x)+b$$\

The problem solved by `liblinear` for `LinearSVC` is a equivant form of primal problem:\
$$\min_{w,b}\frac{1}{2}w^Tw+C\sum_{i=1}\max(0,y_i(w^T\phi(x_i)+b))$$\
which does not involve inner products between samples, and therefore cannot apply kernel tricks\

$$C=\frac{1}{\text{alpha}}$$
"""

# ╔═╡ 30975066-eb17-4682-bda1-fd77f9d44734
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 9088c5d0-f07c-404f-8e08-2c8585f86a2f
import Random

# ╔═╡ fc9d6f80-70c1-4518-b9db-29c6e3e2d280
import Statistics

# ╔═╡ 6a0019bf-ab70-4095-ac52-7606e401d9a5
import LinearAlgebra

# ╔═╡ b8447232-6b1d-4b3a-8d11-a9a2da1cf31b
JuTools = ingredients("../tools.jl").JuTools

# ╔═╡ e4380307-4abe-4a20-9127-e608d1a5fa6e
md"""
### Linear SVM Implementation (Hinge Loss)

[Reference](https://stackoverflow.com/questions/48804198/soft-margin-in-linear-support-vector-machine-using-python)\

Cost function is:\
$$J=\frac{1}{2}w^Tw+\frac{C}{N}\sum^N_{i=1}\max\Big(0,1-y_i(w^T\phi(x_i)+b)\Big)$$

Gradient function for $w$ is:\
$$\frac{\partial J}{\partial w}=w + \frac{C}{N}\sum^N_{i=1}\begin{cases}
0 & y_i(w^T\phi(x_i)+b) \geq 1\\
-y_i\phi(x_i) & \text{otherwise}
\end{cases}$$

Gradient function for $b$ is:\
$$\frac{\partial J}{\partial b}=\frac{C}{N}\sum^N_{i=1}\begin{cases}
0 & y_i(w^T\phi(x_i)+b) \geq 1\\
-y_i & \text{otherwise}
\end{cases}$$
"""

# ╔═╡ 6841f5fe-45ad-4f71-84d5-368d196a40ab
# define a struct to store weights
# this should be returned by a training function
# alpha should be treated as constant
mutable struct WeightsLinearSVM
    C::AbstractFloat
    w::Array{T} where T<:AbstractFloat
    b::AbstractFloat
end

# ╔═╡ a093820c-20a0-4af1-98b3-ea92ad37bc77
# define cost function for linear SVM
# assum Y_data is {-1, 1}
function cost(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, weights::WeightsLinearSVM)::AbstractFloat
    @assert ndims(Y_data) == ndims(weights.w) == 1
    @assert size(X_data) == (size(Y_data)[1], size(weights.w)[1])
    loss_w = 0.5 * (weights.w' * weights.w)
    loss_inner = 1.0 .- Y_data .* vec(X_data * weights.w .+ weights.b)
    loss_inner .= map(m->max(0.0,m), loss_inner)
    loss = loss_w + weights.C * sum(loss_inner) / size(X_data)[1]
    return loss
end

# ╔═╡ 8c3fa961-d88d-4e28-98e9-18b8804c568e
X_data, Y_data = JuTools.data_generate_linear_2d()

# ╔═╡ 6341ae48-9cda-4019-bfe6-443caf8788aa
Y_data .= Y_data .* 2.0 .- 1.0 # convert from {0,1} to {-1,1}

# ╔═╡ 500703d2-e315-44c7-8d6c-c36f40b815df
X_train, X_test, Y_train, Y_test = JuTools.split_data(X_data, Y_data)

# ╔═╡ 7531e45c-95fd-4664-8c37-3364aae7020b
@show size(X_train)

# ╔═╡ 747fabf5-fd7b-40cc-b5da-b15f6925cff0
@show size(X_test)

# ╔═╡ 9504574b-2281-4f8a-83d5-3eea38acec53
@show size(Y_train)

# ╔═╡ c496efb2-5a50-4131-bd65-e9025a3273a4
@show size(Y_test)

# ╔═╡ 46e9b187-d445-4efe-b4f2-ad824e5ec832
weight_test = WeightsLinearSVM(1.0, Random.randn(size(X_data)[2]), Random.randn())

# ╔═╡ 8adbc614-9440-42d3-9672-f3b816e114a4
# define the learning function (gradient descent)
function learn!(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, weights::WeightsLinearSVM, alpha::AbstractFloat)
    @assert ndims(Y_data) == ndims(weights.w) == 1
    @assert size(X_data) == (size(Y_data)[1], size(weights.w)[1])
    # compute deciding feature
    decide = (Y_data .* (X_data * weights.w .+ weights.b)) .< 1 # (? < 1) will be 1, otherwise 0
    # update w
    gradient_w = weights.w .+ (weights.C / size(X_data)[1]) .* vec(-(Y_data .* decide)' * X_data)
    gradient_w .= gradient_w .* alpha
    weights.w .= weights.w .- gradient_w
    # update b
    gradient_b = (weights.C / size(X_data)[1]) * sum(-(Y_data .* decide))
    gradient_b *= alpha
    weights.b = weights.b - gradient_b
    return nothing
end

# ╔═╡ 633040d0-651b-427f-a5a8-83124bbea43b
# define prediction function
function predict_proba(X_predict::Array{T} where T<:Number, weights::WeightsLinearSVM)::Array
    @assert ndims(X_predict) == 2
    @assert size(X_predict)[2] == size(weights.w)[1]
    prediction = vec(X_predict * weights.w .+ weights.b)
    return prediction
end

# ╔═╡ ea07e163-d2a0-4447-9239-64d1aee8b3a0
# output prediction is in {-1, 1}
function predict(X_predict::Array{T} where T<:Number, weights::WeightsLinearSVM)::Array
    @assert ndims(X_predict) == 2
    @assert size(X_predict)[2] == size(weights.w)[1]
    prediction = vec(X_predict * weights.w .+ weights.b)
    prediction .= map(m -> m >= 0 ? 1.0 : -1.0, prediction)
    return prediction
end

# ╔═╡ 77bd6c88-335f-46a9-a46b-3d3596d54bf9
md"""
### SVM (with various kernels)

[Reference](https://en.wikipedia.org/wiki/Quadratic_programming)\

For solving SVM minimization problem:\
$$\begin{align}
\min_\alpha\frac{1}{2}\alpha^TQ\alpha-e^T\alpha\\
\text{subject to }y^T\alpha=0,\\
0\leq\alpha_i\leq C,i=1,...,n
\end{align}$$\

where $e$ is vector of all ones, $Q$ is nxn positive semidefinite matrix, $Q_{ij}\equiv y_iy_jK(x_i,x_j)$, and $K(x_i,x_j)=\phi(x_i)^T\phi(x_j)$ is the kernel
$\alpha_i$ are dual coefficients, upper-bounded by $C$\

We'll be using a Quatratic Programming technique: [Sequential Minimal Optimization](https://en.wikipedia.org/wiki/Sequential_minimal_optimization)(SMO)\

My implementation is inspired from [this blog](https://jonchar.net/notebooks/SVM/), whose code is originally from [this paper](https://www.researchgate.net/publication/234786663_Fast_Training_of_Support_Vector_Machines_Using_Sequential_Minimal_Optimization)\

The problem can also be written into (objective function):\
$$\begin{align}
\max_\alpha e^T\alpha-\frac{1}{2}\alpha^TQ\alpha\\
\text{subject to }y^T\alpha=0,\\
0\leq\alpha_i\leq C,i=1,...,n
\end{align}$$\

The kernel functions that I'm going to implement are (using sklearn names):
* Linear Kernel: $\langle x_i,x_j \rangle$
* Polynomial Kernel: $(\gamma\langle x_i,x_j \rangle + r)^d$ ($d$ is degree, and $r$ is coeficient, $\gamma$ is a parameter)
* Rbf Kernel: $\exp(-\gamma\|x_i-x_j\|^2)$, where $\gamma$ is a parameter
* Sigmoid Kernel: $\tanh(\gamma\langle x_i,x_j \rangle + r)$, where $r$ and $\gamma$ are parameters
"""

# ╔═╡ 111304b2-bfbb-4875-838a-93c153d452be
# define a struct to store information
mutable struct WeightsSVM
    C::AbstractFloat                # constraint
    b::AbstractFloat                # threshold
    gamma::AbstractFloat            # parameter used for polynomial, rbf, and sigmoid kernels
    r::AbstractFloat                # parameter used for polynomial, and sigmoid kernels
    d::AbstractFloat                # parameter used for polynomial kernel
    tol_alpha::AbstractFloat        # tolerance for alpha
    tol_error::AbstractFloat        # tolerance for error
    alpha::Array{T} where T<:Number # alpha array
    error::Array{T} where T<:Number # array for error cache
    kernel::String                  # kernel function name
end

# ╔═╡ 3f6728b3-fb75-4783-ab6e-96ae6301f7cb
# linear kernel
function kernel_linear(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = X1 * X2'
    return result
end

# ╔═╡ e3b0cd51-75bc-41aa-96eb-251b376a183d
# polynomial kernel
function kernel_polynomial(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number;
        d::AbstractFloat=1.0, r::AbstractFloat=0.0, gamma::AbstractFloat=1.0)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = (gamma .* (X1 * X2') .+ r) .^ d
    return result
end

# ╔═╡ 4a0ee69a-6151-423b-9227-65fa4e9f2b5e
# rbf kernel
function kernel_rbf(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number; gamma::AbstractFloat=1.0)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = (sum(X1 .^ 2, dims=2) * ones(size(X2)[1])') .+ (ones(size(X1)[1]) * sum(X2 .^ 2, dims=2)') .- 2.0 .* (X1 * X2')
    result .= broadcast(m->max(0.0, m), result) # ignore very small negative outputs, due to precision
    result .= sqrt.(result)
    result .= (-gamma) .* result
    result .= exp.(result)
    return result
end

# ╔═╡ 85939fbf-de18-46d8-9708-dec8a85f49dd
# sigmoid kernel
function kernel_sigmoid(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number; gamma::AbstractFloat=1.0, r::AbstractFloat=0.0)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = gamma .* (X1 * X2') .+ r
    result .= tanh.(gamma)
    return result
end

# ╔═╡ 99d0a23e-e06b-4c13-89a7-62aeebf7c559
# because we call it cost function, we will use the original formula
function cost(X_data::Array{T} where T<:Number, weights::WeightsSVM)::AbstractFloat
    @assert ndims(X_data) == ndims(weights.alpha) + 1 == 2
    @assert size(X_data)[1] == size(weights.alpha)[1]
    result = nothing
    if weights.kernel == "linear"
        result = kernel_linear(X_data, X_data)
    elseif weights.kernel == "polynomial"
        result = kernel_polynomial(X_data, X_data, d=weights.d, r=weights.r, gamma=weights.gamma)
    elseif weights.kernel == "rbf"
        result = kernel_rbf(X_data, X_data, gamma=weights.gamma)
    elseif weights.kernel == "sigmoid"
        result = kernel_sigmoid(X_data, X_data, gamma=weights.gamma, r=weights.r)
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    result = 0.5 * (weights.alpha' * result * weights.alpha) - sum(weights.alpha)
    return result
end

# ╔═╡ 50e8d679-d4fc-4ed8-9868-b59cc4151d2c
cost(X_data, Y_data, weight_test)

# ╔═╡ ece26477-6e53-4536-a6a8-61947cb2815e
# define learning each step function
# update weights in place, and return num steps
function learn_step!(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number,
        weights::WeightsSVM, id1::Integer, id2::Integer)::Integer
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert size(X_data)[1] == size(weights.alpha)[1]
    @assert size(weights.error) == size(weights.alpha)
    @assert id1 >= 1
    @assert id2 >= 1
    # if choosing same alpha, skip
    if id1 == id2
        return 0
    end
    # prepare data
    alpha1 = weights.alpha[id1]
    alpha2 = weights.alpha[id2]
    Y1 = Y_data[id1]
    Y2 = Y_data[id2]
    error1 = weights.error[id1]
    error2 = weights.error[id2]
    # compute L & H
    L = nothing
    H = nothing
    if Y1 != Y2
        L = max(0.0, alpha2 - alpha1)
        H = min(weights.C, weights.C + alpha2 - alpha1)
    else
        L = max(0.0, alpha1 + alpha2 - weights.C)
        H = min(weights.C, alpha1 + alpha2)
    end
    if L == H
        return 0
    end
    # compute kernel results and 2nd derivative eta
    k11 = nothing
    k12 = nothing
    k22 = nothing
    n_features = size(X_data)[2]
    X_id1 = reshape(X_data[id1, :], (1, n_features))
    X_id2 = reshape(X_data[id2, :], (1, n_features))
    if weights.kernel == "linear"
        k11 = kernel_linear(X_id1, X_id1)[1]
        k12 = kernel_linear(X_id1, X_id2)[1]
        k22 = kernel_linear(X_id2, X_id2)[1]
    elseif weights.kernel == "polynomial"
        k11 = kernel_polynomial(X_id1, X_id1, d=weights.d, r=weights.r, gamma=weights.gamma)[1]
        k12 = kernel_polynomial(X_id1, X_id2, d=weights.d, r=weights.r, gamma=weights.gamma)[1]
        k22 = kernel_polynomial(X_id2, X_id2, d=weights.d, r=weights.r, gamma=weights.gamma)[1]
    elseif weights.kernel == "rbf"
        k11 = kernel_rbf(X_id1, X_id1, gamma=weights.gamma)[1]
        k12 = kernel_rbf(X_id1, X_id2, gamma=weights.gamma)[1]
        k22 = kernel_rbf(X_id2, X_id2, gamma=weights.gamma)[1]
    elseif weights.kernel == "sigmoid"
        k11 = kernel_sigmoid(X_id1, X_id1, gamma=weights.gamma, r=weights.r)[1]
        k12 = kernel_sigmoid(X_id1, X_id2, gamma=weights.gamma, r=weights.r)[1]
        k22 = kernel_sigmoid(X_id2, X_id2, gamma=weights.gamma, r=weights.r)[1]
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    eta = 2 * k12 - k11 - k22
    # compute new alpha2 (a2)
    a2 = nothing
    if eta < 0.0
        a2 = alpha2 - Y2 * (error1 - error2) / eta
        a2 = min(a2, H)
        a2 = max(a2, L)
    else
        weights.alpha[id2] = L
        Lobj = -cost(X_data, weights)
        weights.alpha[id2] = H
        Hobj = -cost(X_data, weights)
        weights.alpha[id2] = alpha2
        if Lobj > (Hobj + weights.tol_alpha)
            a2 = L
        elseif Lobj < (Hobj - weights.tol_alpha)
            a2 = H
        else
            a2 = alpha2
        end
    end
    # push to 0 or C
    if a2 < 1e-8
        a2 = 0.0
    elseif a2 > (weights.C - 1e-8)
        a2 = weights.C
    end
    # skip if cannot be optimized
    if abs(a2 - alpha2) < weights.tol_alpha * (a2 + alpha2 + weights.tol_alpha)
        return 0
    end
    # compute new alpha1 (a1)
    a1 = alpha1 + (Y1 * Y2) * (alpha2 - a2)
    if a1 < 0.0
        a2 += (Y1 * Y2) * a1
        a1 = 0.0
    elseif a1 > weights.C
        a2 += (Y1 * Y2) * (a1 - weights.C)
        a1 = weights.C
    end
    # update threshold
    b1 = error1 + Y1 * (a1 - alpha1) * k11 + Y2 * (a2 - alpha2) * k12 + weights.b
    b2 = error2 + Y1 * (a1 - alpha1) * k12 + Y2 * (a2 - alpha2) * k22 + weights.b
    b_new = nothing
    if 0 < a1 < weights.C
        b_new = b1
    elseif 0 < a2 < weights.C
        b_new = b2
    else
        b_new = (b1 + b2) * 0.5
    end
    # update error cache
    non_optimized_ids = [i for i in 1:size(X_data)[1] if (i != id1 && i != id2 && (0 < weights.alpha[i] < weights.C))]
    kerr1 = nothing
    kerr2 = nothing
    if weights.kernel == "linear"
        kerr1 = vec(kernel_linear(X_id1, X_data[non_optimized_ids, :]))
        kerr2 = vec(kernel_linear(X_id2, X_data[non_optimized_ids, :]))
    elseif weights.kernel == "polynomial"
        kerr1 = vec(kernel_polynomial(X_id1, X_data[non_optimized_ids, :], d=weights.d, r=weights.r, gamma=weights.gamma))
        kerr2 = vec(kernel_polynomial(X_id2, X_data[non_optimized_ids, :], d=weights.d, r=weights.r, gamma=weights.gamma))
    elseif weights.kernel == "rbf"
        kerr1 = vec(kernel_rbf(X_id1, X_data[non_optimized_ids, :], gamma=weights.gamma))
        kerr2 = vec(kernel_rbf(X_id2, X_data[non_optimized_ids, :], gamma=weights.gamma))
    elseif weights.kernel == "sigmoid"
        kerr1 = vec(kernel_sigmoid(X_id1, X_data[non_optimized_ids, :], gamma=weights.gamma, r=weights.r))
        kerr2 = vec(kernel_sigmoid(X_id2, X_data[non_optimized_ids, :], gamma=weights.gamma, r=weights.r))
    end
    weights.error[non_optimized_ids] .= weights.error[non_optimized_ids] .+
        ((Y1*(a1-alpha1)) .* kerr1) .+ ((Y2*(a2-alpha2)) .* kerr2) .+ (b_new - weights.b)
    weights.error[id1] = 0.0
    weights.error[id2] = 0.0
    # update alpha and b
    weights.b = b_new
    weights.alpha[id1] = a1
    weights.alpha[id2] = a2
    return 1
end

# ╔═╡ 59c861e3-8c4d-450d-b727-c59cdf17a84f
# now define the learning function
function learn!(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number,
        weights::WeightsSVM, id::Integer; verbose::Bool=false)::Integer
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert size(X_data)[1] == size(weights.alpha)[1]
    @assert size(weights.error) == size(weights.alpha)
    @assert id >= 1
    Y = Y_data[id]
    alpha = weights.alpha[id]
    error = weights.error[id]
    r = error * Y
    if ((r < -weights.tol_error) && (alpha < weights.C)) || ((r > weights.tol_error) && (alpha > 0))
        alpha_target = [i for (i, m) in enumerate(weights.alpha) if (0.0 < m < weights.C)]
        # try argmax E1 - E2
        new_id = -1
        tmax = 0
        if verbose
            println("Trying argmax(abs(E1 - E2))")
        end
        for i in alpha_target
            tmp = abs(error - weights.error[i])
            if(tmp > tmax)
                tmax = tmp
                new_id = i
            end
        end
        if new_id >= 1
            step = learn_step!(X_data, Y_data, weights, id, new_id)
            if step > 0
                return step
            end
        end
        # loop non-bound alphas, randomly
        if verbose
            println("Trying random non-bound alphas")
        end
        for new_id in alpha_target[Random.randperm(length(alpha_target))]
            step = learn_step!(X_data, Y_data, weights, id, new_id)
            if step > 0
                return step
            end
        end
        # else loop all alphas, randomly
        if verbose
            println("Trying random remaining alphas")
        end
        for new_id in Random.randperm(length(weights.alpha))
            if new_id in alpha_target
                continue # skip the alpha ids that already looked at
            end
            step = learn_step!(X_data, Y_data, weights, id, new_id)
            if step > 0
                return step
            end
        end
    end
    return 0
end

# ╔═╡ 882b95f5-6261-42d4-91b2-16a7ce9c1aa4
# implement predict functions
function predict_proba(X_predict::Array{T} where T<:Number, X_data::Array{T} where T<:Number,
        Y_data::Array{T} where T<:Number, weights::WeightsSVM)::Array
    @assert ndims(X_predict) == ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_predict)[2] == size(X_data)[2]
    result = nothing
    if weights.kernel == "linear"
        result = kernel_linear(X_data, X_predict)
    elseif weights.kernel == "polynomial"
        result = kernel_polynomial(X_data, X_predict, d=weights.d, r=weights.r, gamma=weights.gamma)
    elseif weights.kernel == "rbf"
        result = kernel_rbf(X_data, X_predict, gamma=weights.gamma)
    elseif weights.kernel == "sigmoid"
        result = kernel_sigmoid(X_data, X_predict, gamma=weights.gamma, r=weights.r)
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    prediction = vec((weights.alpha .* Y_data)' * result)
    return prediction
end

# ╔═╡ dde248cd-6fac-4c59-877d-58e5b8f06685
function predict(X_predict::Array{T} where T<:Number, X_data::Array{T} where T<:Number,
        Y_data::Array{T} where T<:Number, weights::WeightsSVM)::Array
    @assert ndims(X_predict) == ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_predict)[2] == size(X_data)[2]
    result = nothing
    if weights.kernel == "linear"
        result = kernel_linear(X_data, X_predict)
    elseif weights.kernel == "polynomial"
        result = kernel_polynomial(X_data, X_predict, d=weights.d, r=weights.r, gamma=weights.gamma)
    elseif weights.kernel == "rbf"
        result = kernel_rbf(X_data, X_predict, gamma=weights.gamma)
    elseif weights.kernel == "sigmoid"
        result = kernel_sigmoid(X_data, X_predict, gamma=weights.gamma, r=weights.r)
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    prediction = vec((weights.alpha .* Y_data)' * result)
    prediction .= map(m -> m >= 0 ? 1.0 : -1.0, prediction)
    return prediction
end

# ╔═╡ aafac217-757c-442a-88b1-17e8895f703a
# training function for linear SVM
# assume Y_data is in {-1, 1}
# this function is similar to the training function for Logistic Regression (Both are gradient descent)
function train_linear(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, C::AbstractFloat;
        learning_rate::AbstractFloat=0.1, max_iter::Integer=1000, n_iter_no_change::Integer=5, tol::AbstractFloat=0.001,
        verbose::Bool=false, shuffle::Bool=true, early_stop::Bool=true)::WeightsLinearSVM
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert max_iter >= 0
    @assert n_iter_no_change >= 0
    @assert tol >= 0
    X_data = Float64.(X_data)
    Y_data = Float64.(Y_data)
    if shuffle
        JuTools.shuffle_data!(X_data, Y_data)
    end
    # is it better to use zero weights than normal weights ?
    weights = WeightsLinearSVM(C, Random.randn(size(X_data)[2]), Random.randn())
    best_cost = nothing
    n_cost_no_change = n_iter_no_change
    for i in 1:max_iter
        if n_cost_no_change <= 0 && early_stop
            break
        end
        learn!(X_data, Y_data, weights, learning_rate)
        new_cost = cost(X_data, Y_data, weights)
        if verbose
            acc = JuTools.compute_accuracy(predict(X_data, weights), Y_data)
            println("Iter: $i")
            println("Cost = $new_cost")
            println("Accuracy = $acc")
            println()
        end
        if early_stop
            if best_cost === nothing || isnan(best_cost)
                best_cost = new_cost
            else
                if new_cost > best_cost - tol
                    n_cost_no_change -= 1
                else
                    best_cost = min(new_cost, best_cost)
                    n_cost_no_change = n_iter_no_change
                end
            end
        end
    end
    return weights
end

# ╔═╡ ec256a57-ea44-4f73-af5b-f5e09ba2fec1
weights = train_linear(X_train, Y_train, 1.0, learning_rate=0.005, max_iter=20, tol=0.001, verbose=true)

# ╔═╡ fb710b22-79b0-477c-b2a6-f85fe8b4271a
JuTools.compute_accuracy(predict(X_test, weights), Y_test)

# ╔═╡ 47846484-1824-41b3-b26d-e4e3953fd079
# finally implement the training function
function train(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, C::AbstractFloat;
        tol_alpha::AbstractFloat=0.001, tol_error::AbstractFloat=0.001, kernel::String="rbf", gamma::String="scale",
        degree::AbstractFloat=1.0, coef::AbstractFloat=0.0, verbose::Bool=false)::WeightsSVM
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    X_data = Float64.(X_data)
    Y_data = Float64.(Y_data)
    # gamma is computed the same way sklearn does
    gamma_num = nothing
    if gamma == "scale"
        gamma_num = 1.0 / (size(X_data)[2] * Statistics.var(X_data))
    elseif gamma == "auto"
        gamma_num = 1.0 / size(X_data)[2]
    else
        throw(ArgumentError("Error: gamma $gamma is not recognized, possible values are 'scale' and 'auto'"))
    end
    weights = WeightsSVM(C, 0.0, gamma_num, coef, degree, tol_alpha, tol_error, Float64.(zeros(size(X_data)[1])), -copy(Y_data), kernel)
    num_changed = 0
    examine_all = true
    total_steps = 0
    while (num_changed > 0) || examine_all
        num_changed = 0
        if examine_all
            if verbose
                println("Scanning all training data")
            end
            for i in 1:size(X_data)[1]
                step = learn!(X_data, Y_data, weights, i, verbose=verbose)
                num_changed += step
                if step > 0 && verbose
                    obj = -cost(X_data, weights)
                    println("1 step further, objective = $obj")
                end
            end
        else
            if verbose
                println("Scanning data whose alpha is not at limit")
            end
            alpha_target = [i for (i, m) in enumerate(weights.alpha) if (m != 0.0 && m != weights.C)]
            for i in alpha_target
                step = learn!(X_data, Y_data, weights, i, verbose=verbose)
                num_changed += step
                if step > 0 && verbose
                    obj = -cost(X_data, weights)
                    println("1 step further, objective = $obj")
                end
            end
        end
        if examine_all
            examine_all = false
        elseif num_changed <= 0
            examine_all = true
        end
        total_steps += num_changed
    end
    if verbose
        println("Training Complete\nTotal steps: $total_steps")
    end
    return weights
end

# ╔═╡ fdd56685-18d1-47b5-9b9c-27090422f250
weights_svm = train(X_train, Y_train, 10.0, kernel="rbf", verbose=false, gamma="auto", coef=0.0, degree=2.0)

# ╔═╡ 1b11d952-a35f-4778-8a7f-779eb72f0436
JuTools.compute_accuracy(predict(X_test, X_train, Y_train, weights_svm), Y_test)

# ╔═╡ Cell order:
# ╟─6a41e9d0-c58c-11eb-078e-c36bf22a38d2
# ╟─9cfc2328-0e4f-4a86-944c-31c0c133fe9a
# ╟─30975066-eb17-4682-bda1-fd77f9d44734
# ╠═9088c5d0-f07c-404f-8e08-2c8585f86a2f
# ╠═fc9d6f80-70c1-4518-b9db-29c6e3e2d280
# ╠═6a0019bf-ab70-4095-ac52-7606e401d9a5
# ╠═b8447232-6b1d-4b3a-8d11-a9a2da1cf31b
# ╟─e4380307-4abe-4a20-9127-e608d1a5fa6e
# ╠═6841f5fe-45ad-4f71-84d5-368d196a40ab
# ╠═a093820c-20a0-4af1-98b3-ea92ad37bc77
# ╠═8c3fa961-d88d-4e28-98e9-18b8804c568e
# ╠═6341ae48-9cda-4019-bfe6-443caf8788aa
# ╠═500703d2-e315-44c7-8d6c-c36f40b815df
# ╠═7531e45c-95fd-4664-8c37-3364aae7020b
# ╠═747fabf5-fd7b-40cc-b5da-b15f6925cff0
# ╠═9504574b-2281-4f8a-83d5-3eea38acec53
# ╠═c496efb2-5a50-4131-bd65-e9025a3273a4
# ╠═46e9b187-d445-4efe-b4f2-ad824e5ec832
# ╠═50e8d679-d4fc-4ed8-9868-b59cc4151d2c
# ╠═8adbc614-9440-42d3-9672-f3b816e114a4
# ╠═633040d0-651b-427f-a5a8-83124bbea43b
# ╠═ea07e163-d2a0-4447-9239-64d1aee8b3a0
# ╠═aafac217-757c-442a-88b1-17e8895f703a
# ╠═ec256a57-ea44-4f73-af5b-f5e09ba2fec1
# ╠═fb710b22-79b0-477c-b2a6-f85fe8b4271a
# ╟─77bd6c88-335f-46a9-a46b-3d3596d54bf9
# ╠═111304b2-bfbb-4875-838a-93c153d452be
# ╠═3f6728b3-fb75-4783-ab6e-96ae6301f7cb
# ╠═e3b0cd51-75bc-41aa-96eb-251b376a183d
# ╠═4a0ee69a-6151-423b-9227-65fa4e9f2b5e
# ╠═85939fbf-de18-46d8-9708-dec8a85f49dd
# ╠═99d0a23e-e06b-4c13-89a7-62aeebf7c559
# ╠═ece26477-6e53-4536-a6a8-61947cb2815e
# ╠═59c861e3-8c4d-450d-b727-c59cdf17a84f
# ╠═882b95f5-6261-42d4-91b2-16a7ce9c1aa4
# ╠═dde248cd-6fac-4c59-877d-58e5b8f06685
# ╠═47846484-1824-41b3-b26d-e4e3953fd079
# ╠═fdd56685-18d1-47b5-9b9c-27090422f250
# ╠═1b11d952-a35f-4778-8a7f-779eb72f0436

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
module HMM

#=
Implementation of a Hidden Markov Model (HMM) with a Gaussian emission distribution

Code in part adapted from:

    Zucchini, W., MacDonald, I. L., & Langrock, R. (2017). Hidden Markov models for
    time series: an introduction using R. CRC press.

TODO: add multinomial distribution
TODO: add forecasting function (p.246)
=#

using Random, Plots, Logging, LinearAlgebra, ProgressMeter, Distributions
using Logging;

"""
    simulate_HMM(M::Int64, T::Int64, Γ::Array{Float64}, δ::Array{Float64}, μ::Array{Float64}, σ::Array{Float64})::Tuple{Array{Float64}, Array{Float64}}

Simulate a dataset of a mixture of M gaussians for T timesteps. This dataset will follow a Markov process.

# Arguments
- `M::Int64`: Number of latent states.
- `T::Array{Int64}`: Length of the observed sequence of data.
- `Γ::Array{Float64}`: Transition probability matrix. Array of dimensions (M x M).
- `δ::Array{Float64}`: Initial transition probabilities. Array of dimensions (M x 1).
- `μ::Array{Float64}`: Array of dimensions (M x 1) containing the means for each latent state.
- `σ::Array{Float64}`: Array of dimensions (M x 1) containing the variances for each latent state.

# Returns
An array of dimensions (∑N x M) containing the data points.

# Examples
```julia-repl
julia>
```
"""
function simulate(M::Int64, T::Int64, Γ::Array{Float64}, μ::Array{Float64}, σ::Array{Float64})::Tuple{Array{Float64}, Array{Int64}}
    # Assert dimensions
    @assert M == size(σ)[1] "Number of variances not equal to the number of latent states ..."
    @assert M == size(μ)[1] "Number of means not equal to the number of latent states ..."
    @assert M == size(Γ)[1] == size(Γ)[2] "Transition probability matrix should be of dimensions M x M ..."
    #@assert M == size(δ)[1] "Initial transition probability array should be of dimensions M x 1"
    # TODO: TPM may not contain zero entries
    # Compute the initial distribution
    δ = (Matrix(I, M, M) .- Γ .+ 1) \ ones(M)
    # Populate X and Z
    X = zeros((T))
    Z = zeros(Int64, (T))
    # Draw from latent states
    Z[1] = Multinomial(1, δ) |>
        x -> rand(x, 1) |>
        x -> argmax(x)[1]
    # Draw first X-value
    X[1] = Normal(μ[Z[1]], σ[Z[1]]) |>
            x -> rand(x, 1)[1]
    for t ∈ 2:T
        # Draw from distribution
        Z[t] = Γ[Z[t-1],:] |>
            x -> Multinomial(1, x) |>
            x -> rand(x, 1) |>
            x -> argmax(x)[1]
        # Draw X-value
        X[t,1] = Normal(μ[Z[t,1]], σ[Z[t,1]]) |>
            x -> rand(x, 1)[1]
    end;
    # Return X and Z
    return X, Z
end;

# Utility function to retrieve the normalized probability density of X given parameters on Z
function normalized_pdf(X::Array{Float64}, μ::Array{Float64}, σ::Array{Float64})::Array{Float64}
    #=
    Compute the probabilities of X for each component distribution given its parameters
    :param X: univariate input data. Array of size T x 1
    :param μ: component distribution means. Array of size M x 1
    :param σ: component distribution standard deviations. Array of size M x 1
    :return: Ψ, which is an array of size T x M, where Ψ[:,m] contains the likelihood of X
             given the parameters of component distribution m
    =#
    @debug "Computing likelihood of X ..."
    # Get dimensions
    T = size(X)[1]
    M = size(μ)[1]
    # Results matrix
    Ψ = zeros((T, M));
    # Compute the likelihood of X given parameters on Z
    for m ∈ 1:M
        Ψ[:,m] = Normal(μ[m], σ[m]) |>
            dist -> pdf.(dist, X)
    end;
    # Normalize across rows s.t. each row sums to unity
    Ψ ./= sum(Ψ, dims=2)
    # Return psi
    @debug "Finished computing likelihood of X ..."
    return Ψ
end;

# Forward algorithm
# This function uses logged probabilities to reduce the chance of underflow issues
function forward_algorithm(Γ::Array{Float64}, Ψ::Array{Float64})::Array{Float64}
    #=
    Compute forward step of the forward-backward algorithm
    :param Ψ: Likelihood of input data X given the M component distributions. Array of size T x M. See also
               normalized_pdf()
    :param Γ: transition probabily matrix of size m x m, where m is the number of hidden states
    :param Ψ: emission distribution probabilities. In the case of gaussian emission distributions, this is:

                p(X | Zk) = N(X | μk, Σk) ∀ k

              We assume that we know these and that this matrix Ψ ∈ R^{m x T}

    :param δ: initial distribution matrix of size 1 x m
    :return: matrix Α ∈ R^{T x m} with logged probabilities computed in the forward step
    =#
    @debug "Computing forward probabilities ..."
    # Shapes
    T = size(Ψ)[1]
    M = size(Γ)[1]
    # Compute the initial distribution
    δ = (Matrix(I, M, M) .- Γ .+ 1) \ ones(M)
    # Set up matrix α
    α = zeros((T, M))
    # Forward step using the initial distribution
    κ = δ .* Ψ[1, :]
    ω = sum(κ)[1]
    η = log(ω)
    # Normalize κ by ω (turning it into a valid distribution)
    κ /= ω
    # Add to alpha
    α[1, :] = log.(κ) .+ η
    # Repeat for 2:T
    for t ∈ 2:T
        # Forward step
        κ = (κ' * Γ)' .* Ψ[t,:]
        # Normalize
        ω = sum(κ)[1]
        κ /= ω
        # Add to total
        η += log(ω)
        # Set alpha
        α[t,:] = log.(κ) .+ η
    end;
    # Return alpha
    @debug "Completed forward step ..."
    return α
end;

# Backward algorithm
function backward_algorithm(Γ::Array{Float64}, Ψ::Array{Float64})::Array{Float64}
    #=
    Compute the backward step of the forward-backward algorithm
    :param Ψ: Likelihood of input data X given the M component distributions. Array of size T x M. See also
               normalized_pdf()
    :param Γ: transition probabily matrix of size m x m, where m is the number of hidden states
    :param Ψ: emission distribution probabilities.
    :return: matrix Β ∈ R^{T x m} with logged probabilities computed in the backward step
    =#
    @debug "Computing backward step ...."
    # Shapes
    T = size(Ψ)[1]
    M = size(Γ)[1]
    # Set up matrix α
    β = zeros((T, M))
    # setting ln(beta(T)) = ln(1) = [0] * M
    β[T,:] = zeros((1, M))
    # Compute backward probabilities at T
    κ = repeat([1/M], outer=[M])
    η = log(M)
    # Loop from T-1 to 1
    for t ∈ reverse(1:(T-1))
        κ = Γ * (Ψ[t+1,:] .* κ)
        # Add to beta
        β[t,:] = @. η + log(κ)
        # Sum across kappa and normalize
        ω = sum(κ)[1]
        κ /= ω
        η += log(ω)
    end;
    # Return beta
    @debug "Completed backward step ..."
    return β
end;

# Wrapper to compute forward-backward probabilities using maximum likelihood
function forward_backward_algorithm_EM(Γ::Array{Float64}, Ψ::Array{Float64})::Tuple{Array{Float64}, Array{Float64}}
    #=
    Compute forward-backward algorithm
    =#
    # Compute forward & backward probabilities
    α = forward_algorithm(Γ, Ψ)
    β = backward_algorithm(Γ, Ψ)
    # Return both
    return α, β
end;

# E-step (Baum-Welch algorithm)
function baum_welch_E_step(α::Array{Float64}, β::Array{Float64}, Ψ::Array{Float64}, Γ::Array{Float64})::Tuple{Array{Float64}, Float64}
    #=
    Compute the expectation step of the Baum-Welch algorithm.
    =#
    T = size(Ψ)[1]
    M = size(Γ)[1]
    # Compute log-likelihood
    s1 = α[T, :]
    s2 = s1[argmax(s1)]
    LL = s2 + log(sum(exp.(s1 .- s2)))
    @debug "Log-likelihood of the forward-backward probabilities is $LL ..."
    # For each state, compute next values
    Γ_next = zeros(size(Γ))
    @debug "Size of Γ is $(size(Γ))"
    @debug "E-step: number of states is $M ..."
    for i ∈ 1:M
        for j ∈ 1:M
            inside = @. α[1:(T-1), i] + Ψ[2:T,j] + β[2:T, j] - LL
            @debug "Computed unnormalized transition probabilities going from state $i to state $j ..."
            Γ_next[i,j] = Γ[i,j] * sum(exp.(inside))
        end;
    end;
    # Return gamma next
    return Γ_next, LL
end;

# M-step (Baum-Welch algorithm)
function baum_welch_M_step(α::Array{Float64}, β::Array{Float64}, X::Array{Float64}, LL::Float64, M::Int64)
    #=
    perform the M-step of the Baum-Welch algorithm
    =#
    # Shapes
    T = size(X)[1]
    # Allocate
    μ_next = zeros((M))
    σ_next = zeros((M))
    @debug "M-step: Number of states M is $M ..."
    for i ∈ 1:M
        # Update mean, sd
        p_i = @. exp(α[:,i] + β[:, i] - LL)
        @debug "Computed mixing proportion for state $i ..."
        μ_next[i] = sum(p_i .* X) / sum(p_i)
        @debug "Computed μ_$i ..."
        X_centered = X .- μ_next[i]
        σ_next[i] = sqrt((p_i .* X_centered)' * X_centered / sum(p_i))
        @debug "Computed σ_$i ..."
    end;
    # Return parameters
    return μ_next, σ_next
end;

# Optimization of HMM parameters
# Modified to reduce the chance of underflow errors.
function optim(Γ::Array{Float64}, X::Array{Float64}, μ::Array{Float64}, σ::Array{Float64}; iterations = 150, tol=1e-20)::Tuple{Array{Float64}, Array{Float64}, Array{Float64}, Float64, Float64, Float64}
    #=
    Optimize the HMM parameters to estimate the TPM and the emission distribution parameters
    :param Γ: initial transition probabily matrix of size m x m, where m is the number of hidden states
    :param X: observed data. Sequence of datapoints of length T (T x 1).
    :param μ: intial mean vector of length M x 1.
    :param σ: initial vector containing standard deviations. This array is of length M x 1.
    :return: Tuple containing updated Γ, μ, and σ as well as the log-likelihood, the AIC and the BIC
    =#
    # Get dimensions
    M = size(Γ)[1]
    T = size(X)[1]
    LL = AIC = BIC = NaN
    # For each iteration, estimate parameters
    for n ∈ 1:iterations
        # Compute the likelihood of the data given the parameters
        #  Ψ ∈ R^{T, M}
        Ψ = normalized_pdf(X, μ, σ)
        # Get forward-backward probabilities
        α, β = forward_backward_algorithm_EM(Γ, Ψ)
        # Compute the E-step (Baum-Welch)
        Γ_next, LL = baum_welch_E_step(α, β, Ψ, Γ)
        # Compute the M-step (Baum-Welch)
        μ_next, σ_next = baum_welch_M_step(α, β, X, LL, M)
        # Normalize new TPM
        Γ_next ./= sum(Γ_next, dims=1)'
        # Compute measures of fit
        params = M*M+M-1
        AIC = 2*(LL-params)
        BIC = -2*LL*params*log(T)
        # Compute stopping criterion
        # (sum of difference in parameter change)
        ϵ = (abs.(μ.-μ_next) |> sum, abs.(σ.-σ_next) |> sum, abs.(Γ.-Γ_next) |> sum) |>
            sum
        @debug "Epsilon has value $ϵ on iteration $n with log-likelihood $LL..."
        # Set new parameters to old ones
        Γ = Γ_next
        μ = μ_next
        σ = σ_next
        # If less than tolerance, break
        if ϵ < tol
            @debug "Stopping condition reached after $n iterations. Exiting now ..."
            return Γ, μ, σ, LL, AIC, BIC
        end;
    end;
    @debug "EM algorithm did not converge ..."
    # Return best parameters and fit statistics
    return Γ, μ, σ, LL, AIC, BIC
end;

# Viterbi algorithm
# Used for decoding the most likely set of state sequences across T
function viterbi_algorithm(Γ::Array{Float64}, Ψ::Array{Float64})::Array{Int64}
    #=
    Compute the maximizing latent state sequence
    =#
    # Dimensions
    T = size(Ψ)[1]
    M = size(Γ)[1]
    # Compute the initial distribution
    δ = (Matrix(I, M, M) .- Γ .+ 1) \ ones(M)
    # State sequences
    sequences = zeros(Int64, (T-1, M))
    # Viterbi probabilities
    Ω = zeros((T, M))
    # First time step
    Ω[1,:] = [log(Ψ[1,m]) + log(δ[m]) for m ∈ 1:M]
    @debug "Forward pass (computing set of likely sequences) ...)"
    # Loop through time steps
    for t ∈ 2:T
        if t % 100 == 0
            @debug "Computing sequences at time step $t ..."
        end;
        # For each state, conditioning on the last state, do ...
        for m ∈ 1:M
            # Compute for each state the probability of ending up at state m ...
            prob = zeros((M))
            for j ∈ 1:M
                prob[j] = Ω[t-1, j] + log(Ψ[t, m]) + log(Γ[j, m])
            end;
            # Get argmax
            sequences[t-1, m] = argmax(prob)
            Ω[t, m] = prob[sequences[t-1,m]]
        end;
    end;
    @debug "Finished forward pass ..."
    states = zeros(Int64, (T))
    states[T] = argmax(Ω[T, :])
    # Backward pass
    for t ∈ reverse(1:(T-1))
        states[t] = sequences[t, states[t+1]]
    end;
    # Return the state sequence, log-likelihood
    return states
end;

# Initialize parameters
function initialize_parameters(M::Int64)::Tuple{Array{Float64},Array{Float64}, Array{Float64}}
    #=
    Initialize the TPM and parameters of the distributions
    =#
    # Gamma
    Γ = ones((M,M)) ./= M
    μ = Uniform(0, 1) |>
        x -> sort(rand(x, M))
    @debug "μ initialized as $μ ..."
    # Standard deviations
    σ = ones((M))
    @debug "σ initialized as $σ ..."
    # Return
    return Γ, μ, σ
end;

# HMM function
function fit(X::Array{Float64}, M::Int64; epochs=3, iterations=100, tol=1e-6)
    #=
    Fit a Gaussian Hidden Markov Model to a dataset
    =#
    # Assertions
    # M > 1
    # Get dimensions
    T = size(X)[1]
    # Allocate
    Γ_best = zeros((M,M))
    μ_best = zeros((M))
    σ_best = zeros((M))
    LL_best = NaN
    AIC_best = NaN
    BIC_best = NaN
    # For each epoch
    @showprogress "Fitting HMM ..." for epoch ∈ 1:epochs
        @debug "Fitting HMM on epoch $epoch ..."
        # Initialize parameters
        Γ0, μ0, σ0 = initialize_parameters(M)
        # Fit HMM
        Γ1, μ1, σ1, LL, AIC, BIC = optim(Γ0, X, μ0, σ0; iterations=iterations, tol=tol)
        @debug "Log-likelihood is $LL ..."
        # If loss is better, save
        if epoch == 1
            Γ_best[:] = Γ1[:]
            μ_best[:] = μ1[:]
            σ_best[:] = σ1[:]
            LL_best = LL
            AIC_best = AIC
            BIC_best = BIC
        elseif LL > LL_best
            @debug "Found log-likelihood $LL that improves on best log-likelihood $LL_best (improvement of $(((-1*LL)-(-1*LL_best)) / (-1*LL_best) |> x-> x*100 |> x->round(x, digits=5))%)..."
            Γ_best[:] = Γ1[:]
            μ_best[:] = μ1[:]
            σ_best[:] = σ1[:]
            LL_best = LL
            AIC_best = AIC
            BIC_best = BIC
        else
            continue
        end;
    end;
    # Predicted state sequence
    Ψ = normalized_pdf(X, μ_best, σ_best) # Likelihood given the best parameters
    S = viterbi_algorithm(Γ_best, Ψ)
    # Return best parameters, fit statistics and the predicted sequences
    return (Γ_best, μ_best, σ_best), (LL_best, AIC_best, BIC_best), S
end;

### End module
end;

ext: jl


text:
struct JsFunction{T <: Union{Nothing, NodeValue}} <: JsObjectType
    ref::NodeObject
    this::T
end
value(::Type{JsFunction}, v::NapiValue; this=nothing) = JsFunction(
    NodeObject(v), isnothing(this) ? nothing : node_value(this)
)

(func::ValueTypes)(
    args...;
    recv=nothing, pass_copy=false,
    result=RESULT_VALUE,
    context=current_context()
) = with_result(result) do
    if pass_copy
        args = deepcopy.(args)
    end
    recv = if isnothing(recv)
        if typeof(func) <: JsFunction
            this = getfield(func, :this)
            isnothing(this) ? get_global(context) : this
        else
            get_global(context)
        end
    else
        recv
    end
    argc = length(args)
    p_ca = COPY_ARRAY[]
    if pass_copy
        COPY_ARRAY[] = true
    end
    argv = argc == 0 ? C_NULL : NapiValue.(collect(args))
    COPY_ARRAY[] = p_ca
    @napi_call napi_call_function(
        recv::NapiValue, func::NapiValue, argc::Csize_t, argv::Ptr{NapiValue}
    )::NapiValue
end

function call_function(
    func_ptr::Ptr{Cvoid},
    args_ptr::Ptr{Cvoid},
    argc::Csize_t,
    _recv::Ptr{Cvoid}
)
    func = get_reference(func_ptr)
    args_ptr = Ptr{NapiValue}(args_ptr)
    args = [value(unsafe_load(args_ptr, i)) for i in 1:argc]
    try
        NapiValue(func[](args...))
    catch e
        node_throw(e)
        nothing
    end
end

function napi_value(f::Function; name=nothing)
    func_ptr = pointer(reference(f))
    name = isnothing(name) ? C_NULL : name
    nv = @napi_call create_function(func_ptr::Ptr{Cvoid}, name::Cstring)::NapiValue
    add_finalizer!(nv, dereference, func_ptr)
    nv
end

# Iterators
# napi_value(v::AbstractIterator)
value(::Type{JsIterator}, v::NapiValue) = JsIterator(NodeObject(v))
function Base.iterate(v::NapiValue, state = nothing; result=RESULT_VALUE)
    iterator = isnothing(state) ? v[_JS_ITERATOR_SYMBOL](result=result) : state
    state = iterator.next(result=result)
    state.done && return nothing
    ret = get(state, "value", nothing; result=result)
    ret, iterator
end
Base.iterate(v::ValueTypes, state = nothing) = @with_scope iterate(NapiValue(v), state)


text:
using ModelingToolkit, OrdinaryDiffEq, CardiovascularModels

@variables t
D = Differential(t)

@named vessel1 = ElasticVessel(Ees = 5e6, Vd = 0.0)
@named resistor = Resistor(R = 10e6)
@named vessel2 = ElasticVessel(Ees = 1e7, Vd = 0.0)

eqs_con = [
connect(vessel1.out, resistor.in),
connect(resistor.out, vessel2.in),
vessel1.in.Q ~ 0.,
vessel2.out.Q ~ 0.
]

volume_start = [
    vessel1.V => 1e-4,
    vessel2.V => 2e-4

]

time_span = (0.0, 5.0)

@named connect_sys = ODESystem(eqs_con, t; systems = [vessel1, vessel2, resistor])
problem = ODEProblem(structural_simplify(connect_sys), volume_start, time_span, [])

sol = solve(problem, Tsit5(), dtmax = 0.01, reltol = 1e-6)

# Plots.plot(sol)


text:
module Multivectors

using ComputedFieldTypes
using LinearAlgebra
using StaticArrays

using ..Defs
using ..Forms

export Multivector
@computed struct Multivector{D,M,T}
    elts::SVector{numelts(Val(D), Val(M)),T}
end

@generated function numelts(::Val{D}, ::Val{M}) where {D,M}
    D::Int
    @assert D >= 0
    M::Unsigned
    @assert count_ones(~zero(M)) >= D
    return count_ones(M)
end

# Constructor without explicit type
Multivector{D,M}(elts::SVector{N,T}) where {D,M,N,T} = Multivector{D,M,T}(elts)

# Constructor with added computed type (which must match)
function Multivector{D,M,T,X}(args...) where {D,M,T,X}
    return Multivector{D,M,T}(args...)::Multivector{D,M,T,X}
end

const IteratorTypes = Union{Base.Generator,Iterators.Flatten}
function Multivector{D,M,T}(gen::IteratorTypes) where {D,M,T}
    N = length(Multivector{D,M})
    return Multivector{D,M,T}(SVector{N,T}(elts))
end
function Multivector{D,M}(gen::IteratorTypes) where {D,M}
    @assert IteratorEltype(typeof(gen)) == HasEltype()
    N = length(Multivector{D,M})
    T = eltype(gen)
    return Multivector{D,M,T}(SVector{N,T}(elts))
end

function Multivector{D,M,T}(tup::Tuple) where {D,M,T}
    N = length(Multivector{D,M})
    return Multivector{D,M,T}(SVector{N,T}(tup))
end
function Multivector{D,M}(tup::Tuple) where {D,M}
    N = length(Multivector{D,M})
    return Multivector{D,M}(SVector{N}(tup))
end
function Multivector{D,M}(tup::Tuple{}) where {D,M}
    return error("Cannot create Multivector from emtpy tuple without specifying type")
end

# Conversions
function Multivector{D,M,T}(f::Multivector{D,M}) where {D,M,T}
    N = length(Multivector{D,M})
    return Multivector{D,M,T}(SVector{N,T}(f.elts))
end

################################################################################

# I/O

function Base.show(io::IO, x::Multivector{D,M,T}) where {D,M,T}
    print(io, "$T⟦")
    for n in 1:length(x)
        n > 1 && print(io, ",")
        print(io, x[n])
    end
    print(io, "⟧{D=$D,M=0b$(string(M; base = 2))}")
    return nothing
end

################################################################################

# Iterate over masks

setbit(bits::Unsigned, n) = bits | (one(bits) << n)
getbit(bits::Unsigned, n) = (bits & (1 << n)) != 0
nbits(D) = setbit(zero(UInt64), D)

struct BinaryTerm{I1,I2}
    i1::I1
    i2::I2
end
@generated function binary_algorithm(::Val{D}, ::Val{M1}, ::Val{M2}) where {D,M1,M2}
    D::Int
    @assert D >= 0
    M1::Unsigned
    @assert M1 < nbits(nbits(D))
    M2::Unsigned
    @assert M2 < nbits(nbits(D))
    M = M1 | M2
    r = BinaryTerm[]
    i1 = i2 = i = 0
    for n in 0:(nbits(D) - 1)
        i1 += getbit(M1, n)
        i2 += getbit(M2, n)
        i += getbit(M, n)
        if getbit(M, n)
            i1′ = getbit(M1, n) ? i1 : nothing
            i2′ = getbit(M2, n) ? i2 : nothing
            push!(r, BinaryTerm(i1′, i2′))
        end
    end
    @assert i1 == count_ones(M1)
    @assert i2 == count_ones(M2)
    @assert i == count_ones(M)
    return Tuple(r)
end

################################################################################

# Comparisons

@inline @inbounds function eval_eq_term(term::BinaryTerm, x1, x2)
    x1i1 = term.i1 === nothing ? zero(eltype(x1)) : x1[term.i1]
    x2i2 = term.i2 === nothing ? zero(eltype(x2)) : x2[term.i2]
    return x1i1 == x2i2
end
function Base.:(==)(x1::Multivector{D,M1,T1}, x2::Multivector{D,M2,T2}) where {D,M1,M2,T1,T2}
    algorithm = binary_algorithm(Val(D), Val(M1), Val(M2))
    # return all(map(term -> eval_eq_term(term, x1, x2), algorithm))
    # We want left-to-right evaluation without shortcuts to enable
    # SIMD vectorization
    return mapfoldl(term -> eval_eq_term(term, x1, x2), &, algorithm; init=true)
end
# function Base.:(==)(x1::Multivector{D}, x2::Multivector{D}) where {D}
#     return all(iszero, x1 - x2)
# end
function Base.isequal(x1::Multivector{D,M}, x2::Multivector{D,M}) where {D,M}
    return isequal(x1.elts, x2.elts)
end
Base.hash(x1::Multivector, h::UInt) = hash(hash(x1.elts, h), UInt(0x60671433))

################################################################################

# Multivectors are collections

Base.eltype(::Type{<:Multivector{D,M,T}}) where {D,M,T} = T
Base.firstindex(::Type{<:Multivector}) = 1
Base.firstindex(x::Multivector) = firstindex(typeof(x))
Base.iterate(x::Multivector, state...) = iterate(x.elts, state...)
function Base.lastindex(::Type{<:Multivector{D,M}}) where {D,M}
    return length(Multivector{D,M})
end
Base.lastindex(x::Multivector) = lastindex(typeof(x))
Base.length(::Type{<:Multivector{D,M}}) where {D,M} = count_ones(M)
Base.length(x::Multivector) = length(typeof(x))

Base.getindex(x::Multivector, ind::Integer) = x.elts[ind]
function Base.setindex(x::Multivector{D,M}, val, ind::Integer) where {D,M}
    return Multivector{D,M}(Base.setindex(x.elts, val, ind))
end

function Base.map(f, x::Multivector{D,M}, ys::Multivector{D,M}...) where {D,M}
    return Multivector{D,M}(map(f, x.elts, map(y -> y.elts, ys)...))
end
function Base.reduce(f, x::Multivector{D,M}, ys::Multivector{D,M}...) where {D,M}
    return reduce(f, x.elts, map(y -> y.elts, ys)...)
end

################################################################################

# Multivectors form a vector space

function Base.rand(::Type{<:Multivector{D,M,T}}) where {D,M,T}
    N = length(Multivector{D,M})
    return Multivector{D,M}(rand(SVector{N,T}))
end
function Base.zero(::Type{<:Multivector{D,M,T}}) where {D,M,T}
    N = length(Multivector{D,M})
    return Multivector{D,M}(zero(SVector{N,T}))
end
function Base.zero(::Type{<:Multivector{D,M}}) where {D,M}
    return zero(Multivector{D,M,Float64})
end
Base.zero(::Type{<:Multivector{D}}) where {D} = zero(Multivector{D,0})
Base.zero(x::Multivector) = zero(typeof(x))

function Defs.unit(::Type{<:Multivector{D,M,T}}, inds::SVector{N,<:Integer}) where {D,M,T,N}
    M::Unsigned
    n = lst2lin(Val(D), Val(M), inds)
    @assert n != 0
    return Multivector{D,M}((one(T),))
end
function Defs.unit(::Type{<:Multivector{D,M}}, inds::SVector{N,<:Integer}) where {D,M,N}
    return unit(Multivector{D,M,Float64}, inds)
end
Defs.unit(F::Type{<:Multivector}, inds::Tuple{}) = unit(F, SVector{0,Int}())
Defs.unit(F::Type{<:Multivector}, inds::Tuple) = unit(F, SVector(inds))
Defs.unit(F::Type{<:Multivector}, inds::Integer...) = unit(F, inds)

Base.:+(x::Multivector{D,M}) where {D,M} = Multivector{D,M}(+x.elts)
Base.:-(x::Multivector{D,M}) where {D,M} = Multivector{D,M}(-x.elts)

@inline @inbounds function eval_add_term(term::BinaryTerm, x1, x2)
    term.i1 === nothing && return x2[term.i2]
    term.i2 === nothing && return x1[term.i1]
    return x1[term.i1] + x2[term.i2]
end
function Base.:+(x1::Multivector{D,M1}, x2::Multivector{D,M2}) where {D,M1,M2}
    M = M1 | M2
    M == 0 && return Multivector{D,M,typeof(zero(eltype(x1)) + zero(eltype(x2)))}(())
    algorithm = binary_algorithm(Val(D), Val(M1), Val(M2))
    return Multivector{D,M}(map(term -> eval_add_term(term, x1, x2), algorithm))
end
@inline @inbounds function eval_sub_term(term::BinaryTerm, x1, x2)
    term.i1 === nothing && return -x2[term.i2]
    term.i2 === nothing && return x1[term.i1]
    return x1[term.i1] - x2[term.i2]
end
function Base.:-(x1::Multivector{D,M1}, x2::Multivector{D,M2}) where {D,M1,M2}
    M = M1 | M2
    M == 0 && return Multivector{D,M,typeof(zero(eltype(x1)) - zero(eltype(x2)))}(())
    algorithm = binary_algorithm(Val(D), Val(M1), Val(M2))
    return Multivector{D,M}(map(term -> eval_sub_term(term, x1, x2), algorithm))
end

Base.:*(x::Multivector{D,M}, a) where {D,M} = Multivector{D,M}(x.elts * a)
Base.:/(x::Multivector{D,M}, a) where {D,M} = Multivector{D,M}(x.elts / a)
Base.:*(a, x::Multivector{D,M}) where {D,M} = Multivector{D,M}(a * x.elts)
Base.:\(a, x::Multivector{D,M}) where {D,M} = Multivector{D,M}(a \ x.elts)

################################################################################

# Internal representation

const bits2uint = Forms.bits2uint
const uint2bits = Forms.uint2bits

"""
all bit indices for dimension `D` with stored elements `M`
"""
function dim_bitindices_slow(D::Int, M::Unsigned)
    bitindices = Vector{SVector{D,Bool}}()
    for ibits in CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> 1, D))
        bits = SVector{D,Bool}(Tuple(ibits))
        if getbit(M, bits2uint(bits))
            push!(bitindices, bits)
        end
    end
    return bitindices
end
@generated function dim_bitindices(::Val{D}, ::Val{M}) where {D,M}
    return dim_bitindices_slow(D, M)
end

"""
all linear indices for dimension `D` with stored elements `M`
"""
function dim_linindices_slow(D::Int, M::Unsigned)
    bitindices = dim_bitindices(Val(D), Val(M))
    linindices = zeros(Int, nbits(D))
    for (lin, bits) in enumerate(bitindices)
        ubits = bits2uint(bits)
        linindices[ubits + 1] = lin
    end
    return linindices
end
@generated function dim_linindices(::Val{D}, ::Val{M}) where {D,M}
    return dim_linindices_slow(D, M)
end

"""
convert a linear index to a bit index
"""
function lin2bit(::Val{D}, ::Val{M}, lin::Int) where {D,M}
    D::Int
    M::Unsigned
    return dim_bitindices(Val(D), Val(M))[lin]::SVector{D,Bool}
end

"""
Convert a bit index to a linear index. Return 0 for failure.
"""
function bit2lin(::Val{D}, ::Val{M}, bits::SVector{D,Bool}) where {D,M}
    D::Int
    M::Unsigned
    return dim_linindices(Val(D), Val(M))[bits2uint(bits) + 1]::Int
end

"""
Convert a bit index to an index list
"""
function bit2lst(::Val{D}, bits::SVector{D,Bool}) where {D}
    D::Int
    list = Int[]
    for d in 1:D
        if bits[d]
            push!(list, d)
        end
    end
    return list::Vector{Int}
end

"""
Convert an index list to a bit index
"""
function lst2bit(::Val{D}, list::AbstractVector{<:Integer}) where {D}
    D::Int
    @assert issorted(list)
    bits = zero(SVector{D,Bool})
    for ind in list
        bits = Base.setindex(bits, true, ind)
    end
    return bits::SVector{D,Bool}
end

function lin2lst(::Val{D}, ::Val{M}, lin::Int) where {D,M}
    return bit2lst(Val(D), lin2bit(Val(D), Val(M), lin))
end

function lst2lin(::Val{D}, ::Val{M}, list::AbstractVector{<:Integer}) where {D,M}
    return bit2lin(Val(D), Val(M), lst2bit(Val(D), list))
end

################################################################################

# Multivector algebra

# one

@inline function Base.one(::Type{<:Multivector{D,M,T}}) where {D,M,T}
    return unit(Multivector{D,M,T})
end
@inline function Base.one(::Type{<:Multivector{D,M}}) where {D,M}
    return one(Multivector{D,M,Float64})
end
@inline function Base.one(::Type{<:Multivector{D}}) where {D}
    return one(Multivector{D,UInt64(1)})
end
@inline Base.one(x::Multivector) = one(typeof(x))

# reverse (~)

@generated function reverse_algorithm(::Val{D}, ::Val{M1}) where {D,M1}
    res = Bool[]
    N1 = numelts(Val(D), Val(M1))
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(M1), n1)
        R1 = count(bits1)
        s = isodd((R1 - 1) * R1 ÷ 2)
        push!(res, s)
    end
    return SVector{length(res),Bool}(res)
end
function Base.reverse(x::Multivector{D,M}) where {D,M}
    M == 0 && return x
    algorithm = reverse_algorithm(Val(D), Val(M))
    return Multivector{D,M}(map((s, x) -> s ? -x : x, algorithm, x.elts))
end
@inline Base.:~(x::Multivector) = reverse(x)

# hodge (⋆, \\star)

@generated function hodge_mask(::Val{D}, ::Val{M1}) where {D,M1}
    N1 = numelts(Val(D), Val(M1))
    M = zero(M1)
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bitsr = .~bits1
        M = setbit(M, bits2uint(bitsr))
    end
    return M
end
@generated function hodge_algorithm(::Val{D}, ::Val{M1}) where {D,M1}
    N1 = numelts(Val(D), Val(M1))
    M = hodge_mask(Val(D), Val(M1))
    N = numelts(Val(D), Val(M))
    elts = Any[nothing for n in 1:N]
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bitsr = .~bits1
        _, parity = sort_perm([bit2lst(Val(D), bits1)
                               bit2lst(Val(D), bitsr)])
        s = isodd(parity)
        ind = bit2lin(Val(D), Val(M), bitsr)
        elts[ind] = (s, n1)
    end
    @assert !any(==(nothing), elts)
    return SVector{N,Tuple{Bool,Int}}(elts)
end
@inline @inbounds function eval_hodge_term(term::Tuple{Bool,Int}, x1)
    s, i = term
    return bitsign(s) * x1[i]
end
function Forms.hodge(x1::Multivector{D,M1}) where {D,M1}
    @assert 0 <= D
    M1::Unsigned
    M = hodge_mask(Val(D), Val(M1))
    M == 0 && return zero(Multivector{D,M,eltype(x1)})
    algorithm = hodge_algorithm(Val(D), Val(M1))::SVector
    return Multivector{D,M}(map(term -> eval_hodge_term(term, x1), algorithm))
end

# invhodge (inv(⋆), inv(\\star))

@generated function invhodge_mask(::Val{D}, ::Val{M1}) where {D,M1}
    N1 = numelts(Val(D), Val(M1))
    M = zero(M1)
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bitsr = .~bits1
        M = setbit(M, bits2uint(bitsr))
    end
    return M
end
@generated function invhodge_algorithm(::Val{D}, ::Val{M1}) where {D,M1}
    N1 = numelts(Val(D), Val(M1))
    M = invhodge_mask(Val(D), Val(M1))
    N = numelts(Val(D), Val(M))
    elts = Any[nothing for n in 1:N]
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bitsr = .~bits1
        # The "opposite" parity as `hodge`
        _, parity = sort_perm([bit2lst(Val(D), bitsr)
                               bit2lst(Val(D), bits1)])
        s = isodd(parity)
        ind = bit2lin(Val(D), Val(M), bitsr)
        elts[ind] = (s, n1)
    end
    @assert !any(==(nothing), elts)
    return SVector{N,Tuple{Bool,Int}}(elts)
end
@inline @inbounds function eval_invhodge_term(term::Tuple{Bool,Int}, x1)
    s, i = term
    return bitsign(s) * x1[i]
end
function Forms.invhodge(x1::Multivector{D,M1}) where {D,M1}
    @assert 0 <= D
    M1::Unsigned
    M = invhodge_mask(Val(D), Val(M1))
    M == 0 && return zero(Multivector{D,M,eltype(x1)})
    algorithm = invhodge_algorithm(Val(D), Val(M1))::SVector
    return Multivector{D,M}(map(term -> eval_invhodge_term(term, x1), algorithm))
end

# conj

@inline function Base.conj(x::Multivector{D,M}) where {D,M}
    return Multivector{D,M}(conj.(x.elts))
end

# wedge (∧, \\wedge)

struct WedgeTerm
    sign::Bool
    n1::Int
    n2::Int
end
struct WedgeElt{N}
    terms::NTuple{N,WedgeTerm}
end
@generated function wedge_mask(::Val{D}, ::Val{M1}, ::Val{M2}) where {D,M1,M2}
    D::Int
    M1::Unsigned
    M2::Unsigned
    N1 = numelts(Val(D), Val(M1))
    N2 = numelts(Val(D), Val(M2))
    M = zero(typeof(M1 | M2))
    for n1 in 1:N1, n2 in 1:N2
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bits2 = lin2bit(Val(D), Val(M2), n2)
        if !any(bits1 .& bits2)
            bitsr = bits1 .| bits2
            M = setbit(M, bits2uint(bitsr))
        end
    end
    return M
end
@generated function wedge_algorithm(::Val{D}, ::Val{M1}, ::Val{M2}) where {D,M1,M2}
    D::Int
    M1::Unsigned
    M2::Unsigned
    M = wedge_mask(Val(D), Val(M1), Val(M2))
    N1 = numelts(Val(D), Val(M1))
    N2 = numelts(Val(D), Val(M2))
    N = numelts(Val(D), Val(M))
    elts = [WedgeTerm[] for n in 1:N]
    for n1 in 1:N1, n2 in 1:N2
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bits2 = lin2bit(Val(D), Val(M2), n2)
        if !any(bits1 .& bits2)
            bitsr = bits1 .| bits2
            @assert getbit(M, bits2uint(bitsr))
            _, parity = sort_perm([bit2lst(Val(D), bits1)
                                   bit2lst(Val(D), bits2)])
            s = isodd(parity)
            ind = bit2lin(Val(D), Val(M), bitsr)
            push!(elts[ind], WedgeTerm(s, n1, n2))
        end
    end
    elts = map(elt -> WedgeElt{length(elt)}(Tuple(elt)), elts)
    return Tuple(elts)
end
@inline @inbounds function eval_wedge_elt(elt::WedgeElt{N}, x1, x2) where {N}
    U = typeof(one(eltype(x1)) * one(eltype(x2)))
    N == 0 && return zero(U)
    term = elt.terms[1]
    r = bitsign(term.sign) * x1[term.n1] * x2[term.n2]
    for n in 2:N
        term = elt.terms[n]
        r += bitsign(term.sign) * x1[term.n1] * x2[term.n2]
    end
    return r
end
function Forms.wedge(x1::Multivector{D,M1}, x2::Multivector{D,M2}) where {D,M1,M2}
    @assert 0 <= D
    M1::Unsigned
    M2::Unsigned
    M = wedge_mask(Val(D), Val(M1), Val(M2))
    M == 0 && return zero(Multivector{D,M,typeof(one(eltype(x1)) * one(eltype(x2)))})
    algorithm = wedge_algorithm(Val(D), Val(M1), Val(M2))::Tuple
    return Multivector{D,M}(map(elt -> eval_wedge_elt(elt, x1, x2), algorithm))
end
@inline Forms.wedge(x::Multivector) = x
@inline function Forms.wedge(x1::Multivector, x2::Multivector, x3s::Multivector...)
    return wedge(wedge(x1, x2), x3s...)
end
@inline function Forms.wedge(xs::SVector{0,<:Multivector{D,M,T}}) where {D,M,T}
    return one(Multivector{D,one(M),T})
end
@inline Forms.wedge(xs::SVector{1,<:Multivector}) = xs[1]
@inline function Forms.wedge(xs::SVector{N,<:Multivector}) where {N}
    return ∧(pop(xs)) ∧ last(xs)
end

# vee (∨, \\vee)

@inline Forms.vee(x1::Multivector, x2::Multivector) = inv(⋆)(⋆x1 ∧ ⋆x2)
@inline Forms.vee(x::Multivector) = x
@inline vee′(x1::Multivector) = x1
@inline vee′(x1::Multivector, x2::Multivector) = x1 ∧ ⋆x2
@inline function vee′(x1::Multivector, x2::Multivector, x3s::Multivector...)
    return vee′(x1 ∧ ⋆x2, x3s...)
end
@inline function Forms.vee(x1::Multivector, x2::Multivector, x3s::Multivector...)
    # vee(x1 ∨ x2, x3s...)
    return inv(⋆)(vee′(⋆x1, x2, x3s...))
end
@inline function Forms.vee(xs::SVector{0,<:Multivector{D,M,T}}) where {D,M,T}
    return ⋆one(Multivector{D,one(M),T})
end
@inline Forms.vee(xs::SVector{1,<:Multivector}) = xs[1]
@inline function vee′(xs::SVector{0,<:Multivector{D,M,T}}) where {D,M,T}
    return one(Multivector{D,one(M),T})
end
@inline vee′(xs::SVector{1,<:Multivector}) = ⋆xs[1]
@inline vee′(xs::SVector{N,<:Multivector}) where {N} = vee′(pop(xs)) ∧ ⋆last(xs)
@inline function Forms.vee(xs::SVector{N,<:Multivector}) where {N}
    # vee(pop(xs)) ∨ last(xs)
    return inv(⋆)(vee′(xs))
end

# dot (⋅, \\cdot)

@inline LinearAlgebra.dot(x1::Multivector, x2::Multivector) = x1 ∨ ⋆x2

@inline function Forms.norm2(x::Multivector{D,M}) where {D,M}
    M == 0 && return zero(eltype(x))
    return (x ⋅ x)[1]
end
@inline LinearAlgebra.norm(x::Multivector) = sqrt(norm2(x))

# cross (×, \\times)

@inline LinearAlgebra.cross(x1::Multivector, x2::Multivector) = ⋆(x1 ∧ x2)

# geometric product (*)

struct MulTerm
    sign::Bool
    n1::Int
    n2::Int
end
struct MulElt{N}
    terms::NTuple{N,MulTerm}
end
@generated function mul_mask(::Val{D}, ::Val{M1}, ::Val{M2}) where {D,M1,M2}
    D::Int
    M1::Unsigned
    M2::Unsigned
    N1 = numelts(Val(D), Val(M1))
    N2 = numelts(Val(D), Val(M2))
    M = zero(typeof(M1 | M2))
    for n1 in 1:N1, n2 in 1:N2
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bits2 = lin2bit(Val(D), Val(M2), n2)
        bitsr = (bits1 .| bits2) .& .~(bits1 .& bits2)
        M = setbit(M, bits2uint(bitsr))
    end
    return M
end
@generated function mul_algorithm(::Val{D}, ::Val{M1}, ::Val{M2}) where {D,M1,M2}
    D::Int
    M1::Unsigned
    M2::Unsigned
    M = mul_mask(Val(D), Val(M1), Val(M2))
    N1 = numelts(Val(D), Val(M1))
    N2 = numelts(Val(D), Val(M2))
    N = numelts(Val(D), Val(M))
    elts = [MulTerm[] for n in 1:N]
    for n1 in 1:N1, n2 in 1:N2
        bits1 = lin2bit(Val(D), Val(M1), n1)
        bits2 = lin2bit(Val(D), Val(M2), n2)
        bitsr = (bits1 .| bits2) .& .~(bits1 .& bits2)
        @assert getbit(M, bits2uint(bitsr))
        _, parity = sort_perm([bit2lst(Val(D), bits1)
                               bit2lst(Val(D), bits2)])
        s = isodd(parity)
        ind = bit2lin(Val(D), Val(M), bitsr)
        push!(elts[ind], MulTerm(s, n1, n2))
    end
    elts = map(elt -> MulElt{length(elt)}(Tuple(elt)), elts)
    return Tuple(elts)
end
@inline @inbounds function eval_mul_elt(elt::MulElt{N}, x1, x2) where {N}
    U = typeof(one(eltype(x1)) * one(eltype(x2)))
    N == 0 && return zero(U)
    term = elt.terms[1]
    r = bitsign(term.sign) * x1[term.n1] * x2[term.n2]
    for n in 2:N
        term = elt.terms[n]
        r += bitsign(term.sign) * x1[term.n1] * x2[term.n2]
    end
    return r
end
function Base.:*(x1::Multivector{D,M1}, x2::Multivector{D,M2}) where {D,M1,M2}
    @assert 0 <= D
    M1::Unsigned
    M2::Unsigned
    M = mul_mask(Val(D), Val(M1), Val(M2))
    M == 0 && return zero(Multivector{D,M,typeof(one(eltype(x1)) * one(eltype(x2)))})
    algorithm = mul_algorithm(Val(D), Val(M1), Val(M2))::Tuple
    return Multivector{D,M}(map(elt -> eval_mul_elt(elt, x1, x2), algorithm))
end
@inline Base.:*(x::Multivector) = x
@inline function Base.:*(x1::Multivector, x2::Multivector, x3s::Multivector...)
    return *(x1 * x2, x3s...)
end
@inline function Base.:*(xs::SVector{0,<:Multivector{D,M,T}}) where {D,M,T}
    return one(Multivector{D,one(M),T})
end
@inline Base.:*(xs::SVector{1,<:Multivector}) = xs[1]
@inline Base.:*(xs::SVector{N,<:Multivector}) where {N} = *(pop(xs)) * last(xs)

end


text:
using Distributions
using SpinGlassNetworks

function ramp(t::T, τ::T, α::T, pi::T, pf::T) where T <: Real
    p = (pf + pi) + (pf - pi) * tanh(α * (2.0 * t / τ - 1.0))
    p / 2.0
end

@testset "Coherent Ising Machine simulator for chimera instances." begin
    L = 128

    ig = ising_graph("$(@__DIR__)/instances/chimera_droplets/$(L)power/001.txt")

    scale = 0.3
    noise = Normal(0.1, 0.3)

    x0 = zeros(L)
    sat = 1.0
    time = 1000.
    pi, pf, α = -15.0, 0.0, 2.0
    momentum = 0.4

    pump = [ramp(t, time, α, pi, pf) for t ∈ 0:time]

    opo = OpticalOscillators{Float64}(ig, scale, noise)
    dyn = OPODynamics{Float64}(x0, sat, pump, momentum)

    pump_nmfa = [ramp(t, time, α, 10., 1.) for t ∈ 0:time]

    opo_nmfa = OpticalOscillators{Float64}(ig, scale, noise)
    dyn_nmfa = OPODynamics{Float64}(x0, sat, pump_nmfa, momentum)

    N = 500
    states = Vector{Vector{Int}}(undef, N)
    states_nmfa = copy(states)
    Threads.@threads for i ∈ 1:N
        states[i] = evolve_optical_oscillators(opo, dyn)
        states_nmfa[i] = noisy_mean_field_annealing(opo_nmfa, dyn_nmfa)
    end

    en = minimum(energy(states, ig))
    en_nmfa = minimum(energy(states_nmfa, ig))

    en_exact = -210.933334
    @testset "Energy found is at least negative and within the bounds" begin
        @test  en_exact <= en < 0.0
        @test en_exact <= en_nmfa < 0.0
    end
end


text:
import CLBLAS

const clblas = CLBLAS
clblas.setup()

using OpenCL: cl

device, ctx, queue = clblas.get_next_compute_context()

N = Unsigned(5)
data = Vector{Float32}(5)
data[1] = 1.1
data[2] = 2.22
data[3] = 3.333
data[4] = 4.4444
data[5] = 5.55555

bufX = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=data)
bufAsum = cl.Buffer(Float32, ctx, :w, length(bufX))
scratchBuff = cl.Buffer(Float32, ctx, :rw, length(bufX))

offx = Unsigned(0)
incx = cl.cl_int(1)
ncq = cl.cl_uint(1)
clq = queue.id
newl = cl.cl_uint(0)

event = cl.UserEvent(ctx, retain=true)
ptrEvent = [event.id]

clblas.clblasSasum( N, bufAsum.id, 0, bufX.id, 0, incx, scratchBuff.id,
                                ncq, [clq], newl, C_NULL, ptrEvent);

cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent)

result = Vector{Float32}(length(1))
cl.enqueue_read_buffer(queue, bufAsum, result, Unsigned(0), nothing, true)

expected = data[1:end] |> sum
println(result)
println("------------")
println(expected)

try
    @assert(isapprox(norm(expected - result[1]), zero(Float32)))
    info("success!")
finally
    clblas.teardown()
end


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

"""
    eta_quotient(g, prec=10)

Return the eta-quotient relative to the given collection of integers
`g` as a power series up to precision `prec`.

The input g is an array of pairs [[t_1,r_1], [t_2,r_2], ..., [t_s,r_s]]
where each t_j is a positive integer and each r_j a nonnegative integer.
An error is thrown if the sum of t_j*r_j for j in {1, ..., s} does not
equal 24, as the function only applies to eta-quotients that are cusp
forms. The output is a cusp form of integral weight.

# Arguments
- `g::Array{Array}`: an array of pairs, collection of integers
- `prec::Integer=10`: precision of the output

# Examples still missing
"""
function eta_quotient(g::Array{Array{Int,1},1}, prec::Int=10)

   R, q = PolynomialRing(ZZ, "q")
   poly = R(1)		#product
   sum = 0 		#error handling
   s = length(g)
   for i in 1:s
      t = g[i][1]
      r = g[i][2]
      poly = poly * eta_qexp(r, prec, q^t)	#check prec
      sum += t*r
   end

   if sum != 24
      error("sum of t_j*r_j for all j must equal 24 if it is a cusp form")
   end

   poly = q * poly
   power_series = poly_to_power_series(poly, prec)

   return power_series

end


text:
# ---
# title: 1104. Path In Zigzag Labelled Binary Tree
# id: problem1104
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Tree
# link: <https://leetcode.com/problems/path-in-zigzag-labelled-binary-tree/description/>
# hidden: true
# ---
#
# In an infinite binary tree where every node has two children, the nodes are
# labelled in row order.
#
# In the odd numbered rows (ie., the first, third, fifth,...), the labelling is
# left to right, while in the even numbered rows (second, fourth, sixth,...),
# the labelling is right to left.
#
# ![](https://assets.leetcode.com/uploads/2019/06/24/tree.png)
#
# Given the `label` of a node in this tree, return the labels in the path from
# the root of the tree to the node with that `label`.
#
#
#
# **Example 1:**
#
#
#
#     Input: label = 14
#     Output: [1,3,4,14]
#
#
# **Example 2:**
#
#
#
#     Input: label = 26
#     Output: [1,2,6,10,26]
#
#
#
#
# **Constraints:**
#
#   * `1 <= label <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =========================================================================================
"""
    tess_get_param(
        inst::TessInst,
        name::AbstractString,
        ::Type{T}
    )::Union{T, Nothing} where T<:Integer

Retrieve an integer parameter from the Tesseract engine.  Returns `nothing` if the value
could not be retrieved.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | inst      |         | The instance to read the parameter from.
| R | name      |         | The name of the parameter to read.
| R | ::Type{T} |         | The type to return.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_get_param(instance, "edges_min_nonhole", Int)
12
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_set_param`](@ref)
"""
function tess_get_param(
            inst::TessInst,
            name::AbstractString,
            ::Type{T}
        )::Union{T, Nothing} where T<:Integer

    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    local value = Ref{Cint}(0)
    local result = ccall(
        (:TessBaseAPIGetIntVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Ptr{Cint}),
        inst,
        name,
        value
    )

    if result == 0
        return nothing
    end

    return value[]
end

# =========================================================================================
"""
    tess_get_param(
        inst::TessInst,
        name::AbstractString,
        ::Type{Bool}
    )::Union{Bool, Nothing}

Retrieve a boolean parameter from the Tesseract engine.  Returns `nothing` if the value
could not be retrieved.

__Arguments:__

| T | Name         | Default | Description
|:--| :----------- | :------ | :----------
| R | inst         |         | The instance to read the variable from.
| R | name         |         | The name of the variable to read.
| R | ::Type{Bool} |         | Identifies that you want a boolean value.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_get_param(instance, "edges_debug", Bool)
false
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_set_param`](@ref)
"""
function tess_get_param(
            inst::TessInst,
            name::AbstractString,
            ::Type{Bool}
        )::Union{Bool, Nothing}

    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    local value = Ref{Cint}(0)
    local result = ccall(
        (:TessBaseAPIGetBoolVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Ptr{Cint}),
        inst,
        name,
        value
    )

    if result == 0
        return nothing
    end

    return value[] == 1
end

# =========================================================================================
"""
    tess_get_param(
        inst::TessInst,
        name::AbstractString,
        ::Type{Float64}
    )::Union{Float64, Nothing}

Retrieve a Float parameter from the Tesseract engine.  Returns `nothing` if the value
could not be retrieved.

__Arguments:__

| T | Name            | Default | Description
|:--| :-------------- | :------ | :----------
| R | inst            |         | The instance to read the variable from.
| R | name            |         | The name of the parameter to read.
| R | ::Type{Float64} |         | Identifies that you want a Float64 value.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_get_param(instance, "classify_min_slope", Float64)
0.414214
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_set_param`](@ref)
"""
function tess_get_param(
            inst::TessInst,
            name::AbstractString,
            ::Type{Float64}
        )::Union{Float64, Nothing}

    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    local value = Ref{Cdouble}(0.0)
    local result = ccall(
        (:TessBaseAPIGetDoubleVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Ptr{Cdouble}),
        inst,
        name,
        value
    )

    if result == 0
        return nothing
    end

    return value[]
end

# =========================================================================================
"""
    tess_get_param(
        inst::TessInst,
        name::AbstractString,
        ::Type{String}
    )::Union{String, Nothing}

Retrieve a String parameter from the Tesseract engine.  Returns `nothing` if the value
could not be retrieved.

__Arguments:__

| T | Name            | Default | Description
|:--| :-------------- | :------ | :----------
| R | inst            |         | The instance to read the variable from.
| R | name            |         | The name of the parameter to read.
| R | ::Type{Float64} |         | Identifies that you want a Float64 value.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_get_param(instance, "page_separator", String)
"\\f"
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_set_param`](@ref)
"""
function tess_get_param(
            inst::TessInst,
            name::AbstractString,
            ::Type{String}
        )::Union{String, Nothing}

    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    local result = ccall(
        (:TessBaseAPIGetStringVariable, TESSERACT),
        Cstring,
        (Ptr{Cvoid}, Cstring),
        inst,
        name
    )

    if result == C_NULL
        return nothing
    end

    return unsafe_string(result)
end


text:
using Simulacrum, Simulacrum.Wip
using HyperGraphs, Symbolics
using Catalyst

@variables t X(t) Y(t) Z(t)
che1 = ChemicalHyperEdge([X], [Y], 2)
che2 = ChemicalHyperEdge([X, Y], [Z], 2)
che3 = ChemicalHyperEdge(SpeciesSet([X], [3]), SpeciesSet{Num}(), 1)
chx = ChemicalHyperGraph([che1, che2, che3])

# full dynamics on chemical hyperedges, purely from assumptions and topology
ode_rate_law(che::ChemicalHyperEdge{Num}) = rn_assumptions(che)
ode_rate_law(chx::ChemicalHyperGraph{Num}) = ode_rate_law.(hyperedges(chx))

## propensities ##

# mass action
@test isequal(mass_action(che1), 2*X)
@test isequal(mass_action(che2), 2*X*Y)
@test isequal(mass_action(che3), X)

# power law
@test isequal(power_law(che1), [X])
@test isequal(power_law(che2), [X, Y])
@test isequal(power_law(che3), [X^3])

# combinatoric power law
@test isequal(combinatoric_power_law(che1), [X])
@test isequal(combinatoric_power_law(che2), [X, Y])
@test isequal(combinatoric_power_law(che3), [(X^3)/6])

## rate laws ##

# ODE rate law
@test isequal(ode_rate_law(che1), 2*X)
@test isequal(ode_rate_law(che2), 2*X*Y)
@test isequal(ode_rate_law(che3), (X^3)/6)

# full equations
D = Differential(:t)
ref_eqs = [ Equation(D(X), - 2*X - 2*X*Y - (X^3)/2),
            Equation(D(Y), 2*X - 2*X*Y),
            Equation(D(Z), 2*X*Y)]
@test isequal(make_equations(chx, ode_rate_law), ref_eqs)

# zeroth order reactions
@test isequal(ode_rate_law(ChemicalHyperEdge(Num[], [X], 3)), 3)
@test isequal(ode_rate_law(ChemicalHyperEdge([X], Num[], 3)), 3X)

## benchmarking against ModelingToolkit to make sure the equations are correct ##

chx_eqs = Simulacrum.get_value(make_equations(chx, ode_rate_law))
r1 = Reaction(2, [X], [Y])
r2 = Reaction(2, [X, Y], [Z])
r3 = Reaction(1, [X], nothing, [3], nothing)
@named rsys = ReactionSystem([r1, r2, r3], t, [X, Y, Z], [])
sys = convert(ODESystem, rsys)
mtk_eqs = equations(sys)
chx_eqs_rhss = [eq.rhs for eq in chx_eqs]
mtk_eqs_rhss = [eq.rhs for eq in mtk_eqs]
@test isequal(chx_eqs_rhss, mtk_eqs_rhss) # for some reason the left hand sides aren't equal


text:
type SysState
    X::Array{Float64,1}
    P::Array{Float64,2}
end

function c2d!(M1::Array{Float64,2},
              M2::Array{Float64,2},
              F::Array{Float64,2},
              G::Array{Float64,2},
              Qc::Array{Float64,2},
              dt::Float64;
              expmfnc::Function=expm)
    # more accurate explicit matrix exponential approach --
    # computationally efficient approaches can be done here to lighten the load

    # M1 = [[F';G']';zeros(12,21)]
    # M1 = zeros(21,21)
    M1[1:9,1:9] = F
    M1[1:9,10:21] = G
    M1[10:21,1:21] = 0.0

    Md1 = expmfnc(M1*dt)
    Phi = Md1[1:9,1:9]
    Gamma = Md1[1:9,10:21]

    #M2 = [[-F';(G*Qc*G')']';[zeros(9,9);F]']
    # M2 = zeros(18,18)
    M2[1:9,1:9] = -F
    M2[1:9,10:18] = (G*Qc*G')
    M2[10:18,10:18] = F'
    M2[10:18,1:9] = 0.0

    Md2 = expmfnc(M2*dt)
    Qd = Phi*(Md2[1:9,10:18])
    return Phi, Gamma, Qd
end

# cpp: ReturnTypePhiGammaQd c2d(const Eigen::MatrixXd &F, const Eigen::MatrixXd &G, const Eigen::MatrixXd &Qc, const double &dt) {
function c2d(F::Array{Float64,2}, G::Array{Float64,2}, Qc::Array{Float64,2}, dt::Float64;
            expmfnc::Function=expm)
    # more accurate explicit matrix exponential approach --
    # computationally efficient approaches can be done here to lighten the load

    # M1 = [[F';G']';zeros(12,21)] # easy concat
    M1 = zeros(21,21)
    # Eigen has M1.topLeftCorner(9,9) = F;
    M1[1:9,1:9] = F
    M1[1:9,10:21] = G

    Md1 = expmfnc(M1*dt)
    # Eigen has Phi = Md1.topLeftCorner(9,9);
    Phi = Md1[1:9,1:9]
    Gamma = Md1[1:9,10:21]

    #M2 = [[-F';(G*Qc*G')']';[zeros(9,9);F]'] # easy concat
    M2 = zeros(18,18)
    M2[1:9,1:9] = -F
    M2[1:9,10:18] = (G*Qc*G')
    M2[10:18,10:18] = F'

    Md2 = expmfnc(M2*dt)
    Qd = Phi*(Md2[1:9,10:18])
    return Phi, Gamma, Qd
end

function filterpropagate!(st::SysState, wQb::Quaternion, dt::Float64, Qcont::Array{Float64,2})
    F, G = FG(wQb)
    Phi, Gamma, Qd = c2d(F, G, Qcont, dt)
    # ultimately move M1 and M2 declaration all the way up to above the iteration loop to allow memory reuse

    # estimated error state
    st.X = Phi*st.X # + Gamma*noiseinput, which we can only approximate
    # noise increases uncertainty (error covariance)
    st.P = Phi*st.P*Phi' + Qd
    nothing
end

function filterupdate!(st::SysState, r::Array{Float64,1}, H::Array{Float64,2}, R::Array{Float64,2})
    # K = P H' (A.inv) = P (H'*Ainv) = P (A' \ H)'
    At = (H*st.P*(H') + R)'
    K = st.P * ((At \ H)')
    # inverse way
    #K = st.P*H'*((H*st.P*(H') + R) \ eye(size(R,1))) # inverse
    st.X += K*r
    KH = K*H
    P = (eye(size(KH,1)) - KH)*st.P
    st.P = P
    nothing
end


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

"""
    CrossValidation(k, [shuffle])

Compare estimation solvers using k-fold cross-validation.
Optionally shuffle the data (default to true).

    CrossValidation(partitioner)

Compare estimation solvers using cross-validation by
splitting the data with a `partitioner`. This method
is a generalization of k-fold cross-validation, which
uses a [`UniformPartitioner`](@ref) to split the data.

The result of the comparison stores the errors for each
variable of the problem.

## Examples

Compare `solver₁` and `solver₂` on a `problem` with variable
`:var` using 10 folds. Plot error distribution:

```julia
julia> result = compare([solver₁, solver₂], problem, CrossValidation(10))
julia> plot(result, bins=50)
```
"""
struct CrossValidation{P} <: AbstractErrorEstimator
  partitioner::P
end

CrossValidation(k::Int, shuffle::Bool) =
  CrossValidation(UniformPartitioner(k, shuffle))
CrossValidation(k::Int) = CrossValidation(k, true)
CrossValidation() = CrossValidation(10)

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::CrossValidation)
  # retrieve problem info
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  partitioner = eestimator.partitioner

  # folds for cross-validation
  folds  = subsets(partition(sdata, partitioner))
  nfolds = length(folds)

  solutions = pmap(1:nfolds) do k
    # source and target indices
    sinds = [ind for i in vcat(1:k-1, k+1:nfolds) for ind in folds[i]]
    tinds = folds[k]

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)

    # setup and solve sub-problem
    subproblem = LearningProblem(train, hold, task(problem))
    solve(subproblem, solver)
  end

  result = pmap(ovars) do var
    𝔏 = defaultloss(sdata[1,var])
    losses = map(1:nfolds) do k
      hold = view(sdata, folds[k])
      ŷ = solutions[k][var]
      y = hold[var]
      𝔏(ŷ, y)
    end
    var => mean(losses)
  end

  Dict(result)
end

function Base.error(solver::AbstractEstimationSolver,
                    problem::EstimationProblem,
                    eestimator::CrossValidation)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)
  partitioner = eestimator.partitioner

  # folds for cross-validation
  folds  = subsets(partition(pdata, partitioner))
  nfolds = length(folds)

  result = []
  for (var, V) in variables(problem)
    # mappings from data to domain locations
    varmap = Dict(datloc => domloc for (domloc, datloc) in datamap(problem, var))

    # k-fold validation loop
    losses = pmap(1:nfolds) do k
      # training and holdout set
      train = [ind for i in vcat(1:k-1, k+1:nfolds) for ind in folds[i]]
      hold  = folds[k]

      # discard indices that were filtered out by mapping strategy (e.g. missing values)
      train = filter(in(keys(varmap)), train)
      hold  = filter(in(keys(varmap)), hold)

      if isempty(train)
        missing
      else
        # copy data to their locations in domain
        mapper = CopyMapper([varmap[ind] for ind in train])

        # setup and solve sub-problem
        subproblem = EstimationProblem(view(pdata, train), pdomain,
                                       var, mapper=mapper)
        solution = solve(subproblem, solver)

        # get true holdout values
        y = pdata[hold,var]

        # get solver estimate at holdout locations
        μ, _ = solution[var]
        ŷ = [μ[varmap[loc]] for loc in hold]

        # loss for fold
        mean((ŷ .- y).^2)
      end
    end

    push!(result, var => mean(skipmissing(losses)))
  end

  Dict(result)
end


text:
#=
utilities for benchmark scripts
=#

function setup_benchmark_dataframe()
    perf = DataFrames.DataFrame(
        example = String[],
        model_type = String[],
        inst_set = String[],
        inst_num = Int[],
        inst_data = Tuple[],
        extender = String[],
        real_T = Type{<:Real}[],
        solver = String[],
        solver_options = Tuple[],
        nondefaults = Tuple[],
        script_status = String[],
        n = Int[],
        p = Int[],
        q = Int[],
        nu = Float64[],
        cone_types = Vector{String}[],
        num_cones = Int[],
        max_q = Int[],
        status = String[],
        solve_time = Float64[],
        iters = Int[],
        primal_obj = Float64[],
        dual_obj = Float64[],
        rel_obj_diff = Float64[],
        compl = Float64[],
        x_viol = Float64[],
        y_viol = Float64[],
        z_viol = Float64[],
        time_rescale = Float64[],
        time_initx = Float64[],
        time_inity = Float64[],
        time_loadsys = Float64[],
        time_upsys = Float64[],
        time_upfact = Float64[],
        time_uprhs = Float64[],
        time_getdir = Float64[],
        time_search = Float64[],
        setup_time = Float64[],
        check_time = Float64[],
        total_time = Float64[],
        )
    DataFrames.allowmissing!(perf, 11:DataFrames.ncol(perf))
    DataFrames.allowmissing!(perf, [:extender, :solver_options, :nondefaults])
    return perf
end

function get_extender(inst::Tuple, ::Type{<:ExampleInstanceJuMP})
    if length(inst) > 1
        ext = inst[2]
        !isnothing(ext) && return string(ext)
    end
    return missing
end

get_extender(inst::Tuple, ::Type{<:ExampleInstance}) = missing

function get_nondefaults(inst::Tuple, ::Type{<:ExampleInstanceJuMP})
    if length(inst) > 2
        @assert length(inst) == 3
        return inst[3]
    end
    return missing
end

function get_nondefaults(inst::Tuple, ::Type{<:ExampleInstance})
    if length(inst) > 1
        @assert length(inst) == 2
        return inst[2]
    end
    return missing
end

function write_perf(
    perf::DataFrames.DataFrame,
    results_path::Union{String, Nothing},
    new_perf::NamedTuple,
    )
    push!(perf, new_perf, cols = :subset)
    if !isnothing(results_path)
        CSV.write(results_path, perf[end:end, :], transform =
            (col, val) -> something(val, missing), append = true)
    end
    return
end

function run_instance_set(
    inst_subset::Vector,
    ex_type_T::Type{<:ExampleInstance},
    info_perf::NamedTuple,
    new_default_options::NamedTuple,
    script_verbose::Bool,
    perf::DataFrames.DataFrame,
    results_path::Union{String, Nothing} = nothing,
    )
    for (inst_num, inst) in enumerate(inst_subset)
        extender_name = get_extender(inst, ex_type_T)
        nondefaults = get_nondefaults(inst, ex_type_T)
        test_info = "inst $inst_num: $(inst[1])"
        if !ismissing(extender_name)
            test_info *= " $extender_name"
        end
        if !ismissing(nondefaults)
            test_info *= " $nondefaults"
        end

        @testset "$test_info" begin
            println(test_info, " ...")

            total_time = @elapsed run_perf = run_instance(ex_type_T, inst...,
                default_options = new_default_options, verbose = script_verbose)

            new_perf = (; info_perf..., run_perf..., total_time, inst_num,
                :solver => "Hypatia", :inst_data => inst[1],
                :extender => extender_name, :nondefaults => nondefaults)
            write_perf(perf, results_path, new_perf)

            @printf("%8.2e seconds\n", total_time)
        end
    end
    return
end


text:
function complex_test()
    c1 = (2 + 3im)
    c2 = (4 + 5im)
    c3 = (c1 + c2)
    @assert(c3 == (6 + 8im))
    c4 = (c1 + 3)
    @assert(c4 == (5 + 3im))
    c5 = (c1 + 4im)
    @assert(c5 == (2 + 7im))
    c6 = (c3 - 2.3)
    @assert(c6 == (3.7 + 8im))
end

function main()
    complex_test()
    println(join(["OK"], " "))
end

main()


text:
input = joinpath(@__DIR__, "input")
lines = readlines(input)

m = hcat(map(line -> map(n -> parse(Int, n), collect(line)), lines)...)

function neighbors(m, idx)
    x, y = Tuple(idx)
    return filter(c -> checkbounds(Bool, m, c), [
        CartesianIndex(x + 1, y),
        CartesianIndex(x - 1, y),
        CartesianIndex(x, y + 1),
        CartesianIndex(x, y - 1)
    ])
end

function lowpoint(m, idx)
    ns = neighbors(m, idx)
    return all(i -> m[i] > m[idx], ns) ? idx : nothing
end

function basins(m, idx)
    basin = Set()
    prospects = [idx]
    while length(prospects) > 0
        curr = pop!(prospects)
        ch = m[curr]
        ns = neighbors(m, curr)
        for n in ns
            if m[n] > ch && m[n] < 9 && !(n in basin)
                push!(prospects, n)
            end
        end
        push!(basin, curr)
    end
    return length(basin)
end

lowpoints = filter(!isnothing, map(idx -> lowpoint(m, idx), CartesianIndices(m)))
p1 = sum([m[n] + 1 for n in lowpoints])
p2 = prod(sort(map(idx -> basins(m, idx), lowpoints), lt=!isless)[1:3])

println("-----------------------------------------------------------------------")
println("smoke basin -- part one :: $p1")
println("smoke basin -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 524)
@assert(p2 == 1235430)


text:
using VectorPartitions
using Test

@testset "VectorPartitions.jl" begin
using Combinatorics

function count_partition_iterator(n,k)
  P = all_vector_partitions(n,k)
  x = zeros(Int,n)
  count = 0
  for partition in P
    x .= partition
    count += 1
  end
  return count
end

function count_partitions_raw(n,k)
  P = all_vector_partitions(n,k)
  x = zeros(Int,n)
  next = iterate(P)
  count = 0
  while next !== nothing
      (x, state) = next
      next = iterate(P, state)
      count += 1
  end
  return count
end

function count_partitions_stateful(n,k)
  P = all_vector_partitions(n,k)
  next = iterate(P)
  if next !== nothing
    (_,state) = next
  else
    return nothing
  end
  count = 0
  while !state.done
    increment_vector_partition_state!(P,state)
    count += 1
  end
  return count
end

# some spot checks that we get the right number of partitions
@test count_partition_iterator(5,2) == stirlings2(5,2)
@test count_partitions_raw(5,2) == stirlings2(5,2)
@test count_partitions_stateful(5,2) == stirlings2(5,2)

@test count_partition_iterator(5,3) == stirlings2(5,3)
@test count_partitions_raw(5,3) == stirlings2(5,3)
@test count_partitions_stateful(5,3) == stirlings2(5,3)

# do we produce the right list of partitions in this case?
three_two_partitions = [[1,1,2],[1,2,1],[1,2,2]]
for (p1, p2) in zip(all_vector_partitions(3,2),three_two_partitions)
  @test all(p1 .== p2)
end

# test increment_vector_partition_state!
I = all_vector_partitions(7,3)
_, state = iterate(I)
state.partition .= [1, 1, 2, 3, 3, 3, 3]
state.upper_bound .= [1, 1, 2, 3, 3, 3, 3]
state.first_index .= [1, 3, 4]
increment_vector_partition_state!(I,state)
@test all(state.partition .== [1,2,1,3,1,1,1])
@test all(state.upper_bound .== [1,2,2,3,3,3,3])
@test all(state.first_index .== [1,2,4])

increment_vector_partition_state!(I,state)
@test all(state.partition .== [1,2,2,3,1,1,1])
@test all(state.upper_bound .== [1,2,2,3,3,3,3])
@test all(state.first_index .== [1,2,4])

# test increment_surjective_monotone_sequence!
seq = [1, 2, 2, 3, 4]
step_locations = [1,2,4,5]
stop = VectorPartitions.increment_surjective_monotone_sequence!(seq,4,step_locations)
@test all(seq .== [1, 2, 3, 3, 4])
@test all(step_locations .== [1,2,3,5])
@test stop == false

seq = [1,2,3,4,4]
step_locations = [1,2,3,4]
stop = VectorPartitions.increment_surjective_monotone_sequence!(seq,4,step_locations)
@test all(seq .== [1, 2, 3, 4, 4])
@test all(step_locations .== [1,2,3,4])
@test stop == true

end


text:
f = open("jawiki-country.txt", "r")

for line in readlines(f)
  if ismatch(r"^\|.+\s=\s",line)
    m = match(r"^\|(.+)\s=\s(.+)",line)
    println(m.captures[1], "\t",replace(m.captures[2], r"'{2,}", ""))
  end
end


text:
import ApproxFun: recA, recB, recC, recα, recβ, recγ
import ApproxFun: jacobirecA, jacobirecB, jacobirecC, jacobirecα, jacobirecβ, jacobirecγ

struct ComplexPlane <: Domain{ComplexF64} end

Base.in(x, d::ComplexPlane) = true

const ℂ = ComplexPlane()

canonicaldomain(d::ComplexPlane) = d
Base.reverse(C::ComplexPlane) = C

Base.intersect(a::ComplexPlane,b::ComplexPlane) = a
Base.union(a::ComplexPlane,b::ComplexPlane) = a

export JacobiQ, LegendreQ, WeightedJacobiQ

struct JacobiQ{D<:Domain,T} <: Space{D,T}
    a::T
    b::T
    domain::D
end
LegendreQ(domain) = JacobiQ(0.,0.,domain)
LegendreQ() = LegendreQ(ChebyshevInterval())
JacobiQ(a,b,d::Domain) = JacobiQ(promote(a,b)...,d)
JacobiQ(a,b,d) = JacobiQ(a,b,Domain(d))
JacobiQ(a,b) = JacobiQ(a,b,ChebyshevInterval())

domain(::JacobiQ) = ComplexPlane()

Base.promote_rule(::Type{JacobiQ{D,T}},::Type{JacobiQ{D,V}}) where {T,V,D} =
    JacobiQ{D,promote_type(T,V)}
convert(::Type{JacobiQ{D,T}},J::JacobiQ{D,V}) where {T,V,D} =
    JacobiQ{D,T}(J.a,J.b,J.domain)

const WeightedJacobiQ{D,T} = JacobiQWeight{JacobiQ{D,T},D}

WeightedJacobiQ(α,β,d::Domain) = JacobiQWeight(α,β,JacobiQ(β,α,d))
WeightedJacobiQ(α,β) = JacobiQWeight(α,β,JacobiQ(β,α))

spacescompatible(a::JacobiQ,b::JacobiQ) = a.a==b.a && a.b==b.b && a.domain == b.domain

function canonicalspace(S::JacobiQ)
    #if isapproxinteger(S.a+0.5) && isapproxinteger(S.b+0.5)
    #    Chebyshev(domain(S))
    #else
        # return space with parameters in (-1,0.]
        JacobiQ(mod(S.a,-1),mod(S.b,-1),S.domain)
    #end
end

for (REC,JREC) in ((:recα,:jacobirecα),(:recβ,:jacobirecβ),(:recγ,:jacobirecγ),
                   (:recA,:jacobirecA),(:recB,:jacobirecB),(:recC,:jacobirecC))
    @eval $REC(::Type{T},sp::JacobiQ,k) where {T} = $JREC(T,sp.a,sp.b,k)
end

stieltjes(f::Fun{<:Jacobi}) = Fun(LegendreQ(domain(f)),2coefficients(f,Legendre(domain(f))))
stieltjes(f::Fun{<:Chebyshev}) = Fun(LegendreQ(domain(f)),2coefficients(f,Legendre(domain(f))))

function stieltjes(f::Fun{<:JacobiWeight})
    # Jacobi parameters need to transform to:
    α,β = f.space.α,f.space.β
    Fun(WeightedJacobiQ(β,α,domain(f)),2coefficients(f,WeightedJacobi(β,α,domain(f))))
end

function evaluate(f::AbstractVector,S::JacobiQ{<:IntervalOrSegment},x)
    isinf(x) && return zero(promote_type(typeof(x), eltype(f)))
    stieltjesintervalrecurrence(S,f,mobius(S.domain,x))./2jacobiQweight(S.b,S.a,mobius(S.domain,x))
end
function evaluate(f::AbstractVector,S::JacobiQ{<:Curve},z::Number)
    isinf(z) && return zero(promote_type(typeof(z), eltype(f)))
    sum(evaluate(f,setcanonicaldomain(S),complexroots(S.domain.curve-z)))
end

# restrict to
function istieltjes(f::Fun{<:JacobiQ})
    # Jacobi parameters need to transform to:
    a,b = f.space.b,f.space.a
    d = f.space.domain
    a == b == 0 && return Fun(Legendre(d), coefficients(f)/2)
    Fun(WeightedJacobiQ(b,a,d),coefficients(g)/2)
end

stieltjes(f::Fun{<:ZeroSpace},z) = zero(promote_type(cfstype(f),typeof(z)))

istieltjes(f::Fun{<:ConstantSpace}) = Fun(ZeroSpace(),cfstype(f)[])
function istieltjes(f::Fun{<:SumSpace})
    is = istieltjes.(components(f))
    filter!(s -> !(space(s) isa ZeroSpace), is)
    Fun(is,PiecewiseSpace)
end

istieltjes(f::Fun{<:SumSpace{<:Tuple{<:Any,<:ConstantSpace}}}) =
    istieltjes(component(f,1))

istieltjes(f::Fun{<:ArraySpace}) = Fun(istieltjes.(Array(f)))

function union_rule(A::JacobiQ, B::JacobiQ)
    if A.domain == B.domain
        Jacobi(min(A.b,B.b),min(A.a,B.a),A.domain)
    else
        NoSpace()
    end
end
function maxspace_rule(A::JacobiQ,B::JacobiQ)
    A.domain == B.domain && return JacobiQ(max(A.b,B.b),max(A.a,B.a),A.domain)
    NoSpace()
end