text:
# test using cross-validation and a search algorithm to optimise the value of a hyper-parameter for scale of factors
using RatPack
using Random
using Distributions
using PlotlyJS

# generate large simulated dataset, and run base cross-validation on it, and optimise hyper-parameters
include("test_cv_3.jl")

#######################################
# optimise the scale factor, given home-ground advantages
r = [0.0, 400.0]
ε = 1.0e-0
score(f) = score_ratings( srule, test_competitions,
                         UpdateIterate( rule=UpdateEloF(;r0 = r0, K = 30.0, θ=θ, factor_scale=f ), batch_size=batch_size ),
                         update_ratings(UpdateIterate( rule=UpdateEloF(;r0 = r0, K = 30.0, θ=θ, factor_scale=f ), batch_size=batch_size ), initial_ratings, input_competitions)
                         )
(x_star2, y_star2, error_estimate2, n_evaluations2) = linesearch( "Golden", x->score(x); range=r, ε=ε)

data1 = scatter(;x=f_s, y=score_f, mode="lines+markers", name="score as a function of f_s")
data2 = scatter(;x=[0, maximum(score_f)], y=[s2, s2], mode="lines", name="upper bound")
data3 = scatter(;x=[x_star2, x_star2], y=[y_star2, y_star2], mode="markers", marker=attr(symbol="diamond-dot", size=12), name="maximum")

data = [data1, data2, data3]
layout = Layout(;
                xaxis=attr(title="f_s" ),
                yaxis=attr(title="(normalised) Brier score" ),
                )
p2 = PlotlyJS.plot( data, layout )


text:
export CoNLL
module CoNLL

function read(f, path::String)
    doc = []
    sent = []
    lines = open(readlines, path)
    for line in lines
        line = chomp(line)
        if length(line) == 0
            length(sent) > 0 && push!(doc, sent)
            sent = []
        else
            items = Vector{String}(split(line,'\t'))
            push!(sent, f(items))
        end
    end
    length(sent) > 0 && push!(doc, sent)
    T = typeof(doc[1][1])
    Vector{Vector{T}}(doc)
end
read(path::String) = read(identity, path)

end

ext: jl

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
# ---
# title: 153. Find Minimum in Rotated Sorted Array
# id: problem153
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/find-minimum-in-rotated-sorted-array/description/>
# hidden: true
# ---
#
# Suppose an array of length `n` sorted in ascending order is **rotated**
# between `1` and `n` times. For example, the array `nums = [0,1,2,4,5,6,7]`
# might become:
#
#   * `[4,5,6,7,0,1,2]` if it was rotated `4` times.
#   * `[0,1,2,4,5,6,7]` if it was rotated `7` times.
#
# Notice that **rotating** an array `[a[0], a[1], a[2], ..., a[n-1]]` 1 time
# results in the array `[a[n-1], a[0], a[1], a[2], ..., a[n-2]]`.
#
# Given the sorted rotated array `nums`, return _the minimum element of this
# array_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [3,4,5,1,2]
#     Output: 1
#     Explanation: The original array was [1,2,3,4,5] rotated 3 times.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [4,5,6,7,0,1,2]
#     Output: 0
#     Explanation: The original array was [0,1,2,4,5,6,7] and it was rotated 4 times.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [11,13,15,17]
#     Output: 11
#     Explanation: The original array was [11,13,15,17] and it was rotated 4 times.
#
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `1 <= n <= 5000`
#   * `-5000 <= nums[i] <= 5000`
#   * All the integers of `nums` are **unique**.
#   * `nums` is sorted and rotated between `1` and `n` times.
#
#
## @lc code=start
using LeetCode

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


text:
using Test

@testset "NoiseCC" begin

    @test 1==1
    @test 2==2

end

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


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

@doc raw"""A node selector represents the union of the results of one or more label queries over a set of nodes; that is, it represents the OR of the selectors represented by the node selector terms.

    IoK8sApiCoreV1NodeSelector(;
        nodeSelectorTerms=nothing,
    )

    - nodeSelectorTerms::Vector{IoK8sApiCoreV1NodeSelectorTerm} : Required. A list of node selector terms. The terms are ORed.
"""
mutable struct IoK8sApiCoreV1NodeSelector <: SwaggerModel
    nodeSelectorTerms::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1NodeSelectorTerm} } # spec name: nodeSelectorTerms

    function IoK8sApiCoreV1NodeSelector(;nodeSelectorTerms=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeSelector, Symbol("nodeSelectorTerms"), nodeSelectorTerms)
        setfield!(o, Symbol("nodeSelectorTerms"), nodeSelectorTerms)
        o
    end
end # type IoK8sApiCoreV1NodeSelector

const _property_map_IoK8sApiCoreV1NodeSelector = Dict{Symbol,Symbol}(Symbol("nodeSelectorTerms")=>Symbol("nodeSelectorTerms"))
const _property_types_IoK8sApiCoreV1NodeSelector = Dict{Symbol,String}(Symbol("nodeSelectorTerms")=>"Vector{IoK8sApiCoreV1NodeSelectorTerm}")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeSelector }) = collect(keys(_property_map_IoK8sApiCoreV1NodeSelector))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeSelector }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeSelector[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeSelector }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeSelector[property_name]

function check_required(o::IoK8sApiCoreV1NodeSelector)
    (getproperty(o, Symbol("nodeSelectorTerms")) === nothing) && (return false)
    true
end

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


text:
"""
$(TYPEDEF)

Hypograph of perspective function of sum-log cone of dimension `dim`.

    $(FUNCTIONNAME){T}(dim::Int, use_dual::Bool = false)
"""
mutable struct HypoPerLog{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    ϕ::T
    ζ::T
    tempw::Vector{T}

    function HypoPerLog{T}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real}
        @assert dim >= 3
        cone = new{T}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        return cone
    end
end

function setup_extra_data!(cone::HypoPerLog{T}) where {T <: Real}
    d = cone.dim - 2
    cone.tempw = zeros(T, d)
    return cone
end

get_nu(cone::HypoPerLog) = cone.dim

function set_initial_point!(arr::AbstractVector, cone::HypoPerLog)
    (arr[1], arr[2], w) = get_central_ray_hypoperlog(cone.dim - 2)
    @views arr[3:end] .= w
    return arr
