text:
## sample.jl : collection of methods to add samples to an Estimator
#
# Methods to add samples to the Estimator.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

function update_samples(estimator::Estimator, n_opt)
    for τ in keys(n_opt)
        n_due = n_opt[τ] - nb_of_samples(estimator, τ)
        n_due > 0 && sample!(estimator, τ, n_due)
    end
end

function update_samples(estimator::Estimator{<:U}, n_opt)
    for τ in keys(n_opt)
        n_opt[τ] > 0 && sample!(estimator, τ, n_opt[τ])
    end
end

function sample!(estimator::Estimator{I}, index::Index, n::Integer) where I
    warm_up = !has_samples_at_index(estimator, index)
    estimator[:verbose] && print_sample!_header(estimator, index, n, warm_up)

    Istart = I <: U ? sum(values(nb_of_samples(estimator))) + 1 : nb_of_samples(estimator, index) + 1 # count zero level if U
    Iend = Istart + n - 1
    t = @elapsed parallel_sample!(estimator, index, Istart, Iend)

    add_to_nb_of_samples(estimator, index, n)
    add_to_total_work(estimator, index, n)
    add_to_total_time(estimator, index, t)

    estimator[:verbose] && print_sample!_footer()
end

#
# MC
#
function parallel_sample!(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index, Istart::Integer, Iend::Integer)

    m = estimator[:nb_of_uncertainties](index)
    f(i) = estimator.sample_function(index, transform.(view(distributions(estimator), 1:m), rand(m)))
    all_workers = workers()
    worker_idcs = 1:min(estimator[:nb_of_workers](index), nworkers())
    pool = CachingPool(all_workers[worker_idcs])
    #batch_size = ceil(Int, (Iend - Istart + 1)/length(worker_idcs))
    batch_size = 1
    retry_delays = ExponentialBackOff(n = 3)
    S = pmap(f, pool, Istart:Iend, batch_size = batch_size, retry_delays = retry_delays)
    clear!(pool)

    s_diff = first.(S)
    s = last.(S)

    append_samples!(estimator, index, s_diff, s)
end

function append_samples!(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index, s_diff, s)
    for n_qoi in 1:estimator[:nb_of_qoi]
        append_samples_diff!(estimator, n_qoi, index, getindex.(s_diff, n_qoi))
        append_samples!(estimator, n_qoi, index, getindex.(s, n_qoi))
    end
end

function append_samples!(estimator::Estimator{<:U, <:MC}, index::Index, s_diff, s)
    for idx in zero(index):index
        for n_qoi in 1:estimator[:nb_of_qoi]
            append_samples_diff!(estimator, n_qoi, idx, getindex.(getindex.(s_diff, Ref(idx + one(idx))), n_qoi))
            append_samples!(estimator, n_qoi, idx, getindex.(getindex.(s, Ref(idx + one(idx))), n_qoi))
            idx == zero(idx) && append_to_accumulator!(estimator, n_qoi, zeros(length(s_diff)))
            add_to_accumulator!(estimator, n_qoi, 1 / Prob(estimator, idx) * getindex.(getindex.(s_diff, Ref(idx + one(idx))), n_qoi))
        end
    end
end

#
# QMC
#
function parallel_sample!(estimator::Estimator{I, <:QMC}, index::Index, Istart::Integer, Iend::Integer) where I<:AbstractIndexSet

    m = estimator[:nb_of_uncertainties](index)
    x(i) = begin
        _index = I <: U ? zero(index) : index # take points from level 0 if unbiased
        point = get_point(generator(estimator, index, i[1]), i[2] - 1)
        length(point) < m && append!(point, rand(m - length(point)))
        view(point, 1:m)
    end
    f(i) = estimator.sample_function(index, transform.(view(distributions(estimator), 1:m), x(i)))
    all_workers = workers()
    worker_idcs = 1:min(estimator[:nb_of_workers](index), nworkers())
    pool = CachingPool(all_workers[worker_idcs])
    #batch_size = ceil(Int, (Iend - Istart + 1)/length(worker_idcs))
    batch_size = 1
    retry_delays = ExponentialBackOff(n = 3)
    S = pmap(f, pool, Iterators.product(1:estimator[:nb_of_shifts](index), Istart:Iend), batch_size = batch_size, retry_delays = retry_delays)
    clear!(pool)

    s_diff = first.(S)
    s = last.(S)

    append_samples!(estimator, index, s_diff, s)
end

function append_samples!(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index, s_diff, s)
    for n_shift in 1:estimator[:nb_of_shifts](index)
        for n_qoi in 1:estimator[:nb_of_qoi]
            append_samples_diff!(estimator, n_qoi, n_shift, index, getindex.(view(s_diff, n_shift, :), n_qoi))
            append_samples!(estimator, n_qoi, n_shift, index, getindex.(view(s, n_shift, :), n_qoi))
        end
    end
end

function append_samples!(estimator::Estimator{<:U, <:QMC}, index::Index, s_diff, s)
    for idx in zero(index):index
        for n_shift in 1:estimator[:nb_of_shifts](index)
            for n_qoi in 1:estimator[:nb_of_qoi]
                append_samples_diff!(estimator, n_qoi, n_shift, idx, getindex.(getindex.(view(s_diff, n_shift, :), Ref(idx + one(idx))), n_qoi))
                append_samples!(estimator, n_qoi, n_shift, idx, getindex.(getindex.(view(s, n_shift, :), Ref(idx + one(idx))), n_qoi))
                idx == zero(idx) && append_to_accumulator!(estimator, n_qoi, n_shift, zeros(size(s_diff, 2)))
                add_to_accumulator!(estimator, n_qoi, n_shift, 1 / Prob(estimator, idx) * getindex.(getindex.(view(s_diff, n_shift, :), Ref(idx + one(idx))), n_qoi))
            end
        end
    end
end


text:
"""
    NegativeMonomials(n) = [ x-> x^-k for k = 1:n ]
"""
struct NegativeMonomials <: AbstractBasis
    n::Int
end

Base.length(B::NegativeMonomials) = B.n

# approxpoints(B::NegativeMonomials) = TODO
# approxtransform(B::NegativeMonomials,f) = TODO

function iterate_basis(B::NegativeMonomials, x)
    length(B) == 0 && return nothing
    p = inv(x)
    return p,(2,p)
end
function iterate_basis(B::NegativeMonomials, x, (i,p))
    i > length(B) && return nothing
    p = x\p
    return p,(i+1,p)
end

evaluate_basis(B::NegativeMonomials,i,x) = x^-i

ext: jl


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

include("model_CapabilityInformation.jl")
include("model_CheckNameAvailabilityParameters.jl")
include("model_DataLakeStoreAccountListResult.jl")
include("model_DataLakeStoreAccountPropertiesBasic.jl")
include("model_DataLakeStoreAccountUpdateParameters.jl")
include("model_DataLakeStoreFirewallRuleListResult.jl")
include("model_DataLakeStoreTrustedIdProviderListResult.jl")
include("model_EncryptionConfig.jl")
include("model_EncryptionIdentity.jl")
include("model_ErrorDetails.jl")
include("model_FirewallRuleProperties.jl")
include("model_KeyVaultMetaInfo.jl")
include("model_NameAvailabilityInformation.jl")
include("model_Operation.jl")
include("model_OperationDisplay.jl")
include("model_OperationListResult.jl")
include("model_Resource.jl")
include("model_SubResource.jl")
include("model_TrustedIdProviderProperties.jl")
include("model_UpdateDataLakeStoreAccountProperties.jl")
include("model_UpdateEncryptionConfig.jl")
include("model_UpdateFirewallRuleParameters.jl")
include("model_UpdateFirewallRuleProperties.jl")
include("model_UpdateKeyVaultMetaInfo.jl")
include("model_UpdateTrustedIdProviderParameters.jl")
include("model_UpdateTrustedIdProviderProperties.jl")
include("model_DataLakeStoreAccount.jl")
include("model_DataLakeStoreAccountBasic.jl")
include("model_DataLakeStoreAccountProperties.jl")
include("model_FirewallRule.jl")
include("model_TrustedIdProvider.jl")


text:
# ---
# title: 49. Group Anagrams
# id: problem49
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, String
# link: <https://leetcode.com/problems/group-anagrams/description/>
# hidden: true
# ---
#
# Given an array of strings `strs`, group **the anagrams** together. You can
# return the answer in **any order**.
#
# An **Anagram** is a word or phrase formed by rearranging the letters of a
# different word or phrase, typically using all the original letters exactly
# once.
#
#
#
# **Example 1:**
#
#
#
#     Input: strs = ["eat","tea","tan","ate","nat","bat"]
#     Output: [["bat"],["nat","tan"],["ate","eat","tea"]]
#
#
# **Example 2:**
#
#
#
#     Input: strs = [""]
#     Output: [[""]]
#
#
# **Example 3:**
#
#
#
#     Input: strs = ["a"]
#     Output: [["a"]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= strs.length <= 104`
#   * `0 <= strs[i].length <= 100`
#   * `strs[i]` consists of lower-case English letters.
#
#
## @lc code=start
using LeetCode

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


text:
module FluxExperimental

using CUDA
using FileIO
using Flux
using Functors
using LinearAlgebra
using Requires
using Serialization
using Statistics
using Wandb
using Zygote

using Flux: hasaffine, ones32, zeros32, _isactive
using Functors: isleaf, children, _default_walk, functor
import NNlibCUDA: batchnorm, ∇batchnorm

function __init__()
    @require ChemistryFeaturization = "6c925690-434a-421d-aea7-51398c5b007a" begin
        include("graph_data.jl")
        # Atomic Graph Net Layers
        export batch_graph_data, BatchedAtomicGraph
    end
end

on_gpu(arr::CuArray) = true
on_gpu(arr::AbstractArray) = false

"""
    on_gpu(layer)

Checks if the parameters and states of the layer is in GPU memory. If one is in GPU and other is not,
then a warning is displayed. If there are no parameters and states then false is returned by default.
"""
on_gpu(layer) = on_gpu(destructure_parameters_states(layer)[1:2]...)

function on_gpu(p, s)
    if length(p) > 0 && length(s) > 0
        p_device = on_gpu(p)
        s_device = on_gpu(s)
        if !(p_device == s_device)
            @warn "Layer parameters and states are not on the same device"
        end
        return p_device && s_device
    elseif length(p) == 0
        return on_gpu(s)
    elseif length(s) == 0
        return on_gpu(p)
    else
        return false
    end
end

abstract type AbstractFluxLayer end

function Base.show(io::IO, l::AbstractFluxLayer)
    p, s, _ = destructure_parameters_states(l)
    device = on_gpu(p, s) ? "GPU" : "CPU"
    return print(io, string(typeof(l).name.name), "() ", string(length(p)), " Trainable Parameters & ",
                 string(length(s)), " States & Device = ", device)
end

"""
    debug_backward_pass(msg::String)

Prints the `msg` during the forward pass. During the backwards pass prints `∇(msg)`. Also prints
the time it took to reach the backward pass since the forward pass happended.
"""
@inline function debug_backward_pass(msg::String)
    start_time = Zygote.@ignore time()
    return Zygote.hook(Δ -> begin
                           println("∇(" * msg * ") -- Took $(time() - start_time) s")
                           return Δ
                       end, println(msg))
end

"""
    enable_fast_mode!()

Enable CUDA.FAST_MATH mode. Note this will reduce floating point precision but we rarely care about that.
"""
function enable_fast_mode!()
    CUDA.math_mode!(CUDA.FAST_MATH)
    return nothing
end

# Needed for destructure
include("layers/fchain.jl")

include("destructure.jl")
include("saving.jl")
include("utils.jl")

# Useful Layer Implementations
include("layers/norm_utils.jl")
include("layers/functional_wrappers.jl")
include("layers/weight_norm.jl")
include("layers/spectral_norm.jl")
include("layers/dropout.jl")
include("layers/utils.jl")
include("layers/agn.jl")
include("layers/group_norm.jl")
include("layers/batch_norm.jl")

# Extending the logging capabilities of wandb
include("wandb.jl")

# Common Utilities
export destructure_parameters_states, destructure_parameters
export enable_fast_mode!
export save_flux_model, load_flux_model
export debug_backward_pass
# Layers
export conv1x1, conv3x3, conv5x5, conv1x1_norm, conv3x3_norm, conv5x5_norm, conv_norm, downsample_module,
       upsample_module
export FChain
export WeightNorm, SpectralNorm
export VariationalHiddenDropout, update_is_variational_hidden_dropout_mask_reset_allowed
export AGNConv, AGNMaxPool, AGNMeanPool
export GroupNormV2, BatchNormV2
export ReshapeLayer, FlattenLayer, SelectDim, NoOpLayer
# Logging
export ParameterStateGradientWatcher

end


text:
using Test

using LinearAlgebra
using StaticArrays
using SatelliteAnalysis
import SatelliteToolbox: date_to_jd, sun_position_i

@testset "Beta angle" verbose = true begin
    include("./beta_angle.jl")
end

@testset "Lighting condition" verbose = true begin
    include("./lighting_condition.jl")
end

@testset "Miscellaneous" verbose = true begin
    include("./misc.jl")
