text:
using YaoBase

export ChainBlock, chain

"""
    ChainBlock{D} <: CompositeBlock{D}

`ChainBlock` is a basic construct tool to create
user defined blocks horizontically. It is a `Vector`
like composite type.
"""
struct ChainBlock{D} <: CompositeBlock{D}
    n::Int
    blocks::Vector{AbstractBlock{D}}
    function ChainBlock(n::Int, blocks::Vector{<:AbstractBlock{D}}) where {D}
        _check_block_sizes(blocks, n)
        return new{D}(n, blocks)
    end
end

ChainBlock(blocks::Vector{<:AbstractBlock{D}}) where {D} = ChainBlock(_check_block_sizes(blocks), blocks)
ChainBlock(blocks::AbstractBlock{D}...) where {D} =
    ChainBlock(collect(AbstractBlock{D}, blocks))

nqudits(c::ChainBlock) = c.n

"""
    chain(blocks...)

Return a [`ChainBlock`](@ref) which chains a list of blocks with same
[`nqudits`](@ref). If there is lazy evaluated
block in `blocks`, chain can infer the number of qudits and create an
instance itself.
"""
chain(blocks::AbstractBlock{D}...) where {D} = ChainBlock(blocks...)
function chain(blocks::Union{AbstractBlock{D},Function}...) where {D}
    blocks = filter(x->x isa AbstractBlock, blocks)
    return chain(map(x -> parse_block(_check_block_sizes(blocks...), x), blocks)...)
end

function chain(list::Vector{<:AbstractVector{D}}) where D
    return ChainBlock(list)
end

# if not all matrix block, try to put the number of qudits.
chain(n::Int, blocks...) = chain(map(x -> parse_block(n, x), blocks)...)
chain(n::Int, itr) = isempty(itr) ? chain(n) : chain(map(x -> parse_block(n, x), itr)...)
# disambiguity
# NOTE: we use parse_block here to make sure the behaviour are the same
chain(n::Int, it::Pair) = chain(n, parse_block(n, it))
chain(n::Int, f::Function) = chain(n, parse_block(n, f))
function chain(n::Int, block::AbstractBlock{D}) where D
    if n != nqudits(block)
        throw(QubitMismatchError("number of qudits mismatch: expect $n, got $(nqudits(block))"))
    end
    return ChainBlock(n, AbstractBlock{D}[block])
end
chain(blocks::Function...) = @λ(n -> chain(n, blocks...))
chain(it) = chain(it...) # forward iterator to vargs, so we could dispatch based on types
chain(it::Pair) = error("got $it, do you mean put($it)?")
chain(blocks...) = @λ(n -> chain(n, blocks))

"""
    chain(n)

Return an empty [`ChainBlock`](@ref) which can be used like a list of blocks.
"""
chain(n::Int; nlevel=2) = ChainBlock(n::Int, AbstractBlock{nlevel}[])

"""
    chain()

Return an lambda `n->chain(n)`.
"""
chain() = @λ(n -> chain(n))

subblocks(c::ChainBlock) = c.blocks
occupied_locs(c::ChainBlock) =
    Tuple(unique(Iterators.flatten(occupied_locs(b) for b in subblocks(c))))

chsubblocks(pb::ChainBlock{D}, blocks::Vector{<:AbstractBlock{D}}) where {D} =
    length(blocks) == 0 ? ChainBlock(pb.n, AbstractBlock{D}[]) : ChainBlock(pb.n, blocks)
chsubblocks(pb::ChainBlock, it) = chain(it...)

function mat(::Type{T}, c::ChainBlock{D}) where {T,D}
    if isempty(c.blocks)
        return IMatrix{D^nqudits(c),T}()
    else
        return prod(x -> mat(T, x), Iterators.reverse(c.blocks))
    end
end

function _apply!(r::AbstractRegister, c::ChainBlock)
    for each in c.blocks
        _apply!(r, each)
    end
    return r
end

cache_key(c::ChainBlock) = Tuple(cache_key(each) for each in c.blocks)

function Base.:(==)(lhs::ChainBlock{D}, rhs::ChainBlock{D}) where {D}
    nqudits(lhs) == nqudits(rhs) &&
    (length(lhs.blocks) == length(rhs.blocks)) &&
    all(lhs.blocks .== rhs.blocks)
end

Base.copy(c::ChainBlock{D}) where {D} = ChainBlock(c.n, copy(c.blocks))
Base.similar(c::ChainBlock{D}) where {D} = ChainBlock(c.n, empty!(similar(c.blocks)))
Base.getindex(c::ChainBlock, index) = getindex(c.blocks, index)
Base.getindex(c::ChainBlock, index::Union{UnitRange,Vector}) =
    ChainBlock(c.n, getindex(c.blocks, index))
Base.setindex!(c::ChainBlock{D}, val::AbstractBlock{D}, index::Integer) where {D} =
    (_check_block_sizes(c, val); setindex!(c.blocks, val, index); c)
Base.insert!(c::ChainBlock{D}, index::Integer, val::AbstractBlock{D}) where {D} =
    (_check_block_sizes(c, val); insert!(c.blocks, index, val); c)
Base.adjoint(blk::ChainBlock{D}) where {D} =
    ChainBlock(blk.n, AbstractBlock{D}[adjoint(b) for b in reverse(subblocks(blk))])
Base.lastindex(c::ChainBlock) = lastindex(c.blocks)
## Iterate contained blocks
Base.iterate(c::ChainBlock, st = 1) = iterate(c.blocks, st)
Base.length(c::ChainBlock) = length(c.blocks)
Base.eltype(c::ChainBlock) = eltype(c.blocks)
Base.eachindex(c::ChainBlock) = eachindex(c.blocks)
Base.popfirst!(c::ChainBlock) = popfirst!(c.blocks)
Base.pop!(c::ChainBlock) = pop!(c.blocks)
Base.push!(c::ChainBlock{D}, m::AbstractBlock{D}) where {D} = (_check_block_sizes(c, m); push!(c.blocks, m); c)
Base.push!(c::ChainBlock{D}, f::Function) where {D} = (push!(c.blocks, f(c.n)); c)
Base.pushfirst!(c::ChainBlock{D}, m::AbstractBlock{D}) where {D} =
    (_check_block_sizes(c, m); pushfirst!(c.blocks, m); c)
Base.pushfirst!(c::ChainBlock{D}, f::Function) where {D} = (pushfirst!(c.blocks, f(c.n)); c)
Base.append!(c::ChainBlock{D}, list::Vector{<:AbstractBlock{D}}) where {D} =
    (_check_block_sizes(c, list...); append!(c.blocks, list); c)
Base.append!(c1::ChainBlock{D}, c2::ChainBlock{D}) where {D} =
    (_check_block_sizes(c1, c2); append!(c1.blocks, c2.blocks); c1)
Base.prepend!(c1::ChainBlock{D}, list::Vector{<:AbstractBlock{D}}) where {D} =
    (_check_block_sizes(c1, list...); prepend!(c1.blocks, list); c1)
Base.prepend!(c1::ChainBlock{D}, c2::ChainBlock{D}) where {D} =
    (_check_block_sizes(c1, c2); prepend!(c1.blocks, c2.blocks); c1)

YaoBase.isunitary(c::ChainBlock) = all(isunitary, c.blocks) || isunitary(mat(c))
YaoBase.isreflexive(c::ChainBlock) =
    (iscommute(c.blocks...) && all(isreflexive, c.blocks)) || isreflexive(mat(c))
YaoBase.ishermitian(c::ChainBlock) =
    (all(isreflexive, c.blocks) && iscommute(c.blocks...)) || isreflexive(mat(c))


text:
# stacks

type Stack{S}   # S is the type of the internal dequeue instance
    store::S
end

Stack{T}(ty::Type{T}) = Stack(Deque{T}())
Stack{T}(ty::Type{T}, blksize::Integer) = Stack(Deque{T}(blksize))

isempty(s::Stack) = isempty(s.store)
length(s::Stack) = length(s.store)

top(s::Stack) = back(s.store)

function push!(s::Stack, x)
    push!(s.store, x)
    s
end

#returns a collection that can be used in a for loop
function iter{T}(s::Stack{Deque{T}})
    a = T[]
    for i in s.store
        unshift!(a,i)
    end
    return a
end

pop!(s::Stack) = pop!(s.store)

ext: jl


text:
# # How to write _X_ in parallel.

using FLoops

using LiterateTest                                                     #src
using Test                                                             #src

# ## In-place mutation
#
# Mutable containers can be allocated in the `init` expressions
# (`zeros(3)` in the example below):

@test begin
    local ys  # hide
    @floop for x in 1:10
        xs = [x, 2x, 3x]
        @reduce() do (ys = zeros(3); xs)
            ys .+= xs
        end
    end
    ys
end == [55, 110, 165]

# Mutating objects allocated in the `init` expressions is not data
# race because each basecase "owns" such mutable objects.  However, it
# is incorrect to mutate objects created outside `init` expressions.
#
# See also: [What is the difference of `@reduce` and `@init` to the approach
# using `state[threadid()]`?](@ref faq-state-threadid)
#
# !!! note
#
#     Technically, it is correct to mutate objects in the loop body if
#     the objects are protected by a lock.  However, it means that the
#     code block protected by the lock can only be executed by a
#     single task.  For efficient data parallel loops, it is highly
#     recommended to use **non**-thread-safe data collection (i.e., no
#     lock) and construct the `@reduce` block that efficiently merge
#     two mutable objects.

# ### INCORRECT EXAMPLE

# This example has data race because the array `ys0` is shared across
# all base cases and mutated in parallel.

ys0 = zeros(3)
@dedent let
    @floop for x in 1:10
        xs = [x, 2x, 3x]
        @reduce() do (ys = ys0; xs)
            ys .+= xs
        end
    end
end

# ## [Data race-free reuse of mutable objects using private variables](@id private-variables)
#
# To avoid allocation for each iteration, it is useful to pre-allocate mutable
# objects and reuse them. We can use [`@init`](@ref) macro to do this in a
# data race-free ("thread-safe") manner:

@test begin
    local ys  # hide
    @floop for x in 1:10
        @init xs = Vector{typeof(x)}(undef, 3)
        xs .= (x, 2x, 3x)
        @reduce() do (ys = zeros(3); xs)
            ys .+= xs
        end
    end
    ys
end == [55, 110, 165]

# See also: [What is the difference of `@reduce` and `@init` to the approach
# using `state[threadid()]`?](@ref faq-state-threadid)

# ## Efficient and reproducible usage patterns of random number generators
#
# Julia's default random number generator (RNG) is data race-free for invoking
# from multiple threads; i.e., calls like `randn()` have well-defined
# behaviors. However, for the performance and reproducibility, it is useful to
# directly creating the RNGs. A convenient approach to this is to use a
# [private variable](@ref private-variables):

using Random

MersenneTwister()  # the first invocation of `MersenneTwister` is not data race-free

@dedent let
    @floop for _ in 1:10
        @init rng = MersenneTwister()
        @reduce(s += rand(rng))
    end
end

# The above approach may work well for exploratory purposes. However, it has a
# problem that the computation is not reproducible and each invocation of
# `MersenneTwister` requires an I/O (reading `/dev/urandom`). These problems
# can be solved by, for example, using `randjump` function. First, let us
# construct `ntasks` RNGs to be used.

using Future

ntasks = Threads.nthreads()  # the number of base cases
rngs = [MersenneTwister(123456789)]
let rng = rngs[end]
    for _ in 2:ntasks
        rng = Future.randjump(rng, big(10)^20)
        push!(rngs, rng)
    end
end

# This list of RNGs can be used with some input array by manually partitioning
# the input into `ntasks` chunks:

@dedent let
    xs = 1:10  # input
    chunks = Iterators.partition(xs, cld(length(xs), length(rngs)))
    @floop ThreadedEx(basesize = 1) for (rng, chnk) in zip(rngs, chunks)
        y = 0
        for _ in chnk
            y += rand(rng)
        end
        @reduce(s += y)
    end
end

# Note that the above pattern can also be used with `@threads for` loop.
#
# Another approach is to use a counter-based RNG as illustrated in
# [Monte-Carlo π · FoldsCUDA](https://juliafolds.github.io/FoldsCUDA.jl/dev/examples/monte_carlo_pi/).
# This approach works both on CPU and GPU.


text:
struct AssociativeMap <: Mapping
    amap::Dict
    inv_amap::Dict
end

function AssociativeMap(amap::Dict{S,T}) where {S,T}
    f(x) = Set((x,))
    if !(S <: Set)
        if !(T <: Set)
            amap = Dict((f.(keys(amap)) .=> f.(values(amap))))
            AssociativeMap(amap, invert(amap))
        else
            amap = Dict((f.(keys(amap)) .=> values(amap)))
            AssociativeMap(amap, invert(amap))
        end
    else
        if !(T <: Set)
            amap = Dict((keys(amap) .=> f.(values(amap))))
            AssociativeMap(amap, invert(amap))
        else
            amap = Dict((keys(amap) .=> values(amap)))
            AssociativeMap(amap, invert(amap))
        end
    end
end

# AssociativeMap(f::Function, d) = AssociativeMap(Dict([(x, f(x)) for x in d]))

dom(m::AssociativeMap) = Set(collect(keys(m.amap)))
image(m::AssociativeMap) = dom(inverse(m))
inverse(m::AssociativeMap) = AssociativeMap(m.inv_amap, m.amap)
inv(m::AssociativeMap) = inverse(m)
applymap(m::AssociativeMap, x) = m.amap[Set((x,))]

function invert(D::Dict)
    invD = Dict()
    for k in keys(D)
        if D[k] in keys(invD)
            push!(invD[D[k]], k)
        else
            invD[D[k]] = [k]
        end
    end
    return invD
end

# Base.ComposedFunction{AssociativeMap, AssociativeMap} <: Core.Function
# Base.ComposedFunction{AssociativeMap, Function} <: Core.Function
# Base.ComposedFunction{Function, AssociativeMap} <: Core.Function