end

function update_feas(cone::HypoPerLog{T}) where {T <: Real}
    @assert !cone.feas_updated
    v = cone.point[2]
    @views w = cone.point[3:end]

    if v > eps(T) && all(>(eps(T)), w)
        u = cone.point[1]
        cone.ϕ = sum(log(wi / v) for wi in w)
        cone.ζ = v * cone.ϕ - u
        cone.is_feas = (cone.ζ > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::HypoPerLog{T}) where {T <: Real}
    u = cone.dual_point[1]
    @views w = cone.dual_point[3:end]

    if all(>(eps(T)), w) && u < -eps(T)
        v = cone.dual_point[2]
        sumlog = sum(log(w_i / -u) for w_i in w)
        return (v - u * (sumlog + length(w)) > eps(T))
    end

    return false
end

function update_grad(cone::HypoPerLog)
    @assert cone.is_feas
    u = cone.point[1]
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    ζ = cone.ζ
    g = cone.grad

    g[1] = inv(cone.ζ)
    g[2] = -(cone.ϕ - d) / ζ - inv(v)
    vζi1 = -1 - v / ζ
    @. g[3:end] = vζi1 / w

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::HypoPerLog)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    v = cone.point[2]
    @views w = cone.point[3:end]
    H = cone.hess.data
    g = cone.grad
    ζ = cone.ζ
    wivζi = cone.tempw
    d = length(w)
    σζi = (cone.ϕ - d) / ζ
    vζi = v / ζ
    @. wivζi = vζi / w

    # u, v
    H[1, 1] = ζ^-2
    H[1, 2] = -σζi / ζ
    H[2, 2] = v^-2 + abs2(σζi) + d / ζ / v

    # u, v, w
    vζi2 = -vζi / ζ
    c1 = ((cone.ϕ - d) * vζi - 1) / ζ
    @. H[1, 3:end] = vζi2 / w
    @. H[2, 3:end] = c1 / w

    # w
    @inbounds for j in eachindex(wivζi)
        j2 = 2 + j
        wivζij = wivζi[j]
        for i in 1:j
            H[2 + i, j2] = wivζi[i] * wivζij
        end
        H[j2, j2] -= g[j2] / w[j]
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::HypoPerLog,
    )
    v = cone.point[2]
    @views w = cone.point[3:end]
    ζ = cone.ζ
    d = length(w)
    σ = cone.ϕ - d
    rwi = cone.tempw
    vζi1 = v / ζ + 1

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @. @views rwi = arr[3:end, j] / w

        qζi = q / ζ
        c0 = sum(rwi) / ζ
        # ∇ϕ[r] = v * c0
        c1 = (v * c0 - p / ζ + σ * qζi) / ζ
        c3 = c1 * v - qζi
        prod[1, j] = -c1
        prod[2, j] = c1 * σ - c0 + (qζi * d + q / v) / v
        @. prod[3:end, j] = (c3 + vζi1 * rwi) / w
    end

    return prod
end

function update_inv_hess(cone::HypoPerLog)
    @assert cone.grad_updated
    @assert !cone.inv_hess_updated
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    Hi = cone.inv_hess.data
    ζ = cone.ζ
    ϕ = cone.ϕ
    ζv = ζ + v
    ζζvi = ζ / ζv
    c3 = v / (ζv + d * v)
    c0 = ϕ - d * ζζvi
    c2 = v * c3
    c4 = c2 * ζv
    c1 = v * ζζvi + c0 * c2

    Hi[1, 1] = abs2(v * ϕ) + ζ * (ζ + d * v) - d * abs2(ζ + v * ϕ) * c3
    Hi[1, 2] = c0 * c4
    Hi[2, 2] = c4

    @. Hi[1, 3:end] = c1 * w
    @. Hi[2, 3:end] = c2 * w

    @inbounds for j in eachindex(w)
        j2 = 2 + j
        Hi[j2, j2] += abs2(w[j])
    end
    @views mul!(Hi[3:end, 3:end], w, w', c2 / ζv, ζζvi)

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::HypoPerLog,
    )
    @assert cone.grad_updated
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    ζ = cone.ζ
    ϕ = cone.ϕ
    ζv = ζ + v
    ζζvi = ζ / ζv
    c3 = v / (ζv + d * v)
    c0 = ϕ - d * ζζvi
    c4 = v * c3 * ζv
    c6 = abs2(v * ϕ) + ζ * (ζ + d * v) - d * abs2(ζ + v * ϕ) * c3
    c7 = c4 * c0
    c8 = c7 + v * ζ
    rw = cone.tempw

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @. @views rw = arr[3:end, j] * w

        c1 = sum(rw) / ζv
        c5 = c0 * p + q + c1
        c2 = v * (ζζvi * p + c3 * c5)
        prod[1, j] = c6 * p + c7 * q + c8 * c1
        prod[2, j] = c4 * c5
        @. prod[3:end, j] = (c2 + ζζvi * rw) * w
    end

    return prod
end

function dder3(cone::HypoPerLog{T}, dir::AbstractVector{T}) where {T <: Real}
    @assert cone.grad_updated
    v = cone.point[2]
    @views w = cone.point[3:end]
    dder3 = cone.dder3
    p = dir[1]
    q = dir[2]
    ζ = cone.ζ
    d = length(w)
    σ = cone.ϕ - d
    viq = q / v
    viq2 = abs2(viq)
    rwi = cone.tempw
    vζi = v / ζ
    vζi1 = vζi + 1

    @. @views rwi = dir[3:end] / w
    c0 = sum(rwi)
    c7 = sum(abs2, rwi)
    ζiχ = (-p + σ * q + c0 * v) / ζ
    c4 = (viq * (-viq * d + 2 * c0) - c7) / ζ / 2
    c1 = (abs2(ζiχ) - v * c4) / ζ
    c3 = -(ζiχ + viq) / ζ
    c5 = c3 * q + vζi * viq2
    c6 = -2 * vζi * viq - c3 * v
    c8 = c5 + c1 * v

    dder3[1] = -c1
    dder3[2] = c1 * σ + (viq2 - (d * c5 + c6 * c0 + vζi * c7)) / v - c4
    @. dder3[3:end] = (c8 + rwi * (c6 + vζi1 * rwi)) / w

    return dder3