end


text:
using Genie, Genie.Router, Genie.Renderer.Json, Genie.Requests
using HTTP
import Genie.Router: route
import Genie.Renderer.Json: json
using ArgParse

# Add the command line argument to optionally specify the port desired
s = ArgParseSettings()
@add_arg_table s begin
    "--port"
        help = "an option with an argument to specify the Port number"
        arg_type = Int
end

# Parse arguments
parsed_args = parse_args(ARGS, s)
port_num = parsed_args["port"]

Genie.config.run_as_server = true

# Parse script: replace Julia code arguments with the actual value inputs from Tableau.
function parse_script(message)

  script = message["script"]
  num_args = length(message["data"])

  # Replace the script _arg* placeholders with the actual values. Return parsed script.
  for index in 1:num_args
    arg = "_arg" * string(index)
    arg_value = message["data"][arg]
    script = replace(script, arg => arg_value)
  end

  script

end

### Set up /info endpoint for Tableau connection.
route("/info") do
  json("""
  {"description": "", "creation_time": "0", "state_path": "", "server_version": "0.8.7", "name": "Julia", "versions": {"v1": {"features": {}}}}
  """)
end

### Evaluate the Julia script with the data provided by Tableau. Return output in JSON format.
route("/evaluate", method = POST) do

  # Parse json from Tableau into Julia Dictionary.
  message = jsonpayload()

  # Parse the script's arguments with their actual values.
  parsed_script = parse_script(message)

  # Evaluate the parsed script. Return in JSON format to Tableau for interpretability.
  output = eval(Meta.parse(parsed_script))
  json(output)

end

# Can change 8011 port to port of interest
Genie.startup(port_num, async = false)

max_stars_repo_name: ApogeeIntegration/Tableau-Utilities

max_issues_repo_name: ApogeeIntegration/Tableau-Utilities

max_forks_repo_name: ApogeeIntegration/Tableau-Utilities


text:
abstract AbstractList{T}

immutable EmptyList{T} <: AbstractList{T} end
EmptyList() = EmptyList{Any}()

immutable PersistentList{T} <: AbstractList{T}
    head::T
    tail::AbstractList{T}
    length::Int
    PersistentList(head::T, tail::AbstractList{T}, length) =
        new(head, tail, length)
    PersistentList() = EmptyList{Any}()
    function PersistentList(v)
        v = reverse(v)
        list = EmptyList{T}()
        for el in v
            list = cons(el, list)
        end
        list
    end
end
PersistentList(itr) = PersistentList{eltype(itr)}(itr)

head(::EmptyList) = throw(BoundsError())
tail(::EmptyList) = throw(BoundsError())
head(l::PersistentList) = l.head
tail(l::PersistentList) = l.tail

Base.first(l::AbstractList) = head(l)

Base.length(::EmptyList) = 0
Base.length(l::PersistentList) = l.length

Base.isempty(::EmptyList) = true
Base.isempty(::PersistentList)      = false

cons{T}(val, ::EmptyList{T}) = PersistentList{T}(val, EmptyList{T}(), 1)
cons{T}(val, l::PersistentList{T})  = PersistentList{T}(val, l, length(l) + 1)
..(val, l::AbstractList) = cons(val, l)

Base.isequal(::EmptyList, ::EmptyList) = true
Base.isequal(l1::PersistentList, l2::PersistentList) =
    isequal(head(l1), head(l2)) && isequal(tail(l1), tail(l2))
==(::EmptyList, ::EmptyList) = true
==(l1::PersistentList, l2::PersistentList) =
    head(l1) == head(l2) && tail(l1) == tail(l2)

Base.start(l::AbstractList) = l
Base.done(::AbstractList, ::EmptyList) = true
Base.done(::AbstractList, ::PersistentList)      = false
Base.next(::AbstractList, l::PersistentList) = (head(l), tail(l))

Base.isequal(a::AbstractArray, l::PersistentList) = isequal(l, a)
Base.isequal(l::PersistentList, a::AbstractArray) =
    isequal(length(l), length(a)) && all((el) -> el[1] == el[2], zip(l, a))
==(a::AbstractArray, l::PersistentList) = isequal(l, a)
==(l::PersistentList, a::AbstractArray) = isequal(l, a)

Base.map(f::(@compat Union{Function, DataType}), e::EmptyList) = e
Base.map(f::(@compat Union{Function, DataType}), l::PersistentList) = cons(f(head(l)), map(f, tail(l)))

Base.reverse(e::EmptyList) = e
function Base.reverse{T}(l::PersistentList{T})
    reversed = EmptyList{T}()
    for val in l
        reversed = val..reversed
    end
    reversed
end

@compat Base.show(io::IO, ::MIME"text/plain", ::EmptyList) = print(io, "()")
@compat function Base.show{T}(io::IO, ::MIME"text/plain", l::PersistentList{T})
    print(io, "$T($(head(l))")
    for val in tail(l)
        print(io, ", $val")
    end
    print(io, ")")
end


text:
struct WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, T<:Real}
    Hwedge::Vector{T}
    Hv::Vector{T}
    Hextra::Vector{T}
    zₐpowers::Vector{Complex{T}}
    zᵧpowers::Vector{Complex{T}}
    z::Vector{Complex{T}}
    function WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, T}() where {ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, T<:Real}
        if ℓₘᵢₙ < 0
            throw(DomainError("ℓₘᵢₙ = $ℓₘᵢₙ",
                    "This function is only valid for non-negative values of ℓ."
            ))
        end
        if ℓₘₐₓ < 0
            throw(DomainError("ℓₘₐₓ = $ℓₘₐₓ",
                    "This function is only valid for non-negative values of ℓ."
            ))
        end
        if ℓₘᵢₙ > ℓₘₐₓ
            throw(DomainError("ℓₘᵢₙ = $ℓₘᵢₙ",
                    "The input ℓₘᵢₙ must be no greater than the input ℓₘₐₓ=$ℓₘₐₓ."
            ))
        end
        new{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, T}(
            zeros(T, WignerHsize(m′ₘₐₓ, ℓₘₐₓ)),
            zeros(T, (ℓₘₐₓ + 1)^2),
            zeros(T, ℓₘₐₓ + 2),
            zeros(Complex{T}, ℓₘₐₓ + 1),
            zeros(Complex{T}, ℓₘₐₓ + 1),
            zeros(Complex{T}, 3)
        )
    end
end

"""
   WignerMatrixCalculator(ℓmin, ℓmax, [m′max], [T])

Construct a `WignerMatrixCalculator` to use in calculating ``𝔇``, ``d``, or
``{}_{s}Y_{ℓ,m}``.

If `m′max` is given, values may be calculated with `ℓ ∈ ℓmin:ℓmax` and `m ∈
-ℓ:ℓ` but only `m′ ∈ ℓmin:m′max.  This will reduce the storage and computation
costs.  Note that if your goal is just to compute spin-weighted spherical
harmonics up to some spin `s`, you should use `m′max=abs(s)`.

If `T<:Real` is given, all calculations will be done using that type; the
default is `Float64`.

See also:
  * [`D!`](@ref) to compute the full ``𝔇`` matrix
  * [`d!`](@ref) to compute the small ``d`` matrix
  * [`Y!`](@ref) to compute spin-weighted spherical harmonic values

"""
WignerMatrixCalculator(ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ=typemax(Int), T=Float64) = WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, min(ℓₘₐₓ, abs(m′ₘₐₓ)), T}()

# Access type parameters
ℓₘᵢₙ(::Type{<:WignerMatrixCalculator{ℓₘᵢₙW, ℓₘₐₓ, m′ₘₐₓ, T}}) where {ℓₘᵢₙW, ℓₘₐₓ, m′ₘₐₓ, T<:Real} = ℓₘᵢₙW
ℓₘₐₓ(::Type{<:WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓW, m′ₘₐₓ, T}}) where {ℓₘᵢₙ, ℓₘₐₓW, m′ₘₐₓ, T<:Real} = ℓₘₐₓW
m′ₘₐₓ(::Type{<:WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓW, T}}) where {ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓW, T<:Real} = m′ₘₐₓW
T(::Type{<:WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, TW}}) where {ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, TW<:Real} = TW
ℓₘᵢₙ(::WignerMatrixCalculator{ℓₘᵢₙW, ℓₘₐₓ, m′ₘₐₓ, T}) where {ℓₘᵢₙW, ℓₘₐₓ, m′ₘₐₓ, T<:Real} = ℓₘᵢₙW
ℓₘₐₓ(::WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓW, m′ₘₐₓ, T}) where {ℓₘᵢₙ, ℓₘₐₓW, m′ₘₐₓ, T<:Real} = ℓₘₐₓW
m′ₘₐₓ(::WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓW, T}) where {ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓW, T<:Real} = m′ₘₐₓW
T(::WignerMatrixCalculator{ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, TW}) where {ℓₘᵢₙ, ℓₘₐₓ, m′ₘₐₓ, TW<:Real} = TW
# ℓₘᵢₙ(w) = ℓₘᵢₙ(typeof(w))
# ℓₘₐₓ(w) = ℓₘₐₓ(typeof(w))
# m′ₘₐₓ(w) = m′ₘₐₓ(typeof(w))
# T(w) = T(typeof(w))

# Access pre-computed array sizes
WignerHsize(::Type{W}) where {W<:WignerMatrixCalculator} = WignerHsize(m′ₘₐₓ(W), ℓₘₐₓ(W))
Wignerdsize(::Type{W}) where {W<:WignerMatrixCalculator} = Wignerdsize(ℓₘᵢₙ(W), m′ₘₐₓ(W), ℓₘₐₓ(W))
WignerDsize(::Type{W}) where {W<:WignerMatrixCalculator} = WignerDsize(ℓₘᵢₙ(W), m′ₘₐₓ(W), ℓₘₐₓ(W))
Ysize(::Type{W}) where {W<:WignerMatrixCalculator} = Ysize(ℓₘᵢₙ(W), ℓₘₐₓ(W))
WignerHsize(w::WignerMatrixCalculator) = WignerHsize(typeof(w))
Wignerdsize(w::WignerMatrixCalculator) = Wignerdsize(typeof(w))
WignerDsize(w::WignerMatrixCalculator) = WignerDsize(typeof(w))
Ysize(w::WignerMatrixCalculator) = Ysize(typeof(w))

# Access pre-computed helper data
function a(::Type{W}) where {W<:WignerMatrixCalculator}
    [√((n+1+m) * (n+1-m) / T(W)((2n+1)*(2n+3))) for n in 0:ℓₘₐₓ(W)+1 for m in 0:n]
end
a(w) = a(typeof(w))
function b(::Type{W}) where {W<:WignerMatrixCalculator}
    [(m<0 ? -1 : 1) * √((n-m-1) * (n-m) / T(W)((2n-1)*(2n+1))) for n in 0:ℓₘₐₓ(W)+1 for m in -n:n]
end
b(w) = b(typeof(w))
function d(::Type{W}) where {W<:WignerMatrixCalculator}
    [(m<0 ? -1 : 1) * (√T(W)((n-m) * (n+m+1))) / 2 for n in 0:ℓₘₐₓ(W)+1 for m in -n:n]
end
d(w) = d(typeof(w))

function Base.show(io::IO, w::WignerMatrixCalculator)
    print(io, typeof(w))
end


text:
using ForwardDiff
using LinearAlgebra
using StaticArrays

norm²(x) = sum( z -> z^2, x)
# this is a bit faster than norm (but don't worry, norm()^2 is fine)

function linesearch(f, Δx, m², x, α, β; t_min = 1e-5)
	t = 1.0
	fx = f(x)
	z = x + t*Δx
	while t > t_min
		z .= x + t*Δx
		if f(z) <= fx + α * t * m²
			return t
		end
		t = t*β
	end

	return NaN64
end

using ForwardDiff: GradientConfig, Chunk

function gradientdecent(f, x0, α, β, ε;
							max_steps = 1e4, traj = true,
							cfg = GradientConfig(f, x0, Chunk{length(x0)}()))
	x = copy(x0)
	i = 0
	∇f = similar(x0)

	if traj
		X = [x0]
	end

	Δx = -∇f
	#@show typeof(Δx) typeof(∇f)
	while i < max_steps
		ForwardDiff.gradient!(∇f, f, x, cfg)
		Δx .= -∇f

		m² = norm²(∇f)
		if m² < ε;	break  end

		t = linesearch(f, Δx, m², x, α, β)
		x .= x + t * Δx

		if isnan(t)
			@error "Gradient decent method did not find a solution."
			return missing
		end

		if traj; push!(X, copy(x)) end

		i += 1
	end

	if i == max_steps
		@error "Gradient decent method did not find a solution."
		return missing
	else
		return traj ? (; min = x, traj = reduce(hcat,X)) : x
	end
end

begin
	f(x) = exp(x[1]+3x[2]-0.1) + exp(x[1]-3x[2]-0.1) + exp(-x[1]-0.1)
	α = 0.25
	β = 0.5
	ε = 1e-15
	x0 = [-1.8, 1.2]

	sol = gradientdecent(f, x0, α, β, ε; traj = true)
	size(sol.traj)
end