# function Base.insert!(m::AssociativeMap, newpair)
#     k = newpair[1]
#     v = newpair[2]
#     if k in keys(m.amap)
#         m.amap[k] = m.amap[k] ∪ Set((v,))
#     else
#         m.amap[k] = v
#     end
#     if m.amap[k] in keys(m.inv_amap)
#         m.inv_amap[m.amap[k]] = m.amap[k] ∪ Set((v,))
#     else
#         m.inv_amap[m.amap[k]] = Set((k,))
#     end
#     return m
# end


text:
module MyInterpolations
function my_lin_interp(grid, vals)
    function func(x)
        if x < grid[1]
            print("error")
        end
        if x > grid[end]
            print("error")
        end
        if grid[1] <= x < grid[end]
        i = searchsortedlast(grid, x)
        x in grid[i]:grid[i+1]
        return vals[i] + (vals[i+1] - vals[i])*(x - grid[i])/(grid[i+1] - grid[i])
        end
        if x == grid[end]
        return vals[end]
        end
    end
    return func
end
    export my_lin_interp
end


text:
abstract Zwierzę

type Drapieżnik <: Zwierzę
  Pozycja :: Array{Int64, 1}
  Drapieżnik(x, y) = new([y, x])
end

type Ofiara <: Zwierzę
  Pozycja :: Array{Int64, 1}
  Ofiara(x, y) = new([y, x])
end

type Pustak <: Zwierzę
  Pustak() = new()
end

type Świat
  Mapa :: Array{Zwierzę, 2}
  Świat(n) = new(fill(Pustak(), n, n))
end

function DodajZwierzę{T <: Zwierzę}(MójŚwiat, x, y, Zwierz::Type{T})
  if(typeof(MójŚwiat.Mapa[y, x]) != Pustak)
    error("Zajęte")
  end
  MójŚwiat.Mapa[y, x] = Zwierz(x, y)
end

function Odległość{T<: Zwierzę, W <: Zwierzę}(Zwierzę1::T, Zwierzę2::W)
  return abs(Zwierzę1.Pozycja[1] - Zwierzę2.Pozycja[1]) + abs(Zwierzę1.Pozycja[1] - Zwierzę2.Pozycja[1])
end

function Usuń(MójŚwiat::Świat, Zwierzę1::Zwierzę)
  MójŚwiat.Mapa[Zwierzę1.Pozycja[1], Zwierzę1.Pozycja[2]] = Pustak()
end

function Interakcja(MójŚwiat::Świat, Zwierzę1::Drapieżnik, Zwierzę2::Ofiara)
  Usuń(MójŚwiat, Zwierzę2)
end

function Uciekaj(MójŚwiat::Świat, Zwierzę1::Ofiara)
  for i in 1:size(MójŚwiat.Mapa)[1]
    for j in 1:size(MójŚwiat.Mapa)[2]
      if(typeof(MójŚwiat.Mapa[i, j]) == Pustak)
        MójŚwiat.Mapa[i, j] = Zwierzę1
        Zwierzę1.Pozycja = [i, j]
        return
      end
    end
  end
end

function Interakcja(MójŚwiat::Świat, Zwierzę1::Ofiara, Zwierzę2::Drapieżnik)
  Usuń(MójŚwiat, Zwierzę1)
  Uciekaj(MójŚwiat::Świat, Zwierzę1::Ofiara)
end

function Interakcja(MójŚwiat::Świat, Zwierzę1::Drapieżnik, Zwierzę2::Drapieżnik)
  return "WRrrr"
end

function Interakcja(MójŚwiat::Świat, Zwierzę1::Ofiara, Zwierzę2::Ofiara)
  return "Beeee"
end


text:
using Visualize
resolution = (500, 500)
w = GLWindow.create_glcontext(resolution = resolution)

proj = orthographicprojection(SimpleRectangle(0, 0, resolution...), -10_000f0, 10_000f0)
uniforms = Uniforms(
    proj,
    Vec4f0(1, 0, 0, 1),
    Vec4f0(1, 0, 1, 1)
)

atlas = get_texture_atlas(w);

text = Text(
    GLSprite2[],
    Char[],
    atlas,
    [0],
    Visualize.defaultfont(),
    45f0,
    0f0,
    Vec4f0(0f0, 0f0, 0f0, 1f0),
    Point2f0(0),
    1.5
);

print(text, 'c')
print(text, " lololol")
text.color = Vec4f0(1f0, 0f0, 0f0, 1f0)
print(text, " Muashahahaha@")

uniform_buff = UniformBuffer(uniforms);
vbo = VertexArray(text.data, face_type = Face{1, OffsetInteger{1, GLint}});

draw_particles = GLRasterizer(
    vbo, (uniform_buff, atlas.images),
    vertex_main, fragment_main;
    geometry_main = geometry_main
)
glDisable(GL_DEPTH_TEST)
glClearColor(1, 1, 1, 0)
GLAbstraction.enabletransparency()

while isopen(w)
    GLWindow.poll_glfw()
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    draw_particles(vbo, (uniform_buff, atlas.images))
    GLWindow.swapbuffers(w)
end
GLFW.DestroyWindow(w)


text:
"""
Dash port of Shiny faithful example:

https://shiny.rstudio.com/gallery/faithful.html

Note: several features have been removed since figure_factory
does not exist for plotly.jl. This is a heavily simplified
version of the original R shiny implementation.
"""

using Dash, DashBootstrapComponents, DashCoreComponents, DashHtmlComponents
using CSV, DataFrames, HTTP

# import plotly.figure_factory as ff

# data_path = "https://cdn.opensource.faculty.ai/old-faithful/data.csv";
# DATA = DataFrame(CSV.File(HTTP.get(data_path).body));

DATA = DataFrame(CSV.File("data.csv"));

app = dash(external_stylesheets = [dbc_themes.BOOTSTRAP]);

dropdown = dbc_formgroup([
    dbc_label("Number of bins in histogram (approximate):"),
    dcc_dropdown(
        id = "dropdown",
        options = [Dict("label" => n, "value" => n) for n in [10, 20, 35, 50]],
        value = 20,
    ),
])

app.layout = dbc_container([
    html_h1("Old Faithful eruption data"),
    html_hr(),
    dbc_row([dbc_col(dropdown)]),
    html_br(),
    dcc_graph(id = "graph"),
])

callback!(app, Output("graph", "figure"), Input("dropdown", "value")) do dropdown_value
    fig_data = [(x = DATA[:, "eruptions"], nbinsx = dropdown_value, type = "histogram")]

    fig_layout = (
        title = "Geyser eruption duration",
        showlegend = false,
        xaxis = (title = "Duration (minutes)",),
        yaxis = (title = "Density",),
    )
    return (data = fig_data, layout = fig_layout)

end;

run_server(app, "0.0.0.0", 8888, debug = true);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
function gradient(H::LocalHamiltonian, Ψρs::FiniteCMPSData, HL = nothing, HR = nothing;
                    gradindices = 1:length(nodes(Ψρs[1].Q)),
                    left_boundary = :fixed,
                    right_boundary = :fixed,
                    kwargs...)

    @assert first(gradindices) == 1
    @assert last(gradindices) == length(nodes(Ψρs[1].Q))
    @assert issorted(gradindices)

    Ψ, ρL, ρR = Ψρs
    if isnothing(HL)
        HL, = leftenv(H, Ψρs; kwargs...)
    end
    if isnothing(HR)
        HR, = rightenv(H, Ψρs; kwargs...)
    end
    a, b = domain(Ψ)
    Q, Rs, vL, vR = Ψ

    Z = dot(vR, ρL(b)*vR)
    E = dot(vR, HL(b)*vR)/Z

    gradQ, gradRs, grad∂Rs = _gradient(H, Ψρs, HL, HR, Z, E)
    ∇Q, ∇Rs = _project(gradQ, gradRs, grad∂Rs; gradindices = gradindices)

    ∇vL = left_boundary == :fixed ? zero(vL) : (HR(a)*vL - E*ρR(a)*vL)/Z
    ∇vR = right_boundary == :fixed ? zero(vR) : (HL(b)*vR - E*ρL(b)*vR)/Z

    return FiniteCMPSTangent(Ψ, ∇Q, ∇Rs, ∇vL, ∇vR, gradindices)
end

# Compute the gradient as a continuous function, represented as Piecewise{TaylorSeries}
function _gradient(H::LocalHamiltonian, Ψρs::FiniteCMPSData, HL, HR, Z, E)
    Ψ, ρL, ρR = Ψρs
    Q, Rs, vL, vR = Ψ

    # gradQ = ∑(coeff * localgradientQ)  +  HL*ρR + ρL*HR
    gradQ = zero(ρL)
    for (coeff, op) in zip(coefficients(H.h), operators(H.h))
        if op isa ContainsDifferentiatedCreation
            if coeff isa Number
                axpy!(coeff/Z, localgradientQ(op, Q, Rs, ρL, ρR), gradQ)
            else
                mul!(gradQ, coeff, localgradientQ(op, Q, Rs, ρL, ρR), 1/Z, 1)
            end
        end
    end
    mul!(gradQ, HL, ρR, 1/Z, 1)
    mul!(gradQ, ρL, HR, 1/Z, 1)
    mul!(gradQ, ρL, ρR, -E/Z, 1)

    # gradR = ∑(coeff * localgradientR)  +  HL*R*ρR + ρL*R*HR
    # however, we treat R and ∂R as independent variables at first in computing the gradient
    gradRs = map(R->zero(ρL), Rs)
    for (coeff, op) in zip(coefficients(H.h), operators(H.h))
        if coeff isa Number
            axpy!.(coeff/Z, localgradientRs(op, Q, Rs, ρL, ρR), gradRs)
        else
            mul!.(gradRs, (coeff,), localgradientRs(op, Q, Rs, ρL, ρR), 1/Z, 1)
        end
    end
    RρRs = Rs .* (ρR,)
    mul!.(gradRs, (HL,), RρRs, 1/Z, 1)
    mul!.(gradRs, (ρL,), Rs .* (HR,), 1/Z, 1)
    mul!.(gradRs, (ρL,), RρRs, -E/Z, 1)

    grad∂Rs = map(R->zero(ρL), Rs)
    for (coeff, op) in zip(coefficients(H.h), operators(H.h))
        if op isa ContainsDifferentiatedCreation
            if coeff isa Number
                axpy!.(coeff/Z, localgradient∂Rs(op, Q, Rs, ρL, ρR), grad∂Rs)
            else
                mul!.(grad∂Rs, (coeff,), localgradient∂Rs(op, Q, Rs, ρL, ρR), 1/Z, 1)
            end
        end
    end

    return gradQ, gradRs, grad∂Rs
end

function Qmetric(grid)
    N = length(grid)
    dv = zeros(N)
    ev = zeros(N-1)
    @inbounds for i = 1:N
        dv[i] = (grid[min(N, i+1)] - grid[max(i-1, 1)])/3
        ev[i] = (grid[i+1] - grid[i])/6
    end
    return SymTridiagonal(dv, ev)
end

function Rmetric(grid)
    g = Qmetric(grid)
    g.dv[1] = g.dv[end] = 1
    g.ev[1] = g.ev[end] = 0
    return g
end


text:
export PC, PCCreate, PCSetFromOptions, KSPGetPC, KSPSetPC, PCSetType, PCGetType, PCFactorSetUseInPlace, PCFactorGetUseInPlace, PCSetReusePreconditioner, PCGetReusePreconditioner, PCFactorSetAllowDiagonalFill, PCFactorGetAllowDiagonalFill, PCFactorSetLevels, PCFactorGetLevels, PCSetReusePreconditioner, PCGetReusePreconditioner, PCBJacobiGetSubKSP, PCFactorSetFill, PCJacobiSetType, PCJacobiGetType

# Shell PC
export PCShellSetApply, PCShellSetApplyTranspose, PCShellSetSetUp,
       PCShellSetContext, PCShellGetContext

# developer PC interface
export PCApply, PCSetUp, PCApplyTranspose, PCApplyTransposeExists

# preconditioner contex
# the KSP object creates the PC contex, so we don't provide a constructor
type PC
  pobj::Ptr{Void}
end

global const PC_NULL = PC(C_NULL)

function PCCreate(comm::MPI.Comm)
  arg2 = Ref{Ptr{Void}}()
  ccall((:PCCreate,petsc),PetscErrorCode,(comm_type,Ptr{Ptr{Void}}),comm,arg2)

  return PC(arg2[])
end

function PetscDestroy(arg1::PC)
  if arg1.pobj != C_NULL
    ccall((:PCDestroy,petsc),PetscErrorCode,(Ptr{Ptr{Void}},),&arg1.pobj)
  end
end

function PCSetFromOptions(arg1::PC)
    ccall((:PCSetFromOptions,petsc),PetscErrorCode,(Ptr{Void},),arg1.pobj)
end

function KSPGetPC(ksp::KSP)
    arr = Array(Ptr{Void}, 1)
    ccall((:KSPGetPC,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),ksp.pobj, arr)
    return PC(arr[1])
end

function KSPSetPC(ksp::KSP, pc::PC)
    ccall((:KSPSetPC,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),ksp.pobj, pc.pobj)
    return nothing
end

function PCSetType(pc::PC, pctype::PCType)
    ccall((:PCSetType,petsc),PetscErrorCode,(Ptr{Void},Cstring), pc.pobj, pctype)
end

function PCGetType(pc::PC)
    arr = Array(Ptr{UInt8}, 1)
    ccall((:PCGetType,petsc),PetscErrorCode,(Ptr{Void},Ptr{Ptr{UInt8}}), pc.pobj, arr)
    return bytestring(arr[1])
end

function PCShellSetApply(arg1::PC,arg2::Ptr{Void})
    ccall((:PCShellSetApply,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),arg1.pobj,arg2)
end

function PCShellSetApplyTranspose(arg1::PC,arg2::Ptr{Void})
    ccall((:PCShellSetApplyTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),arg1.pobj,arg2)
end

function PCShellSetSetUp(arg1::PC,arg2::Ptr{Void})
    ccall((:PCShellSetSetUp,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),arg1.pobj,arg2)
end

function PCShellSetContext(arg1::PC,arg2::Ptr{Void})
    ccall((:PCShellSetContext,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),arg1.pobj,arg2)
end

