text:
using LinearAlgebra, StatsBase, Random, LaTeXStrings, Plots; pyplot()
Random.seed!(0)

L = 10
p0, p1 = 1/2, 3/4
beta = 0.75
pExplore(t) = t^-0.2
alpha(t) = t^-0.2
T = 10^6

function QlearnSim(kappa)
    P0 = diagm(1=>fill(p0,L-1)) + diagm(-1=>fill(1-p0,L-1))
    P0[1,1], P0[L,L] = 1 - p0, p0

    P1 = diagm(1=>fill(p1,L-1)) + diagm(-1=>fill(1-p1,L-1))
    P1[1,1], P1[L,L] = 1 - p1, p1

    R0 = collect(1:L)
    R1 = R0 .- kappa

    nextState(s,a) =
    	a == 0 ? sample(1:L,weights(P0[s,:])) : sample(1:L,weights(P1[s,:]))

    Q = zeros(L,2)
    s = 1
    optimalAction(s) = Q[s,1] >= Q[s,2] ? 0 : 1
    for t in 1:T
        if rand() < pExplore(t)
            a = rand([0,1])
        else
            a = optimalAction(s)
        end
        sNew = nextState(s,a)
        r = a == 0 ? R0[sNew] : R1[sNew]
        Q[s,a+1]=(1-alpha(t))*Q[s,a+1]+alpha(t)*(r+beta*max(Q[sNew,1],Q[sNew,2]))
        s = sNew
    end
    [optimalAction(s) for s in 1:L]
end

kappaGrid = 0.0:0.1:2.0
policyMap = zeros(L,length(kappaGrid))

for (i,kappa) in enumerate(kappaGrid)
    policyMap[:,i] = QlearnSim(kappa)
end

heatmap(policyMap, fill=cgrad([:blue, :red]),
        xticks=(0:1:21, -0.1:0.1:2), yticks=(0:L, 0:L),
        xlabel=L"\kappa", ylabel="State", colorbar_entry=false)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# let
#     W = 2
#     for T ∈ [Float32, Float64, Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64]
#         Wmax = max_vector_width(T)
#         while W ≤ Wmax

#             W += W
#         end
#     end
# end
# function broadcast_str(W::Int, typ::String)
#     vtyp = "<$W x $typ>"
#     """
#         %ie = insertelement $vtyp undef, $typ %0, i32 0
#         %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
#         ret $vtyp %v
#     """
# end
@inline vzero(::Val{1}, ::Type{T}) where {T<:NativeTypes} = zero(T)
@inline vzero(::StaticInt{1}, ::Type{T}) where {T<:NativeTypes} = zero(T)
@generated function vzero(::Union{Val{W},StaticInt{W}}, ::Type{T}) where {W,T<:NativeTypes}
    # isone(W) && return Expr(:block, Expr(:meta,:inline), Expr(:call, :zero, T))
    typ = LLVM_TYPES[T]
    instrs = "ret <$W x $typ> zeroinitializer"
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($instrs, _Vec{$W,$T}, Tuple{}))
    end
end
@generated function vbroadcast(::Union{Val{W},StaticInt{W}}, s::_T) where {W,_T<:NativeTypes}
    isone(W) && return :s
    if _T <: Integer && sizeof(_T) * W > REGISTER_SIZE
        T = pick_integer(W, sizeof(_T))
        if _T <: Unsigned
            T = unsigned(T)
        end
        # ssym = :(s % $T)
        ssym = :(convert($T, s))
    else
        T = _T
        ssym = :s
    end
    typ = LLVM_TYPES[T]
    vtyp = vtype(W, typ)
    instrs = """
        %ie = insertelement $vtyp undef, $typ %0, i32 0
        %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
        ret $vtyp %v
    """
    quote
        $(Expr(:meta,:pure,:inline))
        Vec(llvmcall($instrs, _Vec{$W,$T}, Tuple{$T}, $ssym))
    end
end

# for T ∈ [Float32,Float64,Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64]#, Float16]
#     maxW = pick_vector_width(T)
#     typ = LLVM_TYPES[T]
#     W = 2
#     while W ≤ maxW
#         instrs = "ret <$W x $typ> zeroinitializer"
#         @eval @inline vzero(::Val{$W}, ::Type{$T}) = Vec(llvmcall($instrs, Vec{$W,$T}, Tuple{}, ))
#         instrs = broadcast_str(W, typ)
#         # vtyp = "<$W x $typ>"
#         # instrs = """
#         # %ie = insertelement $vtyp undef, $typ %0, i32 0
#         # %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
#         # ret $vtyp %v
#         # """
#         @eval Base.@pure @inline vbroadcast(::Val{$W}, s::$T) = Vec(llvmcall($instrs, Vec{$W,$T}, Tuple{$T}, s))
#         W += W
#     end
# end

# @generated function vbroadcast(::Val{W}, s::Ptr{T}) where {W, T}
#     typ = JULIAPOINTERTYPE
#     instrs = broadcast_str(W, typ)
#     quote
#         $(Expr(:meta,:inline))
#         Vec(llvmcall( $instrs, _Vec{$W,Ptr{$T}}, Tuple{Ptr{$T}}, s ))
#     end
# end
# @generated function vbroadcast(::Val{W}, s::T) where {W, T <: NativeTypes}
#     typ = LLVM_TYPES[T]
#     instrs = broadcast_str(W, typ)
#     quote
#         $(Expr(:meta,:inline))
#         Vec(llvmcall( $instrs, _Vec{$W,$T}, Tuple{$T}, s))
#     end
# end
@generated function vbroadcast(::Union{Val{W},StaticInt{W}}, ptr::Ptr{T}) where {W, T}
    isone(W) && return Expr(:block, Expr(:meta, :inline), :(vload(ptr)))
    typ = LLVM_TYPES[T]
    ptyp = JULIAPOINTERTYPE
    vtyp = "<$W x $typ>"
    alignment = Base.datatype_alignment(T)
    instrs = """
        %ptr = inttoptr $ptyp %0 to $typ*
        %res = load $typ, $typ* %ptr, align $alignment
        %ie = insertelement $vtyp undef, $typ %res, i32 0
        %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
        ret $vtyp %v
    """
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall( $instrs, _Vec{$W,$T}, Tuple{Ptr{$T}}, ptr ))
    end
end
# @generated function Base.zero(::Type{Vec{W,T}}) where {W,T}
#     typ = LLVM_TYPES[T]
#     instrs = "ret <$W x $typ> zeroinitializer"
#     quote
#         $(Expr(:meta,:inline))
#         Vec(llvmcall($instrs, _Vec{$W,$T}, Tuple{}, ))
#     end
# end
@inline vbroadcast(::Union{Val{W},StaticInt{W}}, v::AbstractSIMDVector{W}) where {W} = v
@inline Vec{W,T}(v::Vec{W,T}) where {W,T} = v
# @inline vbroadcast(::Val{1}, s::T) where {T <: NativeTypes} = s
# @inline vbroadcast(::Val{1}, s::Ptr{T}) where {T <: NativeTypes} = s
@inline vzero(::Union{Val{W},StaticInt{W}}, ::Type{T}) where {W,T} = zero(Vec{W,T})
@inline Base.zero(::Type{Vec{W,T}}) where {W,T} = vzero(Val{W}(), T)
@inline Base.zero(::Vec{W,T}) where {W,T} = zero(Vec{W,T})
@inline Base.one(::Vec{W,T}) where {W,T} = vbroadcast(Val{W}(), one(T))

@inline Base.one(::Type{Vec{W,T}}) where {W,T} = vbroadcast(Val{W}(), one(T))
@inline Base.oneunit(::Type{Vec{W,T}}) where {W,T} = vbroadcast(Val{W}(), one(T))
@inline vzero(::Type{T}) where {T<:Number} = zero(T)
@inline vzero() = vzero(pick_vector_width_val(Float64), Float64)
# @inline sveczero(::Type{T}) where {T} = Svec(vzero(pick_vector_width_val(T)))
# @inline sveczero() = Svec(vzero(pick_vector_width_val(Float64)))

@inline Vec{W,T}(s::Integer) where {W,T<:Integer} = vbroadcast(Val{W}(), s % T)
@inline Vec{W,T}(s::Real) where {W,T} = vbroadcast(Val{W}(), T(s))
@inline Vec{W}(s::T) where {W,T<:NativeTypes} = vbroadcast(Val{W}(), s)
@inline Vec(s::T) where {T<:NativeTypes} = vbroadcast(pick_vector_width_val(T), s)

@generated function Base.zero(::Type{VecUnroll{N,W,T,V}}) where {N,W,T,V}
    t = Expr(:tuple); foreach(_ -> push!(t.args, :(zero(Vec{$W,$T}))), 0:N)
    Expr(:block, Expr(:meta, :inline), :(VecUnroll($t)))
end
@inline Base.zero(::VecUnroll{N,W,T,V}) where {N,W,T,V} = zero(VecUnroll{N,W,T,V})

@generated function VecUnroll{N,W,T,V}(x::S) where {N,W,T,V<:AbstractSIMDVector{W,T},S<:Real}
    t = Expr(:tuple)
    for n ∈ 0:N
        push!(t.args, :(convert($V, x)))
    end
    Expr(:block, Expr(:meta,:inline), :(VecUnroll($t)))
end
@inline VecUnroll{N,W,T}(x::NativeTypesV) where {N,W,T} = VecUnroll{N,W,T,Vec{W,T}}(x)

# @inline vbroadcast(::Union{Val{W},StaticInt{W}}, ::Type{T}, s::T) where {W,T} = vbroadcast(Val{W}(), s)
# @generated function vbroadcast(::Union{Val{W},StaticInt{W}}, ::Type{T}, s::S) where {W,T,S}
#     ex = if sizeof(T) < sizeof(S)
#         vbroadcast(Val{W}(), promote_type(T, S)(s))
#     else

#     end
#     Expr(:block, Expr(:meta, :inline), ex)
# end

ext: jl


text:
using Test

@testset "TinyModia examples with simulation" begin

include("SimpleFilters.jl")

include("FilterCircuit.jl")

include("CauerLowPassFilter.jl")

include("MotorControl.jl")

include("ServoSystem.jl")

end


text:
module Interpreter

export interpret
abstract type Node end

struct Add <: Node end
struct Sub <: Node end
struct Div <: Node end
struct Mul <: Node end

Base.show(io::IO, ::Add) = print(io, "+ₛ")
Base.show(io::IO, ::Sub) = print(io, "-ₛ")
Base.show(io::IO, ::Div) = print(io, "/ₛ")
Base.show(io::IO, ::Mul) = print(io, "*ₛ")

# const Value = Union{Term, Float64}

struct Term
  head::Node
  x::Union{Term, Float64}
  y::Union{Term, Float64}
end

a +ₛ b = Term(Add(), a, b)
a -ₛ b = Term(Sub(), a, b)
a *ₛ b = Term(Mul(), a, b)
a /ₛ b = Term(Div(), a, b)

Base.show(io::IO, term::Term) = print(io, term.x, " ", term.head, " ", term.y)

function interpret(t::Term)
  if t.head == Add()
    return interpret(t.x) + interpret(t.y)
  elseif t.head == Sub()
    return interpret(t.x) - interpret(t.y)
  elseif t.head == Div()
    return interpret(t.x) / interpret(t.y)
  elseif t.head == Mul()
    return interpret(t.x) * interpret(t.y)
  end
end

interpret(t::Float64) = t

function test()
  @show term = 15.3 *ₛ 12.2 +ₛ 41.8 /ₛ 5.2
  interpret(term)
end

end


text:
# # Mesh with implicit geometry

# In this tutorial, we will learn
#
#    -  How to create a mesh of hexahedra.
#    -  How to create a mesh using an implicit definition of the geometry.

# The geometry is a ramp in the form of 3/4 part of a circle.

# The tutorial will produce files for mesh visualization in the
# [Paraview](https://www.paraview.org/) format (VTK). One can display this
# information by loading the file with `paraview.exe`. When the tutorial is
# executed in `mybinder.org`, the graphics file needs to be downloaded to your
# desktop, and then visualized locally.

# We shall use static arrays for the connectivities and the coordinates.
using StaticArrays

# We start with a parametric domain: a rectangular block, with dimensions of
# different physical units. The extent of the domain: the first
# coordinate is angular. The next two coordinates are in the plane generating
# section.
L, W, H = 1.5*pi, 3.0, 1.5
R = 4.0

# This is the number of element edges along each dimension.
M, N, K = 13, 4, 2

# Now we define the indexing of the vertices. The linear index converts the
# triple of indexes along the three dimensions into a single number. the
# dimensions correspond to the number of nodes along direction.
linix = LinearIndices((M+1, N+1, K+1))
# The serial-number function converts the triple of indexes into a single linear
# index.
sn(m, n, k) = linix[m, n, k]

# Now we create the connectivities of the individual hexahedral cells. The
# numbering is given by the topologically regular character of the block domain.
C = SVector{8, Int}[]
for k in 1:K, n in 1:N, m in 1:M
    c = (
        sn(m, n, k), sn(m+1, n, k), sn(m+1, n+1, k), sn(m, n+1, k),
        sn(m, n, k+1), sn(m+1, n, k+1), sn(m+1, n+1, k+1), sn(m, n+1, k+1)
        )
    push!(C, SVector{8, Int}(c))
end

# The locations of the vertices is also perfectly regular within the parametric
# domain. The nodes are distributed uniformly along each direction. the location
# in the physical space is computed from the parametric coordinates: the first
# coordinate is used as an angle, and `R` is the interior radius of the circular
# ramp.
loc = SVector{3, Float64}[]
for k in 1:K+1, n in 1:N+1, m in 1:M+1
    a = (m-1)/M*L
    x, y, z = (R+(n-1)/N*W) * cos(a), (R+(n-1)/N*W) * sin(a), (k-1)/K*H * (1 + a)
    push!(loc, SVector{3, Float64}(x, y, z))
end

# At this point we can create the shape collections.
using MeshCore: ShapeColl, H8, P1
# This is the collection of hexahedral shapes.
elements = ShapeColl(H8, length(C), "elements")
# And this is the collection of the vertices.
vertices = ShapeColl(P1, length(loc), "vertices")

# The geometry is attached is an attribute. Here we use vector attribute which
# actually stores three coordinates per vertex.
using MeshCore: VecAttrib
vertices.attributes["geom"] = VecAttrib(loc)
using MeshCore: IncRel
connectivities = IncRel(elements, vertices, C)

# We now write out a postprocessing file with the elements. Note that we pass
# the `connectivities` incidence relation to get the elements stored in the
# file.
using MeshSteward: vtkwrite
vtkwrite("block", connectivities)

# Now we will define another shape collection. This one will refer to geometry
# defined implicitly: through a function that returns the coordinates  of a
# vertex based on its serial number. Remember that serial numbers are stored in
# the connectivities.

vertices2 = ShapeColl(P1, prod((M+1, N+1, K+1)), "vertices2")

# The Cartesian indices will allow us to map the linear index to a Cartesian
# triple of indices.
cart = CartesianIndices((M+1, N+1, K+1))
using MeshCore: FunAttrib

# Note that this function defines the same geometry in physical space as the
# loop above. However, the geometry is not stored, it is returned as a triple of
# numbers.
access_location(i) = let
    ci  = cart[i]
    m, n, k = ci[1], ci[2], ci[3]
    a = (m-1)/M*L
    x, y, z = (R+(n-1)/N*W) * cos(a), (R+(n-1)/N*W) * sin(a), (k-1)/K*H * (1 + a)
    SVector{3, Float64}(x, y, z)
end

# The geometry attribute is defined by providing access to the locations of the
# vertices through the function `access_location`.
vertices2.attributes["geom"] = FunAttrib(0.0, prod((M+1, N+1, K+1)), access_location);
connectivities = IncRel(elements, vertices2, C)
vtkwrite("block2", connectivities)

# Now both data sets have been written out. Load them into
# [Paraview](https://www.paraview.org/) and it should be possible to verify that
# those two meshes sets overlap perfectly.

true


text:
using Printf
using Random
using LinearAlgebra
using Distributed
using MAT
@everywhere using DistributedArrays
@everywhere using RCAM
@everywhere using Random
@everywhere Random.seed!(123)

#Load the data
fid = matopen("../data/X.mat")
d = read(fid)
X = d["X"]

function pMNtest(X)
# Choltest is just a wrapper for this script

# Params
r = 4 #Rank
alt = 5 #Iterations
perc = 0.8 #Percent of missing values
s1,s2 = size(X)

# Subsample
ind = randperm(s1*s2)
inds = ind[1:Int(round(perc*s1*s2))]
b = copy(X)
b[inds] .= 0.

#Gen Noise
N = rand(Float64,size(b))
noise = 0.5*norm(vec(X))*N/norm(vec(N))
noise[inds] .= 0.

eta = norm(vec(noise))/sqrt(size(X)[1])

# Interpolate
ownerL, ownerR = dclrMN(b+noise, eta, r, alt)

# Gather L and R results
L = fetch(@spawnat ownerL Main.L)
R = fetch(@spawnat ownerR Main.R)

return L,R
end

snr(raw,interp) = -20*log10(norm(interp-raw)/norm(raw))