@elapsed  let
	f(x) = exp(x[1]+3x[2]-0.1) + exp(x[1]-3x[2]-0.1) + exp(-x[1]-0.1)
	α = 0.25
	β = 0.5
	ε = 1e-15
	x0 = @MVector[-1.8, 1.2]
	for i in 1:1000
		gradientdecent(f, x0, 0.25, 0.5, ε; traj = false)
	end
end

using GLMakie
let
	fig = Figure()
	ax = Axis(fig[1,1], title = "Trajectory of gradient decent method")

	lines!(sol.traj[1,:], sol.traj[2,:], linewidth = 2, color = :black)

	xs = LinRange(-2.5,0.5,100)
	ys = LinRange(-1.5,1.5,100)
	z = [f((x,y)) for x in xs, y in ys]
	heatmap!(xs, ys, z, colormap = :heat)
	contour!(xs, ys, z, levels = 100 )

	fig
end

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
#######################################################################
#=
Section of code dealing with unconstrained minimization problems.
Some functions are reused in the equality constrained section
=#
#######################################################################

#################
#Gradient Descent
#################

#Use Calculus.jl library to compute the gradient numerically
function Compute_Gradient(_f, _x)
    return Calculus.gradient(g -> _f(g),_x)
end

#Use Calculus.jl library to compute the hessian numerically
function Compute_Hessian(_f, _x)
    return Calculus.hessian(h -> _f(h),_x)
end

#Choose t using backgracking line search (p. 464, alg. 9.2)
function Back_Line_Search(_x, _f, _fx, _∇f, _Δx, _α, _β, _n)

    _t = 1

    while _f((_x + _t*_Δx)) > _fx + _α*_t*transpose(_∇f)*_Δx
        _t *= _β
    end

    return _t
end

#Main solution method for Unconstrained Gradient Descent Method (p. 466, alg 9.3)
function Unconstrained_Grad_Descent(_f, _x, _n, _α, _β, _ϵ, maxIt, verbose)
    it = 0;
    t = 1
    ∇f = Compute_Gradient(_f, _x)
    normGrad = norm(∇f)

    val = _f(_x)

    #Stopping conditions: norm(∇f) <= η or iterations exceeds maximum iteration tolerance
    while normGrad > _ϵ && it <= maxIt
        ∇f = Compute_Gradient(_f, _x) #Evaluate gradient at current point
        normGrad = norm(∇f) #Compute euclidean norm using LinearAlgebra library

        #Perform backgracking line search to pick t
        bls = Back_Line_Search(_x, _f, val, ∇f, -∇f, _α, _β, _n)
        t = bls[1]

        #update step
        _x -= t*∇f

        val = _f(_x)

        if(verbose)
            print("\n\nIteration: ", it)
            print("\nGradient Norm = ", normGrad)
        #    print("\nDescent Direction = ", -∇fVal)
            print("\nt = ", t)
        #    print("\nx = ",_x)
            print("\nf(x) = ", val)
        end

        it += 1
    end

    if it > maxIt
        print("\n\nMax Iterations Reached!")
    end

    return _x, val
end

#################
#Newton Method
#################

#Computes Newton step (p. 487, alg. 9.5)
function Newton_Step(_∇f, _∇2f, _x, _n)
    return _∇2f\-_∇f
end

#Computes Newton decrement (p. 487, alg. 9.5)
function Newton_Decrement(_∇2f, _Δxnt)
    return transpose(_Δxnt)*_∇2f*_Δxnt
end

#Main solver method for Unconstrained Newton Method (p. 487, alg. 9.5)
function Unconstrained_Newton(_f, _x, _n, _α, _β, _ϵ, maxIt, verbose)

    λ2 = 1 #prime λ
    it = 0
    val = _f(_x)

    while λ2 > _ϵ && it <= maxIt

        #Objective, Gradient, and Hessian
        val = _f(_x)
        ∇f = Compute_Gradient(_f, _x)
        ∇2f = Compute_Hessian(_f, _x)

        #Newton step
        Δxnt = Newton_Step(∇f, ∇2f, _x, _n)

        #Newton decrement
        λ2 = Newton_Decrement(∇2f, Δxnt)

        #Choose t with backtracking line search
        t = Back_Line_Search(_x, _f, val, ∇f, Δxnt, _α, _β, _n)

        #Print data
        if verbose
            print("\n\nIteration: ", it)
        #    print("\nNewton step = ", Δxnt)
            print("\nλ^2 = ", λ2)
            print("\nt = ", t)
        #    print("\nx = ", _x)
            print("\nf(x) = ", val)
        end

        #Update x
        _x += t*Δxnt

        it+=1
    end

    if it > maxIt
        print("\n\nMax Iterations Reached!")
    end

    return _x, val
end

#=
This function combines the consistent convergent rate of gradient descent
with the fast local convergence of Newton's Method by using gradient
descent at first, and then switching to Newton's method when the solution
gets within η > ϵ. Then, Newton's method gets the solution to within ϵ
=#
function Solve_UMP(_f, _x, _n, _α, _β, _η, _ϵ, maxIt, verbose)
    sol = Unconstrained_Grad_Descent(_f, _x, _n, _α, _β, _η, maxIt, verbose)

    print("\n\nCurrent solution is ", sol[2] ,"\nSwitching to Newton's Method\n")
    #use output point of gradient descent solution as starting point for Newton's Method
    return Unconstrained_Newton(_f, sol[1], _n, _α, _β, _ϵ, maxIt, verbose)
end

#######################################################################
#=
Section of code dealing with equality constrained minimization problems.
=#
#######################################################################

#################
#Feasible Start NM
#################

#Computes Newton Step (p. 526 (10.11))
function Compute_Newton_Step(_f, _∇f, _∇2f, _A, _At, _x, _n, _m)
    return ([_∇2f _At; _A zeros(_m,_m)]  \ [-_∇f; zeros(_m)])[1:_n]
end

#Compute Newton Decrement (p. 527 (10.12))
function Compute_Newton_Decrement(_Δxn, _∇2f)
    return abs(transpose(_Δxn) * _∇2f * _Δxn)
end

#Main solver function for Feasible Start Newton Method (p. 528, alg. 10.1)
function Newton_FS(_f, _x, _A, _n, _m, _α, _β, _ϵ, maxIt, verbose)

    obj = Float64
    λ2 = 10 #prime λ with a big number to ensure first iteration is run
    it = 0
    _At = transpose(_A)

    while λ2 / 2 > _ϵ && it <= maxIt

        obj = _f(_x)
        ∇f = Compute_Gradient(_f, _x)
        ∇2f = Compute_Hessian(_f, _x)
        Δxn = Compute_Newton_Step(_f, ∇f, ∇2f, _A, _At, _x, _n, _m)

        λ2 = Compute_Newton_Decrement(Δxn, ∇2f) #returns λ^2

        t = Back_Line_Search(_x, _f, obj, ∇f, Δxn, _α, _β, _n)

        #Update x
        _x += t*Δxn

        if verbose
            print("\n\nIteration ", it)
            print("\nλ^2 =         ",λ2)
            print("\nObjective =   ", obj)
        end

        it+=1
    end

    return _x, obj
end

#################
#Infeasible Start NM
#################

# This function computes the dual. (p. 531)
function Compute_Dual(_∇f, _At)
    return _At \ -_∇f
end

#This function computes the residual see page 533 (10.21)
function Compute_Residual(_∇f, _At, _v, _vr)
    return [_∇f + _At*_v; _vr]
end

#This function computes the residual just like Compute_Residual but in a way that we can use it for backgracking line search.
function R(_x, _v, _At, _A, _b, _f)
    return [Compute_Gradient(_f, _x) + _At*_v; _A*_x - _b]
end

#This function computes the Newton and Dual Step. See page 533 (10.22)
function Compute_Steps(_∇f, _∇2f, _x, _A, _At, _vr, _p, _v, _res)
    return [_∇2f _At; _A zeros(_p,_p)] \ -_res
end

#Backgracking line search for infeasible start method. (p. 534, Alg. 10.2)
function BLS_IF(_x, _v, _A, _At, _b, _vr, _f, _∇f, _Δxnt, _Δvnt, _α, _β, _nrm)
    _t = 1

    while norm(R(_x + _t*_Δxnt, _v + _t*_Δvnt, _At, _A, _b, _f)) > (1 - _α*_t)* _nrm
        _t *= _β
    end

    return _t
end

#Main solver function
function Newton_IFS(_f, _x, _A, _b, _n, _m, _α, _β, _ϵ, maxIt, verbose)

    #Initialize necessary values for first iteration
    At = transpose(_A)
    vr = _A*_x - _b
    ∇f = Compute_Gradient(_f, _x)
    _v = Compute_Dual(∇f, At)
    it = 0
    nrm = 100

    while  (norm(vr) > _ϵ || nrm > _ϵ) && it < maxIt

        it += 1

        ∇f = Compute_Gradient(_f, _x)
        ∇2f = Compute_Hessian(_f, _x)
        _v = Compute_Dual(∇f, At)

        vr = _A*_x - _b
        res = Compute_Residual(∇f, At, _v, vr)

        Δ = Compute_Steps(∇f, ∇2f, _x, _A, At, vr, _m, _v, res)
        Δxnt = Δ[1:_n]
        Δvnt = Δ[_n+1:_n+_m]

        nrm = norm(res)
        t = BLS_IF(_x, _v, _A, At, _b, vr, _f, ∇f, Δxnt, Δvnt, _α, _β, nrm)

        _x += t*Δxnt
        _v += t*Δvnt

        if verbose
            print("\n\nIteration: ", it)
            print("\nt = ", t)
            print("\nnorm(r) = ", nrm)
    #        print("\nx = ", x)
            print("\nObjective = ", _f(_x))
    #        print("\nΔx = ", Δxnt)
        end

    end

    return _x, _f(_x)
end

max_stars_repo_name: ihawn/EqualityConstrainedMinimization

max_issues_repo_name: ihawn/EqualityConstrainedMinimization

max_forks_repo_name: ihawn/EqualityConstrainedMinimization


text:
using DReal

ctx = global_context()
X = Var(Float64, 0.0, 1.0)
exgold = (X < 0.5) | (X > 0.9)
notex = (!)(ctx, exgold)

ex = exgold
push_ctx!(ctx)
lb_constraint = (>=)(ctx, X, 0.0)
DReal.add!(ctx, lb_constraint)
ub_constraint = (<=)(ctx, X, 1.0)
DReal.add!(ctx, ub_constraint)
DReal.add!(ctx, ex)
result = is_satisfiable(ctx)
pop_ctx!(ctx)

ex = notex
push_ctx!(ctx)
lb_constraint = (>=)(ctx, X, 0.0)
DReal.add!(ctx, lb_constraint)
ub_constraint = (<=)(ctx, X, 1.0)
DReal.add!(ctx, ub_constraint)
DReal.add!(ctx, ex)
result = is_satisfiable(ctx)
pop_ctx!(ctx)

ex = exgold
push_ctx!(ctx)
lb_constraint = (>=)(ctx, X, 0.0)
DReal.add!(ctx, lb_constraint)
ub_constraint = (<=)(ctx, X, 0.5)
DReal.add!(ctx, ub_constraint)
DReal.add!(ctx, ex)
result = is_satisfiable(ctx)
pop_ctx!(ctx)

ex = notex
push_ctx!(ctx)
lb_constraint = (>=)(ctx, X, 0.0)
DReal.add!(ctx, lb_constraint)
ub_constraint = (<=)(ctx, X, 0.5)
DReal.add!(ctx, ub_constraint)
DReal.add!(ctx, ex)

DReal.is_satisfiable()


text:
# Solves a basic eigenvalue problem defined through different NEP types through NLEIGS

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using NonlinearEigenproblems.RKHelper
using Test
using SparseArrays

struct CustomNLEIGSNEP <: NEP
    n::Int  # problem size; this is the only required field in a custom NEP type when used with NLEIGS
end

nep_types_test_n = 2
nep_types_test_B = Vector{Matrix{Float64}}([[1 3; 5 6], [3 4; 6 6]])
nep_types_test_C = Vector{Matrix{Float64}}([[1 0; 0 1]])

# implement a few methods used by the solver
import NonlinearEigenproblems.NEPCore.compute_Mder
import NonlinearEigenproblems.NEPCore.compute_Mlincomb
import Base.size
nep_types_test_pep = PEP([nep_types_test_B; nep_types_test_C])
compute_Mder(::CustomNLEIGSNEP, λ::Number) = compute_Mder(nep_types_test_pep, λ)
compute_Mlincomb(::CustomNLEIGSNEP, λ::Number, x::AbstractVecOrMat) = compute_Mlincomb(nep_types_test_pep, λ, x)
size(::CustomNLEIGSNEP, _) = nep_types_test_n