end

# see analysis in
# https://github.com/lkapelevich/HypatiaSupplements.jl/tree/master/centralpoints
function get_central_ray_hypoperlog(d::Int)
    if d <= 10
        # lookup points where x = f'(x)
        return central_rays_hypoperlog[d, :]
    end
    # use nonlinear fit for higher dimensions
    x = inv(d)
    if d <= 70
        u = 4.657876 * x ^ 2 - 3.116192 * x + 0.000647
        v = 0.424682 * x + 0.553392
        w = 0.760412 * x + 1.001795
    else
        u = -3.011166 * x - 0.000122
        v = 0.395308 * x + 0.553955
        w = 0.837545 * x + 1.000024
    end
    return [u, v, w]
end

const central_rays_hypoperlog = [
    -0.827838387  0.805102007  1.290927686
    -0.689607388  0.724605082  1.224617936
    -0.584372665  0.68128058  1.182421942
    -0.503499342  0.65448622  1.153053152
    -0.440285893  0.636444224  1.131466926
    -0.389979809  0.623569352  1.114979519
    -0.349255921  0.613978276  1.102013921
    -0.315769104  0.606589839  1.091577908
    -0.287837744  0.600745284  1.083013
    -0.264242734  0.596019009  1.075868782
    ]


text:
using Ensemble
using Test: @testset

include("TestStats.jl")
include("TestEnsembleSampler.jl")
include("TestEnsembleGibbs.jl")
include("TestEnsembleNest.jl")
include("TestEnsembleKombine.jl")
include("TestEnsemblePTSampler.jl")

using .TestStats
using .TestEnsembleSampler
using .TestEnsembleGibbs
using .TestEnsembleNest
using .TestEnsembleKombine
using .TestEnsemblePTSampler

@testset "Ensemble.jl package tests" begin
    TestStats.testall()
    TestEnsembleSampler.testall()
    TestEnsembleGibbs.testall()
    TestEnsembleNest.testall()
    TestEnsembleKombine.testall()
    TestEnsemblePTSampler.testall()
end


text:
# simple program to test MCL clustering

using Test
using Clustering

@testset "MCL" begin

Random.seed!(34568)

@testset "nonsquare matrices not supported" begin
    @test_throws DimensionMismatch mcl(zeros(Float64, 4, 3))
end

# initialize adjacency matrix of a weighted graph
nodes = [:bat, :bit, :cat, :fit, :hat, :hit]
edges = Tuple{Symbol, Symbol, Float64}[(:cat, :hat, 0.2), (:hat, :bat, 0.16),
                                      (:bat, :cat, 1.0), (:bat, :bit, 0.125),
                                      (:bit, :fit, 0.25), (:fit, :hit, 0.5),
                                      (:hit, :bit, 0.16)]
adj_matrix = zeros(Float64, length(nodes), length(nodes))
for edge in edges
    n1 = findfirst(isequal(edge[1]), nodes)
    n2 = findfirst(isequal(edge[2]), nodes)
    adj_matrix[n1, n2] = adj_matrix[n2, n1] = edge[3]
end
@assert issymmetric(adj_matrix)

@testset "fractional inflation param (1.8)" begin
    res = mcl(adj_matrix, display=:none, inflation=1.8)
    @test isa(res, MCLResult)
    local k = length(res.counts)
    # @show k
    @test k == 2
    @test all(a -> 1 <= a <= k, res.assignments)
    @test length(res.assignments) == length(nodes)
    @test length(res.counts) == k
    local c
    for c in 1:k
        @test count(==(c), res.assignments) == res.counts[c]
    end
    @test res.nunassigned == 0
    @test res.assignments == [1, 2, 1, 2, 1, 2]
end

@testset "integer inflation param (2)" begin
    res = mcl(adj_matrix, display=:none, inflation=2)
    @test isa(res, MCLResult)
    @test length(res.assignments) == length(nodes)
    @test res.nunassigned == 0
end

@testset "test non-integral expansion" begin
    # should not raise an exception
    res = mcl(adj_matrix, display=:none, inflation=1.5, expansion=1.5, save_final_matrix=true)
    @test isa(res, MCLResult)
    @test length(res.assignments) == length(nodes)
    @test size(res.mcl_adj) == size(adj_matrix) # test that the matrix is returned
end

@testset "allow_singles option" begin
    res = mcl(diagm(0 => [1.0, 1.0]), display=:none, allow_singles=true)
    @test length(res.counts) == 2
    @test res.assignments == [1, 2]
    @test res.counts == [1, 1]
    @test res.nunassigned == 0

    res = mcl(diagm(0 => [1.0, 1.0]), display=:none, allow_singles=false)
    @test length(res.counts) == 0
    @test res.assignments == [0, 0]
    @test res.nunassigned == 2
end

@testset "sparse input matrix" begin
    res = mcl(sparse(adj_matrix), display=:none, expansion=2)
    @test isa(res, MCLResult)
    @test length(res.assignments) == length(nodes)
    @test res.nunassigned == 0
    @test eltype(res.mcl_adj) === Float64

    # fractional powers not supported for sparse matrices
    @test_broken mcl(sparse(adj_matrix), display=:none, expansion=2.1)
end

@testset "use Float32 input" begin
    res = mcl(convert(Matrix{Float32},adj_matrix), display=:none, expansion=2)
    @test isa(res, MCLResult)
    @test eltype(res.mcl_adj) === Float32
end

end


text:
const _DIRTY_TOKEN = "dirty"
const _CLEAN_TOKEN = "clean"
function _ready_to_push(repodir)
    !_check_gitdir(repodir) && return true
    out = _run("git -C $(repodir) diff-index --quiet HEAD && echo $(_CLEAN_TOKEN) || echo $(_DIRTY_TOKEN)";
        ignorestatus = false, verbose = false
    )
    return out == _CLEAN_TOKEN
end


text:
using ModelingToolkit, StaticArrays, LinearAlgebra, SparseArrays
using DiffEqBase
using Test

canonequal(a, b) = isequal(simplify(a), simplify(b))

# Calculus
@parameters t σ ρ β
@variables x y z
@test isequal((Differential(z) * Differential(y) * Differential(x))(t),
              Differential(z)(Differential(y)(Differential(x)(t))))