function PCShellGetContext(arg1::PC)
    arg2 = Ref{Ptr{Void}}()
    ccall((:PCShellGetContext,petsc),PetscErrorCode,(Ptr{Void},Ref{Ptr{Void}}),arg1.pobj,arg2)

    return arg2[]
end

function PCFactorSetUseInPlace(pc::PC, arg2::PetscBool)
    ccall((:PCFactorSetUseInPlace,petsc),PetscErrorCode,(Ptr{Void},PetscBool),pc.pobj, arg2)
end

function PCFactorGetUseInPlace(pc::PC)
    arr = Array(PetscBool, 1)
    ccall((:PCFactorGetUseInPlace,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscBool}), pc.pobj, arr)
    return arr[1]
end

function PCSetReusePreconditioner(pc::PC,arg2::PetscBool)
    ccall((:PCSetReusePreconditioner,petsc),PetscErrorCode,(Ptr{Void},PetscBool), pc.pobj, arg2)
end

function PCGetReusePreconditioner(pc::PC)
    arr = Array(PetscBool, 1)
    ccall((:PCGetReusePreconditioner,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscBool}), pc.pobj, arr)
    return arr[1]
end

function PCFactorSetAllowDiagonalFill(pc::PC,arg2::PetscBool)
    ccall((:PCFactorSetAllowDiagonalFill,petsc),PetscErrorCode,(Ptr{Void},PetscBool), pc.pobj, arg2)
end

function PCFactorGetAllowDiagonalFill(pc::PC)
   arr = Array(PetscBool, 1)
    ccall((:PCFactorGetAllowDiagonalFill,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscBool}), pc.pobj, arr)
    return arr[1]
end

function PCFactorSetLevels(pc::PC,arg2::PetscInt)
    ccall((:PCFactorSetLevels,petsc),PetscErrorCode,(Ptr{Void}, PetscInt), pc.pobj, arg2)
end

function PCFactorGetLevels(pc::PC)
    arr = Array(PetscInt, 1)
    ccall((:PCFactorGetLevels,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt}),pc.pobj, arr)
    return arr[1]
end

function PCSetReusePreconditioner(pc::PC, arg2::PetscBool)
    ccall((:PCSetReusePreconditioner,petsc),PetscErrorCode,(Ptr{Void}, PetscBool), pc.pobj, arg2)
end

function PCGetReusePreconditioner(pc::PC)
    arr = Array(PetscBool, 1)
    ccall((:PCGetReusePreconditioner,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscBool}), pc.pobj, arr)
    return arr[1]
end

function PCBJacobiGetSubKSP(pc::PC)
    n_local_arr = Array(PetscInt, 1)
    first_local = Array(PetscInt, 1)
    ksp_ptrarr = Array(Ptr{Ptr{Void}}, 1)
    ccall((:PCBJacobiGetSubKSP,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt},Ptr{Ptr{Ptr{Void}}}), pc.pobj, n_local_arr, first_local, ksp_ptrarr)

    n_local = n_local_arr[1]
    ksp_ptrarr2 = pointer_to_array(ksp_ptrarr[1], n_local)

    ksp_arr = Array(KSP, n_local)
    for i=1:n_local
      ksp_arr[i] = KSP(ksp_ptrarr2[i])
    end

    return n_local, first_local[1], ksp_arr
end

function PCFactorSetFill(pc::PC, fill::PetscReal)
    ccall((:PCFactorSetFill,petsc),PetscErrorCode,(Ptr{Void}, PetscReal), pc.pobj, fill)
end

function PCJacobiSetType(pc::PC, jacobitype::PCJacobiType)
    ccall((:PCJacobiSetType,petsc),PetscErrorCode,(Ptr{Void}, PCJacobiType), pc.pobj, jacobitype)
end

function PCJacobiGetType(pc::PC)
    arr = Array(PCJacobiType, 1)
    ccall((:PCJacobiGetType,petsc),PetscErrorCode,(Ptr{Void},Ptr{PCJacobiType}),pc.pobj, arr)
    return arr[1]
end

function PCApply(arg1::PC,arg2::PetscVec,arg3::PetscVec)
  println("pc.pobj = ", arg1.pobj)
    ccall((:PCApply,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}),arg1.pobj ,arg2.pobj, arg3.pobj)
end

function PCApplyTranspose(arg1::PC,arg2::PetscVec,arg3::PetscVec)
    ccall((:PCApplyTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}),arg1.pobj ,arg2.pobj, arg3.pobj)
end

function PCApplyTransposeExists(arg1::PC)
    arg2 = Ref{PetscBool}()
    ccall((:PCApplyTransposeExists,petsc),PetscErrorCode,(PC,Ptr{PetscBool}),arg1,arg2)

    return arg2[] == 1
end

function PCSetUp(arg1::PC)
    ccall((:PCSetUp,petsc),PetscErrorCode,(Ptr{Void},),arg1.pobj)
end


text:
@testset verbose=true "d" begin

    @testset "Compare d expressions ($T)" for T in [BigFloat, Float64, Float32]
        # This just compares the two versions of the `d` function from test_utilities
        # to ensure that later tests that use those functions are reliable
        tol = ifelse(T === BigFloat, 100, 1) * 30eps(T)
        for β in βrange(T)
            expiβ = cis(β)
            for ℓₘₐₓ in 0:2  # 2 is the max explicitly coded ℓ
                for m′ₘₐₓ in 0:ℓₘₐₓ
                    for n in 0:ℓₘₐₓ
                        for m′ in -min(n, m′ₘₐₓ):min(n, m′ₘₐₓ)
                            for m in -n:n
                                d_expl = ExplicitWignerMatrices.d_explicit(n, m′, m, expiβ)
                                d_form = ExplicitWignerMatrices.d_formula(n, m′, m, expiβ)
                                @test d_expl ≈ d_form atol=tol rtol=tol
                            end
                        end
                    end
                end
            end
        end
    end

    @testset "Compare d to formulaic d ($T)" for T in [BigFloat, Float64, Float32]
        # Now, we're ready to check that d_{n}^{m′,m}(β) matches the expected values
        # for a range of β values
        tol = ifelse(T === BigFloat, 100, 1) * 30eps(T)
        for β in βrange(T)
            expiβ = cis(β)
            for ℓₘₐₓ in 0:4
                d, H_rec_coeffs = dprep(ℓₘₐₓ, T)
                d!(d, expiβ, ℓₘₐₓ, H_rec_coeffs)
                for n in 0:ℓₘₐₓ
                    for m′ in -n:n
                        for m in -n:n
                            d_formula = ExplicitWignerMatrices.d_formula(n, m′, m, expiβ)
                            d_recurrence = d[WignerDindex(n, m′, m)]
                            @test d_formula ≈ d_recurrence atol=tol rtol=tol
                        end
                    end
                end
            end
        end
    end

    @testset "Test d signatures ($T)" for T in [BigFloat, Float64, Float32]
        # 1 d!(𝔡, expiβ::Complex{T}, ℓₘₐₓ) where {T<:Real}
        # 2 d!(𝔡, β::T, ℓₘₐₓ, (a,b,d)) where {T<:Real}
        # 3 d!(𝔡, β::T, ℓₘₐₓ) where {T<:Real}
        # 4 d(expiβ::Complex{T}, ℓₘₐₓ) where {T<:Real}
        # 5 d(β::Real, ℓₘₐₓ)
        ℓₘₐₓ = 8
        for β in βrange(T)
            expiβ = cis(β)
            dA = d(β, ℓₘₐₓ)
            dB = d(expiβ, ℓₘₐₓ)
            @test array_equal(dA, dB)
            dB .= 0
            d!(dB, β, ℓₘₐₓ)
            @test array_equal(dA, dB)
            dB .= 0
            d!(dB, expiβ, ℓₘₐₓ)
            @test array_equal(dA, dB)
            dB .= 0
            H_rec_coeffs = H_recursion_coefficients(ℓₘₐₓ, T)
            d!(dB, β, ℓₘₐₓ, H_rec_coeffs)
            @test array_equal(dA, dB)
        end
    end

end


text:
module VortexDistributions

using Reexport
@reexport using Test
@reexport using JLD2
@reexport using Parameters
@reexport using Interpolations
Core.eval(Main, :(using Interpolations))
@reexport using SpecialFunctions
@reexport using LinearAlgebra
@reexport using ToeplitzMatrices
@reexport using SparseArrays
@reexport using FFTW
@reexport using FileIO
@reexport using ProgressMeter

const Λ = 0.8249
export Field, Torus, Sphere
export Vortex, CoreShape, Ansatz, Exact, ScalarVortex
export PointVortex, rawData, uniform, randcharge
export randPointVortex, randScalarVortex, randVortex
export vortex!, findvortices, Thetad, periodic_dipole!
export foundNear, randVortexField
export Basis, Oscillator, hermite, hermite_polar
export index, spectrum, qnumbers, filter, slowpolar, init_polar, polar, filterH

# export ψi, ψa
# export findwhere, findvortices_jumps, findvortices_grid, findvortices_interp
# export removeedgevortices, corezoom
# export phasejumps, phasejumps!, unwrap, unwrap!
# export gpecore_exact, chebdif, getChebDMatrix, getChebD2Matrix

Interpolations.interpolate(A::Array{Complex{Float64},2}) = interpolate(real(A))+im*interpolate(imag(A))

include("types.jl")
include("pointvortex.jl")
include("detection.jl")
include("creation.jl")
include("analysis.jl")

@load joinpath(@__DIR__,"exactcore.jld2") ψi
@load joinpath(@__DIR__,"ansatzcore.jld2") ψa

end # module


text:
using Clang.Generators
using open62541_jll

cd(@__DIR__)

include_dir = joinpath(open62541_jll.artifact_dir, "include") |> normpath
open62541_header = joinpath(include_dir, "open62541.h") |> normpath
@assert isfile(open62541_header)