function nleigs_nep_types()
    n = nep_types_test_n
    B = nep_types_test_B
    C = nep_types_test_C
    f = [λ -> λ^2]
    Σ = [-10.0-2im, 10-2im, 10+2im, -10+2im]

    # define and solve the same problem in many different ways
    problems = [
        ("SPMF_NEP", SPMF_NEP([B; C], [λ -> λ^0; λ -> λ; λ -> λ^2])),
        ("PEP", PEP([B; C])),
        ("PEP + SPMF", SumNEP(PEP(B), SPMF_NEP(C, f))),
        ("PEP + LowRankFactorizedNEP", SumNEP(PEP(B), LowRankFactorizedNEP([LowRankMatrixAndFunction(sparse(C[1]), f[1])]))),
        ("Custom NEP type", CustomNLEIGSNEP(n))]

    for problem in problems
        @bench @testset "$(problem[1])" begin
            lambda, X = nleigs(problem[2], Σ, maxit=10, v=ones(n).+0im, blksize=5)
            verify_lambdas(4, problem[2], lambda, X)
        end
    end
end

@testset "NLEIGS: NEP Types" begin
    nleigs_nep_types()
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function apply_smoothing!(mesh::StructuredMesh{1}, alpha, alpha_tmp, dg, cache)
  # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha
  # Copy alpha values such that smoothing is indpedenent of the element access order
  alpha_tmp .= alpha

  # So far, alpha smoothing doesn't work for non-periodic initial conditions for structured meshes.
  @assert isperiodic(mesh) "alpha smoothing for structured meshes works only with periodic initial conditions so far"

  # Loop over elements, because there is no interface container
  for element in eachelement(dg,cache)
    # Get neighboring element ids
    left  = cache.elements.left_neighbors[1, element]

    # Apply smoothing
    alpha[left]     = max(alpha_tmp[left],    0.5 * alpha_tmp[element], alpha[left])
    alpha[element]  = max(alpha_tmp[element], 0.5 * alpha_tmp[left],    alpha[element])
  end
end

end # @muladd


text:
addprocs(4)

using ProgressMeter
using ParallelProgressMeter

@parallelprogress for i in 1:50
    sleep(0.1)
end

@parallelprogress 0.5 "Testing..." for i in 1:50
    sleep(0.1)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module PeriodsTest

using Dates
using Test

@testset "basic arithmetic" begin
    @test -Dates.Year(1) == Dates.Year(-1)
    @test Dates.Year(1) > Dates.Year(0)
    @test (Dates.Year(1) < Dates.Year(0)) == false
    @test Dates.Year(1) == Dates.Year(1)
    @test Dates.Year(1) != 1
    @test Dates.Year(1) + Dates.Year(1) == Dates.Year(2)
    @test Dates.Year(1) - Dates.Year(1) == zero(Dates.Year)
    @test 1 == one(Dates.Year)
    @test_throws MethodError Dates.Year(1) * Dates.Year(1) == Dates.Year(1)
    t = Dates.Year(1)
    t2 = Dates.Year(2)
    @test ([t, t, t, t, t] .+ Dates.Year(1)) == ([t2, t2, t2, t2, t2])
    @test (Dates.Year(1) .+ [t, t, t, t, t]) == ([t2, t2, t2, t2, t2])
    @test ([t2, t2, t2, t2, t2] .- Dates.Year(1)) == ([t, t, t, t, t])
    @test_throws MethodError ([t, t, t, t, t] .* Dates.Year(1)) == ([t, t, t, t, t])
    @test ([t, t, t, t, t] * 1) == ([t, t, t, t, t])
    @test ([t, t, t, t, t] .% t2) == ([t, t, t, t, t])
    @test div.([t, t, t, t, t], Dates.Year(1)) == ([1, 1, 1, 1, 1])
    @test mod.([t, t, t, t, t], Dates.Year(2)) == ([t, t, t, t, t])
    @test [t, t, t] / t2 == [0.5, 0.5, 0.5]
    @test abs(-t) == t
    @test sign(t) == sign(t2) == 1
    @test sign(-t) == sign(-t2) == -1
    @test sign(Dates.Year(0)) == 0
end
@testset "div/mod/gcd/lcm/rem" begin
    @test Dates.Year(10) % Dates.Year(4) == Dates.Year(2)
    @test gcd(Dates.Year(10), Dates.Year(4)) == Dates.Year(2)
    @test lcm(Dates.Year(10), Dates.Year(4)) == Dates.Year(20)
    @test div(Dates.Year(10), Dates.Year(3)) == 3
    @test div(Dates.Year(10), Dates.Year(4)) == 2
    @test div(Dates.Year(10), 4) == Dates.Year(2)
    @test Dates.Year(10) / Dates.Year(4) == 2.5

    @test mod(Dates.Year(10), Dates.Year(4)) == Dates.Year(2)
    @test mod(Dates.Year(-10), Dates.Year(4)) == Dates.Year(2)
    @test mod(Dates.Year(10), 4) == Dates.Year(2)
    @test mod(Dates.Year(-10), 4) == Dates.Year(2)

    @test rem(Dates.Year(10), Dates.Year(4)) == Dates.Year(2)
    @test rem(Dates.Year(-10), Dates.Year(4)) == Dates.Year(-2)
    @test rem(Dates.Year(10), 4) == Dates.Year(2)
    @test rem(Dates.Year(-10), 4) == Dates.Year(-2)
end