L,R = pMNtest(X)
println(snr(X,L*R'))


text:
push!(LOAD_PATH, joinpath(dirname(@__FILE__), "../../.."))
push!(LOAD_PATH, dirname(@__FILE__))
ENV["USE_CUDA"] = true
workspace()
using Merlin

function bench()
  a =  Array{Float32,2}[]
  push!(a, zeros(Float32,100,100))
  push!(a, zeros(Float32,100,100))
  for i = 1:1000
    Merlin.concat(a, 1)
    #cat(1,a...)
  end
end

@time bench()

include("token.jl")
include("model.jl")
include("train.jl")
path = "C:/Users/hshindo/Dropbox/tagging"

@time train(path)

max_stars_repo_name: hshindo/POSTagging-benchmark

max_issues_repo_name: hshindo/POSTagging-benchmark

max_forks_repo_name: hshindo/POSTagging-benchmark


text:
using Lux, NNlib, Random, Statistics, Zygote

include("../utils.jl")

rng = Random.default_rng()
Random.seed!(rng, 0)

@testset "BatchNorm" begin
    let m = BatchNorm(2), x = [1.0f0 3.0f0 5.0f0
                               2.0f0 4.0f0 6.0f0]
        ps, st = Lux.setup(rng, m)

        @test Lux.parameterlength(m) == 2 * 2

        @test ps.bias == [0, 0]  # init_bias(2)
        @test ps.scale == [1, 1]  # init_scale(2)

        y, st_ = pullback(m, x, ps, st)[1]
        @test isapprox(y, [-1.22474 0 1.22474; -1.22474 0 1.22474], atol=1.0e-5)
        # julia> x
        #  2×3 Array{Float64,2}:
        #  1.0  3.0  5.0
        #  2.0  4.0  6.0
        #
        # mean of batch will be
        #  (1. + 3. + 5.) / 3 = 3
        #  (2. + 4. + 6.) / 3 = 4
        #
        # ∴ update rule with momentum:
        #  .1 * 3 + 0 = .3
        #  .1 * 4 + 0 = .4
        @test st_.running_mean ≈ reshape([0.3, 0.4], 2, 1)

        # julia> .1 .* var(x, dims = 2, corrected=false) .* (3 / 2).+ .9 .* [1., 1.]
        # 2×1 Array{Float64,2}:
        #  1.3
        #  1.3
        @test st_.running_var ≈
              0.1 .* var(x; dims=2, corrected=false) .* (3 / 2) .+ 0.9 .* [1.0, 1.0]

        st_ = Lux.testmode(st_)
        x′ = m(x, ps, st_)[1]
        @test isapprox(x′[1], (1 .- 0.3) / sqrt(1.3), atol=1.0e-5)

        @inferred m(x, ps, st)

        run_JET_tests(m, x, ps, st)

        test_gradient_correctness_fdm((x, ps) -> sum(first(m(x, ps, st))), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)
    end

    let m = BatchNorm(2; track_stats=false), x = [1.0f0 3.0f0 5.0f0; 2.0f0 4.0f0 6.0f0]
        ps, st = Lux.setup(rng, m)
        @inferred m(x, ps, st)
        run_JET_tests(m, x, ps, st)

        test_gradient_correctness_fdm((x, ps) -> sum(first(m(x, ps, st))), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)
    end

    # with activation function
    let m = BatchNorm(2, sigmoid), x = [1.0f0 3.0f0 5.0f0
                                        2.0f0 4.0f0 6.0f0]
        ps, st = Lux.setup(rng, m)
        st = Lux.testmode(st)
        y, st_ = m(x, ps, st)
        @test isapprox(y,
                       sigmoid.((x .- st_.running_mean) ./
                                sqrt.(st_.running_var .+ m.epsilon)), atol=1.0e-7)
        @inferred m(x, ps, st)
        run_JET_tests(m, x, ps, st)

        test_gradient_correctness_fdm((x, ps) -> sum(first(m(x, ps, st))), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)
    end

    let m = BatchNorm(2), x = reshape(Float32.(1:6), 3, 2, 1)
        ps, st = Lux.setup(rng, m)
        st = Lux.trainmode(st)
        @test_throws AssertionError m(x, ps, st)[1]
    end

    let m = BatchNorm(32), x = randn(Float32, 416, 416, 32, 1)
        ps, st = Lux.setup(rng, m)
        st = Lux.testmode(st)
        m(x, ps, st)
        @test (@allocated m(x, ps, st)) < 100_000_000
        @inferred m(x, ps, st)
        run_JET_tests(m, x, ps, st)
    end
end

@testset "GroupNorm" begin
    # begin tests
    squeeze(x) = dropdims(x; dims=tuple(findall(size(x) .== 1)...)) # To remove all singular dimensions

    let m = GroupNorm(4, 2; track_stats=true), sizes = (3, 4, 2),
        x = reshape(collect(1:prod(sizes)), sizes)

        @test Lux.parameterlength(m) == 2 * 4
        x = Float32.(x)
        ps, st = Lux.setup(rng, m)
        @test ps.bias == [0, 0, 0, 0]   # init_bias(32)
        @test ps.scale == [1, 1, 1, 1]  # init_scale(32)

        y, st_ = pullback(m, x, ps, st)[1]

        # julia> x
        # [:, :, 1]  =
        # 1.0  4.0  7.0  10.0
        # 2.0  5.0  8.0  11.0
        # 3.0  6.0  9.0  12.0
        #
        # [:, :, 2] =
        # 13.0  16.0  19.0  22.0
        # 14.0  17.0  20.0  23.0
        # 15.0  18.0  21.0  24.0
        #
        # mean will be
        # (1. + 2. + 3. + 4. + 5. + 6.) / 6 = 3.5
        # (7. + 8. + 9. + 10. + 11. + 12.) / 6 = 9.5
        #
        # (13. + 14. + 15. + 16. + 17. + 18.) / 6 = 15.5
        # (19. + 20. + 21. + 22. + 23. + 24.) / 6 = 21.5
        #
        # mean =
        # 3.5   15.5
        # 9.5   21.5
        #
        # ∴ update rule with momentum:
        # (1. - .1) * 0 + .1 * (3.5 + 15.5) / 2 = 0.95
        # (1. - .1) * 0 + .1 * (9.5 + 21.5) / 2 = 1.55
        @test st_.running_mean ≈ [0.95, 1.55]
        n = prod(size(x)) ÷ m.groups ÷ size(x)[end]
        corr = n / (n - 1)
        z = reshape(x, 3, 2, 2, 2)
        variance = var(z; dims=(1, 2), corrected=false)
        @test st_.running_var ≈ 0.1 * corr * vec(mean(variance; dims=4)) .+ 0.9 * 1

        st__ = Lux.testmode(st_)
        y, st__ = m(x, ps, st__)
        out = (z .- reshape(st_.running_mean, 1, 1, 2, 1)) ./
              sqrt.(reshape(st_.running_var, 1, 1, 2, 1) .+ 1.0f-5)
        @test y≈reshape(out, size(x)) atol=1.0e-5

        @inferred m(x, ps, st)
        run_JET_tests(m, x, ps, st)
        test_gradient_correctness_fdm(ps -> sum(first(m(x, ps, st))), ps; atol=1.0f-3,
                                      rtol=1.0f-3)
    end
end


text:
@testset "temporal" begin

    hod = HoD()
    @test hod isa Transform
    @test cardinality(hod) == OneToOne()

    @testset "Basic" begin
        x = collect(DateTime(2020, 1, 1, 9, 0):Hour(1):DateTime(2020, 5, 7, 9, 0))
        # Expected result is an hour a day starting and ending on the 9th hour inclusive,
        # with 126 full days in the middle
        expected = [9:23..., repeat(0:23, 126)..., 0:9...]

        @test FeatureTransforms.apply(x, hod) == expected
        @test hod(x) == expected

        @testset "StepRange" begin
            x = DateTime(2020, 1, 1, 9, 0):Hour(1):DateTime(2020, 5, 7, 9, 0)
            @test FeatureTransforms.apply(x, hod) == expected
            @test hod(x) == expected
        end

        @testset "DST" begin
            x = ZonedDateTime(2020, 3, 7, 9, 0, tz"America/New_York"):Hour(1):ZonedDateTime(2020, 3, 8, 9, 0, tz"America/New_York")

            # expected result skips the DST transition hour of 2
            expected_dst = [9:23..., 0, 1, 3:9...]

            @test FeatureTransforms.apply(x, hod) == expected_dst
            @test hod(x) == expected_dst
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MagicVLSI"
version = v"8.3.160"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Keno/magic.git", "8ab1743be8f6e0d7ab95e0c2e86688379fca1aad")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/magic/
apk add tcsh
apk add tcl
update_configure_scripts
cd scripts
autoreconf
cd ..
./configure --with-wish=wish8.6 --prefix=${prefix} --build=${MACHTYPE} --host=${target} LDFLAGS="-Wl,-rpath-link=/opt/${target}/${target}/lib64"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(p -> (Sys.islinux(p) && libc(p) != "musl") || Sys.isfreebsd(p), supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("tclmagic", :tclmagic, "lib/magic/tcl")
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Tcl_jll", uuid="9fc9b8dd-f768-5d0d-bb0b-e65961e00eb6"))
    Dependency(PackageSpec(name="GLU_jll", uuid="bd17208b-e95e-5925-bf81-e2f59b3e5c61"))
    Dependency(PackageSpec(name="Cairo_jll"))
    Dependency(PackageSpec(name="Tk_jll"))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Testing objects and Layout" begin
    layout(3,1,"GA")
    obj1=BasisBlade(2)
    answer1=e1
    @test obj1==answer1
    obj2=[BasisBlade(1),BasisBlade(2),BasisBlade(3),BasisBlade(4),BasisBlade(5)]
    answer2=[id,e1,e2,e3,e4]
    @test obj2==answer2
    obj3=BasisBlade(16)
    answer3=e1234
    @test obj3==answer3

    layout(5,0,"Projective")
    obj1=BasisBlade(2)
    answer1=e1
    @test obj1==answer1
    obj2=[BasisBlade(1),BasisBlade(2),BasisBlade(3),BasisBlade(4),BasisBlade(5),BasisBlade(6) ]
    answer2=[id,e1,e2,e3,e4,e5]
    @test obj2==answer2
    obj3=BasisBlade(32)
    answer3=e12345
    @test obj3==answer3

    layout(3,1,"Conformal")
    obj1=BasisBlade(2)
    answer1=e1
    @test obj1==answer1
    obj2=[BasisBlade(1),BasisBlade(2),BasisBlade(3),BasisBlade(4),BasisBlade(5)]
    answer2=[id,e1,e2,ep,en]
    @test obj2==answer2
    obj3=BasisBlade(16)
    answer3=e12pn
    @test obj3==answer3
    obj4=e0
    obj5=e∞
    @test obj4==-0.5*ep+0.5*en
    @test obj5==1.0*ep+1.0*en

end


text:
using Mustache
using Base.Test

tpl = mt"a:{{x}} b:{{{y}}}"

x, y = "ex", "why"
d = Dict("x"=>"ex", "y"=>"why")
type ThrowAway
    x
    y
end

@test render(tpl, Main) == "a:ex b:why"
@test render(tpl, d) == "a:ex b:why"
@test render(tpl, ThrowAway(x,y)) == "a:ex b:why"

## triple quoted
tpl = mt"""a:{{x}} b:{{y}}"""

@test render(tpl, Main) == "a:ex b:why"
@test render(tpl, d) == "a:ex b:why"
@test render(tpl, ThrowAway(x,y)) == "a:ex b:why"

## conditional
tpl = "{{#b}}this doesn't show{{/b}}{{#a}}this does show{{/a}}"
@test render(tpl, Dict("a" => 1)) == "this does show"

## dict using symbols
d = Dict(:a => x, :b => y)
tpl = "a:{{:a}} b:{{:b}}"
@test render(tpl, d) == "a:ex b:why"

## keyword args
@test render(tpl, a="ex", b="why") == "a:ex b:why"

## unicode
module TMP
α = 1
β = 2
end
@test render("{{α}} + 1 = {{β}}", TMP) == "1 + 1 = 2"
@test render("{{:β}}", β=2) == "2"
@test render("α - {{:β}}", β=2) == "α - 2"
@test render("{{:α}}", α="β") == "β"

## {{.}} test
tpl = mt"{{#:vec}}{{.}} {{/:vec}}"
@test render(tpl, vec=[1,2,3]) == "1 2 3 "

## test of function, see http://mustache.github.io/mustache.5.html (Lambdas)

tpl = mt"""{{#wrapped}}
  {{name}} is awesome.
{{/wrapped}}
"""

d = Dict()
d["name"] =  "Willy"
d["wrapped"] = function()
    function(text, render)
        "" * render(text) * ""
    end
    end

@test Mustache.render(tpl, d) == "\n  Willy is awesome.\n\n" #?? extra \n??

## Test of using Dict in {{#}}/{{/}} things
tpl = mt"{{#:d}}{{x}} and {{y}}{{/:d}}"
d = Dict(); d["x"] = "salt"; d["y"] = "pepper"
@test Mustache.render(tpl, d=d) == "salt and pepper"


text:
"""
Description

      The Watchmaker.plot() method uses the Winston plotting package to plot one or more time series
      on a single plot

Usage

      using MarketData
      plot(cl, op)
      plot(percentchange(CAT["Close"]), title="Caterpillar daily returns")

Method Signature(s)

      plot(ts::TimeSeries.TimeArray...; title::ASCIIString="")

Details

      The Watchmaker.plot() method takes a variable number of TimeSeries.TimeArray objects (maximum of six) and plots them
      on a single plot. The color scheme is based on the default_colors constant that is exported in this package.

References

      https://github.com/nolta/Winston.jl

See Also

      uplot()

"""

function  plot(ts::TimeSeries.TimeArray...; title::ASCIIString="")
    # stop the madness at 6
    length(ts) > 6 ? error("maximum of 6 time series at a time") : n = length(ts)

    # the legend names are taken from the series meta data
    legend = Array(ASCIIString, n)

    # the last time series is outside the loop so just fill it in here
    isa(ts[n].meta, AbstractString) ? legend[n] = ts[n].meta : legend[n] = string(n)

    for t in 1:n- 1
        # fill out the legend
        isa(ts[t].meta, AbstractString) ? legend[t] = ts[t].meta : legend[t] = string(t)

        # plot
        Winston.timeplot(ts[t].timestamp, ts[t].values[:,1], color=default_colors[t])
        Winston.hold(true)
    end

    p = Winston.timeplot(ts[n].timestamp, ts[n].values[:,1], color=default_colors[n], title=title)
    Winston.hold(false)

    for s in 1:n -1
        lab = Winston.PlotLabel(0.1, 1-0.05s, legend[s], color=default_colors[s])
        Winston.add(p, lab)
    end

    lab = Winston.PlotLabel(0.1, 1-0.05n, legend[n], color=default_colors[n])
    return Winston.add(p, lab)

end


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

include("model_ErrorResponse.jl")
include("model_InfoField.jl")
include("model_MeterInfo.jl")
include("model_OfferTermInfo.jl")
include("model_RateCardQueryParameters.jl")
include("model_ResourceRateCardInfo.jl")
include("model_UsageAggregation.jl")
include("model_UsageAggregationListResult.jl")
include("model_UsageSample.jl")
include("model_MonetaryCommitment.jl")
include("model_MonetaryCredit.jl")
include("model_RecurringCharge.jl")


text:
# This list was obtained through AWK with Knitro 11.0 (and then pared down) with the following two commands:
# grep "#define" include/ktr.h | awk '{ printf("const %s = Int32(%s)\n",$2,$3) }'
# grep "#  define" include/ktr.h | awk '{ printf("const %s = Int32(%s)\n",$3,$4) }'
# For version up until 10.3, use following commands :
# grep "#define" knitro.h | awk '{ printf("const %s = Int32(%s)\n",$2,$3) }'
# grep "#  define" knitro.h | awk '{ printf("const %s = Int32(%s)\n",$3,$4) }'

export
    KTR_INFBOUND,
    KTR_PARAMTYPE_INTEGER,
    KTR_PARAMTYPE_FLOAT,
    KTR_PARAMTYPE_STRING,
    KTR_OBJGOAL_MINIMIZE,
    KTR_OBJGOAL_MAXIMIZE,
    KTR_OBJTYPE_CONSTANT,
    KTR_OBJTYPE_GENERAL,
    KTR_OBJTYPE_LINEAR,
    KTR_OBJTYPE_QUADRATIC,
    KTR_CONTYPE_GENERAL,
    KTR_CONTYPE_LINEAR,
    KTR_CONTYPE_QUADRATIC,
    KTR_RESTYPE_GENERAL,
    KTR_RESTYPE_LINEAR,
    KTR_VARTYPE_CONTINUOUS,
    KTR_VARTYPE_INTEGER,
    KTR_VARTYPE_BINARY,
    KTR_FNTYPE_UNCERTAIN,
    KTR_FNTYPE_CONVEX,
    KTR_FNTYPE_NONCONVEX,
    KTR_LINEARVAR_NO,
    KTR_LINEARVAR_YES,
    KTR_RC_BEGINEND,
    KTR_RC_EVALFC,
    KTR_RC_EVALGA,
    KTR_RC_EVALH,
    KTR_RC_EVALX0,
    KTR_RC_FINISHED,
    KTR_RC_NEWPOINT,
    KTR_RC_EVALHV,
    KTR_RC_EVALH_NO_F,
    KTR_RC_EVALHV_NO_F,
    KTR_RC_NODE,
    KTR_RC_MSPROCESS,
    KTR_RC_OPTIMAL_OR_SATISFACTORY,
    KTR_RC_OPTIMAL,
    KTR_RC_NEAR_OPT,
    KTR_RC_FEAS_XTOL,
    KTR_RC_FEAS_NO_IMPROVE,
    KTR_RC_FEAS_FTOL,
    KTR_RC_INFEASIBLE,
    KTR_RC_INFEAS_XTOL,
    KTR_RC_INFEAS_NO_IMPROVE,
    KTR_RC_INFEAS_MULTISTART,
    KTR_RC_INFEAS_CON_BOUNDS,
    KTR_RC_INFEAS_VAR_BOUNDS,
    KTR_RC_UNBOUNDED,
    KTR_RC_ITER_LIMIT_FEAS,
    KTR_RC_ITER_LIMIT,
    KTR_RC_TIME_LIMIT_FEAS,
    KTR_RC_TIME_LIMIT,
    KTR_RC_FEVAL_LIMIT_FEAS,
    KTR_RC_MIP_EXH_FEAS,
    KTR_RC_MIP_EXH,
    KTR_RC_MIP_TERM_FEAS,
    KTR_RC_MIP_FEAS_TERM,
    KTR_RC_MIP_SOLVE_LIMIT_FEAS,
    KTR_RC_MIP_SOLVE_LIMIT,
    KTR_RC_MIP_NODE_LIMIT_FEAS,
    KTR_RC_MIP_NODE_LIMIT,
    KTR_RC_ITER_LIMIT_INFEAS,
    KTR_RC_TIME_LIMIT_INFEAS,
    KTR_RC_FEVAL_LIMIT_INFEAS,
    KTR_RC_MIP_EXH_INFEAS,
    KTR_RC_MIP_SOLVE_LIMIT_INFEAS,
    KTR_RC_MIP_NODE_LIMIT_INFEAS,
    KTR_RC_CALLBACK_ERR,
    KTR_RC_LP_SOLVER_ERR,
    KTR_RC_EVAL_ERR,
    KTR_RC_OUT_OF_MEMORY,
    KTR_RC_USER_TERMINATION,
    KTR_RC_OPEN_FILE_ERR,
    KTR_RC_BAD_N_OR_F,
    KTR_RC_BAD_CONSTRAINT,
    KTR_RC_BAD_JACOBIAN,
    KTR_RC_BAD_HESSIAN,
    KTR_RC_BAD_CON_INDEX,
    KTR_RC_BAD_JAC_INDEX,
    KTR_RC_BAD_HESS_INDEX,
    KTR_RC_BAD_CON_BOUNDS,
    KTR_RC_BAD_VAR_BOUNDS,
    KTR_RC_ILLEGAL_CALL,
    KTR_RC_BAD_KCPTR,
    KTR_RC_NULL_POINTER,
    KTR_RC_BAD_INIT_VALUE,
    KTR_RC_NEWPOINT_HALT,
    KTR_RC_BAD_LICENSE,
    KTR_RC_BAD_PARAMINPUT,
    KTR_RC_LINEAR_SOLVER_ERR,
    KTR_RC_DERIV_CHECK_FAILED,
    KTR_RC_DERIV_CHECK_TERMINATE,
    KTR_RC_OVERFLOW_ERR,
    KTR_RC_INTERNAL_ERROR,
    KTR_PARAM_NEWPOINT,
    KTR_PARAM_HONORBNDS,
    KTR_PARAM_ALGORITHM,
    KTR_PARAM_ALG,
    KTR_PARAM_BAR_MURULE,
    KTR_PARAM_BARRULE,
    KTR_PARAM_BAR_FEASIBLE,
    KTR_PARAM_FEASIBLE,
    KTR_PARAM_GRADOPT,
    KTR_PARAM_HESSOPT,
    KTR_PARAM_BAR_INITPT,
    KTR_PARAM_ACT_LPSOLVER,
    KTR_PARAM_LPSOLVER,
    KTR_PARAM_CG_MAXIT,
    KTR_PARAM_MAXCGIT,
    KTR_PARAM_MAXIT,
    KTR_PARAM_OUTLEV,
    KTR_PARAM_OUTMODE,
    KTR_PARAM_SCALE,
    KTR_PARAM_SHIFTINIT,
    KTR_PARAM_SOC,
    KTR_PARAM_DELTA,
    KTR_PARAM_BAR_FEASMODETOL,
    KTR_PARAM_FEASMODETOL,
    KTR_PARAM_FEASTOL,
    KTR_PARAM_FEASTOLABS,
    KTR_PARAM_MAXTIMECPU,
    KTR_PARAM_BAR_INITMU,
    KTR_PARAM_MU,
    KTR_PARAM_OBJRANGE,
    KTR_PARAM_OPTTOL,
    KTR_PARAM_OPTTOLABS,
    KTR_PARAM_LINSOLVER_PIVOTTOL,
    KTR_PARAM_PIVOT,
    KTR_PARAM_XTOL,
    KTR_PARAM_DEBUG,
    KTR_PARAM_MULTISTART,
    KTR_PARAM_MSENABLE,
    KTR_PARAM_MSMAXSOLVES,
    KTR_PARAM_MSMAXBNDRANGE,
    KTR_PARAM_MSMAXTIMECPU,
    KTR_PARAM_MSMAXTIMEREAL,
    KTR_PARAM_LMSIZE,
    KTR_PARAM_BAR_MAXCROSSIT,
    KTR_PARAM_MAXCROSSIT,
    KTR_PARAM_MAXTIMEREAL,
    KTR_PARAM_CG_PRECOND,
    KTR_PARAM_BLASOPTION,
    KTR_PARAM_BAR_MAXREFACTOR,
    KTR_PARAM_BAR_MAXBACKTRACK,
    KTR_PARAM_LINESEARCH_MAXTRIALS,
    KTR_PARAM_BLASOPTIONLIB,
    KTR_PARAM_OUTAPPEND,
    KTR_PARAM_OUTDIR,
    KTR_PARAM_CPLEXLIB,
    KTR_PARAM_BAR_PENRULE,
    KTR_PARAM_BAR_PENCONS,
    KTR_PARAM_MSNUMTOSAVE,
    KTR_PARAM_MSSAVETOL,
    KTR_PARAM_PRESOLVEDEBUG,
    KTR_PARAM_MSTERMINATE,
    KTR_PARAM_MSSTARTPTRANGE,
    KTR_PARAM_INFEASTOL,
    KTR_PARAM_LINSOLVER,
    KTR_PARAM_BAR_DIRECTINTERVAL,
    KTR_PARAM_PRESOLVE,
    KTR_PARAM_PRESOLVE_TOL,
    KTR_PARAM_BAR_SWITCHRULE,
    KTR_PARAM_HESSIAN_NO_F,
    KTR_PARAM_MA_TERMINATE,
    KTR_PARAM_MA_MAXTIMECPU,
    KTR_PARAM_MA_MAXTIMEREAL,
    KTR_PARAM_MSSEED,
    KTR_PARAM_MA_OUTSUB,
    KTR_PARAM_MS_OUTSUB,
    KTR_PARAM_XPRESSLIB,
    KTR_PARAM_TUNER,
    KTR_PARAM_TUNER_OPTIONSFILE,
    KTR_PARAM_TUNER_MAXTIMECPU,
    KTR_PARAM_TUNER_MAXTIMEREAL,
    KTR_PARAM_TUNER_OUTSUB,
    KTR_PARAM_TUNER_TERMINATE,
    KTR_PARAM_LINSOLVER_OOC,
    KTR_PARAM_BAR_RELAXCONS,
    KTR_PARAM_MSDETERMINISTIC,
    KTR_PARAM_BAR_REFINEMENT,
    KTR_PARAM_DERIVCHECK,
    KTR_PARAM_DERIVCHECK_TYPE,
    KTR_PARAM_DERIVCHECK_TOL,
    KTR_PARAM_LINSOLVER_INEXACT,
    KTR_PARAM_LINSOLVER_INEXACTTOL,
    KTR_PARAM_MAXFEVALS,
    KTR_PARAM_FSTOPVAL,
    KTR_PARAM_DATACHECK,
    KTR_PARAM_DERIVCHECK_TERMINATE,
    KTR_PARAM_BAR_WATCHDOG,
    KTR_PARAM_FTOL,
    KTR_PARAM_FTOL_ITERS,
    KTR_PARAM_ACT_QPALG,
    KTR_PARAM_BAR_INITPI_MPEC,
    KTR_PARAM_XTOL_ITERS,
    KTR_PARAM_LINESEARCH,
    KTR_PARAM_OUT_CSVINFO,
    KTR_PARAM_INITPENALTY,
    KTR_PARAM_ACT_LPFEASTOL,
    KTR_PARAM_CG_STOPTOL,
    KTR_PARAM_RESTARTS,
    KTR_PARAM_RESTARTS_MAXIT,
    KTR_PARAM_BAR_SLACKBOUNDPUSH,
    KTR_PARAM_CG_PMEM,
    KTR_PARAM_BAR_SWITCHOBJ,
    KTR_PARAM_OUTNAME,
    KTR_PARAM_OUT_CSVNAME,
    KTR_PARAM_ACT_PARAMETRIC,
    KTR_PARAM_ACT_LPDUMPMPS,
    KTR_PARAM_ACT_LPALG,
    KTR_PARAM_ACT_LPPRESOLVE,
    KTR_PARAM_ACT_LPPENALTY,
    KTR_PARAM_CONVEX_TUNE,
    KTR_PARAM_OUT_HINTS,
    KTR_PARAM_MIP_METHOD,
    KTR_PARAM_MIP_BRANCHRULE,
    KTR_PARAM_MIP_SELECTRULE,
    KTR_PARAM_MIP_INTGAPABS,
    KTR_PARAM_MIP_INTGAPREL,
    KTR_PARAM_MIP_MAXTIMECPU,
    KTR_PARAM_MIP_MAXTIMEREAL,
    KTR_PARAM_MIP_MAXSOLVES,
    KTR_PARAM_MIP_INTEGERTOL,
    KTR_PARAM_MIP_OUTLEVEL,
    KTR_PARAM_MIP_OUTINTERVAL,
    KTR_PARAM_MIP_OUTSUB,
    KTR_PARAM_MIP_DEBUG,
    KTR_PARAM_MIP_IMPLICATNS,
    KTR_PARAM_MIP_GUB_BRANCH,
    KTR_PARAM_MIP_KNAPSACK,
    KTR_PARAM_MIP_ROUNDING,
    KTR_PARAM_MIP_ROOTALG,
    KTR_PARAM_MIP_LPALG,
    KTR_PARAM_MIP_TERMINATE,
    KTR_PARAM_MIP_MAXNODES,
    KTR_PARAM_MIP_HEURISTIC,
    KTR_PARAM_MIP_HEUR_MAXIT,
    KTR_PARAM_MIP_HEUR_MAXTIMECPU,
    KTR_PARAM_MIP_HEUR_MAXTIMEREAL,
    KTR_PARAM_MIP_PSEUDOINIT,
    KTR_PARAM_MIP_STRONG_MAXIT,
    KTR_PARAM_MIP_STRONG_CANDLIM,
    KTR_PARAM_MIP_STRONG_LEVEL,
    KTR_PARAM_MIP_INTVAR_STRATEGY,
    KTR_PARAM_MIP_RELAXABLE,
    KTR_PARAM_MIP_NODEALG,
    KTR_PARAM_MIP_HEUR_TERMINATE,
    KTR_PARAM_MIP_SELECTDIR,
    KTR_PARAM_PAR_NUMTHREADS,
    KTR_PARAM_PAR_CONCURRENT_EVALS,
    KTR_PARAM_PAR_BLASNUMTHREADS,
    KTR_PARAM_PAR_LSNUMTHREADS,
    KTR_PARAM_PAR_MSNUMTHREADS,
    KTR_NEWPOINT_NONE,
    KTR_NEWPOINT_SAVEONE,
    KTR_NEWPOINT_SAVEALL,
    KTR_NEWPOINT_USER,
    KTR_HONORBNDS_NO,
    KTR_HONORBNDS_ALWAYS,
    KTR_HONORBNDS_INITPT,
    KTR_ALG_AUTOMATIC,
    KTR_ALG_AUTO,
    KTR_ALG_BAR_DIRECT,
    KTR_ALG_BAR_CG,
    KTR_ALG_ACT_CG,
    KTR_ALG_IPDIRECT,
    KTR_ALG_IPCG,
    KTR_ALG_ACTIVE,
    KTR_ALG_ACT_SQP,
    KTR_ALG_MULTI,
    KTR_BAR_MURULE_AUTOMATIC,
    KTR_BAR_MURULE_AUTO,
    KTR_BAR_MURULE_MONOTONE,
    KTR_BAR_MURULE_ADAPTIVE,
    KTR_BAR_MURULE_PROBING,
    KTR_BAR_MURULE_DAMPMPC,
    KTR_BAR_MURULE_FULLMPC,
    KTR_BAR_MURULE_QUALITY,
    KTR_BAR_FEASIBLE_NO,
    KTR_BAR_FEASIBLE_STAY,
    KTR_BAR_FEASIBLE_GET,
    KTR_BAR_FEASIBLE_GET_STAY,
    KTR_FEASIBLE_NO,
    KTR_FEASIBLE_ALWAYS,
    KTR_GRADOPT_EXACT,
    KTR_GRADOPT_FORWARD,
    KTR_GRADOPT_CENTRAL,
    KTR_GRADOPT_USER_FORWARD,
    KTR_GRADOPT_USER_CENTRAL,
    KTR_HESSOPT_AUTO,
    KTR_HESSOPT_EXACT,
    KTR_HESSOPT_BFGS,
    KTR_HESSOPT_SR1,
    KTR_HESSOPT_FINITE_DIFF,
    KTR_HESSOPT_PRODUCT_FINDIFF,
    KTR_HESSOPT_PRODUCT,
    KTR_HESSOPT_LBFGS,
    KTR_HESSOPT_GAUSS_NEWTON,
    KTR_BAR_INITPT_AUTO,
    KTR_BAR_INITPT_YES,
    KTR_BAR_INITPT_STRAT1,
    KTR_BAR_INITPT_NO,
    KTR_BAR_INITPT_STRAT2,
    KTR_BAR_INITPT_STRAT3,
    KTR_ACT_LPSOLVER_INTERNAL,
    KTR_ACT_LPSOLVER_CPLEX,
    KTR_ACT_LPSOLVER_XPRESS,
    KTR_LP_INTERNAL,
    KTR_LP_CPLEX,
    KTR_LP_XPRESS,
    KTR_OUTLEV_NONE,
    KTR_OUTLEV_SUMMARY,
    KTR_OUTLEV_MAJORIT10,
    KTR_OUTLEV_ITER_10,
    KTR_OUTLEV_MAJORIT,
    KTR_OUTLEV_ITER,
    KTR_OUTLEV_ALLIT,
    KTR_OUTLEV_ITER_VERBOSE,
    KTR_OUTLEV_ALLIT_X,
    KTR_OUTLEV_ITER_X,
    KTR_OUTLEV_ALL,
    KTR_OUTMODE_SCREEN,
    KTR_OUTMODE_FILE,
    KTR_OUTMODE_BOTH,
    KTR_SCALE_NEVER,
    KTR_SCALE_NO,
    KTR_SCALE_ALLOW,
    KTR_SCALE_USER_INTERNAL,
    KTR_SCALE_USER_NONE,
    KTR_SCALE_INTERNAL,
    KTR_SOC_NO,
    KTR_SOC_MAYBE,
    KTR_SOC_YES,
    KTR_DEBUG_NONE,
    KTR_DEBUG_PROBLEM,
    KTR_DEBUG_EXECUTION,
    KTR_MULTISTART_NO,
    KTR_MULTISTART_YES,
    KTR_CG_PRECOND_NONE,
    KTR_CG_PRECOND_CHOL,
    KTR_BLASOPTION_KNITRO,
    KTR_BLASOPTION_INTEL,
    KTR_BLASOPTION_DYNAMIC,
    KTR_OUTAPPEND_NO,
    KTR_OUTAPPEND_YES,
    KTR_BAR_PENRULE_AUTO,
    KTR_BAR_PENRULE_SINGLE,
    KTR_BAR_PENRULE_FLEX,
    KTR_BAR_PENCONS_AUTO,
    KTR_BAR_PENCONS_NONE,
    KTR_BAR_PENCONS_ALL,
    KTR_BAR_PENCONS_EQUALITIES,
    KTR_BAR_PENCONS_INFEAS,
    KTR_PRESOLVEDBG_NONE,
    KTR_PRESOLVEDBG_BASIC,
    KTR_PRESOLVEDBG_VERBOSE,
    KTR_MSTERMINATE_MAXSOLVES,
    KTR_MSTERMINATE_OPTIMAL,
    KTR_MSTERMINATE_FEASIBLE,
    KTR_MSTERMINATE_ANY,
    KTR_LINSOLVER_AUTO,
    KTR_LINSOLVER_INTERNAL,
    KTR_LINSOLVER_HYBRID,
    KTR_LINSOLVER_DENSEQR,
    KTR_LINSOLVER_MA27,
    KTR_LINSOLVER_MA57,
    KTR_LINSOLVER_MKLPARDISO,
    KTR_LINSOLVER_MA97,
    KTR_LINSOLVER_MA86,
    KTR_PRESOLVE_NONE,
    KTR_PRESOLVE_BASIC,
    KTR_PRESOLVE_ADVANCED,
    KTR_BAR_SWITCHRULE_AUTO,
    KTR_BAR_SWITCHRULE_NEVER,
    KTR_BAR_SWITCHRULE_LEVEL1,
    KTR_BAR_SWITCHRULE_LEVEL2,
    KTR_HESSIAN_NO_F_FORBID,
    KTR_HESSIAN_NO_F_ALLOW,
    KTR_MA_TERMINATE_ALL,
    KTR_MA_TERMINATE_OPTIMAL,
    KTR_MA_TERMINATE_FEASIBLE,
    KTR_MA_TERMINATE_ANY,
    KTR_MA_OUTSUB_NONE,
    KTR_MA_OUTSUB_YES,
    KTR_MS_OUTSUB_NONE,
    KTR_MS_OUTSUB_YES,
    KTR_TUNER_OFF,
    KTR_TUNER_ON,
    KTR_TUNER_OUTSUB_NONE,
    KTR_TUNER_OUTSUB_YES,
    KTR_TUNER_OUTSUB_SUMMARY,
    KTR_TUNER_OUTSUB_ALL,
    KTR_TUNER_TERMINATE_ALL,
    KTR_TUNER_TERMINATE_OPTIMAL,
    KTR_TUNER_TERMINATE_FEASIBLE,
    KTR_TUNER_TERMINATE_ANY,
    KTR_LINSOLVER_OOC_NO,
    KTR_LINSOLVER_OOC_MAYBE,
    KTR_LINSOLVER_OOC_YES,
    KTR_BAR_RELAXCONS_NONE,
    KTR_BAR_RELAXCONS_EQS,
    KTR_BAR_RELAXCONS_INEQS,
    KTR_BAR_RELAXCONS_ALL,
    KTR_MSDETERMINISTIC_NO,
    KTR_MSDETERMINISTIC_YES,
    KTR_BAR_REFINEMENT_NO,
    KTR_BAR_REFINEMENT_YES,
    KTR_DERIVCHECK_NONE,
    KTR_DERIVCHECK_FIRST,
    KTR_DERIVCHECK_SECOND,
    KTR_DERIVCHECK_ALL,
    KTR_DERIVCHECK_FORWARD,
    KTR_DERIVCHECK_CENTRAL,
    KTR_LINSOLVER_INEXACT_NO,
    KTR_LINSOLVER_INEXACT_YES,
    KTR_DATACHECK_NO,
    KTR_DATACHECK_YES,
    KTR_DERIVCHECK_STOPERROR,
    KTR_DERIVCHECK_STOPALWAYS,
    KTR_BAR_WATCHDOG_NO,
    KTR_BAR_WATCHDOG_YES,
    KTR_ACT_QPALG_AUTO,
    KTR_ACT_QPALG_BAR_DIRECT,
    KTR_ACT_QPALG_BAR_CG,
    KTR_ACT_QPALG_ACT_CG,
    KTR_LINESEARCH_AUTO,
    KTR_LINESEARCH_BACKTRACK,
    KTR_LINESEARCH_INTERPOLATE,
    KTR_OUT_CSVINFO_NO,
    KTR_OUT_CSVINFO_YES,
    KTR_BAR_SWITCHOBJ_NONE,
    KTR_BAR_SWITCHOBJ_SCALARPROX,
    KTR_BAR_SWITCHOBJ_DIAGPROX,
    KTR_ACT_PARAMETRIC_NO,
    KTR_ACT_PARAMETRIC_MAYBE,
    KTR_ACT_PARAMETRIC_YES,
    KTR_ACT_LPDUMPMPS_NO,
    KTR_ACT_LPDUMPMPS_YES,
    KTR_ACT_LPALG_DEFAULT,
    KTR_ACT_LPALG_PRIMAL,
    KTR_ACT_LPALG_DUAL,
    KTR_ACT_LPALG_BARRIER,
    KTR_ACT_LPPRESOLVE_OFF,
    KTR_ACT_LPPRESOLVE_ON,
    KTR_ACT_LPPENALTY_AUTO,
    KTR_ACT_LPPENALTY_ALL,
    KTR_ACT_LPPENALTY_NONLINEAR,
    KTR_ACT_LPPENALTY_DYNAMIC,
    KTR_CONVEX_TUNE_NO,
    KTR_CONVEX_TUNE_YES,
    KTR_OUT_HINTS_NO,
    KTR_OUT_HINTS_YES,
    KTR_MIP_METHOD_AUTO,
    KTR_MIP_METHOD_BB,
    KTR_MIP_METHOD_HQG,
    KTR_MIP_METHOD_MISQP,
    KTR_MIP_BRANCH_AUTO,
    KTR_MIP_BRANCH_MOSTFRAC,
    KTR_MIP_BRANCH_PSEUDOCOST,
    KTR_MIP_BRANCH_STRONG,
    KTR_MIP_SEL_AUTO,
    KTR_MIP_SEL_DEPTHFIRST,
    KTR_MIP_SEL_BESTBOUND,
    KTR_MIP_SEL_COMBO_1,
    KTR_MIP_OUTLEVEL_NONE,
    KTR_MIP_OUTLEVEL_ITERS,
    KTR_MIP_OUTLEVEL_ITERSTIME,
    KTR_MIP_OUTLEVEL_ROOT,
    KTR_MIP_OUTSUB_NONE,
    KTR_MIP_OUTSUB_YES,
    KTR_MIP_OUTSUB_YESPROB,
    KTR_MIP_DEBUG_NONE,
    KTR_MIP_DEBUG_ALL,
    KTR_MIP_IMPLICATNS_NO,
    KTR_MIP_IMPLICATNS_YES,
    KTR_MIP_GUB_BRANCH_NO,
    KTR_MIP_GUB_BRANCH_YES,
    KTR_MIP_KNAPSACK_NO,
    KTR_MIP_KNAPSACK_INEQ,
    KTR_MIP_KNAPSACK_INEQ_EQ,
    KTR_MIP_ROUND_AUTO,
    KTR_MIP_ROUND_NONE,
    KTR_MIP_ROUND_HEURISTIC,
    KTR_MIP_ROUND_NLP_SOME,
    KTR_MIP_ROUND_NLP_ALWAYS,
    KTR_MIP_ROOTALG_AUTO,
    KTR_MIP_ROOTALG_BAR_DIRECT,
    KTR_MIP_ROOTALG_BAR_CG,
    KTR_MIP_ROOTALG_ACT_CG,
    KTR_MIP_ROOTALG_ACT_SQP,
    KTR_MIP_ROOTALG_MULTI,
    KTR_MIP_LPALG_AUTO,
    KTR_MIP_LPALG_BAR_DIRECT,
    KTR_MIP_LPALG_BAR_CG,
    KTR_MIP_LPALG_ACT_CG,
    KTR_MIP_TERMINATE_OPTIMAL,
    KTR_MIP_TERMINATE_FEASIBLE,
    KTR_MIP_HEURISTIC_AUTO,
    KTR_MIP_HEURISTIC_NONE,
    KTR_MIP_HEURISTIC_FEASPUMP,
    KTR_MIP_HEURISTIC_MPEC,
    KTR_MIP_PSEUDOINIT_AUTO,
    KTR_MIP_PSEUDOINIT_AVE,
    KTR_MIP_PSEUDOINIT_STRONG,
    KTR_MIP_INTVAR_STRATEGY_NONE,
    KTR_MIP_INTVAR_STRATEGY_RELAX,
    KTR_MIP_INTVAR_STRATEGY_MPEC,
    KTR_MIP_RELAXABLE_NONE,
    KTR_MIP_RELAXABLE_ALL,
    KTR_MIP_NODEALG_AUTO,
    KTR_MIP_NODEALG_BAR_DIRECT,
    KTR_MIP_NODEALG_BAR_CG,
    KTR_MIP_NODEALG_ACT_CG,
    KTR_MIP_NODEALG_ACT_SQP,
    KTR_MIP_NODEALG_MULTI,
    KTR_MIP_HEUR_TERMINATE_FEASIBLE,
    KTR_MIP_HEUR_TERMINATE_LIMIT,
    KTR_MIP_SELECTDIR_DOWN,
    KTR_MIP_SELECTDIR_UP,
    KTR_PAR_CONCURRENT_EVALS_NO,
    KTR_PAR_CONCURRENT_EVALS_YES,
    KTR_PAR_MSNUMTHREADS_AUTO

const KTR_INFBOUND = Float64(1.0e20)
const KTR_PARAMTYPE_INTEGER = Int32(0)
const KTR_PARAMTYPE_FLOAT = Int32(1)
const KTR_PARAMTYPE_STRING = Int32(2)
const KTR_OBJGOAL_MINIMIZE = Int32(0)
const KTR_OBJGOAL_MAXIMIZE = Int32(1)
const KTR_OBJTYPE_CONSTANT = Int32(-1)
const KTR_OBJTYPE_GENERAL = Int32(0)
const KTR_OBJTYPE_LINEAR = Int32(1)
const KTR_OBJTYPE_QUADRATIC = Int32(2)
const KTR_CONTYPE_GENERAL = Int32(0)
const KTR_CONTYPE_LINEAR = Int32(1)
const KTR_CONTYPE_QUADRATIC = Int32(2)
const KTR_RESTYPE_GENERAL = Int32(0)
const KTR_RESTYPE_LINEAR = Int32(1)
const KTR_VARTYPE_CONTINUOUS = Int32(0)
const KTR_VARTYPE_INTEGER = Int32(1)
const KTR_VARTYPE_BINARY = Int32(2)
const KTR_FNTYPE_UNCERTAIN = Int32(0)
const KTR_FNTYPE_CONVEX = Int32(1)
const KTR_FNTYPE_NONCONVEX = Int32(2)
const KTR_LINEARVAR_NO = Int32(0)
const KTR_LINEARVAR_YES = Int32(1)
const KTR_RC_BEGINEND = Int32(0)
const KTR_RC_EVALFC = Int32(1)
const KTR_RC_EVALGA = Int32(2)
const KTR_RC_EVALH = Int32(3)
const KTR_RC_EVALX0 = Int32(4)
const KTR_RC_FINISHED = Int32(5)
const KTR_RC_NEWPOINT = Int32(6)
const KTR_RC_EVALHV = Int32(7)
const KTR_RC_EVALH_NO_F = Int32(8)
const KTR_RC_EVALHV_NO_F = Int32(9)
const KTR_RC_NODE = Int32(10)
const KTR_RC_MSPROCESS = Int32(11)
const KTR_RC_OPTIMAL_OR_SATISFACTORY = Int32(0)
const KTR_RC_OPTIMAL = Int32(0)
const KTR_RC_NEAR_OPT = Int32(-100)
const KTR_RC_FEAS_XTOL = Int32(-101)
const KTR_RC_FEAS_NO_IMPROVE = Int32(-102)
const KTR_RC_FEAS_FTOL = Int32(-103)
const KTR_RC_INFEASIBLE = Int32(-200)
const KTR_RC_INFEAS_XTOL = Int32(-201)
const KTR_RC_INFEAS_NO_IMPROVE = Int32(-202)
const KTR_RC_INFEAS_MULTISTART = Int32(-203)
const KTR_RC_INFEAS_CON_BOUNDS = Int32(-204)
const KTR_RC_INFEAS_VAR_BOUNDS = Int32(-205)
const KTR_RC_UNBOUNDED = Int32(-300)
const KTR_RC_ITER_LIMIT_FEAS = Int32(-400)
const KTR_RC_ITER_LIMIT = Int32(-400)
const KTR_RC_TIME_LIMIT_FEAS = Int32(-401)
const KTR_RC_TIME_LIMIT = Int32(-401)
const KTR_RC_FEVAL_LIMIT_FEAS = Int32(-402)
const KTR_RC_MIP_EXH_FEAS = Int32(-403)
const KTR_RC_MIP_EXH = Int32(-403)
const KTR_RC_MIP_TERM_FEAS = Int32(-404)
const KTR_RC_MIP_FEAS_TERM = Int32(-404)
const KTR_RC_MIP_SOLVE_LIMIT_FEAS = Int32(-405)
const KTR_RC_MIP_SOLVE_LIMIT = Int32(-405)
const KTR_RC_MIP_NODE_LIMIT_FEAS = Int32(-406)
const KTR_RC_MIP_NODE_LIMIT = Int32(-406)
const KTR_RC_ITER_LIMIT_INFEAS = Int32(-410)
const KTR_RC_TIME_LIMIT_INFEAS = Int32(-411)
const KTR_RC_FEVAL_LIMIT_INFEAS = Int32(-412)
const KTR_RC_MIP_EXH_INFEAS = Int32(-413)
const KTR_RC_MIP_SOLVE_LIMIT_INFEAS = Int32(-415)
const KTR_RC_MIP_NODE_LIMIT_INFEAS = Int32(-416)
const KTR_RC_CALLBACK_ERR = Int32(-500)
const KTR_RC_LP_SOLVER_ERR = Int32(-501)
const KTR_RC_EVAL_ERR = Int32(-502)
const KTR_RC_OUT_OF_MEMORY = Int32(-503)
const KTR_RC_USER_TERMINATION = Int32(-504)
const KTR_RC_OPEN_FILE_ERR = Int32(-505)
const KTR_RC_BAD_N_OR_F = Int32(-506)
const KTR_RC_BAD_CONSTRAINT = Int32(-507)
const KTR_RC_BAD_JACOBIAN = Int32(-508)
const KTR_RC_BAD_HESSIAN = Int32(-509)
const KTR_RC_BAD_CON_INDEX = Int32(-510)
const KTR_RC_BAD_JAC_INDEX = Int32(-511)
const KTR_RC_BAD_HESS_INDEX = Int32(-512)
const KTR_RC_BAD_CON_BOUNDS = Int32(-513)
const KTR_RC_BAD_VAR_BOUNDS = Int32(-514)
const KTR_RC_ILLEGAL_CALL = Int32(-515)
const KTR_RC_BAD_KCPTR = Int32(-516)
const KTR_RC_NULL_POINTER = Int32(-517)
const KTR_RC_BAD_INIT_VALUE = Int32(-518)
const KTR_RC_NEWPOINT_HALT = Int32(-519)
const KTR_RC_BAD_LICENSE = Int32(-520)
const KTR_RC_BAD_PARAMINPUT = Int32(-521)
const KTR_RC_LINEAR_SOLVER_ERR = Int32(-522)
const KTR_RC_DERIV_CHECK_FAILED = Int32(-523)
const KTR_RC_DERIV_CHECK_TERMINATE = Int32(-524)
const KTR_RC_OVERFLOW_ERR = Int32(-525)
const KTR_RC_INTERNAL_ERROR = Int32(-600)
const KTR_PARAM_NEWPOINT = Int32(1001)
const KTR_PARAM_HONORBNDS = Int32(1002)
const KTR_PARAM_ALGORITHM = Int32(1003)
const KTR_PARAM_ALG = Int32(1003)
const KTR_PARAM_BAR_MURULE = Int32(1004)
const KTR_PARAM_BARRULE = Int32(1004)
const KTR_PARAM_BAR_FEASIBLE = Int32(1006)
const KTR_PARAM_FEASIBLE = Int32(1006)
const KTR_PARAM_GRADOPT = Int32(1007)
const KTR_PARAM_HESSOPT = Int32(1008)
const KTR_PARAM_BAR_INITPT = Int32(1009)
const KTR_PARAM_ACT_LPSOLVER = Int32(1012)
const KTR_PARAM_LPSOLVER = Int32(1012)
const KTR_PARAM_CG_MAXIT = Int32(1013)
const KTR_PARAM_MAXCGIT = Int32(1013)
const KTR_PARAM_MAXIT = Int32(1014)
const KTR_PARAM_OUTLEV = Int32(1015)
const KTR_PARAM_OUTMODE = Int32(1016)
const KTR_PARAM_SCALE = Int32(1017)
const KTR_PARAM_SHIFTINIT = Int32(1018)
const KTR_PARAM_SOC = Int32(1019)
const KTR_PARAM_DELTA = Int32(1020)
const KTR_PARAM_BAR_FEASMODETOL = Int32(1021)
const KTR_PARAM_FEASMODETOL = Int32(1021)
const KTR_PARAM_FEASTOL = Int32(1022)
const KTR_PARAM_FEASTOLABS = Int32(1023)
const KTR_PARAM_MAXTIMECPU = Int32(1024)
const KTR_PARAM_BAR_INITMU = Int32(1025)
const KTR_PARAM_MU = Int32(1025)
const KTR_PARAM_OBJRANGE = Int32(1026)
const KTR_PARAM_OPTTOL = Int32(1027)
const KTR_PARAM_OPTTOLABS = Int32(1028)
const KTR_PARAM_LINSOLVER_PIVOTTOL = Int32(1029)
const KTR_PARAM_PIVOT = Int32(1029)
const KTR_PARAM_XTOL = Int32(1030)
const KTR_PARAM_DEBUG = Int32(1031)
const KTR_PARAM_MULTISTART = Int32(1033)
const KTR_PARAM_MSENABLE = Int32(1033)
const KTR_PARAM_MSMAXSOLVES = Int32(1034)
const KTR_PARAM_MSMAXBNDRANGE = Int32(1035)
const KTR_PARAM_MSMAXTIMECPU = Int32(1036)
const KTR_PARAM_MSMAXTIMEREAL = Int32(1037)
const KTR_PARAM_LMSIZE = Int32(1038)
const KTR_PARAM_BAR_MAXCROSSIT = Int32(1039)
const KTR_PARAM_MAXCROSSIT = Int32(1039)
const KTR_PARAM_MAXTIMEREAL = Int32(1040)
const KTR_PARAM_CG_PRECOND = Int32(1041)
const KTR_PARAM_BLASOPTION = Int32(1042)
const KTR_PARAM_BAR_MAXREFACTOR = Int32(1043)
const KTR_PARAM_BAR_MAXBACKTRACK = Int32(1044)
const KTR_PARAM_LINESEARCH_MAXTRIALS = Int32(1044)
const KTR_PARAM_BLASOPTIONLIB = Int32(1045)
const KTR_PARAM_OUTAPPEND = Int32(1046)
const KTR_PARAM_OUTDIR = Int32(1047)
const KTR_PARAM_CPLEXLIB = Int32(1048)
const KTR_PARAM_BAR_PENRULE = Int32(1049)
const KTR_PARAM_BAR_PENCONS = Int32(1050)
const KTR_PARAM_MSNUMTOSAVE = Int32(1051)
const KTR_PARAM_MSSAVETOL = Int32(1052)
const KTR_PARAM_PRESOLVEDEBUG = Int32(1053)
const KTR_PARAM_MSTERMINATE = Int32(1054)
const KTR_PARAM_MSSTARTPTRANGE = Int32(1055)
const KTR_PARAM_INFEASTOL = Int32(1056)
const KTR_PARAM_LINSOLVER = Int32(1057)
const KTR_PARAM_BAR_DIRECTINTERVAL = Int32(1058)
const KTR_PARAM_PRESOLVE = Int32(1059)
const KTR_PARAM_PRESOLVE_TOL = Int32(1060)
const KTR_PARAM_BAR_SWITCHRULE = Int32(1061)
const KTR_PARAM_HESSIAN_NO_F = Int32(1062)
const KTR_PARAM_MA_TERMINATE = Int32(1063)
const KTR_PARAM_MA_MAXTIMECPU = Int32(1064)
const KTR_PARAM_MA_MAXTIMEREAL = Int32(1065)
const KTR_PARAM_MSSEED = Int32(1066)
const KTR_PARAM_MA_OUTSUB = Int32(1067)
const KTR_PARAM_MS_OUTSUB = Int32(1068)
const KTR_PARAM_XPRESSLIB = Int32(1069)
const KTR_PARAM_TUNER = Int32(1070)
const KTR_PARAM_TUNER_OPTIONSFILE = Int32(1071)
const KTR_PARAM_TUNER_MAXTIMECPU = Int32(1072)
const KTR_PARAM_TUNER_MAXTIMEREAL = Int32(1073)
const KTR_PARAM_TUNER_OUTSUB = Int32(1074)
const KTR_PARAM_TUNER_TERMINATE = Int32(1075)
const KTR_PARAM_LINSOLVER_OOC = Int32(1076)
const KTR_PARAM_BAR_RELAXCONS = Int32(1077)
const KTR_PARAM_MSDETERMINISTIC = Int32(1078)
const KTR_PARAM_BAR_REFINEMENT = Int32(1079)
const KTR_PARAM_DERIVCHECK = Int32(1080)
const KTR_PARAM_DERIVCHECK_TYPE = Int32(1081)
const KTR_PARAM_DERIVCHECK_TOL = Int32(1082)
const KTR_PARAM_LINSOLVER_INEXACT = Int32(1083)
const KTR_PARAM_LINSOLVER_INEXACTTOL = Int32(1084)
const KTR_PARAM_MAXFEVALS = Int32(1085)
const KTR_PARAM_FSTOPVAL = Int32(1086)
const KTR_PARAM_DATACHECK = Int32(1087)
const KTR_PARAM_DERIVCHECK_TERMINATE = Int32(1088)
const KTR_PARAM_BAR_WATCHDOG = Int32(1089)
const KTR_PARAM_FTOL = Int32(1090)
const KTR_PARAM_FTOL_ITERS = Int32(1091)
const KTR_PARAM_ACT_QPALG = Int32(1092)
const KTR_PARAM_BAR_INITPI_MPEC = Int32(1093)
const KTR_PARAM_XTOL_ITERS = Int32(1094)
const KTR_PARAM_LINESEARCH = Int32(1095)
const KTR_PARAM_OUT_CSVINFO = Int32(1096)
const KTR_PARAM_INITPENALTY = Int32(1097)
const KTR_PARAM_ACT_LPFEASTOL = Int32(1098)
const KTR_PARAM_CG_STOPTOL = Int32(1099)
const KTR_PARAM_RESTARTS = Int32(1100)
const KTR_PARAM_RESTARTS_MAXIT = Int32(1101)
const KTR_PARAM_BAR_SLACKBOUNDPUSH = Int32(1102)
const KTR_PARAM_CG_PMEM = Int32(1103)
const KTR_PARAM_BAR_SWITCHOBJ = Int32(1104)
const KTR_PARAM_OUTNAME = Int32(1105)
const KTR_PARAM_OUT_CSVNAME = Int32(1106)
const KTR_PARAM_ACT_PARAMETRIC = Int32(1107)
const KTR_PARAM_ACT_LPDUMPMPS = Int32(1108)
const KTR_PARAM_ACT_LPALG = Int32(1109)
const KTR_PARAM_ACT_LPPRESOLVE = Int32(1110)
const KTR_PARAM_ACT_LPPENALTY = Int32(1111)
const KTR_PARAM_CONVEX_TUNE = Int32(1114)
const KTR_PARAM_OUT_HINTS = Int32(1115)
const KTR_PARAM_MIP_METHOD = Int32(2001)
const KTR_PARAM_MIP_BRANCHRULE = Int32(2002)
const KTR_PARAM_MIP_SELECTRULE = Int32(2003)
const KTR_PARAM_MIP_INTGAPABS = Int32(2004)
const KTR_PARAM_MIP_INTGAPREL = Int32(2005)
const KTR_PARAM_MIP_MAXTIMECPU = Int32(2006)
const KTR_PARAM_MIP_MAXTIMEREAL = Int32(2007)
const KTR_PARAM_MIP_MAXSOLVES = Int32(2008)
const KTR_PARAM_MIP_INTEGERTOL = Int32(2009)
const KTR_PARAM_MIP_OUTLEVEL = Int32(2010)
const KTR_PARAM_MIP_OUTINTERVAL = Int32(2011)
const KTR_PARAM_MIP_OUTSUB = Int32(2012)
const KTR_PARAM_MIP_DEBUG = Int32(2013)
const KTR_PARAM_MIP_IMPLICATNS = Int32(2014)
const KTR_PARAM_MIP_GUB_BRANCH = Int32(2015)
const KTR_PARAM_MIP_KNAPSACK = Int32(2016)
const KTR_PARAM_MIP_ROUNDING = Int32(2017)
const KTR_PARAM_MIP_ROOTALG = Int32(2018)
const KTR_PARAM_MIP_LPALG = Int32(2019)
const KTR_PARAM_MIP_TERMINATE = Int32(2020)
const KTR_PARAM_MIP_MAXNODES = Int32(2021)
const KTR_PARAM_MIP_HEURISTIC = Int32(2022)
const KTR_PARAM_MIP_HEUR_MAXIT = Int32(2023)
const KTR_PARAM_MIP_HEUR_MAXTIMECPU = Int32(2024)
const KTR_PARAM_MIP_HEUR_MAXTIMEREAL = Int32(2025)
const KTR_PARAM_MIP_PSEUDOINIT = Int32(2026)
const KTR_PARAM_MIP_STRONG_MAXIT = Int32(2027)
const KTR_PARAM_MIP_STRONG_CANDLIM = Int32(2028)
const KTR_PARAM_MIP_STRONG_LEVEL = Int32(2029)
const KTR_PARAM_MIP_INTVAR_STRATEGY = Int32(2030)
const KTR_PARAM_MIP_RELAXABLE = Int32(2031)
const KTR_PARAM_MIP_NODEALG = Int32(2032)
const KTR_PARAM_MIP_HEUR_TERMINATE = Int32(2033)
const KTR_PARAM_MIP_SELECTDIR = Int32(2034)
const KTR_PARAM_PAR_NUMTHREADS = Int32(3001)
const KTR_PARAM_PAR_CONCURRENT_EVALS = Int32(3002)
const KTR_PARAM_PAR_BLASNUMTHREADS = Int32(3003)
const KTR_PARAM_PAR_LSNUMTHREADS = Int32(3004)
const KTR_PARAM_PAR_MSNUMTHREADS = Int32(3005)

const KTR_NEWPOINT_NONE = Int32(0)
const KTR_NEWPOINT_SAVEONE = Int32(1)
const KTR_NEWPOINT_SAVEALL = Int32(2)
const KTR_NEWPOINT_USER = Int32(3)
const KTR_HONORBNDS_NO = Int32(0)
const KTR_HONORBNDS_ALWAYS = Int32(1)
const KTR_HONORBNDS_INITPT = Int32(2)
const KTR_ALG_AUTOMATIC = Int32(0)
const KTR_ALG_AUTO = Int32(0)
const KTR_ALG_BAR_DIRECT = Int32(1)
const KTR_ALG_BAR_CG = Int32(2)
const KTR_ALG_ACT_CG = Int32(3)
const KTR_ALG_IPDIRECT = Int32(1)
const KTR_ALG_IPCG = Int32(2)
const KTR_ALG_ACTIVE = Int32(3)
const KTR_ALG_ACT_SQP = Int32(4)
const KTR_ALG_MULTI = Int32(5)
const KTR_BAR_MURULE_AUTOMATIC = Int32(0)
const KTR_BAR_MURULE_AUTO = Int32(0)
const KTR_BAR_MURULE_MONOTONE = Int32(1)
const KTR_BAR_MURULE_ADAPTIVE = Int32(2)
const KTR_BAR_MURULE_PROBING = Int32(3)
const KTR_BAR_MURULE_DAMPMPC = Int32(4)
const KTR_BAR_MURULE_FULLMPC = Int32(5)
const KTR_BAR_MURULE_QUALITY = Int32(6)
const KTR_BAR_FEASIBLE_NO = Int32(0)
const KTR_BAR_FEASIBLE_STAY = Int32(1)
const KTR_BAR_FEASIBLE_GET = Int32(2)
const KTR_BAR_FEASIBLE_GET_STAY = Int32(3)
const KTR_FEASIBLE_NO = Int32(0)
const KTR_FEASIBLE_ALWAYS = Int32(1)
const KTR_GRADOPT_EXACT = Int32(1)
const KTR_GRADOPT_FORWARD = Int32(2)
const KTR_GRADOPT_CENTRAL = Int32(3)
const KTR_GRADOPT_USER_FORWARD = Int32(4)
const KTR_GRADOPT_USER_CENTRAL = Int32(5)
const KTR_HESSOPT_AUTO = Int32(0)
const KTR_HESSOPT_EXACT = Int32(1)
const KTR_HESSOPT_BFGS = Int32(2)
const KTR_HESSOPT_SR1 = Int32(3)
const KTR_HESSOPT_FINITE_DIFF = Int32(4)
const KTR_HESSOPT_PRODUCT_FINDIFF = Int32(4)
const KTR_HESSOPT_PRODUCT = Int32(5)
const KTR_HESSOPT_LBFGS = Int32(6)
const KTR_HESSOPT_GAUSS_NEWTON = Int32(7)
const KTR_BAR_INITPT_AUTO = Int32(0)
const KTR_BAR_INITPT_YES = Int32(1)
const KTR_BAR_INITPT_STRAT1 = Int32(1)
const KTR_BAR_INITPT_NO = Int32(2)
const KTR_BAR_INITPT_STRAT2 = Int32(2)
const KTR_BAR_INITPT_STRAT3 = Int32(3)
const KTR_ACT_LPSOLVER_INTERNAL = Int32(1)
const KTR_ACT_LPSOLVER_CPLEX = Int32(2)
const KTR_ACT_LPSOLVER_XPRESS = Int32(3)
const KTR_LP_INTERNAL = Int32(1)
const KTR_LP_CPLEX = Int32(2)
const KTR_LP_XPRESS = Int32(3)
const KTR_OUTLEV_NONE = Int32(0)
const KTR_OUTLEV_SUMMARY = Int32(1)
const KTR_OUTLEV_MAJORIT10 = Int32(2)
const KTR_OUTLEV_ITER_10 = Int32(2)
const KTR_OUTLEV_MAJORIT = Int32(3)
const KTR_OUTLEV_ITER = Int32(3)
const KTR_OUTLEV_ALLIT = Int32(4)
const KTR_OUTLEV_ITER_VERBOSE = Int32(4)
const KTR_OUTLEV_ALLIT_X = Int32(5)
const KTR_OUTLEV_ITER_X = Int32(5)
const KTR_OUTLEV_ALL = Int32(6)
const KTR_OUTMODE_SCREEN = Int32(0)
const KTR_OUTMODE_FILE = Int32(1)
const KTR_OUTMODE_BOTH = Int32(2)
const KTR_SCALE_NEVER = Int32(0)
const KTR_SCALE_NO = Int32(0)
const KTR_SCALE_ALLOW = Int32(1)
const KTR_SCALE_USER_INTERNAL = Int32(1)
const KTR_SCALE_USER_NONE = Int32(2)
const KTR_SCALE_INTERNAL = Int32(3)
const KTR_SOC_NO = Int32(0)
const KTR_SOC_MAYBE = Int32(1)
const KTR_SOC_YES = Int32(2)
const KTR_DEBUG_NONE = Int32(0)
const KTR_DEBUG_PROBLEM = Int32(1)
const KTR_DEBUG_EXECUTION = Int32(2)
const KTR_MULTISTART_NO = Int32(0)
const KTR_MULTISTART_YES = Int32(1)
const KTR_CG_PRECOND_NONE = Int32(0)
const KTR_CG_PRECOND_CHOL = Int32(1)
const KTR_BLASOPTION_KNITRO = Int32(0)
const KTR_BLASOPTION_INTEL = Int32(1)
const KTR_BLASOPTION_DYNAMIC = Int32(2)
const KTR_OUTAPPEND_NO = Int32(0)
const KTR_OUTAPPEND_YES = Int32(1)
const KTR_BAR_PENRULE_AUTO = Int32(0)
const KTR_BAR_PENRULE_SINGLE = Int32(1)
const KTR_BAR_PENRULE_FLEX = Int32(2)
const KTR_BAR_PENCONS_AUTO = Int32(0)
const KTR_BAR_PENCONS_NONE = Int32(1)
const KTR_BAR_PENCONS_ALL = Int32(2)
const KTR_BAR_PENCONS_EQUALITIES = Int32(3)
const KTR_BAR_PENCONS_INFEAS = Int32(4)
const KTR_PRESOLVEDBG_NONE = Int32(0)
const KTR_PRESOLVEDBG_BASIC = Int32(1)
const KTR_PRESOLVEDBG_VERBOSE = Int32(2)
const KTR_MSTERMINATE_MAXSOLVES = Int32(0)
const KTR_MSTERMINATE_OPTIMAL = Int32(1)
const KTR_MSTERMINATE_FEASIBLE = Int32(2)
const KTR_MSTERMINATE_ANY = Int32(3)
const KTR_LINSOLVER_AUTO = Int32(0)
const KTR_LINSOLVER_INTERNAL = Int32(1)
const KTR_LINSOLVER_HYBRID = Int32(2)
const KTR_LINSOLVER_DENSEQR = Int32(3)
const KTR_LINSOLVER_MA27 = Int32(4)
const KTR_LINSOLVER_MA57 = Int32(5)
const KTR_LINSOLVER_MKLPARDISO = Int32(6)
const KTR_LINSOLVER_MA97 = Int32(7)
const KTR_LINSOLVER_MA86 = Int32(8)
const KTR_PRESOLVE_NONE = Int32(0)
const KTR_PRESOLVE_BASIC = Int32(1)
const KTR_PRESOLVE_ADVANCED = Int32(2)
const KTR_BAR_SWITCHRULE_AUTO = Int32(0)
const KTR_BAR_SWITCHRULE_NEVER = Int32(1)
const KTR_BAR_SWITCHRULE_LEVEL1 = Int32(2)
const KTR_BAR_SWITCHRULE_LEVEL2 = Int32(3)
const KTR_HESSIAN_NO_F_FORBID = Int32(0)
const KTR_HESSIAN_NO_F_ALLOW = Int32(1)
const KTR_MA_TERMINATE_ALL = Int32(0)
const KTR_MA_TERMINATE_OPTIMAL = Int32(1)
const KTR_MA_TERMINATE_FEASIBLE = Int32(2)
const KTR_MA_TERMINATE_ANY = Int32(3)
const KTR_MA_OUTSUB_NONE = Int32(0)
const KTR_MA_OUTSUB_YES = Int32(1)
const KTR_MS_OUTSUB_NONE = Int32(0)
const KTR_MS_OUTSUB_YES = Int32(1)
const KTR_TUNER_OFF = Int32(0)
const KTR_TUNER_ON = Int32(1)
const KTR_TUNER_OUTSUB_NONE = Int32(0)
const KTR_TUNER_OUTSUB_YES = Int32(1)
const KTR_TUNER_OUTSUB_SUMMARY = Int32(1)
const KTR_TUNER_OUTSUB_ALL = Int32(2)
const KTR_TUNER_TERMINATE_ALL = Int32(0)
const KTR_TUNER_TERMINATE_OPTIMAL = Int32(1)
const KTR_TUNER_TERMINATE_FEASIBLE = Int32(2)
const KTR_TUNER_TERMINATE_ANY = Int32(3)
const KTR_LINSOLVER_OOC_NO = Int32(0)
const KTR_LINSOLVER_OOC_MAYBE = Int32(1)
const KTR_LINSOLVER_OOC_YES = Int32(2)
const KTR_BAR_RELAXCONS_NONE = Int32(0)
const KTR_BAR_RELAXCONS_EQS = Int32(1)
const KTR_BAR_RELAXCONS_INEQS = Int32(2)
const KTR_BAR_RELAXCONS_ALL = Int32(3)
const KTR_MSDETERMINISTIC_NO = Int32(0)
const KTR_MSDETERMINISTIC_YES = Int32(1)
const KTR_BAR_REFINEMENT_NO = Int32(0)
const KTR_BAR_REFINEMENT_YES = Int32(1)
const KTR_DERIVCHECK_NONE = Int32(0)
const KTR_DERIVCHECK_FIRST = Int32(1)
const KTR_DERIVCHECK_SECOND = Int32(2)
const KTR_DERIVCHECK_ALL = Int32(3)
const KTR_DERIVCHECK_FORWARD = Int32(1)
const KTR_DERIVCHECK_CENTRAL = Int32(2)
const KTR_LINSOLVER_INEXACT_NO = Int32(0)
const KTR_LINSOLVER_INEXACT_YES = Int32(1)
const KTR_DATACHECK_NO = Int32(0)
const KTR_DATACHECK_YES = Int32(1)
const KTR_DERIVCHECK_STOPERROR = Int32(1)
const KTR_DERIVCHECK_STOPALWAYS = Int32(2)
const KTR_BAR_WATCHDOG_NO = Int32(0)
const KTR_BAR_WATCHDOG_YES = Int32(1)
const KTR_ACT_QPALG_AUTO = Int32(0)
const KTR_ACT_QPALG_BAR_DIRECT = Int32(1)
const KTR_ACT_QPALG_BAR_CG = Int32(2)
const KTR_ACT_QPALG_ACT_CG = Int32(3)
const KTR_LINESEARCH_AUTO = Int32(0)
const KTR_LINESEARCH_BACKTRACK = Int32(1)
const KTR_LINESEARCH_INTERPOLATE = Int32(2)
const KTR_OUT_CSVINFO_NO = Int32(0)
const KTR_OUT_CSVINFO_YES = Int32(1)
const KTR_BAR_SWITCHOBJ_NONE = Int32(0)
const KTR_BAR_SWITCHOBJ_SCALARPROX = Int32(1)
const KTR_BAR_SWITCHOBJ_DIAGPROX = Int32(2)
const KTR_ACT_PARAMETRIC_NO = Int32(0)
const KTR_ACT_PARAMETRIC_MAYBE = Int32(1)
const KTR_ACT_PARAMETRIC_YES = Int32(2)
const KTR_ACT_LPDUMPMPS_NO = Int32(0)
const KTR_ACT_LPDUMPMPS_YES = Int32(1)
const KTR_ACT_LPALG_DEFAULT = Int32(0)
const KTR_ACT_LPALG_PRIMAL = Int32(1)
const KTR_ACT_LPALG_DUAL = Int32(2)
const KTR_ACT_LPALG_BARRIER = Int32(3)
const KTR_ACT_LPPRESOLVE_OFF = Int32(0)
const KTR_ACT_LPPRESOLVE_ON = Int32(1)
const KTR_ACT_LPPENALTY_AUTO = Int32(0)
const KTR_ACT_LPPENALTY_ALL = Int32(1)
const KTR_ACT_LPPENALTY_NONLINEAR = Int32(2)
const KTR_ACT_LPPENALTY_DYNAMIC = Int32(3)
const KTR_CONVEX_TUNE_NO = Int32(0)
const KTR_CONVEX_TUNE_YES = Int32(1)
const KTR_OUT_HINTS_NO = Int32(0)
const KTR_OUT_HINTS_YES = Int32(1)
const KTR_MIP_METHOD_AUTO = Int32(0)
const KTR_MIP_METHOD_BB = Int32(1)
const KTR_MIP_METHOD_HQG = Int32(2)
const KTR_MIP_METHOD_MISQP = Int32(3)
const KTR_MIP_BRANCH_AUTO = Int32(0)
const KTR_MIP_BRANCH_MOSTFRAC = Int32(1)
const KTR_MIP_BRANCH_PSEUDOCOST = Int32(2)
const KTR_MIP_BRANCH_STRONG = Int32(3)
const KTR_MIP_SEL_AUTO = Int32(0)
const KTR_MIP_SEL_DEPTHFIRST = Int32(1)
const KTR_MIP_SEL_BESTBOUND = Int32(2)
const KTR_MIP_SEL_COMBO_1 = Int32(3)
const KTR_MIP_OUTLEVEL_NONE = Int32(0)
const KTR_MIP_OUTLEVEL_ITERS = Int32(1)
const KTR_MIP_OUTLEVEL_ITERSTIME = Int32(2)
const KTR_MIP_OUTLEVEL_ROOT = Int32(3)
const KTR_MIP_OUTSUB_NONE = Int32(0)
const KTR_MIP_OUTSUB_YES = Int32(1)
const KTR_MIP_OUTSUB_YESPROB = Int32(2)
const KTR_MIP_DEBUG_NONE = Int32(0)
const KTR_MIP_DEBUG_ALL = Int32(1)
const KTR_MIP_IMPLICATNS_NO = Int32(0)
const KTR_MIP_IMPLICATNS_YES = Int32(1)
const KTR_MIP_GUB_BRANCH_NO = Int32(0)
const KTR_MIP_GUB_BRANCH_YES = Int32(1)
const KTR_MIP_KNAPSACK_NO = Int32(0)
const KTR_MIP_KNAPSACK_INEQ = Int32(1)
const KTR_MIP_KNAPSACK_INEQ_EQ = Int32(2)
const KTR_MIP_ROUND_AUTO = Int32(0)
const KTR_MIP_ROUND_NONE = Int32(1)
const KTR_MIP_ROUND_HEURISTIC = Int32(2)
const KTR_MIP_ROUND_NLP_SOME = Int32(3)
const KTR_MIP_ROUND_NLP_ALWAYS = Int32(4)
const KTR_MIP_ROOTALG_AUTO = Int32(0)
const KTR_MIP_ROOTALG_BAR_DIRECT = Int32(1)
const KTR_MIP_ROOTALG_BAR_CG = Int32(2)
const KTR_MIP_ROOTALG_ACT_CG = Int32(3)
const KTR_MIP_ROOTALG_ACT_SQP = Int32(4)
const KTR_MIP_ROOTALG_MULTI = Int32(5)
const KTR_MIP_LPALG_AUTO = Int32(0)
const KTR_MIP_LPALG_BAR_DIRECT = Int32(1)
const KTR_MIP_LPALG_BAR_CG = Int32(2)
const KTR_MIP_LPALG_ACT_CG = Int32(3)
const KTR_MIP_TERMINATE_OPTIMAL = Int32(0)
const KTR_MIP_TERMINATE_FEASIBLE = Int32(1)
const KTR_MIP_HEURISTIC_AUTO = Int32(0)
const KTR_MIP_HEURISTIC_NONE = Int32(1)
const KTR_MIP_HEURISTIC_FEASPUMP = Int32(2)
const KTR_MIP_HEURISTIC_MPEC = Int32(3)
const KTR_MIP_PSEUDOINIT_AUTO = Int32(0)
const KTR_MIP_PSEUDOINIT_AVE = Int32(1)
const KTR_MIP_PSEUDOINIT_STRONG = Int32(2)
const KTR_MIP_INTVAR_STRATEGY_NONE = Int32(0)
const KTR_MIP_INTVAR_STRATEGY_RELAX = Int32(1)
const KTR_MIP_INTVAR_STRATEGY_MPEC = Int32(2)
const KTR_MIP_RELAXABLE_NONE = Int32(0)
const KTR_MIP_RELAXABLE_ALL = Int32(1)
const KTR_MIP_NODEALG_AUTO = Int32(0)
const KTR_MIP_NODEALG_BAR_DIRECT = Int32(1)
const KTR_MIP_NODEALG_BAR_CG = Int32(2)
const KTR_MIP_NODEALG_ACT_CG = Int32(3)
const KTR_MIP_NODEALG_ACT_SQP = Int32(4)
const KTR_MIP_NODEALG_MULTI = Int32(5)
const KTR_MIP_HEUR_TERMINATE_FEASIBLE = Int32(1)
const KTR_MIP_HEUR_TERMINATE_LIMIT = Int32(2)
const KTR_MIP_SELECTDIR_DOWN = Int32(0)
const KTR_MIP_SELECTDIR_UP = Int32(1)
const KTR_PAR_CONCURRENT_EVALS_NO = Int32(0)
const KTR_PAR_CONCURRENT_EVALS_YES = Int32(1)
const KTR_PAR_MSNUMTHREADS_AUTO = Int32(0)


text:
# Noise radius illustration using e.g. Poincare section
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

ro = Systems.roessler([0.1,0.2,0.1])
fig, axs = subplots(2, 1; figsize = (0.4*figx, 1.5figy))
plane = (2, 0.0)

err = (1e-4, 1e-12)
εs = 10 .^ (-5:0.5:1)

C = ["C0", "C2"]

for (i, e) ∈ enumerate(err)
    p = poincaresos(
        ro, plane, 10000; Ttr = 100.0,
        rootkw = (xrtol = e, atol = e)
    )

    axs[1].scatter(p[:, 1], p[:, 3]; s = 20, alpha = 0.75/i^2,
    label = "tol = \$10^{$(round(Int, log10(e)))}\$", color = C[i])
    Cs = correlationsum(p, εs)
    y = log.(Cs ./ maximum(Cs))
    axs[2].plot(log.(εs), y; color = C[i])

    # add noise line
    if i == 1
        z = log.(εs)[end-3]
        axs[2].plot([z, z], [-14, y[end-3]]; color = "C1", lw = 2, ls = ":")
        axs[2].text(z+0.2, -9, "\$\\sigma\$"; color = "C1")
    end
end

axs[1].set_xlabel("\$x\$";labelpad = -15)
axs[1].set_ylabel("\$z\$")
axs[1].legend(markerscale = 5, fontsize=26)
axs[1].set_yticks(0:6:18)

axs[2].plot([-9, 0], 0.9 .* [-9, 0]; ls = "--", color = C[2])
axs[2].text(-7, -4, "0.9"; color = C[2])
axs[2].plot([-6, -3], 2.0 .* [-6, -3] .- 2; ls = "--", color = C[1])
axs[2].text(-3, -11, "2"; color = C[1])
axs[2].set_ylabel("\$\\log(C)\$"; labelpad = -10)
axs[2].set_xlabel("\$\\log(\\varepsilon)\$"; labelpad = -12)
axs[2].set_xticks(-12:5:2)
axs[2].set_yticks(-14:5:2)

fig.tight_layout(pad = 0.5)
add_identifiers!(fig)
wsave(plotsdir("5", "fractaldim_noise"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# function to generate a Kitaev colorcode
function generateBondColorcodeKitaev(
        lattice :: L
    ) :: Dict{LB,RGB} where {
        LB, N, U,
        S <: AbstractSite,
        B <: AbstractBond{LB,N},
        L <: AbstractLattice{S,B,U}
    }

    # put the default kitaev labels with default colors into a dictonary
    return Dict{LB,RGB}(
        getDefaultLabelX(LB) => RGB(1.0, 0.0, 0.0),
        getDefaultLabelY(LB) => RGB(0.0, 1.0, 0.0),
        getDefaultLabelZ(LB) => RGB(0.0, 0.0, 1.0)
    )
end

# function to generate a Kitaev colorcode
function generateBondColorcodeRandom(lattice::AbstractLattice)
    # TODO
    # This would be way better:
    # labels = unique(label.(bonds(lattice)))
    # colors = distinguishable_colors(length(labels)) # or a colormap
    # Dict([lb => color for (lb, color) in zip(labels, colors)])

    return Dict([
        lb => rand(RGB) for lb in unique(label.(bonds(lattice)))
    ])
end


text:
module tensorboard
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("allocation_description_pb.jl")
  include("tensor_shape_pb.jl")
  include("types_pb.jl")
  include("resource_handle_pb.jl")
  include("tensor_pb.jl")
  include("attr_value_pb.jl")
  include("api_def_pb.jl")
  include("cluster_pb.jl")
  include("cost_graph_pb.jl")
  include("node_def_pb.jl")
  include("op_def_pb.jl")
  include("function_pb.jl")
  include("versions_pb.jl")
  include("graph_pb.jl")
  include("tensor_description_pb.jl")
  include("step_stats_pb.jl")
  include("debug_pb.jl")
  include("verifier_config_pb.jl")
  include("rewriter_config_pb.jl")
  include("config_pb.jl")
  include("cpp_shape_inference_pb.jl")
  include("summary_pb.jl")
  include("event_pb.jl")
  include("variable_pb.jl")
  include("struct_tb_pb.jl")
  include("trackable_object_graph_pb.jl")
  include("saved_object_graph_pb.jl")
  include("saver_pb.jl")
  include("meta_graph_pb.jl")
  include("tfprof_log_pb.jl")

  include("plugins/custom_scalar/layout_pb.jl")
  include("plugins/text/plugin_data_pb.jl")

  #include("plugins/hparams/hparams.jl")

end


text:
"""
for constructing boundary conditions used in the definition of a Boundary
"""
module BoundaryConditions

export AbstractBC
export AbstractComplexBC
export AbstractRealBC
export noBC
export DirichletBC
export NeumannBC
export FloquetBC
export MatchedBC
export LocalBC
export NonLocalBC
export BCLocality
export BCHermiticity
export HermitianBC
export NonHermitianBC

import ..PRINTED_COLOR_DARK
import ..getside

"""
	AbstractBC{SIDE}
"""
abstract type AbstractBC{SIDE} end
"""
	AbstractRealBC{SIDE} <: AbstractBC{SIDE}
"""
abstract type AbstractRealBC{SIDE} <: AbstractBC{SIDE} end
"""
	AbstractComplexBC{SIDE} <: AbstractBC{SIDE}
"""
abstract type AbstractComplexBC{SIDE} <: AbstractBC{SIDE} end

"""
	getside(::AbstractBC) -> side
"""
getside(::AbstractBC{SIDE}) where SIDE = SIDE

"""
	noBC{SIDE} <: AbstractRealBC{SIDE}
"""
struct noBC{SIDE} <: AbstractRealBC{SIDE} end

"""
	DirichletBC{SIDE} <: AbstractRealBC{SIDE}
"""
struct DirichletBC{SIDE} <: AbstractRealBC{SIDE} end

"""
	NeumannBC{SIDE} <: AbstractRealBC{SIDE}
"""
struct NeumannBC{SIDE} <: AbstractRealBC{SIDE} end

"""
	FloquetBC{SIDE} <: AbstractComplexBC{SIDE}

NOT YET IMPLEMENTED
"""
struct FloquetBC{SIDE} <: AbstractComplexBC{SIDE} end

"""
	MatchedBC{SIDE} <: AbstractComplexBC{SIDE}

fields: `in`, `out`
"""
struct MatchedBC{SIDE} <: AbstractComplexBC{SIDE}
	in::Array{Int,1}
	out::Array{Int,1}
end
"""
	MatchedBC{SIDE}([in, out]) -> mbc

`in::Vector` are input channels specified by integer
`out::Vector` are output channels specified by integer
"""
function MatchedBC{SIDE}(;in=Int[],out=Int[]) where SIDE
	typeof(in)<:Number ? in = [in] : nothing
	typeof(out)<:Number ? out = [out] : nothing
	return MatchedBC{SIDE}(in,out)
end

"""
	LocalBC{SIDE}
"""
struct LocalBC{SIDE} end

"""
	LocalBC{SIDE}
"""
struct NonLocalBC{SIDE} end

BCLocality(::Union{noBC{SIDE},DirichletBC{SIDE},NeumannBC{SIDE}}) where SIDE = LocalBC{SIDE}
"""
	BCLocality(bc) -> locality

get locality trait of a boundary condition
"""
BCLocality(::AbstractBC{SIDE}) where SIDE = NonLocalBC{SIDE}

"""
	HermitianBC{SIDE}
"""
struct HermitianBC{SIDE} end

"""
	NonHermitianBC{SIDE}
"""
struct NonHermitianBC{SIDE} end

BCHermiticity(::Union{noBC{SIDE},DirichletBC{SIDE},NeumannBC{SIDE},FloquetBC{SIDE}}) where SIDE = HermitianBC{SIDE}
"""
	BCHermiticity(bc) -> hermiticity

get hermiticity trait of a boundary condition
"""
BCHermiticity(::AbstractBC{SIDE}) where SIDE = NonHermitianBC{SIDE}

for TBC ∈ (noBC,DirichletBC,NeumannBC,FloquetBC)
	@eval begin
		function Base.show(io::IO,::$TBC{SIDE}) where SIDE
			printstyled(io,$TBC,color=PRINTED_COLOR_DARK)
			!get(io,:compact,false) ? print(io," SIDE=",SIDE) : nothing
		end
	end
end

# Pretty Printing
function Base.show(io::IO,mbc::MatchedBC{SIDE}) where SIDE
	printstyled(io,"MatchedBC",color=PRINTED_COLOR_DARK)
	!get(io,:compact,false) ? print(io," SIDE=",SIDE,";") : nothing
	ioc = IOContext(io,:typeinfo=>Array{Int})
	print(ioc," in:", mbc.in)
	print(ioc,"/out:", mbc.out)
end

Base.conj(bc::AbstractRealBC) = bc

Base.conj(bc::MatchedBC{SIDE}) where SIDE = MatchedBC{SIDE}(in=bc.out,out=bc.in)

end # module


text:
__precompile__(true)
module Transpiler

using StaticArrays, Compat
using Sugar, DataStructures

include("clike/shared.jl")
include("clike/opencl/printing.jl")

function empty_caches!()
    empty_replace_cache!()
    empty_hash_dict!()
    return
end

end # module


text:
export @set, @lens, @set!
using MacroTools

"""
    @set assignment

Return a modified copy of deeply nested objects.

# Example
```jldoctest
julia> using Setfield

julia> struct T;a;b end

julia> t = T(1,2)
T(1, 2)

julia> @set t.a = 5
T(5, 2)

julia> t
T(1, 2)

julia> t = @set t.a = T(2,2)
T(T(2, 2), 2)

julia> @set t.a.b = 3
T(T(2, 3), 2)
```
"""
macro set(ex)
    atset_impl(ex, overwrite=false)
end

"""
    @set! assignment

Shortcut for `obj = @set obj...`.

julia> t = (a=1,)
(a = 1,)

julia> @set! t.a=2
(a = 2,)

julia> t
(a = 2,)
"""
macro set!(ex)
    atset_impl(ex, overwrite=true)
end

is_interpolation(x) = x isa Expr && x.head == :$

function parse_obj_lenses(ex)
    if @capture(ex, front_[indices__])
        obj, frontlens = parse_obj_lenses(front)
        if any(is_interpolation, indices)
            if !all(is_interpolation, indices)
                throw(ArgumentError(string(
                    "Constant and non-constant indexing (i.e., indices",
                    " with and without \$) cannot be mixed.")))
            end
            index = esc(Expr(:tuple, [x.args[1] for x in indices]...))
            lens = :(ConstIndexLens{$index}())
        else
            index = esc(Expr(:tuple, indices...))
            lens = :(IndexLens($index))
        end
    elseif @capture(ex, front_.property_)
        obj, frontlens = parse_obj_lenses(front)
        lens = :(PropertyLens{$(QuoteNode(property))}())
    else
        obj = esc(ex)
        return obj, ()
    end
    obj, tuple(frontlens..., lens)
end

function parse_obj_lens(ex)
    obj, lenses = parse_obj_lenses(ex)
    lens = Expr(:call, :compose, lenses...)
    obj, lens
end

function get_update_op(sym::Symbol)
    s = String(sym)
    if !endswith(s, '=') || isdefined(Base, sym)
        # 'x +=' etc. is actually 'x = x +', and so '+=' isn't defined in Base.
        # '>=' however is a function, and not an assignment operator.
        msg = "Operation $sym doesn't look like an assignment"
        throw(ArgumentError(msg))
    end
    Symbol(s[1:end-1])
end

struct _UpdateOp{OP,V}
    op::OP
    val::V
end
(u::_UpdateOp)(x) = u.op(x, u.val)

function atset_impl(ex::Expr; overwrite::Bool)
    @assert ex.head isa Symbol
    @assert length(ex.args) == 2
    ref, val = ex.args
    obj, lens = parse_obj_lens(ref)
    dst = overwrite ? obj : gensym("_")
    val = esc(val)
    ret = if ex.head == :(=)
        quote
            lens = $lens
            $dst = set($obj, lens, $val)
        end
    else
        op = get_update_op(ex.head)
        f = :(_UpdateOp($op,$val))
        quote
            $dst = modify($f, $obj, $lens)
        end
    end
    ret
end

"""
    @lens

Construct a lens from a field access.

# Example

```jldoctest
julia> using Setfield

julia> struct T;a;b;end

julia> t = T("A1", T(T("A3", "B3"), "B2"))
T("A1", T(T("A3", "B3"), "B2"))

julia> l = @lens _.b.a.b
(@lens _.b.a.b)

julia> get(t, l)
"B3"

julia> set(t, l, 100)
T("A1", T(T("A3", 100), "B2"))

julia> t = ("one", "two")
("one", "two")

julia> set(t, (@lens _[1]), "1")
("1", "two")
```

"""
macro lens(ex)
    obj, lens = parse_obj_lens(ex)
    if obj != esc(:_)
        msg = """Cannot parse lens $ex. Lens expressions must start with @lens _"""
        throw(ArgumentError(msg))
    end
    lens
end

has_atlens_support(::Any) = false
has_atlens_support(::Union{PropertyLens, IndexLens, ConstIndexLens, IdentityLens}) =
    true
has_atlens_support(l::ComposedLens) =
    has_atlens_support(l.outer) && has_atlens_support(l.inner)

print_application(io::IO, l::PropertyLens{field}) where {field} = print(io, ".", field)
print_application(io::IO, l::IndexLens) = print(io, "[", join(repr.(l.indices), ", "), "]")
print_application(io::IO, ::ConstIndexLens{I}) where I =
    print(io, "[", join(string.("\$", I), ", "), "]")
print_application(io::IO, l::IdentityLens) = print(io, "")

function print_application(io::IO, l::ComposedLens)
    print_application(io, l.outer)
    print_application(io, l.inner)
end

function Base.show(io::IO, l::Lens)
    if has_atlens_support(l)
        print_in_atlens(io, l)
    else
        show_generic(io, l)
    end
end

function Base.show(io::IO, l::ComposedLens)
    if has_atlens_support(l)
        print_in_atlens(io, l)
    else
        show(io, l.outer)
        print(io, " ∘ ")
        show(io, l.inner)
    end
end

function print_in_atlens(io, l)
    print(io, "(@lens _")
    print_application(io, l)
    print(io, ')')
end

function show_generic(io::IO, args...)
    types = tuple(typeof(io), Base.Iterators.repeated(Any, length(args))...)
    Types = Tuple{types...}
    invoke(show, Types, io, args...)
end
show_generic(args...) = show_generic(stdout, args...)


text:
using Turing, ReverseDiff, Memoization
using DifferentialEquations
using Plots
using Statistics
using StatsPlots
using DataFrames
using CSV
using NLsolve
Turing.setadbackend(:reversediff)
Turing.setrdcache(true)
Turing.turnprogress(true)
using Distributions
using LinearAlgebra

# here we are simulating the signal of two coupled oscillators
# with amplitudes A1 and A2 of kBT/k and kBT/(k+c) where c is the coupling constant
# D = kBT/gamma
# in terms of the standard parameters of a OU process:
# sigma = sqrt(2*D)
# Theta = D/A

function corr_osc(c,N,delta_t)
    # function that returns two time series that are correlated through a coupling coefficient c
    μ = 0.0 # mean is zero
    σ = sqrt(2) # D=1
    Θ1 = 1.0
    Θ2 = 1.0+abs(c)
    total_time = N*delta_t
    W1 = OrnsteinUhlenbeckProcess(Θ1,μ,σ,0.0,1.0)
    W2 = OrnsteinUhlenbeckProcess(Θ2,μ,σ,0.0,1.0)
    prob1 = NoiseProblem(W1,(0.0,total_time))
    prob2 = NoiseProblem(W2,(0.0,total_time))
    sol1 = solve(prob1;dt=delta_t)
    sol2 = solve(prob2;dt=delta_t)

    # creating the two correlated
    x1 = (sol1.u .+ sol2.u)/2
    if c>0
        x2 = (sol1.u .- sol2.u)/2
    else
        x2 = (sol2.u .- sol1.u)/2
    end
    return x1,x2
end

function calc_fundstats(x)
    return x[1]^2+x[end]^2,sum(x[2:end-1].^2),sum(x[1:end-2].*x[2:end-1])
end

function phi_deriv(F,x,a1ep,a1ss,a1c,a2ep,a2ss,a2c,delta_t,N)
    # x[1] = A1, x[2] = A2, x[3]=D
    A1 = x[1]
    A2 = x[2]
    D = x[3]
    b1 = b(D,A1,delta_t)
    b2 = b(D,A2,delta_t)
    Q1 = q(a1ep,a1ss,a1c,b1)
    Q2 = q(a2ep,a2ss,a2c,b2)
    dQ1 = dqdB(a1ep,a1ss,a1c,b1)
    dQ2 = dqdB(a2ep,a2ss,a2c,b2)
    F[1] = -N*A1^2/2 + A1*Q1/2 + b1*D*delta_t*(A1*b1*(N-1)/(1-b1^2)-dQ1/2)
    F[2] = -N*A2^2/2 + A2*Q2/2 + b2*D*delta_t*(A2*b2*(N-1)/(1-b2^2)-dQ2/2)
    F[3] = (b1*(N-1)/(1-b1^2)-dQ1/A1/2)*b1/A1 + (b2*(N-1)/(1-b2^2)-dQ2/A2/2)*b2/A2
end

@model ou_corr(rn1,rn2,T,delta_t) = begin
    ampl1 ~ Uniform(0.0,5.0)
    ampl2 ~ Uniform(0.0,5.0)
    d ~ Uniform(0.0,5.0)
    b1 = exp(-delta_t*d/ampl1)
    b2 = exp(-delta_t*d/ampl2)

    rn1[1] ~ Normal(0,sqrt(ampl1))
    rn2[1] ~ Normal(0,sqrt(ampl2))

    for i=2:T
        rn1[i] ~ Normal(rn1[i-1]*b1,sqrt(ampl1*(1-b1^2)))
        rn2[i] ~ Normal(rn2[i-1]*b2,sqrt(ampl2*(1-b2^2)))
    end
end

# here we would like to run MCMC on several artificial data sets with varying correlation coefficients
# we should test rho from -0.9 to 0.9

results = DataFrame(rho = Float64[],
                    coupling = Float64[],
                    pearson = Float64[],
                    a1 = Float64[],
                    da1 = Float64[],
                    a2 = Float64[],
                    da2 = Float64[],
                    d = Float64[],
                    dd = Float64[],
                    da1da2 = Float64[],
                    da1dd = Float64[],
                    dasdd = Float64[],
                    a1ep = Float64[],
                    a1ss = Float64[],
                    a1c = Float64[],
                    a2ep = Float64[],
                    a2ss = Float64[],
                    a2c = Float64[],
                    c = Float64[],
                    dc = Float64[])

tries = 500
rho = 0.5
delta_t = 0.3
N = 1000
coupling = 2*abs(rho)/(1-abs(rho))*sign(rho)
println("coupling: ",coupling," rho: ",rho)
# do rho tries times
for i in 1:tries
     x1, x2 = corr_osc(coupling,N,delta_t)
    pearson = Statistics.cor(x1,x2)
    println("i = ",i,"Pearson: ",pearson)

    # lets see whether we can estimate c from the data
    y1 = x1 .+ x2
    y2 = x1 .- x2

    a1ep,a1ss,a1c = calc_fundstats(y1)
    a2ep,a2ss,a2c = calc_fundstats(y2)

    # Ornstein-Uhlenbeck process of two coupled oscillators
    model_fct = ou_corr(y1,y2,length(y1),delta_t)
    chn = sample(model_fct, NUTS(0.65), 5000)
#        Turing.emptyrdcache()

    print(describe(chn))
    println(chn)
    ampl1 = Array(chn[:ampl1])
    ampl2 = Array(chn[:ampl2])
    d = Array(chn[:d])

    a1 = mean(ampl1)
    a2 = mean(ampl2)
    da1 = std(ampl1)
    da2 = std(ampl2)
    da1da2 = cov(ampl1,ampl2)
    da1dd = cov(ampl1,d)
    da2dd = cov(ampl2,d)
    d_mean = mean(d)
    dd = std(d)

    if a1>a2
        c = (ampl1 .- ampl2)./ampl2
    else
        c = (ampl1 .- ampl2)./ampl1
    end
    c_mean = mean(c)
    dc = std(c)

    println("A1,A2: ",a1,",",a2)
    println("c estimate: ",c_mean,"std: ",dc)
    println("cross corr: ",da1da2,da1dd,da2dd)

    push!(results,[rho,coupling,pearson,a1,da1,a2,da2,d_mean,dd,da1da2[1],da1dd[1],da2dd[1],a1ep,a1ss,a1c,a2ep,a2ss,a2c,c_mean,dc])
end

println(results)
CSV.write("correlations1k05.csv",results)


text:
# using Revise
using PiecewiseDeterministicMarkovProcesses,DifferentialEquations, LinearAlgebra, Random
const PDMP = PiecewiseDeterministicMarkovProcesses

function F_eva!(xcdot, xc, xd, parms::Vector{Float64}, t::Float64)
	# vector field used for the continuous variable
	xcdot[1] = -(xc[1] - 1.5)
	nothing
end

function R(x)
	return x^4
end

function R_eva(rate, xc, xd, parms, t::Float64, issum::Bool)
	# rate function
	rate_print = parms[1]
	if issum == false
		if xd[1] == 0
			rate[1] = R(xc[1])
			rate[2] = 0.0
			rate[3] = rate_print
			return 0.0, 4.95 #transition 0->1
		else
			rate[1] = 0.0
			rate[2] = 1.0
			rate[3] = rate_print
			return 0.0, 4.95 #transition 1->0
		end
	else
		if xd[1] == 0
			return R(xc[1]) + rate_print, 5. #transition 0->1
		else
			return 1.0 + rate_print, 5. #transition 1->0
		end
	end
end

function Delta_xc_eva(xc, xd, parms, t::Float64, ind_reaction::Int64)
	# this function return the jump in the continuous component
	if ind_reaction == 2
		xc[1] = 0.0
	end
	return true
end

xc0 = [0.0]
xd0 = [0, 1]

nu_eva = [1 0;-1 0;0 1]
parms = [1.]
tf = 100.

println("--> Case simple chv:")
	Random.seed!(1234)
	problem = PDMP.PDMPProblem(F_eva!,R_eva,Delta_xc_eva,nu_eva, xc0, xd0, parms, (0.0, tf))
	dummy_t =  @time PDMP.solve(problem, CHV(Tsit5()); n_jumps = 200)

println("--> For simulations rejection (Tsit5):")
	Random.seed!(123)
	problem = PDMP.PDMPProblem(F_eva!,R_eva,Delta_xc_eva,nu_eva, xc0, xd0, parms, (0.0, tf))
	result1 =  @time PDMP.solve(problem, Rejection(:lsoda); n_jumps = 200)

println("--> Simulation using save_at to see sampling behaviour")
	nj = 51
	parms = [10.0]
	Random.seed!(1234)
	problem = PDMP.PDMPProblem(F_eva!,R_eva,Delta_xc_eva,nu_eva, xc0, xd0, parms, (0.0, tf))
	result3 =  @time PDMP.solve(problem, CHV(Tsit5()); n_jumps = 200, save_positions = (false,true))


text:
using Octavian, StaticArrays, LinearAlgebra, BenchmarkTools, ProgressMeter

# BLAS.set_num_threads(1)

# For laptops that thermally throttle, you can set the `JULIA_SLEEP_BENCH` environment variable for #seconds to sleep before each `@belapsed`
const SLEEPTIME = parse(Float64, get(ENV, "JULIA_SLEEP_BENCH", "0"))
maybe_sleep() = iszero(SLEEPTIME) || sleep(SLEEPTIME)
# BenchmarkTools.DEFAULT_PARAMETERS.samples = 1_000_000
# BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10

matrix_sizes(x::Integer) = (x,x,x)
matrix_sizes(x::NTuple{3}) = x

const matmulmethodnames = [:SMatrix, :MMatrix, :OctavianStatic, :OctavianDynamic];
function fill_bench_results!(br, lp, (M,K,N), t, i, j)
    name = matmulmethodnames[j]
    br[i,j,1] = t
    gflops = 2e-9M*K*N / t
    br[i,j,2] = gflops
    lp[j+1] = (name, gflops)
    nothing
end

function runbenches(sr, ::Type{T} = Float64) where {T}
    bench_results = Array{Float64}(undef, length(sr), 4, 2)
    p = Progress(length(sr))
    last_perfs = Vector{Tuple{Symbol,Union{Float64,NTuple{3,Int}}}}(undef, size(bench_results,2) + 1)
    for (i,s) ∈ enumerate(sr)
        M, K, N = matrix_sizes(s)
        last_perfs[1] = (:Size, (M,K,N))
        Astatic = @SMatrix rand(T, M, K);
        Bstatic = @SMatrix rand(T, K, N);
        maybe_sleep()
        t = @belapsed $(Ref(Astatic))[] * $(Ref(Bstatic))[]
        fill_bench_results!(bench_results, last_perfs, (M,K,N), t, i, 1)
        Amutable = MArray(Astatic);
        Bmutable = MArray(Bstatic);
        Cmutable = MMatrix{M,N,T}(undef);
        maybe_sleep()
        t = @belapsed mul!($Cmutable, $Amutable, $Bmutable)
        fill_bench_results!(bench_results, last_perfs, (M,K,N), t, i, 2)
        Cmutable0 = copy(Cmutable); Cmutable .= NaN;
        maybe_sleep()
        t = @belapsed matmul!($Cmutable, $Amutable, $Bmutable)
        fill_bench_results!(bench_results, last_perfs, (M,K,N), t, i, 3)
        A = Array(Amutable); B = Array(Bmutable); C = Array(Cmutable);
        maybe_sleep()
        t = @belapsed matmul!($C, $A, $B)
        fill_bench_results!(bench_results, last_perfs, (M,K,N), t, i, 4)
        @assert Array(Cmutable) ≈ Array(Cmutable0) ≈ C
        ProgressMeter.next!(p; showvalues = last_perfs)
    end
    bench_results
end

sizerange = 2:48
br = runbenches(sizerange);
using DataFrames, VegaLite

df = DataFrame(@view(br[:,:,2]));
rename!(df, matmulmethodnames);
df.Size = sizerange

function pick_suffix(desc = "")
    suffix = if Bool(Octavian.VectorizationBase.has_feature(Val(:x86_64_avx512f)))
        "AVX512"
    elseif Bool(Octavian.VectorizationBase.has_feature(Val(:x86_64_avx2)))
        "AVX2"
    elseif Bool(Octavian.VectorizationBase.has_feature(Val(:x86_64_avx)))
        "AVX"
    else
        "REGSIZE$(Octavian.VectorizationBase.register_size())"
    end
    if desc != ""
        suffix *= '_' * desc
    end
    "$(Sys.CPU_NAME)_$suffix"
end

dfs = stack(df, matmulmethodnames, variable_name = :MatMulType, value_name = :GFLOPS);
p = dfs |> @vlplot(:line, x = :Size, y = :GFLOPS, width = 900, height = 600, color = {:MatMulType});
save(joinpath(pkgdir(Octavian), "docs/src/assets/sizedarraybenchmarks_$(pick_suffix()).svg"), p)


text:
abstract type CurveInfoSnippet end
struct EmptyInfo <: CurveInfoSnippet end

"""
    CurveDistance(timepoints::AbstractRange)
    CurveDistance(timepoints::Vector{<:AbstractFloat})
Print @info REPL output when curve distance reaches each of the timepoints::Vector{<:AbstractFloat}

Example
c = CurveDistance(0.1:0.1:2.1)

c can then be used as an argument to Verbose(), which then goes as a keyword argument into evolve.
"""
struct CurveDistance{V <: AbstractFloat} <: CurveInfoSnippet
    timepoints::Vector{V}
end

"""
    HamiltonianResidual(timepoints::AbstractRange)
    HamiltonianResidual(timepoints::Vector{<:AbstractFloat})
Print @info REPL output on the numerical value of dHdu, the u-derivative of the Hamiltonian, at timepoints::Vector{<:AbstractFloat}

Example
c = HamiltonianResidual(0.1:0.1:2.1)

c can then be used as an argument to Verbose(), which then goes as a keyword argument into evolve.
"""
struct HamiltonianResidual{V <: AbstractFloat} <: CurveInfoSnippet
    timepoints::Vector{V}
end
CurveDistance(a::AbstractRange) = CurveDistance(a |> collect)
HamiltonianResidual(a::AbstractRange) = HamiltonianResidual(a |> collect)

"""
Construct as `v = Verbose(a::Vector{T}) where T <: CallbackCallable`

Example:
```v = Verbose([HamiltonianResidual(0.1:2.:10), CurveDistance(0.1:0.1:5)])
evolve(c::MDCProblem, ...; mdc_callback=[v, other_mdc_callbacks])
````
In this case, you would get REPL output as the curve evolved, indicating when the curve reached each distance in `0.1:0.1:5`, and indicating the residual on dHdu, which is ideally zero, at `0.1:0.1:5`.
"""
struct Verbose{T <: CurveInfoSnippet} <: CallbackCallable
    snippets::Vector{T}
end

Verbose() = Verbose([EmptyInfo()])
Verbose(snippet::EmptyInfo) = Verbose()
Verbose(snippet::T) where T <: CurveInfoSnippet = Verbose([snippet])

"""
    ParameterBounds(ids, lbs, ubs)
    - ids are the indices of the model parameters that you want to bound
    - lbs are an array of lower bounds, with length == to indices
    - ubs are...well you can guess.
"""
struct ParameterBounds{I <: Integer,T <: Number} <: AdjustmentCallback
    ids::Vector{I}
    lbs::Vector{T}
    ubs::Vector{T}
end

function (c::CurveDistance)(cp::CurveProblem, u, t, integ)
    @info "curve length is $t"
    nothing
end

function (h::HamiltonianResidual)(c::CurveProblem, u, t, integ)
    x = dHdu_residual(c, u, t, nothing)
    @info "dHdu residual = $x at curve length $t"
end

(e::EmptyInfo)(c, u, t, integ) = nothing

function (v::Verbose)(c::CurveProblem)
    to_call = map(v.snippets) do snippet
        (u, t, _integ) -> snippet(c, u, t, _integ)
    end
    return map(to_call, v.snippets) do each, snippet
        FunctionCallingCallback(each; funcat=snippet.timepoints)
    end
end

"""
    DEPRECATED, use Verbose() instead.
    VerboseOutput(level=:low, times = 0:0.1:1.)
    Callback to give online info on how the solution is going, as the MDCurve evolves. activates at curve distances specified by times
"""
function VerboseOutput(level=:low, times=0:0.1:1.)

        function affect!(integ)
            if level == :low
                @info "curve length is $(integ.t)"
            end
            if level == :medium
                @info "dHdu residual = "
            end
            if level == :high

            end
            return integ
        end
        return PresetTimeCallback(times, affect!)
end

"""
    ParameterBounds(ids::Vector{Integer},lbs::Vector{Number},ubs::Vector{Number})
parameters[ids] must fall within lbs and ubs, where lbs and ubs are Arrays of the same size as ids.
Create hard bounds on the parameter space over which the minimally disruptive curve can trace. Curve evolution terminates if it hits a bound.
"""
    function (p::ParameterBounds)(c::CurveProblem)
        function condition(u, t, integrator)
            tests = u[p.ids]
            any(tests .< p.lbs) && return true
            any(tests .> p.ubs) && return true
            return false
        end
        return DiscreteCallback(condition, terminate!)
    end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export skeleton

"""
    skeleton(;kwargs...)
    skeleton(children::Any;kwargs...)
    skeleton(children_maker::Function;kwargs...)

A Skeleton component.
Indicate content loading state. For more information, see: https://mantine.dev/core/skeleton/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Primary content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `circle` (Bool; optional): If Skeleton is a circle, it's width and border-radius will be equal to height
- `height` (String | Real; optional): Skeleton height
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Radius from theme.radius or number to set border-radius in px
- `style` (Dict; optional): Inline style override
- `visible` (Bool; optional): Should skeleton overlay be displayed
- `width` (String | Real; optional): Skeleton width
"""
function skeleton(; kwargs...)
        available_props = Symbol[:children, :id, :circle, :height, :loading_state, :radius, :style, :visible, :width]
        wild_props = Symbol[]
        return Component("skeleton", "Skeleton", "dash_mantine_components", available_props, wild_props; kwargs...)
end

skeleton(children::Any; kwargs...) = skeleton(;kwargs..., children = children)
skeleton(children_maker::Function; kwargs...) = skeleton(children_maker(); kwargs...)


text:
@testset "norm functions" begin
    @testset "$fnorm(x::Array{$T,$(length(sz))})" for
        fnorm in (
            LinearAlgebra.norm1,
            LinearAlgebra.norm2,
            LinearAlgebra.normInf,
            LinearAlgebra.normMinusInf,
        ),
        T in (Float64, ComplexF64),
        sz in [(3,), (3, 3), (3, 2, 1)]

        x = randn(T, sz)
        # finite differences is unstable if maxabs (minabs) values are not well
        # separated from other values
        if fnorm === LinearAlgebra.normInf
            x[end] = 1000rand(T)
            kwargs = (atol=1e-5, rtol=1e-5)
        elseif fnorm == LinearAlgebra.normMinusInf
            x .*= 1000
            x[end] = rand(T)
            kwargs = (atol=1e-5, rtol=1e-5)
        else
            kwargs = NamedTuple()
        end

        y = fnorm(x)
        x̄ = rand_tangent(x)
        ȳ = rand_tangent(y)

        fnorm === LinearAlgebra.norm2 && @testset "frule" begin
            ẋ = rand_tangent(x)
            frule_test(fnorm, (x, ẋ))
        end
        @testset "rrule" begin
            rrule_test(fnorm, ȳ, (x, x̄); kwargs...)
            x isa Matrix && @testset "$MT" for MT in (Diagonal, UpperTriangular, LowerTriangular)
                rrule_test(fnorm, ȳ, (MT(x), MT(x̄)); kwargs...)
            end
            @test extern(rrule(fnorm, zero(x))[2](ȳ)[2]) ≈ zero(x)
            @test rrule(fnorm, x)[2](Zero())[2] isa Zero
        end
    end
    @testset "norm(x::Array{$T,$(length(sz))})" for
        T in (Float64, ComplexF64),
        sz in [(0,), (3,), (3, 3), (3, 2, 1)]

        x = randn(T, sz)
        y = norm(x)
        ẋ = rand_tangent(x)
        x̄ = rand_tangent(x)
        ȳ = rand_tangent(y)

        @testset "frule" begin
            frule_test(norm, (x, ẋ))
            @test frule((Zero(), Zero()), norm, x)[2] isa Zero
            @test iszero(frule((Zero(), ẋ), norm, zero(x))[2])
        end
        @testset "rrule" begin
            rrule_test(norm, ȳ, (x, x̄))
            x isa Matrix && @testset "$MT" for MT in (Diagonal, UpperTriangular, LowerTriangular)
                rrule_test(norm, ȳ, (MT(x), MT(x̄)))
            end
            @test extern(rrule(norm, zero(x))[2](ȳ)[2]) ≈ zero(x)
            @test rrule(norm, x)[2](Zero())[2] isa Zero
        end
    end
    @testset "$fnorm(x::Array{$T,$(length(sz))}, $p) with size $sz" for
        fnorm in (norm, LinearAlgebra.normp),
        p in (1.0, 2.0, Inf, -Inf, 2.5),
        T in (Float64, ComplexF64),
        sz in (fnorm === norm ? [(0,), (3,), (3, 3), (3, 2, 1)] : [(3,), (3, 3), (3, 2, 1)])

        x = randn(T, sz)
        # finite differences is unstable if maxabs (minabs) values are not well
        # separated from other values
        if p == Inf
            if !isempty(x)
                x[end] = 1000rand(T)
            end
            kwargs = (atol=1e-5, rtol=1e-5)
        elseif p == -Inf
            if !isempty(x)
                x .*= 1000
                x[end] = rand(T)
            end
            kwargs = (atol=1e-5, rtol=1e-5)
        else
            kwargs = NamedTuple()
        end

        y = fnorm(x, p)
        x̄ = rand_tangent(x)
        ȳ = rand_tangent(y)
        p̄ = rand_tangent(p)

        rrule_test(fnorm, ȳ, (x, x̄), (p, p̄); kwargs...)
        x isa Matrix && @testset "$MT" for MT in (Diagonal, UpperTriangular, LowerTriangular)
            rrule_test(fnorm, ȳ, (MT(x), MT(x̄)), (p, p̄); kwargs...)
        end
        @test extern(rrule(fnorm, zero(x), p)[2](ȳ)[2]) ≈ zero(x)
        @test rrule(fnorm, x, p)[2](Zero())[2] isa Zero
    end
    @testset "norm($fdual(::Vector{$T}), p)" for
        T in (Float64, ComplexF64),
        fdual in (adjoint, transpose)
        p = 2.5
        n = 3
        x = fdual(randn(T, n))
        y = norm(x, p)
        x̄ = rand_tangent(x)
        ȳ = rand_tangent(y)
        p̄ = rand_tangent(p)
        rrule_test(norm, ȳ, (x, x̄), (p, p̄))
        @test extern(rrule(norm, x, p)[2](ȳ)[2]) isa typeof(x)
    end
    @testset "norm(x::$T, p)" for T in (Float64, ComplexF64)
        @testset "p = $p" for p in (-1.0, 2.0, 2.5)
            x = randn(T)
            y = norm(x, p)
            ẋ, ṗ = rand_tangent.((x, p))
            x̄, p̄, ȳ = rand_tangent.((x, p, y))
            frule_test(norm, (x, ẋ), (p, ṗ))
            rrule_test(norm, ȳ, (x, x̄), (p, p̄))
            _, back = rrule(norm, x, p)
            @test back(Zero()) == (NO_FIELDS, Zero(), Zero())
        end
        @testset "p = 0" begin
            p = 0.0
            x = randn(T)
            y = norm(x, p)
            ẋ, ṗ = rand_tangent.((x, p))
            x̄, p̄, ȳ = rand_tangent.((x, p, y))
            y_fwd, ẏ = frule((Zero(), ẋ, ṗ), norm, x, p)
            @test y_fwd == y
            @test iszero(ẏ)
            y_rev, back = rrule(norm, x, p)
            @test y_rev == y
            @test back(ȳ) == (NO_FIELDS, zero(x), Zero())
            @test back(Zero()) == (NO_FIELDS, Zero(), Zero())
        end
    end
end

@testset "normalize" begin
    @testset "x::Vector{$T}" for T in (Float64, ComplexF64)
        n = 3
        x = randn(T, n)
        y = normalize(x)
        x̄ = rand_tangent(x)
        ȳ = rand_tangent(y)
        rrule_test(normalize, ȳ, (x, x̄))
        @test rrule(normalize, x)[2](Zero()) === (NO_FIELDS, Zero())
    end
    @testset "x::Vector{$T}, p=$p" for T in (Float64, ComplexF64),
        p in (1.0, 2.0, -Inf, Inf, 2.5) # skip p=0, since FD is unstable
        n = 3
        x = randn(T, n)
        y = normalize(x, p)
        x̄ = rand_tangent(x)
        ȳ = rand_tangent(y)
        p̄ = rand_tangent(p)
        rrule_test(normalize, ȳ, (x, x̄), (p, p̄))
        @test rrule(normalize, x, p)[2](Zero()) === (NO_FIELDS, Zero(), Zero())
    end
end


text:
@doc raw"""
# Exponential Growth

"""
module ExponentialGrowth

    using GeometricEquations

    export exponential_growth_ode

    function exponential_growth_ode_f(t, x, f)
        f[1] = x[1]
        nothing
    end

    function exponential_growth_ode(x₀=[1.0])
        ODE(exponential_growth_ode_f, x₀)
    end

end


text:
export Var
export parameter, zerograd!, isnothing, isparam, gradient!, topsort

"""
    Var

Variable struct.

`Var` contains the following members:
* data
* args
* grad

# Example
```julia
T = Float32
x = Var(rand(T,10,5)) # x.grad is set to `nothing`
x = zerograd(rand(T,10,5)) # x.grad is initialized as zero.
```
"""
mutable struct Var
    data
    f
    args
    grad
end

Var(data, f=nothing, args=()) = Var(data, f, args, nothing)

function parameter(data)
    v = Var(data)
    v.grad = fill!(similar(v.data), 0)
    v
end

function zerograd!(x::Var)
    fill!(x.grad, 0)
    x
end

data(x::Var) = x.data
gradient(x::Var) = x.grad
Base.size(x::Var) = size(x.data)
Base.size(x::Var, i::Int) = size(x.data, i)
Base.length(x::Var) = length(x.data)
Base.ndims(x::Var) = ndims(x.data)
Base.eltype(x::Var) = eltype(x.data)
isnothing(x) = x == nothing

#getdevice(x::Var) = getdevice(x.data)
#getdevice(x::Array) = -1
#getdevice(x::CuArray) = CUDA.getdevice(x)

"""
    isparam(x::Var)

Returns whether `x` is a parameter or not
"""
isparam(x) = isa(x,Var) && !isnothing(x.grad) && isempty(x.args)

"""
    topsort(tops::T...)

Topological sort.
"""
function topsort(top::T) where T
    sorted = T[]
    dict = IdDict{T,T}()
    function visit(x::T)
        haskey(dict,x) && return
        dict[x] = x
        for arg in x.args
            isa(arg,T) && visit(arg)
        end
        push!(sorted, x)
    end
    visit(top)
    sorted
end

"""
    gradient!(top::Var)

Compute gradients.
"""
function gradient!(top::Var)
    sorted = topsort(top)
    if isnothing(top.grad)
        top.grad = fill!(similar(top.data), 1)
    end
    for v in sorted
        if !isempty(v.args) && isnothing(v.grad)
            v.grad = fill!(similar(v.data), 0)
        end
    end
    for i = length(sorted):-1:1
        y = sorted[i]
        isnothing(y.grad) && continue
        isempty(y.args) && continue
        y.f(y, y.args...)
    end
    sorted
end

function tocpu!(x::Var)
end

function tocuda!(x::Var)
end

function configure2!(xs::Vararg{Var})
    return

    if iscpu()
        f = tocpu
    elseif iscuda()
        f = tocuda
    end
    for x in xs
        x.data = f(x.data)
        isnothing(x.grad) || (x.grad = f(x.grad))
    end
end