function write_generated_defs(generated_defs_dir::String, open62541_header::String)
    type_s = """
    # Vector of all UA types; Generated with: type_names = [Symbol("UA_", unsafe_string(unsafe_load(type_ptr).typeName)) for type_ptr in UA_TYPES_PTRS] TODO: Automate this step
    const type_names = Symbol[:UA_Boolean, :UA_SByte, :UA_Byte, :UA_Int16, :UA_UInt16, :UA_Int32, :UA_UInt32, :UA_Int64, :UA_UInt64, :UA_Float, :UA_Double, :UA_String, :UA_DateTime, :UA_Guid, :UA_ByteString, :UA_XmlElement, :UA_NodeId, :UA_ExpandedNodeId, :UA_StatusCode, :UA_QualifiedName, :UA_LocalizedText, :UA_ExtensionObject, :UA_DataValue, :UA_Variant, :UA_DiagnosticInfo, :UA_ViewAttributes, :UA_XVType, :UA_ElementOperand, :UA_VariableAttributes, :UA_EnumValueType, :UA_EventFieldList, :UA_MonitoredItemCreateResult, :UA_EUInformation, :UA_ServerDiagnosticsSummaryDataType, :UA_ContentFilterElementResult, :UA_LiteralOperand, :UA_MessageSecurityMode, :UA_UtcTime, :UA_UserIdentityToken, :UA_X509IdentityToken, :UA_MonitoredItemNotification, :UA_StructureType, :UA_ResponseHeader, :UA_SignatureData, :UA_ModifySubscriptionResponse, :UA_NodeAttributes, :UA_ActivateSessionResponse, :UA_EnumField, :UA_VariableTypeAttributes, :UA_CallMethodResult, :UA_MonitoringMode, :UA_SetMonitoringModeResponse, :UA_BrowseResultMask, :UA_RequestHeader, :UA_MonitoredItemModifyResult, :UA_CloseSecureChannelRequest, :UA_NotificationMessage, :UA_CreateSubscriptionResponse, :UA_EnumDefinition, :UA_AxisScaleEnumeration, :UA_BrowseDirection, :UA_CallMethodRequest, :UA_ReadResponse, :UA_TimestampsToReturn, :UA_NodeClass, :UA_ObjectTypeAttributes, :UA_SecurityTokenRequestType, :UA_CloseSessionResponse, :UA_SetPublishingModeRequest, :UA_IssuedIdentityToken, :UA_DeleteMonitoredItemsResponse, :UA_ApplicationType, :UA_BrowseNextRequest, :UA_ModifySubscriptionRequest, :UA_BrowseDescription, :UA_SignedSoftwareCertificate, :UA_BrowsePathTarget, :UA_WriteResponse, :UA_AddNodesResult, :UA_AddReferencesItem, :UA_DeleteReferencesResponse, :UA_RelativePathElement, :UA_SubscriptionAcknowledgement, :UA_TransferResult, :UA_CreateMonitoredItemsResponse, :UA_DeleteReferencesItem, :UA_WriteValue, :UA_DataTypeAttributes, :UA_TransferSubscriptionsResponse, :UA_AddReferencesResponse, :UA_DeadbandType, :UA_DataChangeTrigger, :UA_BuildInfo, :UA_FilterOperand, :UA_MonitoringParameters, :UA_DoubleComplexNumberType, :UA_DeleteNodesItem, :UA_ReadValueId, :UA_CallRequest, :UA_RelativePath, :UA_DeleteNodesRequest, :UA_MonitoredItemModifyRequest, :UA_UserTokenType, :UA_AggregateConfiguration, :UA_LocaleId, :UA_UnregisterNodesResponse, :UA_ContentFilterResult, :UA_UserTokenPolicy, :UA_DeleteMonitoredItemsRequest, :UA_SetMonitoringModeRequest, :UA_Duration, :UA_ReferenceTypeAttributes, :UA_GetEndpointsRequest, :UA_CloseSecureChannelResponse, :UA_ViewDescription, :UA_SetPublishingModeResponse, :UA_StatusChangeNotification, :UA_StructureField, :UA_NodeAttributesMask, :UA_EventFilterResult, :UA_MonitoredItemCreateRequest, :UA_ComplexNumberType, :UA_Range, :UA_DataChangeNotification, :UA_Argument, :UA_TransferSubscriptionsRequest, :UA_ChannelSecurityToken, :UA_ServerState, :UA_EventNotificationList, :UA_AnonymousIdentityToken, :UA_FilterOperator, :UA_AggregateFilter, :UA_RepublishResponse, :UA_DeleteSubscriptionsResponse, :UA_RegisterNodesRequest, :UA_StructureDefinition, :UA_MethodAttributes, :UA_UserNameIdentityToken, :UA_UnregisterNodesRequest, :UA_OpenSecureChannelResponse, :UA_SetTriggeringResponse, :UA_SimpleAttributeOperand, :UA_RepublishRequest, :UA_RegisterNodesResponse, :UA_ModifyMonitoredItemsResponse, :UA_DeleteSubscriptionsRequest, :UA_RedundancySupport, :UA_BrowsePath, :UA_ObjectAttributes, :UA_PublishRequest, :UA_FindServersRequest, :UA_ReferenceDescription, :UA_CreateSubscriptionRequest, :UA_CallResponse, :UA_DeleteNodesResponse, :UA_ModifyMonitoredItemsRequest, :UA_ServiceFault, :UA_PublishResponse, :UA_CreateMonitoredItemsRequest, :UA_OpenSecureChannelRequest, :UA_CloseSessionRequest, :UA_SetTriggeringRequest, :UA_BrowseResult, :UA_AddReferencesRequest, :UA_AddNodesItem, :UA_ServerStatusDataType, :UA_BrowseNextResponse, :UA_AxisInformation, :UA_ApplicationDescription, :UA_ReadRequest, :UA_ActivateSessionRequest, :UA_BrowsePathResult, :UA_AddNodesRequest, :UA_BrowseRequest, :UA_WriteRequest, :UA_AddNodesResponse, :UA_AttributeOperand, :UA_DataChangeFilter, :UA_EndpointDescription, :UA_DeleteReferencesRequest, :UA_TranslateBrowsePathsToNodeIdsRequest, :UA_FindServersResponse, :UA_CreateSessionRequest, :UA_ContentFilterElement, :UA_TranslateBrowsePathsToNodeIdsResponse, :UA_BrowseResponse, :UA_CreateSessionResponse, :UA_ContentFilter, :UA_GetEndpointsResponse, :UA_EventFilter]

    # Julia types corresponding to the UA types in vector type_names
    const julia_types = DataType[Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64, UA_String, Int64, UA_Guid, UA_String, UA_String, UA_NodeId, UA_ExpandedNodeId, UInt32, UA_QualifiedName, UA_LocalizedText, UA_ExtensionObject, UA_DataValue, UA_Variant, UA_DiagnosticInfo, UA_ViewAttributes, UA_XVType, UA_ElementOperand, UA_VariableAttributes, UA_EnumValueType, UA_EventFieldList, UA_MonitoredItemCreateResult, UA_EUInformation, UA_ServerDiagnosticsSummaryDataType, UA_ContentFilterElementResult, UA_LiteralOperand, UA_MessageSecurityMode, Int64, UA_UserIdentityToken, UA_X509IdentityToken, UA_MonitoredItemNotification, UA_StructureType, UA_ResponseHeader, UA_SignatureData, UA_ModifySubscriptionResponse, UA_NodeAttributes, UA_ActivateSessionResponse, UA_EnumField, UA_VariableTypeAttributes, UA_CallMethodResult, UA_MonitoringMode, UA_SetMonitoringModeResponse, UA_BrowseResultMask, UA_RequestHeader, UA_MonitoredItemModifyResult, UA_CloseSecureChannelRequest, UA_NotificationMessage, UA_CreateSubscriptionResponse, UA_EnumDefinition, UA_AxisScaleEnumeration, UA_BrowseDirection, UA_CallMethodRequest, UA_ReadResponse, UA_TimestampsToReturn, UA_NodeClass, UA_ObjectTypeAttributes, UA_SecurityTokenRequestType, UA_CloseSessionResponse, UA_SetPublishingModeRequest, UA_IssuedIdentityToken, UA_DeleteMonitoredItemsResponse, UA_ApplicationType, UA_BrowseNextRequest, UA_ModifySubscriptionRequest, UA_BrowseDescription, UA_SignedSoftwareCertificate, UA_BrowsePathTarget, UA_WriteResponse, UA_AddNodesResult, UA_AddReferencesItem, UA_DeleteReferencesResponse, UA_RelativePathElement, UA_SubscriptionAcknowledgement, UA_TransferResult, UA_CreateMonitoredItemsResponse, UA_DeleteReferencesItem, UA_WriteValue, UA_DataTypeAttributes, UA_TransferSubscriptionsResponse, UA_AddReferencesResponse, UA_DeadbandType, UA_DataChangeTrigger, UA_BuildInfo, Ptr{Nothing}, UA_MonitoringParameters, UA_DoubleComplexNumberType, UA_DeleteNodesItem, UA_ReadValueId, UA_CallRequest, UA_RelativePath, UA_DeleteNodesRequest, UA_MonitoredItemModifyRequest, UA_UserTokenType, UA_AggregateConfiguration, UA_String, UA_UnregisterNodesResponse, UA_ContentFilterResult, UA_UserTokenPolicy, UA_DeleteMonitoredItemsRequest, UA_SetMonitoringModeRequest, Float64, UA_ReferenceTypeAttributes, UA_GetEndpointsRequest, UA_CloseSecureChannelResponse, UA_ViewDescription, UA_SetPublishingModeResponse, UA_StatusChangeNotification, UA_StructureField, UA_NodeAttributesMask, UA_EventFilterResult, UA_MonitoredItemCreateRequest, UA_ComplexNumberType, UA_Range, UA_DataChangeNotification, UA_Argument, UA_TransferSubscriptionsRequest, UA_ChannelSecurityToken, UA_ServerState, UA_EventNotificationList, UA_AnonymousIdentityToken, UA_FilterOperator, UA_AggregateFilter, UA_RepublishResponse, UA_DeleteSubscriptionsResponse, UA_RegisterNodesRequest, UA_StructureDefinition, UA_MethodAttributes, UA_UserNameIdentityToken, UA_UnregisterNodesRequest, UA_OpenSecureChannelResponse, UA_SetTriggeringResponse, UA_SimpleAttributeOperand, UA_RepublishRequest, UA_RegisterNodesResponse, UA_ModifyMonitoredItemsResponse, UA_DeleteSubscriptionsRequest, UA_RedundancySupport, UA_BrowsePath, UA_ObjectAttributes, UA_PublishRequest, UA_FindServersRequest, UA_ReferenceDescription, UA_CreateSubscriptionRequest, UA_CallResponse, UA_DeleteNodesResponse, UA_ModifyMonitoredItemsRequest, UA_ServiceFault, UA_PublishResponse, UA_CreateMonitoredItemsRequest, UA_OpenSecureChannelRequest, UA_CloseSessionRequest, UA_SetTriggeringRequest, UA_BrowseResult, UA_AddReferencesRequest, UA_AddNodesItem, UA_ServerStatusDataType, UA_BrowseNextResponse, UA_AxisInformation, UA_ApplicationDescription, UA_ReadRequest, UA_ActivateSessionRequest, UA_BrowsePathResult, UA_AddNodesRequest, UA_BrowseRequest, UA_WriteRequest, UA_AddNodesResponse, UA_AttributeOperand, UA_DataChangeFilter, UA_EndpointDescription, UA_DeleteReferencesRequest, UA_TranslateBrowsePathsToNodeIdsRequest, UA_FindServersResponse, UA_CreateSessionRequest, UA_ContentFilterElement, UA_TranslateBrowsePathsToNodeIdsResponse, UA_BrowseResponse, UA_CreateSessionResponse, UA_ContentFilter, UA_GetEndpointsResponse, UA_EventFilter]

    # Vector of types that are ambiguously defined via typedef and are not to be used as default type
    types_ambiguous_denylist = [:UA_Duration, :UA_ByteString, :UA_XmlElement, :UA_LocaleId, :UA_DateTime, :UA_UtcTime, :UA_StatusCode]

    """

    inlined_funcs = """

    # Vector of all inlined function names listed in the amalgamated open62541 header file
    const inlined_funcs = $(extract_inlined_funcs(open62541_header))
    """

    data_UA_Client = """

    # UA_Client_ functions data
    const attributes_UA_Client_Service = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}\s(UA_Client_Service_(\w*))\((?:[\s\S]*?)\)(?:[\s\S]*?)UA_\S*", open62541_header))

    const attributes_UA_Client_read = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}\s(UA_Client_read(\w*)Attribute)\((?:[\s\S]*?,\s*){2}(\S*)", open62541_header))

    const attributes_UA_Client_write = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}\s(UA_Client_write(\w*)Attribute)\((?:[\s\S]*?,\s*){2}const\s(\S*)", open62541_header))

    const attributes_UA_Client_read_async = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}\s(UA_Client_read(\w*)Attribute_async)\([\s\S]+?\)[\s\S]+?{[\s\S]+?__UA_Client_readAttribute_async\s*\([\s\S]+?&UA_TYPES\[([\S]+?)\]", open62541_header))

    const attributes_UA_Client_write_async = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}\s(UA_Client_write(\w*)Attribute_async)\s*\(UA_Client\s*\*client,\s*const\s*UA_NodeId\s*nodeId,\s*const\s*(\S*)", open62541_header))
    """

    data_UA_Server = """

        # UA_Server_ functions data
        const attributes_UA_Server_read = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}\s(UA_Server_read(\w*))\s*\(UA_Server\s*\*server,\s*const\s*UA_NodeId\s*nodeId,\s*(\S*)", open62541_header))

        const attributes_UA_Server_write = $(extract_header_data(r"UA_INLINE[\s\S]{0,50}(UA_Server_write(\w*))\s*\(UA_Server\s*\*server,\s*const\s*UA_NodeId\s*nodeId,\s*const (\S*)", open62541_header))
        """

        # const attributes_UA_Client_read = $(extract_header_data(r"[;}^][\s\S]{0,50}UA_INLINE[\s\S]{0,50}\s(UA_Client_read(\w*)Attribute)\((?:[\s\S]*?,\s*){2}(\S*)", open62541_header))

        # const attributes_UA_Client_write = $(extract_header_data(r"[;}^][\s\S]{0,50}UA_INLINE[\s\S]{0,50}\s(UA_Client_write(\w*)Attribute)\((?:[\s\S]*?,\s*){2}const\s(\S*)", open62541_header))

        # const attributes_UA_Client_read_async = $(extract_header_data(r"[;}^][\s\S]{0,50}UA_INLINE[\s\S]{0,50}\s(UA_Client_read(\w*)Attribute_async)\([\s\S]+?\)[\s\S]+?{[\s\S]+?__UA_Client_readAttribute_async\s*\([\s\S]+?&UA_TYPES\[([\S]+?)\]", open62541_header))

        # const attributes_UA_Client_write_async = $(extract_header_data(r"[;}^][\s\S]{0,50}UA_INLINE[\s\S]{0,50}\s(UA_Client_write(\w*)Attribute_async)\s*\(UA_Client\s*\*client,\s*const\s*UA_NodeId\s*nodeId,\s*const\s*(\S*)", open62541_header))

    open(generated_defs_dir, "w") do f
        write(f, type_s)
        write(f, inlined_funcs)
        write(f, data_UA_Client)
        write(f, data_UA_Server)
    end
end

function extract_inlined_funcs(open62541_header::String)
    regex_inlined = r"UA_INLINE[\s]+(?:[\w\*]+[\s]*[\s\S]){0,3}((?:__)?UA_[\w]+)\("
    inlined_funcs = String[]
    open(open62541_header, "r") do f
        data = read(f, String)
        append!(inlined_funcs, vcat(getfield.(collect(eachmatch(regex_inlined, data)), :captures)...)) # Extract inlined functions from header file
    end
    return inlined_funcs
end

function extract_header_data(regex::Regex, open62541_header::String)
    f = open(open62541_header, "r")
    data = read(f, String)
    close(f)
    all_data = getfield.(collect(eachmatch(regex, data)), :captures) # Extract inlined functions from header file
    return all_data
end

# Write static definitions to file generated_defs.jl
write_generated_defs(joinpath(@__DIR__, "../src/generated_defs.jl"), open62541_header)

# Load options from generator.toml
options = load_options(joinpath(@__DIR__, "generator.toml"))

# Extract all inlined functions and move them to codegen denylist
append!(options["general"]["printer_denylist"], extract_inlined_funcs(open62541_header))

# Add compiler flags
args = get_default_args()
push!(args, "-I$include_dir")
push!(args, "-std=c99")

# Create context
ctx = create_context([open62541_header], args, options)

# Run generator
build!(ctx)


text:
using EconPDEs

struct ConstraintModel

  # utility function
  γ::Float64
  ψ::Float64
  ρ::Float64
  δ::Float64

  # FixedShare
  α::Float64
  # Proportion Fixed Share
  ν::Float64

  # consumption
  μ::Float64
  σ::Float64
end

function ConstraintModel(;γ  = 1.5, ψ = 0.7, ρ = 0.001, δ = 0.02, α = 1.0, ν = 0.01, μ = 0.02, σ = 0.041)
  ConstraintModel(γ , ψ, γ, ψ, ρ, δ, ν, μ, σ)
end

function initialize_stategrid(m::ConstraintModel; n = 200)
  OrderedDict(:x => range(0.0, stop = 1.0, length = n))
end

function initialize_y(m::ConstraintModel, stategrid::OrderedDict)
    x = ones(length(stategrid[:x]))
    OrderedDict(:pA => x, :pB => x)
end

# p = ρ / γ + (1- 1/γ) * (r + α * κ + α σ^2)