y = Dates.Year(1)
m = Dates.Month(1)
w = Dates.Week(1)
d = Dates.Day(1)
h = Dates.Hour(1)
mi = Dates.Minute(1)
s = Dates.Second(1)
ms = Dates.Millisecond(1)
us = Dates.Microsecond(1)
ns = Dates.Nanosecond(1)
emptyperiod = ((y + d) - d) - y
@testset "Period arithmetic" begin
    @test Dates.Year(y) == y
    @test Dates.Month(m) == m
    @test Dates.Week(w) == w
    @test Dates.Day(d) == d
    @test Dates.Hour(h) == h
    @test Dates.Minute(mi) == mi
    @test Dates.Second(s) == s
    @test Dates.Millisecond(ms) == ms
    @test Dates.Microsecond(us) == us
    @test Dates.Nanosecond(ns) == ns
    @test Dates.Year(convert(Int8, 1)) == y
    @test Dates.Year(convert(UInt8, 1)) == y
    @test Dates.Year(convert(Int16, 1)) == y
    @test Dates.Year(convert(UInt16, 1)) == y
    @test Dates.Year(convert(Int32, 1)) == y
    @test Dates.Year(convert(UInt32, 1)) == y
    @test Dates.Year(convert(Int64, 1)) == y
    @test Dates.Year(convert(UInt64, 1)) == y
    @test Dates.Year(convert(Int128, 1)) == y
    @test Dates.Year(convert(UInt128, 1)) == y
    @test Dates.Year(convert(BigInt, 1)) == y
    @test Dates.Year(convert(BigFloat, 1)) == y
    @test Dates.Year(convert(Complex, 1)) == y
    @test Dates.Year(convert(Rational, 1)) == y
    @test Dates.Year(convert(Float16, 1)) == y
    @test Dates.Year(convert(Float32, 1)) == y
    @test Dates.Year(convert(Float64, 1)) == y
    @test y == y
    @test m == m
    @test w == w
    @test d == d
    @test h == h
    @test mi == mi
    @test s == s
    @test ms == ms
    @test us == us
    @test ns == ns
    y2 = Dates.Year(2)
    @test y < y2
    @test y2 > y
    @test y != y2

    @test Dates.Year(Int8(1)) == y
    @test Dates.Year(UInt8(1)) == y
    @test Dates.Year(Int16(1)) == y
    @test Dates.Year(UInt16(1)) == y
    @test Dates.Year(Int(1)) == y
    @test Dates.Year(UInt(1)) == y
    @test Dates.Year(Int64(1)) == y
    @test Dates.Year(UInt64(1)) == y
    @test Dates.Year(UInt128(1)) == y
    @test Dates.Year(UInt128(1)) == y
    @test Dates.Year(big(1)) == y
    @test Dates.Year(BigFloat(1)) == y
    @test Dates.Year(float(1)) == y
    @test Dates.Year(Float32(1)) == y
    @test Dates.Year(Rational(1)) == y
    @test Dates.Year(complex(1)) == y
    @test_throws InexactError Dates.Year(BigFloat(1.2)) == y
    @test_throws InexactError Dates.Year(1.2) == y
    @test_throws InexactError Dates.Year(Float32(1.2)) == y
    @test_throws InexactError Dates.Year(3//4) == y
    @test_throws InexactError Dates.Year(complex(1.2)) == y
    @test_throws InexactError Dates.Year(Float16(1.2)) == y
    @test Dates.Year(true) == y
    @test Dates.Year(false) != y
    @test_throws MethodError Dates.Year(:hey) == y
    @test Dates.Year(real(1)) == y
    @test_throws InexactError Dates.Year(m) == y
    @test_throws MethodError Dates.Year(w) == y
    @test_throws MethodError Dates.Year(d) == y
    @test_throws MethodError Dates.Year(h) == y
    @test_throws MethodError Dates.Year(mi) == y
    @test_throws MethodError Dates.Year(s) == y
    @test_throws MethodError Dates.Year(ms) == y
    @test Dates.Year(Dates.Date(2013, 1, 1)) == Dates.Year(2013)
    @test Dates.Year(Dates.DateTime(2013, 1, 1)) == Dates.Year(2013)
    @test typeof(y + m) <: Dates.CompoundPeriod
    @test typeof(m + y) <: Dates.CompoundPeriod
    @test typeof(y + w) <: Dates.CompoundPeriod
    @test typeof(y + d) <: Dates.CompoundPeriod
    @test typeof(y + h) <: Dates.CompoundPeriod
    @test typeof(y + mi) <: Dates.CompoundPeriod
    @test typeof(y + s) <: Dates.CompoundPeriod
    @test typeof(y + ms) <: Dates.CompoundPeriod
    @test typeof(y + us) <: Dates.CompoundPeriod
    @test typeof(y + ns) <: Dates.CompoundPeriod
    @test y > m
    @test d < w
    @test mi < h
    @test ms < h
    @test ms < mi
    @test us < ms
    @test ns < ms
    @test ns < us
    @test ns < w
    @test us < w
    @test typemax(Dates.Year) == Dates.Year(typemax(Int64))
    @test typemax(Dates.Year) + y == Dates.Year(-9223372036854775808)
    @test typemin(Dates.Year) == Dates.Year(-9223372036854775808)
end
@testset "Period-Real arithmetic" begin
    @test_throws MethodError y + 1 == Dates.Year(2)
    @test_throws MethodError y + true == Dates.Year(2)
    @test_throws InexactError y + Dates.Year(1.2)
    @test y + Dates.Year(1f0) == Dates.Year(2)
    @test y * 4 == Dates.Year(4)
    @test y * 4f0 == Dates.Year(4)
    @test_throws InexactError y * 3//4 == Dates.Year(1)
    @test div(y, 2) == Dates.Year(0)
    @test_throws MethodError div(2, y) == Dates.Year(2)
    @test div(y, y) == 1
    @test y*10 % Dates.Year(5) == Dates.Year(0)
    @test_throws MethodError (y > 3) == false
    @test_throws MethodError (4 < y) == false
    @test 1 != y
    t = [y, y, y, y, y]
    @test t .+ Dates.Year(2) == [Dates.Year(3), Dates.Year(3), Dates.Year(3), Dates.Year(3), Dates.Year(3)]

    let x = Dates.Year(5), y = Dates.Year(2)
        @test div(x, y) * y + rem(x, y) == x
        @test fld(x, y) * y + mod(x, y) == x
    end
end
@testset "Associativity" begin
    dt = Dates.DateTime(2012, 12, 21)
    test = ((((((((dt + y) - m) + w) - d) + h) - mi) + s) - ms)
    @test test == dt + y - m + w - d + h - mi + s - ms
    @test test == y - m + w - d + dt + h - mi + s - ms
    @test test == dt - m + y - d + w - mi + h - ms + s
    @test test == dt + (y - m + w - d + h - mi + s - ms)
    @test test == dt + y - m + w - d + (h - mi + s - ms)
    @test (dt + Dates.Year(4)) + Dates.Day(1) == dt + (Dates.Year(4) + Dates.Day(1))
    @test Dates.Date(2014, 1, 29) + Dates.Month(1) + Dates.Day(1) + Dates.Month(1) + Dates.Day(1) ==
        Dates.Date(2014, 1, 29) + Dates.Day(1) + Dates.Month(1) + Dates.Month(1) + Dates.Day(1)
    @test Dates.Date(2014, 1, 29) + Dates.Month(1) + Dates.Day(1) == Dates.Date(2014, 1, 29) + Dates.Day(1) + Dates.Month(1)
end
@testset "traits" begin
    @test Dates._units(Dates.Year(0)) == " years"
    @test Dates._units(Dates.Year(1)) == " year"
    @test Dates._units(Dates.Year(-1)) == " year"
    @test Dates._units(Dates.Year(2)) == " years"
    @test Dates.string(Dates.Year(0)) == "0 years"
    @test Dates.string(Dates.Year(1)) == "1 year"
    @test Dates.string(Dates.Year(-1)) == "-1 year"
    @test Dates.string(Dates.Year(2)) == "2 years"
    @test zero(Dates.Year) == Dates.Year(0)
    @test zero(Dates.Year(10)) == Dates.Year(0)
    @test zero(Dates.Month) == Dates.Month(0)
    @test zero(Dates.Month(10)) == Dates.Month(0)
    @test zero(Dates.Day) == Dates.Day(0)
    @test zero(Dates.Day(10)) == Dates.Day(0)
    @test zero(Dates.Hour) == Dates.Hour(0)
    @test zero(Dates.Hour(10)) == Dates.Hour(0)
    @test zero(Dates.Minute) == Dates.Minute(0)
    @test zero(Dates.Minute(10)) == Dates.Minute(0)
    @test zero(Dates.Second) == Dates.Second(0)
    @test zero(Dates.Second(10)) == Dates.Second(0)
    @test zero(Dates.Millisecond) == Dates.Millisecond(0)
    @test zero(Dates.Millisecond(10)) == Dates.Millisecond(0)
    @test Dates.Year(-1) < Dates.Year(1)
    @test !(Dates.Year(-1) > Dates.Year(1))
    @test Dates.Year(1) == Dates.Year(1)
    @test Dates.Year(1) != 1
    @test 1 != Dates.Year(1)
    @test Dates.Month(-1) < Dates.Month(1)
    @test !(Dates.Month(-1) > Dates.Month(1))
    @test Dates.Month(1) == Dates.Month(1)
    @test Dates.Day(-1) < Dates.Day(1)
    @test !(Dates.Day(-1) > Dates.Day(1))
    @test Dates.Day(1) == Dates.Day(1)
    @test Dates.Hour(-1) < Dates.Hour(1)
    @test !(Dates.Hour(-1) > Dates.Hour(1))
    @test Dates.Hour(1) == Dates.Hour(1)
    @test Dates.Minute(-1) < Dates.Minute(1)
    @test !(Dates.Minute(-1) > Dates.Minute(1))
    @test Dates.Minute(1) == Dates.Minute(1)
    @test Dates.Second(-1) < Dates.Second(1)
    @test !(Dates.Second(-1) > Dates.Second(1))
    @test Dates.Second(1) == Dates.Second(1)
    @test Dates.Millisecond(-1) < Dates.Millisecond(1)
    @test !(Dates.Millisecond(-1) > Dates.Millisecond(1))
    @test Dates.Millisecond(1) == Dates.Millisecond(1)
    @test_throws MethodError Dates.Year(1) < Dates.Millisecond(1)
    @test_throws MethodError Dates.Millisecond(1) < Dates.Year(1)

    # issue #27076
    @test Dates.Year(1) != Dates.Millisecond(1)
    @test Dates.Millisecond(1) != Dates.Year(1)
end

struct Beat <: Dates.Period
    value::Int64
end

Beat(p::Period) = Beat(Dates.toms(p) ÷ 86400)

@testset "comparisons with new subtypes of Period" begin
    # https://en.wikipedia.org/wiki/Swatch_Internet_Time
    Dates.value(b::Beat) = b.value
    Dates.toms(b::Beat) = Dates.value(b) * 86400
    Dates._units(b::Beat) = " beat" * (abs(Dates.value(b)) == 1 ? "" : "s")
    Base.promote_rule(::Type{Dates.Day}, ::Type{Beat}) = Dates.Millisecond
    Base.convert(::Type{T}, b::Beat) where {T<:Dates.Millisecond} = T(Dates.toms(b))

    @test Beat(1000) == Dates.Day(1)
    @test Beat(1) < Dates.Day(1)
end

@testset "basic properties" begin
    @test Dates.Year("1") == y
    @test Dates.Month("1") == m
    @test Dates.Week("1") == w
    @test Dates.Day("1") == d
    @test Dates.Hour("1") == h
    @test Dates.Minute("1") == mi
    @test Dates.Second("1") == s
    @test Dates.Millisecond("1") == ms
    @test Dates.Microsecond("1") == us
    @test Dates.Nanosecond("1") == ns
    @test_throws ArgumentError Dates.Year("1.0")
    @test Dates.Year(parse(Float64, "1.0")) == y

    dt = Dates.DateTime(2014)
    @test typeof(Dates.Year(dt)) <: Dates.Year
    @test typeof(Dates.Month(dt)) <: Dates.Month
    @test typeof(Dates.Week(dt)) <: Dates.Week
    @test typeof(Dates.Day(dt)) <: Dates.Day
    @test typeof(Dates.Hour(dt)) <: Dates.Hour
    @test typeof(Dates.Minute(dt)) <: Dates.Minute
    @test typeof(Dates.Second(dt)) <: Dates.Second
    @test typeof(Dates.Millisecond(dt)) <: Dates.Millisecond
end
@testset "Default values" begin
    @test Dates.default(Dates.Year) == y
    @test Dates.default(Dates.Month) == m
    @test Dates.default(Dates.Week) == w
    @test Dates.default(Dates.Day) == d
    @test Dates.default(Dates.Hour) == zero(Dates.Hour)
    @test Dates.default(Dates.Minute) == zero(Dates.Minute)
    @test Dates.default(Dates.Second) == zero(Dates.Second)
    @test Dates.default(Dates.Millisecond) == zero(Dates.Millisecond)
    @test Dates.default(Dates.Microsecond) == zero(Dates.Microsecond)
    @test Dates.default(Dates.Nanosecond) == zero(Dates.Nanosecond)
end
@testset "Conversions" begin
    @test Dates.toms(ms) == Dates.value(Dates.Millisecond(ms)) == 1
    @test Dates.toms(s)  == Dates.value(Dates.Millisecond(s)) == 1000
    @test Dates.toms(mi) == Dates.value(Dates.Millisecond(mi)) == 60000
    @test Dates.toms(h)  == Dates.value(Dates.Millisecond(h)) == 3600000
    @test Dates.toms(d)  == Dates.value(Dates.Millisecond(d)) == 86400000
    @test Dates.toms(w)  == Dates.value(Dates.Millisecond(w)) == 604800000

    @test Dates.days(ms) == Dates.days(s) == Dates.days(mi) == Dates.days(h) == 0
    @test Dates.days(Dates.Millisecond(86400000)) == 1
    @test Dates.days(Dates.Second(86400)) == 1
    @test Dates.days(Dates.Minute(1440)) == 1
    @test Dates.days(Dates.Hour(24)) == 1
    @test Dates.days(d) == 1
    @test Dates.days(w) == 7
end
@testset "issue #9214" begin
    @test 2s + (7ms + 1ms) == (2s + 7ms) + 1ms == 1ms + (2s + 7ms) == 1ms + (1s + 7ms) + 1s == 1ms + (2s + 3d + 7ms) + (-3d) == (1ms + (2s + 3d)) + (7ms - 3d) == (1ms + (2s + 3d)) - (3d - 7ms)
    @test 1ms - (2s + 7ms) == -((2s + 7ms) - 1ms) == (-6ms) - 2s
    @test emptyperiod == ((d + y) - y) - d == ((d + y) - d) - y
    @test emptyperiod == 2y + (m - d) + ms - ((m - d) + 2y + ms)
    @test emptyperiod == 0ms
    @test string(emptyperiod) == "empty period"
    @test string(ms + mi + d + m + y + w + h + s + 2y + m) == "3 years, 2 months, 1 week, 1 day, 1 hour, 1 minute, 1 second, 1 millisecond"
    @test 8d - s == 1w + 23h + 59mi + 59s
    @test h + 3mi == 63mi
    @test y - m == 11m
end
@testset "compound periods and types" begin
    # compound periods should avoid automatically converting period types
    @test (d - h).periods == Dates.Period[d, -h]
    @test d - h == 23h
    @test !isequal(d - h, 23h)
    @test isequal(d - h, 2d - 2h - 1d + 1h)
    @test sprint(show, y + m) == string(y + m)
    @test convert(Dates.CompoundPeriod, y) + m == y + m
end
@testset "compound period simplification" begin
    # reduce compound periods into the most basic form
    @test Dates.canonicalize(h - mi).periods == Dates.Period[59mi]
    @test Dates.canonicalize(-h + mi).periods == Dates.Period[-59mi]
    @test Dates.canonicalize(-y + d).periods == Dates.Period[-y, d]
    @test Dates.canonicalize(-y + m - w + d).periods == Dates.Period[-11m, -6d]
    @test Dates.canonicalize(-y + m - w + ms).periods == Dates.Period[-11m, -6d, -23h, -59mi, -59s, -999ms]
    @test Dates.canonicalize(y - m + w - d + h - mi + s - ms).periods == Dates.Period[11m, 6d, 59mi, 999ms]
    @test Dates.canonicalize(-y + m - w + d - h + mi - s + ms).periods == Dates.Period[-11m, -6d, -59mi, -999ms]

    @test Dates.Date(2009, 2, 1) - (Dates.Month(1) + Dates.Day(1)) == Dates.Date(2008, 12, 31)
    @test_throws MethodError (Dates.Month(1) + Dates.Day(1)) - Dates.Date(2009,2,1)
end
@testset "unary ops and vectorized period arithmetic" begin
    pa = [1y 1m 1w 1d; 1h 1mi 1s 1ms]
    cpa = [1y + 1s 1m + 1s 1w + 1s 1d + 1s; 1h + 1s 1mi + 1s 2m + 1s 1s + 1ms]

    @test +pa == pa == -(-pa)
    @test -pa == map(-, pa)
    @test 1y .+ pa == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test (1y + 1m) .+ pa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]
    @test pa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test pa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]

    @test 1y .+ cpa == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test (1y + 1m) .+ cpa == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]
    @test cpa .+ 1y == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test cpa .+ (1y + 1m) == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]

    @test 1y .+ pa == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test (1y + 1m) .+ pa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]
    @test pa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test pa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]

    @test 1y .+ cpa == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test (1y + 1m) .+ cpa == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]
    @test cpa .+ 1y == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test cpa .+ (1y + 1m) == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]

    @test 1y .- pa == [0y 1y-1m 1y-1w 1y-1d; 1y-1h 1y-1mi 1y-1s 1y-1ms]
    @test (1y + 1m) .- pa == [1m 1y 1y + 1m-1w 1y + 1m-1d; 1y + 1m-1h 1y + 1m-1mi 1y + 1m-1s 1y + 1m-1ms]
    @test pa .- (1y + 1m) == [-1m -1y -1y-1m + 1w -1y-1m + 1d; -1y-1m + 1h -1y-1m + 1mi -1y-1m + 1s -1y-1m + 1ms]
    @test pa .- 1y == [0y 1m-1y -1y + 1w -1y + 1d; -1y + 1h -1y + 1mi -1y + 1s -1y + 1ms]

    @test 1y .- cpa == [-1s 1y-1m-1s 1y-1w-1s 1y-1d-1s; 1y-1h-1s 1y-1mi-1s 1y-2m-1s 1y-1ms-1s]
    @test (1y + 1m) .- cpa == [1m-1s 1y-1s 1y + 1m-1w-1s 1y + 1m-1d-1s; 1y + 1m-1h-1s 1y + 1m-1mi-1s 1y-1m-1s 1y + 1m-1s-1ms]
    @test cpa .- 1y == [1s -1y + 1m + 1s -1y + 1w + 1s -1y + 1d + 1s; -1y + 1h + 1s -1y + 1mi + 1s -1y + 2m + 1s -1y + 1ms + 1s]
    @test cpa .- (1y + 1m) == [-1m + 1s -1y + 1s -1y-1m + 1w + 1s -1y-1m + 1d + 1s; -1y-1m + 1h + 1s -1y-1m + 1mi + 1s -1y + 1m + 1s -1y + -1m + 1s + 1ms]

    @test [1y 1m; 1w 1d] + [1h 1mi; 1s 1ms] == [1y + 1h 1m + 1mi; 1w + 1s 1d + 1ms]
    @test [1y 1m; 1w 1d] - [1h 1mi; 1s 1ms] == [1y-1h 1m-1mi; 1w-1s 1d-1ms]
    @test [1y 1m; 1w 1d] - [1h 1mi; 1s 1ms] - [1y-1h 1m-1mi; 1w-1s 1d-1ms] == [emptyperiod emptyperiod; emptyperiod emptyperiod]

    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] + [1h 1mi; 1s 1ms] == [1y + 1h + 1s 1m + 1mi + 1s; 1w + 2s 1d + 1s + 1ms]
    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] - [1h 1mi; 1s 1ms] == [1y-1h + 1s 1m-1mi + 1s; 1w 1d + 1s-1ms]

    @test [1y 1m; 1w 1d] + [1h + 1s 1mi + 1s; 1m + 1s 1s + 1ms] == [1y + 1h + 1s 1m + 1mi + 1s; 1w + 1m + 1s 1d + 1s + 1ms]
    @test [1y 1m; 1w 1d] - [1h + 1s 1mi + 1s; 1m + 1s 1s + 1ms] == [1y-1h-1s 1m-1mi-1s; 1w-1m-1s 1d-1s-1ms]

    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] + [1y + 1h 1y + 1mi; 1y + 1s 1y + 1ms] == [2y + 1h + 1s 1y + 1m + 1mi + 1s; 1y + 1w + 2s 1y + 1d + 1s + 1ms]
    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] - [1y + 1h 1y + 1mi; 1y + 1s 1y + 1ms] == [1s-1h 1m + 1s-1y-1mi; 1w-1y 1d + 1s-1y-1ms]