@test canonequal(ModelingToolkit.derivative(sin(cos(x)), x),
                 -sin(x) * cos(cos(x)))

@register_symbolic no_der(x)
@test canonequal(ModelingToolkit.derivative([sin(cos(x)), hypot(x, no_der(x))], x),
                 [
                     -sin(x) * cos(cos(x)),
                     x / hypot(x, no_der(x)) +
                     no_der(x) * Differential(x)(no_der(x)) / hypot(x, no_der(x)),
                 ])

@register_symbolic intfun(x)::Int
@test ModelingToolkit.symtype(intfun(x)) === Int

eqs = [σ * (y - x),
    x * (ρ - z) - y,
    x * y - β * z]

simpexpr = [:($(*)(σ, $(+)(y, $(*)(-1, x))))
            :($(+)($(*)(x, $(+)(ρ, $(*)(-1, z))), $(*)(-1, y)))
            :($(+)($(*)(x, y), $(*)(-1, z, β)))]

σ, β, ρ = 2 // 3, 3 // 4, 4 // 5
x, y, z = 6 // 7, 7 // 8, 8 // 9
for i in 1:3
    @test eval(ModelingToolkit.toexpr.(eqs)[i]) == eval(simpexpr[i])
    @test eval(ModelingToolkit.toexpr.(eqs)[i]) == eval(simpexpr[i])
end

@parameters t σ ρ β
@variables x y z
∂ = ModelingToolkit.jacobian(eqs, [x, y, z])
for i in 1:3
    ∇ = ModelingToolkit.gradient(eqs[i], [x, y, z])
    @test canonequal(∂[i, :], ∇)
end

@test all(canonequal.(ModelingToolkit.gradient(eqs[1], [x, y, z]), [σ * -1, σ, 0]))
@test all(canonequal.(ModelingToolkit.hessian(eqs[1], [x, y, z]), 0))

du = [x^2, y^3, x^4, sin(y), x + y, x + z^2, z + x, x + y^2 + sin(z)]
reference_jac = sparse(ModelingToolkit.jacobian(du, [x, y, z]))

@test findnz(ModelingToolkit.jacobian_sparsity(du, [x, y, z]))[[1, 2]] ==
      findnz(reference_jac)[[1, 2]]

let
    @variables t x(t) y(t) z(t)
    @test ModelingToolkit.exprs_occur_in([x, y, z], x^2 * y) == [true, true, false]
end

@test isequal(ModelingToolkit.sparsejacobian(du, [x, y, z]), reference_jac)

using ModelingToolkit

rosenbrock(X) =
    sum(1:(length(X) - 1)) do i
        100 * (X[i + 1] - X[i]^2)^2 + (1 - X[i])^2
    end

@variables a, b
X = [a, b]

spoly(x) = simplify(x, expand = true)
rr = rosenbrock(X)

reference_hes = ModelingToolkit.hessian(rr, X)
@test findnz(sparse(reference_hes))[1:2] ==
      findnz(ModelingToolkit.hessian_sparsity(rr, X))[1:2]

sp_hess = ModelingToolkit.sparsehessian(rr, X)
@test findnz(sparse(reference_hes))[1:2] == findnz(sp_hess)[1:2]
@test isequal(map(spoly, findnz(sparse(reference_hes))[3]), map(spoly, findnz(sp_hess)[3]))