function (m::ConstraintModel)(state::NamedTuple, y::NamedTuple)
  γ = m.γ ; ψ = m.ψ ; ρ = m.ρ ; δ = m.δ ; ν = m.ν ; μ = m.μ ; σ = m.σ
  x = state.x
  # pA is wealth / consumption ratio of agent A
  # pB is wealth / consumption ratio of agent B
  pA = y.pA ; pAx = y.pAx ; pAxx = y.pAxx ; pB = y.pB ; pBx = y.pBx ; pBxx = y.pBxx

  # Market price of risk κ
  Γ = 1 / (x / γA + (1 - x) / γB)
  p = x * pA + (1 - x) * pB
  σx = σ * x * (Γ / γA - 1) / (1 + Γ * x * (1 - x) / (γA * γB) * ((1 - γB * ψB) / (ψB - 1) * (pBx / pB) - (1 - γA * ψA) / (ψA - 1) * (pAx / pA)))
  σpA = pAx / pA * σx
  σpB = pBx / pB * σx
  σϕ1 = ϕ1x / ϕ1 * σx
  σϕ2 = ϕ2x / ϕ2 * σx
  κ = Γ * (σ - x * (1 - γA * ψA) / (γA * (ψA - 1)) * σpA - (1 - x) * (1 - γB * ψB) / (γB * (ψB - 1)) * σpB)
  σCA = κ / γA + (1 - γA * ψA) / (γA * (ψA - 1)) * σpA
  σCB = κ / γB + (1 - γB * ψB) / (γB * (ψB - 1)) * σpB

  # Interest rate r
  # A.16 Equation in Garleanu Panageas has a typo
  mcA = κ^2 * (1 + ψA) / (2 * γA) + (1 - ψA * γA) / (γA * (ψA - 1)) * κ * σpA - (1 - γA * ψA) / (2 * γA * (ψA - 1)) * σpA^2
  mcB = κ^2 * (1 + ψB) / (2 * γB) + (1 - ψB * γB) / (γB * (ψB - 1)) * κ * σpB - (1 - γB * ψB) / (2 * γB * (ψB - 1)) * σpB^2
  r =  ρ + 1 / (ψA * x  + ψB * (1 - x))  * (μ - x * mcA - (1 - x) * mcB - δ * ((νA / pA + (1 - νA) / pB) * (ϕ1 + ϕ2) - 1))
  μCA = ψA * (r - ρ) + mcA
  μCB = ψB * (r - ρ) + mcB
  μx = x * (μCA - δ - μ) + δ * νA / pA * (ϕ1 + ϕ2) - σ * σx
  μpA = pAx / pA * μx + 0.5 * pAxx / pA * σx^2
  μpB = pBx / pB * μx + 0.5 * pBxx / pB * σx^2
  μϕ1 = ϕ1x / ϕ1 * μx + 0.5 * ϕ1xx / ϕ1 * σx^2
  μϕ2 = ϕ2x / ϕ2 * μx + 0.5 * ϕ2xx / ϕ2 * σx^2

  # Market Pricing
  pAt = pA * (1 / pA + (μCA - δ) + μpA + σCA * σpA - r - κ * (σpA + σCA))
  pBt = pB * (1 / pB + (μCB - δ) + μpB + σCB * σpB - r - κ * (σpB + σCB))
  ϕ1t = ϕ1 * (B1 * ω / ϕ1 + (μ - δ - δ1) + μϕ1 + σ * σϕ1 - r - κ * (σϕ1 + σ))
  ϕ2t = ϕ2 * (B2 * ω / ϕ2 + (μ - δ - δ2) + μϕ2 + σ * σϕ2 - r - κ * (σϕ2 + σ))

  return (pAt, pBt, ϕ1t, ϕ2t), (μx, ), (μx = μx, p = p, pA = pA, pB = pB, κ = κ, r = r, σx = σx)
end

m = ConstraintModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0)


text:
export setTFSF!

function  setTFSF! end

function (sim::Simulation)(::typeof(setTFSF!), reg::Function, gridtype::GridType= p̂)
	sim(set!,sim.Q,x -> true,reg, gridtype)
    sim.activate_tfsf = true;
end
function (sim::Simulation)(::typeof(setTFSF!),dir::Direction{D},side::Low, thickness::Number,gridtype::GridType= p̂) where {D}
    sim(set!,sim.Q,x -> true,x -> x[D] <= thickness, gridtype)
    sim.activate_tfsf = true;
end

function (sim::Simulation)(::typeof(setTFSF!),dir::Direction{D},side::High, thickness::Number,gridtype::GridType= p̂) where {D}
    sim(set!,sim.Q,x -> true,x -> x[D] >= extent(sim.grid,dir)- thickness, gridtype)
    sim.activate_tfsf = true;
end
function (sim::Simulation)(::typeof(setTFSF!),dir::Direction{D}, thickness::Number,gridtype::GridType= p̂) where {D}
sim(setTFSF!,dir,HIGH, thickness,gridtype)
sim(setTFSF!,dir,LOW, thickness,gridtype)
end

function (sim::Simulation{1})(::typeof(setTFSF!), thickness::Number,gridtype::GridType= p̂)
sim(setTFSF!,x̂, thickness,gridtype)
end
function (sim::Simulation{2})(::typeof(setTFSF!), thickness::Number,gridtype::GridType= p̂)
sim(setTFSF!,x̂, thickness,gridtype)
sim(setTFSF!,ŷ, thickness,gridtype)
end
function (sim::Simulation{3})(::typeof(setTFSF!), thickness::Number,gridtype::GridType= p̂)
sim(setTFSF!,x̂, thickness,gridtype)
sim(setTFSF!,ŷ, thickness,gridtype)
sim(setTFSF!,ẑ, thickness,gridtype)
end


text:
@testset "get_skolem_normal_form" begin
    @customtest String(get_skolem_normal_form(Formula("*{5}*{1}?{2}?{3}*{4}(~P(f(1)) & (P(3) | P(4, f(2, 3))))"))) ==
        "*{5}*{1}*{4}(~P(f(1)) & (P(##2(5, 1)) | P(4, f(##1(5, 1), ##2(5, 1)))))"
    @customtest String(get_skolem_normal_form(Formula("*{x}?{y}*{z}?{w}(~P(a, w) | Q(f(x), y))"))) ==
        "*{x}*{z}(~P(a, ##2(x, z)) | Q(f(x), ##1(x)))"
end


text:
using Test
using TimeMachine, MPI, DiffEqBase

using TimeMachine.Callbacks

MPI.Initialized() || MPI.Init()

mutable struct MyCallback
    initialized::Bool
    calls::Int
    finalized::Bool
end
MyCallback() = MyCallback(false, 0, false)

function Callbacks.initialize!(cb::MyCallback, integrator)
    cb.initialized=true
end
function Callbacks.finalize!(cb::MyCallback, integrator)
    cb.finalized=true
end
function (cb::MyCallback)(integrator)
    cb.calls += 1
end

cb1 = MyCallback()
cb2 = MyCallback()
cb3 = MyCallback()
cb4 = MyCallback()
cb5 = MyCallback()

cbs = CallbackSet(
    EveryXSimulationTime(cb1, 1/4),
    EveryXSimulationTime(cb2, 1/2, atinit=true),
    EveryXSimulationSteps(cb3, 1),
    EveryXSimulationSteps(cb4, 4,  atinit=true),
    EveryXSimulationSteps(_ -> sleep(1/32), 1),
    EveryXWallTimeSeconds(cb5, 0.49, MPI.COMM_SELF)
)

solve(const_prob, LSRKEulerMethod(), dt=1/32, callback=cbs)

@test cb1.initialized
@test cb2.initialized
@test cb3.initialized
@test cb4.initialized
@test cb5.initialized

@test cb1.calls == 4
@test cb2.calls == 3
@test cb3.calls == 32
@test cb4.calls == 9
@test cb5.calls >= 2

if isdefined(DiffEqBase, :finalize!)

    @test cb1.finalized
    @test cb2.finalized
    @test cb3.finalized
    @test cb4.finalized
    @test cb5.finalized
end


text:
# replace CxxStd.StdVector{T} of union to separate contents

const StdVector{T}  = cxxt"std::vector<$T>"
const StdVectorR{T} = cxxt"std::vector<$T>&"

Base.convert(::Type{StdVectorR{T}}, v::StdVector{T}) where T = icxx"(std::vector<$T>&) $v;"
#Base.convert(::Type{StdVectorR{T}}, v::StdVector{T}) where T = icxx"($(StdVectorR{T})) $v;"

#function Base.convert(::Type{cxxt"std::vector<$T>&"}, v::cxxt"std::vector<$T>") where T
#    icxx"(std::vector<$T>&) $v;"
#end

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
using StatsModels

@test StatsModels.Terms(rcopy(R"y ~ x + (1 | g)")) == StatsModels.Terms(@formula y ~ x + (1 | g))
@test StatsModels.Terms(rcopy(R"y ~ x + (1 | g) + (1 | d)")) == StatsModels.Terms(@formula y ~ x + (1 | g) + (1 | d))
@test StatsModels.Terms(rcopy(R"y ~ (a + b) * c")) == StatsModels.Terms(@formula y ~ (a + b) * c)
@test StatsModels.Terms(rcopy(R"y ~ c : (a + b)")) == StatsModels.Terms(@formula y ~ c & (a + b))

# testing association
@test StatsModels.Terms(rcopy(R"y ~ a + b + c")) == StatsModels.Terms(@formula y ~ a + b + c)
@test StatsModels.Terms(rcopy(R"y ~ a + b * c * d")) == StatsModels.Terms(@formula y ~ a + b * c * d)
@test StatsModels.Terms(rcopy(R"y ~ a : b : c")) == StatsModels.Terms(@formula y ~ a & b & c)
@test StatsModels.Terms(rcopy(R"y ~ a + b : c : d")) == StatsModels.Terms(@formula y ~ a + b & c & d)

@test rcopy(rcall(Symbol("=="), R"y ~ x + (1 | g)", RObject(@formula y ~ x + (1 | g))))
@test rcopy(rcall(Symbol("=="), R"y ~ x + (1 | g) + (1 | d)", RObject(@formula y ~ x + (1 | g) + (1 | d))))
@test rcopy(rcall(Symbol("=="), R"y ~ (a + b) * c", RObject(@formula y ~ (a + b) * c)))
@test rcopy(rcall(Symbol("=="), R"y ~ c * (a + b)", RObject(@formula y ~ c * (a + b))))

# testing association
@test rcopy(rcall(Symbol("=="), R"y ~ a + b + c", RObject(@formula y ~ a + b + c)))
@test rcopy(rcall(Symbol("=="), R"y ~ a + b * c * d", RObject(@formula y ~ a + b * c * d)))
@test rcopy(rcall(Symbol("=="), R"y ~ a : b : c", RObject(@formula y ~ a & b & c)))
@test rcopy(rcall(Symbol("=="), R"y ~ d + a : b : c", RObject(@formula y ~ d + a & b & c)))


text:
module MaxHelpingHandNoDashRefillSpring

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/NoDashRefillSpring" NoDashRefillSpring(x::Integer, y::Integer, spriteDirectory::String="objects/MaxHelpingHand/noDashRefillSpring", playerCanUse::Bool=true, ignoreLighting::Bool=false)
@mapdef Entity "MaxHelpingHand/NoDashRefillSpringRight" NoDashRefillSpringRight(x::Integer, y::Integer, spriteDirectory::String="objects/MaxHelpingHand/noDashRefillSpring", ignoreLighting::Bool=false)
@mapdef Entity "MaxHelpingHand/NoDashRefillSpringLeft" NoDashRefillSpringLeft(x::Integer, y::Integer, spriteDirectory::String="objects/MaxHelpingHand/noDashRefillSpring", ignoreLighting::Bool=false)

const placements = Ahorn.PlacementDict(
    "No Dash Refill Spring (Up) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        NoDashRefillSpring
    ),
    "No Dash Refill Spring (Left) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        NoDashRefillSpringRight
    ),
    "No Dash Refill Spring (Right) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        NoDashRefillSpringLeft
    ),
)

function Ahorn.selection(entity::NoDashRefillSpring)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 6, y - 3, 12, 5)
end

function Ahorn.selection(entity::NoDashRefillSpringLeft)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 1, y - 6, 5, 12)
end

function Ahorn.selection(entity::NoDashRefillSpringRight)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 4, y - 6, 5, 12)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashRefillSpring, room::Maple.Room) = Ahorn.drawSprite(ctx, get(entity, "spriteDirectory", "objects/MaxHelpingHand/noDashRefillSpring") * "/00.png", 0, -8)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashRefillSpringLeft, room::Maple.Room) = Ahorn.drawSprite(ctx, get(entity, "spriteDirectory", "objects/MaxHelpingHand/noDashRefillSpring") * "/00.png", 24, 0, rot=pi / 2)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashRefillSpringRight, room::Maple.Room) = Ahorn.drawSprite(ctx, get(entity, "spriteDirectory", "objects/MaxHelpingHand/noDashRefillSpring") * "/00.png", -8, 16, rot=-pi / 2)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function GridLayout(spec::GridLayoutSpec)
    gl = GridLayout(; spec.kwargs...)

    for ((rows, cols, side), content) in spec.content
        if content isa GridLayoutSpec
            content = GridLayout(content)
        end
        gl[rows, cols, side] = content
    end

    gl
end

const PosTuple = Tuple{Indexables, Indexables, Side}
to_valid_pos(pos::Tuple{Indexables, Indexables}) = (pos[1], pos[2], Inner())
to_valid_pos(pos::Tuple{Indexables, Indexables, Side}) = pos

function GridLayoutSpec(content::Vector{<:Pair}; kwargs...)
    spec_content = Pair{PosTuple, Any}[Pair{PosTuple, Any}(to_valid_pos(pos), c) for (pos, c) in content]
    GridLayoutSpec(spec_content, Dict{Symbol, Any}(kwargs))
end


text:
using Distributed

addprocs(12, exeflags="--project")   # create worker processes with current project activated
@everywhere include("parallel_setup.jl")

# simulated annealing parameters
T₀ = 0.5
N = 68
Ns = 36
Nt = 10
tol = 1.0

# bounds and step length
ub = repeat([1.0, repeat([0.5, 0.5, 1.0], 3)...], 6)
lb = repeat([0.01, repeat([0.0, 0.1, 0.01], 3)...], 6) # constrain lb of activation to 0.01
append!(ub, [20_000, 2000, 200_000, 100_000, 20_000, 1000, 200_000, 100_000])
append!(lb, zeros(8))
v = ub .- lb