end
@testset "Equality and hashing between FixedPeriod types" begin
    let types = (Dates.Week, Dates.Day, Dates.Hour, Dates.Minute,
                 Dates.Second, Dates.Millisecond, Dates.Microsecond, Dates.Nanosecond)
        for i in 1:length(types), j in i:length(types), x in (0, 1, 235, -4677, 15250)
            local T, U, y, z
            T = types[i]
            U = types[j]
            y = T(x)
            z = convert(U, y)
            @test y == z
            @test hash(y) == hash(z)
        end
    end
end
@testset "Equality and hashing between OtherPeriod types" begin
    for x in (0, 1, 235, -4677, 15250)
        local x, y, z
        y = Dates.Year(x)
        z = convert(Dates.Month, y)
        @test y == z
        @test hash(y) == hash(z)
    end
end

@testset "#30832" begin
    @test Dates.toms(Dates.Second(1) + Dates.Nanosecond(1)) == 1e3
    @test Dates.tons(Dates.Second(1) + Dates.Nanosecond(1)) == 1e9 + 1
    @test Dates.toms(Dates.Second(1) + Dates.Microsecond(1)) == 1e3
end

end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
"""
    svm(dataset::Dataset, labels::Labels; C::Real = 10.0)

perform suppor vector machines learning on `dataset` and `labels`

returns optimal vector alpha which defines which samples are suppor vectors

- keyword argument `C` defines penalty for misclassified points in soft-margin SVM
"""
function svm(dataset::Dataset, labels::Labels; C::Real = 10.0)
    X = dataset()
    y = labels()

    alpha = dual(X, y, C)

    return alpha
end

"""
    linearkernel(xi::Vector{<:Real}, xj::Vector{<:Real})

simplest svm kernel using a simple dot product
"""
function linearkernel(xi::Vector{<:Real}, xj::Vector{<:Real})
    return dot(xi, xj)
end

"""
    polynomialkernel(xi::Vector{<:Real}, xj::Vector{<:Real}; degree::Integer = 1)

polynomial kernel with variable `degree` of the polynomial
"""
function polynomialkernel(xi::Vector{<:Real}, xj::Vector{<:Real}; degree::Integer = 1)
    return (1 + dot(xi, xj))^degree
end

"""
    gaussiankernel(xi::Vector{<:Real}, xj::Vector{<:Real}; variance::Real = 1.0)

kernel using normal distribution and can be parametrized by `variance`
"""
function gaussiankernel(xi::Vector{<:Real}, xj::Vector{<:Real}; variance::Real = 1.0)
    return exp(-dot(xi - xj, xi - xj) / (2 * variance))
end

"""
    kernel(xi::Vector{<:Real}, xj::Vector{<:Real})

kernel wrapper to select from prepared kernels without the need to alter more complicated functions
"""
function kernel(xi::Vector{<:Real}, xj::Vector{<:Real})
    return linearkernel(xi, xj)
    # return guassiankernel(xi, xj)
    # return polynomialkernel(xi, xj; degree=2)
end

"""
    dual(X::Matrix{<:Real}, y::Vector{<:Integer}, C::Real)

finds optimal alpha by solving the dual formulation using quadratic programming solver
"""
function dual(X::Matrix{<:Real}, y::Vector{<:Integer}, C::Real)
    _, n = size(X)

    qp = Model(OSQP.Optimizer)
    set_silent(qp)

    @variable(qp, alpha[1:n])
    @objective(qp, Max, sum(alpha) - sum([alpha[i] * alpha[j] * y[i] * y[j] * kernel(X[:, i], X[:, j]) for i in 1:n for j in 1:n]) / 2)
    @constraint(qp, dot(alpha, y) == 0)
    @constraint(qp, 0 .<= alpha .<= C)

    optimize!(qp)

    return value.(alpha)
end

"""
    classify(x::Vector{<:Real}, traindata::Dataset, trainlabels::Labels, alpha::Vector{<:Real})

classify one data sample `x` with given `alpha` weights
"""
function classify(x::Vector{<:Real}, traindata::Dataset, trainlabels::Labels, alpha::Vector{<:Real})
    X = traindata()
    y = trainlabels()
    n = length(y)
    f = sum([alpha[i] * y[i] * kernel(X[:, i], x) for i in 1:n])
    return f >= 0 ? 1 : -1
end

"""
    classify(test::Dataset, traindata::Dataset, trainlabels::Labels, alpha::Vector{<:Real})

classify the whole dataset `test` with given `alpha` weights and `traindata`, `trainlabels`
"""
function classify(test::Dataset, traindata::Dataset, trainlabels::Labels, alpha::Vector{<:Real})
    X = test()
    _, n = size(X)
    return Labels(test, [classify(X[:, i], traindata, trainlabels, alpha) for i in 1:n])
end


text:
using Random
import Distributions: ValueSupport, Sampleable
import Base: eltype, rand
import Random: rand!
using Phylo
using Distributions
using Missings
using IterableTables: getiterator
using Unitful
import Unitful: numtype
using Statistics

@inline numtype(::Type{R}) where R <: Real = R
iscontinuous(::Type{N}) where N <: Number = numtype(N) <: AbstractFloat

struct Phylogenetics{T <: AbstractTree} <: ValueSupport end
Base.eltype(::Type{Phylogenetics{T}}) where T <: AbstractTree = T
Base.eltype(::Sampleable{S, P}) where {S, P <: Phylogenetics} = eltype(P)

Random.rand!(s::Sampleable, t::AbstractTree) = rand!(Random.GLOBAL_RNG, s, t)

"""
    Nonultrametric{T <: AbstractTree,
                   SAMP <: Sampleable}(n::Int,
                                       sampleable::SAMP = Exponential())
    Nonultrametric{T <: AbstractTree,
                   SAMP <: Sampleable}(tiplabels::Vector{String},
                                       sampleable::SAMP = Exponential())

The sampler for non-ultrametric phylogenetic trees of size `n` or with
tip labels `tiplabels`. Generate random trees by calling rand().
"""
struct Nonultrametric{T <: AbstractTree,
                      SAMP <: Sampleable{Univariate, Continuous},
                      LenUnits <: Number} <:
                          Sampleable{Univariate, Phylogenetics{T}}
    n::Int
    tiplabels::Vector{String}
    sampleable::SAMP
    leafinfo::Any
    height::LenUnits

    function Nonultrametric{T, SAMP,
                            LenUnits}(n::Int, tiplabels::Vector{String},
                                      sampleable::SAMP, leafinfo,
                                      height::LenUnits) where {T, SAMP,
                                                               LenUnits}
        return new{T, SAMP, LenUnits}(n, tiplabels, sampler(sampleable),
                                      leafinfo, height)
    end

    function Nonultrametric{T, SAMP,
                            LenUnits}(n::Int, sampleable::SAMP,
                                      height::LenUnits) where {T, SAMP,
                                                               LenUnits}
        return new{T, SAMP, LenUnits}(n, ["tip $i" for i in Base.OneTo(n)],
                                      sampler(sampleable), missing, height)
    end

    function Nonultrametric{T, SAMP,
                            LenUnits}(tiplabels::Vector{String},
                                      sampleable::SAMP,
                                      height::LenUnits) where {T, SAMP,
                                                               LenUnits}
        return new{T, SAMP, LenUnits}(length(tiplabels), tiplabels,
                                      sampler(sampleable), missing, height)
    end
end

function Nonultrametric{T}(n::Int, height = 1.0) where T <: AbstractTree
    return Nonultrametric{T, Exponential, typeof(height)}(n, Exponential(),
                                                          height)
end

function Nonultrametric{T}(tiplabels::Vector{String},
                           height = 1.0) where T <: AbstractTree
    return Nonultrametric{T, Exponential, typeof(height)}(tiplabels,
                                                          Exponential(),
                                                          height)
end

function Nonultrametric{T}(leafinfo,
                           height = 1.0) where T <: AbstractTree
    tipnames = unique([info[1] for info in getiterator(leafinfo)])
    return Nonultrametric{T, Exponential,
                          typeof(height)}(length(tipnames), tipnames,
                                          Exponential(), leafinfo, height)
end

Nonultrametric(info::LI, height::U = 1.0) where {LI, U} =
    Nonultrametric{Phylo.LT{OneRoot, LI, U}}(info)

Nonultrametric(n::Int, height::U = 1.0) where U =
    Nonultrametric{Phylo.LTD{OneRoot, U}}(n, height)
Nonultrametric(tiplabels::Vector{String}, height::U = 1.0) where U =
    Nonultrametric{Phylo.LTD{OneRoot, U}}(tiplabels, height)

function rand(rng::AbstractRNG, t::Nonultrametric{T, SAMP, U}) where
    {T, SAMP, U}
    t.n >= 2 || error("A tree must have at least 2 tips")
    if ismissing(t.leafinfo)
        tree = T(t.tiplabels)
    else
        tree = T(t.leafinfo)
    end
    while nroots(tree) > 1
        roots = getroots(tree)
        children = sample(rng, collect(roots), 2, replace=false)
        parent = createnode!(tree)
        createbranch!(tree, parent, children[1],
                      rand(rng, t.sampleable) * t.height)
        createbranch!(tree, parent, children[2],
                      rand(rng, t.sampleable) * t.height)
    end
    scale = t.height / mean(getheight.(tree, getleaves(tree)))
    for branch in collect(getbranches(tree))
        source = src(tree, branch)
        destination = dst(tree, branch)
        len = uconvert(unit(t.height), getlength(tree, branch) * scale)
        deletebranch!(tree, branch)
        createbranch!(tree, source, destination, len)
    end
    return tree
end

"""
    Ultrametric{T <: AbstractTree,
                SAMP <: Sampleable,
                LenUnits <: Number}(n::Int,
                                    sampleable::SAMP = Exponential())
    Ultrametric{T <: AbstractTree,
                SAMP <: Sampleable,
                LenUnits <: Number}(tiplabels::Vector{String},
                                    sampleable::SAMP = Exponential())

The sampler for ultrametric phylogenetic trees of size `n` or with
tip labels `tiplabels`. Generate random trees by calling rand().
"""
struct Ultrametric{T <: AbstractTree,
                   SAMP <: Sampleable{Univariate, Continuous},
                   LenUnits <: Number} <:
                       Sampleable{Univariate, Phylogenetics{T}}
    n::Int
    tiplabels::Vector{String}
    sampleable::SAMP
    leafinfo::Any
    height::LenUnits

    function Ultrametric{T, SAMP, U}(n::Int, tiplabels::Vector{String},
                                     sampleable::SAMP, leafinfo,
                                     height::U = 1.0) where {T, SAMP, U}
        return new{T, SAMP, U}(n, tiplabels, sampler(sampleable),
                               leafinfo, height)
    end

    function Ultrametric{T, SAMP, U}(n::Int, sampleable::SAMP,
                                     height::U = 1.0) where {T, SAMP, U}
        return new{T, SAMP, U}(n, ["tip $i" for i in Base.OneTo(n)],
                               sampler(sampleable), missing, height)
    end

    function Ultrametric{T, SAMP, U}(tiplabels::Vector{String},
                                     sampleable::SAMP,
                                     height::U = 1.0) where
        {T, SAMP, U}
        return new{T, SAMP, U}(length(tiplabels), tiplabels,
                               sampler(sampleable), missing, height)
    end
end

function Ultrametric{T}(n::Int, height = 1.0) where T <: AbstractTree
    return Ultrametric{T, Exponential, typeof(height)}(n, Exponential(), height)
end

function Ultrametric{T}(tiplabels::Vector{String},
                        height = 1.0) where T <: AbstractTree
    return Ultrametric{T, Exponential, typeof(height)}(tiplabels, Exponential(),
                                                       height)
end

function Ultrametric{T}(leafinfo, height = 1.0) where T <: AbstractTree
    tipnames = unique([info[1] for info in getiterator(leafinfo)])
    return Ultrametric{T, Exponential, typeof(height)}(length(tipnames),
                                                       tipnames,
                                                       Exponential(), leafinfo,
                                                       height)
end

Ultrametric(info::LI, height::U = 1.0) where {LI, U} =
    Ultrametric{Phylo.LT{OneRoot, LI, U}}(info)
Ultrametric(n::Int, height::U = 1.0) where U =
    Ultrametric{Phylo.LTD{OneRoot, U}}(n, height)
Ultrametric(tiplabels::Vector{String}, height::U = 1.0) where U =
    Ultrametric{Phylo.LTD{OneRoot, U}}(tiplabels, height)