Joop, Jiip = eval.(ModelingToolkit.build_function(∂, [x, y, z], [σ, ρ, β], t))
J = Joop([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J isa Matrix
J2 = copy(J)
Jiip(J2, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J2 == J

Joop, Jiip = eval.(ModelingToolkit.build_function(vcat(∂, ∂), [x, y, z], [σ, ρ, β], t))
J = Joop([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J isa Matrix
J2 = copy(J)
Jiip(J2, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J2 == J

Joop, Jiip = eval.(ModelingToolkit.build_function(hcat(∂, ∂), [x, y, z], [σ, ρ, β], t))
J = Joop([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J isa Matrix
J2 = copy(J)
Jiip(J2, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J2 == J

∂3 = cat(∂, ∂, dims = 3)
Joop, Jiip = eval.(ModelingToolkit.build_function(∂3, [x, y, z], [σ, ρ, β], t))
J = Joop([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test size(J) == (3, 3, 2)
J2 = copy(J)
Jiip(J2, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J2 == J

s∂ = sparse(∂)
@test nnz(s∂) == 8
Joop, Jiip = eval.(ModelingToolkit.build_function(s∂, [x, y, z], [σ, ρ, β], t,
                                                  linenumbers = true))
J = Joop([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test length(nonzeros(s∂)) == 8
J2 = copy(J)
Jiip(J2, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], 1.0)
@test J2 == J

# Function building

@parameters σ ρ β
@variables x y z
eqs = [σ * (y - x),
    x * (ρ - z) - y,
    x * y - β * z]
f1, f2 = ModelingToolkit.build_function(eqs, [x, y, z], [σ, ρ, β])
f = eval(f1)
out = [1.0, 2, 3]
o1 = f([1.0, 2, 3], [1.0, 2, 3])
f = eval(f2)
f(out, [1.0, 2, 3], [1.0, 2, 3])
@test all(o1 .== out)

function test_worldage()
    @parameters σ ρ β
    @variables x y z
    eqs = [σ * (y - x),
        x * (ρ - z) - y,
        x * y - β * z]
    f, f_iip = ModelingToolkit.build_function(eqs, [x, y, z], [σ, ρ, β];
                                              expression = Val{false})
    out = [1.0, 2, 3]
    o1 = f([1.0, 2, 3], [1.0, 2, 3])
    f_iip(out, [1.0, 2, 3], [1.0, 2, 3])
end
test_worldage()

## No parameters
@variables x y z
eqs = [(y - x)^2,
    x * (x - z) - y,
    x * y - y * z]
f1, f2 = ModelingToolkit.build_function(eqs, [x, y, z])
f = eval(f1)
out = zeros(3)
o1 = f([1.0, 2, 3])
f = eval(f2)
f(out, [1.0, 2, 3])
@test all(out .== o1)

# y ^ -1 test
g = let
    f(x, y) = x / y
    @variables x y
    ex = expand_derivatives(Differential(x)(f(x, y)))
    func_ex = build_function(ex, x, y)
    eval(func_ex)
end

@test g(42, 4) == 1 / 4

function test_worldage()
    @variables x y z
    eqs = [(y - x)^2,
        x * (x - z) - y,
        x * y - y * z]
    f, f_iip = ModelingToolkit.build_function(eqs, [x, y, z]; expression = Val{false})
    out = zeros(3)
    o1 = f([1.0, 2, 3])
    f_iip(out, [1.0, 2, 3])
end
test_worldage()

@test_nowarn muladd(x, y, 0)
@test promote(x, 0) == (x, identity(0))
@test_nowarn [x, y, z]'

let
    @register_symbolic foo(x)
    @variables t
    D = Differential(t)

    @test isequal(expand_derivatives(D(foo(t))), D(foo(t)))
    @test isequal(expand_derivatives(D(sin(t) * foo(t))),
                  cos(t) * foo(t) + sin(t) * D(foo(t)))
end

foo(; kw...) = kw
foo(args...; kw...) = args, kw
pp = :name => :cool_name

@named cool_name = foo()
@test collect(cool_name) == [pp]

@named cool_name = foo(42)
@test cool_name[1] == (42,)
@test collect(cool_name[2]) == [pp]

@named cool_name = foo(42; a = 2)
@test cool_name[1] == (42,)
@test collect(cool_name[2]) == [pp; :a => 2]

@named cool_name = foo(a = 2)
@test collect(cool_name) == [pp; :a => 2]

@named cool_name = foo(; a = 2)
@test collect(cool_name) == [pp; :a => 2]

@named cool_name = foo(name = 2)
@test collect(cool_name) == [:name => 2]

@named cool_name = foo(42; name = 3)
@test cool_name[1] == (42,)
@test collect(cool_name[2]) == [:name => 3]

kwargs = (; name = 3)
@named cool_name = foo(42; kwargs...)
@test cool_name[1] == (42,)
@test collect(cool_name[2]) == [:name => 3]

if VERSION >= v"1.5"
    name = 3
    @named cool_name = foo(42; name)
    @test cool_name[1] == (42,)
    @test collect(cool_name[2]) == [:name => name]
    @named cool_name = foo(; name)
    @test collect(cool_name) == [:name => name]

    ff = 3
    @named cool_name = foo(42; ff)
    @test cool_name[1] == (42,)
    @test collect(cool_name[2]) == [pp; :ff => ff]

    @named cool_name = foo(; ff)
    @test collect(cool_name) == [pp; :ff => ff]
end

foo(i; name) = i, name
@named goo[1:3] = foo(10)
@test isequal(goo, [(10, Symbol(:goo_, i)) for i in 1:3])
@named koo 1:3 i->foo(10i)
@test isequal(koo, [(10i, Symbol(:koo_, i)) for i in 1:3])


text:
using SparseUtils
import SparseArrays
import SparseArrays: SparseMatrixCSC
#import SparseUtils: materialize
import SparseUtils: SparseMatrixCOO
import LinearAlgebra
using Serialization
using Test

let # typeof(sparse_array) = SparseMatrixCSC
    sparse_array = open("sparse_array.dat", "r") do io
        deserialize(io)
    end
    expected_size = (30, 70)
    expected_sum = 1.6277046559059591
    @testset "measurements" begin
        @test size(sparse_array) == expected_size
        @test length(sparse_array) == prod(size(sparse_array))
    end

    @testset "density" begin
        expected_density = 0.008095238095238095
        @test density(sparse_array) == expected_density
    end

    @testset "transpose" begin
        @test transpose(sparse_array) |> copy |> size == (expected_size[2], expected_size[1])
        @test sparse_array |> transpose |> copy |> transpose |> copy == sparse_array
    end

    @testset "nnz" begin
        # nnz defined for columns
        @test sum(map(i -> SparseArrays.nnz(sparse_array, i),  1:size(sparse_array)[2])) == SparseArrays.nnz(sparse_array)
    end

    @testset "construction" begin
        S = SparseArrays.sparse([1], [1], [1], 1, 1; sparsetype=SparseMatrixCOO)
        @test isa(S, SparseMatrixCOO{Int,Int})
        @test size(S) == (1, 1)
        @test length(S) == 1
        S1 = SparseArrays.sparse([5, 7], [2, 1], [1.0, 2.0], 10, 10; sparsetype=SparseMatrixCOO)
        @test size(S1) == (10, 10)
    end

    let sparse_array_coo = SparseMatrixCOO(sparse_array)

        @testset "conversion" begin
            @test SparseMatrixCSC(sparse_array_coo) == sparse_array
        end

        @testset "transpose" begin
            pdcoo = permutedims(sparse_array_coo)
            @test copy(transpose(sparse_array_coo)) == pdcoo
            @test isa(transpose(sparse_array_coo), LinearAlgebra.Transpose)
            @test permutedims(pdcoo) == sparse_array_coo
            @test sparse_array_coo |> transpose |> transpose == sparse_array_coo
        end

        @testset "sum" begin
            @test sum(sparse_array) == sum(sparse_array_coo)
        end

        @testset "prod" begin
            @test prod(sparse_array) == prod(sparse_array_coo)
            zerotoone(x) = x == 0 ? one(x) : x
            @test prod(zerotoone, sparse_array) == prod(zerotoone, sparse_array_coo)
        end
    end

#    @test isapprox(sum(sparse_array), expected_sum)
end


text:
module Combination
using LinearAlgebra

export linear_add,SRSS,CQC,envelop

function linear_add(X)
    sum(X,dims=2)
end

function SRSS(X)
    sqrt.(sum(X.^2,dims=2))
end

function CQC(X)
    row,col=size(X)
    cqc=zeros(row)
    for i in 1:col
        for j in 1:col
            X[:,i].*X[:,j]
        end
    end
    sqrt.(cqc)
end

function envelop(X)
    maxX=maximum(X,dims=2)
    minX=mininum(X,dims=2)
    return maxX,minX
end

struct Comb()
    id::String
    hid::Int
    cases::Array
    factors::Array
    method::String
end

end


text:
struct Param end

struct ParameterValue{T} <: JuMP.AbstractVariable
    value::T
end

_invalid_init_error(msg) = error("Invalid initialization of parameter. " * msg * " not supported.")
function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo, ::Param)
    info.has_lb && _invalid_init_error("Lower bound")
    info.has_ub && _invalid_init_error("Upper bound")
    info.binary && _invalid_init_error("Binary")
    info.integer && _invalid_init_error("Integer")
    info.has_start && _invalid_init_error("Initial value")
    if !info.has_fix
        return ParameterValue(0.0)
    else
        return ParameterValue(info.fixed_value)
    end
end

function JuMP.add_variable(m::JuMP.Model, p::ParameterValue, name::String="")
    pref = add_parameter(m, p.value)
    if !isempty(name)
        JuMP.set_name(pref, name)
    end
    return pref
end


text:
export relu

function relu(x::Array)
    p = Ptr{Cvoid}[0]
    dnnPrimitiveAttributesCreate_F32(p)
    attr = p[1]

    layout = create_layout(x)

    p = Ptr{Cvoid}[0]
    dnnReLUCreateForward_F32(p, attr, layout, Cfloat(0))
    primitive = p[1]

    y = similar(x)
    dnnConversionExecute_F32(primitive, x, y)

    dnnPrimitiveAttributesDestroy_F32(attr)
    dnnLayoutDelete_F32(layout)
    dnnDelete_F32(primitive)
    y
end


text:
module TableCellCounts

import DataFrames
import ScientificTypes
import StatsBase
import Tables

_hello_world() = "hello world"

include("public.jl")

end # module


text:
"""
Type: LS_Stopping (specialization of GenericStopping)
Methods: start!, stop!, update_and_start!, update_and_stop!, fill_in!, reinit!, status

LS_Stopping is designed to handle the stopping criterion of line search problems.
Let f:R→Rⁿ, then h(t) = f(x+td) where x and d are vectors and t is a scalar.
h is such that h:R→R.

Stopping structure for non-linear programming problems using NLPModels.
    Input :
       - pb         : an AbstractNLPModel
       - optimality_check : a stopping criterion through an admissibility function
       - state      : The information relative to the problem, see GenericState
       - (opt) meta : Metadata relative to stopping criterion.
       - (opt) main_stp : Stopping of the main loop in case we consider a Stopping
                          of a subproblem.
                          If not a subproblem, then nothing.

 Note:
 * The pb can be a LineModel defined in SolverTools.jl (https://github.com/JuliaSmoothOptimizers/SolverTools.jl)
 * It is possible to define those stopping criterion in a NLPStopping except NLPStopping
   uses vectors operations. LS_Stopping and it's admissible functions (Armijo and Wolfe are provided with Stopping.jl)
   uses scalar operations.
 * optimality_check(pb, state; kwargs...) -> Float64
   For instance, the armijo condition is: h(t)-h(0)-τ₀*t*h'(0) ⩽ 0
   therefore armijo(h, h_at_t) returns the maximum between h(t)-h(0)-τ₀*t*h'(0) and 0.
 """
mutable struct LS_Stopping <: AbstractStopping
    # problem
    pb :: Any

    # stopping criterion proper to linesearch
    optimality_check :: Function

    # shared information with linesearch and other stopping
    meta :: AbstractStoppingMeta

    # current information on linesearch
    current_state :: LSAtT

    # Stopping of the main problem, or nothing
    main_stp :: Union{AbstractStopping, Nothing}

    function LS_Stopping(pb             :: Any,
                         admissible     :: Function,
                         current_state  :: LSAtT;
                         meta           :: AbstractStoppingMeta = StoppingMeta(),
                         main_stp       :: Union{AbstractStopping, Nothing} = nothing,
                         kwargs...)

        if !(isempty(kwargs))
           meta = StoppingMeta(;kwargs...)
		end

        return new(pb, admissible, meta, current_state, main_stp)
    end

end

"""
_unbounded_problem_check!: If x gets too big it is likely that the problem is unbounded
                           This is a specialized version that takes into account
                           that the problem might be unbounded if the objective function
                           is unbounded from below.

Note: evaluate the objective function is state.ht is void.
"""
function _unbounded_problem_check!(stp  :: LS_Stopping,
                                   x    :: Iterate)

 if stp.current_state.ht == nothing && typeof(stp.pb) <: AbstractNLPModel && typeof(x) <: AbstractVector
     stp.current_state.ht = obj(stp.pb, x)
 end
 f_too_large = stp.current_state.ht != nothing && norm(stp.current_state.ht) >= stp.meta.unbounded_threshold

 stp.meta.unbounded_pb = f_too_large

 return stp
end

"""
_resources_check!: check if the optimization algorithm has exhausted the resources.

If the problem is an AbstractNLPModel check the number of evaluations of f and sum.
"""
function _resources_check!(stp    :: LS_Stopping,
                           x      :: Iterate)

 max_evals = false
 max_f     = false

 if typeof(stp.pb) <: AbstractNLPModel
  max_f = stp.meta.max_f < neval_obj(stp.pb)
  max_evals = stp.meta.max_eval < sum_counters(stp.pb)
 end

 # global limit diagnostic
 stp.meta.resources = max_evals || max_f

 return stp
end

"""
_optimality_check: compute the optimality score.

This is the NLP specialized version that takes into account the structure of the
LS_Stopping where the optimality_check function is an input.
"""
function _optimality_check(stp  :: LS_Stopping; kwargs...)

 optimality = stp.optimality_check(stp.pb, stp.current_state; kwargs...)
 stp.current_state.current_score = optimality

 return optimality
end

################################################################################
# line search admissibility functions
################################################################################
include("ls_admissible_functions.jl")


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

using ImageFiltering
using TestImages
using ImageShow # to display images clearly

c = testimage("c")
c

diff1,diff2 = Kernel.prewitt()
println(1/6) # 0.166666
display(diff1)
display(diff2)

c = testimage("c")
imfilter(c, 6. .* Kernel.prewitt())

diff1,diff2 = Kernel.sobel()
println(1/8) # 0.125
display(diff1)
display(diff2)

c = testimage("c")
imfilter(c, 8. .* Kernel.prewitt())


text:
using Base: JLOptions

using Test: @test, @testset, @test_broken, @test_logs, @test_throws

using Suppressor: @capture_err, @suppress

using SimpleMock

@testset "SimpleMock.jl" begin
    @testset "Mock type" begin
        include("mock_type.jl")
    end
    @testset "mock function" begin
        include("mock_fun.jl")
    end
end


text:
using DReal
import DReal: opensmt_context, opensmt_set_verbosity, opensmt_assert, opensmt_check, opensmt_print_expr, opensmt_mk_context,
              opensmt_mk_bool_var, opensmt_mk_int_var, opensmt_mk_minus, opensmt_mk_leq, opensmt_mk_num_from_string,
              opensmt_push, opensmt_pop, opensmt_reset, opensmt_del_context, opensmt_mk_uminus, opensmt_mk_times_2, opensmt_mk_times_3,
              opensmt_set_precision, opensmt_mk_unbounded_real_var, opensmt_mk_num, opensmt_mk_div, opensmt_mk_sqrt, opensmt_mk_pow,
              opensmt_mk_exp, opensmt_define_ode, opensmt_mk_integral, opensmt_mk_eq, opensmt_mk_geq, opensmt_mk_real_var
using Base.Test
using Compat

ctx = opensmt_mk_context(@compat UInt32(1))
opensmt_set_precision(ctx, 1.0)

# Creating Real variables
x = opensmt_mk_real_var(ctx, "x", -100.0, 100.0)
P = opensmt_mk_real_var(ctx, "P", -100.0, 100.0)
x_0 = opensmt_mk_real_var(ctx, "x_0", -100.0, 100.0)
P_0 = opensmt_mk_real_var(ctx, "P_0", -100.0, 100.0)
x_t = opensmt_mk_real_var(ctx, "x_t", -100.0, 100.0)
P_t = opensmt_mk_real_var(ctx, "P_t", -100.0, 100.0)
time_0 = opensmt_mk_real_var(ctx, "time_0", 0.0, 100.0)

# =======================================================================
# Define ODE (flow)
# d/dt[x] = 1.0
# d/dt[P] = 1.0 / ((2.0 * 3.14159265359) ^ 0.5) * exp (- (x ^ 2) / 2.0)
zero = opensmt_mk_num(ctx, 0.0)
one = opensmt_mk_num(ctx, 1.0)
ten = opensmt_mk_num(ctx, 10.0)
neg_ten = opensmt_mk_uminus(ctx, ten)
two = opensmt_mk_num(ctx, 2.0)
pi  = opensmt_mk_num(ctx, 3.14159265359)
two_pi = opensmt_mk_times_2(ctx, two, pi)
# p_rhs = opensmt_mk_uminus(ctx, x)
p_rhs = opensmt_mk_times_2(
  ctx,
  # 1.0 / sqrt(2 * pi)
  # opensmt_mk_div(ctx, one, opensmt_mk_sqrt(ctx, two_pi)),
  opensmt_mk_div(ctx, one, opensmt_mk_pow(ctx, two_pi, opensmt_mk_num(ctx, 0.5))),
  # exp (- x^2 / 2.0)
  opensmt_mk_exp(
    ctx,
    opensmt_mk_uminus(ctx, opensmt_mk_div(ctx, opensmt_mk_pow(ctx, x, two), two))))
#
vars = [x, P]
rhses = [one, p_rhs]
opensmt_define_ode(ctx, "flow_1", vars, rhses, Cuint(2))

# =======================================================================
## Define initial and final state
opensmt_assert(ctx, opensmt_mk_eq(ctx, x_0, neg_ten))  # x_0 = -10.0
opensmt_assert(ctx, opensmt_mk_eq(ctx, P_0, zero))     # P_0 = 0.0
opensmt_assert(ctx, opensmt_mk_eq(ctx, x_t, ten))      # x_t = 10.0
opensmt_assert(ctx, opensmt_mk_geq(ctx, P_t, zero))    # P_t >= 0.0

# =======================================================================
# Add integral constraint
# [x_t P_t] = integral_0^{time_0} [x_0 P_0] | flow_1

vars_t = [x_t, P_t]
vars_0 = [x_0, P_0]
integral = opensmt_mk_integral(ctx, vars_t, zero, time_0, vars_0, Cuint(2), "flow_1")
opensmt_assert(ctx, integral)
res = opensmt_check(ctx)

println( "res = $(res == -1 ? "unsat" : "sat")")

print( "Deleting context\n" )
opensmt_del_context( ctx )


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

using Markdown
using InteractiveUtils

# ╔═╡ 857d0150-b0ea-11eb-2bcb-b70a5fc0a6b7
md"""
# [Day 25: Let It Snow](https://adventofcode.com/2015/day/25)

Merry Christmas! Santa is booting up his weather machine; looks like you might get a [white Christmas](https://adventofcode.com/2015/day/1) after all.

The weather machine beeps! On the console of the machine is a copy protection message asking you to [enter a code from the instruction manual](https://en.wikipedia.org/wiki/Copy_protection#Early_video_games). Apparently, it refuses to run unless you give it that code. No problem; you'll just look up the code in the--

"Ho ho ho", Santa ponders aloud. "I can't seem to find the manual."

You look up the support number for the manufacturer and give them a call. Good thing, too - that 49th star wasn't going to earn itself.

"Oh, that machine is quite old!", they tell you. "That model went out of support six minutes ago, and we just finished shredding all of the manuals. I bet we can find you the code generation algorithm, though."

After putting you on hold for twenty minutes (your call is very important to them, it reminded you repeatedly), they finally find an engineer that remembers how the code system works.

The codes are printed on an infinite sheet of paper, starting in the top-left corner. The codes are filled in by diagonals: starting with the first row with an empty first box, the codes are filled in diagonally up and to the right. This process repeats until the [infinite paper is covered](https://en.wikipedia.org/wiki/Cantor's_diagonal_argument). So, the first few codes are filled in in this order:

```
   | 1   2   3   4   5   6
---+---+---+---+---+---+---+
 1 |  1   3   6  10  15  21
 2 |  2   5   9  14  20
 3 |  4   8  13  19
 4 |  7  12  18
 5 | 11  17
 6 | 16

```

For example, the 12th code would be written to row `4`, column `2`; the 15th code would be written to row `1`, column `5`.

The voice on the other end of the phone continues with how the codes are actually generated. The first code is `20151125`. After that, each code is generated by taking the previous one, multiplying it by `252533`, and then keeping the remainder from dividing that value by `33554393`.

So, to find the second code (which ends up in row `2`, column `1`), start with the previous value, `20151125`. Multiply it by `252533` to get `5088824049625`. Then, divide that by `33554393`, which leaves a remainder of `31916031`. That remainder is the second code.

"Oh!", says the voice. "It looks like we missed a scrap from one of the manuals. Let me read it to you." You write down his numbers:

```
   |    1         2         3         4         5         6
---+---------+---------+---------+---------+---------+---------+
 1 | 20151125  18749137  17289845  30943339  10071777  33511524
 2 | 31916031  21629792  16929656   7726640  15514188   4041754
 3 | 16080970   8057251   1601130   7981243  11661866  16474243
 4 | 24592653  32451966  21345942   9380097  10600672  31527494
 5 |    77061  17552253  28094349   6899651   9250759  31663883
 6 | 33071741   6796745  25397450  24659492   1534922  27995004

```

"Now remember", the voice continues, "that's not even all of the first few numbers; for example, you're missing the one at 7,1 that would come before 6,2. But, it should be enough to let your-- oh, it's time for lunch! Bye!" The call disconnects.

Santa looks nervous. Your puzzle input contains the message on the machine's console. **What code do you give the machine?**
"""

# ╔═╡ 36909b22-41bb-4b02-b2d7-a662da7ffbf3
puzzle_input = "To continue, please consult the code grid in the manual.  Enter the code at row 2978, column 3083."

# ╔═╡ Cell order:
# ╟─857d0150-b0ea-11eb-2bcb-b70a5fc0a6b7
# ╟─36909b22-41bb-4b02-b2d7-a662da7ffbf3


text:
module Utils

Base.Experimental.@optlevel 3

export linspace, half
export yaml_repr, struct2dict
export ORI

import IntervalSets: width
using IntervalSets

const ORI{T} = Interval{:closed,:open,T}
ORI(i::Interval{L,R,T}) where {L,R,T} = ORI{T}(i)

function _atype end
function half end
function width end

width(xs::AbstractVector{T}) where T = width(Interval(extrema(xs)...))
width(xs::AbstractRange{T}) where T = abs(last(xs) - first(xs))
half(xs::AbstractVector{T}) where T = width(xs) / 2
half(i::Interval{L,R,T}) where {L,R,T} = width(i) / 2

_atype(::Type{T}) where {T <: AbstractArray} = T
_atype(a::AbstractArray{T}) where T = typeof(a)

@inline function linspace(start::T, stop::U, len::Integer) where {
    T <: Number,
    U <: Number,
}
    range(start, stop, length = len)
end

@inline function linspace(::Type{T}, start::X, stop::Y, len::Integer) where {
    T <: Number,
    X <: Number,
    Y <: Number,
}
    range(T(start), T(stop), length = len)
end

@inline function linspace(
    i::ClosedInterval{T},
    len::Integer,
) where {T <: Number}
    range(i, len)
end

@inline function linspace(
    ::Type{T},
    i::ClosedInterval{U},
    len::Integer,
) where {T <: Number,U <: Number}
    range(convert(ClosedInterval{T}, i), len)
end

@inline function linspace(
    i::ORI{T},
    len::Integer,
) where {T <: Number}
    range(i, len)
end

@inline function linspace(
    ::Type{T},
    i::ORI{U},
    len::Integer,
) where {T <: Number,U <: Number}
    range(convert(ORI{T}, i), len)
end

yaml_repr(obj) = obj
yaml_repr(ntup::NamedTuple) = struct2dict(ntup)

function struct2dict(obj)
    Dict(map(propertynames(obj)) do p
        field = getproperty(obj, p)
        p => yaml_repr(field)
    end)
end

@inline function _compute_radius(x::T, y::T, d::T)::T where {T <: Real}
    x == zero(T) && return abs(y) * d + one(T)
    y == zero(T) && return abs(x) * d + one(T)
    return √(x^2 + y^2) * d + one(T)
end

@inline function _compute_radius(x::T, y::T)::T where {T <: Real}
    x == 0 && return abs(y) + one(T)
    y == 0 && return abs(x) + one(T)
    return √(x^2 + y^2) + one(T)
end

@inline function _compute_angle(x::T, y::T, d::T)::T where {T <: Real}
    θ = atan(y, x)
    θ < 0 && return (θ + 2π) * d + one(T)
    θ * d + one(T)
end

@inline function _wrap_angle(θ::T, n::Integer, nh::Real)::T where {T <: Real}
    θ >= nh && return one(T)
    θ > n && return T(n)
    θ
end

end # module


text:
using MPI
using P4est_wrapper
using Test

# Skip tests if library was not properly loaded
if P4est_wrapper.P4EST_FOUND

    nprocs = min(Sys.CPU_THREADS, 2)
    dir = dirname(@__FILE__)
    julia = Base.julia_cmd()
    mpiexec = MPI.mpiexec_path
    @testset "Test SC bindings" begin
        cmd = `$mpiexec -n $nprocs --allow-run-as-root --tag-output $(julia) $(joinpath(dir, "test_sc_bindings.jl"))`
        @show cmd
        run(cmd)
    end
    @testset "Test P4EST bindings" begin
        cmd = `$mpiexec -n $nprocs --allow-run-as-root --tag-output $(julia) $(joinpath(dir, "test_p4est_bindings.jl"))`
        @show cmd
        run(cmd)
    end
    @testset "Test P8EST bindings" begin
        cmd = `$mpiexec -n $nprocs --allow-run-as-root --tag-output $(julia) $(joinpath(dir, "test_p8est_bindings.jl"))`
        @show cmd
        run(cmd)
    end
    @testset "Example 1" begin
        cmd = `$mpiexec -n $nprocs --allow-run-as-root --tag-output $(julia) $(joinpath(dir, "example_1.jl"))`
        @show cmd
        run(cmd)
    end
    @testset "Example 2" begin
        cmd = `$mpiexec -n $nprocs --allow-run-as-root --tag-output $(julia) $(joinpath(dir, "example_2.jl"))`
        @show cmd
        run(cmd)
    end
    @testset "Example 3" begin
        cmd = `$mpiexec -n $nprocs --allow-run-as-root --tag-output $(julia) $(joinpath(dir, "example_3.jl"))`
        @show cmd
        run(cmd)
    end
end