# initial guess
randx(lb, ub) = [rand(lb:0.001:ub) for (lb, ub) in zip(lb, ub)]
x₀ = [rand(lb:0.001:ub) for (lb, ub) in zip(lb, ub)] #=vcat(map(x -> inputs.activation_parameters[x], [:he, :ke, :ae, :hf, :kf, :af])...) =#
x₀ = [vcat(map(x -> inputs.activation_parameters[x], [:he, :ke, :ae, :hf, :kf, :af])...)..., p.k1, p.k2, p.k3, p.k4, p.k5, p.k6, p.k7, p.k8]
f₀ = objective(x₀)

# initial structs
current = State(f=f₀, x=x₀, v=v, T=T₀)
result = Result(fopt=f₀, xopt=x₀)
options = Options(func=objective, N=N, Ns=Ns, Nt=Nt, lb=lb, ub=ub, tol=tol, print_status=true)
@time sa!(current, result, options)

open("optimisations/matching/results.csv", "a") do io
    writedlm(io, [inputs.initial_conditions[:vcmx] result.fopt result.xopt...], ',')
end

# # optimisation loop
# submax_velocity = true
# # set initial speed
# @everywhere VCMX = 9.5704
# while submax_velocity
#     # set velocity and reinitialise initial conditions and prob
#     @everywhere begin
#         VCMX += 0.1
#         inputs.initial_conditions[:vcmx] += 0.1
#         p, u₀ = set_values(inputs)
#         prob = ODEProblem(eom, u₀, tspan, p)
#     end

#     println("Rerunning optimisation with vcmx = $VCMX")

#     # (re)initialise sa
#     x₀ = result.xopt    # use previous optimal x as starting point for next optimisation
#     f₀ = objective(x₀)
#     current = State(f = f₀, x = x₀, v = v, T = T₀)
#     result = Result(fopt = f₀, xopt = x₀)
#     options = Options(func = objective, N = N, Ns = Ns, Nt = Nt, lb = lb, ub = ub, tol = tol, print_status = true)

#     # run sa
#     @time sa!(current, result, options)

#     # write results to file
# open("optimisations/matching/results.csv", "a") do io
# writedlm(io, [VCMX result.fopt result.xopt...], ',')
# end
#     submax_velocity = false

#     # check if conditions met to be able to run at VCMX
#     popt, u₀opt = updateParameters(p, result.xopt, u₀)
#     sol1, sol2 = simulate(result.xopt)
#     # animate_model(sol)
#     # submax_velocity = abs(swing_time(sol) - TSW) < 0.001 && abs(step_velocity(sol) - VCMX) < 0.01
# end


text:
"""
    assigndist_binomial(n, p; trunc_lower = -Inf, trunc_upper = Inf)

Assign parameters to a binomial distribution with `n` trials and probability `p`
of success in individual trials, optionally truncating the distribution.
"""
function assigndist_binomial(n, p; trunc_lower = -Inf, trunc_upper = Inf)
    truncated(Binomial(n, p), trunc_lower, trunc_upper)
end

export assigndist_binomial


text:
"""
    UniGrid(m::Int)

where `m` is the number of points on each dimension
Adaptive uniform grid based on [1]

[1] Moreno-Muñoz, P., Artés-Rodríguez, A. & Álvarez, M. A. Continual Multi-task Gaussian Processes. (2019).
"""
struct UniGrid <: OnIPSA
    m::Int # Number of points per dimension
end

Base.show(io::IO, Z::UniGrid) = print(io, "Uniform grid with side length $(Z.m).")

function initZ(
    ::AbstractRNG,
    alg::UniGrid,
    X::Union{AbstractVector{<:Real},AbstractVector{<:AbstractVector{<:Real}}};
    kwargs...,
)
    ndim = length(first(X)) # Take the dimensionality
    bounds = [extrema(x -> getindex(x, i), X) for i in 1:ndim]
    Z = map(bounds) do lims
        LinRange(lims..., alg.m)
    end
    return Z
end

function updateZ!(
    ::AbstractRNG, Z::AbstractVector, alg::UniGrid, X::AbstractVector; kwargs...
)
    ndim = length(Z)
    new_bounds = [extrema(x -> getindex(x, i), X) for i in 1:ndim]
    map!(Z, Z, new_bounds) do Z_d, new_b
        x_start = min(Z_d.start, new_b[1]) # Find the new limits
        x_stop = max(Z_d.stop, new_b[2])
        return LinRange(x_start, x_stop, alg.m) # readapt bounds
    end
    return Z
end


text:
export lanbpro

inclue("lanczos.jl")
# inclue("cg.jl")
# inclue("cgls.jl")
# inclue("lsqr.jl")


text:
# 1D Grid examples
# ===============
using ExtendableGrids

# ## Interval from vector
#
function interval_from_vector()
    X=collect(0:0.05:1)
    grid=simplexgrid(X)
end
# ![](interval_from_vector.svg)

#
# ##Interval with local refinement
#
function interval_localref()

    XLeft=geomspace(0.0,0.5,0.1, 0.01)
    XRight=geomspace(0.5,1.0,0.01,0.1)
    X=glue(XLeft, XRight)
    grid=simplexgrid(X)
end
# ![](interval_localref.svg)

#
# ## Interval with  multiple regions
#
function interval_multiregion()

    X=collect(0:0.05:1)
    grid=simplexgrid(X)
    cellmask!(grid,[0.0],[0.5],3)
    bfacemask!(grid,[0.5], [0.5],4)
    grid
end
# ![](interval_multiregion.svg)
#
# ## Multiple regions and subgrid
#
function interval_subgrid()
    X=collect(0:0.01:1)
    grid=simplexgrid(X)
    bfacemask!(grid,[0.5],[0.5],3)
    cellmask!(grid,[0.0],[0.25],2)
    cellmask!(grid,[0.20],[0.5],3)
    subgrid(grid,[2,3])
end
# ![](interval_subgrid.svg)


text:
"""
    Sorting networks for 1,2,..16 and 17,18..25 values are given.
    These sorting networks should minimize conditional swaps.
    The first sixteen are known to minimize conditional swaps.

    Functions are a sequence of internally parallelizable blocks.
    Functions are written using empty lines as block delmiters.

    `swapsort(x1::T, .., xn::T)::NTuple{n,T}`
    `swapsort(tup::NTuple{n,T})::NTuple{n,T}`
"""
module SortingNetworks

export swapsort, ExchangeSort

struct ExchangeSortAlg  <: Base.Sort.Algorithm end
const ExchangeSort  = ExchangeSortAlg()

include("swapsort.jl")
include("swapsort_17to25.jl")

end # module


text:
# Definitions of physical unit types for dispatch. At the moment this is a
# "dummy" set of definitions that just sets SI units to 1, but this may be
# swapped out for actual units if I have time.

using BasicUnits

const Mass = Real
const Distance = Real
const Radius = Distance
const Temperature = Real
const MassLuminosity = Real
const Luminosity = Real
const Flux = Real
const Pressure = Real
const Opacity = Real
const Acceleration = Real
const Density = Real
const Gravity = Acceleration
const Dimensionless = Real
const OpticalDepth = Dimensionless

#= definitions unused for now
using SIUnits
using SIUnits.ShortUnits
import SIUnits.unit, SIUnits.SIQuantity

# const Unit{T} = SIQuantity{T,...}
const Mass{T} = SIQuantity{T,0,1,0,0,0,0,0,0,0}
const Distance{T} = SIQuantity{T,1,0,0,0,0,0,0,0,0}
const Radius = Distance
const Temperature{T} = SIQuantity{T,0,0,0,0,1,0,0,0,0}
const MassLuminosity{T} = SIQuantity{T,2,0,-3,0,0,0,0,0,0}
const Luminosity{T} = SIQuantity{T,2,1,-3,0,0,0,0,0,0}
const Flux{T} = SIQuantity{T,0,1,-3,0,0,0,0,0,0}
const Pressure{T} = SIQuantity{T,-1,1,-2,0,0,0,0,0,0}
const Opacity{T} = SIQuantity{T,2,-1,0,0,0,0,0,0,0}
const Acceleration{T} = SIQuantity{T,1,0,-2,0,0,0,0,0,0}
const Density{T} = SIQuantity{T,-3,1,0,0,0,0,0,0,0}
const Gravity = Acceleration
const Dimensionless = Real
const OpticalDepth = Dimensionless
=#


text:
using AFMEphys
using Test

testdat = makeFileList(pwd() * "\\example", ".asc")
println(pwd())

@testset "AFMEphys.jl" begin
    @test testdat[1] == pwd() * "\\example\\test.asc"
    @test length(testdat) == 1
end


text:
abstract type AbstractReservesFormulation <: AbstractServiceFormulation end
struct RangeReserve <: AbstractReservesFormulation end
############################### Reserve Variables` #########################################
"""
This function add the variables for reserves to the model
"""
function activeservice_variables!(
    psi_container::PSIContainer,
    service::SR,
    contributing_devices::Vector{<:PSY.Device},
) where {SR <: PSY.Reserve}
    add_variable(
        psi_container,
        contributing_devices,
        variable_name(PSY.get_name(service), SR),
        false;
        lb_value = d -> 0,
    )
    return
end

################################## Reserve Requirement Constraint ##########################
# This function can be generalized later for any constraint of type Sum(req_var) >= requirement,
# it will only need to be specific to the names and get forecast string.
function service_requirement_constraint!(
    psi_container::PSIContainer,
    service::SR,
    model::ServiceModel{SR, RangeReserve},
) where {SR <: PSY.Reserve}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)
    initial_time = model_initial_time(psi_container)
    @debug initial_time
    time_steps = model_time_steps(psi_container)
    name = PSY.get_name(service)
    constraint = get_constraint(psi_container, constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(psi_container, variable_name(name, SR))

    if use_forecast_data
        ts_vector = TS.values(PSY.get_data(PSY.get_forecast(
            PSY.Deterministic,
            service,
            initial_time,
            "get_requirement",
            length(time_steps),
        )))
    else
        ts_vector = ones(time_steps[end])
    end

    requirement = PSY.get_requirement(service)
    if parameters
        param = get_parameter_array(
            psi_container,
            UpdateRef{SR}(SERVICE_REQUIREMENT, "get_requirement"),
        )
        for t in time_steps
            param[name, t] =
                PJ.add_parameter(psi_container.JuMPmodel, ts_vector[t] * requirement)
            constraint[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                sum(reserve_variable[:, t]) >= param[name, t]
            )
        end
    else
        for t in time_steps
            constraint[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                sum(reserve_variable[:, t]) >= ts_vector[t] * requirement
            )
        end
    end
    return
end

function modify_device_model!(
    devices_template::Dict{Symbol, DeviceModel},
    service_model::ServiceModel{<:PSY.Reserve, RangeReserve},
    contributing_devices::Vector{<:PSY.Device},
)
    device_types = unique(typeof.(contributing_devices))
    for dt in device_types
        for (device_model_name, device_model) in devices_template
            # add message here when it exists
            device_model.device_type != dt && continue
            service_model in device_model.services && continue
            push!(device_model.services, service_model)
        end
    end

    return
end

function include_service!(constraint_data::DeviceTimeSeries, services, SM::ServiceModel)
    range_data = constraint_data.range
    isnothing(range_data) && return
    include_service!(range_data, services, SM)
    return
end

function include_service!(
    constraint_data::DeviceRange,
    services,
    ::ServiceModel{SR, <:AbstractReservesFormulation},
) where {SR <: PSY.Reserve{PSY.ReserveUp}}
    for (ix, service) in enumerate(services)
        push!(
            constraint_data.additional_terms_ub,
            constraint_name(PSY.get_name(service), SR),
        )
    end
    return
end

function include_service!(
    constraint_data::DeviceRange,
    services,
    ::ServiceModel{SR, <:AbstractReservesFormulation},
) where {SR <: PSY.Reserve{PSY.ReserveDown}}
    for (ix, service) in enumerate(services)
        push!(
            constraint_data.additional_terms_lb,
            constraint_name(PSY.get_name(service), SR),
        )
    end
    return
end

function add_device_services!(
    constraint_data::RangeConstraintsData,
    device::D,
    model::DeviceModel,
) where {D <: PSY.Device}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.service_type)
            services =
                (s for s in PSY.get_services(device) if isa(s, service_model.service_type))
            @assert !isempty(services)
            include_service!(constraint_data, services, service_model)
        end
    end
    return
end

function add_device_services!(
    constraint_data_in::RangeConstraintsData,
    constraint_data_out::RangeConstraintsData,
    device::D,
    model::DeviceModel{D, <:AbstractStorageFormulation},
) where {D <: PSY.Storage}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.service_type)
            services =
                (s for s in PSY.get_services(device) if isa(s, service_model.service_type))
            @assert !isempty(services)
            if service_model.service_type <: PSY.Reserve{PSY.ReserveDown}
                for service in services
                    push!(
                        constraint_data_in.additional_terms_ub,
                        constraint_name(PSY.get_name(service), service_model.service_type),
                    )
                end
            elseif service_model.service_type <: PSY.Reserve{PSY.ReserveUp}
                for service in services
                    push!(
                        constraint_data_out.additional_terms_ub,
                        constraint_name(PSY.get_name(service), service_model.service_type),
                    )
                end
            end
        end
    end
    return
end


text:
using LoopVectorization, OffsetArrays, Test