function rand(rng::AbstractRNG, t::Ultrametric{T, SAMP}) where {T, SAMP}
    t.n >= 2 || error("A tree must have at least 2 tips")
    if ismissing(t.leafinfo)
        tree = T(t.tiplabels)
    else
        tree = T(t.leafinfo)
    end
    depth = zero(rand(rng, t.sampleable)) * t.height
    leaves = getleaves(tree)
    while nroots(tree) > 1
        roots = getroots(tree)
        tocoalesce = collect(roots)
        coalescers = sample(rng, tocoalesce, 2, replace=false)
        parent = createnode!(tree)
        depth += rand(rng, t.sampleable) * 2.0 / length(tocoalesce) * t.height
        d1 = getheight(tree, first(nodefuture(tree, coalescers[1]) ∩
                                   leaves))
        d2 = getheight(tree, first(nodefuture(tree, coalescers[2]) ∩
                                   leaves))
        createbranch!(tree, parent, coalescers[1], depth - d1)
        createbranch!(tree, parent, coalescers[2], depth - d2)
    end
    scale = t.height / getheight(tree, first(getleaves(tree)))
    for branch in collect(getbranches(tree))
        source = src(tree, branch)
        destination = dst(tree, branch)
        len = uconvert(unit(t.height), getlength(tree, branch) * scale)
        deletebranch!(tree, branch)
        createbranch!(tree, source, destination, len)
    end
    return tree
end

rand(s::S, treenames::AbstractVector{LABEL}) where
{LABEL, TREE <: AbstractTree,
 S <: Sampleable{Univariate, Phylogenetics{TREE}}} =
     rand(Random.GLOBAL_RNG, s, treenames)

function rand(rng::AbstractRNG, s::S, treenames::AbstractVector{LABEL}) where
    {LABEL, TREE <: AbstractTree{OneTree},
     S <: Sampleable{Univariate, Phylogenetics{TREE}}}
    trees = Dict{eltype(treenames), TREE}()
    for name in treenames
        trees[name] = rand(rng, s)
    end
    return TreeSet(trees)
end

rand(rng::AbstractRNG, s::S, n::Tuple{Int}) where
    {TREE <: AbstractTree{OneTree},
     S <: Sampleable{Univariate, Phylogenetics{TREE}}} =
    rand(rng, s, 1:n[1])

abstract type EvolvedTrait{T <: AbstractTree} <: ValueSupport end
Base.eltype(::Type{EvolvedTrait{T}}) where T <: AbstractTree = T
Base.eltype(::Sampleable{S, P}) where {S, P <: EvolvedTrait} = eltype(P)

"""
    BrownianTrait{T <: AbstractTree, N <: Number}

A continuous trait evolved on a phylogenetic tree. This is a
`Sampleable` type, so a random trait can be created using
`rand()`. The trait to be evolved can be any continuous numeric type,
including `Unitful` types for instance, and in the simplest case is
determined by the third argument to the constructor `start`:

function BrownianTrait(tree::AbstractTree, trait::String, start::Number = 0.0;
                       σ² = missing, σ = missing, f::Function = identity)

Note that when `Unitful` is being used, either here or in branch
lengths, `σ`/`σ²` keyword argument units must be appropriate. The final keyword argument, `f`, is a function to transform the evolved gaussian trait into its true value. By default this is the identity function, but can, for instance, be `abs` to force a positive value on the trait, or more complex functions as required, such as a transformation to turn a continuous variable into a discrete trait
"""
struct BrownianTrait{T <: AbstractTree, N <: Number} <:
    Sampleable{Univariate, EvolvedTrait{T}}
    tree::T
    trait::String
    start::N
    _rand::Function
    f::Function
end

function BrownianTrait(tree::T, trait::String, start::N = 0.0;
                       σ² = missing, σ = missing, f::Function = identity) where
    {T <: AbstractTree, N <: Number}
    iscontinuous(N) ||
        throw(TypeError("Type $N must be continuous for a Gaussian Trait"))
    if ismissing(σ)
        σ = sqrt(σ²)
    end

    f ≢ identity && f(start) ≉ start &&
        @warn "Note that the third argument (the starting state) for trait '$trait' is untransformed - transformed version is $(f(start))"

    if ismissing(σ)
        dimension(N) ≡
            dimension(sqrt(getlength(tree, first(getbranches(tree))))) ||
            throw(DimensionMismatch("Dimensions of start, σ[²] and branch lengths must combine correctly if using Unitful"))
        return BrownianTrait{T, N}(tree, trait, start,
                                   ((rng::AbstractRNG, start::N, length) ->
                                    start + randn(rng, typeof(one(start))) *
                                    sqrt(length)), f)
    else
        dimension(N) ≡
            dimension(σ*sqrt(getlength(tree,
                                       first(getbranches(tree))))) ||
                                           throw(DimensionMismatch("Dimensions of start, σ[²] and branch lengths must combine correctly if using Unitful"))
        return BrownianTrait{T, N}(tree, trait, start,
                                   ((rng::AbstractRNG, start::N, length) ->
                                    start +
                                    σ * randn(rng, typeof(one(start))) *
                                    sqrt(length)), f)
    end
end

function rand!(rng::AbstractRNG,
               bm::BrownianTrait{TREE, N},
               tree::TREE) where {TREE <: AbstractTree, N <: Number}
    trait = Dict{nodetype(TREE), N}()
    use_dict = (bm.f ≢ identity)
    for node in traversal(tree, preorder)
        if isroot(tree, node)
            if use_dict
                trait[node] = bm.start
                setnodedata!(tree, node, bm.trait, bm.f(bm.start))
            else
                setnodedata!(tree, node, bm.trait, bm.start)
            end
        else
            inb = getinbound(tree, node)
            prt = src(tree, inb)
            previous = use_dict ? trait[prt] :
                getnodedata(tree, prt, bm.trait)
            value = bm._rand(rng, previous, N(getlength(tree, inb)))
            if use_dict
                trait[node] = value
                setnodedata!(tree, node, bm.trait, bm.f(value))
            else
                setnodedata!(tree, node, bm.trait, value)
            end
        end
    end
    return tree
end

function rand(rng::AbstractRNG,
              bm::BrownianTrait{TREE, N}) where {TREE <: AbstractTree,
                                                 N <: Number}
    untrait = Dict{nodetype(TREE), N}()
    traitbyname = Dict{nodenametype(TREE), typeof(bm.f(bm.start))}()
    for node in traversal(bm.tree, preorder)
        if isroot(bm.tree, node)
            untrait[node] = bm.start
            traitbyname[getnodename(bm.tree, node)] = bm.f(bm.start)
        else
            inb = getinbound(bm.tree, node)
            prt = src(bm.tree, inb)
            previous = untrait[prt]
            value = bm._rand(rng, previous, N(getlength(bm.tree, inb)))
            untrait[node] = value
            traitbyname[getnodename(bm.tree, node)] = bm.f(value)
        end
    end
    return traitbyname
end

abstract type AbstractDiscreteTrait{T <: AbstractTree, E <: Enum} <:
    Sampleable{Univariate, EvolvedTrait{T}} end

"""
    DiscreteTrait{T <: AbstractTree, E <: Enum}

A discrete trait evolved on a phylogenetic tree. This is a
`Sampleable` type, so a random trait can be created using
`rand(dt)`. The trait to be evolved must be an Enum (generally created
using `@enum`), and is the second argument to the constructor:

function DiscreteTrait(tree::AbstractTree, ttype::Type{<:Enum},
                       transition_matrix::AbstractMatrix{Float64},
                       trait::String = "\$ttype")

The transition matrix holds transition rates from row to column (so
row sums must be zero), and the transition probabilities in a branch
are calculated as `exp(transition_matrix .* branch_length)`.
"""
struct DiscreteTrait{T <: AbstractTree, E <: Enum, N <: Number} <:
    Sampleable{Univariate, EvolvedTrait{T}}
    tree::T
    transition_matrix::Matrix{N}
    trait::String
end

function DiscreteTrait(tree::TREE, ttype::Type{TRAIT},
                       transition_matrix::AbstractMatrix{N},
                       trait::String = "$ttype") where
    {TREE <: AbstractTree, TRAIT <: Enum, N <: Number}
    all(size(transition_matrix) .== length(instances(TRAIT))) ||
        error("Transition matrix size must match number of traits in $TRAIT")
    dimension(first(transition_matrix)) ==
        NoDims / dimension(getlength(tree, first(getbranches(tree)))) ||
        throw(DimensionMismatch("Dimensions of transition_matrix and branch lengths must cancel if using Unitful units"))
    dimension(getlength(tree, first(getbranches(tree))) *
              first(transition_matrix)) == NoDims ||
              throw(DimensionMismatch("Dimensions of transition_matrix and branch lengths must cancel correctly if using Unitful"))
    all(sum(transition_matrix .* getlength(tree, first(getbranches(tree))),
            dims=2) .+ 1 .≈ 1) ||
                error("Row sums of transition matrix must be zero")

    return DiscreteTrait{TREE, TRAIT, N}(tree, transition_matrix, trait)
end

function enum_rand(rng::AbstractRNG, current::TRAIT,
                   mat::AbstractMatrix{Float64}) where TRAIT <: Enum
    traits = instances(TRAIT)
    row = findfirst(x -> x == current, traits)
    p = @view mat[row, :]
    tot = length(traits)

    rp = 1.0  # remaining total probability
    i = 1
    while i < tot
        @inbounds pi = p[i]
        if pi < rp
            if rand(rng) < pi / rp
                return traits[i]
            end
            rp -= pi
        else
            return traits[i]
        end
        i += 1
    end

    return traits[tot]
end

function rand!(rng::AbstractRNG, dt::DiscreteTrait{TREE, TRAIT},
               tree::TREE) where {TREE <: AbstractTree, TRAIT <: Enum}
    traits = instances(TRAIT)
    num = length(traits)
    for node in traversal(tree, preorder)
        if isroot(tree, node)
            setnodedata!(tree, node, dt.trait,
                         sample(rng, collect(traits)))
        else
            inb = getinbound(tree, node)
            prt = src(tree, inb)
            previous = getnodedata(tree, prt, dt.trait)
            value = enum_rand(rng, previous,
                              exp(uconvert.(NoUnits,
                                            getlength(tree, inb) .*
                                            dt.transition_matrix)))
            setnodedata!(tree, node, dt.trait, value)
        end
    end
    return tree
end

"""
    SymmetricDiscreteTrait{T <: AbstractTree, E <: Enum}

The simplest possible discrete trait evolved on a phylogenetic
tree. This is a `Sampleable` type, so a random trait can be created
using `rand(sdt)`. The trait to be evolved must be an Enum (generally
created using `@enum`), and is the second argument to the constructor:

function DiscreteTrait(tree::AbstractTree, ttype::Type{<:Enum},
                       transition_rate::Number,
                       trait::String = "\$ttype")

The transition matrix holds transition rates from row to column (so
row sums must be zero), and the transition probabilities in a branch
are calculated as `exp(transition_matrix .* branch_length)`.
"""
struct SymmetricDiscreteTrait{T <: AbstractTree, E <: Enum, N <: Number} <:
    Sampleable{Univariate, EvolvedTrait{T}}
    tree::T
    transition_rate::N
    trait::String
end

function SymmetricDiscreteTrait(tree::TREE, ttype::Type{TRAIT},
                                transition_rate::N,
                                trait::String = "$ttype") where
    {TREE <: AbstractTree, TRAIT <: Enum, N <: Number}
    dimension(getlength(tree, first(getbranches(tree))) *
              transition_rate) == NoDims ||
              throw(DimensionMismatch("Dimensions of transition_matrix and branch lengths must cancel correctly if using Unitful"))

    return SymmetricDiscreteTrait{TREE, TRAIT, N}(tree, transition_rate, trait)
end

function enum_rand(rng::AbstractRNG, current::TRAIT,
                   p_stay::Float64) where TRAIT <: Enum
    if rand(rng) < p_stay
        return current
    end

    traits = instances(TRAIT)
    num = length(traits)
    remain = num - 1  # remaining total options
    from = findfirst(x -> x == current, traits)
    passed = false
    i = 1
    while i < num
        passed |= (from == i)
        if remain > 1
            if rand(rng) < inv(remain)
                return traits[i+passed]
            end
            i += 1
            remain -= 1
        else
            return traits[i+passed]
        end
    end
end

function rand!(rng::AbstractRNG, dt::SymmetricDiscreteTrait{TREE, TRAIT},
               tree::TREE) where {TREE <: AbstractTree, TRAIT <: Enum}
    traits = instances(TRAIT)
    num = length(traits)
    frac = inv(num)
    for node in traversal(tree, preorder)
        if isroot(tree, node)
            setnodedata!(tree, node, dt.trait,
                         sample(rng, collect(instances(TRAIT))))
        else
            inb = getinbound(tree, node)
            prt = src(tree, inb)
            previous = getnodedata(tree, prt, dt.trait)
            value = enum_rand(rng, previous,
                              frac + (1.0 - frac) *
                              exp(-num * getlength(tree, inb) *
                                  dt.transition_rate))
            setnodedata!(tree, node, dt.trait, value)
        end
    end
    return tree
end


text:
"""
A parametric node of phylogenetic tree
"""
struct Node
  in::Vector{Int64}
  out::Vector{Int64}

  function Node()
    new(Int64[], Int64[])
  end
end

"""
A directed parametric branch connecting two nodes of phylogenetic tree
"""
struct Branch
  source::Int64
  target::Int64
  length::Float64

  function Branch(source::Int64, target::Int64, length::Float64)
    if length < 0.
      error("Branch length is < 0.0")
    end
    new(source, target, length)
  end

  function Branch(source::Int64, target::Int64)
    new(source, target, NaN)
  end
end

"""
Parametric phylogenetic tree object
"""
struct Tree
  nodes::Dict{Int64, Node}
  branches::Dict{Int64, Branch}
  height::Float64

  Tree() = new(Dict{Int64, Node}(), Dict{Int64, Branch}(), 0.0)
  Tree(x::Float64) = new(Dict{Int64, Node}(), Dict{Int64, Branch}(), x)
end


text:
export saveactuatorconfig, saveexperimentmatrix

using HDF5

"""
`saveactuatorconfig(h5, dev::AbstractActuator)`

Saves an actuator configuration to a HDF5 file.
"""
function saveactuatorconfig(h5, dev::AbstractActuator)

    dtype = string(typeof(dev))
    dname = devname(dev)
    h5[dname] = [string(a) for a in axesnames(dev)]
    attributes(h5[dname])["devtype"] = dtype
    return
end

function saveactuatorconfig(h5, devs::AbstractVector{<:AbstractActuator})
    for actuator in devs
        saveactuatorconfig(h5, actuator)
    end
    return
end

"""
`saveactuatorconfig(h5, devs::ActuatorSet)`

Saves an [`ActuatorSet`](@ref) configuration to a HDF5 file. It basically
saves the configuration of each individual actuator in the set.
"""
function saveactuatorconfig(h5, devs::ActuatorSet)

    dtype = string(typeof(devs))
    dname = devname(devs)
    axes = axesnames(devs)
    g = create_group(h5, dname)
    attributes(g)["axes"] = axes
    attributes(g)["type"] = dtype
    for actuator in devs.actuators
        saveactuatorconfig(g, actuator)
    end

    return

end

"""
`saveexperimentmatrix(h5, pts::ExperimentMatrix, name="points")`

Saves the experimental matrix points to an HDF5 file.
"""
function saveexperimentmatrix(h5, pts::ExperimentMatrix, name="points")

    params = matrixparams(pts)
    dtype = string(typeof(pts))
    nparams = length(params)
    g = create_group(h5, name)

    g["points"] = pts.pts
    attributes(g)["type"] = dtype
    attributes(g)["params"] = params
    attributes(g)["nparams"] = nparams

    return
end

"""
`saveexperimentmatrix(h5, pts::CartesianExperimentMatrix, name="points")`

Saves the `CartesianExperimentMatrix` points to an HDF5 file. This specific
method is useful when the points is the cartesian product of several axis
and you want to know each axis besides the coordinates of every point.
"""
function saveexperimentmatrix(h5, pts::CartesianExperimentMatrix, name="points")

    dtype = "CartesianExperimentMatrix"
    params = matrixparams(pts)
    nparams = length(params)

    g = create_group(h5, name)

    g["points"] = pts.pts

    attributes(g)["params"] = params
    attributes(g)["nparams"] = params
    attributes(g)["type"] = dtype

    for i in 1:nparams
        nn = numstring(i, 4)
        g[nn] = pts.axes[i]
        attributes(g[nn])["param"] = params[i]
    end

    return

end

function saveexperimentmatrix(h5, pts::ExperimentMatrixProduct, name="points")

    dtype = "ExperimentMatrixProduct"
    params = matrixparams(pts)
    nparams = length(params)

    g = create_group(h5, name)

    g["points"] = experimentpoints(pts)
    g["params"] = params
    g["nparams"] = nparams

    attributes(g)["params"] = params
    attributes(g)["nparams"] = params
    attributes(g)["type"] = dtype
    npoints = length(pts.points)
    attributes(g)["npoints"] = npoints
    for i in 1:npoints
        nn = numstring(i, 4)
        saveexperimentmatrix(g, pts.points[i], nn)
    end

    return

end


text:
# struct TemperedModel <: AbstractPPL.AbstractProbabilisticProgram
#     model :: DynamicPPL.Model
#     β     :: AbstractFloat
# end

function make_tempered_model end


text:
# [ref](https://opensourc.es/blog/random-forest)

mutable struct Node

    feature    :: Int64
    comp_value :: Float64
    data_idxs  :: Vector{Int64}
    mean       :: Float64
    left       :: Union{Nothing, Node}
    right      :: Union{Nothing, Node}

    Node() = new()

end

mutable struct Tree
    root   :: Node
    Tree() = new()
end

function create_random_forest(feature_matrix, train_ys, ntrees)

    forest             = Vector{Tree}()
    tc                 =  1
    total_nfeatures    = size(feature_matrix)[1]
    nfeatures_per_tree = 5
    ncols_per_tree     = length(train_ys)÷2

    while tc < ntrees

        feature_idx = sample(1:total_nfeatures,
                             nfeatures_per_tree;
                             replace=false)

        cols = sample(1:length(train_ys), ncols_per_tree; replace=false)

        tree = create_random_tree(feature_matrix,
                                  feature_idx,
                                  cols,
                                  train_ys,
                                  total_nfeatures)

        if tree.root.left_child != nothing
            push!(forest,tree)
            tc += 1
        else
            println("Couldn't create tree!!")
        end
        flush(stdout)
    end

    forest

end

function create_random_tree(glob_feature_matrix, ...)

    tree = Tree()
    @views feature_matrix = glob_feature_matrix[feature_idx, cols]
    root, queue = create_root_node!(...)

    while length(queue) > 0
        node = popfirst!(queue)
        left_node, right_node = compute_node!(...)
        if left_node != nothing
            queue_compute_nodes!(...)
        end
    end
    tree.root = root
    tree

end


text:
using Random

struct RandBitStringIP <: AbstractIP
    n::Integer
end

generate(ip::RandBitStringIP) = bitrand(ip.n)

max_stars_repo_name: DTUComputeTONIA/SDRLS


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# 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.

@with_kw mutable struct DiminishingStepParameters{T <: AbstractFloat}
    γ₀::T = 0.1
    η::T = 1.0
end

"""
    DiminishingStep

Functor object for using a constant step size in a quasigradient algorithm. Create by supplying a [`Constant`](@ref) object through `step` to `QuasiGradient.Optimizer` or by setting the [`StepSize`](@ref) attribute.

...
# Parameters
- `γ₀::AbstractFloat = 0.1`: Nominal step
- `η::AbstractFloat = 1.0`: Diminishing factor
...
"""
struct DiminishingStep{T <: AbstractFloat} <: AbstractStep
    parameters::DiminishingStepParameters{T}

    function DiminishingStep(::Type{T}; kw...) where T <: AbstractFloat
        return new{T}(DiminishingStepParameters{T}(; kw...))
    end
end
function step(step::DiminishingStep, k::Integer, ::Float64, ::AbstractVector, ::AbstractVector)
    @unpack γ₀, η = step.parameters
    return γ₀ / (1 + η * k)
end

# API
# ------------------------------------------------------------
"""
    Diminishing

Factory object for [`DiminishingStep`](@ref). Pass to `step` in `Quasigradient.Optimizer` or set the [`StepSize`](@ref) attribute. See ?DiminishingStep for parameter descriptions.

"""
struct Diminishing <: AbstractStepSize
    parameters::DiminishingStepParameters{Float64}
end
DiminishingStepParame(γ::AbstractFloat) = Constant(DiminishingStepParameters(; γ = Float64(γ)))
Diminishing(; kw...) = Diminishing(DiminishingStepParameters(; kw...))

function (stepsize::Diminishing)(::Type{T}) where T <: AbstractFloat
    return DiminishingStep(T; type2dict(stepsize.parameters)...)
end

function str(::Diminishing)
    return "diminishing step size"
end


text:
using MittagLefflerFunctions
import RationalApproximations: locate_extrema_idx, aaa_v2, minimax
using PyPlot

α = 0.75
β = 1.0
N = 4
Equad = MLQuad(α, β, N, :hyperbola)
E_ref = MLQuad(α, β, 14, :hyperbola)

function f(t)
    x = (1-t) / (1+t)
    return E_ref(-x)
end

t = collect(range(-1, 1, length=1000))
max_m = 2N+1
tol = 1e-14
F = f.(t)
r, err = aaa_v2(F, t, max_m, tol)
#m = length(r)
m = 2N
resid = F - r[m].(t)
idx = locate_extrema_idx(resid, m)

figure(1)
plot(t, resid, t[idx], resid[idx], "o")
xlabel(L"t")
grid(true)
title("Error in AAA approximation with m = $m")

max_iterations = 5
clamp = :no
opt_r, zmin, zmax = minimax(f, t[idx], r[m].supp_pt, (-1.0,1.0),
                             max_iterations, clamp)
opt_resid = F - opt_r[end].(t)

figure(2)
x = (1 .- t)./(1 .+ t)
Equad_resid = Equad.(-x) - F
plot(t[1:end-1], Equad_resid[1:end-1])
grid(true)
xlabel(L"$t$")
title("Error in quadrature approximation with N = $N")

figure(3)
plot(t, opt_resid, "-")
legend()
xlabel(L"t")
title("Error in best rational approximation with m = $m")
grid(true)

figure(4)
loglog(x[1:end-1], abs.(Equad_resid[1:end-1]), "-", label="N = $N")
loglog(x[1:end-1], abs.(opt_resid[1:end-1]), "--", label="m = $m")
#axis([-1.1, 1.1, 5e-11, 5e-5])
legend()
xlabel(L"t")
grid(true)
savefig("figure14.pdf")


text:
leftgauge(Ψ::UniformCMPS, args...; kwargs...) = leftgauge!(copy(Ψ), args...; kwargs...)
rightgauge(Ψ::UniformCMPS, args...; kwargs...) = rightgauge!(copy(Ψ), args...; kwargs...)

function leftgauge!(Ψ::UniformCMPS, C₀ = one(Ψ.Q);
                    maxreorth = 10,
                    eigalg = Arnoldi(; krylovdim = min(64, length(Ψ.Q[0]))),
                    kwargs...)

    U = one(Ψ.Q)
    if Ψ.gauge == :l
        η = norm(LeftTransfer(Ψ)(U))
        info = ConvergenceInfo(1, nothing, η, 0, 1)
        return Ψ, zero(scalartype(Ψ)), U, info
    end
    tol = eigalg.tol
    ρL, λ, info = leftenv(Ψ, C₀'*C₀; eigalg = eigalg, kwargs...)
    D, V = eigen!(Hermitian(ρL[]))
    Dsqrt = sqrt.(max.(D, defaulttol(D)))
    _, C = qr!(Diagonal(Dsqrt)*V')
    C ./= norm(C)
    CL = UpperTriangular(C)

    Ψ.Q = Constant(rdiv!(CL*Ψ.Q[], CL))
    Qdiag = view(Ψ.Q[], diagind(Ψ.Q[]))
    Qdiag .-= λ
    Ψ.Rs = map(R->Constant(rdiv!(CL*R[], CL)), Ψ.Rs)
    numreorth = 0
    η = norm(LeftTransfer(Ψ)(U))
    numiter = info.numiter
    numops = info.numops
    while η > tol
        ρL, dλ, info = leftenv(Ψ, U; eigalg = eigalg, kwargs...)
        numiter += info.numiter
        numops += info.numops
        D, V = eigen!(Hermitian(ρL[]))
        Dsqrt = sqrt.(max.(D, defaulttol(D)))
        _, dC = qr!(Diagonal(Dsqrt)*V')
        dCL = UpperTriangular(dC)

        Ψ.Q = Constant(rdiv!(dCL*Ψ.Q[], dCL))
        Qdiag = view(Ψ.Q[], diagind(Ψ.Q[]))
        Qdiag .-= dλ
        Ψ.Rs = map(R->Constant(rdiv!(dCL*R[], dCL)), Ψ.Rs)
        λ += dλ
        C = lmul!(dCL, C)
        C ./= norm(C)
        η = norm(LeftTransfer(Ψ)(U))
        numreorth += 1
        if numreorth == maxreorth
            break
        end
    end
    converged = Int(η < tol)
    if converged > 0
        Ψ.gauge = :l
    end
    return Ψ, λ, Constant(C), ConvergenceInfo(converged, nothing, η, numiter, numops)
end

function rightgauge!(Ψ::UniformCMPS, C₀ = one(Ψ.Q); kwargs...)
    if Ψ.gauge == :r
        U = one(Ψ.Q)
        η = norm(RightTransfer(Ψ)(U))
        info = ConvergenceInfo(1, nothing, η, 0, 1)
        return Ψ, zero(scalartype(Ψ)), U, info
    end
    _adjoint!(Ψ.Q[])
    foreach(Ψ.Rs) do R
        _adjoint!(R[])
    end
    Ψ.gauge = :n
    _, λ, C, info = leftgauge!(Ψ, C₀'; kwargs...)
    _adjoint!(Ψ.Q[])
    foreach(Ψ.Rs) do R
        _adjoint!(R[])
    end
    _adjoint!(C[])
    if Ψ.gauge == :l
        Ψ.gauge = :r
    end
    return Ψ, λ, C, info
end
#
# function leftcanonical!(Ψ::UniformCMPS, C₀ = one(Ψ.Q); maxreorth = 10, kwargs...)
#     Ψ, λ, C, infoL = leftgauge!(Ψ, C₀; maxreorth = maxreorth, kwargs...)
#     ρR, infoR = rightenv(Ψ; kwargs...)
#