@testset "copy" begin

  function copyavx1!(x, y)
    @turbo for i ∈ eachindex(x)
      x[i] = y[i]
    end
  end
  function copy_avx1!(x, y)
    @_avx for i ∈ eachindex(x)
      x[i] = y[i]
    end
  end
  function copyavx2!(x, y)
    @turbo for i ∈ eachindex(x)
      yᵢ = y[i]
      x[i] = yᵢ
    end
  end
  function copy_avx2!(x, y)
    @_avx for i ∈ eachindex(x)
      yᵢ = y[i]
      x[i] = yᵢ
    end
  end
  function offset_copy!(A, B)
    @inbounds for i = 1:size(A, 1), j = 1:size(B, 2)
      A[i, j+2] = B[i, j]
    end
  end
  function offset_copyavx1!(A, B)
    @turbo for i = 1:size(A, 1), j = 1:size(B, 2)
      A[i, j+2] = B[i, j]
    end
  end
  function offset_copy_avx1!(A, B)
    @_avx for i = 1:size(A, 1), j = 1:size(B, 2)
      @inbounds A[i, j+2] = B[i, j]
    end
  end
  function offset_copyavx2!(A, B)
    @turbo for i = 1:size(A, 1), j = 1:size(B, 2)
      Bᵢⱼ = B[i, j]
      A[i, j+2] = Bᵢⱼ
    end
  end
  function offset_copy_avx2!(A, B)
    @_avx for i = 1:size(A, 1), j = 1:size(B, 2)
      Bᵢⱼ = B[i, j]
      A[i, j+2] = Bᵢⱼ
    end
  end
  function make2point3avx!(x)
    a = 1.742416161578685
    b = 1.5
    @turbo for i ∈ eachindex(x)
      x[i] = a^b
    end
  end
  function make2point3_avx!(x)
    @_avx for i ∈ eachindex(x)
      x[i] = 2.3
    end
  end
  function make23avx!(x)
    @turbo for i ∈ eachindex(x)
      @inbounds x[i] = 23
    end
  end
  function make23_avx!(x)
    @_avx for i ∈ eachindex(x)
      x[i] = 23
    end
  end
  function myfillavx!(x, a)
    @turbo for i ∈ eachindex(x)
      x[i] = a
    end
  end
  function myfill_avx!(x, a)
    @_avx for i ∈ eachindex(x)
      x[i] = a
    end
  end
  function issue_256!(A)
    @turbo for i = 1:size(A, 1)
      A[i, (x = 1,).x] = 0
    end
    A
  end

  function reversecopy1!(B, A)
    for i in eachindex(A)
      B[i] = A[11-i]
    end
    B
  end
  function reversecopy1avx!(B, A)
    @turbo for i in eachindex(A)
      B[i] = A[11-i]
    end
    B
  end
  function reversecopy2!(B, A)
    for i in eachindex(B)
      B[i] = A[-i]
    end
    B
  end
  function reversecopy2avx!(B, A)
    @turbo for i in eachindex(B)
      B[i] = A[-i]
    end
    B
  end
  function reversecopy3avx!(B, A)
    n = length(A)
    @turbo for i = 1:n
      B[i] = A[n+1-i]
    end
    B
  end
  function copy3!(B, A)
    @assert (length(B) ≥ 3) && (length(A) ≥ 3)
    @turbo for i = 1:3
      B[i] = A[i]
    end
    B
  end
  function copyselfdot!(s, x)
    m = zero(eltype(x))
    @turbo for i ∈ 1:2
      sᵢ = x[i]
      s[i] = sᵢ
      m += sᵢ * sᵢ
    end
    m
  end
  function scattercopyavx!(H, a, j)
    @turbo for i ∈ eachindex(j), k ∈ eachindex(a)
      H[j[i], k] = a[k]
    end
    H
  end
  function scattercopy!(H, a, j)
    @inbounds for i ∈ eachindex(j), k ∈ eachindex(a)
      H[j[i], k] = a[k]
    end
    H
  end

  function issue279!(output, input, idx)
    @turbo for j in axes(output, 2), i in axes(output, 1)
      output[i, j, idx] = input[1, 1, i, j, idx]
    end
    output
  end

  function issue279_plain!(output, input, idx)
    for j in axes(output, 2), i in axes(output, 1)
      output[i, j, idx] = input[1, 1, i, j, idx]
    end
    output
  end
  function issue280!(dest, src)
    @turbo for i in indices((dest, src), (2, 2))
      dest[1, i] = src[2, i]
      dest[2, i] = src[1, i]
    end
  end

  function fill_with_3rd!(dest, src)
    p = 1
    @turbo for i ∈ eachindex(dest)
      dest[i] = src[3*p]
    end
  end
  function fill_with_r1c3!(dest, src::AbstractMatrix)
    p = 1
    @turbo for i ∈ eachindex(dest)
      dest[i] = src[p, 3*p]
    end
  end
  function collect_turbo(N, ::Type{T}) where {T}
    l = Vector{T}(undef, N)
    @turbo for i ∈ 1:length(l)
      l[i] = i
    end
    return l
  end
  for T ∈ (Float32, Float64, Int32, Int64)
    @show T, @__LINE__
    R = T <: Integer ? (-T(100):T(100)) : T
    x = rand(R, 237)
    q1 = similar(x)
    q2 = similar(x)

    fill!(q2, -999999)
    copyavx1!(q2, x)
    @test x == q2
    fill!(q2, -999999)
    copy_avx1!(q2, x)
    @test x == q2
    fill!(q2, -999999)
    copyavx2!(q2, x)
    @test x == q2
    fill!(q2, -999999)
    copy_avx2!(q2, x)
    @test x == q2
    fill!(q2, -999999)
    @turbo q2 .= x
    @test x == q2
    @test all(iszero, issue_256!(reshape(x, (length(x), 1))))

    B = rand(R, 79, 83)
    A1 = zeros(T, 79, 85)
    A2 = zeros(T, 79, 85)
    offset_copy!(A1, B)
    fill!(A2, 0)
    offset_copyavx1!(A2, B)
    @test A1 == A2
    fill!(A2, 0)
    offset_copyavx2!(A2, B)
    @test A1 == A2
    fill!(A2, 0)
    offset_copy_avx1!(A2, B)
    @test A1 == A2
    fill!(A2, 0)
    offset_copy_avx2!(A2, B)
    @test A1 == A2
    fill_with_3rd!(A1, B)
    @test all(==(B[3]), A1)
    fill_with_3rd!(A1, eachindex(B))
    @test all(==(3), A1)

    fill_with_r1c3!(A1, B)
    @test all(==(B[1, 3]), A1)
    fill_with_r1c3!(A1, LinearIndices(B))# check args fails
    @test all(==(LinearIndices(B)[1, 3]), A1)

    a = rand(R)
    myfillavx!(x, a)
    fill!(q2, a)
    @test x == q2
    a = rand(R)
    myfill_avx!(x, a)
    fill!(q2, a)
    @test x == q2
    a = rand(R)
    myfill_avx!(x, a)
    fill!(q2, a)
    @test x == q2
    a = rand(R)
    myfillavx!(x, a)
    fill!(q2, a)
    @test x == q2
    q2 .= 23
    fill!(q1, -99999)
    make23_avx!(q1)
    @test q2 == q1
    fill!(q1, -99999)
    make23avx!(q1)
    @test q2 == q1
    if T <: Union{Float32,Float64}
      make2point3avx!(x)
      fill!(q2, 2.3)
      @test x == q2
      fill!(x, -999999)
      make2point3_avx!(x)
      @test x == q2
    end
    a = rand(R)
    @turbo x .= a
    fill!(q2, a)
    @test x == q2
    a = rand(R)
    @turbo x .= a
    fill!(q2, a)
    @test x == q2

    @test reversecopy1!(zeros(T, 10), collect(1:10)) ==
          reversecopy1avx!(zeros(T, 10), collect(1:10))
    @test reversecopy2!(zeros(T, 10), OffsetArray(collect(1:10), -10:-1)) ==
          reversecopy2avx!(zeros(T, 10), OffsetArray(collect(1:10), -10:-1))
    @test 177:-1:1 == reversecopy3avx!(zeros(T, 177), collect(1:177))

    x = rand(R, 3)
    y = similar(x)
    @test copy3!(y, x) == x
    fill!(y, 0)
    @test copyselfdot!(y, x) ≈ x[1]^2 + x[2]^2
    @test view(x, 1:2) == view(y, 1:2)

    H0 = zeros(T, 10, 10)
    H1 = zeros(T, 10, 10)
    j = [1, 2, 5, 8]
    a = rand(R, 10)
    @test scattercopyavx!(H0, a, j) == scattercopy!(H1, a, j)

    input = rand(R, 2, 2, 5, 5, 1)
    output = Array{T}(undef, size(input)[3:end]...)
    output_plain = similar(output)

    @test issue279!(output, input, 1) ≈ issue279_plain!(output_plain, input, 1)

    src = rand(R, 2, 17)
    dest = similar(src)
    issue280!(dest, src)
    @test dest ≈ vcat(view(src, 2, :)', view(src, 1, :)')
    if VERSION ≥ v"1.6"
      src2 = reinterpret(reshape, T, Vector{Tuple{T,T}}(undef, 17))
      src2 .= src
      dest2 = reinterpret(reshape, T, Vector{Tuple{T,T}}(undef, 17))
      issue280!(dest2, src2)
      @test dest2 ≈ vcat(view(src, 2, :)', view(src, 1, :)')
    end

    rng = 10:20
    arr1 = zeros(30)
    arr2 = zeros(30)
    ifirst = 5
    @turbo for i in rng
      arr1[i-ifirst] += 1
    end
    @. arr2[rng-ifirst] += 1
    @test arr1 == arr2

    @test collect_turbo(77, T) == T.(1:77)
  end
end


text:
@testset "1105.filling-bookcase-shelves.jl" begin
    @test min_height_shelves([[1,1],[2,3],[2,3],[1,1],[1,1],[1,1],[1,2]], 4) == 6
    @test min_height_shelves([[1,3],[2,4],[3,2]], 6) == 4
    @test min_height_shelves([[17, 20], [9, 15], [4, 9], [20, 20], [11, 9], [15, 5], [4, 15], [7, 16], [8, 1], [3, 19], [11, 1], [17, 15], [20, 16], [7, 20], [11, 2], [5, 18], [1, 20], [16, 15], [10, 7], [16, 4]], 20) == 196
end


text:
# Rank-based correlations
#
# - Spearman's correlation
# - Kendall's correlation
#

#######################################
#
#   Spearman correlation
#
#######################################

"""
    corspearman(x, y=x)

Compute Spearman's rank correlation coefficient. If `x` and `y` are vectors, the
output is a float, otherwise it's a matrix corresponding to the pairwise correlations
of the columns of `x` and `y`.
"""
function corspearman(x::RealVector, y::RealVector)
    n = length(x)
    n == length(y) || throw(DimensionMismatch("vectors must have same length"))
    (any(isnan, x) || any(isnan, y)) && return NaN
    return cor(tiedrank(x), tiedrank(y))
end

function corspearman(X::RealMatrix, y::RealVector)
    size(X, 1) == length(y) ||
        throw(DimensionMismatch("X and y have inconsistent dimensions"))
    n = size(X, 2)
    C = Matrix{Float64}(I, n, 1)
    any(isnan, y) && return fill!(C, NaN)
    yrank = tiedrank(y)
    for j = 1:n
        Xj = view(X, :, j)
        if any(isnan, Xj)
            C[j,1] = NaN
        else
            Xjrank = tiedrank(Xj)
            C[j,1] = cor(Xjrank, yrank)
        end
    end
    return C
end

function corspearman(x::RealVector, Y::RealMatrix)
    size(Y, 1) == length(x) ||
        throw(DimensionMismatch("x and Y have inconsistent dimensions"))
    n = size(Y, 2)
    C = Matrix{Float64}(I, 1, n)
    any(isnan, x) && return fill!(C, NaN)
    xrank = tiedrank(x)
    for j = 1:n
        Yj = view(Y, :, j)
        if any(isnan, Yj)
            C[1,j] = NaN
        else
            Yjrank = tiedrank(Yj)
            C[1,j] = cor(xrank, Yjrank)
        end
    end
    return C
end

function corspearman(X::RealMatrix)
    n = size(X, 2)
    C = Matrix{Float64}(I, n, n)
    anynan = Vector{Bool}(undef, n)
    for j = 1:n
        Xj = view(X, :, j)
        anynan[j] = any(isnan, Xj)
        if anynan[j]
            C[:,j] .= NaN
            C[j,:] .= NaN
            C[j,j] = 1
            continue
        end
        Xjrank = tiedrank(Xj)
        for i = 1:(j-1)
            Xi = view(X, :, i)
            if anynan[i]
                C[i,j] = C[j,i] = NaN
            else
                Xirank = tiedrank(Xi)
                C[i,j] = C[j,i] = cor(Xjrank, Xirank)
            end
        end
    end
    return C
end

function corspearman(X::RealMatrix, Y::RealMatrix)
    size(X, 1) == size(Y, 1) ||
        throw(ArgumentError("number of rows in each array must match"))
    nr = size(X, 2)
    nc = size(Y, 2)
    C = Matrix{Float64}(undef, nr, nc)
    for j = 1:nr
        Xj = view(X, :, j)
        if any(isnan, Xj)
            C[j,:] .= NaN
            continue
        end
        Xjrank = tiedrank(Xj)
        for i = 1:nc
            Yi = view(Y, :, i)
            if any(isnan, Yi)
                C[j,i] = NaN
            else
                Yirank = tiedrank(Yi)
                C[j,i] = cor(Xjrank, Yirank)
            end
        end
    end
    return C
end

#######################################
#
#   Kendall correlation
#
#######################################

# Knight, William R. “A Computer Method for Calculating Kendall's Tau with Ungrouped Data.”
# Journal of the American Statistical Association, vol. 61, no. 314, 1966, pp. 436–439.
# JSTOR, www.jstor.org/stable/2282833.
function corkendall!(x::RealVector, y::RealVector, permx::AbstractVector{<:Integer}=sortperm(x))
    if any(isnan, x) || any(isnan, y) return NaN end
    n = length(x)
    if n != length(y) error("Vectors must have same length") end

    # Initial sorting
    permute!(x, permx)
    permute!(y, permx)

    # Use widen to avoid overflows on both 32bit and 64bit
    npairs = div(widen(n) * (n - 1), 2)
    ntiesx = ndoubleties = nswaps = widen(0)
    k = 0

    @inbounds for i = 2:n
        if x[i - 1] == x[i]
            k += 1
        elseif k > 0
            # Sort the corresponding chunk of y, so the rows of hcat(x,y) are
            # sorted first on x, then (where x values are tied) on y. Hence
            # double ties can be counted by calling countties.
            sort!(view(y, (i - k - 1):(i - 1)))
            ntiesx += div(widen(k) * (k + 1), 2) # Must use wide integers here
            ndoubleties += countties(y,  i - k - 1, i - 1)
            k = 0
        end
    end
    if k > 0
        sort!(view(y, (n - k):n))
        ntiesx += div(widen(k) * (k + 1), 2)
        ndoubleties += countties(y, n - k, n)
    end

    nswaps = merge_sort!(y, 1, n)
    ntiesy = countties(y, 1, n)

    # Calls to float below prevent possible overflow errors when
    # length(x) exceeds 77_936 (32 bit) or 5_107_605_667 (64 bit)
    (npairs + ndoubleties - ntiesx - ntiesy - 2 * nswaps) /
     sqrt(float(npairs - ntiesx) * float(npairs - ntiesy))
end

"""
    corkendall(x, y=x)

Compute Kendall's rank correlation coefficient, τ. `x` and `y` must both be either
matrices or vectors.
"""
corkendall(x::RealVector, y::RealVector) = corkendall!(copy(x), copy(y))

function corkendall(X::RealMatrix, y::RealVector)
    permy = sortperm(y)
    return([corkendall!(copy(y), X[:,i], permy) for i in 1:size(X, 2)])
end

function corkendall(x::RealVector, Y::RealMatrix)
    n = size(Y, 2)
    permx = sortperm(x)
    return(reshape([corkendall!(copy(x), Y[:,i], permx) for i in 1:n], 1, n))
end

function corkendall(X::RealMatrix)
    n = size(X, 2)
    C = Matrix{Float64}(I, n, n)
    for j = 2:n
        permx = sortperm(X[:,j])
        for i = 1:j - 1
            C[j,i] = corkendall!(X[:,j], X[:,i], permx)
            C[i,j] = C[j,i]
        end
    end
    return C
end

function corkendall(X::RealMatrix, Y::RealMatrix)
    nr = size(X, 2)
    nc = size(Y, 2)
    C = Matrix{Float64}(undef, nr, nc)
    for j = 1:nr
        permx = sortperm(X[:,j])
        for i = 1:nc
            C[j,i] = corkendall!(X[:,j], Y[:,i], permx)
        end
    end
    return C
end

# Auxilliary functions for Kendall's rank correlation

"""
    countties(x::RealVector, lo::Integer, hi::Integer)

Return the number of ties within `x[lo:hi]`. Assumes `x` is sorted.
"""
function countties(x::AbstractVector, lo::Integer, hi::Integer)
    # Use of widen below prevents possible overflow errors when
    # length(x) exceeds 2^16 (32 bit) or 2^32 (64 bit)
    thistiecount = result = widen(0)
    checkbounds(x, lo:hi)
    @inbounds for i = (lo + 1):hi
        if x[i] == x[i - 1]
            thistiecount += 1
        elseif thistiecount > 0
            result += div(thistiecount * (thistiecount + 1), 2)
            thistiecount = widen(0)
        end
    end

    if thistiecount > 0
        result += div(thistiecount * (thistiecount + 1), 2)
    end
    result
end

# Tests appear to show that a value of 64 is optimal,
# but note that the equivalent constant in base/sort.jl is 20.
const SMALL_THRESHOLD = 64

# merge_sort! copied from Julia Base
# (commit 28330a2fef4d9d149ba0fd3ffa06347b50067647, dated 20 Sep 2020)
"""
    merge_sort!(v::AbstractVector, lo::Integer, hi::Integer, t::AbstractVector=similar(v, 0))

Mutates `v` by sorting elements `x[lo:hi]` using the merge sort algorithm.
This method is a copy-paste-edit of sort! in base/sort.jl, amended to return the bubblesort distance.
"""
function merge_sort!(v::AbstractVector, lo::Integer, hi::Integer, t::AbstractVector=similar(v, 0))
    # Use of widen below prevents possible overflow errors when
    # length(v) exceeds 2^16 (32 bit) or 2^32 (64 bit)
    nswaps = widen(0)
    @inbounds if lo < hi
        hi - lo <= SMALL_THRESHOLD && return insertion_sort!(v, lo, hi)

        m = midpoint(lo, hi)
        (length(t) < m - lo + 1) && resize!(t, m - lo + 1)

        nswaps = merge_sort!(v, lo,  m, t)
        nswaps += merge_sort!(v, m + 1, hi, t)

        i, j = 1, lo
        while j <= m
            t[i] = v[j]
            i += 1
            j += 1
        end

        i, k = 1, lo
        while k < j <= hi
            if v[j] < t[i]
                v[k] = v[j]
                j += 1
                nswaps += m - lo + 1 - (i - 1)
            else
                v[k] = t[i]
                i += 1
            end
            k += 1
        end
        while k < j
            v[k] = t[i]
            k += 1
            i += 1
        end
    end
    return nswaps
end

# insertion_sort! and midpoint copied from Julia Base
# (commit 28330a2fef4d9d149ba0fd3ffa06347b50067647, dated 20 Sep 2020)
midpoint(lo::T, hi::T) where T <: Integer = lo + ((hi - lo) >>> 0x01)
midpoint(lo::Integer, hi::Integer) = midpoint(promote(lo, hi)...)

"""
    insertion_sort!(v::AbstractVector, lo::Integer, hi::Integer)

Mutates `v` by sorting elements `x[lo:hi]` using the insertion sort algorithm.
This method is a copy-paste-edit of sort! in base/sort.jl, amended to return the bubblesort distance.
"""
function insertion_sort!(v::AbstractVector, lo::Integer, hi::Integer)
    if lo == hi return widen(0) end
    nswaps = widen(0)
    @inbounds for i = lo + 1:hi
        j = i
        x = v[i]
        while j > lo
            if x < v[j - 1]
                nswaps += 1
                v[j] = v[j - 1]
                j -= 1
                continue
            end
            break
        end
        v[j] = x
    end
    return nswaps
end


text:
using DataDeps

register(DataDep(
  "mushroom",
  "http://archive.ics.uci.edu/ml/datasets/Mushroom",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/mushroom/agaricus-lepiota.data",
  "e65d082030501a3ebcbcd7c9f7c71aa9d28fdfff463bf4cf4716a3fe13ac360e",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:23,
        categoric_indices=2:23,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@inline function sinh_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return sinh(y)-sinh(x)-(y-x)*cosh(x)
end
@inline function sinh_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (y-x)*sinh(x)
end
@inline function cv_sinh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return sinh(x),cosh(x)
  elseif (xU<=zero(T))
    return line_seg(x,xL,sinh(xL),xU,sinh(xU)),dline_seg(x,xL,sinh(xL),xU,sinh(xU),cosh(x))
  else
    try
      p = newton(xL,xL,zero(T),sinh_env,sinh_envd,xU,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(T),sinh_env,xU,zero(T))
      end
    end
    if (x>p)
      return sinh(x),cosh(x)
    else
      return line_seg(x,p,sinh(p),xU,sinh(xU)),dline_seg(x,p,sinh(p),xU,sinh(xU),cosh(x))
    end
  end
end
@inline function cc_sinh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,sinh(xL),xU,sinh(xU)),dline_seg(x,xL,sinh(xL),xU,sinh(xU),cosh(x))
  elseif (xU<=zero(T))
    return sinh(x),cosh(x)
  else
    try
      p = newton(xU,zero(T),xU,sinh_env,sinh_envd,xL,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(T),xU,sinh_env,xL,zero(T))
      end
    end
    if (x>p)
      return line_seg(x,xL,sinh(xL),p,sinh(p)),dline_seg(x,xL,sinh(xL),p,sinh(p),cosh(x))
    else
      return sinh(x),cosh(x)
    end
  end
end
@inline function sinh(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_sinh(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_sinh(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_sinh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_sinh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_sinh(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_sinh(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, sinh(x.Intv),x.cnst, x.IntvBox, x.xref)
end

@inline function asinh_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (asinh(z)-asinh(x))/(z-x)-one(T)/sqrt(one(T)+x^2)
end
@inline function asinh_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (asinh(z)-asinh(x))/(z-x)^2+x/(x^2+one(T))^(3/2)-one(T)/((z-x)*sqrt(x^2+one(T)))
end
@inline function cv_asinh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,asinh(xL),xU,asinh(xU)),dline_seg(x,xL,asinh(xL),xU,asinh(xU),one(T)/sqrt(x^2+one(T)))
  elseif (xU<=zero(T))
    return asinh(x),one(T)/sqrt(x^2+one(T))
  else
    try
      p = newton(xL/2,xL,zero(x),asinh_env,asinh_envd,xL,xU)
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(x),asinh_env,xL,xU)
      end
    end
    if (x<=p)
      return asinh(x),one(T)/sqrt(x^2+one(T))
    else
      return line_seg(x,p,asinh(p),xU,asinh(xU)),dline_seg(x,p,asinh(p),xU,asinh(xU),one(T)/sqrt(x^2+one(T)))
    end
  end
end
@inline function cc_asinh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p = zero(T)
  if (xL>=zero(T))
    return asinh(x),one(T)/sqrt(x^2+one(T))
  elseif (xU<=zero(T))
    return line_seg(x,xL,asinh(xL),xU,asinh(xU)),dline_seg(x,xL,asinh(xL),xU,asinh(xU),one(T)/sqrt(x^2+one(T)))
  else
    try
      p = newton(xU/2,zero(x),xU,asinh_env,asinh_envd,xL,xU)
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(x),xU,asinh_env,xL,xU)
      end
    end
    if (x<=p)
      return line_seg(x,xL,asinh(xL),p,asinh(p)),dline_seg(x,xL,asinh(xL),p,asinh(p),one(T)/sqrt(x^2+one(T)))
    else
      return asinh(x),one(T)/sqrt(x^2+one(T))
    end
  end
end
@inline function asinh(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_asinh(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_asinh(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_asinh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_asinh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_asinh(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_asinh(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, asinh(x.Intv),x.cnst, x.IntvBox, x.xref)
end

@inline function tanh_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (tanh(y)-tanh(x))/(one(T)-tanh(x)^2)-y+x
end
@inline function tanh_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return -((one(T)-tanh(x)^2)^(-2))*(tanh(y)-tanh(x))
end
@inline function cv_tanh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,tanh(xL),xU,tanh(xU)),dline_seg(x,xL,tanh(xL),xU,tanh(xU),sech(x)^2)
  elseif (xU<=zero(T))
    return tanh(x),sech(x)^2
  else
    try
      p = newton(xL,xL,zero(T),tanh_env,tanh_envd,xU,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(T),tanh_env,xU,zero(T))
      end
    end
    if (x<=p)
      return tanh(x),sech(x)^2
    else
      return line_seg(x,p,tanh(p),xU,tanh(xU)),dline_seg(x,p,tanh(p),xU,tanh(xU),sech(x)^2)
    end
  end
end
@inline function cc_tanh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return tanh(x),sech(x)^2
  elseif (xU<=zero(T))
    return line_seg(x,xL,tanh(xL),xU,tanh(xU)),dline_seg(x,xL,tanh(xL),xU,tanh(xU),sech(x)^2)
  else
    try
      p = newton(xU/2,zero(T),xU,tanh_env,tanh_envd,xL,xU)
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(T),xU,tanh_env,xL,zero(T))
      end
    end
    if (x<=p)
      return line_seg(x,xL,tanh(xL),p,tanh(p)),dline_seg(x,xL,tanh(xL),p,tanh(p),sech(x)^2)
    else
      return tanh(x),sech(x)^2
    end
  end
end
@inline function tanh(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_tanh(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_tanh(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_tanh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_tanh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_tanh(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_tanh(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, tanh(x.Intv),x.cnst, x.IntvBox, x.xref)
end

@inline function atanh_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (x-y)-(one(T)-x^2)*(atan(x)-atan(y))
end
@inline function atanh_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return one(T)+2*x*(atan(x)-atan(y))
end
@inline function cv_atanh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return atanh(x),one(T)/(one(T)-x^2)
  elseif (xU<=zero(T))
    return line_seg(x,xL,atanh(xL),xU,atanh(xU)),dline_seg(x,xL,atanh(xL),xU,atanh(xU),one(T)/(one(T)-x^2))
  else
    try
      p = newton(xU,zero(T),xU,atanh_env,atanh_envd,xL,xU)
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(T),atanh_env,xU,zero(T))
      end
    end
    if (x>p)
      return atanh(x),one(T)/(one(T)-x^2)
    else
      return line_seg(x,p,atanh(p),xU,atanh(xU)),dline_seg(x,p,atanh(p),xU,atanh(xU),one(T)/(one(T)-x^2))
    end
  end
end
@inline function cc_atanh(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,atanh(xL),xU,atanh(xU)),dline_seg(x,xL,atanh(xL),xU,atanh(xU),one(T)/(one(T)-x^2))
  elseif (xU<=zero(T))
    return atanh(x),one(T)/(one(T)-x^2)
  else
    try
      p = newton(xL,xL,zero(T),atanh_env,atanh_envd,xL,xU)
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(T),xU,atanh_env,xL,zero(T))
      end
    end
    if (x>p)
      return line_seg(x,xU,atanh(xU),p,atanh(p)),dline_seg(x,xU,atanh(xU),p,atanh(p),one(T)/(one(T)-x^2))
    else
      return atanh(x),one(T)/(one(T)-x^2)
    end
  end
end
@inline function atanh(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_atanh(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_atanh(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_atanh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_atanh(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_atanh(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_atanh(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, atanh(x.Intv),x.cnst, x.IntvBox, x.xref)
end


text:
using ParallelAccelerator

@acc function blurxy()
    w = 64
    h = 128
    inp = zeros(w,h)
    inp = map(i -> i[1],enumerate(inp))

    blur = map(i -> (inp[max(i[1]-1,1)]+inp[i[1]]+inp[min(i[1]+1,length(inp))])/3,enumerate(inp))
    return blur
end

function main()

    tic()
    blur_out = blurxy()
    time = toq()
    #println("result = ", blur_out)
    println("SELFTIMED ", time)

    #benchmark()

end

function benchmark()
  maxt = 0
  mint = Inf
  totalt = 0
  for i in 1:100
    tic()
    centroids_out = kmeans()
    t = toq()
    println("finished ",i)
    if t > maxt
      maxt = t
    end
    if t < mint
      mint = t
    end
    totalt += t
  end
  println("SELFTIMED AVERAGE ", totalt/100.0)
  println("SELFTIMED MINIMUM", mint)
  println("SELFTIMED MAXIMUM", maxt)
end

main()

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu